# Investigation Progresses and Applications of Fractional Derivative Model in Geotechnical Engineering

Authors: Jinxing Lai; Sheng Mao; Junling Qiu; Haobo Fan; Qian Zhang; Zhinan Hu; Jianxun Chen
Journal: Mathematical Problems in Engineering
 (2016)

---

## Abstract

Over the past couple of decades, as a new mathematical tool for addressing a number of tough problems, fractional calculus has been gaining a continually increasing interest in diverse scientific fields, including geotechnical engineering due primarily to geotechnical rheology phenomenon. Unlike the classical constitutive models in which simulation analysis gradually fails to meet the reasonable accuracy of requirement, the fractional derivative models have shown the merits of hereditary phenomena with long memory. Additionally, it is traced that the fractional derivative model is one of the most effective and accurate approaches to describe the rheology phenomenon. In relation to this, an overview aimed first at model structure and parameter determination in combination with application cases based on fractional calculus was provided. Furthermore, this review paper shed light on the practical application aspects of deformation analysis of circular tunnel, rheological settlement of subgrade, and relevant loess researches subjected to the achievements acquired in geotechnical engineering. Finally, concluding remarks and important future investigation directions were pointed out.


---

## Body

## 1. Introduction

In practical engineering, massive unstable failure is closely related to the rheological characteristic of rock and soil. Recently, owing to guarantee of the long-term safety and stability of geotechnical engineering, an increasing attention has been paid by relevant scholars to the investigations on the rheology, whose core content is constitutive model [1, 2]. Generally, the classical constitutive model for rock and soil can be classified into three distinct approaches, namely, empirical and semiempirical models, rheology theory based models, and visco-elastic-plastic component models. The empirical and semiempirical models, which are based on experiments, were endowed with a small number of parameters. On the other hand, the rheology theory based models, which focus on the mechanical response, thus reveal the properties of rheology at the microscale. The visco-elastic-plastic component model, as a common approach, which is composed of standard elements such as the Hooke spring and the Newtonian dashpot, has the advantage of flexible description of different rheology deformation with visual concept and clear physical significance [3–6]. Nevertheless, the empirical and semiempirical models are restricted by the duration of experiment, whose results can only describe the stable stage of creep. Meanwhile, the visco-elastic-plastic component models require more elements to establish constitutive equations so as to meet precision requirement [7], causing the problems of many parameters, complicated constitutive relation, and difficulty in determination of parameters, resulting in inconvenience for actual engineering. On the whole, because of its natural defects, the classical models have difficulty in characterizing accurately the whole process curve of creep of rock and soil. Specifically, at the primary stage of creep, its theoretical values deviate much from the experimental results, especially if the accelerated stage appears subsequently [8, 9].As a result, more and more investigators employ fractional calculus to the establishment of constitutive models when obtaining a model with a small number of parameters and accurate description is the target point [10–12]. In fact, during last decades, the field of fractional calculus has attracted increasing interest of researchers in very diverse scientific fields including mathematics, chemistry, engineering technology, anomalous diffusion, and even medical machinery and robot control [13–19]. Indeed, fractional calculus which is defined in the form of calculus with concise equation has been found to be an accurate mathematical tool for solving the difficulty during physical and mechanical modeling. As for the rheology, which is featured by time dependency, the fractional derivative model can be a perfect approach for describing rheology phenomenon, in particular to the properties of hereditary phenomena with long memory [20, 21]. Over the years, after devoting many continuous efforts toward the study on the application of fractional calculus by numerous researchers, the fractional derivative model has gained much deep development. Specific to the parameters, which are the key point in the process of application and to determine the feasibility of model, it can be obtained through inversion based on experiment data, numerical calculation, and software fitting. As a distinct approach, fractional derivative model is applied to describe the nonlinear viscoelastic property, deformation property, shear contraction and dilatation, and damage growth during accelerated stage of creep in geotechnical engineering [9, 22, 23]. For more details on the creep phenomena, fractional derivative model is regarded as an effective approach that can not only well reflect the process of nonlinear gradual change at primary stage but also present accurately at steady stage and accelerated stage under high stress. However, due to its pseudodifferential operator, the storage memory property of fractional derivative model brings certain difficulty to numerical calculation despite the fact that the rheological process curve is beautifully presented [24].In this paper, we systematically investigated current research in the field of fractional derivative model structure and acquisition of parameters where the resulting models had been incorporated into creep analysis. We began with an introduction to relevant concepts and developments of fractional calculus, providing an up-to-date review of key constitutive models which had been widely applied in the field of geotechnical engineering. Following this we comprehensively reviewed the state of the art subjected to practical applications of fractional derivative model, including three major aspects of deformation analysis of circular tunnel, rheological settlement of the subgrade, and relevant loess researches.


## 2. Model Structure Research

### 2.1. Presentation and Development of Fractional Calculus

The fractional calculus has been verified, that is, a branch of mathematics which investigates the property of the differentiation and integration of any order. Its embryonic form can be dated back to September 30, 1695, when Leibniz replied to the famous question: “What does the derivativednf/dxn for f(x) mean if n=1/2?” proposed by L’ Hospital [25]. Due to its complicated calculation and undefined physical significance, the focus at early stage is mainly on research in pure mathematics, and for a long period of time, the research on fractional calculus has attracted little attention from scholars in natural science and engineering technology field. In 1936, Gemant firstly applied fractional derivatives to problems in elasticity [26]. Until 1982, Mandelbrot [27] who was the founder of fractal theory pointed out that there was a self-similarity phenomenon between integer and fraction, and many fractal dimensions existed in nature and many technological sciences. Subsequently, researchers found relevant phenomenon or process of fractal geometry, power law phenomenon, and memory process, which could establish a close relationship with fractional calculus, thus leading to an explosion of research activities on fractional calculus. The contribution of Shestopal and Goss [28] was firstly applied to fractional derivative in Maxwell and Kelvin model for analysis of creep buckling. Bagley and Torvik [29] conducted research on three-dimensional constitutive relation of fractional derivative, limitation of thermodynamics on model parameters, and finite element method. Through utilizing the linear standard system of fractional derivative, Zhang [30] established the constitutive relation of nitrile rubber and proved that it was equivalent with Rabotnov model. Metzler and Nonnenmacher [31] took Zener model based on fractional order to explore the abnormal relaxation of polymeric materials. Yin et al. [32] firstly considered fractional derivative element as soft-matter element while building the constitutive model of soil rheology phenomena, representing the material between an ideal fluid and an ideal solid.Up to now, there are many different definitions of fractional calculus [33–35], among which basic ideology is to look at a fractional derivative as the inverse operation of a fractional integral. It is usually chosen as the Riemann-Liouville form:(1)d-γftdt-γ=Dt0t-γft=∫t0tt-ξγ-1Γγfξdξ,t∈t0,b,where f(t) is integrable function on the interval t0,b with the order γ>0, n-1<γ≤n (n is integer) and Γ(γ) is Euler’s Gamma function, which is defined by(2)Γγ=∫0∞e-ttγ-1dt.Correspondingly, the Riemann-Liouville fractional derivative is given by(3)dγftdtγ=Dt0tγft=dndtnDt0t-n-γft.Laplace transformation formula of fractional differential can be expressed as(4)LD0t-γft,p=p-γf¯p,γ≥0,LD0tγft,p=pγf¯p,γ≤0,where f(t) can be integrated near t=0 with 0≤γ≤1 and f¯(p) is Laplace transformation of f(t).


### 2.2. Research Achievements of Model Structure

Just as fractional calculus is a generalization of the commonly used integer-order differentiation and integration, fractional derivative model continues classical model’s feature of being concise and easy to understand, which not only can make up the shortcomings of classical model but also has better description effects. Its essence is to replace a Newtonian dashpot in classical model with the fractional derivative Abel dashpot and thus obtain fractional derivative model through series-parallel connection with other several basic elements. In geotechnical engineering, the functional equations of fractional model are mainly utilized to perform the process curve of rheology phenomena, with an attempt to adjust deformation rate and control the stress (strain) [36, 37].


#### 2.2.1. Fractional Derivative Abel Dashpot Research

There are two kinds of most basic kernel operators: Abel kernel and Rabotnov kernel. The fractional derivative element established based on Rabotnov kernel has the advantage of describing rheological curve with horizontal asymptote, but its equations are quite complicated, resulting in hardly acquired parameters. For a comparison, using Abel kernel to investigate creep phenomenon, the basic ideology is to research the function forms after derivation of creep compliance for modified Abel dashpot. Its relevant model has fewer parameters and stronger adaptability, which is demonstrated by the increasing number of papers and special issues in journals [38–43].Function equation of Abel kernel [36] is written as follows:(5)ϕt=λt-γΓ1-γ.(
1) Abel Dashpot with Constant Coefficient. The constitutive equation of Abel dashpot with constant coefficient shown in Figure 1 is defined as follows:(6)σt=ηγdγεtdtγ0≤γ≤1,where ηγ is the viscosity coefficient. In a special case of γ=1, the constitutive equation (6) is usually applied to describe the Newton dashpot, being equivalent to an ideal fluid with ηγ=η, while a spring in the case of γ=0 is equivalent to an ideal solid with η0=E. Therefore, Abel dashpot with constant coefficient bridges the gap between an ideal fluid and an ideal solid material response by interpolating between the Newton dashpot and a spring.Figure 1
Abel dashpot with constant coefficient [32].Under the circumstance of constant stress in (6), namely, σ(t)=σ, the element denotes the creep strain characteristic. Conducting the fractional integral calculation of (6) based on the Riemann-Liouville operator, we can gain(7)εt=σηγtγΓ1+γ0≤γ≤1.(
2) Abel Dashpot with Variable Coefficient. During the process of rheology, in particular at the accelerated stage, the cracking and damage growth at the microscale must be considered; hence, the viscosity coefficient of rock or soil is no longer remaining unchanged, which should be described by introducing the damage variable [2, 44, 45]. As a result, the constitutive relation of Abel dashpot with variable coefficient shown in Figure 2 deduces that (8)σt=ηγ1-Dωdγεtdtγ0<γ≤1,where Dω is the damage variable defined by 0≤Dω=1-e-αt<1 and ηγ(1-Dω)=ηγe-αt indicates the variable viscosity coefficient; α is related to properties of materials.Figure 2
Abel dashpot with variable coefficient [2].


