#### Abstract

Using a self-designed hydraulic impact drilling test-bed and rock core drill, six groups of cylindrical granite specimens (93 mm dia. × 200 mm) containing central axial holes formed either by impact or nonimpact drilling methods were tested in uniaxial compression to failure on an Instron 1346 universal testing machine to investigate their mechanics and damage properties. The longitudinal acoustic wave velocities were measured before testing. The rock specimens were grouped according to the method of drilling the central hole (impact load exerted by different impact power and different frequencies for an approximately identical impact power, or nonimpact drilling). In this study, a statistical constitutive damage model based on Weibull distribution was used to calculate the degree of rock damage after drilling center holes. The experimental curves were measured to analyze the damage evolution process and the radius of rock damage. These indicate that rock damage increased with the increase of impact power and decreased with increasing impact frequency at constant impact power. This was also verified by the measured longitudinal wave velocity in all rock specimens. These results have significance for guiding the design of composite rock drilling tools that are dedicated to improving rock-breaking efficiency.

#### 1. Introduction

The search for techniques of enhancing rock damage to increase rock-breaking efficiency is a key issue in geotechnical engineering activities, such as mining, tunneling, and general underground excavation, and also in research and development of new rock drills and tunneling equipment. Many wide-ranging studies in China and globally have been conducted on these issues, leading to the rapid development of a focus on rock-damage mechanics in geotechnical engineering [1, 2].

The rock-breaking method of drilling and blasting is widely used in rock excavation around the world because of its low cost, low blasting energy, and ready acceptance by blasting engineers [3–5]. In this method, a high-frequency and high-velocity hammer imparts a force to the rock upon impact by generating a stress wave sufficient to cause the rock to fracture or be damaged. Studies of the impact load on rock mechanics properties have examined the effect of different strain rates [6–8]. Others have investigated the effects of strain rate and confining pressure on rock strength, elastic modulus, and Poisson’s ratio [9–14]. Others have proposed a strength theory that comprehensively takes into account the effects of strain rate and confining pressure on rock strength [15–17]. Yang [18] studied the deformation, peak strength, and crack damage behavior of hollow sandstone specimens with different hole diameters subjected to different confining pressures in triaxial compression experiments. Yang et al. [19] studied jointed specimens to assess the effect of the jointing on peak strength and Young’s modulus, using the damage mechanics method. Li et al. [20] studied the damage evolution process in porous sandstone under loading and proposed a multiscale damage model, which showed that the strain rate and confining pressure had an effect on the strength and mechanical properties of the rock to some extent.

However, during the excavation process, the rock drill causes damage to the rock due to the impact load that it produces. Therefore, it is of great significance to carry out in-depth research on the rock damage extent and mechanical properties under impact loading.

In the process of tunneling, for example, the rocks surrounding the tunnel sustains a certain amount of damage from the impact loads imposed by the process of blasting and the redistribution of stress after excavation [21–24]. Su et al. [25] studied the damage evolution of double free faces subjected to different dynamic disturbances. Saini and Shafei [26] investigated the accuracy of four constitutive models in practical applications to assess the performance of a concrete structure under low-velocity impact loading. Wang et al. [27] established a multiparameter classification combining stress, energy, and damage to predict deterioration in a model of rock failure. Jiang et al. [28] and Wang et al. [29] established a constitutive model of the damage to a rock mass subjected to the impact load of blasting and investigated the damage to the rock surrounding an excavation on the basis of damage mechanics and other related theories. Based on the continuous damage theory, Xu et al. [30] established a thermomechanical coupling damage constitutive model of rock using the Hoek–Brown strength criterion and further investigated the relationship between the confining pressure and damage evolution rate. Golshani et al. and Meglis et al. studied the excavation damage zone and the rock damage surrounding tunnels using on-site ultrasonic testing [31–34], borehole camera technology, and numerical analysis software, among other techniques. Huang et al. [35] used finite element modeling to investigate the excavation-damaged zone using the established damage model. Souissi et al. [36] used numerical simulation analysis to explore the damage zone around the drill bit during the hard-rock drilling process. Zhao et al. and Su et al. investigated crack propagation in rock specimens containing a drilled hole and subjected to a compressive load [37, 38] and explored the way rock strength varied with radial stress using an acoustic emission method. Huang and Yang [39] analyzed the influence of flaw angle and confining pressure on the mechanics properties of rock when subjected to triaxial compression.

