Source: https://mat.caminos.upm.es/wiki/Nivel_piezom%C3%A9trico_en_un_acuifero_confinado._Grupo_10-B
Timestamp: 2018-12-17 06:53:49+00:00

Document:
Nivel piezométrico en un acuifero confinado. Grupo 10-B - MateWiki
Nivel piezométrico en un acuifero confinado. Grupo 10-B
Título Nivel piezométrico en un acuifero confinado. Grupo 10-B
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)
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.
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);
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: 4 Gráfica del comportamiento del nivel piezométrico
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);
46 	Ua=ca(t);
47 	Ub=cb(t);
48 	sol=[Ua;sol;Ub]';
49 	[Mx,Mt]=meshgrid(x,t);
50 	surf(Mx,Mt,sol)
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
5 Resolución del problema y 6 Gráfica del nivel piezométrico para tiempos grandes
1 clear all; clc; clf;
3 % DATOS DEL PROBLEMA
5 % DATOS ESCALARES
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
13 % FUNCIONES
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');
28 %	DISCRETIZACIÓN DEL INTERVALO [a,b]
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
35 % CONDICIONES DE CONTORNO
37 % Dos dirichlet D-D --> sistema de N-1 ecuaciones
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);
46 % CONDICIONES INICIALES
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
52 % RESOLUCIÓN POR EL MÉTODO DEL TRAPECIO
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));
65 % Añadimos condiciones Dirichlet
66 Ua=ca(t);
67 Ub=cb(t);
68 solT=[Ua;solT;Ub];
70 % RESOLUCIÓN POR EL MÉTODO DE EULER IMPLICITO
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);
80 % Añadimos condiciones Dirichlet
81 Ua=ca(t);
82 Ub=cb(t);
83 solEI=[Ua;solEI;Ub];
85 % POSTPROCESO
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;
92 [Mt,Mx]=meshgrid(t,x);
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');
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;
103 [Mt,Mx]=meshgrid(t,x);
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 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
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
16 % INTRODUCIMOS LAS FUNCIONES
18 h1=@(p)0*p+1;                                                   % coeficiente de la derivada de H respecto p
19 g=@(p)0*p;                                                      % término independiente
22 % CREAMOS EL VECTOR POSICIÓN
24 x=linspace(p0,pN,N+1);
25 x=x';
26 xx=x(2:end-1);
28 % CREAMOS NUETRO PROBLEMA DE CONTORNO
30 % CÁLCULO DE B
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;
36 % CÁLCULO DE K
38 K=(2*diag(ones(1,N-1))-diag(ones(1,N-2),-1)-diag(ones(1,N-2),1))/k^2;
40 % CÁLCULO DE G
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
46 % CÁLCULO DE H
48 H=(K+B)\G;
49 H=[ca;H;cb];                                                    % introducimos la condición
51 % CREMOS Y DONDE VAMOS A SUSTITUIR LOS VALORES DEL BUCLE
53 y=zeros(length(t),length(H));
54 y(1,:)=H;
56 % HACEMOS QUE H SEA MATRIZ
58  for i=1:length(t)-1
59      H=(K-B)\G;
60      H=[ca;H;cb];
61      y(i+1,:)=H;
62  end
64  [Mt,Mx]=meshgrid(t,x);
66 % DIBUJAMOS
68 surf(Mx,Mt,Y)
8 Estimación de la capacidad de recuperación de acuífero
1 % APARTADO 8
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
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
12 % CÁLCULO DE LAS FUNCIONES
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
20 % CÁLCULO DEL VECTOR POSICIÓN
22 x=linspace(p0,pN,N+1);
23 x=x';
24 xx=x(1:end-1);
25 l=length(xx);
27 % CREAMOS NUESTRO PROBLEMA DE VALOR INICIAL
29 % CÁLCULO DE B
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
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;
41 % CREAMOS UNA MATRIZ NULA PARA PODER IR SUSTITUYENDO DATOS
43 y=zeros(N,M+1);
44 y(:,1)=H0;
46 % MÉTODO DE HEUN
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
60 % DEFINICIÓN DE LA CONDICIÓN DIRICHLET PARA METERLA EN LA MATRIZ SOLUCIÓN
62 Hb=cb(t);
64 % CREAMOS LA MATRIZ SOLUCIÓN
65 % introducimos la condición Dirichlet
67 Y=[y;Hb];
69 % CREAMOS EL MALLADO
71 [Mt,Mx]=meshgrid(t,x);
73 % DIBUJAMOS
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
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
Obtenido de «https://mat.caminos.upm.es/w/index.php?title=Nivel_piezométrico_en_un_acuifero_confinado._Grupo_10-B&oldid=39596»
Esta página fue modificada por última vez el 7 may 2017 a las 21:47.

References: Resolución 
 resolución 
 Resolución 
 RESOLUCIÓN 
 RESOLUCIÓN 
 resolución