Nivel piezométrico en un acuífero confinado. Grupo 8-B

De MateWiki
Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Nivel piezométrico en acuífero confinado. Grupo 8-B
Asignatura Ecuaciones Diferenciales
Curso Curso 2016-17
Autores José Luis Abia Pascual, Jorge Crespo Jornet, Ignacio Puente Dans, Mario Lamparero de Lucas, Víctor Valenzuela Collado, David Espinosa Aller, Ignacio Sastre Jiménez
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Se entinede como nivel piezométrico a la altura de la superficie libre de agua sobre el nivel del mar, en los acuíferos libres. En los confinados, es la altura que alcanzaría el agua en el interior de un sondeo hasta equilibrarse con la presión atmosférica. Sobre el acuífero se construye un pozo de sección circular y radio [math]\rho _{0} [/math] ,provocando que el nivel piezométrico cambie. Para realizar el estudio utilizaremos la ecuación de conservación de la masa y la ley de Darcy, mediante la que modelaremos el flujo de agua a través de un medio poroso, definiendo que el flujo de agua [math] q [/math] es proporcional a la diferencia de presión.Si consideramos que nos encontramos en un medio poroso homogéno e isótropo, podremos realizar una buena aproximación del comportamiento del agua mediante la Ley de Darcy.

Dichas ecuaciones son:

La ecuación de conservación de la masa: [math] S * \frac{∂h}{∂t} + div(q) = 0 [/math]

La Ley de Darcy: [math] q = −K*∇h [/math]

Croquis del acuífero confinado

1 Cálculo del Laplaciano y ecuación diferencial en coordenadas polares

Si operamos con las dos ecuaciones iniciales tenemos [math] \frac{ \partial h }{ \partial t } - D · \Delta h = 0 [/math] , siendo[math] \quad \rho \gt \rho _{0} [/math], [math] \quad θ\in (0,2\pi ) \quad t\gt0 [/math]

tomando [math] D= \frac{k}{s}[/math] como la conductibilidad hidráulica.

Partimos de la funcion [math] h(\rho,\theta)[/math]

Hacemos el gradiente de la función y luego la divergencia de este y así obtenemos el Laplaciano: : [math] \Delta h=(\frac{\partial ^2 h}{\partial \rho^2}+ \frac{1}{\rho}·\frac{\partial h}{\partial \rho}+\frac{\partial^2 h}{\partial \theta^2}) [/math]

Combinamos el Laplaciano con la ecuación obtenida de operar las dos ecuaciones iniciales y obtenemos: : [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]

Para obener una solución radial, [math]h[/math] ha de depender únicamente de [math]ρ[/math], y así tenemos que la la derivada con respecto a θ sea igual a 0, obteniendo así la ecuación final: : [math]\frac{\partial h}{\partial t}- D·(\frac{\partial ^2 h}{\partial \rho^2}+ \frac{1}{\rho}·\frac{\partial h}{\partial \rho}) = 0[/math]

2 Condiciones de contorno y sistema de ecuaciones

Las condiciones de contorno impuestas (tipo Dirichlet), son una altura de pozo constante [math]h_{p}[/math] y [math]\rho\in (\rho_{0},20)[/math]. Las condiciones serán: : [math]h(\rho _{0},t) = h _{\rho}\quad (1)[/math]: [math] h(20,t) = h _{0}\quad (2)[/math]

Para que se tenga una solución única, tenemos que poner una última condición para [math]t = 0[/math]: : 

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

Nuestro problema se representa en el siguiente sistema: \[\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.\]

3 Resolución del problema por el método de diferencias finitas

La aproximación de las derivadas por diferencias finitas desempeña un papel central en los métodos de diferencias finitas del análisis numérico para la resolución de ecuaciones diferenciales.

Tomando las condiciones iniciales del enunciado, resolvemos el problema por el método de diferencias finitas y aplicando el método del trapecio en tiempo, siendo [math]\Delta \rho = 0.1 [/math] , [math]\Delta t = \Delta \rho [/math] , [math] D= 0.01[/math] y [math]\rho _{0} = 0.5 [/math], de tal forma el sistema será: \[\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)=34 \quad \quad h(20,t)=45 & \\ h(\rho,0)=45 & \end{matrix}\right.\] 

