Coupled Hydromechanical Model of Two-Phase Fluid Flow in Deformable Porous Media

A model of solid-water-air coupling in triphasic mixtures is compared with solid-water coupling in biphasic mixtures with an application to partially saturated porous media. Based on thermodynamics, the mathematical framework governing the behavior of a partially saturated soil is derived using balance equations, and the numerical implementation and drainage tests of a soil column are carried out to validate the obtained formulations.The role of the air phase in the hydro-mechanical behavior of triphasic mixtures can be analyzed from the interactions among multiple phases for the constitutive behavior of a solid skeleton, and the triphasicmixturemodel can be applied in geotechnical engineering problems, such as CO 2 sequestration and air storage in aquifers.


Introduction
For the features of the hydromechanical behavior of multiphase porous media regarding solid-water-gas interactions, the soil model has two kinds of constituent volume fractions in geotechnical engineering problems.One type is when the air phase is considered active in the voids of a threephase soil as a gas phase, and the other is to assume the gas phase as a vacuum in a soil mixture.In soil mechanics, "saturated" means the mixture of soil particles and water, but other fields, such as earth science and multi-phase media, define the mixture concept differently, as shown in Figure 1.
Depending on the role of the air phase, some soil has occluded air bubbles in triphasic mixture, and another one has continuous air bubbles in biphasic mixture.The air phase plays a role in the air-water interface of partially saturated soil or contractile skin of hydromechanical behavior of deformable soil as a four-phase system [1].When performing stress analysis on an element, Fredlund and Rahardjo [1] assumed that partially saturated soil can be visualized as a mixture with soil particles and a contractile skin that approaches equilibrium under applied stress gradients and air and water phases that flow under applied stress gradients.The contractile skin acts like a thin rubber membrane, pulling the particles together when the pore water pressure gradually becomes negative during shrinkage-type experiments involving the drying of a small soil specimen as it is exposed to atmosphere.
Recently, in order to observe coupled solid-water-air phenomenon in detail, many researchers have formulated conservation equations for a three-phase system consisting of a solid and two immiscible fluids: liquid and gas [2][3][4][5][6][7].These researchers also derived expressions for the effective stress tensor in multi-phase porous media exhibiting two porosity scales, micro-and macroporosity, during the course of loading.Figure 2 illustrates the concept of partially saturated porous media with double porosity.The three-phase mixture is composed of porous media consisting of solid and fluids in the micropores and consisting of only fluid constituents, liquid, and gas, in the macropores, as shown in Figure 2.
The developments presented in the literature regarding mixtures of solid, liquid, and gas have been proposed as a constitutive framework of a coupled model to simulate water and air flow in deformable soil.Finite element analysis of partially saturated soil is generally treated as a biphasic (i.e., two field phases) mixture state in geotechnical engineering.Equations regarding the volume and mass of a mixture are defined as the mathematical relations [8,9].The volume of a mixture is V = V  + V  + V  , and the corresponding total mass is  =   +   +   .Similarly, for the  phase,   =   V  (nearly homogeneous), where  = , ,  and   is the true mass density of the  phase.The volume fraction occupied by the  phase is given by   = V  /V, and thus, for water, air, and solid where the porosity  = (V  + V  )/V =   +   .If a material is homogeneous,   = V  /V, whereas if it is heterogeneous, the volume fraction at a material point   = dV  /dV for a differential volume of the mixture.The partial mass density of the  phase is given by   =     , and thus where  = /V is the total mass density of the mixture.As a general notation, phase designations in the superscript form (e.g.,   ) pertain to average or partial quantities, and those in the superscript form with  (e.g.,   ) to intrinsic or real quantities.Based on the current configuration of the mixture (for small strains, theoretically not different from the reference or current configurations), the mass balance equations describe the motions of the water and air phases relative to the motion of the solid phase.In addition, mixture theory assumes that the three phases are smeared together at a spatial position in the current configuration, thus limiting the theory to a continuum representation of a partially saturated soil; it is large enough length scale that the soil behaves as a continuum.
In order to solve for three unknowns (u,   ,   ) using three equations given in (3)-( 5), based on the formulation of Coussy [8], Borja [10], and de Boer [9], we can write the balance of linear momentum for a triphasic mixture and the Figure 2: Schematic description of a mixture with double porosity [7].
balance of the mass of the solid, water phase, and air phase as [11] div  + g = 0, where the total stress is written in terms of the effective stress   (positive in tension) as  =   −   1 + (  −   )1 [10,12] and the degree of saturation, , is defined in the classical form of [13] where   is the effective degree of saturation,   is the residual degree of saturation,   and   are the pore water pressure and the pore air pressure (positive in compression), respectively, and , , and  are the soil water characteristic curve parameters (SWCC).g is the gravity acceleration vector, ṽ =   (v  −v  ) is the Darcy seepage velocity of water where v  is the true velocity of water, and v  is the velocity of the solid skeleton.The constitutive equations are the linear isotropic elasticity for the effective stress and the generalized Darcy's law for the Darcy seepage velocity of water, written, respectively, as, where  is two constituents (water and air), the fourth-order elastic modulus tensor is c  = 1 ⊗ 1 + 2I, the Lame parameters are  and ,  = sym∇u is the symmetric small strain tensor, and   and   are the partially saturated hydraulic conductivity function written as (1) water flow in a partially saturated porous medium [13] (2) air flow in a porous medium [8] where the material property  is the intrinsic permeability of the soil skeleton and the function of the porosity ,   is the relative permeabilities related to, respectively, water and air, and   is the dynamic water and air viscosity.The density of air   is approximately 1.2 kg/m 3 at sea level and at 20 ∘ C and the air bulk modulus   is 10 5 Pa at a constant temperature.The relative permeabilities   and   as a function of  are given in Figure 3.
It can be shown that () =  2 (), where  2 is a parameter of dimension area (m 2 ) and () by the Kozeny-Carman relation (pore space formed by regular packing of spheres) [8] for representing the porosity dependence of hydraulic conductivity.The porosity  is a function of the volumetric strain of the solid skeleton.