#### 2.2.2. Application of Fractional Derivative Model

In terms of the issues in geotechnical engineering, fractional derivative models are usually composed of basic element, constant coefficient Abel dashpot, and variable coefficient Abel dashpot in series and parallel, which now have evolved from simple models such as fractional Maxwell and Kelvin models to the fractional Nishihara model accounting for the efforts researchers devoted. In contrast, the simple models are connected in series or in parallel with the constant coefficient Abel dashpot and spring element, while the three elements’ solid model in series with variable coefficient dashpot leads to the fractional Nishihara model. Additionally, three-dimensional constitutive equation is derived from one-dimensional constitutive equation of fractional derivative models. The commonly used models are shown in Figure3.Figure 3
The fractional derivative element model [36].
(a)
The Maxwell model
(b)
The Kelvin model
(c)
The three elements’ solid model (Kelvin-Voigt)
(d)
The Burgers model
(e)
The Nishihara modelNote that different rheological conditions require different models. In Figure3, (a), (b), (c), and (d) belong to viscoelastic model, while (e) belongs to viscoelastoplastic model. In greater detail, (a), (b), and (c) are all suitable for describing creep that tends to stability finally and (c) shows the highest fitting precision. For comparison, (d) has advantage of describing the creep with decay and flow period; (e) is regarded as special one that can approach the three stages with initial stage, stable stage, and accelerated stage.Taking the relevant investigations on soft soil as examples showing the evolution of model structure, Sun and Zhang [46] adopted the fractional derivative Kelvin model, as shown in Figure 3(b), to simulate the creep laws of two stages of soft soil in Pearl River Delta, in which equation form was simple and unified compared with classical model. Additionally, fewer parameters needed to be adjusted during the process of calculation. In addition, from the results shown in Figure 4, the classical Kelvin model brings a large deviation in comparison with the experimental data near the inflection point of fitted curve, and its offset increases with the time. It can be predicted that accordingly, in a very long time, the classical Kelvin model is seriously out of accordance with the experimental data.Figure 4
Curves, respectively, simulated by fractional Kelvin model and Kelvin model (σ=37.5 kPa) [46].Based on fractional calculus theory, Li and Yue [47] proposed a five-element nonlinear rheological constitutive model, which reserved a Newton dashpot, as shown in Figure 5, for Hong Kong marine clay and Shanghai mucky clay under different stress levels. When σ0<σs, the right side of the viscoplastic dashpot does not participate in the work; thus the model is simplified as fractional derivative three-element solid model. While σ0>σs, the viscoplastic dashpot participates in the work; the model can better characterize the three stages of soft soil under high stress with high fitting accuracy, in which relevant coefficient of the whole curve compared with experimental results reaches 0.9696, shown in Figure 6. In this regard, the model has wide applicability in view simultaneously of a small number of parameters which are easy to determine.Figure 5
Fractional order derivative five-element model [47].Figure 6
Test data and fitting curve of five-element fractional derivative model (σ=243 kPa) [47].To research creep properties of Zhanjiang soft clay, He et al. [48] built the modified Burgers model, and then the expression of three-dimensional equation was derived based on hereditary creep theory. Consider(9)εij=sij2G+σm3kδij+sij2G∫0tκtdt,where sij is the deviator stress tensor and G indicates the shear modulus. σm is the confining pressure and k expresses the bulk modulus of deformation. In a special case of i=j, δij=1, while δij=0 in another case.


## 2.1. Presentation and Development of Fractional Calculus

The fractional calculus has been verified, that is, a branch of mathematics which investigates the property of the differentiation and integration of any order. Its embryonic form can be dated back to September 30, 1695, when Leibniz replied to the famous question: “What does the derivativednf/dxn for f(x) mean if n=1/2?” proposed by L’ Hospital [25]. Due to its complicated calculation and undefined physical significance, the focus at early stage is mainly on research in pure mathematics, and for a long period of time, the research on fractional calculus has attracted little attention from scholars in natural science and engineering technology field. In 1936, Gemant firstly applied fractional derivatives to problems in elasticity [26]. Until 1982, Mandelbrot [27] who was the founder of fractal theory pointed out that there was a self-similarity phenomenon between integer and fraction, and many fractal dimensions existed in nature and many technological sciences. Subsequently, researchers found relevant phenomenon or process of fractal geometry, power law phenomenon, and memory process, which could establish a close relationship with fractional calculus, thus leading to an explosion of research activities on fractional calculus. The contribution of Shestopal and Goss [28] was firstly applied to fractional derivative in Maxwell and Kelvin model for analysis of creep buckling. Bagley and Torvik [29] conducted research on three-dimensional constitutive relation of fractional derivative, limitation of thermodynamics on model parameters, and finite element method. Through utilizing the linear standard system of fractional derivative, Zhang [30] established the constitutive relation of nitrile rubber and proved that it was equivalent with Rabotnov model. Metzler and Nonnenmacher [31] took Zener model based on fractional order to explore the abnormal relaxation of polymeric materials. Yin et al. [32] firstly considered fractional derivative element as soft-matter element while building the constitutive model of soil rheology phenomena, representing the material between an ideal fluid and an ideal solid.Up to now, there are many different definitions of fractional calculus [33–35], among which basic ideology is to look at a fractional derivative as the inverse operation of a fractional integral. It is usually chosen as the Riemann-Liouville form:(1)d-γftdt-γ=Dt0t-γft=∫t0tt-ξγ-1Γγfξdξ,t∈t0,b,where f(t) is integrable function on the interval t0,b with the order γ>0, n-1<γ≤n (n is integer) and Γ(γ) is Euler’s Gamma function, which is defined by(2)Γγ=∫0∞e-ttγ-1dt.Correspondingly, the Riemann-Liouville fractional derivative is given by(3)dγftdtγ=Dt0tγft=dndtnDt0t-n-γft.Laplace transformation formula of fractional differential can be expressed as(4)LD0t-γft,p=p-γf¯p,γ≥0,LD0tγft,p=pγf¯p,γ≤0,where f(t) can be integrated near t=0 with 0≤γ≤1 and f¯(p) is Laplace transformation of f(t).


## 2.2. Research Achievements of Model Structure

Just as fractional calculus is a generalization of the commonly used integer-order differentiation and integration, fractional derivative model continues classical model’s feature of being concise and easy to understand, which not only can make up the shortcomings of classical model but also has better description effects. Its essence is to replace a Newtonian dashpot in classical model with the fractional derivative Abel dashpot and thus obtain fractional derivative model through series-parallel connection with other several basic elements. In geotechnical engineering, the functional equations of fractional model are mainly utilized to perform the process curve of rheology phenomena, with an attempt to adjust deformation rate and control the stress (strain) [36, 37].


### 2.2.1. Fractional Derivative Abel Dashpot Research

There are two kinds of most basic kernel operators: Abel kernel and Rabotnov kernel. The fractional derivative element established based on Rabotnov kernel has the advantage of describing rheological curve with horizontal asymptote, but its equations are quite complicated, resulting in hardly acquired parameters. For a comparison, using Abel kernel to investigate creep phenomenon, the basic ideology is to research the function forms after derivation of creep compliance for modified Abel dashpot. Its relevant model has fewer parameters and stronger adaptability, which is demonstrated by the increasing number of papers and special issues in journals [38–43].Function equation of Abel kernel [36] is written as follows:(5)ϕt=λt-γΓ1-γ.(
1) Abel Dashpot with Constant Coefficient. The constitutive equation of Abel dashpot with constant coefficient shown in Figure 1 is defined as follows:(6)σt=ηγdγεtdtγ0≤γ≤1,where ηγ is the viscosity coefficient. In a special case of γ=1, the constitutive equation (6) is usually applied to describe the Newton dashpot, being equivalent to an ideal fluid with ηγ=η, while a spring in the case of γ=0 is equivalent to an ideal solid with η0=E. Therefore, Abel dashpot with constant coefficient bridges the gap between an ideal fluid and an ideal solid material response by interpolating between the Newton dashpot and a spring.Figure 1
Abel dashpot with constant coefficient [32].Under the circumstance of constant stress in (6), namely, σ(t)=σ, the element denotes the creep strain characteristic. Conducting the fractional integral calculation of (6) based on the Riemann-Liouville operator, we can gain(7)εt=σηγtγΓ1+γ0≤γ≤1.(
2) Abel Dashpot with Variable Coefficient. During the process of rheology, in particular at the accelerated stage, the cracking and damage growth at the microscale must be considered; hence, the viscosity coefficient of rock or soil is no longer remaining unchanged, which should be described by introducing the damage variable [2, 44, 45]. As a result, the constitutive relation of Abel dashpot with variable coefficient shown in Figure 2 deduces that (8)σt=ηγ1-Dωdγεtdtγ0<γ≤1,where Dω is the damage variable defined by 0≤Dω=1-e-αt<1 and ηγ(1-Dω)=ηγe-αt indicates the variable viscosity coefficient; α is related to properties of materials.Figure 2
Abel dashpot with variable coefficient [2].


