Abstract

Loess has a unique structure and water sensitivity, and the immersion of loess leads to many tunnel lining problems in shallowly buried tunnels. Based on a tunnel in Gansu Province in China, two failure glide planes of a shallowly buried loess tunnel and their immersion modes are summarized. Finite element calculation of the structural Duncan-Chang constitutive model is realized via the secondary development of finite element software, by which the loads on the secondary lining are calculated and verified in comparison with measured results. The load characteristics of the secondary lining are studied. The load evolution is closely related to the immersion position and scope of the load. After the loess near the failure glide plane of the arch foot is flooded, the load on the arch foot sharply increases. As the immersion expands, the maximum load moves from the arch end up to the hance. After the loess near the failure glide plane of the hance is flooded, the load on the hance decreases slightly. The stability of the overlying loess decreases gradually, which causes the loads on the vault and arch shoulder to rapidly increase. Additionally, the load distribution characteristics on the secondary lining are summarized.

1. Introduction

With the rapid development and construction of traffic infrastructure in Northwest China, tunnels through the loess area are continually increasing in number and expanding in scale. The special physical and mechanical properties of loess, such as its porosity, macroporosity, soluble salt content, and collapsibility, lead to significant changes in the strength and deformation of loess after water immersion, which greatly impact the stability and safety of existing loess tunnel linings, and a series of problems with the lining structure gradually occur. In loess areas, tunnel defects are closely related to the loess structure and its water sensitivity, and both of these factors should be considered when simulating the deterioration of flooded loess in loess tunnels.

Loess has significant structural characteristics, which is an important reason for the loss of strength and stability of loess under the influence of water. The Duncan-Zhang constitutive relationship is close to the stress-strain relationship of loess [13]. The structural Duncan-Chang constitutive model, which is considering the loess structural influence, can further improve the applicability of Duncan-Chang constitutive model. Studies on loess structure and soil structural constitutive models have been conducted. Xie et al. [4] proposed a structural parameter that comprehensively reflects soil structure. Li et al. [5, 6] established a nonlinear damage constitutive equation for unsaturated loess considering structural effects. Luo et al. [7] and Shao et al. [8] proposed a structural constitutive relation and suggested a humidification constitutive model considering the structure of unsaturated loess. Zhao et al. [9] established an elastic-plastic constitutive model of strain softening and strain hardening of structural loess by a stress-controlled triaxial test. Feng [10] introduced the soil structural parameter into the stress-strain relationship of undisturbed loess and established a structural Duncan-Chang constitutive model of unsaturated undisturbed loess.

Studies on loess structural constitutive models have been widely conducted, which has laid a solid foundation for the application of loess structural constitutive models in loess tunnel engineering. The structural Duncan-Chang constitutive model was shown to accurately reflect the stress-strain relationship of loess based on a comparison with experimental data.

Scholars have performed many studies on the influence of the deterioration of loess on tunnel structure through theoretical analyses and model tests. Reznik [11] and Lai et al. [12, 13] analyzed the stress and deformation characteristics of tunnel structures caused by different water immersion scopes and discussed the mechanism of lining cracking. Jeong et al. [14] analyzed the changes in the shear strength and compressive modulus of loess with water content variations and summarized the influence of law of loess water content on tunnel deformation. Shao et al. [15] divided collapsible loess into environmental grades, which contributes to tunnel design. The causes of loess tunnel lining degradation can include farmland irrigation water, indistinct division of surrounding rock levels, stress evolution of surrounding rock, and lining construction quality [1619]. Liu et al. [20] studied the causes, sequence, regularity, and characteristics of the lining cracking process in loess highway tunnels under the flooding of the surrounding rock. Inokuma and Inano [21] asserted that large loads and displacements would directly affect the stress and deformation of the tunnel structure and its operation safety and that the incidence of tunnel diseases was not directly related to the duration of operation.

In previous studies, the specific immersion water location and immersion scope were not reported, and the constitutive relation used in numerical simulation does not reflect the characteristics of loess. Additionally, the loads on the secondary lining under different immersion cases have not been reported.

Based on a project involving a loess tunnel in Gansu Province in China, this article summarizes four types of failure glide planes and their water immersion modes. Then, the formula for loess structure parameters is obtained using neural network fitting, and the finite element calculation of the structural Duncan-Chang constitutive model is realized through the secondary development of finite element software. The reliability of secondary development is verified by numerical simulation of triaxial tests. The loading characteristics of the secondary lining of shallowly buried loess tunnels under four immersion models are discussed.

