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:

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


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

 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 
52 	x=sol(:,2);
53         plot(x,t)
54 	xlabel('Nivel Piezometrico')
55 	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

  1 clear all; clc; clf;
  2 
  3 % DATOS DEL PROBLEMA
  4 
  5 % DATOS ESCALARES
  6 
  7 a=input('Extremo izquierdo a: ');       % Extremo izquierdo
  8 b=input('Extremo derecho b: ');         % Extremo derecho 
  9 t0=input('Instante inicial: ');         % Instante inicial
 10 tM=input('Instante final: ');           % Instante final
 11 D=input('Difusividad hidraúlica: ');    % Difusividad hidráulica
 12 
 13 % FUNCIONES
 14 
 15 % FUNCIÓN G: Función término independiente (x será vector y t escalar)
 16 g=input('Forma analitica función termino independiente: ','s');
 17 g=inline(g,'x','t');
 18 % FUNCIÓN CA: Función f(t) de la condición de contorno en x=a (t vector)
 19 ca=input('Función f(t) de la condicion de contorno en x=a: ','s');
 20 ca=inline(ca,'t'); 
 21 % FUNCIÓN CB: función f(t) de la condición de contorno en x=b (t vector)
 22 cb=input('Función f(t) de la condicion de contorno en x=b: ','s');
 23 cb=inline(cb,'t');
 24 % FUNCIÓN U0: Forma analítica f(x) función valor inicial (x vector)
 25 u0=input('Forma analitica u0(x) función valor inicial: ','s');
 26 u0=inline(u0,'x');
 27 
 28 %	DISCRETIZACIÓN DEL INTERVALO [a,b]
 29 
 30 h=input('Anchura de paso espacial: ');      % Paso espacial
 31 N=round((b-a)/h);                           % Subintervalos en x
 32 x=linspace(a,b,N+1);                        % Creación del vector x                       
 33 x=x';                                       % El vector x tiene N+1 elementos
 34 
 35 % CONDICIONES DE CONTORNO
 36 
 37 % Dos dirichlet D-D --> sistema de N-1 ecuaciones
 38 
 39 % Tomamos elementos del segundo al penúltimo (N-1 componentes)
 40 xx=x(2:end-1);
 41 % K tendrá (N-1 x N-1)
 42 K=D/h^2*(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1)); 
 43 % Valor inicial
 44 U0=u0(xx);
 45         
 46 % CONDICIONES INICIALES
 47 
 48 k=input('introduzca el paso: '); % paso en t
 49 M=round((tM-t0)/k); % numero de subintervalos en t
 50 t=linspace(t0,tM,M+1); % vector de los tiempos
 51 
 52 % RESOLUCIÓN POR EL MÉTODO DEL TRAPECIO
 53 
 54 solT=zeros(N-1,M+1);                % N-1 (D-D) // M+1 (debido a los tiempos)
 55 solT(:,1)=U0;                       % Primera columna matriz solucion= valor inicial
 56 for i=1:M                           % El bucle obtiene las diferentes columnas
 57 G1T=g(xx,t(i));
 58 G1T(1)=G1T(1)+D*ca(t(i))/h^2;  
 59 G1T(end)=G1T(end)+D*cb(t(i))/h^2;
 60 G2T=g(xx,t(i+1));
 61 G2T(1)=G2T(1)+D*ca(t(i+1))/h^2;  
 62 G2T(end)=G2T(end)+D*cb(t(i+1))/h^2;
 63 solT(:,i+1)=(eye(size(K))+k/2*K)\(solT(:,i)+k/2*(-K*solT(:,i)+G1T+G2T));
 64 end
 65 % Añadimos condiciones Dirichlet
 66 Ua=ca(t);
 67 Ub=cb(t);
 68 solT=[Ua;solT;Ub];
 69  
 70 % RESOLUCIÓN POR EL MÉTODO DE EULER IMPLICITO
 71 
 72 solEI=zeros(N-1,M+1);                % N-1 (D-D) // M+1 (debido a los tiempos)
 73 solEI(:,1)=U0;                       % Primera columna matriz solucion= valor inicial
 74 for i=1:M                            % El bucle obtiene las diferentes columnas
 75 G2EI=g(xx,t(i+1));
 76 G2EI(1)=G2EI(1)+D*ca(t(i+1))/h^2; 
 77 G2EI(end)=G2EI(end)+D*cb(t(i+1))/h^2;
 78 solEI(:,i+1)=(eye(size(K))+k*K)\(solEI(:,i)+k*G2EI);
 79 end
 80 % Añadimos condiciones Dirichlet
 81 Ua=ca(t);
 82 Ub=cb(t);
 83 solEI=[Ua;solEI;Ub];
 84 
 85 % POSTPROCESO
 86 
 87 % SOLUCIÓN APARTADO 5
 88 % introduciendo los datos: 
 89 % escalares: a=0.5; b=20; t0=0; tM=2; D=0.01; k=0.1; h=k;
 90 % funciones: g=0.*x + 0.*t; ca=34+0.*t; cb= 45+0.*t; u0=45+0.*x;
 91 
 92 [Mt,Mx]=meshgrid(t,x);
 93 
 94 subplot(2,2,1), surf(Mt,Mx,solT), title('MÉTODO TRAPECIO');
 95 subplot(2,2,2), surf(Mt,Mx,solEI), title('MÉTODO EULER IMPLICITO');
 96 
 97 % SOLUCIÓN APARTADO 6
 98 % introduciendo los datos: 
 99 % escalares: a=0.5; b=20; t0=0; tM=5000; D=0.01; k=200; h=0.1;