### 2.2.2. Application of Fractional Derivative Model

In terms of the issues in geotechnical engineering, fractional derivative models are usually composed of basic element, constant coefficient Abel dashpot, and variable coefficient Abel dashpot in series and parallel, which now have evolved from simple models such as fractional Maxwell and Kelvin models to the fractional Nishihara model accounting for the efforts researchers devoted. In contrast, the simple models are connected in series or in parallel with the constant coefficient Abel dashpot and spring element, while the three elements’ solid model in series with variable coefficient dashpot leads to the fractional Nishihara model. Additionally, three-dimensional constitutive equation is derived from one-dimensional constitutive equation of fractional derivative models. The commonly used models are shown in Figure3.Figure 3
The fractional derivative element model [36].
(a)
The Maxwell model
(b)
The Kelvin model
(c)
The three elements’ solid model (Kelvin-Voigt)
(d)
The Burgers model
(e)
The Nishihara modelNote that different rheological conditions require different models. In Figure3, (a), (b), (c), and (d) belong to viscoelastic model, while (e) belongs to viscoelastoplastic model. In greater detail, (a), (b), and (c) are all suitable for describing creep that tends to stability finally and (c) shows the highest fitting precision. For comparison, (d) has advantage of describing the creep with decay and flow period; (e) is regarded as special one that can approach the three stages with initial stage, stable stage, and accelerated stage.Taking the relevant investigations on soft soil as examples showing the evolution of model structure, Sun and Zhang [46] adopted the fractional derivative Kelvin model, as shown in Figure 3(b), to simulate the creep laws of two stages of soft soil in Pearl River Delta, in which equation form was simple and unified compared with classical model. Additionally, fewer parameters needed to be adjusted during the process of calculation. In addition, from the results shown in Figure 4, the classical Kelvin model brings a large deviation in comparison with the experimental data near the inflection point of fitted curve, and its offset increases with the time. It can be predicted that accordingly, in a very long time, the classical Kelvin model is seriously out of accordance with the experimental data.Figure 4
Curves, respectively, simulated by fractional Kelvin model and Kelvin model (σ=37.5 kPa) [46].Based on fractional calculus theory, Li and Yue [47] proposed a five-element nonlinear rheological constitutive model, which reserved a Newton dashpot, as shown in Figure 5, for Hong Kong marine clay and Shanghai mucky clay under different stress levels. When σ0<σs, the right side of the viscoplastic dashpot does not participate in the work; thus the model is simplified as fractional derivative three-element solid model. While σ0>σs, the viscoplastic dashpot participates in the work; the model can better characterize the three stages of soft soil under high stress with high fitting accuracy, in which relevant coefficient of the whole curve compared with experimental results reaches 0.9696, shown in Figure 6. In this regard, the model has wide applicability in view simultaneously of a small number of parameters which are easy to determine.Figure 5
Fractional order derivative five-element model [47].Figure 6
Test data and fitting curve of five-element fractional derivative model (σ=243 kPa) [47].To research creep properties of Zhanjiang soft clay, He et al. [48] built the modified Burgers model, and then the expression of three-dimensional equation was derived based on hereditary creep theory. Consider(9)εij=sij2G+σm3kδij+sij2G∫0tκtdt,where sij is the deviator stress tensor and G indicates the shear modulus. σm is the confining pressure and k expresses the bulk modulus of deformation. In a special case of i=j, δij=1, while δij=0 in another case.


## 2.2.1. Fractional Derivative Abel Dashpot Research

There are two kinds of most basic kernel operators: Abel kernel and Rabotnov kernel. The fractional derivative element established based on Rabotnov kernel has the advantage of describing rheological curve with horizontal asymptote, but its equations are quite complicated, resulting in hardly acquired parameters. For a comparison, using Abel kernel to investigate creep phenomenon, the basic ideology is to research the function forms after derivation of creep compliance for modified Abel dashpot. Its relevant model has fewer parameters and stronger adaptability, which is demonstrated by the increasing number of papers and special issues in journals [38–43].Function equation of Abel kernel [36] is written as follows:(5)ϕt=λt-γΓ1-γ.(
1) Abel Dashpot with Constant Coefficient. The constitutive equation of Abel dashpot with constant coefficient shown in Figure 1 is defined as follows:(6)σt=ηγdγεtdtγ0≤γ≤1,where ηγ is the viscosity coefficient. In a special case of γ=1, the constitutive equation (6) is usually applied to describe the Newton dashpot, being equivalent to an ideal fluid with ηγ=η, while a spring in the case of γ=0 is equivalent to an ideal solid with η0=E. Therefore, Abel dashpot with constant coefficient bridges the gap between an ideal fluid and an ideal solid material response by interpolating between the Newton dashpot and a spring.Figure 1
Abel dashpot with constant coefficient [32].Under the circumstance of constant stress in (6), namely, σ(t)=σ, the element denotes the creep strain characteristic. Conducting the fractional integral calculation of (6) based on the Riemann-Liouville operator, we can gain(7)εt=σηγtγΓ1+γ0≤γ≤1.(
2) Abel Dashpot with Variable Coefficient. During the process of rheology, in particular at the accelerated stage, the cracking and damage growth at the microscale must be considered; hence, the viscosity coefficient of rock or soil is no longer remaining unchanged, which should be described by introducing the damage variable [2, 44, 45]. As a result, the constitutive relation of Abel dashpot with variable coefficient shown in Figure 2 deduces that (8)σt=ηγ1-Dωdγεtdtγ0<γ≤1,where Dω is the damage variable defined by 0≤Dω=1-e-αt<1 and ηγ(1-Dω)=ηγe-αt indicates the variable viscosity coefficient; α is related to properties of materials.Figure 2
Abel dashpot with variable coefficient [2].


## 2.2.2. Application of Fractional Derivative Model

In terms of the issues in geotechnical engineering, fractional derivative models are usually composed of basic element, constant coefficient Abel dashpot, and variable coefficient Abel dashpot in series and parallel, which now have evolved from simple models such as fractional Maxwell and Kelvin models to the fractional Nishihara model accounting for the efforts researchers devoted. In contrast, the simple models are connected in series or in parallel with the constant coefficient Abel dashpot and spring element, while the three elements’ solid model in series with variable coefficient dashpot leads to the fractional Nishihara model. Additionally, three-dimensional constitutive equation is derived from one-dimensional constitutive equation of fractional derivative models. The commonly used models are shown in Figure3.Figure 3
The fractional derivative element model [36].
(a)
The Maxwell model
(b)
The Kelvin model
(c)
The three elements’ solid model (Kelvin-Voigt)
(d)
The Burgers model
(e)
The Nishihara modelNote that different rheological conditions require different models. In Figure3, (a), (b), (c), and (d) belong to viscoelastic model, while (e) belongs to viscoelastoplastic model. In greater detail, (a), (b), and (c) are all suitable for describing creep that tends to stability finally and (c) shows the highest fitting precision. For comparison, (d) has advantage of describing the creep with decay and flow period; (e) is regarded as special one that can approach the three stages with initial stage, stable stage, and accelerated stage.Taking the relevant investigations on soft soil as examples showing the evolution of model structure, Sun and Zhang [46] adopted the fractional derivative Kelvin model, as shown in Figure 3(b), to simulate the creep laws of two stages of soft soil in Pearl River Delta, in which equation form was simple and unified compared with classical model. Additionally, fewer parameters needed to be adjusted during the process of calculation. In addition, from the results shown in Figure 4, the classical Kelvin model brings a large deviation in comparison with the experimental data near the inflection point of fitted curve, and its offset increases with the time. It can be predicted that accordingly, in a very long time, the classical Kelvin model is seriously out of accordance with the experimental data.Figure 4
Curves, respectively, simulated by fractional Kelvin model and Kelvin model (σ=37.5 kPa) [46].Based on fractional calculus theory, Li and Yue [47] proposed a five-element nonlinear rheological constitutive model, which reserved a Newton dashpot, as shown in Figure 5, for Hong Kong marine clay and Shanghai mucky clay under different stress levels. When σ0<σs, the right side of the viscoplastic dashpot does not participate in the work; thus the model is simplified as fractional derivative three-element solid model. While σ0>σs, the viscoplastic dashpot participates in the work; the model can better characterize the three stages of soft soil under high stress with high fitting accuracy, in which relevant coefficient of the whole curve compared with experimental results reaches 0.9696, shown in Figure 6. In this regard, the model has wide applicability in view simultaneously of a small number of parameters which are easy to determine.Figure 5
Fractional order derivative five-element model [47].Figure 6
Test data and fitting curve of five-element fractional derivative model (σ=243 kPa) [47].To research creep properties of Zhanjiang soft clay, He et al. [48] built the modified Burgers model, and then the expression of three-dimensional equation was derived based on hereditary creep theory. Consider(9)εij=sij2G+σm3kδij+sij2G∫0tκtdt,where sij is the deviator stress tensor and G indicates the shear modulus. σm is the confining pressure and k expresses the bulk modulus of deformation. In a special case of i=j, δij=1, while δij=0 in another case.


