Nivel piezométrico G5

De MateWiki
Revisión del 23:21 19 may 2014 de Emilio Valero (Discusión | contribuciones) (→‎Capacidad de recuperación en acuífero)

Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Nivel piezométrico G5
Asignatura Ecuaciones Diferenciales
Curso Curso 2013-14
Autores Francisco Durán Muñoz, Javier Bosch Martínez, Manuel Umbert Martín, Miguel Ángel García García, Emilio Valero Muñoz-Rojas,
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura

Definimos nivel piezométrico como la altura que alcanzaría el agua al realizar un sondeo en un punto de un acuífero confinado. Este valor depende de la presión a la que esté el propio acuifero.

Si construimos sobre el acuífero confinado un pozo circular de radio [math]\rho _{0}[/math], el nivel piezométrico varía. Para que el problema sea más sencillo utilizaremos coordenadas cilíndricas. Para poder conocer la variación del nivel piezométrico nos apoyaremos en la ecuación de la conservación de la masa y la ley de Darcy: 

[math] S ·\frac{ \partial h }{ \partial t } + div q = 0[/math]

[math] q = - k ·\nabla h [/math]


La ley de Darcy establece que "el flujo de agua q a través de un medio poroso es proporcional a la diferencia de presión, que a su vez se puede escribir en términos del gradiente del nivel piezométrico en cada punto". La constante K se deduce experimentalmente para cada material y se conoce como la conductividad hidráulica o permeabilidad. La constante S en la ley de conservación de la masa se conoce como almacenamiento específico y se interpreta como la cantidad de agua que libera el acuífero al descender el nivel piezométrico en una unidad, por unidad de volumen. Combinando las ecuaciones de conservación de la masa con la ley de Darcy, obtenemos la ecuación:


[math] \frac{ \partial h }{ \partial t } - D · \Delta h = 0, \quad \rho \gt \rho _{0} \quad θ\in (0,2\pi ) \quad t\gt0 \quad (1) [/math]

Donde [math] D= \frac{k}{s}[/math] se conoce como difusividad hidráulica.


Pozo

1 Obtención del Laplaciano y ecuación diferencial en polares

El Laplaciano es la divergencia del gradiente. Al aplicarle ambos operadores a nuestra función, nos da que:

[math]\Delta h(\rho,\theta)[/math] = [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]

y aplicándolo a la fórmula (1) nos quedaría:

[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], [math]\rho \gt[/math] [math]\rho _{0} [/math]

Finalmente, al decirnos que h debe depender solamente de [math]\rho[/math],es decir h=h([math]\rho[/math]), se deduce que [math]h _{θ}[/math]=0 y por lo tanto [math]h _{θθ}[/math]=0, por lo que nuestra ecuación final será:

