Elastoplastic Model Framework for Saturated Soils Subjected to a Freeze–Thaw Cycle Based on Generalized Plasticity Theory

The failures of soil slopes during the construction of high-speed railway caused by the soil after the freeze–thaw (F–T) cycle and the subsequent threat to construction safety are critical issues. An appropriate constitutive model for soils accurately describing the deformation characteristics of soil slopes after the F–T cycle is very important. Few constitutive models of soils incorporate the F–T cycle, and the associated flow rule has always been employed in previous models, which results in an overestimation of the deformation of soil exposed to the F–T cycle. Generalized plasticity theory is widely used to predict the performance of geotechnical materials and is especially well adapted to deal with this type of generalized cyclic loading (such as a freeze–thaw cycle), and it overcomes the shortcomings of the associated flow rule that causes larger shear deformation. To this end, an elastoplastic model framework based on generalized plasticity theory with double yield surfaces for saturated soils subjected to F–T cycles was developed. Two types of plastic deformation mechanisms, i.e., plastic volumetric compression and plastic shear, were considered in this elastoplastic model. It was found that this model can accurately predict the mechanical behavior and deformation characteristics of saturated soils after F–T cycles.


Introduction
High-speed railways in China have developed rapidly over the last decade [1]. By the end of 2021, the total length of the high-speed railway systems in China is expected to be more than 40,000 km. A high-speed railway has been constructed in a typical seasonally frozen region of Northeast China [2,3]. The freeze-thaw (F-T) cycle is one of the most common physical weathering processes in these regions [4] and significantly affects the mechanical behavior of soils [5][6][7][8][9][10][11][12][13]. Some studies have reported that F-T cycles induce cutting slope failures ( Figure 1) [14] during the construction of high-speed railway. An appropriate constitutive model for soils that accurately describes the deformation characteristics of soil slope after the F-T cycle is very important.
Recently, many studies have been carried out to analyze the effect of wet/dry cycles on the deformation behavior and constitutive stress-strain of soil [15][16][17][18][19]. Furthermore, a number of researchers, such as Li et al. 2019 [20], Li and Yang, 2018 [21], Guo et al. 2021 [22], and Li et al. 2020 [23], have studied the elastoplastic model of unsaturated soils, incorporating different influencing factors. However, the literature about the constitutive model of saturated soils considering the effect of the F-T cycle is relatively sparse. Furthermore, the associated flow rule was usually adopted in the constitutive model of soils [24,25]. The sparse. Furthermore, the associated flow rule was usually adopted in the constitutive model of soils [24,25]. The traditional associated flow rule always induces larger shear deformation [26,27], which overestimates the deformation of soil slopes exposed to the F-T cycle.
Generalized plasticity theory is widely used to predict the performances of geotechnical materials [28,29]. The effect of freezing-thawing cycles on the behavior of soils is a matter of great practical importance. Generalized plasticity is especially well adapted to deal with this type of generalized cyclic loading [29]. The generalized plasticity theory does not obey the traditional assumption of plastic potential function, nor does it comply with the associated plastic flow rule [28]. It overcomes the shortcomings of the associated flow rule causing larger shear deformation. As a result, the use of the framework of generalized plasticity theory is particularly useful to evaluate the mechanical behavior of soil. However, few models for soils subjected to the F-T cycle have been based on the generalized plasticity theory.
In this study, an elastoplastic model for saturated soils subjected to F-T cycles under the framework of generalized plasticity theory with double yield surfaces was developed. The derivation of the model and the estimation of the model parameters are presented in detail. The triaxial test results of the saturated soils were predicted by the proposed model, and finally, the effectiveness of the model was verified using the test data.

An Elastoplastic Model Framework for Saturated Soils in Generalized Plasticity Theory Incorporating the Effects of the Freeze-Thaw Cycle
An elastoplastic model framework for saturated soils subjected to F-T cycles was developed using the generalized plasticity theory. Moreover, a detailed introduction to the generalized plasticity theory is presented in the Appendix.

