Diferencia entre revisiones de «Nivel piezométrico - Grupo 19»
(→Gráfica del nivel piezométrico para tiempos grandes) |
(→Gráfica del nivel piezométrico para tiempos grandes) |
||
| Línea 143: | Línea 143: | ||
surf(RR,TT,sol) | surf(RR,TT,sol) | ||
}} | }} | ||
| + | Gráfica: | ||
==Valor estacionario== | ==Valor estacionario== | ||
Revisión del 16:46 19 may 2014
| 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).
Ley de Darcy particularizada para nuestro problema: :[math]{q}=-{K}\cdot\nabla(h(ρ))[/math]
Ecuación de la conservación de la masa: [math] S ·\frac{ \partial h }{ \partial t } + div q = 0[/math]
Contenido
- 1 Cálculo del Laplaciano y ecuación diferencial en coordenadas polares
- 2 Condiciones de contorno y sistema de ecuaciones
- 3 Resolución del problema
- 4 Gráfica del comportamiento del nivel piezométrico
- 5 Resolución del problema
- 6 Gráfica del nivel piezométrico para tiempos grandes
- 7 Valor estacionario
- 8 Estimación de la capacidad de recuperación de acuífero
- 9 Método de Fourier
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
2 Condiciones de contorno y sistema de ecuaciones
3 Resolución del problema
3.1 Método de las diferencias finitas
3.2 Método del trapecio
4 Gráfica del comportamiento del nivel piezométrico
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)Grafica:
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:
6 Gráfica del nivel piezométrico para tiempos grandes
%Datos
r0=1;rN=20;D=10^(-2);T=5000;
%Dsicretizacion(para tiempos grandes, cambian los pasos de tiempo y el
%t final)
dt=200;dr=0.1;
%Discretizacion
N=(rN-r0)/dr;
%Vectores radio y tiempo
r=(r0+dr):dr:(rN-dr);
t0=0;tN=T;
t=t0:dt:T;
%Preparacion del método
H=ones(1,length(r));H=40*H;H=H';
h0=35;hN=40;
sol(1,:)=[h0 H' hN];
%Interaciones
for k=1:length(t)-1
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;
H=(eye(N-1)+((K-L1).*(dt/2)))\(H+((dt/2).*(-(K-L1)*H+F+F)));
sol(k+1,:)=[h0 H' hN];
end
%Dibujamos
rf=r0:dr:rN;
[RR,TT]=meshgrid(rf,t);
surf(RR,TT,sol)Gráfica:
