An Improved Numerical Simulation Approach for the Failure of Rock Bolts Subjected to Tensile Load in Deep Roadway

Our goal was to develop an effective research tool for roadways with significant deformations supported by rock bolts. The improved numerical simulation approach is constructed through additional development of FLAC3D. The aim is to modify the shortcoming that the original model in FLAC3D regards the plastic tensile strain of any arbitrary rock bolt element node as the rupture discrimination criterion. The FISH programming language is adopted to conduct the secondary development and to embed the revised model into the main program of FLAC3D. Taking an actual mining roadway as the simulation object, two simulation schemes adopting the newly improved approach and the original method were conducted, respectively. The results show that (1) the PILE element that constitutes the rock bolt-free section with the maximum elongation rate ruptures after modification, while the rock bolt tendon elongation rate reaches beyond the predefined tensile rupture elongation rate; (2) the modified model in which the rock bolt is mainly subjected to tension realises the tensile rupture phenomenon at the end of the rock bolt-free section and the rock bolt at the junction between the free section and the anchoring section; and (3) only four rock bolts that are in the roadway sides showed rupture in the modified model, and all rock bolts showed rupture in the original model. The tensile failure of the rock bolt led that the modified model scheme is closer to the actual. Compared with the modified model, in the original model, deformation of the surrounding rock masses is severe. This is resulted by the rupture of all rock bolts in the original model. The analysis shows that the improved numerical simulation approach is much more reliable for large deformation roadway behavior with rock bolt support.


Introduction
Rock bolts are widely used in civil and mining engineering to stabilize underground excavations [1][2][3][4][5]. They are installed to reinforce fractured rock mass by resisting dilation or shear movement along the fractures [6][7][8]. Several rock bolts can create a reinforced zone to improve the self-supporting capacity of the rock mass [9,10]. Currently, mining depths for both metal ore and coal have increased beyond 1000 m throughout the world, including mines in Sweden, Canada, Australia, and South Africa [11]. Certain metal mines even reached depths greater than 3000 m [12,13]. By early 2015, nearly 60 coal mines in China had mining depths that exceeded 1000 m (including 7 coal mines deeper than 1200 m); deep coal mining is now imperative [14,15]. Compared with shallow coal roadways, the stability of roadways in deep mines is difficult to control, because of the larger deformation and the larger mining-induced stress. This made it difficult to install an effective support system [16,17]. Although high-strength bolt (and cables) supports were used, under the effect of soft strata, high ground stress, underground water, mining, and even rock burst [18][19][20][21], the tensile-failure problems of support structures are still very common. Some typical cases are shown in Figure 1 (e.g., in the Xiaojiawa mining area in Shanxi, China).
The tensile-failure problems of rock bolts are complicated. Furthermore, the theoretical analysis method is difficult to use in solving this issue, and the model test method cannot be widely used because of its high cost [22,23]. This requires new studies to properly consider the failure in tension of those rock bolts caused by excessive deformation (elongation). Currently, numerical simulation is an important and widely used way to study rock bolt failure mechanism and to improve support technology [24]. However, little effort has been spent on the support structures themselves, such as the arch supports, rock bolts, or anchor cables [25][26][27][28][29]. Consequently, there are still obvious deficiencies in the simulation of the mechanical behavior of the support structures; these deficiencies are discussed in the following: (1) Nowadays, based on the basic characters of the elongation of rock bolt subjected to tensile loading to the final rupture, the rock bolt elements usually cannot fully show the character of the tensile rupture when numerical simulation is conducted on the rock bolt reinforcement in soils and rock masses (2) In FLAC3D, due to the influence of the algorithm, the imbalanced force exists in the calculating process of each step. This leads to the consequence that the difference between the imbalanced forces of each element nodes in the rock bolt-free section is relatively larger. For example, in the practical anchoring engineering, the rock bolt that is adjacent to the bearing plate and the element node in the junction between the free section and the anchoring section are sub-jected to relatively larger imbalanced force. As for the rest of the element nodes, they are subjected to relatively small imbalanced force. This leads to the result that the difference on the deformation of each structural element in the rock bolt-free section is relatively large. If the plastic tensile strain of a structural element that has the maximum deformation firstly arrives at the tensile rupture strain, it will then show tensile rupture. However, meanwhile, the plastic tensile strain of the rest of the structural elements in the rock bolt-free section is all smaller than the tensile rupture strain. Therefore, from the whole perspective of the rock bolt-free section, this rock bolt shows the tensile rupture prematurely. This is apparently contrary to the practical situation, making the calculating accuracy and precision decrease. Consequently, this restricts the application in several aspects In the above research, when the rock bolt element in the numerical simulation is applied, cableSELs were incorporated proverbially. However, few researches were conducted on the PILE element. Additionally, the tensile rupture strain of the PILE element is set as the rupture elongation rate of the rock bolts. When the local section of the rock bolt tendon shows the tensile rupture, the whole elongation rate of the rock bolt-free section is smaller than the predefined rupture elongation rate. The rock bolt will rupture prematurely. The simulation precision and accuracy decrease largely.
Therefore, the research work in this paper is mainly based on the second development platform of FLAC3D. The previously mentioned problems were solved through modification of the pileSELs, respectively. Moreover, an improved numerical simulation approach was then established, which took the modifications as the core. Finally, the 2 Geofluids approach was verified through a simulation analysis of a typical real engineering case.

