Abstract

For a pile foundation design, the value of proportional coefficient m to define the soil horizontal resist force is a significant parameter. However, different geological conditions and experimental environment have led to different m values. In this paper, an in situ test is firstly carried out on the horizontal bearing capacity of large-size precast square-piles. The piles deformation is then derived by using the optimization method from the measured data. Secondly a back analysis model is established to calculate the m value by using the simplex method, which reveals the evolution rule of the value of proportional coefficient m. Results show that the horizontal bearing characteristics of precast piles depend on the interaction force of piles and soils. The action mechanism of the soils around the piles is gradually developed with the increase in the concrete content. The horizontal critical load and the Eigenvalue of horizontal bearing capacity increased by 16.7% and 20%, respectively. It is also seen that the higher the content of the cement-soil around the piles and the longer the pile length, the bigger the m value obtained. The variation of the proportional coefficient m with the horizontal displacement of pile top is defined by three stages: rapid decaying stage, slow decaying stage, and balanced stage, respectively. The inverse analysis method on the proposed m value can accurately reflect the actual working state of piles and soils. In the depth of 3~18m in the west of Ji'nan, the range of m value is recommended as 4~6.58 MN·. When Δ takes 12mm, the values of m are consistent with the result from the back analysis. In summary, the obtained m value can be effectively used to guide the design of enclosure structure in the super deep foundation pit in the Yellow River alluvial stratum.

1. Introduction

The thick alluvial clay stratum in Ji’nan belongs to the quaternary system with its major content silt, silty soil, and clay. The upper stratum is loose and underconsolidated with poor physical mechanics properties, low bearing capacity, creep, and thixotropic property. The vein network channel of groundwater extends in all directions and a field deep excavation shows distinct spatial effect. The construction of metro stations has exacerbated the tension of urban underground space, making the characteristics of “deep, big, tight, and near” of foundation pit more obvious. The safety of foundation pit excavation is facing great challenges. Due to the constant exploitation of super high-rise buildings and urban underground space, super large foundation pit engineering is constantly emerging. Design of the foundation pit enclosure structure with the combination of the support structure and the main structure emerges at a historic moment. However, few systematic researches on large-size precast piles as principal enclosure structure have been conducted in deep and large foundation pit. The interaction mechanism between prefabricated components and cast-in-situ components is complicated and the design backup of critical supporting parameters is insufficient. Therefore, it is of a vital significance to accurately determine the design parameters of precast piles in deep foundation pit for the construction safety of geotechnical engineering in the alluvial stratum of the Yellow River.

A pile foundation is supported by a stable soil layer or embedded in bedrock through liquefiable soil layers. In a liquefied and seismic subsidence of the shallow soil layer caused by the earthquake, a pile foundation still has sufficient bearing capacity to ensure the stability of the foundation without excessive deformation [13]. However, determination of the horizontal and the vertical bearing property of the pile foundation is a complicated problem under pile-soil interaction. The bearing characteristic covers two aspects: the ultimate horizontal bearing capacity of pile controlled by the strength of soil around piles or the allowable horizontal bearing capacity of pile determined by the allowable deformation of piles [4]. The calculation method of the bearing capacity of a single pile mainly includes p-y curve method [57], analytical method [810], finite element simulation [1115], and artificial intelligence method [16]. Moreover, experimental method is also an effective way to determine the horizontal bearing capacity of a pile foundation. Huang et al. analyzed the horizontal bearing characteristics of bored piles by horizontal static loading tests [17]. Haque et al. used horizontal static loading test to measure the pile strain and analyzed the distribution of pile load transfer under static loading [18]. Zhao et al. used the indoor model test to study the bearing characteristics of a single pile under the combination of vertical and horizontal loadings [19]. The researches mentioned above show that the scale effect of horizontal bearing capacity of pile foundation in laboratory tests cannot be neglected. Also it is difficult for a model test to insure the reliability. In addition the accuracy of numerical simulation depends on the accuracy of test parameters. Comparatively, in situ test of piles is more effective for testing the bearing capacity of pile foundation.

