# Research on Rock Creep Characteristics Based on the Fractional Calculus Meshless Method

Authors: Gang Peng; Zhanqing Chen; Jiarui Chen
Journal: Advances in Civil Engineering
 (2018)

---

## Abstract

The application of fractional calculus in the rheological problems has been widely accepted. In this study, the constitutive relationship of the generalized Kelvin model based on fractional calculus was studied, and the meshless method was introduced so as to derive a new meshless algorithm formula based on the fractional calculus of the generalized Kelvin model. By using the MTS815.02 hydraulic servo rock mechanics test system, the creep test of mudstones is carried out, and the related data of the creep process were obtained. Based on the generalized Kelvin model of fractional calculus, the related creep parameters of the argillaceous sandstone under compression were fitted. The results showed that the solution of the generalized Kelvin model based on fractional calculus was greatly consistent with the numerical method solution. Meanwhile, the meshless algorithm based on fractional calculus had a favorable stability and accuracy.


---

## Body

## 1. Introduction

The rheological properties of rock are very important to the safety and long-term stability of rock engineering. However, it is difficult to test rheological properties of rock in engineering practice under a limited engineering environment. Therefore, the laboratory rheological test is always the principal method to study the rheological properties of rock. Suitable rheological constitutive models have been established to study the creep characteristics of rock by many scholars. Rheological models based on elastic elements and viscous elements have been proposed, such as Kelvin Model, Burgers Model, and Nishihara Model, and improved models of those models have also been investigated extensively [1–5]. Some scholars have carried out compression experiments on the different stress conditions of frozen rock and established a theoretical model of low-temperature rock creep, which accurately reflects the creep properties of frozen rock [6–10]. Some scholars have used the unified creep model and its corresponding creep equation to study the creep properties of rock under the condition of graded loading [11], which is more accurate in response to the increasing stress in the rock mass engineering. A modified Kelvin model was used to describe attenuation creep and steady creep [12].Fractional calculus that is a branch of mathematical analysis has also been introduced into the study of rheological properties of rocks [13–16], and a large amount of achievements have been obtained. Some scholars have proposed an Abel dashpot on the basis of fractional order, further extended the Nishihara model, and derived a new constitutive model of viscoelastic fractional order [17, 18]. In order to describe the rheological constitutive model of salt rock, the fractional calculus theory was introduced to construct the rheological constitutive model of salt rock [19].In the large-scale rock engineering, the deformation of rock is difficult to measure directly due to the complex working conditions, e.g., water seepage [20, 21], high surrounding pressure [22], and dynamic disturbance [23, 24]. The use of a finite element is an effective method for the simulation. In order to improve the accuracy of numerical simulation, some scholars studied the viscoelastic finite element method based on fractional calculus [25]. However, crack propagation, rock fracture, and excessive deformation in the process of rock rheology lead to the finite element mesh deformation and failure [26, 27]. In order to solve these problems, mesh reconstruction is needed to generate new grid information, but this has a great impact on the accuracy and speed of computation. The meshless method is a new method of the numerical analysis. Based on the node information and the geometric boundary information of the computational domain, it can avoid the mesh constraints partially or thoroughly. The process of forming approximate functions in the meshless method is based on a series of discrete nodes. The problems mentioned above can be solved easily by adding the required nodes, instead of reconstructing the mesh. This method gets away from the concept of the unit, which avoids the above weakness of the finite element [28].According to the above advantages, the basic equations of the generalized Kelvin model were studied based on the elementary theory of fractional calculus in this paper. Based on the data obtained from the compression creep test of the argillaceous sandstone, the meshless method was introduced to study the rock creep, whose formula was derived from the generalized Kelvin model based on fractional calculus. The results showed that the algorithm had a good stability and accuracy.


## 2. Meshless Method for Generalized Kelvin Model Based on Fractional Calculus

