Sidebar
This is an old revision of the document!
Table of Contents
HMIC
Description
2D hydraulic microscopic law for solid elements.
Can be parallelised in ELEMB (at the macro-scale) or in the perturbation loop (at the micro-scale).
The model
This law replaces the macroscopic fluid law, by considering a complete hydraulic microstructure, made of dominant horizontal bedding planes, vertical bridging planes and matrix blocks. Under the assumption of the spatial repeticion of the microstructure over the distance $w$, a Representative Element Volume (REV) is built, including fractures and tubes whose behaviours are governed by constitutive laws. Fluid pressures and fluxes are computed at the microscopic scale in that hydraulic network. This way, the law is used for water seepage, air seepage, diffusion and advection (coupled) under non-linear analysis in 2D porous media. Effects of mechanics on the flow are implicitely integrated into the microscale model by means of hydro-mechanical couplings.
Mass balance equation for water
\[ \underbrace{\frac{\partial}{\partial t} (\rho_s . n . S_{r,w}) + div(\rho_w \vec{q_l})}_{\text{Liquid water}} = 0 \]
Liquid water flow
From Darcy's law, the advective component of the liquid water flow respectively reads for a fracture and a tube:
\[
\vec{q_l} = - \frac{k_{r_w}}{\mu_w}\frac{1}{A}\kappa\left[ \vec{grad}(p_w) + g \rho_w \vec{grad}(y)\right]\]
where
\[
\kappa
= {}
\begin{cases}
-\frac{h_b^2}{12}h_b \cdot w, fracture\\
-\pi \frac{D^4}{128}, tube\\
\end{cases}
\]
\[
k_{r_w}
=
\begin{cases}
\frac{S_{r}^{*^2}}{2}(3-S_{r}^{*}), fracture\\
S_{r}^{*^2}, tube
\end{cases}
\]
Liquid state equations
- Density $\rho_w$: \[\rho_w (T, p_w) = \rho_{wo}\left[ 1+\frac{p_w-p_{w0}}{\chi_w} \right]\]
- Intrinsic permeability $k_w$:
Depending on the water saturation degree $S_w$ : $k_{r,w} = f(S_w)$ avec $k_{w,eff} = k_f k_{r,w}$ - Saturation degree $S_w$:
Depending on succion $s = p_a - p_w : S_w = f(s)$
Mass balance equation for air
\[\frac{\partial}{\partial t} (\rho_a . n . S_{r,g}) + div(\rho_a \vec{q_g}) + div(\vec{i_a}) = 0\]
Dry air flow
From Darcy's law, the advective component of the liquid water flow respectively reads for a fracture and a tube:
\[
\vec{q_g} = - \frac{k_{r_g}}{\mu_g}\frac{1}{A}\kappa\left[ \vec{grad}(p_g) + g \rho_g \vec{grad}(y)\right]\]
where
\[
\kappa
= {}
\begin{cases}
-\frac{h_b^2}{12}h_b \cdot w, fracture\\
-\pi \frac{D^4}{128}, tube\\
\end{cases}
\]
\[
k_{r_g}
=
\begin{cases}
(1-S_{r}^*)^3, fracture\\
(1-S_{r}^*)^2, tube
\end{cases}
\]
From Fick's law, the diffusive component of the dissolved air flow respectively reads for a fracture and a tube: \[ \vec{i}_a = - n S_{r,g} \tau D \rho_g \vec{grad} (\omega_a) \]
where $\omega_a = \rho_a/\rho_g$.
Dry gas state equations
- Density $\rho_a$ :
Assumption of classical ideal gas equation of state: \[\rho_a (T, p_a) = \rho_{a,0}\frac{p_a}{p_{a,0}}\frac{T_0}{T} \] - Perméabilité intrinsèque $k_g$:
Depending on the saturation degree $S_g$ : $k_{r,g} = f(S_g)$ avec $k_{g,effectif} = k_{g, intrinsic}k_{a,w}$ - Gaseous saturation degree $S_g$:
Depending on suction $s = p_g - p_w$
$S_g = 1-S_w$
Files
Prepro: LWAVAT.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) | |
|---|---|
| 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 |
Integer parameters
| 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$ | |
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) | |
|---|---|
| PERME | soil isotropic intrinsic permeability ($k_f$) |
Real parameters
| 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
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) |
