# New Maxwell Creep Model Based on Fractional and Elastic-Plastic Elements

Authors: Hao Tang; Dongpo Wang; Zhao Duan
Journal: Advances in Civil Engineering
 (2020)

---

## Abstract

Creep models are mainly used to describe the rheological behaviour of geotechnical materials. An important research focus for studying creep in geotechnical materials is the development of a model with few parameters and good simulation performance. Hence, in this study, by replacing the Newtonian dashpot and spring in the classical Maxwell model with fractional and elastic-plastic elements, a new Maxwell creep model based on fractional derivatives and continuum damage mechanics was developed. One- and three-dimensional (1D/3D) creep equations of the new Maxwell creep model were derived. The 1D creep equation of the new model was used to fit existing creep data of rock salt, and the 3D creep equation was used to fit the creep data of remolded loess. The model curves matched the creep data very well, showing considerably higher accuracy than other models. Furthermore, a sensitivity study was carried out, showing the effects of the fractional derivative orderβ and exponent α on the creep strain of rock salt. This new model is simple with few parameters and can effectively simulate the complete creep behaviour of geotechnical materials.


---

## Body

## 1. Introduction

Creep models are abstract methods for expressing the creep characteristics of geotechnical materials. The creep behaviour of geotechnical materials can be described using differential, integral, or empirical creep equations, and their creep characteristics can be represented using a theoretical creep curve approximating the real creep curve. As a differential creep model, the component model has been applied and expanded by many researchers because of its simple structure and clear physical meaning (e.g., Maxwell, Kelvin, Burgers, and Nishihara models). It is well known that the deformation of geotechnical materials during creep tests can be divided into four stages: instantaneous deformation, transient creep, steady creep, and accelerated creep. These different deformation stages have their own characteristics, and it is difficult to build a model that can simulate all four deformation stages. Therefore, many scholars are looking for new methods to improve the component model for representing all the deformation stages.Fractional calculus involves derivative and integration operators that are applied to fractional orders. As one of the generalisations of classical calculus, fractional calculus is used widely in various fields of science and engineering. Because of the long history dependence, or the so-called memory effect, fractional operators are powerful tools for modelling the viscoelastic behaviour of materials, particularly for building a time-dependent constitutive model. Because of the pioneering work of Scott-Blair [1], who proposed a fractional element analogous to the classical Newtonian dashpot, many creep models based on fractional calculus have been developed in recent years. Rogers [2], Bagly and Torvik [3, 4], and Koeller [5] have performed remarkable studies on creep models using fractional calculus. Welch et al. [6] proposed a four-parameter creep model to characterise the viscoelastic creep of polymeric materials using fractional calculus. Zhou et al. [7, 8] proposed a fractional derivative model to describe the creep behaviour of rock salt. Yin et al. [9] established a variable-order fractional creep model to describe the evolution of the mechanical properties of soft soils. Kang et al. [10] applied fractional calculus to build a nonlinear creep model for coal. Celauro et al. [11] proposed a fractional derivative model to describe the viscoelastic behaviour of asphalt mixtures. In addition, some scholars [12–17] have developed relevant creep models using fractional calculus. Although fractional calculus is widely adopted in developing creep models, a key limitation is its inability to describe the accelerated creep phase of geotechnical materials; in terms of the mathematical description of the obtained curves, fractional calculus cannot suitably describe the accelerated creep stage of rocks.Therefore, using a continuum damage mechanics approach, some creep models incorporate a damage factor to describe the accelerated creep phase. Moreover, the method of introducing a damage factor provides a clear physical meaning in characterising creep damage. Mazotti and Savoia [18] proposed a creep damage model for concrete, considering the combined effect of nonlinear viscous strain evolution and crack nucleation and propagation at high stress levels. Verstrynge et al. [19] proposed an improved Burges model by introducing a viscous damage parameter, which reflected the accumulation of damage over a long period in a masonry structure. Ma et al. [20] proposed a three-dimensional (3D) creep damage model combined with the generalised Hoek–Brown model and validated the same with a series of strength and creep tests for bedded rock salt. Cao et al. [21] proposed an elastic-plastic damage cell model by introducing the damage theory and law of damage proposed by Kachanov [22], which can describe the nonlinear or accelerated rock creep process. In addition, many scholars [23–26] have developed creep models for rock by applying continuum damage mechanics. However, despite recent improvements in the creep models for geotechnical materials, there is no simple or universal model that can fully describe their viscoelastic-plastic creep behaviour. Models with simple structures cannot effectively simulate the complete creep behaviour of geotechnical materials; moreover, many of the effective models have more parameters. Therefore, developing a model with few parameters and better simulation performance is the key goal in the study of creep constitutive models.The Maxwell model is the simplest creep model with only two elements and two parameters. However, its mathematical equation can only express linear material creep. Therefore, the original Maxwell model requires improvements. In recent years, some scholars [27–32] have applied fractional calculus to improve the original Maxwell model. These improved Maxwell models can approximately describe the partial nonlinear creep deformation of geotechnical materials. However, these improved Maxwell models are unable to describe the accelerated creep phase effectively.Therefore, in this paper, we present a new Maxwell creep model based on fractional and elastic-plastic elements to characterise the creep behaviour of rock salt and remolded loess under different stress states. The advantages of the model were identified through validation and sensitivity analysis.


## 2. New Maxwell Creep Model

### 2.1. Fractional Element

