Yiran Yang, Xingping Lai, Tao Luo, Kekuo Yuan and Feng Cui

Study on the viscoelastic–viscoplastic model of layered siltstone using creep test and RBF neural network

De Gruyter | 2021

Abstract

Creep is a fundamental time-dependent property of rock. As one of the main surrounding rocks of underground engineering, layered siltstone is governed by creep to a great extent because of special structure. Based on the structural characteristics of layered siltstone, a viscoelastic–viscoplastic model was proposed to simulate and present its creep property. To verify the accuracy of the model, governing equation of the viscoelastic–viscoplastic model was introduced into finite element difference program to simulate a series of creep tests of layered siltstone. Meanwhile, creep tests on layered siltstone were conducted. Numerical simulation results of the viscoelastic–viscoplastic model were compared with creep test data. Mean relative error of creep test data and numerical simulation result was 0.41%. Combined with Lyapunov function, the radial basis function (RBF) neural network trained with creep test data was adopted. Mean relative error of creep test data and RBF neural network data was 0.57%. The results further showed high accuracy and stability of RBF neural network and the viscoelastic–viscoplastic model.

1 Introduction

Safety construction is one of the main tasks to ensure underground engineering. The continuous creep of disturbed surrounding rock can lead to large deformation and even space closure of underground chamber, which must be taken into account in the design and stability analysis of underground engineering. Creep is one of the time-dependent properties of rock. To obtain creep characteristics and instability mechanism of rock, scholars are devoted to relevant research fields. Numerous models applicable to their research objectives have been established. The classical viscoelastic models can be constructed by a series of springs, sliding blocks, and damping elements in series or in parallel. The constitutive law of this kind of model is to establish the variable calculation relation between strain and stress and define creep of rock through analytical expressions [1,2,3,4,5,6]. With the deepening of researches, characteristics of rock in different mechanical environment are discovered, and the classical viscoelastic model shows its limitations. In this case, a series of rheological models have been proposed through in-depth studies and theoretical breakthroughs [7,8,9,10,11], which can simulate creep behaviors of rock and have a wider range of application [12].

Firme et al. [13] used the multi-mechanism deformation creep model (MD model) to explore the creep behavior of Brazilian rock salt, and the results showed good consistency with the creep test data. Miura et al. [14] proposed a creep failure prediction model for hard rock based on microscopic mechanics, and it was believed that the growth of subcritical cracks and the interaction between cracks were the governing factors of creep failure. Leite et al. [15] constructed a rock creep model based on the power law of transient creep and determined the short-term creep parameters of rock. Bozzano et al. [16] constructed a multi-model analysis method based on morphological evolution simulation to study the creep behavior and failure mechanism of rock slope. Brantut et al. [17] studied the microscopic mechanics of brittle creep and proposed a microscopic mechanical model that can describe brittle creep of rock under triaxial stress.

Results of early creep tests have shown that creep cannot be completely recovered even under minimal stress [18,19,20]. Therefore, relevant scholars assumed that the first, second, and third stages of creep curve were all stages of partial reversible strain, and based on this, they proposed models with different characteristics to evaluate the long-term stability of rock [21,22,23,24]. In addition, many scholars have carried out creep experiments on the viscoplastic behavior of rocks, and many beneficial results have been obtained [25,26].

Liu et al. [27] studied the creep behavior and characteristics of saturated rocks under high stress in uniaxial and cyclic loading modes, providing a basis for the deformation control and disaster prevention of saturated rock. Dubey and Gairola [28] studied the influence of structural anisotropy on creep behavior by experimental means and believed that structural anisotropy has a strong control effect on the instantaneous strain, transient strain, steady strain, and accelerated strain of rock salt. Zivaljevic and Tomanovic [29] used the uniaxial creep test method to analyze the creep characteristics and behavior of marl and focused on the influence of pre-consolidation level and loading time on the creep parameters. Pellet and Fabre [30] carried out static, quasi-static, and cyclic creep experiments on sedimentary rock, and results showed that the content of clay particles in sedimentary rock had a significant influence on the creep behavior. Rahimi and Hosseini [31] carried out a triaxial creep experiment on rock salt and studied the influence of confining pressure, deviating stress, and strain rate on the creep behavior of rock salt. Results showed that the strain rate increased with increment of deviating stress and confining pressure. Grgic and Amitrano [32] studied the influence of water saturation on rock creep and explained the important role of microfracture in the creep process by analyzing monitoring data of strain and acoustic emission.

With the enrichment of machine learning theory and the deep optimization of neural network, artificial intelligence is gradually applied to the study of mechanical properties and simulation of large-scale engineering rock body [33,34,35]. Deep learning has become a major part of machine learning in the past few decades because of the development of computing techniques and huge date collection. At present, deep learning has entered a hot stage in the study of mechanical properties of various materials [36], and with the advent of new methods [37,38], the deep learning has a broader application prospect. Mahdaviara et al. [39] adopted a state-of-the-art machine learning algorithm to estimate permeability of carbonate reservoirs. Saghafi et al. [40] established adaptive network-based fuzzy inference system and genetic programming models for the accurate assessment of reservoir oil formation volume factor.

In the past, researches focused on the creep law and mechanism of soft rocks and soils, but a few researches studied creep property of hard rock with layered structures. With increase of depth, hard rock gradually becomes main surrounding rock in the underground engineering [41]. In particular, the surrounding rock with layered structure such as layered siltstone dominates difficulty and progress of construction. Therefore, it is necessary to study the creep property of layered siltstone to provide guidance for the construction projects [42,43]. In this study, we proposed a viscoelastic–viscoplastic model to simulate and present the creep property of layered siltstone. Meanwhile, the radial basis function (RBF) neural network method was adopted, and creep tests on siltstone were conducted. Based on the creep test data, error analysis was conducted to verify scientific nature and accuracy of the proposed model.