Soils around piles stabilize pile body and resist the external lateral force acting on piles. Taking the soils around piles as an elastic medium, and assuming piles as an elastic beam which supports the soils around piles, the soil resistance force is determined by the linear elastic reaction method to derive the elastic resistance property of the soils. The deformation differential equation of the elastic foundation beam is as follows [24]:where EI is the flexural rigidity of the enclosure structure, E is the elastic modulus of piles, I is the section moment of inertia of the pile, y is the lateral displacement of the enclosure structure, z is the excavating depth, ea(z) is the active earth pressure at the excavating depth z, hn is the excavating depth at the step n, and m is proportional coefficient of horizontal resistance of foundation soils.

The value of m is the proportional coefficient that describes the horizontal resistance coefficient with the depth. When the pile section type, pile stiffness, pile length, and depth of pile bottom are determined, the m value is entirely dependent on the characteristics of the soils around piles. Due to the unreasonable design of row piles accidents frequently occurred in large and deep foundation pits. One of the key factors is the inappropriate selection of the ratio coefficient of the horizontal resistance of foundation soil (m value). Researches indicate that the value of m is a multifactor and regional parameter which relates to soil properties, the pile bearing property under static load, dynamic load, or multicycle grade load, and displacement, etc. [25, 26]. For a long time, researchers have proposed several methods to determine the value of m. Fan et al. studied the proportional coefficient m in muddy soft rock foundation and provided a reference value through indoor horizontal static load tests and numerical simulations [27]. Huang et al. considered calculating the value of m of the composite pile based on the design value of horizontal bearing capacity and the corresponding deformation [28]. Xu et al. determined the proportional coefficient in vertical elastic subgrade beam by using back analysis software Ucode and finite elements software Abaqus combined with the measured deformation of retaining walls [29].

The above researches provide important reference for design of row piles in deep foundation pit. However, the pile foundation bearing capacity and the design parameters m value present obvious regional features —the soil characteristics are complex and different soils have different parameters. As the water level decreased in the Yellow River alluvium for recent years, the shallow soil has better performance in self-stability and generates microconsolidation. The veins network of groundwater extends in all directions, which are different from the properties of ordinary clays and soft soils. The excess pore water is commonly observed. Given the above condition, a new pile-wall superimposition enclosure structure is adopted to support the deep foundation pit in this area (see Figure 1). The two factors, complexity of regional stratum conditions and the new support method, make some results not clear, such as horizontal bearing characteristics of precast piles, and the selection rule of the m value. This is an urgent problem for the construction of extra-large deep foundation pit in alluvial stratum of Yellow River. In this study firstly a variable unique research of large-size precast square-piles by in situ test is conducted to analyze the bearing characteristics of the new pile-wall superimposition structure. Then a calculation model is established to determine the m value based on the engineering survey data and the optimization method. The obtained results can provide theoretical and technical guidance for the enclosure structure design of the deep foundation pit in the alluvial stratum and similar stratum.

2. Field Test

2.1. Parameters of Stratum and Soil

The testing site is located at the Yan Mazhuang Station, Metro line R1, Ji’nan. The design depth of foundation pit is 16~18m and the maximum length of precast piles is 26.6m with 700mm in diameter. The main enclosure structure adopts cement-soil cast-in-place piles with 1100mm in diameter where precast piles with 700×700 mm rectangle section are inserted. The standard interval between piles is 1.5m. The waterproof curtain adopts triple-pipe jet grouting piles with 800mm in diameter and 550 mm in spacing. Enclosure piles are inserted into the cement-soil cast-in-place pile and the length of the cement-soil pile is as long as the enclosure piles. The alluvial stratum is the alluvial-diluvial geomorphic units of the Yellow River and Xiaoqing River. The covering layer is the quaternary alluvial-diluvial strata which is mainly composed of silty clay, silty clay and clay, or sandy soil and pebbles in some regions. The base bearing soil layer is compressible and mainly clay layer and the local is silt. Generally the foundation soils are stable and homogeneous. Hydrogeologic conditions and the distribution of precast piles are shown in Figure 2 and the relevant parameters of each layer are shown in Table 1.

2.2. Test Equipment

