Table of Contents

WAVAT2

Description

Water - air seepage - thermal coupled – vapour diffusion 2nd gradient constitutive law for solid elements\\

The law definition and typical values of parameters for clay materials can be found in Corman (2024)1).

The model

This law is only used for water seepage - air seepage-thermal coupled and vapour diffusion for non linear analysis using the 2nd gradient porous media.

Conservation de la masse d’eau (Liquide + vapeur)

\[ \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 \]

Ecoulement du liquide et de la vapeur

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.

Equations d’état du liquide

  1. Viscosité dynamique $\mu_w$: \[\mu_w (T) = \mu_{w,0} - \alpha_w^T \mu_{w0}(T-T_0)\]
  2. Masse volumique $\rho_w$: \[\rho_w (T, p_w) = \rho_{wo}\left[ 1+\frac{p_w-p_{w0}}{\chi_w} - \beta_w^T (T-T_0) \right]\]
  3. Perméabilité intrinsèque $k_w$:
    Dépendance avec le degré de saturation $S_w$ : $k_{r,w} = f(S_w)$ avec $k_{w,eff} = k_f k_{r,w}$
  4. Degré de saturation $S_w$:
    Dépendance avec la succion $s = p_a - p_w : S_w = f(s)$

Conservation de la masse de l’air sec

\[\frac{\partial}{\partial t} (\rho_a . n . S_{r,g}) + div(\rho_a \vec{q_g}) + div(\vec{i_a}) = 0\]

Ecoulement de la phase gazeuse et l’air sec

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.

Equation d’état de l’air sec

  1. Viscosité dynamique $\mu_a$ : dépendance avec la température \[\mu_a (T) = \mu_{a,0} - \alpha_a^T \mu_{a0}(T-T_0)\]
  2. Masse volumique $\rho_a$ :
    Hypothèse : L’air est considéré comme un gaz parfait. \[\rho_a (T, p_a) = \rho_{a,0}\frac{p_a}{p_{a,0}}\frac{T_0}{T} \]
  3. Perméabilité intrinsèque $k_g$:
    Dépendance avec le degré de saturation $S_g$ : $k_{r,g} = f(S_g)$ avec $k_{g,effectif} = k_{g, intrinsic}k_{a,w}$
  4. Degré de saturation en liquide $S_g$:
    Dépendance avec la succion $s = p_g - p_w$
    $S_g = 1-S_w$ d’où voir partie liquide

Conservation en volume de la chaleur

\[\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 \]

Quantité de chaleur emmagasinée par unité de volume

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 \]

Transfert de la chaleur par unité de volume

\[ \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} \]

Equations d’état

  1. $\rho_w, \vec{f_w}, S_w$: Voir partie eau
  2. $\rho_a, \vec{f_a}, S_a$ : Voir partie air
  3. Les conductivités thermiques $\Gamma_w,\Gamma_a$ et $\Gamma_s$ : \[\begin{array}{l}\Gamma_w(T) = \Gamma_{w,0} + \Gamma_{w,0} \gamma_w^T (T-T_0)\\ \Gamma_a(T) = \Gamma_{a,0} + \Gamma_{a,0} \gamma_a^T (T-T_0)\\ \Gamma_s(T) = \Gamma_{s,0} + \Gamma_{s,0} \gamma_s^T (T-T_0) \end{array}\]
  4. Les chaleurs spécifiques $c_{p,w},c_{p,a}$ et $c_{p,s}$ :\[\begin{array}{l} c_{p,w}(T) = c_{p,w0} + c_{p,w0} H_w^T (T-T_0) \\ c_{p,a}(T) = c_{p,a0} + c_{p,a0} H_a^T (T-T_0)\\ c_{p,s}(T) = c_{p,s0} + c_{p,s0} H_s^T (T-T_0)\end{array}\]

Files

Prepro: LWAVAT2.F

Availability

Plane stress state NO
Plane strain state YES
Axisymmetric state YES
3D state YES
Generalized plane state NO

Input file

Parameters defining the type of constitutive law

Line 1 (2I5, 60A1)
ILLaw number
ITYPE 182
COMMENT Any comment (up to 60 characters) that will be reproduced on the output listing

Integer parameters

Line 1 (18I5)
IANI= 0, isotropic permeability case
≠ 0, anisotropic permeability case
IKWformulation index for $k_w$
IKAformulation index for $k_a$
ISRWformulation index for $S_W$
ITHERMformulation 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$

Real parameters: permeabilities definition

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)
PERMEsoil isotropic intrinsic permeability ($k_f$)

Real parameters