From a mathematical perspective, the constitutive equations of spring and dashpot elements can be described as follows:(1)σt=Ed0εtdt0,σt=ηd1εtdt1.As shown in equation (1), the stress of an elastic solid is proportional to the zeroth-order derivative of the strain, whereas the stress of a Newtonian fluid is proportional to the first-order derivative of the strain. Based on the mathematical transformation relations mentioned above, a general constitutive equation for a viscoelastic body may be written as(2)σt=ξdβεtdtβ,0≤β≤1.Equation (2) is known as the Scott-Blair [1] fractional element, and it is commonly referred to as the Abel dashpot. The symbolic representation of this element is shown in Figure 1. In equation (2), σ(t) and ε(t) are the stress and strain, respectively; ξ is the viscosity coefficient of the material, with units of (stress·time)β; t is the time; β is the fractional order, with 0 ≤ β ≤ 0. Equation (2) degrades to show the constitutive relationship of a spring when β = 0 and that of a dashpot when β = 1. However, when 0 < β < 1, the fractional element can exhibit the characteristics of both a spring and dashpot.Figure 1
Fractional element.Fractional calculus can be defined in many ways, such as Riemann–Liouville and Caputo derivatives [33]. The Riemann–Liouville fractional integration of function f(t) with order β is defined as(3)Dt−β0ft=d−βftdt−β=∫t0tt−τβ−1Γβfτdτ,and the fractional derivative of function f(t) with order β is defined as(4)Dtβ0ft=dβftdtβ=dndtnDt−n−βt0ft,where β > 0 and n − 1 < β ≤ n (n is a positive integer), and Γ(x) is the gamma function defined as(5)Γx=∫0∞e−ttx−1dt,Rex>0.If equation (2) is considered in terms of its physical interpretation, it can be used to characterise the mechanical properties of geotechnical materials ranging between an ideal solid and ideal fluid during creep. Thus, if the stress σ(t) is a constant, the fractional element will describe the change in the creep behaviour of the geotechnical materials. The creep equation of materials can be obtained by employing the Riemann–Liouville fractional calculus theory as follows:(6)εt=σ0ξ⋅tβΓ1+β,0<β<1,where σ0 denotes the constant stress level. Two parameters are included in equation (6), namely, the viscosity coefficient ξ and fractional order β. ξ Reflects viscosity of the materials, and β reflects the variation laws of the creep curves of materials. Substituting σ(t) = 10 MPa and ξ = 2000 MPa·hβ into equation (6), a series of creep curves with different β values were obtained, as shown in Figure 2.Figure 2
Fractional element creep curves for a range of fractional orders.As observed in Figure2, when the stress remains constant, the creep strain increases gradually as β increases. When 0 < β < 1, the morphological characteristics of the creep curves demonstrate that the fractional element can be used to describe the transient and steady creep phase of rock.


### 2.2. Elastic-Plastic Element

The deformation of geotechnical materials during a creep test can be divided into four phases: instantaneous deformation, transient creep, steady creep, and accelerated creep. Based on previous research experience, many classical component models can simulate the first three deformation stages reliably. However, many of the proposed models have limitations in terms of simulating the accelerated creep stage. Therefore, it is always challenging to simulate the accelerated creep stage; many researchers have tried to address this problem by proposing various nonlinear elements. Through comparison and analysis, a type of elastic-plastic element based on the damage rate change proposed by Cao et al. [21] was identified for succinct and accurate simulation of the accelerated creep stage in this study. Based on the damage variable proposed by Kachanov [22], considering the critical stress (σs) of the geotechnical material during creep, the new damage variable improved by Cao et al. [21] can be written as(7)D=0,σ<σS,1−1−ttFα,σ≥σS,where tF is the time from initial deformation to final failure and α is the material constant. Substituting equation (7) into the creep constitutive relation of the elastic element, the creep equation of the elastic-plastic element can be expressed as(8)ε=σE,σ<σS,σE1−ttF−α,σ≥σS.The constitutive equation of the elastic-plastic element is effective to a certain extent for simulating the creep characteristics of geotechnical materials. Whenσ < σs, the elastic-plastic element degenerates into a spring element, and when σ ≥ σs and the time is tF, the element can describe the accelerated creep phases of geotechnical materials.


### 2.3. New Maxwell Creep Model

As mentioned previously, the deformation rock during a creep test can be divided into four phases: instantaneous deformation, transient creep, steady creep, and accelerated creep. Therefore, an excellent creep model must have the ability to simulate the intact creep behaviour of rock, especially in the accelerated creep phase. To achieve this objective, a new creep model based on a fractional element and an elastic-plastic element placed in series was developed, as shown in Figure3. This model is referred to as the new Maxwell creep model.Figure 3
New Maxwell creep model.Based on the connection characteristics of the elements, in the one-dimensional (1D) creep state, this new Maxwell creep model satisfies the following relationship:(9)ε=σE+σξ⋅tβΓ1+β,σ<σs,0<β<1,ε=σE1−t/tFα+σξ⋅tβΓ1+β,σ≥σs,0<β<1,where E and ξ are the elastic modulus and viscosity coefficient of loess, respectively; α is a constant pertaining to the properties of the rock; and β is the fractional order.Equation (9) shows that the new Maxwell creep model will degenerate into spring and fractional elements in series when the loading stress is lower than the critical stress. Thus, the new model can be used to describe instantaneous deformation, transient creep, and steady creep of rock. When the loading stress is larger than the critical stress, this model can be used to describe the entire deformation process, including accelerated creep. This model has the advantage of using only four parameters, thus resulting in a simple structure.In the 3D stress state, the stress tensorσij can be decomposed into the stress sphere tensor σm and partial tensor Sij. Similarly, the strain tensor εij can be decomposed into the strain sphere tensor εm and partial tensor eij. Their relationship can be expressed as follows:(10)σij=Sij+δijσm,εij=eij+δijεm,where δij is the Kronecker symbol. According to the generalised Hook’s law, the 3D constitutive equation for a Hookean body can be written as(11)εije=Sij2G+σm3Kδij,where εije, G, and K represent the strain tensor, shear modulus, and volume modulus of the Hookean body, respectively. Assuming that the material rheology is mainly governed by shear deformation, the 3D constitutive relationship of a viscoelastic body (fractional element) can be expressed as(12)eijve=Sij2ξ′⋅tβΓ1+β,where eijve and ξ′ are the partial strain tensor and shear modulus of the viscoelastic body, respectively.Based on equations (11) and (12), the 3D creep equation of the new Maxwell creep model can be written as(13)εij=Sij2G+σm3Kδij+Sij2ξ′⋅tβΓ1+β,Sij<σs,0<β<1,εij=Sij2G∗+σm3Kδij+Sij2ξ′⋅tβΓ1+β,Sij≥σs,0<β<1,where G∗=G1−D.When the triaxial compression creep stress state isσ2 = σ3 and constant, with all levels of load σ1 being constant, then(14)σm=13σ1+2σ3,S11=σ1−σm=23σ1−σ3.and thus equation (13) can be rewritten as(15)εij=σ1−σ33G+σ1+2σ39K+σ1−σ33ξ′⋅tβΓ1+β,σ1−σ3<σs,0<β<1,εij=σ1−σ33G1−t/tFα+σ1+2σ39K+σ1−σ33ξ′⋅tβΓ1+β,σ1−σ3≥σs,0<β<1.