100 % funciones: g=0.*x + 0.*t; ca=34+0.*t; cb= 45+0.*t; u0=45+0.*x;
101 
102 
103 [Mt,Mx]=meshgrid(t,x);
104 
105 subplot(2,2,1), plot3(Mt,Mx,solT); title('MÉTODO TRAPECIO'); 
106 subplot(2,2,2), plot3(Mt,Mx,solEI); title('MÉTODO EULER IMPLICITO');
107 subplot(2,2,3), surf(Mt,Mx,solT); title('MÉTODO TRAPECIO'); 
108 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\]


 1 % PETICIÓN DE DATOS
 2 
 3 ca=input('introduce la condicion inicial: ');                   % condición dirichlet en a (= 34)
 4 cb=input('introduce la condición final: ');                     % condición dirichtlet en b (=45)
 5 p0=input('introduce el valor inicial de p: ');                  % valor incial de p (= 0,5)
 6 pN=input('introduce el valor final de p: ');                    % valor final de p (=20)
 7 D=input('introduce el valor de la difusividad hidráulica: ');   % valor de la difusividad (=0,01)
 8 k=input('introduce el paso espacial: ');                        % valor del paso espacial (= 0,1)
 9 h=input('introduce el valor del paso temporal: ');              % valor del paso temporal (=200)
10 t0=input('introduce el valor para t0: ');                       % valor inicial del tiempo (= 0)
11 tN=input('introduce el valor para tN: ');                       % valor final del tiempo (= 5000)
12 p=(p0+k):k:(pN-k);                                              % vector p
13 t=t0:h:tN;                                                      % vector de tiempos
14 N=round((pN-p0)/k);                                             % numero de elementos
15 
16 % INTRODUCIMOS LAS FUNCIONES
17 
18 h1=@(p)0*p+1;                                                   % coeficiente de la derivada de H respecto p
19 g=@(p)0*p;                                                      % término independiente
20 
21 
22 % CREAMOS EL VECTOR POSICIÓN
23 
24 x=linspace(p0,pN,N+1);
25 x=x';
26 xx=x(2:end-1);
27 
28 % CREAMOS NUETRO PROBLEMA DE CONTORNO
29 
30 % CÁLCULO DE B
31 
32 H1=h1(xx);
33 B=(diag(H1(1:end-1),1)-diag(H1(2:end),-1))/(2*k);
34 B=diag(1./p)*B;
35 
36 % CÁLCULO DE K
37 
38 K=(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1))/k^2;
39 
40 % CÁLCULO DE G
41 
42 G=g(xx);
43 G(1)=G(1)+ca/k^2+ca/(2*k*p(1));                                 % condicion dirichlet
44 G(end)=G(end)+cb/k^2-cb/(2*k*p(end));                           % condicion dirichlet
45 
46 % CÁLCULO DE H
47 
48 H=(K+B)\G;
49 H=[ca;H;cb];                                                    % introducimos la condición
50 
51 % CREMOS Y DONDE VAMOS A SUSTITUIR LOS VALORES DEL BUCLE
52 
53 y=zeros(length(t),length(H));
54 y(1,:)=H;
55 
56 % HACEMOS QUE H SEA MATRIZ
57 
58  for i=1:length(t)-1
59      H=(K-B)\G;
60      H=[ca;H;cb];
61      y(i+1,:)=H;
62  end
63 
64  [Mt,Mx]=meshgrid(t,x);
65 
66 % DIBUJAMOS
67 
68 surf(Mx,Mt,Y)