2 Theoretical study on viscoelastic–viscoplastic model

Constitutive relation is an important mean to study the mechanical behavior of rock. The breakthrough in constitutive relation theory is the key to establish constitutive model. In this study, we focus on the rheological property of rock, which is a nonlinear function of time. Therefore, it is necessary to distinguish the aging part and the non-aging part in the process of modeling [44,45]. The viscoelastic–viscoplastic model consists of three components, namely, elastoplastic component (EP), viscoplastic component (VP), and viscoelastic component (VE), as shown in Figure 1. The EP consists of a spring and a friction plate. Parameters of the spring are represented by an elastic modulus Ke, and starting stress of the friction plate is represented by σy. The VP is composed of a viscous element and a friction plate element in parallel. Viscosity coefficient of viscous element is represented by ηv, and starting stress of the friction plate is represented by σp. The VE is composed of a series of viscous elements and elastic elements in parallel and series, which can simulate the layered creep characteristic, for example, progressive failure, of layered siltstone. In the VE, viscosity coefficient is ηi (i = 1, 2,..., m), and elasticity coefficient of spring element is Ki (i = 1, 2,..., m). The three components are combined in a series to construct the viscoelastic–viscoplastic model. Many studies have shown that most rock have hardening creep behavior when under stress. The hardening model shows initial hardening envelope after plastic yield starts, which leads to the plastic strain even under low stress. Based on the structural characteristics of layered siltstone, the viscoelastic–viscoplastic model was constructed. The proposed model can simulate elastoplastic characteristic of instantaneous deformation and elastoplastic viscosity deviation of hysteresis deformation; in addition, this model can be well used to simulate stratified deformation and damage characteristics of layered siltstone. ̇

Figure 1 The viscoelastic–viscoplastic model.

Figure 1

The viscoelastic–viscoplastic model.

Total strain of the model is divided into three parts, which can be presented as:

(1) ε ( t ) = ε ep ( t ) + ε ve ( t ) + ε vp ( t ) ,
where ε ep( t) is the elastic strain and irreversible plastic strain of the EP, ε ve( t) is the elastic strain and viscous strain of the VE, and ε vp( t) is the irreversible viscoplastic strain of the VP.

The εep(t) and the εve(t) can be presented as:

(2) ε ep ( t ) = ε k e ( t ) + ε y p ( t ) ,
(3) ε ve ( t ) = ε i e ( t ) + ε i v ( t ) ,
where ε k e ( t ) is the elastic strain of the EP, ε y p ( t ) is the plastic strain generated by friction plate in the EP, ε i e ( t ) is the instantaneous elastic strain of Maxwell cell i, and ε i v ( t ) is the viscous damping strain of Maxwell cell i.

According to the intrinsic property of each component, the stress–strain relationship can be presented as:

(4) σ k = K e ε k e ( t ) ,
(5) σ i = K i ε i e ( t ) ,
(6) σ i v = η i d d t ε i v ( t ) = η i ε i v ( t ) ,
where ε i v ( t ) is the viscous strain rate of Maxwell cell i.

As three components of the viscoelastic–viscoplastic model are in series, total stress of each component is equal. Total stress can be determined by elements in the VE, which can be presented as:

(7) σ ( t ) = i = 1 m K i [ ε ( t ) ε vp ( t ) ε ep ( t ) ε i v ( t ) ] .

Equation (7) indicates that total stress of VE is obtained by stress superposition of each parallel element. In addition, because elastic element and the viscous element in each parallel element are in series, stress of each parallel element is determined by elastic element or viscous element.

Substituting equation (2) into equation (7), and according to equation (4), σ(t) can be obtained as:

(8) σ ( t ) = K e i = 1 m K i K e + i = 1 m K i ε ( t ) ε vp ( t ) ε p ( t ) ε i v ( t ) .

For VE, elastic stress of the parallel element i is equal to viscous stress. Therefore, we have:

(9) σ i v = K i ε ( t ) ε vp ( t ) ε ep ( t ) ε i v ( t ) = η i ε i v ( t ) .

Substituting equation (2) into equation (9), we can get:

(10) σ i v = K i ε ( t ) ε vp ( t ) ε p ( t ) ε i v ( t ) ε k e ( t ) = η i ε i v ( t ) .

As equation (4) is equal to equation (7), it can be written as:

(11) ε k e ( t ) = ε ( t ) ε vp ( t ) ε p ( t ) ε i v ( t ) K e i = 1 m K i + 1 .

By substituting equation (11) into equation (10), it can be obtained that:

(12) η i K i ε i v ( t ) = K e K e + i = 1 m K i ε ( t ) ε vp ( t ) ε p ( t ) ε i v ( t ) .

The below can be obtained after further calculation:

(13) η i K e + i = 1 m K i K i K e ε i v ( t ) + ε i v ( t ) = ε ( t ) ε vp ( t ) ε y p ( t ) .

Assuming τi > 0, we can get:

(14) τ i = η i K e + i = 1 m K i K i K e 0 ,
where τ i is the relaxation time of the ith cell.

Equation (13) can be simplified as:

(15) ε i v ( t ) + 1 τ i ε i v ( t ) = 1 τ i ε ( t ) ε vp ( t ) ε y p ( t ) .

In equation (15), the Colonnetti criterion is adopted. In addition, plastic strain can be transformed into applied distortion strain, which is a known value. Accordingly, the differential solution of equation (15) can be obtained as:

(16) ε i v ( t ) = ε ( t ) ε vp ( t ) ε y p ( t ) t e t t ¯ τ i ε ( t ¯ ) ε vp ( t ¯ ) ε y p ( t ¯ ) d t ¯ .

Substituting equation (16) into equation (8), we can get:

(17) σ ( t ) = t G ( t t ¯ ) ε ( t ¯ ) ε vp ( t ¯ ) ε y p ( t ¯ ) d t ¯ ,
where G( t) is the relaxation function, which can be presented as:
(18) G ( t ) = K e i = 1 m K i K e + i = 1 m K i e t τ i .

Equation (17) shows the relationship between total stress and total strain. The total strain is usually a known value, and the viscoplastic strain rate ε vp ( t ) and elastoplastic strain rate ε y p ( t ¯ ) are unknown values that need to be determined according to the plasticity theory.

Viscoplastic strain is obtained by integrating viscoplastic strain rate within plastic history, and the formula is:

(19) ε vp ( t ) = t ε vp ( t ¯ ) d t ¯ .

According to the plasticity theory, the viscoplastic activation function can be presented as:

(20) Φ ( σ ) = σ σ p .

According to the plasticity theory, when Φ < 0, stress of VP is less than the start value of friction plate σ p , and so the component does not generate plastic strain. When Φ > 0, stress of the VP reaches or exceeds the start value of friction plate, and at this point, viscoplastic strain follows the plastic flow law:

(21) ε vp ( t ) = Φ( σ ) η v Φ( σ ) σ = Φ( σ ) η v sign( σ ),
where Φ( σ ) = ( Φ( σ ) + Φ( σ ) ) / 2 .

The viscoplastic strain rate ε vp ( t ) is determined by equations (19)–(21). To obtain ε y p ( t ¯ ) , the following assumptions are proposed:

  1. (1)

    Stress of the friction element in the EP does not exceed σ y .

    (22) f σ = σ σ y .

    This formula is also known as the elastoplastic yield condition.

  2. (2)

    Assuming σ σ y < 0 , then ε y p ( t ¯ ) = 0 .

  3. (3)

    If σ = σ y > 0 , then ε y p ( t ¯ ) = γ 0 . If σ = σ y < 0 , then ε y p ( t ¯ ) = γ < 0 .

Thus, the elastoplastic flow law can be obtained as:

(23) ε y p ( t ¯ ) = γ sign( σ ) .

Partial derivative of stress on equation (22) can be obtained as follows:

(24) f ̇ = δ f δ σ K e ε ̇ ε y p ( t ¯ ) = δ f δ σ K e ε ̇ γ δ f δ σ K e sign( σ ) .

Since we know that

(25) σ σ = sign( σ ) δ f δ σ = sign( σ ),
and [ sign( σ ) ] 2 = 1 . Therefore, we have:
(26) f ̇ = 0 γ = ε ̇  sign( σ ) .

Substituting equation (26) into equation (23), when f ( σ ) = 0 , f ( σ ) = 0 , we can get:

(27) ε y p ( t ¯ ) = ε ̇ .

Equations (17)–(27) jointly determine and present the viscoelastic–viscoplastic model. The calculation process follows traditional discrete integral method, whereas evolution relationship of change rate of viscoplasticity is solved by the Euler difference method. In particular, when the viscoelastic–viscoplastic model is used for numerical calculation, model parameter identification has been done to provide basis for comparative analysis and verification.

3 Construction and stability proof of RBF neural network

3.1 Introduction to RBF neural network

Neural network has a certain advantage in model research. At present, scholars in related fields have carried out in-depth theoretical research [46,47], of which RBF neural network is a forward network with good performance, which can be presented in the following form [48]:

(28) Y N = f n n Z = i = 1 N ω i s i ( Z ) = W T S ( Z ) .

Structure of the RBF neural network is similar to that of multilayer forward network [49], which is generally composed of input layer, hidden layer, and output layer (Figure 2). The activation function is R i ( x ) = exp x c i 2 2 σ i 2 .

Figure 2 The structure of RBF neural network.

Figure 2

The structure of RBF neural network.

3.2 Construction of the RBF neural network and its stability proof

The RBF neural network has the ability to approximate complex nonlinear functions through simple linear composite mapping, which does not require specific mathematical models. Construction and stability of the RBF neural network are proved as follows.

According to universal approximation property of the neural network, if f(x) is a continuous smooth function about x(k), then an ideal neural network W*TS(x(k)) satisfies the following equation:

(29) f ( x ) = W T S ( x ) + ε ( k ) ,
where W = [ w 1 , w 2 , , w l ] T R l is the ideal neural network weight vector, S ( x ) R l is the activation function of neural network, and ε( k) is an arbitrarily small approximation error.

Suppose that there is an ideal constitutive model:

(30) y ( k + 1 ) = f ( x ) ,
where f( x) is the nonlinear function of the ideal constitutive model, which is approximated by an ideal neural network: f ( x ) = W T S ( x ) + ε ( k ) .

As it is hard to idealize nonlinear function of the constitutive model, a mirror model is constructed as follows:

(31) y ˆ ( k + 1 ) = f ˆ ( x ) = W ˆ T ( k ) S ( x ) ,
where y ˆ ( k ) is the output value of the mirror model, f ˆ ( x ) is the nonlinear function of the mirror model, and W ˆ ( k ) is the estimated value of the ideal neural network weight vector W *.

The model error is:

(32) z ( k + 1 ) = y ˆ ( k + 1 ) y ( k + 1 ) = W ˜ T ( k ) S ( x ) ε ( k ) ,
where z( k) is the model error, and W ˜ ( k ) = W ˆ ( k ) W is the estimation error vector of neural network.

The neural network update rate is:

(33) W ˆ ( k + 1 ) = W ˆ ( k ) r ( S ( x ) z ( k + 1 ) + δ W ˆ ( k ) ) .

It is necessary to select appropriate learning parameters r and δ so that the mirror model can approach the ideal model. The weight W ˆ ( k ) of the neural network is converged to its ideal value W*. The proof process is as follows:

Based on a Lyapunov function method, we can get:

(34) V ( k ) = z 2 ( k ) + 1 r W ˜ T ( k ) W ˜ ( k ) .