## 2.1. Fractional Element

From a mathematical perspective, the constitutive equations of spring and dashpot elements can be described as follows:(1)σt=Ed0εtdt0,σt=ηd1εtdt1.As shown in equation (1), the stress of an elastic solid is proportional to the zeroth-order derivative of the strain, whereas the stress of a Newtonian fluid is proportional to the first-order derivative of the strain. Based on the mathematical transformation relations mentioned above, a general constitutive equation for a viscoelastic body may be written as(2)σt=ξdβεtdtβ,0≤β≤1.Equation (2) is known as the Scott-Blair [1] fractional element, and it is commonly referred to as the Abel dashpot. The symbolic representation of this element is shown in Figure 1. In equation (2), σ(t) and ε(t) are the stress and strain, respectively; ξ is the viscosity coefficient of the material, with units of (stress·time)β; t is the time; β is the fractional order, with 0 ≤ β ≤ 0. Equation (2) degrades to show the constitutive relationship of a spring when β = 0 and that of a dashpot when β = 1. However, when 0 < β < 1, the fractional element can exhibit the characteristics of both a spring and dashpot.Figure 1
Fractional element.Fractional calculus can be defined in many ways, such as Riemann–Liouville and Caputo derivatives [33]. The Riemann–Liouville fractional integration of function f(t) with order β is defined as(3)Dt−β0ft=d−βftdt−β=∫t0tt−τβ−1Γβfτdτ,and the fractional derivative of function f(t) with order β is defined as(4)Dtβ0ft=dβftdtβ=dndtnDt−n−βt0ft,where β > 0 and n − 1 < β ≤ n (n is a positive integer), and Γ(x) is the gamma function defined as(5)Γx=∫0∞e−ttx−1dt,Rex>0.If equation (2) is considered in terms of its physical interpretation, it can be used to characterise the mechanical properties of geotechnical materials ranging between an ideal solid and ideal fluid during creep. Thus, if the stress σ(t) is a constant, the fractional element will describe the change in the creep behaviour of the geotechnical materials. The creep equation of materials can be obtained by employing the Riemann–Liouville fractional calculus theory as follows:(6)εt=σ0ξ⋅tβΓ1+β,0<β<1,where σ0 denotes the constant stress level. Two parameters are included in equation (6), namely, the viscosity coefficient ξ and fractional order β. ξ Reflects viscosity of the materials, and β reflects the variation laws of the creep curves of materials. Substituting σ(t) = 10 MPa and ξ = 2000 MPa·hβ into equation (6), a series of creep curves with different β values were obtained, as shown in Figure 2.Figure 2
Fractional element creep curves for a range of fractional orders.As observed in Figure2, when the stress remains constant, the creep strain increases gradually as β increases. When 0 < β < 1, the morphological characteristics of the creep curves demonstrate that the fractional element can be used to describe the transient and steady creep phase of rock.


## 2.2. Elastic-Plastic Element

The deformation of geotechnical materials during a creep test can be divided into four phases: instantaneous deformation, transient creep, steady creep, and accelerated creep. Based on previous research experience, many classical component models can simulate the first three deformation stages reliably. However, many of the proposed models have limitations in terms of simulating the accelerated creep stage. Therefore, it is always challenging to simulate the accelerated creep stage; many researchers have tried to address this problem by proposing various nonlinear elements. Through comparison and analysis, a type of elastic-plastic element based on the damage rate change proposed by Cao et al. [21] was identified for succinct and accurate simulation of the accelerated creep stage in this study. Based on the damage variable proposed by Kachanov [22], considering the critical stress (σs) of the geotechnical material during creep, the new damage variable improved by Cao et al. [21] can be written as(7)D=0,σ<σS,1−1−ttFα,σ≥σS,where tF is the time from initial deformation to final failure and α is the material constant. Substituting equation (7) into the creep constitutive relation of the elastic element, the creep equation of the elastic-plastic element can be expressed as(8)ε=σE,σ<σS,σE1−ttF−α,σ≥σS.The constitutive equation of the elastic-plastic element is effective to a certain extent for simulating the creep characteristics of geotechnical materials. Whenσ < σs, the elastic-plastic element degenerates into a spring element, and when σ ≥ σs and the time is tF, the element can describe the accelerated creep phases of geotechnical materials.


## 2.3. New Maxwell Creep Model

As mentioned previously, the deformation rock during a creep test can be divided into four phases: instantaneous deformation, transient creep, steady creep, and accelerated creep. Therefore, an excellent creep model must have the ability to simulate the intact creep behaviour of rock, especially in the accelerated creep phase. To achieve this objective, a new creep model based on a fractional element and an elastic-plastic element placed in series was developed, as shown in Figure3. This model is referred to as the new Maxwell creep model.Figure 3
New Maxwell creep model.Based on the connection characteristics of the elements, in the one-dimensional (1D) creep state, this new Maxwell creep model satisfies the following relationship:(9)ε=σE+σξ⋅tβΓ1+β,σ<σs,0<β<1,ε=σE1−t/tFα+σξ⋅tβΓ1+β,σ≥σs,0<β<1,where E and ξ are the elastic modulus and viscosity coefficient of loess, respectively; α is a constant pertaining to the properties of the rock; and β is the fractional order.Equation (9) shows that the new Maxwell creep model will degenerate into spring and fractional elements in series when the loading stress is lower than the critical stress. Thus, the new model can be used to describe instantaneous deformation, transient creep, and steady creep of rock. When the loading stress is larger than the critical stress, this model can be used to describe the entire deformation process, including accelerated creep. This model has the advantage of using only four parameters, thus resulting in a simple structure.In the 3D stress state, the stress tensorσij can be decomposed into the stress sphere tensor σm and partial tensor Sij. Similarly, the strain tensor εij can be decomposed into the strain sphere tensor εm and partial tensor eij. Their relationship can be expressed as follows:(10)σij=Sij+δijσm,εij=eij+δijεm,where δij is the Kronecker symbol. According to the generalised Hook’s law, the 3D constitutive equation for a Hookean body can be written as(11)εije=Sij2G+σm3Kδij,where εije, G, and K represent the strain tensor, shear modulus, and volume modulus of the Hookean body, respectively. Assuming that the material rheology is mainly governed by shear deformation, the 3D constitutive relationship of a viscoelastic body (fractional element) can be expressed as(12)eijve=Sij2ξ′⋅tβΓ1+β,where eijve and ξ′ are the partial strain tensor and shear modulus of the viscoelastic body, respectively.Based on equations (11) and (12), the 3D creep equation of the new Maxwell creep model can be written as(13)εij=Sij2G+σm3Kδij+Sij2ξ′⋅tβΓ1+β,Sij<σs,0<β<1,εij=Sij2G∗+σm3Kδij+Sij2ξ′⋅tβΓ1+β,Sij≥σs,0<β<1,where G∗=G1−D.When the triaxial compression creep stress state isσ2 = σ3 and constant, with all levels of load σ1 being constant, then(14)σm=13σ1+2σ3,S11=σ1−σm=23σ1−σ3.and thus equation (13) can be rewritten as(15)εij=σ1−σ33G+σ1+2σ39K+σ1−σ33ξ′⋅tβΓ1+β,σ1−σ3<σs,0<β<1,εij=σ1−σ33G1−t/tFα+σ1+2σ39K+σ1−σ33ξ′⋅tβΓ1+β,σ1−σ3≥σs,0<β<1.


## 3. Analysis of New Maxwell Creep Model in 1D Creep State

### 3.1. Data Preparation

Based on previous research results [34], typical axial creep experimental data for rock salt in the Changshan Salt Mine (Zigong City, Sichuan Province, China) were used to verify the proposed new Maxwell creep model. Figure 4 shows the axial creep curves for rock salt under loading stresses of 9.93, 14.41, and 14.72 MPa. The creep curve for a loading stress of 14.41 MPa exhibits the three typical creep phases and thus represents the typical intact creep behaviour of rock. The creep curve for a loading stress of 9.93 MPa is characteristic of the case when the loading stress is lower than the yield stress and that for a loading stress of 14.72 MPa is characteristic of the case when the loading stress is larger than the yield stress.Figure 4
Axial creep curves for Changshan rock salt.


### 3.2. Parameter Determination

We used the 1D creep equation of the new Maxwell model to fit the experimental data for rock salt shown above. The parameters in Equation (9) can be determined using the Levenberg–Marquardt method, which is a nonlinear least-squares fitting method, and the results are presented in Table 1. The efficacy of the proposed creep model in fitting the rock salt data is shown in Figure 5. To further analyse the simulation capability of the proposed creep model, we compared it to the Burgers and Nishihara models, as shown in Figure 5; the parameters of the Burgers and Nishihara models are also given in Table 1.Table 1
Parameters determined through fitting analysis based on creep tests of rock salt.
New Maxwell creep model
σ (MPa)
E (GPa)
ξ (GPa·hβ)
β
α
σ
s (MPa)
t
F (h)
R
2
9.93
0.875
4.544
0.260
13.5
0.918
14.41
0.590
3.662
0.252
0.191
13.5
1056
0.988
14.72
0.415
3.904
0.241
0.110
13.5
385
0.978
Burgers model
σ (MPa)
E
H (GPa)
E
K (GPa)
η
K (GPa·h)
η (GPa·h)
R
2
9.93
0.860
0.931
3446.914
63.999
0.993
14.41
0.568
1.271
541.794
30.426
0.968
14.72
0.410
2.222
18.242
221.010
0.901
Nishihara model
σ (MPa)
E
H (GPa)
E
K (GPa)
η
K (GPa·h)
η
S (GPa·h)
σ
s (MPa)
R
2
9.93
0.829
0.820
79.527
13.5
0.975
14.41
0.568
1.271
30.425
34.214
13.5
0.968
14.72
0.410
2.222
18.242
18.317
13.5
0.901Figure 5
Comparison between theoretical and experimental curves of (a) new Maxwell model and Burgers model and (b) new Maxwell model and Nishihara model.
(a)
(b)From the analysis results presented in Table1, it is seen that the coefficient of determination R2 > 0.9, indicating that the new Maxwell creep model proposed herein is in good agreement with the experimental data of rock salt. Moreover, Figure 5 shows that the new Maxwell creep model has good consistency with the experimental data, which further proves that the proposed model is effective and reliable. In addition, Figure 5 shows that the new Maxwell creep model provides better simulation of the creep properties of rock salt than the Burgers and Nishihara models when the loading stress is larger than the critical stress, especially for the accelerated creep phase; further, when the loading stress is lower than the critical stress, the new Maxwell creep model can also show good simulation ability, like that of the Burgers model. Overall, the new Maxwell creep model proposed in this work is capable of describing the creep properties of rock, and its structure is very simple, containing only four parameters; therefore, it has high potential for practical applications in rock engineering.


### 3.3. Sensitivity Studies Using the New Maxwell Creep Model

It is clear that the new Maxwell creep model is advantageous in simulating the accelerated creep stage of rock salt under a 1D stress state. These advantages depend on just a few parameters used in this model, as shown in Equation (9), among which the fractional derivative order β and exponent α are the most important ones. To obtain a better understanding regarding the effects of these parameters, sensitivity analyses were carried out. Substituting σ = 14.41 MPa, σs = 13.5 MPa, E = 0.590 GPa, ξ = 3.662 GPa·hβ, tF = 1056 h, and α = 0.191 (data from Table 1) into equation (9) and letting the fractional derivative order β change from 0.2 to 0.3 with an interval of 0.02, we obtained a series of curves, as shown in Figure 6; these indicate that the creep strain and creep rate in the steady creep stage mainly depend on the fractional derivative order. Small increments in the fractional derivative order can cause the creep curve to rise in the steady creep stage, thus leading to a greater steady creep strain and higher steady creep strain rate.Figure 6
Sensitivity of the creep strain to fractional derivative orderβ (σ = 14.41 MPa, σs = 13.5 MPa, E = 0.590 GPa, ξ = 3.662 GPa·hβ, tF = 1056 h, and α = 0.191).By changing exponentα and keeping the other parameters constant, a series of curves were obtained (Figure 7). It is clear that the strain in the accelerated creep stage and creep rate increase with exponent α. Changing α values can realise simulation of the accelerated creep stage with different creep rates.Figure 7
Sensitivity of the creep strain to exponentα (σ = 14.41 MPa, σs = 13.5 MPa, E = 0.590 GPa, ξ = 3.662 GPa·hβ, tF = 1056 h, and β = 0.252).


## 3.1. Data Preparation

Based on previous research results [34], typical axial creep experimental data for rock salt in the Changshan Salt Mine (Zigong City, Sichuan Province, China) were used to verify the proposed new Maxwell creep model. Figure 4 shows the axial creep curves for rock salt under loading stresses of 9.93, 14.41, and 14.72 MPa. The creep curve for a loading stress of 14.41 MPa exhibits the three typical creep phases and thus represents the typical intact creep behaviour of rock. The creep curve for a loading stress of 9.93 MPa is characteristic of the case when the loading stress is lower than the yield stress and that for a loading stress of 14.72 MPa is characteristic of the case when the loading stress is larger than the yield stress.Figure 4
Axial creep curves for Changshan rock salt.


## 3.2. Parameter Determination

We used the 1D creep equation of the new Maxwell model to fit the experimental data for rock salt shown above. The parameters in Equation (9) can be determined using the Levenberg–Marquardt method, which is a nonlinear least-squares fitting method, and the results are presented in Table 1. The efficacy of the proposed creep model in fitting the rock salt data is shown in Figure 5. To further analyse the simulation capability of the proposed creep model, we compared it to the Burgers and Nishihara models, as shown in Figure 5; the parameters of the Burgers and Nishihara models are also given in Table 1.Table 1
Parameters determined through fitting analysis based on creep tests of rock salt.
New Maxwell creep model
σ (MPa)
E (GPa)
ξ (GPa·hβ)
β
α
σ
s (MPa)
t
F (h)
R
2
9.93
0.875
4.544
0.260
13.5
0.918
14.41
0.590
3.662
0.252
0.191
13.5
1056
0.988
14.72
0.415
3.904
0.241
0.110
13.5
385
0.978
Burgers model
σ (MPa)
E
H (GPa)
E
K (GPa)
η
K (GPa·h)
η (GPa·h)
R
2
9.93
0.860
0.931
3446.914
63.999
0.993
14.41
0.568
1.271
541.794
30.426
0.968
14.72
0.410
2.222
18.242
221.010
0.901
Nishihara model
σ (MPa)
E
H (GPa)
E
K (GPa)
η
K (GPa·h)
η
S (GPa·h)
σ
s (MPa)
R
2
9.93
0.829
0.820
79.527
13.5
0.975
14.41
0.568
1.271
30.425
34.214
13.5
0.968
14.72
0.410
2.222
18.242
18.317
13.5
0.901Figure 5
Comparison between theoretical and experimental curves of (a) new Maxwell model and Burgers model and (b) new Maxwell model and Nishihara model.
(a)
(b)From the analysis results presented in Table1, it is seen that the coefficient of determination R2 > 0.9, indicating that the new Maxwell creep model proposed herein is in good agreement with the experimental data of rock salt. Moreover, Figure 5 shows that the new Maxwell creep model has good consistency with the experimental data, which further proves that the proposed model is effective and reliable. In addition, Figure 5 shows that the new Maxwell creep model provides better simulation of the creep properties of rock salt than the Burgers and Nishihara models when the loading stress is larger than the critical stress, especially for the accelerated creep phase; further, when the loading stress is lower than the critical stress, the new Maxwell creep model can also show good simulation ability, like that of the Burgers model. Overall, the new Maxwell creep model proposed in this work is capable of describing the creep properties of rock, and its structure is very simple, containing only four parameters; therefore, it has high potential for practical applications in rock engineering.


## 3.3. Sensitivity Studies Using the New Maxwell Creep Model

It is clear that the new Maxwell creep model is advantageous in simulating the accelerated creep stage of rock salt under a 1D stress state. These advantages depend on just a few parameters used in this model, as shown in Equation (9), among which the fractional derivative order β and exponent α are the most important ones. To obtain a better understanding regarding the effects of these parameters, sensitivity analyses were carried out. Substituting σ = 14.41 MPa, σs = 13.5 MPa, E = 0.590 GPa, ξ = 3.662 GPa·hβ, tF = 1056 h, and α = 0.191 (data from Table 1) into equation (9) and letting the fractional derivative order β change from 0.2 to 0.3 with an interval of 0.02, we obtained a series of curves, as shown in Figure 6; these indicate that the creep strain and creep rate in the steady creep stage mainly depend on the fractional derivative order. Small increments in the fractional derivative order can cause the creep curve to rise in the steady creep stage, thus leading to a greater steady creep strain and higher steady creep strain rate.Figure 6
Sensitivity of the creep strain to fractional derivative orderβ (σ = 14.41 MPa, σs = 13.5 MPa, E = 0.590 GPa, ξ = 3.662 GPa·hβ, tF = 1056 h, and α = 0.191).By changing exponentα and keeping the other parameters constant, a series of curves were obtained (Figure 7). It is clear that the strain in the accelerated creep stage and creep rate increase with exponent α. Changing α values can realise simulation of the accelerated creep stage with different creep rates.Figure 7
Sensitivity of the creep strain to exponentα (σ = 14.41 MPa, σs = 13.5 MPa, E = 0.590 GPa, ξ = 3.662 GPa·hβ, tF = 1056 h, and β = 0.252).


## 4. Analysis of New Maxwell Creep Model in 3D Creep State

### 4.1. Test Conditions and Results

To further explore the ability of the proposed model to simulate the creep behaviour of geotechnical materials under a 3D stress state, the creep results of remolded loess under triaxial conditions were obtained. The loess samples were collected from the L6 loess layer in the new district of Yan’an city; it belongs to the quaternary middle Pleistocene loess, also referred to as theQ2 loess. The undisturbed loess was ground, dried, and sieved through a 2-mm screen to have a moisture content of 10%; it was kept for 24 h without allowing evaporation to allow the water to be uniformly distributed in the loess. By controlling the compactness levels (k) at 0.84, 0.89, 0.95, and 0.99, the amount of loess required for forming remolded loess under different compactness levels was prepared. The prepared loess was filled into a cylindrical mold to prepare remolded loess specimens and then compressed as per set compression rates using a sample-making machine. The prepared remolded loess specimens were cylinders with diameters and heights of 39.8 and 80 mm, respectively.The SLB-1A stress-strain controlled test machine was used for triaxial creep shear measurements, as shown in Figure8. This instrument can be stress-loaded through a pneumatic pump, so the set stress can reach a predetermined value in a short time. Before performing the creep tests of remolded loess, the ultimate deviatoric stress values of remolded loess specimens were determined via triaxial consolidated undrained compression tests at a fixed confining pressure of 100 kPa. Under the same confining pressure, the specimens were first consolidated for 24 h; then, indoor creep tests using remolded loess specimens with different compactness levels were conducted under different axial deviatoric stress levels. The axial deviatoric stress was applied step by step using stress levels of R = 20, 40, 60, 80, 90, 95, and 98% (R denotes the ratio of the deviatoric stress to the ultimate stress), and the axial pressure exerted on the specimens during loading was kept constant with a precision of ±5 kPa for 12 h.Figure 8
Triaxial creep shear test apparatus: (a) instrument details; (b) creep test state.A cylindrical specimen was prepared and wrapped on the outside with a rubber film with a tubular mold. Filter paper and permeable stone were placed at both ends of the sample, and the lower end of the sample rubber band was tied to the base and the upper end to the top. Next, the pressure chamber cover was installed and turned on the water switch to fill the pressure chamber with water. At this point, we opened the drain valve and applied the required confining pressure of 100 kPa through stress loading; this pressure was maintained for 24 h to allow the specimen to consolidate. We observed changes in the pore water pressure during consolidation, and we then closed the drain valve and applied axial deviatoric stresses were applied in a stepwise manner at the foregoing levels and recorded axial deformation data using an external strain gauge during creep.Based on the results of the above creep tests, creep curves at various compactness levels and deviatoric stresses were drawn, as shown in Figure9. It is evident that axial creep deformations under various compactness levels are roughly similar; the primary, secondary, and tertiary creep stages can be observed under partial deviatoric stress. The specimens at each step during application of the deviatoric stress showed the primary creep stage where the strain rate gradually decreased. After the primary creep stage, the specimens immediately entered the longer secondary creep stage; in this stage, the strain rate tended to be zero under low deviatoric stress; however, under high deviatoric stress, the creep strain increased linearly with time when the strain rate was greater than zero and constant. Next, under a particular high deviatoric stress, the tertiary creep stage developed with time following the secondary creep stage and the axial strain and strain rate significantly increased at this stage.Figure 9
Axial strain versus time at different compactness levels: (a)k = 0.84; (b) k = 0.89; (c) k = 0.95; (d) k = 0.99.
(a)
(b)
(c)
(d)


