Numerical Implementation of a Hydro-Mechanical Coupling Constitutive Model for Unsaturated Soil Considering the Effect of Micro-Pore Structure

: An unsaturated soil constitutive model considering the influence of microscopic pore structure can more accurately describe the hydraulic–mechanical behavior of unsaturated soil, but its numerical implementation is more complicated. Based on the fully implicit Euler backward integration algorithm, the ABAQUS software is used to develop the established hydro-mechanical coupling constitutive model for unsaturated soil, considering the influence of micro-pore structure, and a new User-defined Material Mechanical Behavior (UMAT) subroutine is established to realize the numerical application of the proposed model. The developed numerical program is used to simulate the drying/wetting cycle process of the standard triaxial specimen. The simulation results are basically consistent with those calculated by the Fortran program, which verifies the rationality of the developed numerical program. the Fortran program of the model is into ABAQUS to create a UMAT subroutine for secondary development to realize the application of the model. The drying–wetting cycle test of a triaxial sample model is used to demonstrate the simulation performance of the developed numerical program. The results indicate that the developed program can reasonably describe the distribution of pore water pressure and void ratio during the wetting–drying process of a sample model. The prediction results of matric suction and void ratio during the drying–wetting process from the developed program are basically consistent with the prediction results from the Fortran program, indicating that the developed program can effectively describe the unsaturated soil hydro-mechani-cal coupling characteristics during a drying–wetting cycle.


Introduction
The hydro-mechanical coupling characteristics of unsaturated soils are closely related to many geotechnical engineering problems, such as the collapse of unsaturated soil slopes, the uneven settlement of unsaturated subgrade soil, and the movement of landfill leachate. Different from saturated soil, unsaturated soil is a three-phase mixture of solid, pore water and pore air, and its hydraulic-mechanical properties are complex. It is of great significance to establish a constitutive model to accurately describe the hydro-mechanical coupling characteristics of unsaturated soils.
Many scholars have conducted research on unsaturated hydro-mechanical coupling constitutive models and obtained valuable results: Wheeler et al. [1] found that the change in the degree of saturation r S will have an important effect on the unsaturated soil mechanical properties in the case of the same matric suction and mean net stress. At the same time, the effective stress p * and modified suction s * are used as stress variables, the conjugate variable volumetric strain v ε and the degree of saturation r S are used as strain variables, a coupled elastic-plastic model (the Glasgow coupled model, GCM) of the stress-strain behavior and liquid phase lag of unsaturated soil under isotropic stress state was proposed, and the coupling mechanism of the deformation characteristics and water retention characteristics was analyzed. Marti et al. [2] extended Wheeler's model above to a three-dimensional stress state and elaborated several coupling mechanisms under a three-dimensional stress state. Sun et al. [3] established a three-dimensional unsaturated soil hydro-mechanical coupling elastoplastic model. This model considers the effect of suction on the hydraulic and mechanical properties of unsaturated soil, the influence of the r S on the stress-strain relationship and strength, and the effect of soil deformation on the soil-water characteristic curve (SWCC). Sheng et al. [4] first discussed coupled hydro-mechanical models from the perspective of thermodynamics and proposed a coupled hydro-mechanical model that can smoothly transition from an unsaturated state to a saturated state. In the framework of the constitutive model of continuous soil mechanics, Liu et al. [5] established an unsaturated hydro-mechanical coupling constitutive model: the generalized effective stress, modified suction and gas pressure are used as state variables. Zhou et al. [6] established a new unsaturated hydro-mechanical coupling constitutive model by using the degree of capillary saturation and the effective inter-particle stress. Li et al. [7] introduced the uploading surface to establish the hydro-mechanical coupling constitutive model of overconsolidated, unsaturated soil. Mahmoodabadi et al. [8] established an unsaturated hydro-mechanical coupling constitutive model by redefining the hysteresis and elastic shear strain components of the Sheng, Fredlund, Gens (SFG) model. Dat et al. [9] established a new unsaturated hydro-mechanical coupling constitutive model based on the generic thermodynamics.
Most of the existing constitutive models do not consider the influence of pore structure on the unsaturated soil's hydro-mechanical coupling characteristics. Some scholars [10,11] found that the pore structure of unsaturated soil will change under loading or during a drying/wetting cycle, which is more obvious for soils with pore structures at two scales. Compared with intra-aggregate pores, the inter-aggregate pores in a soil with two pore structures are more sensitive to the interaction of water and force. Based on the GCM model proposed by Wheeler, a modified unsaturated hydro-mechanical coupling constitutive model, considering the influence of microscopic pore structure (the modified Glasgow coupled model, MGCM), was established for the isotropic state by Cai et al. [12]. In this model, the influence of macro-and micro-pores on the hydro-mechanical characteristics is considered by replacing r S with e S .
When using the constitutive model to solve the problem-that is, the boundary value problem-it is necessary to use a suitable integration algorithm to numerically implement the model at the Gaussian point scale. According to different integration modes, the integral algorithm of the unsaturated soil constitutive equation can be divided into implicit integrals and explicit integrals. The explicit constitutive integral algorithm adopts the forward Euler integral method to calculate the yield surface and the plastic potential gradient by using the stress at the starting point of each incremental step, which can be directly given by the constitutive matrix. Generally, an incremental step in the display algorithm will be divided into smaller substeps, and appropriate error control technology is used to ensure that the error is within the allowable range. Sheng et al. [13], Sánchez et al. [14], Sołowski et al. [15] and Zhou et al. [16] adopted an explicit algorithm to analyze constitutive models of unsaturated soil. For the implicit algorithm, the yield surface and plastic potential gradient are calculated by using the stress of the unknown stress point in the current step, and the constitutive matrix is also unknown. Therefore, numerical iteration is required to solve the constitutive equation. The implicit algorithm has a faster convergence speed and produces a more accurate calculation than an explicit algorithm. Zhang et al. [17,18], Hoyos et al. [19], Liu et al. [20], Sołowski et al. [21], Enrico [22], Barbara [23] and Ma [24] adopted an implicit algorithm to analyze unsaturated soil. Hu et al. [25] and Sołowski et al. [26] used the two algorithms mentioned above to numerically implement a constitutive model of unsaturated soil and then compared and analyzed the performance of the two algorithms. At present, there are few reports on the numerical implementation of the constitutive model, considering the influence of the unsaturated soils' micro-pore structure. Considering the influence of the micro-pore structure of unsaturated soils is of great significance for improving the accuracy of the numerical analysis of related engineering problems.
In this paper, the established MGCM [8] is summarized first. Then, based on Euler's backward implicit integration algorithm, the Fortran calculation program of the proposed model is imported into ABAQUS for secondary development. The numerical application of the model is implemented. The developed numerical program is used to simulate the wetting-drying-wetting paths test of a triaxial sample, and the rationality of the developed numerical program is verified by comparison with the calculation results of the Fortran program.