Then the one-step forward difference is:

(35) Δ V ( k ) = V ( k + 1 ) V ( k ) = z 2 ( k + 1 ) + 1 r W ˜ T ( k + 1 ) W ˜ ( k + 1 ) z 2 ( k ) 1 r W ˜ T ( k ) W ˜ ( k ) .

Further calculation of equation (33) can be obtained as:

(36) W ˜ ( k + 1 ) = W ˜ ( k ) r ( S ( x ) z ( k + 1 ) + δ W ˆ ( k ) ) .

Substituting equation (36) into equation (35), we can get:

(37) Δ V ( k ) = z 2 ( k + 1 ) z 2 ( k ) 1 r W ˜ T ( k ) W ˜ ( k ) + ( W ˜ ( k ) r ( S ( x ) z ( k + 1 ) + δ W ˆ ( k ) ) T 1 r ( W ˜ ( k ) r ( S ( x ) z ( k + 1 ) + δ W ˆ ( k ) ) ) .

Further calculation of equation (37) can be obtained as:

(38) Δ V ( k ) = z 2 ( k + 1 ) z 2 ( k ) 2 W ˜ ( k ) ( S ( x ) z ( k + 1 ) + δ W ˆ ( k ) ) + r S T ( x ) S ( x ) z 2 ( k + 1 ) + 2 r δ W ˆ T ( k ) S ( x ) z ( k + 1 ) ) + r δ 2 W ˆ T ( k ) W ˆ ( k ) .

Equation (39) can be obtained from equation (32) as follows:

(39) z 2 ( k + 1 ) + z ( k + 1 ) ε ( k ) = W ˜ T ( k ) S ( x ) z ( k + 1 ) .

Substituting equation (39) into equation (38), we can get:

(40) Δ V ( k ) = z 2 ( k + 1 ) z 2 ( k ) 2 z ( k + 1 ) ε ( k ) 2 W ˜ T ( k ) δ W ˆ ( k ) + r S T ( x ) S ( x ) z 2 ( k + 1 ) + 2 r δ W ˆ T ( k ) S ( x ) z ( k + 1 ) + r δ 2 W ˆ T ( k ) W ˆ ( k ) .

Based on the neural network activation function and Young inequality, equation (41) can be obtained (derivations of equations (4144) are added into Appendix):

(41) S T ( x ) S ( x ) l ,
(42) 2 z ( k + 1 ) ε ( k ) r z 2 ( k + 1 ) + 1 r ε 2 ( k ) ,
(43) 2 r δ W ˆ T ( k ) S ( x ) z ( k + 1 ) r l z 2 ( k + 1 ) + r δ 2 W ˆ T ( k ) W ˆ ( k ) ,
(44) 2 δ W ˜ T ( k ) W ˆ ( k ) = δ W T W δ W ˜ T ( k ) W ˜ ( k ) W ˆ T ( k ) W ˆ ( k ) .

Substituting equations (41)–(44) into equation (40), we can get:

(45) Δ V ( k ) z 2 ( k + 1 ) z 2 ( k ) + r z 2 ( k + 1 ) + 1 r ε 2 ( k ) + 2 r l z 2 ( k + 1 ) + 2 r δ 2 W ˆ T ( k ) W ˆ ( k ) + δ W T W δ W ˜ T ( k ) W ˜ ( k ) W ˆ T ( k ) W ˆ ( k ) z 2 ( k ) ( 1 r 2 r l ) z 2 ( k + 1 ) δ W ˜ T ( k ) W ˜ ( k ) ( 1 2 r δ 2 ) W ˆ T ( k ) W ˆ ( k ) + 1 r ε 2 ( k ) + δ W T W .

Selecting appropriate learning parameters r and δ satisfies the following conditions:

(46) 0 < r < 1 1 + 2 l , 0 < δ < 1 2 r .

Then, equation (45) can be converted as follows:

(47) Δ V ( k ) z 2 ( k ) + 1 r ε 2 ( k ) + δ W T W .

When z ( k ) > 1 r ε 2 ( k ) + δ W 2 , Δ V ( k ) < 0 . The error z(k) converges to μ , μ > 1 r ε 2 k + δ W 2 . According to the universal approximation of the neural network, ε(k) can be arbitrarily small, and μ can also be arbitrarily small by selecting a small learning parameter. This shows that the mirror model error z(k) can converge to small neighborhood of the origin, that is, the mirror model can approach the ideal model.

The proof process of the neural network weight W ˆ ( k ) converges to the vicinity of its ideal value W* is as follows:

(48) W ˜ ( k + 1 ) = ( I r δ I r S ( x ) S T ( x ( k ) ) ) W ˜ ( k ) + r ( S ( x ) ε ( k ) δ W ( k ) ) .

Based on equation (46), I r δ I r S ( x ) S T ( x ) < 1 can be obtained, where ε(k) is an arbitrarily small value, and r , S ( x ) are bound, so r S ( x ) ε 1 ( k ) is also arbitrarily small. Supposing that there is a small correction parameter, then r δ W ( k ) is also an arbitrarily small value. Therefore, equation (48) can be rewritten as:

(49) W ˜ ( k + 1 ) = A ( k ) W ˜ ( k ) + d ( k ) ,
where A( k) is the symmetric matrix and A ( k ) < 1 . d( k) is the small perturbation term.

The minimum form of equation (49) is:

(50) W ˜ ( k + 1 ) = A ( k ) W ˜ ( k ) .

According to the property of A(k), the state transition matrix ϕ ( k 1 , k 0 ) = k = k 0 k = k 1 1 A ( k ) < 1 in equation (50), W ˜ ( k ) converges to the origin.

According to the Lyapunov inverse theorem, there are positive definite symmetric matrices P and Q, which satisfy AT(k) PA(k) + P = −Q.

