# Elastic-Viscoplastic Constitutive Model of Soil under Cyclic Loading

Authors: Yuan Ma; Hua Liang; Yongxue You; Weiguang Zhang; Limin Guo; Jianwei Fan; Zhihui Cao
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

The creep problems are often involved in soil under cyclic loading, and its behaviors of soil under cyclic loading sparks many arguments in current research field. To propose one new model to demonstrate these creep behaviors of soil under cyclic loading, the cyclic loading was simplified equivalently, and the elastic-viscoplastic model (EVPM) for soil under cyclic loading was established based on the Bingham model. The yield criterion for soil under cyclic loading with constant amplitude was proposed following the simplified load. A constitutive equation based on the EVPM was constructed by using the flow criterion related to the yield criterion. Meanwhile, the parameters of EVPM were identified and discussed. In addition, the case analysis of the EVPM was also performed. The results indicate that the stable and destructive creep behaviors of soil under cyclic loading could be well described by the recommended EVPM, and the obtained parameters in the model exhibited a clear regularity with the increase of dynamic stress amplitude. Besides, the established model could be selected to predict the stable and destructive creep behavior of soil under cyclic loading.


---

## Body

## 1. Introduction

In geotechnical engineering, the building foundations often involves creep behaviors of soil under cyclic loading, such as earth dams under the conditions of cyclical changes in water level of reservoir and subgrade soil under traffic load and tunnels under high-speed railway load. The strength and creep characteristics of soil under cyclic loading are critical to the safety of engineering [1]. Under high-cycle stress, with the increase of loading time, the creep rate of soil becomes bigger; thus, the excessive deformation is caused. Meanwhile, the creep behavior will also occur under low-cycle stress. Therefore, the research concerning the creep characteristics of soil under cyclic loading has certain theoretical significance for the long-term stability evaluation as well as design of geotechnical engineering.The creep behaviors of soil under cyclic loading are divided into stable, critical, and destructive creep (Figure1) in accordance with the εN−N curves of soil under cyclic loading [2–9]. When the peak deviatoric stress σ1max−σ3 of cyclic loading is less than the critical dynamic stress σcr of soil, the creep rate soil is faster at the initial stage of loading; as the loading time is prolonged, the strain does not increase but stabilizes; in this case, the creep behavior of soil belongs to stable creep. When the peak deviatoric stress σ1max−σ3 of cyclic loading is near the critical stress value σcr of soil, the strain of soil sometimes increases and sometimes decreases with the increase of cyclic times N. In this case, the creep behavior of soil belongs to critical creep. When σ1max−σ3>σcr, the creep rate of soil increases until the soil is destroyed with the increase of N; at this time, the creep behavior of soil belongs to destructive creep behavior. Sometimes, the destructive creep behaviors show three creep stages of obvious constant speed, deceleration, and acceleration (Figure 2).Figure 1
Typical creep curves of soil under cyclic loading.Figure 2
Three-stage creep curve of soil under cyclic loading.At present, many scholars have carried out a lot of research on the models which describe the strain accumulation curves of soil shown in Figures1 and 2. Many achievements have been made in these studies. There are two main aspects in the study of soil plastic cumulative model under cyclic loading, i.e., theoretical models and mathematical empirical models. The theoretical models mainly include plastic accumulation energy model [10], rheology theory model [11–14], modified Cambridge model [15], and HCA model [16]. Mathematical empirical models include stress-strain hyperbolic model [17], exponential model [18], modified hyperbolic model on the basis of hyperbolic model and exponential model [19], grey theory model [20], and semiempirical and semitheoretical model [21].Although many constitutive models for soil under cyclic loading have been developed, many existing models cannot well reflect the plastic creep characteristics of soil under cyclic loading. Therefore, the relevant research still needs to be strengthened. In this paper, according to the existing theory, the cyclic loading was simplified equivalently. A dynamic yield criterion of soil was proposed based on the simplified load. A new constitutive equation based on the improved Bingham model was established by using the flow criterion related to the dynamic yield criterion. Meanwhile, the creep constitutive equation was obtained. The parameters of the model were identified and application analysis was also conducted on the proposed model.


## 2. Creep Theory of Soil under Cyclic Loading

### 2.1. Cyclic Loading and Simplified Load

In engineering, the most frequently encountered types of cyclic loading are sine-wave, cosine-wave, and rectangle- and triangular-wave cyclic loading. However, the rectangular- and triangular-wave cyclic loading can be transformed into sine-wave or cosine-wave cyclic loading by Fourier transformation. Therefore, in this study, the sine-wave cyclic loading equation (1) is taken as an example to study the creep behavior of soil under cyclic loading:(1)σt=σav+σd2sinωt,where σav is the vibration center of cyclic loading, that is, the average value of cyclic loading, σd is the stress double amplitude of cyclic loading, ω is the angular frequency of cyclic loading, ω=2πf (f is loading frequency), and t is the time of loading.For the sake of simplicity of the study [14], the cyclic loading (see equation (1)) acting on the soil was simplified into equation (2). However, the confining pressure was not considered in the simplified load:(2)σeq=σavexpσmax−σcrσcrf,where σmax is the upper limit of cyclic stress and σcr is the critical dynamic stress of soil.In accordance with equation (2), when σmax>σcr, the simplified load σeq is larger than the average value σav of cyclic loading; in this case, the destructive creep of soil will appear. When σmax<σcr and σeq<σav, in this case, the stable creep of soil will appear. For equation (2), the equivalent stress σeq increases with increasing frequency f and σmax of cyclic loading acting on soil. Although equation (2) was treated as a static form load on the surface, the properties of cyclic loading in essence can be reflected by it.When we consider the confining pressure, the problem is studied under the state of principal stress. It is assumed that the direction of cyclic loading is consistent with the direction of the maximal principal stressσ1; then, the cyclic loading under confining pressure can be written as(3)σ1t=pc+σst+σd2sinωt,where pc is confining pressure and σst is static deviatoric stress.Based on the idea of equation (2), when we take the influence of confining pressure on creep behavior into consideration, the deviatoric stress of soil subjected to the cyclic loading shown in equation (3) can be equivalently simplified as(4)σ1eq−σ3=σavexpσ1max−σ3−σcrσcrgf,where σav is the average value of the deviatoric stress of cyclic loading, namely, the vibration center value of the load σ1t−σ3, σ1max is the upper limit of cyclic loading, σ1max=pc+σst+σd/2, σcr is the critical strength of soil under confining pressure, only the destructive creep behavior appears, will the maximum stress value determined by the σ1t−σ3−εN curve of soil under cyclic loading occur, σ3 is the confining pressure acing on soil, and gf is a function of frequency.In equation (4), 0<σ1max−σ3−σcr/σcr<1; if equation (4) can be used to reflect the characteristic of cyclic loading, then there must be 0≤gf<1. If the gf value is too large, the simplified stress will exceed σ1max; as a result, the simplified stress is not consistent with the actual situation. For this reason, the frequency function in equation (4) is determined by(5)gf=f10n,where n is an integer and its value is determined according to the following procedure.With relation0≤gf<1, if cyclic loading is with different frequencies in the experiment when testing, the parameter n should satisfy(6)0≤fmax10n<1,where fmax is the maximum frequency value of cyclic loading used during testing, and it is determined by(7)fmax=maxf1,f2,…,fn,where f1,f2,…,fn are the frequency values of cyclic load used during the test.From equation (6), it is well known that n≥0 and n≥lgfmax. When lgfmax<0, n value is 0. When lgfmax>0, n rounds up to the nearest multiple of lgfmax.According to equation (4), it is obvious that when the peak deviatoric stress of cyclic loading σ1max−σ3>σcr is larger than σcr, the equivalent deviatoric stress σ1eq−σ3>σav, the destructive creep of soil will appear, and when σ1max−σ3<σcr, the equivalent deviatoric stress σ1eq−σ3<σav, and the destructive creep of soil will not appear.


### 2.2. Dynamic Yield Criterion for Soil under Cyclic Loading

#### 2.2.1. Establishment of Dynamic Yield Criterion