Layout of the horizontal static loading test is shown in Figure 3. In the test the horizontal force is applied by the lifting jack. A spherical hinge seat is placed at the contact between the jack and the test pile to ensure that the jack force can pass horizontally through the axis of the pile body. The horizontal displacement of the pile is measured by large range dial gauges. Two dial gauges are installed on the plane surface and 50cm above the plane. Lower gauge is used to measure the horizontal displacement of the pile on the ground, while the upper one is used to measure the horizontal displacement on the top of the pile. The rotation angle of the pile body above the ground can be obtained according to the ratio of the displacement difference between the two tables and the distance of the two tables. The base pile for fixing the dial gauges is set in the opposite direction of the test pile, while the net distance to the test pile is less than one time of the pile diameter. With the excavation of the foundation pit, the stress balance is destroyed and the superposition enclosure structure is mainly deformed towards the inside of the foundation pit. Therefore, horizontal loads are applied perpendicularly to the outside of the pit in the test. The horizontal counterforce is developed from the consolidated soils. The horizontal thrust is applied by the hydraulic jack with a spherical joint, while the horizontal action line coincides with the top elevation (± 0.00m) of the pile foundation.

2.3. Test Method

According to the requirements of Technical Code for Testing of Building Foundation Piles (JGJ106-2014) [30], one-way multicycle graded loading method is adopted in the test. The loading capacity of each grade is 1/10 of the design bearing capacity of pile foundation. After each grade of load is applied, the load remains unchanged for 4 minutes before the horizontal displacement is measured. After data recording, unloading is conducted to zero and the load is maintained for 2 minutes, when the residual deformation is recorded. A load cycle is completed according to the above steps. Deformation detection of the first grade loading is completed after five cycles.

2.4. Test Pile Type and Basic Parameters

The pile tests are based on the principle of unique variable. Three piles are selected to conduct the experiment, in which two are equipped into one group. Three groups are adopted in total in the experiment. Specific parameters are listed in Table 2, where the PCSP is the abbreviation of precast concrete square-pile

3. Test Results

3.1. Relation between Displacement of Piles and Horizontal Load

Figure 4 shows the relationship between the pile top displacement Y0 and the horizontal load H. The displacement trajectory of PCSP No.2 and PCSP No.3 is the same before loading to 500kN. Due to the breakage of PCSP No.2, the loading process is stopped and the maximum horizontal displacement is 11.33mm.

Compared with the displacement curves of PCSP No.1 and PCSP No.3, the displacement of pile top increases gradually during loading process and the curve shows a parabolic trend. Under the same stress level, the displacement of PCSP No.1 pile with lower cement content is more than that of PCSP No.3 from the beginning to the end and the displacement difference increases gradually. Finally, the differential value of the maximum displacement is 6.61mm. When the piles are subjected to the lateral load, the cemented soils around piles play important roles for stabilizing piles and resisting the external force transfer from piles.

In order to more intuitively reveal the effect of horizontal loads on the deformation of retaining structures, the ratio of horizontal load change per displacement unit (ΔY H) over the horizontal increment is summarized, as shown in Figure 5. The ordinate represents the characterization of displacement grade under loading, that is, the displacement variation amplitude under a unit force. The abscissa represents the applied horizontal load in the test.

From the pile displacement grade curve of PCSP No.3 (Figure 5), it is seen that in the range of 100kN to 400kN, the variation of the displacement grade is not large. The maximum displacement grade observed after the applied horizontal load is more than 400kN. Similar to that of PCSP No.3, the displacement grade of the PCSP No.2 increases suddenly and tends to increase with a higher rate when loading is beyond 400kN. On the contrary, PCSP No.1 without the cement-soil mixture has a significant increase in the grade curve. The grade curve has the same trend in the loading interval between 160kN and 480kN. After loading to 480kN, the slope of the grade curve suddenly increases, which indicates that the use of cement-soil enhances the bearing capacity of the pile in a certain range.

3.2. Analysis of the Relationship between the Displacement (Y0) and lgt

Figures 6(a), 6(b), and 6(c) show the Y0-lgt relationship curves of PCSP Nos.1, 2, and 3 in static loading tests, respectively. It is seen that, when loading is applied to 400kN, 400kN, and 480kN, respectively, the displacement is obviously higher than the displacement difference mentioned in the previous section. According to Technical Code for Testing of Building Foundation Piles [30], the horizontal critical loads of PCSP Nos.1, 2, and 3 can be judged as 400kN, 400kN, and 480kN, respectively, and the Eigenvalue of horizontal bearing capacity Rh is 300kN, 300kN, and 360kN, accordingly. It is seen that the maintenance effect of soils around the piles is improved by 8.9%, while the horizontal critical load and the Eigenvalue of horizontal bearing capacity increase by 16.7% and 20%, respectively.

