Introduction

Roughly 75% of the Earth's surface is covered by metamorphic and sedimentary rocks, which often play a significant role in rock engineering projects (Li et al. 2021). These include various applications such as shallow civil projects, deep mining, waste disposal, and oil and gas projects (Chiarelli et al. 2003; Corkum and Martin 2007; Han et al. 2020; Li et al. 2022; Ma et al. 2018; Meier et al. 2014; Stjern et al. 2003; Wang et al. 2024; Yin et al. 2020; Zhang et al. 2013; Zheng et al. 2019). In those projects, layered rocks are very common such as shale, slate, schist, sandstone, and coal (Gholami and Rasouli 2013). The presence of weak bedding planes in these rocks can significantly impact engineering safety, motivating numerous researchers investigate the failure mechanisms of the rock mass with weakness planes (Heng et al. 2021; Li et al. 2020a, 2021; Shen et al. 2021; Simanjuntak et al. 2016; Vu et al. 2012; Wang et al. 2024, 2018; Yang et al. 2023). In terms of experimental studies, scholars extracted and prepared specimens with bedding planes through two methodologies. The first one obtained sample directly through coring in the field (Corkum and Martin 2007; Duan et al. 2022; Gholami and Rasouli 2013; Kou et al. 2022; Li et al. 2020a; Saeidi et al. 2014, 2013; Tan et al. 2014). The other one was lab-based through casting the layered rocks with synthetic material followed by coring the samples from different angles to achieve various bedding angles (Aliabadian et al. 2019; Li et al. 2021; Shen et al. 2021; Yang et al. 2019; Yao et al. 2023; Yin and Yang 2020). Once appropriate samples are obtained, uniaxial or triaxial compression tests (Duan et al. 2022; Li et al. 2021; Saeidi et al. 2014, 2013; Shen et al. 2021; Tien et al. 2006; Yang et al. 2019; Yin and Yang 2020), Brazilian splitting tests (Aliabadian et al. 2019; Li et al. 2020a; Ma et al. 2018; Tan et al. 2014), three-point bending tests (Lei et al. 2021; Zuo et al. 2020), torsion shear tests (Togashi et al. 2018), and quasi-static loading tests (Yang et al. 2023) are normally conducted to study the influences of various factors on the mechanical properties of the rock including specimen size, bedding angle, bedding spacing, bedding morphology, and confining pressure. On the flip side, the experiments are always costly and time consuming and as such, numerical simulations including particle flow code numerical simulation methods (PFC) (Huan et al. 2019), discrete element numerical simulation methods (DEM) (Li et al. 2020b; Park and Min 2015; Shang et al. 2018), and realistic failure process analysis (RFPA) (Wang et al. 2012) have been employed to study the deformation and failure characteristics of layered rocks and their impact on tunnel engineering. In addition, for a deeper understanding of the deformation and failure patterns of layered rocks to better guide engineering practice, the constitutive relationships of layered rocks and the mechanical responses of layered rock masses with opening have also received significant attentions (Heng et al. 2021; Kou et al. 2022; Li et al. 2020a, 2021; Shen et al. 2021; Simanjuntak et al. 2016; Togashi et al. 2018; Vu et al. 2012; Wang et al. 2018; Zhang et al. 2017; Zhang and Sun 2011; Zhao et al. 2022). The analytical models would require five fundamental elastic mechanical parameters including the elastic modulus and Poisson's ratio in the transverse isotropic plane, as well as the elastic modulus, Poisson's ratio, and shear modulus perpendicular to the transverse isotropic plane (Chiarelli et al. 2003; Li et al. 2020a, 2021; Wang et al. 2018). To obtain these five parameters, it's necessary to measure the relevant mechanical parameters from at least two different bedding angles for calculation (Cho et al. 2012). Also, strain gauges need to be installed on the specimens during the testing process to measure strains at different locations. To date, many researchers have developed constitutive models for investigating the influence of bedding angles and confining pressures on the mechanical properties, deformation, and failure mechanisms of layered rocks (Kou et al. 2022; Shen et al. 2021; Wang et al. 2018). The common limitation is these studies are only qualitative based. Moreover, the process of deriving constitutive relationships for layered rocks is cumbersome, especially in terms of the time and effort required for installing strain gauges on the surface of specimens.

Therefore, this study conducted uniaxial and triaxial compression tests on two groups of rocks with different bedding angles under various confining pressures. The experimental results were then analyzed through orthogonal matrix to determine the weight of the influences of various factors including confining pressure, lithology, and bedding angles on the deformation of layered rocks. Consequently, a coupled statistical damage constitutive model for layered rock was proposed, considering the combined effects of confining pressure, lithology, and bedding angles. The model theoretical results were compared with the experimental results, validating the rationality of the constitutive model. The new findings of this study contribute to a quantitative understanding of various influential factors affecting the deformation and failure mechanisms in layered rocks and hence provide guidance for rock engineering projects.

Samples preparation and experimental methods

Preparation of two groups briquettes contained bedding