Weak Form and Coupled Finite Element Formulation
It is assumed that the whole domain of the body  is partially saturated.Applying the method of weighted residuals [11,15], the coupled weak form for a triphasic mixture is written as where w is the weighting function for the displacement u, t is the applied traction,  and  are the weighting function for the pore water pressure (  ) and pore air pressure (  ), respectively, and   is the positive inward water seepage on the boundary Γ  .Assuming a mixed finite element formulation as indicated by the quadrilateral elements in the example mesh in Figure 4, the discretized displacement u ℎ is interpolated biquadratically, and the pore water pressure ( ℎ  ) and pore air pressure ( ℎ  ) bilinearly [15].Introducing the shape functions and expressing in matrix form, the coupled nonlinear finite element form is written as where A  el =1 is the element assembly operator for element  over the number of elements  el , c  is the element nodal weighting function values for  ℎ ,   and   are the element nodal weighting function values for  ℎ and  ℎ , respectively (both of which are arbitrary except where they are zero at the boundaries with essential boundary conditions), d  is the element nodal displacement vector, and   and   are the element nodal pore water and air pressure vector.
After applying the boundary conditions and assembling the finite element equations, the matrix form is obtained  as the monolithically coupled nonlinear first-order ordinary differential equation to solve for D as where, where C is the combination of the damping matrix and the stiffness matrix of the d.o.f.vector time derivative and K is the stiffness matrix.Then, the location matrix (LM) can be used to assemble the individual 26 × 26 and 26 × 3 contributions to the global "damping" matrix C, the stiffness matrix K, and the forcing vector F, and the matrix form uses generalized trapezoidal integration to solve transient equations [15].For consolidation analysis, the generalized trapezoidal rule is used to integrate transient problems by FE coupled balance of mass and linear momentum equations at time  +1 and derived from difference formulas for D +1 and V +1 ,  where velocity V +1 is Ḋ( +1 ) and  is the time integration parameter: The form of ( 14) allows us to consider a semi-implicit integration scheme of a linear form which is written as