Gu et al. [40] proposed a damage constitutive model that takes the crack compaction effect into account based on accumulated acoustic emission (AE) counts. Jin et al. [41] proposed a damage constitutive model of coal-rock using AE experiments during coal-rock uniaxial compression. Ma and An and Hamdi et al. analyzed rock fracturing and damage caused by blasting using numerical simulation software and a numerical algorithm based on image analysis [42, 43]. Wang et al. studied the influence of confining pressure on the failure mode of surrounding rock using a numerical analysis method [44]. With the aid of Irazu geomechanical simulation software, Vazaios et al. simulated brittle behavior and spall fracture mechanisms in a hard, blocky rock mass using the finite-discrete element method [45].

As seen in the literature, various research methods—theoretical analysis, numerical simulation, and laboratory tests—have been reported for the study of rock damage under different loading conditions; however, a quantitative analysis of the rock damage range under different impact loads still needs further exploration. At present, relatively few studies have been conducted on the influence of rock drill impact power and frequency on rock damage, so research in this field will provide important reference values for practical engineering contexts.

Rock damage induced by a combination of dynamic and static loading in the process of tunneling was investigated by Zou et al. and Tao et al. [46–48], using the split Hopkinson pressure bar (SHPB) apparatus to study the compressive strength and elastic modulus of rock specimens with boreholes, and the properties of the damage induced in the surrounding rock. Improved SHPB apparatus was used in studies of the dynamic and static strength properties and dynamic fracture behavior of sandstone and other sedimentary rocks [49–52]. Taheri et al. [53] studied the mechanical properties of brown coal under different loading conditions to assess mine stability. Wang et al. [54, 55] investigated the creep properties of salt rock subjected to triaxial compression. Wu et al. [56] studied the nonlinear characteristics of salt rock at the creep stage, from which they proposed a constitutive model of creep damage in salt rock, and verified the validity of the model.

According to the modern mining manual [57], the blast-hole spacing in open pit mining operations is generally 30 × blast-hole diameter. To increase the spacing and therefore reduce the number of blast holes to save mining energy consumption, it is feasible to increase the blast-hole diameter. Based on this assumption, the principle of composite rock drilling shown in Figure 1 is proposed. To determine the appropriate diameter, drilling experiments on rock specimens in this study adopted granite as the rock material and high-frequency hydraulic impact rock drilling. By combining theoretical and experimental results, the damage evolution and radius of the damage in the surrounding granite when subjected to a range of high-frequency impact loads were analyzed to provide a scientifically based guide for the design of a composite rock drill (Figure 1).

#### 2. Study of Rock Damage Zone due to Impact Loading

##### 2.1. Analysis of Stress in Rock Induced by Impact Loading

The stress distribution in rock subjected to the impact of the drill bit was considered to approximate Boussinnesq’s solution [58]. Figure 2 shows a penetration force *P* acting in the plane *Z* = 0. By using the method of elasticity mechanics, the components of the stress due to *P* are obtained in cylindrical coordinates.

The stress components are given bywhere *σ*_{r}, *σ*_{θ}, and *σ*_{z} are the radial normal stress, circumferential normal stress, and the axial normal stress, respectively; *P* is the penetration force; *r* is the radius; *Z* is the depth of the rock unit to the point of penetration force; *μ* is Poisson’s ratio; *τ*_{rz} is the radial shear stress; and *τ*_{zr} is the axial shear stress.