Due to the disturbances caused by tectonic movements and underground excavations, coring samples in the field is extremely challenging. Furthermore, the majority of the rocks surrounding roadway in coal mine are coal, which has a lower overall strength, making sampling more challenging. Previous studies have suggested the briquettes might be an alternative option (Jasinge et al. 2011; Meng et al. 2021; Zhao et al. 2021). Therefore, this study focuses on investigating the deformation and failure mechanism of briquettes subject to the coupling effects of stress, lithology, and bedding. Here, stress represents confining pressure, lithology represents the briquettes with different proportions, and bedding represents the bedding angle. Cement with a grade of C42.5 is used as a binder. Coal powder with a diameter of 0.075 mm, cement, and purified water are mixed in proportions of 10:1:1.2 and 10:2:2.4, respectively, to create briquettes representing different lithologies. These two briquette ratios are named Group A samples and Group B samples, respectively. As shown in Fig. 1, after thoroughly mixing the coal powder, cement, and purified water, the mixture is compacted in three stages using pressure plates set at angles of 0°, 10°, 20°, 30°, and 45° and a compaction machine. This process is conducted to create briquettes specimens with bedding angles of 0°, 10°, 20°, 30°, and 45°, each having a height of 100 mm and a diameter of 50 mm (Fig. 2). The noteworthy point is that all specimens have a fixed bedding spacing of 20 mm. After the pouring is completed, the specimens are stabilized using the compaction machine for 20 min before demolding. The briquettes are then cured for 28 days to reach the desired strength.

Fig. 1
figure 1

Specimens preparation and experimental procedure flowchart

Full size image
Fig. 2
figure 2

Preparation schematic of briquettes with different bedding angles

Full size image

Taking specimen with 0° and 30° bedding angles, the sample preparation process is illustrated in Fig. 2. To cast the 0° bedding sample, the first layer is initially poured into the mold with height 130 mm followed by compaction by the plate at 0°. Then, layers 2 and 3 were casted in sequence before left for curing for 20 min. The specimen is then ready for test after demolding upon the completion of curing. As shown in Fig. 2, the specimens in this study consist of three layers. To ensure that the compaction stress on each specimen is the same, after compacting each layer, the exposed length of the pressure plates is measured using a ruler. It is ensured that the exposed length of the pressure plates after compacting the third layer of each specimen is consistent, and then this pressure is maintained for 20 min. In this case, the height of the mold is the same, and the compression amount of the pressure plate is also the same, thus ensuring that the molding pressure applied by the compactor on each specimen is the same. The similar process was applied to all other samples with various bedding angles. The only difference in the sample preparation was the oriented angle of the compaction plate set up according to the bedding angle for each specimen as demonstrated in Fig. 2.

Experimental methods and microscopic analysis of samples

As shown in the Fig. 1, the X-ray diffraction (XRD) and Scanning Electron Microscope (SEM) analyses are employed to examine the composition and internal microstructure of the two groups different briquettes. Finally, the Geotechnical Consulting & Testing Systems (GCTS RTR-1000) was employed to carry out the triaxial tests for two groups of samples under various confining pressures including 0 MPa, 0.5 MPa, 1.0 MPa, 1.5 MPa, and 2.0 MPa. The experimental program is listed in Table 1. The strain rate was set constantly at 0.05%/min. In Table 1, A0-1 represents the confining pressure of 0 MPa for the specimen with a bedding angle of 0° in Group A. Similarly, A10-5 represents a confining pressure of 2.0 MPa for the specimen with a bedding angle of 10° in Group A, and so forth.

Table 1 Physical characteristics of specimens and experimental methods
Full size table

The mineral composition and specific content of Group A samples and B samples were obtained through XRD. Figure 3 indicates both Group A samples and Group B samples contain quartz, kaolinite, calcite, gypsum, and analcime, but in varying proportions. The quartz content in Group A samples is 24%, kaolinite content is 31.7%, calcite content is 35.7%, gypsum content is 4.2%, and analcime content is 4.4%. On the other hand, Group B samples has 23.8% quartz, 25.3% kaolinite, 37.3% calcite, 8.5% gypsum, and 5.1% analcime. Among these, kaolinite is clay minerals, and the clay mineral content in Group A samples is 31.7%, while in Group B samples, it is 25.3%. Quartz and calcite can increase the strength of rocks, and gypsum can increase the strength and brittleness of rocks. Although the quartz and calcite content in Group A samples and Group B samples is similar, the higher cement content in Group B samples leads to a higher gypsum content compared to Group A samples. Based on this analysis, due to the same composition but different ratios, Groups A and B samples ultimately exhibit different physical and mechanical properties.

Fig. 3
figure 3

Quantitative XRD analysis of specimens

Full size image

SEM were used to observe the internal microstructures of Group A samples and B samples. Due to the consistent uniformity of the specimens, the molded Group A and Group B specimens were, respectively, crushed into blocks with a volume smaller than 1 cm3. Subsequently, SEM tests were conducted on them to observe the differences in the internal microstructure between the two groups of specimens. As depicted in Fig. 4, observations were taken at magnifications of 500×, 1000×, and 5000× at both specimens. It is evident that Group A samples possesses larger grain sizes and more pores than Group B samples indicating that Group B samples are denser and more cohesive than Group A samples. This can be attributed to the significant higher content of the gypsum in Group B samples (8.5%) being over twice of that in Group A samples (4.2%). Gypsum can interact with other minerals to fill the pores within the rock, consequently reduce pore size and connectivity and enhance the physical compactness of the rock. This further corroborates that Groups A and B samples have different physical and mechanical properties. These findings provide the basis for investigating the deformation and failure mechanisms of coal under the coupling effects of stress, lithology, and bedding.

Fig. 4
figure 4

SEM results of different specimens with different magnification factors

Full size image

Experimental results and analysis

Mechanical properties and failure modes

Figures 5 and 6 depict the stress–strain curves of Group A samples and Group B samples under different confining pressures. It can be observed that both Group A samples and Group B samples undergo elastic, yielding, and failure stages in sequence regardless of confining pressures, whereas the compaction stage is not significant. This arises because both Group A samples and Group B samples were pre-compacted using a compaction machine during their preparation, leading to an overall lower porosity. Furthermore, it is evident from Figs. 5 and 6 that both Group A samples and Group B samples exhibit brittle failure and the yielding stage became less phenomenal at lower confining pressures. In addition, the brittle failure of Group B samples is more pronounced than those of Group A samples at each confining pressure level. This could be attributed to the higher gypsum content in Group B samples (8.5%) being over twice of that in Group A samples (4.2%).

