The mechanism of decrease in bearing capacity and performance evaluation of the yielding bolt in coal mine

With the increase of the coal mining depth, the surrounding rock shows characteristics of high stress and large deformation. However, the traditional bolt could not adapt to this surrounding rock environment and become invalid due to low elongation. In order to solve the problem, the new type of yielding bolt has been designed. It could release the deformation energy of the surrounding rock through its yielding pipe. However, few studies have been reported on the mechanism in bearing capacity decrease of the yielding pipe. Furthermore, the evaluation of the yielding bolt was not fully understood. Under such a background, a numerical model of the yielding model was established. According to the numerical simulation, the mechanism of the decrease in the bearing capacity of the yielding bolt is that the buckling occurs at different positions of the yielding pipe, making the material stiffness matrix turn negative, and the bearing capacity decreases significantly. Also, the performance of the yielding bolt could be evaluated from the use efficiency of the yielding pipe, the plastic strain characteristics of the components of the yielding bolt, and the energy absorption law. The evaluation results indicated that: (a) the deviation of the yielding pipe's Mises stresses was within the range of 1.0%‐4.0%. The stresses were distributed evenly, and the material use efficiency was high; (b) only the yielding pipe occurred the plastic failure during the yielding process, and other components were at the elastic stages. After the yielding, the yielding bolt was still effective for support. (c) The yielding pipe could absorb energy smoothly, which is helpful for the stability of the yielding bolt.

constraints, the roadway surrounding rock's deformation will become more serious. 9 In order to improve the bolt's adaptability to large surrounding rock deformation, extensive researches have been carried out by scholars both at home and abroad, [10][11][12][13] and many extensible bolts have been invented, which could be divided into the bolts with extensible rod body and the bolts with extensible structural elements according to the operating principles. 14 The former has simple structure, manufacturing, and installation, but they have higher requirements for the materials and have a low level of initial prestressing force. The latter is featured with adjustable expansion, constant resistance, yieldable capacity, stable, and reliable operation, and large elongation.. 15,16 With such advantages, the bolts with extensible structural elements have been widely applied. [17][18][19] At present, the popular type is the new yielding bolt developed by Jennmar, US, in 2006. Figure 1 shows its structure. It constitutes the nut, yielding pipe, tray, and high-strength rod. The yielding pipe is made of a seamless steel pipe and is the core component. The scale of deformation of the yielding pipe is generally 25-35 mm.
During the yielding process, the excess deformation energy of the surrounding rock could be released through the deformation of the yielding pipe. After the yielding, the yielding bolt has stronger constraints for surrounding rocks and could effectively control the roadway deformation, thus ensuring the stability of the roadway.
Currently, the studies about the yielding bolts mainly focus on the yielding mechanism, coupling relationships between the yielding bolt and the surrounding rock, and key parameters of the yielding bolt. (a) Yielding mechanism. Liu et al 20 established a mechanical model of the yielding bolt based on Hooke's law and obtained that the stress increase rate in the rod body was less than that of the conventional bolt due to deformation and energy absorption of the yielding pipe. Therefore, the stress of the yielding bolt rod was designed significantly less than that of the conventional bolt under the same surrounding rock deformation. (b) Coupling relationships between the yielding bolt and the surrounding rock.
Lian et al 21 discussed the coupling relationships between the yielding bolt characteristic parameters (the yielding point and the maximum yielding distance), bolt design parameters (row and space, diameter, length, and prestressing force), and the surrounding rock; Zhang et al 22 systematically analyzed the variations of surrounding rock stress and displacement before and after the support of the yielding bolt by FLAC 3D software. (c) The key parameters of the yielding bolt. Based on the uniaxial compression test, Yang et al 23 studied the basic mechanical parameters, such as the yielding point, the stability of the yielding point, and the maximum yielding distance, which could provide the basis for the selection of the yielding bolt; Shan et al 24 divided the yielding bolt's load-displacement curve into the elastic stage, yielding stage, and plastic strengthening stage. Then, according to the energy constitutive equations at the three stages, the computational formula of the yielding point of the yielding bolt was given.
However, the current studies only present the relationships between the load and displacement for yielding bolt during this process and do not explain the mechanism of bearing capacity reduction during the yielding process. Additionally, the performance evaluation of the yielding bolt has not been studied thoroughly. Under such background, a numerical model of a yielding bolt was established by Nastran software in this paper. Then, the mechanism of bearing capacity reduction of yielding bolt was studied. Finally, the performance of the yielding bolt was evaluated from the use efficiency of the yielding pipe, plastic strain characteristics, and the energy absorption laws of the yielding bolt.

