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

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 as CO₂ sequestration and air storage in aquifers.

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 Figure 1.

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 (𝑠), water (𝑤), and air (𝑎). Equations regarding the volume and mass of a mixture are defined as the mathematical relations [8, 9]. The volume of a mixture is 𝑣=𝑣_𝑠+𝑣_𝑤+𝑣_𝑎, and the corresponding total mass is 𝑚=𝑚_𝑠+𝑚_𝑤+𝑚_𝑎. Similarly, for the 𝛼 phase, 𝑚_𝛼=𝜌^𝛼𝑅𝑣_𝛼 (nearly homogeneous), where 𝛼=𝑠,𝑤,𝑎 and 𝜌^𝛼𝑅 is the true mass density of the 𝛼 phase. The volume fraction occupied by the 𝛼 phase is given by 𝑛^𝛼=𝑣_𝛼/𝑣, and thus, for water, air, and solid 𝑛^𝑠+𝑛=𝑛^𝑠+𝑛^𝑤+𝑛^𝑎=1, (1) where the porosity 𝑛=(𝑣_𝑤+𝑣_𝑎)/𝑣=𝑛^𝑤+𝑛^𝑎. If a material is homogeneous, 𝑛^𝛼=𝑣_𝛼/𝑣, whereas if it is heterogeneous, the volume fraction at a material point 𝑛^𝛼=d𝑣_𝛼/d𝑣 for a differential volume of the mixture. The partial mass density of the 𝛼 phase is given by 𝜌^𝛼=𝑛^𝛼𝜌^𝛼𝑅, and thus 𝜌^𝑠+𝜌^𝑤+𝜌^𝑎=𝜌, (2) where 𝜌=𝑚/𝑣 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 (𝐮,𝑝_𝑤,𝑝_𝑎) 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] 𝑛div𝝈+𝜌𝐠=𝟎,𝜕𝑆𝑝̇𝜕𝑠^𝑎−𝑝̇_𝑤+𝑆div𝐯=−div𝐯^𝑤,−𝑛𝜕𝑆+𝑛()𝜕𝑠1−𝑆𝐾_𝑎𝑝̇_𝑎+𝑛𝜕𝑆𝑝̇𝜕𝑠^𝑤+()11−𝑆div𝐯=−𝜌^𝑎𝑅div𝜌^𝑎𝑅𝐯^𝑎, (3)(4)(5) where the total stress is written in terms of the effective stress 𝝈^′ (positive in tension) as 𝝈=𝝈^′−𝑝_𝑎𝟏+𝑆(𝑝_𝑎−𝑝_𝑤)𝟏 [10, 12] and the degree of saturation, 𝑆, is defined in the classical form of [13] 𝑆=𝑆_𝑟+𝑆_𝑒1−𝑆_𝑟,𝑆_𝑒=1+(𝛼𝑠)^𝑛−𝑚, (6) 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). 𝐠 is the gravity acceleration vector, 𝐯^𝑤=𝑛^𝑤(𝐯_𝑤−𝐯_𝑠) is the Darcy seepage velocity of water where 𝐯_𝑤 is the true velocity of water, and 𝐯_𝑠 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) where 𝛼 is two constituents (water and air), the fourth-order elastic modulus tensor is 𝐜^𝑒=𝜆𝟏⊗𝟏+2𝜇𝐈, the Lame parameters are 𝜆 and 𝜇, 𝝐=sym𝛁𝐮 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] 𝐯^𝑤=𝑘_𝑤(𝑛,𝑆)−𝜕𝑝_𝑤𝜕𝐱+𝜌^𝑤𝑅𝐠,𝑘_𝑤()=𝜘(𝑛)𝑛,𝑆𝜂_𝑤𝑘_𝑟𝑤(𝑆),𝑘_𝑟𝑤(𝑆)=√𝑆1−(1−𝑆^1/𝑚)^𝑚2, (8)(2) air flow in a porous medium [8] 𝐯^𝑎=𝑘_𝑎(𝑛,𝑆)−𝜕𝑝_𝑎𝜕𝐱+𝜌^𝑎𝑅𝐠,𝑘_𝑎()=𝜘(𝑛)𝑛,𝑆𝜂_𝑎𝑘_𝑟𝑎(𝑆),𝑘_𝑟𝑎(𝑆)=√1−𝑆1−𝑆^1/𝑚2𝑚,𝜌^𝑎𝑅𝑝=exp^𝑎𝐾_𝑎(𝑛)=𝑛,𝛿³1−𝑛², (9)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.2kg/m³ at sea level and at 20^∘C and the air bulk modulus 𝐾_𝑎 is 10⁵ Pa at a constant temperature. The relative permeabilities 𝑘_𝑟𝑤 and 𝑘_𝑟𝑎 as a function of 𝑆 are given in Figure 3.

