1 Introduction

The stability of underground excavations directly influences the safety of underground operators, inducing a threat on the production of mining activities. To improve the integrity of rock masses and maintain the safety of underground openings, rock support and rock reinforcement methods are always used (Kang 2014; Wu et al. 2017; Zhang et al. 2020).

Windsor (1997) indicated that rock reinforcement represents the techniques or devices that are used to improve the rock mass properties from within rock masses, such as rock bolts and ground anchors. In this paper, the authors concentrated on studying the axial pulling behaviour of fully encapsulated rock bolts.

Since fully encapsulated rock bolts were used in mining engineering, it saw a wide application of rock bolts in roadways, entries, chambers and stopes (He et al. 2020; Wang et al. 2020; Zhang et al. 2020a). The installed rock bolts, the bonding grout in the borehole and the rock masses constitute the rock reinforcement system. However, the rock masses have a tendency to move towards the excavated area. Then, relative displacement occurs between the rock masses and the rock bolt. This leads to occurring of shear stresses in the grout annulus. As a result, the force can be transferred among rock bolts and rock masses (Chen et al. 2016b, 2019).

To evaluate the loading performance among rock bolts and rock masses, numerical simulation is becoming more widely used. A special finite element was developed by Mitri and Rajaie (1990). Their study indicated that with the rock reinforcement elements, the stability and safety level of underground excavations was largely improved. Ivanovic et al. (2003) developed a lumped parameter model to evaluate the performance of rock bolts. This numerical model was successfully confirmed with experimental work. Moosavi and Grayeli (2006) numerically modelled the bond failure in the rock reinforcement system and coupled it into the discontinuous deformation analysis (DDA) algorithm. A numerical pulling test was conducted to evaluate the maximum bearing capacity and the load distribution of the bolt. However, the bolt and the blocks were assumed to be two-dimensional. Deb and Das (2011) used the enriched finite element method (EFEM) to study the mechanical performance of fully encapsulated rock bolts. Based on the EFEM, the displacement and the stress in the rock bolt tendon can be simulated. Credibility of this numerical model was successfully validated. Li et al. (2012) used the commercial software FLAC3D to simulate the performance of rock bolts. The structural element of cable was used. More importantly, the input parameters regarding the cable were successfully calibrated through back analysis Nemcik et al. (2014). revised the structural element in the two-dimensional program FLAC to simulate the failure propagation in a rock bolting system. The yielding behaviour of the interface between the structural element and the mesh can be set manually. Then, the interfacial yielding behaviour can be simulated. However, in this case, users should input the shear coupling data manually, making the simulation complicated. Amarnath and Sitharam (2016) used the pile element in FLAC3D to simulate the pulling behaviour of rock bolts. Experimental field tests were performed in a rock slope. Then, the modelling quality was certified. Ma et al. (2016) studied the shear coupling behaviour of the interface between the bolt and the grout and developed a numerical model to simulate the performance of rock bolts. It was found that the shear stress distributed non-linearly along the embedment length. Bahrani and Hadjigeorgiou (2017) compared the performance of rock bolts simulated with the structural elements of cable and rockbolt in UDEC. The results indicated that compared with the rockbolt element, the cable element is simpler in terms of the input parameters. Moreover, those input parameters can be acquired from the laboratory pulling experiments. Chang et al. (2017) used the finite element program to simulate the grout cracking and debonding of rock bolts. Experimental short embedment pulling tests were used to validate the numerical model. Yan et al. (2018) studied the rock bolting system with numerical modelling in FLAC3D. In the modelling process, axial failure of the bolt itself was investigated. Results revealed that once elongation of the rock bolt was higher than the critical limit, rupture of the rock tendon would occur, leading to failure of the rock bolting system. Yokota et al. (2019) adopted the DDA method to study the behaviour of the interface between a rock bolt and the surrounding grout. The results showed that the rock bolt rib angle had a significant effect on the performance of rock bolts. Ho et al. (2019) investigated the bond failure mechanism in the anchorage body. It was found that when the rock bolt was subjected to the high confining stress condition, the bolt/grout interface failed by shearing.