%trabajo Ecus Pregunta 3
%ecuacion Ht-D*(Hrr+Hr/r)=0
%D=K/s=conductividad hidraulica/almacenamiento especifico=difusividad
%hidraulica
%h(r0=0.5,t)=ha=hp=34 %condicion dirichlet
%h(rf=20,t)=hb=h0=45 %condicion dirichlet
%h(r,t=0)=h0=45 %condicion inicial
%introducimos los datos
r0=0.5; rf=20; D=10^(-2); hr=0.1; ht=hr; %dt=ht dr=hr
t0=0; tf=2;
r=(r0+hr):hr:(rf-hr);
%Introducimos las funciones
h1=@(r)0*r+1; %coeficiente de Hr
g=@(r,t)0*r*t; %coeficiente independiente
ha=@(t)0*t+34; %funcion de la condicion dirichlet
hb=@(t)0*t+45; %funcion de la condicion dirichlet
h0=@(r)0*r+45; %funcion de la condicion inicial
%como tenemos h,calculamos N
N=round((rf-r0)/hr);
%definimos el vector posicion
x=linspace(r0,rf,N+1);
x=x';
xx=x(2:end-1);
%Transformamos lo que tenemos en un PVI forma H'=-D*(K+B)*H+G; H(0)=h0 de la C.I;
%Calculamos B
H1=h1(xx);
B=(diag(H1(1:end-1),1)-diag(H1(2:end),-1))./(2*hr);
B=diag(1./r).*B;
%calculamos K
K=(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1))/hr^2;
%calculamos valor inicial H(0)=ho
H0=h0(xx);
%calculamos el termino independiente G como PVI
%calculamos un vector de tiempos
M=round((tf-t0)/ht);
%calculamos el vector de tiempos
t=linspace(t0,tf,M+1);
%Generamos una matriz de 0
y=zeros(N-1,M+1);
y(:,1)=H0;
%aplicamos el metodo de heun
for i=1:M       
    G=g(xx,t(i));  
    G(1)=G(1)+D*ha(t(i))/hr^2;          %dirichlet
    G(end)=G(end)+D*hb(t(i))/hr^2;      %dirichlet
    
    Gi=g(xx,t(i)+hr);       
    Gi(1)=Gi(1)+D*ha(t(i))/hr^2;          %dirichlet
    Gi(end)=Gi(end)+D*hb(t(i))/hr^2;      %dirichlet
    
    k1=-D*(K+B)*y(:,i)+G;
    k2=-D*(K+B)*y(:,i)+Gi+k1*ht;
    y(:,i+1)=y(:,i)+ht*(k1+k2)/2;
end
%definimos las condiciones dirichlet
Ha=ha(t);
Hb=hb(t);
%formulamos la solucion
Y=[Ha;y;Hb];
%hacemos el grafico de la solucion en 3d
%mallamos 
[Mt,Mx]=meshgrid(t,x);
surf(Mx,Mt,Y)
shading flat
title('Evolucion del nivel piezometrico por Heun')
xlabel('Eje Posicion');ylabel('Eje Tiempo');
zlabel('Nivel piezomatrizo');
Gráfica solución por diferencias finitas

4 Gráfica del comportamiento del nivel piezométrico

En este caso tendremos este sistema: \[\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)=34 \quad \quad h(20,t)=45 & \\ h(\rho,0)=45 & \end{matrix}\right.\]

Siendo [math] D= 0.01[/math] y [math]\rho _{0} = 0.5 [/math] 