2. Project Background

2.1. Project Overview

This article studies a loess highway tunnel in Gansu Province, which is located on the Loess Plateau and is covered by thick loess. The strata lithology around the tunnel can be divided into three types from top to bottom: Quaternary Upper Pleistocene loess Q3, Quaternary Lower Pleistocene loess Q2 and Upper Tertiary Pliocene mudstone N2.

The natural bulk density of Q2 loess is 20.26 kN/m3, the natural water content is 21.72%, the liquid limit is 28.34%, the plastic limit is 20.12%, the cohesive force is 28.67 kPa, the internal friction angle is 28°28′, and the natural porosity ratio is 0.63.

The natural bulk density of Q3 loess is 15.96 kN/m3, the natural water content is 15.66%, the liquid limit is 28.54%, the plastic limit is 19.99%, the cohesive force is 45.29 kPa, the internal friction angle is 28°03′, and the natural porosity ratio is 0.98.

The monitoring section of the secondary lining is located in the Quaternary Upper Pleistocene loess Q3, with a burial depth of 30 m. The mountain surface of this section is an agricultural irrigation area, as shown in Figure 1. The secondary lining is made of C25 concrete with a thickness of 40 cm, whereas the primary support is shotcrete with a thickness of 25 cm. The inverted arch is made of C20 concrete with a thickness of 25 cm, as shown in Figure 2.

Two surface cracks are found at the surface of the shallow burial area of the tunnel site, as shown in Figure 3. The cracks extend basically parallel to the longitudinal direction of the tunnel and are distributed on the left and right sides approximately 23 m from the tunnel axis. The surface is land irrigated by flooding in spring and winter every year, and the irrigation water aggravates surface cracks.

During the excavation of shallowly buried loess tunnels, the overlying loess moves downward. It was observed that the loess at the arch foot cracked first, the columnar joints propagated due to the settlement of the loess, and the loess cracks gradually propagated upward toward the surface. Finally, the failure glide plane in the surrounding rock of the shallowly buried tunnel appeared, as shown in Figure 4. The failure glide planes were also found to extend from the arch foot to the surface at an inclination of approximately 70°, which is consistent with a previous report [22]. Additionally, by means of physical detection and pit exploration, the failure glide planes in the surrounding rock were found to extend from the arch waist to the surface at an inclination of 64.3° [23].

During rainfall and agricultural irrigation, the failure glide planes in the surrounding rock provide the best channels for water infiltration. The loess near the failure glide plane was softened by the water immersion, and its deformation increased. The downward displacement of the loess overlying the tunnel increased obviously, and the secondary lining load increased significantly, which negatively impacts the safe operation of the tunnel.

2.2. Failure Glide Planes and Immersion Modes of the Surrounding Rock

The surrounding rock failure glide plane of the loess tunnel widely existed in the shallow-buried loess tunnel. The failure glide plane is one of the main causes of the collapse of the shallow tunnel, the load increasing on lining, and structure damage during the construction period [24, 25]. The influence of the loess immersion did not consider the immersion of failure glide plane in the previous studies. The influence of loess immersion on the loads of the secondary lining in shallow tunnel cannot be ignored. Furthermore, the secondary lining load is an important basis for the prevention of lining structure diseases, and it is of great significance for the maintenance and repair of tunnels.

Under different tunnel excavation methods, two typical types of failure glide planes occur in the surrounding rock of shallowly buried loess tunnels: one type extends from the arch foot to the surface at a certain inclination angle along both sides of the tunnel (full section excavation), whereas the other type extends from the hance to the surface at a certain inclination angle along both sides of the tunnel (bench excavation) [13, 20]. Corresponding to these two types of failure glide planes, there are four modes of water immersion. The two failure glide planes in the shallowly buried tunnel and its water immersion modes are shown in Figure 5.

3. Secondary Development of Structural Duncan-Chang Constitutive Model

3.1. Expression and Formula Fitting of Loess Structural Parameter

Because the parameters of the nonlinear model can be determined by conventional geotechnical tests, the nonlinear constitutive models (the modified Cambridge constitutive model and the Duncan-Chang constitutive model) are the most widely used in engineering applications [26].

