A fully coupled hydro-mechanical model for the analysis of partially saturated multiphase geo-structures using the bounding surface plasticity

. This paper presents a fully coupled plasticity model for the rigorous analysis of fluid flow and solid skeleton deformation in variably saturated porous media. The flow deformation behaviour of unsaturated soils is described through the simultaneous solution of governing partial differential equations, including the equilibrium equation and the continuity equations of pore fluids. The coupling between solid and fluid phases is established based on the effective stress concept with the effective stress parameter as a function of suction. The elasto-plastic behaviour of geomaterials is formulated using the bounding surface theory. The fully coupled flow deformation framework along with the UNSW bounding surface constitutive law are implemented into COMSOL Multiphysics commercial software to simulate the results of several laboratory tests conducted on soils with different degrees of saturation under various drainage conditions and loading paths. The close agreement between the numerical solutions and experimental data from the literature indicates the good performance of this new elastic-plastic unsaturated model and emphasises its capability to capture the different characteristic features of response of soils with a variety of saturation degrees and over-consolidation ratios.


Introduction
Experimental and theoretical investigations of the hydro-mechanical behaviour of unsaturated geomaterials have been a subject of great interest amongst geotechnical researchers in the past few decades [e.g., [1][2][3][4][5]. As a multiphase porous system, the response of a partially saturated soil subjected to the mechanical and hydraulic loadings is generally governed by the deformation of solid skeleton, the simultaneous flow of water and air through the porous medium, and the coupling between the fluid flow and soil deformation. The interrelations between the suction and the effective stress, and the evolution of soil water characteristic curve with the matrix deformation, are additional mechanisms contributing to the complexity of the hydro-mechanical response in unsaturated soils.
The physical processes involved in the flow deformation problems are mathematically described by a set of coupled partial differential equations, including the force equilibrium equation and the mass conservation law. Due to the complexity and the highly non-linear nature of the governing equations, the analytical solutions are difficult, if not almost impossible, to obtain, and the numerical simulations are frequently employed as an effective tool for the hydromechanical analysis of partially saturated deformable porous media. COMSOL Multiphysics [6] provides a platform for tackling the complex coupling systems in an efficient and reliable manner. It renders coupled modelling capabilities, suitable for the numerical modelling of the hydro-mechanical problems in partially saturated geo-structures available to practising engineers.
This paper presents a fully coupled framework based on the theory of bounding surface plasticity to predict the flow deformation response of unsaturated porous materials. After introduction, the governing partial differential equations, the main concepts of the UNSW bounding surface model, and its implementation in COMSOL are briefly described. The capability of the model implemented is investigated via comparing the simulation results with the experimental data available in the literature.

Governing equations
The flow and deformation phenomena in partially saturated geo-materials can mathematically be expressed by the following set of partial differential equations [5]: where is the elasto-plastic stiffness matrix, u is the soil skeleton displacement vector, and are the pore water and air pressures, respectively, is the Kronecker delta vector, F is the body force per unit volume, ⁄ , is the saturation degree, n is the porosity, is the volumetric strain, and are the relative permeabilities of pore water and air, respectively, k is the intrinsic permeability of soil, and are the dynamic viscosities of water and air, respectively, g is the gravitational acceleration vector, and are the mass densities of water and air, respectively, ⁄ , is the matric suction, , , and are the water and air compressibility coefficients, respectively, and and denote the volumetric contents of water and air phases, respectively.
In this formulation, the effective stress representation by Khalili and Khabbaz [7], in which the effective stress parameter ( ) depends on suction, is adopted: (4) where 1 for . for (5) is the effective stress tensor, is the net stress tensor, is the total external stress tensor, and is equal to air entry value ( ) for the drying path and to air expulsion value ( ) for the wetting path.
In addition, the soil water characteristic curve is defined using Brooks and Corey [8] equation: 1 for for (6) where 1 ⁄ is the effective saturation degree, is the saturation degree, is the residual saturation degree, and is the pore size distribution index.

Bounding surface plasticity model
The elastic-plastic stiffness matrix ( ) is formed based on the bounding surface plasticity theory. Specifically, the constitutive law, referred to as the UNSW bounding surface plasticity model, is adopted to simulate the stress-strain behaviour of cohesive and granular soils subjected to complex monotonic and cyclic loading paths [e.g., 5,9]. The model permits a smooth transition between elastic and elasto-plastic responses and captures accumulation of plastic strains in the so-called "elastic region" [5]. Within this context, the shape of the bounding surface is described by: in which ̅ , and ̅ are the mean effective stress, deviatoric stress and Lode angle on the bounding surface, ̅ is the isotropic pre-consolidation pressure which varies with suction and plastic volumetric strain ( ) [5], is the slope of critical state line (CSL) in the plane, and and are material constants controlling the shape of bounding surface.
The critical state in the plane is represented by a line passing through the origin, with the slope of . In the ln plane, the critical state is defined as a line parallel to the limiting isotropic compression line (LICL), with a constant shift along the recompression line. The isotropic compression line for unsaturated soils is given by: where is the specific volume on the LICL, is the intercept of LICL at 1 kPa, and is the slope of LICL.
The loading surface assumes the same shape as the bounding surface, with the centre of homology located on the origin of stress coordinate system, for first time loading, and moving to the last point of stress reversal, for unloading/reloading conditions. The current stress state always lies on the loading surface and its image on the bounding surface can be found using a mapping rule [10], such that the normal to the loading surface at stress point and to the bounding surface at image point are the same.
The plastic potential determines the direction of plastic strain increment vector and can be stated as: , , , ̅ 1 for 1 (9) , , , ̅ ln for 1 (10) where , and are the mean effective stress, deviatoric stress and Lode angle at the current stress point, is a dummy parameter controlling the size of the plastic potential, ̅ is the loading direction multiplier, and is a material constant.
The hardening modulus at the current stress includes two terms: where ℎ is the plastic modulus at the image point on the bounding surface, and ℎ is some arbitrary modulus at the stress point on the loading surface, depending on the distance between the current stress and the image point. These hardening moduli are defined as: where , ̂ defines the size of loading surface, 1 2 , is the specific volume at , is the specific volume on the critical state line corresponding to , ⁄ , and is a material parameter. Accordingly, the elastic-plastic stiffness matrix is expressed in the following general form as: (14) where is the elastic stiffness matrix, ⁄ ‖ ⁄ ‖ ⁄ is the unit normal vector to the bounding surface at the image point, and ⁄ ‖ ⁄ ‖ ⁄ denotes the unit normal vector to the plastic potential at the stress point.

Implementation in COMSOL
The numerical solution of the governing differential equations is obtained through implementing the elastoplastic model proposed into COMSOL software. To this end, the physical processes involved in the hydromechanical analysis of partially saturated soils, namely: the mechanical deformation, and water and air flow in the porous medium, are accounted for through adopting the relevant modules and physics interfaces in COMSOL. They include the generic "Solid Mechanics" and "Darcy's Law" physics interfaces. These interfaces are then coupled manually in order to solve a simultaneous system of equations.
Moreover, in order to add the UNSW bounding surface plasticity model into COMSOL, a dynamic link library (DLL) for this constitutive law is created in C language, which is accessed using the External Material feature.
The fully coupled flow deformation framework developed in COMSOL, along with the constitutive model implemented, constitute a robust computational tool for the hydro-mechanical analysis of variably saturated multiphase geo-structures based on the plasticity theory. The intention is to make available the use of advanced numerical models for unsaturated soils for practical applications.

Model application
In this section, the fully coupled plastic model proposed is validated through simulating the results of experimental tests conducted under a range of saturation and drainage conditions.

Drained triaxial tests on Weald clay
Model performance in capturing the response of saturated soils during drained loading is examined using the results of compression triaxial tests on normally consolidated and heavily over-consolidated samples of Weald clay [11]. The initial conditions for samples used in the simulations, including the initial mean effective pressure and the initial void ratio, were: 207 and 0.632 for the normally consolidated sample, and 34.5 and 0.617 for the overconsolidated sample. The material parameters used in the analyses are summarized in Table 1.
Comparisons between the model results and experimental data are presented in Figures 1 and 2. It is seen that the proposed framework can reproduce the drained strength and volumetric behaviours of both normally consolidated and over-consolidated clays with acceptable accuracy.

Undrained triaxial tests on Cardiff kaolin clay
To demonstrate the capability of proposed flow deformation plasticity model to capture the response of saturated soils subjected to undrained loading, the results of the triaxial compression tests on a soft kaolin clay reported in [12] are analysed. The tests involved samples with three different over-consolidation ratios equal to 1, 2 and 5, with the material parameters presented in  Model simulation results in terms of the effective stress path, and the deviatoric stress evolution with the axial strain are shown in Figures 3 and 4. The good agreement between model predictions and experimental data for all three cases emphasises the ability of developed numerical model to describe the undrained behaviour of clays.

Unsaturated triaxial compression tests on Bourke silt
The application of the model to the analysis of variably saturated soils is indicated using the results of a series of suction controlled triaxial tests performed on samples of silt from Bourke region of New South Wales, Australia [13]. These tests involve: 1) pre-consolidating the samples to an isotropic stress of 200 kPa, 2) unloading them to the initial net stress of 50 kPa, and then 3) performing drained triaxial compression tests at constant suction values of 100 and 300 kPa. Table 3 summarises the corresponding bounding surface and unsaturated model parameters for Bourke silt. Furthermore, the variations of two material parameters and λ(s) with suction are presented in Table 4. The results obtained from the numerical modelling of this set of tests along with the experimental data are displayed in Figures 5 and 6. It is observed that the presented model can predict both deviatoric and volumetric responses of the soil with a reasonable accuracy compared to the test data. As expected, the suction effect has led to an increase in the soil strength and more volumetric strain under the applied deviatoric load.

Conclusions
In this paper, a bounding surface plasticity framework for the hydro-mechanical analysis of flow deformation phenomena in deformable variably saturated porous media is presented. The model is formulated based on the effective stress principle, with the effective stress parameter defined as a function of suction. The processes involved in this fully coupled phenomena (including the geo-mechanical deformation, and coupled water and air flow) were described mathematically using the mass and momentum balance equations. The interaction between the various processes in the system is established through the effective stress equation and the water retention curve.
To capture the elasto-plastic deformation behaviour of geo-materials, the UNSW bounding surface plasticity is adopted. The stiffening effect of suction on the response of the solid skeleton is accounted for in the definition of isotropic compression line, and consequently, the parameter controlling the size of bounding surface. The model is implemented in the commercially available COMSOL software for practical applications. The model implemented provides a robust computational tool for analysing coupled flow deformation problems in unsaturated porous media within an elastic-plastic setting The performance of the model was investigated through comparing the numerical simulations with experimental data from the literature. Good agreement between the results highlights the capabilities of the model in capturing the different aspects of the behaviour of soils for a range of saturation states and loading conditions.