Creep Mechanical Properties and Constitutive Model of Hard Brittle Mud Shale

Abstract

Mud shale is characterized by low strength and strong swelling, with rheological effects and deformation caused by drilling fluid and formation water. Establishing a rheological model to characterize the deformation characteristics is key to solving the problem of wellbore stability. The influence of moisture content on rock strength and creep mechanical properties were studied by means of water absorption, uniaxial compression, and creep tests. The tests showed that with the increase in moisture content, the elastic modulus and strength of hard brittle mud shale decreased. Further, under the same load, the instantaneous strain increased with increasing moisture content. Meanwhile, under various loading stresses, rock creep exhibited non-linear characteristics, which can be divided into three different creep stages: attenuated creep, stable creep, and accelerated creep. Starting with a non-linear viscous dashpot, and then introducing aging degradation and water-bearing weakening effects, based on the water-bearing creep characteristics of hard brittle shale as well as the modeling ideas of the classic component combination model, a new improved creep model based on the Nishihara model was established to describe the characteristics of the accelerated creep stage of hard brittle mud shale with various moisture contents. Subsequently, the Levenberg–Marquardt non-linear, least-squares method was adopted to invert the creep parameters. The results showed that the simulated creep curves achieved by employing the new creep model were consistent with the experimental results, thereby confirming the ability of the new non-linear creep model to provide a theoretical reference for the study of wellbore stability of hard brittle mud shale.

Introduction

In contrast to traditional mud shale which is characterized by highly hydratable expansive properties, the mineral composition of hard brittle mud shale comprises brittle and clay minerals, with the latter being present in relatively small amounts. The clay mineralogy is dominated by illite with non-swelling ability, containing no smectite, only a small amount of a mixed-layer illite-smectite, and a relatively small cation exchange capacity (CEC). Therefore, although hard brittle mud shale is characterized by a certain brittleness and hydration, the hydration performance is not similar to that of traditional expansive mud shale. Moreover, the hydration reaction is not severe. Hard brittle mud shale is a type of fine-grained heterogeneous clay rock, with its mechanical properties determined by the type of clay mineral. Hard brittle shale formed by such clay minerals is among the primary reservoirs of shale gas and oil.

Borehole-wall instability in hard brittle mud shale containing clay minerals is critical. Few previous studies have been conducted on the relationship between clay minerals and the hydration damage to mechanical properties as well as the creep characteristics of hard brittle mud shale. Because of percolation, diffusion, and exchange of ions and water molecules, the ion concentration among the clay particles changes, which causes a change in the potential between the plates. The repulsive force of the electrical double layer changes, therefore, the original mechanical balance between clay particles is broken, and the macro-mechanical properties are affected easily by hydration (Deng et al., 2006; Mao et al., 2010; Zhang et al., 2000). In the petroleum industry, borehole-wall instability is a complex problem in the drilling process. It can be divided roughly into: broken body, plastic body, and mud-shale instabilities, with the third accounting for >90% of the difficulties (Masoud et al., 2016; Mohiuddin & Khan, 2007). Some researchers believe that wellbore stability is the same as mud-shale stability. After excavation, however, the mud-shale wall rock is under stress for an extended period and exhibits the characteristics of slow deformation over time, i.e. creep deformation (Abdulaziz et al., 2021; Al Ajmi & Zimmerman, 2009; Zou et al., 2021). Moreover, the mud-shale wall rock is affected by both the underground fluid environment and the chemical action of the inflow fluid erosion, which has a rheological effect. The creep deformation of the rock is increased significantly by the effect of water on the rock microstructure and mechanical parameters. In addition, mud-shale hydration expansion to produce creep deformation results in the damage to the wall of the bore hole, including shrinkage as well as borehole instability, sticking, squeezing, and destruction of the casing after cementing, which results in considerable economic losses (Abolfazl & Yaser, 2019; Liu et al., 2019; Liu & Zhu, 2018). Study of the creep characteristics of hard brittle mud shale can, therefore, aid in establishing the rheological law of hard brittle mud shale containing water while providing more comprehensive and reliable data for the maintenance of stable bore-hole walls in hard brittle mud shale.

The influence of water on the mechanical properties of rocks has always been a focus of research.