## 3. Determination of Model Parameters

Simplicity and accuracy are two significant benchmarks to measure the applicability of constitutive model [39–41]. Lewandowski and Chorążyczewski [42] pointed out that an important problem, connected with the fractional rheological models, is an estimation of the model parameters from experimental data. To some extent, the key to accurately apply the fractional derivative model to practical engineering lies in the precision of fitting parameters [42, 43].Literature [25] investigated that the creep equations of different models have unified form given by(10)εt=Jtσ,where J(t) refers to creep compliance, which have different forms. Their forms [49–52] and the number of parameters are summarized in Table 1.Table 1
Characteristics of the commonly used rheological model of fractional derivative.
Fractional derivative model
Creep compliance (J)
Number of parameters
Kelvin
1
q
0
1
-
E
γ
t, t≥0 Eγt=1+∑n=1∞(-1)n(t/η)γnΓ(1+γn)
3
Three elements’ solid (Kelvin-Voigt)
1
E
0
+
1
η
γ
∑
n
=
0
∞
-
E
1
/
η
γ
n
t
γ
1
+
n
γ
1
+
n
Γ
1
+
n
γ
4
Burgers
1
E
1
1
-
1
-
E
1
+
E
2
E
2
∑
n
=
0
∞
-
1
n
E
2
η
2
n
+
1
t
γ
1
+
n
Γ
γ
1
+
n
+
1
+
1
η
1
t
γ
Γ
r
+
1
5
Nishihara
1
E
0
+
1
E
1
1
-
∑
n
=
0
∞
-
1
n
t
E
1
/
η
1
γ
n
Γ
1
+
γ
n
,
σ
<
σ
s 1E0+1E11-∑n=0∞-1ntE1/η1γnΓ1+γn+1η2γ2etγ2,σ>σs
6


### 3.1. Parameters Identification with Method of Numerical Analysis

Generally, the numerical methods for determining the parameters of rheological model mainly include regression inversion method, least square method, and decomposition method based on rheological curve. Nevertheless, all these optimized algorithms have some restricted conditions when dealing with relevant issues; that is, the optimized function of constitutive relation should be differentiable. Along with the rapid development of numerical calculation and analysis, the estimation of the fractional derivative model has been brought to a whole new level.


#### 3.1.1. Ant Colony Algorithm

Specific to the interval variables of parameters in viscoelastic Voigt constitutive model that was the fractional derivative three elements’ solid model, Zhao [53] firstly established the solving models for parameters, including the semianalytical finite element model based on time domain homogenous field and the difference finite element model on the basis of regional heterogeneous field. And then using optimized algorithm depended on ant colony theory, as shown in Figure 7; the single and combined identification of constitutive parameters all was realized.Figure 7
Ant colony algorithm [53].


#### 3.1.2. Genetic Algorithm

In order to address the uniaxial creep of artificially frozen soil in the area of Huainan-Huaibei, Wang [36] adopted fractional derivative Kelvin model to simulate the creep characteristics of the first and second stages and then obtained the model parameters through global optimization based on genetic algorithm. As shown in Figure 8, fractional derivative Kelvin model has a higher accuracy of simulation than the classical Kelvin model. In contrast to test data, the former has the relevant coefficient of more than 0.99.Figure 8
Calculation value of frozen soil and the uniaxial creep test value contrast [36].
(a)
Fractional derivative Kelvin model
(b)
The Kelvin model


#### 3.1.3. Particle Swarm Optimization Algorithm

Chen et al. [54] applied particle swarm optimization for the purpose of identifying the creep parameters of generalized Kelvin model, which was built with the experimental data on steady stage of artificially frozen soil under different loading stress; thus material parameters of generalized Kelvin model were obtained.Similarly, Niu [55] built the fractional derivative Nishihara model intended to describe the deformation of frozen soil at accelerating rheology stage; furthermore, particle swarm optimization was applied to optimize the parameters of fractional derivative creep model. The parameters using particle swarm optimization algorithm were shown in Table 2.Table 2
Particle swarm parameters [55].
Populations
Close to populations
C
1
C
2
Maximum acceleration
80
15
3.0
3.0
15Special to the experiment conditions, sandy clay had the loading values of0.7σs, respectively, at the temperature of −5°, −10°, and −15°. The obtained parameters of the fractional derivative Nishihara model were shown in Table 3. As a result, the relevant coefficient for fitting model parameters through particle swarm optimization algorithm was more than 0.98, which highlighted the unique advantage of particle swarm optimization algorithm in application of optimizing the fractional derivative Nishihara model for describing accelerating rheology of frozen soil.Table 3
Parameters of particle swarm fractional order derivative Nishihara model [55].
Soil
Temperature
Loading coefficient
σ
E
1
η
2
γ
σ
s
Correlation coefficient
Clay
−5
0.3
1.05
0.06
—
0.1116
—
0.9869
0.5
1.75
0.28
—
0.72
—
0.9990
−10
0.3
1.5
0.24
—
0.0068
0.9943
0.5
2.5
1.16
0.0198
—
0.9947
−15
0.3
2.1
0.65
—
0.3119
—
0.9966
0.5
3.5
0.18
—
0.3901
—
0.9992
−5
0.7
2.45
0.000646
0.000631
1.0494
2.37
0.9919
−10
0.7
3.5
0.11854
1.0638
0.00811
3.25
0.9873
−15
0.7
4.9
0.01117
0.00548
1.139
4.74
0.9927


#### 3.1.4. Intelligence Algorithm

Simulated annealing (SA) algorithm [56, 57] was incorporated into standard particle swarm optimization (PSO) algorithm by Xu et al. [58], which could minimize the error caused by the limitations of PSO algorithm. The new method obtained was called SAPSO algorithm, as shown in Figure 9. Since PSO algorithm only accepts the optimal solution, point A in Figure 9 will fall into the local optimal solution on the left during the calculating, resulting in the global optimal solution being usually neglected. Subjected to the SA algorithm, point A can be only accepted with a certain probability. Therefore, point B as the global optimal solution has opportunity to be obtained by getting out of the local optimal solution.Figure 9
Sketches of SAPSO algorithm [58].After optimization, rheological curve decomposition-particle swarm-simulated annealing algorithm is established to inverse the multiparameter constitutive model of rock rheology, leading to effectively solving the initial value dependence of nonlinear problem and gaining global optimal solution simultaneously. Additionally, we can obtain from Figure10 that SAPSO algorithm shows rapid changes with fitness value in the first 16 iterative steps and then reaches a high convergence precision.Figure 10
Relationship between fitness and iteration times [58].


### 3.2. Parameter Solving Based on Software

Based on Monte Carlo, Gen et al. [59] took the minimum error of experimental data as the criterion and then made the experimental curves of nine samples fitted with standard fractional derivative linear model, thus determining the relaxation time and coefficient of viscosity of Malan loess. Meanwhile, Sun and Zhang [46] and Li and Yue [47] by means of MATLAB performed effective identification for parameters aspects of fractional derivative Kelvin model under the condition of low stress of soft soil and fractional derivative five-element model under the condition of high stress, which was featured by simplicity and wide applicability.He [60] performed parameter fitting for fractional derivative Maxwell model through 1stOpt mathematical software shown in Figure 11. As a result, the changes law of rheological model parameters for Q2 loess in cases of different water content was gained, as shown in Table 4.Table 4
Parameters fitting results of uniaxial rheological model ofQ2 loess in Yan’an [60].
ω (%)
E
0 (105 kPa)
η (105 kPa⋅min)
γ
5
3.029
57.093
σ
-
0.5813
0.0109
σ
-
0.3369
7
1.361
11.068
σ
-
0.3529
0.0337
σ
-
0.1172
10
1.109
16.234
σ
-
0.3711
0.0496
σ
-
0.0964
12
0.787
3.1993
σ
-
0.1921
0.0265
σ
-
0.2053
15
0.577
3.166
σ
-
0.3786
0.0122
σ
-
0.3588
18
0.419
4.9773
σ
-
0.3832
0.026
σ
-
0.1839
20
0.311
2.5315
σ
-
0.2093
0.0424
σ
-
0.1495
22
0.229
1.2174
σ
-
0.1585
0.0238
σ
-
0.2741Figure 11
Interface of 1stOpt [60].On the basis of parameter identification, hence, constitutive equation of steady creep model for Yan’anQ2 loess under the condition of uniaxial stress was obtained. Consider(11)εt=Jsetσ=1E0+1ηtγ∫0∞e-ttγdtσ,σ<σs,where σs is determined by uniaxial creep experiment.


### 3.3. Parameter Analysis

#### 3.3.1. Parameter Sensitivity Analysis

In order to get an adequate understanding of the effects of parametersγ and α on the time-dependent strain and stress, Zhou et al. [61] applied programmable rheometer to carry out uniaxial rheological experiments for salt rock sample of well Wangchu-1 of Jianghan oilfield. Subsequently, based on the results of experiment, fitting analysis of parameters γ and α on rheological strain was conducted and its effects were discussed simultaneously. Figure 12(a) shows that rheological value and rheological speed greatly depend on fractional derivative order, and, the higher the fractional derivative order, the higher the creep strain, which reflects the flexibility of fractional derivate model in description of rheology phenomena. Moreover, in the special case of α=0, rheological procedures of salt rock never appear to be obvious accelerated stage as shown in Figure 12(b), namely, the third stage, which only appears in the case of α>0. What is more, the accelerated stage is more likely to occur as the values of α increase. The results obtained mutually verify the conclusions of Xu et al. [58], Liu and Zhang [62], Yang and Fang [63], and He et al. [64].Figure 12
Sensitivity of the creep strain to the fractional derivative order in the case [61].
(a)
Parameterγ
(b)
Parameterα