%introducimos los datos
r0=0.5; rf=20; D=10^(-2); hr=0.1; ht=hr;
t0=0; tf=2;
r=(r0+hr):hr:(rf-hr);
%Introducimos las funciones
h1=@(r)0*r+1; %coeficiente de Hr
g=@(r,t)0*r*t; %coeficiente independiente
ha=@(t)0*t+34; %funcion de la condicion dirichlet
hb=@(t)0*t+45; %funcion de la condicion dirichlet
h0=@(r)0*r+45; %funcion de la condicion inicial
%como tenemos h,calculamos N
N=round((rf-r0)/hr);
%definimos el vector posicion
x=linspace(r0,rf,N+1);
x=x';
xx=x(2:end-1);
%Transformamos lo que tenemos en un PVI forma H'=-D*(K+B)*H+G; H(0)=h0 de la C.I;
%Calculamos B
H1=h1(xx);
B=(diag(H1(1:end-1),1)-diag(H1(2:end),-1))./(2*hr);
B=diag(1./r).*B;
%calculamos K
K=(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1))/hr^2;
%calculamos valor inicial H(0)=ho
H0=h0(xx);
%calculamos el termino independiente G como PVI
%calculamos un vector de tiempos
M=round((tf-t0)/ht);
%calculamos el vector de tiempos
t=linspace(t0,tf,M+1);
%Generamos una matriz de 0
y=zeros(N-1,M+1);
y(:,1)=H0;
%aplicamos el metodo de heun
for i=1:M       
    G=g(xx,t(i));  
    G(1)=G(1)+D*ha(t(i))/hr^2;          %dirichlet
    G(end)=G(end)+D*hb(t(i))/hr^2;      %dirichlet
    
    Gi=g(xx,t(i)+hr);       
    Gi(1)=Gi(1)+D*ha(t(i))/hr^2;          %dirichlet
    Gi(end)=Gi(end)+D*hb(t(i))/hr^2;      %dirichlet
    
    k1=-D*(K+B)*y(:,i)+G;
    k2=-D*(K+B)*y(:,i)+Gi+k1*ht;
    y(:,i+1)=y(:,i)+ht*(k1+k2)/2;
end
%definimos las condiciones dirichlet
Ha=ha(t);
Hb=hb(t);
%formulamos la solucion
Y=[Ha;y;Hb];
%hacemos el grafico de la solucion en 2d
s=2;
s=round((s-r0)/hr+1);
z=Y(s,:);
plot(t,z)
title('Comportamiento del nivel piezometrico para r=2')
ylabel('Nivel piezometrico')
xlabel('Eje Tiempo')
Gráfica del comportamiento del nivel piezométrico

5 Resolución por método del trapecio y el método de Euler implícito

De nuevo, el sistema será: \[\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)=34 \quad \quad h(20,t)=45 & \\ h(\rho,0)=45 & \end{matrix}\right.\]

Siendo [math] D= 0.01[/math] y [math]\rho _{0} = 0.5 [/math] 

%trabajo Ecus Pregunta 5
%ecuacion Ht-D*(Hrr+Hr/r)=0
%D=K/s=conductividad hidraulica/almacenamiento especifico=difusividad
%hidraulica
%h(r0=0.5,t)=ha=hp=34 %condicion dirichlet
%h(rf=20,t)=hb=h0=45 %condicion dirichlet
%h(r,t=0)=h0=45 %condicion inicial
%introducimos los datos
r0=0.5; rf=20; D=10^(-2); hr=0.1; ht=hr;
t0=0; tf=2;
r=(r0+hr):hr:(rf-hr);
%Introducimos las funciones
h1=@(r)0*r+1; %coeficiente de Hr
g=@(r,t)0*r*t; %coeficiente independiente
ha=@(t)0*t+34; %funcion de la condicion dirichlet
hb=@(t)0*t+45; %funcion de la condicion dirichlet
h0=@(r)0*r+45; %funcion de la condicion inicial
%como tenemos h,calculamos N
N=round((rf-r0)/hr);
%definimos el vector posicion
x=linspace(r0,rf,N+1);
x=x';
xx=x(2:end-1);
%Transformamos lo que tenemos en un PVI forma H'=-D*(K+B)*H+G; H(0)=h0 de la C.I;
%Calculamos B
H1=h1(xx);
B=(diag(H1(1:end-1),1)-diag(H1(2:end),-1))./(2*hr);
B=diag(1./r).*B;
%calculamos K
K=(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1))/hr^2;
%calculamos valor inicial H(0)=ho
H0=h0(xx);
%calculamos el termino independiente G como PVI
%calculamos un vector de tiempos
M=round((tf-t0)/ht);
%calculamos el vector de tiempos
t=linspace(t0,tf,M+1);
%Generamos la matriz de 0
y=zeros(N-1,M+1);
y(:,1)=H0;
%Generamos la matriz de 0
y1=zeros(N-1,M+1);
y1(:,1)=H0;
%aplicamos el metodo de Euler implicito
for i=1:M       
    G=g(xx,t(i+1));        
    G(1)=G(1)+D*ha(t(i+1))/hr^2;          %dirichlet
    G(end)=G(end)+D*hb(t(i+1))/hr^2;      %dirichlet
    
    y(:,i+1)=(eye(size(K))+ht*D*(B+K))\(y(:,i)+ht*G);    