Fig. 5
figure 5

Deviator stress–strain curves of Group A samples under different confining pressures

Full size image
Fig. 6
figure 6

Deviator stress–strain curves of Group B samples under different confining pressures

Full size image

Under the coupled effects of stress, lithology, and bedding, specimens exhibit different failure characteristics. As shown in Fig. 7, from left to right are typical failure photographs of Group A samples and Group B samples with bedding angles of 0°, 10°, 20°, 30°, and 45°, respectively. Under uniaxial compression, the specimens of both Group A samples and Group B samples show predominantly axial extension fractures with little shear fracture, indicating tensile failure across the intact rock components and bedding planes. It's worth noting that both Group A samples and Group B samples demonstrate composite tension-shear failure at a bedding angle of 45°, involving fractures crossing both bedding and parallel to bedding planes, under uniaxial compression. Upon applying confining pressure, both specimens of Group A samples and Group B samples with the 0°, 10°, 20°, and 30° initially show shear-tension composite failure, with increasing confining pressure emphasizing the role of shear failure. With further increases in confining pressure, the expansion of tension cracks within the specimen is suppressed, transitioning to shear failure. And the 45° samples exhibits shear slip-tension composite failure. Specifically, at 0.5 MPa confining pressure, tension cracks are still observed in the 0°, 10°, 20°, and 30° bedding specimens of both Group A samples and Group B samples, primarily accompanied by shear cracks. At a bedding angle of 45°, tension cracks connecting two weak bedding planes become visible, signifying shear slip-tension composite failure dominated by shear sliding along the bedding planes. When the confining pressure reaches 1.0 MPa, tension cracks are almost absent in both Group A samples and Group B samples. In the 45° bedding specimens, tension cracks appear between the two bedding planes, and showing distinct sliding marks. This indicates shear slip-tension composite failure predominantly driven by bedding plane sliding. Like the 1.0 MPa confining pressure, at 1.5 MPa and 2.0 MPa confining pressures, the 0°, 10°, 20°, and 30° bedding specimens of both Group A samples and Group B samples primarily display shear failure with secondary shear cracks along the main shear crack path. The 45° bedding specimens continue to exhibit shear slip-tension composite failure dominated by shear sliding along the bedding planes. In addition, it can also be observed from the Fig. 7 that the direction of crack propagation exhibits two forms when encountering bedding planes. One is where the crack maintains its original direction and directly crosses the bedding plane. The other is where the crack changes direction as it crosses the bedding plane or continues in the original direction or changes direction after extending along the bedding plane. This phenomenon becomes increasingly apparent with the increase in bedding plane angles, and is more pronounced under low confining pressure.

Fig. 7
figure 7

Failure modes of specimens with varying bedding angles under different confining pressures (the white dashed line represents the bedding planes position, and the yellow dashed line represents cracks)

Full size image

In general, the 0°, 10°, 20°, and 30° bedding specimens primarily exhibit tension-shear composite failure under uniaxial compression. However, under triaxial compression, they display shear-tension composite failure dominated by shear, with increasing confining pressure progressively limiting the sliding of weak bedding planes. Consequently, after the confining pressure increases to 1.0 MPa, the failure mode shifts to shear failure, and the shear failure angle decreases as the confining pressure increases. On the other hand, the 45° bedding specimens exhibit a composite tension-shear failure involving both parallel to and crossing the bedding planes under uniaxial compression. Under triaxial compression, tension failure appears between the two bedding planes, characterized by evident slip displacement marks on the bedding planes. The overall behavior is dominated by shear slip-tension composite failure driven by sliding along the bedding planes. Furthermore, due to the higher gypsum content in Group B samples (8.5%) compared to Group A samples (4.2%), Group B samples displays more pronounced brittle failure characteristics at higher confining pressures (Fig. 7d and e) compared to Group A samples.

Summarizing the crack types observed in the experimental failure diagrams in Fig. 7, three crack types are illustrated in Fig. 8. Specifically, the first type includes tensile cracks and tensile wing cracks, occurring in specimens subjected to uniaxial compression at 0°, 10°, 20°, and 30°. The second type consists of shear cracks, predominantly appearing in specimens subjected to relatively high confining pressures (1.0 MPa, 1.5 MPa, and 2.0 MPa) at 0°, 10°, 20°, and 30°, with shear cracks also observed in the 30° specimen under 0.5 MPa confining pressure. The third type is a combination of tensile and shear cracks, exhibiting varied forms in specimens under different confining pressures and bedding angles. For instance, the 45° specimen displays tensile cracks in the upper and lower parts during uniaxial compression, with shear cracks along the bedding planes in the central part. In triaxial compression, the upper and lower parts exhibit shear cracks along the bedding planes, while tensile cracks occur between the bedding planes. The 0°, 10°, and 20° specimens show tensile cracks in the upper part and shear cracks in the central and lower parts under 0.5 MPa confining pressure. Additionally, it is noteworthy that when cracks encounter bedding planes, they manifest in two ways: either directly penetrating through and continuing to propagate along the original crack direction or altering the propagation direction after crossing the bedding planes, followed by continued propagation along the altered crack direction. As indicated by the red circles in Fig. 8, the latter mode of crossing bedding planes becomes more pronounced with increasing bedding angles. From Figs. 7 and 8, it can be observed that confining pressure and bedding angle have varying degrees of influence on the crack types and failure modes of layered rocks. Next, the specific effects of confining pressure and bedding angle on the mechanical parameters of rocks will be analyzed.

