Table of Contents Author Guidelines Submit a Manuscript
Mathematical Problems in Engineering
Volume 2013, Article ID 589452, 8 pages
Research Article

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

Department of Civil Engineering, Chonbuk National University, Jeonju 561-756, Republic of Korea

Received 3 April 2013; Revised 13 August 2013; Accepted 15 August 2013

Academic Editor: Piermarco Cannarsa

Copyright © 2013 You-Seong Kim and Jaehong Kim. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.


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 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 [27]. 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 where the porosity . If a material is homogeneous, , whereas if it is heterogeneous, the volume fraction at a material point for a differential volume of the mixture. The partial mass density of the phase is given by , and thus 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] 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] 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, where is two constituents (water and air), the fourth-order elastic modulus tensor is , the Lame parameters are and , 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 at sea level and at C and the air bulk modulus is  Pa at a constant temperature. The relative permeabilities and as a function of are given in Figure 3.

Figure 3: Relative permeabilities (a) and (b) plotted against when varying from to .

It can be shown that , where is a parameter of dimension area () 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 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 where is the element assembly operator for element over the number of elements , 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 where, 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 and 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 and derived from difference formulas for and , where velocity is 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

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 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: . 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 69. 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.

Figure 6: Comparison of the numerical results and experimental data with pore water pressure and matric suction.
Figure 7: Comparison of numerical results of air pressure, .
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 and the bulk modulus of the solid grains ( MPa) and the water ( 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 .

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


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.


  1. D. G. Fredlund and H. Rahardjo, Soil Mechanics for Unsaturated Soils, John Wiley & Sons, New York, NY, USA, 1993.
  2. B. A. Schrefler and Z. X. Zhan, “A fully coupled model for water flow and airflow in deformable porous media,” Water Resources Research, vol. 29, no. 1, pp. 155–167, 1993. View at Publisher · View at Google Scholar · View at Scopus
  3. N. Khalili and S. Valliappan, “Unified theory of flow and deformation in double porous media,” European Journal of Mechanics A, vol. 15, no. 2, pp. 321–336, 1996. View at Google Scholar · View at Scopus
  4. N. Khalili, R. Witt, L. Laloui, L. Vulliet, and A. Koliji, “Effective stress in double porous media with two immiscible fluids,” Geophysical Research Letters, vol. 32, no. 15, Article ID L15309, 2005. View at Publisher · View at Google Scholar · View at Scopus
  5. W. G. Gray and B. A. Schrefler, “Thermodynamic approach to effective stress in partially saturated porous media,” European Journal of Mechanics A, vol. 20, no. 4, pp. 521–538, 2001. View at Publisher · View at Google Scholar · View at Scopus
  6. W. G. Gray and B. A. Schrefler, “Analysis of the solid phase stress tensor in multiphase porous media,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 31, no. 4, pp. 541–581, 2007. View at Publisher · View at Google Scholar · View at Scopus
  7. R. I. Borja and A. Koliji, “On the effective stress in unsaturated porous continua with double porosity,” Journal of the Mechanics and Physics of Solids, vol. 57, no. 8, pp. 1182–1193, 2009. View at Publisher · View at Google Scholar · View at Scopus
  8. O. Coussy, Poromechanics, pp. 45–51, 157–168, John Wiley & Sons, New York, NY, USA, 2004.
  9. R. de Boer, Trends in Continuum Mechanics of Porous Media: Theory and Applications of Transport in Porous Media, Springer, New York, NY, USA, 2005.
  10. R. I. Borja, “Cam-Clay plasticity, part V: a mathematical framework for three-phase deformation and strain localization analyses of partially saturated porous media,” Computer Methods in Applied Mechanics and Engineering, vol. 193, pp. 5301–5338, 2004. View at Publisher · View at Google Scholar · View at Scopus
  11. J. Kim, Plasticity modeling and coupled finite element analysis fo partially-saturated soils [Ph.D. thesis], University of Colorado at Boulder, Boulder, Colo, USA, 2010.
  12. R. W. Lewis and B. A. Schrefler, The Finite Element Method in the Deformation and Consolidation of Porous Media, pp. 6–20, chapter 2, John Wiley & Sons, New York, NY, USA, 1987.
  13. M. T. van Genuchten, “Closed-form equation for predicting the hydraulic conductivity of unsaturated soils,” Soil Science Society of America Journal, vol. 44, no. 5, pp. 35–53, 1980. View at Google Scholar · View at Scopus
  14. L. Laloui, G. Klubertanz, and L. Vulliet, “Solid-liquid-air coupling in multiphase porous media,” International Journal for Numerical and Analytical Methods in Geomechanics, vol. 27, no. 3, pp. 183–206, 2003. View at Publisher · View at Google Scholar · View at Scopus
  15. T. J. Hughes, The Finite Element Method, pp. 1–51, 57–75, Prentice-Hall, Upper Saddle River, NJ, USA, 1987.
  16. A. C. Liakopoulos, Transient flow through unsaturated porous media [Ph.D. thesis], University of California, Berkeley, Calif, USA, 1965.
  17. B. A. Schrefler and R. Scotta, “A fully coupled dynamic model for two-phase fluid flow in deformable porous media,” Computer Methods in Applied Mechanics and Engineering, vol. 190, no. 24-25, pp. 3223–3246, 2001. View at Publisher · View at Google Scholar · View at Scopus
  18. D. Gawin, L. Simoni, and B. A. Schrefler, “Numerical model for hydro-mechanical behaviour in deformable porous media: a benchmark problem,” in Proceedings of the 9th International Conference on Computer Methods and Advances in Geomechanics, pp. 1143–1148, Wuhan, China, November 1997.
  19. R. N. Brooks and A. T. Corey, “Properties of porous media affecting fluid flow,” Journal of Irrigation Draining Division, vol. 92, pp. 61–88, 1966. View at Google Scholar