end
%aplicamos el metodo del Trapecio
for i=1:M     
    G=g(xx,t(i));  
    G(1)=G(1)+D*ha(t(i))/hr^2;          %dirichlet
    G(end)=G(end)+D*hb(t(i))/hr^2;      %dirichlet
    
    Gi=g(xx,t(i+1));       
    Gi(1)=Gi(1)+D*ha(t(i+1))/hr^2;          %dirichlet
    Gi(end)=Gi(end)+D*hb(t(i+1))/hr^2;      %dirichlet
    
    y1(:,i+1)=(eye(size(K))+ht*D*(B+K)/2)\(y1(:,i)+ht/2*(Gi-D*(B+K)*y1(:,i)+G));
end
%definimos las condiciones dirichlet
Ha=ha(t);
Hb=hb(t);
%formulamos la solucion Euler Implicito
Y=[Ha;y;Hb];
%formulamos la solucion Trapecio
Y1=[Ha;y1;Hb];
%hacemos el grafico de la solucion en 3d
%mallamos Euler Implicito
[Mt,Mx]=meshgrid(t,x);
surf(Mx,Mt,Y)
shading flat
xlabel('Eje Posicion');ylabel('Eje Tiempo');
zlabel('Nivel piezomatrizo');
title('Evolucion del nivel piezometrico por Euler Implicito');
%mallamos Trapecio
[Mt,Mx]=meshgrid(t,x);
figure
surf(Mx,Mt,Y1)
shading flat
xlabel('Eje Posicion');ylabel('Eje Tiempo');
zlabel('Nivel piezomatrizo');
title('Evolucion del nivel piezometrico por Trapecio');
Gráfica solución según método de Euler implícito
Gráfica solución según método del trapecio

6 Gráfica solución

Sistema: \[\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)=34 \quad \quad h(20,t)=45 & \\ h(\rho,0)=45 & \end{matrix}\right.\]

Siendo [math] D= 0.01[/math] y [math]\rho _{0} = 0.5 [/math] 

%trabajo Ecus Pregunta 6
%ecuacion Ht-D*(Hrr+Hr/r)=0
%D=K/s=conductividad hidraulica/almacenamiento especifico=difusividad
%hidraulica
%h(r0=0.5,t)=ha=hp=34 %condicion dirichlet
%h(rf=20,t)=hb=h0=45 %condicion dirichlet
%h(r,t=0)=h0=45 %condicion inicial
%introducimos los datos
r0=0.5; rf=20; D=10^(-2); hr=0.1; ht=0.1;
t0=0; tf=5000;
r=(r0+hr):hr:(rf-hr);
%Introducimos las funciones
h1=@(r)0*r+1; %coeficiente de Hr
g=@(r,t)0*r*t; %coeficiente independiente
ha=@(t)0*t+34; %funcion de la condicion dirichlet
hb=@(t)0*t+45; %funcion de la condicion dirichlet
h0=@(r)0*r+45; %funcion de la condicion inicial
%como tenemos h,calculamos N
N=round((rf-r0)/hr);
%definimos el vector posicion
x=linspace(r0,rf,N+1);
x=x';
xx=x(2:end-1);
%Transformamos lo que tenemos en un PVI forma H'=-D*(K+B)*H+G; H(0)=h0 de la C.I;
%Calculamos B
H1=h1(xx);
B=(diag(H1(1:end-1),1)-diag(H1(2:end),-1))./(2*hr);
B=diag(1./r).*B;
%calculamos K
K=(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1))/hr^2;
%calculamos valor inicial H(0)=ho
H0=h0(xx);
%calculamos el termino independiente G como PVI
%calculamos un vector de tiempos
M=round((tf-t0)/ht);
%calculamos el vector de tiempos
t=linspace(t0,tf,M+1);
%Generamos una matriz de 0
y=zeros(N-1,M+1);
y(:,1)=H0;
%aplicamos el metodo de heun
for i=1:M       
    G=g(xx,t(i));  
    G(1)=G(1)+D*ha(t(i))/hr^2;          %dirichlet
    G(end)=G(end)+D*hb(t(i))/hr^2;      %dirichlet
    
    Gi=g(xx,t(i)+hr);       
    Gi(1)=Gi(1)+D*ha(t(i))/hr^2;          %dirichlet
    Gi(end)=Gi(end)+D*hb(t(i))/hr^2;      %dirichlet
    
    k1=-D*(K+B)*y(:,i)+G;
    k2=-D*(K+B)*y(:,i)+Gi+k1*ht;
    y(:,i+1)=y(:,i)+ht*(k1+k2)/2;