| THE NUMERICAL METHOD
The traditional theoretical calculation is difficult to adapt the yielding bolt's complex geometric model and nonlinearity due to the contact type and material. Therefore, the numerical computation method was adopted in this paper to analyze the elements, iterative method, and rapid convergence.

| The element description
The hexahedral element has excellent computational efficiency and could be used for simple geometries. For complex geometries, the tetrahedral element shall be used, because the tetrahedral element has strong adaptability. Two kinds of elements need to be described in the numerical model.

| The tetrahedral element
The tetrahedral element is made of four nodes, and each node has three displacements. Figure 2 shows the element's nodes and the displacements of the nodes.
(1) The element's geometry and node description.
The element's nodal displacement array q e and nodal force array P e could be expressed as follows 25,26 : (2) Expression of the element displacement field. This element has four nodes, and the element's nodal displacement has 12 degrees of freedom. Therefore, the displacement field in each direction could set four undetermined coefficients. According to the number of nodes and the basic principle of determining displacement mode, the element's displacement mode could be selected, as follows 25,27 : where Ni = 1 6V (ai + bix + ciy + diz), i = 1,2,3,4; V is the volume of the tetrahedron, and ai,bi,ci,di are the coefficients related to the node's geometric position.
The element's strain could be expressed as follows 25,28 : Where B i is: According to the physical equation for spatial problems, the stress field could be expressed as follows 25 : where D is the elastic coefficient matrix of the spatial problem.
(3) The element's stiffness matrix and the equivalent nodal load matrix.
After getting the geometric matrix B, the element's stiffness matrix could be calculated by the formula of stiffness matrix 29 : (4) The element stiffness equation.
The tetrahedral element with four nodes The element's stiffness equation could be obtained by taking the first-order extremum for the nodal displacement q e in the formula for the potential energy of an element, as follows 25 : where K e is the stiffness matrix of the element.

| The hexahedral element
(1) The element's geometry and node description.
As shown in Figure 3, the element consists of eight nodes, and each node has three displacements. The element's nodal displacement has 24 degrees of freedom. The element's displacement array q e and nodal force array P e could be obtained, as follows 25,26 : The element's shape function matrix is as follows 25,28 : The nodal displacements are up to 24, so it is difficult to use the node conditions to directly fix the undetermined coefficients in the displacement mode and the method of the shape function matrix. The element's natural coordinates could be used to get the shape function matrix by Lagrange's interpolation formula.
(2) Expression of the element's strain field.
Based on the geometric equation of the plane problems in elasticity, the element's strain could be expressed as follows 25 : (3) The element's stiffness matrix and the equivalent nodal load matrix.
The element's stress expression could be obtained by the physical equation of the plane problems. Then, the element's potential energy was calculated. Similar to the previous process of deriving other elements, the element stiffness matrix could be obtained, as follows 30 : (4) The element stiffness equation.
The element's stiffness equation could be obtained by taking the first-order extremum for the nodal displacement q e in the formula for the potential energy of an element, as follows 25 :

| The iterative method
In the advanced nonlinear analysis, the whole computational process is divided into multiple time increment steps. Within each time increment, the Newton-Raphson method is used for iterative equilibrium. The equilibrium equation of the Nastran at t + Δt could be expressed as follows 31,32 : where t+Δt K (i−1) is the tangent stiffness matrix based on the solution calculated at the end of the iteration (i − 1) at time t + Δt; t+Δt U (i) is the displacement vector at the end of iteration i at time t + Δt; t+Δt R is the externally applied load vector at time t + Δt; t+Δt F (i) is the consistent nodal force vector at the end of iteration i at time t + Δt.
The main idea of the Newton-Raphson iteration algorithm is to use the method of successive approximation. At each