Numerical Results for Two-Fluid Flow in a Deformable Porous Medium
As previously mentioned, since there is no exact solution for the problem of water and air flow in deformable partially Mathematical Problems in Engineering saturated soils, numerical modeling of the experimental results of the drainage of a sand column is performed.For validation and application of a coupling model of solid, water, and air in partially saturated soils based on the water drainage experiment of a sand column described by Liakopoulos [16], the numerical solutions given by Schrefler and Scotta [17] and Gawin et al. [18] are compared to various results obtained from the coupled model.The mesh of this example, which is composed of a column of 20 nine-node isoparametric Lagrangian elements of equal size, was employed for all numerical simulations.Numerical integration was semiimplicit, and the triphasic model associated with linear elasticity used mesh of 2D plain strain of nine integration gauss points.Figure 5 shows the initial and boundary conditions (left) for numerical analysis and the procedure of the experimental test (right).A soil column test with a column 0.5 m in height is carried out to investigate the approximate value of pore air pressure for the validation of numerical results, even though the simulation used a soil column 1 m in height for the numerical analysis.
The physical experiment consisted of a soil column 1 m in high and a constant flow through the soil column corresponding to a water pressure gradient initially equal to zero initially.At the starting time steps, the water inflow is cut at the top of the soil column and the water is flowed out at the bottom.Air pressure is equal to atmospheric pressure at both the top and bottom of the column, with zero vertical load at the top and no deformation at the bottom and on the lateral walls of the column.The gravity-governed changes in the constituent volume fractions only depend on the soil and water parameters.In the numerical test, with the same properties and boundary conditions implemented by Schrefler and Scotta [17], the coupled model also uses the relationship of Brooks and Corey [19] for the relative permeability of gas pressure and the experimental function of Schrefler and Scotta [17] for the hydraulic properties of the soil, as shown in (16).The material properties used for the numerical test are summarized in Table 1.
One has where   is the dynamic viscosity and   is the relative permeability of the  phase, which depends on the relative saturation   through the experimental relationship   =   (  ),  is the intrinsic permeability, and the respective degrees of saturation   and   sum to one:   +   = 1.Even if the data of the mechanical behavior and the parameters of the Del Monte sand used by Liakopoulos [16] were missing and unpublished, his solutions have been obtained numerically by trial and error techniques.Thus,  is 0.1 and the residual saturation   is 0.06689 for sand [17,18].The triphasic mixture analysis of Schrefler and Scotta [17], which was results of the numerical solutions based on the Liakopoulos [16], are compared to those from coupled code, as shown in Figures 6-9.As no measurement of pore air pressure was made by Liakopoulos [16], the numerical results are plotted and also compared to other results [17,18] in Figure 7.The evolution of air pressure is more sensitive to the analysis method than that of water pressure.The comparison of pore water pressure in Figure 6(a) is similar to that of Schrefler and Scotta [17], but the results (Figure 6(b)) of the coupled code showed suction increases slower than those found by Schrefler and Scotta [17] since the air pressure response from the methods applied is sensitive, as shown in the suction evolution in Figure 7.
Comparing with the two previous numerical results, the air pressure profiles from the coupled model fit closer to that of Gawin et al. [18] than that of Schrefler and Scotta [17].These differences are produced by choosing different sets of governing equations and numerical algorithms.In particular, the averaged density of the mixture  = (1 − )  +   + (1 − )(1 − )  and the bulk modulus of the solid grains (10 6 MPa) and the water (2 × 10 3 MPa) used by Schrefler and Scotta [17] are different from those used by Gawin et al. [18] and the coupled code.Gawin et al. [18] and the coupled code both derived the mass balance equation assuming the bulk moduli (  and   ) are infinite due to the large values.The averaged density of the mixture is  = (1 − )  +   + (1 − )  .
For sand and weathered soil types, experimental tests are performed to investigate the pore air pressure at 5 cm, 10 cm, and 15 cm place from top surface of soil sample.Figure 8 shows that experimental results are similar to those of the coupled code although the air pore pressures are just measured by three sensors at the top portion of the soil column.
As shown in Figure 9, the vertical displacements at the top surface of the soil sample show little difference with time, but the final vertical displacements coincide with those of Schrefler and Scotta [17] and Gawin et al. [18] under identical initial conditions.Because the coupled code has the hydraulic conductivity,   , which is the function of porosity, the vertical displacement of the coupled code deforms little faster than those of other numerical solutions at early time step.

Conclusions
We have implemented a numerical integration algorithm (semi-implicit solution) for solid-water-air coupling finite element formulation using balance equations.Based on Liakopoulos [16]'s experimental results, Gawin et al. [18] and Schrefler and Scotta [17] presented numerical simulations for the behavior and the diffusion of air pressure in a drainage test of a soil column.In this study, the developed coupled finite element model for a deformable partially saturated soil based on linear isotropic elasticity describes the poromechanical behavior of a soil column by linking solid displacement, pore water pressure and air pressure simultaneously.The results of the coupled model approach the simulation of drainage test because it uses partially saturated permeability which is the function of porosity.The numerical results of the coupled model are more similar to those of Gawin et al. [18] rather than to those of Schrefler and Scotta [17] regarding the diffusion and dissipation of air pressure, matric suction, and vertical displacement.The coupled model was validated through comparisons with the literature and through laboratory tests of the drainage of a soil column, and the results of two fluids flow obtained by semi-implicit linear solution also demonstrate the stability of the solution by comparing nonlinear models of Gawin et al. [18] and Schrefler and Scotta [17].

Figure 5 :
Figure 5: Diagram of the numerical analysis and the experimental test.

Figure 6 :
Figure 6: Comparison of the numerical results and experimental data with pore water pressure and matric suction.

Figure 7 :Figure 8 :
Figure 7: Comparison of numerical results of air pressure,   .

Figure 9 :
Figure 9: Displacement at the top surface of a drainage test in a triphasic mixture.

Table 1 :
Soil parameters for triphasic mixture implementation.