Lyapunov function is constructed for equation (49):

(51) V 2 ( k ) = W ˜ T ( k ) P W ˜ ( k ) .

Then the one-step forward difference of equation (51) is:

(52) Δ V 2 ( k ) = V 2 ( k + 1 ) V 2 ( k ) = W ˜ T ( k + 1 ) P W ˜ ( k + 1 ) W ˜ T P W ˜ ( k ) = λ min ( Q ) W ˜ T ( k ) W ˜ ( k ) λ max ( P ) d T ( k ) d ( k ) 2 λ max A T ( k ) P W ˜ T ( k ) W ˜ ( k ) d T ( k ) d ( k ) .

So, W ˆ ( k ) eventually converges to μ 2 :

(53) μ 2 > d max λ min ( Q ) λ max ( A T ( k ) P ) + λ max 2 A T ( k ) P + λ min ( Q ) λ max ( P ) .

As d max d ( k ) max is relatively small, error vector W ˜ ( k ) of the neural network converges to the small neighborhood of the origin. Which means that estimation vector W ˆ ( k ) of the neural network converges to the vicinity of the neural network ideal vector W*.

It can be proved from the above process that the mirror model can approach the ideal model, and the neural network weight W ˆ ( k ) converges to the ideal value W*. Thus, the RBF neural network can be established, and Lyapunov function is constructed to verify the convergence and stability of the model. Training and learning of the RBF neural network is arranged in Section 5.

4 Experimental setup and data processing

4.1 Rock rheological test

The rock rheological test system is used in the creep test (Figure 3). The calibration curve of the test system shows that relationship between effective stress σe and applied stress σa is σe = 0.244σa – 0.511, which meets the requirements of the study. The diameter and height of rock samples are 50.0 and 100.0 mm, respectively. The dip angle between layered siltstone and horizontal plane is 87° approximately, and the number of rock samples is 6. All rock samples are prepared according to ISRM (2014) requirements. During the test, strain gauges and DD1 cantilever strain sensors were used to measure axial strain and radial strain. In the test, six rock samples (samples A–F) were compressed under static load, and the loading stress was set as σ CI + σ CD 2 ( σ CI is the crack initial stress, σ CD is the crack damage fracture stress). The stress control method was used for loading, and the loading rate was controlled at 0.01 kN/s. In the initial, intermediate, and the final stages of the creep test, 5.0 min, 1.0 h, and 5.0 min were used as time intervals to read and record data. The input mechanical parameters were obtained through creep tests (Table 1).

Figure 3 The rock rheological test system.

Figure 3

The rock rheological test system.

Table 1

Mechanical parameters of rock

Parameter γ (kN/m3) φ (°) C (MPa) υ K (MPa) σt (KPa) GM (MPa) ηM (Pa s) Gk (MPa) ηk (Pa s)
Value 24.6 29.0 0.89 0.28 833 36 522 6.37 × 1014 467 1.28 × 1014

Creep curves of part rock samples are consistent with the typical creep curves in morphology (Figure 4). The failure modes of rock samples are mainly sliding failure along structural plane and shear failure through structural plane. The average axial creep strain of rock samples in the first stage is 0.53% H (H is the height of test samples). In the steady creep stage, the average axial strain rate was 0.05% H/h, and the average creep duration was 42.7 h. The axial strain of the rock sample at the time of failure remained at 0.59% H on average.

Figure 4 Creep test results of layered siltstone. (a) Sample A, (b) sample B, (c) sample C, (d) sample D, (e) sample E, (f) sample F.

Figure 4

Creep test results of layered siltstone. (a) Sample A, (b) sample B, (c) sample C, (d) sample D, (e) sample E, (f) sample F.

4.2 Numerical calculation experiment using the viscoelastic–viscoplastic model

To analyze the performance of the viscoelastic–viscoplastic model in simulating layered siltstone, we adopted numerical simulation program based on two-dimensional finite element method. Numerical samples were 100.0 mm in height and 50.0 mm in diameter with a 2:1 ratio of height to diameter. The mechanical parameters assigned to the numerical samples were obtained through creep tests (Table 2). To avoid the influence of stress reflection, 10 diameters of the model were designed from center to boundary. Bottom displacement of the model was constrained. Stress boundary was set to reflect only the gravity load, and grid was arranged radially outward without influence of support.

Table 2

Model calculation input parameters

Parameter ε0 εK μ1 μ2 μ3 μv K1 K2 K3 Ke KM σp (MPa) σy (MPa) ω
Value 0.3 2.1 38 2465 2691 3968 10.6 8.1 5.3 39.3 27.2 15.5 38.7 0.01

4.3 RBF neural network training and output

A three-layer RBF neural network was adopted in this part. At the input layer, load stress σ and layered structure angle α were used as nodes. The selected parameters will make a difference to creep strain, which shows creep behavior and property of layered siltstone. Therefore, the output layer was set as creep strain ε. MATLAB was used to normalize the stress–strain data obtained from creep test of rock samples. Then the processed data were used as learning sample to train the RBF neural network. Error of the training target was set as 10−4. When the training steps reached 57,600, the model tended to reach a stable state (Figure 5). Figure 6 shows the loss value and accuracy vs the training steps on all dataset. The output data of RBF neural network were reversely normalized with min–max normalization method and Z-score normalization method based on the number of parameters.

Figure 5 Error reduction evolution during RBF neural network training process.

Figure 5

Error reduction evolution during RBF neural network training process.

Figure 6 The loss value and accuracy vs the training steps on all dataset.

Figure 6

The loss value and accuracy vs the training steps on all dataset.

5 Comparative analysis of the creep test, viscoelastic–viscoplastic model, and RBF neural network