Gráfico:Grafia7.png


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

 1 % APARTADO 8 
 2 
 3 % INTRODUCCIÓN DE DATOS
 4 % vamos a cambiar unicamente los datos que cambien respecto al apartado
 5 % anterior para que se ejecute todo sobre el mismo
 6 
 7 h=input('introduce el valor del paso temporal: ');       % nuevo valor del paso temporal (= 0.1)
 8 p=p0:k:(pN-k);                                           % no tiene las mismas dimensiones que antes
 9 M=round((tN-t0)/h);                                      % numero de elementos del vector tiempo
10 t=linspace(t0,tN,M+1);                                   % vector de tiempos
11 
12 % CÁLCULO DE LAS FUNCIONES
13 
14 g=@(r,t)0*r;                                             % término independiente
15 hra=@(t)0*t;                                             % condición Neumann
16 cb=@(t)0*t+45;                                           % condición Dirichlet
17 H0=y(end,:);
18 H0=H0(1:N);                                              % condición incial
19 
20 % CÁLCULO DEL VECTOR POSICIÓN
21 
22 x=linspace(p0,pN,N+1);
23 x=x';
24 xx=x(1:end-1);
25 l=length(xx);
26 
27 % CREAMOS NUESTRO PROBLEMA DE VALOR INICIAL
28 
29 % CÁLCULO DE B
30 
31 H1=h1(xx);
32 B=(diag(H1(1:end-1),1)-diag(H1(2:end),-1))/(2*k);
33 B=diag(1./p)*B;
34 B(1,2)=0;                                                % condición Neumann
35 
36 % CÁLCULO DE K
37 
38 K=(2*diag(ones(1,N))-diag(ones(1,N-1),-1)-diag(ones(1,N-1),1))/k^2;
39 K(1,2)=-2/k^2;
40 
41 % CREAMOS UNA MATRIZ NULA PARA PODER IR SUSTITUYENDO DATOS
42 
43 y=zeros(N,M+1);
44 y(:,1)=H0;
45 
46 % MÉTODO DE HEUN
47 
48  for i=1:M       
49      G=g(xx,t(i));  
50      G(1)=G(1)+2*D*hra(t(i))/k;                         % condición Neumann
51      G(end)=G(end)+D*cb(t(i))/k^2;                      % condición Dirichlet     
52      Gi=g(xx,t(i)+k);       
53      Gi(1)=Gi(1)+2*D*hra(t(i))/k;                       % condición Neumann
54      Gi(end)=Gi(end)+D*cb(t(i))/k^2;                    % condición Dirichlet     
55      k1=-D*(K+B)*y(:,i)+G;
56      k2=-D*(K+B)*y(:,i)+Gi+k1*h;
57      y(:,i+1)=y(:,i)+h*(k1+k2)/2;
58  end
59  
60 % DEFINICIÓN DE LA CONDICIÓN DIRICHLET PARA METERLA EN LA MATRIZ SOLUCIÓN
61 
62 Hb=cb(t);
63 
64 % CREAMOS LA MATRIZ SOLUCIÓN
65 % introducimos la condición Dirichlet
66 
67 Y=[y;Hb];
68 
69 % CREAMOS EL MALLADO
70 
71 [Mt,Mx]=meshgrid(t,x);
72 
73 % DIBUJAMOS
74 
75 surf(Mx,Mt,Y)
76 title('Nivel piezométrico en función del t y la X');
77 xlabel('Eje X');ylabel('Eje T');
78 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)

 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