Diferencia entre revisiones de «Nivel piezométrico en acuífero confinado-Grupo 12»

Revisión del 15:46 19 may 2014

El contenido de este artículo plantea la situación del nivel piezométrico, definido como la altura que alcanzaría el agua al realizar un sondeo en un punto de un acuífero confinado entre dos capas de terreno horizontales impermeables al construirse un pozo, de sección circular y radio [math]\rho _{0} [/math] ,provocando que el nivel piezométrico cambie. Las ecuaciones empleadas para el estudio son la ecuación de conservación de la masa junto con la ley de Darcy, que es una ley experimental que modela el flujo de agua en un medio poroso, estableciendo que el flujo de agua [math] q [/math] es proporcional a la diferencia de presión. La ley de Darcy nos proporciona una buena aproximación del comportamiento del agua en un medio poroso siempre que éste sea homogéneo e isótropo.

1 Interpretación

1.1 Hipótesis iniciales

El análisis que sigue es de aplicación al caso de pozos aislados que penetran totalmente el acuífero por lo que el derrame es horizontal. Los supuestos de partida son los siguientes:

• El acuífero es extensivo, homogéneo e isótropo por lo que se cumple la ley de Darcy.
• El acuífero es un medio poroso saturado de agua, ocupando una región infinita y en equilibrio, por tanto [math] h(x,y)= h_{0} [/math] constante para todo [math](x,y)[/math].
• La pérdida de carga de entrada al pozo es despreciable.
• El acuífero se asienta sobre una superficie horizontal.
• El nivel piezométrico antes del bombeo es (casi) horizontal.
• El caudal extraído del acuífero se produce al mismo tiempo que el descenso del nivel piezométrico.
• La sección transversal que atraviesa el agua permanece constante (el espesor del acuífero es constante).
• El flujo de agua tiene simetría radial.

Nivel piezométrico de acuífero confinado en el que se construye un pozo

Un pozo perfora totalmente un acuífero confinado. La situación inicial al bombeo se corresponde a un régimen hidrostático cuyo límite superior es el nivel piezométrico [math] h_\rho [/math]. Una vez que se comienza el bombeo se produce una depresión en dicho nivel de las zonas próximas al pozo, que disminuye gradualmente al aumentar la distancia a éste. La presencia del pozo hace que el nivel piezométrico cambie, definido este nivel, [math] h(x,y)[/math], por la altura que alcanzaría el agua al realizar un sondeo en el punto [math](x,y)[/math],por tanto va a depender únicamente de la distancia a dicha excavación debido a la simetría del problema. El descenso del nivel piezométrico origina un gradiente hidráulico hacia el pozo que hace que el acuífero suministre el caudal q elevado por la bomba. La velocidad del agua es proporcional a dicho gradiente por lo que la superficie libre adopta una forma acampanada cuya pendiente se incrementa en la proximidad al mismo.

1.2 Ecuaciones y parámetros

Si tomamos coordenadas cilíndricas de manera que OZ coincide con el eje de simetría del pozo, [math]\ h = h_\rho [/math] donde [math] \rho [/math] [math] = \sqrt{x^2\; +\; y^2\;} [/math] . Trabajamos por tanto en coordenadas polares en el plano [math](\rho ,\theta) [/math].

La primera ecuación con la que podemos determinar [math]\ h = h_\rho [/math] es la de conservación de la masa:

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

siendo S el almacenamiento específico, que expresa la masa de agua extraída (o almacenada) por unidad de volumen de acuífero cuando desciende el nivel piezométrico [math] h [/math].

Para conocer completamente [math]\ h = h_\rho [/math] hace falta combinar con una segunda ecuación la de la conservación de la masa , que en este caso es la que explica la ley de Darcy:

:[math] q = - K ·\nabla h [/math]

siendo K la permeabilidad o conductividad hidraulica y deducida para cada material experimentalmente. El flujo de agua que provoca un cambio de presión será mayor, cuanto mayor sea K.

Combinando ambas ecuaciones tenemos como resultado:

[math] \frac{ \partial h }{ \partial t } - D · \Delta h = 0 [/math][math]\rho \gt [/math] [math]\rho _{0}[/math] [math] \frac{ \partial h }{ \partial t } - D · \Delta h = 0 [/math] siendo [math] D=\frac{k}{s}[/math] la difusividad hidráulica que la suponemos constante.

