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

Authors: You-Seong Kim; Jaehong Kim
Journal: Mathematical Problems in Engineering
 (2013)

---

## Abstract

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 triphasic mixture model can be
applied in geotechnical engineering problems, such asCO2 sequestration and
air storage in aquifers.


---

## Body

## 1. 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 three-phase 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 Figure1.Figure 1
Schematic description of three-phase mixtures [14].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–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.Figure 2
Schematic description of a mixture with double porosity [7].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.


## 2. Balance Equations and Constitutive Equations

Partially saturated soil which has governing equations based on mixture theory consists of solid (s), water (w), and air (a). Equations regarding the volume and mass of a mixture are defined as the mathematical relations [8, 9]. The volume of a mixture is v=vs+vw+va, and the corresponding total mass is m=ms+mw+ma. Similarly, for the α phase, mα=ραRvα (nearly homogeneous), where α=s,w,a and ραR is the true mass density of the α phase. The volume fraction occupied by the α phase is given by nα=vα/v, and thus, for water, air, and solid
(1)ns+n=ns+nw+na=1,
where the porosity n=(vw+va)/v=nw+na. If a material is homogeneous, nα=vα/v, whereas if it is heterogeneous, the volume fraction at a material point nα=dvα/dv for a differential volume of the mixture. The partial mass density of the α phase is given by ρα=nαραR, and thus
(2)ρs+ρw+ρa=ρ,
where ρ=m/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 R (e.g., ραR) 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,pw,pa) 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 balance of the mass of the solid, water phase, and air phase as [11]
(3)divσ+ρg=0,(4)n∂S∂s(p˙a-p˙w)+Sdivv=-divv~w,(5)[-n∂S∂s+n(1-S)Ka]p˙a+n∂S∂sp˙w+(1-S)divv0000=-1ρaRdiv(ρaRv~a),
where the total stress is written in terms of the effective stress σ′ (positive in tension) as σ=σ′-pa1+S(pa-pw)1 [10, 12] and the degree of saturation, S, is defined in the classical form of [13]
(6)S=Sr+Se(1-Sr),Se=[1+(αs)n]-m,
where Se is the effective degree of saturation, Sr is the residual degree of saturation, pw and pa are the pore water pressure and the pore air pressure (positive in compression), respectively, and α, n, and m are the soil water characteristic curve parameters (SWCC). g is the gravity acceleration vector, v~w=nw(vw-vs) is the Darcy seepage velocity of water where vw is the true velocity of water, and vs 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,
(7)σ′=ce:ϵ,v~α=kα(S,n)[-∇pα+ραRg],
where α is two constituents (water and air), the fourth-order elastic modulus tensor is ce=λ1⊗1+2μI, the Lame parameters are λ and μ, ϵ=sym∇u is the symmetric small strain tensor, and kw and ka are the partially saturated hydraulic conductivity function written as(1)
water flow in a partially saturated porous medium [13]
(8)v~w=kw(n,S)[-∂pw∂x+ρwRg],kw(n,S)=ϰ(n)ηwkrw(S),krw(S)=S(1-(1-S1/m)m)2,(2)
air flow in a porous medium [8]
(9)v~a=ka(n,S)[-∂pa∂x+ρaRg],ka(n,S)=ϰ(n)ηakra(S),kra(S)=1-S(1-S1/m)2m,ρaR=exp[paKa],δ(n)=n31-n2,where the material property ϰ is the intrinsic permeability of the soil skeleton and the function of the porosity n, krα is the relative permeabilities related to, respectively, water and air, and ηα is the dynamic water and air viscosity. The density of air ρaR is approximately 1.2kg/m3 at sea level and at 20°C and the air bulk modulus Ka is 105 Pa at a constant temperature. The relative permeabilities krw and kra as a function of S are given in Figure 3.Relative permeabilitieskrw (a) and kra (b) plotted against S when varying m from 0.1 to 1.
(a)
(b)It can be shown thatϰ(n)=l2δ(n), where l2 is a parameter of dimension area (m2) and δ(n) by the Kozeny-Carman relation (pore space formed by regular packing of spheres) [8] for representing the porosity dependence of hydraulic conductivity. The porosity n is a function of the volumetric strain of the solid skeleton.


## 3. Weak Form and Coupled Finite Element Formulation

