Nivel piezométrico de un acuífero confinado (A5)

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 Obtención de la ecuación

Para poder resolver la ecuación anterior y que el resultado sea único, hemos de establecer condiciones iniciales y de contorno.

Las condiciones de contorno serán dos pues aparece la derivada segunda de la posición: [math] h(\rho_0 , t) = h_p \\ h(20,t) = h_0 [/math]

La condición inicial será única pues sólo aparece una derivada de primer orden en tiempo en nuestra ecuación: [math] h(\rho,0) = h_0 [/math]

De este modo la ecuación resultante será:

[math] \left\{\begin{matrix} \frac{\partial h}{\partial t} - \frac{k}{s} (\frac{1}{\rho} \frac{\partial h}{\partial \rho} + \frac{\partial^2 h}{\partial \rho^2})= 0 & \rho \epsilon (\rho_{0}, 2c) & \theta \epsilon (0, 2\Pi ) & t \gt 0 \\ h(\rho_{0},t) = h_{p}; h(20, t)= h_{0} \\ h(\rho, 0) = h_{0} \end{matrix}\right. [/math]

4 Resolución por método de diferencias finitas

Suponiendo [math]ρ_0 =0.5m[/math], [math]h_0 =45m[/math], [math]hp =34[/math], D=10−2, [math](ρ,0)=h_0[/math]; y tomando ∆ρ = 0.1 y ∆t = ∆ρ. Vamos a resolver la ecuación anterior por el método de diferencias finitas.

Primero tenemos que justificar la interpretación de ∆ρ = 0.1:

Empezamos calculando el gradiente [math]\rho=\frac{\partial\rho}{\partial \rho}\ Finalmente, sustituyendo todo lo anterior en la ecuación, llegamos a: \ltmath\gt \left\{\begin{matrix} \frac{\partial h}{\partial t} - 10^{2} (0.1 \frac{\partial h}{\partial \rho} + \frac{\partial^2 h}{\partial \rho^2})= 0 & \rho \epsilon (0.5, 20) & \theta \epsilon (0, 2\Pi ) & t \gt 0 \\ h(0.5,t) = 34; h(20, t)= 45 \\ h(\rho, 0) = 45 \end{matrix}\right. [/math]

% Escribe aquí tu código

a = 1 + 2;
disp(a);

https://wiki.geogebra.org/es/C%C3%B3digo_LaTeX_para_las_f%C3%B3rmulas_m%C3%A1s_comunes