end
%definimos las condiciones dirichlet
Ha=ha(t);
Hb=hb(t);
%formulamos la solucion
Y=[Ha;y;Hb];
%hacemos el grafico de la solucion en 3d
%mallamos 
[Mt,Mx]=meshgrid(t,x);
surf(Mx,Mt,Y)
title('Evolucion del nivel piezometrico a lo largo de 5000H')
axis tight
shading flat
xlabel('Eje Posicion');ylabel('Eje Tiempo');
zlabel('Nivel piezomatrizo');
Gráfica solución

7 Estado estacionario

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

Siendo [math] D= 0.01[/math] y [math]\rho _{0} = 0.5 [/math] 

%cogemos el t=0, 100, 1000, 5000 del apartado 6
figure
s=0;
s=round((s-t0)/ht+1);
t00=Y(:,s)';
plot(x,t00)
hold on
s=100;
s=round((s-t0)/ht+1);
t00=Y(:,s)';
plot(x,t00)
s=1000;
s=round((s-t0)/ht+1);
t00=Y(:,s)';
plot(x,t00)
s=5000;
s=round((s-t0)/ht+1);
t00=Y(:,s)';
plot(x,t00)


%Problema 7
%ecuacion Hrr+Hr/r=0
%D=K/s=conductividad hidraulica/almacenamiento especifico=difusividad
%hidraulica
%h(r0=0.5)=ha=hp=34 %condicion dirichlet
%h(rf=20)=hb=h0=45 %condicion dirichlet
ha=34; hb=45;
%introducimos los datos
r0=0.5; rf=20; D=10^(-2); hr=0.01; ht=0.1; %dt=ht dr=hr
t0=0; tf=5;
r=(r0+hr):hr:(rf-hr);
%Introducimos las funciones
h1=@(r)0*r+1; %coeficiente de Hr
g=@(r)0*r; %coeficiente independiente
%calculamos el vector de tiempos
t=t0:ht:tf;
%como tenemos h,calculamos N
N=round((rf-r0)/hr);
%definimos el vector posicion
x=linspace(r0,rf,N+1);
x=x';
xx=x(2:end-1);
%Transformamos lo que tenemos en un PC forma (K+B)*H=G
%Calculamos B
H1=h1(xx);
B=(diag(H1(1:end-1),1)-diag(H1(2:end),-1))/(2*hr);
B=diag(1./r)*B;
%calculamos K
K=(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1))/hr^2;
%calculamos G
G=g(xx);
G(1)=G(1)+ha/hr^2+ha/(2*hr*r(1)); %condicion dirichlet
G(end)=G(end)+hb/hr^2-hb/(2*hr*r(end)); %condicion dirichlet
%calculamos H=(K+B)\G y añadimos las condiciones dirichlet
H=(K+B)\G;
H=[ha;H;hb];
%generamos una matriz de 0
y=zeros(length(t),length(H));
y(1,:)=H;
%transformamos H en una matriz
for i=1:length(t)-1
    H=(K-B)\G;
    H=[ha;H;hb];
    y(i+1,:)=H;
end
%dibujamos el valor estacionario
plot(x,y(end,:))
title('cambio de nivel piezometrico')
legend('t en 0','t en 100','t en 1000','t en 5000','Sol estacionaria')
xlabel('Eje Posicion');ylabel('Nivel piezomatrizo');
hold off
Gráfica Nivel piezométrico

8 Capacidad de recuperación del acuífero

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

Siendo [math] D= 0.01[/math] y [math]\rho _{0} = 0.5 [/math] 

%trabajo Ecus Preguntas 7+8
%ecuacion Ht-D*(Hrr+Hr/r)=0
%D=K/s=conductividad hidraulica/almacenamiento especifico=difusividad
%hidraulica
%hr(r0=0.5,t)=ha=hp=34 %condicion neumann
%h(rf=20,t)=hb=h0=45 %condicion dirichlet
%h(r,t=0)=y(1,:) %condicion inicial, procede de P7