To verify scientific accuracy of the viscoelastic–viscoplastic model and stability of the RBF neural network, comparison was made among output data of the viscoelastic–viscoplastic model, RBF neural network, and creep test. We processed the rock strain data from the RBF neural network, creep test, and the viscoelastic–viscoplastic model, and the mean values are 0.643, 0.679, and 0.712, the median values are 0.709, 0.728, and 0.753, respectively. We considered K-fold as the way to split the dataset, and the dataset is split into two sub-datasets, which are training dataset and testing dataset. The latter is not used for training, but for evaluation of the model. Then, to solve the overfitting problem, a portion of the training data is then set aside as validation data to evaluate the training effectiveness of the model. The validation data split from training dataset is an unseen dataset, which can be used to validate the model. R-square and root mean square error for train data, test data, and total database are calculated (Table 3).

Table 3

R-square and root mean square error of data

Parameter Train data Test data Total database
R2 0.951 0.973 0.940
Root mean square error 0.023 0.020 0.036

Figure 7 shows the comparison of rock creep strain predicted by RBF neural network, creep test, and the viscoelastic–viscoplastic model (part of rock samples). The average relative error can be used to evaluate the accuracy of the RBF neural network:

(54) ε = 1 N i = 1 N E i P i E i × 100 ,
where E i is the creep strain of rock samples obtained by creep test, P i is the creep strain of rock obtained by RBF neural network or the viscoelastic–viscoplastic model, and N is the number of data collection points. It can be seen from Figure 7 that average relative error of data from RBF neural network and the viscoelastic–viscoplastic model is 0.57 and 0.41%, respectively, based on the creep test data, which shows a high degree of consistency. This further shows that the RBF neural network has high accuracy and stability and also confirms that the viscoelastic–viscoplastic model is scientific and accurate.

Figure 7 Strain comparison of the RBF neural network, creep test, and the viscoelastic–viscoplastic model. (a) Sample A, (b) sample B, (c) sample C.

Figure 7

Strain comparison of the RBF neural network, creep test, and the viscoelastic–viscoplastic model. (a) Sample A, (b) sample B, (c) sample C.

6 Discussion

In this paper, we studied creep behavior of layered siltstone. Our studies established the viscoelastic–viscoplastic model and trained the RBF neural network to reveal the effect of layered structure on creep behavior of layered siltstone, and the creep test results suggested stability and accuracy of the model and the RBF neural network. Overall, by means of combining constitutive model with RBF neural network and creep test, our studies propose an innovative approach to obtain creep behavior of layered siltstone.

Many researchers devoted themselves to study the mechanical behavior of rock and achieved useful results [13,17,27]. Wang et al. [50] adopted 3D convolutional neural network to realize voxel model of rock samples. Shlyannikov and Tumanov [51] constructed creep damage model for the fracture of the process zone using a stress and ductility-based formulation. In addition, to assess the mechanical behavior of rock salt under different pressure fluctuations, Han et al. [52] built a modified creep model for cyclic characterization of rock salt considering the effects of mean stress, half‑amplitude, and cycle period. These studies used innovative methods to obtain creep behavior and property of rock under different boundary conditions from macro to meso and then to microscale, and part results were applied in practical engineering. Although there are important discoveries revealed by these studies, most of them focus on homogeneous rock. These studies neglect the critical effect of primary structure, for example, layered structure, on creep behavior of rock. Therefore, their stability and accuracy cannot be guaranteed. Our study takes the primary structure of layered siltstone into consideration, part component of the viscoelastic–viscoplastic model is composed of a series of viscous elements and elastic elements in parallel and in series, which can simulate the layered creep characteristic, for example, progressive failure, of layered siltstone. In addition, all training data are obtained from creep test on layered siltstone. Therefore, the established model and RBF neural network can well reflect creep behavior of layered siltstone.

In the future study, we plan to build coupling model of structure effect and frost heave effect to explore creep property of layered siltstone in the low temperature and provide application basis for rock engineering construction in the cold region. In addition, RBF neural network is an alternative way to obtain mechanical behavior. Based on the comparative analysis in Section 5, we ensure that if there are plenty of train sets, its accuracy can be well guaranteed. Accordingly, we need to make progress on a new RBF neural network, which requires less training sets but has higher accuracy in the future study.

As one of the main surrounding rocks of underground engineering, layered siltstone is governed by creep to a great extent because of special structure. Therefore, our study on creep behavior, creep property, and long-term mechanical mechanism has long-term theoretical significance and practical value.

7 Conclusions

  1. (1)

    The viscoelastic–viscoplastic model was established and theoretically analyzed to simulate and present its creep property. Particularly, this model could simulate its progressive failure, elastoplastic characteristic of instantaneous deformation, and elastoplastic viscosity deviation of hysteresis deformation.

  2. (2)

    The RBF neural network was established and optimized using creep test data; then, the Lyapunov function is constructed to prove its stability. Data validation showed that the RBF neural network had good performance in predicting creep property of layered siltstone.

  3. (3)

    The mean relative error of creep test data and the RBF neural network was 0.57%. The mean relative error of creep test data and the viscoelastic–viscoplastic model was 0.41%. This indicated that the RBF neural network had high accuracy and stability and also confirmed scientific nature and accuracy of the viscoelastic–viscoplastic model.

Acknowledgments

The authors are grateful to China Scholarship Council (CSC). A special thanks to Professor Pinnaduwa H. S. W. Kulatilake, Department of Mining and Geological Engineering, The University of Arizona, Tucson, AZ, USA.

References

[1] Golshani A , Oda M , Okui Y , Takemura T , Munkhtogoo E . Numerical simulation of the excavation damaged zone around an opening in brittle rock. Int J Rock Mech Min Sci. 2007 Sept 1;44(6):835–45. Search in Google Scholar

[2] Zhifa Y , Zhiyin W , Luqing Z , Ruiguang Z , Nianxing X . Back-analysis of viscoelastic displacements in a soft rock road tunnel. Int J Rock Mech Min Sci. 2001 April 1;38:331–41. Search in Google Scholar