4. Analysis of the Value of in Thick Alluvial Foundation Soils

Support design of a deep foundation pit in thick alluvium area is always a very complex common problem. The determination of the ratio coefficient m of the horizontal resistance for foundation soils also has great controversy in various areas due to different geological conditions. The initial value of m in each region in China was recommended by Xilin et al. (former Soviet Union), which was then adopted by the Technical Code for Design of Railway, Highway and Urban Road Bridges and Culverts of the former Soviet Union [31]. In 1974, the China research center corrected the recommended value of m by single pile horizontal static loading tests. The Code for Design of Foundation and Foundation of Highway Bridge Culvert (JTJ024-85) [22] adopted the revised value in 1985, and the railway ministry took the same value as the benchmark. Now the various codes and the research contents of scholars are used in this study for comparison and collected in Table 3.

Although the current specifications and guidelines present recommended values for m values, they are only obtained by test statistics. The suggested m values for each type of soils are in a large range, and it is still difficult to accurately determine the values.

4.1. Calculation of m Value Based on the Horizontal Static Loading Tests

The target objective of vertical elastic subgrade beam method is used to solve plane strain problems. In this method, the enclosure structure is assumed as a vertical elastic foundation beam. The effect to the enclosure structure induced by the soil mass below the excavating surface is stimulated by a series of Winkler springs, while the supporting effect to the enclosure structure induced by the lateral bracing (steel shotcrete or concrete support) is stimulated by the elastic mountings. In order to achieve the balance of calculation, the effect of the soils behind the wall acting on the enclosure structure is represented by the Rankine earth pressure.

Due to the soil stratification, the values of m are different in different layers. Therefore, the beam element in the finite element method is adopted. That is, the elastic foundation beam will be divided into some units based on the soil layers along the vertical direction, and the differential equations above of each unit are listed to solve them subsequently.

The method of flat vertical elastic foundation beam for the foundation pit enclosure structure is substantially evolved from the computational method of laterally loaded piles. Therefore, the determination of the value of m is actually based on the results from the horizontal load tests.

According to Technical Code for Testing of Building Foundation Piles (JGJ106-2014) [30], the computational equation of m value and the horizontal deformation coefficient α arewhere H is the horizontal thrust of single pile, x is the horizontal displacement caused by pile under loading, b0 is the calculative width of pile, vx is the displacement coefficient at the pile top, according to the Technical Code for Building Pile Foundations [30], α is calculated by comparing the values of αh and finding the corresponding values in Table 4, and h is the injection depth of pile.

For circular pile, when pile diameter D is less than or equal to 1m, b0=0.9(1.5D+0.5). When pile diameter is more than 1m, b0=0.9(D+1). For precast concrete square-pile, when the sectional dimension B is less than or equal to 1m, b0=1.5B+0.5. When B is more than 1m, b0=B+1.

Determination of flexural rigidity (EI) of piles is explained as follows. In the loading process on the precast concrete square-piles, the cemented soils and the precast core pile work together. Comparing the effects of flexural rigidity, it is seen that the flexural rigidity of core pile is 6.9×105 kN·m2, while the flexural rigidity of cement-soil is 0.32×105 kN·m2, so the /=4.6%. Obviously, the flexural rigidity of cement-soil is relatively tiny, so that it is replaced by the flexural rigidity in this paper.

According to (2), the variation trend of the value of m in different loading grades is obtained, as shown in Figure 7. The relationship between the value of m and horizontal force is resembled in an inverse proportion relation. When loading to 200kN, the m value of PCSP Nos. 1, 2, and 3 is 277.3MN/m4, 157.4 MN/m4, and 100 MN/m4, respectively. With the increase of the horizontal loading, the value of m in each case gradually decreases and eventually tends to be the same. In the case of large displacement of piles, the coefficient of horizontal resistance of the pile from the soils around piles is gradually reduced, and it gradually tends to be a stable value.

It is also found that the values of m vary with the formation and environment. It seems that there is a special relationship between the value of m and the load magnitude in the test, such as the displacement at the pile top. From the test data, the relation between the value of m and the horizontal displacement at the pile top is described as shown in Figure 8.

According to the m value curves of PCSP Nos. 1, 2, and 3, the relation between the coefficient m and the horizontal displacement at the pile top can be divided into three stages: rapid decaying stage, slow decaying stage, and balanced stage, respectively.