Line 1 (5G10.0)
POROSsoil porosity $(= n)$
TORTUsoil tortuosity $(=\tau )$
T0definition temperature $(=T_0)\ \left[K\right]$
PW0definition liquid pression $(=p_{w,0})\ \left[Pa\right]$
PA0definition gaz pression $(=p_{a,0})\ \left[Pa\right]$
Line 2 (7G10.0)
VISCW0liquid dynamic viscosity $(=\mu_{w,0}\ \left[ Pa.s \right]$
ALPHW0liquid dynamic viscosity thermal coefficient $(=\alpha_{w}^T)\ \left[K^{-1}\right]$
RHOW0
UXHIW0liquid compressibility coefficient $(=1/ \chi_{w})\ \left[ Pa^{-1}\right]$
BETAW0liquid thermal expansion coefficient $(=\beta{w}^T\ \left[K^{-1}\right]$
CONW0liquid thermal conductivity $(=\Gamma_{w,0})\ \left[W.m^{-1}.K^{-1}\right]$
GAMW0liquid thermal conductivity coefficient $(=\gamma_{w}^T\ \left[K^{-1}\right]$
Line 3 (3G10.0)
CPW0liquid specific heat $(=c_{p,w0})\ \left[J.kg^{-1}.K^{-1}\right]$
HEATW0liquid specific heat coefficient $(=H_{w}^T)\ \left[K^{-1}\right]$
EMMAGstorage coefficient $(=E_s)\ \left[ P_a^{-1}\right]$
Line 4 (7G10.0)
VISCA0
ALPHW0gaz dynamic viscosity thermal coefficient$(=\alpha_{a}^T)\ \left[K^{-1}\right]$
RHOA0gaz density$(=\rho_{a,0})\ \left[kg.m^{-3}\right]$
CONA0gaz thermal conductivity $(=\Gamma_{a,0})\ \left[W.m^{-1}.K^{-1}\right]$
GAMA0gaz thermal conductivity coefficient $(=\gamma_{a}^T)\ \left[K^{-1}\right]$
CPA0gaz specific heat $(=c_{p,a0})\ \left[J.kg^{-1}.K^{-1}\right]$
HEATA0gaz specific heat coefficient $(=H_{a}^T)\ \left[K^{-1}\right]$
Line 5 (5G10.0)
BETAS0solid thermal expansion coefficient $(=\beta_{s}^T)\ \left[K^{-1}\right]$
CONS0solid thermal conduction$(=\Gamma_{s,0})\ \left[W.m^{-1}.K^{-1}\right]$
GAMS0solid conduction coefficient $(=\gamma_{s}^T)\ \left[K^{-1}\right]$
CPS0solid specific heat $(=c_{p,s0})\ \left[J.kg^{-1}. K^{-1}\right]$
HEATS0solid specific heat coefficient $(=H_{s}^T)\ \left[ K^{-1}\right]$
Line 6 (3G10.0)
CKW11st coefficient of the function $k_{rw}$
CKW22nd coefficient of the function $k_{rw}$
CKW33rd coefficient of the function $k_{rw}$
Line 7 (3G10.0)
CKA11st coefficient of the function $k_{ra}$
CKA22nd coefficient of the function $k_{ra}$
CSR55th coefficient of the function $S_w$
Line 8 (7G10.0/)
CSR11st coefficient of the function $S_w$
CSR22nd coefficient of the function $S_w$
CSR33rd coefficient of the function $S_w$
CSR44th coefficient of the function $S_w$
SRESresidual saturation degree $(=S_{res})$
SRFIELDfield saturation degree $(=S_{r, field})$
AIREVair entry value $\left[Pa\right]$
Line 9 (7G10.0)
CLT11st coefficient of the function $\Gamma_T$
CLT22nd coefficient of the function $\Gamma_T$
CLT33rd coefficient of the function $\Gamma_T$
CLT44th coefficient of the function $\Gamma_T$
RHOCcoefficient for enthalpie $\rho C_p$ (if ienth = 1)
RHOC11st coefficient for enthalpie $\rho C_p$ (if ienth = 2)
RHOC22nd coefficient for enthalpie $\rho C_p$ (if ienth = 2)
Line 10 (4G10.0)
KRMINMinimum value of $kr$
HENRYHenry coefficient
EXPMfor IKRN=1: $m$ Exponent of Kozeni-Karmann formulation
for IKRN=3 : A parameter
EXPNn 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

Stresses

Number of stresses

24

Meaning

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

State variables

Number of state variables

= 26 in 2D cases
= 16 in 3D cases

List of state variables

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)
1)
Corman, G. (2024). Hydro-mechanical modelling of gas transport processes in clay host rocks in the context of a nuclear waste repository. PhD thesis, University of Liège. https://hdl.handle.net/2268/307996