Diferencia entre revisiones de «Nivel piezométrico G5»

Revisión actual del 01:22 23 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

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 ρ=1 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

La diferencia entre los dos métodos de Euler ,implícito y explícito. La diferencia esta en la partición dt. En el primero, la partición dt es la misma que dr , mientras que en el segundo [math] dt \leq \frac{1}{2}dr^{2}. [/math] De esta manera, habrá muchos mas puntos en el método explícito y de ahí que toda la gráfica nos salga de color negro.

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.

En esta gráfica observamos que a diferencia de la primera gráfica, la curva que representa el cambio del nivel piezométrico es mucho más pronunciada. Esto es debido a que el intervalo de tiempo es muchísimo más grande ( T=5000) y que dt=1, con lo que el programa dibujará soluciones cada 1 intervalos de tiempo. También podemos observar que el pico debido a la discontinuidad inicial en la primera gráfica desaparece en esta. Esto se debe a que t=0 está muy muy próximo a t=5 o t= 20 con lo que el pico directamente pasa a ser una pequeña curva que con el paso del tiempo se pronuncia más.

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] el nivel piezométrico pasa del estado inicial al estacionario directamente. La distancia entre la ecuación original y el nivel estacionario para diferentes valores de t es:

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:

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

Estas condiciones de frontera no son homogéneas con lo que habrá que cambiar el sistema de manera que satisfaga unas condiciones de frontera homogéneas. Operando nos queda que el nuevo sistema será de la forma:

\[\left\{\begin{matrix}\ \frac{\partial \tilde{h}}{\partial t}- D·(\frac{\partial ^2 \tilde{h}}{\partial \rho^2}+ \frac{1}{\rho}·\frac{\partial \tilde{h}}{\partial \rho})=0 , & \\ \tilde{h}_{\rho}(0,t)=0 \quad \tilde{h}(20,t)=0, & \\ \tilde{h}(\rho ,0)=f(\rho)-40 & \end{matrix}\right.\]

Decimos que [math] \tilde{h}(ρ,t)=X(ρ).T(t) [/math] ya que la solución dependerá tanto de ρ como de t. Por variables separadas obtenemos que:

\[\left\{\begin{matrix}\ X′′(ρ)+ \frac{X′(ρ)}{ρ} + λ X(ρ)=0 & \\ X(20)=0 \quad \quad X′(0)=0 & \\ & \end{matrix}\right.\]

Este es el problema de autovalores asociado. Sabemos que la solución a la ecuación de la forma:

[math] r^{2} {\varphi}''(r) + r{\varphi}'(r))+ r^{2}\varphi(r)=0 \quad h _{0}=40m., \quad \varphi(0)=1 \quad {\varphi}'(0)=0 \quad [/math]

es conocida como función de Bessel de primera especie, [math] J_{0}(r) [/math]. Las autofunciones son de la forma [math] \varphi_{k}(ρ)=J_{0}(ρ\sqrt{λ})[/math] Sustituimos las autofunciones en la ecuación que pasa a ser [math] λ J_{0}''(ρ\sqrt{λ} + \sqrt{λ} \frac{J_{0}'(ρ\sqrt{λ})}{ρ} + λJ_{0}(ρ\sqrt{λ})=0 [/math]. Las condiciones de contorno pasan a ser [math] J_{0}(20\sqrt{λ})=0 [/math] y [math] J_{0}'(0)=0 [/math] Al comparar las ecucaciones sacamos que [math] r=ρ\sqrt{λ} [/math] .Las autofunciones satisfacen el problema de autovalores, satisface X'(0)=0 y X(20)=0. Imponemos la seguna condición y usamos los 5 primeros ceros de la función para calcular con la relación [math] λ= \frac{r^{2}}{\rho^{2}} [/math] los 5 primeros autovalores. Estos 5 primeros autovalores son:

\[\left\{\begin{matrix}\ λ_{1}=0,0145 & \\ λ_{2}=0,0762 & \\ λ_{3}=0,1872 & \\ λ_{4}=0,3476 & \\ λ_{5}=0,5573 & \end{matrix}\right.\]

%apartado 9
x1=0.0145;
x2=0.0762;
x3=0.1872;
x4=0.3476;
x5=0.5575;
r=0:0.01:20;
y1=besselj(0,r.*(sqrt(x1)));
y2=besselj(0,r.*(sqrt(x2)));
y3=besselj(0,r.*(sqrt(x3)));
y4=besselj(0,r.*(sqrt(x4)));
y5=besselj(0,r.*(sqrt(x5)));
hold on
plot(r,y1,'r')
plot(r,y2,'g')
plot(r,y3,'c')
plot(r,y4,'k')
plot(r,y5)
hold off

Gráfica con los 5 autovalores