Fig. 8
figure 8

Crack types of specimens with varying bedding angles under different confining pressures

Full size image

Mechanical parameter analysis

Figure 9 illustrates the influence of bedding angle on the peak deviatoric stress of Group A samples and Group B samples. From Fig. 9, it can be observed that the peak deviatoric stress of both Group A samples and Group B samples decreases as the bedding angle increases. When the bedding angle is 0°, there is no induced shear stress along the weak bedding plane. As the weak bedding plane angle increases from 0° to 45°, the induced shear stress component along the weak bedding plane increases, consequently leading to a reduction in the overall shear strength of the rock. In addition, the peak deviatoric stresses of both Group A samples and Group B samples increases with the confining pressure. The largest increase occurs when the confining pressure changes from 0 to 0.5 MPa, and the rate of increase slows down when it goes from 0.5 to 2.0 MPa. This phenomenon can be explained by two main factors. On one hand, it is known that when confining pressure increases, a greater stress is required to trigger internal fracture and sliding within the rock according to the relationship between the maximum principal stress and confining pressure in the Mohr–Coulomb criterion. On the other hand, the increase in confining pressure leads to an increase in the contact area between the grains or particles inside the rock, enhancing the frictional forces among the particles, thereby increasing the compressive strength of the rock. Consequently, as confining pressure continues to increase, the rate of increase in the overall compressive strength of the rock slows down.

Fig. 9
figure 9

Relationship between deviatoric stress and bedding angle in specimen

Full size image

Figure 10 illustrates the impact of bedding angle and confining pressure on the elastic modulus of Group A samples and Group B samples. From Fig. 10, it can be observed that the elastic modulus of both Group A samples and Group B samples increases with the bedding angle. When the bedding angle is 0°, the stress loading direction is perpendicular to the bedding planes, which is essentially the same as that of the specimen without bedding. The difference lies in the fact that the strength of the bedding planes has lower strength compared to the rock components, and the overall strength of the rock is controlled by the weaker bedding planes. This leads to greater axial deformation. With an increasing angle of bedding planes, the stress acting on the specimen is influenced by the presence of weak bedding planes, resulting in a decrease in the axial stress component. Furthermore, when the same axial stress is applied, a greater bedding plane angle leads to a smaller axial stress component. This ultimately results in a lower stress driving axial deformation as the bedding plane angle increases, manifesting as a higher elastic modulus. The experimental results observed in this study, showing a gradual increase in elastic modulus with an increasing angle of bedding planes, are consistent with Shen et al. (2021). Additionally, as the confining pressure increases, the elastic modulus also increases gradually, following a similar trend to the increase in peak deviatoric stress. Similarly, when the confining pressure changes from 0 to 0.5 MPa, there is a larger change in the elastic modulus, followed by a gradual decrease in the rate of increase with further increments in confining pressure. This is attributed to the effect of increasing confining pressure on gradually compacting the internal pore structure of the specimen, enhancing its ability to resist deformation, and thus leading to an increase in elastic modulus.

Fig. 10
figure 10

The relationship between elastic modulus and bedding angle in specimen

Full size image

Figure 11 shows the variations in cohesion and internal friction angle for Group A samples and Group B samples at different bedding angles. It can be observed that both cohesion and internal friction angle gradually decrease with increasing bedding angle. This is because as the bedding angle increases, the contact area between particles on bedding planes gradually increases, leading to an enhanced sliding action between bedding planes. This makes the relative motion of particles inside the sample easier and requires them to overcome additional sliding resistance to resist shear failure, which results in a decrease in both cohesion and internal friction angle.

Fig. 11
figure 11

The cohesion and internal friction angle of samples at different bedding angles

Full size image

Based on the comprehensive analysis above, it can be concluded that confining pressure, lithology, and bedding angle are three factors that exert varying degrees of influence on the deformation and failure characteristics as well as the physical–mechanical parameters of specimens. Subsequently, the specific impact weights of stress, lithology, and bedding on the deformation and failure of rock will be determined using the orthogonal matrix analysis method.

Constitutive model

Orthogonal matrix analysis of influential factors

In this study, stress, lithology, and bedding are selected as input parameters, while peak differential stress, elastic modulus, and peak strain of rock are chosen as responses. Following the rules of orthogonal experiment table design, as shown in Table 2, the orthogonal experiment table can be developed. In Table 2, the meanings of k1, k2, k3, k4, and k5 correspond to the average values of the experimental indicators at a particular level of each factor. For instance, when visually analyzing the data for differential stress in Table 2 , k1 under the stress factor represents the average differential stress at the first level of the stress factor, while k5 represents the average differential stress at the fifth level of the stress factor. Then, use the orthogonal matrix analysis method to conduct orthogonal matrix analysis on the mechanical properties of rocks under the influence of multiple factors and multiple indicators. The orthogonal matrix analysis method, as well as specific derivation process of the weight matrices of the stress, lithology, and bedding, can be found in Appendix A.

Table 2 Results of orthogonal test
Full size table

The specific impact weights of stress, lithology, and bedding on rock can be summarized by Eq. (1). Specifically, the influence weights of stress, lithology, and bedding on rock mechanics properties are 40.99%, 23.87%, and 35.14%, respectively. Such impact weighting matrices of stress, lithology, and bedding were used to develop a constitutive model for characterizing the failure process of rocks subject to triaxial loading condition in next section.