SELs Used to Simulate Rock Bolts in FLAC3D
FLAC3D is the numerical calculating software that is applicable for the large deformation of the geotechnical engineering [30]. In FLAC3D, there are primarily two types of SELs that can be used to model a rock bolt, i.e., cable element and pile element [30,31]. Those two SEL types describe an abstract interaction between rock zones and bolts instead of physical analysis. CableSELs integrated in FLAC3D can be imagined as a long steel bar, which can only carry the axial tensile and axial compressive loads [2,31]. Except the axial tensile and compressive loads, pileSELs can carry a bending moment induced from rock movement normal to the bolt axis. Thus, the pileSEL can replicate the real behavior of a bolt subjected to combined shear and axial loads, which is more common in field conditions [31][32][33]. It is the pileSEL, which can relatively accurately simulate the tensile loading along the longitudinal direction, lateral shearing, and the anchoring effect between rock bolts and surrounding rock masses. Meanwhile, it can realise the tensile rupture function of rock bolts through defining the tensile rupture strain. When the pileSEL realises the tensile rupture function of rock bolts, its rupture function is based on the plastic tensile strain of any arbitrary element nodes. It means that when the plastic tensile strain of elements is larger than or equal to the tensile rupture strain, this element shows rupture [30].

The Current pileSEL Model and Problems in FLAC3D
2.1.1. The Basic Components and Theory of the pileSEL. The single pileSEL is composed of two structural nodes and a structural component [30]. Between the nodes and the entity elements, there exists the coupled spring-slider model along the shear direction and normal direction. This can relatively better reflect the transfer of the load and moment between those two. Therefore, it can better simulate the mutual effect between it and the entity zones that are connected with it along the normal direction and the axial direction. Therefore, the pileSEL can better reflect the phenomenon of the tension along the axial direction and normal direction. Its mechanical model is shown in Figure 2. It shows that the pileSEL and the entity zones transfer the load and bending moment through the spring-slider model along the shearing direction and normal direction. Figure 2 mainly involves the geometric and mechanical parameters including the radius of the rock bolt tendon, the normal stiffness, the longitudinal stiffness, the adhesive force of the grout along the normal direction, the friction angle, the stiffness, the adhesive force of the grout along the shearing direction, the friction angle, the stiffness, and the drilled borehole diameter. The nodes of the pileSEL are used to store the data such as the displacement and the unbalanced force of grid points [30]. The line between two structural nodes is called as the structural component, which is used to store the data such as the axial force, shear force, bending moment, and length of the structural elements. To conveniently describe the force and bending moment in the structural component, each structural element should have its local coordinate system, as shown in Figure 3. The axis of x is the longitudinal direction, pointing from node 1 to node 2. The axes of y and z are shearing that are perpendicular to the longitudinal direction. The rectangular coordinate system is composed of xyz, yielding the right-hand rule. Each structural node has 6 degrees of freedom: three translational degrees of freedom and three rotational degrees of freedom. Each structural element has a unique number (CID), which is used to assign properties and locate structural elements. Through combining multiple elements, a full-length rock bolt can be simulated. Additionally, Figure 3 shows the rules for the positive direction of the force and bending moment of the structural elements in the local coordinate system. Specifically, F x indicates the axial force. F y and F z indicate the shear force. M x , M y , and M z indicate the bending moment [30].
The reinforcing function of the pileSELs is that the difference on the tangential displacement of two ends of the pileSELs leads to the longitudinal deformation of the rock bolt. Consequently, axial force is generated in the rock bolt. Furthermore, the grout relied and acted on the surrounding rock masses. Additionally, the difference on the normal displacement of two ends of the pileSEL leads to the bending moment and shear force in the rock bolts. Meantime, the grout is acted on to act on the surrounding rock masses. Finally, the balance is realised.