Modified Glasgow Coupled Model (MGCM) Overview
Considering the interaction between the changes in the microscopic pore structure and the soil behavior of unsaturated soil with two pore structures, Cai. et al. [12] introduced the e S into the p * in the GCM. Cai et al. established an unsaturated hydro-mechanical coupling constitutive model considering the influence of microscopic pore structure for an isotropic stress state and extended the model to the three-dimensional stress space. The model is described as follows.

Effective Degree of Saturation Considering the Influence of the Microscopic Pore Structure.
The e S can be defined as follows to consider the effect of the microscopic pore structure [27]: where m r S is the r S of the intra-aggregate pores, which is numerically equal to the residual degree of saturation res S of the unsaturated soil. The e S essentially reflects the degree to which the macropores are filled by water, and its value range is (0,1). For a soil with a single pore structure, there is no difference in the saturations between the large and small pores-that is, m r =0 S , m e r = S S . For soil with two pore structures, the pores at each scale account for a certain proportion of the total pore volume. With the increase in clay content, the degree of micro-pore saturation increases gradually, the large pores are gradually replaced by small ones, and the saturation limit is reached. At this time, Using the e S instead of the r S to establish the model can distinguish the contribution of different types of pores to the soil deformation.

Constitutive Variables
Based on the energy equation, taking into account the conjugation of the stress and strain variables into energy, the effective stress [28] p * considering the micro-pores and the modified suction s * are used as the two stress variables: where p is the mean net stress, p is the total stress, s is the matric suction, and n is the porosity. Regardless of whether the suction is zero, as long as the sample is saturated, p * can always return to the effective stress expression of the saturated soil.
From Equations (2) and (3), the stress variable increments can be obtained as follows: where v dε and e dS are the conjugate strain increments corresponding to the p * and s * , respectively.

