Abstract
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.
1. Introduction
Rock bolts are widely used in civil and mining engineering to stabilize underground excavations [1–5]. They are installed to reinforce fractured rock mass by resisting dilation or shear movement along the fractures [6–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–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).
(a)
(b)
(c)
(d)
(e)
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–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 subjected 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 approach was verified through a simulation analysis of a typical real engineering case.
2. 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–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].
2.1. 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 is the longitudinal direction, pointing from node 1 to node 2. The axes of and are shearing that are perpendicular to the longitudinal direction. The rectangular coordinate system is composed of , 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, indicates the axial force. and indicate the shear force. , , and 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.
2.1.2. 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 is the axial force in the rock bolt; is the tensile strength of the rock bolt; is the cross area of the rock bolt; is the elastic modulus; is the elastic strain; is the plastic tensile strain; tfstrain is the tensile failure strain.
The elastic perfectly plastic mechanical model of the pileSEL tendon can be classified as two stages, namely, (1)elastic stage (section )
In the elastic stage, there is a linear proportional relationship between the axial force and deformation. The slope is the axial stiffness of the pileSEL, namely, where is the axial stiffness of the rock bolt tendon; is the cross-sectional area of the rock bolt tendon; is the elastic modulus of the rock bolt tendon material; and is the length of the structural component. For each parameter, the international standard unit is used. (2)plastic stage (after the point )
When the axial force of the pileSEL reaches the tensile yield force , the rock bolt tendon enters the yield state. The axial force will not rise again and the longitudinal stiffness 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 . The axial force will not increase again, and the axial stiffness 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 () can be recorded. Its calculating expression is shown in the following equation: where is the total tensile yielding strain of the pileSEL; is the diameter of the structural element; is the length of the structural element; is the average rotational angle of the structural element; 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.1.3. The Functional Character of the pileSEL
From the above analysis, it can be known that the pileSELs in FLAC3D have the following characters when it is used to simulate rock bolts: (1)The pileSEL is a one-dimensional linear component, which can better simulate the axial tension and compression effect between rock bolts and surrounding rock masses(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 bolt-free 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 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 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 ), which is the same as the original model(2)plastic stage (section ), which is same as the original model(3)rupture stage (after the point ), in which the axial force in the rock bolt decreases to zero
The axial stiffness decreases to zero and the increment of the axial force 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 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 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.
3.2. 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 , 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 , 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 ), the calculation ended(2)If the procedure has not reached convergence, it then enters the rupture prejudgement module (i)From the first rock bolt (), first, summing is conducted on the plastic tensile strain of each structural elements of the rock bolt-free section (the number changes from to ). Judging is conducted to check whether 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 is true, judging is conducted to calculate the elongation rate of the rock bolt-free section, (iii)If the condition is false, the judgement is moved to the next rock bolt element with the number of (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 in the free section, (the elongation rate of the rock bolt-free section is the difference between the length of the rock bolt-free section and the initial length of the rock bolt-free section . Then, the acquired result is divided by ) 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 is set as a value that is larger than 0. However, it is quite close to 0. In this paper, 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 . Then, the calculation enters the rock bolt (). If the condition is false, calculation enters directly the rock bolt ()(ii)Judging whether the rock bolt number is larger than ( 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.
3.2.2. Specific Methods
The above theory and the modification action are realised through the embedded programming language FISH in FLAC3D. Through the internally embedded functions and commands in FLAC3D, the length of each structural element in the free section and the elongation rate of the rock bolt free section can be monitored. Furthermore, assignment is conducted on the tensile failure element of the elements. Through the FISH functions sp_tfstr(sp_pnt) and sp_length(sp_pnt) in 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.
3.3. Mechanical Matching Test of Rock Bolt Reinforcement Accessories
Laboratory tests were conducted to analyze the mechanical matching characters of the rock bolt accessories under the axial loading condition [35–38]. The rock bolt performance test platform consists of chucks, a displacement measurement, and a testing data acquisition apparatus to simulate the actual working condition of rock bolts. A rock bolt with a full length of 2400 mm and an anchored length of 800 mm is used in the test, during which the rock bolt rod, anchored section, nut, and plate installed on the test platform are synchronously stressed by the mechanical gripping force exerted by a fixed chuck on the test platform. The matching between the rock bolt rod, anchored section, and its components is tested, and the axial force and elongation characteristics of the rock bolt under tensile load conditions are then obtained. A tensile force loads on the mobile chuck, and a force sensor is installed between the chuck and the faceplate to obtain the axial force. In the pull-out test, the axial force of rock bolts and displacement of mobile chuck are recorded. This section presents the results of the comparison between a pull-out test for a 2400 mm long, 22 mm diameter, rock bolt encapsulated in concrete and the simulated model in FLAC3D. The Young modulus of the steel bolt was 210 GPa. The axial force and displacement at the loaded end of the rock bolt were recorded during the pull-out test, as shown in 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.
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.
4. 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 section as the judgement criterion to distinguish the tensile rupture type. This section is based on the engineering geological background. Taking the rock bolt support roadway case in the Xiaojiawa coal mine as the engineering background, the approach is verified through comparative analysis of the results of the original model scheme and the improved approach scheme (modified model scheme).
4.1. General Situation of the Engineering Case
The Xiaojiawa coal mine is located in Lvliang, Shanxi Province, China, with a production of 88,000,000 tons per year (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 semicircle 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 . 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 . The rock bolt reinforcement and the concrete spraying are used to support the roadway. The dimension of the rock bolt is 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.
4.2. Modelling and Parameters
4.2.1. Establishment of the Model and Determination of the Material Parameters
The modified scheme adopts the improved approach, and the original scheme adopts the built-in models. The main differences in the two schemes are the modelling methods and the parameters of rock bolts. Under the rock bolt reinforcement parameters in Figure 11(a), tests and analysis are conducted between the original model and the modified model.
The in situ stress in the mining area was measured using borehole relief methods with CSIRO cells. The in situ stress measurement shows that the vertical stress is 22 MPa, the horizontal principal stresses are 26.6 MPa and 17.6 MPa, respectively. The shape of the roadway is a semicircle with two straight walls, and the size is the same as the one in practice. Rock bolt supports were used to control the roadway. The rock bolts were 22 mm in diameter and 2.4 m in length. The pileSELs are used to simulate rock bolts. To increase the simulation precision, each rock bolt is divided into 24 segments. For each segment, the length is 0.1 m. To replicate filed conditions, the dimension of the established test model is (). A vertical load was applied to the model top to simulate overburden (depth of cover is 750 m), while roller boundaries were placed on the other five sides of the model. Normal displacement constraints are applied on all of the four surrounding surfaces, and all displacement constrains are applied on the bottom surface. There is no displacement constraint on the top surface, but the compensating ground stress of 20 MPa is applied. The gravity was set as 10 m/s2. The lateral pressure coefficient was 1.0. For the sandy mudstones and the coal, the strain softening yielding criterion is used. For the poststone and the concrete spraying layer, the Mohr-Coulomb yielding criterion is used. The SHELL element is used to simulate the bearing plate and the hexahedron zones are used to simulate the sprayed layer.
This model is mainly aimed at the reinforcement effect of the modified model applied to engineering examples, so this paper draws lessons from the results of relevant scholars for the value of rock layer parameters, as shown below [39–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 . And the elastic modulus is Pa. The Poisson ratio is 0.3. The inertia moment along the axis is m4, and the inertia moment along the axis is m4. The polar inertia moment is m4. 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.
4.3. Result Analysis
4.3.1. 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 . 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 elements of rock bolts when the modified model and the original model arrive at the calculation balanced. 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.
(a) Modified model
(b) Original model
(a) Modified model
(b) Original model
(a) Modified model
(b) Original model
(a) Modified model
(b) Original model
(a) Modified model
(b) Original model
(a) Modified model
(b) Original model
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 14th 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 15th rock bolts in the modified model and the second and 14th rock bolts in the original model occurs in the end of the free section that is adjacent to the plate. Figures 14 and 15 show that the maximum horizontal and vertical displacement of the surrounding rock masses in the modified model is 395 mm and 298 mm. However, in the original model, the maximum horizontal and vertical displacement of the surrounding rock masses in the original model is 428 mm and 332 mm. Compared with the modified model, the original model shows severe deformation, which is resulted by all rupture of the rock bolts in the original model. Figures 16 and 17 show that due to the difference on the roadway deformation, the load in the sprayed layer is different. The maximum horizontal stress of the sprayed layer is concentrated on the top of the roadway. In the modified model, the maximum horizontal stress is 13.655 MPa while in the original model, the maximum horizontal stress is 13.664 MPa, which is 0.01 MPa larger than the modified model. However, the maximum vertical stress of the sprayed layer is concentrated in the bottom of the roadway arch. In the modified model, the maximum vertical stress is 13.20 MPa while in the original model, the maximum vertical stress is 13.36 MPa, which is 0.16 MPa larger than the modified model.
Based on 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.
4.3.2. 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.
(1) Analysis on the Rupture Process of the Rock Bolt No. 8 at the Arch Top. The rock bolt No. 8 at the arch top is composed of 24 pileSELs. The CID number is ranged from 529 to 552, in which the number of the free section is ranged from 529 to 542 and the number of the anchoring section is ranged from 543 to 552. Under the critical tensile state and when the calculating balance is reached, for the rock bolt No. 8 in the arch top, the elongation of each element in the rock bolt is shown in Figure 18. Under the calculating condition in the modified model and the original model, the axial load and elongation of the dangerous structural element () and the variation law of the elongation of the whole free section belonging to the rock bolt No. 8 with the calculating timestep are shown in Figure 19.
Figure 18 shows that under the critical tensile rupture state, the elongation of the structural element () 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 () and the structural element () 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 () 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 () 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 is the dangerous element and the structural element of 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 () 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 () 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 bolt-free 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 () 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 them is 59.4 mm. Figure 21 shows that in the modified model, the structural element () 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 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 and the elongation of the structural element of 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 ruptures, the elongation velocity increases suddenly. The ultimate elongation of the element is 233.5 mm. After the structural element of 86 ruptures, the elongation of the structural element () 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 ruptures when the calculation reaches 749 timesteps. The axial force decreases suddenly from 297 kN to 0. When the structural element of 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 ruptures, the elongation velocity increases suddenly. The ultimate elongation of the element is 292.9 mm. However, after the structural element of ruptures, the elongation of the structural element of 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.
(3) Analysis on the Rupture Process of the Rock Bolt No. 15 at the Roadway Side Bottom. In the modified model and the original model, the rock bolt No. 1 and No. 15 at the bottom of the roadway side rupture. The rock bolt No. 15 is taken as an example, to analyze the rupture process of the rock bolts at the bottom of the roadway side. The rock bolt 15% is composed of 24 pileSELs, and the CID number is ranged from 697 to 720. Among them, the number of the free section ranges from 697 to 710 and the number of the anchoring section ranges from 711 to 720. Under the critical tension state and when the calculation reaches balance, for the rock bolt No. 15, the elongation of each element in the rock bolt in those two models is shown in Figure 20. Under the calculating condition in the modified model and the original model, the structural element of is the dangerous structural element. The structural element of is less dangerous. The axial force and distribution of the dangerous and less dangerous structural elements and the variation law between the whole elongation of the free section of the rock bolt No. 15 and the calculating timesteps are shown in Figure 22.
From Figure 23, it can be known that in the original model, the elongation of the structural element ( 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 bolt-free 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 structural element () 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 () 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 section is relatively small. Figure 22 shows that in the modified model, the structural element of ruptures when the calculation reaches 18092 timesteps. The axial load in the structural element decreases from 297 kN suddenly to 0. When the calculation reaches 18118 timesteps, the axial force in the structural element of decreases gradually. With the calculating timestep increasing, the axial force decreases gradually to 0. Consequently, the rupture process of the rock bolt-free section from the local to the whole is realised. The elongation of the structural element () and the structural element () increases with the calculating timestep. Both of them increase rapidly. The increasing velocity of the front one is relatively larger. When the calculation reaches the critical rupture state, the elongation of two structural elements reaches 63.5 mm and 37.2 mm, respectively. Compared with the latter one, the front one is larger by 26.3 mm. After the structural element of mm ruptures, its elongation velocity increases suddenly. The ultimate elongation of the element is 100.9 mm. However, after the structural element of ruptures, the length of the structural element () shrinks, decreasing by 0.6 mm. The ultimate elongation of the element is 36.6 mm. After the rock bolt ruptures, the elongation velocity of the rock bolt-free section increases suddenly. The transition effect of the curve is pretty obvious. When the calculating reaches balance, the ultimate elongation of the rock bolt-free section is 172 mm. In the original model, when calculating to 1187 timesteps, the structural element of ruptures. The axial force decreases from 297 kN to 0 suddenly. After that, the axial force in the structural element of decreases gradually. With the calculating timestep increasing, the axial force is close to 0. After the structural element of ruptures, its elongation velocity increases suddenly. The ultimate elongation of the element is 226 mm. However, after the structural element of ruptures, the elongation of the structural element () decreases by 0.4 mm. The ultimate elongation of the element is 1.5 mm. The ultimate elongation of the rock bolt-free section is 232.9 mm, which is 60.9 mm larger than the modified model.
Based on the analysis on the tensile rupture process of rock bolt No. 8 at the arch top and rock bolt No. 14 and No. 15 at the roadway side, it can be acquired that in the original model, when the elongation of the rock bolt-free section ranges from 21.3 mm to 26.8 mm, tensile rupture occurs. At this time, the elongation of the free section is only 15.2%-19.1% of the ultimate elongation. Its mechanical response action is apparently different from the practice. However, for the modified model, when the rupture of the tendon occurs, the elongation of the rock bolt-free section equals the ultimate elongation, which is 140 mm. The mechanical 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.
5. Discussion
(1)The rock bolt failure modes including tensile breakage, rod-anchorage interface failure, and rock-anchorage 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, lattice-bolt union support, and pile-cable union support in a large deformation slope or foundation pit engineering
6. 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(4)The modified element and the realising method extend the application range of FLAC3D. The simulation ability is enhanced and the calculating precision is improved. However, due to the self-limitation of the pileSEL, FLAC3D, and the complexity of the anchoring failure, the accurate simulation of the anchoring failure still needs to be further studied
Data Availability
Data is available in this paper or through emailing the corresponding author.
Conflicts of Interest
The authors declare no conflict of interest.
Authors’ Contributions
R.W. conceived the research and wrote the paper. J.B. directed the paper structure. S.Y. analyzed the data and contributed to numerical simulation. Y.B. and G.D. participated in collecting the data and language editing. All authors have read and approved the final manuscript.
Acknowledgments
This research has been supported by the National Natural Science Foundation of China (grant no. 51927807), a Project Funded by the Priority Academic Program Development of Jiangsu Higher Education Institutions (PAPD), National Natural Science Foundation of China (General Program: 52074239) and Independent Research Project of State Key Laboratory of Coal Resources and Safe Mining, CUMT (SKLCRSM18X08).