Analysis of mineral composition

Previous studies (Liu et al., Reference Liu, Ma and Huang2011; Yan et al., Reference Yan, Cui and Gen2016) showed that the mineral components of the hard brittle mud shale of the Longmaxi Formation are primarily quartz, feldspar, carbonate, clay minerals, and pyrite. The proportions of the various minerals in the Longmaxi Formation in areas of the Sichuan Basin (Fig. 4) indicate that the clay minerals and quartz are the primary components of the shale matrix. The mineral compositions at various depths and in different areas vary, however, and such variation is the main cause for differences among the various shales. The amount of quartz present influences significantly the brittleness of the shale, and the clay minerals influence significantly the hydration expansion of the shale.

Fig. 4 Comparison of mineral composition distribution of hard and brittle shale of Longmaxi Formation in different areas of the Sichuan Basin. a Eastern Sichuan; b Western Sichuan; c South Sichuan

Results revealed (Fig. 5) that the mineral components of hard brittle shale were divided into brittle minerals and clay minerals, which is consistent with previous studies. The brittle minerals are quartz, feldspar, carbonate, and a small amount of pyrite with particle sizes of 102 and 1–20 μm for quartz and pyrite, respectively. The clay minerals were illite, chlorite, and mixed-layer illite-smectite with particle sizes averaging <7.6 μm. Montmorillonite and kaolinite were not detected. The clay-mineral content was large, up to 37.3%, while that of quartz was 26.8%, and the brittleness index was 0.549 (Fig. 5a), which is characterized as a typical hard brittle shale. Illite and mixed-layer illite-smectite are the primary clay minerals with a cumulative content of 77.6% (Fig. 5b). Therefore, the expansion capacity of the shale in water was weak (Table 2).

Fig. 5 Mineral composition: a whole-rock mineral analysis; b composition of clay minerals

Table 2 Creep damage corresponding to different moisture contents

Influence of clay-mineral microstructure on hydration swelling

The microstructure of a rock can reveal its macroscopic characteristics. Studying the microstructural characteristics of hard brittle mud shale can help to explain the mechanical behavior of rocks during the hydration process. The typical microstructural characteristics of hard brittle mud shale were obtained by scanning electron microscopy (Fig. 6), which revealed that the clay minerals are closely cemented with non-clay minerals, with a compact structure and high degree of compaction. Several microcracks, micropores, and other fractures developed in the shale, and the clay minerals tended to be distributed and oriented within these spaces. However, the morphology of the microcracks was diverse, with some being connected to micropores. The microcracks enable water to be pushed from the drilling fluid into the rock, whereas the micropores provide a wider area for hydration.

Fig. 6 Scanning electron microscope images of hard brittle mud shale: a clay minerals and non-clay minerals are densely cemented; b clay minerals tend to be oriented; c micro-cracks; d interlayer micro-cracks; e micropores; f micro holes

Influence of physical properties of clay minerals on hydration swelling

(1) Specific surface area: The chemical activity of a rock is related closely to the specific surface area, and hard brittle mud shale with a large specific surface area can undergo rapid chemical and physicochemical reactions with external fluids that intrude into the formation. Experiments performed illustrated that the specific surface area of the hard brittle mud shale was 7.23 cm²/g (Fig. 7b), while the average pore size was 6.578 nm (Fig. 7c), exerting a high capillary force and promoting the hydration reaction. Owing to the crystal layer surface of clay minerals being composed primarily of oxygen atoms or hydroxyl radicals, the larger the specific surface area of rocks, the better the adsorption of water molecules. This indicates that the greater the chemical activity of rocks, the more likely the hydration reaction is to occur. The development of large pore sizes in hard brittle mud shale is related to the clay mineral content, and the specific surface area and pore volume of large pores increase with increase in the clay-mineral content. Further, the specific surface area reflects the degree of hydration of hard brittle mud shale to a certain extent.

Fig. 7 a Cation exchange capacity, b specific surface area, c pore size