Yield Curves
The LC, SI and SD yield curves for isotropic stress states are considered approximately straight lines, as shown in Figure 1. LC is the loading collapse yield curve in * * -p s plane, SI is the suction increase yield curve and SD is the suction decrease yield curve. Only elastic strains occur in the rectangular region enclosed by the three yield curves.

Hardening Law
Under hydro-mechanical coupling, the evolution law of each yield curve-namely, the hardening law-can be obtained through the interaction of solid and liquid phases. Assume that the coupling mechanism can be expressed as the following three cases.
Plastic volumetric strain occurs only when the stress point reaches the LC yield curve, causing the coupling of the SI and SD curves to move up. 2 k is the coupling parameter, and the coupled movements equations of the SI and SD yield curves are given as follows: The plastic r S decreases when the stress point reaches the SI yield curve, causing the coupling of the SD curve to move up and the LC curve to move out. The plastic r S increases when the stress point reaches the SD yield curve, causing the coupling of the SI The coupling relationship can be given by: In general, the movement of the LC yield curve is determined by its own direct movement and the coupled movement of the SI or SD curves. The hardening equation of the solid phase is given by: where λ is the saturated normal compression line slope and w Similarly, the liquid phase hardening equation is given by:

Flow Rule
One or more plastic mechanisms are activated during the yield process. The flow rule for the LC, SD and SI yield curves corresponds to the associated flow rule. We assume that only plastic volumetric strain occurs when yielding on the LC yield curve but causes no changes in the Similarly, we assume that only plastic p e S occurs when yielding on the SI or SD yield curves but causes no changes in the plastic volumetric strain. The associated flow rule applied on the SD or SI yield curve is given as follows:

Stress-Strain Relationship
The incremental expressions of the elastic strain are as follows: where ν is the specific volume: According to the solid-phase hardening equation and liquid-phase hardening equation defined in the MGCM, the plastic strain increment relationship of unsaturated soil can be derived: The total strain increment equation is as follows: ( ) ( )

Fortran Program Prediction and Validation
Cai et al. have verified the predictive ability of the model under isotropic compression conditions, and the predictive ability of the proposed model during a drying-wetting cycle will be verified in this paper. Chen [29] carried out a number of constant mean net stress drying-wetting tests of Wuhan Zaoyang expansive soil by using a computer-controlled triaxial stress-path testing system inherited from Zhan [30]. The triaxial cell is a Bishop and Wesley stress-path cell for testing up to 100 mm diameter cylindrical specimen. The change in suction is achieved by varying pore water pressure. Three sets of drying-wetting test data were selected for model validation.

Solid-Phase and Liquid-Phase Yield Equations
To develop the UMAT subroutine, the stress state of unsaturated soil must be considered. Under natural conditions, most soils will be subjected to effective or net stress, so it is necessary to define a solid-phase yield equation in the * -p q stress plane. In this paper, the modified Cambridge Clay model [31] is defined as the yield equation: where q is the deviator stress, and M is the critical state line slope. This paper mainly considers the soil deformation behavior caused by water-force coupling during wetting-drying processes, so it is also necessary to define a liquid-phase yield equation in the * * -s p stress space. In this paper, Equation (8)

Matrix Form Stress-Strain Relation
To better realize secondary development in the ABAQUS finite element platform, the modified suction is regarded as a stress component, and the e S is regarded as a strain component. Three positive stress components, three shear components and the modified suction component are combined into a 1 × 7 generalized stress vector. Three positive strain components, three shear components and the effective degree of saturation component are combined to form a 1 × 7 generalized strain vector. Therefore, the stress and strain tensors are given as the following vectors: ( ) determined by the plastic mechanism of the loading process. The plasticity mechanism is divided into yielding on LC only, yielding on SI/SD only and yielding on LC and SI/SD at the same time.