Because of the structural behavior and porosity of soil, the plastic volume and shear strain will occur when soil is subjected to cyclic loading. The creep behavior of soil under cyclic loading is related not only to its stress state but also to the characteristics of cyclic loading. Consequently, the dynamic yield function of soil under cyclic loading can be established in accordance with the equivalent load as follows:(8)FI1eq,J2eq,ϕcq,ϕeq,σav=αeqI1eq+2J2eq−Keq−13σav,where J2eq is the second invariant of equivalent stress deviator, J2eq = 1/6σ1eq−σ22+σ2−σ32+σ3−σ1eq2, I1eq is first invariant of equivalent stress tensor, αeqI1eq reflects the volumetric strain under cyclic loading, that is, the expansibility of soil, αeq is the equivalent influence coefficient of volume strain of soil, and Keq is the equivalent strength eigenvalue of soil under cyclic loading.The influence of hydrostatic pressure on the plastic volumetric strain of geomaterials is considered byαI1eq of this criterion, and the influence of stress deviator on the shear strain of soil is reflected by J2.The parametersαeq and Keq in equation (8) is taken as the equivalent form of α and K of Druker–Prager criterion [22], i.e.,(9)αeq=sinϕeq33+sin2ϕeq,Keq=3ceqcosϕeq3+sin2ϕeq,where ϕeq is equivalent internal friction angle of soil under cyclic loading and ceq is the equivalent cohesion of soil under cyclic loading.


#### 2.2.2. The Definition and Meaning ofφeq and ceq

φeq and ceq are two dynamic indicators of soil under cyclic loading, that is, two strength indexes for soil resistance to destruction, which can be calculated in accordance with the Mohr–Coulomb criterion. The calculation process is as follows. Calculate the equivalent stress σ1eq of soil in the critical state according to equation (4) and then equivalent the Mohr–Coulomb diagram plotted, as shown in Figure 3. Finally, the φeq and ceq values can be obtained from Figure 3 and also be solved by equation (10) according to the test data under two different stress states:(10)sinϕeq=σ1eq−σ3σ1eq+σ3+2ceqcotϕeq.Figure 3
Equivalent Mohr–Coulomb diagram.


## 2.1. Cyclic Loading and Simplified Load

In engineering, the most frequently encountered types of cyclic loading are sine-wave, cosine-wave, and rectangle- and triangular-wave cyclic loading. However, the rectangular- and triangular-wave cyclic loading can be transformed into sine-wave or cosine-wave cyclic loading by Fourier transformation. Therefore, in this study, the sine-wave cyclic loading equation (1) is taken as an example to study the creep behavior of soil under cyclic loading:(1)σt=σav+σd2sinωt,where σav is the vibration center of cyclic loading, that is, the average value of cyclic loading, σd is the stress double amplitude of cyclic loading, ω is the angular frequency of cyclic loading, ω=2πf (f is loading frequency), and t is the time of loading.For the sake of simplicity of the study [14], the cyclic loading (see equation (1)) acting on the soil was simplified into equation (2). However, the confining pressure was not considered in the simplified load:(2)σeq=σavexpσmax−σcrσcrf,where σmax is the upper limit of cyclic stress and σcr is the critical dynamic stress of soil.In accordance with equation (2), when σmax>σcr, the simplified load σeq is larger than the average value σav of cyclic loading; in this case, the destructive creep of soil will appear. When σmax<σcr and σeq<σav, in this case, the stable creep of soil will appear. For equation (2), the equivalent stress σeq increases with increasing frequency f and σmax of cyclic loading acting on soil. Although equation (2) was treated as a static form load on the surface, the properties of cyclic loading in essence can be reflected by it.When we consider the confining pressure, the problem is studied under the state of principal stress. It is assumed that the direction of cyclic loading is consistent with the direction of the maximal principal stressσ1; then, the cyclic loading under confining pressure can be written as(3)σ1t=pc+σst+σd2sinωt,where pc is confining pressure and σst is static deviatoric stress.Based on the idea of equation (2), when we take the influence of confining pressure on creep behavior into consideration, the deviatoric stress of soil subjected to the cyclic loading shown in equation (3) can be equivalently simplified as(4)σ1eq−σ3=σavexpσ1max−σ3−σcrσcrgf,where σav is the average value of the deviatoric stress of cyclic loading, namely, the vibration center value of the load σ1t−σ3, σ1max is the upper limit of cyclic loading, σ1max=pc+σst+σd/2, σcr is the critical strength of soil under confining pressure, only the destructive creep behavior appears, will the maximum stress value determined by the σ1t−σ3−εN curve of soil under cyclic loading occur, σ3 is the confining pressure acing on soil, and gf is a function of frequency.In equation (4), 0<σ1max−σ3−σcr/σcr<1; if equation (4) can be used to reflect the characteristic of cyclic loading, then there must be 0≤gf<1. If the gf value is too large, the simplified stress will exceed σ1max; as a result, the simplified stress is not consistent with the actual situation. For this reason, the frequency function in equation (4) is determined by(5)gf=f10n,where n is an integer and its value is determined according to the following procedure.With relation0≤gf<1, if cyclic loading is with different frequencies in the experiment when testing, the parameter n should satisfy(6)0≤fmax10n<1,where fmax is the maximum frequency value of cyclic loading used during testing, and it is determined by(7)fmax=maxf1,f2,…,fn,where f1,f2,…,fn are the frequency values of cyclic load used during the test.From equation (6), it is well known that n≥0 and n≥lgfmax. When lgfmax<0, n value is 0. When lgfmax>0, n rounds up to the nearest multiple of lgfmax.According to equation (4), it is obvious that when the peak deviatoric stress of cyclic loading σ1max−σ3>σcr is larger than σcr, the equivalent deviatoric stress σ1eq−σ3>σav, the destructive creep of soil will appear, and when σ1max−σ3<σcr, the equivalent deviatoric stress σ1eq−σ3<σav, and the destructive creep of soil will not appear.


## 2.2. Dynamic Yield Criterion for Soil under Cyclic Loading

### 2.2.1. Establishment of Dynamic Yield Criterion

Because of the structural behavior and porosity of soil, the plastic volume and shear strain will occur when soil is subjected to cyclic loading. The creep behavior of soil under cyclic loading is related not only to its stress state but also to the characteristics of cyclic loading. Consequently, the dynamic yield function of soil under cyclic loading can be established in accordance with the equivalent load as follows:(8)FI1eq,J2eq,ϕcq,ϕeq,σav=αeqI1eq+2J2eq−Keq−13σav,where J2eq is the second invariant of equivalent stress deviator, J2eq = 1/6σ1eq−σ22+σ2−σ32+σ3−σ1eq2, I1eq is first invariant of equivalent stress tensor, αeqI1eq reflects the volumetric strain under cyclic loading, that is, the expansibility of soil, αeq is the equivalent influence coefficient of volume strain of soil, and Keq is the equivalent strength eigenvalue of soil under cyclic loading.The influence of hydrostatic pressure on the plastic volumetric strain of geomaterials is considered byαI1eq of this criterion, and the influence of stress deviator on the shear strain of soil is reflected by J2.The parametersαeq and Keq in equation (8) is taken as the equivalent form of α and K of Druker–Prager criterion [22], i.e.,(9)αeq=sinϕeq33+sin2ϕeq,Keq=3ceqcosϕeq3+sin2ϕeq,where ϕeq is equivalent internal friction angle of soil under cyclic loading and ceq is the equivalent cohesion of soil under cyclic loading.


### 2.2.2. The Definition and Meaning ofφeq and ceq

φeq and ceq are two dynamic indicators of soil under cyclic loading, that is, two strength indexes for soil resistance to destruction, which can be calculated in accordance with the Mohr–Coulomb criterion. The calculation process is as follows. Calculate the equivalent stress σ1eq of soil in the critical state according to equation (4) and then equivalent the Mohr–Coulomb diagram plotted, as shown in Figure 3. Finally, the φeq and ceq values can be obtained from Figure 3 and also be solved by equation (10) according to the test data under two different stress states:(10)sinϕeq=σ1eq−σ3σ1eq+σ3+2ceqcotϕeq.Figure 3
Equivalent Mohr–Coulomb diagram.