%Problema 7
%ecuacion Hrr+Hr/r=0
%D=K/s=conductividad hidraulica/almacenamiento especifico=difusividad
%hidraulica
%h(r0=0.5)=ha=hp=34 %condicion dirichlet
%h(rf=20)=hb=h0=45 %condicion dirichlet
ha=34; hb=45;
%introducimos los datos
r0=0.5; rf=20; D=10^(-2); hr=0.1; ht=100; 
t0=0; tf=5000;
r=(r0+hr):hr:(rf-hr);
%Introducimos las funciones
h1=@(r)0*r+1; %coeficiente de Hr
g=@(r)0*r; %coeficiente independiente
%calculamos el vector de tiempos
t=t0:ht:tf;
%como tenemos h,calculamos N
N=round((rf-r0)/hr);
%definimos el vector posicion
x=linspace(r0,rf,N+1);
x=x';
xx=x(2:end-1);
%Transformamos lo que tenemos en un PC forma (K+B)*H=G
%Calculamos B
H1=h1(xx);
B=(diag(H1(1:end-1),1)-diag(H1(2:end),-1))/(2*hr);
B=diag(1./r)*B;
%calculamos K
K=(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1))/hr^2;
%calculamos G
G=g(xx);
G(1)=G(1)+ha/hr^2+ha/(2*hr*r(1)); %condicion dirichlet
G(end)=G(end)+hb/hr^2-hb/(2*hr*r(end)); %condicion dirichlet
%calculamos H=(K+B)\G y añadimos las condiciones dirichlet
H=(K+B)\G;
H=[ha;H;hb];
%generamos una matriz de 0
y=zeros(length(t),length(H));
y(1,:)=H;
%transformamos H en una matriz
for i=1:length(t)-1
    H=(K-B)\G;
    H=[ha;H;hb];
    y(i+1,:)=H;
end

%Pregunta 8
%hr(r0=0.5,t)=hra=hp=0 %condicion neumann
%h(rf=20,t)=hb=h0=45 %condicion dirichlet
%h(r,t=0)=y(1,:) %condicion inicial, procede de P7
hr=0.1; ht=0.1;
%introducimos los datos
r=r0:hr:(rf-hr);
%Introducimos las funciones
g=@(r,t)0*r; %coeficiente independiente
hra=@(t)0*t; %funcion de la condicion neumann
hb=@(t)0*t+45; %funcion de la condicion dirichlet
H0=y(end,:);
%como tenemos h,calculamos N
N=round((rf-r0)/hr);
%definimos el vector posicion
x=linspace(r0,rf,N+1);
x=x';
xx=x(1:end-1);
l=length(xx);
%Transformamos lo que tenemos en un PVI forma H'=-D*(K+B)*H+G; H(0)=h0 de la C.I;
%Calculamos B
H1=h1(xx);
B=(diag(H1(1:end-1),1)-diag(H1(2:end),-1))/(2*hr);
B=diag(1./r)*B;
B(1,2)=0;%condición neumann
%calculamos K
K=(2*diag(ones(1,N))-diag(ones(1,N-1),-1)-diag(ones(1,N-1),1))/hr^2;
K(1,2)=-2/hr^2;
%Valor inicial
H0=H0(1:N);
%calculamos el termino independiente G como PVI
%calculamos un vector de tiempos
M=round((tf-t0)/ht);
%calculamos el vector de tiempos
t=linspace(t0,tf,M+1);
%Generamos una matriz de 0
y=zeros(N,M+1);
y(:,1)=H0;
%aplicamos el metodo de heun
for i=1:M       
    G=g(xx,t(i));  
    G(1)=G(1)+2*D*hra(t(i))/hr;          %neumann
    G(end)=G(end)+D*hb(t(i))/hr^2;      %dirichlet
    
    Gi=g(xx,t(i)+hr);       
    Gi(1)=Gi(1)+2*D*hra(t(i))/hr;          %neumann
    Gi(end)=Gi(end)+D*hb(t(i))/hr^2;      %dirichlet
    
    k1=-D*(K+B)*y(:,i)+G;
    k2=-D*(K+B)*y(:,i)+Gi+k1*ht;
    y(:,i+1)=y(:,i)+ht*(k1+k2)/2;
