Diferencia entre revisiones de «Nivel piezométrico G5»

Revisión del 13:15 18 may 2014

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

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

5 Método de Euler

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

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=(H)+dt.*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=200;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)



%apartado8

dt1=100;

T1=2000;

t01=0;tN1=T;

t1=t01:dt1:T1;

h01=40;hN1=40;

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

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

figure(2)

[RR1,TT1]=meshgrid(rf,t1);

surf(RR1,TT1,sol1)

Nivel piezométrico despreciando [math]h _{t}[/math]

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.

9 Método de Fourier