## 2.2.1. Establishment of Dynamic Yield Criterion

Because of the structural behavior and porosity of soil, the plastic volume and shear strain will occur when soil is subjected to cyclic loading. The creep behavior of soil under cyclic loading is related not only to its stress state but also to the characteristics of cyclic loading. Consequently, the dynamic yield function of soil under cyclic loading can be established in accordance with the equivalent load as follows:(8)FI1eq,J2eq,ϕcq,ϕeq,σav=αeqI1eq+2J2eq−Keq−13σav,where J2eq is the second invariant of equivalent stress deviator, J2eq = 1/6σ1eq−σ22+σ2−σ32+σ3−σ1eq2, I1eq is first invariant of equivalent stress tensor, αeqI1eq reflects the volumetric strain under cyclic loading, that is, the expansibility of soil, αeq is the equivalent influence coefficient of volume strain of soil, and Keq is the equivalent strength eigenvalue of soil under cyclic loading.The influence of hydrostatic pressure on the plastic volumetric strain of geomaterials is considered byαI1eq of this criterion, and the influence of stress deviator on the shear strain of soil is reflected by J2.The parametersαeq and Keq in equation (8) is taken as the equivalent form of α and K of Druker–Prager criterion [22], i.e.,(9)αeq=sinϕeq33+sin2ϕeq,Keq=3ceqcosϕeq3+sin2ϕeq,where ϕeq is equivalent internal friction angle of soil under cyclic loading and ceq is the equivalent cohesion of soil under cyclic loading.


## 2.2.2. The Definition and Meaning ofφeq and ceq

φeq and ceq are two dynamic indicators of soil under cyclic loading, that is, two strength indexes for soil resistance to destruction, which can be calculated in accordance with the Mohr–Coulomb criterion. The calculation process is as follows. Calculate the equivalent stress σ1eq of soil in the critical state according to equation (4) and then equivalent the Mohr–Coulomb diagram plotted, as shown in Figure 3. Finally, the φeq and ceq values can be obtained from Figure 3 and also be solved by equation (10) according to the test data under two different stress states:(10)sinϕeq=σ1eq−σ3σ1eq+σ3+2ceqcotϕeq.Figure 3
Equivalent Mohr–Coulomb diagram.


## 3. Elastic-Viscoplastic Constitutive Equation for Soil under Cyclic Loading

### 3.1. Elastic-Viscoplastic Model Based on the Bingham Model

The Bingham model is composed of elastic element and ideal viscoplastic body in series, which cannot reflect the nonlinear creep behavior characteristics of geotechnical materials [23]. In order to describe the nonlinear creep behavior of soil under cyclic loading, a modified Bingham constitutive model (see Figure 4) is proposed by replacing the constant elastic and viscous elements in the Bingham model with nonlinear elements. A three-dimensional elasto-viscoplastic fatigue model (EVPM) for soil was established on the basis of the improved Bingham model. In Figure 4, GN is the elastic shear coefficient and ηN is viscoplastic shear coefficient.Figure 4
Improved Bingham fatigue model.


#### 3.1.1. Elastic Element

The relationships between dynamic strength, dynamic shear modulus, and the damping ratio of geomaterials under cyclic loading and the number of cyclesN (or loads time t) were studied by many researchers such as Lai Xilei [24] and Mu Kun [25]. These results show that the dynamic parameters of geotechnical materials decrease with the increase of N (or loading time t), and most of them decay in the form of exponential functions. For this reason, it is assumed that the elastic shear fatigue coefficient GN in the model, shown in Figure 4, is(11)GN=G0exp−λNfλ,where G0 is the initial value of shear coefficient of elastic element, and its physical dimension is stress, and λ is a dimensionless parameter, which is used to reflect the change rule of elastic shear coefficient and the strain rate of soil at the stable creep stage, and it can be determined by test data and λ>0.


#### 3.1.2. Viscous Element

Based on fractional calculus theory, the fractional element is defined and the fractional Kelvin–Voigt rheological model is constructed by Guo Jiaqi et al. [26]. It is found that the strain compliance Jt of the fractional element is proportional to the rheological time tβ, in which β is a constant and t is rheology time. Schiessel et al. [27] gave the strain compliance of the fractional unit proportional to tβ by considering the time accumulation of viscous elements in the Kelvin model and constructed a fractional rheological constitutive equation by defining fractional units according to fractional derivative. From the above, it can be assumed that the coefficient ηN of the viscous element in the model, shown in Figure 4, is(12)ηN=η0t1−n=η0Nf1−n,where n is the strain rate parameter for soil destructive creep behavior, n≠1, η0 is an initial value of viscoplastic element, and η0≠0, its physical dimension is stress·time n.


### 3.2. Constitutive Equation for Stable Creep of Soil under Cyclic Loading

When the stable creep behavior occurs, the creep behavior is described by the nonlinear elastic element in the model, as shown in Figure4. At this time, the model is degenerated to the classical nonlinear elastic model. According to the basic theory of elastic mechanics, the state equation of elastic element can be expressed as(13)Sijeq=2GNeijeN,εm=13KNσmeq,where Sijeq is the equivalent stress deviator, Sijeq=σijeq−σmeqδij (δij is Kronecker delta mark), eijeN is the strain tensor of elastic element, σmeq is equivalent average stress, εm is the volumetric strain, and KN is the fatigue bulk modulus of soil, assuming that its variation law is the same as that of GN; accordingly, we have(14)KN=K0exp−λNfλ,where K0 is the initial value of volume coefficient.When the cycle number of load isN, the strain of the elastic element is(15)εijeN=12G0Sijeq+σkkeqδij9K01exp−λN/fλ,where σkkeq is the first invariant of the stress tensor.As indicated in equation (15), we have(1)
When the cycles’ numberN⟶0 of cyclic load σ1t, ΔεijeN⟶1/2G0Sijeq+1/9K0σkkeqδij, i.e., which is the transient elastic strain.(2)
WhenN⟶∞, ΔεijN is directly proportional to the exponential function expλN/fλ. When different λ are chosen, the stable creep behavior of soil is reflected.


### 3.3. Constitutive Equation for Soil Destructive Creep Behavior

When the destructive creep of soil occurs, in the model shown in Figure4, the soil strain is composed of elastic and viscoplastic strain. When the cycle number of cyclic loading is N, soil strain is(16)εijN=εijeN+εijvpN.In accordance with the basic principle of rheological mechanics, the viscoplastic flow rate in equation (16) can be determined by [28].(17)ε˙l,vpt=1ηMφF∂Q∂σ,where ηM is viscoplastic flow coefficient and φF is calculated by equation (18). The main forms of φF are equations (19) and (20), respectively [29]:(18)φF=0,F<0,φF,F>0,(19)φF=Fb⋅b≥1,(20)φF=expF−1.The plastic strain of soil under cyclic loading can be determined by equation (17). With t=N/f and dt=dN/f, the viscoplastic flow rate of viscoplastic bodies in the model shown in Figure 4 can be written as(21)ε˙ijvpt=fε˙ijvpN=1ηNφF∂Q∂σij,where ε˙ijvpt is the viscoplastic flow rate of viscoplastic bodies in the model at time t, e˙ijvpN is the viscoplastic strain rate of the viscoplastic body in the model shown in Figure 4 when the number of cycles is N, F is the dynamic yield function of soil, and Q is the dynamic plastic potential function of soil.φF is calculated by equation (19) in this article, and taking b as 1 in processes of calculating, then equation (17) is expressed as(22)φF=F.If the flow rule related to yield criterion is adopted, thenF=Q. When the destructive creep behavior occurs, the strain of viscoplastic body in the model can be obtained by integrating equation (22) as follows:(23)εijvpN=1η0nF∂F∂σijNfn.When the destructive creep behavior occurs, it is assumed that the main deformation is plastic deformation. Without considering the volume strain and instantaneous elastic strain, the constitutive equation of elastic-viscoplastic fatigue based on the model shown in Figure4 can be obtained from equations (15) and (23) as follows:(24)ΔεijpN=12G0exp−λN/fλSijeq+1η0nF∂F∂σijNfn−12G0Sijeq.WhenN≠0 and the different parameters G0, λ, and η0 are chosen, equation (24) can be used to describe the destructive creep characteristics of soil under the condition of (σ1max−σ3) > σer. In addition, equation (24) can also be used to describe the three creep stages of soil under cyclic loading.