Longitudinal Mechanical
Model of the pileSEL Tendon and Realising of the Tensile Rupture. The material of the rock bolt tendon is usually steel. Tension and compression along the longitudinal direction can be regarded as the elastic perfectly plastic material. Therefore, when the pileSEL is used to simulate rock bolts, for the constitutive law of the longitudinal tension of the rock bolt tendon, the elastic perfectly plastic model is used, as shown in Figure 4. Its mechanical expression is shown below: where F a is the axial force in the rock bolt; F t max is the tensile strength of the rock bolt; A is the cross area of the rock bolt; E is the elastic modulus; ε is the elastic strain; ε pl is the plastic tensile strain; tfstrain is the tensile failure strain. The elastic perfectly plastic mechanical model of the pile-SEL tendon can be classified as two stages, namely, (1) elastic stage (section OA) In the elastic stage, there is a linear proportional relationship between the axial force and deformation. The slope is the axial stiffness K of the pileSEL, namely, where K is the axial stiffness of the rock bolt tendon; A is the cross-sectional area of the rock bolt tendon; E is the elastic modulus of the rock bolt tendon material; and L is the length of the structural component. For each parameter, the international standard unit is used.
(2) plastic stage (after the point A) When the axial force of the pileSEL reaches the tensile yield force F t max , the rock bolt tendon enters the yield state. The axial force will not rise again and the longitudinal stiff-ness K decreases to 0. However, at this time, the rock bolt tendon element cannot elongate infinitely. When the tendon elongation reaches the ultimate value, the tendon shows the necking phenomenon. Then, the tensile rupture will occur. Meantime, the axial force in the rock bolt tendon decreases to 0. Attention should be paid that when the rock bolt tendon is subjected to compression state, the rock bolt tendon enters the yielding state once the axial force in the pileSEL reaches the compressive yielding force F c max . The axial force will not increase again, and the axial stiffness K decreases to 0. The rock bolt element can be stretched and compressed infinitely. When the pileSEL reaches the yielding strength, decreasing of the external load makes the axial force in the element show the decreasing tendency. The loading path is shown with the dash arrows in Figure 4. It should be mentioned that in FLAC3D, it is not necessary to compulsorily define the yielding parameters. The above contents describe the tension and compression situation of rock bolts. In fact, in the tensile rupture failure mode of rock bolts, the tensile situation is mainly considered.
When the pileSEL is used to simulate the rock bolts, the plastic tensile yielding strain of each structural nodes of the pileSEL (ε pl ) can be recorded. Its calculating expression is shown in the following equation: where ε pl is the total tensile yielding strain of the pileSEL; d is the diameter of the structural element; L is the length of the structural element; θ is the average rotational angle of the structural element; ε ax pl is the axial strain of the structural element. From Equation (3), it can be seen that the total tensile yielding strain is the sum of the axial strain and the bending strain. In the numerical calculating process, the users can define the tensile rupture strain (tfstrain). When the plastic tensile strain of the structural elements is more than or equal to the tensile rupture strain, the force and the bending moment suddenly decrease to 0. At this time, the element fails.    (2) The grout has the mechanical character of tangential strain softening. It indicates that the tangential cohesion and the tangential internal friction angle can be weakened according to the tangential displacement of the grout. The mechanical character is closer to the reality (3) Through defining the tensile rupture strain, the pileSEL can realise the tensile rupture. However, when the tensile rupture strain is set as the rupture elongation rate of the rock bolt, the rock bolt tendon usually shows the premature rupture. When the rupture occurs, the whole elongation rate of the rock bolt-free section is far smaller than the predefined elongation rate, largely decreasing the simulation precision and accuracy 2.1.4. Shortcomings in the pileSEL Model. The pileSEL can realise the tensile rupture function of rock bolts. However, its rupture function is based on the judgement criterion of the plastic tensile strain of any element node. Therefore, it cannot realise the tensile rupture function of the rock bolt tendon from the perspective of the whole elongation rate of the free section. When the rock bolt rupture elongation rate is defined as the tensile rupture strain of the element, the rock bolt usually ruptures prematurely. When the rock bolt rupture occurs, the whole elongation rate of the rock bolt-free section is far smaller than the predefined rupture elongation rate, largely decreasing the simulation precision and accuracy.
3. Modification of the pileSEL and the Realising Process with Programming 3.1. The Current pileSEL Model and Its Problems in FLAC3D.
In engineering practices, when the rock bolt tendon generates large tensile deformation, rupture will occur once it is beyond the ultimate elongation rate. At this time, the axial force in the rupture position suddenly decreases to 0. As for the rest of the positions, the rock bolt will show different mechanical states according to the different situations. Li et al. adopted the cableSEL, in which the elongation rate of each element is regarded as the judgement criterion [29]. When the elongation rate is more than or equal to the predefined allowable elongation value, the cableSEL will show the tensile rupture. This judgement criterion is similar to the criterion that the plastic tensile strain of any element in the rock bolt pileSEL is regarded as the judgement criterion. The rock bolt tendon will prematurely rupture. Aiming at this shortcoming, Yang et al. proposed the judgement criterion that the whole elongation rate of the rock bolt-free section is regarded as the tensile rupture judgement criterion, based on the cableSEL [34]. This criterion is more in accordance with the reality. However, the rupture position is set to occur in the ending section in the rock bolt-free section that is adjacent to the bearing plate. The practical elongation rate of each element in the rock bolt-free section under the critical rupture state is not considered to judge the rupture position of the rock bolt-free section. Overall, the current secondary development process of the rock bolt numerical simulation still has certain shortcomings. Based on the abovementioned shortcomings that occur in the secondary development of the rock bolt elements, this paper incorporates the tensile failure rupture judgement criterion, based on the current pileSEL. It means that the whole elongation rate of the rock bolt-free section is regarded as the judgement criterion. If the elongation rate of the rock boltfree section is more than or equal to the rock bolt rupture elongation rate, modification is conducted on the mechanical model of the structural element that has the largest elongation rate in the free section. The modified mechanical constitutive model of the established pileSEL tendon is shown in Figure 5, being expressed with the following equation: where F ai is the axial force in the structural element that has the largest elongation rate in the rock bolt-free section; δ is the elongation rate of the rock bolt-free section (the percentage between the elongation rate of the rock bolt-free section and the initial length of the free section); and δ max is the tensile rupture elongation rate of the rock bolt.
Combining with the statement in the above context, in reality, when the tensile rupture failure of the rock bolts is considered, the tensile situation should be analyzed. Therefore, in Figure 5, detailed modification is stated for the tensile rupture model of the pileSEL. For the pileSEL shown in Figure 5, the modified mechanical model of the tensile rupture of the rock bolt tendon mainly includes three stages, namely, (1) elastic stage (section OA), which is the same as the original model (2) plastic stage (section AB), which is same as the original model (3) rupture stage (after the point B), in which the axial force in the rock bolt F decreases to zero The axial stiffness K decreases to zero and the increment of the axial force △F decreases to zero.
Compared with the original model, the differences are that when the rock bolt enters the yielding state and the axial elongation reaches the defined rupture judgement criterion, the rock bolt enters the rupture state. The mentioned rupture state is that the axial force of the rock bolt suddenly decreases to zero and becomes constant. Additionally, the axial stiffness K decreases to 0, indicating that the rock bolt can deform arbitrarily along the axial direction while the axial force is constant to 0. Besides the abovementioned differences, there is no modification on the rest of the pileSELs. It should be 5 Geofluids noted that the abovementioned constitutive law uses the length of the rock bolt elements or the variation of the length. In fact, it is equivalent with using the elongation rate to describe. This is because the rupture length is correspondent with the rupture elongation rate.
From Figure 5 and Equation (4), it can be acquired that the modified tensile rupture mechanical model can not only realise regarding the elongation rate of the rock bolt-free section as the judgement criterion for the rupture of the rock bolt tendon. Moreover, it can judge the rupture position of the rock bolt-free section based on the elongation rate of each element in the free section under the critical rupture state. Then, the macro rupture effect from the local to the whole can be further realised.