The literatures above made a significant contribution in revealing the load transfer mechanism of rock bolts. Moreover, it was beneficial for improving the optimisation of rock bolts and safety of mining production. Although numerical simulation has its advantages in solving complex geotechnical problems, it still has certain shortcomings. For instance, in FLAC3D, the cableSELs are always used to simulate rock bolts and cable bolts. However, in the original constitutive law of the cableSELs, the grout annulus was treated as an elastic perfectly-plastic material. Due to this shortcoming, the interface shear stress between the rock bolt and the grout material would not decay after arriving at the interface shear strength (Fig. 1). Consequently, the softening behaviour of the interface between the bolt and the grout was not described properly.

Fig. 1
figure 1

Shear behaviour of the interface between the bolt and the grout in the original cableSELs, after Itasca (2016)

Full size image

Considering this issue, in this paper, the authors used a two-stage shear coupling model to depict the shear behaviour of the interface between the bolt and the grout. Then, an approach was proposed to incorporate this two-stage shear coupling model into FLAC3D to modify the current cableSELs. After it, a numerical pulling test was conducted with the modified cableSELs. The numerical pulling results were validated with experimental results, confirming the credibility of the modified cableSELs. Last, based on the modified cableSELs, the influence of bolt diameter and grout shear strength on rock bolt performance was studied. When pulling force increased, the shear stress distribution at the bolt/grout interface along the bolt was investigated.

2 The shear coupling model for the interface

The purpose of the shear coupling model was to determine the interface shear stress between the rock bolt and grout material. In this study, specifically, a two-stage model was used (Pan and Wu 2014). Initially, as the shear slippage increased, the interface shear stress between the rock bolt and the grout material rose until the interface shear strength. After that, the shear stress decreased exponentially, which is depicted with Eqs. (1)–(2) and Fig. 2:

$$\tau = \frac{{\tau_{\max } }}{{\delta_{0} }}\delta \, \left( {\delta < \delta_{0} } \right)$$
(1)
$$\tau = \frac{{\tau_{\max } }}{{e^{{c\left( {\delta - \delta_{0} } \right)}} }} \, \left( {\delta \ge \delta_{0} } \right)$$
(2)

where, τ is the interface shear stress between the rock bolt and the grout material (Pa); τmax is the interface shear strength (Pa); δ is the relative shear slippage between the rock bolt and the grout material (m); δ0 is the relative shear slippage between the rock bolt and the grout material when the interface shear strength was reached (m); c is the coefficient for the plastic stage (m−1).

Fig. 2
figure 2

Two-stage shear coupling model

Full size image

Compared with Fig. 1, it shows that for the two-stage shear coupling model, the linear elastic behaviour of the interface between the bolt and the grout is equal to the original constitutive law of the cableSELs. However, after δ0, the original constitutive law of the cableSELs assumed that the shear stress of the interface between the bolt and the grout was constant. Compared with that, in the two-stage shear coupling model, the interface shear stress between the rock bolt and the grout material reduced. With this method, the post-failure behaviour of the interface between the rock bolt and the grout material can be simulated.

3 Incorporating the shear coupling model into FLAC3D

FLAC3D is the numerical software developed by the Itasca Company. Due to its convenience and high efficiency, FLAC3D has been widely used in solving geotechnical problems. In the software, structural elements, namely the cable, the pile, the shell, are created to simulate the support and reinforcement that are applied on soils and rock masses. Specifically, the cableSELs are always used to simulate the rock reinforcement tendons, such as rock bolts and cable bolts (Itasca 2016). In this study, the FLAC3D 5.01 was used.

The cableSELs are straight finite elements with two nodes. For each node, there is one axially oriented translational degree-of-freedom. Most importantly, the spring-slider is designed to simulate the bonding performance of the interface between the cable and the mesh. The maximum shear force that the grout annulus can bear is composed of two parts: the chemical adhesion and the friction. Moreover, the maximum shear force per unit length of the cable/mesh interface can be calculated with Eq. (3) (Itasca 2016):

$$\frac{{F_{s}^{\max } }}{L} = c_{\text{g}} + \sigma_{\text{m}} p_{\text{g}} \tan \varphi_{\text{g}}$$
(3)

