#### Abstract

Artificial frozen soil is a kind of typical creep material, and the frozen clay under the unloading stress paths of high-confining pressure conforms to the improved the Zienkiewicz–Pande parabola-type yield criterion, and the Mohr–Coulomb yield function can describe the shear yield surface of artificial frozen clay under low-confining pressure. Based on the results of triaxial creep and shear tests for artificial frozen soil, the viscoplastic damage variable and evolution rule of artificial frozen clay were obtained by using the theory of viscoelastic-plastic mechanics and damage mechanics. An improved Zienkiewicz–Pande parabola-type yield criterion was used instead of a linear Newton body to obtain a coupled constitutive model of viscoelastic-plastic damage in the frozen soil under the unloading stress paths and to derive the coupling flexibility matrix for viscoelastic and viscoplastic damage. A finite element program of artificial frozen soil considering creep damage was written in the Visual Fortran 6.6A environment and embedded into the nonlinear finite element software ADINA as a user subroutine. The results of numerical simulation and laboratory testing were identical, with a maximum error of no more than 4.8%. This work shows that it is reasonable to describe the creep constitutive model of frozen soil with the viscoelastic-plastic-coupled constitutive model.

#### 1. Introduction

Andersland and Douglas [1] successively proposed the creep theory of frozen soil. A series of constitutive models describing the creep yield process of frozen soil were presented, and a nonlinear yield criterion for frozen soil was proposed by Ladanyi [2, 3]. In the design of an artificial frozen wall during this period, attention was paid to the creep properties of frozen soil. It is found that the relationship between the uniaxial compressive strength of frozen soil and its temperature conforms to the power function [4], and the power function is applicable to both conditions of compression and tension [5–7].

Frozen soil has special static-dynamic characteristics, so it is very important to study its basic mechanical properties [8]. Many scholars have conducted experiments, numerical simulations, and constitutive parameter inversions to study the creep damaged characteristics of frozen soil. However, there have been relatively few studies of the creep damaged characteristics of deep artificial frozen soil under unloading stress paths, especially the numerical calculation of the creep damaged model of artificial frozen clay under these paths [9, 10]. The damage variable of artificial frozen soil and its law of damage evolution can be derived based on the triaxial creep test, and its creep constitutive behavior can be studied using the finite element analysis.

#### 2. Damage Constitutive Model of Artificial Frozen Soil Based on Viscoelastic-Plastic Theory

##### 2.1. Creep Damaged Constitutive Model of Artificial Frozen Soil

Creep tests performed under unloading stress paths of artificial frozen soil show the following [11]:(1)The artificial frozen soil has transient elastic strain under the unloading stress paths, and the instantaneous elastic strain is about 1.0%, less than 10% of the total creep deformation generally.(2)The creep test of deep artificial frozen soil under the lower deviatoric stress shows two stages of unsteady creep and steady creep.(3)The accelerating creep stage occurs when the stress level is above a certain critical value and will end after about 3 h.

The Nishihara model can fully reflect the viscoelastic-plastic constitutive model in geotechnical creep material. Therefore, this model was selected according to the experimental results of artificial frozen soil under unloading stress paths and as the basis for improving the instantaneous elastic deformation and viscoelastic-plastic items, to apply the complex creep deformation law of deep artificial frozen soil. According to the rule of shear and creep deformation of artificial frozen soil under unloading stress path, the modified Zienkiewicz–Pande parabola-type yield surface was analyzed and used instead of the Saint-Venant body model. The damage variable was added for a new viscoelastic-plastic constitutive model suitable for deep artificial frozen soil under unloading stress paths of high-confining pressure, as shown in Figure 1.

For calculation, according to classical rheological finite element theory, the deep artificial frozen soil was regarded as a nonlinear continuum medium. The total strain includes a nonlinear viscoelastic strain component and a nonlinear viscoplastic damage strain component . The nonlinear viscoelastic strain includes transient elastic strain and viscoelastic creep, and the nonlinear viscoplastic damage strain includes viscoplastic creep and damage creep. Therefore, the total strain of the coupled constitutive model under unloading stress paths can be expressed as follows:

The nonlinear viscoelastic strain component and the nonlinear viscoplastic damage strain component are solved below. For numerical calculation, the nonlinear ADINA finite element program was used to simulate the mechanical process of deep shaft excavation and unloading construction. In this program, the secondary development of the constitutive model requires the element flexibility matrix of each component. Therefore, the element flexibility matrix of the nonlinear viscoelastic components and the element flexibility matrix of nonlinear viscoplastic damage were derived below.

##### 2.2. Element Flexibility Matrix of Nonlinear Viscoelastic Body

This section presents the viscoelastic-plastic constitutive model of artificial frozen soil under unloading stress paths of high-confining pressure, based on the test results of artificial frozen soil, in which the model of rheological elements includes the viscoelastic deformation property. The expression of the nonlinear viscoelastic strain component can be obtained directly according to the principle of rheological mechanics:where is the viscoelastic modulus of deep artificial frozen soil (MPa); is the transient modulus (MPa); is the viscosity coefficient (MPa·h); is the load time (h).

The viscoelastic properties of deep artificial frozen soil can be treated as exhibiting instantaneous elasticity in the numerical calculation. By comparing with the generalized Hooke’s law, the viscoelastic flexibility matrix can be given directly bywhere is the equivalent nonlinear viscoelastic modulus and is the sum of the instantaneous elastic modulus and viscoelastic modulus. Based on elastic-plastic mechanics, the nonlinear viscoelastic modulus can be directly given by the following formula:

##### 2.3. Element Flexibility Matrix of the Nonlinear Viscoplastic Body of Artificial Frozen Soil

The last item in the viscoelastic-plastic constitutive model of artificial frozen soil under unloading stress paths of high-confining pressure is the term for nonlinear viscoelastic-plastic damage, as shown in Figure 1. The expression of the nonlinear viscoplastic strain component can be obtained directly from the rheomechanics principle:where is the parabolic yield function of artificial frozen soil, as determined by (6). For frozen soil, ; is the viscosity coefficient (MPa·h); and is load time (h):where is the internal friction angle of artificial frozen soil under unloading stress paths (°), determined by the triaxial test, and is the cohesive force (MPa), determined by triaxial test.

In order to improve the calculation efficiency and reliability, the implicit integration method was used to calculate the nonlinear viscoplastic flexibility matrix of a large geotechnical underground project. The basic idea was to determine the strain increment generated in a time increment interval by assuming the displacements and nodal stresses in the frozen wall at the moment of time :where

According to the theory of geotechnical rheological mechanics, the viscoplastic strain rate can be expressed as follows:wherewhere , , and are the coefficient matrixes of the following constitutive model:

Substitution of the above formulas into (9) yields

Since , , and coefficient matrixes are symmetrical matrixes, the transposed matrix of and coefficients in are the same, as shown in the following formula:

Therefore, the first item in (9) can be expressed as follows:

Similarly, since , , and coefficient matrixes are symmetrical matrixes, the second item in (9) can be expressed as follows:

The specific expressions of the nonlinear viscoplastic element flexibility matrix of deep artificial frozen soil under unloading stress paths are given below. It is assumed that

Therefore, the specific nonlinear viscoplastic element flexibility matrix of the artificial frozen soil under unloading stress paths can be expressed in 3D space as follows:

From (18), the nonlinear viscoplastic element flexibility matrix contains three parameters for frozen soil, that is, , which can be obtained through the triaxial creep test.

##### 2.4. Element Flexibility Matrix of Artificial Frozen Soil Considering Creep Damage

