The stress of the bolt is shown in Figure 2. T ( z ) shows the axial force of the bolt; τ ( z ) shows the shear stress on the surface of the bolt body; S b ' represents the interface between the anchorage section and the non-anchorage section; and S b shows the arbitrary section in the anchorage section. The axial force T ( z ) of the bolt is constant in the non-anchorage section and its value is the same as the working resistance T 0 . In the anchorage section, it is an unknown function about z. The shear stress τ ( z ) is 0 in the non-anchorage section and an unknown function about z in the anchorage section.

The bolt body from the end of the anchor to the section S b is taken as the research object. In the cylindrical coordinate system, the direction of shear stress in Figure 2 is negative. According to the balance condition of the force system, the relationship between the axial force T ( z ) and shear stress of τ ( z ) can be obtained as follows: T ( z ) = T 0 + 2 π r b ∫ 0 z τ ( z ) d z . (1)

In Eq. 1, r b represents the radius of the bolt. The elastic modulus of the bolt is much larger than that of the surrounding rock and the cross-section of the bolt is much smaller than that of the roadway. Therefore, it can be assumed that the normal stress σ ( z ) on the cross-section of the bolt is evenly distributed. Then, the relationship between the axial force T ( z ) and normal stress σ ( z ) is T ( z ) = π ⋅ ( r b ) 2 ⋅ σ ( z ) . (2)

According to the field observation data and the design requirements of the bolt support, the bolts are in the elastic state during roadway driving. Therefore, the relationship between the axial normal stress σ ( z ) and the axial strain ε ( z ) of the bolt is σ ( z ) = E ⋅ ε ( z ) . (3)

In Eq. 3, E represents the elastic modulus of the bolt. Substituting Eq. 3 into Eq. 2, we can get T ( z ) = π E ⋅ ( r b ) 2 ⋅ ε ( z ) . (4)

By substituting Eq. 4 into Eq. 1, the relationship between the bolt surface shear stress τ ( z ) and the bolt axial strain ε ( z ) is 2 π r b ⋅ ∫ 0 z τ ( z ) d z − π E ⋅ ( r b ) 2 ⋅ ε ( z ) + T 0 = 0. (5)

The shear stress and strain of the surrounding rock at the anchorage surface are the same as the shear stress and axial strain of the bolt. Then, the shear stress and strain of the surrounding rock at the anchorage surface can be substituted into Eq. 5. By solving Eq. 5, the expressions of axial force and shear stress of the bolt are derived.

Assuming that the extrusion force of the bearing plate on the surrounding rock is uniformly distributed and there is no shear load on the surrounding rock, the force of the bearing plate on the surrounding rock is shown in Figure 3:

In Figure 3, the calculation method of p t and r t are as follows: r t = S t π ,       p t = T 0 S t . (6)

In Eq. 6, S t represents the area of the bearing plate. When the Love displacement function is taken in the form as shown in Eq. 7, the exact solution of the mechanical problem shown in Figure 3 can be obtained. Φ r = ∫ 0 z ( ϕ 0 + z ⋅ ϕ 3 ) d z (7)