#### 3.3.2. Special Case

Introducing the Mittag-Leffler function, the equations based on the fractional derivative Nishihara model deduce that(12)
ε
t
=
σ
E
0
+
σ
E
1
-
σ
E
1
E
γ
,
1
-
E
1
η
1
γ
t
γ
,
σ
≤
σ
s
,
σ
E
0
+
σ
E
1
-
σ
E
1
E
γ
,
1
-
E
1
η
1
γ
t
γ
+
σ
-
σ
s
η
2
γ
t
γ
E
1,1
+
γ
α
t
,
σ
≥
σ
s
,where Ex,y(z) is Mittag-Leffler function with the case of E1,1(z)=ez and E1,2=ez-1/z.In a special case ofγ=1 and α=0, (12) is simplified as(13)
ε
t
=
σ
E
0
+
σ
E
1
1
-
e
-
E
1
/
η
1
t
,
σ
≤
σ
s
,
σ
E
0
+
σ
E
1
1
-
e
-
E
1
/
η
1
t
+
σ
-
σ
s
η
2
t
,
σ
≥
σ
s
.Note that (13) refers to the constitutive relation of the classical Nishihara model. The classical Nishihara model is a special case of the fractional derivative Nishihara model accordingly [60]. The conclusion is mutually verified with the results of Yin et al. [32], He et al. [48], and Xu et al. [58].


## 3.1. Parameters Identification with Method of Numerical Analysis

Generally, the numerical methods for determining the parameters of rheological model mainly include regression inversion method, least square method, and decomposition method based on rheological curve. Nevertheless, all these optimized algorithms have some restricted conditions when dealing with relevant issues; that is, the optimized function of constitutive relation should be differentiable. Along with the rapid development of numerical calculation and analysis, the estimation of the fractional derivative model has been brought to a whole new level.


### 3.1.1. Ant Colony Algorithm

Specific to the interval variables of parameters in viscoelastic Voigt constitutive model that was the fractional derivative three elements’ solid model, Zhao [53] firstly established the solving models for parameters, including the semianalytical finite element model based on time domain homogenous field and the difference finite element model on the basis of regional heterogeneous field. And then using optimized algorithm depended on ant colony theory, as shown in Figure 7; the single and combined identification of constitutive parameters all was realized.Figure 7
Ant colony algorithm [53].


### 3.1.2. Genetic Algorithm

In order to address the uniaxial creep of artificially frozen soil in the area of Huainan-Huaibei, Wang [36] adopted fractional derivative Kelvin model to simulate the creep characteristics of the first and second stages and then obtained the model parameters through global optimization based on genetic algorithm. As shown in Figure 8, fractional derivative Kelvin model has a higher accuracy of simulation than the classical Kelvin model. In contrast to test data, the former has the relevant coefficient of more than 0.99.Figure 8
Calculation value of frozen soil and the uniaxial creep test value contrast [36].
(a)
Fractional derivative Kelvin model
(b)
The Kelvin model


### 3.1.3. Particle Swarm Optimization Algorithm

Chen et al. [54] applied particle swarm optimization for the purpose of identifying the creep parameters of generalized Kelvin model, which was built with the experimental data on steady stage of artificially frozen soil under different loading stress; thus material parameters of generalized Kelvin model were obtained.Similarly, Niu [55] built the fractional derivative Nishihara model intended to describe the deformation of frozen soil at accelerating rheology stage; furthermore, particle swarm optimization was applied to optimize the parameters of fractional derivative creep model. The parameters using particle swarm optimization algorithm were shown in Table 2.Table 2
Particle swarm parameters [55].
Populations
Close to populations
C
1
C
2
Maximum acceleration
80
15
3.0
3.0
15Special to the experiment conditions, sandy clay had the loading values of0.7σs, respectively, at the temperature of −5°, −10°, and −15°. The obtained parameters of the fractional derivative Nishihara model were shown in Table 3. As a result, the relevant coefficient for fitting model parameters through particle swarm optimization algorithm was more than 0.98, which highlighted the unique advantage of particle swarm optimization algorithm in application of optimizing the fractional derivative Nishihara model for describing accelerating rheology of frozen soil.Table 3
Parameters of particle swarm fractional order derivative Nishihara model [55].
Soil
Temperature
Loading coefficient
σ
E
1
η
2
γ
σ
s
Correlation coefficient
Clay
−5
0.3
1.05
0.06
—
0.1116
—
0.9869
0.5
1.75
0.28
—
0.72
—
0.9990
−10
0.3
1.5
0.24
—
0.0068
0.9943
0.5
2.5
1.16
0.0198
—
0.9947
−15
0.3
2.1
0.65
—
0.3119
—
0.9966
0.5
3.5
0.18
—
0.3901
—
0.9992
−5
0.7
2.45
0.000646
0.000631
1.0494
2.37
0.9919
−10
0.7
3.5
0.11854
1.0638
0.00811
3.25
0.9873
−15
0.7
4.9
0.01117
0.00548
1.139
4.74
0.9927


### 3.1.4. Intelligence Algorithm

Simulated annealing (SA) algorithm [56, 57] was incorporated into standard particle swarm optimization (PSO) algorithm by Xu et al. [58], which could minimize the error caused by the limitations of PSO algorithm. The new method obtained was called SAPSO algorithm, as shown in Figure 9. Since PSO algorithm only accepts the optimal solution, point A in Figure 9 will fall into the local optimal solution on the left during the calculating, resulting in the global optimal solution being usually neglected. Subjected to the SA algorithm, point A can be only accepted with a certain probability. Therefore, point B as the global optimal solution has opportunity to be obtained by getting out of the local optimal solution.Figure 9
Sketches of SAPSO algorithm [58].After optimization, rheological curve decomposition-particle swarm-simulated annealing algorithm is established to inverse the multiparameter constitutive model of rock rheology, leading to effectively solving the initial value dependence of nonlinear problem and gaining global optimal solution simultaneously. Additionally, we can obtain from Figure10 that SAPSO algorithm shows rapid changes with fitness value in the first 16 iterative steps and then reaches a high convergence precision.Figure 10
Relationship between fitness and iteration times [58].


## 3.1.1. Ant Colony Algorithm

Specific to the interval variables of parameters in viscoelastic Voigt constitutive model that was the fractional derivative three elements’ solid model, Zhao [53] firstly established the solving models for parameters, including the semianalytical finite element model based on time domain homogenous field and the difference finite element model on the basis of regional heterogeneous field. And then using optimized algorithm depended on ant colony theory, as shown in Figure 7; the single and combined identification of constitutive parameters all was realized.Figure 7
Ant colony algorithm [53].


## 3.1.2. Genetic Algorithm

In order to address the uniaxial creep of artificially frozen soil in the area of Huainan-Huaibei, Wang [36] adopted fractional derivative Kelvin model to simulate the creep characteristics of the first and second stages and then obtained the model parameters through global optimization based on genetic algorithm. As shown in Figure 8, fractional derivative Kelvin model has a higher accuracy of simulation than the classical Kelvin model. In contrast to test data, the former has the relevant coefficient of more than 0.99.Figure 8
Calculation value of frozen soil and the uniaxial creep test value contrast [36].
(a)
Fractional derivative Kelvin model
(b)
The Kelvin model


## 3.1.3. Particle Swarm Optimization Algorithm

Chen et al. [54] applied particle swarm optimization for the purpose of identifying the creep parameters of generalized Kelvin model, which was built with the experimental data on steady stage of artificially frozen soil under different loading stress; thus material parameters of generalized Kelvin model were obtained.Similarly, Niu [55] built the fractional derivative Nishihara model intended to describe the deformation of frozen soil at accelerating rheology stage; furthermore, particle swarm optimization was applied to optimize the parameters of fractional derivative creep model. The parameters using particle swarm optimization algorithm were shown in Table 2.Table 2
Particle swarm parameters [55].
Populations
Close to populations
C
1
C
2
Maximum acceleration
80
15
3.0
3.0
15Special to the experiment conditions, sandy clay had the loading values of0.7σs, respectively, at the temperature of −5°, −10°, and −15°. The obtained parameters of the fractional derivative Nishihara model were shown in Table 3. As a result, the relevant coefficient for fitting model parameters through particle swarm optimization algorithm was more than 0.98, which highlighted the unique advantage of particle swarm optimization algorithm in application of optimizing the fractional derivative Nishihara model for describing accelerating rheology of frozen soil.Table 3
Parameters of particle swarm fractional order derivative Nishihara model [55].
Soil
Temperature
Loading coefficient
σ
E
1
η
2
γ
σ
s
Correlation coefficient
Clay
−5
0.3
1.05
0.06
—
0.1116
—
0.9869
0.5
1.75
0.28
—
0.72
—
0.9990
−10
0.3
1.5
0.24
—
0.0068
0.9943
0.5
2.5
1.16
0.0198
—
0.9947
−15
0.3
2.1
0.65
—
0.3119
—
0.9966
0.5
3.5
0.18
—
0.3901
—
0.9992
−5
0.7
2.45
0.000646
0.000631
1.0494
2.37
0.9919
−10
0.7
3.5
0.11854
1.0638
0.00811
3.25
0.9873
−15
0.7
4.9
0.01117
0.00548
1.139
4.74
0.9927


