A finite strain elastoplastic constitutive model for unsaturated soils incorporating mechanisms of compaction and hydraulic collapse

Although the deformation of unsaturated soils has usually been described based on simple infinitesimal theory, simulation methods based on the rational framework of finite strain theory are attracting attention especially when solving geotechnical problems such as slope failure induced by heavy rain in which large a deformation is expected. The purpose of this study is to reformulate an existing constitutive model for unsaturated soils (Kikumoto et al., 2010) on the basis of finite strain theory. The proposed model is based on a critical state soil model, modified Cam-clay, implementing a hyperelastic model and a bilogarithmic lnv-lnP’ (v, specific volume; P’, effective mean Kirchhoff stress) relation for a finite strain. The model is incorporated with a soil water characteristic curve based on the van Genuchten model (1990) modified to be able to consider the effect of deformation of solid matrices. The key points of this model in describing the characteristics of unsaturated soils are as follows: (1) the movement of the normal consolidation line in lnv-lnP’ resulted from the degree of saturation (Q, deviatoric Kirchhoff stress), and (2) the effect of specific volume on a water retention curve. Applicability of the model is shown through element simulations of compaction and successive soaking behavior.


Introduction
Although the stress-strain relationship of unsaturated soils has usually been described based on simple infinitesimal theory, application of the framework of finite strain theory to unsaturated soil mechanics is attracting attention (e.g.[1,2]), especially when predicting geotechnical issues in which large deformation of the ground (such as a failure of slope or embankment owing to heavy rain or earthquake) is expected.
The constitutive framework of finite strain elastoplasticity has been developed [3][4][5][6] based on the multiplicative decomposition of the deformation gradient [7].For geomaterials, several researchers [8][9][10] have proposed stress-update algorithms based on return mapping in principal space [3,11,12].Borja and Tamagnini [8] developed the infinitesimal version of the modified Cam-Clay model for finite strain theory based on a multiplicative decomposition of the deformation gradient ( ep  F F F ).In their formulation, a yield function is defined as a function of the mean Kirchhoff stress (P) and deviatoric Kirchhoff stress (Q) instead of mean Cauchy stress (p) and deviatoric Cauchy stress (q), respectively.In addition, for modeling the behavior of isotropic consolidation of soil, a bilogarithmic lnv-lnP (v, specific volume) relation [13] is applied for some significant advantages [14,15].For unsaturated soil, Song and Borja [1,2] developed a mathematical framework for coupled solid-deformation/fluid-diffusion in a finite strain range.
In this study, we aim to develop a method that can suitably predict the entire life of an embankment from its construction process to the failure stage owing to heavy rain, for instance.In order to achieve this, an infinitesimal constitutive model for unsaturated soils incorporating compaction and collapse mechanisms (Kikumoto et al. [16]) is extended to a finite-strain model that can capture the behavior of unsaturated soils at a large strain level.In the model, the shear strength of unsaturated soil is controlled by the movement of the normal consolidation line in the direction of the specific volume (v) axis with a variation in the degree of saturation (S r ), and a soil water characteristic curve incorporating the effect of changes in the void ratio (e) is adopted.In this paper, we present an outline of a constitutive model for unsaturated soils based on finite strain theory.The applicability of the model is finally discussed through element simulations of compaction and the soaking behavior of unsaturated soils.

Outline of a constitutive model for unsaturated soils
In this section, a model for unsaturated soils incorporating compaction and collapse mechanisms under a finite strain range is derived.This model consists of two parts, namely, an elastoplastic stress-strain relationship taking account of the effect of the degree of saturation, and an advanced soil water retention curve model considering effect of volumetric deformation.

Elastoplastic constitutive model
We first derive a stress-strain relationship for unsaturated soils based on the effective Kirchhoff stress: where J  τ σ (J, Jacobian; σ , Cauchy stress tensor) is the Kirchhoff stress tensor, S r is degree of saturation, and f Ju   ( f u , pore fluid pressure) is the Kirchhoff pore fluid pressure.I is the second-order identity tensor.Subscripts w and a denote water and air, respectively.