(2) Cation exchange capacity: CEC refers to the number of exchangeable cations on the clay-mineral surface, reflecting the strength of the water-adsorption capacity of the rock sample. The larger the CEC, the better is the surface hydration of the mud shale. The experimental results showed that the CEC of the rock samples was relatively high, ranging from 150 to 350 mmol/kg (Fig. 7a). The CEC is related to the degree of dispersion of clay and increases with an increase in the dispersion degree of clay. The results indicate that hard brittle mud shale has certain ability to hydrate, to swell, and to disperse.

The effect of moisture content on the mechanical properties of hard brittle mud shale

In the present study, multiple sets of uniaxial compressive strength experiments of hard brittle mud shale with various moisture contents were performed. The experimental results of the uniaxial compressive strength (UCS) of hard brittle mud shale at various moisture contents are shown in Fig. 8. The maximum UCS of hard brittle mud shale under dry conditions was 87.307 MPa, whereas that under saturated condition was 40.52 MPa corresponding to a softening coefficient of 0.464. This implies that the strength of the rock was obviously reduced with the increase in moisture content. Similarly, the elastic modulus (E) decreased with increasing moisture content from 13.735 to 8.923 GPa. Equation 2 was obtained through non-linear fitting, which is the softening equation for UCS and E of hard brittle mud shale with moisture content, where w is the moisture content of the rock:

(2) $\begin{array}{c}{\sigma }_c=46.738\ast exp\left(-2.319\ast w\right)+40.569\\ E=4.796\ast exp\left(-1.367\ast w\right)+8.939\end{array}$

Fig. 8 Relationship between rock-strength parameters and the moisture content of hard brittle mud shale

The effect of step loading on the creep characteristics of mud shale

Instantaneous elastic deformation was observed during stress application at each level (Fig. 9), and the instantaneous deformation of the rock sample with 3.6% moisture content was the largest. Further, the creep curve with 0% moisture content was below all the curves, and failure occurred only after a long creep time and ultimate strain. Moreover, all rock samples exhibited creep characteristics, implying that the deformation increases with time under a certain stress level. At the same load level, the instantaneous strain increased with increasing moisture content, whereas, under different loading stresses, it exhibited different creep types, which were characterized as non-linear. In the case of medium and low stresses, the creep rate of the rock decreased gradually and tended to be constant, i.e. attenuated creep and stable creep. Under high stress, however, the creep rate of the rock increased, which is shown as accelerated creep. Further, under the last load, the rock sample presented an obviously accelerated creep stage, and the creep curve was slightly upwarped, which immediately led to the failure of the rock sample.

Fig. 9 Creep curves of hard brittle mud shale with various moisture contents under step loading

As is evident from the creep experimental curves of hard brittle mud shale samples with various moisture contents (Fig. 9), under the same external loading conditions, the creep deformation of the dry rock sample was the smallest, that of the rock sample with 0.6% water content was slightly larger, and that of the samples with 3.6% water content was the largest. Thus, the creep deformation of mud shale increased with increasing moisture content.

Isochronous viscoelastic stress–strain characteristics

The axial isochronous stress–strain curves of rock samples with various moisture contents under the same confining pressure were obtained by analyzing the experimental creep data. The stress–strain curves at different times were similar in shape and exhibited a cluster of similar curves (Fig. 10). The isochronous stress–strain curve (Fig. 10) of the rock sample is clearly not a straight line; the strain increased with time, while the deformation modulus decreased gradually with an increase in time, indicating that the creep of mud shale is non-linear. Furthermore, with time, the increase in viscous deformation resulted in isochronous curves deviating from linearity with the strain, ε axis, i.e. the x-axis direction. In addition, with an increase in the stress level, σ, the degree of deviation from linearity of the stress–strain isochronous curve increased, indicating that the degree of non-linearity of the stress–strain curve increased with increase in the stress level. As the elastic deformation in viscoelastic deformation was still large at the early stage of creep, the slope of the curve varied significantly. However, when the creep reached the stable stage, the elastic deformation decreased gradually, the slope of the curve decreased with time, and the variation rate decreased as well. Therefore, the damage was considered to be small at this stage.