Elastic Deformations
In this section, the elastic deformations are discussed, and the soil is regarded as isotropic and elastic within the yield surface. The increments in the strain invariants consist of two parts [30]

An Elastoplastic Model Framework for Saturated Soils in Generalized Plasticity Theory Incorporating the Effects of the Freeze-Thaw Cycle
An elastoplastic model framework for saturated soils subjected to F-T cycles was developed using the generalized plasticity theory. Moreover, a detailed introduction to the generalized plasticity theory is presented in the Appendix A.

Elastic Deformations
In this section, the elastic deformations are discussed, and the soil is regarded as isotropic and elastic within the yield surface. The increments in the strain invariants consist of two parts [30], i.e., dε s = dε e s +dε p s (1) where dε s , dε e s , and dε p s are the increment of shear strain, elastic shear strain, and plastic shear strain, respectively; and dε v , dε e v , and dε p v are the increment of volumetric strain, elastic volumetric strain, and plastic volumetric strain, respectively.

Yield Surfaces and Plastic Potential Functions
The deformation of soil slope after the F-T cycle, which includes the shear deformation, compression deformation, or a combination of the two deformations, is complicated. The double yield surfaces proposed by Yin (1988) [31] could reflect two types of plastic deformation mechanisms, namely, plastic volumetric compression and plastic shear for soils, and it is often employed by researchers [24,32] to present the mechanical and deformation characteristics of soils. Therefore, the double yield surfaces proposed by Yin (1988) [19] were used in this paper. Figure 2a shows the two yield surfaces proposed by Yin (1988) [31] in the q-p plane. Point A is the intersection of the two yield surfaces. Referring to the yield surfaces in [31], the yield surfaces of soils subjected to F-T cycling are plotted in Figure 2b. Two yield surfaces divide the q-p plane into four parts [31]: region Z 0 only has elastic deformation, region Z 1 is only related to the first yield surface, region Z 2 is only related to the second yield surface, and the two kinds of plastic deformation exist simultaneously in region Z 3 .

Yield Surfaces and Plastic Potential Functions
The deformation of soil slope after the F-T cycle, which includes the shear deformation, compression deformation, or a combination of the two deformations, is complicated. The double yield surfaces proposed by Yin (1988) [31] could reflect two types of plastic deformation mechanisms, namely, plastic volumetric compression and plastic shear for soils, and it is often employed by researchers [24,32] to present the mechanical and deformation characteristics of soils. Therefore, the double yield surfaces proposed by Yin (1988) [19] were used in this paper. Figure 2a shows the two yield surfaces proposed by Yin (1988) [31] in the q − p′ plane. Point A is the intersection of the two yield surfaces. Referring to the yield surfaces in [31], the yield surfaces of soils subjected to F-T cycling are plotted in Figure 2b. Two yield surfaces divide the q − p′ plane into four parts [31]: region Z0 only has elastic deformation, region Z1 is only related to the first yield surface, region Z2 is only related to the second yield surface, and the two kinds of plastic deformation exist simultaneously in region Z3. In the first yield surface (fv), the loading-collapse (LC) yield surface represents the compression mechanism and is expressed as the following form: In the first yield surface (f v ), the loading-collapse (LC) yield surface represents the compression mechanism and is expressed as the following form: where M 1 , p 0 , and p r are the LC yield surface indexes that are related to the shape of the stress-strain curve of soils, the initial mean effective stress, and the intercept of the failure line on the p axis, respectively. In addition, symbols M i , p i , p r,i , and p 0,i are the parameters under various numbers of freeze-thaw cycles (N FT ), which are the slope of failure line, mean effective stress, intercept of the failure line on the p axis, and initial mean effective stress, respectively. Subscript i represents the number of freeze-thaw cycles (i = 0, 1, 3, 7, and 11), which is the same meaning as N FT . For convenience, the symbol i was used in this section.
A hyperbolic curve is used to predict the relation between plastic volumetric strain (ε p v ) and p 0 [31]: where χ 1 , χ 2 , and p a are the elastoplastic compression indexes and air pressure, respectively. The ε p v was chosen as the hardening parameter. A non-associated plastic flow rule was adopted in this model. Referring to [33], the mean effective stress p is assumed to be the plastic potential function (Q v ) corresponding to the plastic volumetric strain.
The second yield surface (f q ) proposed by Yin (1988) [31] is employed to present the shearing mechanism of the soils: where a and M 2 are the elastoplastic dilation index and shear yield surface index that is related to the slope of failure line (M), respectively. ε p s is the plastic shear strain. Similarly, ε p s is used as the hardening parameter for the second yield surface. Referring to [33], the plastic potential function (Q q ) is assumed to be q, corresponding to the plastic shear strain.

