====== HYPOFE2 ======
===== Description =====

Multiscale law for water-air seepage, pollutant diffusion and advection. Inspired from WAVAT and ADVEC.

Can be parallelized with ELEMB (macroscale) or at the perturbation loop (microscale).

Takes into account the hysteresis in the water retention law when used with FKRSAT. Can also be used with osmotic suction (under development).
==== The model ====
This law is only  used for water seepage - air seepage- pollutant diffusion and advection (coupled with water or gas flows) for non linear analysis in 2D porous media.

=== Mass conservation of water (liquid and vapour) ===
\[
\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})}_{\text{Vapeur}} = 0
\]

=== Liquid and vapour flows ===
Starting from Darcy's law, the liquid water velocity is:
\[
\vec{q_l} = - \frac{k_w}{\mu_w}\left[ \vec{grad}(p_w) + g \; \rho_w \; \vec{grad}(y) \right]\ \text{where}\ k_w = K_w\; \frac{\mu_w}{\rho_w\; g}\left[ m^2\right]
\]

The water vapour only flows in unsaturated pores and depends on the tortuosity of the path:
\[
\vec{i}_v = - n \; S_{r,g} \; \tau D\;  \rho_s \; \vec{grad} \omega_v
\]
Where $\omega_v = \rho_v/\rho_g$ is the dry air mass content in the gaseous mix.


=== Liquid State Equations ===

  - Density: $\rho_w$: \[\rho_w (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)$ with $k_{w,eff} = k_f k_{r,w}$
  - Saturation degree $S_w$: \\ Depending on suction $s = p_g - p_w : S_w = f(s)$

=== Saturation degree equation (with FKRSAT) ===
ISR = 53 Van Genuchten model (ISR=5) with hysteresis implemented.

The main water retention curves (d=drying, w=wetting) are, according to the Van Genuchten model:
\[S_{ed} = S_{res} + (S_{max}-S_{res}) \left[1 + \left(\frac{s}{a_d}\right)^{n_d}\right]^{-m_d}\] 
\[S_{ew} = S_{res} + (S_{max}-S_{res}) \left[1 + \left(\frac{s}{a_w}\right)^{n_w}\right]^{-m_w}\]

The hysteresis is then defined by:
\[\frac{\partial S_{es}}{\partial s} (\text{wetting}) = \left(\frac{s_w}{s}\right)^b\left(\frac{\partial S_{ew}}{\partial s}\right) \text{ with } s_w = a_w \left(S_e^{-1/m_w}-1\right)^{1/n_w}\] 
\[\frac{\partial S_{es}}{\partial s} (\text{drying}) = \left(\frac{s_d}{s}\right)^{-b}\left(\frac{\partial S_{ed}}{\partial s}\right) \text{ with } s_d = a_d \left(S_e^{-1/m_d}-1\right)^{1/n_d}\]

And therefore:
\[S_e^{t+1} = S_e^t + \left(\frac{\partial S_{es}}{\partial s}\right)\times ds\]

The ISR=53 parameters are: CSRW1=$a_d$, CSRW2=$n_d$, CSRW3=$a_w$, CSRW4=$n_w$ and CSRW5=$b$

=== Osmotic suction model ===

TO BE COMPLETED. 

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

=== Dry air and dissolved gas flows ===
Starting from Darcy's law, the gas velocity is:
\[
\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]
\]

The diffusion velocity of dry air is proportional to a density gradient. Using the diffusion theory adapted to porous medium, one writes:
\[
\vec{i}_a = - n S_{r,g} \tau D \rho_g \vec{grad} \omega_a = -\vec{I}_v
\]
Where $\omega_a = \rho_a/\rho_g$ is the dry air mass content inside the gas mix.

=== Gas State Equation ===

  - Density $\rho_a$ :\\ //Hypothesis// : The air is supposed to be a perfect gas. \[\rho_a (p_a) = \rho_{a,0}\frac{p_a}{p_{a,0}} \]
  - Intrinsic Permeability $k_g$: \\ Depending on the saturation degree $S_g$ : $k_{r,g} = f(S_g)$ with $k_{g,effectif} = k_{g, intrinsic}k_{a,w}$
  - Saturation degree $S_g$: \\ Depending on suction $s = p_g - p_w : S_g = f(s) = 1 - S_w$