For deep shaft freezing engineering, the strength and creep characteristics of artificial frozen soil are closely related to the internal damage of frozen wall. Due to the rheology of artificial frozen soil in the deep shaft freezing method, microcracks are produced inside the frozen wall and can develop continuously with time and gradually reduce the elastic modulus, Poisson’s ratio, and the strength of the frozen wall. The artificial viscoplastic creep exhibits obvious nonlinearity, which is reflected by the decrease of the viscoplastic viscosity coefficient. The decrease of the plastic viscosity coefficient with time is mainly caused by microelement damage to the artificial frozen soil. Therefore, the frozen soil microelement strength criterion () is assumed to conform to the Weibull-type probability distribution, and the probability density function can be expressed as follows:where is the Weibull-type probability distribution variable for microelement damage in artificial frozen soil and m and are the Weibull probability distribution parameters, as determined by the creep test for artificial frozen soil.

According to the theory of probability theory, the damage probability of artificial frozen soil, represented by the damage variable can be expressed as follows:

The first stress invariable and second deviatoric stress invariable that are expressed by effective stress and overstress definitions are as follows:

Therefore, the element flexibility matrix of viscoplastic damage coupling of artificial frozen soil under unloading stress paths is as follows:where

When the artificial frozen soil conforms to the improved Zienkiewicz–Pande yield function within the main stress space, the viscoplastic creep begins. As the time of creep increases, microcracks in the frozen soil begin to expand and damage to the frozen soil appears. The damage gradually increases at a growing rate with the increase of deviatoric stress. The frozen soil structure is destroyed when many microcracks are formed, and the artificial frozen soil is completely damaged.

#### 3. Finite Element Method to Solve the Problem of Artificial Frozen Soil for Viscoelastic Damage

The basic equation of the finite element method is the general equilibrium differential equation. At time , this equation is expressed as follows:where is the creep damaged matrix of artificial frozen soil at time . The artificial frozen soil is a kind of viscoelastic-plastic damage material, so the creep matrix is related to the stress state; is the total stress in artificial frozen soil at time ; and is the equivalent nodal load of artificial frozen soil at the time of .

Formula (25) is expressed in the incremental form aswhere is the total stress increment in artificial frozen soil during the period of and is the equivalent nodal load change during the period of .

Calculation of the stress during the period of with creep increment can be done using

The Euler integral method can then be applied:

The Taylor series expansion can then be done for viscoplastic creep damage of artificial frozen soil, to obtain

Substitution of (28) into (29) gives the arithmetic expression of the strain increment for viscoplastic creep damage as follows:

Substitution of (30) into (27) yields:

Therefore, in the incremental form, the equilibrium equation of the artificial frozen soil can be expressed as

Therefore, the expression for stress increment of artificial frozen soil is as follows:

The expression for strain increment of viscoplastic damage in artificial frozen soil is as follows:

Formula (23) shows the flexibility matrix of finite element calculation that uses the viscoelastic-plastic damage coupled constitutive model of artificial frozen soil under unloading stress paths. Numerical calculation can be performed using the nonlinear finite element program ADINA, operating in the Visual Fortran 6.6A secondary development environment to generate the file of the dynamic link library [14]. The viscoelastic-plastic damage constitutive model of artificial frozen soil presented herein is written into an identifiable constitutive program in Fortran language by using the incremental finite element format in (34).

#### 4. Determination and Verification of Parameters of the Constitutive Model of Artificial Frozen Soil

##### 4.1. Determination of Parameters

The coupled constitutive model for creep damage presented herein contains eight parameters: .

###### 4.1.1. Determination of the Elastic Parameters

Based on the shear tests performed at different negative temperatures and axial pressures under the triaxial unloading stress paths, the curve was fitted on the coordinate plane of the first stress invariable and the second stress variable of deviatoric stress, to obtain the strength parameters of artificial frozen soil at different temperatures.

represent the instantaneous Poisson’s ratio and the modulus of elasticity of the artificial frozen soil, respectively. The method to determine the parameters is determined consistent with the determined parameters method of the Duncan–Chang model, a hyperbolic constitutive model considering the nonlinearity of instantaneous deformation [12, 13]. Another important elastic parameter of artificial frozen soil is Poisson’s ratio, which is determined by the ratio of axial strain and lateral strain in triaxial shear tests.

