Nivel piezométrico - Grupo 19

De MateWiki
Saltar a: navegación, buscar
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

Como se puede comprobar en la gráfica, al tener un intervalo tan pequeño, es complicado apreciar las variaciones en el nivel piezométrico con respecto a la variable independiente que es el tiempo.

4 Gráfica del comportamiento del nivel piezométrico

%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
rf=r0:dr:rN;


x=res(:,2);
plot(x,t)
xlabel('Nivel Piezométrico')
ylabel('tiempo')



Comportmiento

Grafica que representa el comportamiento en los puntos para los cuales [math] ρ=2 [/math] en una gráfica bidimensional. Para representarlo, hemos tomado un intervalo de tiempo entre 0 y 2 y la variacion de nivel resultante queda entre los valores 37 y 40.

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

En esta solución se puede apreciar un color más oscuro en la gráfica. Esto es debido a que se nos ha impuesto un paso de tiempo muy pequeño comparado con el intervalo total de tiempo aumentando el mallado, que se representa en negro, oscureciendo la tonalidad final de la gráfica.

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


Al ser el método de euler implícito más cercano a la realidad por tomar el valor en el punto estudiado, es comprensible suponer que ésta gráfica se adecuará más al fenómeno físico estudiado.

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

Tras obtener este resultado del código, advertimos que este resultado no puede darse ya que en el radio 0 no debería cambiar el nivel piezométrico con respecto al tiempo y en una situacion dónde no hay un aporte puntual de agua debería ser estable con respecto al tiempo y al radio, manteniéndose constante a medida que se aleja del pozo resultando casi imperceptible para radios 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 resultante


Una vez empezamos la extracción del agua del pozo, el nivel del agua va disminuyendo en éste, provocando un cambio de presiones intersticiales en los alrededores que provocarán un flujo de agua de las partes mas alejadas del pozo hacia este. Al cabo de un determinado tiempo t, podemos encontrar una situación de equilibrio dónde el agua aportada por el sistema es igual al agua extraída por bombeo en el pozo. En esta situación, podemos decir que el nivel piezométrico no varia con respecto al tiempo, por lo tanto la derivada con respecto al tiempo de la función del nivel piezométrico debe estar igualada a cero, haciendo que el nivel piezométrico permanezca constante y haciéndola variar únicamente con el radio. Haciéndo las simplificaciones explicadas anteriormente, la ecuacion de modelizacion resultante es: :

[math] D·(\frac{\partial ^2 h}{\partial \rho^2}+ \frac{1}{\rho}·\frac{\partial h}{\partial \rho})=0[/math]

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:

Capacidad de recuperación

Para el estudio de la capacidad de recuperacion de un acuifero, hay que, partiendo de una situación dónde hay producido un desequilibrio de presiones instersiticiales como ha sido el caso del apartado anterior, impermeabilizar las paredes del pozo o, en su defecto, detener el bombeo del agua del acuífero. Para ello, la forma más fácil de representarlo en nuestro sistema de ecuaciones es decir que el el flujo para r=r0 es nulo. En este caso, dónde no está siendo drenado, el acuífero tiende a recuperar un nivel piezométrico constante en todos los puntos del sistema. Efectivamente, el valor estacionario tiende a este nivel piezométrico en todos los puntos. Es decir, tiende a recuperar la posición inicial.

9 Sistema de Fourier

Suponemos que el radio del pozo es despreciable y que una vez extraido el agua el pozo se cierra quedando el acuifero de nuevo homogéneo. Tomamos la región [math] ρ\lt20 [/math] . Nuestro objetivo es aproximar las soluciones usando el método de Fourier. Suponemos que la solución sólo depende de [math]ρ[/math] y [math]t[/math], y tomaremos ahora como condición frontera [math] h(20,t)=40 [/math]

Tomamos como problema inicial : : \[\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)=f(\rho) \quad \quad h(20,t)=40, & \\ h_{\rho}(0,t)=0 & \end{matrix}\right.\]

Operamos en el sistema anterior para hacerlo homogéneo ([math]{h}(20,t)=0 [/math] y [math]{h}(\rho ,0)=f(\rho)-4 [/math]), tomamos [math]{h}(ρ,t)=X(ρ).T(t) [/math] y resolviendo por el método de variables separadas obtenemos: : \[\left\{\begin{matrix}\ X′′(ρ)+ \frac{X′(ρ)}{ρ} + λ X(ρ)=0 & \\ X′(0)=0 \quad \quad X(20)=0 & \\ & \end{matrix}\right.\]

Este problema se asimila a uno de Bessel de primera especie cuyas autofunciones son del tipo [math]Ψ(\rho)=J_0·(\rho\sqrt{λ})[/math] y si tomamos [math] ρ=20 [/math] obtenemos los 5 primeros autovalores del problema asociado: :

[math] λ_{1}=0,0145 \quad λ_{2}=0,0762 \quad λ_{3}=0,1872 \quad λ_{4}=0,3476 \quad λ_{5}=0,5573 [/math] 

%Introducimos los autovalores del sistema, que hemos calculado, tomándolos
%de menor a mayor. 
x1=0.0145;
x2=0.0762;
x3=0.1872;
x4=0.3476;
x5=0.5575;
%Vector con r
r=0:0.01:20;
%Aplicamos la funcion besselj
z1=besselj(0,r.*(sqrt(x1)));
z2=besselj(0,r.*(sqrt(x2)));
z3=besselj(0,r.*(sqrt(x3)));
z4=besselj(0,r.*(sqrt(x4)));
z5=besselj(0,r.*(sqrt(x5)));
%Dibujamos las rectas obtenidas en la gráfica
hold on
plot(r,z1)
plot(r,z2)
plot(r,z3)
plot(r,z4)
plot(r,z5)
hold off


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