(a)

(b)

Figure 3 

Relative permeabilities 𝑘𝑟𝑤 (a) and 𝑘𝑟𝑎 (b) plotted against 𝑆 when varying 𝑚 from 0.1 to 1.

It can be shown that 𝜘(𝑛)=𝑙²𝛿(𝑛), where 𝑙² is a parameter of dimension area (m²) 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.

3. 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 𝛁𝐰∶𝝈_𝐵^′−𝑆𝑝_𝑤𝟏−(1−𝑆)𝑝_𝑎𝟏=d𝑣_𝐵𝜌𝐰⋅𝐠d𝑣+_Γ𝑡𝐰⋅𝐭d𝑎,_𝐵𝑤𝜂𝑛𝜕𝑆𝑝̇𝜕𝑠^𝑎−𝑝̇_𝑤d𝑣+_𝐵𝑤𝜂𝑆div𝐯d𝑣−_𝐵𝑤𝛁𝜂⋅𝐯^𝑤=d𝑣_Γ𝑤𝑠𝜂𝑆^𝑤−d𝑠,_𝐵𝑎𝜑𝑛𝜕𝑆𝑝̇𝜕𝑠^𝑎d𝑣+_𝐵𝑎𝜑𝑛(1−𝑆)𝐾_𝑎𝑝̇_𝑎d𝑣+_𝐵𝑎𝜑𝑛𝜕𝑆𝑝̇𝜕𝑠^𝑤+d𝑣_𝐵𝑎𝜑(1−𝑆)div𝐯d𝑣+_𝐵𝑎𝜑𝐾_𝑎𝛁𝑝_𝑎⋅𝐯^𝑎d𝑣−_𝐵𝑎𝛁𝜑⋅𝐯^𝑎=d𝑣_Γ𝑎𝑠𝜑𝑆^𝑎d𝑠, (10) where 𝐰 is the weighting function for the displacement 𝐮, 𝐭 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 𝐮^ℎ is interpolated biquadratically, and the pore water pressure (𝑝^ℎ_𝑤) and pore air pressure (𝑝^ℎ_𝑎) 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 𝐀^𝑛el_𝑒=1𝐜^𝑒𝑇⋅𝐤^{𝑒,𝑑𝑑}⋅𝐝^𝑒−𝐤^{𝑒,𝑑𝜃}𝜽^𝑒,𝜻^𝑒⋅𝜽^𝑒−𝐤^{𝑒,𝑑𝜁}𝜽^𝑒,𝜻^𝑒⋅𝜻^𝑒=𝐟_𝑓^𝑒,𝑑𝐝^𝑒,𝜽^𝑒,𝜻^𝑒+𝐟_𝑡^𝑒,𝑑,𝐀^𝑛el_𝑒=1𝜶^𝑒𝑇⋅𝐤^{𝑒,𝜃𝜁}𝐝^𝑒,𝜽^𝑒,𝜻^𝑒⋅𝜻̇^𝑒−𝐤^{𝑒,𝜃𝜁}𝐝^𝑒,𝜽^𝑒,𝜻^𝑒⋅𝜽̇^𝑒+𝐤^{𝑒,𝜃𝑑}𝜽^𝑒,𝜻^𝑒⋅𝐝̇^𝑒+𝐤^{𝑒,𝜃𝜃}𝐝^𝑒,𝜽^𝑒,𝜻^𝑒⋅𝜽^𝑒=𝐟^𝑒,𝜃𝐝^𝑒,𝜽^𝑒,𝜻^𝑒+𝐟_𝑠^𝑒,𝜃,𝐀^𝑛el_𝑒=1𝜷^𝑒𝑇⋅𝐤^{𝑒,𝜁𝜃}𝐝^𝑒,𝜽^𝑒,𝜻^𝑒⋅𝜻̇^𝑒+𝐤^{𝑒,𝜃𝜁}𝐝^𝑒,𝜽^𝑒,𝜻^𝑒⋅𝜽̇^𝑒+𝐤^{𝑒,𝜁𝑑}𝜽^𝑒,𝜻^𝑒⋅𝐝̇^𝑒+𝐤^{𝑒,𝜁𝜁}𝐝^𝑒,𝜽^𝑒,𝜻^𝑒⋅𝜻^𝑒=𝐟^𝑒,𝜁𝐝^𝑒,𝜽^𝑒,𝜻^𝑒+𝐟_𝑠^𝑒,𝜁, (11) where 𝐀^𝑛el_𝑒=1 is the element assembly operator for element 𝑒 over the number of elements 𝑛_el, 𝐜^𝑒 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), 𝐝^𝑒 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 𝐃 as 𝐂(𝐃)⋅𝐃̇+𝐊(𝐃)⋅𝐃=𝐅(𝐃), (12) where, ̇𝐝̇𝜽̇𝜻̇𝐝𝜽𝜻,𝐂(𝐃)=𝐊𝐃=,𝐃=𝟎𝟎𝟎_𝜃𝑑(𝐃)−𝐊^𝜃𝜁(𝐃)𝐊^𝜃𝜁(𝐃)𝐊^𝜁𝑑(𝐃)𝐊^𝜃𝜁(𝐃)𝐊^𝜁𝜃(𝐃),𝐊(𝐃)=𝐊^𝑑𝑑−𝐊^𝑑𝜃(𝐃)−𝐊^𝑑𝜁(𝐃)𝟎𝐊^𝜃𝜃(𝐃)𝟎𝟎𝟎𝐊^𝜁𝜁(𝐃),𝐅(𝐃)=𝐟_𝑓^𝑒,𝑑(𝐃)+𝐟_𝑡^𝑒,𝑑𝐟^𝑒,𝜃(𝐃)+𝐟_𝑠^𝑒,𝜃𝐟^𝑒,𝜁(𝐃)+𝐟_𝑠^𝑒,𝜁=𝐅^𝑑(𝐃)𝐅^𝜃(𝐃)𝐅^𝜁(𝐃), (13) where 𝐂 is the combination of the damping matrix and the stiffness matrix of the d.o.f. vector time derivative and 𝐊 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 𝐂, the stiffness matrix 𝐊, and the forcing vector 𝐅, 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 𝐃_𝑛+1 and 𝐕_𝑛+1, where velocity 𝐕_𝑛+1 is 𝐃̇(𝑡_𝑛+1) and 𝛼 is the time integration parameter: 𝐂𝐃_𝑛+1⋅𝐕_𝑛+1+𝐊𝐃_𝑛+1⋅𝐃_𝑛+1=𝐅_𝑛+1𝐃_𝑛+1,𝐃_𝑛+1=𝐃_𝑛+Δ𝑡𝐕_𝑛+𝑎,𝐕_𝑛+𝑎=(1−𝛼)𝐕_𝑛+𝛼𝐕_𝑛+1. (14)