## 3.1.4. Intelligence Algorithm

Simulated annealing (SA) algorithm [56, 57] was incorporated into standard particle swarm optimization (PSO) algorithm by Xu et al. [58], which could minimize the error caused by the limitations of PSO algorithm. The new method obtained was called SAPSO algorithm, as shown in Figure 9. Since PSO algorithm only accepts the optimal solution, point A in Figure 9 will fall into the local optimal solution on the left during the calculating, resulting in the global optimal solution being usually neglected. Subjected to the SA algorithm, point A can be only accepted with a certain probability. Therefore, point B as the global optimal solution has opportunity to be obtained by getting out of the local optimal solution.Figure 9
Sketches of SAPSO algorithm [58].After optimization, rheological curve decomposition-particle swarm-simulated annealing algorithm is established to inverse the multiparameter constitutive model of rock rheology, leading to effectively solving the initial value dependence of nonlinear problem and gaining global optimal solution simultaneously. Additionally, we can obtain from Figure10 that SAPSO algorithm shows rapid changes with fitness value in the first 16 iterative steps and then reaches a high convergence precision.Figure 10
Relationship between fitness and iteration times [58].


## 3.2. Parameter Solving Based on Software

Based on Monte Carlo, Gen et al. [59] took the minimum error of experimental data as the criterion and then made the experimental curves of nine samples fitted with standard fractional derivative linear model, thus determining the relaxation time and coefficient of viscosity of Malan loess. Meanwhile, Sun and Zhang [46] and Li and Yue [47] by means of MATLAB performed effective identification for parameters aspects of fractional derivative Kelvin model under the condition of low stress of soft soil and fractional derivative five-element model under the condition of high stress, which was featured by simplicity and wide applicability.He [60] performed parameter fitting for fractional derivative Maxwell model through 1stOpt mathematical software shown in Figure 11. As a result, the changes law of rheological model parameters for Q2 loess in cases of different water content was gained, as shown in Table 4.Table 4
Parameters fitting results of uniaxial rheological model ofQ2 loess in Yan’an [60].
ω (%)
E
0 (105 kPa)
η (105 kPa⋅min)
γ
5
3.029
57.093
σ
-
0.5813
0.0109
σ
-
0.3369
7
1.361
11.068
σ
-
0.3529
0.0337
σ
-
0.1172
10
1.109
16.234
σ
-
0.3711
0.0496
σ
-
0.0964
12
0.787
3.1993
σ
-
0.1921
0.0265
σ
-
0.2053
15
0.577
3.166
σ
-
0.3786
0.0122
σ
-
0.3588
18
0.419
4.9773
σ
-
0.3832
0.026
σ
-
0.1839
20
0.311
2.5315
σ
-
0.2093
0.0424
σ
-
0.1495
22
0.229
1.2174
σ
-
0.1585
0.0238
σ
-
0.2741Figure 11
Interface of 1stOpt [60].On the basis of parameter identification, hence, constitutive equation of steady creep model for Yan’anQ2 loess under the condition of uniaxial stress was obtained. Consider(11)εt=Jsetσ=1E0+1ηtγ∫0∞e-ttγdtσ,σ<σs,where σs is determined by uniaxial creep experiment.


## 3.3. Parameter Analysis

### 3.3.1. Parameter Sensitivity Analysis

In order to get an adequate understanding of the effects of parametersγ and α on the time-dependent strain and stress, Zhou et al. [61] applied programmable rheometer to carry out uniaxial rheological experiments for salt rock sample of well Wangchu-1 of Jianghan oilfield. Subsequently, based on the results of experiment, fitting analysis of parameters γ and α on rheological strain was conducted and its effects were discussed simultaneously. Figure 12(a) shows that rheological value and rheological speed greatly depend on fractional derivative order, and, the higher the fractional derivative order, the higher the creep strain, which reflects the flexibility of fractional derivate model in description of rheology phenomena. Moreover, in the special case of α=0, rheological procedures of salt rock never appear to be obvious accelerated stage as shown in Figure 12(b), namely, the third stage, which only appears in the case of α>0. What is more, the accelerated stage is more likely to occur as the values of α increase. The results obtained mutually verify the conclusions of Xu et al. [58], Liu and Zhang [62], Yang and Fang [63], and He et al. [64].Figure 12
Sensitivity of the creep strain to the fractional derivative order in the case [61].
(a)
Parameterγ
(b)
Parameterα


### 3.3.2. Special Case

Introducing the Mittag-Leffler function, the equations based on the fractional derivative Nishihara model deduce that(12)
ε
t
=
σ
E
0
+
σ
E
1
-
σ
E
1
E
γ
,
1
-
E
1
η
1
γ
t
γ
,
σ
≤
σ
s
,
σ
E
0
+
σ
E
1
-
σ
E
1
E
γ
,
1
-
E
1
η
1
γ
t
γ
+
σ
-
σ
s
η
2
γ
t
γ
E
1,1
+
γ
α
t
,
σ
≥
σ
s
,where Ex,y(z) is Mittag-Leffler function with the case of E1,1(z)=ez and E1,2=ez-1/z.In a special case ofγ=1 and α=0, (12) is simplified as(13)
ε
t
=
σ
E
0
+
σ
E
1
1
-
e
-
E
1
/
η
1
t
,
σ
≤
σ
s
,
σ
E
0
+
σ
E
1
1
-
e
-
E
1
/
η
1
t
+
σ
-
σ
s
η
2
t
,
σ
≥
σ
s
.Note that (13) refers to the constitutive relation of the classical Nishihara model. The classical Nishihara model is a special case of the fractional derivative Nishihara model accordingly [60]. The conclusion is mutually verified with the results of Yin et al. [32], He et al. [48], and Xu et al. [58].


## 3.3.1. Parameter Sensitivity Analysis

In order to get an adequate understanding of the effects of parametersγ and α on the time-dependent strain and stress, Zhou et al. [61] applied programmable rheometer to carry out uniaxial rheological experiments for salt rock sample of well Wangchu-1 of Jianghan oilfield. Subsequently, based on the results of experiment, fitting analysis of parameters γ and α on rheological strain was conducted and its effects were discussed simultaneously. Figure 12(a) shows that rheological value and rheological speed greatly depend on fractional derivative order, and, the higher the fractional derivative order, the higher the creep strain, which reflects the flexibility of fractional derivate model in description of rheology phenomena. Moreover, in the special case of α=0, rheological procedures of salt rock never appear to be obvious accelerated stage as shown in Figure 12(b), namely, the third stage, which only appears in the case of α>0. What is more, the accelerated stage is more likely to occur as the values of α increase. The results obtained mutually verify the conclusions of Xu et al. [58], Liu and Zhang [62], Yang and Fang [63], and He et al. [64].Figure 12
Sensitivity of the creep strain to the fractional derivative order in the case [61].
(a)
Parameterγ
(b)
Parameterα


## 3.3.2. Special Case

Introducing the Mittag-Leffler function, the equations based on the fractional derivative Nishihara model deduce that(12)
ε
t
=
σ
E
0
+
σ
E
1
-
σ
E
1
E
γ
,
1
-
E
1
η
1
γ
t
γ
,
σ
≤
σ
s
,
σ
E
0
+
σ
E
1
-
σ
E
1
E
γ
,
1
-
E
1
η
1
γ
t
γ
+
σ
-
σ
s
η
2
γ
t
γ
E
1,1
+
γ
α
t
,
σ
≥
σ
s
,where Ex,y(z) is Mittag-Leffler function with the case of E1,1(z)=ez and E1,2=ez-1/z.In a special case ofγ=1 and α=0, (12) is simplified as(13)
ε
t
=
σ
E
0
+
σ
E
1
1
-
e
-
E
1
/
η
1
t
,
σ
≤
σ
s
,
σ
E
0
+
σ
E
1
1
-
e
-
E
1
/
η
1
t
+
σ
-
σ
s
η
2
t
,
σ
≥
σ
s
.Note that (13) refers to the constitutive relation of the classical Nishihara model. The classical Nishihara model is a special case of the fractional derivative Nishihara model accordingly [60]. The conclusion is mutually verified with the results of Yin et al. [32], He et al. [48], and Xu et al. [58].


## 4. Engineering Application

### 4.1. Deformation Analysis of Circular Tunnel

