Introduction

Coal mining can cause damage to the roof structure, releasing a large amount of energy stored in strata, and when this energy is transferred to the roadway area it may cause damage to the roadway, or even rockburst1,2,3,4. Relevant researches have shown that the bending strain energy stored in the roof structure is an important factor in inducing rockburst5,6,7. When there are multiple layers of roofs with different conditions above the coal seam, they collectively form a layered composite roof structure. However, there is no reasonable theory to explain the bending deformation and energy storage characteristics of such roof structures. Therefore, it is of great theoretical significance to study the bending deformation and energy storage mechanism of LCRS.

To address the deformation and energy storage laws of Layered Composite Roof Structure, scholars have conducted extensive research via theoretical analysis, laboratory experiments, numerical simulations, and engineering measurements. In theoretical research, many scholars have proposed various theoretical models considering the characteristics of two or even multiple layers of roof breaking based on classical roof breaking models8,9,10,11. These include the double key strata structure12, composite pressure-arch model13, elastic foundation composite beam14,15, and interacting hard roof structures16. To verify and analyze the deformation process before roof breaking, scholars have conducted numerous laboratory experiments on layered roofs, deeply studying the interaction relationship between multiple layers of roofs from perspectives such as deformation field17,18, bearing capacity19,20, fracture characteristics21,22,23, and energy evolution laws24,25,26. In order to more realistically reflect the complexity of engineering geological conditions, scholars have utilized numerical simulation methods to analyze the effects of variables such as roof thickness27,28, rock beam span29,30,31, mechanical parameters32,33, and dynamic stress34 on the deformation and fracture characteristics of Layered Composite Roof Structure. Furthermore, scholars have also conducted on-site detection studies using a variety of approaches, including borehole detection35,36, microseismic monitoring37,38, and roof measurements39,40, among others.

As mentioned above, scholars have studied the deformation and energy characteristics of LCRS from different perspectives, and have achieved valuable research results. However, existing studies have primarily concentrated on qualitative research and have not yet established a mechanical model to describe the bending deformation issues of LCRS, leading to the inability to calculate the energy storage under such conditions. Therefore, this study established a mechanical model for the bending energy storage of LCRS under three conditions, then calculated the expression of bending energy under the above cases. Finally, the model was verified by numerical simulation and laboratory experiments. This research aims to provide a theoretical basis for the prevention and control of rockburst caused by roof fracture.

Basic models and assumptions

Force analysis of composite roof

The key layer is a load-bearing structure composed of a single layer or adjacent multiple hard rock layers, which plays a vital role in controlling the movement of the overlying rock. When the key layer consists of two layers of rock together, this bearing structure can be simplified to a composite beam that will only support the loads of the overlying strata before the key layer fracture, as shown in Fig. 1.

Fig. 1
figure 1

Schematic diagram of composite roof structure.

Full size image

Assuming the load exerted by the overlying rock layer is a concentrated one, this section of the composite beam undergoes transverse bending. Based on the plane section assumption, the normal stress distribution across the cross-section of a rectangular beam during bending deformation is as follows:

$$\sigma =\frac{{My}}{{{I_z}}}$$
(1)

Considering that the upper and lower layers of the structure consist of distinct lithologies, characterized by elastic moduli E1 and E2, and thicknesses h1 and h2, respectively, the formulations for shear force and bending moment within the beam are presented in Eqs. (2) and (3), as illustrated in Fig. 2.

Fig. 2
figure 2

Schematic diagram of force analysis of composite beam.

