Abstract
Open-pit slopes contain numerous nonpenetrating, intermittent joints which maintain stability under blasting operations. The tip dynamic response coefficient (DRC) of parallel cracks in a typical rock mass under combined dynamic and static loading conditions was calculated in this study based on the superposition principle. The dynamic response law of the intermittent joint in the slope under blasting was determined accordingly. The influence of many factors (the disturbance amplitude of dynamic load, the lateral confining pressure, the length of rock bridge, the length between cracks, the staggered distance between cracks, and the crack inclination angle) on the dynamic response was theoretically analyzed as well. The ABAQUS numerical assessments were conducted on simulation models with parallel cracks under combined dynamic and static loading conditions. The results show that a larger dynamic load amplitude and smaller crack inclination angle/confining pressure result in greater Type II dynamic strengthening effect on the crack tip. When the length of the rock bridge between cracks () is smaller than the half length of the crack (), the dynamic strengthening effect at the crack tip weakens gradually with increase in ; when, the strengthening effect is almost unchanged. With the increase in the staggered distance between cracks (), the dynamic strengthening effect of the crack tip weakens at first and then strengthens; the strengthening effect is weakest when ; the crack propagation under combined dynamic and static loading is the most sensitive to the lateral confining pressure () and is the least sensitive to the inclination angle of the cracks (). Theoretical results are validated by comparison with numerical simulation results. Such information regarding the dynamic response law of the parallel cracks in rock masses under dynamic and static loading conditions is conducive to further research on the mesofailure mechanism of open-pit mine jointed rock slopes under blasting operations.
1. Introduction
The fracture failure characteristics of the deep rock mass in an open-pit mine slope are different from those of the shallow rock mass under blasting conditions. Deep rock mass failure is the result of the joint action of high ground stress and explosion load. The role of ground stress, in particular, is very important [1–3]. There are a large number of discontinuous joints in rock slope, so ground stress may play a crucial important role in controlling the strength of the rock mass and the failure mode of the rock slope. Many scholars have explored the mesofailure mechanism of rock masses with joints. Bohloli and De Pater [4], for example, studied the cracking and propagation mechanism of cracks in rock masses under different surrounding stress conditions and found that confining stress has strong influence on the fracturing behaviour of the compacted sand. Tang et al. [5] analyzed the mechanism of crack generation and propagation in a rock mass with intermittent joints by using the RFPA 2D software and demonstrated that the crack growth is stable and stops at some finite crack length with a confining pressure. Zhou et al. [6] performed dynamic laboratory tests on 3D-printed artificial rock samples and found that wing cracks can continuously extend to the sample ends under dynamic compression, whereas they can intermittently extend only a limited distance under static compression. Tang et al. [7, 8] studied the brittle fracture of a rock sample with an angled crack under combined tensile and compressive loading conditions using linear elastic fracture mechanics (LEFM) and found that the high compressive confining pressure and substantial friction can increase the possibility of shear crack extension. Li et al. [9] analyzed the influences of preexisting flaw inclination angle on the cracking processes by means of finite element method (FEM) and nonlinear dynamics method and found that the crack initiation sequences and the overall crack pattern are different under different loading conditions.
There have been many other valuable contributions to the literature. Zhou and Chen [10] studied the parallel crack interactions in rock under unloading conditions based on fracture mechanics and found that a smaller length of rock bridge or staggered distance between the cracks results in more severe and sensitive interactions at the crack tip. Zou et al. [11] performed the research of cracking processes in a natural rock under dynamic loading conditions and found that only the antiwing and shear patches evolve into two symmetrical pairs of shear cracks, which result in the specimen failure under dynamic compression. Li et al. [12] analyzed the effects of preexisting flaws with different flaw angles and lengths on the dynamic mechanical properties and found that the presence of an artificial flaw may change the failure mode of marble under dynamic loading. Zhu et al. [13] performed uniaxial compression to characterize and visualize the mechanical and fracture properties of the 3DP rock and demonstrated that the initial internal voids and cracks dominate the spatial fracture evolution and failure patterns within the rock.
However, there have been few previous studies on the propagation and evolution of parallel cracks and the dynamic response between cracks under dynamic and static loading conditions. Additionally, the research methods on the propagation law of cracks under combined dynamic and static loading are mainly focused on rock specimen testing and numerical simulation while the relevant theoretical quantitative researches are yet lacking. There is no fully feasible, effective, straightforward methods for calculating the stress intensity factor (SIF) at the crack tip of rock masses with multiple cracks. The dynamic response of parallel cracks under the combined action of dynamic and static loads was investigated in this study based on the superposition principle. Variations in the self-defined dynamic response coefficient (DRC) of the crack tip were investigated in detail as per the disturbance amplitude of dynamic load, the lateral confining pressure, the length of rock bridge, the length between cracks, the staggered distance between cracks, and the crack inclination angle. A numerical dynamic/static load crack propagation simulation was conducted in ABAQUS, and the results were compared against the theoretical results to validate them.
2. Stress Intensity Factor (SIF) Calculations
2.1. Fracture Mechanics Models
Three failure forms may occur in the rock mass under the combined action of in situ stress and explosion stress: shear failure, tensile crack failure, and compression-shear failure, where the failure is dominated by both I-II tension-shear mixed mode failure (when the stress acted on the crack surface is tensile stress) and I-II compression-shear mixed mode failure (when the stress acted on the crack surface is compression stress) along the structural plane [14, 15]. Tensile failure of the structural plane is generally caused by a tensile stress wave formed by the explosion stress wave reflected by the free surface. The structural plane in the deep rock mass of the slope is usually closed under compression due to the substantial in situ stress, which leads to compression-shear failure of the internal structural plane under in situ stress and explosion stress. Therefore, structural-plane tensile failure is not discussed here; this study focused on shear failure of the rock mass along the structural plane.
Figure 1 shows several closed parallel nonpenetrating intermittent joints in the slope rock mass. A “key block” with two parallel intermittent joints in the slope was selected as the research object here for the sake of simplicity.
As shown in Figure 1, due to the propagation of blasting stress waves, the key block is subjected not only to the static load of vertical and horizontal in situ stress ( and ) in the near area but also to the horizontal dynamic load () under blasting.
According to the superposition principle of fracture mechanics theory, the total stress field caused by two or more different loading systems near the crack tip can be obtained by the algebra of each SIF on the basis of linear elastic mechanics under the same boundary conditions [16, 17]. Therefore, the SIF of the crack tip in the rock mass stress State A () under the combined action of dynamic and static loads (Figure 2) can be considered equivalent to the sum of the static SIF () and the dynamic SIF (). where and are the SIFs at the crack tip under biaxial static loading and dynamic loading alone, respectively.
2.2. Calculation of SIFs () under Biaxial Compression
The key block in the jointed rock mass of the slope can be regarded as an infinite plate with parallel double cracks [18]. As shown in Figure 3, the vertical action of the plate has a uniformly distributed load, and the transverse action has a uniformly distributed load . The center of the plate contains two parallel cracks; the crack inclination angle is . The lengths of Crack 1 and Crack 2 are both . is the distance between the tips of the two parallel cracks perpendicular to the direction of the cracks, which is referred to here as the “staggered distance” of the cracks. is the distance between the tips of the cracks along the direction of the cracks, which is referred to as the “length of the rock bridge between the cracks.”
According to the stress circle theory in material mechanics, the stress State B of the rock mass falls under biaxial compression [19]. The stress State B under biaxial compression can be equivalent to the case of vertical stress and shear stress acting on the crack surface, as shown in the stress State B1 in Figure 3.
According to the superposition principle, the stress State B1 of the rock mass under biaxial compression can then be decomposed into the stress State B2 and the stress State B3 (Figure 4) [17].
According to the stress-circle theory, the vertical stress and shear stress acting on the crack surface in Figure 3 can be expressed as follows [20, 21]:
For the State B2 with vertical stress, the cracks are closed under compression, and there is only SIFs of Type II () at the crack tips, which can be calculated according to Formula (3). For the State B3 with shear stress, there is also only SIFs of Type II () at the crack tips, which can be calculated according to Formula (4) [22]. where and are SIFs of Type II at the crack tips in the stress States B2 and B3, respectively; and are the geometric correlation coefficients of SIFs related to the length of the cracks (), the length of the rock bridge, and the staggered distance between the cracks, which can reflect the interaction between multiple cracks and can be determined by consulting the appendix or reference [22]. Therefore, according to the superposition sketch shown in Figure 4, the SIFs at the crack tips under biaxial compression can be expressed as follows: where and are SIFs of Type I and Type II at the crack tips in the stress State B, respectively.
2.3. Calculation of SIFs () under Dynamic Loading
The stress State C of the rock mass under dynamic loading can be considered equivalent to the stress State C1, as shown in Figure 5. The stress State C1 can be considered equivalent to the superposition of the stress States C2 and C3, as shown in Figure 6 (); that is, . According to the approximate method proposed by Broberg in analyzing the problem of dynamic crack propagation in an infinite elastic body, the dynamic SIFs can be expressed as a function of static SIFs [23, 24]: where is the dynamic SIF when the crack growth rate is , is the static value of , and is the velocity influence coefficient independent of the geometric shape of the crack. Once has been determined, the dynamic SIF calculation problem can be approximately regarded as a quasistatic SIF problem, that is
According to references [24, 25], the general forms of and can be expressed as follows: where is the velocity of the Rayleigh wave in the rock mass, which is related to the density , elastic modulus , and Poisson’s ratio of the rock mass. can be calculated as follows [25]:
The dynamic SIFs at the crack tips can be obtained as based on Formulas (7) and (8) as follows (assuming that the crack is closed):
It is assumed that the blasting dynamic load on the key block is a simple harmonic wave in the form of a sinusoidal function (Figure 7). The simple harmonic wave can be expressed as where is the amplitude, is the frequency, and is the duration of the disturbance wave; is the number of longitudinal waves and , where is the velocity of the longitudinal wave. Let the simple harmonic wave be in the form of an index, i.e., the imaginary part of . , so the stress intensity of the incident wave in its propagation direction is as follows [26]: where is the amplitude of the disturbance stress intensity, is the velocity of the shear wave, and is the Lamé constant. The disturbance stress of the blasting dynamic load on the end of the “key block” can be expressed as follows:
2.4. Calculation of SIFs () under Dynamic and Static Loading
Substituting Formulas (5) and (10) into Formula (1) yields the following expression of the SIFs at the crack tips under combined dynamic and static loading conditions: where and are Type I and Type II SIFs at the crack tips under dynamic and static loading conditions, respectively.
3. Dynamic Response Analysis
The effects of dynamic loading among multiple cracks on the SIF at the crack tip may be reinforcing, null, or shielding [27–29]. To observe the effects of dynamic load on crack propagation, the ratio of the SIF at the crack tip under combined dynamic and static loading to the SIF at the crack tip under biaxial static loading was defined in this study as the dynamic response coefficient (DRC) : where is the DRC of Type II at the crack tip, which reflects the effect of dynamic loading on crack initiation and propagation at the crack tip. and are Type II SIFs at the crack tips under combined dynamic and static loading condition and biaxial static loading condition, respectively. In addition, the dynamic load has three effects on the crack initiation and propagation of the crack tip: strengthening, null, and shielding. The crack falls into the reinforcing effect zone (strengthening), the null effect zone (null), and the shielding effect zone (shielding) when , , and , respectively.
3.1. Effect of Dynamic Load Amplitude ()
Due to space constraints, only the inner tip of Crack 1 (position (1) in Figure 3) is reported here as an example (similarly below) to analyze the DRC variation of Type II with the dynamic load amplitude (). In order to simplify the calculation, let , , , ( is the half lengths of Crack 1 and Crack 2), , and . According to reference [22], when , and . When , and. Substituting and into Formula (13) yields the dynamic load function () varying with time. And then, the DRCs of Type II () at different times can be obtained by Formula (15). The DRC () calculation results and variation curves of the inner tip of Crack 1 with dynamic load amplitude () obtained are shown in Table 1 and Figure 8.
(a)
(b)
As shown in Figure 8, when the length of the rock bridge between the cracks is constant, the Type II DRCs () of the inner tip of Crack 1 vary continuously as the dynamic load action time increases. The DRCs are greater when the dynamic load amplitude () is larger. The propagation effect of dynamic load on crack tip appears to constantly change over time. The strengthening effect of Type II on the crack tip increases as the amplitude increases, that is, the initiation and propagation of the crack tip is more significantly affected by the dynamic load, and the initiation is more likely to occur at the crack tip in this case.
3.2. Effect of Confining Pressure ()
To observe the DRC variations of Type II of the inner tip of Crack 1 with confining pressure (), let , , , , , and . According to reference [22], when , and . When , and . The calculation results and variation curves of the DRCs of the inner tip of Crack 1 with confining pressure () are shown in Table 2 and Figure 9.
(a)
(b)
As shown in Figure 9, the Type II DRCs of the inner tip of Crack 1 increase with the continuous decrease in confining pressure (); the dynamic strengthening effect on the crack tip becomes greater as decreases. Crack tips in this case are more easily initiated under dynamic loading when confining pressure is high. The existence of confining pressure increases the initiation strength of the crack surface, which hinders the propagation of the crack and makes the crack surface less breakable. Increase in confining pressure also increases the crack initiation strength and weakens the effect of the dynamic load.
The above analysis indicates that as the confining pressure of the key block located far away from the slope surface (“Y” in Figure 1) is larger than that of the key block located closer to the slope surface (“J” in Figure 1), the nonpenetrating joints located in the deep zone of the actual slope are more significantly affected by the dynamic load than those located in the shallow zone. Under these conditions, the crack tip is more likely to crack during blasting.
3.3. Effect of the Rock Bridge Length () and the Staggered Distance () between Cracks
To observe the DRC variations of Type II of the inner tip of Crack 1 with the length ratio of the rock bridge () and staggered distance ratio (), let , , , , , and . The DRC calculation results and variation curves of the inner tip of Crack 1 with and are shown in Table 3 and Figure 10.
(a)
(b)
As shown in Figure 10(a), when the staggered distance () between cracks is constant, the DRCs of Type II at the inner tip of Crack 1 decrease sharply at first and then remain basically unchanged with further increase in the length ratio of the rock bridge between cracks (). When the length of rock bridge between cracks is small (), the dynamic strengthening effect at the crack tip is substantial; as increases, however, the dynamic strengthening effect weakens gradually. When , the dynamic strengthening effect is basically unchanged, and the dynamic response is almost unaffected by the length of the rock bridge between the cracks.
As shown in Figure 10(b), when the rock bridge length between cracks is constant, the DRCs of Type II at the tip of Crack 1 decrease at first and then increase with increase in the staggered distance ratio (). The minimum value is reached when , indicating that the dynamic response strengthening effect of crack tip weakens at first and then strengthens as the staggered distance between cracks increases; the dynamic response strengthening effect is the weakest when .
It appears that the staggered distance between the nonpenetrating joints and the length of the intermediate rock bridge in the actual slope have an important influence on the initiation of the crack tip and the penetration between the cracks during blasting operations. When the length of the rock bridge between the cracks is small and the staggered distance between cracks is large, the crack tip is more strongly affected by the dynamic load and cracks are more likely in general.
3.4. Effect of Crack Inclination Angle (α)
To observe the DRC variations of Type II of the inner tip of Crack 1 with the inclination angle (α), let , , , , , and . The DRC calculation results and variation curves of the inner tip of Crack 1 with inclination angle (α) are shown in Table 4 and Figure 11.
(a)
(b)
As shown in Figure 11, when the length of the rock bridge () between cracks is constant, the DRCs of Type II at the inner tip of Crack 1 decrease gradually as the inclination angle (α) increases. This suggests that the dynamic response strengthening effect of crack tip weakens, and the crack tip is less likely to break under dynamic load as higher inclination angles.
The above analysis indicates that a smaller inclination angle of the nonpenetrating joint in the slope makes the crack tip more significantly affected by the dynamic load, in which case cracks are generally more likely to occur.
4. Numerical Simulation and Analysis
4.1. Numerical Models
To validate the theoretical analyses presented above, numerical simulations of 20 groups of calculation models (Table 1) under dynamic and static loading conditions were carried out in the ABAQUS software. As shown in Figure 12, the models are in size with 20,618 grid elements and 35,324 nodes in total (due to space constraints, only portions of the models are shown here). Elastic modulus () of the model material is 4.2 GPa and Poisson’s ratio () is 0.26. The finite element models adopted the plane stress eight-node quadrilateral elements, and the local refinement is carried out in the area around the cracks. In addition, the singular elements are arranged in the area of the crack tips considering the strong singularity of the crack tips after the crack surfaces are closed, and the cracks are arranged in the middle of the models to prevent the boundary effect from influencing the crack failure. The geometric parameters of the cracks in the models are listed in Table 5.
(a)
(b)
(c)
(d)
4.2. Simulation Scheme
In the simulation, considering the effect of the original rock stress on the crack propagation in the slope, the static load was preapplied to simulate the initial stress field before the dynamic load was applied to the model. The multiload step technique in ABAQUS was used for this purpose. The stress application sequence and path are shown in Figure 13. First, axial pressure () and lateral confining pressure () were loaded to at a rate of 0.1 MPa/s simultaneously under the loading control mode. Then, the lateral confining pressure () was kept constant as axial pressure was loaded to 2 MPa at a rate of 0.1 MPa/s, and the SIFs at the tips of Crack 2 in this state (State A) were calculated accordingly. Finally, the dynamic load in the form of a sinusoidal incident wave with amplitude of 0.5 MPa and frequency of 0.5 Hz was input from the left side of the model until the model was destroyed (according to the material parameters of the models, after many debugging calculations, when the amplitude of the dynamic load is 0.5 MPa, and the frequency is 0.5 HZ, the calculation results can meet the requirements). The SIFs at the inner tip of Crack 2 in each model at different times were also calculated at this time.
4.3. Analysis of the Simulation Results
4.3.1. Crack Propagation Process
Due to the space constraints, only the initiation and propagation processes of cracks in Models 1, 5, 13, and 15 under combined dynamic and static loading are shown here (Figures 14–17).
As shown in Figures 14–17, cracks appear in the inner and outer tips of Crack 1 and Crack 2 and propagate continuously as the dynamic loading process continues. The crack propagation rates in Model 13 are much faster than in Model 1, Model 5, or Model 15. The inner tip of Crack 2 in Model 13 is the first to crack (), followed by that in Model 1 (), Model 5 (), and finally Model 15 (). To this effect, the rock bridge lengths between cracks and the crack inclination angle both exert an important influence on the initiation and propagation of cracks in the model under dynamic loading. When the length of the rock bridge and the inclination angle are relatively small, cracks occur more easily under the dynamic load and propagate at a faster rate. The dynamic strengthening effect of the crack tip is higher in this case, which supports the conclusions discussed in Sections 3.3 and 3.4.
4.3.2. Variation Analysis of Dynamic Response and Dynamic Initiation SIF (DIS)
The SIFs (, ) of the inner tip of Crack 2 in the dynamic initiation state of different models (Figure 12) calculated by ABAQUS are as follows: 0.078 MPa m1/2 (, Model 1), 9.16 MPa m1/2 (, Model 1), 0.092 MPa m1/2 (, Model 5), 9.71 MPa m1/2 (, Model 5), 0.061 MPa m1/2 (, Model 13), 8.95 MPa m1/2 (, Model 13), 0.12 MPa m1/2 (, Model 15), and 9.89 MPa m1/2 (, Model 15). During the initiation of Crack 2, the SIFs of Type II were much larger than the SIFs of Type I. Therefore, the crack failure in the model under dynamic and static loading conditions can be mainly attributed to Type II failure, and the effect of Type II SIFs on crack propagation is separate from the effects of Type I SIFs, as the crack in the model was closed under confining pressure.
In order to validate the dynamic response law of crack propagation obtained by theoretical analysis (Sections 3.1–3.4), the ratio of the dynamic SIF () of the inner tip of Crack 2 at any time () in the dynamic loading process (as calculated in ABAQUS) to the static SIF () of the inner tip of Crack 2 when the axial pressure and the lateral static confining pressure both reach the preset value predynamic-loading was defined as the DRC. . The calculated DRC results when and and the dynamic initiation SIFs of the inner tip of Crack 2 in different models are shown in Table 6 and Figure 18.
(a)
(b)
(c)
(d)
Figures 18(a), 18(c), and 18(d) show that the DRCs of the tip of Crack 2 decrease gradually with increase in rock bridge length, crack inclination angle, and lateral confining pressure when and . The dynamic initiation SIFs increase with increase in rock bridge length, crack inclination angle, and lateral confining pressure as well. To this effect, these three factors all have important effects on crack initiation in the tip. The dynamic response of the crack tip is more pronounced under dynamic loading when these three values are smaller, i.e., cracks occur more easily in the tip.
As shown in Figure 18(b), as the staggered distance between cracks increases, the DRCs of the tip of Crack 2 decrease at first, then increase, and finally tend toward stability. The dynamic initiation SIFs at the tip of Crack 2 increase at first and then decrease in this case. The DRC is minimal and the dynamic initiation SIF is maximal when is around 4 mm, at which point the dynamic response of the crack tip is the weakest under the dynamic load, and the tip is less likely to crack. The dynamic response increases as the staggered distance between cracks increases. When reaches about 10 mm, the DRC of the crack tip is almost unchanged with any further increase in staggered distance; in other words, the dynamic response of the crack tip is no longer sensitive to changes in fault distance once the staggered distance reaches . This numerically simulated dynamic response law is consistent with the theoretical law discussed in Sections 3.1–3.4.
4.3.3. Sensitivity Analysis of Various Factors
Suppose the number of elements of the comparison column and reference column are both . And, and can be expressed as , respectively. Write and in the form of a matrix, as shown in Formula (16) and Formula (17) [30]. where represents the value of the -th variable and represents the value of the dependent variable corresponding to the -th variable.
The interval relative valuing is used to dedimensionalize the matrices of and . And, the calculation formulas are as shown in Formula (18) and Formula (19):
The elements in the dedimensionalized matrices of and are subtracted and taken as absolute values for differentiation, as shown in Formula (20).
Then, the calculation Formula of the correlation coefficient can be expressed as follows: where is the resolution coefficient, which is used to increase the significance of the correlation coefficient of each factor, . In general, the value of is 0.5. Generally, a large number of correlation coefficients are calculated, and the information is relatively scattered, which is not conducive to the correlation degree analysis. Therefore, , the average value of, is taken as the average correlation coefficient to analyze the correlation degree among the factors. where . The closer the average correlation coefficient is to 0, the smaller the correlation degree of the reference column to the comparison column is, and the lower the sensitivity of the reference column is to the comparison column is; the closer the correlation coefficient is to 1, the greater the average correlation degree of the reference column to the comparison column is, and the higher the sensitivity of the reference column to the comparison column is.
Select the change values of the influencing factors in Table 5 and the DRCs of State B and State C in Table 6 as the comparison matrix and the reference matrix, respectively. The average correlation coefficients (and ) of State B and State C can be obtained by substituting them into Formulas (18)–(21):
According to the Grey correlation theory, the sensitivity of DRC to influencing factors is in the following order: . It appears that the crack propagation under combined dynamic and static loading is the most sensitive to the lateral confining pressure (), followed by the length of rock bridge () and the staggered distance (), and the crack propagation is the least sensitive to the inclination angle of the cracks ().
5. Conclusions
The dynamic response laws of parallel cracks under combined dynamic and static loading conditions were analyzed in this study based on the multiple superposition principle and SIF calculations under fracture mechanics theory. A series of ABAQUS numerical simulations were performed to validate the theoretical analysis results. The conclusions can be summarized as follows. (1)The dynamic load amplitude, confining pressure, and crack inclination angle all markedly influence the crack dynamic response law. A larger load amplitude, smaller inclination angle, and smaller confining pressure result in greater dynamic strengthening effect of Type II on the crack tip. That is, crack initiation and propagation in the tip are significantly affected by dynamic loading(2)The length of the rock bridge and the staggered distance between cracks also significantly influence the crack dynamic response law. When the length of the rock bridge between cracks is small (), the dynamic strengthening effect at the crack tip is large. It weakens gradually; however, when , the dynamic strengthening effect remains almost unchanged, and the dynamic response is practically unaffected by the length of the rock bridge between the cracks. The dynamic response strengthening effect is weakest when (3)The crack propagation under combined dynamic and static loading is the most sensitive to the lateral confining pressure (), followed by the length of rock bridge () and the staggered distance () and is the least sensitive to the inclination angle of the cracks ()(4)The ABAQUS simulation produced results in close agreement with the theoretical analysis results, which validates the dynamic response laws of parallel double cracks under combined dynamic and static loading conditions reported in this paper.
Appendix
As shown in Figures 4 and 6, the geometric correlation coefficients of SIFs in Equations (3), (4), (5), (14), and (15) can be obtained as shown in Figures 19 and 20 [22].
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 they have no conflicts of interest.
Acknowledgments
This research was supported by the National Key R&D Programs of China (grant nos. 2016YFC0801602 and 2017YFC1503103).