Hyperelastic model
In infinitesimal strain theory, a hypoelastic formulation is usually employed.This approach, however, causes an irrational dissipation of stored energy especially under conditions of cyclic loading [17].On the other hand, the hyperelastic model guarantees the conservation of stored energy.Therefore, we apply a hyperelastic model [8] with pressure-dependent bulk and shear moduli to describe the elastic response, which was originally proposed by Houlsby [18] in an infinitesimal strain range and extended to finite strain theory by Borja and Tamagnini [8].
The  )

Yield function
Since unsaturated soils having a lower degree of saturation usually exhibit stiffer behavior, the normal consolidation line (NCL) defined in the ln Pln v plane is assumed to move in the direction of the specific volume v axis owing to the variation of degree of saturation S r (Figure 1) in the proposed model.
lny in Equation ( 7) is the separation between the NCL for current, the unsaturated state, and that for the saturated state in logarithmic scale.lny works as a state variable considering the effect of the degree of saturation r S on the strength of the unsaturated soils.lny is assumed to increase as r S decreases: r ln (1 ) S yx  (9) x is a material parameter representing the vertical distance of the state boundary surface for dried and saturated samples in the compression plane.The volumetric movements of NCL and CSL are also considered in a similar way in the previous work [19].In order to take the effect of density or the overconsolidation ratio into account in the proposed model, we also define a logarithmic volumetric distance ln c (Figure 2).Thus, the specific volume v can be written as follows: The state variable c is assumed to decrease with the development of the plastic deformation and finally converges to 1 as follows: ln a c  c  (11) a is a material parameter controlling the effect of the density or overconsolidation ratio, where  is the increment of plastic multiplier.
From Equation (10), the volumetric strain for isotropic consolidation can be written as where subscript 0 denotes the values at initial state.The elastic component of the volumetric strain for pure isotropic loading is derived from Equation (5) as follows: From Equations ( 12) and ( 13), we obtain the plastic volumetric strain for isotropic consolidation as Equation ( 14): Selecting an elliptic-shaped yield surface as the modified Cam-Clay model [20], the plastic volumetric strain owing to shearing, namely dilatancy, is defined as where M is a material constant and a critical state stress ratio, which controls the aspect ratio of the ellipsoidal yield surface.From Equations ( 14) and ( 15), the yield function can be finally written as where    .For the plastic flow rule, we employ a formulation of the exponential approximation [3,12]: where e,tr  is the principal Kirchhoff effective stress, and t

Water retention curve model
In this study, an extended version of a water retention curve model proposed by van Genuchten [21] is applied.This model is capable of considering the effect of volume change on retention characteristics.The model proposed by van Genuchten is first given as a function of the Kirchhoff suction as: where min r S and max r S are the minimum and maximum degrees of saturation, respectively. , n, m are material parameters, and Js  S (s, suction) is the Kirchhoff suction.In order to consider the effect of the density of soils, we propose to replace S with a modified Kirchhoff suction * S as Equations ( 19) and (20).
  Here, ref v is a reference specific volume, and v x is a material parameter for controlling the effect of the volumetric deformation of soil. Figure 4 shows the water retention curves for three kinds of specific volume (v = 1.90, 1.75, and 1.60).It is properly simulated by the model that soil with a higher density has a higher degree of saturation under the same value of suction, which is reported by experimental results [22].

Return mapping
Next, we derive the return mapping procedure based on the elastic principal logarithmic strains [23].In the return mapping of proposed model, e ,1 An   (A = 1, 2, 3), where  β is the principal effective Kirchhoff stress vector.The backward-Euler approximation is applied to Equation (11).
1 n y  can be calculated from Equation (9)   and (19).Equation ( 21) is nonlinear, so the solution x needs to be evaluated iteratively.To this end, we use Newton's method with a local Jacobian matrix ( ) /  r x x until () rx < TOL is satisfied.
In the elastic regime (i.e., unloading), the above procedure is not implemented.This can be judged by the value of the yield function with trial values tr The trial values except for tr

Algorithmic tangent moduli
An algorithmic tangent moduli that is consistent with a formulation based on the principal space is presented here.The use of the algorithmic tangent moduli is necessary to ensure the convergence of iterations.-98.0 Specific volume v 0 1.85

A differentiation of the plastic flow rule (17) in
From Equations ( 2) and ( 23), the stress-strain relation with the Hessian matrix is derived as [4] where  dc is obtained by differentiating the backward-Euler approximation of Equation ( 11) as follows: Inserting Equations ( 23), (24), and (26) into the consistency condition, we obtain where is the vector in the principal space.The algorithmic tangent moduli   These algorithmic tangent moduli correspond with conventional elastoplastic tangent moduli when 0  .