The modified Cambridge model is based on the normal consolidation of saturated remolded soil, and it has few parameters and a clear physical definition. However, this model cannot accurately describe over consolidated clay. The Duncan-Chang model is a nonlinear elastic model based on incremental generalized Hooke’s law. It has clear physical definition and is easy to program, and its parameters could be easy to determine. However, this model does not consider the stress history of soil and does not reflect the dilatancy and compression hardness of geotechnical materials. In addition, both modified Cambridge model and Duncan-Chang model cannot accurately describe the stress-strain relationship of structural soil [27, 28].

Loess is a typical structural soil. Establishing a nonlinear constitutive model considering structural properties is an effective way to study the mechanical properties of loess. Considering the structural changes in the loess stress process, the structural parameters are introduced into the Duncan-Chang constitutive model, which is more in line with the stress-strain changes of structural loess.

3.1.1. Expression of Loess Structural Parameter

The difference between the structural Duncan-Chang constitutive model and the Duncan-Chang constitutive model is that the former considers the structure of the loess. Loess is a special type of soil with skeleton and void characteristics and particle cementation between soil particles. The skeleton and void characteristics can be eliminated by remolding and loading, and the cementation between the particles can be eliminated by water immersion. Considering the strength of soil, the loess structural parameter can be expressed by the principal stress difference of undisturbed soil, saturated undisturbed soil, and remolded soil corresponding to the same strain under the triaxial shear condition, which can be used to determine the structural strength of loess.

Under triaxial shear stress conditions, the ratio of the principal stress difference between remolded loess and undisturbed loess and the ratio of the principal stress difference between undisturbed loess and saturated loess can comprehensively reflect the spatial ordering and the connection characteristics of loess grains. The division of these two ratios is the mathematical expression of the loess structural parameter under triaxial shear conditions. The expression is as follows:where is the principal stress difference corresponding to undisturbed loess, saturated undisturbed loess, and remolded loess when i = o, r, and s, respectively.

This expression implies that the stronger the connection of undisturbed loess, the greater the strength loss of remolded loess and the greater the loess structural parameter. Moreover, the greater the structural damage of the flooded loess, the greater the strength loss of undisturbed saturated loess and the greater the loess structural parameter.

3.1.2. Formula Fitting of Loess Structural Parameter

The formula fitting of the loess structural parameter was conducted using the BP neural network model of MATLAB, and the input layer, hidden layer, and output layer were applied. In the input layer, the loess structural parameter and its corresponding factors need to be input, and there are 1114 groups of data; thus, the input vector is a 4 × 1114 matrix. All data are quoted from the study by Xie et al. [4]. The loess structural parameter is correlated with four factors: water content ω, dry density ρd, principal strain ε1, and confining pressure σ3. In the hidden layer, the number of nodes is 5, which is the largest average number of nodes. The output layer of this model is the value of the loess structural parameter . The hyperbolic tangent function tansig is taken as the transfer function among neurons, and the linear function is taken as the transfer function in the output layer.

The expression of the BP network in this model is shown below:where net is the generated BP network objects; is the network input vector, which has 1114 group test levels; T is the network target vector, which has 1114 loess structural parameters ; and trainlm is the Levenberg-Marquardt optimized algorithm.

Finally, the similarity coefficient R2 of this formula is as high as 0.9222, and the formula of the loess structural parameter is as follows:where y1 is the normalized value of water content ω, y2 is the normalized value of principal strain ε1, y3 is the normalized value of the surrounding rock σ3, and y4 is the normalized value of dry density ρd.

3.2. Expression of the Structural Duncan-Chang Constitutive Model

Loess has different structural stress-strain curves under different confining pressures, as shown in Figure 6, and these structural stress-strain curves are the hyperbolic. Similarly, the stress-strain curve of loess also obeys the hyperbola. Based on the derivation process of the Duncan-Zhang constitutive model, the structural Duncan-Zhang constitutive model can be obtained. The hyperbolic fitting can be expressed aswhere a and b are testing constants; σ1 and σ3 are the first and third principal stresses, respectively; ε1 is the principal strain; and m is the loess structural parameter.