En un caso general , la [math]h[/math] dependería de [math]\rho[/math] y [math]\theta[/math]. Por la definición del Laplaciano, esto quedaría como:[math]\Delta h(\rho,\theta)[/math]. Esta ecuación también se puede obtener si hallamos la divergencia del gradiente de h. Lo primero para hallar q es el gradiente de h (formula gradiente), por tanto q es un campo vectorial. Si se quiere introducir en la ecuación de conservación de la masa, podríamos sustituir el valor de q obtenido dentro de la fórmula que nos obligaría a realizar la divergencia de dicho campo vectorial. Con esto obtendríamos la misma ecuación. Un vez hecho esto para [math]\rho \gt[/math] [math]\rho _{0}[/math] la ecuación diferencial sería:

:[math] \frac{ \partial h }{ \partial t } - D·(\frac{\partial ^2 h}{\partial \rho^2}+ \frac{1}{\rho}·\frac{\partial h}{\partial \rho}+\frac{\partial^2 h}{\partial \theta^2})= 0[/math]

Como queremos imponer que nuestra solución sea radial, obligamos a que [math]h[/math] dependa únicamente de [math]\rho[/math] , por tanto su derivada con respecto a [math]\theta[/math] es igual a cero. La segunda por consiguiente también. Teniendo en cuenta todo esto, obtenemos la siguiente ecuación:

:[math]\frac{\partial h}{\partial t}- D·(\frac{\partial ^2 h}{\partial \rho^2}+ \frac{1}{\rho}·\frac{\partial h}{\partial \rho}) = 0[/math](1)

Ciñéndonos a una región finita, [math]\rho\in (\rho_{0},20)[/math], y considerando que la altura del pozo se mantiene constante [math]h_{p}[/math], la solución de nuestra ecuación (1) requiere dos condiciones de contorno que para este caso, se han propuesto de tipo Dirichlet (valor de la variable preescrito):

[math]h(\rho _{0},t) = h _{\rho} \quad \quad h(20,t) = h _{0}[/math]

Debido al carácter transitorio del problema se necesita una condición inicial, que implica al tiempo dado en un instante inicial [math]t = 0[/math] y que haga referencia al estado del nivel piezométrico [math]h_{0}[/math], antes de ejecutar la excavación del pozo, es decir [math]h(\rho,0) = h _{0}[/math]. Así para completar el sistema imponemos esa condición inicial para que nuestro problema tenga una única solución. Teniendo en cuenta todo lo anterior el sistema será el siguiente:

\[\left\{\begin{matrix}\ \frac{\partial h}{\partial t}- D·(\frac{\partial ^2 h}{\partial \rho^2}+ \frac{1}{\rho}·\frac{\partial h}{\partial \rho})=0 , & \\ h(\rho _{0},t)=h _{\rho} \quad \quad h(20,t)=h _{0} & \\ h(\rho,0)=h _{0} & \end{matrix}\right.\]

2 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}= 1[/math] m, [math]h_{0}= 40[/math] m, [math]h_{p}= 35[/math], [math] D = 0.1[/math], [math]h(\rho,0)= h_{0}[/math] y el tiempo medido en horas.

2.1 Método del trapecio

Plantearemos el problema primero por el método del trapecio, tomando un [math]\Delta \rho = 0.1 [/math] y un [math]\Delta t = \Delta \rho [/math]

2.1.1 Programación en Matlab

