#### Abstract

In recent years, terrorist attacks and regional conflicts have increased and research on the resistance of building structures has been more highlighted. Q235B is a commonly used steel type in building structures, and there is an urgent need for research on its dynamic material properties. In this paper, Paul constitutive model was calibrated based on the material properties of Q235B steel. Based on the obtained results, it was found that the original model was not able to take into account temperature-softening effect. Therefore, a temperature-softening parameter was added into the temperature term of the original Paul model, and the modified M-Paul constitutive model was developed. Then, numerical prediction results obtained from Q235B steel penetration experiments using J-C (Johnson–Cook) and M-Paul constitutive models were compared. It was observed that J-C constitutive model significantly overestimated ballistic limit while estimations obtained from M-Paul constitutive model were closer to experimental results. This observation verified the correctness and effectiveness of M-Paul model and established a solid foundation for relevant research works on the impact resistance of building structures.

#### 1. Introduction

In recent years, terrorist attacks and regional conflicts have increased and research on the resistance of building structures has been more highlighted. Under dynamic loadings, the stress-strain relationships of materials become more complicated. Therefore, the selection of a reasonable and effective material constitutive relationship plays an important role in obtaining accurate simulation results. The Johnson–Cook (J-C) constitutive model [1], which takes into account the effects of strain hardening, strain rate hardening, and temperature softening, has been successfully applied for various metal materials. The Zerilli–Armstrong (Z-A) constitutive model [2] has been employed for different metal materials such as FCC and BCC, but it was not applicable for the investigation of the mechanical behaviors of metal materials at high temperatures (*T* > 0.6 Tm) and low strain rates. The Khan–Liang (K-L) constitutive model [3] was found to be more suitable for the evaluation of the mechanical behaviors of several BCC metal materials with high strains (15%), strain rates (10^{−6}–10^{−4} s^{−1}), and temperatures (77–600°F). The Lin–Wagoner (L-W) constitutive model [4] was successfully applied for gapless steel and 310 stainless steel and took into account strain rate and temperature effects. The Rusinek–Klepaczko (R-K) constitutive model [5] was suitable for evaluating the mechanical behaviors of low-carbon steels with high strain rates (10^{−4} − 5 × 10^{3} s^{−1}) and temperatures (213–393 K). Yang [6] developed a new elastic-peak plastic-softening-fracture constitutive model for borehole expansion problems. The Paul [7] constitutive model successfully characterized the mechanical behaviors of metal materials with high strain rates and temperatures. It was shown that this model was more applicable for DP600 grade steel than J-C, L-W, and K-L models at high strain rates.

Yang [8] investigated the sensitivity of the strain rate of S690 structural steel and developed a new Cowper–Symonds model to predict the dynamic increase factor and used the J-C model to fit true stress-strain curves. Amitava [9] predicted the thermal deformation behavior of aluminum 5083 by using J-C, Z-A, and strain-compensated Arrhenius models. Forni [10] employed J-C and Cowper–Symonds models to evaluate the strain rate behavior of S355 structural steel. Lin [11, 12] focused on the mechanical properties of Q235B steel with different strain rates, temperatures, and stress triaxialities and modified softening term in the J-C model. Xiao [13, 14] conducted Taylor impact tests on aluminum alloy 2024-T351 and studied the effect of Lode parameters on its fracture behavior. Li [15] also carried out Taylor impact tests on aluminum alloy 6061-T6511H and performed numerical simulations using Abaqus software based on the modified J-C constitutive model, modified J-C fracture criterion, and Lou fracture criterion. Chen [16, 17] conducted Taylor impact and target plate penetration tests to verify the accuracy of the J-C material model.

It can be seen in the above paragraphs that significant research works have been conducted on the constitutive relationships of materials and a large number of effective constitutive relations have been widely used in different engineering fields. Previous research on the properties of metal materials has been mostly based on J-C and modified J-C constitutive models. But, it is not entirely suitable for expressing the dynamic mechanical behaviors of structural steel and some better constitutive models should be found necessarily. At present, research on the dynamic mechanical behaviors of materials mostly adopts a combination of experiments and numerical simulations. Numerical simulations are frequently used, and their accuracies are verified by Taylor impact or target penetration tests. They are also commonly used to verify the accuracies of constitutive models and fracture criterion.