## 3.1. Elastic-Viscoplastic Model Based on the Bingham Model

The Bingham model is composed of elastic element and ideal viscoplastic body in series, which cannot reflect the nonlinear creep behavior characteristics of geotechnical materials [23]. In order to describe the nonlinear creep behavior of soil under cyclic loading, a modified Bingham constitutive model (see Figure 4) is proposed by replacing the constant elastic and viscous elements in the Bingham model with nonlinear elements. A three-dimensional elasto-viscoplastic fatigue model (EVPM) for soil was established on the basis of the improved Bingham model. In Figure 4, GN is the elastic shear coefficient and ηN is viscoplastic shear coefficient.Figure 4
Improved Bingham fatigue model.


### 3.1.1. Elastic Element

The relationships between dynamic strength, dynamic shear modulus, and the damping ratio of geomaterials under cyclic loading and the number of cyclesN (or loads time t) were studied by many researchers such as Lai Xilei [24] and Mu Kun [25]. These results show that the dynamic parameters of geotechnical materials decrease with the increase of N (or loading time t), and most of them decay in the form of exponential functions. For this reason, it is assumed that the elastic shear fatigue coefficient GN in the model, shown in Figure 4, is(11)GN=G0exp−λNfλ,where G0 is the initial value of shear coefficient of elastic element, and its physical dimension is stress, and λ is a dimensionless parameter, which is used to reflect the change rule of elastic shear coefficient and the strain rate of soil at the stable creep stage, and it can be determined by test data and λ>0.


### 3.1.2. Viscous Element

Based on fractional calculus theory, the fractional element is defined and the fractional Kelvin–Voigt rheological model is constructed by Guo Jiaqi et al. [26]. It is found that the strain compliance Jt of the fractional element is proportional to the rheological time tβ, in which β is a constant and t is rheology time. Schiessel et al. [27] gave the strain compliance of the fractional unit proportional to tβ by considering the time accumulation of viscous elements in the Kelvin model and constructed a fractional rheological constitutive equation by defining fractional units according to fractional derivative. From the above, it can be assumed that the coefficient ηN of the viscous element in the model, shown in Figure 4, is(12)ηN=η0t1−n=η0Nf1−n,where n is the strain rate parameter for soil destructive creep behavior, n≠1, η0 is an initial value of viscoplastic element, and η0≠0, its physical dimension is stress·time n.


## 3.1.1. Elastic Element

The relationships between dynamic strength, dynamic shear modulus, and the damping ratio of geomaterials under cyclic loading and the number of cyclesN (or loads time t) were studied by many researchers such as Lai Xilei [24] and Mu Kun [25]. These results show that the dynamic parameters of geotechnical materials decrease with the increase of N (or loading time t), and most of them decay in the form of exponential functions. For this reason, it is assumed that the elastic shear fatigue coefficient GN in the model, shown in Figure 4, is(11)GN=G0exp−λNfλ,where G0 is the initial value of shear coefficient of elastic element, and its physical dimension is stress, and λ is a dimensionless parameter, which is used to reflect the change rule of elastic shear coefficient and the strain rate of soil at the stable creep stage, and it can be determined by test data and λ>0.


## 3.1.2. Viscous Element

Based on fractional calculus theory, the fractional element is defined and the fractional Kelvin–Voigt rheological model is constructed by Guo Jiaqi et al. [26]. It is found that the strain compliance Jt of the fractional element is proportional to the rheological time tβ, in which β is a constant and t is rheology time. Schiessel et al. [27] gave the strain compliance of the fractional unit proportional to tβ by considering the time accumulation of viscous elements in the Kelvin model and constructed a fractional rheological constitutive equation by defining fractional units according to fractional derivative. From the above, it can be assumed that the coefficient ηN of the viscous element in the model, shown in Figure 4, is(12)ηN=η0t1−n=η0Nf1−n,where n is the strain rate parameter for soil destructive creep behavior, n≠1, η0 is an initial value of viscoplastic element, and η0≠0, its physical dimension is stress·time n.


## 3.2. Constitutive Equation for Stable Creep of Soil under Cyclic Loading

When the stable creep behavior occurs, the creep behavior is described by the nonlinear elastic element in the model, as shown in Figure4. At this time, the model is degenerated to the classical nonlinear elastic model. According to the basic theory of elastic mechanics, the state equation of elastic element can be expressed as(13)Sijeq=2GNeijeN,εm=13KNσmeq,where Sijeq is the equivalent stress deviator, Sijeq=σijeq−σmeqδij (δij is Kronecker delta mark), eijeN is the strain tensor of elastic element, σmeq is equivalent average stress, εm is the volumetric strain, and KN is the fatigue bulk modulus of soil, assuming that its variation law is the same as that of GN; accordingly, we have(14)KN=K0exp−λNfλ,where K0 is the initial value of volume coefficient.When the cycle number of load isN, the strain of the elastic element is(15)εijeN=12G0Sijeq+σkkeqδij9K01exp−λN/fλ,where σkkeq is the first invariant of the stress tensor.As indicated in equation (15), we have(1)
When the cycles’ numberN⟶0 of cyclic load σ1t, ΔεijeN⟶1/2G0Sijeq+1/9K0σkkeqδij, i.e., which is the transient elastic strain.(2)
WhenN⟶∞, ΔεijN is directly proportional to the exponential function expλN/fλ. When different λ are chosen, the stable creep behavior of soil is reflected.


## 3.3. Constitutive Equation for Soil Destructive Creep Behavior

When the destructive creep of soil occurs, in the model shown in Figure4, the soil strain is composed of elastic and viscoplastic strain. When the cycle number of cyclic loading is N, soil strain is(16)εijN=εijeN+εijvpN.In accordance with the basic principle of rheological mechanics, the viscoplastic flow rate in equation (16) can be determined by [28].(17)ε˙l,vpt=1ηMφF∂Q∂σ,where ηM is viscoplastic flow coefficient and φF is calculated by equation (18). The main forms of φF are equations (19) and (20), respectively [29]:(18)φF=0,F<0,φF,F>0,(19)φF=Fb⋅b≥1,(20)φF=expF−1.The plastic strain of soil under cyclic loading can be determined by equation (17). With t=N/f and dt=dN/f, the viscoplastic flow rate of viscoplastic bodies in the model shown in Figure 4 can be written as(21)ε˙ijvpt=fε˙ijvpN=1ηNφF∂Q∂σij,where ε˙ijvpt is the viscoplastic flow rate of viscoplastic bodies in the model at time t, e˙ijvpN is the viscoplastic strain rate of the viscoplastic body in the model shown in Figure 4 when the number of cycles is N, F is the dynamic yield function of soil, and Q is the dynamic plastic potential function of soil.φF is calculated by equation (19) in this article, and taking b as 1 in processes of calculating, then equation (17) is expressed as(22)φF=F.If the flow rule related to yield criterion is adopted, thenF=Q. When the destructive creep behavior occurs, the strain of viscoplastic body in the model can be obtained by integrating equation (22) as follows:(23)εijvpN=1η0nF∂F∂σijNfn.When the destructive creep behavior occurs, it is assumed that the main deformation is plastic deformation. Without considering the volume strain and instantaneous elastic strain, the constitutive equation of elastic-viscoplastic fatigue based on the model shown in Figure4 can be obtained from equations (15) and (23) as follows:(24)ΔεijpN=12G0exp−λN/fλSijeq+1η0nF∂F∂σijNfn−12G0Sijeq.WhenN≠0 and the different parameters G0, λ, and η0 are chosen, equation (24) can be used to describe the destructive creep characteristics of soil under the condition of (σ1max−σ3) > σer. In addition, equation (24) can also be used to describe the three creep stages of soil under cyclic loading.


## 4. Constitutive Equation of Soil under Triaxial Compression

