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

Revisión actual del 09:58 23 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, \qquad \; \rho \gt \rho _{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]. Esta ecuación también se puede obtener si hallamos la divergencia del gradiente de h, es decir, el Laplaciano. Por la definición del Laplaciano, esto quedaría como:

[math]\Delta h(\rho,\theta) = (\frac{\partial ^2 h}{\partial \rho^2}+ \frac{1}{\rho}·\frac{\partial h}{\partial \rho}+\frac{\partial^2 h}{\partial \theta^2}) [/math]

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 \quad(1)[/math]

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.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 .

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], dado lugar a la siguiente resolución numérica:

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)));
sol1(1,:)=[40;h0;40]'; 
h=h0;
for n=1:M
  h=(((eye(N)-(dt*D)*(K-A)))*h+(dt)*(F));
  sol1(n+1,:)=[35 h' 40];
end

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

subplot(1,2,1);
surf(xx,tt1,sol1)
title('Nivel piezométrico en función del tiempo y la distancia');
ylabel('Tiempo');
xlabel('Distancia');
zlabel('Nivel piezométrico')

%Vario el valor de T
L=20;
T=4;
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)));
sol2(1,:)=[40;h0;40]'; 
h=h0;
for n=1:M
  h=(((eye(N)-(dt*D)*(K-A)))*h+(dt)*(F));
  sol2(n+1,:)=[35 h' 40];
end

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

subplot(1,2,2);
surf(xx,tt2,sol2)
title('Nivel piezométrico en función del tiempo y la distancia');
ylabel('Tiempo');
xlabel('Distancia');
zlabel('Nivel piezométrico')

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

La gráfica que obtenemos es la siguiente

En el primer caso, hemos tomado un tiempo [math] T = 4 [/math](gráfica derecha). Por lo que podemos comprobar, para un tiempo tan pequeño no se aprecia con precisión y claridad el desarrollo y la variación del nivel piezométrico en función del tiempo y la distancia. En el caso de la izquierda, hemos incrementado nuestro valor de [math] T [/math] a [math] T = 500 [/math], por lo que se comprueba como va variando dicho nivel de una manera más clara al tener un paso de tiempo mayor. Como consecuencia nuestra gráfica aparece de color negro, ya que el mallado es muy pequeño debido a que nuestra [math] dx [/math] que nos viene impuesta es demasiado pequeña y el tiempo [math] T [/math] aplicado es muy grande.

El comportamiento de este nivel piezométrico en los puntos para los cuales [math] \rho = 2[/math] en una gráfica bidimensional (nivel piezométrico/tiempo), queda representado en la siguiente gráfica, para un tiempo [math] T = 5000 [/math]. En el programa únicamente cambiaríamos el valor del tiempo, y quitaríamos su comparación con otra posible gráfica de tiempo menor como hemos hecho en el caso anterior.

Podemos observar como evoluciona dicho punto ([math] \rho = 2[/math]) en el tiempo, convirtiéndose en un problema de una única variable. En un inicio éste se encuentra a la condición inicial a la que se encuentra todo el acuífero, es decir, [math] h_{0} = 40 [/math]. Cabe destacar como tras la ejecución del pozo el nivel de agua cae debido al fuerte desnivel causado por la apertura del mismo. A medida que el nivel de agua disminuye, el gradiente ocasionado por éste es menor, ya que también ha disminuido la diferencia entre el nivel piezométrico de nuestro punto y el nivel considerado estable al cual se llega al cabo de cierto tiempo, y en el cual el nivel piezométrico se equilibra y deja de variar. Esto se manifiesta en la figura anterior con una pendiente más o menos pronunciada en su zona inicial que más tarde conforme pasa el tiempo va suavizándose, hasta el momento en el que se convierte en horizontal.

2.2 Método de Euler

Tomando los mismos datos iniciales anteriores, realizamos el problema con el método de euler explícito e implícito respectivamente.

2.2.1 Euler Explícito

A continuación se muestra el programa usado para la resolución numérica de nuestro problema mediante el método de Euler explícito, utilizando los mismo datos iniciales que en el problema resuelto por el método del trapecio:

%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');

La gráfica resultante es:

Esta aproximación no logra ser tan precisa como llega a ser el método del trapecio, por lo que posteriormente realizaremos el método de Euler implícito que mejora la aproximación de la gráfica anterior.

2.2.2 Euler Implícito

En el caso de Euler implícito nuestro programa es:

%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');

Como resultado obtenemos la gráfica siguiente:

La figura como resultado de este método se asemeja mejor a la situación de nuestro problema a estudiar, debido a que hace un refinamiento en la aproximación de Euler explícito, tomando un promedio entre las pendientes, pareciendo este modo más razonable para la obtención de un valor más preciso.

2.3 Comparación de resultados

Realizamos una comparación entre los tres métodos anteriores mediante el siguiente código Matlab:

clear all
clc
L=20;
T=200;
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)+0.5*(dt*D)*(K-A))\(((eye(N)-0.5*(dt*D)*(K-A)))*h+(dt/2)*(F+F));
  sol(n+1,:)=[35 h' 40];
end

[xx,tt]=meshgrid(x,t);
hold on
plot(sol(2001,:),'*');
title('Nivel piezométrico en función del radio');
xlabel('Radio');
ylabel('Nivel piezométrico');
pause
L=20;
T=200;
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);
plot(sol(2001,:),'g -');
pause
L=20;
T=200;
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);
plot(sol(2001,:),'r .');