Full size image
$${F_S}\left( x \right)=\left\{ \begin{gathered} \begin{array}{*{20}{c}} {P/2}&{} \end{array}\begin{array}{*{20}{c}} {}&{(0 \leqslant x<S/2)} \end{array} \hfill \\ - P/2\begin{array}{*{20}{c}} {}&{}&{(S/2 \leqslant x<S)} \end{array} \hfill \\ \end{gathered} \right.$$
(2)
$$M\left( x \right)=\left\{ \begin{gathered} Px/2\begin{array}{*{20}{c}} {}&{} \end{array}\begin{array}{*{20}{c}} {}&{\begin{array}{*{20}{c}} {\begin{array}{*{20}{c}} {}&{(0 \leqslant x<S/2)} \end{array}}&{} \end{array}} \end{array} \hfill \\ PS/2 - Px/2\begin{array}{*{20}{c}} {}&{(S/2 \leqslant x \leqslant S)} \end{array} \hfill \\ \end{gathered} \right.$$
(3)

where S denotes the span of the rock beam.

Based on the plane assumption, it can be deduced that the longitudinal linear strain at any point on the cross-section of the composite beam varies linearly with respect to the height of the section. The strain is expressed as:

$$\varepsilon =\frac{y}{\rho }$$
(4)

Where, ρ is the radius of curvature of the neutral axis.

Neutral axis position criterion

When both materials are within their linear elastic ranges, the bending normal stresses in the upper and lower layers of the beam’s cross-section are calculated using Hooke’s Law, as follows:

$$\left\{ \begin{gathered} {\sigma _1}={E_1}\varepsilon =\frac{{{E_1}y}}{\rho } \hfill \\ {\sigma _2}={E_2}\varepsilon =\frac{{{E_2}y}}{\rho } \hfill \\ \end{gathered} \right.$$
(5)

Based on static equilibrium principles, the resultant force on any cross-section of the beam must equal zero. Therefore, we can derive:

$$\int\limits_{{{A_1}}} {{\sigma _1}{\text{d}}A} +\int\limits_{{{A_2}}} {{\sigma _2}{\text{d}}A} =0$$
(6)

A coordinate system is established with the origin of the Y-axis located at the neutral axis of the composite beam. Assuming that the distance from the interface (between the lower and upper rock beams) to the neutral axis is denoted as h’, the expression for h’ will vary depending on the ___location of the neutral axis, as shown in Fig. 3.

Fig. 3
figure 3

Three conditions for neutral axis position of composite beams.

Full size image

(1)When the neutral axis is located at the interface.

When the neutral axis is positioned precisely at the boundary separating the rock layers, Eq. (6) can be reformulated as follows:

$$\int_{0}^{{{h_1}}} {{E_1}\frac{y}{\rho }} b{\text{d}}y - \int_{0}^{{ - {h_2}}} {{E_2}\frac{y}{\rho }} b{\text{d}}y=0$$
(7)

where b denotes the width of the rock beam.

Under this specific condition, Eq. (7) can be solved accordingly.

$${E_1}{h_1}^{2}={E_2}{h_2}^{2}$$
(8)

Therefore, when Eq. (8) is satisfied, it signifies that the neutral axis of the composite beam aligns precisely with the interface of the rock layers.

(2)When the neutral axis is located in the lower rock beam.

When the neutral axis is positioned within the lower rock beam, Eq. (6) can be reformulated or adjusted to reflect this scenario.

$$\int_{{h^\prime}}^{{h^\prime+{h_1}}} {{E_1}\frac{y}{\rho }} b{\text{d}}y+\int_{0}^{{h^\prime}} {{E_2}\frac{y}{\rho }} b{\text{d}}y - \int_{0}^{{ - \left( {{h_2} - h^\prime} \right)}} {{E_2}\frac{y}{\rho }} b{\text{d}}y=0$$
(9)

Under this specific condition, Eq. (9) can be solved accordingly.

$$h^\prime=\frac{{{E_2}{h_2}^{2} - {E_1}{h_1}^{2}}}{{2\left( {{E_1}{h_1}+{E_2}{h_2}} \right)}}$$
(10)

Considering that h’ > 0, the condition required for the neutral axis of the composite beam to be positioned within the lower rock beam is:

$${E_1}{h_1}^{2}<{E_2}{h_2}^{2}$$
(11)

(3)When the neutral axis is located in the upper rock beam.

When the neutral axis is positioned within the upper rock beam, Eq. (6) can be rewritten as:

$$\int_{0}^{{{h_1} - h^\prime}} {{E_1}\frac{y}{\rho }} b{\text{d}}y - \int_{0}^{{ - h^\prime}} {{E_1}\frac{y}{\rho }} b{\text{d}}y - \int_{{ - h^\prime}}^{{ - \left( {{h_2}+h^\prime} \right)}} {{E_2}\frac{y}{\rho }} b{\text{d}}y=0$$
(12)

Under this condition, the Eq. (12) can be solved:

$$h^\prime=\frac{{{E_1}{h_1}^{2} - {E_2}{h_2}^{2}}}{{2\left( {{E_1}{h_1}+{E_2}{h_2}} \right)}}$$
(13)

Since h’ > 0, the condition required for the neutral axis of the composite beam to be positioned within the upper rock beam is:

$${E_1}{h_1}^{2}>{E_2}{h_2}^{2}$$
(14)

Stress distribution in composite roof

Stress calculation

Based on the methodology outlined in Sect. "Neutral axis position criterion" for determining the neutral axis of a composite beam, the corresponding static equation governing the bending moment of the composite beam can be derived as follows:

$$\int_{{{{\text{A}}_1}}} {{y_1}{\sigma _1}{\text{d}}A} +\int_{{{{\text{A}}_2}}} {{y_2}{\sigma _2}{\text{d}}A} =M\left( x \right)$$
(15)

By substituting Eq. (5) into Eq. (15), we can derive the explicit forms of the static equation for bending moment corresponding to three distinct conditions, each of which is addressed individually.

When the neutral axis is located at the interface of rock layers, the equation is:

$$\int_{0}^{{{h_1}}} {{E_1}\frac{{{y^2}}}{\rho }} b{\text{d}}y+\int_{0}^{{ - {h_2}}} {{E_2}\frac{{{y^2}}}{\rho }} b{\text{d}}y=M\left( x \right)$$
(16)

Solving Eq. (16) gives:

$$\frac{1}{\rho }=\frac{{3 M\left( x \right)}}{{{E_1}b{h_1}^{3}+{E_2}b{h_2}^{3}}}$$
(17)

By substituting Eq. (17) into Eq. (5), we can derive the distribution of normal stress within the beam when the neutral axis is situated at the boundary between rock layers.

$$\left\{ \begin{array}{l}{\sigma _1} = \frac{{ - 3M\left( x \right){E_1}y}}{{{E_1}b{h_1}^3 + {E_2}b{h_2}^3}}\begin{array}{*{20}{c}}{}&{\left( {0 < y \le {h_1}} \right)}\end{array}\\{\sigma _2} = \frac{{ - 3M\left( x \right){E_2}y}}{{{E_1}b{h_1}^3 + {E_2}b{h_2}^3}}\begin{array}{*{20}{c}}{}&{\left( { - {h_2} \le y < 0} \right)}\end{array}\end{array} \right.$$
(18)

When the neutral axis is positioned within the lower rock beam, the static equation governing the bending moment is expressed as:

$$\int_{{h^\prime}}^{{h^\prime+{h_1}}} {{E_1}\frac{{{y^2}}}{\rho }} b{\text{d}}y+\int_{0}^{{h^\prime}} {{E_2}\frac{{{y^2}}}{\rho }} b{\text{d}}y+\int_{{ - \left( {{h_2} - h^\prime} \right)}}^{0} {{E_2}\frac{{{y^2}}}{\rho }} b{\text{d}}y=M\left( x \right)$$
(19)

Solving Eq. (19) gives:

$$\frac{1}{\rho }=\frac{{3 M\left( x \right)}}{{{E_1}b\left( {{h_1}^{3}+3{h_1}h^{{\prime}{2}}+3{h_1}^{2}h^\prime} \right)+{E_2}b\left( {{h_2}^{3}+3{h_2}h^{\prime {2}} - 3{h_2}^{2}h^{\prime}} \right)}}$$
(20)

By substituting Eq. (20) into Eq. (5), we can determine the distribution of normal stress in the beam when the neutral axis is situated within the lower rock beam.

$$\left\{ \begin{gathered} \sigma _{1} = \frac{{ - 3M\left( x \right)E_{1} y}}{{E_{1} b\left( {h_{1} ^{3} + 3h_{1} h^{{{\prime }2}} + 3h_{1} ^{2} h^{\prime }} \right) + E_{2} b\left( {h_{2} ^{3} + 3h_{2} h^{{{\prime }2}} - 3h_{2} ^{2} h^{\prime }} \right)}}\begin{array}{*{20}c} {} & {\left( {h^{\prime } < y \le h^{\prime } + h_{1} } \right)} \\ \end{array} \hfill \\ \sigma _{2} = \frac{{ - 3M\left( x \right)E_{2} y}}{{E_{1} b\left( {h_{1} ^{3} + 3h_{1} h^{{{\prime }2}} + 3h_{1} ^{2} h^{\prime }} \right) + E_{2} b\left( {h_{2} ^{3} + 3h_{2} h^{{{\prime }2}} - 3h_{2} ^{2} h^{\prime }} \right)}}\begin{array}{*{20}c} {} & {\left( {h^{\prime } - h_{2} \le y < h^{\prime }} \right)} \\ \end{array} \hfill \\ \end{gathered} \right.$$
(21)

When the neutral axis is positioned within the upper rock beam, the corresponding static equation that describes the bending moment is:

$$\int_{0}^{{{h_1} - h^\prime}} {{E_1}\frac{{{y^2}}}{\rho }} b{\text{d}}y+\int_{{ - h^\prime}}^{0} {{E_1}\frac{{{y^2}}}{\rho }} b{\text{d}}y+\int_{{ - \left( {{h_2}+h^\prime} \right)}}^{{ - h^\prime}} {{E_2}\frac{{{y^2}}}{\rho }} b{\text{d}}y=M\left( x \right)$$
(22)

Solving Eq. (22) gives:

$$\frac{1}{\rho } = \frac{{3M\left( x \right)}}{{E_{1} b\left( {h_{1} ^{3} + 3h_{1} h^{{{\prime }2}} - 3h_{1} ^{2} h^{\prime }} \right) + E_{2} b\left( {h_{2} ^{3} + 3h_{2} h^{{{\prime }2}} + 3h_{2} ^{2} h^{\prime }} \right)}}$$
(23)

By substituting Eq. (23) into Eq. (5), we can ascertain the distribution of normal stress within the beam when the neutral axis is positioned within the upper rock beam.

$$\left\{ \begin{gathered} \sigma _{1} = \frac{{ - 3M\left( x \right)E_{1} y}}{{E_{1} b\left( {h_{1} ^{3} + 3h_{1} h^{{{\prime }2}} - 3h_{1} ^{2} h^{\prime }} \right) + E_{2} b\left( {h_{2} ^{3} + 3h_{2} h^{{{\prime }2}} + 3h_{2} ^{2} h^{\prime }} \right)}}\begin{array}{*{20}c} {} & {\left( {-h^{\prime } < y \le h_1 - h^{\prime} } \right)} \\ \end{array} \hfill \\ \sigma _{2} = \frac{{ - 3M\left( x \right)E_{2} y}}{{E_{1} b\left( {h_{1} ^{3} + 3h_{1} h^{{{\prime }2}} - 3h_{1} ^{2} h^{\prime }} \right) + E_{2} b\left( {h_{2} ^{3} + 3h_{2} h^{{{\prime }2}} + 3h_{2} ^{2} h^{\prime }} \right)}}\begin{array}{*{20}c} {} & {\left( {-h^{\prime } - h_{2} \le y < -h^{\prime }} \right)} \\ \end{array} \hfill \\ \end{gathered} \right.$$
(24)

Equations (18), (21), and (24) represent the formulations that describe the stress distribution in composite beams under three distinct conditions.

Ultimate load calculation

According to the first strength theory, a point in a beam will undergo failure when the stress at that point satisfies Eq. (25).

$$\sigma \geqslant \left[ {{\sigma _b}} \right]$$
(25)

When the neutral axis of the composite beam is positioned at the interface between rock layers, it is evident that the point where the tensile stress reaches its maximum is located at y=-h2, corresponding to point C in Fig. 4a. At this point, the tensile stress is given by:

$${\sigma _{y= - {h_2}}}=\frac{{3PS{E_2}{h_2}}}{{4\left( {{E_1}b{h_1}^{3}+{E_2}b{h_2}^{3}} \right)}}$$
(26)

Similarly, when the neutral axis is situated in the lower rock beam, the tensile stress at point C (y = h’-h2) in Fig. 4b is:

$${\sigma _{y=h^\prime - {h_2}}}=\frac{{3PS{E_2}\left( {{h_2} - h^\prime} \right)}}{{4\left[ {{E_1}b\left( {{h_1}^{3}+3{h_1}h^{\prime{2}}+3{h_1}^{2}h^\prime} \right)+{E_2}b\left( {{h_2}^{3}+3{h_2}h^{\prime{2}} - 3{h_2}^{2}h^\prime} \right)} \right]}}$$
(27)

When the neutral axis is positioned within the upper rock beam, the positions of maximum tensile stress in the upper and lower rock beams are respectively at y=-h’ and y=-h’-h2. In other words, points B and C in Fig. 4c signify potential failure locations within the composite beam. Under these conditions, the tensile stresses at these two points are calculated as follows:

$$\left\{ \begin{gathered} {\sigma _{y= - h^\prime}}=\frac{{3PS{E_2}h^\prime}}{{4\left[ {{E_1}b\left( {{h_1}^{3}+3{h_1}h^{\prime{2}} - 3{h_1}^{2}h^\prime} \right)+{E_2}b\left( {{h_2}^{3}+3{h_2}h^{\prime{2}}+3{h_2}^{2}h^\prime} \right)} \right]}} \hfill \\ {\sigma _{y= - h^\prime - {h_2}}}=\frac{{3PS{E_2}\left( {h^\prime{\text{+}}{h_2}} \right)}}{{4\left[ {{E_1}b\left( {{h_1}^{3}+3{h_1}h^{\prime{2}} - 3{h_1}^{2}h^\prime} \right)+{E_2}b\left( {{h_2}^{3}+3{h_2}h^{\prime{2}}+3{h_2}^{2}h^{\prime}} \right)} \right]}} \hfill \\ \end{gathered} \right.$$
(28)
Fig. 4
figure 4

Stress distribution within composite beams under three different conditions.

Full size image

Equation (26) to (28) represent the ultimate forces at potential failure points in composite beams that are in an elastic state under three different conditions. By substituting Eq. (25) into Equations (26) to (28), we can determine the loads at which the composite beams reach their elastic limit state under these three conditions.

When the neutral axis is positioned at the interface between rock layers:

$${P_{\hbox{max} }}=\frac{{4{E_1}b{h_1}^{3}+{E_2}b{h_2}^{3}}}{{3{E_2}{h_2}S}} \cdot \left[ {{\sigma _{{{\text{b}}_2}}}} \right]$$
(29)

When the neutral axis is situated within the lower rock beam:

$${P_{\hbox{max} }}=\frac{{4{E_1}b\left( {{h_1}^{3}+3{h_1}h^{\prime{2}}+3{h_1}^{2}h^{\prime}} \right)+{E_2}b\left( {{h_2}^{3}+3{h_2}h^{\prime{2}} - 3{h_2}^{2}h^{\prime}} \right)}}{{3{E_2}\left( {h^\prime - {h_2}} \right)S}} \cdot \left[ {{\sigma _{{{\text{b}}_2}}}} \right]$$
(30)

When the neutral axis is located within the upper rock beam:

$${P_{\hbox{max} }}=\hbox{min} \left\{ \begin{gathered} \frac{{4{E_1}b\left( {{h_1}^{3}+3{h_1}h^{\prime{2}} - 3{h_1}^{2}h^{\prime}} \right)+{E_2}b\left( {{h_2}^{3}+3{h_2}h^{\prime{2}}+3{h_2}^{2}h^{\prime}} \right)}}{{3{E_2}h^{\prime}S}} \cdot \left[ {{\sigma _{{{\text{b}}_1}}}} \right], \hfill \\ \frac{{4{E_1}b\left( {{h_1}^{3}+3{h_1}h^{\prime{2}} - 3{h_1}^{2}h^{\prime}} \right)+{E_2}b\left( {{h_2}^{3}+3{h_2}h^{\prime{2}}+3{h_2}^{2}h^{\prime}} \right)}}{{3{E_2}\left( {h^{\prime}+{h_2}} \right)S}} \cdot \left[ {{\sigma _{{{\text{b}}_2}}}} \right] \hfill \\ \end{gathered} \right\}$$
(31)

Since the layered composite roof structure serves as a key layer in the strata, it bears the load of the overlying rock layers. Therefore, there is an upper limit to the load applied on the composite beam, which should be equal to the self-weight stress of the overlying rock layers within the control range of the key layer. This is represented by Eq. (32).

$${P_{\lim }}=\sum\limits_{{i=1}}^{n} {{\gamma _i}} {H_i}$$
(32)

Where, n represents the number of rock layers controlled above the key layer; γi is the unit weight of the i-th rock layer; Hi is the thickness of the i-th rock layer.

Therefore, the fracture state of the composite roof structure can be determined through the use of Equations (29) to (32). When Plim<Pmax, the load from the overlying rock layers is insufficient to cause the composite roof structure to reach its elastic limit state, and thus, fracture will not occur. When the composite roof undergoes fracture, it indicates that the load imparted by the overlying rock layers has pushed the composite roof to its elastic limit state. Therefore, Plim≥Pmax is a sufficient but not necessary condition for the fracture of the composite roof structure.

Mechanical model of bending energy storage of layered composite roof structure

The composite beam accumulates elastic strain energy when subjected to bending loads, and this stored energy can be mathematically expressed as:

$${U_{\varepsilon b}}=\iiint\limits_{{{{\text{V}}_1}}} {\frac{{{\sigma _1}^{2}}}{{2{E_1}}}{\text{d}}V}+\iiint\limits_{{{{\text{V}}_2}}} {\frac{{{\sigma _2}^{2}}}{{2{E_2}}}{\text{d}}V}$$
(33)

Based on the stress expression derived in Sect. "Stress calculation", the total energy stored in the composite beam can be calculated under three specific conditions.

When the neutral axis is positioned at the interface between rock layers:

$${U_{\varepsilon b}}=\frac{{bS}}{{2{E_1}}}\int\limits_{0}^{{{h_1}}} {\int\limits_{0}^{S} {{{\left[ {\frac{{3 M\left( x \right){E_1}y}}{{{E_1}b{h_1}^{3}+{E_2}b{h_2}^{3}}}} \right]}^2}} } {\text{d}}x{\text{d}}y+\frac{{bS}}{{2{E_2}}}\int\limits_{0}^{{{h_2}}} {\int\limits_{0}^{S} {{{\left[ {\frac{{3 M\left( x \right){E_2}y}}{{{E_1}b{h_1}^{3}+{E_2}b{h_2}^{3}}}} \right]}^2}} } {\text{d}}x{\text{d}}y$$
(34)

When the neutral axis is situated within the lower rock beam:

$$\begin{gathered} {U_{\varepsilon b}}=\frac{{bS}}{{2{E_1}}}\int\limits_{0}^{{{h_1}}} {\int\limits_{0}^{S} {{{\left[ {\frac{{3 M\left( x \right){E_1}y}}{{{E_1}b\left( {{h_1}^{3}+3{h_1}h^{\prime {2}}+3{h_1}^{2}h\prime} \right)+{E_2}b\left( {{h_2}^{3}+3{h_2}h^{\prime {2}} - 3{h_2}^{2}h^{\prime}} \right)}}} \right]}^2}} } {\text{d}}x{\text{d}}y+ \hfill \\ \begin{array}{*{20}{c}} {^{{_{{}}}}}&{_{{{}_{{{}_{{\mathop {}\nolimits_{{}}^{{}} }}^{{}}}}}}^{{}}} \end{array}\frac{{bS}}{{2{E_2}}}\int\limits_{0}^{{{h_2}}} {\int\limits_{0}^{S} {{{\left[ {\frac{{3 M\left( x \right){E_2}y}}{{{E_1}b\left( {{h_1}^{3}+3{h_1}h{\prime 2}+3{h_1}^{2}h^{\prime}} \right)+{E_2}b\left( {{h_2}^{3}+3{h_2}h^{\prime {2}} - 3{h_2}^{2}h^{\prime}} \right)}}} \right]}^2}} } {\text{d}}x{\text{d}}y \hfill \\ \end{gathered}$$
(35)

When the neutral axis is located within the upper rock beam:

$$\begin{gathered} {U_{\varepsilon b}}=\frac{{bS}}{{2{E_1}}}\int\limits_{0}^{{{h_1}}} {\int\limits_{0}^{S} {{{\left[ {\frac{{3 M\left( x \right){E_1}y}}{{{E_1}b\left( {{h_1}^{3}+3{h_1}h^{\prime {2}}-3{h_1}^{2}h\prime} \right)+{E_2}b\left( {{h_2}^{3}+3{h_2}h^{\prime {2}} + 3{h_2}^{2}h^{\prime}} \right)}}} \right]}^2}} } {\text{d}}x{\text{d}}y+ \hfill \\ \begin{array}{*{20}{c}} {^{{_{{}}}}}&{_{{{}_{{{}_{{\mathop {}\nolimits_{{}}^{{}} }}^{{}}}}}}^{{}}} \end{array}\frac{{bS}}{{2{E_2}}}\int\limits_{0}^{{{h_2}}} {\int\limits_{0}^{S} {{{\left[ {\frac{{3 M\left( x \right){E_2}y}}{{{E_1}b\left( {{h_1}^{3}+3{h_1}h{\prime 2}-3{h_1}^{2}h^{\prime}} \right)+{E_2}b\left( {{h_2}^{3}+3{h_2}h^{\prime {2}} + 3{h_2}^{2}h^{\prime}} \right)}}} \right]}^2}} } {\text{d}}x{\text{d}}y \hfill \\ \end{gathered}$$
(36)

In order to simplify the expression of the model, let’s define the ratio of the thickness between the upper and lower rock beams as λ, and the ratio of their elastic moduli as η. Mathematically, this can be expressed as:

$$\lambda =\frac{{{h_2}}}{{{h_1}}}$$
(37)
$$\eta =\frac{{{E_2}}}{{{E_1}}}$$
(38)

By substituting Eqs. (3), (37), and (38) into Equations (34) to (36), a simplified mechanical model for the bending energy storage of the Layered Composite Rock Structure (LCRS) is obtained.

When the neutral axis is positioned at the interface between rock layers:

$${U_{\varepsilon b}}=\frac{{3{S^3}{{\left( {1+\lambda } \right)}^3}\left( {1+\eta } \right){P^2}}}{{64b\left( {1+{\lambda ^3}\eta } \right){E_1}{h_1}^{3}}}$$
(39)

When the neutral axis is situated within the lower rock beam:

$${U_{\varepsilon b}}=\frac{{3{S^3}\left( {1+\eta } \right){P^2}}}{{64b{E_1}{h_1}^{3}{{\left[ {1+{\lambda ^3}\eta - \frac{{3{{\left( {1 - {\lambda ^2}\eta } \right)}^2}}}{{4\left( {1+\lambda \eta } \right)}}} \right]}^2}}}$$
(40)

When the neutral axis is located within the upper rock beam:

$${U_{\varepsilon b}}=\frac{{3{S^3}\left( {1+\eta } \right){P^2}}}{{64b{E_1}{h_1}^{3}{{\left[ {1+{\lambda ^3}\eta +\frac{{9{{\left( {1 - {\lambda ^2}\eta } \right)}^2}}}{{4\left( {1+\lambda \eta } \right)}}} \right]}^2}}}$$
(41)

Equation (39) to (41) represent the bending energy storage models of the layered composite roof structure under three different conditions. The significance of these models lies in their ability to determine the stored elastic strain energy of composite beams, based on their intrinsic properties and stress state. Since the denominators of Eqs. (40) and (41) are always greater than 0, Eq. (40) is always greater than Eq. (41) when all variables take the same values. This indicates that under the same conditions of lithology, thickness, and loading, more energy will be stored when the neutral axis is located in the lower rock beam.

Numerical and experimental verification of model accuracy

Numerical simulation

To verify the accuracy of the theoretical model established above, a numerical model of FLAC3D containing a double-layered rock beam was established, as shown in Fig. 5. The elastic model was used in the numerical validation and the variables studied were the height and the elastic modulus of the rock beams, i.e., h and E. Five sets of numerical simulation schemes were developed, as shown in Table 1.

Fig. 5
figure 5

Schematic diagram of the numerical model.

Full size image
Table 1 Numerical simulation scheme.
Full size table

Figures 6 and 7 depict the stress distribution diagram and strain energy density distribution diagram, respectively, for both the theoretical model and the numerical simulation. When compared under identical conditions, the distribution patterns of stress and strain energy density in both the theoretical model and the numerical simulation results exhibit a remarkable consistency, featuring comparable peak stress values. This alignment underscores the high precision of the stress calculation equation employed within the theoretical model.

Fig. 6
figure 6

Stress distribution of composite beam.

Full size image

When the elastic moduli of the upper and lower rock beams are equivalent, the composite beam can be regarded as a unified thick beam, with the neutral axis positioned at the geometric centroid. Consequently, the distribution of tensile and compressive stress regions within the composite beam exhibits perfect symmetry, and the strain energy density is also symmetrical about the interface between the rock layers, as illustrated in Figs. 6a and b and 7a, and 7b.

When the thicknesses of the upper and lower beams are identical, and the elastic modulus of the upper beam exceeds that of the lower beam, the neutral axis of the composite beam shifts to within the upper rock beam. Consequently, the tensile stress region in the beam becomes larger than the compressive stress region, and the strain energy density in the lower beam is higher than that in the upper beam, as depicted in Figs. 6c and d and 7c and d. If the elastic modulus of the upper beam is smaller than that of the lower beam, the neutral axis of the composite beam shifts to within the lower rock beam. Consequently, the range of compressive stress in the beam becomes greater than that of tensile stress, and the strain energy density in the upper beam exceeds that in the lower beam, as illustrated in Fig. 6e and f, and Fig. 7e and f.

When the thickness and elastic modulus of the upper and lower beams are different, the determination of the neutral axis’s position necessitates the use of the method described in Sect. "Neutral axis position criterion", as exemplified in Figs. 6g and j and 7g and j.

Fig. 7
figure 7

Strain energy density distribution of composite beam.

Full size image

Experimental analysis

To validate the rationality of the energy calculation equation from the view of experiment, a three-point bending test of the composite beam was carried out. The composite beam specimen consists of limestone and fine sandstone, taken from the roof area of a coal mine located in the west of Shandong Province, and its basic mechanical parameters are shown in Table 2. First, large blocks of rock were transported to the laboratory and then cut into specimens of 200 mm in length, 40 mm in width and different heights. Finally, the specimens were bonded together with gypsum to form rock beam specimens of different combination of types, and the numbering and dimensions of the composited specimens are shown in Table 3.

Table 2 Basic physical and mechanical parameters of composite beams.
Full size table
Table 3 The numbering and dimensions of the composited specimens.
Full size table

The test system consists of a loading system and a digital speckle correlation measurement system (DSCM), as shown in Fig. 8. During the test, the rock beam span was set to be 150 mm, and displacement loading was applied with a loading rate of 0.05 mm/min. The DSCM system is used to obtain the specimen deformation field, which is consequently utilized for energy field calculations. A CCD camera was employed for image acquisition with an effective pixel resolution of 2048 × 1536 and a sampling frequency of 5 Hz. Meanwhile, the specimen observation surface was artificially treated with scattering before the laboratory test.

Fig. 8
figure 8

Three-point bending test system.

Full size image

Based on DSCM data, the energy evolution characteristics during the deformation and failure process of the rock beam have been calculated, with the strain energy density and strain energy being derived from Eqs. (42) and (43), respectively.

$$U=\frac{E}{2}{\left( {{\varepsilon _1}^{2}+{\varepsilon _2}^{2} - 2\mu {\varepsilon _1}{\varepsilon _2}^{2}} \right)^2}$$
(42)
$$W=\iiint\limits_{{{{\text{V}}_1}}} {{U_1}{\text{d}}V}+\iiint\limits_{{{{\text{V}}_2}}} {{U_2}{\text{d}}V}$$
(43)

Where, E denotes the elastic modulus of the rock sample; µ represents the Poisson’s ratio of the rock sample; ε1 and ε2 are the principal strains on the surface of the rock sample.

Figure 9 presents the evolution of strain energy for different rock beam specimens. As shown in the figure, the strain energy evolution during the loading of composite beams can be broadly classified into two types: J-type and S-type. This may be related to the sequence of fracture in the rock beams. When the composite beam contains a thick and hard layer, the low-strength rock beam fractures first, leading to a continuous increase in the overall strain energy. Subsequently, the fracture of the high-strength beam causes the overall strain energy to reach its peak, as shown in Fig. 9a and d, and f. In contrast, for composite beams without a thick and hard layer, both layers of rock beams fracture almost simultaneously during loading, resulting in a sudden increase in the overall strain energy to its peak, as illustrated in Fig. 9b and c, and e.

Fig. 9
figure 9

The evolution of strain energy for different rock beam specimens.

Full size image

Based on the bending energy storage models of composite beams derived in Sect. "Mechanical model of bending energy storage of layered composite roof structure", Table 4 presents the theoretical calculations and experimental results for the strain energy stored in the composite beams during the six tests. According to the theoretical calculations, the bending strain energy stored prior to the fracture of the composite beams is closely linked to the presence of thick and hard layers. Specifically, composite beams with thick and hard layers store more energy before fracture compared to those without such layers. Notably, the errors between five out of the six test results and the theoretical model proposed in this study are below 7%, while the error of the remaining result is 11.8%. This further validates the accuracy of the theoretical model.

Table 4 Theoretical calculation of energy storage for composite beams.
Full size table

Conclusion

(1) Considering the locations of the neutral layer which is located at the upper rock beam, lower rock beam and interface of rock layers, the bending strain energy calculation of the LCRS has been modeled. The stress distribution, ultimate load and strain energy calculation equations for the LCRS in the above three states are given.

(2) The distribution of bending strain energy is closely related to the ___location of the neutral axis and the presence of a thick hard layer. When the neutral axis is located in the lower part of the rock beam, more bending strain energy is stored, and LCRS with a thick hard layer can stores more energy before fracture than LCRS without a thick hard layer.

(3) The accuracy of the model is verified through numerical simulations and three-point bending tests. And the error between the five sets of test results and the theoretical model calculations is less than 7%. However, the work does not consider the influence of interfacial shear interaction. Therefore, it is necessary to consider the above factor for further research.