This is an old revision of the document!
Water - air seepage - thermal coupled – vapour diffusion 2D/3D constitutive law for solid elements
This law is only used for water seepage - air seepage-thermal coupled and vapour diffusion for non linear analysis in 2D/3D/2nd gradient porous media.
\[ \underbrace{\frac{\partial}{\partial t} (\rho_s . n . S_{r,w}) + div(\rho_w \vec{q_l})}_{\text{Liquide}} + \underbrace{\frac{\partial}{\partial t} (\rho_v . n . S_{r,g}) + div(\rho_v \vec{q_g})}_{Vapeur} = 0 \]
En partant de l’équation de Darcy, la vitesse du liquide (Volume de fluide par unité de surface de sol) est donnée par : \[ \vec{q_l} = - \frac{k_w}{\mu_w}\left[ \vec{grad}(p_w) + g \rho_w \vec{grad}(y)\right]\ \text{où}\ k_w = K_w \frac{\mu_w}{\rho_w g}\left[ m^2\right] \]
La vapeur se déplace uniquement dans les pores occupés par la phase gazeuse et le trajet qu'elle parcourt dépend de la tortuosité : \[ \vec{i}_v = - n S_{r,g} \tau D \rho_s \vec{grad} \omega_v \] Où $\omega_v = \rho_v/\rho_g$ est la teneur massique d’air sec dans le mélange gazeux.
\[\frac{\partial}{\partial t} (\rho_a . n . S_{r,g}) + div(\rho_a \vec{q_g}) + div(\vec{i_a}) = 0\]
En partant de l’équation de Darcy , la vitesse de la phase gazeuse est donnée par : \[ \vec{q_g} = - \frac{k_g}{\mu_g}\left[ \vec{grad}(p_g) + g \rho_g \vec{grad}(y)\right]\ \text{où}\ k_g = K_g \frac{\mu_g}{\rho_g g}\left[ m^2\right] \] La vitesse de diffusion de l'air sec est liée à un gradient de la masse volumique de l'air. En utilisant la théorie de la diffusion adaptée au cas des milieux poreux, nous pouvons écrire : \[ \vec{i}_a = - n S_{r,g} \tau D \rho_g \vec{grad} \omega_a = -\vec{I}_v \] Où $\omega_a = \rho_a/\rho_g$ est la teneur massique d’air sec dans le mélange gazeux.
\[\dot{S_T} + \dot{E}_{H_2O}^{w\rightarrow v}L + div(\vec{V}_T) - Q_T = 0\]
Où $S_T$ représente la quantité de chaleur emmagasinée, $L$ la chaleur latente de vaporisation, $\vec{V}_T$ le flux de chaleur et $Q_T$ une source de chaleur en volume.
Cette dernière équation peut être transformée en utilisant l'équation de bilan de la vapeur d'eau: \[ \dot{S_T} + \dot{S}_vL + div(\vec{V}_T) + div(\vec{V}_v) L - Q_T = 0 \]
La quantité d'enthalpie du système est définie comme la somme des contributions de chaque espèce du système: \[ S_T = H_m = \sum_i H_i = \sum_i \theta_i \rho_i c_{p,i} (T-T_0) \] Les contributions de chaque composante à l’enthalpie du système s'exprime selon : \[ \begin{array}{l} H_w = n.S_{r,w} \rho_w c_{p,w} (T-T_0) \\ H_a = n(1-S_{r,g})\rho_ac_{p,a} (T-T_0)\\ H_s = (1-n)\rho_s c_{p,s} (T-T_0)\\ H_v = n(1-S_{r,g})\rho_v c_{p,v} (T-T_0) \end{array} \] Un dernier terme, lié à la vaporisation de l’eau, contribue également à l'emmagasinement de chaleur et dépend de la quantité de vapeur et la chaleur latente de vaporisation. \[ H_{Lat} = nS_{r,g} \rho_{v} L \]
\[ \vec{V_T} + \vec{V_v}L = \underbrace{- \Gamma_m \vec{\nabla}T}_{conduction} + \underbrace{\left[c_{p,w}\rho_w \vec{q}_l + c_{p,a}(\vec{i}_a + \rho_a \vec{q}_g) + c_{pv}(\vec{i}_v+\rho_v\vec{q}_g)\right](T-T_0)}_{convection} + \underbrace{\left[\vec{i}_v + \rho_v \vec{q}_g\right] L}_{Latente} \]
Prepro: LWAVAT.F
| Plane stress state | NO |
| Plane strain state | YES |
| Axisymmetric state | YES |
| 3D state | YES |
| Generalized plane state | NO |
| Line 1 (2I5, 60A1) | |
|---|---|
| IL | Law number |
| ITYPE | 171 (= 174 in LOI2 for 3D state) |
| COMMENT | Any comment (up to 60 characters) that will be reproduced on the output listing |
| Line 1 (18I5) | |
|---|---|
| IANI | = 0, isotropic permeability case |
| ≠ 0, anisotropic permeability case | |
| IKW | formulation index for $k_w$ |
| IKA | formulation index for $k_a$ |
| ISRW | formulation index for $S_W$ |
| ITHERM | formulation index for $\Gamma_T$ |
| IVAP | = 0 if no vapour diffusion in the problem |
| = 1 else | |
| IFORM | = 1 tangent formulation \[ \left\{\begin{array}{l}f_{we} = \dot{M}_w = (\dot{\varepsilon}_v S_w + n S_w\frac{\dot{\rho}_w}{\chi_w} + n \dot{S}_w) \rho_w \\ f_{ae} = \dot{M}_a = (\dot{\varepsilon}_V S_a + n S_a\frac{\dot{p}_a}{p_a} + n \dot{S}_a) \rho_a \end{array}\right. \] |
| = 0 secant formulation \[ \left\{\begin{array}{l} f_{we} = \dot{M}_w = (n^BS_w^B\rho_w^B - n^AS_w^A\rho_w^A)/\Delta t \\ f_{ae} = \dot{M}_a = (n^BS_a^B\rho_a^B - n^AS_a^A\rho_a^A)/\Delta t \end{array}\right. \] | |
| ICONV | = 0 if no convectif term in the heat transport problem |
| = 1 else | |
| ITEMOIN | = 0 if analytic matrix (can be used only if IVAP = 0 if no vapour diffusion in the problem) |
| = 1 if semi-analytic matrix (can be used in all the problems) | |
| IKRN | = 1 if Kozeni-Karmann formulation |
| = 2, GDR Momas relation $K=f(n)$ | |
| = 3 Coupling permeability-deformation $K=f(\varepsilon_n)$ (only in 2D) | |
| IGAS | = 0 if gas is air |
| = 1 if gas is $H_2$ | |
| = 2 if gas is $N_2$ | |
| = 3 if gas is Ar | |
| = 4 if gas is He | |
| = 5 if gas is $CO_2$ | |
| IENTH | = 0 if we define $\rho$ and $C_p$ for each constituent \[\left\{\begin{array}{l}H_w = N.S_{r,w} \rho_w c_{p,w} (T-T_0) \\ H_v = n(1-S_{r,w})\rho_v c_{p,v} (T-T_0)\\ H_a = n(1-S_{r,w})\rho_ac_{p,a} (T-T_0)\\ H_{a-d} = n S_{r,w} H\rho_ac_{p,a} (T-T_0)\\ H_s = (1-n)\rho_s c_{p,s} (T-T_0) \end{array}\right.\] |
| = 1 if we define $\rho C_p$ equivalent for the medium and constant \[H_m = \rho C_p (T-T_0)\] | |
| = 2 if we define $\rho C_p$ equivalent for the medium depending on temperature: \[\left\{\begin{array}{l} \rho C_p = RHOC1\ \text{if }\ T>T_u \\ \rho C_p = RHOC2\ \text{if }\ T<T_f \\ \rho C_p = \frac{RHOC1-RHOC2}{T_u-T_f}(T-T_u) + RHOC1\ \text{if }\ T_f\leq T \leq T_u\\ \end{array} \right.\] \[T_u = CLT3\] \[T_f = CLT4\] | |
| IANITH | = 0, isotropic conductivity case |
| = 1, anisotropic conductivity case | |
| IVISCW | = 0, $\mu_w = \mu_{w,0}\left( 1-ALPHW0(T-T_0)) \right)$ |
| = 1, $\mu_w = 0.6612(T-229)^{-1.562}$ | |
| IXHIW | = 0, constant water compressibility |
| = 1, $\chi_w = \chi_w + \frac{H}{p_a}$ ( $p_a$ is partial pressure of air and H is Henry coefficient) | |
| IDIFF | = 0, with diffusion of dissolved air |
| ≠ 0, divisor (integer becomes real) of diffusion coefficient of dissolved air | |
| ISTRUCT | = 0, constant permeability |
| = 1, permeability depends on microstructure evolution | |
| ICOAL | = 0, solid conductivity: $\Gamma_s = \Gamma_{s,0}(1+GAMS0(T-T_0))$ |
| = 1, solid conductivity: $\Gamma_s = \Gamma_{s,0}GAMS0(T-T_0)^3$ | |
The permeability $k_f$ is an \underline{intrinsic} permeability ($\left[L^2\right]$) $\boxed{ \begin{array}{l} k_{f, intrinsic} = K_f \frac{\mu_f}{\rho_f g}\\ \left[ L^2 \right] = \left[ LT^{-1} \right] \frac{ \left[ ML^{-1}T^{-1}\right]}{\left[ML^{-3}\right]\left[LT^{-2}\right]} \end{array}}$
| If IANI ≠ 0 (4G10.0) - Repeated IANI times | |
|---|---|
| PERME(I) | soil anisotropic int. permeability ($k_f$) in the direction I |
| COSX(I) | director cosinus of the direction I (in 3d state) |
| COSY(I) | |
| COSZ(I) | |
Permeabilities in different directions can be input ( $I_{max}= 10$ ).The effect of these permeabilities will be summed.
| If IANI = 0 (1G10.0) | |
|---|---|
| PERME | soil isotropic intrinsic permeability ($k_f$) |
| Line 1 (5G10.0) | |
|---|---|
| POROS | soil porosity $(= n)$ |
| TORTU | soil tortuosity $(=\tau )$ |
| T0 | definition temperature $(=T_0)\ \left[K\right]$ |
| PW0 | definition liquid pression $(=p_{w,0})\ \left[Pa\right]$ |
| PA0 | definition gaz pression $(=p_{a,0})\ \left[Pa\right]$ |
| Line 2 (7G10.0) | |
| VISCW0 | liquid dynamic viscosity $(=\mu_{w,0}\ \left[ Pa.s \right]$ |
| ALPHW0 | liquid dynamic viscosity thermal coefficient $(=\alpha_{w}^T)\ \left[K^{-1}\right]$ |
| RHOW0 | |
| UXHIW0 | liquid compressibility coefficient $(=1/ \chi_{w})\ \left[ Pa^{-1}\right]$ |
| BETAW0 | liquid thermal expansion coefficient $(=\beta{w}^T\ \left[K^{-1}\right]$ |
| CONW0 | liquid thermal conductivity $(=\Gamma_{w,0})\ \left[W.m^{-1}.K^{-1}\right]$ |
| GAMW0 | liquid thermal conductivity coefficient $(=\gamma_{w}^T\ \left[K^{-1}\right]$ |
| Line 3 (3G10.0) | |
| CPW0 | liquid specific heat $(=c_{p,w0})\ \left[J.kg^{-1}.K^{-1}\right]$ |
| HEATW0 | liquid specific heat coefficient $(=H_{w}^T)\ \left[K^{-1}\right]$ |
| EMMAG | storage coefficient $(=E_s)\ \left[ P_a^{-1}\right]$ |
| Line 4 (7G10.0) | |
| VISCA0 | |
| ALPHW0 | gaz dynamic viscosity thermal coefficient$(=\alpha_{a}^T)\ \left[K^{-1}\right]$ |
| RHOA0 | gaz density$(=\rho_{a,0})\ \left[kg.m^{-3}\right]$ |
| CONA0 | gaz thermal conductivity $(=\Gamma_{a,0})\ \left[W.m^{-1}.K^{-1}\right]$ |
| GAMA0 | gaz thermal conductivity coefficient $(=\gamma_{a}^T)\ \left[K^{-1}\right]$ |
| CPA0 | gaz specific heat $(=c_{p,a0})\ \left[J.kg^{-1}.K^{-1}\right]$ |
| HEATA0 | gaz specific heat coefficient $(=H_{a}^T)\ \left[K^{-1}\right]$ |
| Line 5 (5G10.0) | |
| BETAS0 | solid thermal expansion coefficient $(=\beta_{s}^T)\ \left[K^{-1}\right]$ |
| CONS0 | solid thermal conduction$(=\Gamma_{s,0})\ \left[W.m^{-1}.K^{-1}\right]$ |
| GAMS0 | solid conduction coefficient $(=\gamma_{s}^T)\ \left[K^{-1}\right]$ |
| CPS0 | solid specific heat $(=c_{p,s0})\ \left[J.kg^{-1}. K^{-1}\right]$ |
| HEATS0 | solid specific heat coefficient $(=H_{s}^T)\ \left[ K^{-1}\right]$ |
| Line 6 (3G10.0) | |
| CKW1 | 1st coefficient of the function $k_{rw}$ |
| CKW2 | 2nd coefficient of the function $k_{rw}$ |
| CKW3 | 3rd coefficient of the function $k_{rw}$ |
| Line 7 (3G10.0) | |
| CKA1 | 1st coefficient of the function $k_{ra}$ |
| CKA2 | 2nd coefficient of the function $k_{ra}$ |
| CSR5 | 5th coefficient of the function $S_w$ |
| Line 8 (7G10.0/) | |
| CSR1 | 1st coefficient of the function $S_w$ |
| CSR2 | 2nd coefficient of the function $S_w$ |
| CSR3 | 3rd coefficient of the function $S_w$ |
| CSR4 | 4th coefficient of the function $S_w$ |
| SRES | residual saturation degree $(=S_{res})$ |
| SRFIELD | field saturation degree $(=S_{r, field})$ |
| AIREV | air entry value $\left[Pa\right]$ |
| Line 9 (7G10.0) | |
| CLT1 | 1st coefficient of the function $\Gamma_T$ |
| CLT2 | 2nd coefficient of the function $\Gamma_T$ |
| CLT3 | 3rd coefficient of the function $\Gamma_T$ |
| CLT4 | 4th coefficient of the function $\Gamma_T$ |
| RHOC | coefficient for enthalpie $\rho C_p$ (if ienth = 1) |
| RHOC1 | 1st coefficient for enthalpie $\rho C_p$ (if ienth = 2) |
| RHOC2 | 2nd coefficient for enthalpie $\rho C_p$ (if ienth = 2) |
| Line 10 (4G10.0) | |
| KRMIN | Minimum value of $kr$ |
| HENRY | Henry coefficient |
| EXPM | for IKRN=1: $m$ Exponent of Kozeni-Karmann formulation |
| for IKRN=3 : A parameter | |
| EXPN | n Exponent of Kozeni-Karmann formulation |
| Line 11 - Only if IANITH ≠ 0 (4G10.0) - Repeat IANITH times | |
| CONDUC(I) | soil anisotropic conductivity in the direction I |
| COSX(I) | director cosinus of the direction I (in 3d state) |
| COSY(I) | |
| COSZ(I) | |
| Thermal conductivities in different directions can be input (Imax = 10). The effect of these conductivities will be summed. In that case of anisotropic conductivity, conductivities remain constants for each direction: coefficients ITHERM, CLT1, CLT2, CLT3, CT4, CONS0, CONW0, and CONA0 are so not used in that case | |
| If IANITH = 0: nothing | |
| Isotropic thermal conductivity is already defined by preceding coefficients ITHERM, CLT1, CLT2, CLT3, CT4, CONS0, CONW0, CONA0 … | |
Following empirical formulations for describing the evolution of the relative permeability, the thermal conductivity and saturation with the suction are possible: see Appendix 8.
For any suction lower than air entry value (AIREV), the saturation is equal to SRFIELD value.
Kozeny Karman formulation:
\[K = C_0 \frac{n^{EXPN}}{(1-n)^{EXPM}}\]
$C_0$ is computed automatically from $C_0 = K_0 \frac{(1-n_0)^{EXPM}}{(n_0)^{EXPn}}$
GDR Momas formulation:
\[
\frac{k}{k_0} = 1+EXPM\left[ \phi - \phi_0\right]^{EXPN}\ \text{où}\ EXPM = 2.10^{12}\ \text{et}\ EXPN = 3
\]
Coupling permeability-deformation formulation: (only in 2D)
\[
K_{ij} = \sum_n K_n^0 (1+A\varepsilon_n^T)^3\beta_{ij} (\alpha_n)
\]
$\varepsilon_n^T$ tensile deformation
$\alpha_n$ crack orientation with horizontal
$A$ parameter of the crack
24
In 2D state :
| SIG(1) | liquid velocity in the X direction $(=f_{wx})$ | |
| SIG(2) | liquid velocity in the Y direction $(=f_{wy})$ | |
| SIG(3) | liquid velocity stored $(=f_{we})$ | |
| SIG(4) | none | |
| SIG(5) | gas total velocity in the X direction $(=f_{ax})$ | gas advection + gas diffusion + dissolved gas advection + dissolved gas diffusion |
| SIG(6) | gas total velocity in the Y direction $(=f_{ay})$ | |
| SIG(7) | gas total velocity stored $(=f_{ae})$ | |
| SIG(8) | none | |
| SIG(9) | conductive heat flow in the X direction $(=f_{tx})$ | |
| SIG(10) | conductive heat flow in the Y direction $(=f_{ty})$ | |
| SIG(11) | energy accumulated by heat capacity $(=f_{te})$ | |
| SIG(12) | none | |
| SIG(13) | Water vapour velocity in the X direction $(=f_{vx})$ | |
| SIG(14) | Water vapour velocity in the Y direction $(=f_{vy})$ | |
| SIG(15) | Water vapour stored $(=f_{ve})$ | |
| SIG(16) | none | |
| SIG(17) | dissolved gas advection and diffusion velocity in the X direction | |
| SIG(18) | dissolved gas advection and diffusion velocity in the Y direction | |
| SIG(19) | dissolved gas advection and diffusion velocity stored | |
| SIG(20) | none | |
| SIG(21) | dissolved gas diffusion velocity in the X direction | |
| SIG(22) | dissolved gas diffusion velocity in the Y direction | |
| SIG(23) | dissolved gas and diffusion velocity stored | |
| SIG(24) | none | |
In 3D state :
| SIG(1) | liquid velocity in the X direction $(=f_{wx})$ | |
| SIG(2) | liquid velocity in the Y direction $(=f_{wy})$ | |
| SIG(3) | liquid velocity in the Z direction $(=f_{wz})$ | |
| SIG(4) | liquid velocity stored $(=f_{we})$ | |
| SIG(5) | gas total velocity in the X direction $(=f_{ax})$ | gas advection + gas diffusion + dissolved gas advection + dissolved gas diffusion |
| SIG(6) | gas total velocity in the Y direction $(=f_{ay})$ | |
| SIG(7) | gas total velocity in the Z direction $(=f_{az})$ | |
| SIG(8) | gas total velocity stored $(=f_{az})$ | |
| SIG(9) | conductive heat flow in the X direction $(=f_{tx})$ | |
| SIG(10) | conductive heat flow in the Y direction $(=f_{ty})$ | |
| SIG(11) | conductive heat flow in the Z direction $(=f_{tz})$ | |
| SIG(12) | energy accumulated by heat capacity $(=f_{te})$ | |
| SIG(13) | Water vapour velocity in the X direction $(=f_{yx})$ | |
| SIG(14) | Water vapour velocity in the Y direction $(=f_{yy})$ | |
| SIG(15) | Water vapour velocity in the Z direction $(=f_{yz})$ | |
| SIG(16) | Water vapour stored $(=f_{ye})$ | |
| SIG(17) | dissolved gas advection and diffusion velocity in the X direction | |
| SIG(18) | dissolved gas advection and diffusion velocity in the Y direction | |
| SIG(19) | dissolved gas advection and diffusion velocity in the Z direction | |
| SIG(20) | dissolved gas advection and diffusion velocity stored | |
| SIG(21) | dissolved gas diffusion velocity in the X direction | |
| SIG(22) | dissolved gas diffusion velocity in the Y direction | |
| SIG(23) | dissolved gas diffusion velocity in the Z direction | |
| SIG(24) | dissolved gas and diffusion velocity stored | |
= 26 in 2D cases
= 16 in 3D cases
| Q(1) | water relative permeability $(=k_{rw})$ |
| Q(2) | air relative permeability $(=k_{ra})$ |
| Q(3) | Soil porosity (= n) |
| Q(4) | Soil saturation degree $(=S_w)$ |
| Q(5) | Suction $(=p_c = p_a-p_w)$ |
| Q(6) | water specific mass $(=\rho_w)$ |
| Q(7) | air specific mass $(=\rho_a)$ |
| Q(8) | “Pe number” = convective effect / conductive effect \[= \frac{\rho_f . c_f . T . \vec{q}}{\Gamma_{av} . \vec{grad} (T)}\] |
| Q(9) | Water content (=w) |
| Q(10) | Vapour specific mass $(=\rho_v)$ |
| Q(11) | Vapour pressure $(=p_v)$ |
| Q(12) | Relative humidity $(=H_r)$ |
| Q(13) | Liquid water mass per unit soil volume |
| Q(14) | Dry air mass per unit soil volume |
| Q(15) | Vapour mass per unit soil volume |
| Q(16) | Intrinsic permeability |
| Q(17) | Gas soil saturation degree $(=S_g)$ |
| Q(18) | $\alpha (H_2, N_2, …)$ partial pressure $(=p_a^g = p^g - p_{H_2O}^g = \text{gas pressure-vapour pressure})$ |
| Q(19) | Area associated to one integration point |
| Q(20) | Dissolved air concentration $=\frac{\rho_{a-d}}{\rho_w + \rho_{a-d}} = \frac{H_a \rho_a}{\rho_w + H_a \rho_a}$ |
| Q(21) | $K_{xx}$ (or zero if IANI = 0) |
| Q(22) | $K_{yy}$ (or zero if IANI = 0) |
| Q(23) | $K_{xy}$ (or zero if IANI = 0) |
| Q(24) | $\varepsilon_1$ |
| Q(25) | $\varepsilon_2$ |
| Q(26) | $\alpha$ (= angle between principal stress and horizontal) |