F I G U R E 3 The hexahedral element with eight nodes
x y z increment step of load, the obtained displacement value is used to get the displacement-dependent elastic-plastic matrix values. Next, the linear computation will be conducted. Finally, the set accuracy could be realized by iterating over the difference between the load value through repeatedly adjusting and calculating and the set load value. Therefore, the total load is divided into a series of load segments for the iteration of the linear equation. The main steps of the method are as follows 33 : Step 1: To divide the external load F into a series of load segment, F (1) ,F (2) ,F (3) , … F (i) ; Step 2: The multistep cyclic iteration is carried out in each segment until the convergence occurs in the segment. Formula (15) is used for each step's iterative calculation; Step 3: To conduct the cyclic iterations of all load segments and accumulate the results, as shown in Figure 4.
In the N-R iterative method, the tangent stiffness matrix should be reformed at each time. After that, the inversion algorithm is necessary, which needs a huge calculation force. If the initial tangent stiffness matrix is always used and remained unchanged, the amount of calculation could be significantly reduced. This method is named the modified N-R method, as shown in Figure 5.

| Linear search
The yielding process of the bolt involves many nonlinear problems, such as the materials' large deformation, multiple contact sets, and sudden change of the boundary conditions. The linear search could solve the aforementioned nonlinear problems by controlling displacement increments in iteration during the numerical calculation, thus improving the convergence of the model and the calculation efficiency. The specific incremental displacement vector in iteration could be obtained by the following formula 34 : where t+Δt U (i) is the displacement vector at the end of iteration i at time t + Δt. t+Δt U (i) is the displacement vector at the end of the iteration (i − 1) at time t + Δt; ΔU (i) is the incremental displacement vector in iteration i; and β (i) is a scaling factor obtained from a line search in the direction ΔU (i) in order to reduce out-of-balance residuals, according to the following criterion. 35 The tangent stiffness matrix needs to be recalculated every time where STOL is a line search convergence tolerance, by default is 0.5, and t+Δ F (i) is calculated using the total displacement vector t+Δt U (i) .
The following bounds also govern the magnitude of β.

| Automatic-Time-Stepping (ATS) method
Automatic-Time-Stepping method is one kind of computing methods to achieve fast convergence by increasing or decreasing the size of the incremental step. If there is no convergence in the numerical calculation of the previous incremental step, the ATS method will continue to reduce the size of the incremental step until the iteration convergence. If the convergence is calculated in the previous incremental step, the ATS method will continue to increase the size of the incremental step. 36,37 In order to choose proper incremental steps in the numerical calculation, the maximum and minimum incremental steps of the ATS method are given. Generally, a tripling of incremental step is set as the maximum incremental step, and the 0.1 times of incremental step is set as the minimum incremental step. 36 The yielding bolt's gasket and yielding pipe belong to metal materials. They may occur plastic deformation during the yielding process, so it is necessary to select an appropriate elastic-plastic constitutive model. According to the stress-strain relationships, metal materials exist two linear elastoplastic constitutive model, multilinear constitutive model, and the elastoplastic constitutive model expressed by self-defining function. In the field of hardening criterion, the kinematic hardening model, the isotropic hardening model, and mixed hardening model exist. Due to the isotropic materials of the components of the yielding bolt, the classical bilinear isotropic hardening constitutive model was applied in the simulation. 38 In this model, the stressstrain relationship curve of the metal materials was simplified into a double line, as shown in Figure 6. Two linear segments were used to simulate the constitutive relations of the elastoplastic materials. The model's stress-strain curve was drawn. OA represents the linear elasticity of the material; the slope E is the elasticity modulus. When the stress exceeds the initial yield stress σ y0 , the plastic strain occurs. AB represents the plasticity of the materials, and the slope E T refers to the tangent modulus, and 0 < E T < E. The stress-strain relationship could be expressed as follows 39 : where E is the steel's elasticity modulus; E T is the hardening slope; ε y0 is the yield strain.

