# Construction of 3D Creep Model of Landslide Slip-Surface Soil and Secondary Development Based on FLAC3D

Authors: Li Wang; Shimei Wang; Gao Li; Lin Wang
Journal: Advances in Civil Engineering
 (2020)

---

## Abstract

Landslide deformations have predictable rheological properties, but, when using finite element analysis, only the elastoplastic constitutive equation of soil is considered. Soil creep properties are not considered, often rendering inaccurate results. In this paper, the three-dimensional state equation of the extended Burgers model is derived based on the creep test results of slip-surface soil and the Mohr–Coulomb elastoplastic model. This paper uses FLAC3D, a finite-element software platform, for secondary development and to verify the accuracy of the model and program. A strength reduction method that considers rheological properties is proposed. A numerical simulation of rheological properties is conducted on a landslide and compared with conventional viscoelastic constitutive results. It can be found that the landslide displacement calculation is significantly smaller when the rheological properties are not considered. The stability coefficient of the landslide calculated by the strength reduction method, considering the rheological properties, is smaller than the coefficient calculated without considering rheological properties.


---

## Body

## 1. Introduction

The understanding of the influence of water-level changes on the deformation of reservoir landslides has developed and changed over time [1, 2]. For example, the deformation of the Bazimen and Shuping landslides in the China Three Gorges reservoir area showed obvious rheological properties [3, 4], such as the stepwise Global Positioning System displacement of the surface with a decrease of the reservoir water level. Even when the reservoir and groundwater levels were stable, its deformation continued to increase. When using a general finite element analysis, it is impractical to consider only the elastoplastic constitutive model to simulate the landslide problem. Instead, the rheological model should be applied to landslide deformation and long-term stability prediction.Wang et al. [5] paid much attention to the time effect of landslide deformation, believing that the instability of landslides was closely related to the viscoelastic properties of rock masses. Several other authors [6–11] noted significant creep effects in landslide deformation and proposed indicators of reactive creep properties using predictive models. Rutter and Green [12] conducted creep monitoring on the Mam Tor landslide and established a nonlinear empirical model of pore-water pressure and creep. Other long-term deformation predictions for landslide bodies have been studied [13–17], in which predictions of landslides were based on mathematical models. However, the mechanical mechanism of creep was not considered. Furthermore, the meaning of model parameters used was not explained. Chang et al. [18] used the Mohr–Coulomb and anisotropic models to simulate and explain the causes of slope change and creep behaviours. Wang et al. [19] provided model parameters for numerical analyses of landslide deformation based on the results of slip-zone soil rheological tests. To consider the creep properties of rock and soil using finite-element simulation, current geotechnical finite-element software contains partial creep calculation modules. For example, ABAQUS software provides related creep models, including the time-hardening creep model, the strain-hardening creep variable model, and the Singh–Mitchell creep model. FLAC software also provides related creep calculation modules. However, these models have limitations, and the respective creep models are relatively small, presenting gaps between the actual states of the rock and soil. Most numerical analysis models apply an empirical creep model, causing the meaning of parameters to remain unclear. Furthermore, the creep models in commercial software have limited adaptability. With these tools, some scholars had achieved decent results in rheological model parameter fitting. However, the commonly fitted parameters are those of a 1-dimensional (1D) rheological model. These models do not reflect 3-dimensional (3D) soil characteristics, making it difficult to directly apply to engineering problems.Based on the triaxial creep test data of a landslide slip-surface soil, this paper expands the Burgers creep model and derives a 3D creep model that is more suitable for rheological landslide properties. This paper also solves the parameters of the 3D model. The extended Burgers model is thus deduced into a 3D finite difference form combined with the Mohr–Coulomb elastoplastic model. Thus, the secondary development functions of FLAC3D software are used to write a user-defined constitutive model based on C++, and the correctness of the custom model is tested. Then, the custom model is applied to the rheological analysis of the landslide. Furthermore, this paper proposes a strength reduction method that considers the rheological properties of soils, which is then applied to the calculation and compared with the conventional finite element strength reduction method.


## 2. Triaxial Creep Test of Slip-Surface Soil

In this study, a creep test was performed on slip-surface soil using a triaxial creep apparatus, using the detailed test steps provided in a past study [20].The confining pressure system controlled by air pressure in the triaxial creep apparatus kept the pressure constant and exerted axial stress on the soil. This method, combined with the conventional gravity-loading method and the selected varied weight loading method, can keep the axial stress at a constant value for a long time. If changes to the cross-sectional area of sample soil are ignored with time and loading, constant axial stress can be ensured. The triaxial rheological instrument is shown in Figure 1.Figure 1
Triaxial rheological apparatus.A test soil specimen was taken from the deep slip surface of the Dayanyu landslide, located downstream of the Shuibuya Dam in the Qingjiang River, Hubei Province, China. The basic physical and mechanical properties of soil samples are tested in accordance with China geotechnical engineering test code (SL237-1999) [21]. These basic properties are shown in Table 1.Table 1
Basic physical parameters of slide soil.
Specific gravityGs
Moisture contentw (%)
Densityρ (g/cm3)
Dry densityρd (g/cm3)
Void ratioe
Liquid limitwL (%)
Plastic limitwP (%)
Cohesionc (kPa)
Internal friction angleΦ (°)
1
23.5
2.00
1.70
0.30
41.5
17
14.0
22.0In the triaxial creep test, a constant load (i.e., stress) was applied to the tested soil sample, and the strain values of the soil at different times under the same stress were recorded. The experiment was divided into three groups of confining pressure:σ3 = 100 kPa; σ3 = 200 kPa; and σ3 = 300 kPa. The main steps are described in the following.For each consolidation pressure level, five different deviating stresses were implemented. If the conventional consolidation undrained shear strength wasτcu, the following five kinds of deviator stresses were used: 0.3τcu, 0.5τcu, 0.6τcu, 0.7τcu, and 0.8τcu.Creep test data ofσ3 at 100 kPa, 200 kPa, and 300 kPa were selected. The creep test curves and the stress—strain—time relationship diagrams are drawn as shown in Figure 2. The creep process generally includes an instantaneous deformation stage, an attenuation deformation stage, and an isovelocity deformation stage. Because the soil sample in our case was a highly plastic slip-surface soil, the plastic deformation was large, and there was no accelerated creep stage in the test range. The creep of each stage increased considerably with the increase of axial stress level.Figure 2
Creep curve under different confining pressures: (a) 100 kPa; (b) 100 kPa (first 40 min); (c) 200 kPa; (d) 200 kPa (first 40 min); (e) 300 kPa; (f) 300 kPa (first 40 min).
(a)
(b)
(c)
(d)
(e)
(f)


## 3. 3D Creep Model Construction of Slip-Surface Soil

### 3.1. Establishment of the 3D Burgers Creep Equation

The Burgers model is a rheological model that is consistent with the rheological properties of the slip surface and can be fitted with parameters using experimental data. It can reflect the transient and attenuation stages of creep deformation. The prolonged burner creep model is actually a Maxwell body in a series havingn Kelvin bodies, as shown in Figure 3, whose creep constitutive relationship is as follows [22]:(1)εt=σ^ctηm+σ^cEm+σ^c∑i=1n1−exp−Eiηit.Figure 3
Extended Burgers creep model (E0, E1, and E2 are the elastic moduli, and η0, η1, and η2 are the viscosity coefficients).In practice, the soil is placed in a 3D stress state. If the 1D creep model is not applicable to determine the stress and strain of the soil, it becomes necessary to derive a universally applicable 3D stress state creep equation using a few basic assumptions [21]. The soil is set as an ideal viscoelastic body. The shape change caused by soil rheology is mainly caused by deviating stress, and ’Poisson’s ratio of soil remains unchanged. The 3D form of the extended burgers model is thus obtained as follows [21]:(2)εt=σm3K+Sij21G0+1η0t+∑i=121Gi1−exp−Giηit,where σm is the spherical stress tensor, sij is the deviator stress tensor, and Gis the shear modulus. When time t = 0, triaxial rheological experiments are carried out on the specimens. Furthermore, the applied axial stress, σ1, the confining pressure, σ3, and the specimens are assumed to be continuously isotropic. Thus, the creep stress tensor state can be expressed as follows:(3)σ11t=σ1,σ22=σ33=σ3,σ12=σ13=σ23=0,where subscripts 1, 2, and 3 represent three orthogonal coordinate axes of the Cartesian coordinate system and coordinate axis 1 represents the axial direction of the test piece. Therefore,(4)σkk=σ1+2σ3,S11=σ11−13σkk=23σ1−σ3,because σm=σkk/3=1/3σ11+σ22+σ33=1/3σ1+σ2+σ3. Using the above formula, the extended Burgers (M-2K) model axial creep equation in the unsaturated triaxial rheological test can be obtained:(5)ε1t=ε11=e11+σkk9K=σΙΙ9K+σΙ31G1+1η1t+∑i=1n1Gi1−exp−Giηit,where σΙ=σ1−σ3,σΙΙ=σ1+2σ3, and e11 is the corresponding partial strain of S11. The other parameters have the same meaning as before.


### 3.2. Determining the 3D Creep Parameters of the Extended Burgers Model

A set of known experimental data were fitted using least squares. The function is as follows:(6)y=fX,b1,b2,…,bm.Because the Burgers creep model has negative exponential terms, a set of initial values can be initially obtained. The Matlablsqcurvefit nonlinear regression tool was used to obtain the model parameters. This brings us to our first hypothesis.HypothesisA=σΙΙ/9K+σΙ/3G0, B=σΙ/3η0, C=σΙ/3G1, D=G1/η1, E=σΙ/3G2, and F=G2/η2. Therefore, the creep equation (5) of the Burgers model can be rewritten as(7)ε1t=A+Bt+C1−e−Dt+E1−e−Ft.Lett = 0. The right-hand end simplification of equation (7) is equal to A, and the parameter A can be determined using creep data of the instantaneous deformation stage. If the sample creeps long enough, the last two parameters, C1−e−Dt+E1−e−Ft, tend to be constant: C+E. When the initial value is selected, C=EandD=Fcan be approximated. From the decay-creep phase, equation (7) can be simplified as to(8)ε1t=A+2C1−e−Dt.The constant term in the above formula is moved to the left side, and logarithm is taken to obtain a logarithmic coordinate system. The linear equation then becomes(9)lnε1t−A−2C=ln−2C−Dt,where A is determined and D is the slope of the line. Thus, the parameter C can be determined. From this, a set of initial values can also be determined. Then lsqcurvefit nonlinear regression tool in MATLAB software is used to solve the stable values of A, B, C, and D. The bulk modulus K=E/31−2μ, the shear modulus G=E/21+μ, and the μvalue is 0.4, which allows us to get A=σΙΙ/9K+σΙ/3G0=15σ1−12σ3/14G0. According to the above method, the 3D creep parameters of the Burgers model, having a confining pressure of 100 kPa, 200 kPa, and 300 kPa, can be obtained. Limited by article length, only the parameters associated with the confining pressure 300 kPa are shown in Table 2.Table 2
Three-axis rheological experiment fitted according to the improved Burgers model (confining pressure: 300 kPa).
Deviatoric stressq
Shear modulusG0
Shear modulusG1
Shear modulusG2
Viscosity coefficientη0
Viscosity coefficientη1
Viscosity coefficientη2
Bulk modulusK
Correlation coefficientR2
86.3
4.081E + 04
2.520E + 04
2.520E + 04
5.364e + 08
4.352E + 06
4.352E + 06
1.904E + 05
0.9924
141.6
7.301E + 03
1.403E + 04
1.403E + 04
5.643e + 08
6.031E + 05
6.0318E + 05
3.407E + 04
0.9576
168.6
5.507E + 03
7.225E + 03
6.783E + 03
5.243e + 08
4.042E + 05
2.420E + 04
2.569E + 04
0.9765
196.1
4.113E + 03
2.097E + 03
1.949E + 03
4.211e + 08
1.753E + 04
1.932E + 03
1.919E + 04
0.9462
211.8
2.466E + 03
1.582E + 03
1.521E + 03
2.794e + 08
9.722E + 03
1.934E + 03
1.151E + 04
0.9534


## 3.1. Establishment of the 3D Burgers Creep Equation

The Burgers model is a rheological model that is consistent with the rheological properties of the slip surface and can be fitted with parameters using experimental data. It can reflect the transient and attenuation stages of creep deformation. The prolonged burner creep model is actually a Maxwell body in a series havingn Kelvin bodies, as shown in Figure 3, whose creep constitutive relationship is as follows [22]:(1)εt=σ^ctηm+σ^cEm+σ^c∑i=1n1−exp−Eiηit.Figure 3
Extended Burgers creep model (E0, E1, and E2 are the elastic moduli, and η0, η1, and η2 are the viscosity coefficients).In practice, the soil is placed in a 3D stress state. If the 1D creep model is not applicable to determine the stress and strain of the soil, it becomes necessary to derive a universally applicable 3D stress state creep equation using a few basic assumptions [21]. The soil is set as an ideal viscoelastic body. The shape change caused by soil rheology is mainly caused by deviating stress, and ’Poisson’s ratio of soil remains unchanged. The 3D form of the extended burgers model is thus obtained as follows [21]:(2)εt=σm3K+Sij21G0+1η0t+∑i=121Gi1−exp−Giηit,where σm is the spherical stress tensor, sij is the deviator stress tensor, and Gis the shear modulus. When time t = 0, triaxial rheological experiments are carried out on the specimens. Furthermore, the applied axial stress, σ1, the confining pressure, σ3, and the specimens are assumed to be continuously isotropic. Thus, the creep stress tensor state can be expressed as follows:(3)σ11t=σ1,σ22=σ33=σ3,σ12=σ13=σ23=0,where subscripts 1, 2, and 3 represent three orthogonal coordinate axes of the Cartesian coordinate system and coordinate axis 1 represents the axial direction of the test piece. Therefore,(4)σkk=σ1+2σ3,S11=σ11−13σkk=23σ1−σ3,because σm=σkk/3=1/3σ11+σ22+σ33=1/3σ1+σ2+σ3. Using the above formula, the extended Burgers (M-2K) model axial creep equation in the unsaturated triaxial rheological test can be obtained:(5)ε1t=ε11=e11+σkk9K=σΙΙ9K+σΙ31G1+1η1t+∑i=1n1Gi1−exp−Giηit,where σΙ=σ1−σ3,σΙΙ=σ1+2σ3, and e11 is the corresponding partial strain of S11. The other parameters have the same meaning as before.


## 3.2. Determining the 3D Creep Parameters of the Extended Burgers Model

A set of known experimental data were fitted using least squares. The function is as follows:(6)y=fX,b1,b2,…,bm.Because the Burgers creep model has negative exponential terms, a set of initial values can be initially obtained. The Matlablsqcurvefit nonlinear regression tool was used to obtain the model parameters. This brings us to our first hypothesis.HypothesisA=σΙΙ/9K+σΙ/3G0, B=σΙ/3η0, C=σΙ/3G1, D=G1/η1, E=σΙ/3G2, and F=G2/η2. Therefore, the creep equation (5) of the Burgers model can be rewritten as(7)ε1t=A+Bt+C1−e−Dt+E1−e−Ft.Lett = 0. The right-hand end simplification of equation (7) is equal to A, and the parameter A can be determined using creep data of the instantaneous deformation stage. If the sample creeps long enough, the last two parameters, C1−e−Dt+E1−e−Ft, tend to be constant: C+E. When the initial value is selected, C=EandD=Fcan be approximated. From the decay-creep phase, equation (7) can be simplified as to(8)ε1t=A+2C1−e−Dt.The constant term in the above formula is moved to the left side, and logarithm is taken to obtain a logarithmic coordinate system. The linear equation then becomes(9)lnε1t−A−2C=ln−2C−Dt,where A is determined and D is the slope of the line. Thus, the parameter C can be determined. From this, a set of initial values can also be determined. Then lsqcurvefit nonlinear regression tool in MATLAB software is used to solve the stable values of A, B, C, and D. The bulk modulus K=E/31−2μ, the shear modulus G=E/21+μ, and the μvalue is 0.4, which allows us to get A=σΙΙ/9K+σΙ/3G0=15σ1−12σ3/14G0. According to the above method, the 3D creep parameters of the Burgers model, having a confining pressure of 100 kPa, 200 kPa, and 300 kPa, can be obtained. Limited by article length, only the parameters associated with the confining pressure 300 kPa are shown in Table 2.Table 2
Three-axis rheological experiment fitted according to the improved Burgers model (confining pressure: 300 kPa).
Deviatoric stressq
Shear modulusG0
Shear modulusG1
Shear modulusG2
Viscosity coefficientη0
Viscosity coefficientη1
Viscosity coefficientη2
Bulk modulusK
Correlation coefficientR2
86.3
4.081E + 04
2.520E + 04
2.520E + 04
5.364e + 08
4.352E + 06
4.352E + 06
1.904E + 05
0.9924
141.6
7.301E + 03
1.403E + 04
1.403E + 04
5.643e + 08
6.031E + 05
6.0318E + 05
3.407E + 04
0.9576
168.6
5.507E + 03
7.225E + 03
6.783E + 03
5.243e + 08
4.042E + 05
2.420E + 04
2.569E + 04
0.9765
196.1
4.113E + 03
2.097E + 03
1.949E + 03
4.211e + 08
1.753E + 04
1.932E + 03
1.919E + 04
0.9462
211.8
2.466E + 03
1.582E + 03
1.521E + 03
2.794e + 08
9.722E + 03
1.934E + 03
1.151E + 04
0.9534


## 4. Secondary Development of FLAC3D Based on the Slip-Surface Soil Rheological Model

The 3D Fast Lagrangian Analysis of Continua (FLAC3D) tool is a 3D explicit finite difference program developed by ITASCA, USA, which is applicable to the mechanical analysis of civil engineering problems. It can better simulate the mechanical behaviour of material failures and plastic flows. This paper uses the platform developed by FLAC3D to realise the development of the extended Burgers viscoelasticity model, verifying the model’s correctness.


### 4.1. Finite Difference Form of the Extended Burgers Viscoelastic Model

The extended Burgers viscoelastic model comprises two Kelvin bodies and one Maxwell body in a series. In the following formula, the superscriptN represents the new value, and the superscript O represents the old value. The Kelvin body strain rate is(10)ε˙1=σ1β1,and the viscous force is(11)σ1=σ−E1ε1.In the kelvin volume,ε˙1=σ−E1ε1/β1, from equation (20), can be obtained from equations (10) and (11). They are sorted as follows:(12)ε1N−ε1oΔt=σ¯−E1ε¯1β1=σN+σo−E1ε1N+ε1o2β1.Organize the equation as(13)ε1N=ε1o+σN+σo−E1ε1N+ε1oΔt2β1.Obtain the following formula in the same way:(14)ε2N=ε2o+σN+σo−E2ε2N+ε2oΔt2β2.For the Maxwell body,(15)ε˙3=ε˙m=σ˙E3+σ¯β3.The transformation of the Maxwell body gives us(16)εmN=εmo+σN−σoE3+σN+σo2β3⋅Δt.According to the above formula, the total strain increment is(17)Δε=εN−εo=ε1N−ε1o+ε2N−ε2o+ε3N−ε3o.According to equation (13), the following can be obtained:(18)1+E1Δt2β1ε1N=1−E1Δt2β1ε1o+σN+σoΔt2β1,which can be simplified as follows:(19)ε1N=1ABε1o+σN+σo⋅Δt2β1.Then, obtain the following formula in the same way:(20)1+E2Δt2β2ε2N=1−E2Δt2β2ε2o+σN+σoΔt2β2,which is simplified as follows:(21)ε2N=1CDε2o+σN+σo⋅Δt2β2.According to (13), (14), (16), and (17), the following can be obtained:(22)σN=1XεN−εo+Yσo−BA−1ε1o−DC−1ε2o.Therefore, the 1D M-2k model (i.e., the finite difference form of the improved Burgers viscoelastic model) is obtained:(23)ε1N=1ABε2o+σN+σo⋅Δt2β2,ε2N=1CDε2o+σN+σo⋅Δt2β2,σN=1XεN−εo+Yσo−BA−1ε1o−DC−1ε2o,where A=1+E1Δt/2β1, B=1−E1Δt/2β1, C=1+E2Δt/2β2, and D=1−E2Δt/2β2. X=1/A⋅Δt/2β1+1/C⋅Δt/2β2+1/E3+Δt/2β3. Y=1/E3−1/A⋅Δt/2β1−1/C⋅Δt/2β2−Δt/2β3.According to the method derived from the 1D to the 3D rheological model, the finite difference form of 3D model can be obtained as follows:(24)eij1N=1abeij1o+SijN+Sijo⋅Δt4η1,eij2N=1cdeij2o+SijN+Sijo⋅Δt4η2,SijN=1xΔeij+ySijo−ba−1eij1o−dc−1eij2o,Δσm=KΔεkk,σijN=SijN+13σkkN+Δσmδij,where Δeij=Δεij−1/3Δεijδij, a=1+G1Δt/2η1, b=1−G1Δt/2η1, c=1+G2Δt/2η2, and d=1−G2Δt/2η2. x=1/A⋅Δt/4η1+1/C⋅Δt/4η2+1/2G3+Δt/4η3. y=1/2G3−1/A⋅Δt/4η1−1/C⋅Δt/4η2−Δt/4η3.The Mohr–Coulomb strain rate is as follows:(25)e˙ijp=λ∗∂g∂σij−13e˙volpδij,e˙volp=λ∗∂g∂σ11+∂g∂σ22+∂g∂σ33.The Mohr–Coulomb yield criterion isf = 0, and the formula on the principal axis is the shear yield:(26)f=σ1−σ3Nϕ+2CNϕ.The tension yield is as follows:(27)f=σt−σ3,where C and ϕ are the cohesion and friction angle of the material:(28)Nϕ=1−sinϕ1+sinϕ,where σt is the tension strength. The potential function, g, has the following form, where the shear damage is(29)g=σ1−σ3Nψ.The tensile damage is as follows:(30)g=−σ3,where Ψ is the material dilatancy angle and Nψ=1+sinψ1−sinψ. λ is determined by the yield condition, f = 0, and is not equal to zero in the plastic flow phase.The stick-elastic increment form is adopted in FLAC3D, and the new trial stress component is determined by equation (24). The principal stress components are calculated and classified according to the trial stress, and the yield function is obtained. If f  ≥  0, it is considered that the trial stress is a new stress. If f < 0, it is considered that plastic flow occurs. The plastic strain increment is calculated, and the value of the trial stress must be corrected by the plastic strain increment before the application of new stress. The method of correcting the trial stress is described in the following.The stress is defined according to the following partial component:(31)σ1N=σ^1N−α1Δε1p+α2Δε2p+Δε3p,σ2N=σ^2N−α1Δε2p+α2Δε1p+Δε3p,σ3N=σ^3N−α1Δε3p+α2Δε1p+Δε2p,where α1=K+2/3x and α2=K+2/3x. For shear yielding, the stress can be obtained as(32)σ1N=σ^1N−λα1−α2Nϕ,σ2N=σ^2N−λα21−Nϕ,σ3N=σ^3N−λα2−α1Nϕ,where λ=σ^1N−σ^3NNϕ+2CNϕ/α1−α2Nϕ−α2−α1NϕNϕ. For tensile yielding, the stress can be obtained as follows:(33)σ1N=σ^1N+λα2,σ2N=σ^2N+λα2,σ3N=σ^3N+λα1,where λ=σt−σ^3N/α1.


### 4.2. Program Verification

To verify the correctness of the custom model writing, model parameters were selected as the result of the fitting the test data under the conditions of a confining pressure of 300 kPa and an axial deviatoric stress of 168.6 kPa. The cylindrical sample model is shown in Figure4. The sample diameter is 61.8 mm, the height is 120 mm, the time unit is minutes, and the density unit is g/mm3. The bottom surface is constrained, and the upper surface is subjected to a deviatoric stress of 39.2 kPa and a confining pressure of 100 kPa. It can be seen from Figure 5 that the extended Burgers viscoelastic-plastic model fits well with the experimental results, illustrating the rationality of the custom model.Figure 4
Calculation model of the specimen.Figure 5
Comparison of theoretical and experimental values of the extended Burgers visco-elastoplastic model.


## 4.1. Finite Difference Form of the Extended Burgers Viscoelastic Model

The extended Burgers viscoelastic model comprises two Kelvin bodies and one Maxwell body in a series. In the following formula, the superscriptN represents the new value, and the superscript O represents the old value. The Kelvin body strain rate is(10)ε˙1=σ1β1,and the viscous force is(11)σ1=σ−E1ε1.In the kelvin volume,ε˙1=σ−E1ε1/β1, from equation (20), can be obtained from equations (10) and (11). They are sorted as follows:(12)ε1N−ε1oΔt=σ¯−E1ε¯1β1=σN+σo−E1ε1N+ε1o2β1.Organize the equation as(13)ε1N=ε1o+σN+σo−E1ε1N+ε1oΔt2β1.Obtain the following formula in the same way:(14)ε2N=ε2o+σN+σo−E2ε2N+ε2oΔt2β2.For the Maxwell body,(15)ε˙3=ε˙m=σ˙E3+σ¯β3.The transformation of the Maxwell body gives us(16)εmN=εmo+σN−σoE3+σN+σo2β3⋅Δt.According to the above formula, the total strain increment is(17)Δε=εN−εo=ε1N−ε1o+ε2N−ε2o+ε3N−ε3o.According to equation (13), the following can be obtained:(18)1+E1Δt2β1ε1N=1−E1Δt2β1ε1o+σN+σoΔt2β1,which can be simplified as follows:(19)ε1N=1ABε1o+σN+σo⋅Δt2β1.Then, obtain the following formula in the same way:(20)1+E2Δt2β2ε2N=1−E2Δt2β2ε2o+σN+σoΔt2β2,which is simplified as follows:(21)ε2N=1CDε2o+σN+σo⋅Δt2β2.According to (13), (14), (16), and (17), the following can be obtained:(22)σN=1XεN−εo+Yσo−BA−1ε1o−DC−1ε2o.Therefore, the 1D M-2k model (i.e., the finite difference form of the improved Burgers viscoelastic model) is obtained:(23)ε1N=1ABε2o+σN+σo⋅Δt2β2,ε2N=1CDε2o+σN+σo⋅Δt2β2,σN=1XεN−εo+Yσo−BA−1ε1o−DC−1ε2o,where A=1+E1Δt/2β1, B=1−E1Δt/2β1, C=1+E2Δt/2β2, and D=1−E2Δt/2β2. X=1/A⋅Δt/2β1+1/C⋅Δt/2β2+1/E3+Δt/2β3. Y=1/E3−1/A⋅Δt/2β1−1/C⋅Δt/2β2−Δt/2β3.According to the method derived from the 1D to the 3D rheological model, the finite difference form of 3D model can be obtained as follows:(24)eij1N=1abeij1o+SijN+Sijo⋅Δt4η1,eij2N=1cdeij2o+SijN+Sijo⋅Δt4η2,SijN=1xΔeij+ySijo−ba−1eij1o−dc−1eij2o,Δσm=KΔεkk,σijN=SijN+13σkkN+Δσmδij,where Δeij=Δεij−1/3Δεijδij, a=1+G1Δt/2η1, b=1−G1Δt/2η1, c=1+G2Δt/2η2, and d=1−G2Δt/2η2. x=1/A⋅Δt/4η1+1/C⋅Δt/4η2+1/2G3+Δt/4η3. y=1/2G3−1/A⋅Δt/4η1−1/C⋅Δt/4η2−Δt/4η3.The Mohr–Coulomb strain rate is as follows:(25)e˙ijp=λ∗∂g∂σij−13e˙volpδij,e˙volp=λ∗∂g∂σ11+∂g∂σ22+∂g∂σ33.The Mohr–Coulomb yield criterion isf = 0, and the formula on the principal axis is the shear yield:(26)f=σ1−σ3Nϕ+2CNϕ.The tension yield is as follows:(27)f=σt−σ3,where C and ϕ are the cohesion and friction angle of the material:(28)Nϕ=1−sinϕ1+sinϕ,where σt is the tension strength. The potential function, g, has the following form, where the shear damage is(29)g=σ1−σ3Nψ.The tensile damage is as follows:(30)g=−σ3,where Ψ is the material dilatancy angle and Nψ=1+sinψ1−sinψ. λ is determined by the yield condition, f = 0, and is not equal to zero in the plastic flow phase.The stick-elastic increment form is adopted in FLAC3D, and the new trial stress component is determined by equation (24). The principal stress components are calculated and classified according to the trial stress, and the yield function is obtained. If f  ≥  0, it is considered that the trial stress is a new stress. If f < 0, it is considered that plastic flow occurs. The plastic strain increment is calculated, and the value of the trial stress must be corrected by the plastic strain increment before the application of new stress. The method of correcting the trial stress is described in the following.The stress is defined according to the following partial component:(31)σ1N=σ^1N−α1Δε1p+α2Δε2p+Δε3p,σ2N=σ^2N−α1Δε2p+α2Δε1p+Δε3p,σ3N=σ^3N−α1Δε3p+α2Δε1p+Δε2p,where α1=K+2/3x and α2=K+2/3x. For shear yielding, the stress can be obtained as(32)σ1N=σ^1N−λα1−α2Nϕ,σ2N=σ^2N−λα21−Nϕ,σ3N=σ^3N−λα2−α1Nϕ,where λ=σ^1N−σ^3NNϕ+2CNϕ/α1−α2Nϕ−α2−α1NϕNϕ. For tensile yielding, the stress can be obtained as follows:(33)σ1N=σ^1N+λα2,σ2N=σ^2N+λα2,σ3N=σ^3N+λα1,where λ=σt−σ^3N/α1.


## 4.2. Program Verification

To verify the correctness of the custom model writing, model parameters were selected as the result of the fitting the test data under the conditions of a confining pressure of 300 kPa and an axial deviatoric stress of 168.6 kPa. The cylindrical sample model is shown in Figure4. The sample diameter is 61.8 mm, the height is 120 mm, the time unit is minutes, and the density unit is g/mm3. The bottom surface is constrained, and the upper surface is subjected to a deviatoric stress of 39.2 kPa and a confining pressure of 100 kPa. It can be seen from Figure 5 that the extended Burgers viscoelastic-plastic model fits well with the experimental results, illustrating the rationality of the custom model.Figure 4
Calculation model of the specimen.Figure 5
Comparison of theoretical and experimental values of the extended Burgers visco-elastoplastic model.


## 5. Engineering Application

By using the modified improved Burgers viscoelastic-plastic model in FLAC3D, the rheological properties of a landslide can be analysed, and the deformation development trend and long-term stability can be predicted and compared with the calculation results of the elastoplastic model without considering the rheological properties.


### 5.1. Landslide Overview

The Qingjiang Shuibuya Water Control Project is located in Badong County, Hubei Province, and is an important water conservancy project in the Qingjiang cascade development. The Dayanyu landslide body is located on the left bank of the Qingjiang River at the lower reaches of the dam site. The landslide is 1.1 km from the dam axis. It was originally a microsand bedrock landslide, and it evolved into a main landslide body and with similar bodies on both sides via local disintegration. The complex main sliding body is elongated, and the upper sliding body is slow and gentle. The secondary sliding bodies on both sides have fan shapes. The total area of the landslide is 0.196 km2, the thickness of the sliding body is 25–40 m, the thickest part is 64.8 m, and the total volume is ∼5.88 × 106 m3.Because the landslide is in the area along the left bank’s antiscouring wall of the dam spillway scouring pit, its stability is directly affected by construction and spillway flooding and fogging, which is a great threat to the safe operation of the engineering building. Thus, it has always been the focus of attention.


### 5.2. Model Construction

The landslide consists of three parts: the slip mass, the slip surface, and the slip bed. The landslide material mesh division is shown in Figure6. In the figure, red indicates the slip mass, green indicates the slip surface, and blue indicates the slip bed. The calculation model divides 1,680 nodes and 811 units. The direction along the landslide is the x-axis, the height direction is the y-axis, and the unit thickness is along the z-axis. The landslide is 950 m along the x direction and 450 m along the y direction.Figure 6
Zoning map of landslide material.


### 5.3. Material Calculation Parameter

In the rheological calculation, the improved Burgers viscoelastic model is used for slip mass, the slip surface, and the slip bed. The material calculation parameters are obtained by analogy and trial calculation according to the fitting parameters of the triaxial creep test data of a landslide, as shown in Table3.Table 3
Improved parameter table of the Burgers visco-elastoplastic model.
Materials
Shear modulusG0 (kPa)
Shear modulusG1 (kPa)
Shear modulusG2 (kPa)
Viscosity coefficientη0 (kPa)
Viscosity coefficientη1 (kPa)
Viscosity coefficientη2 (kPa)
Cohesionc (kPa)
Internal friction angleφ (°)
Elastic modulusE (kPa)
Unit weightγ (kN/m3)
Tensile strengthσt (kPa)
Slip mass
1.659E + 05
8.721E + 04
7.666E + 04
6.063E + 09
3.059E + 06
3.268E + 06
18
25
7.745E + 05
28
0
Slip surface
9.221E + 03
4.845E + 03
4.259E + 03
3.368E + 08
1.699E + 05
1.815E + 05
14
22
4.303E + 04
20
0
Slip bed
1.659E + 07
8.721E + 06
7.666E + 06
6.063E + 11
3.059E + 08
3.268E + 08
32
580
7.745E + 07
40
3500


### 5.4. Boundary Conditions and Load

The model boundary constraint requires the front and rear sides of the landslide to be fixed in thex and z directions. The two sides perpendicular to the z-axis are fixed in the x and z directions. The nodes on the four faces can be deformed along the y direction, but there is no deformation in other directions. Fixed constraints exist in three directions. The duration of the landslide rheology is 300 days. Only the dead-weight load is considered in the calculation. When the rheological calculation is performed, the initial stress state of the landslide is calculated first. Then, the rheological calculations have been started.


### 5.5. Calculation Results

The results of the landslide rheological calculation show displacement and stress field distributions. Only the displacement field was used as the calculation result. When analysing the rheology, the displacement field distribution maps on days 1, 10, 20, 100, 200, and 300 (Figure7 for the M-2K model) were selected. We compared them with the displacement field without considering the rheological properties (Figure 8) and conducted comparative analyses. During the rheological analysis, the horizontal displacement changed the most during the first 20 days. On the first day, the maximum horizontal displacement was 24 mm. On the 10th day, it was 150 mm. On the 20th day, it was 160 mm. After 20 days, the variation range became smaller, and, by the 100th day, it was 160 mm. Notably, the 160 mm displacement range was larger than it was 20 days previously. It was 170 mm on day 200 and 170 mm on day 300. The range was slightly larger than the range at 200 days, but the change of the displacement isoline was small, and there was nearly no change. The results calculated by the Mohr–Coulomb model were obviously smaller than those of the rheological calculation, and the maximum horizontal displacement was only 80 mm. It can be seen that, if the special rheological effects of the soil body are not taken into account in the calculation, the resulting displacement field will be significantly smaller than the actual situation. This result also shows the correctness and accuracy of the secondary development of the model.Figure 7
Horizontal displacement contour map (extended M-2K model). (a) Day 1. (b) Day 10. (c) Day 20. (d) Day 100. (e) Day 200. (f) Day 300.
(a)
(b)
(c)
(d)
(e)
(f)Figure 8
Elastoplastic calculation of horizontal displacement contour map.Ten key points were selected on the landslide. Five were on the slip body, and the remaining five points were on the slip surface. The distribution of the strategic points on the landslide is shown in Figure9. During the rheological calculation, the horizontal displacement of these ten points is shown as a function of time (Figure 10). The creep curves of the key points on the sliding zone and the sliding mass were comparable to those obtained by the triaxial creep test. The creep curves of the 10 key points experienced two stages, including attenuation. During the creep and isotropic creep phases, the decay creep phase started at the beginning of the 20th day. The creep curve has no accelerated creep stage, so the landslide is stable.Figure 9
Schematic of the distribution of key points in the landslide.Figure 10
Curve of the horizontal displacement of the key points over time in (a) the slip surface and (b) the slip body.
(a)
(b)


## 5.1. Landslide Overview

The Qingjiang Shuibuya Water Control Project is located in Badong County, Hubei Province, and is an important water conservancy project in the Qingjiang cascade development. The Dayanyu landslide body is located on the left bank of the Qingjiang River at the lower reaches of the dam site. The landslide is 1.1 km from the dam axis. It was originally a microsand bedrock landslide, and it evolved into a main landslide body and with similar bodies on both sides via local disintegration. The complex main sliding body is elongated, and the upper sliding body is slow and gentle. The secondary sliding bodies on both sides have fan shapes. The total area of the landslide is 0.196 km2, the thickness of the sliding body is 25–40 m, the thickest part is 64.8 m, and the total volume is ∼5.88 × 106 m3.Because the landslide is in the area along the left bank’s antiscouring wall of the dam spillway scouring pit, its stability is directly affected by construction and spillway flooding and fogging, which is a great threat to the safe operation of the engineering building. Thus, it has always been the focus of attention.


## 5.2. Model Construction

The landslide consists of three parts: the slip mass, the slip surface, and the slip bed. The landslide material mesh division is shown in Figure6. In the figure, red indicates the slip mass, green indicates the slip surface, and blue indicates the slip bed. The calculation model divides 1,680 nodes and 811 units. The direction along the landslide is the x-axis, the height direction is the y-axis, and the unit thickness is along the z-axis. The landslide is 950 m along the x direction and 450 m along the y direction.Figure 6
Zoning map of landslide material.


## 5.3. Material Calculation Parameter

In the rheological calculation, the improved Burgers viscoelastic model is used for slip mass, the slip surface, and the slip bed. The material calculation parameters are obtained by analogy and trial calculation according to the fitting parameters of the triaxial creep test data of a landslide, as shown in Table3.Table 3
Improved parameter table of the Burgers visco-elastoplastic model.
Materials
Shear modulusG0 (kPa)
Shear modulusG1 (kPa)
Shear modulusG2 (kPa)
Viscosity coefficientη0 (kPa)
Viscosity coefficientη1 (kPa)
Viscosity coefficientη2 (kPa)
Cohesionc (kPa)
Internal friction angleφ (°)
Elastic modulusE (kPa)
Unit weightγ (kN/m3)
Tensile strengthσt (kPa)
Slip mass
1.659E + 05
8.721E + 04
7.666E + 04
6.063E + 09
3.059E + 06
3.268E + 06
18
25
7.745E + 05
28
0
Slip surface
9.221E + 03
4.845E + 03
4.259E + 03
3.368E + 08
1.699E + 05
1.815E + 05
14
22
4.303E + 04
20
0
Slip bed
1.659E + 07
8.721E + 06
7.666E + 06
6.063E + 11
3.059E + 08
3.268E + 08
32
580
7.745E + 07
40
3500


## 5.4. Boundary Conditions and Load

The model boundary constraint requires the front and rear sides of the landslide to be fixed in thex and z directions. The two sides perpendicular to the z-axis are fixed in the x and z directions. The nodes on the four faces can be deformed along the y direction, but there is no deformation in other directions. Fixed constraints exist in three directions. The duration of the landslide rheology is 300 days. Only the dead-weight load is considered in the calculation. When the rheological calculation is performed, the initial stress state of the landslide is calculated first. Then, the rheological calculations have been started.


## 5.5. Calculation Results

The results of the landslide rheological calculation show displacement and stress field distributions. Only the displacement field was used as the calculation result. When analysing the rheology, the displacement field distribution maps on days 1, 10, 20, 100, 200, and 300 (Figure7 for the M-2K model) were selected. We compared them with the displacement field without considering the rheological properties (Figure 8) and conducted comparative analyses. During the rheological analysis, the horizontal displacement changed the most during the first 20 days. On the first day, the maximum horizontal displacement was 24 mm. On the 10th day, it was 150 mm. On the 20th day, it was 160 mm. After 20 days, the variation range became smaller, and, by the 100th day, it was 160 mm. Notably, the 160 mm displacement range was larger than it was 20 days previously. It was 170 mm on day 200 and 170 mm on day 300. The range was slightly larger than the range at 200 days, but the change of the displacement isoline was small, and there was nearly no change. The results calculated by the Mohr–Coulomb model were obviously smaller than those of the rheological calculation, and the maximum horizontal displacement was only 80 mm. It can be seen that, if the special rheological effects of the soil body are not taken into account in the calculation, the resulting displacement field will be significantly smaller than the actual situation. This result also shows the correctness and accuracy of the secondary development of the model.Figure 7
Horizontal displacement contour map (extended M-2K model). (a) Day 1. (b) Day 10. (c) Day 20. (d) Day 100. (e) Day 200. (f) Day 300.
(a)
(b)
(c)
(d)
(e)
(f)Figure 8
Elastoplastic calculation of horizontal displacement contour map.Ten key points were selected on the landslide. Five were on the slip body, and the remaining five points were on the slip surface. The distribution of the strategic points on the landslide is shown in Figure9. During the rheological calculation, the horizontal displacement of these ten points is shown as a function of time (Figure 10). The creep curves of the key points on the sliding zone and the sliding mass were comparable to those obtained by the triaxial creep test. The creep curves of the 10 key points experienced two stages, including attenuation. During the creep and isotropic creep phases, the decay creep phase started at the beginning of the 20th day. The creep curve has no accelerated creep stage, so the landslide is stable.Figure 9
Schematic of the distribution of key points in the landslide.Figure 10
Curve of the horizontal displacement of the key points over time in (a) the slip surface and (b) the slip body.
(a)
(b)


## 6. Analysis of Landslide Stability Based on the Rheological Calculation

The rigid body limit equilibrium method is the conventional method for determining the stability factor of a landslide, which is a simplified calculation method. However, it cannot perform stress-strain analysis. To carry out stress and strain analyses and determine the stability coefficient of the landslide, the method of determining the stability coefficient of the landslide by the strength reduction method is applied in the numerical analysis method. However, the strength reduction method is only used for elastoplastic analysis. The strength reduction method in rheological analysis is used to determine the stability coefficient of the landslide.The so-called strength reduction results from gradually narrowing the shear strength parameters (i.e., cohesion and internal friction angle) of the slope rock and soil in the ideal elastoplastic finite difference calculation until it reaches its ultimate failure state. Then, the program automatically measures the elastoplasticity. The numerical results demonstrated that the sliding surface of the slope was destroyed and that the safety factor of the strength of the slope was obtained. The sliding surface was a plastic strain shear zone where the elastic strain and displacement were abrupt.The strength reduction stability factor is expressed asf=τ/τ′, where τ is the initial shear strength of the rock and soil material and τ′ is the shear strength when the slope is brought to the limit state after reduction. The finite-difference strength reduction method uses the Mohr–Coulomb strength yield criterion, in which the strength reduction process is described by(34)τ′=τf=c+σtgφf=cf+σtgφf=c′+σtgφ′.Then,(35)c′=cf,tgφ′=tgφf.The finite difference strength reduction method reduces the values ofc and φ according to equation (35). This strength reduction form is in agreement with the definition of the stability coefficient of the traditional limit equilibrium strip method for slope stability analysis. Only the limit equilibrium method of traditional slope stability analysis presupposes a sliding surface and calculates the stability coefficient according to the balance of force (moment). The stability coefficient is described as the sliding resistance (moment) and sliding force (moment) of the sliding surface ratio.The stability coefficient of a landslide was calculated by the finite difference strength reduction method to be 1.24, and the calculation result is shown in Figure11. The slip surface of the landslide changed abruptly along the shear-strain rate. The slip surface calculated by the strength reduction method was roughly consistent with the slip surface of the actual landslide. At the leading edge of the landslide, the shear strain rate reached the maximum value. The stability coefficient is shown in the figure, where the factor-of-safety value was 1.24.Figure 11
Calculation results of elastoplastic strength reduction method.Extant research has reflected many achievements in the calculation principle of the elastoplastic strength reduction method, the determination of slip surface, the criterion of instability, and the improvement of calculation precision. These methods have been widely used in engineering practice. The cohesive force,C, and the internal friction in the improved Burgers viscoelastic model were correspondingly reduced, and the creep curves of the key points on each of the reduction factors are shown in Figure 12. According to the stage experienced by the creep curve, it was determined whether or not the landslide was unstable according to the deformation rate of the constant velocity creep stage. However, it is difficult to unify the critical deformation rate of different landslides, and the critical failure deformation rate of this landslide was not clear. Thus, this paper did not rely on this method. Instead, based on the fact that the key points on the slide zone did not experience a stage of decay creep, they instead directly experienced the stage of constant speed creep or accelerated creep to determine stability.Figure 12
Creep displacement diagram of key points of slip-surface soil with different reduction factors: (a) 1.02; (b) 1.09; (c) 1.13; (d) 1.14; (e) 1.15; (g) 1.16. (f) Creep curve for the first 20 days when the strength reduction factor = 1.15. (h) Creep curve for the first 20 days when the strength reduction factor = 1.16.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)It can be seen from Figure12 that, when the strength reduction factor was between 1.02 and 1.14, the creep curves of the five key points on the slip zone were similar to the creep curves of the triaxial creep test of the landslide slip zone. When the reduction factor was 1.15 and 1.16, the creep curve of the key points on the sliding zone showed no decay creep phase.Under different reduction factors, the maximum displacement and maximum deformation rates are shown in Table4 to have occurred at strategic point 5. When the reduction factor was 1.15, the creep curve had no decay crawl stage. The first 20 days reflected a constant velocity creep curve, and the deformation rate reached 34.8 mm/d, which is several hundred times the stable creep. On the 20th day, the maximum horizontal displacement was 695.1 mm. On the 300th day, the landslide was destroyed when the maximum horizontal displacement reached 6,500 mm. When the reduction factor was 1.16, the first 20 days reflected the accelerated creep phase, and the landslide was obviously destroyed when the maximum deformation rate reached 275.2 mm/d. When the reduction factor was 1.14, the landslide was stable and in a critical equilibrium state.Table 4
Maximum displacement and maximum deformation rates under different reduction factors.
Reduction factor
Creep curve
Maximum horizontal displacement (elastoplasticity)
Maximum horizontal displacement (viscoelasticity)
Maximum deformation rate
Day 20
Day 300
1.02
Attenuation, constant creep
154.2
179.0
110.0
0.0887
1.09
Attenuation, constant creep
199.2
231.8
138.6
0.116
1.13
Attenuation, constant creep
281.9
347.0
175.5
0.232
1.14
Attenuation, constant creep
347.7
449.0
226.4
0.362
1.15
Constant creep
695.1
6500.0
276.1
34.8
1.16
Accelerated creep
1522.9
67800.0
319.5
275.2
Displacement unit: mm. Deformation rate unit: mm/d.The horizontal displacement and deformation rates of fundamental point 5, having the largest horizontal displacement on the sliding belt, were selected as the research object. The variation law of the displacement and deformation rates having the reduction factor was examined. As shown in Figures13 and 14, the reduction factor was 1.14. The inflection point at which the displacement and deformation rate was abrupt, in combination with the above, and the stability coefficient of the landslide was 1.14 when considering rheological properties.Figure 13
Horizontal displacements of the 300th day with a key point of 5 under different reduction factors.Figure 14
Horizontal deformation rate at the 300th day with a key point of 5 under different reduction factors.


## 7. Conclusions

Based on the triaxial creep test data of a landslide slip-surface soil, the rheological characteristics of the landslide were analysed. The FLAC3D user-defined constitutive model was written with C++, and the improved model was applied to the landslide rheological analysis. This paper proposed a 3D creep equation derivation method suitable for conventional triaxial tests. For the creep test data and creep curve of a landslide slip zone, a 3D rheological model suitable for the rheological properties of the landslide slip surface was constructed, and the corresponding parameters were determined.Using the FLAC3D secondary development function, based on the extended Burgers rheological model, the FLAC3D user-defined constitutive model was compiled. Distribution analysis considered rheological effects, and the practical deformation of the landslide was considered without rheological effects. The results showed that the landslide displacement calculation was significantly smaller when the rheological effects were not considered.The method of determining the stability coefficient using the elastoplastic strength reduction method was applied to the rheological calculation. Cohesionc and internal friction angle φ in the rheological parameters were reduced by this method. Then, the landslide was used in the case of each reduction factor. The rheological analysis was then carried out. When considering the rheological properties, the landslide stability coefficient was 1.14, and the elastoplastic strength reduction method obtained a stability coefficient of 1.24. Considering that the rheological stability coefficient was smaller than when the rheology was not considered, it indicated geotechnical material. The rheological properties had an adverse effect on the stability of the slope. The rheological properties of the landslide should therefore be considered when determining the stability factor of the landslide.



---

*Source: 2694651-2020-01-10.xml*