USERMAT, which is a subroutine of the constitutive model in ANSYS, is programmed by user programmable features (UPFs) in FORTRAN to achieve the structural Duncan-Chang constitutive model in ANSYS. Based on the given strain increment, the stress increment is calculated to achieve the stress updating (). The Jacobian matrix () is determined by iteration.

The elastic parameter in USERMAT is updated by the subroutine getLam. The mathematical expressions for the elastic modulus and bulk modulus in the loess structural Duncan-Chang constitutive model are as follows:

The different levels of structural stress determine whether the unit is loaded or unloaded, and S is limited from 0 to 0.95.

During the loading process, the structural tangent elastic modulus is given as follows:where k and n are the coefficient and index of the structural elastic modulus, respectively; pa is the standard atmospheric pressure; Rsf is the structural damage ratio; Cs is the loess structural cohesive force; and is the loess structural internal friction angle.

During the unloading process, the structural unloading modulus is given as follows:where Kur and nur are the coefficient and index of the structural unloading modulus, respectively.

The formula of structural bulk modulus Bs is as follows:

Poisson’s ratio for loess is derived from the tangent bulk modulus Bs and the tangent elastic modulus Est. The formula is as follows:

Poisson’s ratios of ideal elastic materials vary from 0 to 0.49; thus, the tangent bulk modulus Bs is 0.33–17 times the value of Est in the program.

The structural Duncan-Chang constitutive model has a few parameters and clear physical meaning, and all parameters can be determined by triaxial shear tests. This constitutive model can describe the stress-strain relationship of loess under complex stress paths.

3.3. Verification of the Structural Duncan-Chang Constitutive Model

The developed structural Duncan-Chang constitutive model is verified by comparing the numerical simulation results from triaxial tests with the calculation results of test constants. The triaxial tests are simulated in Solid185, which is the latest technology unit in ANSYS. The size of the model is φ39.1 × 80 mm, which is the same as the standard specimen in triaxial tests. A vertical displacement constraint is applied to the bottom of this model. The tops of the model loads are graded with an increment of 10 kPa, and confining pressures of 200 kPa, 300 kPa, and 400 kPa are adopted. The model parameters and testing constants are shown in Tables 1 and 2, respectively. The dry density of the loess sample is 1.26 g/cm3, and its water content is 17.8%.

Figure 7 shows the comparison of theoretical solution and numerical solutions of the triaxial test. The theoretical solutions are calculated according to the expression formula of the structural Duncan-Chang constitutive model under the given parameters, and the numerical solutions are calculated according to the numerical simulation of the triaxial test using the structural Duncan-Chang constitutive model. The comparison shows that the simulation solutions are in close agreement with the theoretical solutions, showing that the developed structural Duncan-Chang constitutive model in ANSYS is reasonable and the results are reliable. Additionally, a previous study has proved that the theoretical solutions fit the triaxial test results quite well [7]; therefore, the numerical simulation of the triaxial test using the structural Duncan-Chang constitutive model in ANSYS corresponds with the triaxial test results.

4. Loads on Secondary Lining in a Flooded Loess Tunnel

4.1. Element Selection

Loess is simulated by the structural Duncan-Chang constitutive model in the numerical model. Because the subroutine USERMAT in ANSYS does not support traditional units, plane182 is selected, which is the latest technology unit and has properties such as plasticity, superelasticity, stress stiffening, large deformation, and large strain. Geotextiles and waterproof boards are inserted between the secondary lining and primary support, and the interaction between these elements is represented by the contact unit conta172 and target unit targe169 [29, 30]. Conta172 is set on the inner surface of the primary support, and targe169 is set on the outer surface of the secondary lining.

4.2. Boundary Conditions and Model Parameters

The accuracy of the solution is closely related to the model boundary. Under the same conditions, the bigger the model boundary, the closer the solution will be to the real value. When the boundary range reaches a certain value, expanding the boundary range provides no obvious improvement on the accuracy of the solution. To minimize the adverse effects of the boundary constraint on the calculation results, the boundaries of the calculation model are set as follows: the horizontal width of the surrounding rock on both sides of the tunnel is 35 m; the bottom boundary is 35 m from the tunnel bottom; the depth of the tunnel is 30 m, which is the actual buried depth in monitoring section; and the whole model is 82 m wide and 75 m high. There are horizontal and vertical displacement constraints on the bottom and horizontal displacement constraints on both sides of the model. The top of the model is a free surface. The calculation model is shown in Figure 8. Under surface water infiltration, shallowly buried loess tunnels have three sections immersed in water: the surface, the loess near the failure glide planes of the hance, and the loess near the failure glide planes of the arch foot.