| Nylon materials
Nylon, the first synthetic fiber in the world, is one of the common materials used to substitute metal materials such as steel, iron, and copper. For nylon materials, the strain energy polynomial in formula (20) could be used to analyze the experimental data by a least-squares regression method. The multinomial coefficient A ij could be calculated.
where ̇I 1 and ̇I 2 are first-order invariant and the second-order invariant of the shear deformation strain tensor except the volume deformation, respectively. A 01 , A 02 , A 10 , A 11 , A 20 is the multinomial coefficient relative to the shear deformation.
There was material nonlinearity, contact nonlinearity, and geometric nonlinearity during the compression process, so the advanced nonlinear solver was used in this paper. Under the condition that Na ≤ 3 required by the solver, the tensile The stress-strain curve of the bilinear isotropic hardening constitutive model testing data, shear test data, or pure shear test data were used to get five shear deformation coefficients A ij , as shown in Table 1.

| Contact model
There were multi-group direct contacts among the nylon gasket, metal gasket, yielding pipe, and tray. Besides, during the yielding process, the yielding pipe could contact with itself. For the above contact behaviors, the description of the contact rule is the basis to guarantee the reliability of the numerical simulation.
The following condition should be satisfied during the contact process 40 : where g is the contact clearance distance, and λ is the contact reaction force.
The constraint function method, the Lagrange segmentation method, and the stiffness target method are proposed to describe the contact behaviors. Among them, the constraint function method is commonly used to judge and calculate the distance between contact surfaces, whose expression is as follows 41 : where w(g, ) is the friction constraint function; N and T are the friction constraints along the direction and the tangential direction, respectively; and v is the sliding velocity.
Except N and T , other parameters were automatically calculated by the software. The empirical values N and T were given by Nastran, which were 1.0 × 10 −12 and 1.0 × 10 −3 , respectively. 40,42 Except for the constraint function, the reasonable description of friction behaviors is also very important. Nastran has given 12 friction calculation models, including Static/ Dynamic, vs Sliding Velocity, and vs Contact Force, to calculate friction coefficients related to velocity and contact force under the variations of motion state. In this paper, considering the process of the associated components from static motion to dynamic motion during the compression process, the static/dynamic method is chosen to calculate the friction coefficients. [43][44][45] The static friction coefficient of 0.2 and the dynamic friction coefficient of 0.15 could be obtained.

| Model establishment
The uniaxial compression was directly conducted in the experiment to analyze the load-displacement relationships and obtain the yielding characteristic parameters. However, the mutual contacts among the yielding pipe, tray, and gasket were neglected in the experiment. Meanwhile, it was hard to monitor the whole yielding pipe's stress variation characteristics during the yielding process in the experiment. In order to accurately get the mechanism of the decrease in bearing capacity and evaluate the performance of the yielding bolt, the numerical model was established, as shown in Figure 7. Firstly, the commands of stretch/rotate/cut in three-dimensional parametric modeling software Pro/E were applied to build the 3D part models of the tray, the yielding pipe, and gaskets. Secondly, according to the mutual positional relationships of the components, the 3D components were assembled and saved to igs file. 46,47 The files were then imported into the numerical modeling software for geometry clean up and idealization. Since the model has a symmetrical structure, a half model was adopted for the numerical simulation to reduce the calculated amount and save the calculation costs. In the numerical software, the tetrahedron and hexahedron elements were utilized for mesh generation of the solid model. When simulating the compression among the nuts, gaskets, and trays due to the extension of the bolt, a rigid plane was applied to replace the bolt and the nut. The rigid surfaces and reference points were defined in the numerical modeling, respectively, and the mesh generation was carried out for the rigid surfaces. 48 A total of 32 688 tetrahedron elements, 1956 hexahedral elements, and 1296 rigid elements were discretized in the model. The bilinear isotropic hardening model was used for the tray, the yielding pipe, and metal gaskets in the model. The constitutive equation based on the strain energy function in Section 3.1.2 was adopted for the nylon gasket. The fixed constraint was added at the bottom of the model, and the symmetry constrain was used in the model. To simulate the mutual contacts among the components, eight sets of contact were defined. 49 Table 2 shows the specific material parameters. The different component deformation of the yielding bolt shows the movement law during the yielding process. As shown in Figure 8, the component deformations under various displacements were extracted, and the displacement values were given. Figure 8 indicates that: During the yielding process, the nylon gasket, metal gasket, and the yielding pipe gradually move downward while the tray and the metal gasket remain basically unchanged, because the tray has a large thickness and high intensity.