where, Fmax s is the maximum shear force of the cable/mesh interface (N); L is the embedment length (m); cg is the grout cohesive force per unit length (N/m); σm is the effective confining stress (Pa); pg is the grout exposed perimeter (m); and φg is the grout friction angle (°).

The cableSELs are composed of segments and nodes. In this study, it was assumed that all nodes in the cableSELs have the same properties. Therefore, for each node, Eq. (4) can be acquired:

$$\frac{{F_{si}^{\max } }}{{l_{i} }} = c_{\text{g}} + \sigma_{\text{m}} p_{\text{g}} \tan \varphi_{\text{g}}$$
(4)

where, Fsimax is the maximum shear force of the cable/mesh interface at node i (N); li is the corresponding segment length for node i (m).

By multiplying (πD)−1, the following equation was acquired:

$$\frac{{F_{si}^{\max } }}{{\pi Dl_{i} }} = \frac{{c_{\text{g}} + \sigma_{\text{m}} p_{\text{g}} \tan \varphi_{\text{g}} }}{\pi D}$$
(5)

Since the full length of the cableSELs is divided into a number of segments, the shear stress along each segment can be assumed uniform (Benmokrane et al. 1995). Then, Eq. (5) can be rewritten as:

$$\tau_{i} = \frac{{c_{\text{g}} + \sigma_{\text{m}} p_{\text{g}} \tan \varphi_{\text{g}} }}{\pi D}$$
(6)

where, τi is the shear stress at node i (Pa).

In this paper, the shear coupling behaviour of the interface between the bolt and the grout is independent of the grout friction angle. Therefore, substituting φg = 0 into Eq. (6) leads to the following equation (Ma et al. 2014):

$$\tau_{i} = \frac{{c_{\text{g}} }}{\pi D}$$
(7)

By coupling Eqs. (1), (2) and (7) together, the following equation was acquired:

$$c_{\text{g}} = \frac{{\tau_{\max } \pi D}}{{\delta_{0} }}\delta \, \left( {\delta < \delta_{0} } \right)$$
(8)
$$c_{\text{g}} = \frac{{\tau_{\max } \pi D}}{{e^{{c\left( {\delta - \delta_{0} } \right)}} }} \, \left( {\delta \ge \delta_{0} } \right)$$
(9)

Since cg in Eq. (8) is the input parameter for the cableSELs, the two-stage shear coupling model can be incorporated into cableSELs. Consequently, the current constitutive law of the cableSELs can be modified with the two-stage shear coupling model.

To technically realise this, a FISH function was created. A flow chart is used to depict the logic of the created FISH function, as shown in Fig. 3.

Fig. 3
figure 3

Logic of the created FISH function

Full size image

After this FISH function is executed, the head pointer of the cableSELs is extracted and assigned to the address variable of a. Then, logical judgement is conducted to check whether the address variable of a is null. If the address variable of a is not null, the displacement in the grout at the address variable of a under the large strain calculating mode is extracted and assigned to the shear slippage (δ). Then, δ is compared with δ0. If δ is smaller than δ0, δ is substituted to Eq. (8) to calculate the corresponding shear stress. Otherwise, δ is substituted to Eq. (9) to calculate the corresponding shear stress. Following that, the calculated shear stress is assigned to the grout cohesive force per unit length at the address variable of a. Last, the next address of the structural element is acquired and assigned to the address variable of a. This FISH function ends when the address variable of a is null.

To guarantee that the shear stress of the interface between the bolt and the grout can be varied automatically with the shear slippage, this FISH function is executed in every timestep in the calculation process in FLAC3D.

4 Validation of this numerical model

To validate the credibility of this numerical approach, a numerical pulling test was conducted. To further confirm whether the numerical pulling results were appropriate, experimental pulling results were introduced for comparison and validation.

4.1 Pulling scenario