Fig. 10. Isochronous stress–strain curves of rock samples with different moisture contents. a w = 0%; b w = 1.8%; c w = 2.6%; d w = 3.6%

The isochronous stress–strain curve of the rock sample with 0.6% moisture content (Fig. 10a) showed that the stress–strain relationship can be regarded as a linear relationship when the stress level is low, whereas a non-linear relationship is observed when the stress level is high, with the non-linearity of the stress–strain curve becoming more obvious with time. The isochronous stress–strain curve of the rock sample with 1.8% moisture content (Fig. 10b) was a cluster of curves. With an increase in time, the curve deviated to the strain axis, i.e. the stress–strain curve softened with increasing time. This is because water entered the rock and reacted with the clay minerals physically and chemically, weakening the overall internal structure of the rock (Liu et al., Reference Liu, Meng and Li2010; Yang et al., Reference Yang, Wang and Li2007). The rock sample with 2.6% moisture content (Fig. 10c) was similar to that of the rock sample with 1.8% moisture content (Fig. 10b). With an increase in time, the isochronous stress–strain curve was closer to the axial strain axis, i.e. the stress–strain curve exhibited softening, implying that the physical and chemical action of water on rocks changed their composition and structure and weakened their overall internal structure.

To study the effect of different moisture contents on viscosity in the creep process, this study defined the average value of the isochronous curve slope in the stable creep stage as a new viscous modulus, E _v. Its value is equal to the average slope of the isochronous curve in the stable creep stage under the same moisture content; this value is used to describe the evolution law of viscous strain during the entire creep process. Based on the viscoelastic modulus calculated for various moisture contents, the viscoelastic modulus softening curve was obtained via fitting (Fig. 11). As shown in Fig. 11, the value of the viscoelastic modulus, E _v, decreased non-linearly with increasing moisture content. The softening equation of E _v obtained by fitting was:

(3) $E_v=33.796\ast exp\left(-1.287\ast w\right)+23.656$

Fig. 11. Viscoelastic modulus softening curve obtained by fitting

Nishihara model

The Nishihara mechanical model is composed of Hooke, Kelvin, and ideal viscoplastic bodies in series (Nishihara, Reference Nishihara1952), as shown in Fig. 12. It can reflect comprehensively the basic elasticity, viscoelasticity, and elastic-plastic properties of rock and describe the attenuated, stable, and accelerated creep of rock. Attenuated, stable, and accelerated creep refer to the stage at which the creep rate decreases, remains constant, and increases, respectively. This is an advantage of the widely used Nishihara model.

Fig. 12. Nishihara mechanical model composed of elastic, viscoelastic, and viscoplastic elements

Because of the existence of plastic elements, when σ < σ _s , the Nishihara model is equivalent to a generalized Kelvin model; however, when σ ≥ σ _s , the model deformation is similar to that of the Burgers model if the plastic resistance part of σ _s is removed. Thus, the constitutive equation of Nishihara’s model is:

(4) $\left(\begin{array}{c}\frac{{\eta }_1}{E_0}\dot{\sigma }+\left(1+\frac{E_1}{E_0}\right)\sigma ={\eta }_1\dot{\epsilon }+E_1\epsilon \left(\sigma <{\sigma }_s\right)\\ \ddot{\sigma }+\left(\frac{E_1}{{\eta }_1}+\frac{E_1}{{\eta }_2}+\frac{E_0}{{\eta }_1}\right)\dot{\sigma }+\frac{E_0{E}_1}{{\eta }_0{\eta }_1}\left(\sigma -{\sigma }_s\right)=E_1\ddot{\epsilon }+\frac{E_0{E}_1}{{\eta }_1}\dot{\epsilon }\left(\sigma \ge {\sigma }_s\right)\end{array}\right)$

The creep equation is further obtained as:

(5) $\left(\begin{array}{c}\epsilon =\frac{\sigma }{E_0}+\frac{\sigma }{E_1}\left(1-e^{-\frac{E_1}{{\eta }_1}t}\right)\left(\sigma <{\sigma }_s\right)\\ \epsilon =\frac{\sigma }{E_0}+\frac{\sigma }{E_1}\left(1-e^{-\frac{E_1}{{\eta }_1}t}\right)+\frac{\sigma -{\sigma }_s}{{\eta }_2}t\left(\sigma \ge {\sigma }_s\right)\end{array}\right)$