### 4.2. Parameter Determination

Considering the creep test results of remolded loess samples with a compactness of 0.95 as an example, we used the 3D creep equation of the new Maxwell model to fit the experimental data. The parameters in equation (15) can also be determined by the Levenberg–Marquardt method, and the results are presented in Table 2. The efficacy of the proposed creep model in fitting the data of remolded loess is shown in Figure 10. To further analyse the ability of the model to describe the accelerated creep stage of remolded loess under a 3D stress state, we selected the creep data corresponding to the last stress loading and compared it to the Burgers and Nishihara models, as shown in Figure 11; the parameters of the Burgers and Nishihara models are given in Table 2.Table 2
Parameters determined through fitting the analysis based on creep test results of remolded loess.
New Maxwell creep model
σ
1 – σ3 (kPa)
G (kPa)
K (kPa)
ξ′(kPa·hβ)
β
α
σ
s (kPa)
t
F (h)
R
2
122
29856.761
742.362
0.018
0.041
0.974
244
209.192
91.334
0.003
0.097
0.975
366
5832.610
42.649
0.001
0.181
0.954
488
90.682
345.728
0.0008
0.332
0.978
550
2702.403
33.217
0.005
0.456
0.785
490
6.3
0.998
Burgers model
σ
1 – σ3 (kPa)
G
0 (kPa)
K (kPa)
G
k (kPa·h)
η
k (kPa·h)
η (kPa·h)
R
2
550
105.281
48.865
0.235
25516.903
431.671
0.767
Nishihara model
σ
1 – σ3 (kPa)
G
0 (kPa)
K (kPa)
G
k (kPa·h)
η
k (kPa·h)
η
s (kPa·h)
σ
s (kPa)
R
2
550
80.301
71.612
0.155
28047.865
47.019
490
0.797Figure 10
Comparison between theoretical and experimental curves.Figure 11
Comparison between theoretical and experimental curves of (a) new Maxwell model and Burgers model and (b) new Maxwell model and Nishihara model.
(a)
(b)As observed in Figure10 and Table 2, the calculated results using the new Maxwell model are in good agreement with the test results obtained under deviatoric stresses of 122, 244, 366, 488, and 550 kPa. The correlation coefficient squares (R2) of the five curves are 0.974, 0.975, 0.954, 0.978, and 0.996, respectively, indicating that the prediction precision of the new Maxwell creep model is high when simulating the 3D creep properties of remolded loess. According to Figure 11 and Table 2, the accuracies of the Nishihara and Burgers models when fitting the creep test data corresponding to the last stress loading of remolded loess are obviously lower than those of the new Maxwell model. According to the analysis presented above, the new Maxwell creep model is capable of describing the complete creep behaviour of remolded loess, especially the accelerated creep phase, under a 3D stress state.


## 4.1. Test Conditions and Results

To further explore the ability of the proposed model to simulate the creep behaviour of geotechnical materials under a 3D stress state, the creep results of remolded loess under triaxial conditions were obtained. The loess samples were collected from the L6 loess layer in the new district of Yan’an city; it belongs to the quaternary middle Pleistocene loess, also referred to as theQ2 loess. The undisturbed loess was ground, dried, and sieved through a 2-mm screen to have a moisture content of 10%; it was kept for 24 h without allowing evaporation to allow the water to be uniformly distributed in the loess. By controlling the compactness levels (k) at 0.84, 0.89, 0.95, and 0.99, the amount of loess required for forming remolded loess under different compactness levels was prepared. The prepared loess was filled into a cylindrical mold to prepare remolded loess specimens and then compressed as per set compression rates using a sample-making machine. The prepared remolded loess specimens were cylinders with diameters and heights of 39.8 and 80 mm, respectively.The SLB-1A stress-strain controlled test machine was used for triaxial creep shear measurements, as shown in Figure8. This instrument can be stress-loaded through a pneumatic pump, so the set stress can reach a predetermined value in a short time. Before performing the creep tests of remolded loess, the ultimate deviatoric stress values of remolded loess specimens were determined via triaxial consolidated undrained compression tests at a fixed confining pressure of 100 kPa. Under the same confining pressure, the specimens were first consolidated for 24 h; then, indoor creep tests using remolded loess specimens with different compactness levels were conducted under different axial deviatoric stress levels. The axial deviatoric stress was applied step by step using stress levels of R = 20, 40, 60, 80, 90, 95, and 98% (R denotes the ratio of the deviatoric stress to the ultimate stress), and the axial pressure exerted on the specimens during loading was kept constant with a precision of ±5 kPa for 12 h.Figure 8
Triaxial creep shear test apparatus: (a) instrument details; (b) creep test state.A cylindrical specimen was prepared and wrapped on the outside with a rubber film with a tubular mold. Filter paper and permeable stone were placed at both ends of the sample, and the lower end of the sample rubber band was tied to the base and the upper end to the top. Next, the pressure chamber cover was installed and turned on the water switch to fill the pressure chamber with water. At this point, we opened the drain valve and applied the required confining pressure of 100 kPa through stress loading; this pressure was maintained for 24 h to allow the specimen to consolidate. We observed changes in the pore water pressure during consolidation, and we then closed the drain valve and applied axial deviatoric stresses were applied in a stepwise manner at the foregoing levels and recorded axial deformation data using an external strain gauge during creep.Based on the results of the above creep tests, creep curves at various compactness levels and deviatoric stresses were drawn, as shown in Figure9. It is evident that axial creep deformations under various compactness levels are roughly similar; the primary, secondary, and tertiary creep stages can be observed under partial deviatoric stress. The specimens at each step during application of the deviatoric stress showed the primary creep stage where the strain rate gradually decreased. After the primary creep stage, the specimens immediately entered the longer secondary creep stage; in this stage, the strain rate tended to be zero under low deviatoric stress; however, under high deviatoric stress, the creep strain increased linearly with time when the strain rate was greater than zero and constant. Next, under a particular high deviatoric stress, the tertiary creep stage developed with time following the secondary creep stage and the axial strain and strain rate significantly increased at this stage.Figure 9
Axial strain versus time at different compactness levels: (a)k = 0.84; (b) k = 0.89; (c) k = 0.95; (d) k = 0.99.
(a)
(b)
(c)
(d)


