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

De MateWiki
Saltar a: navegación, buscar

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, usando el Método de Heun para realizar la aproximación en tiempo.

Primero, sustituyendo los datos anteriores en la ecuación, llegamos a:

[math] \left\{\begin{matrix} \frac{\partial h}{\partial t} - 10^{2} (\frac{1}{\rho} \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]

Aplicando el método de las diferencias finitas llegamos al siguiente código: 

%TRPRIMERA FASE. ESPACIO

a=0.5; b=20; %Condiciones de contorno
drho=0.1; N=(b-a)/drho; %División del intervalo

rho=a+drho:drho:b-drho; %Sólo nodos interiores

%Creamos las matrices

L=diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
A=diag(1./rho);
L=10^(-2)*(1/(2*drho))*A*L;

K=-2*diag(ones(1,N-1))+diag(ones(1,N-2),1)+diag(ones(1,N-2),-1);
K=10^(-2)*(1/drho^2)*K; 

F=zeros(N-1,1);
F(1,1)=10^(-2)*(-17/(rho(1)*drho)+34/drho^2);
F(N-1,1)=10^(-2)*(45/(rho(N-1)*2*drho)+45/drho^2);

%SEGUDA FASE. TIEMPO

H=ones(1,length(rho)); H=H*45; H=H';
dt=drho; T=10; t=0:dt:T;
h0=34; hN=45;
sol=[h0 H' hN];

for j=1:length(t)-1 
    
    %Método de Heun
    k1=F+(K+L)*H;
    k2=F+(K+L)*(H+k1*dt);
    H=H+(dt/2)*(k1+k2);
 
    sol(j+1,:)=[h0 H' hN]; %Salvamos las solución, incluimos los extremos
end

rho1=a:drho:b; %Incluimos los nodos exteriores

[rr,tt]=meshgrid(rho1,t);
surf(rr,tt,sol); %Dibujamos la solución
title('Nivel Piezométrico. Método Diferencias finitas con Heun en tiempo.')
xlabel('Distancia radial al centro del pozo.')
ylabel('Tiempo.')
zlabel('Nivel piezométrico.')


Mediante este código obtenemos la siguiente gráfica:

Gráfica Nivel Piezométrico aproximada con el Método de Diferencias Finitas con Heun en tiempo.



















Ahora dibujaremos la gráfica mediante los métodos de Euler implícito y trapecio y compararemos los resultados obtenidos.


Sobre el programa anterior cambiaremos las líneas correspondientes al método de Heun.


  • Método de Euler implícito:

[math] \left\{\begin{matrix} y_{0}\\ y_{n+1}=y_{n} + hf(t_{n+1},y_{n+1})\\ \end{matrix}\right. [/math] 


for j=1:length(t)-1 
    
    %Método de Euler implícito
    
    H=(eye(N-1)-dt*(K+L))\(H+(F*dt)); 
 
    sol(j+1,:)=[h0 H' hN]; %Salvamos las solución, incluimos los extremos
end

rho1=a:drho:b; %Incluimos los nodos exteriores

[rr,tt]=meshgrid(rho1,t);
surf(rr,tt,sol); %Dibujamos la solución
title('Nivel Piezométrico. Método Diferencias finitas con Euler en tiempo.')
xlabel('Distancia radial al centro del pozo.')
ylabel('Tiempo.')
zlabel('Nivel piezométrico.')


  • Método del Trapecio:

[math] \left\{\begin{matrix} y_{0}\\ y_{n+1}=y_{n} + \frac{h}{2}f((t_{n},y_{n}+f(t_{n+1},y_{n+1}))\\ \end{matrix}\right. [/math] 

for j=1:length(t)-1 
    
    %Método del trapecio
    k1=F+(K+L)*H;
    k2=F+(K+L)*(H+k1*dt); 
    H=(eye(N-1)-(dt/2)*(K+L))\(H+dt*F+(dt/2)*(K+L)*H);
 
    sol(j+1,:)=[h0 H' hN]; %Salvamos las solución, incluimos los extremos
end

rho1=a:drho:b; %Incluimos los nodos exteriores

[rr,tt]=meshgrid(rho1,t);
surf(rr,tt,sol); %Dibujamos la solución
title('Nivel Piezométrico. Método Diferencias finitas con trapecio en tiempo.')
xlabel('Distancia radial al centro del pozo.')
ylabel('Tiempo.')
zlabel('Nivel piezométrico.')



Como podemos ver a continuación las tres gráficas resultantes son prácticamente iguales por lo que los tres métodos hacen una aproximación similar.


Gráfica Nivel Piezométrico aproximada con el Método de Diferencias Finitas con Trapecio en tiempo.
Gráfica Nivel Piezométrico aproximada con el Método de Diferencias Finitas con Euler Implícito en tiempo.































































5 Comportamiento del nivel piezométrico con [math] \rho [/math] constante

Ahora supondremos que [math] \rho [/math] es constante y que su valor es 2 y dibujaremos el comportamiento del nivel piezométrico en una gráfica 2D nivel piezométrico-tiempo.

Nuevamente utilizando el programa anterior con el método de Heun:

Gráfica Nivel Piezométrico-tiempo para [math] \rho =2 [/math] para T=10.
Gráfica Nivel Piezométrico-tiempo para [math] \rho =2 [/math] para T=1000.































































%PRIMERA FASE. ESPACIO

a=0.5; b=20; %Condiciones de contorno
drho=0.1; N=(b-a)/drho; %División del intervalo

rho=a+drho:drho:b-drho; %Sólo nodos interiores

%Creamos las matrices

L=diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
A=diag(1./rho);
L=10^(-2)*(1/(2*drho))*A*L;

K=-2*diag(ones(1,N-1))+diag(ones(1,N-2),1)+diag(ones(1,N-2),-1);
K=10^(-2)*(1/drho^2)*K; 

F=zeros(N-1,1);
F(1,1)=10^(-2)*(-17/(rho(1)*drho)+34/drho^2);
F(N-1,1)=10^(-2)*(45/(rho(N-1)*2*drho)+45/drho^2);

%SEGUDA FASE. TIEMPO

H=ones(1,length(rho)); H=H*45; H=H';
dt=drho; 
T=10 
%Comprobamos también la solución para un tiempo mayor 
%T=1000
t=0:dt:T;
h0=34; hN=45;
sol=[h0 H' hN];

for j=1:length(t)-1 
    
    %Método de Heun
    k1=F+(K+L)*H;
    k2=F+(K+L)*(H+k1*dt);
    H=H+(dt/2)*(k1+k2);
 
    sol(j+1,:)=[h0 H' hN]; %Salvamos las solución, incluimos los extremos
end

rho1=a:drho:b; %Incluimos los nodos exteriores

m=sol(:,rho==2); %Vector que incluye la solución en tiempo para rho=2

plot(t,m);

title('Comportamiento del nivel piezométrico en el tiempo para p=2')
xlabel('Tiempo')
ylabel('Nivel piezométrico.')



La primera gráfica muestra que para tiempos pequeños la variación del nivel piezométrico es muy pequeña (del orden de 0,10/0,15 m), aún así como observamos en la segunda gráfica a medida que avanzamos en el tiempo la variación del nivel disminuye aún más ya que la curva está cada vez menos inclinada.

6 Variación del nivel piezométrico en tiempos grandes

Ahora vamos a comprobar que la variación del nivel piezométrico para tiempos elevados (T=5000) es mínima. Para ello dibujaremos las soluciones cada 100 horas de en el intervalo temporal [1; 5000], en una misma gráfica.

Nivel piezométrico para tiempos grandes





































%PRIMERA FASE. ESPACIO

a=0.5; b=20; %Condiciones de contorno
drho=0.1; N=(b-a)/drho; %División del intervalo

rho=a+drho:drho:b-drho; %Sólo nodos interiores

%Creamos las matrices

L=diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
A=diag(1./rho);
L=10^(-2)*(1/(2*drho))*A*L;

K=-2*diag(ones(1,N-1))+diag(ones(1,N-2),1)+diag(ones(1,N-2),-1);
K=10^(-2)*(1/drho^2)*K; 

F=zeros(N-1,1);
F(1,1)=10^(-2)*(-17/(rho(1)*drho)+34/drho^2);
F(N-1,1)=10^(-2)*(45/(rho(N-1)*2*drho)+45/drho^2);

%SEGUDA FASE. TIEMPO

H=ones(1,length(rho)); H=H*45; H=H';
dt=0.1; T=5000; t=0:dt:T;
h0=34; hN=45;
sol=[h0 H' hN];

figure(1);
hold on;


for j=1:length(t)-1 
    
    %Método de Heun
    k1=F+(K+L)*H;
    k2=F+(K+L)*(H+k1*dt);
    H=H+(dt/2)*(k1+k2);
 
    sol(j+1,:)=[h0 H' hN]; %Salvamos las solución, incluimos los extremos
    
 if rem(t(j),100)==0; 
    
    rho1=a:drho:b; %Incluimos los nodos exteriores
  
  plot(rho1,sol(j,:)) %pintamos la solución cada 100 horas
  
  title('Variación del nivel piezométrico para tiempos grandes')
  pause(1)

end

end



Como se puede observar en la gráfica según vamos avanzando en el tiempo las curvas están más juntas lo que quiere decir que la variación del nivel es cada vez menor aproximándose a un régimen estacionario.

Nivel piezométrico estacionario

7 Problema estacionario

Como hemos visto en el apartado anterior para tiempos grandes el nivel piezométrico se mantiene practicamente constante a lo largo del tiempo, esto es debido a que se establece un equilibrio entre el caudal de agua extraida del acuifero y la migración de aguas desde los puntos más alejados del pozo hacia este. Por lo tanto se verifica que: 

[math]\frac{\partial h}{\partial t}=0[/math]

Y por lo tanto el nivel piezométrico se puede aproximar por:

[math] \left\{\begin{matrix} - \frac{k}{s} (\frac{1}{\rho} \frac{\partial h}{\partial \rho} + \frac{\partial^2 h}{\partial \rho^2})= 0 & \rho \epsilon (1, 20) & \theta \epsilon (0, 2\Pi ) & 0 \lt t \lt 5000 \\ h(1,t) = 35; h(20, t)= 40 \\ h(\rho, 0) = 40 \end{matrix}\right. [/math] Resolvemos por diferencias finitas: 

%datos
    t=5000;
    l=20;
    D=0.01;
    hp=35;
    ho=40;
    %discretización
    s=0.1;
    x0=1;
    x=x0:s:l;
    xint=x0+s:s:l-s;
    N=length(x)-2; 
    %construimos las matrices para el método de las diferencias finitas
    K=(1/s^2)*(diag(2*ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1));
    A=(1/(s*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/s^2)-(hp/(2*s)*1/xint(1)));
    F(N)=D*((ho/s^2)+ho/(s*2*xint(N)));
    H(1,:)=[40;h0;40]'; 
    h=h0;
    for n=1:M
      h=(D)*(K-A)\((F)); %utilizamos el comando backslash para resolver el sistema de ecuaciones
      H(n+1,:)=[hp h' ho];
    end
    %mallado y representación
    [X,T]=meshgrid(x,t);
    surf(X,T,H)
    xlabel('Distancia (metros)');
    ylabel('Tiempo (horas)');
    zlabel('Nivel piezométrico (metros)')

Efectivamente, como podemos comprobar en la gráfica, la solución es independiente del tiempo salvo para tiempos pequeños en los cuales no queda muy bien definida

También adjuntamos un programa que calcula el error entre la aproximación del problema considerándolo estacionario y la aproximación original para ciertos instantes: 

%TRPRIMERA FASE. ESPACIO
 
a=0.5; b=20; %Condiciones de contorno
drho=0.1; N=(b-a)/drho; %División del intervalo
 
rho=a+drho:drho:b-drho; %Sólo nodos interiores
 
%Creamos las matrices
 
L=diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
A=diag(1./rho);
L=10^(-2)*(1/(2*drho))*A*L;
 
K=-2*diag(ones(1,N-1))+diag(ones(1,N-2),1)+diag(ones(1,N-2),-1);
K=10^(-2)*(1/drho^2)*K; 
 
F=zeros(N-1,1);
F(1,1)=10^(-2)*(-17/(rho(1)*drho)+34/drho^2);
F(N-1,1)=10^(-2)*(45/(rho(N-1)*2*drho)+45/drho^2);
 
%SEGUDA FASE. TIEMPO
 
H=ones(1,length(rho)); H=H*45; H=H';
dt=1; T=5000; t=0:dt:T;
h0=34; hN=45;
sol=[h0 H' hN];
 
for j=1:length(t)-1 
 
    %Método de Heun
    k1=F+(K+L)*H;
    k2=F+(K+L)*(H+k1*dt);
    H=H+(dt/2)*(k1+k2);
 
    sol(j+1,:)=[h0 H' hN]; %Salvamos las solución, incluimos los extremos
end



%CALCULAMOS SOLUCIÓN ESTACIONARIA

%PRIMERA FASE. ESPACIO
 
a=0.5; b=20; %Condiciones de contorno
drho=0.1; N=(b-a)/drho; %División del intervalo
 
rho=a+drho:drho:b-drho; %Sólo nodos interiores
 
%Creamos las matrices
 
L=diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
A=diag(1./rho);
L=10^(-2)*(1/(2*drho))*A*L;
 
K=-2*diag(ones(1,N-1))+diag(ones(1,N-2),1)+diag(ones(1,N-2),-1);
K=10^(-2)*(1/drho^2)*K; 
 
F=zeros(N-1,1);
F(1,1)=10^(-2)*(-17/(rho(1)*drho)+34/drho^2);
F(N-1,1)=10^(-2)*(45/(rho(N-1)*2*drho)+45/drho^2);
 
%SEGUDA FASE. TIEMPO
 
H=ones(1,length(rho)); H=H*45; H=H';
dt=1; T=5000; t=0:dt:T;
h0=34; hN=45;
solest=[h0 H' hN];
 
for j=1:length(t)-1 
 
    H=-(L+K)\F;
 
    solest(j+1,:)=[h0 H' hN]; %Salvamos las solución, incluimos los extremos
end

e0=max(abs(sol(1,:)-solest(1,:)))
e100=max(abs(sol(101,:)-solest(101,:)))
e1000=max(abs(sol(1001,:)-solest(1001,:)))
e5000=max(abs(sol(5001,:)-solest(5001,:)))

Error para t=0: 0.
error para t=100:5.6443e+67
error para t=1000: 0
error para t=5000: 0

8 Capacidad de recuperación del acuífero

En este apartado queremos calcular una ley que nos describa el comportamiento del acuífero e el caso de que impermeabilizásemos las paredes del pozo partiendo del estado estacionario anterior, esto significaría que el fujo de agua en los bordes del pozo sería nulo y el nivel piezométrico iría variando hasta llegar a otro estado estacionario. El estado estacionarioque tenderá a alcanzar el fluido será el inicial, ya que al impermeabilizar las paredes volvemos al caso del acuífero confinado y por ello es nivel piezométrico quedará constante a lo largo del tiempo y de la extensión del acuífero. La condición frontera a imponer será que el flujo de agua en r0 será 0, ya que estará impermeabilizado. Esto lo podemos plasmar diciendo que la parcial del nivel piezométrico con respecto a ro en r0 es 0, ya que un aislamiento impide el flujo de líquido, y si el flujo de líquido es nulo el nivel piezométrico no varía. [math]\frac{\partial h}{\partial \rho}=0[/math]

Con la condición que acabamos de proponer resolvemos numéricamente con Matlab:

Recuperación del acuífero tras impermeabilizar





















clear all

%REPETIMOS EL PROGRAMA ANTERIOR PARA SACAR EL VALOR ESTACIONARIO
 
a=0.5; b=20; %Condiciones de contorno
drho=0.1; N=(b-a)/drho; %División del intervalo
rho=a:drho:b-drho; %Sólo nodos interiores
L=diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
L(1,2)=0;
A=diag(1./rho);
L=10^(-2)*(1/(2*drho))*A*L;
K=-2*diag(ones(1,N))+diag(ones(1,N-1),1)+diag(ones(1,N-1),-1);
K(1,2)=2;
K=10^(-2)*(1/drho^2)*K; 
F=zeros(N,1);
F(1,1)=0; %Ahora esto es cero
F(N,1)=10^(-2)*(45/(rho(N)*2*drho)+45/drho^2);
H=ones(1,length(rho)); H=H*45; H=H';
dt=200; T=5000; t=0:dt:T;
h0=34; hN=45;
sol=[h0 H' hN];
 
for j=1:length(t)-1 
 
    H=-(L+K)\F;
    
    MAT=[h0 H' hN];
 
    sol(j+1,:)=[h0 H' hN]; %Salvamos las solución, incluimos los extremos
end



%AHORA RESOLVEMOS PARA EL CASO DE LA RECUPERACIÓN

%PRIMERA FASE. ESPACIO
 
a=0.5; b=20; %Condiciones de contorno
drho=0.1; N=(b-a)/drho; %División del intervalo
 
rho=a:drho:b-drho; %Sólo nodos interiores
 
%Creamos las matrices
 
L=diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
L(1,2)=0;
A=diag(1./rho);
L=10^(-2)*(1/(2*drho))*A*L;
 
K=-2*diag(ones(1,N))+diag(ones(1,N-1),1)+diag(ones(1,N-1),-1);
K(1,2)=2;
K=10^(-2)*(1/drho^2)*K; 
 
F=zeros(N,1);
F(1)=0; %Ahora esto es cero
F(N)=10^(-2)*(45/(rho(N)*2*drho)+45/drho^2);
 
%SEGUDA FASE. TIEMPO
 
H=ones(1,length(rho)); 
dt=200; T=5000; t=0:dt:T;
hN=45;
sol=[MAT(1:length(rho)) hN];

for j=1:length(t)-1 
    
    H=(-(L+K))\F;
 
    sol(j+1,:)=[H' hN]; %Salvamos las solución, incluimos los extremos
end
 
rho1=a:drho:b; %Incluimos los nodos exteriores
 
[rr,tt]=meshgrid(rho1,t);
surf(rr,tt,sol); %Dibujamos la solución
title('Recuperación del acuífero. Condición estacionaria.')
xlabel('Distancia radial al centro del pozo.')
ylabel('Tiempo.')
zlabel('Nivel piezométrico.')

Este programa nos representa la evolución del nivel piezométrico del acuífero tras ser impermeabilizado, partiendo del estado estacionario descrito anteriormente hasta llegar al estado inicial con h constante.


También construimos un programa aparte que calcula, con un error del 5%, cuánto dura esta recuperación del nivel piezométrico: 

clear all

%PRIMERA FASE. ESPACIO
 
a=0.5; b=20; %Condiciones de contorno
drho=0.1; N=(b-a)/drho; %División del intervalo
 
rho=a:drho:b-drho; %Sólo nodos interiores
 
%Creamos las matrices
 
L=diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
L(1,2)=0;
A=diag(1./rho);
L=10^(-2)*(1/(2*drho))*A*L;
 
K=-2*diag(ones(1,N))+diag(ones(1,N-1),1)+diag(ones(1,N-1),-1);
K(1,2)=2;
K=10^(-2)*(1/drho^2)*K; 
 
F=zeros(N,1);
F(1,1)=0; %Ahora esto es cero
F(N,1)=10^(-2)*(45/(rho(N)*2*drho)+45/drho^2);
 
%SEGUDA FASE. TIEMPO
 
H=ones(1,length(rho)); H=H*45; H=H';
dt=1; T=5000; t=0:dt:T;
h0=34; hN=45;
sol=[h0 H' hN];
 
for j=1:length(t)-1 
 
    H=-(L+K)\F;
    
    MAT=[h0 H' hN];
 
    sol(j+1,:)=[h0 H' hN]; %Salvamos las solución, incluimos los extremos
end
 

%PRIMERA FASE. ESPACIO
 
a=0.5; b=20; %Condiciones de contorno
drho=0.1; N=(b-a)/drho; %División del intervalo
 
rho=a:drho:b-drho; %Sólo nodos interiores
 
%Creamos las matrices
 
L=diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);
L(1,2)=0;
A=diag(1./rho);
L=10^(-2)*(1/(2*drho))*A*L;
 
K=-2*diag(ones(1,N))+diag(ones(1,N-1),1)+diag(ones(1,N-1),-1);
K(1,2)=2;
K=10^(-2)*(1/drho^2)*K; 
 
F=zeros(N,1);
F(1,1)=0; %Ahora esto es cero
F(N,1)=10^(-2)*(45/(rho(N)*2*drho)+45/drho^2);
 
%SEGUDA FASE. TIEMPO
 
H=ones(1,length(rho)); 
dt=1; T=5000; t=0:dt:T;
hN=45;
sol=[MAT(1:length(rho)) hN];

for j=1:length(t)-1 
    
    H=-(L+K)\F;
 
    sol(j+1,:)=[H' hN]; %Salvamos las solución, incluimos los extremos
    
    error=max(abs(sol(j+1,:)-hN*ones(1,N+1)));
    
    tiempo=j; %Contador de tiempo, cada bucle una hora
    if error <=0.05
        break
    end
end

tiempo %Saca en pantalla el tiempo que tarda en alcanzar el régimen estacionario

El resultado es de 1 hora. *Nota: los programas son reiterativos y definen las mismas variables por estar diseñados para compilar por separado.

9 Método de Fourier

Ahora suponemos que el radio del pozo es despreciable y que una vez extraido el agua el pozo se cierra quedando el acuifero de nuevo homogéneo. En este caso tomaremos la región ρ < 20. Nuestro objetivo es aproximar las soluciones usando el método de Fourier. De nuevo supondremos que la solución sólo depende de ρ y t, y tomaremos ahora como condición frontera h(20, t) = 40. Vamos a considerar el sistema de ecuaciones para h(ρ,t)=h(ρ,t)-40: 

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

Obtenemos el problema de autovalores asociado al sistema anterior con el método de separación de variables: 

[math] \left\{\begin{matrix}
D (\frac{X'(\rho)}{\rho} + X''(\rho))+ λX(\rho)=0 \\
X'(0) = 0;   X(20)= 0 \\ 

\end{matrix}\right. [/math]
derecha

Sabemos que los problemas de esta forma tienen como solución la función de Bessel de primera especie que forma una familia de autofunciones: [math]Ψ(\rho)=J_0·(\rho\sqrt{λ})[/math]. Estas se anulan cuando [math](\rho\sqrt{λ})[/math] se iguala a 2.4048, 5.5201, 8.6537, 11.7915 y 14.9309, siendo estos los 5 primeros ceros de la función de Bessel. De esta manera y tomando ρ=20 obtenemos los 5 primeros autovalores: [math](20\sqrt{λ})=2.4048[/math]Haciendo la operación anterior para los 5 ceros obtenemos los siguientes autovalores: 0.0145, 0.0762, 0.1872, 0.3476, 0.5573. Ahora vamos a dibujar las 5 primeras autofunciones utilizando la función besselj: 

a1=0.0145;
a2=0.0762;
a3=0.1872;
a4=0.3476;
a5=0.5575;
r=0:0.02:20;
F1=besselj(0,r.*(sqrt(a1))); 
F2=besselj(0,r.*(sqrt(a2)));
F3=besselj(0,r.*(sqrt(a3)));
F4=besselj(0,r.*(sqrt(a4)));
F5=besselj(0,r.*(sqrt(a5)));
hold on
plot(r,F1)
plot(r,F2,'r')
plot(r,F3,'g')	
plot(r,F4,'c')
plot(r,F5,'y')
hold off
derecha

Aproximamos ahora la solución de la ecuación tomando como dato inicial es[math] h (\rho,0)= -log \frac{|\rho + 0.01|}{20.01} [/math] y usando 5 términos del desarrollo de Fourier, teniendo en cuenta que las autofunciones son ortogonales respecto al producto escalar: [math] (f(\rho),g(\rho)) = \int_{0}^{20} \rho f(\rho) g(\rho) \cdot dx [/math] 





























x0=0; L=20;
c0=[2.4048,5.5201,8.6537,11.7915,14.9309]; %Vector de ceros de la función de Bessel
A=(c0./20).^2; %Vector de autovalores
K=5; %Número de términos del desarrollo de Fourier
x=x0:0.1:L;
N=length(x)-1;
J=zeros(N+1,K);
z=zeros(N+1,K);
t0=0; tf=500;
t=t0:0.1:tf;
[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); %Condición inicial
   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
surf(xx,tt,h)
title('Aproximación de la solución')
Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Nivel_piezométrico_de_un_acuífero_confinado_(A5)&oldid=30535»