In this paper, the Paul constitutive model was calibrated by several material tests. It was found that the model could not take into account temperature-softening effect for the Q235B steel test specimen. Therefore, temperature-softening parameter was added into the temperature term of the model, and the modified M-Paul constitutive model was proposed. Moreover, it was observed in the test that flat-body projectile-penetrated Q235B steel target plate and the damage form as well as the ballistic limit of the target plate are obtained, which are subjected to the projectile penetration. The numerical model was developed using Abaqus software, and M-Paul, modified J-C, and J-C constitutive models were used in the simulation. It was observed that the J-C constitutive model significantly overestimated ballistic limit, while estimations obtained from the M-Paul constitutive model were closer to experimental results, and the accuracy of the proposed material model was verified.

#### 2. M-Paul Constitutive Model and Parameters Calibration

The original Paul [7] constitutive model considered the effect of plastic strain, strain rates, and temperature on the flow stress, and the homologous temperature was taken into account for easy calculation; the expression is given aswhere *B*, *C*, and *β* are the strain-hardening coefficients, *A* and *H* are the strain rate sensitivity coefficients, and *k* and *G* are the thermal-softening coefficients. *σ*_{0} is the yield strength and *ε*_{eq} is the equivalent plastic strain. The dimensionless plastic strain rate defined as , is the reference strain rate, and is the strain rate. And, , where *T* is the material temperature, *T*_{r} is a reference temperature, and *T*_{m} is the melting point of the considered material. Parameters *A* and *k* represent the yield strength item, parameters *B*, *C*, and *β* represent the strain hardening item, and parameters *H* and *G* represent the plastic deformation item.

In this paper, the mechanical behavior of Q235B steel was described by the Paul constitutive model; therefore, the calibration of unknown material parameters of the Paul constitutive model depends on all kinds of material test results in literature [12].

##### 2.1. Calibration of Quasistatic Item

The corresponding quasistatic item of the Paul constitutive model is shown in Equation (2); the model represents the material mechanical behavior at the reference strain rate (quasistatic) and temperature (room temperature). In this paper, the parameters B, C, and *β* were determined by fitting equivalent plastic stress-equivalent plastic strain curve obtained by tensile tests on the round bar specimens:

Because the conversion formula from engineering stress to true stress is no longer applicable after necking [14], the true stress-true strain curve before necking is used to fit using Origin to obtain *B* *=* 583.1 MPa, *C* *=* 454.3 MPa, and *β* *=* 1.

After that, the numerical model was generated by using Abaqus/Standard, and user-defined material mechanical behavior (UMAT) of the Paul model was developed. Then, the initial value of *B*, *C*, and *β* was entered in the UMAT, and the finite element simulation of the smooth round bar tensile test was performed. The load-displacement curves of finite element simulation were output and compared with the test results, as shown in Figure 1. It is noticeable that the yield plateau of finite element simulation load-displacement curve is not significant, so the test curve after yield point will be compared to the finite element simulation result.

In order to obtain the more accurate parameter value, this paper refers to the numerical optimization process in literature [18], namely, the Abaqus and Excel components were integrated in Isight5.9–2 to optimize the parameters. The Isight optimization gives the following parameters: *B* *=* 210.1 MPa, *C* *=* 740.7 MPa, and *β* *=* 1.409. The load-displacement curve is obtained through running the Abaqus finite element model again with the above parameter, as shown in Figure 1. By comparing the three curves, it can be seen that the optimized Isight curve containing the postnecking data points is more similar to the test curve compared with the Origin fitting curve which only use the prenecking data. Therefore, it is more reasonable to adopt the parameters optimized by Isight, so the quasistatic parameters of the Paul constitutive model is *B* *=* 210.1 MPa, *C* *=* 740.7 MPa, and *β* *=* 1.409.

##### 2.2. Calibration of Strain Rate Item

At reference temperature (room temperature), the Paul constitutive model is expressed as

At the yield point (*ε*_{eq} = 0), the relationship between yield strength and strain rates is shown as follows

The coefficient *A* can be determined by Equation (4) or fitting test data. The yield stress of Q235B steel at reference temperature and different strain rates is summarized in Figure 2. Equation (4) is used to fit the test data giving *A* = 0.0317.