The dynamic triaxial compression test is an indoor test method for studying the creep behavior of soil. The rationality of the proposed model can be verified by the dynamic triaxial test data. Therefore, the constitutive equation of soil under triaxial compression is obtained in view of equations (15) and (24).According to equation (4), under triaxial cyclic compression, the equivalent stress state of soil subjected to(25)σijeq=σavexpσ1max−σ3−σcrσcrf+σ3000σ3000σ3,where σ3 is the confining pressure acting on the specimen.According to equation (8), the dynamic yield function under triaxial condition is(26)FI1eq,J2eq,ϕcq,ϕeq,σav=αeqσ1eq+2σ3+13σ1eq−σ3−Keq+13σ1eq−σ3−σav.By plugging equations (25) and (26) into equations (15) and (24), the constitutive equations of soil under triaxial compression are obtained. So, the constitutive equations of soil under triaxial compression are shown in equations (27) and (28), respectively:(27)ε11N=σ1eq−σ33G0exp−λN/fλ+σ1eq+2σ39K0exp−λN/fλ,(28)Δε11pN=σ1eq−σ33G0exp−λN/fλ+α+232σ1eq+2αeq2+23αeq−43σ3−Keq+13σavαeq+231η0nNfn−σ1eq−σ33G0.


## 5. Adaptability Verification and Example Analysis of EVPM

### 5.1. Parameter Determination and Adaptability Verification of EVPFM

#### 5.1.1. Determination of Parameters

1stop is leading software for nonlinear curve fitting and comprehensive optimization analysis in the world. It has been widely recognized in the field of nonlinear regression, fitting nonlinear curve, and parameter estimation of complex engineering models. In this paper, the software was used to solve the parameters of the proposed model.In equation (27), there are 3 parameters: G0, K0, and λ. In equation (28), there are 4 parameters: G0, λ, η0, and n. In order to verify the adaptability of the model, the PO (Powell Optimization) method and SA (Simulated Annealing) method in 1stop software are used to calculate the parameters of the constitutive equation for destructive and stable creep of soil under cyclic loading, respectively.The process of solving model parameters and fitting test data is as follows. According to equation (4), the equivalent stress under different stress amplitudes was obtained. When the stable creep behavior occurs, equation (27) and the corresponding experimental data were input into the software. When the destruction creep behavior occurs, the constitutive equation (28) and the corresponding test data were input, and then, the corresponding optimization algorithm for parameter calculation was selected. Finally, the curves and parameters of the proposed model could be obtained by running the program.


#### 5.1.2. Calculation of EVPM Parameters for Laterite Mixed with Cement

Kong et al. [30] used the British GDS dynamic triaxial test system to carry out cyclic test on laterite mixed with cement under sine-wave periodic load with frequency f=5Hz; the stress level is D=σd/2τf (τf is the static shear strength) when testing. According to the test, the critical dynamic stress amplitudes σdc of laterite mixed with cement under confining pressure pc of 25, 50, and 100 kPa are 216, 274, and 339 kPa respectively; the static shear strength τf is 720 kPa, 917 kPa, and 1130 kPa, respectively. The test conditions are shown in Table 1. The vibration center of the axial dynamic stress is σ=pc+σd, and the axial dynamic stress is(29)σ1t=σ3+σd+σd2sinϖt.Table 1
Test conditions [30].
Confining pressure (kPa)Dynamic stress amplitudeσd (kPa)25254565951502252452605050100150175200285310—100100200300330360370380390According to equation (29), it is well known that σ1max=σ3+σd+σd/2, and the average value of the deviatoric stress under cyclic loading is σav=σd.Under 25, 50, and 100 kPa confining pressure, when the critical state occurs, the average value of cyclic loading is 216, 274, and 339 kPa; then, the critical strengthσcr=σdc+σdc/2 are 314, 411, and 508.5 kPa, respectively. The corresponding σ1max at critical state are 457, 598, and 778 kPa, respectively.n in equation (4) takes 1 here. The critical equivalent stress σ1eq of soil under 25, 50, and 100 kPa confining pressure is 243.54, 324.78, and 401.82 kPa, respectively.According to equation (10), the equivalent internal friction angle ϕeq=21.10∘ and the equivalent cohesion ceq=78.42kPa of laterite mixed with cement can be obtained by selecting the equivalent stress state under 25 and 100 kPa confining pressure. αeq=0.12 and Keq=76.73kPa can be obtained from equation (9).According to the obtained parameters, the constitutive equation (28) of laterite mixed with cement can be further simplified to(30)Δε11pN=σ1eq−σ33G0exp−λN/fλ+σ1eq+2σ39K0exp−λN/fλ+1.62σ1eq−1.17σ3−1.2776.73+13σav1η0nNfn.According to the data from Table1 and equation (4), the equivalent deviator stress σ1eq−σ3 under 25 and 100 kPa confining pressure can be calculated; the results of equivalent deviator stresses under different dynamic stress amplitudes σd are shown in Table 2.Table 2
The equivalent deviator stress of cyclic loading.
Confining pressure (kPa)Equivalentσ1eq−σ3 (kPa)2516.1030.4745.8172.52130.14234.18266.74293.1510070.46162.12282.53326.72371.33387.42403.50422.48It can be seen from Table2 that if σ1max−σ3>σcr of cyclic loading, the equivalent deviatoric stress σ1eq−σ3>σav obtained from equation (2). On the contrary, if σ1max−σ3<σcr, σ1eq−σ3<σav. In addition, the equivalent stress σ1eq−σ3 obtained was near σav of cyclic loading, indicating that the equivalent of equation (4) is reasonable and the effect of cyclic loading on the soil can be reflected by it.In this research, the obtained parameters of the model of soil under the confining pressure of 100 kPa are shown in Table3, and the comparison between model curves and test curves is shown in Figure 5.Table 3
EVPM parameters of laterite mixed with cement under 100 kPa confining pressure [30].
σd (kPa)G0 (kPa)K0 (kPa)λη0 (kPa·sn)nR1007.113 × 10153.969 × 1040.1840.8702005.735 × 10152.918 × 1040.1760.8923005.472 × 10153.000 × 1040.1790.8933304.960 × 10152.510 × 1040.1720.8663603.348 × 1030.1972.210 × 107−0.0030.9993704.158 × 1030.1856.793 × 1061.3580.9983803.702 × 1030.1571.917 × 1061.3400.9963902.679 × 1030.1166.035 × 1051.3570.999Figure 5
Fitting the curves of laterite mixed with cement [30] under the confining pressure of 100 kPa by EVPFM.From Figure5 and Table 3, it can be seen that the elastic-viscoplastic constitutive model based on the Bingham model can describe the stable and destructive creep behavior of soil under cyclic loading. The fitting precision of various deformation curves is good. The fitting correlation coefficient is above 0.870 and up to 0.999 under high and low dynamic stress. Therefore, the constitutive equation established in this paper can be used to reflect the various creep behavior of soil under cyclic loading.


#### 5.1.3. Analysis of EVPM Parameters of Laterite Mixed with Cement

Under the confining pressure of 100 kPa, the obtained model curve is shown in Figure5, and the fitted parameters are shown in Table 3. By analyzing the model parameters in Table 3, it can be found that when stable creep behavior occurs, the model parameters λ, K0, and G0 decayed with the increase of σ1eq−σ3 (Figures 6–8) and when destructive creep behavior occurs, λ, η0 and G0 decreased with the increase of σ1eq−σ3 (Figures 9–11). When the upper limit of the cyclic loading is near the critical dynamic stress, n was close to 0, and when the cyclic load increased, n was close to a fixed value. R in Table 3 is correlation coefficient.Figure 6
Relationship betweenG0 and σ1eq−σ3 for stable creep.Figure 7
Relationship betweenK0 and σ1eq−σ3 for stable creep.Figure 8
Relationship betweenλ and σ1eq−σ3 for stable creep.Figure 9
Relationship betweenλ and σ1eq−σ3 for destructive creep.Figure 10
Relationship betweenη0 and σ1eq−σ3 for destructive creep.Figure 11
Relationship betweenG0 and σ1eq−σ3 for destructive creep.