where E ₀ is the elastic modulus of the Hooke body, E ₁ and η₁ are the Kelvin elastic modulus and viscosity coefficient, respectively, η₂ is the viscosity coefficient of the ideal viscoplastic body, σ _s is the yield stress of the rock, σ is the load stress, and t is the loading time.

As is evident from Eq. 5, when the stress is less than the yield stress (σ < σ _s ), it degenerates into the generalized Kelvin model. However, when the creep time is sufficiently large, its expression is

(6) $\epsilon =\underset{t\to \infty }{lim}\left(\frac{\sigma }{E_0}+\frac{\sigma }{E_1}\right)$

When the stress is greater than the yield stress (σ ≥ σ _s ) and the time is sufficiently large, the creep expression of the Nishihara model is:

(7) $\epsilon =\underset{t\to \infty }{lim}\left(\frac{\sigma }{E_0}+\frac{\sigma }{E_1}\right)+\frac{\sigma -{\sigma }_s}{{\eta }_2}t$

The strain-time curves before and after stress yielding are shown in Fig. 13.

Fig. 13. Nishihara model creep curve. a σ < σ _s . b σ > σ _s

Using Eqs 6 and 7, the variation in the Nishihara creep model curve before and after the yield stress was determined. When the axial stress is less than the yield stress, the creep rate in the stable creep stage is 0. In contrast, when the load stress is greater than the yield stress, the stable creep rate is the ratio of the load stress minus the yield stress to the viscosity coefficient in the ideal viscoplastic body.

Creep constitutive model considering damage

According to the non-linear creep characteristics of hard brittle mud shale with various moisture contents, a creep constitutive model with aging and hydration damage of rock in a water-bearing state can be established, which can reflect attenuation creep, stable creep, and accelerated creep. As evident from the test results, the creep property of water-bearing rock is related to the moisture content, stress state, and loading time. Based on the Nishihara model, a creep constitutive model that considers aging and hydration damage is proposed here.

Damage variables

Damage is a microstructural change that causes deterioration of solid material properties under a certain load and environment. This type of microstructural change results in the destruction of solid materials to a certain extent; thus, damage to materials under general conditions can be considered as the process of damage accumulation (Zhang et al., Reference Zhang, Yang and Ren2003). The damage variable is used primarily to characterize quantitatively the irreversible changes in the material structure. However, the damage variable cannot be measured directly and accurately; thus, it must be determined reasonably with the aid of certain intermediate variables (Hu et al., Reference Hu, Li and Zhao2002; Liu et al., Reference Liu, He and Zhao2006). In the present study, the damage variables were defined based on the elastic modulus, and the aging damage variables, D _t , of the rock are as follows:

(8) $D_t=1-\frac{E_t}{E_0}$

where E ₀ is the elastic modulus at the beginning and E _t is the elastic modulus at any time, t.

The softening law of the elastic modulus in the instantaneous loading stage and that of the viscosity modulus in the creep process were obtained in the previous section. According to the basic definition of the damage theory, the damage variable, D, is expressed as the ratio of the material defect area to the total effective bearing area of the material:

(9) $D=\frac{S_0-S_w}{S_0}$

where S ₀ is the effective bearing area without damage and S _w is the effective bearing area in the damaged state. Therefore, the damage variable of the strength of the water-adsorbing hard brittle mud shale is:

(10) $D_w=1-\frac{E_w}{E_0}$

Combined with the fitted equation of the elastic modulus varying with moisture content, the damage equation of the elastic modulus can be rewritten as:

(11) $D_w=0.349-0.349\ast exp\left(-1.367\ast w\right)$

