Abstract
The load transfer control equations under bolt-surrounding rock interaction are established on the basis of classical beam theory and the trilinear shear slip model. The axial stress and transverse shear force distributions of the anchorage body are obtained by solving the equations. The equivalent forces obtained by the transverse force and axial shear stress of the bolts are applied to rock mass elements to simulate the support effect. A new dynamic algorithm for bolts is proposed in considering of the axial and transverse deformation of the anchorage body. The rationality of the algorithm is verified by comparing with laboratory pullout and shear tests of bolts. A dynamic time-history case study of underground caverns is conducted using this algorithm. Results indicate that (1) the algorithm may reflect the stress and deformation characteristics of bolts during an earthquake; (2) for the antiseismic support effect of the surrounding rock at fault, the bolt algorithm in this study is more valid than the algorithm that considered only the axial deformation of bolts; (3) in the support force of the bolt to the surrounding rock, transverse force is the key to limit fault dislocation and reduce the dynamic damage of the rock at fault.
1. Introduction
Southwest China is rich in water resources and has several hydropower stations. Restricted by geological conditions, these hydropower stations are mainly built across underground rock. The region is located at plate boundaries with a high probability of earthquakes; therefore, the antiseismic performance of underground caverns is directly related to the operational safety of hydropower stations [1]. As a flexible seismic reinforcement measure, bolt support is widely used in mining, water conservancy, and transportation industries due to its economic, simple, and reliable advantages. A commonly used bolt type in underground engineering is fully grouted rock bolt. Rock mass is strengthened by the rigidity and strength of the bolts, and the deformation of the surrounding rock is limited by transferring shear stress at the interface between the anchorage body and the surrounding rock [2]. Freeman [3] first proposed neutral point theory, which effectively describes the axial stress of a bolt in rock mass and has become important for analyzing bolt-surrounding rock interaction. Recent research on axial mechanical transform mechanisms between the anchorage body and rock mass is mainly based on pullout tests [4–6] and theoretical models [7–9].
The abovementioned studies only considered the effect of axial stress on bolts. However, in an actual rock mass, numerous faults and joint surfaces exist. The bolts not only reinforce the rock mass in the axial direction but also effectively stop or delay the displacement of the joint surface in the transverse direction [10]. Therefore, bolt mechanisms that consider transverse deformation require further study. The effects of material, dimension, inclination, prestress, bolt profile, rock strength, and joint roughness on the shear resistance of bolts have been studied through tests [11–16]. Pellet and Egger [17] determined the relationships between axial force and deformation and transverse shear force and bolt deformation. Jalalifar and Aziz [18] studied the force and deformation distribution of a shear-stressed bolt in the elastic and plastic stages and showed the relationship between the plastic hinge point and transverse deformation of the bolt. Ge and Liu [19] proposed a theoretical formula for estimating the shear strength of bolted rock joints. Zhang and Liu [20] comprehensively considered the axial and transverse deformability of bolts and proposed a deformation model of bolts under transverse shear load.
Several studies have presented static and dynamic calculations for bolt data using numerical methods in underground caverns [2, 21–23]. However, these studies have considered only the axial effect and neglected the transverse effect of bolts on the surrounding rock. Research on transverse deformation and force of bolts is mainly based on small-scale bolt shear tests and analytical methods. The mechanical characteristics of bolts considering axial and transverse deformation and their support effect on the surrounding rock should be studied further. Therefore, the current study takes underground caverns as research background, considers the axial and transverse deformation of bolts, and establishes the governing equations of load transfer between bolts and the surrounding rock. A new analytical algorithm for a fully grouted bolt, which is suitable for seismic time-history analysis under complex conditions, such as joints and fault, is proposed. The rationality of the algorithm is verified, and a 3D finite element method (FEM) model of underground caverns of a hydropower station is used to study the force and deformation regularity and support effect of bolts under seismic load.
2. Analytical Model for Bolt-Rock Interaction
In this study, the implicit bolt element, which is buried in rock mass element and has no effect on the division of rock mass is used for bolt simulation [24]. This study is mainly based on the following assumptions: (1) the bolt is made of a homogeneous elastic material; (2) the perfect combination between the bolt and mortar causes them to deform together; (3) the slip failure of the bolt only occurs at the interface between mortar and rock; (4) the transverse deformation of anchorage body is consistent with that of the rock.
2.1. Transverse Deformation and Force of Bolt
The mutual displacement between the joint surface of the rock causes the axial displacement and transverse displacement of the bolt, resulting in a similar deformation of the bolt. The mechanical model is shown in Figure 1 [25]. Point is the intersection of the bolt and the joint surface. The rotation angle at point is 0. The shear force and rotation angle are 0, and the bending moment is the largest at point [18]. is called plastic hinge segment, and is the plastic hinge point [17]. and are, respectively, the axial force and transverse shear force generated at point when the bolt is transversely deformed.
Segment in Figure 1 is simplified to a static beam, and the force analysis of this segment is conducted on the basis of classical beam theory, as shown in Figure 2. is the reaction force concentration of the rock on the compression side of the bolt, and its distribution along the bolt is parabolic [10], i.e.,where is the length of , and the expression of is semiempirical, [18]; is the yield strength of the bolt; and is the compressive strength of the surrounding rock. is the diameter of anchorage body, , where is the diameter of bolt, and is the thickness of mortar.
The shear force at point is 0; thus, the force balance analysis of is performed as
The bending moment of is
On the basis of the differential equation of beam deflection curve, the approximate deflection curve expression of can be obtained by the integral of equation (3) as
The rotation angle and deflection at point are 0; hence, equation (4) can be transformed as
Transverse displacement at point of the bolt can be obtained by substituting into equation (5) as
The transverse load-displacement-controlled equation between the bolt and rock can be obtained by combining equations (2) and (6) intowhere is the equivalent elastic modulus of the anchorage body; and are the elastic moduli of the bolt and mortar, respectively; , where , , and are the cross-sectional area of the anchorage body, bolt, and mortar, respectively; is the moment of inertia of the anchorage body, .
2.2. Axial Deformation and Stress of Bolt
2.2.1. Bond Slip Model for Interface
The trilinear shear slip model [26] shown in Figure 3 is adopted to describe the relationship between shear stress and shear displacement at the interface. Shear displacement is the relative displacement between the anchorage body and surrounding rock, i.e.,where and are the respective axial displacements of the surrounding rock and anchorage body at .
(1) Elastic Stage 1. The anchoring interface is in a linear-elastic state, and the expression iswhere is the shear stiffness of the interface during the elastic stage, which is calculated by the following equation [27]:where is the shear stiffness of the rock, which can be valued in accordance with different rock masses; is 5–10 and 1.5–3 GPa/m for hard and soft rock, respectively [1]. is the shear stiffness of mortar, which is calculated by the following equation [28]:where is the shear modulus of mortar; generally, .
(2) Plastic Softening Stage 2. Shear stress decreases when the relative displacement increases, and the expression iswhere is the shear stiffness of the interface during the plastic softening stage; is the peak shear stress; is the relative displacement corresponding to . , where and are the cohesive force and internal frictional angle of the interface, respectively; and is the hydrostatic confining pressure of the rock element perpendicular to the anchoring surface [22].
(3) Residual Strength Stage 3. The interface is completely damaged, and only residual shear stress exists, namely,
2.2.2. Load Transfer Differential Equation
A microsegment of the anchorage body, which has a coordinate of and length of , is considered for study. A 1D local coordinate system is established along the axial direction of the anchorage body and has a positive direction from the anchor head to the deep rock. The force and deformation of this anchorage body and rock are shown in Figure 4.
The balanced differential equation of the microsegment is
The following equation can be derived using equation (14):
Equations (9) and (12) are derived and, respectively, combined with equation (8), i.e.,where is the axial force of the anchorage body at ; is the shear stress of the anchoring interface; and are the strains corresponding to and , respectively.
In the transverse deformation of the bolt (shown in Figure 1), the axial displacement of the bolt in point is composed of the deformation caused by and , which can be derived from Hooke’s law and the formula of maximum bending moment normal stress [10], i.e.,
The microsegment shown in Figure 4 is supposed to be a section of a transverse deformation bolt, and corresponds to point in Figure 1. The following equation can be deduced from equation (18) via the geometric method:where is the bending section modulus of the anchorage body; ; is the transverse force of the bolt at ; ; is the transverse displacement of the bolt at .
The axial controlled equation for the interaction between the rock and anchorage body can be obtained by combining equations (16) and (17) with equations (15) and (19). The unified expression iswhere ; ; and . is replaced with in the elastic stage and with in the plastic softening stage.
3. Solution of the Controlled Equation
After the excavation of underground caverns, the displacement of the surrounding rock is a series of element node displacements in the directions of , , and ; the calculated and are discrete values. The functional expression of is difficult to obtain through the analytic method. Thus, the finite difference method is adopted to solve equation (20), thereby obtaining the numerical solution of .
Bolt length is denoted as and equally divided into segments along the axial direction. Each segment length is . The starting point of the bolt near the free surface is , and the end point is . The segmented points are numbered in order as 1, 2, …, . In accordance with the finite difference principle, the corresponding derivatives are replaced with the difference quotient of segmented points, and equation (20) is transformed aswhere , , and are the values of , , and at the -th segmented point, respectively. The range of in equation (21) is . When and , the following boundary conditions are needed:where is the prestress applied to the anchor head; if no prestress exists, .
and are taken, and equation (21) can be rewritten in a matrix form aswhere
Equation (23) is a tridiagonal linear equation system for element unknown variable , which can be solved using the Chase method. When the interface is in different stages, the values of parameters in equation (23) are different.(1)Elastic stage 1:(2)Plastic softening stage 2:(3)Residual strength stage 3:where and is the transverse displacement at the -th segmented point of the bolt.
The expression of can be obtained by transforming equation(14) into and are taken, and equation (28) is discretized by the finite difference method. The values of and at different segmented points can be calculated by
4. Dynamic Calculation Theory with Bolt Support
With regard to the dynamic finite element calculation of underground caverns, the differential equation of motion by the Lagrange method is established as [29].where and are the lumped mass and damping matrices of the model nodes, respectively; and are the acceleration and velocity of the nodes, respectively; and are the external and internal forces of the nodes, respectively.
4.1. Mechanism of Bolt Support
The mechanical restraint of the bolt to the surrounding rock is reflected by its support reaction force . The dynamic calculation differential equation under bolt support is obtained by substituting into equation (31) as
The axial shear stresses at the interface can be equivalently transformed into concentrated loads at the segmented points. The loads are exerted on the rock elements in an opposite direction, which can simulate the axial support effect of bolts. The shear stresses are supposed to be linearly distributed on the anchoring interface, and the equivalent load of each segmented point is expressed as
The three components of in the global coordinate system of the model are given aswhere , , and are the directional cosines of the bolt.
The plastic hinge segment at each segmented point must be calculated in accordance with the of the rock element where the segmented point is located. at the segmented point of bolt is calculated using equation (7), which can be exerted on rock elements in an opposite direction to simulate the transverse support effect of bolts. The three components of in the global coordinate system are given aswhere , , and are the directional cosines of transverse deformation at segmented points of the bolt.
In accordance with interpolation theory of shape function, the support reaction force exerted on the rock by the -th segmented point of the bolt is expressed aswhere is the value of the shape function at the -th node of the rock element at the segmented point of the bolt. and are the components of and in the directions of , , and in the global coordinate system, respectively.
4.2. Calculation Procedure
(1)The seismic response results with bolt support in the -th time step are obtained, and the free deformation of the surrounding rock without support in the -th time step can be obtained by solving equation (31). On the basis of rock displacement, transverse displacement and axial displacement at the segmented points of the bolt are obtained using interpolation theory of finite element shape function. The and of bolts are, respectively, calculated by equations (7), (23), and (29), and the support force in the transverse and axial directions is deduced.(2)The support force is exerted on the surrounding rock in the -th time step. The model is resolved with equation (32), the displacement of the surrounding rock in the -th time step with bolt supports is obtained, and this displacement is the real deformations of rock with bolts support during an earthquake.
The detailed calculation procedure is shown in Figure 5.
5. Verification of the Numerical Model of the Bolt
5.1. Example 1
Zhao et al. [30] conducted a bolt pullout test and studied the distribution of axial force and shear stress along the axial direction of the bolt under different pullout forces. In the three test groups, the tensile force is 80, 160, and 240 kN, respectively; the physical and mechanical parameters of the materials are shown in Table 1.
During the test, the bolt was not subjected to transverse force, and could be considered approximately equal to of the bolt. Therefore, equation (19) can be rewritten as
Equation (37) is combined with equation (20), and the differential equation is expressed as
In the test, the bolt was buried in the fixed concrete, i.e., . Therefore, equation (38) can be written in a matrix form aswhere
The anchoring interface is supposed to be in the elastic stage, and . On the basis of the parameters of rock, mortar, and bolt, can be calculated using equations (10) and (11); . The comparison of numerical values with the test values of the axial force and stress of the anchoring interface are shown in Figure 6. The numerical values of axial force and shear stress are nearly consistent with those of the test values, indicating that the bolt algorithm proposed in this study can accurately simulate bolt force in a pullout test.
(a)
(b)
5.2. Example 2
A transverse shear test of bolts in jointed rock mass [31] is conducted to study the mechanical behavior of bolts. The schematic of the test model is shown in Figure 7 in which the transverse shear force is 102, 128, and 154 N, respectively. The distribution of the axial strain of bolts is obtained by strain gauges 1–3, which are fixed on the bolt, and strain gauge 2 is close to the joint surface. The parameters of the materials are provided in Table 2.
The bolt is divided into four equal segments and five segmented points to calculate the axial strain of the bolt using the proposed algorithm. The positive direction is from strain gauge 1 to strain gauge 3, and the segmented point 3 is located on the joint surface. The value of plastic hinge section at segment point 3 is obtained. With the value of , the transverse displacement at the joint surface can be calculated using (7). Supposing the anchoring interface is in the elastic stage, the transverse displacement in segmented point 3 is and the transverse displacements at the other four segmented points are 0. The mortar is not considered in this test; therefore, the expression of can be rewritten as . On the basis of the material parameters, is taken, and the other parameters are consistent with the test.
Axial force – at the segmented points can be obtained by solving (23), and axial force and transverse shear force are substituted into (18) to obtain the axial displacements – at the segmented points. Axial strains – of the bolt at segmented points are obtained by geometric relations. The comparison of numerical values with the test values of the axial strain of the bolt is shown in Figure 8.
The distribution curves of the axial strain of the bolt, which are obtained from the test and calculation, show a trend of “large in the middle and small at both sides.” The closer the bolt is to the joint, the greater its axial strain. The axial strains of the bolt at the segmented points are obtained by the numerical method in this study and consistent with the strains measured by the test, indicating that the proposed algorithm can accurately simulate bolt deformation in a shear test.
6. Engineering Case Study
6.1. Engineering Profile
Figure 9 shows the FEM model of 3# of an underground powerhouse. The model is mainly composed of the main powerhouse, main transformer cavern, and surrounding rock. The sizes of the main powerhouse and main transformer cavern are 34.5 m × 30.0 m × 84.2 m (length × width × height) and 34.5 m × 18.5 m × 34.2 m, respectively. The model consists of 25,839 nodes and 22,710 hexahedron elements (eight nodes). The maximum mesh size is 10 m which meets the required dynamic calculation accuracy. The calculation ranges of the , , and directions are 289.3, 47.6, and 294.1 m, respectively. The buried depth of the underground powerhouse is approximately 150 m. The region where the project is located is earthquake prone, with a basic seismic intensity of VII. The initial geostress field is obtained via stress inversion of the measured points. The disturbing stress field after the excavation of the underground caverns is considered the initial condition for dynamic calculation. The blue elements are rock 1, the red elements are rock 2, and the green elements at the top arch of the cavern are fault. The middle section of the model is selected as the typical section, in which six monitoring points on the surrounding rock of the main powerhouse and main transfer cavern are taken as research objects (shown in Figure 10) to study the stress and deformation of bolts and the stress and displacement response of the surrounding rock.
(a)
(b)
6.2. Calculation Conditions
The elastic-plastic damage constitutive relation based on the Mohr–Coulomb criterion is adopted for surrounding rocks, and the scalar damage coefficient is expressed as [29].where is the damage coefficient of each principal stress direction, is the plastic bias strain in the -th principal strain direction, and is the material damage constant.
The bolts in the main caverns are arranged alternately in accordance with , . The physical and mechanical parameters of materials are shown in Table 3. A few bolts at the top arch of the caverns span two different rock mass materials, rock mass 1, and fault. The interfaces between all bolts and the surrounding rock are supposed to be , , , and to simplify the calculation. The input acceleration and displacement time history (shown in Figure 11) are synthesized in accordance with the seismic ground motion criterion of 5% exceedance probability in 50 years. The seismic wave is input from the bottom of the model. The direction seismic waves are shown in Figure 11, and the direction seismic waves are 2/3 that of the direction. The viscoelastic artificial boundary is applied to the bottom of the model, the top is the free surface, and the free field artificial boundary is applied at the other sides. The bolt dynamic algorithm proposed in this study is embedded into a dynamic FEM numerical simulation platform [29] to study seismic calculation with support bolts.
(a)
(b)
6.3. Seismic Response of Underground Caverns under Different Support Conditions
The influence of bolts on the antiseismic reinforcement of the surrounding rock is studied under three conditions: (1) without bolt, (2) considering only the axial deformation of bolt (the axial bolt algorithm), and (3) considering both the axial and transverse deformations of the bolt.
In condition (2), the transverse deformation of the bolt is not considered; thus, transverse force can be considered to be approximately 0, and can be considered approximately equal to the length of the bolt. Therefore, control equation (20) can be transformed into equation (38), which can be rewritten in a matrix form aswherein stage 1, ; in stage 2, . The other parameters are consistent with those in Section 3.
6.3.1. Analysis on the Axial Force of Bolts
The distribution regularities of the normal and shear stress of bolts A and B under conditions (2) and (3) conform to neutral point theory. The maximal normal stress appears at the neutral point, whereas the shear stress at both sides of the bolts reaches the maximum, as shown in Figure 12. Under condition (3), the maximal normal stress of bolt A increases considerably compared with that of condition (2), whereas the maximal normal stress of the bolt B increases slightly. Bolt A at fault has a large transverse deformation, whereas the transverse deformation of bolt B at the sidewall is small. The axial and transverse deformation of the bolt both contribute to the axial force for the algorithm in this study. Therefore, using the proposed algorithm can improve the support effect of bolts on the surrounding rock in the axial direction.
(a)
(b)
6.3.2. Analysis on the Transverse Force of Bolts
Using the bolt algorithm proposed in this study, Figure 13 shows the distribution regularity of the transverse deformation and force of bolts A, B, C, and D after an earthquake. Bolts A and D located at fault have large transverse deformation and force near the interface between the fault and rock 1, and bolts B and C at the sidewall have almost no transverse deformation and force. A transverse deformation of bolts can produce transverse force, which is the key to restraint dislocation between fault and rock mass. The bolts that pass through faults or joints undergo large transverse deformation at joints under seismic loads. Therefore, the axial bolt algorithm cannot fully reflect the support effect in a complex rock mass environment. The dynamic bolt algorithm, which considers both axial and transverse deformations, is necessary for practical engineering.
(a)
(b)
6.3.3. Displacement Response of Surrounding Rock
In this study, the relative displacement of monitoring points with the corresponding center points (E and F) of the floor on the typical section in the direction is taken to reflect the deformation law of surrounding rock. Figure 14(a) shows the relative displacement time history of points A and B with the corresponding center point E under three different conditions. The three time-history curves of point A are different in amplitude, whereas the three curves of point B are similar in waveform and wave distribution. The relative displacement peak of monitoring points (A, B, C, and D) is shown in Figure 14(b). The relative displacement peak of points A and D under condition (3) is lower than that under conditions (1) and (2), whereas the relative displacement peak of points B and C at the sidewall has minimal difference under the three conditions. The bolt algorithm has an inconsiderable influence on the numerical results of sidewall deformation, but considering the transverse deformation and force of bolts can effectively reduce the deformation of fault at the top arch.
(a)
(b)
6.3.4. Stress Response of Surrounding Rock
The tensile strength of rock is much smaller than its compressive strength; thus, tension fracture is a common failure form of rock material. The maximum principal stress time histories of point A under three conditions are taken for analysis (as shown in Figure 15(a)). The maximum principal stresses are mainly compressive stresses, and the fluctuation amplitudes of these curves are small in the period of 0–1.8 s. From 1.8 s to 7.2 s, the curves fluctuate violently, with their fluctuation ranging from −0.46 MPa to 0.82 MPa. In the period of 7.2–20 s, the fluctuation tends to be steady. A comparison of Figures 11(a) and 15(a) shows that the time when the maximum principal stress at point A occurs under the three conditions is consistent with the time when maximum acceleration occurs. The maximum principal stress time-history curves are similar to the input acceleration curve, indicating that the stress response of the surrounding rock is affected by the input seismic wave. The stress responses are greatly different for point A under three conditions. The fluctuations of maximum principal stress are most intense under condition (1), and the fluctuations are lowest under condition (3). The effects of bolt support on the peak values of maximum principal stresses are shown in Figure 15(b). The peak values of the maximum principal stress decrease considerably when considering bolt support. The peaks of the maximum principal stress of the top arches (points A and D) under condition (2) are larger than those under condition (3), whereas the peaks of the maximum principal stress of the sidewalls (points B and C) are not much different under conditions (2) and (3). Hence, the stress response of fault can be reduced effectively by using the bolt algorithm proposed in this study.
(a)
(b)
6.3.5. Damage Coefficient of Surrounding Rock
The damage coefficient of surrounding rock can directly reflect its seismic damage degree. The distributions of damage coefficient on the typical section under three conditions after the earthquake are plotted in Figure 16. The damage zone without bolt support is mainly distributed at the top arches and sidewalls of the main powerhouse and main transformer cavern and the intersection of the main powerhouse with tailrace tunnel, diversion tunnel, and main electrical wire hall (as shown in Figure 16(a)). The damage coefficient for majority of the area is approximately 0.5. Damage is most serious at the top of the caverns, which are crossed by fault, and the maximum damage coefficient exceeds 0.9. When considering the axial support of bolts, the damage area and degree of surrounding rock are considerably reduced, while the damage coefficient of sidewalls and the intersections of caverns are reduced to 0.3. However, the damage coefficient of the top arches remains large; the maximum of which is approximately 0.8 (as shown in Figure 16(b)). Compared with Figure 16(b), Figure 16(c) depicts that the damage degree, which is obtained by the proposed bolt algorithm, at the top arches of the caverns is considerably reduced. The damage coefficient for majority of the top arches is approximately 0.3, and the maximum damage coefficient does not exceed 0.5. Therefore, the transverse support of bolts can effectively reduce the damage degree of the surrounding rock at the fault under seismic load.
(a)
(b)
(c)
7. Conclusions
On the basis of classical beam theory and the trilinear shear slip model, a new dynamic algorithm that considers the axial and transverse deformations of bolts is proposed. The algorithm is applied to the seismic support analysis of an underground project, and the following conclusions are drawn:(1)The proposed dynamic algorithm fully considers the characteristic of slippage and yield of interface between the anchorage body and surrounding rock and the transverse deformation characteristic of bolts. The algorithm can effectively simulate the stress and deformation of bolts located at faults or joints under seismic load. The rationality and feasibility of the dynamic algorithm are verified through comparison using pullout data and a shear test of bolts.(2)The distribution of the axial normal and shear stresses of bolts obtained by the proposed algorithm is similar to the results calculated by the axial bolt algorithm; both of which fit in with neutral point theory. However, the proposed algorithm is more effective in restraining the surrounding rock in the axial direction of the bolts due to the contribution of the transverse deformation to the axial force. The bolts located at the fault have large transverse deformation, whereas the bolts located at the rock without joint or fault have no transverse deformation. This result indicates that the transverse deformation of bolts in complex rock environment conditions during an earthquake should be considered.(3)Bolt support can effectively reduce the directional deformation of the fault at the top arches of the caverns, and exerts a minimal effect on the sidewalls. The directional deformation of the top arches under earthquake calculated by the proposed bolt algorithm is lower than that calculated by the axial bolt algorithm, which indicates that the transverse force of bolts affects the restriction of the mutual slip of rock masses at the fault. After considering bolt support, the maximum principal stress of the surrounding rock is considerably reduced under seismic load, and the reduction extent of the maximum principal stress of the fault calculated by the proposed bolt algorithm is greater than that calculated by the axial bolt algorithm. Considering bolt support can effectively reduce the damage degree of the entire underground caverns, and the reduction extent of damage for the surrounding rock at the fault calculated by the proposed bolt algorithm is more remarkable than that calculated by the axial bolt algorithm.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This study was supported by the National Key Basic Research Program of China (973 program) (Grant no. 2015CB057904) and the Major Program of the National Natural Science Foundation of China (Grant nos. 51579191 and 51279136). These supports are greatly acknowledged and appreciated.