The generalized Kelvin model is presented in Figure1, consisting of a Kelvin body and a Hooke body in series. It reflects the gradual change in deformation with time.Figure 1
Generalized Kelvin model.k
1, η1, and ε1 are the elastic modulus, viscosity coefficient, and strain of the Kelvin body, respectively; k2 and ε2 are the elastic modulus and dependent variable, respectively. The fractional form of the stress σ and strain ε relations in the model can be expressed by component parameters:(1)εt=ε1t+ε2t,σ=k1ε1t+η1dαε1tdtα,σ=k2ε2t.Equation (2) can be obtained by Equation (1):(2)εt=ε1t+ε2t,dαε1tdtα+k1η1ε1t=ση1,ε2t=σk2.According to the basic definition of the Riemann–Liouville fractional calculus, the fractional order integral of Riemann–Liouville [29] is given as follows.For any pluralα>0, Equation (3) can be obtained as(3)Ia+aft=1Γα∫atfτdτt−τ1−αt>a,α>0,where Γα=∫0∞tα−1e−tdt is a gamma function.The Riemann–Liouville fractional calculus has strong singularity, which is not easy for engineering and physical modeling. Italy geophysicist Caputo proposed a fractional differential definition with weak singularity [30], which can be expressed as follows:(4)When0<α<1,Dtαft=1Γ1−α∫atf′τdτt−τt>a.Equation (2), Equation (3), and Equation (4) are combined. Set a=k1/η1 and b=σ/η1, considering an initial condition ε1=0 when t=0, Equation (5) can be obtained:(5)dαε1tdtα+aε1t=b.Equation (5) is transformed by Laplace as(6)sαε1s+aε1s=bs.After transferring, it is expressed as(7)ε1s=bssα+a.Equation (7) is transformed by Laplace as(8)ε1t=b∫0tt−sα−1Eα,α−at−sαds,where Eα,αz=∑n=0∞zn/Γnα+α, and then the Mittag-Leffler equation is obtained [29].Solving Equation (8), a,b are brought in:(9)ε1t=ση1∑n=0∞−1nΓαn+α+1η1tαkn+1.Substitute Equation (9) into Equation (2), which can be written as(10)εt=σk2+ση1∑n=0∞−1nΓαn+α+1η1tαk1n+1.In the computational domainΩ, the integral treatment of Equation (10) is carried out:(11)∫ΩεtdV=∫Ωσk2+ση1∑n=0∞−1nΓαn+α+1η1tαk1n+1dV.Then the equation can be changed as(12)Ktut=P,where(13)Kt=∑kIJt,(14)KIJt=∫VBITDk2BJltdV,(15)lt=11+k2/η1∑n=0∞−1n/Γαn+α+1η1tα/k1n+1,(16)P=∫VΦITfdV+∫sσΦITTdV,where Kt is the total stiffness matrix of the equivalent meshless method, KIJt is the stiffness matrix between the meshless node I and the node J, lt is the fractional order form of the generalized Kelvin model, ut is the node displacement vector, and P is the node load vector.


## 3. Experimental Study on Compressive Creep Properties of Rock

Specimens were derived from the unweathered and nonvibratory argillaceous sandstone collected from coal mines. A standard cylindrical sample was produced with a diameter of 55 mm and a height of 100 mm, and the processed rock samples would be maintained the state of natural moisture. During the experiment, indoor temperature and humidity conditions were kept constant as much as possible to reduce the influence of external factors on the experimental results.The three-axis compression rheology test was carried out in the State Key Laboratory for Geomechanics and Deep Underground Engineering, China University of Mining and Technology. The MTS815.02 electrohydraulic servo rock mechanics test system was used, as shown in Figure2. The test system has the advantages of high loading control precision and automatic real-time acquisition of test results. It is equipped with three sets of independent closed-loop servo control systems, such as axial compression, confining pressure, and pore water pressure. It has five basic test functions, namely, uniaxial compression test, conventional triaxial compression test, true triaxial compression test, pore water pressure test, and water penetration test [31–33].Figure 2
MTS815.02 electrohydraulic servo rock mechanics test system.For rock rheological tests, there are usually two loading methods, namely, separate loading method and multistage loading method. The multistage loading method was adopted in this test. The rheological tests were carried out with a confining pressure of 0 MPa and 10 MPa, respectively, the axial pressure of the specimen was changed, and the multistage loading was carried out. The test was terminated when the creep deformation at all levels of load tended to be stable or had a constant steady-state creeping state, and then it was transferred to the next level of loading. The failure samples are shown in Figure3. No RS01 was tested under confining pressure 0 MPa, and No RS02 was tested under confining pressure 10 MPa. The results of the experiment are presented in Figures 4 and 5.Figure 3
Failure specimens after the test.Figure 4
Creep curve: the confining pressure is 0 MPa.Figure 5
Creep curve: the confining pressure is 10 MPa.From the results of the creep test, it is shown that, at low stress levels (when the load is less than the long-term strength of the rock mass), the deformation of the specimen consists of two parts: instantaneous deformation and stable creep deformation. According to the constitutive relation of the generalized Kelvin model of fractional calculus, there are 4 unknown parameters in the model, which are two elastic modulik1, η1 and two viscosity coefficients k2, α. Since the constitutive equations of the model are complex, it is difficult to obtain them directly from the experimental results. According to Equation (10), the parameters of the model can be solved by the least square method, namely,(17)εLSt,k1,k2,η1,α=∑i=1Nεi−εti2,where εi and ti are the strain and corresponding time of the experimental test, respectively, and N is the number of data. The fitting of parameters can be obtained by general nonlinear least squares fitting. When x is small, the Mittag-Leffler function converges faster, but when x is large, the series converges slowly and is difficult to calculate. Therefore, the asymptotic approximation of the Mittag-Leffler function is introduced.(18)Eαx≅−∑n=1∞x−nΓ1−αn.With the method in the literature [19], the results of various parameters can be obtained.Under different stress conditions, the rheological parameters of rock are not constant. Theoretically, there is a corresponding function relationship between the rheological parameters of rock and the stress level. However, because of the complexity in rock engineering practice and many uncertain factors in the rock mechanics test, it is not necessary to emphasize the difference of the rheological parameters under different stresses so as to improve the calculation precision. Therefore, the rheological parameters of rock under different confining pressure are assumed to be constant, and the final rheological parameters of rock can be determined according to the statistical average of the rheological parameters under various loads. The final result of parameter fitting based on the experiment is shown in Table1.Table 1
Parameters determined by fitting analysis based on creep tests.
Confining pressure (MPa)
k
1
GPa
η
1
10
14
Pa
⋅
s
k
2
GPa
α
0
4.43
2.42
6.24
0.90
10
4.64
2.06
5.60
0.91
Average value
4.54
2.24
5.92
0.91