=== Balance Equation of Pollutant ===
\[\frac{\partial}{\partial x_i} (v_i^p) = 0\]

=== Pollutant flows ===
\[ v_i^p = v_i^{advection} + v_i^{diffusion+dispersion} = C_M v_i^{w/g} - D \frac{\partial C_m}{\partial x_i} \]\\
With $C_M$ and $C_m$ [-] the concentration in pollutant at the macroscale and subscale, respectively. $v_i^{w/g}$ is the water or gas velocity obtained from Darcy's law and $D$ [m$^2$/s] is the diffusion and dispersion coefficient.
==== Files ====
Prepro: LHYPOFE2.F \\
Lagamine: HYPOFE2.F \\
===== Availability =====
|Plane stress state| NO |
|Plane strain state| YES |
|Axisymmetric state| NO|
|3D state| YES |
|Generalized plane state| NO |
===== Input file =====
==== Parameters defining the type of constitutive law ====
^ Line 1 (2I5, 60A1)^^
|IL|Law number|
|ITYPE| 629 |
|COMMENT| Any comment (up to 60 characters) that will be reproduced on the output listing|
==== Integer parameters ====
^ Line 1 (3I10,2G10.0) ^^
|NLAWFEM2|Number of constitutive laws at the subscale|
|KFLU|Number of DOF at the microscale: 1 = $P_w$, 2 = $P_w+C$, 3 = $P_w+P_g$, 4 = $P_w+C+P_g$ with $C$ the concentration in pollutant|
|MITER|Maximum number of iterations at the subscale|
|CNORM|Norm for the solver of the subscale|
|FACONV|Units of conversion of the RVE (it has a size of 1*FACONV[-])|



==== Real parameters ====
^ Line 1 (3E10.2,2G10.0) ^^
|VISCW0|Liquid dynamic viscosity $(=\mu_{w,0})\ \left[ Pa.s \right]$|
|RHOW0|Liquid density $(=\rho_{w,0})\ \left[ kg.m^{-3}\right]$|
|UXHIW|Liquid compressibility coefficient $(=1/ \chi_{w})\ \left[ Pa^{-1}\right]$|
|PW0|Initial water pressure $\left[ Pa\right]$|
|T0|Initial temperature $\left[ K\right]$|
^ Line 2 (1G10.0) ^^
|CPINI|Initial pollutant concentration $\left[ -\right]$|
^ Line 3 (3E10.2,2G10.0) ^^
|VISCA0|Gas dynamic viscosity $(=\mu_{a,0})\ \left[Pa.s \right]$|
|RHOA0|Gaz density $(=\rho_{a,0})\ \left[kg.m^{-3}\right]$|
|PMGAS|Gas molar mass $[g/mol]$|
|PG0|Initial gas pressure $\left[ Pa\right]$|
|PHENRY|Henry coefficient|
^ Line 4 (4I10) ^^
|IVAP|= 1 for vapour, = 0 if liquid water only|
|IGAS|= 0 for air, =1 for $H_2$, =2 for $N_2$, = 3 for $Ar$, = 4 for $He$, = 5 for $CO_2$, = 6 for $CH_4$|
|IOSMOTIC|= 0 to neglect osmotic suction, = 1 for osmotic suction with Van't Hoff model, = 2 for osmotic suction with Kelvin (water activity) and Pitzer model|
|IDIFF|= 0 for the pollutant to diffuse through water, = 1 through gas|
^ Line 5 (4I10) ^^
|ISR|Retention curve (=53 for Van Genuchten with hysteresis)|
|IKW|Water relative permeability curve (=7 for Van Genuchten)|
|IKA|Gas relative permeability curve (=6 for Van Genuchten)|
|N_SUBINCR|Number of additional multiplicator for the number of subincrement in the hysteresis model|
^ Line 6 (3G10.0)^^
|CKW1|First parameter of IKW|
|CKW2|Second paremeter of IKW|
|CKW3|Third parameter of IKW|
^ Line 7 (2G10.0)^^
|CKA1|First parameter of IKA|
|CKA2|Second parameter of IKA|
^ Line 8 (5G10.0)^^
|CSR1|First parameter of ISR|
|CSR2|Second parameter of ISR|
|CSR3|Third parameter of ISR|
|CSR4|Fourth parameter of ISR|
|CSR5|Fifth parameter of ISR|
^ Line 9 (5G10.0)^^
|SRES|Residual saturation degree $(=S_{res})$|
|SRFIELD|Field saturation degree $(=S_{r, field})$|
|AIREV|Air entry pressure $\left[Pa\right]$|
|AKRMIN|Minimum value of relative permeabikity|
|SRINI|Initial saturation degree|