(a) When the horizontal displacement of the pile head Y0 ranges from 0 to 5mm, the value of m is on the first rapid decaying stage.

(b) When the horizontal displacement of the pile head Y0 ranges from 5 to 15mm, the value of m is on the second slow decaying stage.

(c) When the horizontal displacement of the pile head Y0 is more than 15mm, the value of m is on the third slow decaying stage, balanced stage.

It is found that the values of m for three groups of piles are eventually maintained at the section between 10 and 30 MN· along with the increase in the horizontal displacement at the pile top. On the one hand, the value of m of the precast pile with different variables varies greatly in the first stage. The relative values of the second and third phases decreased rapidly, tending to be in a uniform stable interval. On the other hand, although the piles are under different circumstances, they vary following the same trend. The three change stages of m value and the displacement at the pile top referred to in this paper all apply to the different circumstances for both precast piles and bored piles.

4.2. Displacement Back Analysis of the Value of m

The principle of the displacement back analysis calculation method adopts the displacement obtained from the actual measurement of the project, and the original stress or mechanical parameters are back-derived according to the selected calculation model. This method was put forward in 1971 by K. T. Kavanagh and R. W. Clough [32]. Due to the maneuverability and reliability of the displacement measurement, in recent years the displacement back analysis method is applied widely in the field of geotechnical engineering [3337]. The displacement inversion of m value is based on the compiling of the acquisition information of pile displacement. The ‘m’ method and the boundary conditions of the mechanical model of the enclosure structure are optimally selected to serve the objective function. The difference between the theoretical value and the measured value is minimized by the constant correction on a given initial value. Finally, optimization technique is used to derive the resistance coefficient of each soil layer. In order to make the calculated value mostly close to all the measured values, the square sum of deviations is required to be the smallest; thus the objective function can be described aswhere M is the construction step, N is the test number, and and are the calculated value of horizontal displacement and the measured displacement of piles based on the test measuring point i in the construction step j, respectively. The measured values are obtained by using the inclinometer along the precast piles.

In (4), the soil layer parameters m1, m2, m3 mn minimized by the objective function are the parameters to be sought, and they can be solved by the simplex optimization method [3840]. Specifically, a group of initial values used to develop the initial simplex according to certain rules is provided first. And then the optimal value is obtained by searching the step length. The accuracy of the convergence is controlled by the step length.

As the proportional coefficient is quite different between the upper and lower limits of the soil layers, the m value of various types of soil mass has great limitations. The following constraints are introduced:where and are the upper and lower bound of values for the foundation soil in the layer i. Combining (4) and (5), a set of nonlinear programming problems is obtained, e.g., the problem of function optimizing under constraint conditions. This method can not only fully exploit the characteristics of the back analysis algorithm, but also make the obtained m values meet the requirements. The flowchart of back analysis is shown in Figure 9.

4.3. A Calculating Example

In this paper, the deep foundation pit in the Yanmazhuang west station of Ji’nan rail transit line R1 is taken as an example. The deep excavation adopts the cut-and-cover method with 16.5m in excavation depth. The enclosure structure is composed of 700×700mm precast square-piles and two steel braces and a concrete support. The 400×400mm precast core columns are used in the middle of the enclosure structure. After the excavation of the foundation pit, the encased concrete forms the permanent structural columns. The excavation sequence of foundation pit is shown in Table 5.

The back analysis is based on the measured displacement of the inclinometer pipe of the support structure under different working conditions. The initial value of the proportionality coefficient m is substituted into the calculation process of back analysis as shown in Figure 9 for iterative computations until a reasonable m value of the foundation soil is obtained. Figure 10 shows the influence curve of the iterations for the optimal value of the inversion. It is seen that the convergence of the inversion process meets the requirements.

Figure 11 shows the comparison between the deformation of the enclosure structure calculated from the initial m value and inverse m value, and the measured deformation when excavating to the bottom of the foundation pit. Results show that the measured deformation is much larger than that of the deformation calculated from the initial m value. The maximum horizontal displacement difference is more than 4.2mm. The measured data is 66% higher than the calculated data. It is believed the m value of the foundation soil is obviously larger. After the iteration of the back analysis, the range of m value for each soil layer in this area is in the range of 4~8 MN·m−4. It has been found that the deformation obtained by the inversion of m values is consistent with the measured one.