| The load-displacement curve
The rigid surface was used in the numerical simulation. The movement of the rigid surface was controlled by the reference points whose load-displacement curve was extracted through the output command, as shown in Figure 9(A). To explain the load-displacement relationships, Figure 9(A) gives Mises stress of the structural components during the yielding process, and Figure 9(B) shows the components' buckling and contact behaviors.
According to the load-displacement curve, the whole yielding process could be divided into five stages, as shown in the Figure 9(A).
1. At stage 1, when the displacement was 0-0.3 mm, the tray and gaskets were at the elastic stage, most of the yielding pipe was also at the elastic stage, and the load and the displacement were a linear relationship. When the displacement was 0.3 mm, the load of half model was 79.9 kN. Therefore, the load of the whole model is 159.8 kN. When the displacement was between 0.3 mm and 3.3 mm, the tray and the yielding pipe occurred plastic strain. The stiffness matrix of the materials was reduced significantly, and the load increased slowly with the increase of the displacement. When the displacement was 3.3 mm, the load of half model was 80.7 kN. Therefore, a load of the whole model (yielding point) was 161.4 kN. 2. At stage 2, when the displacement was 3.3-23.25 mm, the load decreased with the increase of the displacement. When the displacement was 3.3 mm, the load of half model was 43.9 kN. Therefore, the load of the whole model is 97.8 kN. Figure 9(B) indicates that the reason for the decrease in the load is the buckling in the middle and on the bottom of the yielding pipe, thus making the material stiffness matrix turn negative, and the bearing capacity decrease. 3. At stage 3, when the displacement was 23.25-27 mm, the load increased with the increase of the displacement. When the displacement was 27 mm, the load of half model was 50.9 kN. Therefore, the load of the whole model was 101.8 kN. Figure 9(B) indicates that the contact between the middle position of the yielding pipe and the gasket changes the contact status of the yielding pipe, which increases the bearing capacity.

Name
Poisson's ratio μ 4. At stage 4, when the displacement was 27-30.9 mm, the load reduced with the increase of the displacement. When the displacement was 30.9 mm, the load of half model was 44.1 kN. Therefore, the load of the whole model was 88.2 kN. As shown in Figure 9(B), buckling occurs on the top of the yielding pipe, causing the material stiffness matrix negative, and leads to the decrease of the bearing capacity. 5. At stage 5, when the displacement was 30.9-32.1 mm, the load increased with the increase of the displacement. Figure 9(B) shows the multiple contacts between the top of the yielding pipe, the metal gasket, and the self-contact of the inner cavity of the yielding pipe.

Hardening rate E T /GPa
The analyses mentioned above indicate that the mechanism of the yielding bolt is that the buckling gradually occurs at different positions of the yielding bolt, causing the material stiffness matrix to become negative and the bearing capacity to decrease obviously.

| Efficiency of the yielding pipe
The yielding pipe is the yielding component of the bolt. During the yielding process, an essential indicator of material use efficiency is whether the stress distribution of the yielding pipe is even or not. According to analyses in Section 4.1, the yielding process mainly occurred within the range of 0-30 mm, so Mises stresses of the components within this range were extracted, as shown in Figure 10.
The Figure 10 indicates that: When the displacement was 10.8 mm, the Mises stresses of 28 nodes at the outer boundary of the yielding pipe were extracted, as shown in Figure 11. The stresses are 568-606 MPa with an average value of 584.1 MPa, and the relative standard deviation is 1.5%.
As shown in Figure 12, in order to further analyze the use efficiency of the yielding pipe, the element's Mises stresses were continuously extracted after the mesh generation of the finite model, and the data were exported.
Application Programming Interface (API) was used to realize the above steps. The Mises stresses of 13 077 tetrahedron elements were counted. Their ranges, average values, and relative standard deviations were given, as shown in Table 3. The deviation is within the range of 1.0%-4.0% and stresses distributed evenly in the yielding pipe. When the relative standard deviation is <10%, the material use efficiency is high. Therefore, the material selected has high efficiency.

