Ver el código de «Nivel piezométrico (ECUS 7)»


{{ Beta }}

INTRODUCCIÓN
Sobre un acuífero confinado se construye un pozo de sección circular y radio ρ0 para extraer agua. La presencia de este pozo hace que el nivel piezómetro (altura que alcanzaría el agua si hiciéramos un sondeo y que depende de la presión a la que se encuentra el agua) cambie. 
Dada la simetría (cilíndrica) del problema podemos suponer que h(x, y) sólo va a depender de la distancia al pozo. Si tomamos coordenadas cilíndricas de forma que el eje OZ coincida con el eje de simetría del pozo, entonces h=h(ρ) donde ρ=√(x^2+y^2 ). Trabajaremos por tanto en coordenadas polares en el plano (ρ, θ). Las ecuaciones que permiten conocer h(ρ) son:
-	La ecuación de conservación de la masa     \[ S •\frac{ \partial h }{ \partial t } + div q = 0\]
-	La ley de Darcy \[ q = - K •\nabla h \]
Que es una ley experimental que modela el flujo de agua en un medio poroso. La ley de Darcy establece que el flujo de agua que 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. Cuanto mayor es la constante K mayor es el flujo de agua provocado por un cambio de presión.  
La ley de Darcy proporciona una buena aproximación del comportamiento del agua en un medio poroso siempre que éste sea homogéneo e isótropo. 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:

\[ \frac{ \partial h }{ \partial t } - D • \Delta h = 0, \qquad \; \rho > \rho _{0}\] siendo \( D=\frac{k}{s} \)
 Donde \( D=\frac{k}{s} \) es la difusividad hidráulica que la suponemos constante.

'''1 Cálculo del Laplaciano y ecuación diferencial en coordenadas polares'''
En un caso general , la \(h\) dependería de \(\rho\) y \(\theta\). Esta ecuación también se puede obtener si hallamos la divergencia del gradiente de h, es decir, el Laplaciano. Por la definición del Laplaciano, esto quedaría como\[\Delta h(\rho,\theta) = (\frac{\partial ^2 h}{\partial \rho^2}+ \frac{1}{\rho}•\frac{\partial h}{\partial \rho}+\frac{\partial^2 h}{\partial \theta^2}) \]

Lo primero para hallar q es el gradiente de h (formula gradiente), por tanto q es un campo vectorial. Si se quiere introducir en la ecuación de conservación de la masa, podríamos sustituir el valor de q obtenido dentro de la fórmula que nos obligaría a realizar la divergencia de dicho campo vectorial. Con esto obtendríamos la misma ecuación. Un vez hecho esto para \(\rho >\) \(\rho _{0}\) la ecuación diferencial sería:

\[ \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\]

Como queremos imponer que nuestra solución sea radial, obligamos a que \(h\) dependa únicamente de \(\rho\) , por tanto su derivada con respecto a \(\theta\) es igual a cero. La segunda por consiguiente también. Teniendo en cuenta todo esto, obtenemos la siguiente ecuación::

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


Ciñéndonos a una región finita, \(\rho\in (\rho_{0},20)\), y considerando que la altura del pozo se mantiene constante \(h_{p}\), la solución de nuestra ecuación (1) requiere dos condiciones de contorno que para este caso, se han propuesto de tipo Dirichlet (valor de la variable preescrito): \[h(\rho _{0},t) = h _{\rho} \quad \quad h(20,t) = h _{0}\]

Debido al carácter transitorio del problema se necesita una condición inicial, que implica al tiempo dado en un instante inicial \(t = 0\) y que haga referencia al estado del nivel piezométrico \(h_{0}\), antes de ejecutar la excavación del pozo, es decir \(h(\rho,0) = h _{0}\). Así para completar el sistema imponemos esa condición inicial para que nuestro problema tenga una única solución. Teniendo en cuenta todo lo anterior el sistema será el siguiente:

\[\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.\]

'''2 Condiciones de contorno y sistema de ecuaciones'''
Las condiciones de contorno impuestas (tipo Dirichlet), son una altura de pozo constante \(h_{p}\) y \(\rho\in (\rho_{0},20)\). Las condiciones serán: \[h(\rho _{0},t) = h _{\rho}\quad (1)\]\[ h(20,t) = h _{0}\quad (2)\]