[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 Sistema completo de ecuaciones

Las condiciones de frontera serán: 

[math]h(\rho _{0},t) = h _{\rho0} \quad \quad h(20,t) = h _{0}[/math]

Para que el problema tenga una única solución, nos falta una condición inicial del tipo [math]h(\rho,0) = h _{i}(\rho)[/math]. El problema quedará:

\[\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 _{i}(\rho) , & \end{matrix}\right.\]

3 Resolución del problema por diferencias finitas y método del trapecio

Para obtener resultados concretos al problema planteado en el apartado 2, le damos valores numéricos a los parámetros y resolvemos el problema por el método de diferencias finitas aplicando el método del trapecio. [math] \rho _{0} =1m., \quad h _{0}=40m., \quad h _{\rho}=35, \quad D=-10^{-2}, \quad h(\rho,0) = h _{0}[/math] 

%apartado 3 trabajo2 ecuaciones

r0=1;rN=20;D=10^(-2);

dt=0.1;dr=dt;

N=(rN-r0)/dr;T=2;

r=(r0+dr):dr:(rN-dr);

t0=0;tN=T;

t=t0:dt:T;

H=ones(1,length(r));H=40*H;H=H';

h0=35;hN=40;

sol(1,:)=[h0 H' hN];

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

rf=r0:dr:rN;

figure(1)

[RR,TT]=meshgrid(rf,t);

surf(RR,TT,sol)



%pintamos ahora el comportamiento de los puntos con r=2 a lo largo del

%tiempo

x=sol(:,2);

figure (2)
Grafica metodo separación de variables


Observamos en la gráfica como el nivel piezométrico va variando a lo largo del tiempo. La variación es mínima ya que t final es igual a 2. Podemos observar que en la situación inicial hay un salto de discontinuidad debido a que la condición inicial es h=40 y la condición de frontera correspondiente a ρ=0 es h=35, que es la altura que alcanza el agua al extraerla del pozo. Podemos observar que el "pico" que hay al principio se va curvando, lo que representa su variación a lo largo del tiempo como hemos dicho anteriormente.

4 Dibujo del comportamiento del nivel piezométrico

Aprovechamos el programa anterior para ver el comportamiento del nivel piezométrico en una gráfica 2D nivel piezométrico/tiempo en el punto ρ=2.

Grafica nivel piezomético/tiempo

Aquí podemos observar con mucha mas precisión la deformación mencionada en el apartado 3 ya que en un punto concreto observamos la variación, y no como antes que la variación era en cada punto del terreno

5 Método de Euler

Hacemos lo mismo que en el apartado 3 pero en vez de aplicar el método del trapecio, aplicamos los métodos de Euler implícito y explícito. 

%apartado 5 de trabajo euler implicito

clear all

r0=1;rN=20;D=10^(-2);

dt=0.1;dr=dt;

N=(rN-r0)/dr;T=2;

r=(r0+dr):dr:(rN-dr);

t0=0;tN=T;

t=t0:dt:T;

H=ones(1,length(r));H=40*H;H=H';

h0=35;hN=40;

sol(1,:)=[h0 H' hN];

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)));

    sol(k+1,:)=[h0 H' hN];

end

rf=r0:dr:rN;

[RR,TT]=meshgrid(rf,t);

surf(RR,TT,sol)


%apartado 5 trabajo euler explicito
clear all
r0=1;rN=20;D=10^(-2);
dr=0.1;dt=0.5*dr^2;
N=(rN-r0)/dr;T=2;
r=(r0+dr):dr:(rN-dr);
t0=0;tN=T;
t=t0:dt:T;
H=ones(1,length(r));H=40*H;H=H';
h0=35;hN=40;
sol(1,:)=[h0 H' hN];
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=(H)+dt.*(-(K-L)*H+F);
    sol(k+1,:)=[h0 H' hN];
end
rf=r0:dr:rN;
[RR,TT]=meshgrid(rf,t);
surf(RR,TT,sol)
Método de Euler implicito
Método de Euler explicito

6 Nivel piezométrico en intervalos de tiempo grandes

%apartado 6
clear all

r0=1;rN=20;D=10^(-2);
dt=1;dr=0.1;
N=(rN-r0)/dr;T=5000;
r=(r0+dr):dr:(rN-dr);
t0=0;tN=T;
t=t0:dt:T;
H=ones(1,length(r));H=40*H;H=H';
h0=35;hN=40;
sol(1,:)=[h0 H' hN];
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
rf=r0:dr:rN;
figure(1)
[RR,TT]=meshgrid(rf,t);
surf(RR,TT,sol)
Nivel piezométrico con soluciones cada 200 horas.

7 Estado estacionario

El nivel piezométrico varía desde que se empieza a extraer agua del pozo en t=0 hasta un determinado tiempo t a partir del cual el nivel piezométrico se mantiene constante a lo largo del tiempo, siempre que se mantengan las mismas condiciones de permeabilidad y se continúe con la extracción de agua. Por tanto, una vez que se ha llegado a dicha t donde el nivel no varía con el paso del tiempo, se puede despreciar la variación de h respecto del tiempo, [math]h _{t}[/math] ,en la ecuación hallada anteriormente.

Como consecuencia, se considerará que el nivel piezométrico solamente variará respecto de [math]\rho[/math]:

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


%apartado 7 trabajo

clear all

r0=1;rN=20;D=10^(-2);

dt=0.1;dr=dt;

N=(rN-r0)/dr;T=2;

r=(r0+dr):dr:(rN-dr);

t0=0;tN=T;

t=t0:dt:T;

H=ones(1,length(r));H=40*H;H=H';

h0=35;hN=40;

sol(1,:)=[h0 H' hN];

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=(((K-L1)))\(F);

    sol(k+1,:)=[h0 H' hN];

end

rf=r0:dr:rN;

figure(1)

[RR,TT]=meshgrid(rf,t);

surf(RR,TT,sol)
Nivel piezométrico despreciando [math]h _{t}[/math]

Al despreciar [math]h _{t}[/math] existe un pequeño error de distancia entre la solución que hemos adoptado como estacionaria y la solución al problema u(ρ,t).

Para diferentes valores de t, ésta distancia varía:

Para t=0, la distancia es de 4,8150 m.

Para t=100, la distancia es de 2,4080 m.

Para t=1000, la distancia es de o,96910 m.

Para t=5000, la distancia es de 1,0918*10^-4 m.

8 Capacidad de recuperación en acuífero

El acuífero recuperará el nivel piezométrico inicial de valor [math]h _{0}[/math] al impermeabilizarse el pozo, ya que el agua no se filtrará y el pozo ya no influirá en este nivel pi, por lo tanto con el paso del tiempo, evolucionará de nuestro estado estacionario del cuál partimos ahora, hasta estabilizarse en el nivel inicial de nuevo. 

%apartado8
r0=1;rN=20;r1=r0:dr:(rN-dr);
N=(rN-r0)/dr;
dt1=0.5;
T1=5000;
t01=0;tN1=T;
t1=t01:dt1:T1;
hN1=40;
P=sol(length(t),:);P(length(t))=[];P=P';H1=[P' hN1];
sol1(1,:)=H1';
K1=2.*diag(ones(1,N))-diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);K1(1,2)=-2;
K1=K1./(dr^2);K1=D*K1;
L2=diag(zeros(1,N))+diag(ones(1,N-1),1)-diag(ones(1,N-1),-1);L2(1,2)=0;
L2=L2./(2*dr);L2=D*L2;
L2=(diag(1./r1))*L2;
F1=zeros(length(r1),1);
F1(N)=D*((hN1/(dr^2))+((1/r1(N)))*(hN1/(2*dr)));
for k=1:length(t1)-1
    P=(eye(N)+((K1-L2).*(dt1/2)))\(P+((dt1/2).*(-(K1-L2)*P+F1+F1)));
    sol1(k+1,:)=[P' hN];
end
figure(2)
[RR1,TT1]=meshgrid(rf,t1);
surf(RR1,TT1,sol1)

for j=1:T1
    x=max(abs(sol1(j,:)-40))/40;
    if x<0.05
        break
    end
end


Nivel piezométrico impermeabilizando las paredes del pozo.

Una vez impermeabilizadas las paredes del pozo, el nivel piezométrico volverá a recuperar su nivel inicial y volverá a ser estacionario en 2712 horas aproximadamente (con un error de 5%).

9 Método de Fourier

En este caso despreciamos el radio del pozo y lo cerramos cuando ya hemos extraído todo el agua de él. De esta manera el acuífero queda de nuevo homogéneo, como antes de empezar el problema. Tomamos la región ρ<20 y usamos para este apartado el método de Fourier. Ahora tomamos como condición de frontera h(ρ,t)=40. Partimos del problema inicial:

Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Nivel_piezométrico_G5&oldid=13216»