| Plastic strain of the components
The plastic strain refers to the deformation that occurs after the yielding of the materials. It is one of the essential indexes to reflect material damage. After the yielding bolt suffers from larger loads, the yielding pipe will play the yielding role. Therefore, it is necessary to avoid the plastic strain of other components. Figure 13 shows the plastic strain of the components during the yielding process.
The Figure 13 indicates that: which is less than the bolt's yield load. It is proved that the yielding pipe occurs failure, and other components are at the elastic stage during the yielding process. 3. After yielding, the other components except the yielding pipe have an ability of constraint for surrounding rock. Therefore, the performance of the compression bolt is promising.

| Energy absorption law
The yielding bolt will absorb a certain amount of energy during the yielding process. The amount and rate of energy absorbed during this process could reflect the yielding effect of the bolt to a certain degree. Figure 14 shows the energy density of the yielding bolt, as well as the amounts of the energy absorbed with the displacements of 0-30 mm by the API during the yielding process. The Figure 14 shows that: 1. The amounts of the energy absorbed during the yielding process increased linearly, indicating that the yielding bolt could absorb the energy continually, which is helpful for the stability of yielding bolt. 2. When the displacement is 30 mm, the total amount of the absorbed energy is 3800 J.

| ENGINEERING APPLICATION
Tashan Coal Mine is a coal mine with an annual output of ten million tons. The gob-side roadway will be affected by the cantilever beam structure of the goaf in 8214 working face, which is prone to be large deformation. Therefore, the on-site roadway is supported by yielding bolts. The yield strength was 600 MPa, and the tensile strength was 800 MPa. The bolt had a diameter of 22 mm and a length of 2.4 m. The bolt space was 0.8 × 0.8 m, and the pretightening force was 100 kN. The high-intensity tray with a size of 130 × 130 × 10 mm was used. The yield bolt has a height of 40 mm, an inner diameter of 25 mm, and an outer diameter of 30 mm. As shown in Figure 15, the deformation characteristics and force for the yielding pipes are obtained from the field.
The Figure 15 shows that: 1. The deformation characteristics of the field and numerical simulation shows a similar trend. Both are bulging out of the middle of the yielding pipe. According to the numerical simulation, it can be concluded that the efficiency of the material of the yielding pipe is high. 2. During the mining of the workface, the force of the yielding pipe increased gradually, and the maximum force is 178 kN. In the numerical simulation, the yield point of the yielding pipe is 161.4 kN, and the error between the numerical model and the field measurement is about 9.3%, which verifies the reliability of the numerical simulation.

| CONCLUSIONS
In this paper, based on the analysis of the numerical method, the yielding bolt's numerical model was established. The conclusions are as follows: 1. The mechanism of the yielding bolt is that the buckling occurs gradually at different positions of the yielding pipe, making the bearing capacity decrease significantly. From the numerical simulation, the yielding point and yielding displacement are 161.4 kN and 30.1 mm, respectively. 2. The performance of the yielding bolt could be evaluated from the material use efficiency of the yielding pipe, the plastic strain characteristics of the components, and the energy absorption laws. During the yielding process, the relative standard deviation of Mises stress is 4%, and only the yielding pipe occurred the plastic failure. The yielding pipe absorbed energy stably, which indicated that the performance of the yield bolt is excellent. 3. The deformation of the yielding bolt in the numerical simulation is almost the same as field measurement results. The error of the yielding point between the numerical model and the field measurement is 9.3%. The accuracy verifies the reliability of the numerical simulation and shows an admiring performance of the yielding bolt.

ACKNOWLEDGMENTS
This work was supported by the Chinese Scholarship Council (CSC); State Key Research Development Program of China (2018YFC0604500); and LiaoNing Revitalization Talents Program (XLYC1807219). The authors gratefully acknowledge the financial support from the organization.

CONFLICT OF INTEREST
The authors declare no competing interests.