Para que se tenga una solución única, tenemos que poner una última condición para \(t = 0\): \[h(\rho,0) = h _{0}\] 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'''

{{matlab|codigo=
1)	clear
2)	closeall
3)	a=0.5;
4)	b=20;
5)	t0=0;
6)	tn=2;
7)	q=0.01;
8)	g='0.*x+0.*t';
9)	g=inline(g,'x','t');
10)	ca='34+0.*t';
11)	ca=inline(ca,'t');
12)	cb='45+0.*t';
13)	cb=inline(cb,'t');
14)	u0='45+0.*x';
15)	u0=inline(u0,'x');
16)	h=0.1;
17)	N=round((b-a)/h);
18)	x=linspace(a,b,N+1);
19)	x=x';
20)	xx=x(2:end-1);
21)	%K derivada parcial respecto ro de segundo orden
22)	K=q/(h^2).*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
23)	U0=u0(xx);
24)	%B derivada parcial respecto ro de primer orden
25)	a1='1./x';
26)	a1=inline(a1,'x');
27)	A1=a1(xx);
28)	B=q/(2*h).*(diag(A1(1:end-1),1)-diag(A1(2:end),-1));
29)	k=h;
30)	M=round((tn-t0)/k);
31)	t=linspace(t0,tn,M+1);
32)	sol=zeros(N-1,M+1);
33)	sol(:,1)=U0;
34)	%heun
35)	for i=1:M
36)	    G1=g(xx,t(i));
37)	    G1(1)=G1(1)+q.*ca(t(i))/(h^2);
38)	    G1(end)=G1(end)+q.*cb(t(i))/(h^2);
39)	G2=g(xx,t(i)+k/2);
40)	    G2(1)=G2(1)+q*ca(t(i)+k/2)/h^2;
41)	    G2(end)=G2(end)+q*cb(t(i)+k/2)/h^2;
42)	K1=-(K+B)*sol(:,i)+G1;
43)	    K2=-(K+B)*(sol(:,i)+1/2*K1*k)+G2;
44)	    sol(:,i+1)=sol(:,i)+k/2*(K1+K2);
45)	end
46)	Ua=ca(t);
47)	Ub=cb(t);
48)	sol=[Ua;sol;Ub]';
49)	[Mx,Mt]=meshgrid(x,t);
50)	surf(Mx,Mt,sol)
}}

Gráfica
[[Archivo:Figuraapartado3.jpg]]
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'''

{{matlab|codigo=

1)	clear
2)	closeall
3)	a=0.5;
4)	b=20;
5)	t0=0;
6)	tn=2;
7)	q=0.01;
8)	g='0.*x+0.*t';
9)	g=inline(g,'x','t');
10)	ca='34+0.*t';
11)	ca=inline(ca,'t');
12)	cb='45+0.*t';
13)	cb=inline(cb,'t');
14)	u0='45+0.*x';
15)	u0=inline(u0,'x');
16)	h=0.1;
17)	N=round((b-a)/h);
18)	x=linspace(a,b,N+1);
19)	x=x';
20)	xx=x(2:end-1);
21)	%K derivada parcial respecto ro de segundo orden
22)	K=q/(h^2).*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
23)	U0=u0(xx);
24)	%B derivada parcial respecto ro de primer orden
25)	a1='1./x';
26)	a1=inline(a1,'x');
27)	A1=a1(xx);
28)	B=q/(2*h).*(diag(A1(1:end-1),1)-diag(A1(2:end),-1));
29)	k=h;
30)	M=round((tn-t0)/k);
31)	t=linspace(t0,tn,M+1);
32)	sol=zeros(N-1,M+1);
33)	sol(:,1)=U0;
34)	%heun
35)	for i=1:M
36)	    G1=g(xx,t(i));
37)	    G1(1)=G1(1)+q.*ca(t(i))/(h^2);
38)	    G1(end)=G1(end)+q.*cb(t(i))/(h^2);
39)	G2=g(xx,t(i)+k/2);
40)	    G2(1)=G2(1)+q*ca(t(i)+k/2)/h^2;
41)	    G2(end)=G2(end)+q*cb(t(i)+k/2)/h^2;
42)	K1=-(K+B)*sol(:,i)+G1;
43)	    K2=-(K+B)*(sol(:,i)+1/2*K1*k)+G2;
44)	    sol(:,i+1)=sol(:,i)+k/2*(K1+K2);
45)	end
46)	Ua=ca(t);
47)	Ub=cb(t);
48)	sol=[Ua;sol;Ub]';
49)	[Mx,Mt]=meshgrid(x,t);
50)	surf(Mx,Mt,sol)

51)	x=sol(:,2);
52)	plot(x,t)
53)	xlabel('Nivel Piezometrico')
54)	ylabel('tiempo')
}}
Gráfica
[[Archivo:figuraapartado4.jpg]]

'''5 Resolución del problema y 6 Gráfica del nivel piezométrico para tiempos grandes'''

{{matlab|codigo=
clear all; clc; clf;

% DATOS DEL PROBLEMA

% DATOS ESCALARES

a=input('Extremo izquierdo a: ');       % Extremo izquierdo
b=input('Extremo derecho b: ');         % Extremo derecho 
t0=input('Instante inicial: ');         % Instante inicial
tM=input('Instante final: ');           % Instante final
D=input('Difusividad hidraúlica: ');    % Difusividad hidráulica

% FUNCIONES

% FUNCIÓN G: Función término independiente (x será vector y t escalar)
g=input('Forma analitica función termino independiente: ','s');
g=inline(g,'x','t');
% FUNCIÓN CA: Función f(t) de la condición de contorno en x=a (t vector)
ca=input('Función f(t) de la condicion de contorno en x=a: ','s');
ca=inline(ca,'t'); 
% FUNCIÓN CB: función f(t) de la condición de contorno en x=b (t vector)
cb=input('Función f(t) de la condicion de contorno en x=b: ','s');
cb=inline(cb,'t');
% FUNCIÓN U0: Forma analítica f(x) función valor inicial (x vector)
u0=input('Forma analitica u0(x) función valor inicial: ','s');
u0=inline(u0,'x');

%	DISCRETIZACIÓN DEL INTERVALO [a,b]

h=input('Anchura de paso espacial: ');      % Paso espacial
N=round((b-a)/h);                           % Subintervalos en x
x=linspace(a,b,N+1);                        % Creación del vector x                       
x=x';                                       % El vector x tiene N+1 elementos

% CONDICIONES DE CONTORNO

% Dos dirichlet D-D --> sistema de N-1 ecuaciones

% Tomamos elementos del segundo al penúltimo (N-1 componentes)
xx=x(2:end-1);
% K tendrá (N-1 x N-1)
K=D/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1)); 
% Valor inicial
U0=u0(xx);
        
% CONDICIONES INICIALES

k=input('introduzca el paso: '); % paso en t
M=round((tM-t0)/k); % numero de subintervalos en t
t=linspace(t0,tM,M+1); % vector de los tiempos

% RESOLUCIÓN POR EL MÉTODO DEL TRAPECIO

solT=zeros(N-1,M+1);                % N-1 (D-D) // M+1 (debido a los tiempos)
solT(:,1)=U0;                       % Primera columna matriz solucion= valor inicial
for i=1:M                           % El bucle obtiene las diferentes columnas
G1T=g(xx,t(i));
G1T(1)=G1T(1)+D*ca(t(i))/h^2;  
G1T(end)=G1T(end)+D*cb(t(i))/h^2;
G2T=g(xx,t(i+1));
G2T(1)=G2T(1)+D*ca(t(i+1))/h^2;  
G2T(end)=G2T(end)+D*cb(t(i+1))/h^2;
solT(:,i+1)=(eye(size(K))+k/2*K)\(solT(:,i)+k/2*(-K*solT(:,i)+G1T+G2T));
end
% Añadimos condiciones Dirichlet
Ua=ca(t);
Ub=cb(t);
solT=[Ua;solT;Ub];
 
% RESOLUCIÓN POR EL MÉTODO DE EULER IMPLICITO

solEI=zeros(N-1,M+1);                % N-1 (D-D) // M+1 (debido a los tiempos)
solEI(:,1)=U0;                       % Primera columna matriz solucion= valor inicial
for i=1:M                            % El bucle obtiene las diferentes columnas
G2EI=g(xx,t(i+1));
G2EI(1)=G2EI(1)+D*ca(t(i+1))/h^2; 
G2EI(end)=G2EI(end)+D*cb(t(i+1))/h^2;
solEI(:,i+1)=(eye(size(K))+k*K)\(solEI(:,i)+k*G2EI);
end
% Añadimos condiciones Dirichlet
Ua=ca(t);
Ub=cb(t);
solEI=[Ua;solEI;Ub];

% POSTPROCESO

% SOLUCIÓN APARTADO 5
% introduciendo los datos: 
% escalares: a=0.5; b=20; t0=0; tM=2; D=0.01; k=0.1; h=k;
% funciones: g=0.*x + 0.*t; ca=34+0.*t; cb= 45+0.*t; u0=45+0.*x;

[Mt,Mx]=meshgrid(t,x);

subplot(2,2,1), surf(Mt,Mx,solT), title('MÉTODO TRAPECIO');
subplot(2,2,2), surf(Mt,Mx,solEI), title('MÉTODO EULER IMPLICITO');

% SOLUCIÓN APARTADO 6
% introduciendo los datos: 
% escalares: a=0.5; b=20; t0=0; tM=5000; D=0.01; k=200; h=0.1;
% funciones: g=0.*x + 0.*t; ca=34+0.*t; cb= 45+0.*t; u0=45+0.*x;

hold on
surf(Mt,Mx,SolT)
surf(Mt,Mx,SolEI)
hold off
}}

Gráfica
[[Archivo:Grafica5.png]]

'''7 Valor estacionario'''
Para tiempos grandes la solución apenas varía en tiempo y el término ht se puede despreciar. El nivel piezométrico toma un valor estacionario en cada punto. ¿Qué ecuaciones crees que debería verificar este estado estacionario? Resolver numéricamente ese estado estacionario usando el método de diferencias finitas. Calcular la distancia entre esa solución estacionaria y la solución u (ρ, t) para t = 0, 100, 1000, 5000

Una vez empezamos la extracción del agua del pozo, el nivel del agua va disminuyendo, provocando un cambio de presiones intersticiales en los alrededores que provocarán un flujo de agua de las partes más alejadas. 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 varía 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 ecuación de modelización resultante es:
 \[ D·(\frac{\partial ^2 h}{\partial \rho^2}+ \frac{1}{\rho}·\frac{\partial h}{\partial \rho})=0\]


{{matlab|codigo=
1)	%Datos 
2)p0=1;pN=20;D=10^(-2);
3)%Discretizacion
4) dt=0.1;dr=dt;
5)%Vectores tiempo y radio
6)N=(pN-p0)/dr;T=2;
7)r=(p0+dr):dr:(pN-dr);
8)t0=0;tN=T;
9)t=t0:dt:T;
10)%Inicializacion     
11)H=ones(1,length(r));H=40*H;H=H';
12)h0P=35;hN=40;
13)sol(1,:)=[h0P H' hN];
14)%Aplicacion
15)G=zeros(length(r),1);
16)G(1)=D*((h0P/(dr^2))-((1/r(1)))*(h0P/(2*dr)));
17)G(N-1)=D*((hN/(dr^2))+((1/r(N-1)))*(hN/(2*dr)));
18)%creo la matriz K para pasar el sistema a forma matricial y resolverlo
19)K=2.*diag(ones(1,N-1))-diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
20)K=K./(dr^2);K=D*K;
21)L=diag(zeros(1,N-1))+diag(ones(1,N-2),1)-diag(ones(1,N-2),-1);
22)L=L./(2*dr);L=D*L;
23)L1=(diag(1./r))*L;
24)%Tnteracciones
25)for k=1:length(t)-1
26)H=(((K-L1)))\(F);

27)sol(k+1,:)=[h0P H' hN];
28)end
29)%Dibujamos solucion
30)pf=p0:dr:pN;
31)[PP,TT]=meshgrid(pf,t);
32)surf(PP,TT,sol)
}}

Gráfico:
[[Archivo:grafica7.png]]


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

{{matlab|codigo=

1)%Datos
2)dt1=100;T1=2000;t01=0;tN1=T;h01=40;hN1=40;
3)%Vector tiempo
4)t1=t01:dt1:T1;
5)%Inicializacion
6)P=sol(21,:);
7)P(1)=[];
8)P(19)=[];
9)P=P';
10)H1=[h0P'hN];
11)sol1(1,:)=H1';
12)F1=zeros(length(r),1);F1(1)=D*((h01/(dr^2))-((1/r(1)))*(h01/(2*dr)));
13)F1(N-1)=D*((hN1/(dr^2))+((1/r(N-1)))*(hN1/(2*dr)));
14)%Interacciones
15)fork=1:length(t1)-1
16)P=(eye(N-1)+((K-L1).*(dt1/2)))\(P+((dt1/2).*(-(K-L1)*P+F1+F1)));
17)sol1(k+1,:)=[h0P'hN];
18)end
19)[RR1,TT1]=meshgrid(rf,t1);
20)surf(RR1,TT1,sol1)
}}

'''9 Sistema de Fourier'''

Vamos ahora a suponer que el radio del pozo es despreciable y que una vez extraído el agua el pozo se cierra quedando el acuífero de nuevo homogéneo. En este caso tomaremos la región ρ < 20. De nuevo supondremos que la solución sólo depende de ρ y t, y tomaremos ahora como condición frontera h (20, t) = 40. Aproximar la solución de la ecuación cuando el dato inicial es h˜(ρ, 0) = log |ρ+0.01| 20.01 .
(El método de resolución está pendiente de estudio)

{{matlab|codigo=

1 %Introducimos los autovalores del sistema, que hemos calculado, tomándolos
 2 %de menor a mayor. 
 3 x1=;
 4 x2=;
 5 x3=;
 6 x4=;
 7 x5=;
 8 %Vector con r
 9 r=0:0.01:20;
10 %Aplicamos la funcion besselj
11 z1=besselj(0,r.*(sqrt(x1)));
12 z2=besselj(0,r.*(sqrt(x2)));
13 z3=besselj(0,r.*(sqrt(x3)));
14 z4=besselj(0,r.*(sqrt(x4)));
15 z5=besselj(0,r.*(sqrt(x5)));
16 %Dibujamos las rectas obtenidas en la gráfica
17 hold on
18 plot(r,z1)
19 plot(r,z2)
20 plot(r,z3)
21 plot(r,z4)
22 plot(r,z5)
23 hold off

}}