%ht-D*(hxx+hx/x)=0
%h(0,t)=hp
%h(L,t)=ho
%h(x,0)=ho
clear all
clc
L=20;
T=500;
D=0.01;
hp=35;
ho=40;
dx=1/10;
x0=1;
x=x0:dx:L;
xint=x0+dx:dx:L-dx;
N=length(x)-2; 
K=(1/dx^2)*(diag(2*ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1));
A=(1/(dx*2))*(diag(ones(1,N-1),1)-diag(ones(1,N-1),-1));
A=diag(1./xint)*A;
h0=(40*ones(N,1));
dt=dx; 
t=0:dt:T;
M=length(t)-1;
F=(zeros(N,1));
F(1)=D*((hp/dx^2)-(hp/(2*dx)*1/xint(1)));
F(N)=D*((ho/dx^2)+ho/(dx*2*xint(N)));
sol(1,:)=[40;h0;40]'; 
h=h0;
for n=1:M
  h=(((eye(N)-(dt*D)*(K-A)))*h+(dt)*(F));
  sol(n+1,:)=[35 h' 40];
end

[xx,tt]=meshgrid(x,t);

surf(xx,tt,sol)
title('Nivel piezométrico en función del tiempo y la distancia');
ylabel('Tiempo');
xlabel('Distancia');
zlabel('Nivel piezométrico')
pause

plot(sol(:,21));
title('Nivel piezométrico en función del tiempo para radio 2');
xlabel('Tiempo');
ylabel('Nivel piezométrico');

2.1.2 Gráfica

En lo referente al metodo trapezoidal podemos observar dos mejoras respecto al otro par de aproximaciones (que se veran a continuacion). Se trata de un metodo implicito, lo cual hace reducir su error y su aproximacion por la tangente no tiene en cuenta tanta cantidad de area externa a la funcion, lo cual ofrece una perspectiva mas realista del proceso natural. En la grafica se puede observar que a tiempo cero, el nivel piezometrico se mantiene constante independientemente de la distancia al pozo. Eso podria interpretarse como que en el instante inicial no se pueden percibir claras diferencias entre el estado natural y estado tras la excavacion. A medida que avanzamos en la linea del tiempo, podemos observar que a radios lejanos el nivel piezometrico se mantiene en la constante inicial. No asi el de aquellos puntos mas cercanos a la excavacion, los cuales van perdiendo altura gradualmente siguiendo una funcion similar a la ecuacion del calor (?) hasta que la curva final puede aproximarse a una parabola cuyo vertice se encuentra al nivel del agua dentro del pozo.

2.2 Euler explícito

2.2.1 Programacion en Matlab

%ht-D*(hxx+hx/x)=0
%h(0,t)=hp
%h(L,t)=ho
%h(x,0)=ho
clear all
clc
L=20;
T=50;
D=0.01;
hp=35;
ho=40;
dx=1/10;
x0=1;
x=x0:dx:L;
xint=x0+dx:dx:L-dx;
N=length(x)-2; 
K=(1/dx^2)*(diag(2*ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1));
A=(1/(dx*2))*(diag(ones(1,N-1),1)-diag(ones(1,N-1),-1));
A=diag(1./xint)*A;
h0=(40*ones(N,1));
dt=dx; 
t=0:dt:T;
M=length(t)-1;
F=(zeros(N,1));
F(1)=D*((hp/dx^2)-(hp/(2*dx)*1/xint(1)));
F(N)=D*((ho/dx^2)+ho/(dx*2*xint(N)));
sol(1,:)=[40;h0;40]'; 
h=h0;
for n=1:M
  h=(((eye(N)-(dt*D)*(K-A)))*h+(dt)*(F));
  sol(n+1,:)=[35 h' 40];
end

[xx,tt]=meshgrid(x,t);

surf(xx,tt,sol)
title('Nivel piezométrico en función del tiempo y la distancia');
ylabel('Tiempo');
xlabel('Distancia');
zlabel('Nivel piezométrico')
pause

plot(sol(:,21));
title('Nivel piezométrico en función del tiempo para radio 2');
xlabel('Tiempo');
ylabel('Nivel piezométrico');

2.2.2 Grafica

2.3 Euler implícito

2.3.1 Programacion en Matlab

%ht-D*(hxx+hx/x)=0
%h(0,t)=hp
%h(L,t)=ho
%h(x,0)=ho
clear all
clc
L=20;
T=500;
D=0.01;
hp=35;
ho=40;
dx=1/10;
x0=1;
x=x0:dx:L;
xint=x0+dx:dx:L-dx;
N=length(x)-2; 
K=(1/dx^2)*(diag(2*ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1));
A=(1/(dx*2))*(diag(ones(1,N-1),1)-diag(ones(1,N-1),-1));
A=diag(1./xint)*A;
h0=(40*ones(N,1));
dt=dx; 
t=0:dt:T;
M=length(t)-1;
F=(zeros(N,1));
F(1)=D*((hp/dx^2)-(hp/(2*dx)*1/xint(1)));
F(N)=D*((ho/dx^2)+ho/(dx*2*xint(N)));
sol(1,:)=[40;h0;40]'; 
h=h0;
for n=1:M
  h=(eye(N)+(dt*D)*(K-A))\(h+(dt)*(F));
  sol(n+1,:)=[35 h' 40];
end

[xx,tt]=meshgrid(x,t);

surf(xx,tt,sol)
title('Nivel piezométrico en función del tiempo y la distancia');
ylabel('Tiempo');
xlabel('Distancia');
zlabel('Nivel piezométrico')
pause

plot(sol(:,21));
title('Nivel piezométrico en función del tiempo para radio 2');
xlabel('Tiempo');
ylabel('Nivel piezométrico');

2.3.2 Grafica

2.4 Comparación de resultados

Como el euler implicito tiene en cuenta el resultado de ese mismo punto para calcular su posicion a lo largo del tiempo, la aproximacion en la grafica respecto al proceso natural se acoplara en mejor medida que el euler explicito. La simplicidad del razonamiento de euler permite la aplicacion de los metodos implicitos con mayor comodidad que cualquier otro metodo iterativo. Se podra observar una mayor poligonalizacion de la curva en cada instante en la grafica del euler explicito debido a que su aproximacion por la tangente se realiza de forma mas burda que en el anterior.

3 Estado estacionario

El nivel piezométrico a lo largo del tiempo va bajando a los puntos mas cercanos a la excavación del pozo, hasta un momento hipotético en el cual la altura piezométrica del punto más cercano a la excavación debe coincidir con la altura del nivel del agua.

Esta hipótesis, solo tomaría lugar en ese supuesto caso de tiempo infinito pero para la modelización, las pequeñas variaciones de altura piezométrica en el punto mas cercano tras el paso de mucho tiempo son ínfimos y despreciables. De ahí que pueda suponerse que en un tiempo lejano las alturas piezométricas de los puntos del terreno no varíen y se genere una situación estacionaria. Las ecuaciones que debe verificar son las de la parábola cuyo vértice se encuentre a la altura del nivel del agua dentro del pozo y al radio infinito tienda a la altura piezométrica inicial en el terreno.

Al ser el caso de un estado estacionario nuestra función no varia en función del tiempo, por tanto implica que: :[math]\frac{\partial h}{\partial t}=0 [/math]

Por consiguiente, la ecuación resultante es: :[math] D·(\frac{\partial ^2 h}{\partial \rho^2}+ \frac{1}{\rho}·\frac{\partial h}{\partial \rho})=0[/math]

3.1 Programación en Matlab

clear all
clc
L=20;
T=5000;
D=0.01;
hp=35;
ho=40;
dx=1/10;
x0=1;
x=x0:dx:L;
xint=x0+dx:dx:L-dx;
N=length(x)-2; 
K=(1/dx^2)*(diag(2*ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1));
A=(1/(dx*2))*(diag(ones(1,N-1),1)-diag(ones(1,N-1),-1));
A=diag(1./xint)*A;
h0=(40*ones(N,1));
dt=200; 
t=0:dt:T;
M=length(t)-1;
F=(zeros(N,1));
F(1)=D*((hp/dx^2)-(hp/(2*dx)*1/xint(1)));
F(N)=D*((ho/dx^2)+ho/(dx*2*xint(N)));
sol(1,:)=[40;h0;40]'; 
h=h0;
for n=1:M
  h=(D)*(K-A)\((F));
  sol(n+1,:)=[35 h' 40];
end
[xx,tt]=meshgrid(x,t);
surf(xx,tt,sol)
title('Nivel piezométrico en función del tiempo y la distancia');
ylabel('Tiempo');
xlabel('Distancia');
zlabel('Nivel piezométrico')
pause
plot(sol(:,21));
title('Nivel piezométrico en función del tiempo para radio 2');
xlabel('Tiempo');
ylabel('Nivel piezométrico');

3.2 Gráfica

Tras haber propuesto el modelo para la solución estacionaria llegamos a la siguiente gráfica, en la que se puede comprobar que el paso del tiempo no influye en el valor del nivel piezométrico del punto, sólo la distancia. Sin embargo se observa, que para valores de tiempo corto, el hecho de que la solución y nuestro estado inicial no concuerden, nos lleva a un tramo en el cual si se aprecian variaciones en función del tiempo, dado que existe una transición bastante breve entre el estado inicial y la solución estacionaria.

4 Método de Fourier

1. Para [math] ħ(\rho,t) = h(\rho,t) - 40 [/math] el sistema será:

[math] \frac{ \partial h }{ \partial t } - D · \Delta h = 0 [/math][math]\rho \gt [/math] [math]\rho _{0}[/math] [math] \frac{ \partial h }{ \partial t } - D · \Delta h = 0 [/math] siendo [math] D=\frac{k}{s}[/math]