As the traditional issue of underground engineering, the stress of surrounding rock mass is redistributed due to tunnel excavation, resulting in occasionally a large delayed deformation that might lead to a delayed failure of structure and rock mass collapse, especially for the tunnel in soft rock [65, 66]. Therefore, the deformation mechanism of surrounding rock should be investigated in the perspective of rheology behavior, since the property of associated time effects.Liu and Zhang [62] proposed to take advantage of fractional derivative Kelvin model to investigate the deformation properties of horizontal circular tunnel [67, 68]. With the analysis of the result, it was shown that fractional derivative model could perform the relaxation properties of viscoelastic rock mass. In addition, different viscoelastic rock could be precisely described through adjusting fractional order, which was proved to the conclusions as mentioned previously. The displacement and strain components of circular tunnel based on fractional Kelvin derivative model deduced that(14)ut=a2σ02Gr1-Ea-tτγ,εrt=-a2σ02Gr21-Ea-tτγ,εθt=a2σ02Gr21-Ea-tτγ.Chen [69] took the effects of tunnel face into account, namely, the spatial effect, as shown in Figure 13; at the same time the region within quadruple tunnel radius from tunnel face was defined as the near section, contrarily as the far section. Using fractional derivative Maxwell model and classic Burgers model, respectively, comparative analyses on the displacement changed with time aspects after circular tunnel excavation and support were performed [70], with an attempt to research the deformation characteristics during two stages of decay creep and stable creep of surrounding rock.Figure 13
Tunnel excavation.When we analyze the case of tunnel without support, taking the displacement calculation of tunnel without sidewall as an example, the creep displacement of its far section region is investigated [71], that is, neglecting the influence of tunnel face with λ=1. Thus the viscoelastic solution for sidewall displacement in far section region based on fractional derivative Maxwell model is obtained:(15)Urt=σ0b22r1G+1ηtγΓ1+γ,0≤γ≤1.The parameters of both Burgers model and fractional derivative Maxwell model are firstly confirmed with particle swarm optimization, and then the displacement values of tunnel sidewall in different periods are gained, depicted in Figure14. With the duration time of 1000 days, Figure 14(a) denotes that the description of fractional derivative Maxwell model is essentially equivalent to classical model. But when the duration time comes to 10000 days, the calculation results of two models are close to being consistent only in a certain period of time, rather than the whole procedure, shown in Figure 14(b). In greater detail, it indicates that, compared with the fractional derivative Maxwell model, the calculation results change greatly with time especially in later stage based on Burgers model.Figure 14
The creep displacement of tunnel sidewall without support [69].
(a)
0–1000 days
(b)
0–10000 daysUnder the circumstance of support measures, as shown in Figure15, the impacts of the interaction between tunnel face and support and surrounding rock and support need to be considered simultaneously [72]. It is assumed that time t0 has passed when support structure is performed, in the case of t<t0, representing the state before the support structure is applied, while t=t0 represents the state of performing the support structure. Additionally, the hypothesis is proposed that the process of performing the support structure is completed instantaneously, indicating that no displacement is produced during the installation of support. Thus the viscoelastic solution for sidewall displacement in far section region (λ=1) based on fractional derivative Maxwell model is obtained:(16)Ubt=σ0bKs+2G∗t1+Ks2G,where the stiffness of supporting structure can be estimated by the empirical formula proposed by Brady and Brown [73]:(17)Ks=E′b2-a21+ν′1-2ν′b2+a2.Figure 15
Schematic diagram of support.Subsequently, the displacement values of tunnel sidewall were depicted in Figure16 with the same methods as mentioned previously. It can be known that the change trends of sidewall displacement based on the fractional derivative Maxwell model are basically consistent with Burgers model in the range of 0~1010 days. Furthermore, the deformation results of two kinds of models tend to be convergent with the same final value, caused by the equation limt→∞⁡G(t)=limt→∞⁡1/J(t)=0. Moreover, through the comparative analysis of Figures 14 and 16, it is shown that the displacement value of tunnel sidewall decreases significantly and then tends to be convergent after working of support structure.Figure 16
The creep displacement of tunnel sidewall after supporting [69].


### 4.2. Rheological Settlement of Subgrade

Mentioned in traditional consolidation theory, the soil skeleton is usually regarded as elastomer with pore water pressure [74]. In relation to this, the additional stress formula of embankment is established with an assumption that foundation is considered as an ideal linear elasticity [75, 76]. However, in practical road engineering, the state of embankment foundation is more than extremely complicated with mechanical property of distinct viscoelasticity, which mainly depends on structure composition, temperature, load frequency and load history, and so forth. At the consolidation time scale, processes such as settlement of subgrade can be approached by creep deformation, subjected to a viscous rheology.On the basis of viscoelastic theory and in combination with the characterization of subgrade fill, Wang [77] utilized the fractional derivative Kelvin-Voigt model and simplified the rheological settlement problem of embankment with high fill to be a plane strain problem. Meanwhile, the EDJ-1 type device from Nanjing Soil Instrument Factory was selected, as shown in Figure 17, devoted to undertake shear creep test on compacted soil. The stress-strain relationship was confirmed by the rheological constitutive model subsequently, which was established through the curve cluster of consolidation creep under step loading as shown in Figure 18.Figure 17
The main test device [77].Figure 18
Sketch of curve cluster [77].In greater detail, taking the corresponding load stressσ0 and σ0+Δσ as example, the curves of consolidation creep were obtained from the experimental data treated by linear interpolation. Furthermore, through subtracting correspondingly the curve values of the load stress σ0+Δσ from σ0, the consolidation creep law at a certain depth resulted from the fact that additional load Δσ could be calculated; thus parameters were determined through parameter identification subsequently [78]. In terms of practical subgrade, according to the method as mentioned above, consolidation creep test was performed to each layer of stratified subgrade [79–81], following acquirement of the parameters of model for each layer. Hence, the total settlement of foundation within compressed layer was given by [77](18)S=∑i=1nΔSi=q∑i=1nΔηix,z1E0i+1E1i1-Eγi,1-αiγitτγiHi,where Hi refers to the layer thickness of foundation and its determination method is the same as the layer-wise summation method.


### 4.3. Relevant Loess Researches

Through replacing Newton dashpot in classical model with Abel dashpot, Gen et al. [59] established fractional derivative standard linear model, and its analytical expression of shear stress relaxation was obtained. Furthermore, the test for direct shear relaxation was conducted under the conditions of different positive pressure for Malan loess with different water contents, resulting in the conclusion that the relaxation features of Malan loess could be represented effectively by the proposed model subjected to fitting verification.Taking water damage variableDω of Q2 loess into account, Tang et al. [82] obtained the equations of water damage evolution in terms of various creep parameters accordingly. Through introducing Abel dashpot with variable coefficient, the nonlinear creep model of Q2 loess was established; thus its constitutive equation for creep behavior was given by [82](19)
ε
t
,
ω
=
σ
E
H
∗
+
σ
η
1
∗
t
γ
1
Γ
1
+
γ
1
,
σ
≤
σ
s
,
σ
E
H
∗
+
σ
η
1
∗
t
γ
1
Γ
1
+
γ
1
+
σ
-
σ
s
η
2
t
γ
2
Γ
1
+
γ
2
,
σ
≥
σ
s
,where EH∗=EH[-DεH(ω)], η1∗=η1[1-DεH(ω)], 0<γ1<1, γ2>1.In addition, to verify the correctness and accuracy of the obtained constitutive equation, a series of triaxial creep tests toQ2 loess under different states of water content was undertaken. Meanwhile, based on the test data, the regression analysis for parameters estimation was performed using 1stOpt mathematical software. In the case of water content of 8%, five group curves of model theoretical values and experimental values were depicted in Figure 19, representing different procedure. For a comparison, it exhibits that the established model for Q2 loess is able to characterize effectively each stage of the instantaneous deformation, decay creep, steady creep, and accelerated creep, thus having a bright future for application.Figure 19
Test values and theoretic values for creep of loess sample (moisture content 8%) [80].


## 4.1. Deformation Analysis of Circular Tunnel

As the traditional issue of underground engineering, the stress of surrounding rock mass is redistributed due to tunnel excavation, resulting in occasionally a large delayed deformation that might lead to a delayed failure of structure and rock mass collapse, especially for the tunnel in soft rock [65, 66]. Therefore, the deformation mechanism of surrounding rock should be investigated in the perspective of rheology behavior, since the property of associated time effects.Liu and Zhang [62] proposed to take advantage of fractional derivative Kelvin model to investigate the deformation properties of horizontal circular tunnel [67, 68]. With the analysis of the result, it was shown that fractional derivative model could perform the relaxation properties of viscoelastic rock mass. In addition, different viscoelastic rock could be precisely described through adjusting fractional order, which was proved to the conclusions as mentioned previously. The displacement and strain components of circular tunnel based on fractional Kelvin derivative model deduced that(14)ut=a2σ02Gr1-Ea-tτγ,εrt=-a2σ02Gr21-Ea-tτγ,εθt=a2σ02Gr21-Ea-tτγ.Chen [69] took the effects of tunnel face into account, namely, the spatial effect, as shown in Figure 13; at the same time the region within quadruple tunnel radius from tunnel face was defined as the near section, contrarily as the far section. Using fractional derivative Maxwell model and classic Burgers model, respectively, comparative analyses on the displacement changed with time aspects after circular tunnel excavation and support were performed [70], with an attempt to research the deformation characteristics during two stages of decay creep and stable creep of surrounding rock.Figure 13
Tunnel excavation.When we analyze the case of tunnel without support, taking the displacement calculation of tunnel without sidewall as an example, the creep displacement of its far section region is investigated [71], that is, neglecting the influence of tunnel face with λ=1. Thus the viscoelastic solution for sidewall displacement in far section region based on fractional derivative Maxwell model is obtained:(15)Urt=σ0b22r1G+1ηtγΓ1+γ,0≤γ≤1.The parameters of both Burgers model and fractional derivative Maxwell model are firstly confirmed with particle swarm optimization, and then the displacement values of tunnel sidewall in different periods are gained, depicted in Figure14. With the duration time of 1000 days, Figure 14(a) denotes that the description of fractional derivative Maxwell model is essentially equivalent to classical model. But when the duration time comes to 10000 days, the calculation results of two models are close to being consistent only in a certain period of time, rather than the whole procedure, shown in Figure 14(b). In greater detail, it indicates that, compared with the fractional derivative Maxwell model, the calculation results change greatly with time especially in later stage based on Burgers model.Figure 14
The creep displacement of tunnel sidewall without support [69].
(a)
0–1000 days
(b)
0–10000 daysUnder the circumstance of support measures, as shown in Figure15, the impacts of the interaction between tunnel face and support and surrounding rock and support need to be considered simultaneously [72]. It is assumed that time t0 has passed when support structure is performed, in the case of t<t0, representing the state before the support structure is applied, while t=t0 represents the state of performing the support structure. Additionally, the hypothesis is proposed that the process of performing the support structure is completed instantaneously, indicating that no displacement is produced during the installation of support. Thus the viscoelastic solution for sidewall displacement in far section region (λ=1) based on fractional derivative Maxwell model is obtained:(16)Ubt=σ0bKs+2G∗t1+Ks2G,where the stiffness of supporting structure can be estimated by the empirical formula proposed by Brady and Brown [73]:(17)Ks=E′b2-a21+ν′1-2ν′b2+a2.Figure 15
Schematic diagram of support.Subsequently, the displacement values of tunnel sidewall were depicted in Figure16 with the same methods as mentioned previously. It can be known that the change trends of sidewall displacement based on the fractional derivative Maxwell model are basically consistent with Burgers model in the range of 0~1010 days. Furthermore, the deformation results of two kinds of models tend to be convergent with the same final value, caused by the equation limt→∞⁡G(t)=limt→∞⁡1/J(t)=0. Moreover, through the comparative analysis of Figures 14 and 16, it is shown that the displacement value of tunnel sidewall decreases significantly and then tends to be convergent after working of support structure.Figure 16
The creep displacement of tunnel sidewall after supporting [69].