Las gráficas al ser idénticas debido al pequeño intervalo de tiempo no se diferencien unas de otras, aunque vemos que no permiten que descienda el nivel piezométrico. Necesitariamos un intervalo mayor para poder ver mejor las soluciones, ya que al poner las tres en una misma gráfica todas ellas se superponen.

Por otro lado, podemos comprobar que para tiempos grandes el nivel piezométrico apenas cambia en el tiempo llegando a lo que posteriormente definiremos como estado estacionario.

Su programación es:

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=1; 
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
hold on
plot(sol(1,:));
plot(sol(501,:),'r');
plot(sol(1501,:),'c');
plot(sol(2001,:),'g');
plot(sol(2501,:),'m');
plot(sol(3001,:),'k');
plot(sol(3501,:))
plot(sol(4001,:),'r');
plot(sol(4501,:),'y');
plot(sol(5001,:),'c');
title('Nivel piezométrico en función del radio para distintos tiempos');
xlabel('Radio');
ylabel('Nivel piezométrico');
hold off

Para realizar esta gráfica se han tomado plots desde [math]t = 1[/math] (para matlab no existe el 0) hasta [math] t = 5001[/math] (al no existir 0, el 5000 se convierte en 5001) tomando intervalos de 500 horas entre plot y plot, debido a que tomando un intervalo menor, saldría un número desproporcionado de curvas, que prácticamente no permitirían visualizar correctamente el gráfico. Aun teniendo "solo" 11 curvas, se observa perfectamente como el nivel piezométrico va bajando conforme aumenta el tiempo, pero cada vez baja menos, hasta que llega a un valor cuasiestacionario. El motivo por el que no se observa correctamente la curva de [math]t = 1[/math] es que en este momento el programa interpreta que el nivel piezométrico vale 40 para todo el radio.

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.

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 que se debe verificar en este estado es: :[math] D·(\frac{\partial ^2 h}{\partial \rho^2}+ \frac{1}{\rho}·\frac{\partial h}{\partial \rho})=0[/math]

Resolvemos numéricamente el problema mediante diferencias finitas :

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');

La gráfica resultante es:

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.

Las distancias entre la solución inicial original y la estacionaria para distintos valores de t, 0 100 1000 y 50000, son:

3.1 Recuperación del acuífero

Suponemos que partimos del estado estacionario del apartado anterior e impermeabilizamos las paredes del pozo, o bien dejamos de bombear agua, para conseguir averiguar cual es la capacidad de recuperación del acuífero. En aquellas zonas donde ya dejamos de drenar agua, el acuífero comienza a recuperar un nivel piezométrico constante en todo punto, es decir, tiende a recuperar su valor inicial y por tanto volverá a la situación previa de la realización del pozo.

El código de Matlab que manifiesta esta situación es:

clear all
clc
AAAAAA=[35,35.1858735714553,35.3555777156503,35.5116930348834,35.6562262263997,35.7907684349353,35.9166025197330,36.0347778133697,36.1461635186139,36.2514876719172,36.3513661148673,36.4463243983211,36.5368145916572,36.6232283559116,36.7059072347965,36.7851508450155,36.8612234602207,36.9343593523642,37.0047671615986,37.0726334992756,37.1381259400480,37.2013955232753,37.2625788572152,37.3217998993502,37.3791714708652,37.4347965515233,37.4887693920635,37.5411764741244,37.5920973420963,37.6416053268661,37.6897681778812,37.7366486171181,37.7823048262509,37.8267908764520,37.8701571087410,37.9124504715510,37.9537148211565,37.9939911897550,38.0333180252913,38.0717314065193,38.1092652363055,38.1459514157593,38.1818200014237,38.2168993474649,38.2512162345415,38.2847959868217,38.3176625784293,38.3498387304413,38.3813459994237,38.4122048583733,38.4424347708315,38.4720542588465,38.5010809653860,38.5295317117315,38.5574225503283,38.5847688135163,38.6115851585154,38.6378856090069,38.6636835936105,38.6889919815315,38.7138231156201,38.7381888430646,38.7621005439167,38.7855691576278,38.8086052077595,38.8312188250142,38.8534197687206,38.8752174468947,38.8966209349858,38.9176389934098,38.9382800839600,38.9585523851787,38.9784638067684,38.9980220031096,39.0172343859517,39.0361081363345,39.0546502157951,39.0728673769089,39.0907661732121,39.1083529685458,39.1256339458617,39.1426151155250,39.1593023231478,39.1757012569826,39.1918174549059,39.2076563110160,39.2232230818715,39.2385228923912,39.2535607414378,39.2683415071032,39.2828699517154,39.2971507265824,39.3111883764892,39.3249873439631,39.3385519733194,39.3518865145026,39.3649951267321,39.3778818819665,39.3905507681940,39.4030056925614,39.4152504843482,39.4272888977964,39.4391246148030,39.4507612474825,39.4622023406073,39.4734513739318,39.4845117644072,39.4953868682917,39.5060799831629,39.5165943498368,39.5269331541987,39.5370995289496,39.5470965552740,39.5569272644322,39.5665946392808,39.5761016157256,39.5854510841102,39.5946458905431,39.6036888381670,39.6125826883729,39.6213301619616,39.6299339402552,39.6383966661617,39.6467209451942,39.6549093464470,39.6629644035313,39.6708886154717,39.6786844475664,39.6863543322114,39.6939006696915,39.7013258289400,39.7086321482670,39.7158219360594,39.7228974714532,39.7298610049796,39.7367147591856,39.7434609292312,39.7501016834638,39.7566391639704,39.7630754871091,39.7694127440220,39.7756530011270,39.7817983005940,39.7878506608016,39.7938120767794,39.7996845206331,39.8054699419551,39.8111702682209,39.8167874051717,39.8223232371834,39.8277796276239,39.8331584191975,39.8384614342781,39.8436904752319,39.8488473247285,39.8539337460427,39.8589514833463,39.8639022619902,39.8687877887785,39.8736097522326,39.8783698228486,39.8830696533453,39.8877108789055,39.8922951174097,39.8968239696624,39.9012990196119,39.9057218345633,39.9100939653852,39.9144169467103,39.9186922971300,39.9229215193828,39.9271061005382,39.9312475121738,39.9353472105483,39.9394066367688,39.9434272169529,39.9474103623873,39.9513574696797,39.9552699209077,39.9591490837628,39.9629963116898,39.9668129440220,39.9706003061124,39.9743597094609,39.9780924518371,39.9817998173996,39.9854830768115,39.9891434873513,39.9927822930212,39.9964007246514,40];
L=20;
T=50000;%valor con el que sepamos que se va a llegar al nivel estacionario
D=0.01;
hp=35;
ho=40;
dx=1/10;
x0=1;
x=x0:dx:L;
xint=x0:dx:L-dx;
N=length(x)-1; 
K=(1/dx^2)*(diag(2*ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1));
K(1,2)=-2/dx^2;
A=(1/(dx*2))*(diag(ones(1,N-1),1)-diag(ones(1,N-1),-1));
A=diag(1./xint)*A;
A(1,2)=0;
h0=(AAAAAA(1:N)');
dt=1; 
t=0:dt:T;
M=length(t)-1;
F=(zeros(N,1));
F(1)=D*(-(ho/(2*dx)*1/xint(1)))*0;
F(N)=D*((ho/dx^2)+ho/(dx*2*xint(N)));
sol(1,:)=[h0;40]'; 
h=h0;
for n=1:M
  h=(eye(N)+(dt*D)*(K-A))\(h+(dt)*(F));
  sol(n+1,:)=[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')

La gráfica resultante es la siguiente:

Aquí el nivel en 0 empieza a crecer porque el agua no fluye hacia el pozo.

Para averiguar cuanto tarda nuestro nivel piezométrico en alcanzar el estado estacionario para un error del 5% lo haremos con el siguiente código en Matlab:

clear all
clc
AAAAAA=[35,35.1858735714553,35.3555777156503,35.5116930348834,35.6562262263997,35.7907684349353,35.9166025197330,36.0347778133697,36.1461635186139,36.2514876719172,36.3513661148673,36.4463243983211,36.5368145916572,36.6232283559116,36.7059072347965,36.7851508450155,36.8612234602207,36.9343593523642,37.0047671615986,37.0726334992756,37.1381259400480,37.2013955232753,37.2625788572152,37.3217998993502,37.3791714708652,37.4347965515233,37.4887693920635,37.5411764741244,37.5920973420963,37.6416053268661,37.6897681778812,37.7366486171181,37.7823048262509,37.8267908764520,37.8701571087410,37.9124504715510,37.9537148211565,37.9939911897550,38.0333180252913,38.0717314065193,38.1092652363055,38.1459514157593,38.1818200014237,38.2168993474649,38.2512162345415,38.2847959868217,38.3176625784293,38.3498387304413,38.3813459994237,38.4122048583733,38.4424347708315,38.4720542588465,38.5010809653860,38.5295317117315,38.5574225503283,38.5847688135163,38.6115851585154,38.6378856090069,38.6636835936105,38.6889919815315,38.7138231156201,38.7381888430646,38.7621005439167,38.7855691576278,38.8086052077595,38.8312188250142,38.8534197687206,38.8752174468947,38.8966209349858,38.9176389934098,38.9382800839600,38.9585523851787,38.9784638067684,38.9980220031096,39.0172343859517,39.0361081363345,39.0546502157951,39.0728673769089,39.0907661732121,39.1083529685458,39.1256339458617,39.1426151155250,39.1593023231478,39.1757012569826,39.1918174549059,39.2076563110160,39.2232230818715,39.2385228923912,39.2535607414378,39.2683415071032,39.2828699517154,39.2971507265824,39.3111883764892,39.3249873439631,39.3385519733194,39.3518865145026,39.3649951267321,39.3778818819665,39.3905507681940,39.4030056925614,39.4152504843482,39.4272888977964,39.4391246148030,39.4507612474825,39.4622023406073,39.4734513739318,39.4845117644072,39.4953868682917,39.5060799831629,39.5165943498368,39.5269331541987,39.5370995289496,39.5470965552740,39.5569272644322,39.5665946392808,39.5761016157256,39.5854510841102,39.5946458905431,39.6036888381670,39.6125826883729,39.6213301619616,39.6299339402552,39.6383966661617,39.6467209451942,39.6549093464470,39.6629644035313,39.6708886154717,39.6786844475664,39.6863543322114,39.6939006696915,39.7013258289400,39.7086321482670,39.7158219360594,39.7228974714532,39.7298610049796,39.7367147591856,39.7434609292312,39.7501016834638,39.7566391639704,39.7630754871091,39.7694127440220,39.7756530011270,39.7817983005940,39.7878506608016,39.7938120767794,39.7996845206331,39.8054699419551,39.8111702682209,39.8167874051717,39.8223232371834,39.8277796276239,39.8331584191975,39.8384614342781,39.8436904752319,39.8488473247285,39.8539337460427,39.8589514833463,39.8639022619902,39.8687877887785,39.8736097522326,39.8783698228486,39.8830696533453,39.8877108789055,39.8922951174097,39.8968239696624,39.9012990196119,39.9057218345633,39.9100939653852,39.9144169467103,39.9186922971300,39.9229215193828,39.9271061005382,39.9312475121738,39.9353472105483,39.9394066367688,39.9434272169529,39.9474103623873,39.9513574696797,39.9552699209077,39.9591490837628,39.9629963116898,39.9668129440220,39.9706003061124,39.9743597094609,39.9780924518371,39.9817998173996,39.9854830768115,39.9891434873513,39.9927822930212,39.9964007246514,40];
L=20;
T=50000;%valor con el que sepamos que se va a llegar al nivel estacionario
D=0.01;
hp=35;
ho=40;
dx=1/10;
x0=1;
x=x0:dx:L;
xint=x0:dx:L-dx;
N=length(x)-1; 
K=(1/dx^2)*(diag(2*ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1));
K(1,2)=-2/dx^2;
A=(1/(dx*2))*(diag(ones(1,N-1),1)-diag(ones(1,N-1),-1));
A=diag(1./xint)*A;
A(1,2)=0;
h0=(AAAAAA(1:N)');
dt=1; 
t=0:dt:T;
M=length(t)-1;
F=(zeros(N,1));
F(1)=D*(-(ho/(2*dx)*1/xint(1)))*0;
F(N)=D*((ho/dx^2)+ho/(dx*2*xint(N)));
sol(1,:)=[h0;40]'; 
h=h0;
for n=1:M
  h=(eye(N)+(dt*D)*(K-A))\(h+(dt)*(F));
  sol(n+1,:)=[h' 40];
  Dis=max(abs((sol(n+1,:)-40*ones(1,N+1))/40));
  i=n; %contador de tiempo
  if Dis<=0.05 %distancia entre soluciones <=5%
      break %una vez alcanzado el valor estacionario con un error <=5% se rompe el bucle
  end
end
i %tiempo necesario

Lo que podemos concluir en que el tiempo de recuperación para un error del 5% son 830 horas.

4 Método de Fourier

Mediante el método de Fourier realizamos una aproximación alternativa a las anteriormente realizadas. Suponiendo que hemos extraido el agua, el pozo comienza a cerrarse quedando un acuífero nuevamente homogéneo. Para esto, consideramos despreciable el radio del pozo, diremos que [math]\rho_{0}=0[/math]. Para [math] ħ(\rho,t) = h(\rho,t) - 40 [/math], siendo la condición [math] h(20,t)=40 [/math] el sistema será:

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

Para que se cumpla la condición de diferenciabilidad en [math]\rho_{0}=0[/math] se tendrá que imponer el requesito [math]\rho'(0)=0[/math] ya que no puede haber ningún pico en dicho punto.

Teniendo en cuenta lo anterior el problema de autovalores asociado obtenido es:

[math]\frac{\partial ħ}{\partial t}- D·(\frac{\partial ^2 ħ}{\partial \rho^2}+ \frac{1}{\rho}·\frac{\partial ħ}{\partial \rho})=0 \ltbr /\gt \qquad \; ħ_{t}-D·(ħ_{\rho\rho}-\frac{1}{\rho}·ħ_{\rho})=0[/math], es decir, [math]T'(t)· R(\rho) - D·(R''(\rho)·T(t)+ \frac{1}{\rho}· R'(\rho)·T(t))=0[/math] [math] \qquad \; [/math] [math]T'(t)·R(\rho)=T(t)·D·(R''(\rho)+\frac{1}{\rho}·R(\rho))[/math] [math] \qquad [/math] [math]\frac{T'(t)}{T(t)}= - λ =\frac{D}{R(\rho)}·(R''(\rho)+\frac{1}{\rho}·R'(\rho))[/math] \[\left\{\begin{matrix}\ T'(t)+ λ·T(t)= 0 & \\ \quad R(\rho)+ \frac{1}{\rho}·R'(\rho)+ λ· \frac{R(\rho)}{D}=0& \end{matrix}\right.\]

Como se puede observar nos queda similar a la función ofrecida en el enunciado, con condiciones de frontera homogéneas, y las soluciones de este problema tendrán la forma de función de Bessel de primera especie que forma una familia de autofunciones: [math]R_{k}(\rho)=J_{0}·(\rho\sqrt{λ})[/math] , que se anulan cuando [math](\rho\sqrt{λ})[/math] se iguala a 2.4048, 5.5201, 8.6537, 11.7915 o 14.9309, con lo cual, y para [math] \rho = 20[/math] podemos obtener los 5 primeros autovalores de este problema de contorno asociado.

Como se puede observar nos queda similar a la función de Bessel que tendrá como autofunciones [math]R_{k}(\rho)=J_{0}·(\rho\sqrt{λ})[/math], siendo las autofunciones:

• [math]20·(\sqrt{λ_{1}})=2,4048 \qquad λ_{1}=0,0145[/math]
• [math]20·(\sqrt{λ_{2}})=5,5201 \qquad λ_{2}=0,0762[/math]
• [math]20·(\sqrt{λ_{3}})=8,6537 \qquad λ_{3}=0,1872[/math]
• [math]20·(\sqrt{λ_{4}})=11,7915 \qquad λ_{4}=0,3480[/math]
• [math]20·(\sqrt{λ_{5}})=14,9309 \qquad λ_{5}=0,5573[/math]

4.1 Resolución numérica

Para obtener la gráfica con las primeras cinco autofunciones, utilizaremos la función besselj de Matlab:

clear all
clc
a1=0.0145;
a2=0.0762;
a3=0.1872;
a4=0.3476;
a5=0.5575;
r=0:0.02:20;
R1=besselj(0,r.*(sqrt(a1))); %Ecuacion de Bessel para el primer cero
R2=besselj(0,r.*(sqrt(a2)));
R3=besselj(0,r.*(sqrt(a3)));
R4=besselj(0,r.*(sqrt(a4)));
R5=besselj(0,r.*(sqrt(a5)));
hold on
plot(r,R1)
plot(r,R2,'g')
plot(r,R3,'y')
plot(r,R4,'r')
plot(r,R5,'c')
hold off

Obtenemos la siguiente gráfica:

Por lo que podemos observar en la gráfica es que la función de Bessel es una función oscilante(como las funciones seno o coseno) que decaen proporcionalmente a [math]1/\sqrt{x}[/math].

Por último, aproximamos nuestra ecuación cuando el dato inicial es [math] h (\rho,0)= -log \frac{|\rho + 0.01|}{20.01} [/math] usando cinco términos del desarrollo de Fourier. De la serie de Fourier de [math] h(\rho,0)= -log \frac{|\rho + 0.01|}{20.01} [/math], obtenemos los coeficientes de Fourier. En este caso al tratarse de funciones de Bessel, las autofunciones son ortogonales respecto al producto escalar, es decir [math] (f(\rho),g(\rho)) = \int_{0}^{20} \rho f(\rho) g(\rho) \cdot dx [/math]

x0=0; L=20;
ceros=[2.4048,5.5201,8.6537,11.7915,14.9309]; %Ceros de la función de Bessel
A=(ceros./20).^2; %Autovalores
K=5; %Número de autovalores y autofunciones buscadas
x=x0:0.1:L;
N=length(x)-1;
J=zeros(N+1,K);
z=zeros(N+1,K);
t0=0; T=500;
t=t0:0.1:T;
[xx,tt]=meshgrid(x,t);
h=0;
for k=1:K
   z(:,k)=(x)'*sqrt(A(k));
   J(:,k)=besselj(0,z(:,k));
   h0=-log((abs(x+0.1))/20.01); %Condicion iniciales
   c(k)=trapz(x,x'.*h0'.*J(:,k))/trapz(x,x'.*J(:,k).*J(:,k)); %Coeficientes de Fourier
   h=h+c(k)*exp(-A(k)*t)'*J(:,k)';
end
%Representación gráfica de la solución obtenida
surf(xx,tt,h)