==== Subscale parameters ====
To be repeated as many time as NLAWFEM2.
^ Line 1 (2I5) ^^
|ILAW2|Number of the subscale constitutive law (=1:NLAWFEM2)|
|ITYPE2|Type of subscale law (=1 for Hydraulic pollutant microscale law)|
^ Line 2 (4G10.0) ^^
|POROS|Material porosity ($=n$)|
|PERMINT|Material intrinsic permeability ($=k_{int}$) $[m^2]$|
|DIFFC|Material diffusion coefficient of the pollutant ($D_{app}$) $[m^2/s]$|
|TORTU|Material tortuosity ($=\tau$)|

===== Stresses =====
==== Number of stresses ====
28
==== Meaning ====
__In 2D state :__
|SIG(1)|$\sigma_x$ (unused)|
|SIG(2)|$\sigma_y$ (unused)|
|SIG(3)|$\sigma_{xy}$ (unused)|
|SIG(4)|$\sigma_z$ (unused)|
|SIG(5)|Homogenised liquid flow along $x$ $(=f_{wx})$|
|SIG(6)|Homogenised liquid flow along $y$ $(=f_{wy})$|
|SIG(7)|Homogenised liquid flow stored $(=f_{we})$|
|SIG(8)|Homogenised mean flow of the pollutant along $x$ $(=(f_{px,a}+f_{px,b})/2)$|
|SIG(9)|Homogenised mean flow of the pollutant along $y$ $(=(f_{py,a}+f_{py,b})/2)$|
|SIG(10)|Homogenised pollutant flow stored (takes advection into account) $(=f_{pe})$|
|SIG(11)|Homogenised diffusive flow of the pollutant along $x$ for the current step $(=f_{px,b})$|
|SIG(12)|Homogenised diffusive flow of the pollutant along $y$ for the current step $(=f_{py,b})$|
|SIG(13)|Homogenised gas flow along $x$ $(=f_{gx})$|
|SIG(14)|Homogenised gas flow along $y$ $(=f_{gy})$|
|SIG(15)|Homogenised gas flow stored $(=f_{ge})$|
|SIG(16)|Advective flow of dissolved gas along $x$ (unused)|
|SIG(17)|Advective flow of dissolved gas along $y$ (unused)|
|SIG(18)|Vapour flow along $x$ $(=f_{vx})$|
|SIG(19)|Vapour flow along $y$ $(=f_{vy})$|
|SIG(20)|Vapour flow stored $(=f_{ve})$|
|SIG(21)|Unused|
|SIG(22)|Unused|
|SIG(23)|Unused|
|SIG(24)|Unused|
|SIG(25)|Unused|
|SIG(26)|Unused|
|SIG(27)|Unused|
|SIG(28)|Unused|

===== State variables =====
==== Number of state variables ====
11 + 5*(Number of Subscale Nodes)\\
/!\ The state variables vector also contains the following information for each subscale node: $X$, $Y$, $P_w$, $C$, $P_g$
==== List of state variables ====
|Q(1)|Liquid water mass at the RVE|
|Q(2)|Pollutant mass at the RVE|
|Q(3)|Gaseous air mass at the RVE|
|Q(4)|Homogenised macroscale porosity|
|Q(5)|Water saturation degree|
|Q(6)|Homogenised water relative permeability|
|Q(7)|Homogenised gas relative permeability|
|Q(8)|Homogenised macroscale tortuosity|
|Q(9)|Vapour mass at the RVE|
|Q(10)|Homogenised total suction $(= p_g - p_w + osmotic)$|
|Q(11)|Homogenised osmotic suction $(= osmotic)$|
|Q(12 + (i-1)*5)|$X_i$|
|Q(12 + (i-1)*5 +1)|$Y_i$|
|Q(12 + (i-1)*5 +2)|$P_{w,i}$|
|Q(12 + (i-1)*5 +3)|$C_i$|
|Q(12 + (i-1)*5 +4)|$P_{g,i}$|