Specifically, axial loading tests on fully encapsulated rock bolts were performed by Liu et al. (2011). The tested tendon had a diameter of 42 mm. In the testing process, those rock bolts were encapsulated in boreholes that have a diameter of 150 mm. The embedment length is 3000 mm. Cement-based grout with a strength of 27.1 MPa and elastic modulus of 26 GPa was poured into the borehole to bond the rock bolt with surrounding rock. To measure the axial load variation along the tendon, strain gauges were adhered on the tendon surface. The distance between strain gages was 150 mm. Additionally, the force and displacement of the rock bolt external end were also measured.

In this paper, this pulling process was simulated with FLAC3D. Specifically, a numerical block with a length of 3000 mm was established. The length, width and height of the numerical model were 3 m, 0.4 m and 0.4 m. Along the axes of X, Y and Z, the number of zones was 10, 10 and 30. The whole numerical model was composed of 3000 zones and 3751 grid points. In the centre, the rock bolt was simulated with the installed cableSELs. The cableSELs are composed of 30 segments and 31 nodes. The endpoint of the rock bolt was located at the grid point. However, the endpoint of the rock bolt was coupled to the grid point via the grout, following the shear behaviour of the bolt/grout interface. Therefore, the endpoint of the rock bolt was not slaved on the grid point. Under the axial loading, the endpoint of the rock bolt can still move along the axial direction of the rock bolt.

The geometry of the numerical block and the cableSELs is shown in Fig. 4. In this figure, the blue line represents the segment of the cableSELs and the green sphere represents the nodes of cableSELs.

Fig. 4
figure 4

Numerical pulling scenario for the rock bolt

Full size image

4.2 Input parameters in the numerical model

The input parameters for the shear coupling model include the shear strength of the interface between the bolt and the grout (τmax), the shear slippage when the shear strength of the interface between the bolt and the grout is reached (δ0) and the coefficient for the plastic stage (c). These three parameters determined the bond-slip behaviour of the bolt/grout interface and can be acquired from conducting pulling tests on short encapsulated rock bolts. Specifically, Benmokrane et al. (1995) mentioned that once the rock bolt encapsulation was short, the interface shear stress can be regarded uniform. Therefore, through conducting pulling tests on short encapsulated rock bolts, the bond-slip behaviour of the bolt/grout interface can be acquired.

For this scenario, Liu et al. (2011) conducted short encapsulation pulling tests on rock bolts and acquired the bond-slip behaviour of the bolt/grout interface. Based on their results, the input parameters tabulated in Table 1 were used in this study.

Table 1 Input parameters for the shear coupling model
Full size table

It is admitted that the characteristics of the interface between the bolt and grout is dependent on the bolt length for fully grouted bolts. This may be more apparent if the grouted length was long enough. If the grouted length was long enough, the grouting quality at different locations along the rock bolt may be different.

However, in this study, the purpose was to modify the original structural elements of cable in FLAC3D with the two-stage shear coupling model. And in this study, a relatively short length of 3 m was used. Therefore, it was assumed that in this study, the characteristics of the interface between the bolt and grout was uniform along the bolt length.

As for the diameter and cross-sectional area of the rock bolt, they can be acquired from the corresponding experimental pulling scenario. The Young’s modulus of the rock bolt, can be referred from the handbook or acquired from the rock bolt manufacturer. The grout shear stiffness can be calculated from the input shear coupling model, following Eq. (10):

$$k_{\text{g}} = \frac{{\tau_{\max } \pi D}}{{\delta_{0} }}$$
(10)

where, kg is the grout shearing stiffness (N/m/m). Therefore, once the input shear coupling model was determined. The shear stiffness of the grouting material can be computed. For instance, in this scenario, substituting τmax and δ0 in Table 1 into Eq. (10), leaded to the result that the grout shear stiffness was 2.375 × 109 N/m/m.

Then, there was only one input parameter left, which was the number of segments that constituted a rock bolt. In this paper, for each rock bolt, the number of segments was set as 30. Later, a parametric study will be conducted to demonstrate the influence of segments on the performance of rock bolt. Input parameters of cableSELs are summarised in Table 2.

Table 2 Input parameters for the cableSELs
Full size table

4.3 Comparison with experimental results

In this numerical simulation, the Mohr–Coulomb model was used as the constitutive law for the surrounding rock. To simulate the pulling process of the rock bolt, an initial pulling velocity of 1 × 10–6 m per step was applied at the loaded end of the cableSELs, as shown in Fig. 5. This pulling velocity kept constant during the whole calculation process.