Realising of the Modified Tensile Rupture Mechanical
Model of the pileSEL with Programming 3.2.1. Theory and Realising Procedures with Programming. The procedures to realise the modified tensile rupture mechanical model of the pileSEL is shown in Figure 6. When the model is calculated to step m, it is taken as an example, to illustrate the modification procedures of the pileSEL. After the model parameters are determined (at this time, input of the parameters for the pileSEL still follows the rule of the current FLAC3D), it enters the main calculating procedure of FLAC3D to perform the calculation.
(1) When calculating to step m, the maximum unbalanced force ratio is regarded as the judgement criterion. First, it is needed to judge whether the procedure has reached convergence. If it is (namely, the maximum unbalanced force ratio is less than 1:0 × 10 −5 ), the calculation ended (2) If the procedure has not reached convergence, it then enters the rupture prejudgement module (i) From the first rock bolt (i = 1), first, summing is conducted on the plastic tensile strain of each structural elements of the rock bolt-free section (the number changes from i 1 to i j ). Judging is conducted to check whether T i equals 0 (through the numerical calculating, it can be acquired that the default value of the plastic tensile strain of the structural elements of rock bolts is 0, meaning that the rock bolt element is in the intact state) (ii) If the condition T i = 0 is true, judging is conducted to calculate the elongation rate of the rock bolt-free section, δ i (iii) If the condition T i = 0 is false, the judgement is moved to the next rock bolt element with the number of (i + 1).
The abovementioned rupture prejudgement module aims at improving the calculating efficiency. It has no influence on the calculating timesteps and results.
(3) The modification module of the pileSEL (i) Judging the elongation rate of the rock bolt i in the free section, δ i (the elongation rate of the rock bolt-free section is the difference between the length of the rock bolt-free section L i and the initial length of the rock bolt-free section L f . Then, the acquired result is divided by L f ) is larger or equal to the rock bolt rupture elongation rate. If the condition is true, the pileSEL that has the largest elongation rate in the free section should be found. The tensile rupture strain of this structural element tfstrainðixÞ is set as a value that is larger than 0. However, it is quite close to 0. In this paper, 1:0 × 10 −10 is used. At this time, the force and bending moment in this element suddenly decrease to 0, realising the tensile rupture function of the rock bolt i. Then, the calculation enters the rock bolt (i + 1). If the condition is false, calculation enters directly the rock bolt (i + 1) (ii) Judging whether the rock bolt number is larger than n (n is the total number of rock bolt elements). If the condition is true, calculation enters the next step. If the condition is false, the rock bolt calculating procedure in the previous step is repeated The core of this procedure is that for the pileSEL, when the elongation rate of the rock bolt-free section is smaller than the predefined rupture elongation rate, the tensile failure strain of the pileSEL changes to the default value of 0. At this time, the tensile failure strain has no restricting effect on the structural elements. However, when the elongation rate of the rock bolt-free section is larger or equal to the predefined rupture elongation rate, the tensile failure strain of the structural element that has the largest elongation rate in the free section is set to the positive value that is extremely sensitive and infinitely close to 0. Then, the tensile rupture of the rock bolts can be realised.  Table 1, the plastic tensile strain and length of each element in the free section in each calculating timestep can be monitored. Through mathematical calculation, prejudgement can be conducted to check whether the rupture occurs and acquire the elongation rate of the rock bolt-free section. When the elongation rate of the free section is larger or equal to the manually defined rock bolt rupture elongation rate, the command line of IF-ENDIF can be used to find the structural element that has the largest elongation rate. With the function of s_cid(sp_pnt), the CID number of this element can be acquired. Then, the command "sel pile property tfstrain" can be used to assign properties for the structural element that has the largest elongation rate in the free section. Then, the sudden tensile rupture of the rock bolt can be realised.  Figure 7. The load-displacement curves for the pile element from the FLAC3D model with and without using the FISH modification were compared with the experimental test, as shown in Figure 8. The numerical model results for the modified pileSEL model in FLAC3D closely match the laboratory test results once the test bolt reaches its maximum tensile strength. The original model formulation instead is not able to represent the actual response of the test rock bolt.