Implicit Integration Algorithm
In each iteration, after obtaining the displacement increment u Δ , suction increment s Δ and degree of saturation increment r S Δ , the steps to determine the new stress state are as follows.
(1) Elastic trial calculation The geometric equation is used to calculate the increment of the solid deformation, and the trial value of the stress is calculated according to the pure elastic stress-strain relationship.
is the stress value at the end of the last iteration step.
(2) Yield determination According to the stress test value, the yield function value is calculated: , the test stress value is outside the yield surface, and this state is not allowed. A stress value outside the yield surface should be corrected to coincide with the yield surface. According to the actual yield situation, the plasticity correction should be carried out using Equation (18).

Model Validation
In this section, we use the UMAT subroutine of the MGCM to simulate drying-wetting cycle tests of unsaturated soil triaxial samples for the model verification. By comparing the pore water pressure and void ratio changes in the simulation process with the prediction results of the Fortran program, the practicability of the UMAT subroutine is verified.
The sample is a triaxial sample with a diameter of 40 mm and a height of 80 mm. The cell type is selected as a three-dimensional pore water pressure element with 8 nodes (C3D8P). Figure 3 shows a schematic diagram of the triaxial sample model. The sample material parameters are as follows: The initial hardening parameter of effective stress is 0 =200 kPa p * , and the permeability coefficient of the sample model is 3.6 × 10 − 5 m/s. In this experiment, we set the displacement boundary condition at the initial analysis step to constrain the axial degree of freedom of the bottom surface. The pore water pressure boundary condition is set to a uniform pore water pressure of −100 kPa, and its corresponding r S is 0.1.
The 0 e is 1.186, the three initial normal stresses are −100 kPa, and the three initial shear stresses are zero (tension is positive in ABAQUS). In the consolidation analysis step, an equivalent pressure load of 100 kPa is set to balance the initial stress. During the dryingwetting cycle, a negative pore water pressure (matric suction) is applied on the top of the sample by using the Amplitude function in ABAQUS. The distribution laws of the pore water pressure and void ratio in the initial equilibrium state, wetting process and drying process are shown in Figures 4-6.  The initial equilibrium process of the sample model is similar to the triaxial consolidation process of triaxial samples of the same size under 100 kPa of confining pressure. According to the pore water pressure distribution result and void ratio distribution obtained in the equilibrium analysis step, the applied pressure boundary, displacement boundary, pore water pressure boundary and initial stress are well balanced. During the wetting process, the pore water pressure gradually decreases from the bottom of the sample model upward, and the r S also decreases. At the same time, due to wetting, the matric suction inside the sample model gradually decreases, which causes the effective stress to decrease, so the specific volume of the sample model also increases. During the drying process, the pore water pressure gradually increases from the bottom of the sample model upward, and the r S also increases. At the same time, due to drying, the internal matric suction of the sample model gradually increases, which causes the effective stress to increase, so the specific volume of the sample model continues to decrease.
The relationship between the void ratio and matric suction during the drying-wetting cycle is shown in Figure 7. A comparison between the UMAT subroutine and the Fortran program calculation results is also shown. The calculation results of the UMAT subroutine and Fortran program are basically the same: the sample model produces elastic expansion deformation during the wetting process, and drying process causes the plastic deformation. Therefore, the UMAT subroutine can accurately simulate the unsaturated soil hydro-mechanical coupling characteristics during a drying-wetting cycle.

Conclusions
In this article, we summarize the established unsaturated hydro-mechanical coupling model considering the effect of microscopic pore structure and the drying-wetting test data are used for the model verification. Based on Euler's backward implicit integration 10  UMAT subroutine predictions Fortran program predictions algorithm, the Fortran program of the model is imported into ABAQUS to create a UMAT subroutine for secondary development to realize the application of the model. The drying-wetting cycle test of a triaxial sample model is used to demonstrate the simulation performance of the developed numerical program. The results indicate that the developed program can reasonably describe the distribution of pore water pressure and void ratio during the wetting-drying process of a sample model. The prediction results of matric suction and void ratio during the drying-wetting process from the developed program are basically consistent with the prediction results from the Fortran program, indicating that the developed program can effectively describe the unsaturated soil hydro-mechanical coupling characteristics during a drying-wetting cycle.