Abstract
Determination of the required supporting pressure is the premise of tunnel support design. Only when the support design meets the requirements can the tunnel be safe and stable during construction and operation. This paper focuses on a shallow tunnel in layered rock strata and proposes a method for predicting the required supporting pressure. In this method, the 2D and 3D failure mechanisms are constructed, respectively. The analytical solutions of the supporting pressure corresponding to the two cases are derived on the basis of upper bound theorem and Hoek–Brown failure criterion. Then, the proposed method is validated by comparing with the results of existing research studies. Furthermore, a shallow tunnel in two-layer rock strata is chosen to illustrate the difference between the two solutions. The comparison shows that the supporting pressure in the 2D case is greater than that in the 3D case in general and it tends to be conservative for tunnel design. Conversely, the 3D solution may help to reduce the support cost. Furthermore, the change laws of the supporting pressure and failure range corresponding to varying parameters are obtained. These results may practically provide theoretical references for tunnel support design in layered rock strata.
1. Introduction
Now, as problems including heavy ground traffic and shortage of land resources have become increasingly acute in cities, construction of tunnels or other underground works has become a mainstream trend of urban development in the world. According to the burial depth, tunnels fall into two categories: shallow tunnels and deep tunnels. As for shallow tunnels, due to thin overlying rock strata, tunnel excavation has a great effect on ground disturbance; subsidence failure of surrounding rock masses is obvious. Especially when the geological conditions of tunneling strata are poor or in adverse conditions such as underground water, fault fracture zone, soft sand layer, and additional ground load, stability control over surrounding rock masses becomes more difficult. Once the support is unreasonable, disasters such as collapse are liable to occur, thus bringing a huge impact on engineering construction and operation safety.
As a hotspot, the researchers have carried out numerous research studies on the failure mechanism of shallow tunnels. In general, the following three research approaches can be employed: numerical simulation, theoretical analysis, and experimental method. For example, Karakouzian et al. [1] and Karami et al. [2] utilized the finite element method to investigate the effects of overburden height and internal transient pressure on the hydraulic fracturing of a concrete-lined pressure tunnel. The upper bound method is a classical theoretical analysis method. It has gained great recognitions and extensive applications. In this method, a kinematically admissible velocity field for tunnel failure is required to be built in advance, and the limit load can be solved on the basis of the energy balance principle. By doing so, the complex and tedious calculation can be simplified effectively. Then, a close-to-actual failure mechanism can be obtained. For example, Leca and Dormieux [3] proposed three failure mechanisms of the tunnel face based on the movement of rigid conical blocks and obtained three upper bound solutions of the face pressure for a shallow tunnel driven in a frictional material. Lee and Nam [4] incorporated the influence of seepage force into the upper bound limit analysis of the tunnel face stability. Yamamoto et al. [5] focused on the stability of a square tunnel in cohesive-frictional soils subjected to surcharge loading and proposed the upper bound solutions for the ultimate surcharge loading by applying numerical limit analysis techniques. Mollon et al. [6] proposed two new continuous velocity fields for both collapse and blowout of an air-pressurized tunnel face in a purely cohesive soil and obtained the expressions of the collapse pressure using the upper bound method. Huang and Song [7] investigated the undrained stability of a plane strain tunnel heading in cohesive soil based on a multi-rigid-block mechanism. Zhang et al. [8] investigated the 3D failure characteristics of shallow circular tunnel faces in cohesive-frictional soils. Han et al. [9] and Li et al. [10, 11] investigated the influence of nonhomogeneity and anisotropy of soil masses on the face stability of shallow tunnels. In the above-mentioned works, the soils or rocks in the strata are all regarded as the traditional Mohr–Coulomb materials.
In view of the nonlinear characteristics of the strength envelops when the soils or rocks fail [12–15], the nonlinear failure criteria may be more suitable to describe the failure of the tunnel’s surrounding rock masses. Fraldi and Guarracino [16, 17] initially introduced the nonlinear Hoek–Brown failure criterion to analyze the roof collapse mechanisms for deep tunnels based on upper bound theorem. Recently, they [18] proposed a general characterization of tunnels depth and analyzed the collapse characteristics of intermediate tunnels. Based on their method, Huang and Yang [19] and Zhang et al. [20] conducted upper bound limit analysis on the roof collapse of deep circular tunnels. Furthermore, Yang and Huang [21], Huang et al. [22], and Guan et al. [23] proposed the three-dimensional collapse mechanisms of deep cavities in homogeneous Hoek–Brown rock masses. Qin et al. [24] and Yang et al. [25, 26] incorporated the nonhomogeneity and stratification characteristics of soil masses or rock masses into two-dimensional collapse analysis of deep tunnels. Similarly, they [27–30] also incorporated this influence into the three-dimensional mechanisms of deep tunnels or cavities. In the field of shallow tunnel collapse, Yang and Huang [31] proposed the collapse mechanism for a shallow circular tunnel. Yang and Li [32] conducted the upper bound limit analysis on roof collapse for a shallow tunnel in two-layered rock strata. Wang et al. [33, 34] proposed two kinds of collapse mechanisms for a shallow tunnel by incorporating the effects of changing groundwater table and pore water pressure. Lyu and Zeng [35] investigated the collapse of shallow tunnels in inclined rock stratum.
It is worth noticing that all the abovementioned research studies mainly focused on the prediction of tunnel collapse scope and mode. Specifically, in these works, the objective function constructed on the basis of the internal energy dissipation rate and the work rates done by external forces is the total energy dissipation rate. Moreover, the investigation on the shallow tunnel is relatively less on the whole. It should be noted that, due to smaller thickness of upper covering rock masses, the shallow tunnel is difficult to form an effective bearing arch structure inside the roof surrounding rocks and to liberate the self-bearing function of rock masses. Thus the determination of the required supporting pressure is of great significance to ensure the stability of shallow tunnels. Huang et al. [36] and Jiang et al. [37] once constructed the objective function of the supporting pressure and derived the analytical expressions of the required supporting pressure, but the tunnel’s surrounding rock masses were assumed to be homogeneous and the proposed mechanisms were only limited to 2D case. The required roof supporting pressure for a shallow tunnel in layered rock strata and in 3D failure case needs further investigation. Consequently, based on the preceding research works, a shallow tunnel with arbitrary cross section in layered rock strata is chosen for study in this paper, and the 2D and 3D failure cases are considered simultaneously. The required supporting pressure corresponding to the 2D and 3D failure mechanisms is determined on the basis of the upper bound theorem. The effect of the thickness of varying rock strata, rock strength parameters, the additional ground load, etc. on the supporting pressure and failure range are also analyzed. The research results may serve as theoretical guidance for tunnel support design and construction in layered rock strata.
2. 2D and 3D Roof Failure Mechanisms for a Shallow Tunnel in Layered Rock Strata
A tunnel is a long and narrow underground structure. Therefore, it can be regarded as a plane strain problem to analyze the stress and deformation of the tunnel’s surrounding rock masses. However, in actual engineering (Figure 1), tunnel collapse usually involves a typical 3D failure problem; the roof failure range is limited along the length of the tunnel. In addition, the thickness of the overlying rock strata is small for shallow tunnels; the roof failure often extends to the Earth’s surface; as a result, a collapse arch bearing structure cannot be formed inside the surrounding rock masses; this is also an important reason why stability control over shallow tunnels is difficult. According to the abovementioned features, we propose the 2D and 3D failure mechanisms for a shallow rectangular tunnel in layered rock strata, respectively, as shown in Figure 2.
(a)
(b)
(a)
(b)
Specifically, in Figure 2, the tunnel burial depth is H; the roof includes n layers of rock masses; the thickness of each layer is (i = 1, 2, … n). In Figure 2(a), we assume that roof failure occurs in the xoy plane; the corresponding failure curve (i = 1, 2, … n) consists of n segments: , , … for layer 1, layer 2, ... layer n, respectively. Accordingly, in Figure 2(b), we assume that the rock masses within the failure range constitute a 3D axisymmetric rotating body, which may be formed by rotating the curve in Figure 2(a) around the y axis within the xoy plane; the corresponding 3D failure surface equation for layer i is . Meanwhile, the roof failure is also affected by both the additional ground load and the supporting pressure .
Furthermore, according to the existing indoor rock mechanics test results, in the plane and the plane, the strength envelope of the rock should be nonlinear rather than the linear relation used in the traditional Mohr–Coulomb strength criterion; a nonlinear strength criterion should be more suitable to describe the nonlinear failure characteristic of the tunnel’s surrounding rock masses. Consequently, the widely used Hoek–Brown strength criterion [12, 13] is introduced in this paper. In the Mohr plane , it can be expressed as
In equation (1), and are dimensionless empirical parameters related to properties of rock masses; and refer to the shear stress and normal stress at the rock fracture surface; and and refer to the compressive strength and tensile strength of the rock mass. Meanwhile, we assume that the roof rock masses are ideal rigid-plastic; failure of the rock masses meets the Hoek–Brown strength criterion and its associated flow rules. Accordingly, we may set the corresponding yield function and the plastic potential function of rock failure to be equal; then, according to the potential theory, the plastic strain rates produced by rock failure in layer i may be solved by the following equation:where refers to the plasticity factor; refers to the plasticity strain rate component; and and refers to the stress component. The yield function corresponding to the Hoek–Brown strength criterion can be expressed as follows:
3. 2D Limit Analysis on the Required Supporting Pressure for a Shallow Tunnel
According to Chen’s [38] theory, if the upper bound method is used to solve the required supporting pressure for a shallow tunnel, a kinematically admissible velocity field satisfying the deformation coordination needs to be established in advance; after the internal energy dissipation rate and external force power corresponding to the velocity field are obtained, the virtual work-rate equation can be constructed accordingly. Specifically, according to the failure mechanism shown in Figure 2(a), we may assume that the rock masses within the failure range are collapsing downward at the speed of under the force of gravity and additional ground load. Conversely, the rock masses in the area that is not collapsing stay static. In addition, the rock masses are assumed to be ideal rigid-plastic. The collapsing rock masses and the surrounding static rock masses can be considered to be rigid, while the rock masses at the failure surfaces are in the plastic flow state. Consequently, internal energy dissipation occurs only at the roof failure surfaces.
Accordingly, the failure curve in layer i is a velocity discontinuity line. We may obtain the plastic shear strain rate and plastic normal strain rate at the failure surface by substituting equation (3) into equation (2). Thus, the corresponding internal energy dissipation rate in layer i can be expressed as follows:where refers to the thickness of the failure surface in layer i.
Considering the symmetry of the 2D failure mechanism, half of the failure range is selected for study. By integrating equation (4) along the n failure surfaces in all rock layers, we obtain the equation for the total rate of internal energy dissipation:where (i = 1, 2, …, n) refers to the half width of the failure range in layer i.
The work rate done by the weight of the collapsing rock masses iswhere is the unit weight of rock layer i.
The work rate produced by the additional ground load is
The work rate produced by the supporting pressure is
Based on equations (6)–(8), we can obtain the total work rate done by external forces. Then, by using the virtual work-rate principle, we obtain
Substituting equations (5)–(8) into equation (9) yields
According to equation (10), the required roof supporting pressure for a shallow tunnel can be expressed aswhere is
Equation (11) is a functional of x and . According to the upper bound theorem, in all potential roof failure mechanisms, the corresponding to the real mechanism shall enable equation (11) to obtain the maximum value. This is a typical variational problem. The following Euler–Lagrange equation can be used to obtain the optimal upper bound solution:
Substituting equation (12) into equation (13) results in
By integrating equation (14), we obtainwhere is the constant to be determined. According to the progressive failure characteristics of the roof’s surrounding rock masses, the failure surface should be first formed in the bottom, continue to develop upward, and eventually extend to the Earth’s surface. Similarly, in layer i, rock failure also develops upward from the bottom of this layer. In addition, the layer i is internally homogeneous. If the rock failure in layer i does not extend to the upper layer, this is very similar to the failure characteristics in homogeneous formation. Furthermore, according to the symmetry of the 2D failure mechanism in Figure 2(a) and the failure characteristics in homogeneous rock media [16, 17], the first derivative of in layer i is set equal to 0 at the point ; then parameter can be determined as 0 and equation (15) can be expressed aswhere .
By integrating equation (16), we obtain
Furthermore, by substituting equations (16) and (17) into equation (12), we obtain
Equation (18) is substituted into equation (11) to obtain the required roof supporting pressure:
Note that the required supporting pressure involved in equation (19) is a function related to the width of the failure range and increases with the increase of . To ensure the safety of engineering design, the worst case of roof failure shall be considered. Specifically, we may set the bottom width of the failure range to be equal to the tunnel width; in this case, the tunnel requires the maximum supporting pressure, which can be obtained bywhere (half width of the tunnel) shall meet . In equation (20), and (i = 1, 2, 3, … n) are still unknown. According to the geometrical relationships in Figure 2(a), the failure curve meets
By substituting equation (17) into equation (21), we obtain the integration constant :
Furthermore, the following equations can be derived from equation (21).
Based on equations (21)–(23), the width and constant can be solved. Then the required supporting pressure and the failure range for the shallow tunnel in the 2D failure case can be determined accordingly.
4. 3D Limit Analysis on the Required Supporting Pressure for a Shallow Tunnel
As shown in Figure 2(b), in the case of 3D roof failure for a shallow tunnel, the specific derivation process of the required supporting pressure is in accordance with that in Section 3. Similarly, the rotational collapsing rock masses are assumed to be a rigid block with the downward velocity of . The failure surface in layer i is approximately conical. The plastic strain rates can also be solved by using equations (2) and (3). Consequently, by integrating equation (4) along the 3D failure surfaces, we obtain the equation for calculating the total rate of internal energy dissipation:where (i = 1, 2, …, n) refers to the radius of the rotational failure block in layer i, and it is similar to the half width of the failure range in the 2D case.
The work rate done by the weight of collapsing rock masses is
The work rate produced by the additional ground load is
The work rate produced by the supporting pressure is
Based on equations (24)–(27), the virtual work-rate principle can be utilized to obtain the following equation:
Then, the required supporting pressure in the 3D case can be solved aswhere is
Similarly, to obtain the optimal upper bound solution, we can also use the Euler–Lagrange equation in equation (13). Thus, by substituting equation (30) into equation (13), we obtain
By integrating equation (31), we obtainwhere is an unknown constant.
Note that it is very complicated and difficult to solve equation (32) to obtain the analytical expression of . To facilitate engineering design, the handling method of equation (15) in Section 3 can be utilized to solve this problem. Specifically, the failure characteristics in homogeneous strata and the symmetry of the 3D rotational failure surfaces shall be taken into consideration. The first derivative of can also be assumed to be zero at the point ; then the constant here can be determined as zero. Furthermore, equation (32) can be expressed aswhere .
Furthermore, by integrating equation (33), we obtainwhere is another unknown constant.
Based on equation (34), the equation of the 3D failure surface can be expressed as
Substituting equations (33) and (34) into equation (30) results in
Furthermore, substituting equation (36) into equation (29) results in
Similarly, to ensure the safety of engineering design, according to the design principle for the worst case in Section 3, we set the bottom width of the failure block to be equal to the tunnel width. Then the required supporting pressure corresponding to the 3D mechanism in Figure 2(b) can be expressed as
The unknown parameters can be solved by using the geometrical relationships in Figure 2(b). Particularly, the curve (i = 1, 2, 3, … n) meets the following equations.
By substituting equation (34) into equation (39), we obtain
Furthermore, the following equations can also be derived from equation (39):
By combining equations (39)–(41), we can determine the constant and the radius of the 3D failure block. Then the required supporting pressure and the failure range of the shallow tunnel in the 3D case can be determined.
5. Results and Discussion
5.1. Comparison with Existing Research Studies
In order to predict the required supporting pressure for shallow tunnels, Huang et al. [36] proposed the failure mechanism of a rectangular tunnel in homogeneous rock media and derived a 2D analytical solution of supporting pressure. Accordingly, in order to validate the effectiveness of the proposed method in this paper, we set the strength parameters of all rock layers in Figure 2(a) to be the same value; now this mechanism can be turned into the mechanism by Huang et al. [36]. Furthermore, Figure 3 lists the comparison results between the required supporting pressure calculated by this paper and that calculated by Huang et al. [36]. The comparison involves the following parameters: and , , , , , , and . It can be seen from Figure 3 that the required supporting pressure increases with the burial depth and span of a shallow tunnel. Meanwhile, the results calculated by this paper are slightly smaller than those calculated by Huang et al. [36], but the maximum difference is not more than 8%. Thus, the proposed method in this paper is validated to be effective.
5.2. Comparison of the Required Supporting Pressure in 2D and 3D Mechanisms
According to the above theoretical derivation, we obtain the analytical solutions of the required supporting pressure in the 2D and 3D cases, respectively, which can provide theoretical reference for tunnel support design in practice. However, the difference between these two solutions remains elusive. To solve this problem, a shallow rectangular tunnel in two-layer rock strata is introduced as an example, and a comparative analysis is performed in this section. Figure 4 lists the required supporting pressure in the 2D and 3D mechanisms corresponding to varying tunnel burial depths, wherein the half width of the tunnel is 4 m, 5 m, and 6 m, respectively. The additional ground load is 20 kPa. The strength parameters of the upper and lower rock layers are as follows: , , , , ; , , , and .
It can be seen from Figure 4 that the calculation results corresponding to the 2D mechanism are greater than those corresponding to the 3D mechanism. In addition, when the tunnel burial depth is small or as the span increases, the difference between the 2D results and the 3D results decreases. This shows that the 2D failure mechanism based on the traditional plane strain hypothesis tends to overestimate the required supporting pressure and is conservative in tunnel support design. On the contrary, it is relatively economic to use the 3D failure mechanism to design the roof support. Therefore, in actual engineering, on the premise that the tunnel roof safety is guaranteed, a better support cost can be obtained with the 3D mechanism proposed in this paper. This can provide certain theoretical guidance for support design of shallow tunnels in layered strata.
5.3. Comparative Analysis for a Two-Layer Tunnel Roof
In order to better guide tunnel support design, in-depth analysis on effect laws of varying parameters on the roof supporting pressure is an issue of great significance. This is because the tunnel support design can be more accurate in actual engineering only after the key factors influencing the supporting pressure are identified, thus avoiding waste of support costs and guaranteeing the safety of the support structure. In this section, the supporting pressure corresponding to the 3D failure mechanism in two roof rock layers is chosen for analysis. Figure 5 shows the change laws of the supporting pressure under varying parameters, wherein the tunnel burial depth is ; the span is ; the thickness of the upper layer is ; the thickness of the corresponding lower layer is ; and the additional ground load is . The strength parameters of the upper and lower layers are as follows: , , , , ; , , , , and . Note that the other parameters are fixed when the change of a parameter is analyzed.
(a)
(b)
(c)
(d)
(e)
(f)
(g)
(h)
As can be seen in Figure 5, when the tunnel burial depth is certain, the overall strength of the roof decreases and the required supporting pressure increases as the thickness of the upper rock strata increases; this is because the strength of the lower rock layer is greater than the strength of the upper rock layer. Furthermore, the required roof supporting pressure is positively correlated with the additional ground load, empirical parameter B, and unit weight; it is negatively correlated with empirical parameter A, compressive strength, and tensile strength of rock mass. Thus, in actual engineering, when the field geological conditions are certain, the adverse effects of the additional ground load should be minimized as far as possible. Meanwhile, the advance grouting technology should be adopted to improve the strength parameters of the tunnel’s surrounding rock masses [39–43]. These are effective measures for reducing the required supporting pressure and support costs.
Furthermore, we plot the failure curves corresponding to the varying parameters in Figure 5, as shown in Figure 6. As can be seen in Figure 6, when the thickness of the lower layer increases, the failure range tends to decrease as the required supporting pressure decreases. As the parameter A and the compressive strength of the upper and lower layers increase, the failure range also decreases. However, as the tunnel span, parameter B, and unit weight increase, the failure range increases. In addition, the effect of parameter B, rock compressive strength, and unit weight on the failure range is less significant, while the effect of parameter A and tunnel span is more significant. Priority should be given to parameter A and tunnel span in actual engineering design.
(a)
(b)
(c)
(d)
(e)
(f)
5.4. Comparative Analysis for a Nonhomogeneous Tunnel Roof along Depth
In practical engineering, due to the influence of long-term compression from the overlying strata, the strata that the tunnel passes through may present obvious nonhomogeneous characteristics along depth. This means the rock parameters vary with depth. For example, Fraldi et al. [18] considered the profundity of the excavation and the variability of the rock mass parameters and proposed a general characterization of tunnels depth. In their research works, the unit weight , parameter , and rock compressive strength are defined as the functions of the depth and tend to increase along depth. This nonhomogeneity can exert significant influence on roof collapse and safety design of tunnel support. With reference to [18], in order to incorporate this influence into the 2D and 3D layered collapse mechanisms proposed in Figure 2, we can equally divide the tunnel roof into n layers, as shown in Figure 7. Moreover, the parameters , , and in Layer 1 and those in Layer n are assumed to obey the following relationships:where , , and are nondimensional coefficients characterizing the nonhomogeneity of roof rock masses and can be obtained through laboratory test or field test. Meanwhile, for convenience, we further assume that the rock parameters in varying layers are variable and increase linearly with the depth, while the rock mass in one layer is considered locally homogeneous.
Then in Layer i, the corresponding parameters , , and can be expressed as
Furthermore, by utilizing equation (20) and equation (38), the 2D and 3D analytical solutions of the required supporting pressure for a nonhomogeneous tunnel roof can be obtained. In order to clarify the influence of the nonhomogeneity of rock masses, the supporting pressures corresponding to varying coefficients , , and under the 3D case are calculated and are listed in Figures 8–10, wherein the tunnel burial depth is 5 m, and the roof is divided into five layers. The rock parameters in Layer 1 are , , , , and . When one of the parameters varies in varying layers, the other parameters remain constant.
(a)
(b)
(a)
(b)
(a)
(b)
As can be seen in Figures 8–10, the nonhomogeneity of the rock parameters can significantly affect the required supporting pressure. When the coefficients , , and are equal to 1, the rock parameters do not vary along the depth, and the tunnel roof is homogeneous. As the coefficient increases, the weight of the collapsed rock masses increases and the required supporting pressure tends to increase accordingly. Conversely, as the coefficients and increase, the whole strength of the tunnel roof increases and the required supporting pressure tends to decrease, wherein the coefficient has the most significant influence. For example, compared to , , and , the required supporting pressure corresponding to , , and changes by 34.5%, 11.3%, and 1.38%, respectively. Thus, in tunnel support design, the nonhomogeneity of rock unit weight should be paid much care in real sceneries.
6. Conclusions
This paper constructs the 2D and 3D failure mechanisms for a shallow tunnel in layered strata. The analytical solutions of the required supporting pressure corresponding to 2D and 3D cases are derived with upper bound method. Furthermore, the effectiveness of the proposed method is validated by comparison with existing research works. The difference between the two solutions corresponding to the 2D and 3D mechanisms is clarified. The effect laws of varying parameters on the roof supporting pressure and failure range are obtained. This paper mainly draws the following conclusions:(1)The required supporting pressure corresponding to the 2D mechanism is greater than that corresponding to the 3D mechanism. In actual engineering, on the premise that the tunnel roof safety is guaranteed, a better support cost can be obtained with the 3D failure mechanism proposed in this paper.(2)The required supporting pressure is positively correlated with the tunnel burial depth, span, additional ground load, empirical parameter B, and unit weight; it is negatively correlated with the thickness of the rock layer with higher strength, empirical parameter A, rock compressive strength, and tensile strength. In actual engineering, the adverse effects of the additional ground load should be minimized as far as possible or the advance grouting technology should be adopted to improve the strength parameters of surrounding rock masses. These are effective measures for reducing the required supporting pressure and support costs.(3)The tunnel failure range corresponding to the required supporting pressure is negatively correlated with the thickness of the rock layer with higher strength, parameter A, and compressive strength; it is positively correlated with the tunnel span, parameter B, and unit weight. The effect of parameter B, compressive strength, and unit weight on the failure range is less significant, while the effect of parameter A and tunnel span is more significant. Priority should be given to parameter A and tunnel span in actual engineering design.(4)The proposed method in this paper is employed to predict the required supporting pressure for a nonhomogeneous tunnel roof along depth. The results show that the required supporting pressure is positively correlated with the nonhomogeneity coefficient of rock unit weight and is negatively correlated with the nonhomogeneity coefficients of parameter A and compressive strength.
Data Availability
The data used to support the findings of this study are included within the article.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
The authors would like to acknowledge the financial support from the National Natural Science Foundation of China (Nos. 51704177 and 51809159), a Project of Shandong Province Higher Educational Science and Technology Program (No. J16LG04), Shandong Co-Innovation Center for Disaster Prevention and Mitigation of Civil Structures (No. XTP201911), and the Doctoral Research Fund of Shandong Jianzhu University (No. XNBS1501).