The impact of microstructure and clay-mineral composition on mechanical properties was studied by Shi (2011) using scanning electron microscopy (SEM). The influence of microscopic geological characteristics, i.e. the mineral composition of hard brittle mud shale on wellbore stability was studied by Zhao et al. (2007) using SEM, X-ray-diffraction (XRD) analysis, and mercury injection. The relationship between clay-mineral composition and the mechanical strength of hard brittle shale was established by Amine et al. (2018) to evaluate the ‘fracability’ of the rock for shale-gas exploration purposes. Tests were conducted by Hawkins and McConnell (1992) on 35 groups of sandstones and those authors found a negative exponential relationship between moisture content and strength, i.e. the strength decreased gradually with increasing moisture content. Frolova et al. (2021) studied the influence of water on the time-dependent deformation characteristics of granite and concluded that the time-dependent deformation of dry granite increased significantly when it was exposed to water. Further, through experimental studies, Fereidooni (2016) observed that the creep failure time of saturated granite was shortened by three orders of magnitude compared to that of the natural state. Aydan et al. (2014) studied the influence of water saturation on the strength, elastic modulus, crack-initiation stress, and damage-stress threshold of siltstone. Zhou et al. (2005) conducted tests on soft rocks with different water-absorption times and found that the quantitative characterization relationship of compressive strength, tensile strength, and shear-strength changes followed the rule of exponential change. Ju and Huang (2016) studied the creep characteristics of red mudstone with various water contents under triaxial compression, and the results showed that the rates of deformation during the initial, steady, and accelerated creep stages of red mudstone increased with increasing water content. Although many studies have been conducted on the influence of water on the mechanical properties of rocks, they have mostly been focused on the strength and creep tests under the two extreme conditions of dry or saturated conditions, and particularly less on the mechanical properties of creep of hard brittle mud shale composed mainly of illite and montmorillonite under various water contents.

Several indoor creep experiments have been conducted to study the variation law of stress–strain relationship with time during the creep process and to explore the mechanism of creep failure (Deng et al., 2016; Liu et al., 2020; Tomanovic et al., 2006). The theories used to study rock creep constitutive models include the empirical element model, the viscoplastic model, and the damage model. A fourth model, the component theoretical model, describes the creep behavior by combining basic components (including the Hooke or elastic body (H), the Newton or viscosity body (N), and Saint-Venant or plastic body (St. V)) in series or parallel. Currently, the Kelvin, Burgers, and Nishihara models are used widely. The Nishihara model comprises a Hooke body combined with one of the others to form either a viscoelastic body (H/N) or a viscoplastic body (H/St. V) (Cao et al., 2012; Qi, 2008), which offers the advantage of a more comprehensive expression of creep deformation and is frequently used in creep experiments and numerical simulations, more often than the Kelvin and Burgers models. Viscoelastic body (H/N) refers to those wherein the Hooke and Newton bodies are in series, whereas viscoplastic bodies (H/St. V) are those wherein the Newton and Saint-Venant bodies are in parallel. Singh and Mitchell (1968) used the Nishihara model to conduct a viscoelastoplastic rheological analysis of deep plagioclase, and the fitted curves were consistent with the creep and relaxation curves obtained via the test. Grgic and Amitrano (2009) regarded the creep characteristics of coal as a viscoelastoplastic body that adheres to the Nishihara model, and obtained a model that can predict the deformation of a coal pillar. In terms of water-damage creep, Chen et al. (2019) introduced transient elastic damage and long-term creep-damage variables into the Burgers model, and established a creep constitutive equation considering the water damage effect. Peng et al. (2020) conducted rheological tests on silty mudstone with various moisture contents under various confining pressures and analyzed the influence of moisture content on its creep variation and long-term strength. Guessous et al. (1987) obtained experimentally the elastic modulus and viscosity coefficient of rock with gradually increasing moisture content. The research above on the Nishihara model and moisture content creep focused primarily on the steady-state creep stage; focus on the accelerated creep stage only is insufficient. Moreover, the traditional Nishihara model cannot reflect the characteristics of non-linear, accelerated creep.

Based on the analysis above, the purpose of the current study was, therefore, to consider hard brittle mud shale from the Longmaxi Formation in the western Sichuan Basin as an example for conducting conventional uniaxial compression strength and creep tests in a continuous water environment for specimens with various moisture contents. Subsequently, based on the creep mechanics test, the objective was to improve the Nishihara model by introducing the effects of aging degradation and water reduction, starting from the perspective of a non-linear viscous dashpot. A new component model, thus, was envisioned in which the creep acceleration behavior of hard brittle mud shale could be predicted and verified. The results of this study are expected to provide theoretical guidance for shale wall rock and safe production under complex underground groundwater and complex in situ stress conditions.