Fig. 5
figure 5

Applying a pulling velocity on the loaded end of the cableSELs

Full size image

The out-of-balance force along the pulling direction of the loaded end was recorded and regarded as the pulling load. Also, the displacement along the pulling direction of the loaded end was recorded and regarded as the pulling displacement.

After those properties were set, the numerical pulling test was initiated. During this numerical pulling process, the two-stage shear coupling model was incorporated into cableSELs to revise the current constitutive law. The relationship between the pulling load and displacement was acquired and compared with experimental pulling results, as shown in Fig. 6. It is apparent that there was a favourable agreement between results. This confirms that the modified cableSELs can simulate the load transfer behaviour of fully encapsulated rock bolts.

Fig. 6
figure 6

Comparison between the numerical pulling result and the experimental pulling result acquired by Liu et al. (2011)

Full size image

A numerical pulling test with the original cableSELs was also conducted and the result is plotted in Fig. 6. It can be seen that there was a large difference between the numerical simulation results acquired with the original cableSELs and the experimental results. In fact, the original cableSELs overestimated the load bearing capacity of rock bolts.

To further examine whether this two-stage shear coupling model has been incorporated into cableSELs, the shear stress and shear slippage at nodes along the embedment length were monitored. Then, the relationship between the shear stress and shear slippage at nodes can be acquired. Specifically, the results acquired from Nodes 2, 4 and 6 were taken as an example. Node 2 represents the node that was adjacent to the loaded end while Node 6 represents the node that was relatively far away from the loaded end. As for Node 4, it was located between Nodes 2 and 6.

Attention should be paid that in the cableSELs in FLAC3D, the shear force per unit length of the nodes is used, and the unit for it is N/m. To create the relationship between the shear force per unit length of the nodes, and the shear stress, Eq. (11) was used (Ma et al. 2014).

$$\tau_{i} = \frac{{f_{si} }}{\pi D}$$
(11)

where, fsi is the shear force per unit length at node i, N/m.

It can be seen that Eq. (11) was similar to Eq. (7). In fact, Eq. (7) depicted the relationship between the shear stress at the bolt/grout interface and the grout cohesive force per unit length in the input stage. And Eq. (11) depicted the relationship between the shear stress at the bolt/grout interface and the exported shear force per unit length at nodes in the output stage. Therefore, the outputted shear force per unit length at nodes basically equalled the inputted grout cohesive force per unit length.

Figure 7 shows the relationship between the shear stress and shear slippage at nodes along the modified cableSELs. It can be seen that the shear stress versus shear slippage relationship of those three nodes was consistent. Specifically, when the shear slippage of the node was smaller than 0.15 mm, the shear stress generally increased linearly to the peak. After it, the interface shear stress started reducing with an exponential trend.

Fig. 7
figure 7

The relationship between shear stress and shear slippage at nodes in the modified cableSELs

Full size image

In the numerical pulling test conducted with the original cableSELs, the shear stress versus shear slippage relationship of those three nodes was also recorded, as shown in Fig. 8. In the original cableSELs, after the shear stress of nodes reached the shear strength, it remained constant. There was no decreasing of the shear stress for the bolt/grout interface, which was not a true reflection of the reality. This is also the reason why the original cableSELs overestimated the load bearing capacity of rock bolts.

Fig. 8
figure 8

The relationship between shear stress and shear slippage at nodes in the original cableSELs

Full size image

Then, it can be concluded that the original constitutive law of cableSELs has been successfully modified with the two-stage shear coupling model. More importantly, the modified cableSELs can simulate the post-failure behaviour of the interface between the rock bolt and the grout material.

It was also found that in the modified cableSELs, the shear stress state of those three nodes were different. Specifically, in Fig. 7, when the shear stress at Node 2 decreased to 0.48 MPa, the shear stress at Node 4 was 0.7 MPa and the shear stress at Node 6 was 1 MPa. The results also revealed that at the same time, the shear stress state at different nodes was various even though the constitutive law of those nodes were same.