[3] Kontogianni V , Psimoulis P , Stiros S . What is the contribution of time-dependent deformation in tunnel convergence? Eng Geol. 2006 Feb 1;82:264–7. Search in Google Scholar

[4] Fahimifar A , Tehrani F , Hedayat A , Vakilzadeh A . Analytical solution for the excavation of circular tunnels in a visco-elastic Burger’s material under hydrostatic stress field. Tunn Undergr Sp Technol. 2010 July 1;25:297–304. Search in Google Scholar

[5] Cui F , Yang Y , Lai X , Jia C , Shan P . Experimental study on the effect of advancing speed and stoping time on the energy release of overburden in an upward mining coal working face with a hard roof. Sustainability. 2020;12(1):37. Search in Google Scholar

[6] Cui F , Zhang T , Lai X , Cao J , Shan P . Study on the evolution law of overburden breaking angle under repeated mining and the application of roof pressure relief. Energies. 2019;12(23):4513. Search in Google Scholar

[7] Li X , Qi C , Zhang P . A micro-macro confined compressive fatigue creep failure model in brittle solids. Int J Fatigue. 2019 Sept 1;130:105278. Search in Google Scholar

[8] Haghighat E , Rassouli F , Zoback M , Juanes R . A viscoplastic model of creep in shale. Geophysics. 2020 Sept 2;85:1–47. Search in Google Scholar

[9] Barral M , Chatzigeorgiou G , Meraghni F , Léon R . Homogenization using modified Mori-Tanaka and TFA framework for elastoplastic-viscoelastic-viscoplastic composites: theory and numerical validation. Int J Plast. 2019 Dec 1;127:1–22. Search in Google Scholar

[10] Lim H , Choi H , Zhu F-Y , Webbe T , Yun G . Multiscale damage plasticity modeling and inverse characterization for particulate composites. Mech Mater. 2020 Aug 1;149:103564. Search in Google Scholar

[11] Lin YC , Huang J , Li H-B , Chen D-D . Phase transformation and constitutive models of a hot compressed TC18 titanium alloy in the α + β regime. Vacuum. 2018 Nov 1;157:83–91. Search in Google Scholar

[12] Guan Z , Jiang Y , Tanabashi Y , Huang H . A new rheological model and its application in mountain tunnelling. Tunn Undergr Sp Technol. 2007 Aug 10;23:2008. Search in Google Scholar

[13] Firme P , Roehl D , Romanel C . An assessment of the creep behaviour of Brazilian salt rocks using the multi-mechanism deformation model. Acta Geotech. 2016 Mar 22;11:1445. Search in Google Scholar

[14] Miura K , Okui Y , Horii H . Micromechanics-based prediction of creep failure of hard rock for long-term safety of high-level radioactive waste disposal system. Mech Mater. 2003 Mar 1;35:587–601. Search in Google Scholar

[15] Leite MH , Ladanyi B , Gill D . Determination of creep parameters of rock salt by means of an in situ sharp cone test. Int J Rock Mech Min Sci Geomech Abs. 1993 June 1;30:219–32. Search in Google Scholar

[16] Bozzano F , Della Seta M , Martino S . Time-dependent evolution of rock slopes by a multi-modelling approach. Geomorphology. 2016 April 1;263:113–31. Search in Google Scholar

[17] Brantut N , Baud P , Heap MJ , Meredith P . Micromechanics of brittle creep in rocks. J Geophys Res (Solid Earth). 2012 Aug 1;117:8412. Search in Google Scholar

[18] Cristescu N . Time effects in rock mechanics. Society for experimental mechanics (SEM) – annual conference and exposition on experimental and applied mechanics (vol 2). 2009 Jan 1. Search in Google Scholar

[19] Taheri R , Pak A . Casing failure in salt rock: numerical investigation of its causes. Rock Mech Rock Eng. 2020 June 3;53:3903–18. Search in Google Scholar

[20] Taheri R , Pak A , Shad S , Mehrgini B , Razifar M . Investigation of rock salt layer creep and its effects on casing collapse. Int J Min Sci Technol. 2020 May 1;30:357–65. Search in Google Scholar

[21] Malan DF . Time-dependent behaviour of deep level tabular excavations in hard rock. Rock Mech Rock Eng. 1999 May 1;32(2):123–55. Search in Google Scholar

[22] Parsapour D , Fahimifar A . Semi-analytical solution for time-dependent deformations in swelling rocks around circular tunnels. Geosci J. 2016 Aug 1;20:517–28. Search in Google Scholar

[23] Tran-Manh H , Sulem J , Subrin D , Billaux D . Anisotropic time-dependent modeling of tunnel excavation in squeezing ground. Rock Mech Rock Eng. 2015 Mar 4;48:2301. Search in Google Scholar

[24] Isotton G , Teatini P , Ferronato M , Janna C , Spiezia N , Mantica S , et al. Robust numerical implementation of a 3D rate-dependent model for reservoir geomechanical simulations. Int J Numer Anal Methods Geomech. 2019 Sept 1;43:2752. Search in Google Scholar

[25] Li Y , Xia C . Time-dependent tests on intact rocks in uniaxial compression. Int J Rock Mech Min Sci. 2000 April 1;37:467–75. Search in Google Scholar

[26] Zhang H , Wang Z , Zheng Y , Duan P , Ding S . Study on tri-axial creep experiment and constitutive relation of different rock salt. Saf Sci. 2012 April 1;50:801. Search in Google Scholar

[27] Liu L , Wang G-m , Chen J-h , Yang S . Creep experiment and rheological model of deep saturated rock. Trans Nonferrous Met Soc China. 2013 Feb 1;23:478–83. Search in Google Scholar