Elastoplastic Stress-Strain Relations
The incremental stress-strain relationship of soils in p − q plane can be obtained [28]: where dλ v and dλ q are the plastic factors (see Appendix A). The elastic stiffness matrix, [D e ], can be determined by Equation (13): Based on Equations (9), (11), (12), (A12), and (A13), the following expressions can be obtained [28]: (15) where A v and A q are the plastic parameters. The related parameters are given by Materials 2021, 14, 6485 {dσ}= dp dq (22) {dε}= dε v dε q The elastoplastic incremental stress-strain relationship can be given by According to [28], the elastoplastic stiffness matrix [D ep ] takes the following form: where w 1 and w 2 are given by where t 1 , t 2 , t 3 , t 4 , t 5 , and t 6 can be obtained by Equations (A16)-(A21) (see Appendix A).

Determination of Parameters and Model Validation
The parameters employed in the model fall into three categories: elastic parameters, parameters related to the compressive mechanism, and parameters related to the shear mechanism. The model was employed to predict the test data for soils after the F-T cycle.

Elastic Parameters
The elastic parameters (G and K) can be determined from triaxial compressive loadingunloading-reloading testing. Moreover, G and K can also be obtained from Equations (5) and (6), whereas the value of Poisson's ratio (v) equal to 0.3 is assumed in this model [34,35].

Parameters Related to a Compressive Mechanism
The LC yield surface index, M 1 , which is related to the shape of stress-strain curve, is calculated through the following equation [31]: where β is the ratio of the volumetric strain (ε v ) to axial strain (ε a ) when the stress level equals 75%; the mean value was used in the model under different confining pressure (σ c ) and various N FT values. M is the slope of the failure line, and it is also related to N FT . The value of p r can be obtained from the failure line under various N FT values. The elastoplastic compression indexes (χ 1 and χ 2 ) have been reported by Yin (1988) [31] and will not be repeated here. However, a material parameter Γ was added in the equation Yin 1988 [31]) to represent the material characteristics (it is related to volumetric compression deformation, Γ = 2.5-3.0).

Parameters Related to a Shear Mechanism
The shear yield surface index, M 2 , which is related to the shape of the q-deviatoric strain (ε s ) curve, can be estimated by the following equation [31]: where R f is the failure ratio, which is defined by (30) where q f and q ult are the failure strength and the asymptotic value of deviatoric stress for the q − ε s curve, respectively. The elastoplastic dilation index (a) reflects the proportion of the strain induced by the shear yield surface in terms of the total strain. It can be estimated by the equation proposed by Yin (1988) (a = 0.25 − 0.15d, d is the slope of the ε v − ε a curve at a stress level between 75% and 95%).
In summary, two elastic parameters (G and K), four parameters (M 1 , p r , χ 1 , and χ 2 ) related to the compressive mechanism, and two parameters (M 2 and a) related to the shear mechanism were used in the proposed model. These eight parameters change along N FT and they are all functions of N FT . These parameters can be obtained by test results under different N FT values and empirical formulas.