Materials and Methods

Samples

The Sichuan Basin is located west of the Yangtze Platform, with the Mizangshan uplift-Dabashan fold belt in the north, the Emeishan-Liangshan thrust belt in the south, the Longmenshan orogenic belt in the west, and the Hunan-Guizhian-Hubei thrust belt in the east (Fig. 1). During the Longmaxi Period of the Early Silurian, the margins of the Mid–Upper Yangtze continental blocks were in the process of compressional fold orogeny, and areas such as the Western Sichuan-Kang Dian ancient land, Hannan, and Middle Sichuan uplift expanded increasingly. Subsequently, the Xuefeng uplift and the Miaoling uplift were connected to form the Yunnan-Guizhou-Guangxi uplift belt, producing a tectonic pattern of multiple uplifts and bends. Thus, the Longmaxi period of the Early Silurian was the most intense compression period in southern China. Simultaneously, under the influence of global transgression, a set of black, hard, and brittle mud shales with a wide distribution and large thickness was formed during the Longmaxi period (Fig. 1).

Fig. 1 Regional geological map and stratigraphic sequence of the Longmaxi Formation in western Sichuan, China

The hard brittle mud shale examined in this study originated from the Longmaxi Formation in the western Sichuan Basin, as shown in Fig. 2. It is a typical hard brittle mud shale with stiff texture, gray and black-gray color, and parallel bedding. However, when drilling cores using the hydraulic drilling method, the rock samples are easily disconnected along the bedding plane, rendering the achievement of the required length for the experiment a challenge; thus, preparing rock samples is extremely difficult. The rock was frozen and then drilled with kerosene, with the drilling rate controlled during the drilling process to obtain a standard rock sample that satisfied the requirements of the experiment. Thereafter, according to the standards recommended by the International Society of Rock Mechanics, the core was processed into a cylindrical standard sample with dimensions of φ25 mm×50 mm.

Fig. 2 Samples of hard brittle mud shale from Longmaxi Formation: a black, no bedding, 2558.0 m; b light black, bedding, 2562.3 m; c gray, with rhythmic silt bedding, 2568.5 m; d gray, silt bedding, 2573.6 m; e gray black, argillaceous bedding, 2578.5 m; f gray, silt bedding, 2582.6 m; g gray, filled with lamellar calcite, 2588.7 m; h gray black, bedding, 2596.5 m

Test Methods

Determination of mineral composition

X-ray diffraction analysis of the rock and clay-mineral components of samples was performed employing the D8 DISCOVER (Bruker, Billerica, Massachusetts, USA) X-ray diffractometer at a laboratory temperature of 22°C and relative humidity (RH) of 30%. Quantitative analysis of the mineral composition was based on SY/T5163-2010 (Zhou et al., 2017). The clay minerals were extracted using the suspension method, and 10 g of the sample was poured into 50 mL of distilled water, stirred evenly, and maintained for 6–8 h. The clay minerals were suspended in water while the other substances settled. Thereafter, distilled water was extracted and centrifuged to obtain the clay minerals. The XRD of the clay mineral components was performed in three steps: (1) the samples were air-dried for ~12 h, and then the natural air-dried oriented samples (N pieces) were obtained and placed in the diffractometer for analysis; (2) the samples were placed in an incubator with ethylene glycol for more than 7 h (40–50°C); and (3) the sample was heated for 2 h in a muffle furnace at 550°C. XRD patterns were then obtained from each sample thus treated.

Determination of physical and chemical properties

The method for obtaining the cation exchange capacity (CEC) of the hard brittle shale was the industry standard ‘Test Method for Physical and Chemical Properties of Shale (SY/T 5613-2000)’ (Zhou et al., 2017). The rock samples were made into powder with particle sizes of <0.15 mm. Thereafter, the shale slurry was prepared using distilled water, hydrogen peroxide solution, and sulfuric acid solution. Further, the CEC was calculated using methylene blue titration and the specific surface area was measured using a specific surface area analyzer.

Test procedure of uniaxial compression and creep

