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.

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] y [math]\Delta t = \Delta \rho [/math] 

%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);
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;
%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;
    G(end)=G(end)+D*hb(t(i))/hr^2;
    k1=-D*(K+B)*y(:,i)+G;
    k2=-D*(K+B)*y(:,i)+k1*ht+G;
    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)
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

%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 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;
    G(end)=G(end)+D*hb(t(i))/hr^2;
    k1=-D*(K+B)*y(:,i)+G;
    k2=-D*(K+B)*y(:,i)+k1*ht+G;
    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
z=Y(2,:);
plot(z,t)
xlabel('Nivel piezometrico')
ylabel('Eje Ttiempo')
Gráfica del comportamiento del nivel piezométrico

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

%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 2 matrice de 0
y=zeros(N-1,M+1);
y(:,1)=H0;
y1=zeros(N-1,M+1);
y1(:,1)=H0;
%aplicamos el metodo de Euler implicito
for i=1:M
    G=g(xx,t(i));
    G(1)=G(1)+D*ha(t(i))/hr^2;
    G(end)=G(end)+D*hb(t(i))/hr^2;
    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;
    G(end)=G(end)+D*hb(t(i))/hr^2;
    y1(:,i+1)=(eye(size(K))+ht*D*(B+K)/2)\(y1(:,i)+ht/2*(-D*(B+K)*y1(:,i)+2*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('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('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

%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;
    G(end)=G(end)+D*hb(t(i))/hr^2;
    k1=-D*(K+B)*y(:,i)+G;
    k2=-D*(K+B)*y(:,i)+k1*ht+G;
    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)
axis tight
shading flat
xlabel('Eje Posicion');ylabel('Eje Tiempo');
zlabel('Nivel piezomatrizo');
Gráfica solución

7 Estado estacionario

8 Capacidad de recuperación del acuífero

9 Aproximación de la solución

%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(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(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*ha(t(i))/hr;%condicion Neymann
    G(end)=G(end)+D*hb(t(i))/hr^2;
    k1=-D*(K+B)*y(:,i)+G;
    k2=-D*(K+B)*y(:,i)+k1*ht+G;
    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
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=39401»