In this model, the secondary lining is constructed of concrete C25 with a thickness of 40 cm, and the primary support is shotcrete with a thickness of 25 cm. The parameters of the lining concrete are shown in Table 3.

According to the geological conditions, the loess in the model can be divided into two layers. The upper layer is Quaternary Upper Pleistocene loess Q3, and the sublayer is Quaternary Lower Pleistocene loess Q2. Table 4 shows the model parameters of the structural Duncan-Chang constitutive of loess Q3 after water immersion and loess Q2 with natural water content.

4.3. Location and Range of Immersion

Large-scale in situ water immersion tests suggest that if loess with a failure glide plane is flooded, the failure glide plane first fills with water, after which water infiltrates vertically and horizontally at its boundary [26]. The water infiltrates downward and outward in loess Q3, and the geometric relationship between infiltration depth Z1 and width Z2 is shown in Figure 9. The calculation formula for infiltration depth and width is shown in equation (8), and the water diffusion angle of loess is a = 31°. It is assumed that the water in the failure glide plane seeps around at the same time. The water infiltration times are 12 h, 24 h, 48 h, and 72 h, and the water infiltration widths are 1.2 m, 1.96 m, 3.3 m, and 4.54 m, respectively. Assuming that the water infiltration widths (2 × Z2) of the failure glide planes are 1 m, 2 m, 3 m, 4 m, and 5 m, the relative water infiltration depths Z1 (surface infiltration depth) are 0.83 m, 1.66 m, 2.49 m, 3.32 m, and 4.16 m, respectively. The horizontal and vertical immersion ranges of the surface and the failure glide planes are shown in Table 5.

Here, is the infiltration time and its unit is h, Z1 is the water infiltration depth, and Z2 is the water infiltration width. The units of both Z1 and Z2 are m.

4.4. Calculation Cases

Considering the immersion modes and the immersion range, the calculation cases are divided into four cases: loess immersion near the failure glide planes of both arch feet (Case 1), loess immersion near the failure glide plane of one arch foot (Case 2), loess immersion near the failure glide planes of both hances (Case 3), and loess immersion near the failure glide plane of one hance (Case 4). Case 1 corresponds to the immersion conditions of the project case. All calculation cases are shown in Table 6 in detail, and the immersion locations in different cases are shown in Figure 8.

5. Calculation Results and Analysis

The load on the secondary lining is the contact pressure between the secondary lining and the primary support, which is monitored during the immersion process. There are 15 monitoring points distributed symmetrically on the outside of the lining at the vault, arch shoulder, hance, arch end, sidewall, arch foot, and inverted arch, as shown in Figure 10.

5.1. Variation of the Loads on the Secondary Lining

Figures 11 and 12 show the curves of loads on the secondary lining in cases 1 and 2, in which the loess near the failure glide plane of the arch foot is flooded. With no immersion (infiltration width is 0 m), loads from the hance to the vault are uniformly distributed, and loads from the sidewall to the inverted arch are less than those on the arch. When the loess near the failure glide plane is flooded, the loads on the secondary lining increase. The loads on the arch foot and hance increase faster than those on the other positions on the immersed side, and the load on the arch foot increases to the greatest extent. The distribution of loads on the arch is shaped like a cat’s ears, which is in accordance with the measured results in previous references [31, 32], and the changes in the load on the vault are not obvious. The load on the immersed side increases more than that on the nonimmersed side.

Loess immersion on one arch foot may cause the load to increase on both sides of the arch foot. The reason for the increased load on the arch foot of the nonimmersed side may be that the increasing loads on the immersed side form a biased pressure and the lining deforms toward the nonimmersed side, leading to an increased stratum resistance on the arch foot of the immersed side. As the immersion range expands, the maximum load on the arch on the immersed side moves upward from the arch end to the hance.

Figures 13 and 14 show the load curves on the secondary lining in Cases 3 and 4, in which the loess near the failure glide plane of the hance is flooded. In general, when the loess is flooded, the loads on the secondary lining increase. The loads are relatively large at three positions on the immersed side: the arch shoulder, the arch end, and the arch foot. The loads on the arch shoulder and arch foot increase obviously. The load on the vault increases continuously, which is different from Case 1 and Case 2, where the loads undergo little change.