At reference temperature, the Paul constitutive model can be expressed aswhere and . Equation (5) is equal to the quasistatic item equation after the parameters *σ*_{0}, *B*, *C*, and *β* were magnified *D*_{A}, *D*_{H}, *D*_{H}, and one times, respectively. In this paper, the Abaqus and Excel components were integrated in Isight to obtain coefficient *H*. In order to reduce workload and computational time, the tensile test at different strain rates was simplified as the quasistatic tensile test to be simulated. Only the parameters are appropriately enlarged, Abaqus/Standard module was still used to solve the problem. The Isight optimization flow is shown in Figure 3.

The optimization programs are as follows: First, coefficient *H* was imported into Excel-1, Excel-2, Excel-3, and Excel-4 to get different strain rates corresponding to *D*_{A} and *D*_{H}, respectively, and *B*, *C*, and *β* were obtained by multiplying the corresponding amplification factor. Second, the *B*, *C*, and *β* were input to four Abaqus component INP files to carry out the numerical simulation. Third, the relevant load and displacement data were output from Abaqus component to Excel-5, Excel-6, Excel-7, and Excel-8, respectively. At the same time, ten character points of test data were input to the above Excel to calculate four objective function values and target parameters, respectively. At last, four target parameters were input to the Excel components to calculate their mean square error and the final target parameter was obtained and input into the optimization component. Then, the optimization algorithm is selected, and the bounds of target parameter and the minimum value of objective function are set. By Isight optimization, *H* = 0.0409.

##### 2.3. Modified Expression and Calibration of Temperature Coefficient

The Paul constitutive model at the reference strain rate is shown as follows:

At the yield point (*ε*_{eq} = 0), the relationship between yield strength and temperature is shown as

The coefficient of *k* can be determined by Equation (7) or by fitting test data. The tensile test yield stress of Q235B steel at the reference strain rate and different temperatures is summarized in Figure 4. Equation (7) is used to fit the test data to get *k* = 2.071. According to Figure 4, it is hard to predict thermal-softening behavior of Q235B steel by the Paul constitutive model, and in order to improve the constitutive model prediction accuracy of material mechanical behavior under different temperatures, the relative temperature items of the Paul constitutive model were modified in this paper. The modified expression is expressed as follows:where *m* is the thermal-softening coefficient in the modified expression. Equation (17) is used to fit the test data to get *k* = 11.066 and *m* = 2.925, as shown in Figure 4. And, the constitutive model expression and coefficient are shown in Table 1.

Similar to the temperature item of parameter *k*, the temperature item of the Paul constitutive model where parameter *G* was located was slightly modified and then was modified as . And, the modified Paul constitutive model (M-Paul constitutive model) is expressed as

At reference strain rates, the M-Paul constitutive model can be expressed aswhere and .

Equation (10) is equal to quasistatic formula in which parameters , *B*, *C*, and *β* were magnified *D*_{k}, *D*_{G}, *D*_{G}, and one times. The above way that Abaqus and Excel components were integrated in Isight was adopted and then *G* = 1.977 and *n* = 1.515 were used as the initial value to be optimized. Finally, the optimization gives *G* = 1.989 and *n* = 1.575. And, the related material parameter of Q235B steel is shown in Table 2.

##### 2.4. Comparative Analysis of Different Constitutive Models Prediction

At reference temperature, the modified J-C constitutive model is similar to the J-C constitutive model. Therefore, the J-C constitutive model parameters of literature [12] and the above Paul constitutive model parameters were put into Equations (11) and (3), respectively. Then, the flow curves of Q235B steel at different strain rates were obtained, as shown in Figure 5:

**(a)**

**(b)**

**(c)**

**(d)**

By comparing the flow curves of Q235B steel at different strain rates, it is found that compared with the J-C constitutive model, the Paul constitutive model predicts the test better and could better reflect the mechanical behavior of Q235B steel at different strain rates. Especially at a high strain rate, the Paul model is more accurate in predicting the test equivalent stress-equivalent plastic strain curve than at the low strain rate, which can better simulate the mechanical behavior of Q235B steel at a high strain rate.

At reference strain rates, the flow curves of Q235B steel at different temperatures can be obtained by introducing the parameters of the modified J-C constitutive model and the M-Paul constitutive model calibrated above into equations (12) and (13) respectively, as shown in Figure 6:

Since the extensometer used in the test cannot monitor the deformation of the specimen at high temperature, the displacement used in processing the data of the tensile test at high temperature is the beam displacement of the test machine [12]. Compared with the displacement of the gauge section output by the extensometer, there are some deviations in the data, which will result in a certain deviation between the flow curve of the constitutive model and the test curve. However, from the prediction of ES low-carbon steel by the Paul constitutive model in literature [7], the Paul model can better characterize the dynamic mechanical properties of low-carbon steel materials. By comparing the flow curves of Q235B steel at different temperatures, it is found that compared with the modified J-C constitutive model, the modified Paul constitutive model can predict the test relatively well and can better reflect the mechanical behavior of Q235B steel at different temperatures.

#### 3. Ballistic Test of Q235B Steel

Q235B target plate material is a 70 mm diameter rod produced by Shanghai Baosteel, which penetrates into a 6 mm thick target plate by the wire-electrode cutting and milling machine, and 12 round holes whose diameters are 3 mm are uniformly processed on the circle 31 mm from the center. The projectile material is hardened 9CrSi, and the average diameter and length of the cylindrical blunt projectile are 5.95 mm and 29.82 mm, as shown in Figure 7.

**(a)**

**(b)**

The tests were performed at a one-state compressed gas gun installed at School of Civil Engineering in Nanyang Institute of Technology; the initial and residual velocities of projectile were obtained by FASTCAM SA-Z high-speed camera. The initial velocity of projectiles is controlled by pressured gas.

The total 14 ballistic tests of Q235B steel target plate were completed. Among them, there were 5 tests in which the trajectory of the projectile was not horizontal before penetrating the target plate, so the number of effective tests is 9, and the results of the effective test are shown in Table3. Among them, B-1 represents the first test of Q235B steel subjected to penetration by blunt projectile and later the same.

As Table 3 shows, the minimum initial velocity of the projectile penetrating the target plate is 283.8 m/s and the maximum impenetrable velocity is 272.48 m/s, so the ballistic limit velocity is between 272.48 m/s and 283.8 m/s. The initial-residual velocity curve of the projectile penetrating the target plate can be fitted using formula (R-I) proposed by Recht and Ipson [19], shown as follows:

In the expression, *V*_{i} is the projectile initial velocity, *V*_{r} is the projectile’s residual velocity after penetrating the target plates, *V*_{bl} is the ballistic limit velocity, and *a* and *p* are the model parameters. According to Equation (14), the method of least square is used to fit the test data to get *V*_{bl} = 280.88 m/s, *a* = 0.53, and *p* = 3.11; the fitting results are shown in Figure 8.