Mechanical Matching
Based on these results, the authors believe that the modified FLAC3D pile element is suitable to model the rock bolt behavior in terms of both load-displacement and axial force distribution. However, important input parameters such as the pile element number and the ultimate extension are determined by trial and error.

Validation of the Numerical Simulation Approach
The improved simulation approach is constructed based on the above modification works, because in the anchorage engineering practice, rock bolts are usually subjected to composite effect including tension, shearing, and bending moment. Under some specific conditions, such as the relatively weak and uniform surrounding rock mass roadways with large deformation, rock bolts are generally installed perpendicular to the roadway surface. The rock bolt tendon is mainly subjected to tension. Under this condition, the rock bolt rupture can be simplified as using the elongation rate of the free The main calculating procedure of FLAC3D Import the parameters of the pileSEL and other models Start Y N Step m

Convergence
End The rupture prejudgement module The modification module of the PILE element   (Figure 9(a)). The working face 110802 is located in the No. 8 coal seam, whose thickness is 3.2 m. The direct roof is mainly composed of sandy mudstone, and the main roof is mainly composed of medium sandstone. A generalized stratigraphic column for the mine is shown in Figure 9(b).
The ventilation roadway in the working face 110802 is chosen as the research object. The layout of the roadway is shown in Figure 9(c). The shape of the roadway is a semicir-cle with straight walls, the excavation radius is 2.5 m, and the height of the wall is 1.8 m. The dimension of the tunnelling cross section of the roadway is 5000 mm × 4300 mm. Concrete with a thickness of 100 mm is used to spray the arch section and the wall section. The net cross-sectional dimension of the roadway is 4800 mm × 4200 mm. The rock bolt reinforcement and the concrete spraying are used to support the roadway. The dimension of the rock bolt is ϕ 22 mm × 2400 mm. The material is the high-strength threaded GFRP rock bolts. At the arch section and the top section of the wall, the rock bolt interval is 800 mm (3-13 rock bolts). At the wall section, the rock bolt interval is 700 mm and the spacing is 800 mm. The mean rock bolt anchoring length is 1000 mm. When the rock bolts in the bottom section are installed, there is an angle between the axial direction of the rock bolt and the horizontal direction, which is 15°. To conveniently describe the rock bolts in the following contexts, from the left corner, along the clockwise direction, the rock bolts are numbered from 1 to 15.
The ventilation roadway is affected by the lateral abutment pressure of the working face 110801 and the advance abutment pressure of the working face 110802, and the deformation of the roadway solid coal rib is more serious. The deformation is described below. From the surface deformation curves of the roadway tunnel, as shown in Figure 10, the deformation speeds are high and the deformations at 120 d are huge. The roof displacement reached 900 mm with an average convergence of 561.8 mm, and the value on the side wall reached almost 500 mm. In the process of mining in this working face, there are a number of broken GFRP anchor rods and local cracking of plastic net. According to the field investigation, the broken fiberglass bolt belongs to tensile breaking type, which is consistent with the low plastic characteristics of fiberglass bolt. The anchor rod still has great strength after reaching the failure load. However, there is no obvious flow plastic stage when the bolt reaches the failure load. The in situ stress in the mining area was measured using borehole relief methods with CSIRO cells. The in situ stress  [39][40][41][42]. In the model, the physical and mechanical parameters of the involved rocks and the sprayed layer are shown in Table 2. The parameters of the strain softening model for the sandy mudstone and coal are shown in Table 3. The material of the SHELL element is isotropic, and the thickness is 10 mm. The elastic modulus is 200 GPa, and the Poisson ratio is 0.3. The ultimate load of the rock bolt is 297 kN, and the rupture elongation rate is 10%. The length of the free section is 1400 mm. Therefore, the ultimate elongation of the free section is 140 mm. The cross-sectional area of the rock bolt simulated with the pileSELs is 3:8 × 10 −4 m 2 . And the elastic modulus is 2:0 × 10 11 Pa. The Poisson ratio is 0.3. The inertia moment along the axis y is 1:15 × 10 −8 m 4 , and the inertia moment along the axis z is 1:15 × 10 −8 m 4 . The polar inertia moment is 2:30 × 10 −8 m 4 . In the modified model, the initial tensile rupture strain of the structural elements of rock bolts maintains its default value. In the original model, the tensile failure strain of the structural elements of rock bolts was set 0.1. The parameters of the grout are shown in Table 4. The excavation of the roadway was finished in one step, and after the excavation, rock bolts are installed immediately. After that, the concrete with a thickness of 100 mm was sprayed. For the rock bolt model, the original model of the pileSEL or the modified model of the pileSEL in FLAC3D is used. The rock bolts are installed uniformly in the middle cross section of the model (0.4 m along the thickness direction). The position, the dimensional parameters, and the modelling effect are shown in Figure 11

