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

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.

Archivo:Https://mat.caminos.upm.es/w/images/f/fc/Dibujol.JPG

Aquí se representa un acuífero con nivel piezométrico constante, de superficie horizontal en el que se instala un pozo que extrae un caudal S, esto hace que baje el nivel del acuífero hasta que este se vuelva a mantener constante.

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 divergenia de este y así obtenemos el Laplaciano: : [math] (\frac{\partial ^2 h}{\partial \rho^2}+ \frac{1}{\rho}·\frac{\partial h}{\partial \rho}+\frac{\partial^2 h}{\partial \theta^2})= 0[/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');