In the two tests that the projectile did not penetrate the target, the projectile was all embedded in the targets, but the results of the two tests were slightly different. In test B-1 of low initial velocity (*V*_{i} = 234.51 m/s) of projectile, the target plate failure model is rear bulge. However, in the test B-8 with high speed (*V*_{i} = 272.48 m*/*s), an obvious bulge found in the back of the target plate is uplifted and basically perforated and the plug was basically formed. Because shear stress is the major contribution in the whole process of deformation, the target plate failure model is shear plugging, as shown in Table 4.

Figure 9 shows the whole process of a target under projectile impact at *V*_{i} = 272.48 m*/*s by using a high-speed camera. The forming process of the plug is marked in red box. After the projectile impact of the steel target plate, the plug is not separated from the target and no real plug is formed.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

In the 7 tests when the initial velocity of the projectile is higher than the ballistic limit, the projectile penetrates the target plate and the failure model is the shear plugging. Moreover, different degrees of bulge occurred near the perforation on the back of the target plate, as shown in Table 4. The work done by plastic deformation was converted into heat, the heat was not out of the dissemination timely, and the projectile and the surface of plug generate high temperature. The special shear plugging failure is named as adiabatic shear failure. Through observation, it was found that a large part of the materials in the punching hole were blue because the high temperature oxidized the materials in this area, which further indicated that the penetration test of the target plate could be used to study the mechanical behavior of the materials at high strain rate and high temperature.

When the initial velocity of projectile is slightly higher than the ballistic limit, the back concave crater shape of the target plate is quite regular. With the improvement of the projectile initial velocity, the uneven arch height began to appear near the perforation on the back of the target plate.

When *V*_{i} = 422.91 m*/*s, the maximum and minimum bulge heights on the back of target plate are 1.93 mm and 3.11 mm, respectively. Therefore, with the increase of the initial velocity of the projectile, the uneven arch height around the perforation on the back of the target plate becomes more and more serious. In addition, the maximum heights of bulge on the back of the target plate in the four tests are 1.78 mm, 2.62 mm, and 3.62 mm, respectively. Thus, the maximum arch height on the back of the target plate increases with the initial velocity of the projectile.

Figure 10 shows the cracking condition of the surface of the plug after the above four tests. It is found that cracks appear on the plug surface, which is similar to the phenomenon in the literature [18]. It was caused by the fact that the tensile strain produced during the forming process of the plug exceeds the fracture strain of the material.

**(a)**

**(b)**

**(c)**

**(d)**

Figure 11 shows the whole process of the target under projectile impact at *V*_{i} = 327.55 *m/s* by a high-speed camera. The forming process of plug is marked in the red box. The failure model of target under projectile impact is shear plugging; a complete plug was formed. In addition, Figure 11 shows that the flight attitude level of the projectile before impact on the target plate and no obvious deformation occurred during the perforation.

**(a)**

**(b)**

**(c)**

**(d)**

**(e)**

**(f)**

#### 4. Computational Model and Numerical Simulation Result

##### 4.1. Material Model

A half of a full-scale 3D finite element model was developed by Abaqus/Explicit, and the modeling process is referred to literature [18], as shown in Figure 12. In the process of numerical simulation, the impact of the projectile will result in large deformation of some elements. When the equivalent plastic strain of an element is greater than or equal to 3, it is assumed that the element loses its bearing capacity due to large deformation and the unit is automatically deleted. Because the time of impact is rather short, the impact progress is considered to be adiabatic, and the temperature rise can be calculated as [16]where *ρ* is the density, *C*_{p} is the specific heat, and *χ* is the fraction of plastic work converted to heat. It is usually accepted that *χ* = 0.9.

M-Paul, modified J-C, and J-C constitutive models were selected for the material strength model, and the modified J-C fracture criteria were adopted for the failure model. Because M-Paul constitutive model, modified J-C constitutive law, and fracture criterion are not provided in Abaqus material base, corresponding VUMAT subroutines should be developed to complete the input of material model by calling it during the numerical simulation. And, the M-Paul constitutive model is expressed as Equation (9), and the relevant material parameters are shown in Table 2.

Lin li et al. [12] revised thermal-softening fraction of J-C constitutive law and fracture criterion, as shown in Equations (16) and (17). The relevant material parameters are shown in Table 5:

The J-C constitutive model is expressed as Equation (18). Except *m* = 0.757, the other parameters are shown in Table 5:

Refer to the thermal-softening item (1 − *FT ^{∗}*

^{m}) of MJ-C constitutive model in literature [12], the thermal softening item of J-C fracture criterion was modified in this paper. The modified fracture criterion is expressed as

The parameters of thermal-softening items *D*_{5} and *D*_{6} in the above fracture criteria were recalibrated, as shown in Figure 13. A curve fit gave *D*_{5} = 42.996 and *D*_{6} = 4.216.

##### 4.2. Numerical Simulation Results

In this paper, the numerical simulation of Q235B steel target plate under blunt projectile impact was performed. The residual kinetic energy of the projectile after perforation is output and converted into the residual velocity of the projectile.

* Case 1. *The numerical simulation of the ballistic test was carried out by the M-Paul constitutive model and modified J-C fracture criterion.

* Case 2. *The numerical simulation of the ballistic test was carried out by the modified J-C constitutive model and fracture criterion.

* Case 3. *The numerical simulation of the ballistic test was carried out by the J-C constitutive model and modified J-C fracture criterion.

In cases 1 and 2, a total of 7 numerical simulations were completed, respectively. The initial-residual velocity and the states of projectiles and targets after impact are shown in Table 6, where C1 represents the numerical simulation under Case 1, and the latter is the same. The ballistic limit ranges from 283.8 m/s to 290 m/s in Case 1, and the ballistic limit ranges from 290 m/s to 303.2 m/s in Case 2.

In Case 3, a total of 8 numerical simulations were completed. The initial-residual velocity, deformation, and failure of projectiles and targets are shown in Table 6; the ballistic limit ranges from 323.3 m/s to 340 m/s. Based on formula R-I, the ballistic limit is fitted by using the method of least squares, given

*V*

_{bl}= 327.2 m/s, which is 16.5% more than the test value of 280.88 m/s. Due to the larger error, the numerical simulation result of Case 3 is ignored (Table 7).

A comparison of ballistic limit between ballistic test and numerical simulation is shown in Figure 8. The initial-residual velocity curve in Case 1 is in close agreement with the test result compared with Case 2. The ballistic limits of cases 1 and 2 that can be fitted by formula R-I are 285.9 m/s and 299.8 m/s, respectively; both are found to be 1.8% and 6.7% more than the test ballistic limit, which is 280.88 m/s. The R-I model constants of ballistic test numerical simulation for Q235B steel target are shown in Table 8.

By observing the numerical simulation process of target plate subjected to projectile impact at different velocities under different working conditions, the test phenomenon is shown in Table 9.

As Table 9 shows, when the initial velocity of projectile is lower than the ballistic limit (the target under projectile impact with low velocity), the stress distribution and target failure model of Case 1 are the same as those of Case 2. It indicates that M-Paul and modified J-C constitutive models have little difference in predicting mechanical behavior of materials at a low strain rate. When the initial velocity of projectile is higher than the ballistic limit, the stress distribution and target failure model of Case 1 are similar to those of Case 2. But, the deformation of the target after impact is different, the height of bulge on both sides of the rear hole is uneven in Case 1, and the occurrence is not obvious in Case 2, as shown in Figure 14. The phenomenon in Case 1 is consistent with the above test results.

In addition, it can be seen from Figure 14 that with the improvement of the initial velocity of the projectile, the arch height on the back of the target plate in Case 1 has increased, which is consistent with the test result, while no similar phenomenon has been observed in Case 2.

The above phenomenon indicates that with the improvement of strain rates, M-Paul and modified J-C constitutive models began to show difference in the prediction of material mechanical behavior. The numerical simulation result of Case 1 is more consistent with the test result; this illustrates that the mechanical behavior of Q235B steel can be better predicted by the M-Paul constitutive model.

**(a)**

**(b)**

**(c)**

**(d)**

#### 5. Conclusions

In this paper, the dynamic constitutive model of Q235B steel was studied systematically through theoretical derivations, numerical optimization, tests, and numerical simulation. The M-Paul model for characterizing Q235B steel mechanical behavior under larger strain and high strain rate and temperature is proposed. The tests and numerical simulations of Q235B steel target plate under blunt projectile impact were performed. The mechanical properties of material under impact was investigated, and the accuracy of M-Paul model was verified. The results are as follows:(1)When the parameters of the Paul constitutive model are calibrated by combing theory with numerical optimization, it was found that thermal-softening item could not reasonably describe thermal-softening behavior of Q235B steel, so the original model was modified to obtain the M-Paul constitutive model. By comparing the prediction of flow curves under different strain rates and temperatures between M-Paul and Johnson–Cook constitutive models, it is found that M-Paul constitutive models can more accurately represent the dynamic mechanical behavior of Q235B steel.(2)The ballistic test of Q235B steel target plate under blunt projectile impact was performed; three failure models, i.e., back bulge, shear plugging failure, and adiabatic shear failure, were concluded. It was found that with the improvement of the initial velocity of projectile, the phenomenon that the bulge of concave crater is jagged becomes more serious and, also, the maximum bulge height is increased. In addition, the tensile strain produced in the forming process of the plug exceeds the fracture strain of the material, and thus, cracks appear on the surface of the plug.(3)The numerical simulation of Q235B steel target plate under blunt projectile impact was performed. By comparing the ballistic limit and the numerical simulation results, it was found that the ballistic limit error obtained by using the J-C constitutive model was very large, and then, compared with the modified J-C constitutive model, M-Paul can more accurately describe the dynamic mechanical property of Q235B steel under complex conditions such as large strain and high strain rates and temperature.

#### Data Availability

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

#### Conflicts of Interest

The authors declare that they have no conflicts of interest.

#### Acknowledgments

The authors wish to acknowledge the support from the National Science Foundation of Heilongjiang Province of China (no. LH2019E060) and the National Nature Science Foundation of China (grant nos. 51578515 and 11502120).