$$Q = k^{\rm T} = \left( {\begin{array}{*{20}c} {Q_{{\text{S}}} } & {Q_{{\text{L}}} } & {Q_{{\text{B}}} } \\ \end{array} } \right) = \left( {0.4099\qquad 0.2387\qquad 0.3514} \right)$$
(1)

where, Q represents the weight matrix that characterizes the specific impact weights of stress, lithology, and bedding on rock deformation and failure. QS, QL, and QB, respectively stand for the impact weights of stress, lithology, and bedding on rock deformation and failure.

Model establishment

Assuming that the internal strength of individual microelements within the rock follows a Weibull distribution, according to statistics, its probability density function is as follows (Shen et al. 2021).

$$P\left( F \right) = \frac{m}{{F_{0} }}\left( {\frac{F}{{F_{0} }}} \right)^{m - 1} \exp \left[ { - \left( {\frac{F}{{F_{0} }}} \right)^{m} } \right]$$
(2)

F represents the random distribution variable of microelement strength within the rock, while m and F0 are parameters of the Weibull distribution. The cumulative damage of microelements within the rock will lead to macroscopic damage and failure. Assuming that at a certain moment and under a certain stress state, the number of microelements that fail within the rock is Nd out of a total of N microelements, the damage variable D can be defined as follows.

$$D = \frac{{N_{{\text{d}}} }}{N}$$
(3)

Based on this, when F ≥ 0, the random distribution variable is within any interval [F, F + dF], and the number of failed rock microelements is NP(x)dx. When F < 0, no damage occurs in the rock and hence the damage variable D is 0. Therefore, when the stress level reaches F, the number of microelements that fail at that point can be represented by the following Eq. (4).

$$N_{{\text{d}}} \left( F \right) = \int_{0}^{F} {NP\left( x \right)} {\text{d}}x = N\int_{0}^{F} {P\left( x \right)} {\text{d}}x = N\left\{ {1 - \exp \left[ { - \left( {\frac{F}{{F_{0} }}} \right)^{m} } \right]} \right\}$$
(4)

When F ≥ 0, the expression for the damage variable D is as follows.

$$D = \frac{{N_{{\text{d}}} }}{N} = \frac{{N\left\{ {1 - \exp \left[ { - \left( {\frac{F}{{F_{0} }}} \right)^{m} } \right]} \right\}}}{N} = 1 - \exp \left[ { - \left( {\frac{F}{{F_{0} }}} \right)^{m} } \right]$$
(5)

Therefore, the complete expression for the damage variable D is as follows.