The RLW-2000 triaxial rheometer controlled by a microcomputer is a multifunction rheometer developed by Changchun Chaoyang Testing Machine Factory (China) which can perform uniaxial and triaxial compression and creep seepage tests of rocks. The control system is an ideal controller with an all-digital servo controller, facilitating closed-loop control of the test force, deformation, displacement, and smooth switching of the three control modes. The rheometer and extensometer are shown in Fig. 3.

Fig. 3 Experimental instruments. a RLW-2000 rheometer; b Extensometer

To explore the effect of the water content on the mechanical parameters of hard brittle mud shale and to determine reasonable confining pressure levels for creep tests, uniaxial compression tests were performed on ten groups of rock samples. The procedure was as follows: the specimens were baked in a drying oven at 105°C for 24 h, cooled to room temperature, selected dried specimens were soaked in water for a certain time, and then the surface was sealed. Conventional uniaxial compressive strength (UCS) tests were then performed on samples with various moisture contents. Displacement-controlled loading was adopted and the loading rate was set to 0.005 mm/s. Based on the various soaking times, ten groups with different moisture contents were established: dry (0%), 1 h (0.46%), 2 h (0.6%), 4 h (0.78%), 18 h (1%), 22 h (1.8%), 26 h (2%), 28 h (2.6%), 32 h (3.4%), and 36 h (3.6%).

The two most commonly used creep loading methods are single-stage loading and stepwise incremental loading. The latter was applied in the current experiment, utilizing the electro-hydraulic servo test system mentioned above. The UCS obtained from that experiment was then used to compare the creep response of the various samples at any given stress state. Considering the low UCS of the saturated sample observed previously (unpublished data of the authors), the initial creep stress was set as 16 MPa. The detailed steps of these experiments were as follows. An axial preload was applied to the specimen and thereafter the load was increased slowly to the rated load (16 MPa). Next, the axial stress was maintained constant, the instantaneous strain was recorded, and sample displacement or deformation was observed continuously. Following the creep process, the next stage was entered, in which the creep time of each stage was maintained for 24 h or until creep failure of the sample occurred. In the initial stage of the experiment, the creep rate changed significantly and the time interval of recording was smaller. After a certain period of time, when the rate of creep strain tended to be stable, it was recorded every 1–2 h. Further, in terms of experimental data processing, owing to the rock creep being non-linear, the data were processed according to Chen’s loading method (Chen & Kang, 1991). However, to obtain the curve of the accelerated failure stage of creep and to avoid brittle failure of the specimen in the process of applying excessive instantaneous load, the stress increment gradient was adjusted at the later stage of loading and reduced to half of the initial design value until the accelerated creep sample was destroyed under the action of the last level of load.

Water adsorption test and determination of moisture content

To obtain the creep mechanics, deformation, and failure characteristics of hard brittle mud shale under various water contents and stress states, water adsorption tests were performed before the strength and creep tests. First, the specimen was baked at 105°C for 24 h, following which the mud shale sample was considered to be in a dry state (moisture content 0%). Following cooling to room temperature, the dry specimen was weighed and immersed immediately in a container filled with distilled water. The specimen was removed every 30 min within the first 10 h of immersion, the surface moisture was wiped off with a wet cloth, and the specimen was weighed on a high-precision (resolution of 0.001 g) balance every hour. Subsequently, the moisture content, w, of the rock after soaking for a certain time was calculated using the following formula:

(1) $w=\frac{\left(M_t-M_0\right)}{M_0}\times 100\%$

where M _t is the mass of the soaked rock specimen at time t and M ₀ is the mass of the dry rock specimen. Five different moisture contents were selected: 0% (dry state), 0.6%, 1.8%, 2.6%, and 3.6% (saturated state) in this experiment, and representative results are given in Table 1.

Table 1 Experimental sample parameters

Results and Discussion

Factors influencing the swelling of hard brittle mud shale

Analysis of mineral composition

Previous studies (Liu et al., 2011; Yan et al., 2016) 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., 2010; Yang et al., 2007). 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

Creep constitutive model considering aging and hydration damage

Nishihara model

The Nishihara mechanical model is composed of Hooke, Kelvin, and ideal viscoplastic bodies in series (Nishihara, 1952), 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., 2003). 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., 2002; Liu et al., 2006). 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, 2010; Zhu et al., 2008), 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