5. Discussion on the Empirical Equation of Value

The thick alluvium stratum of the Yellow River in Ji’nan area is similar to the soft soil stratum with obvious regional characteristics. When m value is applied in the calculation of elastic foundation beam, the recommended values of building codes or other local codes should not be copied directly. At present, experiences have been achieved for the calculation of m value in soft soil and clay areas in southern China. In the absence of any single horizontal static load test and lack of experience, the following empirical calculation method is used according to the Technical Specification for Retaining and Protection of Building Foundation Excavations [41]:where φ is the standard value of internal friction angle (°), c is the standard value of cohesion (kPa), and Δ is the displacement at the bottom of foundation pit, which can be evaluated by the regional experience. 10mm is suggested in the absence of engineering experience.

By taking the typical soil layer, the different frictional angles and cohesiveness of natural quick shear test and consolidation quick shear test are obtained by consulting the geological survey data. The m value is obtained by substituting it into the empirical equation. Results are shown in Table 6.

It is obvious that the value of m obtained by applying the cohesive and the internal friction angle under the natural fast shear test is close to the value obtained from the back analysis. However, the value of m obtained by consolidated quick shear test is significantly larger than that of the back analysis method. The empirical equation deserves a further optimization by the parameters of natural fast shear test.

Figure 12 is the experiential correction curve of (6). As the value of Δ in the equation is difficult to determine, the calculated value of m is quite different from that of the regional experience value. Yang et al. pointed out that the m value of the rock stratum in Guangzhou is obviously lower by using empirical equation [42], while in the local standard foundation pit engineering technical specification of Hubei province, the above empirical equation had multiplied an empirical coefficient. Specifically, the empirical coefficients take 1.0 for general cohesive and sandy soil, 1.8 ~ 2.0 for old cohesive soil and the gravel stratum over medium density, and 0.6 to 0.8 for silt and mucky soil. In view of the blank in the correction of the empirical equation of m value in the alluvial stratum of Yellow River, the experience value of Δ is verified and the obtained results are compared with the value of m obtained from back analysis. It is seen that the value obtained by (6) is larger, which is contrary to that in Guangzhou. When Δ takes 12mm, the values of m are consistent with the result of the back analysis. Therefore, it is more reasonable for the horizontal displacement Δ to take 12mm in this problem.

6. Conclusions

This paper raises the research issue on the horizontal bearing capacity of the precast pile and the key supporting parameters of the deep foundation pit in the Yellow River alluvium. Variable unique analysis for the large-size precast square-piles is conducted by in situ test, and the bearing characteristics of the new pile-wall superimposition structure are explored. A calculation model is established for resistance coefficient m of foundation soil based on the engineering survey data and the optimization method. The conclusions are drawn as follows.

(a) As a new type of precast piles, the increase in cemented soils can strengthen the bearing capacity of a single pile to some extent. Under the joint action of pile and cemented soils, the envelop effect of the soils around piles increase by 8.9%, while the horizontal critical load and Eigenvalue of horizontal bearing capacity increase by 16.7% and 20%, respectively.

(b) According to the test results, the ratio coefficient of the horizontal resist force in foundation soils is not a constant. The higher the intensity of the cemented soils around the piles and the insertion ration, the bigger the value of m. With the increase in the horizontal displacement at the pile top, the value of m undergoes three stages: rapid decaying stage, slow decaying stage, and balanced stage. The change mode of each stage is different.

(c) A back analysis model is established for the m value calculated by using the simplex method, which reflects the m value in real work conditions. In the west area of Jinan city, the alluvial stratum of the Yellow River in the depth of 3~18m is regarded as a normally consolidated viscous stratum, and the range of m value is recommended to be in the range of 4 ~ 6.58 MN·.

(d) Based on the measured data, the displacement of the support structure is predicted by the calculation parameter, which proves that the value of m is more accurate with the internal friction angle and cohesive force under the natural quick shear. When Δ takes 12mm, the values of m are consistent with that from the back analysis.

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 paper is supported by the National Natural Science Foundation of China [Grant nos. 51774196 and 41472280], China Postdoctoral Science Foundation [Grant no. 2016M592221], and SDUST Young Teachers Teaching Talent Training Plan [Grant no. BJRC20160501].