$$D = \left\{ {\begin{array}{*{20}l} {1 - \exp \left[ { - \left( {\frac{F}{{F_{0} }}} \right)^{m} } \right]} \hfill & {F \ge 0} \hfill \\ 0 \hfill & {F < 0} \hfill \\ \end{array} } \right.$$
(6)

Based on the assumption of equivalent strain and the deformation compatibility between the damaged and undamaged parts within the rock, the three-dimensional statistical damage constitutive equation for rock can be derived according to (Cao et al. 2019; Lemaitre 1984).

$$\sigma_{1} = E\varepsilon_{1} \left( {1 - D} \right) + \mu \left( {\sigma_{2} + \sigma_{3} } \right)$$
(7)

By substituting Eqs. (6) into (7), the three-dimensional statistical damage constitutive equation rearranged accordingingly as:

$$\sigma_{1} = \left\{ {\begin{array}{*{20}l} {E\varepsilon_{1} \exp \left[ { - \left( {\frac{F}{{F_{0} }}} \right)^{m} } \right] + \mu \left( {\sigma_{2} + \sigma_{3} } \right)} \hfill & {F \ge 0} \hfill \\ {E\varepsilon_{1} + \mu \left( {\sigma_{2} + \sigma_{3} } \right)} \hfill & {F < 0} \hfill \\ \end{array} } \right.$$
(8)

where, E represents the elastic modulus of the rock, μ is the Poisson's ratio of the rock, ε1 is the axial strain of the rock, and σ3 is the confining pressure during the triaxial testing on rocks, in this case, σ2 = σ3. Therefore, to obtain the constitutive equation for damage softening of rock under three-dimensional conditions, it is necessary to determine the random distribution variable F, as well as the parameters m and F0 of the Weibull distribution.

In the past, researchers typically relied on rock failure criteria to determine the random distribution variable F. The rock failure criteria can be generally expressed according to Eq. (9) (Deng and Gu 2011; Zhao et al. 2017).

$$f\left( \sigma \right) - K_{0} = 0$$
(9)

where, K0 is a constant related to the experiment, often referred to the damage threshold. When the stress level on the rock exceeds the damage threshold, microelements within the rock may yield or fail. Hence, previous researchers often treated f (σ)−K0 as the random distribution variable F (Deng and Gu 2011; Zhao et al. 2017). However, in layered rocks, they are influenced by stress, lithology, and bedding. Considering the impact of the damage threshold, an orthogonal matrix analysis was used to derive a stress-lithology-bedding weight matrix, which was applied to modify this variable, resulting in a new random distribution variable F.

$$F = Q\left[ {f\left( \sigma \right) - K_{0} } \right] = \left( {0.4099\sigma_{3} + 0.2387\sigma_{c} + 0.3514\beta_{b} } \right)\left[ {f\left( \sigma \right) - K_{0} } \right]$$
(10)

In the Eq. (10), Q represents the stress-lithology-bedding weight matrix, where σ3 stands for the confining pressure applied onto the rock, σc represents the uniaxial compressive strength of the intact rock, and βb denotes the bedding angle of the rock.

The expression for the random distribution variable F based on the commonly used rock failure criterion, the Drucker–Prager criterion, was derived as shown in Eq. (11) (Deng and Gu 2011; Zhao et al. 2017).

$$f\left( \sigma \right) - K_{0} = \alpha_{0} I_{1} + \sqrt {J_{2} } - K_{0}$$
(11)

In Eq. (11), I1 represents the first invariant of the stress tensor, J2 stands for the second invariant of the stress deviator, α0 and K0 are parameters related to the properties of the rock. The specific expressions for these four parameters are as follows (Zhao et al. 2017).

$$I_{1} = \sigma_{x}^{ * } + \sigma_{y}^{ * } + \sigma_{z}^{ * } = \sigma_{1}^{ * } + \sigma_{2}^{ * } + \sigma_{3}^{ * }$$
(12)
$$\begin{gathered} J_{2} = \frac{1}{6}\left[ {\left( {\sigma_{1}^{ * } - \sigma_{2}^{ * } } \right)^{2} + \left( {\sigma_{2}^{ * } - \sigma_{3}^{ * } } \right)^{2} + \left( {\sigma_{1}^{ * } - \sigma_{3}^{ * } } \right)^{2} } \right] \\ = \frac{1}{6}\left[ {\left( {\sigma_{x}^{ * } - \sigma_{y}^{ * } } \right)^{2} + \left( {\sigma_{x}^{ * } - \sigma_{z}^{ * } } \right)^{2} + \left( {\sigma_{y}^{ * } - \sigma_{z}^{ * } } \right)^{2} } \right] + \left( {\tau_{xy}^{ * } } \right)^{2} + \left( {\tau_{yz}^{ * } } \right)^{2} + \left( {\tau_{xz}^{ * } } \right)^{2} \\ \end{gathered}$$
(13)
$$\alpha_{0} = \frac{\sin \varphi }{{\sqrt {9 + 3\sin^{2} \varphi } }}$$
(14)
$$K_{0} = \frac{\sqrt 3 c\cos \varphi }{{\sqrt {3 + \sin^{2} \varphi } }}$$
(15)

In Eqs. (12) and (13), the effective stresses σ* 1, σ* 2, and σ* 3 can be obtained from the triaxial test stresses σ1, σ2, and σ3. φ represents the internal friction angle of the rock. According to Hooke's law, the following equations can be obtained.

$$\varepsilon_{1} = \frac{1}{E}\left[ {\sigma_{1}^{ * } - \mu \left( {\sigma_{2}^{ * } + \sigma_{3}^{ * } } \right)} \right]$$
(16)
$$\sigma_{3}^{ * } = \sigma_{2}^{ * } = \frac{{\sigma_{3} }}{1 - D}$$
(17)
$$\sigma_{1}^{ * } = \frac{{\sigma_{1} }}{1 - D}$$
(18)
$$I_{1} = \sigma_{1}^{ * } + 2\sigma_{3}^{ * } = \frac{{\sigma_{1} + 2\sigma_{3} }}{1 - D} = \frac{{\left( {\sigma_{1} + 2\sigma_{3} } \right)\left( {\sigma_{1} - 2\mu \sigma_{3} } \right)}}{{\left( {\sigma_{1} - 2\mu \sigma_{3} } \right)\left( {1 - D} \right)}} = \frac{{\left( {\sigma_{1} + 2\sigma_{3} } \right)E\varepsilon_{1} }}{{\sigma_{1} - 2\mu \sigma_{3} }}$$
(19)
$$\sqrt {J_{2} } = \frac{{\sigma_{1}^{ * } - \sigma_{3}^{ * } }}{\sqrt 3 } = \frac{{\sigma_{1} - \sigma_{3} }}{{\sqrt 3 \left( {1 - D} \right)}} = \frac{{\left( {\sigma_{1} - \sigma_{3} } \right)\left( {\sigma_{1} - 2\mu \sigma_{3} } \right)}}{{\sqrt 3 \left( {\sigma_{1} - 2\mu \sigma_{3} } \right)\left( {1 - D} \right)}} = \frac{{\left( {\sigma_{1} - \sigma_{3} } \right)E\varepsilon_{1} }}{{\sqrt 3 \left( {\sigma_{1} - 2\mu \sigma_{3} } \right)}}$$
(20)

By substituting Eqs. (19) and (20) into Eq. (11), and then inserting Eq. (11) into Eq. (10), the expression for the random distribution variable F under the Drucker–Prager criterion can be obtained.

$$\begin{gathered} F = Q\left( {\alpha_{0} I_{1} + \sqrt {J_{2} } - K_{0} } \right) = \left( {0.4099\sigma_{3} + 0.2387\sigma_{c} + 0.3514\beta_{b} } \right) \\ \left[ {\frac{{\sin \varphi \left( {\sigma_{1} + 2\sigma_{3} } \right)E\varepsilon_{1} }}{{\left( {\sigma_{1} - 2\mu \sigma_{3} } \right)\sqrt {9 + 3\sin^{2} \varphi } }} + \frac{{\left( {\sigma_{1} - \sigma_{3} } \right)E\varepsilon_{1} }}{{\sqrt 3 \left( {\sigma_{1} - 2\mu \sigma_{3} } \right)}} - \frac{\sqrt 3 c\cos \varphi }{{\sqrt {3 + \sin^{2} \varphi } }}} \right] \\ \end{gathered}$$
(21)

Parameters determination

So far, the constitutive model incorporating the rock stress, lithology, and bedding was developed, and now the Weibull distribution parameters, m and F0, should be quantified. Since damage variable D is 0 when F is less than 0, it is necessary to analyze the three-dimensional damage-softening constitutive equation for cases where F is greater than or equal to 0. When the rock exhibits strain-softening characteristics, the stress–strain curve at the peak point (σ1 = σsc, ε1 = εsc) exhibits extremum properties.

$$\frac{{\partial \sigma_{1} }}{{\partial \varepsilon_{1} }}\left| {_{{\sigma_{1} = \sigma_{sc} ,\varepsilon_{1} = \varepsilon_{sc} }} } \right. = 0$$
(22)

Taking the partial derivative of Eq. (7) and substituting the boundary conditions from Eq. (22), the following equation can be obtained.

$$\frac{{\partial \sigma_{1} }}{{\partial \varepsilon_{1} }}\left| {_{{\sigma_{1} = \sigma_{sc} ,\varepsilon_{1} = \varepsilon_{sc} }} } \right. = E\left( {1 - D_{sc} } \right) - E\varepsilon_{sc} \frac{\partial D}{{\partial \varepsilon_{1} }} = 0$$
(23)

Rearranging the Eq. (23) leads to the following equation:

$$\frac{\partial D}{{\partial \varepsilon_{1} }} = \frac{{1 - D_{sc} }}{{\varepsilon_{sc} }}$$
(24)

where Dsc represents the damage value at the peak point.

Similarly, taking the partial derivative of the first equation in Eq. (6), and letting σ1 = σsc, ε1 = εsc, and F = Fsc, the following equation can be obtained.

$$\frac{\partial D}{{\partial \varepsilon_{1} }}\left| {_{{\sigma_{1} = \sigma_{sc} ,\varepsilon_{1} = \varepsilon_{sc} ,F = F_{sc} }} } \right. = - \exp \left[ { - \left( {\frac{{F_{sc} }}{{F_{0} }}} \right)^{m} } \right]\left[ { - \frac{m}{{F_{0} }}\left( {\frac{{F_{sc} }}{{F_{0} }}} \right)^{m - 1} } \right]\frac{\partial F}{{\partial \varepsilon_{1} }}$$
(25)

where Fsc represents the random distributed variable corresponding to the peak point. Furthermore, the Eq. (26) can be obtained by rearranging the terms in the Eq. (6).

$$\exp \left[ { - \left( {\frac{{F_{{{\text{sc}}}} }}{{F_{0} }}} \right)^{m} } \right] = 1 - D_{{{\text{sc}}}}$$
(26)

Rearranging the Eq. (26) gives the Eq. (27).

$$- \left( {\frac{{F_{{{\text{sc}}}} }}{{F_{0} }}} \right)^{m - 1} = \frac{{F_{0} }}{{F_{{{\text{sc}}}} }}\ln \left( {1 - D_{{{\text{sc}}}} } \right)$$
(27)

Simultaneously, taking the partial derivative of Eq. (21), Eq. (28) can be obtained.

$$\frac{\partial F}{{\partial \varepsilon_{1} }}\left| {_{{\sigma_{1} = \sigma_{{{\text{sc}}}} }} } \right. = Q\left[ {\frac{{\alpha_{0} E\left( {\sigma_{{{\text{sc}}}} + 2\sigma_{3} } \right) + \left( {\sigma_{{{\text{sc}}}} - \sigma_{3} } \right)E/\sqrt 3 }}{{\sigma_{{{\text{sc}}}} - 2\mu \sigma_{3} }}} \right]$$
(28)

Substituting Eqs. (26)–(28) into Eq. (25), the following equation can be obtained.

$$\frac{\partial D}{{\partial \varepsilon_{1} }} = - \left( {1 - D_{{{\text{sc}}}} } \right)\frac{m}{{F_{{{\text{sc}}}} }}\ln \left( {1 - D_{{{\text{sc}}}} } \right)Q\left[ {\frac{{\alpha_{0} E\left( {\sigma_{{{\text{sc}}}} + 2\sigma_{3} } \right) + \left( {\sigma_{{{\text{sc}}}} - \sigma_{3} } \right)E/\sqrt 3 }}{{\sigma_{{{\text{sc}}}} - 2\mu \sigma_{3} }}} \right]$$
(29)

Combining Eq. (24) and Eq. (29).

$$\frac{{1 - D_{sc} }}{{\varepsilon_{sc} }} = - \left( {1 - D_{{{\text{sc}}}} } \right)\frac{m}{{F_{{{\text{sc}}}} }}\ln \left( {1 - D_{{{\text{sc}}}} } \right)Q\left[ {\frac{{\alpha_{0} E\left( {\sigma_{{{\text{sc}}}} + 2\sigma_{3} } \right) + \left( {\sigma_{{{\text{sc}}}} - \sigma_{3} } \right)E/\sqrt 3 }}{{\sigma_{{{\text{sc}}}} - 2\mu \sigma_{3} }}} \right]$$
(30)

The expressions for m and F0 can be determined from Eqs. (30) and (26) as follows.

$$m = - \frac{{F_{{{\text{sc}}}} }}{{\ln \left( {1 - D_{{{\text{sc}}}} } \right)QB}}$$
(31)
$$F_{0} = \frac{{F_{{{\text{sc}}}} }}{{\left[ { - \ln \left( {1 - D_{{{\text{sc}}}} } \right)} \right]^{1/m} }}$$
(32)

where,

$$F_{sc} = Q\left[ {\frac{{\sin \varphi \left( {\sigma_{{{\text{sc}}}} + 2\sigma_{3} } \right)E\varepsilon_{{{\text{sc}}}} }}{{\left( {\sigma_{{{\text{sc}}}} - 2\mu \sigma_{3} } \right)\sqrt {9 + 3\sin^{2} \varphi } }} + \frac{{\left( {\sigma_{{{\text{sc}}}} - \sigma_{3} } \right)E\varepsilon_{{{\text{sc}}}} }}{{\sqrt 3 \left( {\sigma_{{{\text{sc}}}} - 2\mu \sigma_{3} } \right)}} - \frac{\sqrt 3 c\cos \varphi }{{\sqrt {3 + \sin^{2} \varphi } }}} \right]$$
(33)
$$D_{{{\text{sc}}}} = 1 - \frac{{\sigma_{{{\text{sc}}}} - 2\mu \sigma_{3} }}{{E\varepsilon_{{{\text{sc}}}} }}$$
(34)
$$B = \frac{{\sin \varphi \left( {\sigma_{{{\text{sc}}}} + 2\sigma_{3} } \right)E\varepsilon_{{{\text{sc}}}} }}{{\left( {\sigma_{{{\text{sc}}}} - 2\mu \sigma_{3} } \right)\sqrt {9 + 3\sin^{2} \varphi } }} + \frac{{\left( {\sigma_{{{\text{sc}}}} - \sigma_{3} } \right)E\varepsilon_{{{\text{sc}}}} }}{{\sqrt 3 \left( {\sigma_{{{\text{sc}}}} - 2\mu \sigma_{3} } \right)}}$$
(35)

Model validation

The constitutive modeling results were validated by the experimental results of Groups A and B samples under different confining pressures in Sect. "Experimental results and analysis". And the parameters used in the constitutive modeling are presented in Table 3. It is worth noting that a certain amount of coal powder, cement, and purified water mixture can be poured into the mold depicted in Fig. 1 all at once. Subsequently, by following the sample preparation method outlined in Fig. 1, briquette samples without bedding planes can be obtained. The uniaxial compressive strength σc for the intact rock can then be determined using the GCTS shown in Fig. 1. And the other parameters in Table 3 can be obtained through the experimental data presented in the Sect. "Experimental results and analysis".

Table 3 All parameters used in constitutive modeling of samples
Full size table

As shown in Figs. 12 and 13, the modeling results agree well with the experimental results confirming the capability of the proposed model for simulating the behavior of the layered rock with various bedding angles under triaxial loading condition. It is worthy to highlight that this model is not only able to capture the pre-peak deformation characteristics of rocks (Groups A and B) at different confining pressures and bedding angles, also able to effectively reflect the sensitivity of peak strength of the rock subject to any changes in stress, lithology, and bedding. Furthermore, this model is able to predict the post-peak strain-softening behavior of rocks under low confining pressures. However, it struggles to accurately capture the strain-hardening behavior after the peak. This limitation is more apparent in the testing data of Group B samples (Fig. 13), where the higher gypsum content leads to a greater likelihood of strain hardening, resulting in the discrepancies between the modeling results and the experimental results. This discrepancy can be attributed to the assumption of equivalent strain used in the model, which fails to account for the strain-hardening characteristics of the rock.

Fig. 12
figure 12

Comparison of experimental and constitutive modeling results for Group A samples

Full size image
Fig. 13
figure 13

Comparison of experimental and constitutive modeling results for Group B samples

Full size image

The constitutive model proposed in this study was further validated with axial stress–strain data obtained from literature (Shen et al. 2021), where triaxial testing subject to 0.5 MPa and 1.5 MPa confining pressures were conducted on interbedded rock samples with bedding angles of 0°, 30°, 45°, and 60°. As shown in Fig. 14, the modeling results agree well with the experimental data reported by Shen et al. (2021). However, there is a slight deviation at post-peak stage. This is understandable as the core interest of this model is for the pre-peak stage whereas the post-peak is less attendant.

Fig. 14
figure 14

Comparison between the constitutive modeling results and experimental data from reference (Shen et al. 2021)

Full size image

Conclusions

Triaxial compressive tests were conducted on two groups of rocks at various bedding angles and confining pressures in this study. The orthogonal matrix analysis was applied to quantify the weightings of the influence of the stress, lithology, and bedding on the strength, elastic modulus and peak strain for the rock. Finally, a new statistical constitutive model was developed to simulate the stress–strain behavior of layered rocks subject to various stresses, lithologies, and bedding angles. The key findings are summarized below.

  1. (1)

    In the preparation of briquettes with different bedding angle in this paper, a method of pouring in layers sequentially and then compacting layer by layer was used, which replicates the lithification process of layered rock formations in actual engineering. The method for preparing briquettes described in this paper has certain guiding significance for the preparation of layered rock samples with different bedding angles.

  2. (2)

    For specimens with bedding angles of 0°, 10°, 20°, and 30°, tensional failure prevails during the uniaxial compressive test. At 45° bedding angle, both tensional and shear sliding failure along the bedding plane are observed. The direction of crack propagation deviates after crossing the bedding plane. Under triaxial compression, shear dominates the failure at 0°, 10°, 20°, and 30° angles whereas at 45°, the predominant mode of failure is a shear sliding-tension combination. When cracks encounter bedding planes, they can either follow the bedding plane or change direction to keep propagating. This phenomenon becomes increasingly noticeable with an increase in the bedding angle and gradually diminishes with higher confining pressure.

  3. (3)

    As a result of multi-variate analysis based on orthogonal matrix approach, the weightings of the influence of stress, lithology, and bedding on rock deformation and failure were determined to be 40.99%, 23.87% and 35.14%, respectively. The new statistical constitutive model proposed in this study is evidenced to successfully simulate the stress–strain behavior of the interbedded layered rocks under both uniaxial and triaxial compressive loading conditions in spite of a minor limitation in capturing the post-peak behavior for the rocks.