In Eq. 7, ϕ 0 and ϕ 3 are displacement functions, and the expression is { ϕ 0 = ( 1 − 2 ν ) p t r t ⋅ ∫ 0 + ∞ J 1 ( r t λ ) ⋅ J 0 ( ρ λ ) ⋅ e − λ z λ 2 d λ , ϕ 3 = − p t r t ⋅ ∫ 0 + ∞ J 1 ( r t λ ) ⋅ J 0 ( ρ λ ) ⋅ e − λ z λ d λ . (8)

In Eq. 8, J 0 () and J 1 () represent the first kind of Bessel functions of order 0 and order 1; v represents Poisson’s ratio of the surrounding rock; e represents the base of the natural logarithm; and λ represents the integral variable.

By substituting the Love displacement function Eq. 7 into Eq. 9, the displacement distribution law of the surrounding rock when it is squeezed by the bearing plate can be obtained. { u ρ = 1 2 G ⋅ ∂ 2 Φ ∂ ρ ∂ z u θ = 0 u z = ν − 1 G ⋅ ∇ 2 Φ + 1 2 G ⋅ ∂ 2 Φ ∂ z 2 (9)

In Eq. 9, Φ represents the Love displacement function and G represents the shear modulus of the surrounding rock. By substituting the displacement component of the surrounding rock obtained from Eq. 9 into the geometric equation, the strain tensor of the surrounding rock under the bearing plate extrusion can be obtained, and then the stress tensor of the surrounding rock can be obtained by Hooke’s law.

Assuming that the stress distribution of roadway surrounding rock is a plane strain problem, the proposed complex function method is used to solve the strain distribution law of rectangular roadways (Muskhelishvili and Noordhoff, 1953; Feng et al., 2014; Tran Manh et al., 2015; Shen et al., 2017; Chang et al., 2020). The stress of a rectangular roadway is shown in Figure 4.

The stress distribution of the surrounding rock of a rectangular roadway can be characterized by two complex functions φ ( ξ ) and ψ ( ξ ) . The form of complex functions φ ( ξ ) and ψ ( ξ ) are { φ ( ξ ) = σ V + σ H 4 ⋅ ω ( ξ ) + φ 0 ( ξ ) , ψ ( ξ ) = ( σ V − σ H ) 2 + 4 τ 2 ⋅ e 2 i α ω ( ξ ) 4 + ψ 0 ( ξ ) . (10)

In Eq. 10, the function ω ( ξ ) represents the conformal mapping function from the outer domain of the rectangular roadway to the unit circle on the complex plane, and the solution method is shown in the literature (Nazem et al., 2015; Yuan et al., 2018; Baddoo and Crowdy, 2019; Badreddine et al., 2019; He et al., 2022). Analytic functions φ 0 ( ξ ) and ψ 0 ( ξ ) satisfy Cauchy–Riemann conditions. The values of variables ξ and α are as follows: { ξ = η e i β , α = arctan ( σ V − σ H + ( σ V − σ H ) 2 + 4 τ 2 2 τ ) . (11)

In Eq. 11, η and β represent the coordinate components of the curvilinear coordinate system determined by the conformal mapping function ω ( ξ ) . The surrounding rock stress of the rectangular roadway is { σ η r + σ β r = 4 Re ( φ ′ ( ξ ) ω ′ ( ξ ) ) , σ η r − σ β r + 2 i τ η β r = 2 ξ 2 η 2 ⋅ [ ω ( ξ ) ¯ ω ′ ( ξ ) ¯ ⋅ ( φ ′ ( ξ ) ω ′ ( ξ ) ) ′ + ψ ′ ( ξ ) ω ′ ( ξ ) ¯ ] . (12)

In Eq. 12, σ η , σ β , and τ η β represent the stress components in the curvilinear coordinate system. According to the method in the reference (He et al., 2022), the stress tensor in the rectangular coordinate system can be obtained. Then, the strain tensor of the surrounding rock can be obtained from the stress tensor.

The load function P b ( z ) is used to equivalent the effect of the bolt on the surrounding rock. The load function is distributed along the axial direction of the bolt. It is an unknown function of variable z in the anchorage section and 0 in the non-anchorage section. The effect of the load on the surrounding rock is the same as that of the bolt on surrounding rock. The force of the load on the surrounding rock is shown in Figure 5.

Take microelements on the bolt for analysis. At this time, the surrounding rock z = δ is under the action of the concentrated load P b ( z ) d z , and when the Love displacement function is in the form shown in Eq. 13, the stress distribution law of the surrounding rock under the concentrated load P b ( z ) d z can be obtained by the following equation. Φ = P b ( δ ) d δ ⋅ F ( z , ρ , δ ) . (13)

In Eq. 13, the function F ( z , r , δ ) is a known function and satisfies the following equation: 8 π ⋅ ( 1 − ν ) ⋅ F ( z , ρ , δ ) = R 1 + R 2 ⋅ ( 8 ν − 8 ν 2 − 1 ) + ( 4 − 8 ν ) ( z − ν z − ν δ ) ⋅ log ( R 2 + z + δ ) − 2 δ z R 2 . (14)

In Eq. 14, R 1 and R 2 are { R 1 = ρ 2 + ( z − δ ) 2 , R 2 = ρ 2 + ( z + δ ) 2 . (15)

Integrate Eq. 13 on the bolt body to obtain the Love displacement function: Φ b = ∫ 0 L P b ( δ ) ⋅ F ( z , ρ , δ ) d δ . (16)

Through Eq. 16, the effect of the distribution load P b ( z ) on the surrounding rock can be obtained, that is, the effect of the bolt on the surrounding rock. Substituting Eq. 16 into Eq. 9 can obtain the displacement component of the surrounding rock under the action of the bolt, and then obtain the stress and strain tensor of the surrounding rock under the action of the bolt as follows: { σ i j b = ∫ 0 L P b ( δ ) ⋅ σ i j 0 d δ , ε i j b = ∫ 0 L P b ( δ ) ⋅ ε i j 0 d δ . (17)

In Eq. 17, the sum of tensors σ i j 0 and ε i j 0 is a known quantity, which is only related to the function F ( z , r , δ ) .

According to the uniqueness theorem of solution in elasticity, there is only one exact solution satisfying the corresponding boundary conditions. The contact surface between the bolt and the surrounding rock can be regarded as a boundary condition. When the surrounding rock is taken as the research object, it has two boundaries. The first is the free surface of the surrounding rock and the second is the contact surface between the bolt and the surrounding rock, that is, the anchorage surface. When the two boundary conditions are consistent, the stress-strain state in the surrounding rock is unique and determined. The boundary condition of the surrounding rock at the free surface is not affected by the bolt. When the anchorage surface does not slide, the stress and strain are continuous on the anchorage surface. The stress and strain on the anchorage surface meet both the mechanical equation of the surrounding rock and the mechanical equation of the bolt. Therefore, the stress and strain of the surrounding rock at the anchorage surface can be substituted into the mechanical equation of the bolt, and the functional equation with the distributed load function P b ( ι ) as the unknown function can be obtained. The solution of the unknown function P b ( ι ) can be obtained by solving this equation. The calculation formulas of the axial force and shear stress of the bolt can be obtained by substituting the obtained function P b ( ι ) into the relevant formulas.

Substituting the values of stress and strain at the anchorage surface into Eq. 5, we can get 2 π r b ⋅ ∫ 0 z τ r z m | r = r b d z − π E ⋅ ( r b ) 2 ⋅ ε z m | r = r b + T 0 = 0. (18)

In Eq. 18, τ r z m and ε z m are { τ r z m = τ r z r + τ r z t + τ r z b , ε z m = ε z r + ε z t + ε z b . (19)

Substitute Eq. 17 into Eq. 18 and simplify it to obtain ∫ 0 L P b ( δ ) ⋅ K ( z , δ ) d δ = f ( z ) (20)

Eq. 20 is the first kind of Fredholm integral equation (Mesgarani and Azari, 2019; Khan et al., 2020), in which the function K ( z , δ ) is called the kernel of the integral equation and the function f ( z ) is called the free term of the integral equation. The kernel function and free term are known functions and their expressions are { K ( z , δ ) = ∫ 0 z τ r z 0 | r = r b d z − E r b 2 ⋅ ε z 0 | r = r b , f ( z ) = E r b 2 ⋅ ( ε z r + ε z t ) | r = r b − T 0 2 π r b − ∫ 0 z ( τ r z r + τ r z t ) | r = r b d z . (21)

According to the relevant theory of integral equation, the outgoing load function P b ( z ) can be solved from Eq. 20. Substituting the load function obtained from the solution into Eq. 19 and then substituting the second equation of Eq. 19 into Eq. 4, the expression of the bolt axial force can be obtained T ( z ) = π E ⋅ ( r b ) 2 ⋅ [ ∫ 0 L P b ( δ ) ε z 0 d δ + ε z t + ε z r ] | r = r b . (22)

Substituting Eq. 22 into Eq. 1, and then differentiating and sorting Eq. 1, the expression of the surface shear stress τ ( z ) can be obtained by τ ( z ) = E r b 2 ⋅ { ∫ 0 L P b ( δ ) ⋅ d ε z 0 d z ⋅ d δ + d d z ( ε z t + ε z r ) } | r = r b . (23)

By substituting the stress variables of the surrounding rock under different engineering conditions into Eq. 18, the distribution curves of the bolt axial force and shear stress under the corresponding engineering conditions can be obtained.

The schematic diagram of a driving roadway is shown in Figure 6. When the roadheader cuts out the complete roadway section, some bolts are installed immediately to support the roadway, and its position is shown at point A of Figure 6. At this time, the surrounding rock is supported by the front coal wall without deformation or the deformation is very small, which can be ignored compared with the deformation of the surrounding rock in the later stage. When it is far away from the coal wall, as shown at point B, the bolt is affected by the deformation of the surrounding rock and the axial force of the bolt changes (Chang et al., 2021).

L AB is used to indicate the influence range of the driving face. It is assumed that at point A, the surrounding rock is not deformed, and at point B, the surrounding rock reaches a stable state. The creep of the surrounding rock and other factors are not considered in this study. According to the numerical simulation results and on-site roadway deformation observation data, the exponential function is used to describe the strain of the surrounding rock between point A and point B. ε i j s = ε i j r ( 1 − e − a 1 L A B ⋅ l ) . (24)

In Eq. 24, the parameter a ₁ reflects the severity of the surrounding rock deformation within the influence range of the driving face; l represents the distance between the bolt and the working face; ε i j r represents the strain tensor of the surrounding rock when it is not affected by the driving face; and ε i j s represents the strain tensor of the surrounding rock within the influence range of the driving face. By substituting ε i j s into Eqs 21, 22, 23, the distribution laws of the axial force and shear stress of bolts with different anchor lengths in a rectangular roadway during roadway driving can be obtained. Reference Chang et al. (2020) gives the stress distribution law of the roadway surrounding rock during working face mining. Under the influence of mining stress of the working face, the vertical stress of the roadway can be expressed as σ V ( 1 ) = H ( 1 s − 1 ) ⋅ k s ⋅ σ V 0 e a 2 ⋅ 1 s − 1 ⋅ ( e a 2 ⋅ 1 − 1 ) + H ( 1 − 1 s ) ⋅ σ V 0 ⋅ [ ( k s − 1 ) ⋅ e b 2 ⋅ ( 1 s − 1 ) + 1 ] . (25)

In Eq. 25, function H( ) represents the Heaviside step function; l represents the distance from the working face; σ V 0 represents the vertical pressure of the original rock; k _s represents the vertical stress concentration coefficient at the peak of the abutment pressure; l _s represents the peak position of the abutment pressure; and coefficients a ₂ and b ₂ are used to describe the change severity of the abutment pressure curve. We can obtain the horizontal stress σ H and shear stress τ through the pressure measurement coefficient λ and shear stress coefficient η . Similar to the influence of the driving face on the bolt, we can obtain the stress distribution of the bolt under the influence of the working face through Eq. 25.

By substituting the relevant data into the bolt mechanical model, the variation curve of the axial force of a full-length anchoring bolt in a rectangular roadway during roadway driving and working face mining can be obtained. The variation curve of the axial force of the bolts A and B with the advance of driving face and working face is shown in Figure 10.

It can be seen from Figure 10 that for the full-length anchoring bolt in the middle of the roof and side of the rectangular roadway, the distribution law of its axial force along the bolt body direction is roughly the same. With the advance of the driving face, the axial force rapidly evolves from a monotonous decreasing trend to the change law of first increasing and then decreasing. But with the advance of the working face, the axial force increases first and then decreases. The axial force change of the bolt in the middle of the roof is gentler than that in the middle of the side. The maximum axial force point of the bolt quickly stabilizes at the neutral point from the orifice position with the advance of the driving face. The neutral point of the bolt in the middle of the roof is 0.75 m away from the orifice, and the neutral point of the bolt in the middle of the side is 0.40 m away from the orifice. The neutral point of the bolt on the roadway side is closer to the roadway surface. The working resistance of the bolt shows a monotonous increasing trend with the advance of the driving face and the working face.

According to Eq. 23, the distribution curve of the shear stress of the bolt body can be obtained. The evolution curve of the shear stress of the bolt body with the advance of the driving face and the working face is shown in Figure 11. In Figure 11, the sub coordinate system is a local magnification of the main coordinate system.

It can be seen from Figure 11 that the evolution law of the shear stress in the middle of the roof and the bolt in the middle of the side is the same with the advance of the driving face and the working face. When there is no neutral point in the bolt, that is, when the bolt is installed in the surrounding rock, the shear stress of the bolt shows a monotonous increasing trend. When there is a neutral point in the bolt, that is, when the bolt is far away from the driving face, the shear stress of the bolt shows a monotonous decreasing trend. The variation range of the shear stress of the bolt in the middle of the side is greater than that of the bolt in the middle of the roof. The shear stress curve of the bolt intersects at one point, that is, the neutral point of the bolt, which indicates that the neutral point position of the bolt does not change during roadway driving and working face mining.

Comparing Figures 10, 11, it can be seen that when the shear stress of the bolt is less than 0, the axial force of the bolt decreases. When the shear stress of the bolt is greater than 0, the axial force of the bolt increases. On both sides of the neutral point of the bolt, the sign of the shear stress of the bolt is different, which is consistent with the neutral point theory.

The mining roadway of 1,131(1) working face is driving for 4.8 m every day. After installing the force measuring bolt, it is observed twice on the first day and once every day, after that, for a total of 7 days. The observed data are shown in Figure 12.

It can be seen from Figure 12 that the axial force distribution law of the bolt measured by the force measuring bolt is basically consistent with the axial force distribution law obtained by theoretical calculations. When the force measuring bolt is installed into the surrounding rock, the axial force of the bolt shows a monotonous decreasing trend. When it is far from the driving face, the axial force of the bolt increases first and then decreases. The maximum value is reached at the neutral point of the bolt. The neutral point of the roof force measuring bolt is between 0.52 and 1.04 m. The neutral point of the side force measuring bolt is about 0.46 m, which is consistent with the theoretical calculation results. The variation range of the axial force of the force measuring bolt in the middle of the roadway side is greater than that of the force measuring bolt in the middle of the roof. The working resistance of the force measuring bolt in the middle of the roadway side is less than that of the force measuring bolt in the middle of the roof, which is consistent with the theoretical calculation results. The correctness of the theoretical calculation results can be verified from the observation data of the force measuring bolt.