Consistent with the concept of defining instantaneous elastic damage, the definition of long-term creep damage is based on the variation in the viscous modulus. The rock was considered to be in a dry state without damage, and its viscous modulus was E _v0. When the moisture content of shale varies, the creep modulus of the rock is attenuated gradually. Subsequently, when the water content reaches saturation, the creep modulus is stabilized. Thus, the creep damage variable can be defined as:

(12) $D_c=1-\frac{E_v}{E_{v0}}$

Combined with the fitting equation of the creep modulus varying with moisture content, the damage equation of the creep modulus becomes:

(13) $D_c=0.588-0.588\ast exp\left(-1.287\ast w\right)$

The creep damage increased with increasing moisture content, and the corresponding damage value was 0 in the dry state. When the hard brittle mud shale sample reached the saturated water content, the damage reached a maximum, approaching 0.588. The damage parameters of the samples with different moisture contents are listed in Table 3.

Table 3 Fitted parameters of non-linear damage creep model

Accelerated creep model

After the rock-creep process reaches the accelerated creep stage, the strain increases non-linearly, and internal crack propagation and damage are often irreversible, viscoplastic, and plastic. Strain growth at this stage is the most important factor controlling rock failure; thus, strain-related parameters can be selected to characterize whether the rock enters the accelerated creep stage. The creep and axial strain rate curves under a load of 56 MPa in the dry state are shown in Fig. 14. The entire creep process is subdivided into three stages: the initial creep lasts for 2 h; this is followed by the stable creep stage; and, starting from 8.4 h, the accelerated creep phase is reached.

Fig. 14. Creep curve and axial strain rate in the dry state

Generally, an elastic body is used to describe the instantaneous elastic behavior of the rock during loading. Under creep load, the deformation property conforms completely to Hooke’s law, considering that the deformation is related only to stress and not time, yielding the corresponding constitutive equation as σ = E ∙ ε. The viscoelastic behavior of rocks is characterized by a viscous body, the mechanics of which conform to the definition of a Newtonian fluid, i.e. stress is proportional to the strain rate, with the constitutive equation $\sigma =\eta \times \frac{d\epsilon }{\mathrm{dt}}=\eta \dot{\epsilon }$ . Therefore, in the process of studying the accelerated creep model, the present study aimed mainly to determine the internal relationship between stress and strain or strain rate and then establish a new element to describe the accelerated creep behavior.

By processing the data of the accelerated creep failure stage in the creep experiment, the strain rate of the hard brittle mud shale in the accelerated creep stage was fitted (Fig. 15), and the strain rate dε/dt increased approximately exponentially with time, which can be fitted with a quadratic function. Further, $d\dot{\epsilon }/\mathrm{dt}$ was derived continually, and finally determined that $\ddot{\epsilon }$ increases linearly with time. Therefore, a new non-linear viscous dashpot model with strain triggering was established to describe the accelerated creep process. The model works only when it begins to enter the accelerated creep stage, which is characterized by a strain greater than ε _N . When the overall strain of the model is less than ε _N , the viscous dashpot does not function, as shown in Fig. 12. Similar to the ideal viscosity definition, a viscous dashpot is defined as one with a stress proportional to $d\ddot{\epsilon }/\mathrm{dt}$ , independent of time, i.e.

(14) $\sigma ={\eta }_N\frac{d\ddot{\epsilon }}{\mathrm{dt}}={\eta }_N\stackrel{⃛}{\epsilon }$

Fig. 15. Variation trend of creep rate in accelerated creep stage

Thus, a new creep model was formed in series with the viscous dashpot model and Nishihara model, referred to here as the improved Nishihara model (Fig. 16). When hard brittle mud shale is dry, the rock is in a state of no damage. Further, when σ < σ _s , the model degenerates into a simple Kelvin model, which can describe the instantaneous deformation and initial creep process of the creep curve. When σ ≥ σ _s and ε < ε _c , the model degenerates into the Nishihara model, which describes the initial creep and steady creep stages in the creep curve. Finally, when σ ≥ σ _s and ε < ε _N , the new viscous dashpot model, taking damage into account, begins to function and can be used to describe the accelerated creep stage in the creep process.

Fig. 16. Improved Nishihara model composed of elastic, viscoelastic, viscoplastic, and viscous dashpot elements

