Diferencia entre revisiones de «Nivel piezométrico - Grupo 19»

De MateWiki
Saltar a: navegación, buscar
(→‎Método de las diferencias finitas)
(→‎Método de las diferencias finitas)
Línea 47: Línea 47:
  
 
===Método de las diferencias finitas===
 
===Método de las 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>
 
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>
 
{{matlab|codigo=
 
{{matlab|codigo=

Revisión del 18:48 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]

Acuifero confinado

Representa un manto acuifero horizontal, homogéneo y de espesor constante, confinado entre dos acuiclusos(capas impermeables) y drenado por un pozo de radio r. Sin la existencia del pozo, el nivel piezométrico en el acuífero era constante. Construido el pozo, se extrae de él un caudal S, que hace bajar el nivel del agua en el pozo hasta un nivel donde se mantiene 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

3.1 Método de las 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] 

%Datos
r0=1;rN=20;D=10^(-2);T=2;
%Discretizacion para las dos variables
dt=0.1;dr=dt;
N=(rN-r0)/dr;
%Vector 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];
%Tnteraciones
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)));
    res(k+1,:)=[h0 H' hN];
end
%el problema por el m�etodo de diferencias finitas y aplicando el m�etodo del trapecio en tiempo
rf=r0:dr:rN;
figure(1)
[RR,TT]=meshgrid(rf,t);
surf(RR,TT,res)

Gráfica:

Diferencias finitas

4 Gráfica del comportamiento del nivel piezométrico

Comportmiento

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)
xlabel('Distancia');
ylabel('Tiempo');
zlabel('Nivel piezométrico')

Grafica:

Euler explícito

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)
xlabel('Distancia');
ylabel('Tiempo');
zlabel('Nivel piezométrico')

Gráfica:

Euler implícito

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)
xlabel('Distancia');
ylabel('Tiempo');
zlabel('Nivel piezométrico')

Gráfica:

tiempos grandes

7 Valor estacionario

%Datos
    r0=1;rN=20;D=10^(-2);
    %Discretizacion
    dt=0.1;dr=dt;
    %Vectores tiempo y radio
    N=(rN-r0)/dr;T=2;
    r=(r0+dr):dr:(rN-dr);
    t0=0;tN=T;
    t=t0:dt:T;
    %Inicializacion     
    H=ones(1,length(r));H=40*H;H=H';
    h0=35;hN=40;
    sol(1,:)=[h0 H' hN];
    %Aplicacion
    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;
    %Tnteracciones
    for k=1:length(t)-1
        H=(((K-L1)))\(F);
        sol(k+1,:)=[h0 H' hN];
    end
    %Dibujamos solucion
    rf=r0:dr:rN;
    [RR,TT]=meshgrid(rf,t);
    surf(RR,TT,sol)

Gráfica:

center

8 Estimación de la capacidad de recuperación de acuífero

%Datos
    dt1=100;T1=2000;t01=0;tN1=T;h01=40;hN1=40;
    %Vector tiempo
    t1=t01:dt1:T1;
    %Inicializacion
    P=sol(21,:);
    P(1)=[];
    P(19)=[];
    P=P';
    H1=[h0 P' hN];
    sol1(1,:)=H1';
    F1=zeros(length(r),1);F1(1)=D*((h01/(dr^2))-((1/r(1)))*(h01/(2*dr)));
    F1(N-1)=D*((hN1/(dr^2))+((1/r(N-1)))*(hN1/(2*dr)));
    %Interacciones
    for k=1:length(t1)-1
        P=(eye(N-1)+((K-L1).*(dt1/2)))\(P+((dt1/2).*(-(K-L1)*P+F1+F1)));
        sol1(k+1,:)=[h0 P' hN];
    end
    [RR1,TT1]=meshgrid(rf,t1);
    surf(RR1,TT1,sol1)

Gráfica:

9 Método de Fourier

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Nivel_piezométrico_-_Grupo_19&oldid=12883»