The loads on the immersed side increase much more than those on the nonimmersed side, except the load on the arch foot on the nonimmersed side, which is larger than that on the immersed side, opposite to Case 2. There are two explanations for this phenomenon: One is the “bias pressure,” which was discussed for Case 2. The other is that the increasing loads on the arch shoulder have a significant effect on the increased resistance of the arch foot. The increasing stratum resistance causes the load on the arch foot to increase. During the immersion process, the maximum load on the arch first moves upward from the arch end to the arch shoulder and then moves back to the arch end.

5.2. Relationship between the Infiltration Width and the Load

Figure 15 shows the load variation curves on the sidewall and arch end in Cases 1 and 2. The loads on the sidewall and arch end decrease at first and then increase, and the load on the sidewall decreases earlier than that on the arch end. The load on the arch end in Case 2 decreases earlier than that in Case 1. In other words, when the immersion is close to the load monitoring point, the load decreases at first and then increases slowly as immersion continues. The reason for this phenomenon may be that the loess softens when the water arrives, and the softened loess releases some of the load on the lining. As immersion continues, the stability of the overlying loess is weakened and the loess stratum deforms downward, which causes the loads to increase.

Figure 16 represents the load variation curves from the hance to the vault in Cases 3 and 4. The loads on the hance, the arch shoulder, and the arch end decrease at first and then increase, except for the load on the arch shoulder in Case 4. The load on the hance decreases first, and then the load on the arch end decreases. In addition, in the course of water immersion, the loads on the monitoring point on the upper positions of the failure glide planes are more affected by the immersion than the loads on the monitoring points on the lower positions of the failure glide planes.

5.3. Comparison between Calculated Results and Measured Results

A comparison between the calculated results and the measured results for the secondary lining loads is shown in Figure 17. The red line indicates the measured loads on the secondary lining in the shallow loess tunnel buried 30 m deep. The blue line indicates calculated loads on the secondary lining in Case 1, in which the infiltration horizontal width is 3 m. The calculated loads and measured loads on the secondary lining are generally symmetrically distributed and basically have the same regularities of distribution. Loads on the sidewalls and arch end are similar.

The maximum loads are located at the arch foot on both sides, and the maximum difference load on the left arch foot is 474 kPa. The reason for this difference may be that the pressure on the primary support is relieved before the construction of the secondary lining, and normally, the pressure relief cannot be simulated well in numerical calculations. Thus, the measured loads on the secondary lining are smaller than the calculated loads.

6. Conclusions

(1)Finite element calculation of the structural Duncan-Zhang constitutive model is realized in ANSYS and proved to be reliable by comparing the numerical results of triaxial tests and the calculation results of testing constants.(2)The variation and distribution of the load on the secondary lining are closely related to the immersion position and range:(a)When the loess near the failure glide planes of the arch foot is immersed (Case 1 and Case 2), the strength of the flooded loess decreases, leading to a sharply increased load on the arch foot. As the immersion expands, the loads on the sidewall and arch end first dissipate partially and then increase gradually. The maximum load on the arch moves upward from the arch end to the hance. The load increase on the immersed side results in bias pressure, and the load on the arch foot on the nonimmersed side increases significantly.(b)When the loess near the failure glide planes of the hance is immersed (Case 3 and Case 4), the load on the hance dissipates slightly. The stability of the overlying loess decreases gradually, causing a general load increase on the upper lining. The loads on the vault and arch shoulder increase rapidly and work downward on the lining. Therefore, the load on the arch foot increases against the stratum resistance. The maximum load on the arch moves upward from the arch end to the arch shoulder and then down to the arch end.(3)The load on the arch foot lining is always large after the surrounding rock is flooded, and the load on the arch foot when the loess is immersion near the arch foot is always larger than that in cases of loess immersion near the hance. When the immersed area is close to the load monitoring point, the load on this point decreases at first and then increases slowly.

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 potential conflicts of interest.

Acknowledgments

This study was supported by the National Natural Science Foundation of China (grant nos. 51908051 and 51978064), the Fundamental Research Funds for the Central Universities, CHD (grant no. 300102289101), and the Traffic Construction Research Funds of Shaanxi Province (grant no. 2016-1-3).