According to the series-parallel relation, the current creep model equation can be obtained through the following theoretical derivation. When ε ≥ ε _N , the rock enters the accelerated creep stage; therefore, the total strain of the improved Nishihara model is:

(15) $\epsilon ={\epsilon }_e+{\epsilon }_{\mathrm{ve}}+{\epsilon }_{\mathrm{vp}}+{\epsilon }_N$

where ε_e, ε_ve, ε_vp, and ε_N are the elastic, viscoelastic, viscoplastic, and accelerated strains, respectively. Its overall creep equation is given by:

(16) $\left(\begin{array}{c}\begin{array}{c}\sigma ={\sigma }_1={\sigma }_2={\sigma }_3={\sigma }_4\\ {\sigma }_e=E_0\left(1-D_t\right)\left(1-D_w\right){\epsilon }_1\end{array}\\ {\sigma }_{\mathrm{ve}}=E_1\left(1-D_t\right)\left(1-D_c\right){\epsilon }_2+{\eta }_1\left(1-D_t\right)\left(1-D_c\right){\dot{\epsilon }}_2\\ \begin{array}{c}{\sigma }_{\mathrm{vp}}={\eta }_2\left(1-D_t\right)\left(1-D_c\right){\dot{\epsilon }}_2\\ {\sigma }_N={\eta }_N{\stackrel{⃛}{\epsilon }}_4\end{array}\end{array}\right)$

When σ ≥ σ_s and ε ≥ ε_N, where σ₁, σ₂, σ₃, σ₄ are the stresses, and when each element in the model is connected in series, σ is the total stress and ε is the total strain. Further, E ₀ and E ₁ are the elastic coefficients of the spring, η₁ and η₂ are the viscosity coefficients of the viscoplastic body, σ_s is the yield stress of creep of hard brittle mud shale, D _w is the hydration damage variable, D _c is the damage variable of the creep viscosity coefficient, and D _t is the aging damage variable.

Results described above indicated that the elastic damage coefficient acted only on the elastic components in the model, whereas the viscosity coefficient acted on the viscous components in the model. By the Laplace transform:

(17) $\epsilon \left(s\right)={\epsilon }_e\left(s\right)+{\epsilon }_{\mathrm{ve}}\left(s\right)+{\epsilon }_{\mathrm{vp}}\left(s\right)+{\epsilon }_N\left(s\right)=\frac{\sigma }{E_0\left(1-D_t\right)\left(1-D_w\right)s}+\frac{\sigma }{\left(E_1\left(1-D_t\right)\left(1-D_c\right)+{\eta }_1\left(1-D_t\right)\left(1-D_c\right)s\right)s}+\frac{\sigma -{\sigma }_s}{{\eta }_2\left(1-D_t\right)\left(1-D_c\right)s^2}+\frac{\sigma }{{\eta }_N{s}^4}$

The creep constitutive equations for different creep stages were obtained using the Laplace inverse transformation. When σ < σ_s, the creep constitutive equation of the damaged generalized Kelvin model can be obtained as

(18) $\sigma +\frac{{\eta }_1^{\ast }}{E_0^{\ast }+E_1^{\ast }}\dot{\sigma }=\frac{E_0^{\ast }{E}_1^{\ast }}{E_0^{\ast }+E_1^{\ast }}\epsilon +\frac{E_0^{\ast }{\eta }_1^{\ast }}{E_0^{\ast }+E_1^{\ast }}\dot{\epsilon }$

where $E_0^{\ast }=E_0\left(1-D_t\right)\left(1-D_w\right),E_1^{\ast }=E_1\left(1-D_t\right)\left(1-D_c\right),{\eta }_1^{\ast }={\eta }_1\left(1-D_t\right)\left(1-D_c\right),{\eta }_2^{\ast }={\eta }_2\left(1-D_t\right)\left(1-D_c\right)$

When σ ≥ σ _s and ε < ε _c , the creep constitutive equation of the damaged Nishihara model is:

(19) $\ddot{\sigma }+\frac{\left(E_0^{\ast }{\eta }_2^{\ast }+E_1^{\ast }{\eta }_2^{\ast }+E_1^{\ast }{\eta }_1^{\ast }\right)}{{\eta }_1^{\ast }{\eta }_2^{\ast }}\dot{\sigma }+\frac{E_0^{\ast }{E}_1^{\ast }}{{\eta }_1^{\ast }{\eta }_2^{\ast }}\sigma =E_1^{\ast }\ddot{\epsilon }+\frac{E_0^{\ast }{E}_1^{\ast }}{{\eta }_1^{\ast }}\dot{\epsilon }$

When σ ≥ σ _s and ε < ε _c , the creep equation of the improved damaged Nishihara model is:

(20) $\epsilon =\frac{\sigma }{E_0\left(1-D_t\right)\left(1-D_c\right)}+\frac{\sigma }{E_1\left(1-D_t\right)\left(1-D_c\right)}\left(1-exp\left(-\frac{E_1\left(1-D_w\right)}{{\eta }_1\left(1-D_c\right)}t\right)\right)+\frac{\sigma -{\sigma }_s}{{\eta }_2\left(1-D_t\right)\left(1-D_c\right)}t+\frac{\sigma }{6{\eta }_N}{\tau }^3$

where τ is accelerated creep time; $\tau =t-t_{\epsilon ={\epsilon }_N}$ .

Determination of parameters and verification of model

Based on the experimental creep data from the Longmaxi Formation hard brittle mud shale with various moisture contents, the time at which creep entered the accelerated creep stage was determined. For example, hard brittle mud shale entered the accelerated creep stage when t = 8.6 h in the dry state, and the corresponding strain was used as the trigger strain value of the viscous dashpot.

First, the damaged viscoelastic element parameters $E_0^{\ast }$ , $E_1^{\ast }$ , ${\eta }_1^{\ast }$ , and ${\eta }_2^{\ast }$ were obtained by fitting the curves of the first two stages. Then, the accelerated creep parameter η _N was obtained by fitting the curves of the accelerated creep stage after substituting the parameters obtained above into the formula. Thereafter, based on the test results, the Levenberg–Marquardt non-linear, least-squares method was adopted in order to obtain the damage creep model parameters under various moisture contents (Xu & Zhang, Reference Xu and Zhang2010; Zhu et al., Reference Zhu, Wang and Wu2008), as shown in Table 3. Comparison of the uniaxial creep experimental curve by step loading with the damaged creep model curve of hard brittle shale with moisture contents of 0, 1.8, 2.6, and 3.6% (Fig. 17) revealed that the experimental curve is consistent with the theoretical curves, demonstrating the accuracy of the model in describing the entire creep process, including the initial creep, stable state creep, and accelerated creep. In addition, it also describes the creep characteristics under various moisture contents, thereby characterizing the creep deformation law of hard brittle mud shale under various moisture contents.

Fig. 17. Comparison of test and simulated creep curves with the established model. a w = 0%; b w = 1.8%; c w = 2.6%; d w = 3.6%

As shown in Table 3, the elastic parameters $E_0^{\ast }$ , $E_1^{\ast }$ , and ${\eta }_1^{\ast }$ obtained by the damage creep model under the same load exhibited little difference. For example, when the stress level was 16 MPa, the range of $E_0^{\ast }$ under various moisture contents was 5.01–6.26 GPa, and the range of $E_1^{\ast }$ was 172.53–209.67 GPa, while ${\eta }_1^{\ast }$ ranged from 775.27 to 877.82 GPa h. Considering the heterogeneity and discreteness of rock materials, the damage creep parameters were statistically analyzed within a reasonable range to achieve the variation rule of each parameter of hard brittle mud shale with various water contents under different creep stress, as shown in Fig. 18. After fitting the mean value of each parameter (Fig. 18), the results showed that E ₀ ^* increased linearly with the creep stress σ, whereas E ₁ ^* and ${\eta }_1^{\ast }$ decreased negatively with σ.

Fig. 18. Relation between creep parameters and stress σ at various moisture contents