To further confirm this, the relationship between the shear stress of those three nodes and the pulling displacement is plotted, as shown in Fig. 9. It shows that with the pulling displacement increasing, the shear stress at Node 2 first reached the shear strength of the interface between the bolt and the grout. At the same time, for Nodes 4 and 6, the shear stress was still rising. With the pulling displacement increasing, the shear stress of the Node 2 decreased while the shear stress at Node 4 was going to reach the shear strength of the interface between the bolt and the grout. When the pulling displacement increased to 0.57 mm, the shear stress at Node 6 finally reached the shear strength of the interface between the bolt and the grout. However, meanwhile, the shear stress at Node 2 and Node 4 already decreased to a much lower value.

Fig. 9
figure 9

The relationship between the shear stress at nodes and the pulling displacement

Full size image

To further confirm the credibility of this numerical model, the load distribution along the rock bolt was recorded and compared with experimental results. The comparison between the load distribution results acquired from the numerical pulling test and the load distribution results acquired from the experimental pulling test when the pulling load was 160 kN, 250 kN and 300 kN, is shown in Fig. 10. The results further confirmed that the modified numerical method can correctly simulate the load distribution state along the fully encapsulated rock bolt. This further validates the credibility of the modified cableSELs, confirming that the modified cableSELs can simulate the performance of fully encapsulated rock bolts properly.

Fig. 10
figure 10

Comparison between the numerical load distribution results and the experimental load distribution results acquired by Liu et al. (2011)

Full size image

4.4 Parametric study on the number of segments

In the input parameters for the cableSELs, there was an empirical parameter, namely the number of segments that constitute a rock bolt. In Sect. 4.2, this parameter was set as 30. Then, a numerical investigation was performed, evaluating the influence of this parameter on the rock bolt performance.

Specifically, the pulling scenario was still same as the Sect. 4.1. However, the number of segments was different, varying from 10 to 50. Then, the axial loading trend of rock bolts was computed. The results were plotted in Fig. 11.

Fig. 11
figure 11

Load–displacement performance of rock bolts when the number of segments was different

Full size image

The results showed that although the number of segments was different, the load–displacement trend of the rock bolts was consistent. This indicated that the number of segments had marginal effect in deciding the load–displacement trend of the rock bolt. With the number of segments increasing, the maximum load bearing capacity of the rock bolt increased slightly. Moreover, once the number of segments was larger than a limit, the maximum load bearing capacity of the rock bolt was almost constant. For instance, in this case, when the number of segments increased from 30 to 50, the maximum load bearing capacity of rock bolts was almost equal.

This parametric study demonstrated that the number of segments had marginal effect in deciding the loading performance of rock bolts. Moreover, once this parameter was beyond a limit, its effect on the loading performance of rock bolts can be neglected. Therefore, when using the modified cableSELs, once the other parameters are acquired, the number of segments can be calibrated from the experimental load–displacement curve.

5 A study with the modified cableSELs

After the modified cableSELs were validated with experimental results, a study was conducted to evaluate the performance of rock bolts.

In fact, the performance of rock bolts was influenced by a number of parameters. For instance, the author of this paper conducted experimental tests and numerical simulation to study the influence of rock mass strength on the performance of fully grouted bolts in the previous research (Chen et al. 2016a, 2018). It was found that the surrounding rock mass strength had an important role in determining the performance of bolts. With the strength of the surrounding rock mass increasing, the shear stress in the grout increased apparently. Consequently, this leaded to increasing of the load bearing capacity of bolts. For this paper, as the preliminary study, the influence of the bolt diameter on the performance of rock bolts was first studied.

5.1 Effect of the bolt diameter

The geometry of the numerical model was consistent with Fig. 4 except that the length of the model was 1.5 m. Specifically, the length, width and height of the numerical model were 1.5 m, 0.4 m and 0.4 m. Along the axes of X, Y and Z, the number of zones was 10, 10 and 30. The whole numerical model was composed of 3000 zones and 3751 grid points.

The bolt diameter changed from 15 to 35 mm with an interval of 10 mm. An embedment length of 1.5 m was used. For the interface between the rock bolt and the grout material, input parameters including τmax, δ0 and c, are summarised in Table 3.