## 4. Experimental Analysis and Numerical Method Verification

According to the parameters determined in Table1, σ=6.35MPa is taken. According to the actual situation of the experiment, the strain curve of 160 h is determined by using Equation (10). The result can be seen as an analytical solution of the experiment.In order to further verify the meshless method of the generalized Kelvin model of fractional calculus, numerical analysis is carried out by using Equation (13). There are many ways to build meshless normal form functions. To facilitate the application of essential boundary conditions, the meshless radial basis interpolation (RPIM) method is used to construct the shape function. In the RPIM formula, due to the use of the Galerkin weak form, the surface force boundary condition has been naturally integrated into the discrete system equations formed during the formulation process. Because the RPIM shape function has the property of the Kronecker δ function, its essential boundary condition can be applied as easily as the finite element method. Since the result determined by using Equation (10) is the strain and the result calculated by using Equation (13) is displacement, the numerical method determining the unique value is converted to strain for ease of comparison.The fractional model solution determined by using Equation (10) and the numerical method solution determined by using Equation (13) are shown in Figure 6, and some results are shown in Table 2.Figure 6
Fractional model solution and numerical method solution.Table 2
Partial data of the fractional model solution and numerical method solution.
Time (h)
0
40
80
120
160
Fractional model (ε10−3)
10.727
22.255
24.047
24.511
24.680
Numerical method (ε10−3)
10.729
22.259
24.052
24.516
24.685As shown in Figure6 and Table 2, the results of the fractional model solution and numerical method solution are very close. The results show that the algorithm based on the generalized Kelvin model of fractional calculus is reliable and stable.A study on creep behavior in the field of rock mechanics and mining science has been performed extensively in the literature, involving theoretical, experimental, and numerical methods. In general, the rheological behavior of rocks is typically characteristic of the viscosity, elasticity, and plasticity. Many creep models have been proposed by researchers to describe the visco-elastoplastic rheological behavior. Although the classical Burgers models are widely adopted in practical engineering due to their simplicity, the fractal derivative model shows the flexibility in the description of the rheological property of soft rock and how to interpret the physical meaning of the fractal derivative order requires further research.


## 5. Conclusion

(1)
In this paper, the common forms of fractional calculus were introduced. The constitutive relation of the generalized Kelvin model based on fractional calculus was studied. This constitutive relation was also introduced into the meshless method, deriving a new formula for numerical calculation.
(2)
Based on the elementary theory of rock rheology, using the MTS815.02 hydraulic servo rock mechanics test system, the creep test of the argillaceous sandstone was carried out, and the related data of the creep process in soft rock were obtained. Based on the generalized Kelvin model of fractional calculus, the related parameters were fitted.
(3)
Using the constitutive relation of the generalized Kelvin model of fractional calculus and the obtained parameters, the analytical curves of argillaceous sandstone creep were solved.
(4)
Using the newly obtained numerical method, the model was calculated and compared with the fractional model results. The results showed that the curves of the two methods were very close with a minor difference. It is shown that the meshless method based on fractional calculus had a good stability and accuracy for the generalized Kelvin model.



---

*Source: 1472840-2018-09-26.xml*