It is assumed that the whole domain of the bodyB is partially saturated. Applying the method of weighted residuals [11, 15], the coupled weak form for a triphasic mixture is written as
(10)∫B‍[∇w:(σ′-Spw1-(1-S)pa1)]dv0000=∫B‍ρw·gdv+∫Γt‍w·tda,∫Bw‍ηn∂S∂s(p˙a-p˙w)dv+∫BwηSdivvdv-∫Bw‍∇η·v~wdv00=∫Γsw‍ηSwds,-∫Ba‍φn∂S∂sp˙adv+∫Ba‍φn(1-S)Kap˙adv+∫Ba‍φn∂S∂sp˙wdv0000+∫Ba‍φ(1-S)divvdv+∫BaφKa∇pa·v~adv-∫Ba‍∇φ·v~adv0000=∫Γsa‍φSads,
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 (pw) and pore air pressure (pa), respectively, and Sw is the positive inward water seepage on the boundary Γs. Assuming a mixed finite element formulation as indicated by the quadrilateral elements in the example mesh in Figure 4, the discretized displacement uh is interpolated biquadratically, and the pore water pressure (pwh) and pore air pressure (pah) bilinearly [15].Figure 4
Finite element mesh (20 elements) for numerical simulation.Introducing the shape functions and expressing in matrix form, the coupled nonlinear finite element form is written as(11)Ae=1nel(ce)T·[ke,dd·de-ke,dθ(θe,ζe)·θe000000000000-ke,dζ(θe,ζe)·ζe000000000000=ffe,d(de,θe,ζe)+fte,d],Ae=1nel(αe)T·[ke,θζ(de,θe,ζe)·ζ˙e00000000000-ke,θζ(de,θe,ζe)·θ˙e00000000000+ke,θd(θe,ζe)·d˙e00000000000+ke,θθ(de,θe,ζe)·θe00000000000=fe,θ(de,θe,ζe)+fse,θ],Ae=1nel(βe)T·[ke,ζθ(de,θe,ζe)·ζ˙e00000000000+ke,θζ(de,θe,ζe)·θ˙e00000000000+ke,ζd(θe,ζe)·d˙e00000000000+ke,ζζ(de,θe,ζe)·ζe00000000000=fe,ζ(de,θe,ζe)+fse,ζ],
where Ae=1nel is the element assembly operator for element e over the number of elements nel, ce is the element nodal weighting function values for ωh, αe and βe are the element nodal weighting function values for ηh and φh, respectively (both of which are arbitrary except where they are zero at the boundaries with essential boundary conditions), de is the element nodal displacement vector, and θe and ζe 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 forD as
(12)C(D)·D˙+K(D)·D=F(D),
where,
(13)D˙=[d˙θ˙ζ˙],D=[dθζ],C(D)=[000Kθd(D)-Kθζ(D)Kθζ(D)Kζd(D)Kθζ(D)Kζθ(D)],K(D)=[Kdd-Kdθ(D)-Kdζ(D)0Kθθ(D)000Kζζ(D)],F(D)=[ffe,d(D)+fte,dfe,θ(D)+fse,θfe,ζ(D)+fse,ζ]=[Fd(D)Fθ(D)Fζ(D)],
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 timetn+1 and derived from difference formulas for Dn+1 and Vn+1, where velocity Vn+1 is D˙(tn+1) and α is the time integration parameter:
(14)C(Dn+1)·Vn+1+K(Dn+1)·Dn+10000=Fn+1(Dn+1),Dn+1=Dn+ΔtVn+a,Vn+a=(1-α)Vn+αVn+1.The form of (14) allows us to consider a semi-implicit integration scheme of a linear form which is written as
(15)C(Dn)·Vn+1+K(Dn)·Dn+1=F(Dn).


## 4. 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 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 semi-implicit, and the triphasic model associated with linear elasticity used mesh of 2D plain strain of nine integration gauss points.Figure5 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.Figure 5
Diagram of the numerical analysis and the experimental test.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.Table 1
Soil parameters for triphasic mixture implementation.
Soil parameters
Value
Young’s modulus,E
1.3 MPa
Poisson ratio,ν
0.4
Solid real density,ρsR
2000 kg/m3
Bulk modulus of air,Ka
0.1 MPa
Water real density,ρwR
1000 kg/m3
Air real density,ρaR
1.2 kg/m3
Initial porosity,n
0.2975
Intrinsic permeability,ϰ
4.5
×
10
-
13 m2
Water viscosity,ηw
1.0
×
10
-
3 Pas
Air viscosity,ηa
1.8
×
10
-
5 PasOne has(16)kra=(1-Se)2(1-Se(2+3λ)/λ),Se=S-Sr1-Sr,krw=1-2.207(1-Se)1.0121,S=1-1.9722×10-11·s2.4279,kαϰηαkrα,α=w,a,
where ηα is the dynamic viscosity and krα is the relative permeability of the α phase, which depends on the relative saturation Sα through the experimental relationship krα=krα(Sα), ϰ is the intrinsic permeability, and the respective degrees of saturation Sw and Sa sum to one: Sw+Sa=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 Sr 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.Comparison of the numerical results and experimental data with pore water pressure and matric suction.
(a)
(b)Comparison of numerical results of air pressure,pa.
(a)
(b)Comparison between the coupled code and the experimental air pressure data.
(a)
(b)Figure 9
Displacement at the top surface of a drainage test in a triphasic mixture.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-n)ρsR+nSρwR+(1-n)(1-S)ρaR and the bulk modulus of the solid grains (106 MPa) and the water (2×103 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 (Ks and Kw) are infinite due to the large values. The averaged density of the mixture is ρ=(1-n)ρsR+nSρwR+n(1-S)ρaR.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. Figure8 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 Figure9, 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, kw, 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.


## 5. 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].



---

*Source: 589452-2013-09-18.xml*