The form of (14) allows us to consider a semi-implicit integration scheme of a linear form which is written as 𝐂𝐃_𝑛⋅𝐕_𝑛+1+𝐊𝐃_𝑛⋅𝐃_𝑛+1=𝐅𝐃_𝑛. (15)

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.

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.

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.

One has 𝑘_𝑟𝑎=1−𝑆_𝑒²1−𝑆_𝑒^{(2+3𝜆)/𝜆},𝑆_𝑒=𝑆−𝑆_𝑟1−𝑆_𝑟,𝑘_𝑟𝑤=1−2.2071−𝑆_𝑒^1.0121,𝑆=1−1.9722×10⁻¹¹⋅𝑠^2.4279,𝑘_𝛼𝜘𝜂_𝛼𝑘_𝑟𝛼,𝛼=𝑤,𝑎, (16) 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.

(a)

(b)

Figure 6 

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

(a)

(b)

Figure 7 

Comparison of numerical results of air pressure, 𝑝𝑎.

(a)

(b)

Figure 8 

Comparison between the coupled code and the experimental air pressure data.

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−𝑛)𝜌^𝑠𝑅+𝑛𝑆𝜌^𝑤𝑅+(1−𝑛)(1−𝑆)𝜌^𝑎𝑅 and the bulk modulus of the solid grains (10⁶ MPa) and the water (2×10³ 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.

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

Acknowledgments

This work was supported by the Energy Efficiency & Resources of the Korea Institute of Energy Technology Evaluation and Planning (KETEP) Grant funded by the Korea Government Ministry of Knowledge Economy (no. 20122020200010) and by research funds of Chonbuk National University in 2013.