### 5.2. Application Analysis

From the variation of model parameters under different dynamic stress amplitudes [31–33], when the confining pressure is 100 kPa, it is known that the value of each parameter varies obviously with σ1eq−σ3. According to the above analysis, it can be assumed that the parameter variation laws under any dynamic pressure and different dynamic stress are consistent with 100 kPa, but the numerical values are inconsistent. Under different confining pressures, the values of each parameter can be determined by a set of tests under specific confining pressures. Consequently, it can be assumed that the model parameters under a certain confining pressure with different dynamic stress amplitudes are determined by equations (31)–(33) when a stable creep behavior occurs:(31)K0=1.17×105σ1eq−σ3−0.259×β1sσ1eq,σ2,σ3,(32)G0=1.770×1016×σ1eq−σ3−0.216×β2sσ1eq,σ2,σ3,(33)λ=0.210×σ1eq−σ3−0.033×β3sσ1eq,σ2,σ3.In equations (31)–(33), β1s,β2s,andβ3s values can be obtained by plugging parameters K0, G0, and λ determined by a group test data from stable creep into equations (31)–(33). Of course, in order to get more accurate values, it can be determined by combining more test data.When destructive creep occurs, the parametersG0, λ, and η0 are determined by the following equations (34)–(36):(34)G0=1.425×1016×σ1eq−σ3−4.840×β1fσ1eq,σ2,σ3,(35)λ=1.159×σ1eq−σ3−3.795×β2fσ1eq,σ2,σ3,(36)η=8.285×1076σ1eq−σ3−28.241×β3fσ1eq,σ2,σ3In equations (34)–(36), β1f,β2f,andβ3f values can be obtained by plugging parameters G0, λ, and η0 determined by a group test data from destructive creep into equations (31)–(33). Of course, in order to get more accurate values, it can be determined by combining a set of test data.When the upper limit of stress is close to the critical dynamic stress, the parametern value tends to 0. Under the other different conditions of σd, it is temporarily determined to be 1.350 in combination with Table 3, but the specific change law of n should be further studied in combination with a large number of experiments.In this paper, when stable creep occurs,β1s=1/3.44, β2s=1/1.92, and β3s=1/0.946 are obtained by plugging the parameter K0, G0, and λ value under the dynamic stress amplitude of 25 kPa and 25 kPa confining pressure into equations (31)–(33). When calculating, because the part result of calculation is slightly smaller, s is taken as 1/4. Finally, the parameters K0, G0, and λ under dynamic stress amplitude of 45, 65, 95, 150, and 225 kPa are calculated in Table 4. The deformation under dynamic stress amplitudes 45, 65, 95, 150, and 225 kPa is predicted by using the calculated parameters. The predicted results are compared with the experimental results, as shown in Figure 12, when stable creep behavior occurs.Table 4
EVPM parameters of laterite mixed with cement under 25 kPa confining pressure [30].
σd (kPa)G0 (kPa)K0 (kPa)λη0 (kPa·sn)nRemark255.045 × 10151.665 × 1040.204Test data454.406 × 10151.413 × 1040.200Calculation654.036 × 10151.270 × 1040.197953.656 × 10151127 × 1040.1931503.219 × 10159.709 × 1030.1902252.836 × 10158313 × 1030.18624532540.2169.160 × 106-0.009Test data26020570.1526.305 × 1051.350CalculationFigure 12
Comparison between model predictive value and test data of laterite mixed with cement [30].When the destructive creep behavior occurs, in this paper, the parameter value ofβ1f,β2f,andβ3f can be obtained by combining the experimental parameter value under confining pressure 25 kPa and dynamic stress amplitude 245 kPa with equations (32)–(34). β1f=0.126, β2f=0.302, and β3f=6.476×10−5 are obtained by plugging the parameter K0, η0, and λ value under the dynamic stress amplitude of 25 kPa and 25 kPa confining pressure into equations (32)–(34). According to these obtained parameters, the creep of soil under cyclic loading with a stress amplitude of 260 kPa is calculated, and the predicted results are compared with the experimental results, as shown in Figure 13.Figure 13
Comparison between model predictive value and test data of laterite mixed with cement [30] when destructive deformation occurs.The calculation results of soil deformation under confining pressure 25 kPa show that the proposed EVPM can predict the deformation of soil under cyclic loading. There is a certain error between the calculated value of the model and the experimental value. The main reason for the error is that the experimental data are few, which leads to the inaccuracy of the changing rules of the parameters and the calculation error. Further research is needed in combination with experiments. However, the calculated results of the model are roughly consistent with the experimental data, which further illustrates the feasibility of the model.


## 5.1. Parameter Determination and Adaptability Verification of EVPFM

### 5.1.1. Determination of Parameters

1stop is leading software for nonlinear curve fitting and comprehensive optimization analysis in the world. It has been widely recognized in the field of nonlinear regression, fitting nonlinear curve, and parameter estimation of complex engineering models. In this paper, the software was used to solve the parameters of the proposed model.In equation (27), there are 3 parameters: G0, K0, and λ. In equation (28), there are 4 parameters: G0, λ, η0, and n. In order to verify the adaptability of the model, the PO (Powell Optimization) method and SA (Simulated Annealing) method in 1stop software are used to calculate the parameters of the constitutive equation for destructive and stable creep of soil under cyclic loading, respectively.The process of solving model parameters and fitting test data is as follows. According to equation (4), the equivalent stress under different stress amplitudes was obtained. When the stable creep behavior occurs, equation (27) and the corresponding experimental data were input into the software. When the destruction creep behavior occurs, the constitutive equation (28) and the corresponding test data were input, and then, the corresponding optimization algorithm for parameter calculation was selected. Finally, the curves and parameters of the proposed model could be obtained by running the program.


### 5.1.2. Calculation of EVPM Parameters for Laterite Mixed with Cement

Kong et al. [30] used the British GDS dynamic triaxial test system to carry out cyclic test on laterite mixed with cement under sine-wave periodic load with frequency f=5Hz; the stress level is D=σd/2τf (τf is the static shear strength) when testing. According to the test, the critical dynamic stress amplitudes σdc of laterite mixed with cement under confining pressure pc of 25, 50, and 100 kPa are 216, 274, and 339 kPa respectively; the static shear strength τf is 720 kPa, 917 kPa, and 1130 kPa, respectively. The test conditions are shown in Table 1. The vibration center of the axial dynamic stress is σ=pc+σd, and the axial dynamic stress is(29)σ1t=σ3+σd+σd2sinϖt.Table 1
Test conditions [30].
Confining pressure (kPa)Dynamic stress amplitudeσd (kPa)25254565951502252452605050100150175200285310—100100200300330360370380390According to equation (29), it is well known that σ1max=σ3+σd+σd/2, and the average value of the deviatoric stress under cyclic loading is σav=σd.Under 25, 50, and 100 kPa confining pressure, when the critical state occurs, the average value of cyclic loading is 216, 274, and 339 kPa; then, the critical strengthσcr=σdc+σdc/2 are 314, 411, and 508.5 kPa, respectively. The corresponding σ1max at critical state are 457, 598, and 778 kPa, respectively.n in equation (4) takes 1 here. The critical equivalent stress σ1eq of soil under 25, 50, and 100 kPa confining pressure is 243.54, 324.78, and 401.82 kPa, respectively.According to equation (10), the equivalent internal friction angle ϕeq=21.10∘ and the equivalent cohesion ceq=78.42kPa of laterite mixed with cement can be obtained by selecting the equivalent stress state under 25 and 100 kPa confining pressure. αeq=0.12 and Keq=76.73kPa can be obtained from equation (9).According to the obtained parameters, the constitutive equation (28) of laterite mixed with cement can be further simplified to(30)Δε11pN=σ1eq−σ33G0exp−λN/fλ+σ1eq+2σ39K0exp−λN/fλ+1.62σ1eq−1.17σ3−1.2776.73+13σav1η0nNfn.According to the data from Table1 and equation (4), the equivalent deviator stress σ1eq−σ3 under 25 and 100 kPa confining pressure can be calculated; the results of equivalent deviator stresses under different dynamic stress amplitudes σd are shown in Table 2.Table 2
The equivalent deviator stress of cyclic loading.
Confining pressure (kPa)Equivalentσ1eq−σ3 (kPa)2516.1030.4745.8172.52130.14234.18266.74293.1510070.46162.12282.53326.72371.33387.42403.50422.48It can be seen from Table2 that if σ1max−σ3>σcr of cyclic loading, the equivalent deviatoric stress σ1eq−σ3>σav obtained from equation (2). On the contrary, if σ1max−σ3<σcr, σ1eq−σ3<σav. In addition, the equivalent stress σ1eq−σ3 obtained was near σav of cyclic loading, indicating that the equivalent of equation (4) is reasonable and the effect of cyclic loading on the soil can be reflected by it.In this research, the obtained parameters of the model of soil under the confining pressure of 100 kPa are shown in Table3, and the comparison between model curves and test curves is shown in Figure 5.Table 3
EVPM parameters of laterite mixed with cement under 100 kPa confining pressure [30].
σd (kPa)G0 (kPa)K0 (kPa)λη0 (kPa·sn)nR1007.113 × 10153.969 × 1040.1840.8702005.735 × 10152.918 × 1040.1760.8923005.472 × 10153.000 × 1040.1790.8933304.960 × 10152.510 × 1040.1720.8663603.348 × 1030.1972.210 × 107−0.0030.9993704.158 × 1030.1856.793 × 1061.3580.9983803.702 × 1030.1571.917 × 1061.3400.9963902.679 × 1030.1166.035 × 1051.3570.999Figure 5
Fitting the curves of laterite mixed with cement [30] under the confining pressure of 100 kPa by EVPFM.From Figure5 and Table 3, it can be seen that the elastic-viscoplastic constitutive model based on the Bingham model can describe the stable and destructive creep behavior of soil under cyclic loading. The fitting precision of various deformation curves is good. The fitting correlation coefficient is above 0.870 and up to 0.999 under high and low dynamic stress. Therefore, the constitutive equation established in this paper can be used to reflect the various creep behavior of soil under cyclic loading.


