Nivel piezométrico en un acuifero confinado. Grupo 10-B

Trabajo realizado por estudiantes
Título	Nivel piezométrico en un acuifero confinado. Grupo 10-B
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2016-17
Autores	Pablo Merayo Vidal(1909) María Guadalupe Aranda Sánchez (2049) Marcos Buitrago Peña (656) Dominique Manuela Cazar Espín (1850) Ignacio del Río Pérez (1693) Javier Martín Salgado(1899) Sara Fernández Gago (1663)
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

1 Introducción

Consideramos un acuífero confinado entre dos capas de terreno impermeables horizontales. Llamaremos h(x, y) al nivel piezométrico del acuífero, definido por la altura que alcanzaría el agua si hiciéramos un sondeo en el punto (x, y). Obviamente este valor depende de la presión a la que se encuentra el agua confinada en el acuífero. Supondremos que el acuífero es un medio poroso saturado de agua que ocupa una región infinita y está en equilibrio de manera que [math]h(x, y) = h_0[/math] constante para todo (x, y) ∈ [math]IR^{2}[/math].

Sobre el acuifero se construye un pozo de sección circular y radio[math] ρ_0[/math] para extraer agua. La presencia del pozo hace que el nivel piezométrico cambie. Dada la simetría del problema podemos suponer que h(x, y) sólo va a depender de la distancia al pozo. Si tomamos coordenadas cilíndricas de forma que el eje OZ coincida con el eje de simetría del pozo, entonces h = h(ρ) donde [math]ρ = \sqrt{x^{2} + y^{2}}[/math] . Trabajaremos por tanto en coordenadas polares en el plano (ρ, θ).

Las ecuaciones que permiten conocer h(ρ) son la ecuación de conservación de la masa:

[math] S = \frac{∂h}{∂t} + div q = 0 [/math]

junto con la ley de Darcy

[math] q = −K∇h [/math]

que es una ley experimental que modela el flujo de agua en un medio poroso. La ley de Darcy establece que el flujo de agua q a través de un medio poroso es proporcional a la diferencia de presión, que a su vez se puede escribir en términos del gradiente del nivel piezométrico en cada punto. La constante K se deduce experimentalmente para cada material y se conoce como la conductividad hidráulica o permeabilidad. Cuanto mayor es la constante K mayor es el flujo de agua provocado por un cambio de presión. La ley de Darcy proporciona una buena aproximación del comportamiento del agua en un medio poroso siempre que este sea homogéneo e isótropo.

La constante S en la ley de conservación de la masa se conoce como almacenamiento específico y se interpreta como la cantidad de agua que libera el acuífero al descender el nivel piezométrico en una unidad, por unidad de volumen.

Combinando las ecuaciones de conservación de la masa con la ley de Darcy, obtenemos la ecuación:

[math] \frac{∂h}{∂t} - D∆h = 0, ρ \gt ρ_0, θ ∈ (0, 2π), t \gt 0, (1) [/math]

donde D = K/S es la difusividad hidráulica que supondremos constante.

2 Deducción de la fórmula

En este apartado buscamos demostrar la ecuación (1) como combinación de la ecuación de conservación de masas y la Ley de Darcy. Por un lado sabemos que el laplaciaco equivale a la divergencia del gradiente:

[math] (\Delta f = \nabla \cdot \nabla f = div \left (\nabla f) \right) [/math]

Por otro lado, sabemos que en coordenadas polares:

[math]x = ρ cos \theta \\ y = ρ sen \theta [/math]

[math]\bar{g_ρ } = cos \theta \bar{i} + sin \theta \bar{j} \\ \bar{g_\theta } = - ρ sin \theta \bar{i} + ρ cos \theta \bar{j} [/math]

[math] G = \left( \begin{array}{cc} 1& 0 \\ 0 & ρ ^{2} \\ \end{array} \right); \text {donde la raiz del determinante es} \sqrt{g} = ρ [/math]

Como indica el enunciado, al sustituir la Ley de Darcy en la ecuación de conservación de la masa obtenemos y dividir entre S, teniendo en cuenta que el divergente del gradiente es el laplaciano como hemos indicado al principio de este apartado, obtenemos la ecuación (1): [math] \frac{∂h}{∂t} - \frac{k}{S}∆h = 0 [/math]

En esta ecuación queda por desarrollar el Laplaciano de h. Partiremos de la función h(ρ,θ), de la cual calcularemos primero el gradiente:

[math] \nabla h = \frac{\partial h}{\partial \rho}\hat{\rho} + \frac{\partial h}{\partial \rho}\hat{\theta} = \frac{\partial h}{\partial \rho}\bar{\rho} + \frac{1}{\rho^{2}} \frac{\partial h}{\partial \rho}\bar{\theta} [/math]

teniendo en cuenta que [math]\bar{g^{i} } = \frac{\bar{g_i } }{ | \bar{g_i } | ^{2} } [/math]

Ahora hacemos el divergente del gradiente anterior y obtenemos el Laplaciano en coordenadas polares:

[math]
 \Delta h = div(∇h) 
= {1 \over \rho} {\partial \over \partial \rho}
  \left( \rho {\partial h \over \partial \rho} \right) 
+ {1 \over \rho^2} {\partial^2 h \over \partial \theta^2}
= {1 \over \rho} {\partial h \over \partial \rho} 
+ {\partial^2 h \over \partial \rho^2}
+ {1 \over \rho^2} {\partial^2 h \over \partial \theta^2}
[/math]

Si consideramos que, debido a la simetría del problema, h sólo depende de su distancia al pozo, entonces obtenemos la función radial h(ρ), que nos permite simplificar el Laplaciano del modo siguiente:

[math]
 \Delta h = div(∇h) 
= {1 \over \rho} {\partial h \over \partial \rho} 
+ {\partial^2 h \over \partial \rho^2}
[/math]

3 Resolución numérica por el método de diferencias finitas

Utilizamos el método de diferencias finitas con el objetivo de aproximar la ecuación en derivadas parciales a un problema discreto cuya resolución conduce a un sistema lineal de ecuaciones en cada paso de tiempo. Este método permite tratar problemas con fronteras irregulares sin tener que definir ecuaciones especiales para cada una de estas. Además, como la ubicación nodal no se ajustar necesariamente a una malla regular, es posible ubicar con precisión las extracciones y discontinuidades de entorno y del propio acuífero. Los resultados de estas simulaciones permitieron demostrar la utilidad de este método para tratar geometrías complejas y distribuciones irregulares de extracciones y definir el patrón natural de flujo y las alteraciones ocasionadas por las extracciones, así como su influencia en los niveles piezométricos.

Suponemos para nuestro problema un [math]\rho_{0}= 0.5[/math] m, [math]h_{0}= 45[/math] m, [math]h_{p}= 34[/math], [math] D = 0.01[/math], [math]h(\rho,0)= h_{0}[/math] y el tiempo medido en horas, que en los siguientes apartados realizaremos aproximaciones del problema expuesto con tres métodos distintos: trapecio, euler explícito e implícito .