Deformation of the Roadway Surrounding Rock Masses
and Analysis on the Rock Bolt Ultimate State. When the software calculation arrives at balance, the maximum unbalanced force ratio is less than 1:0 × 10 −5 . Statistics and analysis are conducted on the calculated results. Figure 12 shows the rock bolt axial load distribution state when the modified model and the original model arrive at the calculation balance. Figure 13 shows the ultimate state of the structural   Figure 14 shows the horizontal displacement contour of the surrounding rock masses when the modified model and the original model arrive at the calculation balance. Figure 15 shows the vertical displacement contour of the surrounding rock masses when the modified model and the original model arrive at the calculation balance. Figure 16 shows the horizontal stress contour of the sprayed layer when the modified model and the original model arrive at the calculation balance. Figure 17 shows the vertical stress contour of the sprayed layer when the modified model and the original model arrive at the calculation balance. Figure 12 shows that in the modified model, 4 rock bolts rupture, namely, rock bolt No. 1, No. 2, No. 14, and No. 15. Furthermore, the ruptured rock bolts are located in the roadway sides. As for the rest of the rock bolts that are not ruptured, the axial load in the free section arrives at 297 kN.
In the original model, all rock bolts rupture. The axial load in the free section is nearly 0. In the anchoring section, the maximum axial load occurs in the rock bolts installed in the roadway side, which is 262 kN. Figure 13 shows that the structural elements that have relatively larger elongation are ruptured elements. It can be known that in the modified model and the original model, rupture of the second and 14 th rock bolt occurs in the junction between the free section and the anchoring section along the free section side. The rupture of the rest of the rock bolts except the first and 15 th rock bolts in the modified model and the second and 14 th rock bolts in the original model occurs in the end of the free section that is adjacent to the plate. Figures 14 and 15 Figures 12-17, it can be known that the horizontal displacement of the surrounding rock masses in the roadway side is relatively larger. At the top, the vertical displacement is relatively larger. Furthermore, compared with the maximum value of the vertical displacement, the maximum value of the horizontal displacement is relatively larger than 100 mm. Deformation of the surrounding rock masses mainly concentrate at the roadway sides. At the bottom of the roadway side, the difference between the surface of the surrounding rock masses and the deep rock layers which correspond with two ends of the rock bolt-free section is relatively apparent. This leads to the consequence that the elongation of the rock bolt-free section is relatively large. In the modified model, only four rock bolts located in the roadway sides rupture, which relatively better reflect the rock bolt rupture position. This is consistent with the engineering practice. However, in the original model, all rock bolts rupture, even though in the position where the difference of the displacement between the surface of the surrounding rock masses and the deep rock layers is not apparent. This is obviously different from the engineering practices, and this cannot truly reflect the rock bolt reinforcing effect. Overall, the results indicate that in the anchoring engineering, using the rock bolt elements in the modified model is more consistent with the engineering practices, compared with the original model.