Simulations
Soaking and compaction behaviors of unsaturated soils are simulated by the proposed model, and the applicability of the model is discussed here.In the simulations presented here, the same sets of material parameters for the water retention curve and elastoplastic model summarized in Table 1 and 2 are adopted.The same set of initial parameters (such as density and the stress state) is also used in the simulations, which are summarized in Table 3.As shown in Table 3, the initial  E 2016 saturated soil in the beginning stage of compression; the degree of saturation increases owing to volumetric compression (even though suction is kept constant); there is rapid compression behavior to the NCL for saturated soil with an asymptotic increase in the degree of saturation to 100%.Such a tendency has been reported in past experimental works by authors such as Wheeler and Sivakumar [24] and Kayadelen [25].
We show a response of the proposed model from soaking.In Figure 6, the soil is soaked under constant net stress (-200, -300, -500, and -800 kPa) after the isotropic compression shown in Figure 5 (case: s = 100 [kPa]).In simulations, the soaking collapse behavior of soils can be seen by decreasing suction s to zero.
The compaction behavior of soils is simulated in Figure 7. Suction s is first increased under a constant net stress ( net 8kPa 9    ) until the prescribed water content w is achieved.Compaction behavior can be regarded as the exclusion of entrapped air without significant drainage of void water.Thus, in this simulation, the total Cauchy stress is increased from -98.0 kPa to a predetermined maximum value (-200, -400, -600, -800, -1000, and -1200 kPa) with constant water content, and returns to -98.0 kPa. Figure 7 shows that the proposed model can simulate the typical compaction behavior of soils.
We validate the proposed model in terms of the shear strength of unsaturated soils.In Figure 8, we simulate triaxial tests under constant suction (drained water and drained air).Suction was first increased to prescribed values (100, 200, and 300 kPa) before shearing.Shearing is then simulated under a constant confining pressure.Figure 7 shows that the proposed model can simulate the difference in the strength against shearing, i.e., unsaturated soil is stiffer than saturated soil.

Conclusion
A model for unsaturated soils that can describe compaction and collapse mechanisms under a finite strain range is presented here.The model is an extended version of an infinitesimal model proposed by Kikumoto et al. [16].Essential concepts of the model in describing the behavior of unsaturated soils are as follows: shifting the NCL in the direction of the specific volume with a change in the degree of saturation, and a dependency of the degree of saturation on specific volume or density.The proposed model is formulated based on Kirchhoff stress invariants in order to predict the large deformation behavior of unsaturated soil based on finite strain theory.It is indicated through the simulations presented in this paper that the proposed model can describe the typical behavior of unsaturated soils such as soaking collapse phenomena, compaction of soils and shearing behavior.

1 lFigure 1 .
Figure 1.Assumed normal consolidation line for unsaturated soils: dependence of the position of NCL on variation of S r .

Figure 2 .
Figure 2. Volumetric behavior of unsaturated soil considering the effects of degree of saturation S r and specific volume v.
specific volume v on the normal consolidation line at ref PP   .l is the compression index representing the slope of the compression line in the ln Pln v relation, which is linked to the compression index l that defines the slope of the compression line in ln p - ln v [8].

Figure 4 .
Figure 4. Soil water characteristic curve for three kinds of specific volumes.

Table 1 .
Material parameters for water retention curveMaximum saturation S is the Hessian matrix. n+1

Figure 5 .
Figure 5. Simulations of isotropic compression of unsaturated soils under constant suction.

Figure 7 .Figure 8 .
Figure 7. Compaction curves for several kinds of applied maximum stresses (final states of compaction simulations) Figure 5 shows the isotropic compression curves of unsaturated soil under constant suction.In this simulation, the suction was first increased to prescribed values (50, 100, and 150 kPa) under constant net stress ( net 8kPa 9    ) and three different values of the degree of saturation are obtained for each case.It is indicated that the proposed model can described the typical behaviors of unsaturated soils where: unsaturated soil is able to stay over the NCL for

Table 2 .
Material parameters for elastoplastic model

Table 3 .
Initial state of soil for simulations 32)