end
%definimos las condiciones dirichlet
Hb=hb(t);
%formulamos la solucion
Y=[y;Hb];
%hacemos el grafico de la solucion en 3d
%mallamos 
[Mt,Mx]=meshgrid(t,x);
surf(Mx,Mt,Y)
shading flat
axis tight
title('Nivel piezométrico en función del tiempo y la distancia');
xlabel('Eje Posicion');ylabel('Eje Tiempo');
zlabel('Nivel piezométrico');
Gráfica de recuperación

9 Aproximación de la solución

El nuevo sistema será: \[\left\{\begin{matrix}\ \frac{\partial ^2 h}{\partial \rho^2}+ \frac{1}{\rho}·\frac{\partial h}{\partial \rho}=0 , & \\ \frac{\partial h}{\partial \rho}(\rho _{0},t)=0 \quad \quad h(20,t)=40 & \\ h(\rho,0)=log(|r+0.01|)/20.01 & \end{matrix}\right.\]

Siendo [math] D= 0.01[/math] y [math]\rho _{0} = 0.5 [/math] 

%trabajo Ecus Pregunta 9
%ecuacion Ht-D*(Hrr+Hr/r)=0
%D=K/s=conductividad hidraulica/almacenamiento especifico=difusividad
%hidraulica
%h'(r0=0.5,t)=ha=hp=0 %condicion neumann
%h(rf=20,t)=hb=h0=40 %condicion dirichlet
%h(r,t=0)=h0=log(abs(r+0.01)/20.01) %condicion inicial
%introducimos los datos
r0=0.5; rf=20; D=10^(-2); hr=0.1; ht=hr; %dt=ht dr=hr
t0=0; tf=2;
r=(r0):hr:(rf-hr);
%Introducimos las funciones
h1=@(r)0*r+1; %coeficiente de Hr
g=@(r,t)0*r; %coeficiente independiente
ha=@(t)0*t; %funcion de la condicion Neumann
hb=@(t)0*t+40; %funcion de la condicion dirichlet
h0=@(r)log(abs(r+0.01)/20.01); %funcion de la condicion inicial
%como tenemos h,calculamos N
N=round((rf-r0)/hr);
%definimos el vector posicion
x=linspace(r0,rf,N+1);
x=x';
xx=x(1:end-1);
%Transformamos lo que tenemos en un PVI forma H'=-D*(K+B)*H+G; H(0)=h0 de la C.I;
%Calculamos B
H1=h1(xx);
B=(diag(H1(1:end-1),1)-diag(H1(2:end),-1))./(2*hr);
B=diag(1./r)*B;
B(1,2)=0;
%calculamos K
K=(2*diag(ones(1,N))-diag(ones(1,N-1),-1)-diag(ones(1,N-1),1))/hr^2;
K(1,2)=-2/hr^2;%condicion Neymann
%calculamos valor inicial H(0)=ho
H0=h0(xx);
%calculamos el termino independiente G como PVI
%calculamos un vector de tiempos
M=round((tf-t0)/ht);
%calculamos el vector de tiempos
t=linspace(t0,tf,M+1);
%Generamos una matriz de 0
y=zeros(N,M+1);
y(:,1)=H0;
%aplicamos el metodo de heun
for i=1:M       
    G=g(xx,t(i));  
    G(1)=G(1)+2*D*ha(t(i))/hr;          %neumann
    G(end)=G(end)+D*hb(t(i))/hr^2;      %dirichlet
    
    Gi=g(xx,t(i)+hr);       
    Gi(1)=Gi(1)+2*D*ha(t(i))/hr;          %neumann
    Gi(end)=Gi(end)+D*hb(t(i))/hr^2;      %dirichlet
    
    k1=-D*(K+B)*y(:,i)+G;
    k2=-D*(K+B)*y(:,i)+Gi+k1*ht;
    y(:,i+1)=y(:,i)+ht*(k1+k2)/2;
end
%definimos la condicion dirichlet
Hb=hb(t);
%formulamos la solucion
Y=[y;Hb];
%hacemos el grafico de la solucion en 3d
%mallamos 
[Mt,Mx]=meshgrid(t,x);
surf(Mx,Mt,Y)
shading flat
axis tight
title('Nuevo acuifero')
xlabel('Eje Posicion');ylabel('Eje Tiempo');
zlabel('Nivel piezomatrizo');
Gráfico de la aproximación de la solución
Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Nivel_piezométrico_en_un_acuífero_confinado._Grupo_8-B&oldid=39580»