### 5.1.3. Analysis of EVPM Parameters of Laterite Mixed with Cement

Under the confining pressure of 100 kPa, the obtained model curve is shown in Figure5, and the fitted parameters are shown in Table 3. By analyzing the model parameters in Table 3, it can be found that when stable creep behavior occurs, the model parameters λ, K0, and G0 decayed with the increase of σ1eq−σ3 (Figures 6–8) and when destructive creep behavior occurs, λ, η0 and G0 decreased with the increase of σ1eq−σ3 (Figures 9–11). When the upper limit of the cyclic loading is near the critical dynamic stress, n was close to 0, and when the cyclic load increased, n was close to a fixed value. R in Table 3 is correlation coefficient.Figure 6
Relationship betweenG0 and σ1eq−σ3 for stable creep.Figure 7
Relationship betweenK0 and σ1eq−σ3 for stable creep.Figure 8
Relationship betweenλ and σ1eq−σ3 for stable creep.Figure 9
Relationship betweenλ and σ1eq−σ3 for destructive creep.Figure 10
Relationship betweenη0 and σ1eq−σ3 for destructive creep.Figure 11
Relationship betweenG0 and σ1eq−σ3 for destructive creep.


## 5.1.1. Determination of Parameters

1stop is leading software for nonlinear curve fitting and comprehensive optimization analysis in the world. It has been widely recognized in the field of nonlinear regression, fitting nonlinear curve, and parameter estimation of complex engineering models. In this paper, the software was used to solve the parameters of the proposed model.In equation (27), there are 3 parameters: G0, K0, and λ. In equation (28), there are 4 parameters: G0, λ, η0, and n. In order to verify the adaptability of the model, the PO (Powell Optimization) method and SA (Simulated Annealing) method in 1stop software are used to calculate the parameters of the constitutive equation for destructive and stable creep of soil under cyclic loading, respectively.The process of solving model parameters and fitting test data is as follows. According to equation (4), the equivalent stress under different stress amplitudes was obtained. When the stable creep behavior occurs, equation (27) and the corresponding experimental data were input into the software. When the destruction creep behavior occurs, the constitutive equation (28) and the corresponding test data were input, and then, the corresponding optimization algorithm for parameter calculation was selected. Finally, the curves and parameters of the proposed model could be obtained by running the program.


## 5.1.2. Calculation of EVPM Parameters for Laterite Mixed with Cement

Kong et al. [30] used the British GDS dynamic triaxial test system to carry out cyclic test on laterite mixed with cement under sine-wave periodic load with frequency f=5Hz; the stress level is D=σd/2τf (τf is the static shear strength) when testing. According to the test, the critical dynamic stress amplitudes σdc of laterite mixed with cement under confining pressure pc of 25, 50, and 100 kPa are 216, 274, and 339 kPa respectively; the static shear strength τf is 720 kPa, 917 kPa, and 1130 kPa, respectively. The test conditions are shown in Table 1. The vibration center of the axial dynamic stress is σ=pc+σd, and the axial dynamic stress is(29)σ1t=σ3+σd+σd2sinϖt.Table 1
Test conditions [30].
Confining pressure (kPa)Dynamic stress amplitudeσd (kPa)25254565951502252452605050100150175200285310—100100200300330360370380390According to equation (29), it is well known that σ1max=σ3+σd+σd/2, and the average value of the deviatoric stress under cyclic loading is σav=σd.Under 25, 50, and 100 kPa confining pressure, when the critical state occurs, the average value of cyclic loading is 216, 274, and 339 kPa; then, the critical strengthσcr=σdc+σdc/2 are 314, 411, and 508.5 kPa, respectively. The corresponding σ1max at critical state are 457, 598, and 778 kPa, respectively.n in equation (4) takes 1 here. The critical equivalent stress σ1eq of soil under 25, 50, and 100 kPa confining pressure is 243.54, 324.78, and 401.82 kPa, respectively.According to equation (10), the equivalent internal friction angle ϕeq=21.10∘ and the equivalent cohesion ceq=78.42kPa of laterite mixed with cement can be obtained by selecting the equivalent stress state under 25 and 100 kPa confining pressure. αeq=0.12 and Keq=76.73kPa can be obtained from equation (9).According to the obtained parameters, the constitutive equation (28) of laterite mixed with cement can be further simplified to(30)Δε11pN=σ1eq−σ33G0exp−λN/fλ+σ1eq+2σ39K0exp−λN/fλ+1.62σ1eq−1.17σ3−1.2776.73+13σav1η0nNfn.According to the data from Table1 and equation (4), the equivalent deviator stress σ1eq−σ3 under 25 and 100 kPa confining pressure can be calculated; the results of equivalent deviator stresses under different dynamic stress amplitudes σd are shown in Table 2.Table 2
The equivalent deviator stress of cyclic loading.
Confining pressure (kPa)Equivalentσ1eq−σ3 (kPa)2516.1030.4745.8172.52130.14234.18266.74293.1510070.46162.12282.53326.72371.33387.42403.50422.48It can be seen from Table2 that if σ1max−σ3>σcr of cyclic loading, the equivalent deviatoric stress σ1eq−σ3>σav obtained from equation (2). On the contrary, if σ1max−σ3<σcr, σ1eq−σ3<σav. In addition, the equivalent stress σ1eq−σ3 obtained was near σav of cyclic loading, indicating that the equivalent of equation (4) is reasonable and the effect of cyclic loading on the soil can be reflected by it.In this research, the obtained parameters of the model of soil under the confining pressure of 100 kPa are shown in Table3, and the comparison between model curves and test curves is shown in Figure 5.Table 3
EVPM parameters of laterite mixed with cement under 100 kPa confining pressure [30].
σd (kPa)G0 (kPa)K0 (kPa)λη0 (kPa·sn)nR1007.113 × 10153.969 × 1040.1840.8702005.735 × 10152.918 × 1040.1760.8923005.472 × 10153.000 × 1040.1790.8933304.960 × 10152.510 × 1040.1720.8663603.348 × 1030.1972.210 × 107−0.0030.9993704.158 × 1030.1856.793 × 1061.3580.9983803.702 × 1030.1571.917 × 1061.3400.9963902.679 × 1030.1166.035 × 1051.3570.999Figure 5
Fitting the curves of laterite mixed with cement [30] under the confining pressure of 100 kPa by EVPFM.From Figure5 and Table 3, it can be seen that the elastic-viscoplastic constitutive model based on the Bingham model can describe the stable and destructive creep behavior of soil under cyclic loading. The fitting precision of various deformation curves is good. The fitting correlation coefficient is above 0.870 and up to 0.999 under high and low dynamic stress. Therefore, the constitutive equation established in this paper can be used to reflect the various creep behavior of soil under cyclic loading.