Comparison and Analysis on the Rupture Failure
Status. Under the calculating condition of the modified model and the original model, rock bolts show the tensile rupture phenomenon. To compare and analyze the difference of the rock bolt tendon rupture process in those two models and validate the engineering applicability and reliability of the modified model, in the calculating process, analysis is conducted on the axial force in the element, the element elongation, and the variation law of the whole elongation of the rock bolt-free section with the calculating timestep, of the dangerous and the less dangerous structural elements including the rock bolt No. 8 at the arch top, the rock bolt No. 14 at the middle of the right roadway side, and the rock bolt No. 15 at the bottom of the right roadway side. Additionally, analysis is conducted on the distribution character of the axial force and elongation in the rock bolt elements when the calculation reaches balance.   Figure 19. Figure 18 shows that under the critical tensile rupture state, the elongation of the structural element (CID = 529) at the end of the rock bolt-free section is maximal, which is 11.4 mm. As for the rest of other positions, the element elongation is all less than 1.4 mm. The total elongation of the rock bolt-free section is 21.3 mm, which is only 15.2% of the ultimate elongation (140 mm). The rock bolt tendon ruptures prematurely. However, when the calculation balance is reached, in the modified model, the elongation of the structural element at the rock bolt-free end (CID = 52) and the structural element (CID = 542) at the junction between the free section and the anchorage section is relatively large, namely, 26.4 mm and 23.2 mm, respectively. However, in the original model, the elongation of the structural element at the end of the free section is maximal, reaching 99.7 mm. As for the rest of the elements, the elongation is less than 1 mm. Figure 19 shows that in the modified model, during the calculation process, the axial force in the structural element (CID = 529) basically maintains at 297 kN. The element elongation increases with the calculating timestep. When the calculating balance is reached, the elongation arrives at 26.4 mm. For the rock bolt No. 8, the ultimate elongation of the free section reaches 68.2 mm, which is only 48.7% of the ultimate elongation. Therefore, the rock bolt does not rupture. Additionally, in the original model, when the calculating timestep reaches 851, the axial force in the structural element (CID = 529) suddenly decreases from 297 kN to 0. The rock bolt ruptures, and after it, the element elongation velocity increases suddenly. The ultimate elongation of the element is 99.7 mm, which is 73.3 mm larger than the modified model. The ultimate elongation of the rock bolt-free section is 104.6 mm, which is 36.4 mm larger than the modified model.
(2) Analysis of the Rupture Process of the Rock Bolt No. 14 at the Roadway Side. In the modified model and the original model, the rock bolt No. 2 and No. 14 which are located in the middle of the roadway side rupture. The rock bolt No. 14 is taken as an example, to analyze the rupture process of the rock bolt located in the middle of the roadway side. The rock bolt No. 14 is composed of 24 pileSELs. The CID number is ranged from 673 to 696, in which the free section is ranged from 673 to 686 and the anchoring section is ranged from 687 to 696. Under the critical tension state and when the calculation reaches balance, in those two different models, the elongation of each element for the rock bolt No. 14 is shown in Figure 20. In the modified model and the original model, under the calculation condition, the structural element of CID = 686 is the dangerous element and the structural element of CID = 673 is the less dangerous element. The axial force and elongation of the dangerous and less dangerous structural element and the variation law between the whole elongation of the rock bolt-free section and the calculating timestep are shown in Figure 21. Figure 20 shows that in the original model, the elongation of the structural element (CID = 686) at the junction between the rock bolt-free section and the anchoring section is maximal, which is 11.3 mm. The elongation of the structural element (CID = 673) at the end of the free section is 8.4 mm. As for the rest positions, the elongation of the elements is all less than 0.8 mm. The total elongation of the rock boltfree section is 26.8 mm, which is only 19.1% of the ultimate elongation or 140 mm. The rock bolt ruptures prematurely. In the modified model, the elongation of the structural element at the junction between the free section and the anchoring section is maximal, which is 78.0 mm. The elongation of the structural element at the end of the free section is smaller than it, which is 44.9 mm. As for the rest position, the elongation of the elements is relatively smaller. The elongation of the whole free section is 140 mm, which is the ultimate elongation of the free section. And it is more consistent with the reality. When the calculation reaches balance, the elongation of the structural elements (CID = 686) at the junction between the free section and the anchoring section in the modified model and the original model is maximal, which is 233.5 mm and 292.9 mm. The relative difference between   Figure 21 shows that in the modified model, the structural element (CID = 686) ruptures when the calculation reaches 4282 timesteps. The axial force in the element decreases suddenly from 297 kN to 0. The axial force in the structural element of CID = 673 decreases when the calculation reaches 4426 steps. With the calculating timesteps increasing, the axial force decreases gradually to 0. Then, the rupture process of the free section from the local to the whole is realised. The elongation of the structural element of CID = 673 and the elongation of the structural element of CID = 686 increase gradually with the calculating timesteps. Both of them increase rapidly. The increasing velocity of the latter one is larger. When the calculating reaches the critical rupture state, the elongation of those two structural elements is 44.9 mm and 78.0 mm, respectively. The latter one is 33.1 mm larger than the previous one. After the structural element of CID = 686 ruptures, the elongation velocity increases suddenly. The ultimate elongation of the element is 233.5 mm. After the structural element of CID = 686 ruptures, the elongation of the structural element (CID = 673) continuously increases by 0.4 mm under the residual axial loading effect. The ultimate elongation of the element is 45.3 mm. After the rock bolt tendon ruptures, the elongation velocity of the rock bolt-free section increases suddenly. The ultimate elongation of the rock bolt-free section is 291.5 mm. In the original model, the structural element of CID = 686 ruptures when the calculation reaches 749 timesteps. The axial force decreases suddenly from 297 kN to 0. When the structural element of CID = 673 calculates to 963 timesteps, the axial force shows vibration. When the calculation reaches 1068 timesteps, the axial force decreases gradually. With the calculating timestep increasing, it gradually reaches 0. After the structural element of CID = 686 ruptures, the elongation velocity increases suddenly.    14 Geofluids The ultimate elongation of the element is 292.9 mm. However, after the structural element of CID = 686 ruptures, the elongation of the structural element of CID = 673 continuously increases by 0.8 mm under the residual axial loading effect. The ultimate elongation of the element is 9.2 mm. The whole elongation of the rock bolt-free section is 305.8 mm ultimately, which is 14.3 mm larger than the modified model.     From Figure 23, it can be known that in the original model, the elongation of the structural element (CID = 697 mm) at the end of the free section is maximal, which is 11.4 mm. As for the rest positions, the element elongation is all less than 2.0 mm. The whole elongation of the rock boltfree section is 23.3 mm, which is only 16.6% of the ultimate elongation. The rock bolt tendon ruptures prematurely. In the modified model, it is same that the elongation of the structural element at the end of the free section is maximal, which is 63.5 mm. Following it is the elongation of the struc-tural element (CID = 710) at the junction between the free section and the anchoring section, which is 37.2 mm. As for the rest free sections, the elongation of the element is less than 9.5 mm. In the anchoring section, the elongation of each element is all less than 0.4 mm. The whole elongation of the free section is 140 mm, which is the ultimate elongation of the free section. And this is in accordance with the reality. When the calculating balance is reached, in the modified model and original model, the elongation of the structural element (CID = 697) at the end of the rock bolt free section is relatively large, which is 100.9 mm and 226.0 mm. The relative difference is 125.1 mm. In the modified model and the original model, the elongation of each element in the anchoring section is close to 0. This is because the anchoring section is anchored in the stable rock strata in the floor and the floor deformation is relatively small. This leads to the consequence that the suffered load of the rock bolt tendon in the anchoring    16 Geofluids section is relatively small. Figure 22 shows that in the modi-  18 Geofluids response action of the rock bolt tendon is consistent with the engineering practice.
In summary, the proposed approach can accurately simulate the actual situations of the failure of rock bolts subjected to tensile load in the field.