## 4.2. Rheological Settlement of Subgrade

Mentioned in traditional consolidation theory, the soil skeleton is usually regarded as elastomer with pore water pressure [74]. In relation to this, the additional stress formula of embankment is established with an assumption that foundation is considered as an ideal linear elasticity [75, 76]. However, in practical road engineering, the state of embankment foundation is more than extremely complicated with mechanical property of distinct viscoelasticity, which mainly depends on structure composition, temperature, load frequency and load history, and so forth. At the consolidation time scale, processes such as settlement of subgrade can be approached by creep deformation, subjected to a viscous rheology.On the basis of viscoelastic theory and in combination with the characterization of subgrade fill, Wang [77] utilized the fractional derivative Kelvin-Voigt model and simplified the rheological settlement problem of embankment with high fill to be a plane strain problem. Meanwhile, the EDJ-1 type device from Nanjing Soil Instrument Factory was selected, as shown in Figure 17, devoted to undertake shear creep test on compacted soil. The stress-strain relationship was confirmed by the rheological constitutive model subsequently, which was established through the curve cluster of consolidation creep under step loading as shown in Figure 18.Figure 17
The main test device [77].Figure 18
Sketch of curve cluster [77].In greater detail, taking the corresponding load stressσ0 and σ0+Δσ as example, the curves of consolidation creep were obtained from the experimental data treated by linear interpolation. Furthermore, through subtracting correspondingly the curve values of the load stress σ0+Δσ from σ0, the consolidation creep law at a certain depth resulted from the fact that additional load Δσ could be calculated; thus parameters were determined through parameter identification subsequently [78]. In terms of practical subgrade, according to the method as mentioned above, consolidation creep test was performed to each layer of stratified subgrade [79–81], following acquirement of the parameters of model for each layer. Hence, the total settlement of foundation within compressed layer was given by [77](18)S=∑i=1nΔSi=q∑i=1nΔηix,z1E0i+1E1i1-Eγi,1-αiγitτγiHi,where Hi refers to the layer thickness of foundation and its determination method is the same as the layer-wise summation method.


## 4.3. Relevant Loess Researches

Through replacing Newton dashpot in classical model with Abel dashpot, Gen et al. [59] established fractional derivative standard linear model, and its analytical expression of shear stress relaxation was obtained. Furthermore, the test for direct shear relaxation was conducted under the conditions of different positive pressure for Malan loess with different water contents, resulting in the conclusion that the relaxation features of Malan loess could be represented effectively by the proposed model subjected to fitting verification.Taking water damage variableDω of Q2 loess into account, Tang et al. [82] obtained the equations of water damage evolution in terms of various creep parameters accordingly. Through introducing Abel dashpot with variable coefficient, the nonlinear creep model of Q2 loess was established; thus its constitutive equation for creep behavior was given by [82](19)
ε
t
,
ω
=
σ
E
H
∗
+
σ
η
1
∗
t
γ
1
Γ
1
+
γ
1
,
σ
≤
σ
s
,
σ
E
H
∗
+
σ
η
1
∗
t
γ
1
Γ
1
+
γ
1
+
σ
-
σ
s
η
2
t
γ
2
Γ
1
+
γ
2
,
σ
≥
σ
s
,where EH∗=EH[-DεH(ω)], η1∗=η1[1-DεH(ω)], 0<γ1<1, γ2>1.In addition, to verify the correctness and accuracy of the obtained constitutive equation, a series of triaxial creep tests toQ2 loess under different states of water content was undertaken. Meanwhile, based on the test data, the regression analysis for parameters estimation was performed using 1stOpt mathematical software. In the case of water content of 8%, five group curves of model theoretical values and experimental values were depicted in Figure 19, representing different procedure. For a comparison, it exhibits that the established model for Q2 loess is able to characterize effectively each stage of the instantaneous deformation, decay creep, steady creep, and accelerated creep, thus having a bright future for application.Figure 19
Test values and theoretic values for creep of loess sample (moisture content 8%) [80].


## 5. Summary and Perspectives

For the purpose of describing the rheology phenomena in geotechnical engineering, this paper has provided an up-to-date review of key constitutive models in combination with model structure and acquisition of parameters. Compared with classical constitutive model, the fractional derivative models have shown the extruding advantage of more concise equation, fewer parameters, and higher precision, which also have been improved, amended, and perfected constantly. Currently, with the significant improvement of numerical methods, the parameter of the fractional derivative model that is rather difficult to determine in the pass time has not been the fetters to hinder the research; thus the application of fractional derivative models in geotechnical engineering has been successful with gradually mature technology. We therefore comprehensively reviewed the state-of-the-art of fractional derivative model for practical application, including three major aspects of deformation analysis of circular tunnel, rheological settlement of the subgrade, and relevant loess researches. To some extent, the rheological models based on fractional calculus are developing toward precision, theorization, design of three dimensions, and so forth and have established a better research framework. However, through the application case analysis and summary mentioned above, we can obtain the following conclusions and perspectives.(1) As a powerful mathematical tool, the fractional derivative models have advantage of describing effectively the whole procedure of rheology phenomenon and are endowed with a small number of parameters. Nevertheless, along with the constant refinement of analysis and design in geotechnical engineering, what is crucial and then is proposed during the process of application is a higher demand on parameters precision. Therefore, how to obtain more precise parameters based on experimental data is the target that investigators assiduously seek. Meanwhile, theoretical research ought to be deepened and tends to clear the physical significance and find out its change laws of some parameters such asγ and η.(2) The application case analysis shows that the description of the fractional derivative model has a preferable applicability, by adjusting the fractional derivative order of fractional derivative dashpot. In other words, the fractional derivative model is one of the most effective ways to describe the rheology, which is more flexible for the dynamic deformation. But we need to know that fractional derivative model has many restrictions in the actual application of engineering; for instance, the tunnel under study is assumed to be horizontally circular. Therefore, actual experiences should be accumulated in the processing of engineering, and on this basis the applicability of fractional derivative model ought to keep on improving. Specifically, monitoring in situ test of tunneling should be strengthened, and then the predictive value for corresponding deformation will be gained in combination with the research on accelerated stage of creep. Moreover, in order to ensure the safety of construction site and avoid the occurrence of accident in the long run, the fractional derivative model should be applied to predict the failure phenomena such as collapse and large deformation of tunnel engineering.(3) According to the complicated problems of geotechnical engineering, the results obtained based on classical model due to its limitation deviate largely from the actual field condition. Obviously, the current finite element software programs, which are based on classical constitutive model, are powerless in numerical simulation. Furthermore, owing to the bad potential of operation and cost-effectiveness, laboratory simulation experiment has difficulty in satisfying engineering demands. Considering the complicated mechanics properties of practical engineering, it is necessary to deepen the investigation on three-dimensional fractional derivative model. Moreover, it can be envisaged that fractional derivative model would be more suitable for the analysis of engineering problems in the coming years, especially if the existing numerical software is equipped with fractional derivative models by conducting secondary development.(4) In the field of geotechnical engineering, there were accumulated numerous constitutive models (including classical models and fractional models) with diversified forms. Hence, it is essential to conclude and summarize the existing identification technologies, and then constitutive model database of geotechnical materials should be established in combination with various identification algorithms, devoted to the convenient choices and efficient applications with clear classification of constitutive models. For instance, the fractional derivative Maxwell model should be labeled as one that can only describe the rheology whose initial value tends toward zero, which is not applicable to the case of compacted soil. Finally, taking our results into consideration, a constitutive model based on fractional derivative, which is concise and is able to describe rheology phenomenon accurately, may be incorporated into model database and may be the best choice in the future.



---

*Source: 9183296-2016-05-04.xml*