[28] Dubey RK , Gairola VK . Influence of structural anisotropy on creep of rocksalt from Simla Himalaya, India: an experimental approach. J Struct Geol. 2008 June 1;30:710–8. Search in Google Scholar

[29] Zivaljevic S , Tomanovic Z . Experimental research of the effects of preconsolidation on the time-dependent deformations – creep of marl. Mech Time-Depend Mater. 2015 Feb 1;19:43. Search in Google Scholar

[30] Pellet F , Fabre G . Damage evaluation with P-wave velocity measurements during uniaxial compression tests on argillaceous rocks. Int J Geomech. 2007 Nov 1;7:431. Search in Google Scholar

[31] Rahimi S , Hosseini M . Laboratory studies of creep behavior on thick-walled hollow cylindrical salt rock specimens. Arab J Geosci. 2014 Sept 24;8:1–9. Search in Google Scholar

[32] Grgic D , Amitrano D . Creep of a porous rock and associated acoustic emission under different hydrous conditions. J Geophys Res Solid Earth. 2009 Sept 1;114:B10201. Search in Google Scholar

[33] Anitescu C , Atroshchenko E , Alajlan N , Rabczuk T . Artificial neural network methods for the solution of second order boundary value problems. Comput Mater Contin. 2019 Jan 1;59:345–59. Search in Google Scholar

[34] Mahdaviara M , Rostami A , Shahbazi K . State-of-the-art modeling permeability of the heterogeneous carbonate oil reservoirs using robust computational approaches. Fuel. 2020 May 15;268:117389. Search in Google Scholar

[35] Chen D-D , Lin YC , Wu F . A design framework for optimizing forming processing parameters based on matrix cellular automaton and neural network-based model predictive control methods. Appl Math Model. 2019 Dec 1;76:918–37. Search in Google Scholar

[36] Hongwei G , Zhuang X , Rabczuk T . A deep collocation method for the bending analysis of kirchhoff plate. Comput Mater Contin. 2019 Jan 1;58:433–56. Search in Google Scholar

[37] Samaniego E , Anitescu C , Goswami S , Nguyen-Thanh VM , Guo H , Hamdia K , et al. An energy approach to the solution of partial differential equations in computational mechanics via machine learning: concepts, implementation and applications. Comput Method Appl Mech Eng. 2020 April 1;362:112790. Search in Google Scholar

[38] Lin YC , Liang Y-J , Chen M-S , Chen X-M . A comparative study on phenomenon and deep belief network models for hot deformation behavior of an Al–Zn–Mg–Cu alloy. Appl Phys A. 2016 Dec 22;123(1):68. Search in Google Scholar

[39] Mahdaviara M , Rostami A , Keivanimehr F , Shahbazi K . Accurate determination of permeability in carbonate reservoirs using gaussian process regression. J Pet Sci Eng. 2021 Jan 1;196:107807. Search in Google Scholar

[40] Saghafi HR , Rostami A , Arabloo M . Evolving new strategies to estimate reservoir oil formation volume factor: smart modeling and correlation development. J Pet Sci Eng. 2019 Oct 1;181:106180. Search in Google Scholar

[41] Yi H , Zhou H , Wang R , Liu D , Ding J . On the relationship between creep strain and permeability of granite: experiment and model investigation. Energies. 2018 Oct 22;11:2859. Search in Google Scholar

[42] Muhammad N , Bresser JHP , Spiers CJ , Peach C . Creep behaviour of bischofite, carnallite and mixed bischofite-carnallite-halite salt rock. Geotectonic Res. 2015 Sept 1;97:15–17. Search in Google Scholar

[43] Hashiba K , Fukui K . Time-dependent behaviors of granite: loading-rate dependence, creep, and relaxation. Rock Mech Rock Eng. 2016 Mar 17;49:2569. Search in Google Scholar

[44] Chen K . Constitutive model of rock triaxial damage based on the rock strength statistics. Int J Damage Mech. 2020 Nov 5;29:105678952092372. Search in Google Scholar

[45] Shi S , Zhang F , Feng D , Tang K . Creep constitutive model for frozen soils based on hardening and damage effects. KSCE J Civ Eng. 2020 Mar 11;24:1146. Search in Google Scholar

[46] Rostami A , Masoudi M , Ghaderi-Ardakani A , Arabloo M , Amani M . Effective thermal conductivity modeling of sandstones: SVM framework analysis. Int J Thermophys. 2016 June 1;37:59. Search in Google Scholar

[47] Lin YC , Li J , Chen M-S , Liu Y-X , Liang Y-J . A deep belief network to predict the hot deformation behavior of a Ni-based superalloy. Neural Comput Appl. 2018 June 1;29(11):1015–23. Search in Google Scholar

[48] Yu Q , Hou Z-S , Bu X , Yu Q . RBFNN-based data-driven predictive iterative learning control for nonaffine nonlinear systems. IEEE Trans Neural Netw Learn Syst. 2019 June 25;31:1–13. Search in Google Scholar

[49] Zhang P , Zhang J , Zhang Z . Design of RBFNN-based adaptive sliding mode control strategy for active rehabilitation robot. IEEE Access. 2020 Aug 24;8:1. Search in Google Scholar

[50] Wang Y , Teng Q , He X , Feng J , Zhang T . CT-image of rock samples super resolution using 3D convolutional neural network. Comput Geosci. 2019 Aug 1;133:104314. Search in Google Scholar

[51] Shlyannikov VN , Tumanov A . The effect of creep damage model formulation on crack path prediction. Frattura Ed Integr Strutt. 2019 April 1;13:77–86. Search in Google Scholar

[52] Han Y , Ma H , Yang C , Zhang N , Daemen J . A modified creep model for cyclic characterization of rock salt considering the effects of the mean stress, half-amplitude and cycle period. Rock Mech Rock Eng. 2020 July 1;53:3223. Search in Google Scholar