Table 3 Input parameters for the bolt diameter effect test
Full size table

Following this, the numerical axial loading test was performed. And the load–displacement relationship of rock bolts with different bolt diameters was acquired, as shown in Fig. 12.

Fig. 12
figure 12

Rock bolt performance when the diameter was different

Full size image

It is apparent that the rock bolt performance was largely dependent on the rock bolt diameter. Under the same numerical pulling scenario, when the rock bolt diameter was 15 mm, the rock bolt had a maximum load bearing capacity of 40 kN and the corresponding pulling displacement was 1.16 mm, showing apparent ductile behaviour. However, when rock bolt diameter increased to 35 mm, the rock bolt reached a maximum load bearing capacity of 150 kN and the corresponding pulling displacement decreased to 0.74 mm, showing brittle behaviour. The results revealed that there was a positive relationship between the rock bolt diameter and the load bearing capacity. Furthermore, the rock bolt was more likely to show the brittle behaviour.

Axial displacement distribution along the rock bolt was studied. Specifically, the rock bolt with a diameter of 25 mm was taken as an example. During the numerical pulling process, the displacement of three different nodes, namely Node 2, Node 4 and Node 6, were monitored. Node 2 was relatively adjacent to the loaded end while Node 6 was relatively far away from the loaded end. Moreover, Fig. 13 shows the axial displacement of those three nodes versus the pulling displacement. It can be seen that with the pulling displacement increasing, the displacement of the nodes along the cableSELs was different. Furthermore, the displacement at node that was adjacent to the loaded end was always larger than the displacement at node that was far away from the loaded end.

Fig. 13
figure 13

Displacement variation at nodes in the cableSELs with the pulling displacement increasing

Full size image

To further confirm this, at the end of this numerical pulling test, the displacement distribution of the rock bolt is plotted, as shown in Fig. 14. It can be seen that the rock bolt displacement distributed non-uniformly. Specifically, the displacement of the rock bolt decreased non-linearly from the loaded end to the other end. The results revealed that when the full encapsulated rock bolt was loaded, along the axial direction of the rock bolt, the displacement of the tendon was various.

Fig. 14
figure 14

Displacement distribution along the embedment length

Full size image

With the modified cableSELs, the shear stress distribution in the grout can also be plotted along the rock bolt length. Specifically, in FLAC3D, the rock bolt was composed of segments and nodes. Then, the bolt tendon interacted with the grout through nodes. The shear stress in the grout was stored in the nodes and it can be exported. For instance, for the rock bolt with a diameter of 25 mm, when the pull-out displacement arrived at 0.3 mm, the shear stress distribution in the grout can be exported, as shown in Fig. 15.

Fig. 15
figure 15

Shear stress distribution in the grout material along the bolt

Full size image

Considering that the displacement of the rock bolt along the embedment length had a significant effect in determining the shear stress of the interface between the bolt and the grout, the shear stress distribution at the interface between the bolt and the grout along the embedment length was analysed. For the numerical pulling test that was conducted on the rock bolt with a diameter of 25 mm, the performance of the rock bolt when the pulling displacement was 0.3 mm, 0.6 mm and 0.9 mm was recorded. The corresponding load was 60 kN, 77 kN and 83 kN respectively, as shown in Fig. 16.

Fig. 16
figure 16

Selecting three points to record and analyse the performance of the rock bolt

Full size image

At these three selected points, the shear stress distribution along the interface between the bolt and the grout is plotted, as shown in Fig. 17.

Fig. 17
figure 17

Distribution of the interface shear stress: a Pulling load of 60 kN; b Pulling load of 77 kN; c Pulling load of 83 kN

Full size image

It can be seen that in the loading process of the rock bolt, the shear stress was not uniform along the full embedment length. When the pulling load was 60 kN, the maximum shear stress occurred at the rock bolt that was adjacent to the loaded end. Then, when the pulling load arrived at 77 kN, the maximum shear stress shifted to the central section of the encapsulated length. When the pulling load increased to 83 kN, the maximum shear stress propagated to the section that was far away from the loaded end.

