Ecuación general de flujo no saturado
El flujo en zonas saturadas indica que cuando se alcanza una cierta zona(llamada zona saturada), la humedad es constante, así no habría nada que modelar en ese caso.
El flujo en zonas saturadas indica que cuando se alcanza una cierta zona(llamada zona saturada), la humedad es constante, así no habría nada que modelar en ese caso.
De acuerdo con la ley de Darcy, para una dimensión de flujo vertical, el volumen de flujo q(cm^3/cm^3/h) puede ser escrito como
$$q = -K \frac{\partial }{\partial z}(h-z) =-K (\frac{\partial h}{\partial z}-1) $$
Donde K es la conductividad hiráulica(cm/h), "h" es la presión de agua en el suelo(relativo a la atmósfera) expresado en cm de agua y "z" es la componente vertical considerado positivo cuando va hacia abajo.
Del paper obtenemos las ecuaciones a usar:
$$ C = \frac{\mathrm{d} \theta}{\mathrm{d} h} $$
$$C(h) \frac{\partial h}{\partial t} = \frac{\partial }{\partial z} [K(h) (\frac{\partial h}{\partial z}-1)]$$
La funcion theta indica la humedad para una profundidad y tiempo dado.
Para el estudio, las condiciones iniciales son definidas como
$$\theta (z, t=0) = 0.1$$
$$\theta (z=0, t) = 0.267$$
La primera condición indica que el suelo va estar recibiendo una fuente de fluido constantemente, siempre en z=0(la superficie) se tendra una humedad de 0.267.
La segunda condición indica que inicialmente el suelo estaba casi seco en la zona no saturada.
Simulación
Luego del tiempo 0 empieza la simulación y del paper de Celia(referenciadas al principio) tenemos las siguientes relaciones:
$$K = K_s \frac{A}{A+ \left | h^{\beta_1} \right |}$$
Donde Ks=34 cm/h, A=1175000, b1 =4.74 y
$$ \theta = \frac{\alpha (\theta_s \theta_r)}{\alpha+ \left | h^{\beta_2} \right |}+\theta_r$$
Donde theta_s =0.287, theta_r = 0.075, alfa = 1611000, b2 = 3.96.
Usando diferencias finitas
$$ C_i^{j+a} \frac{h_i^{j+1} -h_i^j}{\Delta t} = \frac{1}{\Delta z}\left [
K_{i+1/2}^{j+a}\left (\frac{h_{i+1}^{j+a} -h_i^{j+a}}{\Delta z} -1\right )-
K_{i-1/2}^{j+a}\left ( \frac{h_{i}^{j+a} -h_{i-1}^{j+a}}{\Delta z} -1\right ) \right ] $$
Donde "a" es un factor de peso entre 0 y 1, siendo a = 0 el esquema explicito, a = 0.5 el esquema de Crank - Nicolson y a = 1 el esquema implicito.
Ejemplo con 41 nodos de profundidad (0-40)
Si se desea el esquema explicito (a = 0), se puede emplear simplemente el siguiente código:
Los valores de c y de k se aproximan mediante:
$$ C_i^{j+a} = F_1 = (1-a)C_i^j + a C_i^{j+1}$$
$$ K_{i+1/2}^{j+a} = F_2 = (1-a)K_{i+1/2}^j + a K_{i+1/2}^{j+1}= (1-a)\sqrt{K_i^j K_{i+1}^j} +a\sqrt{K_i^{j+1} K_{i+1}^{j+1}}$$
$$ K_{i-1/2}^{j+a} = F_3 = (1-a)K_{i-1/2}^j + a K_{i-1/2}^{j+1}= (1-a)\sqrt{K_{i-1}^j K_{i}^j} +a\sqrt{K_{i-1}^{j+1} K_{i}^{j+1}}$$
Acomodando la ecuación para poder formar un sistema lineal
$$\left [ -aF_3 \frac{\Delta t}{(\Delta z)^2} \right ] h_{i-1}^{j+1}+
\left [ F_1+aF_2 \frac{\Delta t}{(\Delta z)^2}+aF_3 \frac{\Delta t}{(\Delta z)^2} \right ] h_{i}^{j+1}+
\left [ -aF_2 \frac{\Delta t}{(\Delta z)^2} \right ] h_{i+1}^{j+1}$$
$$ =
\left [ (1-a)F_3 \frac{\Delta t}{(\Delta z)^2} \right ] h_{i-1}^{j}+
\left [ F_1-(1-a)F_2 \frac{\Delta t}{(\Delta z)^2}-(1-a)F_3 \frac{\Delta t}{(\Delta z)^2} \right ] h_{i}^{j}+
\left [ (1-a)F_2 \frac{\Delta t}{(\Delta z)^2} \right ] h_{i+1}^{j}+
(F_3 - F_2)\frac{\Delta t}{\Delta z} $$
Tenemos que resolver el siguiente sistema lineal
$$\begin{bmatrix}
F_1^1+aF_2^1 \frac{\Delta t}{(\Delta z)^2}+aF_3^1 \frac{\Delta t}{(\Delta z)^2} & -aF_2^1 \frac{\Delta t}{(\Delta z)^2} & 0 & \cdots & 0\\
-aF_3^2 \frac{\Delta t}{(\Delta z)^2}& F_1^2+aF_2^2 \frac{\Delta t}{(\Delta z)^2}+aF_3^2 \frac{\Delta t}{(\Delta z)^2} & -aF_2^2 \frac{\Delta t}{(\Delta z)^2} & \cdots & 0\\
& & & & \\
\vdots & \ddots & \ddots & \ddots & \vdots \\
0& \cdots & 0 & -aF_3^{39} \frac{\Delta t}{(\Delta z)^2} & F_1^{39}+aF_2^{39} \frac{\Delta t}{(\Delta z)^2}+aF_3^{39} \frac{\Delta t}{(\Delta z)^2}
\end{bmatrix}
\begin{bmatrix}
h_{1}^{j+1}\\
h_{2}^{j+1}\\
h_{3}^{j+1}\\
\vdots \\
h_{39}^{j+1}
\end{bmatrix}$$
$$ = $$
$$\begin{bmatrix}
F_1^1-(1-a)F_2^1 \frac{\Delta t}{(\Delta z)^2}-(1-a)F_3^1 \frac{\Delta t}{(\Delta z)^2} &(1-a)F_2^1 \frac{\Delta t}{(\Delta z)^2} & 0 & \cdots & 0\\
(1-a)F_3^2 \frac{\Delta t}{(\Delta z)^2}& F_1^2-(1-a)F_2^2 \frac{\Delta t}{(\Delta z)^2}-(1-a)F_3^2 \frac{\Delta t}{(\Delta z)^2} & (1-a)F_2^2 \frac{\Delta t}{(\Delta z)^2} & \cdots & 0\\
& & & & \\
\vdots & \ddots & \ddots & \ddots & \vdots \\
0& \cdots & 0 & (1-a)F_3^{39} \frac{\Delta t}{(\Delta z)^2} & F_1^{39}-(1-a)F_2^{39} \frac{\Delta t}{(\Delta z)^2}-(1-a)F_3^{39} \frac{\Delta t}{(\Delta z)^2}
\end{bmatrix}
\begin{bmatrix}
h_{1}^{j}+(F_3^1 - F_2^1)\frac{\Delta t}{\Delta z}+aF_3^1 \frac{\Delta t}{(\Delta z)^2} h_0^{j+1} + (1-a)F_3^1 \frac{\Delta t}{(\Delta z)^2} h_0^{j} \\
h_{2}^{j}\\
h_{3}^{j}\\
\vdots \\
h_{39}^{j}+(F_3^{39} - F_2^{39})\frac{\Delta t}{\Delta z}+aF_3^{39} \frac{\Delta t}{(\Delta z)^2} h_{40}^{j+1} + (1-a)F_3^{39} \frac{\Delta t}{(\Delta z)^2} h_{40}^{j}
\end{bmatrix}$$
Valido para j = 0,1, ... , M , M es un número natural tal que dt = T/M(puede ser igual o redondearse al mínimo entero).
El código general es el siguiente, agrupe los términos por comodidad en t1,t2,t3 y ti(término independiente).
Este ejemplo esta con el esquema de Crank Nickolson (a = 0.5), puedes modificar haciendo varia el "a", "dt", "T", "n", "dz", consultar las referencias para ver los valores de convergencia.
Muestra para T=360, T=720 y T= 0.8 horas = 8 x 360 s
Muestra de los esquemas para dt = 0.4, T=360, n=40 y valores de a entre 0 y 1.
Para dt = 3 , se ve cuando a< 0.5 diverge y cuando es mayor converge.
Continuaré en otro post, espero haberles ayudado ;)
Donde K es la conductividad hiráulica(cm/h), "h" es la presión de agua en el suelo(relativo a la atmósfera) expresado en cm de agua y "z" es la componente vertical considerado positivo cuando va hacia abajo.
Del paper obtenemos las ecuaciones a usar:
$$ C = \frac{\mathrm{d} \theta}{\mathrm{d} h} $$
$$C(h) \frac{\partial h}{\partial t} = \frac{\partial }{\partial z} [K(h) (\frac{\partial h}{\partial z}-1)]$$
La funcion theta indica la humedad para una profundidad y tiempo dado.
Para el estudio, las condiciones iniciales son definidas como
$$\theta (z, t=0) = 0.1$$
$$\theta (z=0, t) = 0.267$$
La primera condición indica que el suelo va estar recibiendo una fuente de fluido constantemente, siempre en z=0(la superficie) se tendra una humedad de 0.267.
La segunda condición indica que inicialmente el suelo estaba casi seco en la zona no saturada.
Simulación
Luego del tiempo 0 empieza la simulación y del paper de Celia(referenciadas al principio) tenemos las siguientes relaciones:
$$K = K_s \frac{A}{A+ \left | h^{\beta_1} \right |}$$
Donde Ks=34 cm/h, A=1175000, b1 =4.74 y
$$ \theta = \frac{\alpha (\theta_s \theta_r)}{\alpha+ \left | h^{\beta_2} \right |}+\theta_r$$
Donde theta_s =0.287, theta_r = 0.075, alfa = 1611000, b2 = 3.96.
Usando diferencias finitas
$$ C_i^{j+a} \frac{h_i^{j+1} -h_i^j}{\Delta t} = \frac{1}{\Delta z}\left [
K_{i+1/2}^{j+a}\left (\frac{h_{i+1}^{j+a} -h_i^{j+a}}{\Delta z} -1\right )-
K_{i-1/2}^{j+a}\left ( \frac{h_{i}^{j+a} -h_{i-1}^{j+a}}{\Delta z} -1\right ) \right ] $$
Donde "a" es un factor de peso entre 0 y 1, siendo a = 0 el esquema explicito, a = 0.5 el esquema de Crank - Nicolson y a = 1 el esquema implicito.
Ejemplo con 41 nodos de profundidad (0-40)
Si se desea el esquema explicito (a = 0), se puede emplear simplemente el siguiente código:
thetaS = 0.287; thetaR = 0.075; alfa = 1.611e6; beta2 = 3.96; tf = @(h) (alfa.*(thetaS-thetaR))./(alfa+(abs(h)).^beta2) + thetaR; Ks = 0.00944; Ar = 1.175e6; beta1 = 4.74; %Ks debe estar en cm/s kf = @(h) Ks.*Ar./(Ar+ abs(h).^beta1); cf =@(h) -alfa*beta2*(thetaS-thetaR)*sign(h).*abs(h).^(beta2-1)./(alfa+abs(h).^beta2).^2; hf =@(t) abs(alfa.*(thetaS-t)./(t-thetaR)).^(1/beta2); z_max = 40; n=40; T=360; dz=z_max/n; dt=0.4; t=0.1*ones(1,n); h = -hf(t); h(1) = -hf(0.267); h(n) = -hf(0.1); k=kf(h); c =cf(h); temp = h; for j=1:dt:T for i=2:n-1 temp(i) = h(i)+dt/dz/c(i)*(sqrt(k(i)*k(i+1))*((h(i+1)-h(i))/dz -1)... - sqrt(k(i-1)*k(i))*((h(i)-h(i-1))/dz -1)); end h = temp; k = kf(h); c = cf(h); end z=0:n-1; plot(z,h,'k');
Los valores de c y de k se aproximan mediante:
$$ C_i^{j+a} = F_1 = (1-a)C_i^j + a C_i^{j+1}$$
$$ K_{i+1/2}^{j+a} = F_2 = (1-a)K_{i+1/2}^j + a K_{i+1/2}^{j+1}= (1-a)\sqrt{K_i^j K_{i+1}^j} +a\sqrt{K_i^{j+1} K_{i+1}^{j+1}}$$
$$ K_{i-1/2}^{j+a} = F_3 = (1-a)K_{i-1/2}^j + a K_{i-1/2}^{j+1}= (1-a)\sqrt{K_{i-1}^j K_{i}^j} +a\sqrt{K_{i-1}^{j+1} K_{i}^{j+1}}$$
Acomodando la ecuación para poder formar un sistema lineal
$$\left [ -aF_3 \frac{\Delta t}{(\Delta z)^2} \right ] h_{i-1}^{j+1}+
\left [ F_1+aF_2 \frac{\Delta t}{(\Delta z)^2}+aF_3 \frac{\Delta t}{(\Delta z)^2} \right ] h_{i}^{j+1}+
\left [ -aF_2 \frac{\Delta t}{(\Delta z)^2} \right ] h_{i+1}^{j+1}$$
$$ =
\left [ (1-a)F_3 \frac{\Delta t}{(\Delta z)^2} \right ] h_{i-1}^{j}+
\left [ F_1-(1-a)F_2 \frac{\Delta t}{(\Delta z)^2}-(1-a)F_3 \frac{\Delta t}{(\Delta z)^2} \right ] h_{i}^{j}+
\left [ (1-a)F_2 \frac{\Delta t}{(\Delta z)^2} \right ] h_{i+1}^{j}+
(F_3 - F_2)\frac{\Delta t}{\Delta z} $$
Tenemos que resolver el siguiente sistema lineal
$$\begin{bmatrix}
F_1^1+aF_2^1 \frac{\Delta t}{(\Delta z)^2}+aF_3^1 \frac{\Delta t}{(\Delta z)^2} & -aF_2^1 \frac{\Delta t}{(\Delta z)^2} & 0 & \cdots & 0\\
-aF_3^2 \frac{\Delta t}{(\Delta z)^2}& F_1^2+aF_2^2 \frac{\Delta t}{(\Delta z)^2}+aF_3^2 \frac{\Delta t}{(\Delta z)^2} & -aF_2^2 \frac{\Delta t}{(\Delta z)^2} & \cdots & 0\\
& & & & \\
\vdots & \ddots & \ddots & \ddots & \vdots \\
0& \cdots & 0 & -aF_3^{39} \frac{\Delta t}{(\Delta z)^2} & F_1^{39}+aF_2^{39} \frac{\Delta t}{(\Delta z)^2}+aF_3^{39} \frac{\Delta t}{(\Delta z)^2}
\end{bmatrix}
\begin{bmatrix}
h_{1}^{j+1}\\
h_{2}^{j+1}\\
h_{3}^{j+1}\\
\vdots \\
h_{39}^{j+1}
\end{bmatrix}$$
$$ = $$
$$\begin{bmatrix}
F_1^1-(1-a)F_2^1 \frac{\Delta t}{(\Delta z)^2}-(1-a)F_3^1 \frac{\Delta t}{(\Delta z)^2} &(1-a)F_2^1 \frac{\Delta t}{(\Delta z)^2} & 0 & \cdots & 0\\
(1-a)F_3^2 \frac{\Delta t}{(\Delta z)^2}& F_1^2-(1-a)F_2^2 \frac{\Delta t}{(\Delta z)^2}-(1-a)F_3^2 \frac{\Delta t}{(\Delta z)^2} & (1-a)F_2^2 \frac{\Delta t}{(\Delta z)^2} & \cdots & 0\\
& & & & \\
\vdots & \ddots & \ddots & \ddots & \vdots \\
0& \cdots & 0 & (1-a)F_3^{39} \frac{\Delta t}{(\Delta z)^2} & F_1^{39}-(1-a)F_2^{39} \frac{\Delta t}{(\Delta z)^2}-(1-a)F_3^{39} \frac{\Delta t}{(\Delta z)^2}
\end{bmatrix}
\begin{bmatrix}
h_{1}^{j}+(F_3^1 - F_2^1)\frac{\Delta t}{\Delta z}+aF_3^1 \frac{\Delta t}{(\Delta z)^2} h_0^{j+1} + (1-a)F_3^1 \frac{\Delta t}{(\Delta z)^2} h_0^{j} \\
h_{2}^{j}\\
h_{3}^{j}\\
\vdots \\
h_{39}^{j}+(F_3^{39} - F_2^{39})\frac{\Delta t}{\Delta z}+aF_3^{39} \frac{\Delta t}{(\Delta z)^2} h_{40}^{j+1} + (1-a)F_3^{39} \frac{\Delta t}{(\Delta z)^2} h_{40}^{j}
\end{bmatrix}$$
Valido para j = 0,1, ... , M , M es un número natural tal que dt = T/M(puede ser igual o redondearse al mínimo entero).
El código general es el siguiente, agrupe los términos por comodidad en t1,t2,t3 y ti(término independiente).
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 | % Solucón Numérica para la ecuacuón de Richards % Infiltración de líquido en medio poroso % % % % % Forma general tic (); thetaS = 0.287; thetaR = 0.075; alfa = 1.611e6; beta2 = 3.96; tf = @(h) (alfa.*(thetaS-thetaR))./(alfa+(abs(h)).^beta2) + thetaR; Ks = 0.00944; Ar = 1.175e6; beta1 = 4.74; %Ks debe estar en cm/s kf = @(h) Ks.*Ar./(Ar+ abs(h).^beta1); cf =@(h) -alfa*beta2*(thetaS-thetaR)*sign(h).*abs(h).^(beta2-1)./(alfa+abs(h).^beta2).^2; hf =@(t) abs(alfa.*(thetaS-t)./(t-thetaR)).^(1/beta2); z_max = 40; n=z_max; T=300; dz=z_max/n; dt=0.4; a = 0.5; u=dt/dz^2; t=0.1*ones(1,n); h = -hf(t); h(1) = -hf(0.267); h(n) = -hf(0.1); k=kf(h); c =cf(h); hj1 = h; cj1=c; kj1 = k; temp = h; for j=1:dt:360 %PREDICTOR for i=2:n-1 temp(i) = h(i)+dt/dz/c(i)*(sqrt(k(i)*k(i+1))*((h(i+1)-h(i))/dz -1)... - sqrt(k(i-1)*k(i))*((h(i)-h(i-1))/dz -1)); end hj1 = temp; kj1 = kf(hj1); cj1 = cf(hj1); %CORRECTOR for i=2:n-1 f1(i-1) = (1-a)*c(i)+a*cj1(i); f2(i-1) = (1-a)*sqrt(k(i)*k(i+1))+a*sqrt(kj1(i)*kj1(i+1)); f3(i-1) = (1-a)*sqrt(k(i)*k(i-1))+a*sqrt(kj1(i)*kj1(i-1)); t1(i-1) = -a*u*f3(i-1); t2(i-1) = f1(i-1)+a*u*f2(i-1)+a*u*f3(i-1); t3(i-1) = -a*u*f2(i-1); t1b(i-1) = (1-a)*u*f3(i-1); t2b(i-1) = f1(i-1)-(1-a)*u*f2(i-1)-(1-a)*u*f3(i-1); t3b(i-1) = (1-a)*u*f2(i-1); ti(i-1) = dt/dz*(f3(i-1)-f2(i-1)); end ti(1) = ti(1) - t1(1)*h(1)+t1b(1)*h(1); ti(n-2) = ti(n-2) - t3(n-2)*h(n)+t3b(n-2)*h(n); B=diag(t1b(2:end),-1)+diag(t2b)+... diag(t3b(1:n-3),+1); A=diag(t1(2:end),-1)+diag(t2)+... diag(t3(1:n-3),+1); b = ti; x0 = h(2:n-1); h(2:n-1) = A\(B*x0'+b'); k = kf(h); c = cf(h); temp = h; end z=0:n-1; elapsed_time = toc (); fprintf('se demoro %4.4f segundos \n',elapsed_time); plot(z,h,'r'); |
Este ejemplo esta con el esquema de Crank Nickolson (a = 0.5), puedes modificar haciendo varia el "a", "dt", "T", "n", "dz", consultar las referencias para ver los valores de convergencia.
Muestra para T=360, T=720 y T= 0.8 horas = 8 x 360 s
Muestra de los esquemas para dt = 0.4, T=360, n=40 y valores de a entre 0 y 1.
Casi no se nota la diferencia.
Para dt = 3 , se ve cuando a< 0.5 diverge y cuando es mayor converge.
Continuaré en otro post, espero haberles ayudado ;)