Discussion
(1) The rock bolt failure modes including tensile breakage, rod-anchorage interface failure, and rockanchorage interface failure have not been considered in the proposed rock bolt simulation method. If the above situations need to be modeled, additional targeted development of FLAC3D needs to be conducted (2) The pileSEL can also be used to simulate the soil nails, anchor bolts, and anchor cables in almost all similar engineering. Thus, the developing framework of the simulation approach in the paper can be directly copied in engineering fields that have issues similar to broken rock bolts, such as the anchor bolt support, soil nail support, anchor cable support, lattice frame beam support, latticebolt union support, and pile-cable union support in a large deformation slope or foundation pit engineering

Conclusions
In this paper, based on the self-embedded rock bolt pileSEL in the FLAC3D, through analyzing the original mechanical model of the pileSEL, modification is conducted on its axial tensile mechanical model. Then, the modified tensile rupture mechanical model is established. The secondary development is conducted, combining with the programming language FISH. The modified model is embedded in the main program of FLAC3D. Validation is conducted to confirm the rationality of the modified tensile rupture model. Furthermore, the engineering applicability and the reliability of the modified mechanical model are analyzed. The main conclusions include the following: (1) The tensile rupture failure criterion of the pileSEL is proposed. The tensile rupture mechanical model in which the whole elongation of the rock bolt-free section is regarded as the judgement criterion is established. When the elongation of the free section is larger than or equal to the rock bolt rupture elongation rate, the structural element that has the largest elongation in the free section shows the tensile rupture. And the axial force decreases suddenly to 0. Then, the sudden tensile rupture of the rock bolt tendon is realised (2) The programming language FISH is adopted to embed the modified model into the main program of FLAC3D. This realised the "utilisation" of the rock bolt tensile rupture failure. The macro rock bolt rupture effect can be realised from local to whole. This is consistent with the mechanical mechanism of the rock bolt tensile rupture failure (3) In an actual mining roadway, a numerical calculating model where the rock bolt is mainly subjected to tensile loading is established. The tensile rupture process of the rock bolt tendon is analyzed. The function of the tensile rupture at the end of the rock bolt-free section and the junction between the free section and the anchoring section is realised. In the original model, all rock bolts rupture. Furthermore, when rupture occurs, the elongation of the rock bolt-free section is far smaller than the predefined rock bolt rupture elongation rate. However, in the modified model, only four rock bolts at the roadway sides show rupture. When rupture occurs, the elongation of the rock bolt-free section is equal to the predefined rupture elongation rate. Based on the stress contour of the surrounding rock masses and the sprayed layer in those two models, it indicates that compared with the modified model, the deformation of the original model is severer compared with the modified model. This is resulted by the fact that in the original model, all rock bolts rupture. Therefore, calculating results of the modified model are more consistent with the engineering practices. This improves the engineering simulation accuracy and precision