The results also revealed that when the pulling load was different, there was a slight difference between the maximum interface shear stress. For instance, with the pulling load rising from 60 to 83 kN, the maximum interface shear stresses were 7.6323 × 104 N/m/m, 7.6152 × 104 N/m/m and 7.7022 × 104 N/m/m. The relative differences between the latter two values and the first value were  − 0.22% and 0.92%. It can be seen that there was a tiny relative difference between them. And this tiny relative difference was probably due to the calculating method of FLAC3D. The calculated results were the approximate solution, which was close to the true solution. However, there may be a difference between the approximate solution and the true solution.

5.2 Effect of the grout shear strength

Then, the effect of grout shear strength on the tensile performance of bolts was studied. The numerical pull-out scenario was basically consistent with that in Sect. 5.1, where a rock bolt with a diameter of 25 mm was fully embedded in the numerical block.

During calculating, the grout shear strength changed from 0.5 to 1.5 MPa. The numerical pull-out tests were then performed. Meanwhile, the loading behaviour of the bolt was monitored. The calculated results are shown in Fig. 18. Apparently, the grout shear strength influenced the rock bolting performance. Larger shear strength of the grouting material leaded to higher peak load of rock bolts. For instance, as the shear strength of the grouting material was 0.5 MPa, the peak force was only 67 kN. However, as the shear strength was 1.5 MPa, the peak force increased to 107 kN. Moreover, with the grout shear strength rising, the displacement where peak force occurred also increased. Specifically, with the grout shear strength increasing from 0.5 to 1.5 MPa, the axial displacement where peak load occurred increased from 0.64 to 1.16 mm. Therefore, improving the grout shear strength leaded to better loading performance of rock bolts.

Fig. 18
figure 18

Performance of bolts when grout shear strength was different

Full size image

6 Application and future work

In this paper, a two-stage shear coupling model was incorporated into the current structural elements in FLAC3D to simulate the bond failure process of the bolt/grout interface. The incorporation was realised with FISH function. Since structural elements are also used in other software, such as UDEC and 3DEC, this two-stage shear coupling model can also be incorporated into the other software.

The practical application of the research result is that since the modified cableSELs have been successfully validated with experimental results, they can be further used to study the interaction between rock bolts and rock mass under the in-situ condition. This is also the future work of the authors in this study. In the future work, more investigation will be conducted to use the modified cableSELs to reinforce roadways and excavations under the in-situ condition. The performance of rock bolts in keeping the stability of underground roadways and excavations will be further studied.

7 Conclusions

The original constitutive law of the cableSELs in FLAC3D assumes that the shear behaviour of the bolt/mesh interface follows an elastic perfectly-plastic model. Once yielding, the interface shear stress did not decline. To overcome this shortcoming, a two-stage shear coupling model was used to simulate the bonding performance of the interface between the rock bolt and the grout material. First, the interface shear stress rose with slippage. At a certain slippage, it reached the strength. Then, in the second stage, there was an exponential decreasing relationship between the interface shear stress and the interface shear slippage. With this method, the post-failure behaviour of the bolt/grout interface was modelled.

This two-stage shear coupling model was incorporated into FLAC3D to revise the current cableSELs. A numerical pulling test was conducted to evaluate the load transfer performance of fully encapsulated rock bolts with the modified cableSELs. Experimental results were adopted to validate the accuracy of the numerical model. The results show that there was a close match between numerical and experimental results. Therefore, the credibility of the modified cableSELs was validated.

Based on the modified cableSELs, the influence of the bolt diameter on the performance of rock bolts was studied. Based on numerical modelling, it was found that the load bearing capacity of rock bolts increased with diameter. Furthermore, under the same pulling scenario, with the bolt diameter increasing, the performance of the rock bolting system was likely to change from the ductile behaviour to the brittle behaviour. During the pulling process, the displacement of the rock bolt was non-uniform along the embedment length. More importantly, from the rock bolt loaded end to the other end of the rock bolt, the axial displacement exponentially reduced. Additionally, with the pulling displacement increasing, the maximum shear stress of the interface between the bolt and the grout mobilised, shifting from the rock bolt loaded end to the other end of the rock bolt.