Model Validation
The influence of the F-T cycle on the mechanical properties of silty clay was investigated by Cui et al. (2015) [36]. The triaxial test results of saturated silty clay exposed to three and five F-T cycles were used to verify the correctness of the model. The volumetric strain was measured under consolidation and shearing in [36]. The volumetric strain induced by consolidation was not given, so the volumetric strain during shearing could not be obtained. Therefore, only the stress-strain curves of the saturated silty clay were employed to verify the model's correctness. The model parameters used in [36] for the third and fifth F-T cycle are shown in Tables 1 and 2. Figure 3 plots the comparison between the experimental data and predicted results. From Figure 3, considering some parameter determination based on the empirical formula, the predicted values are within reasonable range, although there are some differences. In addition, the triaxial test results of bioenzyme-treated saturated soils with 3% bioenzyme from Wen and Wang (2018) [37] were also used to verify the model's correctness. (although the bioenzyme-treated saturated soils are not subjected to the F-T cycle, the validation results can be used as proof that the model framework can be widely used). Model parameters used in [37] are given in Table 3. Form Figure 4, the stress-strain curve and volumetric-deviatoric strain curve of bioenzyme-treated saturated soils were reasonably predicted by the proposed model.

Future Fields of Application
This study presents an elastoplastic model framework for saturated soils subjected to the freeze-thaw cycle and also serves as a simple methodology to establish the constitutive model for soils. Two types of plastic deformation mechanisms, i.e., plastic volumetric compression and plastic shear, were considered in this elastoplastic model. Therefore, this model could be employed to analyze the relatively complex deformations of engineering, such as soil slope, railway roadbed, dam, and building foundation. In addition, the analysis discussed in this paper is mainly used to describe the deformation characteristics of strain hardening for clay soils. In practical applications, this important aspect must be given proper attention. Further study is needed for a more quantitative description of the strain softening of soils.

Conclusions
In this study, an elastoplastic model framework for saturated soils subjected to the freeze-thaw cycle based on the generalized plasticity theory was proposed. The primary conclusions from this work can be given as follows: (1) This model was on the framework of the generalized plasticity theory overcoming the disadvantage of the traditionally associated flow rule, inducing larger shear deformation. The model adequately describes the contractive shear behavior of saturated soils under different freeze-thaw cycles.

Future Fields of Application
This study presents an elastoplastic model framework for saturated soils subjected to the freeze-thaw cycle and also serves as a simple methodology to establish the constitutive model for soils. Two types of plastic deformation mechanisms, i.e., plastic volumetric compression and plastic shear, were considered in this elastoplastic model. Therefore, this model could be employed to analyze the relatively complex deformations of engineering, such as soil slope, railway roadbed, dam, and building foundation. In addition, the analysis discussed in this paper is mainly used to describe the deformation characteristics of strain hardening for clay soils. In practical applications, this important aspect must be given proper attention. Further study is needed for a more quantitative description of the strain softening of soils.

Conclusions
In this study, an elastoplastic model framework for saturated soils subjected to the freeze-thaw cycle based on the generalized plasticity theory was proposed. The primary conclusions from this work can be given as follows: (1) This model was on the framework of the generalized plasticity theory overcoming the disadvantage of the traditionally associated flow rule, inducing larger shear deformation. The model adequately describes the contractive shear behavior of saturated soils under different freeze-thaw cycles.
(2) The capabilities of this model are illustrated by predicting the triaxial tests of silty clay and bioenzyme-treated soils. In this model, the overall elastoplastic model behaviors of saturated soils are described through eight model parameters that are specifically associated with the freeze-thaw cycle.
Further work will focus on the derivation of the elastoplastic constitutive model of soil after the freeze-thaw cycle in a real triaxial stress state and consider unsaturated conditions and more general stress paths.