###### 4.1.2. Determination of the Viscoelastic Parameters

Based on the uniaxial or triaxial shear creep test under unloading stress paths, the relationship between axial viscoelastic creep and time was obtained. The nonlinear fitting of the curve with the Gauss–Newton method according to the principle of least squares was used to obtain three viscoelastic parameters of artificial frozen soil: .

###### 4.1.3. Determination of the Viscoplastic Damage Parameters

Based on the uniaxial or triaxial shear creep test under unloading stress path, the relationship between axial viscoplastic creep and time was obtained. The nonlinear fitting of the curve with the Gauss–Newton method according to the principle of least squares was used to obtain two viscoplastic parameters of artificial frozen soil: .

The creep test for frozen soil samples in the air shaft at Huainan coal mine was performed under triaxial unloading conditions, and the following parameters are obtained for the constitutive model of artificial frozen soil herein: .

##### 4.2. Verification of the Triaxial Test for the Constitutive Model of Artificial Frozen Soil

The creep deformation law of artificial frozen clay under test conditions was obtained through numerical simulation, with the axial stress-time relation curves, shown in Figure 2.

**(a)**

**(b)**

**(c)**

**(d)**

#### 5. Conclusions

The creep damaged constitutive program for viscoelastic-plastic damage of artificial frozen soil under unloading stress paths was developed into a large nonlinear finite element program using ADINA. The following conclusions can be obtained:(1)The triaxial shear and creep deformation laws of deep artificial frozen soil were obtained by simulating stress paths after freezing and unloading to establish the constitutive equation of artificial frozen soil under complex stress paths.(2)The value from the triaxial creep test of deep artificial frozen soil under unloading stress paths was consistent with the simulated value, with a maximum error of 4.8%. A constitutive model can reflect the primary creep, steady creep, or the accelerating creep of artificial frozen soil was proposed. These results demonstrate that it is reasonable to describe the nonlinear creep of artificial frozen soil with the viscoelastic-plastic damage coupled model using the improved Zienkiewicz–Pande parabola yield criterion.(3)A combination of basic rheological elements and yield surface models was used to describe the creep damaged deformation law of artificial frozen soil. This approach allowed development of a viscoelastic-plastic damage coupled constitutive model of deep artificial frozen soil under unloading stress paths. The calculation parameters can be obtained through a simple test, and the model parameters are few and easy to be determined. Overall, this model can be easily and widely applied during artificial freezing engineering.(4)The constitutive model provides an effective research approach to study the mechanical behavior of field frozen soil structure and presents theoretical guiding significance for the long-term stability analysis and the prediction of artificial frozen soil structure.

#### Symbols

: | The total strain |

: | The nonlinear viscoelastic strain component |

: | The nonlinear viscoplastic damage strain component |

: | The viscoelastic modulus (MPa) |

: | The transient modulus (MPa) |

: | The viscosity coefficients (MPa·h) |

: | The viscoelastic flexibility matrix |

: | The equivalent Poisson’s ratio and nonlinear viscoelastic modulus |

: | The parabolic yield function of artificial frozen soil |

: | The internal friction angle (°) |

: | The cohesive force (MPa) |

, , : | Coefficient matrixes of the constitutive model |

: | The Weibull-type probability distribution variable for microelement damage |

: | Weibull probability distribution parameters |

: | The creep damaged matrix of artificial frozen soil at time |

: | The total stress increment in artificial frozen soil during the period of |

: | The equivalent nodal load change during the period of . |

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

The authors gratefully acknowledge the support of the National Natural Science Foundation of China (Grant nos. 41672278, 41271071, and 51504070), the Fujian Provincial Natural Science Foundation Projects (Grant no. 2017Y4001), the National Basic Research Program of China (973 Program) (Grant no. 2012CB026102), and the State Key Laboratory Program of Frozen Soil Engineering (Grant no. SKLFSE201204).