## 4.2. Parameter Determination

Considering the creep test results of remolded loess samples with a compactness of 0.95 as an example, we used the 3D creep equation of the new Maxwell model to fit the experimental data. The parameters in equation (15) can also be determined by the Levenberg–Marquardt method, and the results are presented in Table 2. The efficacy of the proposed creep model in fitting the data of remolded loess is shown in Figure 10. To further analyse the ability of the model to describe the accelerated creep stage of remolded loess under a 3D stress state, we selected the creep data corresponding to the last stress loading and compared it to the Burgers and Nishihara models, as shown in Figure 11; the parameters of the Burgers and Nishihara models are given in Table 2.Table 2
Parameters determined through fitting the analysis based on creep test results of remolded loess.
New Maxwell creep model
σ
1 – σ3 (kPa)
G (kPa)
K (kPa)
ξ′(kPa·hβ)
β
α
σ
s (kPa)
t
F (h)
R
2
122
29856.761
742.362
0.018
0.041
0.974
244
209.192
91.334
0.003
0.097
0.975
366
5832.610
42.649
0.001
0.181
0.954
488
90.682
345.728
0.0008
0.332
0.978
550
2702.403
33.217
0.005
0.456
0.785
490
6.3
0.998
Burgers model
σ
1 – σ3 (kPa)
G
0 (kPa)
K (kPa)
G
k (kPa·h)
η
k (kPa·h)
η (kPa·h)
R
2
550
105.281
48.865
0.235
25516.903
431.671
0.767
Nishihara model
σ
1 – σ3 (kPa)
G
0 (kPa)
K (kPa)
G
k (kPa·h)
η
k (kPa·h)
η
s (kPa·h)
σ
s (kPa)
R
2
550
80.301
71.612
0.155
28047.865
47.019
490
0.797Figure 10
Comparison between theoretical and experimental curves.Figure 11
Comparison between theoretical and experimental curves of (a) new Maxwell model and Burgers model and (b) new Maxwell model and Nishihara model.
(a)
(b)As observed in Figure10 and Table 2, the calculated results using the new Maxwell model are in good agreement with the test results obtained under deviatoric stresses of 122, 244, 366, 488, and 550 kPa. The correlation coefficient squares (R2) of the five curves are 0.974, 0.975, 0.954, 0.978, and 0.996, respectively, indicating that the prediction precision of the new Maxwell creep model is high when simulating the 3D creep properties of remolded loess. According to Figure 11 and Table 2, the accuracies of the Nishihara and Burgers models when fitting the creep test data corresponding to the last stress loading of remolded loess are obviously lower than those of the new Maxwell model. According to the analysis presented above, the new Maxwell creep model is capable of describing the complete creep behaviour of remolded loess, especially the accelerated creep phase, under a 3D stress state.


## 5. Discussion

To describe the time-dependence of the mechanical properties of geotechnical materials in the entire creep process, a new Maxwell creep model based on fractional derivatives and elastic-plastic elements was developed, and creep constitutive equations were established with explanations. Through simulation capability analysis of the model proposed herein, we found that the new model showed good agreement with the experimental creep data of rock salt and remolded loess, especially for the accelerated creep phase. Further, we can observe from Figure5 that the new Maxwell creep model has better accuracy in fitting the experimental data of rock salt under a 1D stress state, especially in the case of the intact creep behaviour, compared with the classical Burgers and Nishihara models. The new Maxwell creep model has only two elements and four parameters, thus making it very simple compared to most other creep models; moreover, by changing parameters β and α, the steady creep and accelerated creep stages with different creep rates can be simulated effectively. Similarly, as shown in Figures 10 and 11, compared with the Burgers and Nishihara models, the new Maxwell creep model also has a better effect in describing the creep characteristics of remolded loess under a 3D stress state. In general, our results demonstrate that the new Maxwell creep model can characterise the creep behaviour of geotechnical materials with better accuracy and fewer parameters than extant models. Therefore, this model can be used extensively in the field of geotechnical engineering, and it can help in developing and clarifying creep constitutive models.


## 6. Conclusions

In this study, we employed fractional and elastic-plastic elements to construct a new Maxwell creep model and derived the corresponding creep functions. The developed model is simple and comprises four parameters and two elements. The fitting results provided by the new Maxwell creep model were also analysed. The results indicate that this model is capable of analysing the creep behaviour of geotechnical materials under different loading conditions. Moreover, it can effectively capture the accelerated creep behaviour of geotechnical materials, which is difficult to observe in most of existing creep models. Additionally, the encouraging results provided by this model can enable other scholars to study the creep models with renewed interest.



---

*Source: 9170706-2020-02-10.xml*