Discussion

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). 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 the clay minerals and increases with an increase in dispersion. It illustrated that hard brittle mud shale had a certain capacity for hydration, swelling, and dispersion.

The Nishihara model can describe the elastic, viscoelastic, and elastic-plastic properties of mud shale. The Nishihara model was adopted by Cao et al. (2012), Hoxha et al. (2005), and Sun et al. (2015) for deep-inclined long amphibolite, coal rock, and gypsum mines, respectively, to perform a viscoelastic-plastic rheological analysis, simultaneously receiving the creep and relaxation deformation curves. However, the weakening effect of water on rock-mechanics parameters was not considered, which cannot reflect the influence of moisture content on creep and characterize the accelerating creep property. Although Song et al. (2015) mentioned that water content would affect creep parameters in the process of studying the creep characteristics of carbonaceous slate, he still adopted the traditional Nishihara model, resulting in the effect of water content on viscoelastic parameters not being reflected. In the current study, based on the Nishihara model and damage theory, aging degradation and the water-bearing weakening effect were introduced, and the viscoelastic parameter degradation characteristics in the creep process were considered comprehensively, which can reflect the hydration damage of hard brittle mud shale. In contrast to other existing studies, the present study conducted a statistical analysis on the elastic parameters of different moisture contents, and found an exponential relationship between the elastic parameters and moisture content. Consequently, the concept of the viscoelastic modulus representing this elastic parameter influenced by the weakening effect of water was proposed to reflect the non-linear variation trend of the elastic modulus with increase in the moisture content (Fig. 8).

In terms of characterizing the stage of accelerated creep, based on an empirical model, Yang et al. (2007) considered the influence of confining pressure and axial pressure on the time-dependent stress–strain relationship and proposed an exponential equation to describe initial creep and stable creep. Although this empirical model can accurately express the rheological properties of specific rock materials, the physical meaning of the model parameters is unclear. Based on the Nishihara model, Qi (2008) introduced a new non-linear element to characterize the accelerated creep characteristics but did not consider the effects of aging degradation and water weakening effects on the viscoelastic parameters. In the present study, the strain rate of hard brittle mud shale in the accelerated creep stage was fitted with the data from the accelerated creep failure stage (Fig. 11). The strain rate dε/dt was varied exponentially with time and the growth relationship was fitted to a quadratic function. Herein, $d\dot{\epsilon }/\mathrm{dt}$ for the creep strain rate was derived continuously, i.e. $\ddot{\epsilon }$ exhibited a linear growth relationship with time. Therefore, a new non-linear sticky pot model with strain triggering was established to describe the accelerated creep process. Moreover, the model operated only when the strain was greater than ε _N and the accelerated creep stage was entered. When the overall strain of the model was less than ε _N , the viscosity pot did not work. Therefore, the constitutive model proposed in this study reflects the non-linear creep characteristics, considering both the aging degradation and water weakening effects and the accelerated creep characteristics.

However, although certain excellent fitted effects have been achieved with the new model, more experiments are required to verify further its popularization and accuracy due to the limited number of samples and regional geology.

Conclusions

The uniaxial compressive strength tests of several groups of hard brittle mud shale with various moisture contents showed that the rock strength had been damaged with increasing moisture content. Similarly, the elastic modulus decreases with increase in the water content.

1. (1) Hard brittle mud shale exhibited creep characteristics. Moreover, under a certain stress level, the deformation increased with time. In addition, under the same load level, the instantaneous strain increased with increasing moisture content, while under various loading stresses, the creep of the rock showed various types of creep, i.e. attenuated, stable, and accelerated creep stages.

2. (2) The creep of hard brittle mud shale exhibited non-linear characteristics. The isochronous stress–strain curve showed that the stress–strain relationship can be regarded as a linear relationship when the stress level is low; however, a non-linear relationship was obvious when the stress level was high, and became more obvious with increase in time.

3. (3) By introducing the aging and hydration damage effects and referring to the modeling ideas of the classical element combination model, starting from the perspective of a non-linear viscous dashpot, the Nishihara model was improved. The creep model obtained by considering the aging and hydration damage can describe the accelerated creep characteristics of hard brittle mud shale. Moreover, the creep curve simulated by the model was consistent with the experimental results, indicating that the model is correct and can be used for actual production.