## 5.1.3. Analysis of EVPM Parameters of Laterite Mixed with Cement

Under the confining pressure of 100 kPa, the obtained model curve is shown in Figure5, and the fitted parameters are shown in Table 3. By analyzing the model parameters in Table 3, it can be found that when stable creep behavior occurs, the model parameters λ, K0, and G0 decayed with the increase of σ1eq−σ3 (Figures 6–8) and when destructive creep behavior occurs, λ, η0 and G0 decreased with the increase of σ1eq−σ3 (Figures 9–11). When the upper limit of the cyclic loading is near the critical dynamic stress, n was close to 0, and when the cyclic load increased, n was close to a fixed value. R in Table 3 is correlation coefficient.Figure 6
Relationship betweenG0 and σ1eq−σ3 for stable creep.Figure 7
Relationship betweenK0 and σ1eq−σ3 for stable creep.Figure 8
Relationship betweenλ and σ1eq−σ3 for stable creep.Figure 9
Relationship betweenλ and σ1eq−σ3 for destructive creep.Figure 10
Relationship betweenη0 and σ1eq−σ3 for destructive creep.Figure 11
Relationship betweenG0 and σ1eq−σ3 for destructive creep.


## 5.2. Application Analysis

From the variation of model parameters under different dynamic stress amplitudes [31–33], when the confining pressure is 100 kPa, it is known that the value of each parameter varies obviously with σ1eq−σ3. According to the above analysis, it can be assumed that the parameter variation laws under any dynamic pressure and different dynamic stress are consistent with 100 kPa, but the numerical values are inconsistent. Under different confining pressures, the values of each parameter can be determined by a set of tests under specific confining pressures. Consequently, it can be assumed that the model parameters under a certain confining pressure with different dynamic stress amplitudes are determined by equations (31)–(33) when a stable creep behavior occurs:(31)K0=1.17×105σ1eq−σ3−0.259×β1sσ1eq,σ2,σ3,(32)G0=1.770×1016×σ1eq−σ3−0.216×β2sσ1eq,σ2,σ3,(33)λ=0.210×σ1eq−σ3−0.033×β3sσ1eq,σ2,σ3.In equations (31)–(33), β1s,β2s,andβ3s values can be obtained by plugging parameters K0, G0, and λ determined by a group test data from stable creep into equations (31)–(33). Of course, in order to get more accurate values, it can be determined by combining more test data.When destructive creep occurs, the parametersG0, λ, and η0 are determined by the following equations (34)–(36):(34)G0=1.425×1016×σ1eq−σ3−4.840×β1fσ1eq,σ2,σ3,(35)λ=1.159×σ1eq−σ3−3.795×β2fσ1eq,σ2,σ3,(36)η=8.285×1076σ1eq−σ3−28.241×β3fσ1eq,σ2,σ3In equations (34)–(36), β1f,β2f,andβ3f values can be obtained by plugging parameters G0, λ, and η0 determined by a group test data from destructive creep into equations (31)–(33). Of course, in order to get more accurate values, it can be determined by combining a set of test data.When the upper limit of stress is close to the critical dynamic stress, the parametern value tends to 0. Under the other different conditions of σd, it is temporarily determined to be 1.350 in combination with Table 3, but the specific change law of n should be further studied in combination with a large number of experiments.In this paper, when stable creep occurs,β1s=1/3.44, β2s=1/1.92, and β3s=1/0.946 are obtained by plugging the parameter K0, G0, and λ value under the dynamic stress amplitude of 25 kPa and 25 kPa confining pressure into equations (31)–(33). When calculating, because the part result of calculation is slightly smaller, s is taken as 1/4. Finally, the parameters K0, G0, and λ under dynamic stress amplitude of 45, 65, 95, 150, and 225 kPa are calculated in Table 4. The deformation under dynamic stress amplitudes 45, 65, 95, 150, and 225 kPa is predicted by using the calculated parameters. The predicted results are compared with the experimental results, as shown in Figure 12, when stable creep behavior occurs.Table 4
EVPM parameters of laterite mixed with cement under 25 kPa confining pressure [30].
σd (kPa)G0 (kPa)K0 (kPa)λη0 (kPa·sn)nRemark255.045 × 10151.665 × 1040.204Test data454.406 × 10151.413 × 1040.200Calculation654.036 × 10151.270 × 1040.197953.656 × 10151127 × 1040.1931503.219 × 10159.709 × 1030.1902252.836 × 10158313 × 1030.18624532540.2169.160 × 106-0.009Test data26020570.1526.305 × 1051.350CalculationFigure 12
Comparison between model predictive value and test data of laterite mixed with cement [30].When the destructive creep behavior occurs, in this paper, the parameter value ofβ1f,β2f,andβ3f can be obtained by combining the experimental parameter value under confining pressure 25 kPa and dynamic stress amplitude 245 kPa with equations (32)–(34). β1f=0.126, β2f=0.302, and β3f=6.476×10−5 are obtained by plugging the parameter K0, η0, and λ value under the dynamic stress amplitude of 25 kPa and 25 kPa confining pressure into equations (32)–(34). According to these obtained parameters, the creep of soil under cyclic loading with a stress amplitude of 260 kPa is calculated, and the predicted results are compared with the experimental results, as shown in Figure 13.Figure 13
Comparison between model predictive value and test data of laterite mixed with cement [30] when destructive deformation occurs.The calculation results of soil deformation under confining pressure 25 kPa show that the proposed EVPM can predict the deformation of soil under cyclic loading. There is a certain error between the calculated value of the model and the experimental value. The main reason for the error is that the experimental data are few, which leads to the inaccuracy of the changing rules of the parameters and the calculation error. Further research is needed in combination with experiments. However, the calculated results of the model are roughly consistent with the experimental data, which further illustrates the feasibility of the model.


## 6. Conclusions

In this paper, the elastic-viscoplastic model of soil under cyclic loading is established based on the Bingham model and the yield criterion of soil under equal amplitude cyclic loading is proposed. The major findings are summarized as follows:(1)
Cyclic loading is simplified equivalently; a dynamic yield criterion of soil under cyclic loading is proposed in view of the simplified load. Based on the improved Bingham model, a EVPFM of soil is established by using the flow criterion related to the dynamic yield criterion. Meanwhile, the parameters of the model are determined and the adaptability of the model is verified.(2)
The stable and destructive creep behavior of soil under cyclic loading can be reflected by the established EVPM, and the relevant parameters in the model can also be easily obtained by numerical methods.(3)
When the peak deviatoric stressσ1max−σ3 of periodic load is smaller than the critical dynamic stress σcr value of soil, the model is a nonlinear model describing the stable creep behavior of the soil. Instead, the model is an EVPM describing the soil destructive creep behavior.(4)
When stable creep behavior occurs, the model has three parameters, namely, the initial value of the elastic shear fatigue coefficientG0, the initial value of the volume fatigue coefficient K0, and the elastic attenuation coefficient λ. All three parameters decay with increasing equivalent deviatoric stress σ1eq−σ3.(5)
When destructive creep behavior occurs, the model has four parameters,G0, λ, η0, and n, respectively. All three parameters decay with increasing equivalent deviatoric stress σ1eq−σ3. The law of change of n needs to be further studied.(6)
Comparing the calculated results of strain under different dynamic stresses with the experimental results, the established model can predict the deformation of soil under cyclic loading.



---

*Source: 1602431-2021-09-01.xml*