Nivel piezométrico - Grupo 19

Trabajo realizado por estudiantes
Título	Nivel piezométrico. Grupo 19-B
Asignatura	Ecuaciones Diferenciales
Curso	Curso 2013-14
Autores	Javier Abad, Jesús Castaño, Ignacio Embid, Javier Pérez
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Se define el concepto de nivel piezométrico como 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. Para poder conocer la variación del nivel piezométrico se utiliza la ecuación de la conservación de la masa y la ley de Darcy.

La Ley de Darcy describe la relación entre la cantidad o la velocidad de flujo del agua, la permeabilidad del acuífero y el gradiente piezométrico (o gradiente hidráulico).

Ecuación de la conservación de la masa: [math] S ·\frac{ \partial h }{ \partial t } + div q = 0[/math]

Ley de Darcy particularizada para nuestro problema: :[math]{q}=-{K}\cdot\nabla(h(ρ))[/math]

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

Cálculo del Laplaciano y ecuación diferencial en coordenadas polares: Divergencia: [math] \nabla\cdot\vec F = \frac{\partial F_x}{\partial x}+ \frac{\partial F_y}{\partial y}+ \frac{\partial F_z}{\partial z} [/math] Si operamos las dos ecuaciones iniciales, tenemos [math] \frac{ \partial h }{ \partial t } - D · \Delta h = 0 [/math]

5 Resolución del problema

5.1 Método de Éuler explícito

%Datos
L=20;T=50;D=0.01;hp=35;ho=40;
%Discretizacion
dx=1/10;dt=dx;
%Vector x
x0=1;
x=x0:dx:L;
xint=x0+dx:dx:L-dx;
%Vector tiempo
t=0:dt:T;
%Aplicar euler explicito
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));
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
figure(1)
[xx,tt]=meshgrid(x,t); 
surf(xx,tt,sol)

5.2 Método de Éuler implícito

%Datos
r0=1;rN=20;D=10^(-2);T=2;
%Discretizacion
dt=0.1;dr=dt;
N=(rN-r0)/dr;
%Vectores radio y tiempo
r=(r0+dr):dr:(rN-dr);
t0=0;tN=T;
t=t0:dt:T;
%Preparacion para resolucion
H=ones(1,length(r));H=40*H;H=H';
h0=35;hN=40;
res(1,:)=[h0 H' hN];
%Aplicamos euler implicito
F=zeros(length(r),1);F(1)=D*((h0/(dr^2))-((1/r(1)))*(h0/(2*dr)));
F(N-1)=D*((hN/(dr^2))+((1/r(N-1)))*(hN/(2*dr)));
K=2.*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
K=K./(dr^2);K=D*K;
L=diag(zeros(1,N-1))+diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
L=L./(2*dr);L=D*L;
L1=(diag(1./r))*L;
for k=1:length(t)-1
    H=(eye(N-1)+((K-L1).*(dt)))\(H+((dt).*(F)));
    res(k+1,:)=[h0 H' hN];
end
rf=r0:dr:rN;
[RR,TT]=meshgrid(rf,t);
surf(RR,TT,res)

Gráfica: