Nivel piezométrico en un acuifero confinado. Grupo 10-B

De MateWiki
Saltar a: navegación, buscar
Trabajo realizado por estudiantes
Título Nivel piezométrico en un acuifero confinado. Grupo 10-B
Asignatura Ecuaciones Diferenciales
Curso Curso 2016-17
Autores Pablo Merayo Vidal(1909) María Guadalupe Aranda Sánchez (2049) Marcos Buitrago Peña (656) Dominique Manuela Cazar Espín (1850) Ignacio del Río Pérez (1693) Javier Martín Salgado(1899) Sara Fernández Gago (1663)
Este artículo ha sido escrito por estudiantes como parte de su evaluación en la asignatura


1 Introducción

Consideramos un acuífero confinado entre dos capas de terreno impermeables horizontales. Llamaremos h(x, y) al nivel piezométrico del acuífero, definido por la altura que alcanzaría el agua si hiciéramos un sondeo en el punto (x, y). Obviamente este valor depende de la presión a la que se encuentra el agua confinada en el acuífero. Supondremos que el acuífero es un medio poroso saturado de agua que ocupa una región infinita y está en equilibrio de manera que [math]h(x, y) = h_0[/math] constante para todo (x, y) ∈ [math]IR^{2}[/math].

Sobre el acuifero se construye un pozo de sección circular y radio[math] ρ_0[/math] para extraer agua. La presencia del pozo hace que el nivel piezométrico cambie. Dada la simetría 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 [math]ρ = \sqrt{x^{2} + y^{2}}[/math] . 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:

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

junto con la ley de Darcy

[math] q = −K∇h [/math]

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 este 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:

[math] \frac{∂h}{∂t} - D∆h = 0, ρ \gt ρ_0, θ ∈ (0, 2π), t \gt 0, (1) [/math]

donde D = K/S es la difusividad hidráulica que supondremos constante.

Acu (2).png

2 Deducción de la fórmula

Si tomamos coordenadas cilíndricas de manera que OZ coincide con el eje de simetría del pozo, [math]\ h = h_\rho [/math] donde [math] \rho [/math] [math] = \sqrt{x^2\; +\; y^2\;} [/math] . Trabajamos por tanto en coordenadas polares en el plano [math](\rho ,\theta) [/math].

La primera ecuación con la que podemos determinar [math]\ h = h_\rho [/math] es la de conservación de la masa:

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

siendo S el almacenamiento específico, que expresa la masa de agua extraída (o almacenada) por unidad de volumen de acuífero cuando desciende el nivel piezométrico [math] h [/math].

Para conocer completamente [math]\ h = h_\rho [/math] hace falta combinar con una segunda ecuación la de la conservación de la masa , que en este caso es la que explica la ley de Darcy:

:[math] q = - K ·\nabla h [/math]

siendo K la permeabilidad o conductividad hidraulica y deducida para cada material experimentalmente. El flujo de agua que provoca un cambio de presión será mayor, cuanto mayor sea K.


Combinando ambas ecuaciones tenemos como resultado:

[math] \frac{ \partial h }{ \partial t } - D · \Delta h = 0, \qquad \; \rho \gt \rho _{0}[/math] siendo [math] D=\frac{k}{s} [/math] la difusividad hidráulica que la suponemos constante.


En un caso general , la [math]h[/math] dependería de [math]\rho[/math] y [math]\theta[/math].En esta ecuación queda por desarrollar el Laplaciano de h. Partiremos de la función h(ρ,θ), de la cual calcularemos primero el gradiente\[ \nabla h = \frac{\partial h}{\partial \rho}\hat{\rho} + \frac{\partial h}{\partial \rho}\hat{\theta} = \frac{\partial h}{\partial \rho}\bar{\rho} + \frac{1}{\rho^{2}} \frac{\partial h}{\partial \rho}\bar{\theta} \] teniendo en cuenta que \(\bar{g^{i} } = \frac{\bar{g_i } }{ | \bar{g_i } | ^{2} } \) Ahora hacemos el divergente del gradiente anterior y obtenemos el Laplaciano en coordenadas polares\[ \Delta h = div(∇h) = {1 \over \rho} {\partial \over \partial \rho} \left( \rho {\partial h \over \partial \rho} \right) + {1 \over \rho^2} {\partial^2 h \over \partial \theta^2} = {1 \over \rho} {\partial h \over \partial \rho} + {\partial^2 h \over \partial \rho^2} + {1 \over \rho^2} {\partial^2 h \over \partial \theta^2} \]

Si consideramos que, debido a la simetría del problema, h sólo depende de su distancia al pozo, entonces obtenemos la función radial h(ρ), que nos permite simplificar el Laplaciano del modo siguiente\[ \Delta h = div(∇h) = {1 \over \rho} {\partial h \over \partial \rho} + {\partial^2 h \over \partial \rho^2} \]


Ciñéndonos a una región finita, [math]\rho\in (\rho_{0},20)[/math], y considerando que la altura del pozo se mantiene constante [math]h_{p}[/math], 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):

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

Debido al carácter transitorio del problema se necesita una condición inicial, que implica al tiempo dado en un instante inicial [math]t = 0[/math] y que haga referencia al estado del nivel piezométrico [math]h_{0}[/math], antes de ejecutar la excavación del pozo, es decir [math]h(\rho,0) = h _{0}[/math]. 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.\]

3 Resolución numérica por el método de diferencias finitas

Utilizamos el método de diferencias finitas con el objetivo de aproximar la ecuación en derivadas parciales a un problema discreto cuya resolución conduce a un sistema lineal de ecuaciones en cada paso de tiempo. Este método permite tratar problemas con fronteras irregulares sin tener que definir ecuaciones especiales para cada una de estas. Además, como la ubicación nodal no se ajustar necesariamente a una malla regular, es posible ubicar con precisión las extracciones y discontinuidades de entorno y del propio acuífero. Los resultados de estas simulaciones permitieron demostrar la utilidad de este método para tratar geometrías complejas y distribuciones irregulares de extracciones y definir el patrón natural de flujo y las alteraciones ocasionadas por las extracciones, así como su influencia en los niveles piezométricos.

Suponiendo \(ρ_0 =0.5m\), \(h_0 =45m\), \(hp =34\), D=10−2, \((ρ,0)=h_0\); y tomando ∆ρ = 0.1 y ∆t = ∆ρ realizaremos aproximaciones del problema expuesto: 

a=0.5;
b=20;
t0=0;
tn=2;
q=0.01;
g='0.*x+0.*t';
g=inline(g,'x','t');
ca='34+0.*t';
ca=inline(ca,'t');
cb='45+0.*t';
cb=inline(cb,'t');
u0='45+0.*x';
u0=inline(u0,'x');
h=0.1;
N=round((b-a)/h);
x=linspace(a,b,N+1);
x=x';
xx=x(2:end-1);
%K derivada parcial respecto ro de segundo orden
K=q/(h^2).*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
U0=u0(xx);
%B derivada parcial respecto ro de primer orden
a1='1./x';
a1=inline(a1,'x');
A1=a1(xx);
B=q/(2*h).*(diag(A1(1:end-1),1)-diag(A1(2:end),-1));
k=h;
M=round((tn-t0)/k);
t=linspace(t0,tn,M+1);
sol=zeros(N-1,M+1);
sol(:,1)=U0;
%heun
for i=1:M
    G1=g(xx,t(i));
    G1(1)=G1(1)+q.*ca(t(i))/(h^2);
    G1(end)=G1(end)+q.*cb(t(i))/(h^2);
    G2=g(xx,t(i)+k/2);
    G2(1)=G2(1)+q*ca(t(i)+k/2)/h^2;
    G2(end)=G2(end)+q*cb(t(i)+k/2)/h^2;
    K1=-(K+B)*sol(:,i)+G1;
    K2=-(K+B)*(sol(:,i)+1/2*K1*k)+G2;
    sol(:,i+1)=sol(:,i)+k/2*(K1+K2);
end
Ua=ca(t);
Ub=cb(t);
sol=[Ua;sol;Ub]';
[Mx,Mt]=meshgrid(x,t);
surf(Mx,Mt,sol)


Mediante este código obtenemos la siguiente gráfica: Figuraapartado3.pdf 4 Gráfica del comportamiento del nivel piezométrico 

clear
	closeall
        a=0.5;
	b=20;
	t0=0;
	tn=2;
	q=0.01;
	g='0.*x+0.*t';
	g=inline(g,'x','t');
	ca='34+0.*t';
	ca=inline(ca,'t');
	cb='45+0.*t';
	cb=inline(cb,'t');
	u0='45+0.*x';
	u0=inline(u0,'x');
	h=0.1;
	N=round((b-a)/h);
	x=linspace(a,b,N+1);
	x=x';
	xx=x(2:end-1);
	%K derivada parcial respecto ro de segundo orden
	K=q/(h^2).*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1));
	U0=u0(xx);
	%B derivada parcial respecto ro de primer orden
	a1='1./x';
	a1=inline(a1,'x');
	A1=a1(xx);
	B=q/(2*h).*(diag(A1(1:end-1),1)-diag(A1(2:end),-1));
	k=h;
	M=round((tn-t0)/k);
	t=linspace(t0,tn,M+1);
	sol=zeros(N-1,M+1);
	sol(:,1)=U0;
	%heun
	for i=1:M
	    G1=g(xx,t(i));
	    G1(1)=G1(1)+q.*ca(t(i))/(h^2);
	    G1(end)=G1(end)+q.*cb(t(i))/(h^2);
	G2=g(xx,t(i)+k/2);
	    G2(1)=G2(1)+q*ca(t(i)+k/2)/h^2;
	    G2(end)=G2(end)+q*cb(t(i)+k/2)/h^2;
	K1=-(K+B)*sol(:,i)+G1;
	    K2=-(K+B)*(sol(:,i)+1/2*K1*k)+G2;
	    sol(:,i+1)=sol(:,i)+k/2*(K1+K2);
	end
	Ua=ca(t);
	Ub=cb(t);
	sol=[Ua;sol;Ub]';
	[Mx,Mt]=meshgrid(x,t);
	surf(Mx,Mt,sol)

	x=sol(:,2);
        plot(x,t)
	xlabel('Nivel Piezometrico')
	ylabel('tiempo')

Gráfica Dibujamos el comportamiento del nivel piezométrico en los puntos para los cuales ρ = 2 en una gráfica 2D nivel piezométrico/tiempo

Aa4 (2).png

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

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;


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

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


Gráfica del apartado 5 Grafica5.png Gráfica del apartado 6 Archivo:PREGUN.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 ecuacion de modelizacion resultante es: 

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



% PETICIÓN DE DATOS

ca=input('introduce la condicion inicial: ');                   % condición dirichlet en a (= 34)
cb=input('introduce la condición final: ');                     % condición dirichtlet en b (=45)
p0=input('introduce el valor inicial de p: ');                  % valor incial de p (= 0,5)
pN=input('introduce el valor final de p: ');                    % valor final de p (=20)
D=input('introduce el valor de la difusividad hidráulica: ');   % valor de la difusividad (=0,01)
k=input('introduce el paso espacial: ');                        % valor del paso espacial (= 0,1)
h=input('introduce el valor del paso temporal: ');              % valor del paso temporal (=200)
t0=input('introduce el valor para t0: ');                       % valor inicial del tiempo (= 0)
tN=input('introduce el valor para tN: ');                       % valor final del tiempo (= 5000)
p=(p0+k):k:(pN-k);                                              % vector p
t=t0:h:tN;                                                      % vector de tiempos
N=round((pN-p0)/k);                                             % numero de elementos

% INTRODUCIMOS LAS FUNCIONES

h1=@(p)0*p+1;                                                   % coeficiente de la derivada de H respecto p
g=@(p)0*p;                                                      % término independiente


% CREAMOS EL VECTOR POSICIÓN

x=linspace(p0,pN,N+1);
x=x';
xx=x(2:end-1);

% CREAMOS NUETRO PROBLEMA DE CONTORNO

% CÁLCULO DE B

H1=h1(xx);
B=(diag(H1(1:end-1),1)-diag(H1(2:end),-1))/(2*k);
B=diag(1./p)*B;

% CÁLCULO DE K

K=(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1))/k^2;

% CÁLCULO DE G

G=g(xx);
G(1)=G(1)+ca/k^2+ca/(2*k*p(1));                                 % condicion dirichlet
G(end)=G(end)+cb/k^2-cb/(2*k*p(end));                           % condicion dirichlet

% CÁLCULO DE H

H=(K+B)\G;
H=[ca;H;cb];                                                    % introducimos la condición

% CREMOS Y DONDE VAMOS A SUSTITUIR LOS VALORES DEL BUCLE

y=zeros(length(t),length(H));
y(1,:)=H;

% HACEMOS QUE H SEA MATRIZ

 for i=1:length(t)-1
     H=(K-B)\G;
     H=[ca;H;cb];
     y(i+1,:)=H;
 end

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

% DIBUJAMOS

surf(Mx,Mt,Y)


Gráfico:Grafia7.png


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

% APARTADO 8 

% INTRODUCCIÓN DE DATOS
% vamos a cambiar unicamente los datos que cambien respecto al apartado
% anterior para que se ejecute todo sobre el mismo

h=input('introduce el valor del paso temporal: ');       % nuevo valor del paso temporal (= 0.1)
p=p0:k:(pN-k);                                           % no tiene las mismas dimensiones que antes
M=round((tN-t0)/h);                                      % numero de elementos del vector tiempo
t=linspace(t0,tN,M+1);                                   % vector de tiempos

% CÁLCULO DE LAS FUNCIONES

g=@(r,t)0*r;                                             % término independiente
hra=@(t)0*t;                                             % condición Neumann
cb=@(t)0*t+45;                                           % condición Dirichlet
H0=y(end,:);
H0=H0(1:N);                                              % condición incial

% CÁLCULO DEL VECTOR POSICIÓN

x=linspace(p0,pN,N+1);
x=x';
xx=x(1:end-1);
l=length(xx);

% CREAMOS NUESTRO PROBLEMA DE VALOR INICIAL

% CÁLCULO DE B

H1=h1(xx);
B=(diag(H1(1:end-1),1)-diag(H1(2:end),-1))/(2*k);
B=diag(1./p)*B;
B(1,2)=0;                                                % condición Neumann

% CÁLCULO DE K

K=(2*diag(ones(1,N))-diag(ones(1,N-1),-1)-diag(ones(1,N-1),1))/k^2;
K(1,2)=-2/k^2;

% CREAMOS UNA MATRIZ NULA PARA PODER IR SUSTITUYENDO DATOS

y=zeros(N,M+1);
y(:,1)=H0;

% MÉTODO DE HEUN

 for i=1:M       
     G=g(xx,t(i));  
     G(1)=G(1)+2*D*hra(t(i))/k;                         % condición Neumann
     G(end)=G(end)+D*cb(t(i))/k^2;                      % condición Dirichlet     
     Gi=g(xx,t(i)+k);       
     Gi(1)=Gi(1)+2*D*hra(t(i))/k;                       % condición Neumann
     Gi(end)=Gi(end)+D*cb(t(i))/k^2;                    % condición Dirichlet     
     k1=-D*(K+B)*y(:,i)+G;
     k2=-D*(K+B)*y(:,i)+Gi+k1*h;
     y(:,i+1)=y(:,i)+h*(k1+k2)/2;
 end
 
% DEFINICIÓN DE LA CONDICIÓN DIRICHLET PARA METERLA EN LA MATRIZ SOLUCIÓN

Hb=cb(t);

% CREAMOS LA MATRIZ SOLUCIÓN
% introducimos la condición Dirichlet

Y=[y;Hb];

% CREAMOS EL MALLADO

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

% DIBUJAMOS

surf(Mx,Mt,Y)
title('Nivel piezométrico en función del t y la X');
xlabel('Eje X');ylabel('Eje T');
zlabel('H');


Gráfica apartado 8 GRAFICA 8.png

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) 

%Introducimos los autovalores del sistema, que hemos calculado, tomándolos
  %de menor a mayor. 
  x1=;
  x2=;
  x3=;
  x4=;
  x5=;
  %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
Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Nivel_piezométrico_en_un_acuifero_confinado._Grupo_10-B&oldid=39596»