By applying impact dynamics [59] to the penetration analysis of impact drilling system composed of piston, drill pipe, and rock block, the penetration differential equation can be obtained as below:where *t* is time; *K* is the penetration coefficient (N/cm); *m* is the wave resistance (kg/s); and *p* is the incident wave (MPa). According to the impact dynamics, when the piston strikes the drill rod, *p* is given bywhere *V* is the terminal velocity of the piston when it strikes impacts the drill rod and *M* is the mass of the piston. By substituting equation (3) into equation (2) and then integrating it, *P* is calculated fromwhere *γ* = *m*^{2}/*MK* is the impact penetration index (dimensionless). According the experimental conditions and selected drilling tools [60], *γ* = 0.18. Applying the maximum principle of differential equations to equation (4), the maximum *P*_{m} is obtained from

For the hydraulic rock drill, the relationship between input and output power is expressed aswhere *P*_{oil} and *Q*_{oil} is the pressure (MPa) and the flow rate (L/min) of hydraulic oil input of the hydraulic rock drill; *η*_{hy} is the working efficiency of the hydraulic rock drill, generally *η*_{hy} = 0.55; *f* is the impact frequency of the piston (Hz); and is the single impact energy of the piston (*J*). Combining equations (5) and (6), the maximum impact load of the hydraulic rock drill during impact drilling, *P*_{m}, is obtained from

##### 2.2. Zone Distribution of the Surrounding Rock and Analysis of Elastic and Plastic Zone Distribution

If the stress in the rock surrounding the drill hole within a certain range exceeds the ultimate rock strength, it will sustain damage, causing cracks or deformations. Once this has occurred, it cannot withstand any further excessive stress. Any such area (e.g., area A in Figure 3) is termed as the “plastic zone”. Beyond the plastic zone, an elastic zone still exists (e.g., area B in Figure 3) and can withstand greater stress. Beyond the elastic zone is the original rock zone (e.g., area C in Figure 3).

As shown in Figure 2, *Z* = 0 during the process of drilling; therefore, equation (1) gives the radial stress *σ*_{r} at the central hole wall under impact load *P*_{m} as

In the plastic zone (*r*_{0} < *r* ≤ *r*_{1}), the rock complies with the Mohr–Coulomb strength criterion for plastic conditions, expressed aswhere *c* and *φ* are the cohesion and internal friction angle; *r*_{1} is the radius in the plastic zone; *σ*_{rp} is the radial normal stress; and *σ*_{θp} is the shear normal stress.

As shown in Figure 4, take one rock element in the plastic zone for stress analysis; according to the equilibrium principle of radial stress and shear stress in a rock unit, equation (10) can be obtained:

Equation (10) is the stress equilibrium equation for axisymmetrical conditions. Combining equation (9) with equation (10), according to boundary condition , the radial normal stress in the plastic zone, *σ*_{rp}, is given by

And the shear normal stress in the plastic zone, *σ*_{θp}, is expressed by

Assume *σ*_{0} is the initial rock stress and *σ*_{r1} is the radial normal stress at the elastic/plastic zone boundary, the radial normal stress *σ*_{re} and the circumferential normal stress *σ*_{θe} in the elastic zone of the surrounding rock (*r* ≥ *r*_{1}) are expressed by

From the above equations, the sum of the stresses in the plastic zone and the elastic zone are obtained, respectively, as follows:

The sum of radial stress and shear stress in the elastic zone is equal to that in the plastic zone at the elastic/plastic zone boundary [61]; thus, the implicit equation for the radius of the plastic zone boundary *r*_{1} is

In impact drilling, the existence of unrecoverable strain also indicates that the rock has sustained a certain amount of damage. Elastoplasticity theory gives the strain in the plastic zone aswhere *ε*_{r} and *ε*_{θ} are the radial and tangential strains in the plastic zone; *u*_{r} and *u*_{θ} are the radial and tangential deformations in the plastic zone; *f*_{a} is the plasticity function (*f*_{a} = 0 in the elastic deformation stage); and *E*, *μ*, and *G* are the elastic modulus, Poisson’s ratio, and shear modulus, respectively.

From equation (16), the volumetric strain is derived as

Integrating equation (17) gives the total deformation in the plastic zone:

Using cylindrical cavity expansion theory to define the initial damage to the surrounding rock due to the impact load, the following is obtained:where *R* is the external diameter of the rock specimen.

#### 3. Rock Constitutive Equation for High-Frequency Impact Load

##### 3.1. Basic Assumptions

(1)A viscous body is not damaged and does no work under approximately static load conditions. The constitutive relationship iswhere *σ*_{b} is the stress of the viscous body and *η* is the viscosity coefficient determined by time-dependent deformation (creep) tests. The properties of the viscous body are reflected by *η*, where 0.1 ≤ *η* ≤ 0.5.(2)The strength of the mass of every microunit of damaged rock obeys the Weibull distribution [62–64]. The probability density function iswhere *F* is the distribution variable of the microunit strength and *m* and *F*_{0} are the Weibull distribution parameters that reflect the mechanical properties of the rock materials.(3)The mass of damaged rock is macroscopically isotropic. It is linearly elastic before damage and loses bearing capacity after rock failure.

##### 3.2. Statistical Damage Variable

Supposing that fractures of rock specimens occur as a result of constant failure of the microunits and that the total number of microunits of the rock specimens is *N* and the number of fractured microunits under certain load *F* is *N*_{f}, the statistical damage variable *D* is defined as

Thus, from equation (22), *N*_{f} is derived:

Combining equations (21)–(23), the statistical damage variable *D* on the basis of Weibull distribution can be deduced as follows:

Equation (24) shows that the damage variable is closely associated with rock microunit stress influenced by the stress state. To illustrate the impact of a complex stress state on rock strength, the failure criterion of rock is taken into consideration. Suppose the general failure criterion of each microunit is

Because the Drucker–Prager yield criterion has a simple parameter form and simultaneously reflects the influences of the volumetric stress, shear stress, and intermediate principal stress on rock strength, it is more widely used than other strength theories to reflect the actual situation. The strength of a rock element is given bywhere *φ* is the angle of internal friction; *I*_{1} is the first invariant of the stress tensor; and *J*_{2} is the second invariant of the stress deviator. Hooke’s law gives the following:

When there is no confining pressure (i.e., the one-dimensional stress state), the minimum principal stress *σ*_{2} _{=} *σ*_{3} _{=} 0 and the strain *ε*_{1} = *ε*. Combining equations (26) to (29), *F* is given by

##### 3.3. Rock Constitutive Model

As the microunit model is a parallel combination, its stress and strain are equal to those suffered by two separate units. Considering the damage *D*_{0} caused by drilling the center hole using high-frequency impact loading and that the strain rate in a uniaxial compression strength test is close to zero under approximately static loading, then the strain equivalence hypothesis gives the constitutive relationship of the rock damage mass as

##### 3.4. Calculation of Model Parameters

The model parameters *m* and *F*_{0} in equation (31) must be further determined to establish the damage constitutive model of the rock specimens. On the basis of the extremum of a multivariate function, the slope is zero at the vertex of the stress-strain curve, that is, the coordinates of the peak of the stress-strain curve are (*σ*_{m}, *ε*_{m}), mathematically expressed as

Substituting the peak value into equation (31) and combining it with equation (32), the damage model parameters *m* and *F*_{0} are

#### 4. Experimental Verification

##### 4.1. Specimen Preparation

The rock specimens for the experiments were selected from identical granite with good integrity and homogeneity and with no obvious macroscopic cracks on the surface. Six rock blocks measuring 42 cm × 28 cm × 30 cm were prepared for the experiments. Parallel 30 mm diameter holes were drilled in five of the blocks (i.e., A, B, C, D, and E) by impact drilling on a purpose-designed and manufactured experiment test bed.

The implementation process of impact drilling is illustrated in Figure 5. At selected frequencies, the impact piston of the hydraulic rock drill strikes the shank adapter with a reciprocal motion (Figure 5(a)), transferring its kinetic energy in the form of a stress wave to the rock by way of the shank adapter, drill pipe, and drill bit, causing the rock to fracture and fragment, forming the hole, and also damaging the rock surrounding the hole. The rotational mechanism causes the drill pipe and drill bit to rotate through the shank adapter, grinding the fractured and fragmented rock into dust that was later removed by compressed air.

**(a)**

**(b)**

**(c)**

**(d)**

The hydraulic rock drill (the core component of the impact drilling test bed) was produced by Dazao Machinery Co., Ltd (Guangzhou, China). Its impact performance parameters are listed in Table 1.

In order to determine the optimum impact parameters, the spare rock block had been impacted by gradually increasing the oil pressure before the holes were drilled. When the oil pressure reached about 10 MPa, the rock had already been fractured. As a result, the maximum oil pressure was set as 7 MPa for this experiment. As hydraulic oil is very sensitive to changes of temperature, and since oil spills inevitably occur in the hydraulic system, it is difficult to maintain the oil pressure at a fixed value. As a result, it is normal in practice that both oil pressure and flow fluctuate within a certain range during actual operation. In the process of drilling the holes in the five experimental rock blocks (i.e., rock block A, B, C, D, and E), the pressure and flow data for the hydraulic rock drill were recorded using a hand-held Hydrotechnik tester (shown in Figure 6) in order to calculate the input power of the hydraulic rock drill at any given time.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

For the final block (block F), a coring drill was used to drill the holes using a nonimpact method, a rotating diamond coring bit to leave the hole after the core was removed. The process is shown in Figure 7.

**(a)**

**(b)**

**(c)**

**(d)**

Figure 8 shows the five rock blocks A, B, C, D, and E that were drilled with center holes using the impact method and the rock block F that was drilled center holes using a nonimpact method. Six groups of rock specimens were then prepared in accordance with the international experimental performance-testing requirements for rock mechanics. The final rock specimens are shown in Figure 8(c).

**(a)**

**(b)**

**(c)**

Because varying degrees of damage was caused to the rock by the impact load at different power settings, the longitudinal wave velocity also varied. To thoroughly explore the damage properties of rock caused by impact load and to correspondingly verify the results of the uniaxial compressive experiments, a HS-YS4A rock acoustic parameter tester was used to measure the longitudinal wave velocity at three positions in two randomly chosen rock specimens in each group before the uniaxial compressive experiment. Then, a servocontrolled MTS 815 triaxial testing system and Instron 1346 universal testing machine were used to measure the physical and mechanical properties of the granite under static load. The test results are listed in Tables 2 and 3.

##### 4.2. Experimental Equipment and Principle

An Instron 1346 universal testing machine, Figure 9 maximum load 2000 kN (static loading) and 1000 kN (dynamic loading); measurement accuracy ± 0.5%) at the Institute of Mechanics Properties, Advanced Research Center, Central South University, Changsha, China, was used to conduct the uniaxial compressive tests on the rock specimens attached with strain gauges. Through the experiment, the stress-strain curve, Poisson’s ratio, and elastic modulus of each rock specimen were obtained.

**(a)**

**(b)**

The physical parameters and related impact drilling parameters are listed in Table 4.

##### 4.3. Experimental Results

The stress-strain curve of every rock specimen was plotted. Figure 10 shows the stress-strain curves for each specimen in each of the six rock specimen groups; the average peak stress of each group is noted on the figures.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

The following findings are based on the stress-strain curves for each rock specimen group shown in Figure 10.

In the initial stage of loading, the stress-strain curves of all the rock specimens almost coincide, that is, each rock specimen has the same initial elastic modulus. Then, the slope of each curve gradually increases as the initial cracks inside each rock specimen start to close, leading to an increase in the elastic modulus. As the peak stress is approaching subsequently, the slope of each curve gradually decreases, which indicates that continued loading contributed to the expansion of cracks and larger deformation inside each rock specimen in the last stage, consequently causing a decrease of elastic modulus and the deterioration of material properties.

The average peak stress in the stress-strain curve of rock specimen group F was clearly higher than for the other five rock specimen groups, notably almost 12% higher than group E. Because the center holes in group F specimens were drilled by a nonimpact method, almost no damage was caused to the surrounding rock. Conversely, the other five specimen groups underwent impact drilling, whose radial stress component caused damage to the surrounding rock and reduced its compressive strength.

The stress-strain curves for groups A, B, and C, with increased impact power of the hydraulic rock drill, that the compressive strength of the rock was on the decline. Equation (7) states that the output impact load *P*_{m} of the drill is proportional to the impact power (i.e., the product of pressure × flow) to the ½ power. In addition, equation (1) states that the radial stress at the wall of the center hole, *σ*_{r}, is proportional to the impact load. Therefore, with the increase of input power, the radial stress on the wall of the center hole increased and the damage to the rock surrounding the center hole also increased, thus lowering the overall compressive strength, as observed.

Although rock specimen groups C, D, and E experienced approximately the same impact power, the compressive strength of the rock increased with higher impact frequency, implying that less damage occurs to the surrounding rock at higher impact frequencies and (approximately) constant input power from the rock drill.

This observation is connected to the intrinsic properties of a hydraulic rock drill, that is, each single impact load (energy) produced by high-frequency impact is less than for low-frequency impacts at any given power. This indicates that the parameters of the damage to the surrounding rock are independent of drill impact frequency, but rather they are determined by impact load. However, drilling speed is closely related to frequency in open-cut mining and tunneling.

In addition, as Figure 11 shows, with increasing impact power, both the uniaxial compressive strength and the longitudinal wave velocity of rock specimens display a similar tendency to decline. This is because the greater impact power results in greater expansion of internal cracks and damage in the rock surrounding the center hole, weakening the uniaxial compressive strength, and increasing the propagation time of the longitudinal waves (i.e., lower wave velocity). Conversely, according to Figure 12, with the increase of impact frequency, both the uniaxial compressive strength and the longitudinal wave velocity were on the rise under basically the same power. This occurred because, with the increase of impact frequency, the force of single impact of the hydraulic rock drill diminished, leading to decreased damage to the surrounding rock and an increase in both uniaxial compressive strength and longitudinal wave velocity.

**(a)**

**(b)**

**(a)**

**(b)**

##### 4.4. Comparison between Experimental and Theoretical Constitutive Relations of Rock under High-Frequency Impact Load

From the damage constitutive model in equation (32) and the physical and mechanical parameters of granite in Table 3, calculation and analysis of the experimental data gave the fitted parameters of the constitutive model of six rock specimens selected randomly from each group (i.e., rock specimens A1, B5, C1, D3, E5, and F6), as shown in Table 5, and the theoretical curves of the damage constitutive model were calculated.

Figure 13 presents a comparison between the six experimental stress-strain curves derived from the uniaxial compressive tests and the theoretical curves derived from the damage constitutive model for those rock specimens.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

According to Figure 13, the experimental curve for each rock specimen is consistent with its theoretical curve, which to some extent shows the rationality of the proposed damage model and reflects well the relationship between the strength and the strain of rock under impact load.

However, during the closure of internal cracks and the consequent linearly elastic stage, the stress value in the theoretical curve is larger than that in the experimental curve. The proposed damage constitutive model fails to reflect the slowly increasing stress at that stage. The theoretical curve implies that the rock is linearly elastic up to peak stress; thus, there is some error in modeling the initial stage of the curve, and it indicates that the model needs further improvement.

##### 4.5. Properties of Damage in Surrounding Rock under High-Frequency Impact Load

According to the abovementioned equations, it is deduced that the surrounding rock damage *D*_{0} is a function of the impact load *P*_{m}, center hole diameter *a*, and radius in the plastic zone, *r*_{1}. The ratio *r*_{1}: *a* is defined in this study as the damage ratio, denoted by *β*. The damage *D*_{0} and damage ratio *β* were calculated from the physical parameters and experimental data.

Figure 14 shows that *D*_{0} and *β* both tend to increase with increasing impact power, whereas for constant impact power, *D*_{0} and *β* tend to decrease with increasing impact frequency.

**(a)**

**(b)**

#### 5. Conclusion

The experimental results indicate that the uniaxial compressive strength diminished in specimens whose central hole had been drilled with greater impact power, and the amount of damage to the rock surrounding the hole was seen to have increased.

At any given impact power, an increase in the impact frequency was associated with an increase in the uniaxial compressive strength of the rock specimen, that is, there was evidently less damage to the rock surrounding the center hole as the impact frequency was increased. Analysis of the longitudinal wave velocity measurements was consistent with the above experimental results.

It was inferred from the experimental data that, of all the selected impact drilling parameters, impact power (impact load) had a decisive effect on the surrounding rock damage, and that impact frequency appeared to affect only the drilling rate. However, further investigation is still needed on how a further increase of impact frequency might affect surrounding rock damage, or even if the impact frequency were to exceed the natural frequency of the rock.

On the basis of analysis and calculation of the experimental data, the most appropriate parameters for a proposed rock damage constitutive model were obtained. Comparisons between the theoretical results of the model and the experimental results verified that the experimental curve was in good accordance with the theoretical curve. Therefore, the constitutive equation can be applied to predict the stress-strain relationship of rock subjected to high-frequency impact loading and has a good potential for engineering application.

#### Data Availability

The data used to support the findings of this study are include within the article.

#### Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

#### Acknowledgments

The research presented in this paper was supported by the National Natural Science Foundation of China (Grant no. 51774340).