3. Constitutive equations

3.1. Conservation equations

This is just a reminder, how to set them up is shown in [R7.01.10] [4].

3.1.1. Mechanical balance

Noting \(\sigma\) the total mechanical stress tensor and \(r\) the homogenized density of the medium, the mechanical balance is written:

\(\text{Div}(\sigma )\text{+}r{F}^{m}\text{=}0\) eq 3.1.1-1

We recall that \(r\) is linked to fluid mass variations by the relationship:

\(r={r}_{0}\text{+}{m}_{w}\text{+}{m}_{\text{ad}}\text{+}{m}_{\text{vp}}\text{+}{m}_{\text{as}}\) eq 3.1.1-2

3.1.2. Conservation of fluid masses

For the fluid the derivative \(\dot{a}=\frac{{d}^{s}a}{\text{dt}}\) is in fact an Eulerian derivative and the equations we write for the fluid include transport terms, even if they may be hidden by the choice of unknowns. The equations for the conservation of fluid masses can then be written as:

\(\{\begin{array}{}\dot{{m}_{w}}+{\dot{m}}_{\text{vp}}+\text{Div}({M}_{w}+{M}_{\text{vp}})=0\\ \dot{{m}_{\text{as}}}+{\dot{m}}_{\text{ad}}\text{+}\text{Div}({M}_{\text{as}}\text{+}{M}_{\text{ad}})=0\end{array}\) eq 3.1.2-1

3.1.3. Conservation of energy: thermal equation

\(\begin{array}{c}{h}_{w}^{m}\dot{{m}_{w}}\text{+}{h}_{\text{ad}}^{m}{\dot{m}}_{\text{ad}}\text{+}{h}_{\text{vp}}^{m}{\dot{m}}_{\text{vp}}+{h}_{\text{as}}^{m}{\dot{m}}_{\text{as}}\text{+}\dot{Q}\text{'}+\text{Div}({h}_{w}^{m}{\mathrm{M}}_{\mathrm{w}}\text{+}{h}_{\text{ad}}^{m}{\mathrm{M}}_{\text{ad}}\text{+}{h}_{\text{vp}}^{m}{\mathrm{M}}_{\text{vp}}\text{+}{h}_{\text{as}}^{m}{\mathrm{M}}_{\text{as}})\text{+}\text{Div}(\mathrm{q})\text{=}\\ ({\mathrm{M}}_{\mathrm{w}}\text{+}{\mathrm{M}}_{\text{ad}}\text{+}{\mathrm{M}}_{\text{vp}}+{\mathrm{M}}_{\text{as}}){F}^{m}\text{+}\Theta \end{array}\) eq 3.1.3-1

3.2. Behavioral equations

3.2.1. Evolution of porosity

3.2.1.1. Isotropic general case

\(d\mathrm{\varphi }=\left(b-\mathrm{\varphi }\right)\left(d{\mathrm{\epsilon }}_{V}-3{\mathrm{\alpha }}_{0}\text{dT}+\frac{{\text{dp}}_{\text{gz}}-{S}_{\text{lq}}{\text{dp}}_{c}}{{K}_{s}}\right)\) eq 3.2.1-1

In this equation, the coefficients \(b\) and \({K}_{s}\) appear. \(b\) is the Biot coefficient and \({K}_{s}\) is the compressibility module of solid grains. If \({K}_{0}\) designates the « drained » compressibility module of the porous medium, we have the relationship:

\(b=1-\frac{{K}_{0}}{{K}_{s}}\) eq 3.2.1-2

Note:

The expression 3.2.1-1 may seem unusual considering the standard definition of the Biot coefficient. This is due to the fact that we use Eulerian porosity \(\phi\) while the usual definition of the Biot coefficient is based on the Lagrangian \(\Phi\) definition of this quantity.

The two definitions are linked by the relationship:

\(\mathrm{\Phi }=\left(1\text{+}{\mathrm{\epsilon }}_{V}\right)\mathrm{\varphi }\) eq 3.2.1-3

In the case of an isothermal evolution for a saturated medium, the variation in Lagrangian porosity is simply proportional to the variation in volume:

\(d\Phi =b\text{.}d{\varepsilon }_{V}\) eq 3.2.1-4

Representing eq 3.2.1-4 in eq 3.2.1-3, we find:

\(\left(1+{\mathrm{\epsilon }}_{V}\right)d\mathrm{\varphi }=\left(b-\mathrm{\varphi }\right)d{\mathrm{\epsilon }}_{V}\)

For \({\mathrm{\varepsilon }}_{V}\) small (hypothesis retained in THM models), we obtain \(d\varphi \text{=}(b\mathrm{-}\varphi )d{\varepsilon }_{V}\), which corresponds well to [éq 3.2.1-1] for the case examined.

3.2.1.2. Transverse isotropic case

\(d\phi \text{=}B\mathrm{:}d\epsilon \text{-}\phi d{\epsilon }_{V}\text{-}3{\alpha }_{\phi }\text{dT}\text{+}\frac{{\text{dp}}_{\text{gz}}\text{-}{S}_{\text{lq}}{\text{dp}}_{c}}{{M}_{\phi }}\) eq 3.2.1-5

In this equation, we see the Biot tensor \(B\), the Biot module of the solid matrix \(\frac{1}{{M}_{\phi }}\) and \({\alpha }_{\phi }\) the differential expansion coefficient appear.

The expression for the differential expansion coefficient is given by the relationship:

\({\alpha }_{\phi }\text{=}\frac{(B\text{-}\phi \delta )\mathrm{:}{\alpha }_{0}}{3}\) eq 3.2.1-6

• \(\delta\) the identity matrix,

• \(B\text{=}({b}_{L}({e}_{L}\otimes {e}_{L}\text{+}{e}_{T}\otimes {e}_{T})\text{+}{b}_{N}({e}_{N}\otimes {e}_{N}))\) the Biot tensor as a function of the average Biot coefficients \({b}_{L}\) and \({b}_{N}\) in the directions N and L of the local orthotropy coordinate system (L, T, N). In the isotropic case \({b}_{L}={b}_{N}=b\),

• \({\alpha }_{0}\text{=}({\alpha }_{L}({e}_{L}\mathrm{\otimes }{e}_{L}\text{+}{e}_{T}\mathrm{\otimes }{e}_{T})\text{+}{\alpha }_{N}({e}_{N}\mathrm{\otimes }{e}_{N}))\) the thermal expansion tensor depends on \({\alpha }_{N}\) and \({\alpha }_{L}\) the mean expansion coefficients of the porous medium and solid grains in the N and L directions of the local orthotropy coordinate system (L, T, N) .In the isotropic case \({\alpha }_{L}\mathrm{=}{\alpha }_{N}\mathrm{=}{\alpha }_{0}\).

Note:

In the isotropic case the expression for the differential expansion coefficient is given by:

\({\alpha }_{\varphi }\text{=}(b\text{-}\varphi ){\alpha }_{0}\)

The expression for the Biot module of the solid matrix is given by the following relationship ([11] and [13]):

\(\frac{1}{{M}_{\varphi }}\text{=}(B\text{-}\varphi \delta )\mathrm{:}{S}_{0}^{S}\mathrm{:}\delta\) eq 3.2.1-7

Note:

In the isotropic case, the expression of the Biot modulus of the solid matrix is given by:

\(\frac{1}{{M}_{\varphi }}\text{=}3(b\text{-}\varphi )(\frac{1\mathrm{-}2{\nu }^{S}}{{E}^{S}})\)

• \({S}_{0}^{S}\) the skeletal flexibility matrix, a function of the Young’s modulus of the solid matrix \({E}^{S}\) and the Poisson’s ratio of the solid matrix \({\nu }^{S}\). Assuming that microscopically the skeleton is isotropic homogeneous [12], the expression of the flexibility matrix is given by:

\({S}_{0}^{S}\text{=}(\begin{array}{cccccc}\frac{1}{{E}^{S}}& \text{-}\frac{{\nu }^{S}}{{E}^{S}}& \text{-}\frac{{\nu }^{S}}{{E}^{S}}& 0& 0& 0\\ \text{-}\frac{{\nu }^{S}}{{E}^{S}}& \frac{1}{{E}^{S}}& \text{-}\frac{{\nu }^{S}}{{E}^{S}}& 0& 0& 0\\ \text{-}\frac{{\nu }^{S}}{{E}^{S}}& \text{-}\frac{{\nu }^{S}}{{E}^{S}}& \frac{1}{{E}^{S}}& 0& 0& 0\\ 0& 0& 0& 2\frac{(1+{\nu }^{S})}{{E}^{S}}& 0& 0\\ 0& 0& 0& 0& 2\frac{(1+{\nu }^{S})}{{E}^{S}}& 0\\ 0& 0& 0& 0& 0& 2\frac{(1+{\nu }^{S})}{{E}^{S}}\end{array})\) eq 3.2.1-8

In practice, it is very difficult to access the microscopic parameters \({E}^{S}\) and \({\nu }^{S}\) characterizing the solid matrix. However, we can avoid knowing these two parameters by deducing the compressibility coefficients of solid grains from the elastic parameters of the porous medium.

Taking into account what was said above, we chose to set the value of \({\nu }^{S}\mathrm{=}\mathrm{0,3}\) arbitrarily during programming (see [12]). The following approach is then adopted in order to determine the flexibility tensor of the solid matrix:

• it is assumed that at the microscopic scale, the solid matrix is homogeneous. In addition, it is considered that the volume deformation of the solid matrix between the initial instant and the final instant (after mechanical loading) is negligible. Thus, the definition of the compressibility module of solid grains is given by:

\({K}_{S}=\frac{{E}^{S}}{3\left(1-2{\nu }^{S}\right)}\)

• the expressions of the components of the Biot tensor, along the axes \(L\) and \(N\) of the local orthotropy coordinate system are given by [12]:

\(\mathrm{\{}\begin{array}{c}{b}_{L}\mathrm{=}1\mathrm{-}\frac{{M}_{11}+{M}_{12}+{M}_{13}}{3{K}_{S}}\\ {b}_{N}\mathrm{=}1\mathrm{-}\frac{2{M}_{13}+{M}_{33}}{3{K}_{S}}\end{array}\)

with:

\({M}_{11}\mathrm{=}\frac{{E}_{L}({E}_{N}\mathrm{-}{E}_{L}{\nu }_{\text{LN}}^{2})}{(1+{\nu }_{\text{LT}})({E}_{N}\mathrm{-}{E}_{N}{\nu }_{\text{LT}}\mathrm{-}2{E}_{L}{\nu }_{\text{LN}}^{2})}\)

\({M}_{12}=\frac{{E}_{L}({E}_{N}{\nu }_{\text{LT}}+{E}_{N}{\nu }_{\text{LN}}{\nu }_{\text{LT}}+{E}_{L}{\nu }_{\text{LN}}^{2})}{(1+{\nu }_{\text{LT}})({E}_{N}-{E}_{N}{\nu }_{\text{LT}}-2{E}_{L}{\nu }_{\text{LN}}^{2})}\)

\({M}_{13}\mathrm{=}\frac{{E}_{L}{E}_{N}{\nu }_{\text{LN}}}{(1+{\nu }_{\text{LT}})({E}_{N}\mathrm{-}{E}_{N}{\nu }_{\text{LT}}\mathrm{-}2{E}_{L}{\nu }_{\text{LN}}^{2})}\)

\({M}_{33}\mathrm{=}\frac{{E}_{L}^{2}(1\mathrm{-}{\nu }_{\text{LT}})}{{E}_{N}\mathrm{-}{E}_{N}{\nu }_{\text{LT}}\mathrm{-}2{E}_{L}{\nu }_{\text{LN}}^{2}}\)

By considering either of the two expressions above (taking into account the isotropy hypothesis made on \({K}_{S}\)), we are then in a position to isolate \({K}_{S}\).

• from the expression for \({K}_{S}\) determined previously, we can have access to the Young’s modulus \({E}^{S}\) of the solid matrix knowing that:

\(\{\begin{array}{c}{\nu }_{S}=\mathrm{0,3}\\ {E}_{S}=3\left(1-2{\nu }^{S}\right){K}^{S}\end{array}\) eq 3.2.1-9

3.2.1.3. Case without mechanical coupling

For purely hydraulic or thermo-hydraulic models without mechanical coupling, it is however possible to vary the porosity via a storage coefficient \({E}_{m}\). The latter then relates the variation in porosity to the variation in liquid pressure such as:

\(d\varphi ={E}_{\mathrm{m.}}{\mathrm{dp}}_{\text{lq}}\)

This coefficient is not taken into account in the case of mechanical models (HM)

3.2.2. Evolution of fluid mass inputs

Using the definition of fluid mass inputs and using purely geometric arguments, we find:

\(\begin{array}{}{m}_{w}={\rho }_{w}(1+{\varepsilon }_{V})\varphi {S}_{\text{lq}}-{\rho }_{w}^{0}{\varphi }^{0}{S}_{\text{lq}}^{0}\\ {m}_{\text{ad}}={\rho }_{\text{ad}}(1+{\varepsilon }_{V})\varphi {S}_{\text{lq}}-{\rho }_{\text{ad}}^{0}{\varphi }^{0}{S}_{\text{lq}}^{0}\\ {m}_{\text{as}}={\rho }_{\text{as}}(1+{\varepsilon }_{V})\varphi (1-{S}_{\text{lq}})-{\rho }_{\text{as}}^{0}{\varphi }^{0}(1-{S}_{\text{lq}}^{0})\\ {m}_{\text{vp}}={\rho }_{\text{vp}}(1+{\varepsilon }_{V})\varphi (1-{S}_{\text{lq}})-{\rho }_{\text{vp}}^{0}{\varphi }^{0}(1-{S}_{\text{lq}}^{0})\end{array}\) eq 3.2.2-1

**Note:*

If we used Lagrangian porosity, the mass inputs would be written as:

\(\begin{array}{}{m}_{w}={\rho }_{w}\varphi {S}_{\text{lq}}-{\rho }_{w}^{0}{\varphi }^{0}{S}_{\text{lq}}^{0}\\ {m}_{\text{ad}}={\rho }_{\text{ad}}\varphi {S}_{\text{lq}}-{\rho }_{\text{ad}}^{0}{\varphi }^{0}{S}_{\text{lq}}^{0}\\ {m}_{\text{as}}={\rho }_{\text{as}}\varphi (1-{S}_{\text{lq}})-{\rho }_{\text{as}}^{0}{\varphi }^{0}(1-{S}_{\text{lq}}^{0})\\ {m}_{\text{vp}}={\rho }_{\text{vp}}\varphi (1-{S}_{\text{lq}})-{\rho }_{\text{vp}}^{0}{\varphi }^{0}(1-{S}_{\text{lq}}^{0})\end{array}\)

As an example, we demonstrate the first relationship in the saturated case \({S}_{\text{lq}}\text{=}1\) (with \({\mathrm{\rho }}_{\text{lq}}\text{=}{\mathrm{\rho }}_{w}\)).

Let it be an elementary domain of porous medium with volume \(O\). We note \({O}_{s}\) the volume occupied by the solid grains and \({O}_{l}\) the volume occupied by the liquid and the gas. We note \({O}^{0}\), \({O}_{{}_{s}}^{0}\), \({O}_{{}_{l}}^{0}\) the same volumes in the initial state. We recall that \({\mathrm{\varepsilon }}_{V}\) notes the variation in volume of the porous medium and we note \({\mathrm{\varepsilon }}_{{V}_{s}}\) the volume variation of solid grains.

By definition, we have: \(\phi \text{=}\frac{{\Omega }_{l}}{\Omega }\)

\({\Omega }_{s}=\Omega -{\Omega }_{l}=\Omega (1-\varphi )={\Omega }_{s}^{0}(1+{\epsilon }_{\text{Vs}})\)

But \(\Omega (1-\varphi )={\Omega }^{0}(1+{\epsilon }_{V})(1-\varphi )\).

From this we deduce:

\({\Omega }^{0}(1+{\varepsilon }_{V})(1-\varphi )={\Omega }_{s}^{0}(1+{\varepsilon }_{\text{Vs}})\)

All you have to do is write \({\Omega }_{s}^{0}={\Omega }^{0}(1-{\varphi }^{0})\) to get:

\({\Omega }^{0}(1+{\varepsilon }_{V})(1-\varphi )={\Omega }^{0}(1-{\varphi }^{0})(1+{\varepsilon }_{\text{Vs}})\)

From which we deduce:

\({\varepsilon }_{\text{Vs}}(1-{\varphi }^{0})={\varepsilon }_{V}(1-\varphi )-(\varphi -{\varphi }^{0})\)

We use the Eulerian definition of homogenized density \(r\text{'}\) (not to be confused with the Lagrangian definition of equation [éq 3.1.1-2]):

\(r\text{'}={\rho }_{(s)}(1-\varphi )+{\rho }_{(\mathrm{lq})}\varphi\)

and the definition of liquid mass intake:

\({r}^{\text{'}}\Omega \text{=}({r}_{0}\text{+}{m}_{\text{lq}}){\Omega }^{0}\)

We get:

\({\rho }_{s}(1-\varphi )\Omega +{\rho }_{\text{lq}}\varphi \Omega ={\rho }_{s}^{0}(1-{\varphi }^{0}){\Omega }^{0}+{\rho }_{\text{lq}}^{0}{\varphi }^{0}{\Omega }^{0}+{m}_{\text{lq}}{\Omega }^{0}\)

or again:

\({\rho }_{s}{\Omega }_{s}+{\rho }_{\text{lq}}\varphi (1\text{+}{e}_{V}){\Omega }^{0}={\rho }_{{}_{s}}^{0}{\Omega }_{s}^{0}+{\rho }_{{}_{\text{lq}}}^{0}{\varphi }^{0}{\Omega }^{0}\text{+}{m}_{\text{lq}}{\Omega }^{0}\)

Using the conservation of the mass of solid grains: \({\rho }_{s}{\Omega }_{s}={\rho }_{s}^{0}{\Omega }_{s}^{0}\) we finally obtain:

\({\rho }_{\text{lq}}\varphi (1\text{+}{\varepsilon }_{V})={\rho }_{{}_{\text{lq}}}^{0}{\varphi }^{0}\text{+}{m}_{\text{lq}}\)

3.2.3. Laws of fluid behavior

3.2.3.1. Liquid

\(\frac{d{\rho }_{w}}{{\rho }_{w}}=\frac{{\text{dp}}_{w}}{{K}_{w}}-3{\alpha }_{w}\text{dT}\) eq 3.2.3.1-1

We see the appearance of the water compressibility module \({K}_{w}\) and its expansion module \({\alpha }_{w}\).

3.2.3.2. Gaz

For the gas behavior equations, we take the ideal gas law:

\(\frac{{p}_{\text{vp}}}{{\mathrm{\rho }}_{\text{vp}}}\text{=}\frac{R}{{M}_{\text{vp}}^{\text{ol}}}T\) eq 3.2.3.2-1

\(\frac{{p}_{\text{as}}}{{\mathrm{\rho }}_{\text{as}}}\text{=}\frac{R}{{M}_{\text{as}}^{\text{ol}}}T\) eq 3.2.3.2-2

We see the molar mass of steam, \({M}_{\text{vp}}^{\text{ol}}\), and that of dry air, \({M}_{\text{as}}^{\text{ol}}\), appear.

3.2.4. Evolution of enthalpies

3.2.4.1. Liquid enthalpy

\({\text{dh}}_{w}^{m}={C}_{w}^{p}\text{dT}+(1-3{\alpha }_{w}T)\frac{{\text{dp}}_{w}}{{\rho }_{w}}\) eq 3.2.4.1-1

We see the appearance of specific heat at the constant pressure of water: \({C}_{w}^{p}\).

By replacing the pressure of the liquid in this expression by its value as a function of the capillary pressure and the gas pressure, we have:

\({\text{dh}}_{w}^{m}=(1-3{\alpha }_{w}T)\frac{{\text{dp}}_{\text{gz}}-{\text{dp}}_{c}-{\text{dp}}_{\text{ad}}}{{\rho }_{w}}+{C}_{w}^{p}\text{dT}\) eq 3.2.4.1-2

Noting \({C}_{\text{ad}}^{p}\) the specific heat at constant pressure of dissolved air, we have:

\({\text{dh}}_{\text{ad}}^{m}\text{=}{C}_{\text{ad}}^{p}\text{dT}\) eq 3.2.4.1-3

3.2.4.2. Enthalpy of gases

\({\text{dh}}_{\text{vp}}^{m}={C}_{\text{vp}}^{p}\text{dT}\) eq 3.2.4.2-1

\({\text{dh}}_{\text{as}}^{m}={C}_{\text{as}}^{p}\text{dT}\) eq 3.2.4.2-2

We see the specific heat at constant pressure of dry air \({C}_{\text{as}}^{p}\) and that of steam \({C}_{\text{vp}}^{p}\) appear.

3.2.4.3. Heat supply excluding fluids

It is quantity \(\mathrm{\delta }{Q}^{\text{'}}\) which represents the heat received by the system excluding the enthalpy supply of fluids.

\(\delta {Q}^{\text{'}}\text{=}({C}_{0}\mathrm{:}{\alpha }_{0}\mathrm{:}\text{d}\epsilon )T\text{+}3{\alpha }_{\text{lq}}^{m}{\text{Tdp}}_{c}\text{-}\left(3{\alpha }_{\text{gz}}^{m}\text{+}3{\alpha }_{\text{lq}}^{m}\right){\text{Tdp}}_{\text{gz}}\text{+}{C}_{\epsilon }^{0}\text{dT}\) Eq 3.2.4.3-1

Several expansion coefficients appear: \({\alpha }_{\text{lq}}^{m},{\alpha }_{\text{gz}}^{m}\). The components of tensor \({\alpha }_{0}\) are material data. In the isotropic case, the components of the tensor are equal \(({\alpha }_{L}\text{=}{\alpha }_{N}\text{=}{\alpha }_{0})\) and correspond both to the expansion coefficients of the porous medium and to those of solid grains (which are necessarily equal in the theory that we present here).

We also see the \({C}_{0}\) tensor, corresponding to the skeleton’s Hooke matrix, appear. This tensor is a function of the mechanical parameters of the porous medium \(({E}_{L},{E}_{N},{\nu }_{\text{LT}},{\nu }_{\text{LN}},{G}_{\text{LN}})\) in the transverse isotropic case or \((E,\nu )\) in the isotropic case.

Note:

In the isotropic case, the expression for the supply of heat excluding fluid is given by:

\(\delta {Q}^{\text{'}}\text{=}(3{K}_{0}{\alpha }_{0}\text{d}{\epsilon }_{V})T\text{+}3{\alpha }_{\text{lq}}^{m}{\text{Tdp}}_{c}\text{-}\left(3{\alpha }_{\text{gz}}^{m}\text{+}3{\alpha }_{\text{lq}}^{m}\right){\text{Tdp}}_{\text{gz}}\text{+}{C}_{\epsilon }^{0}\text{dT}\)

with \({K}_{0}\) the coefficient of compressibility of the porous medium, a function of the elastic parameters of the porous medium \((E,\nu )\)

\({\mathrm{\alpha }}_{\text{lq}}^{m},{\mathrm{\alpha }}_{\text{gz}}^{m}\) are given by relationships:

\({\alpha }_{\text{gz}}^{m}\text{=}(1\text{-}{S}_{\text{lq}}){\alpha }_{\varphi }+\frac{\varphi (1\text{-}{S}_{\text{lq}})}{\mathrm{3T}}\) eq 3.2.4.3-2

\({\alpha }_{\text{lq}}^{m}\text{=}{S}_{\text{lq}}{\alpha }_{\phi }\text{+}{\alpha }_{\text{lq}}\phi {S}_{\text{lq}}\) eq 3.2.4.3-3

We also see the appearance in [éq 3.2.4.3-1] of the specific heat at constant deformation of the porous medium \({C}_{\varepsilon }^{0}\), which depends on the specific heat at constant stress of the porous medium \({C}_{\sigma }^{0}\) by the relationship:

\({C}_{\epsilon }^{0}\text{=}{C}_{\sigma }^{0}\text{-}T({C}_{0}\mathrm{:}{\alpha }_{0})\mathrm{:}{\alpha }_{0}\) eq 3.2.4.3-4

Note:

In the isotropic case, the expression of the specific heat at constant deformation of the porous medium is given by:

\({C}_{\epsilon }^{0}\text{=}{C}_{\sigma }^{0}\text{-}9T{K}_{0}{\alpha }_{0}^{2}\)

\({C}_{\sigma }^{0}\) is given by a law of mixture:

\({C}_{\sigma }^{0}\text{=}(1\text{-}\varphi ){\rho }_{s}{C}_{\sigma }^{s}\text{+}{S}_{\text{lq}}\varphi ({\rho }_{w}{C}_{w}^{p}\text{+}{\rho }_{\text{ad}}{C}_{\text{ad}}^{p})\text{+}(1\text{-}{S}_{\text{lq}})\varphi ({\rho }_{\text{vp}}{C}_{\text{vp}}^{p}\text{+}{\rho }_{\text{as}}{C}_{\text{as}}^{p})\) eq 3.2.4.3-5

where \({C}_{\sigma }^{s}\) represents the specific heat at constant stress of solid grains and \({\rho }_{s}\) represents the density of solid grains. For the calculation of \({\rho }_{s}\), the deformation of solid grains is neglected, so \({\rho }_{s}\) is confused with its initial value \({\rho }_{s}^{0}\), which is in fact calculated from the initial specific mass of the porous medium \({r}_{0}\) by the following formula for mixtures:

\((1\text{-}{\varphi }^{0}){\rho }_{{}_{s}}^{0}\text{=}{r}_{0}\text{-}({\rho }_{w}^{0}\text{+}{\rho }_{\text{ad}}^{0}){S}_{\text{lq}}^{0}{\varphi }^{0}\text{-}(1\text{-}{S}_{\text{lq}}^{0}){\varphi }^{0}({\rho }_{\text{vp}}^{0}\text{+}{\rho }_{\text{as}}^{0})\) eq 3.2.4.3-6

3.2.5. Diffusion laws (complementary laws)

3.2.5.1. Heat distribution

We take Fourrier’s classical law:

\(q\text{=}\text{-}{\lambda }^{T}\text{.}\nabla T\) eq 3.2.5.1-1

where we see the thermal conductivity tensor \({\lambda }^{T}\) appearing.

The conductivity tensor is a function of porosity, saturation, and temperature and is given as the product of three functions plus one constant:

\({\lambda }^{T}\mathrm{=}{\lambda }_{\varphi }^{T}(\phi )\text{.}{\lambda }_{S}^{T}({S}_{\text{lq}})\text{.}{\lambda }_{T}^{T}(T)+{\lambda }_{\text{cte}}^{T}\) eq 3.2.5.1-2

In the isotropic case, we have \({\lambda }^{T}\text{=}{\lambda }^{T}\text{.}1\), \({\lambda }_{T}^{T}(T)\text{=}{\lambda }_{T}^{T}(T)\text{.}1\) and \({\lambda }_{\mathit{cte}}^{T}\text{=}{\lambda }_{\mathit{cte}}^{T}\text{.}1\).

3.2.5.2. Diffusion of fluids

These are Darcy’s laws, to which we add Fick’s law in the presence of steam.

Darcy’s laws are written for gas and for liquid:

\(\frac{{\mathrm{M}}_{\text{lq}}}{{\rho }_{\text{lq}}}\text{=}{\lambda }_{{\text{lq}}_{}}^{H}(\text{-}\mathrm{\nabla }{p}_{\text{lq}}\text{+}{\rho }_{\text{lq}}{\mathrm{F}}^{m})\) eq 3.2.5.2-1

\(\frac{{\mathrm{M}}_{\text{gz}}}{{\rho }_{\text{gz}}}\text{=}{\lambda }_{{\text{gz}}_{}}^{H}(\text{-}\mathrm{\nabla }{p}_{\text{gz}}\text{+}{\rho }_{\text{gz}}{\mathrm{F}}^{m})\) eq 3.2.5.2-2

where we see the hydraulic conductivity tensors \({\lambda }_{{\text{lq}}_{}}^{H}\) and \({\lambda }_{{\text{gz}}_{}}^{H}\) appear for liquid and gas respectively. In isotropic cases, \({\lambda }_{{\text{lq}}_{}}^{H}\text{=}{\lambda }_{{\text{lq}}_{}}^{H}\text{.}1\) and \({\lambda }_{{\text{gz}}_{}}^{H}\text{=}{\lambda }_{{\text{gz}}_{}}^{H}\text{.}1\).

We make the approximation that \(\frac{{M}_{w}}{{\rho }_{w}}={\lambda }_{w}^{H}(-\nabla {p}_{\text{lq}}+{\rho }_{\text{lq}}{F}^{m})\).

Note:

In this expression of Darcy’s law, the differential acceleration of water is overlooked. In the case of very permeable and very porous environments subject to seismic loading, this may constitute a limit.

Diffusion in the gas mixture is given by Fick’s law thanks to the relationship:

\(\frac{{M}_{\text{vp}}}{{\rho }_{\text{vp}}}-\frac{{M}_{\text{as}}}{{\rho }_{\text{as}}}=-\frac{{D}_{\text{vp}}}{{C}_{\text{vp}}(1-{C}_{\text{vp}})}\nabla (\frac{{p}_{\text{vp}}}{{p}_{\text{gz}}})\) eq 3.2.5.2-3

where \({D}_{\text{vp}}\) is the Fick diffusion coefficient of the gas mixture (L2.T-1), \({F}_{\text{vp}}\) is subsequently noted as such that:

\({F}_{\text{vp}}=\frac{{D}_{\text{vp}}}{{C}_{\text{vp}}(1-{C}_{\text{vp}})}\) eq 3.2.5.2-4

And with

\({C}_{\text{vp}}\text{=}\frac{{p}_{\text{vp}}}{{p}_{\text{gz}}}\) eq 3.2.5.2-5

So we have:

\(\frac{{M}_{\text{vp}}}{{\rho }_{\text{vp}}}-\frac{{M}_{\text{as}}}{{\rho }_{\text{as}}}=-{F}_{\text{vp}}\nabla {C}_{\text{vp}}\) eq 3.2.5.2-6

In addition, we have:

\(\frac{{M}_{\text{gz}}}{{\rho }_{\text{gz}}}=(1-{C}_{\text{vp}})\frac{{M}_{\text{as}}}{{\rho }_{\text{as}}}+{C}_{\text{vp}}\frac{{M}_{\text{vp}}}{{\rho }_{\text{vp}}}\) eq 3.2.5.2-7

and:

\({\rho }_{\text{gz}}\text{=}{\rho }_{\text{vp}}\text{+}{\rho }_{\text{as}}\) eq 3.2.5.2-8

For the diffusion of the liquid mixture, the usual writing is as follows:

\({M}_{\text{ad}}-{M}_{w}=-{D}_{\text{ad}}\nabla {\rho }_{\text{ad}}\) eq 3.2.5.2-9

where \({D}_{\text{ad}}\) is the Fick diffusion coefficient of the liquid mixture. In order to keep the writing homogeneous with that of the gaseous mixture, we subsequently note \({F}_{\text{ad}}\) such that:

\({F}_{\text{ad}}\text{=}{D}_{\text{ad}}\) eq 3.2.5.2-10

And the concentration \({C}_{\text{ad}}\) corresponds here to the density of dissolved air:

\({C}_{\text{ad}}\text{=}{\rho }_{\text{ad}}\) eq 3.2.5.2-11

\({M}_{\text{ad}}-{M}_{w}=-{F}_{\text{ad}}\nabla {C}_{\text{ad}}\) eq 3.2.5.2-12

With regard to liquid, it was accepted that the liquid Darcy law applies to the speed of liquid water. There is therefore no need to define an average liquid speed.

\(\frac{{M}_{w}}{{\rho }_{w}}={\lambda }_{\text{lq}}^{H}(-\nabla {p}_{\text{lq}}+{\rho }_{\text{lq}}{F}^{m})\) eq 3.2.5.2-13

and:

\({\rho }_{\text{lq}}\text{=}{\rho }_{w}\text{+}{\rho }_{\text{ad}}\) eq 3.2.5.2-14

By combining these relationships, we then find:

\(\frac{{\mathrm{M}}_{\text{as}}}{{\rho }_{\text{as}}}\text{=}{\lambda }_{\text{gz}}^{H}(\text{-}\mathrm{\nabla }{p}_{\text{gz}}\text{+}{\rho }_{\text{gz}}{\mathrm{F}}^{m})\text{+}{C}_{\text{vp}}{F}_{\text{vp}}\mathrm{\nabla }{C}_{\text{vp}}\) eq 3.2.5.2-15

\(\frac{{\mathrm{M}}_{\text{vp}}}{{\rho }_{\text{vp}}}\text{=}{\lambda }_{{\text{gz}}_{}}^{H}(\text{-}\mathrm{\nabla }{p}_{\text{gz}}\text{+}{\rho }_{\text{gz}}{\mathrm{F}}^{m})\text{-}(1\text{-}{C}_{\text{vp}}){F}_{\text{vp}}\mathrm{\nabla }{C}_{\text{vp}}\) eq 3.2.5.2-16

\(\frac{{\mathrm{M}}_{w}}{{\rho }_{w}}\text{=}{\lambda }_{\text{lq}}^{H}(\text{-}\mathrm{\nabla }{p}_{\text{lq}}\text{+}{\rho }_{\text{lq}}{\mathrm{F}}^{m})\) eq 3.2.5.2-17

\({\mathrm{M}}_{\text{ad}}\text{=}{\rho }_{\text{ad}}{\lambda }_{\text{lq}}^{H}(\text{-}\mathrm{\nabla }{p}_{\text{lq}}\text{+}{\rho }_{\text{lq}}{\mathrm{F}}^{m})\text{-}{F}_{\text{ad}}\mathrm{\nabla }{C}_{\text{ad}}\) eq 3.2.5.2-18

Hydraulic conductivity tensors \({\lambda }_{\text{lq}}^{H}\) and \({\lambda }_{{\text{gz}}_{}}^{H}\) are not directly data and their value is known from the formulas:

\({\lambda }_{\text{lq}}^{H}\text{=}\frac{{K}^{\text{int}}(\varphi )\text{.}{k}_{\text{lq}}^{\text{rel}}({S}_{\text{lq}})}{{\mu }_{w}(T)}\) eq 3.2.5.2-19

\({\lambda }_{{\text{gz}}_{}}^{H}\text{=}\frac{{K}^{\text{int}}(\varphi )\text{.}{k}_{\text{gz}}^{\text{rel}}({S}_{\text{lq}},{p}_{\text{gz}})}{{\mu }_{\text{gz}}(T)}\) eq 3.2.5.2-20

\({K}^{\text{int}}\) is the intrinsic permeability tensor, characteristic of the porous medium and user data, any function of porosity. In the isotropic case, \({K}^{\text{int}}\text{=}{K}^{\text{int}}\text{.}1\);

\({\mu }_{w}\) is the dynamic viscosity of water, characteristic of water and user data, any function of temperature;

\({\mu }_{\text{gz}}\) is the dynamic viscosity of the gas, gas characteristic and user data, any function of the temperature;

\({k}_{\text{lq}}^{\text{rel}}\) is the relative permeability to the liquid, characteristic of the porous medium and user data, any function of liquid saturation;

\({k}_{\text{gz}}^{\text{rel}}\) is the relative permeability to gas, characteristic of the porous medium and user data, any function of liquid saturation and gas pressure.

Note:

The hydraulic conductivities defined here are not necessarily very familiar for soil mechanics, who commonly use permeability \(k\) for saturated environments, which is homogeneous at one speed.

The relationship between tensors \(k\) and \({\lambda }_{\text{lq}}^{H}\) is as follows: \({\lambda }_{\text{lq}}^{H}\text{=}\frac{k}{{\rho }_{w}g}\) where \(g\) is the acceleration of gravity.

The Fick diffusion coefficient of the gas mixture \({F}_{\text{vp}}\) is a characteristic of the porous medium, a user data that is a function of any one of the vapour pressure, the gas pressure, the saturation and the temperature which will be written as a function product of each of these variables: \({F}_{\text{vp}}({P}_{\text{vp}},{P}_{\text{gz}},T,S)\text{=}{f}_{\text{vp}}^{\text{vp}}({P}_{\text{vp}})\text{.}{f}_{\text{vp}}^{\text{gz}}({P}_{\text{gz}})\text{.}{f}_{\text{vp}}^{T}(T)\text{.}{f}_{\text{vp}}^{S}(S)\) the derivatives with respect to the vapour pressure and to the saturation will be neglected. In the same way for the Fick diffusion coefficient of the liquid medium: \({F}_{\text{ad}}({P}_{\text{ad}},{P}_{\text{lq}},T,S)\text{=}{f}_{\text{ad}}^{\text{ad}}({P}_{\text{ad}})\text{.}{f}_{\text{ad}}^{\text{lq}}({P}_{\text{lq}})\text{.}{f}_{\text{ad}}^{T}(T)\text{.}{f}_{\text{ad}}^{S}(S)\), only the derivative as a function of temperature is taken into account.

3.2.6. Water-steam balance

This relationship is essential and has the consequence of reducing the number of pressure unknowns.

We note \({h}_{w}^{m}\) the mass enthalpy of water, \({s}_{w}^{m}\) its entropy and \({g}_{w}^{m}={h}_{w}^{m}-{\text{Ts}}_{w}^{m}\) its free enthalpy.

We note \({h}_{\text{vp}}^{m}\) the mass enthalpy of steam, \({s}_{\text{vp}}^{m}\) its entropy and \({g}_{\text{vp}}^{m}={h}_{\text{vp}}^{m}-{\text{Ts}}_{\text{vp}}^{m}\) its free enthalpy.

The water-steam balance is written as:

\({g}_{\text{vp}}^{m}\text{=}{g}_{w}^{m}\) eq 3.2.6-1

Who gives:

\({h}_{\text{vp}}^{m}-{h}_{w}^{m}=T({s}_{\text{vp}}^{m}-{s}_{w}^{m})\) eq 3.2.6-2

In addition, the definition of free enthalpy tells us that: \(\text{dg}\text{=}\frac{\text{dp}}{\rho }-\text{sdT}\), which, when applied to steam and water, combined with the relationship \({\text{dg}}_{{}_{\text{vp}}}^{m}\text{=}{\text{dg}}_{{}_{w}}^{m}\) and using [éq 3.2.6-2] gives:

\(\frac{{\text{dp}}_{\text{vp}}}{{\rho }_{\text{vp}}}=\frac{{\text{dp}}_{w}}{{\rho }_{w}}+({h}_{\text{vp}}^{m}-{h}_{w}^{m})\frac{\text{dT}}{T}\) eq 3.2.6-3

This relationship can be expressed as a function of capillary pressure and gas pressure:

\({\text{dp}}_{\text{vp}}=\frac{{\rho }_{\text{vp}}}{{\rho }_{w}}({\text{dp}}_{\text{gz}}-{\text{dp}}_{c}-{\text{dp}}_{\text{ad}})+{\rho }_{\text{vp}}({h}_{\text{vp}}^{m}-{h}_{w}^{m})\frac{\text{dT}}{T}\) eq 3.2.6-4

3.2.7. Balance between dry air and dissolved air

Dissolved air is defined via the Henry constant \({K}_{H}\), which relates the molar concentration of dissolved air \({C}_{\text{ad}}^{\text{ol}}\) (\(\mathit{moles}\mathrm{/}{m}^{3}\)) to the pressure of dry air:

\({C}_{\text{ad}}^{\text{ol}}\text{=}\frac{{p}_{\text{as}}}{{K}_{H}}\) eq 3.2.7-1

with \({C}_{\text{ad}}^{\text{ol}}\text{=}\frac{{\rho }_{\text{ad}}}{{M}_{\text{ad}}^{\text{ol}}}\) eq 3.2.7-2

The molar mass of dissolved air, \({M}_{\text{ad}}^{\text{ol}}\), is logically the same as that of dry air \({M}_{\text{as}}^{\text{ol}}\). For dissolved air, we take the ideal gas law:

\(\frac{{p}_{\text{ad}}}{{\mathrm{\rho }}_{\text{ad}}}\text{=}\frac{R}{{M}_{\text{as}}^{\text{ol}}}T\) eq 3.2.7-3

The pressure of dissolved air is therefore connected to that of dry air by:

\({p}_{\text{ad}}\text{=}\frac{{p}_{\text{as}}}{{K}_{H}}\text{RT}\) eq 3.2.7-4

3.2.8. Mechanical behavior

We will write it in differential form:

\(d\sigma \text{=}d\sigma \text{'}+d{\sigma }_{p}\) eq 3.2.8-1

In the general case (except for the “HYDR_TABBAL” law), using a Bishop [10] formulation extended to unsaturated media we write:

\(d{\mathrm{\sigma }}_{p}\text{=}\text{-}B({\mathit{dp}}_{\text{gz}}\text{-}{S}_{\text{lq}}{\mathit{dp}}_{c})\) eq 3.2.8-2

(which corresponds well in the saturated case to the classical formulation of Biot \(d{\mathrm{\sigma }}_{p}\text{=}\text{-}B{\mathit{dp}}_{\text{lq}}\)).

We therefore see the Biot tensor appearing in the formula [éq 3.2.8-2]. In the isotropic case \(B\text{=}b\text{.}1\) and \({\mathrm{\sigma }}_{p}\text{=}{\mathrm{\sigma }}_{p}\text{.}1\).

This formulation is valid for « reasonable » drying ranges, namely relative humidities (or degree of humidity) RH between 50% and 100%. It should be noted that relative humidity (RH) and capillary pressure are linked by Kelvin’s law (cf. section 3.3).

For relative humidity ranges falling below 50%, the model resulting from the thesis of Ginger El Tabbal [14] is recommended. This model is called by the HYDR_TABBAL keyword under RELATION_KIT. In addition to the capillary effects used in the above formulation, the effect of adsorption is taken into account here. In this case, the hydraulic stress is written in the form \(d{\mathrm{\sigma }}_{p}\text{=}\text{-}Bd\mathrm{\pi }\) with:

math:

`begin {array} {c} dmathrm {pi} = {mathit {dp}} _ {text {gz}} - {S} _ {mathit {BJH}}} {mathit {BJH}}} {mathit {dp}}} {mathit {dp}}} _ {text {c}}} {mathit {dp}}} _ {mathit {dp}}} _ {mathit {dp}}} _ {text {c}}} {mathit {dp}}} _ {mathit {dp}}} _ {text {c}}} {mathit {BJH}}) {mathit {dS}}} _ {mathit {BJH}}\ +frac {2} {3} ({A} _ {0}/{mathrm {phi}}}} ^ {phi}}}} _ {mathit {BJH}}/{mathrm {phi}}}/{mathrm {phi}}} ^ {0}) (- ({mathrm {omega}}} _ {mathit {BJH}}/{mathrm {phi}}}/{mathrm {phi}}} {phi}}} ^ {0}) (- ({mathrm {omega}}} _ {mathit {BJH}} {BJH}})mathit {dt} (mathit {HR}) +t (mathit {HR}) d {mathrm {omega}} _ {mathit {BJH}}} ({S}} ({S} _ {mathit {HR}}))mathrm {.}} {mathit {dp}} _ {text {c}} - {text {c}} - {left (frac {partialgamma}} {epsilon} _ {s}}right) |} _ {mu} d {mu} d {mathrm {c}}} - {mathrm {omega}}} _ {mathit {BJH}})end {array}} _ {mu} d {mathrm {omega}}} _ {mathit {omega}}} _ {mathit {BJH}}))end {array} `eq 3.2.8-3

With:math: {S} _ {mathit {BJH}}}, :math: {A} _ {0}, :math: {mathrm {omega}} _ {mathit {BJH}}}}, :math: {BJH}}} , :math: `t (mathit {HR}}), :math: {left (frac {partialgamma}} _ {mathit {BJH}}}}, :math: t (mathit {HR}}) `, :math: `t (mathit {HR}}), :math: {left (frac {partialgamma}} {partialgamma}} {partial}}}, :math: `tepsilon} _ {s}}right) |} _ {mu} `the model input data calculated beforehand. The definition is given here:

\({S}_{\mathit{BJH}}\) Volume fraction of liquid water as a function of RH (non-adsorbed water). This function is calculated using the BJH [16] method.

\({A}_{0}\) Total pore area (or specific surface area of the material) per unit volume m²/m3. It is the product of the dry density (g/m3) of the material by the specific surface area (m2/g) determined by method BET [17] or BJH [16].

\({\mathrm{\omega }}_{\mathit{BJH}}\) Surface fraction of unsaturated pores as a function of relative humidity (RH). This is the ratio between the surface area of the unsaturated pores (per unit volume) determined by method BJH (at a given relative humidity) and the total surface area of the pores \({A}_{0}\).

\(t(\mathit{HR})\) Thickness of the adsorbed water layer. This thickness is evaluated by an empirical relationship called Badmann [15].

math:

{left (frac {partialgamma} {partialgamma} {partial {epsilon} _ {s}}right) |} _ {mu} Shuttleworth term, strongly dependent on the orientation and morphology of the solid surface. This term is recalculated on the basis of drying shrinkage experience.

\({S}_{\mathit{BJH}}\), \({A}_{0}\), \({\mathrm{\omega }}_{\mathit{BJH}}\), are therefore calculated from experimental data by numerical methods called BET or BJH. The methodology for obtaining these parameters as well as the thickness of the adsorbed layer and the Shuttleworth parameter is described in detail in doc U2.04.05.

In the relationship [éq 3.2.8-1] the evolution of the effective stress tensor \({\sigma }^{\text{'}}\) is assumed to depend only on the displacement of the skeleton and on the internal variables \(\alpha\). The usual terms related to thermal deformation are integrated into the calculation of the effective stress:

\(d\sigma \text{'}=f(d\varepsilon -{\alpha }_{0}\text{dT}I,d\alpha )\) eq 3.2.8-4

The reason for this choice is to be able to use any classical thermomechanical law for the calculation of effective stresses, laws which, in most writings, are in accordance with [éq 3.2.8-3].

In practice, thermal deformation is evaluated by the following formula:

(3.1)\[ {\ mathrm {\ epsilon}}} ^ {\ mathit {th}}\ left (T\ right) = {\ mathrm {\ alpha}} _ {0}\ left (T- {T} _ {\ mathit {th}} _ {\ mathit {th}}\ right)\]

The reference temperature is given by THM_INIT/TEMPdans DEFI_MATERIAU.

3.2.9. Details on the terms of two-phase transfers

3.2.9.1. Sorption isotherm

To close the system, there is still a relationship to be written, connecting saturation and pressures. We chose to consider that liquid saturation was any function of capillary pressure, that this function was a characteristic of the porous medium and provided with data by the user.

Since the user can very well provide a function \({S}_{\text{lq}}({p}_{c})\) refined piecewise, and given that the derivative of this function, \(\frac{\partial {S}_{\text{lq}}}{\partial {p}_{c}}\), plays an essential physical role, we chose to ask the user to also provide this curve, leaving it up to him to ensure the consistency of the data thus specified. However, there is the possibility for the user to use an analytical saturation model and its derivative coded « hard » in the source: the Mualem-Van Genuchten model (see next section).

It should be noted that in the present approach, we are talking about a one-to-one relationship between saturation and capillary pressure. It is known that for most porous media, it is not the same relationship that must be used for drying paths and for hydration paths. This is one of the limitations of the current approach.

3.2.9.2. The Mualem - Van Genuchten model

Regarding hydraulic behavior, the user currently has two choices: to enter manually, and in tabulated form, relative permeabilities, saturation laws and their derivatives while ensuring their consistency (keyword HYDR_UTIL or HYDR_TABBAL under RELATION_KIT), or to use a known model programmed in analytical form: the Mualem — Van Genuchten model (keyword HYDR_VGM or HYDR_VGC under), or to use a known model programmed in analytical form: the Mualem — Van Genuchten model (keyword or under), or to use a known model programmed in analytical form: the Mualem — Van Genuchten model (keyword or under). RELATION_KIT).

Note:

There are of course other classical models to describe hydraulic behavior (Brook Corey for example) but they are currently not available in Aster. It is then up to the user to enter them in tabulated form in the command file. The Mualem- Van Genuchten model is a classical model for describing clay materials typical of underground storage problems. For relative gas permeability, it is common to use either a Van Genuchten version or a cubic version.

The Mualem-Van Genuchten model is translated into a capillary saturation/pressure law (Van Genuchten to which an inlet pressure is added) such that:

\({S}_{\text{we}}\text{=}\frac{1}{{\left[1\text{+}{\left(\frac{{P}_{c}-{P}_{e}}{{P}_{r}}\right)}^{n}\right]}^{m}}\)

with \({S}_{\text{we}}=\frac{S-{S}_{\text{wr}}}{1-{S}_{\text{wr}}}\) and \(m=1-\frac{1}{n}\)

\(n\), \({S}_{\text{wr}}\) (residual water saturation) and \({P}_{r}\) are model parameters entered by the user in his command file. \({P}_{e}\) is also a parameter (by default taken equal to 0) corresponding to an inlet pressure (the medium only desaturates if the capillary pressure is greater than this pressure).

Relative permeability to water is then expressed by integrating the prediction model proposed by Mualem (1976) into Van Genuchten’s capillarity model:

\({k}_{r}^{w}=\sqrt{{S}_{\text{we}}}{(1-{(1-{S}_{{\text{we}}^{1/m}})}^{m})}^{2}\)

The permeability to gas is formulated in a similar way in the case of model HYDR_VGM:

\({k}_{r}^{\text{gz}}=\sqrt{(1-{S}_{\text{we}})}{(1-{S}_{{\text{we}}^{1/m}})}^{2m}\)

or by a simple cubic law in the case of HYDR_VGC:

\({k}_{r}^{\text{gz}}={(1-S)}^{3}\)

In code_aster, we do a numerical treatment of this model in the « border » areas (\(S\mathrm{=}1\) or \(S\mathrm{=}0\)). For this we use two additional parameters corresponding to a treatment that is performed on these curves, \(\mathit{Smax}\) and \(\mathit{CSAT}\):

For \(S>\mathit{Smax}\), these curves are interpolated by a polynomial of degree 2 \(\mathit{C1}\) in \(\mathit{Smax}\), so as to avoid having to deal with derivatives of infinite values. Indeed, for \(S\mathrm{=}1\):

\(\frac{\partial {k}_{r}^{w}(S)}{\partial S}\text{=}\infty\)

And for the case HYDR_VGM \(\frac{\partial {k}_{r}^{\text{gz}}(S)}{\partial S}\text{=}\infty\)

To avoid then having to deal with this problem (which in principle has no physical meaning), these functions are replaced from a saturation \(\mathit{Smax}\) by a second-order polynomial \({C}^{1}\) at this point.

For function \({k}_{r}^{w}(S)\), this is as follows:

For \(S\mathrm{=}\mathit{Smax}\), we determine the polynomial \(\mathit{PL}(S)\) such that:

\(\{\begin{array}{}\text{PL}(S\text{max})\text{=}{k}_{r}^{w}(S\text{max})\\ P{L}^{\text{'}}(S\text{max})\text{=}\frac{\partial {k}_{r}^{w}}{\partial S}(S\text{max})\end{array}\) and \(\mathit{PL}(1)\mathrm{=}0\)

For \(S>\mathit{Smax}\), \({k}_{r}^{w}(S)\) is replaced by \(\mathit{PL}(S)\).

And for \({k}_{r}^{\text{gz}}(S)\) in the case HYDR_VGM:

For \(S\mathrm{=}\mathit{Smax}\), we determine the polynomial \(\mathit{PG}(S)\) such that:

\(\{\begin{array}{}\text{PG}(S\text{max})\text{=}{k}_{r}^{\mathrm{gz}}(S\text{max})\\ P{G}^{\text{'}}(S\text{max})\text{=}\frac{\partial {k}_{r}^{\mathrm{gz}}}{\partial S}(S\text{max})\end{array}\) and \(\mathit{PG}(0)\mathrm{=}1\)

For \(S>\mathit{Smax}\), \({k}_{r}^{\text{gz}}(S)\) is replaced by \(\mathit{PG}(S)\).

For suction \(S(\mathit{Pc})\) and for \(\mathit{Pc}<\mathit{Pcmin}\) (with \(S(\mathit{Pcmin})\mathrm{=}\mathit{Smax}\)) we extend the curve \(S(\mathit{Pc})\) by a hyperbola such that the curve is \({C}^{1}\) at this point:

For \(S>\mathit{Smax}\):

\(S(\text{Pc})=1-\frac{A}{B-\text{Pc}}\)

with \(A\) and \(B\) such that the curve is \(\mathit{C1}\) in \(\mathit{Smax}\).

So we do have a decreasing curve that tends to 1 when \(\mathit{Pc}\) tends to \(\mathrm{-}\mathrm{\infty }\). This treatment allows us to manage negative capillary pressures (in these areas the term capillary pressure is abusive, it is implicitly a change in variable making it possible to treat a nearly saturated problem).

\(S(\mathit{Pc})\) is then multiplied by a « safety » coefficient \(\mathit{CSAT}\) so that saturation never reaches 1 (a problem that we don’t know how to deal with). We recommend taking a value of \(\mathit{CSAT}\) very close to 1 (0.999999 for example).

3.2.10. Summary of material characteristics and user data

• the drained Young’s modules \({E}_{L}\) and \({E}_{N}\) as well as the drained Poisson coefficients \({\nu }_{\text{LT}}\)

and \({\nu }_{\text{LN}}\) which allow the calculation of the flexibility matrix of the porous medium. (transverse isotropic case),

• the drained Young’s modulus \(E\), and the drained Poisson’s ratio \(\nu\), which allow the calculation of the flexibility matrix of the porous medium. (isotropic case),

• the average Biot coefficients \({b}_{L}\) and \({b}_{N}\) which allow the calculation of the Biot module of the solid matrix \(\frac{1}{{M}_{\phi }}\) (transverse isotropic case),

• the Biot coefficient \(b\text{=}1\text{-}\frac{{K}_{0}}{{K}_{S}}\) makes it possible to calculate the compressibility module of solid grains \({K}_{S}\), useful in calculating the porosity for the isotropic case,

• the compressibility module of water \({K}_{w}\),

• the expansion coefficient of water \({\alpha }_{w}\),

• the ideal gas constant \(R\),

• the molar mass of steam \({M}_{\text{vp}}^{\text{ol}}\),

• the molar mass of dry air \({M}_{\text{as}}^{\text{ol}}\), (= \({M}_{\text{ad}}^{\text{ol}}\))

• the specific heat at constant pressure of water \({C}_{w}^{p}\),

• the specific heat at constant pressure of dissolved air \({C}_{\text{ad}}^{p}\)

• the specific heat at constant pressure of dry air \({C}_{\text{as}}^{p}\),

• the specific heat at constant pressure of steam \({C}_{\text{vp}}^{p}\),

• the average expansion coefficients of the porous medium \({\alpha }_{L}\) and \({\alpha }_{N}\), which are also those of solid grains. They allow the calculation of \({\alpha }_{\phi }\) (transverse isotropic case),

• the expansion coefficient of the porous medium \({\alpha }_{0}\), which is also that of solid grains. It allows the calculation of \({\alpha }_{\phi }\text{=}(b\text{-}\phi ){\alpha }_{0}\) (isotropic case),

• the specific heat at constant stress of solid grains \({C}_{\sigma }^{s}\),

• the thermal conductivity tensor of solid grains alone, \({\lambda }_{s}^{T}\), any function of temperature,

• the thermal conductivity tensor of the liquid, \({\lambda }_{\text{lq}}^{T}\), any function of the temperature,

• the thermal conductivity tensor of dry air, \({\lambda }_{\text{as}}^{T}\), any function of temperature,

• the Fick diffusion coefficient for the gas mixture, \({F}_{\text{vp}}\), a function of any temperature, gas pressure, vapor pressure, and saturation

• the Fick diffusion coefficient for the liquid mixture, \({F}_{\text{ad}}\), a function of any liquid temperature and pressure, dissolved air pressure, and saturation.

• Henry’s constant \({K}_{H}\) as a function of temperature,

• the intrinsic permeability tensor, \({K}^{\text{int}}\), any function of porosity,

• the dynamic viscosity of water, \({\mu }_{w}\), any function of temperature,

• the dynamic viscosity of the gas, \({\mu }_{\text{gz}}\), any function of the temperature,

• the relative permeability to the liquid, \({k}_{\text{lq}}^{\text{rel}}\), any function of the liquid saturation,

• the relative permeability to gas, \({k}_{\text{gz}}^{\text{rel}}\), a function of any liquid saturation and gas pressure,

• the saturation/capillary pressure relationship, \({S}_{\mathrm{lq}}({p}_{c})\), any function of capillary pressure,

• in general, the initial state is characterized by:

• the initial temperature,

• the initial pressures from which we deduce the initial saturation \({S}_{{}_{\text{lq}}}^{0}({p}_{{}_{c}}^{0})\),

• the initial specific mass of water \({\rho }_{w}^{0}\),

• the initial porosity \({\mathrm{\phi }}^{0}\),

• the initial pressure of steam \({p}_{\text{vp}}^{0}\) from which we deduce the initial density of steam \({\rho }_{\text{vp}}^{0}\),

• the initial pressure of dry air \({p}_{\text{as}}^{0}\) from which we deduce the initial density of dry air \({\rho }_{\text{as}}^{0}\).

• the initial homogenized density of the porous medium \({r}_{0}\),

• the initial enthalpies of water, dissolved air, steam, and dry air.

3.3. The reference state and the initial state

The introduction of initial conditions is very important, especially for enthalpies.

In practice, we can reason by considering that we have three states for fluids:

• the current state,

• the reference state: it is that of fluids in their free state. Very often the pressures of water and air will be taken to mean atmospheric pressure. In this reference state, it can be considered that the enthalpies are zero,

• the initial state: it is important to note that, in the initial state of the porous medium, water is in a hygroscopic state different from that of free water. For the enthalpies of water and steam we should take:

\(\begin{array}{c}{}^{\text{init}}{h}_{w}^{m}=\frac{{p}_{w}^{\text{init}}-{p}_{l}^{\text{ref}}}{{\rho }_{w}}=\frac{{p}_{w}^{\text{init}}-{p}_{\text{atm}}}{{\rho }_{w}}\\ {}^{\text{init}}{h}_{\text{vp}}^{m}=L({T}^{\text{init}})\text{=}\text{chaleur latente de vaporisation}\\ \begin{array}{}{}^{\text{init}}{h}_{\text{as}}^{m}=0\\ {}^{\text{init}}{h}_{\text{ad}}^{m}=0\end{array}\end{array}\)

Notes:

The initial vapour pressure should be taken in accordance with the**other data. Very often, we start from knowing an initial state of humidity. The humidity level is the ratio between the vapour pressure and the saturating vapour pressure at the temperature in question. We then use Kelvin’s law which gives the pressure of the liquid as a function of the vapor pressure, the temperature and the saturated vapor pressure: *:math:`frac{{p}_{w}-{p}_{w}^{text{ref}}}{{rho }_{w}}=frac{R}{{M}_{text{vp}}^{text{ol}}}Ttext{ln}(frac{{p}_{text{vp}}}{{p}_{text{vp}}^{text{sat}}(T)})`*. This relationship is only valid for isothermal evolutions. We emphasize that\({p}_{w}^{\text{ref}}`*corresponds to a state of equilibrium to which corresponds to*:math:`{p}_{\text{vp}}^{\text{sat}}`*, this state of equilibrium corresponds in fact to*:math:`{p}_{w}^{0}\text{=}{p}_{\text{gz}}^{0}\text{=}1`*ATM.For evolutions with temperature variation, knowing a law giving the saturated vapor pressure at temperature*:math:`{T}_{0}`*, for example: * *, for example:* :math:`{p}_{\text{vp}}^{\text{sat}}({T}_{0})={\text{10}}^{(2\text{.}\text{7858}+\frac{{T}_{0}-\text{273}\text{.}5}{\text{31}\text{.}\text{559}+\text{.}\text{1354}({T}_{0}-\text{273}\text{.}5)})}\) , for example: , and a degree of hygrometry \(\text{HR}\) , we deduce the vapor pressure thanks to \({p}_{\text{vp}}({T}_{0})\text{=}\text{HR}{p}_{\text{vp}}^{\text{sat}}({T}_{0})\) .

3.4. Nodal unknowns, initial values, and reference values

Here we are addressing a point that is more related to programming choices than to real formulation aspects. However, we are exposing it because it has important practical consequences. The main unknowns, which are also the values of the degrees of freedom, are noted:

\({\left\{u\right\}}^{\text{ddl}}\text{=}\left\{\begin{array}{c}{u}_{x}\\ {u}_{y}\\ {u}_{z}\\ \text{PRE}{1}^{\text{ddl}}\\ \text{PRE}{2}^{\text{ddl}}\\ {T}^{\text{ddl}}\end{array}\right\}\)

Depending on the model, they can have different meanings:

	LIQU_SATU	LIQU_VAPE	LIQU_GAZ_ATM		GAZ	LIQU_VAPE_GAZ
PRE1	\({p}_{w}\)	\({p}_{w}\)		\({p}_{c}=-{p}_{w}\)	\({p}_{\text{gz}}\)	\({p}_{c}={p}_{\text{gz}}-{p}_{w}\)
PRE2					\({p}_{\text{gz}}\)

	LIQU_GAZ	LIQU_AD_GAZ_VAPE	LIQU_AD_GAZ
PRE1	\({p}_{c}={p}_{\text{gz}}-{p}_{w}\)	\({p}_{c}={p}_{\text{gz}}-{p}_{w}-{p}_{\text{ad}}\)	\({p}_{c}={p}_{\text{gz}}-{p}_{w}-{p}_{\text{ad}}\)
PRE2	\({p}_{\text{gz}}\)	\({p}_{\text{gz}}\)	\({p}_{\text{gz}}\)

	LIQU_GAZ	LIQU_AD_GAZ_VAPE	LIQU_AD_GAZ
PRE1	\({p}_{c}={p}_{\text{gz}}-{p}_{w}\)	\({p}_{c}={p}_{\text{gz}}-{p}_{w}-{p}_{\text{ad}}\)	\({p}_{c}={p}_{\text{gz}}-{p}_{w}-{p}_{\text{ad}}\)
PRE2	\({p}_{\text{gz}}\)	\({p}_{\text{gz}}\)	\({p}_{\text{gz}}\)

We will then define the real pressures and the real temperature by:

\(p\text{=}{p}^{\text{ddl}}\text{+}{p}^{\text{init}}\) for pressures PRE1 and PRE2 and \(T\text{=}{T}^{\text{ddl}}\text{+}{T}^{\text{init}}\) for temperatures, where \({p}^{\text{init}}\) and \({T}^{\text{init}}\) are defined under the THM_INIT keyword in the DEFI_MATERIAU command.

The values written by IMPR_RESU are the values of \({\left\{u\right\}}^{\text{ddl}}\). Boundary conditions are set for \({\left\{u\right\}}^{\text{ddl}}\). The DEPL keyword in the ETAT_INIT factor keyword in the STAT_NON_LINE command defines the initial values of \({\left\{u\right\}}^{\text{ddl}}\). The initial values of the enthalpies, which belong to the generalized constraints, are defined using the keyword SIGM of the keyword factor ETAT_INIT of the command STAT_NON_LINE. The real pressures and the real temperature are used in the laws of behavior, in particular laws of the type \({S}_{\mathrm{lq}}\text{=}f({p}_{c})\) or \(\frac{\text{dp}}{p}\text{=}\frac{d\rho }{\rho }+\frac{\text{dT}}{T}\). The initial values of the densities of steam and dry air are defined from the initial values of the gas and steam pressures (values read under the keyword THM_INIT of the command DEFI_MATERIAU). Note that, for the movements, the decomposition \(u\text{=}{u}^{\text{ddl}}\text{+}{u}^{\text{init}}\) is not done: the keyword THM_INIT of the command DEFI_MATERIAU does not therefore make it possible to define initial movements. The only way to initialize the movements is therefore to give them an initial value using the keyword factor ETAT_INIT from the STAT_NON_LINE command.

3.5. Effective constraints and total constraints. Stress limit conditionsThe partition of stresses into total and effective stresses is written as:

\(\sigma \text{=}\sigma \text{'}\text{+}{\sigma }_{p}\)

\(\sigma\) is the total constraint, i.e. the one that verifies: \(\text{Div}(\sigma )\text{+}r{F}^{m}\text{=}0\)

\({\sigma }^{\text{'}}\) is the effective constraint. For the laws of effective constraints, it checks: \(d{\sigma }^{\text{'}}\text{=}f(d\varepsilon \text{-}{\alpha }_{0}\text{dT},\alpha )\), where \(\varepsilon \text{=}1/2(\nabla u\text{+}{}^{T}\nabla u)\) and \(\alpha\) represent the internal variables. The tensor \({\sigma }_{p}\) is calculated according to hydraulic pressures. The writing adopted is incremental and, if we want the value of \({\sigma }_{p}\) to be consistent with the value \({p}^{\text{init}}\) defined under the keyword THM_INIT, we must initialize the components of \({\sigma }_{p}\) with the keyword SIGM of the keyword factor ETAT_INIT of the command STAT_NON_LINE.

In the result files, we find the effective constraints \({\sigma }^{\text{'}}\) under the component names SIXX… and \({\sigma }_{p}\) under the name SIPXX… The stress boundary conditions are written as total stresses.

3.6. Some numerical values

Here we give some reasonable values for some coefficients. These values are not programmed in*Code_Aster*, they are provided here for information purposes only:

For ideal gases, the following values are retained:

\(R=8\text{.}\text{3144}J\text{.}{K}^{-1}\)

\({M}_{\text{vp}}^{\text{ol}}=\text{18}\text{.}{\text{10}}^{\text{-}3}\text{kg}\text{.}{\text{mole}}^{-1}\)

\({M}_{\text{as}}^{\text{ol}}=\text{28}\text{.}\text{96}{\text{10}}^{\text{-}3}\text{kg}\text{.}{\text{mole}}^{-1}\)

For CO2, the value of the Henry constant at \(20°C\) is:

\({K}_{H}\text{=}\text{3162}\text{Pa}\text{.}{m}^{3}{\text{mole}}^{-1}\)

For liquid water, we have:

\({\rho }_{w}\text{=}\text{1000}\text{kg}/{m}^{3}\)

\({K}_{w}\text{=}\text{2000}\text{MPa}\)

The coefficient of thermal expansion of water is correctly approximated by the formula:

\({\alpha }_{w}=9\text{.}\text{52}{\text{10}}^{\text{-}5}\text{ln}(T\text{-}\text{273})-2\text{.}\text{19}{\text{10}}^{\text{-}4}\) (\({K}^{\text{-}1}\))

The calorific capacities have the following values:

\(\begin{array}{}{C}_{\sigma }^{s}=\text{800}{\text{J.kg}}^{-1}{K}^{-1}\\ {C}_{\text{lq}}^{p}=\text{4180}{\text{J.kg}}^{-1}{K}^{-1}\\ {C}_{\text{vp}}^{p}=\text{1870}{\text{J.kg}}^{-1}{K}^{-1}\\ {C}_{\text{as}}^{p}=\text{1000}{\text{J.kg}}^{-1}{K}^{-1}\end{array}\)

A law of evolution of the latent heat of change of liquid vapor phase is also given:

\(L(T)=\text{2500800}-\text{2443}(T-\text{273}\text{.}\text{15})J/\text{kg}\)