Abstract

The development of underground space is fast because of the lack of space. To build shafts in underwater tunnels, the vertical tunneling method (VTM) was invented in 1974 in China, which can act as a freshwater intake or sewage outlet. During the operation of the VTM, the jacking force is one of the essential factors that draw attention. This paper conducts a numerical study of the jacking force and its influencing factors during the vertical tunneling process. First, based on the finite element software ABAQUS, a numerical model of the vertical tunneling process is established according to the VTM project in Beihai, China. Second, in accordance with the Latin hypercube sampling method and the multivariate significance analysis, the mechanical parameters are determined or back-analyzed. Then, the calculated jacking force of the numerical model is compared with the measured jacking force. It turns out that the changing trend of the jacking force in the numerical model and the measured jacking force is relatively consistent. Finally, the influencing factors of the jacking force, such as elastic modulus and the angle of internal friction, are analyzed based on the numerical model. The results show that the elastic modulus and the angle of internal friction of soil are the main influencing factors of the jacking force. The secondary factors are Poisson’s ratio, static earth pressure coefficient, unit weight, and cohesion.

1. Introduction

Due to the world’s population increase and the lack of space, not only the aboveground space is used efficiently, such as high-speed railways [1], but also the underground space has witnessed a development boom in recent years [2, 3]. To build freshwater intake or sewage outlet shafts in underwater tunnels, the vertical tunneling method (VTM) was invented in 1974 in China [47].

There are several critical construction steps of the VTM. The schematic diagram of the construction process and the on-site construction diagram are shown in Figure 1. First, the ceiling segment and the first shaft section are pre-embedded during the construction of the horizontal tunnel. Second, the second shaft section is welded with the first shaft section. Third, all the ceiling segments, the first and second shaft sections, are jacked upward. Then, the third shaft section, the fourth shaft section, etc., are welded with the presections and jacked upward until the entire shaft sections are jacked. Since the ceiling segment has covered the top of the shaft, the shaft is jacked upward by the compression of the above soil, i.e., soil displacement, not soil removal. At last, the ceiling segment is replaced with a one-way valve by divers underwater, and the shaft is constructed successfully. The detailed construction procedure can be found in the previous research [6].

Yang et al. [8] simulated the construction process of vertical tunneling using a particle flow program. The overlying soil layer’s damage mechanism, scope, and development path during the vertical tunneling process are compared with the laboratory test. Xu et al. [9] used theoretical calculation combined with a three-dimensional numerical method to study the structural deformation of horizontal tunnels induced by vertical tunneling. The effect of the jacking force on the deformation characteristics of segments and joints is discussed.

The implementation of vertical tunneling is achieved by compressing the soil above the shaft to cause soil displacement and then uplifting. The principle is similar to that of the penetration of the static pressure pile. Thus, the uplifting of the shaft can be regarded as a reverse piling process [5]. When using numerical models to simulate the penetration process of the static pressure pile in soil, there are mainly three methods (see Figure 1) [10]. The first is the pressure method, the second is the displacement penetration method, and the third is the cavity expansion method.

The pressure method applies pressure directly to the top of the pile. When the penetration distance is small and pressure is too enormous, the calculation may be abnormal. When the pressure is too small, it cannot penetrate and will cause plastic strain. Therefore, determining the appropriate pressure is the key to the operation of the pressure method.

The displacement penetration method means that the pile is pre-embedded at a certain depth, and then, the pile is directly penetrated down to a smaller depth. Pile penetration is achieved by applying a displacement boundary condition to the top of the pile, and no external force is required.

Cavity expansion theory has been widely applied to in situ soil testing, deep foundations, tunnels, underground excavations in soil and rock, and wellbore instability in the oil industry. The cavity expansion theory was initially applied to metal processing. It was later extended to pressure test analysis, bearing capacity analysis of deep foundations, and soil disturbance analysis due to pile driving in geotechnical engineering. The study of pile installation in soil is a large strain problem that involves strong material and geometric nonlinearities. It can predict piles’ end bearing and shaft capacities in soil and rock [11].

It should be noted that the shaft is built through soil displacement by the VTM (see Figure 2). As a result, one of the essential factors during the VTM construction is the jacking force. If the jacking force is too large, the horizontal tunnel may be affected, and horizontal lining segments may crack, even resulting in lining damage. On the other hand, if the jacking force is too small, the shaft cannot be uplifted successfully. During the vertical tunneling process, the soil above the shaft is squeezed until it cannot be compressed anymore. At this time, soil reaches the ultimate state, the jacking force increases to maximum, and soil slippage appears. After that, the jacking force decreases because the soil’s primary state is disturbed.

The jacking force is influenced by ground conditions, shaft materials, jacking speed, etc. The research on the jacking force during vertical tunneling is insufficient; however, some researchers in China have already researched the changing trend of the jacking force. For example, Wang and Ge [12] believed that the jacking force decreases with an increase in the jacking distance. Chen and Huang [13] thought that the jacking force increases in the first place but then decreases. Wang et al. [6] concluded that the changing trend of the jacking force increases to the peak value when the jacking distance reaches one-third of the total jacking distance and then decreases to the minimum value. Thus, it can be summarized that the conclusions are different so far. Therefore, it is vital to figure out the changing trend of the jacking force and its influencing factors during the vertical tunneling process.

To further study the jacking force and its influencing factors, and based on the finite element software ABAQUS, this paper establishes a numerical model for the vertical tunneling process according to the vertical tunneling project in Beihai, China. According to the Latin hypercube sampling method and the multivariate significance analysis, the mechanical parameters of the numerical model are determined or back-analyzed. Moreover, the calculated jacking force of the model is compared with the measured jacking force, and the influencing factors of the jacking force are analyzed.

2. Structure of Vertical Shaft Sections

The structure of the ceiling segment and the first shaft section are shown in Figures 3 and 4, respectively. The corresponding dimensions are shown in Tables 1 and 2. The structure of the ceiling segment is shown in Figure 3. It consists of 1 round steel plate with 12 Φ24 bolt holes, 1 round pipe, four stiffened panels ①, and four stiffened panels ②. The stiffened panels are installed to support the surface of the ceiling segment and the soil pressure above it. They can strengthen the ceiling segment and keep it stable without too much deformation. Besides, they can reduce the thickness of the ceiling segment, reducing the investment. There is a waterproof rubber sheet made of neoprene that connects the ceiling segment with the first shaft section through twelve Φ 24 M20 bolts. The structure of the first shaft section is much more complicated than the subsequent shaft section, viz., the standard shaft section. The first shaft section consists of twelve rectangular stiffened panels, twelve trapezoidal stiffened panels, one flange plate, and two different round pipes, as shown in Figure 4.

Figures 3 and 4 show the main parts of the vertical shaft sections. It can be seen directly from these two figures that not like the pipe jacking method or the shield method, the soil above the shaft sections during VTM construction is not removed but is compressed instead.

3. Numerical Model

The flowchart of the analysis process is shown in Figure 5.

3.1. Model Dimensions and Materials

To discharge sewage in Tieshangang District, Beihai, China, a sewage treatment plant project was carried out. The site is about 41 kilometers away from Beihai city. There were a total of 21 shaft sections for each shaft, the first shaft section was 50 cm in length, and the remaining sections were 60 cm in length. The shafts buried in the homogeneous poorly graded sand in this project are chosen to be analyzed in this paper. Readers interested in greater details about the Beihai VTM project and the geological conditions can read related information from Wang et al. [6].

The No. 6 shaft in the vertical tunneling project in Beihai is regarded as a benchmark to build the numerical model. The No. 6 shaft only passes through homogeneous poorly graded sand. The soil profile of the poorly graded sand is shown in Figure 6. In Figure 6, is the water content, is the liquid limit, is the plastic limit, γ is the unit weight of the soil, e is the void ratio, k is the hydraulic conductivity, Sr is the saturation, c is the cohesive force, and φ is the internal friction angle.

The outer diameter of the shaft is 0.5 m, and the length is 12.5 m. Owing to the stiffness of the shaft, which is much greater than that of the soil, the soil is set to be an axisymmetric deformable part, whose element type is CAX4R. In contrast, the pile is set to be an axisymmetric discrete rigid part, whose element type is RAX2. The original material parameters and the element type of the poorly graded sand are listed in Table 3. The constitutive model for the soil adopts the Mohr–Coulomb model. The cohesion and the angle of internal friction are 5.3 kPa and 8.7°, respectively.

Wang et al. [4] proposed that the vertical tunneling process can be regarded as a reverse piling process. In detail, the vertical tunneling process can be regarded as compressing the soil upwards with the inverted static pressure pile. According to the numerical simulation method of static pressure piles, the improved displacement penetration method is applied to simulate the construction process of the vertical tunneling process. Only half of the model was built for analysis because of the symmetry of typical modeling. In the real project, the top structure of the shaft is complicated (see Figures 3 and 4). To avoid the large deformation of the soil during the simulation, the ceiling segment on the top of the shaft is simplified into a smoother structure, like a pile shoe structure (Figure 7).

3.2. Mesh, Steps, and Boundary Conditions

In the model, the grid size is uniform along the depth direction. Due to larger soil deformation, a fine mesh was adopted in an area near the pile, while a coarse mesh was adopted in the soil away from the pile. The horizontal mesh size in the 2.5 m near the shaft is 0.3 m to 0.5 m. In comparison, the mesh size in the other region is 0.5 m to 1 m. The vertical mesh size is 0.15 m in all the regions.

There are three steps: the initial step, geostatic step, and penetration step. In the penetration step, the vertical displacement of the shaft reference point is 12.5 m, which is applied to simulate the penetration of the shaft.

The boundary conditions were as follows: the bottom boundary was fixed, and the top boundary was free. The soil movement in the normal direction to the right boundary is prohibited. The vertical soil movement in the left boundary during the penetration step is prohibited. Moreover, the vertical displacement of the reference point of the shaft is set to be 12.5 m in the penetration step, which is 0 in the initial step and geostatic step.

3.3. Modeling of the Shaft

The construction of the shaft can be approximated as reverse pile driving [4]. In the classical finite element calculations, the zipper-type technique is utilized to simulate the driving process of the circular pile in the soil [1417]. In the zipper-type technique, a rigid tube with a diameter of t = 1 mm is pre-embedded in the soil. Before the pile is penetrated, it is frictionless between the tube and the surrounding soil. During pile penetration, the pile slides along the rigid tube, separating the soil from the rigid pipe. Thus, the contact between the penetrating pile and the surrounding soil can be established.

This paper applies the zipper-type technique to simulate the vertical tunneling process. Before shaft driving, a small rigid tube with a radius of 1 mm is set to contact the soil. After the shaft is jacked upward, the shaft can be driven upward along the small rigid tube and separate the tube from the soil. The contact between the small rigid tube and soil and that between the shaft and soil are surface-to-surface contacts. The normal contact is a “hard” contact, and the tangential contact is frictionless.

Note that during the numerical simulation, the construction process of the horizontal tunnel is not considered. The response of the horizontal tunnel during the vertical tunneling process is not considered. Besides, groundwater is not considered.

4. Back Analysis of Mechanical Parameters of Soil

Although Table 3 shows the mechanical parameters of the poorly graded sand in the geotechnical experiment, the parameters still need to be adjusted according to the actual simulation situation. For example, in geotechnical experiments, the angle of internal friction of the poorly graded sand is 8.74°, which is quite different from the value range of 20°–40° for sand in general. The difference may be caused by factors such as artificial errors in the experimental process. Based on this, it is necessary to conduct a salience analysis of the soil mechanical parameters relative to the jacking force and determine the main factors and secondary factors that affect the jacking force, which can provide a theoretical basis for selecting parameters for subsequent numerical analysis.

4.1. Salience Analysis of Parameters Influencing the Jacking Force

Based on the ABAQUS software and the MATLAB program, according to the Latin hypercube sampling (LHS) method and the multivariate analysis of variance (multivariate ANOVA), the salience of each parameter relative to the jacking force is evaluated.

The LHS method is a special multidimensional stratified sampling method. Its main advantage is that the sampling value can fully reflect the overall distribution of random variables and ensure that all probability intervals are covered by sampling points [18]. The sampling process includes four steps.

First, the number of samples, namely Ns, is defined.

Second, each random variable is equally divided into Ns small regions that do not cross each other, and d groups of Ns small region sets are generated. d is the number of data groups composed of different random variables produced by the LHS method. The probability for each group of random variables to be in a small area is 1/Ns.

Third, different Latin hypercube configurations are obtained by randomly selecting small areas in each group without repetition.

Fourth, random samples are taken from the regions selected in the previous step, and these sampling points are from different Latin hypercube sampling results.

When two or more factors impact the dependent variable, the method of multivariate ANOVA can be used to determine whether each factor significantly impacts the dependent variable through hypothesis testing or not. Then, the best combination of dependent variables is found. The individual influence of the influencing factors on the dependent variable is called the main effect, and the common influence of the influencing factors on the dependent variable is called the interaction effect [19].

The basic steps of the method of multivariate ANOVA can be divided into the following three steps [20]:(1)The original hypothesis is proposed based on the dependent and multiple independent variables that need to be analyzed.The original hypothesis of multivariate ANOVA is that when the respective variables are at different levels, the mean of the dependent variable is not significantly different. The independent variable’s main effect and interaction effect are both zero; that is, the independent variable and their interaction effect have no significant influence on the dependent variable.(2)Decomposition of the independent variable variance is as follows:Two independent variables, A and B, are taken as an example, and the dependent variables are decomposed intoHere, SST is the total variance of the dependent variable, SSA and SSB are the variances caused by the independent effects of the independent variables A and B, respectively, SSAB is the variance caused by the interaction between independent variables A and B, and SSE is the variance caused by random factors. In this paper, the influence of random factors is not considered, that is, SSE = 0.In formula (1), SST is defined asHere, is the number of levels of the ith independent variable, o is the number of levels of the jth independent variable, xijg is the th sample’s value under the ith level of the independent variable A and the jth level of the independent variable B, nij is the number of samples under the ith level of the independent variable A and the j-h level of the independent variable B, and is the mean of the dependent variable.SSA is satisfied as follows:Here, is the mean of the dependent variable at the ith level of the independent variable A.SSB is satisfied as follows:Here, is the mean of the dependent variable at the jth level of the independent variable B.(3)The companion probability is calculated, the significance level α is given, and the final decision is made in this step.

The companion probability of the test statistic is compared with the given significance level α.

In the fixed-effect pattern, if the companion probability corresponding to the independent variable A is less than or equal to the given significance level α, null hypothesis will be rejected. That is, the observed population means are significantly different under different levels of the independent variable A. On the contrary, if the companion probability corresponding to the independent variable A is greater than the given significance level α, null hypothesis will be accepted. That is, there is no significant difference between the observed population means under different levels of the independent variable A, and different levels of the independent variable A have no significant effects on the dependent variable.

It should be noted that the influence of random factors is not considered, and only the main effect is considered for this paper. Moreover, the significance level α is set to 0.05.

Based on the LHS method and multivariate ANOVA, 20 numerical experiments are carried out to study the influence of multiple geotechnical parameters on the jacking force.

Regarding implementing the LHS method, Khaled [21] modified the lhsdesign function in MATLAB and obtained the lhsdesign-modified function. When the value range of the basic variable is known, the lhsdesign-modified function can be used to generate the Latin hypercube sampling array within the specified interval in MATLAB.

Under the condition that the value range of each independent variable is known, we first use the lhsdesign-modified function to obtain the Latin hypercube sampling array corresponding to each variable, i.e., each geotechnical parameter. Then, only the main effect is considered, and 20 sets of the geomechanical parameters obtained by the LHS method are input into Abaqus, different models are established, and the maximum jacking force of each model is calculated. After that, the Latin hypercube sampling array and the maximum jacking force corresponding to each array are set as the input data. The multivariate ANOVA is carried out on the influencing factors of the jacking force by the means of MATLAB.

4.2. Multivariate Significance Analysis of the Jacking Force

Based on the geotechnical parameters of the vertical tunneling project in Beihai, the value ranges of the geotechnical parameters of the project are shown in Table 4.

In detail, the elastic modulus of soil E is 0.5 times∼5 times the compressive modulus of the poorly graded sand layer, i.e., Es = 5.42 MPa. According to the suggested values of the static earth pressure coefficients and Poisson’s ratios for various kinds of soils (see Table 5), the value range of static earth pressure coefficients is 0.18∼0.33 and the range of Poisson’s ratio equals that of the gravel soil and sandy soil, that is, 0.15∼0.25.

Direct shear testing is an easy and cost-effective way to figure out soil properties [2326]. Based on the results of in situ tests, as well as indoor geotechnical tests including the direct shear tests, Zhong [27] evaluated the geotechnical engineering conditions of the sedimentary strata in Beihai and determined the design parameters of each geotechnical stratum (see Table 6). According to the geotechnical design parameters of Zhong [27] and the indoor geotechnical test results in the vertical tunneling project in Beihai in this paper, the cohesion range of the poorly graded sand is 0 kPa∼10 kPa and the angle of internal friction range is 5°∼35°.

4.3. Inversion Results of the Model Parameters

Twenty groups of LHS parameters corresponding to six geotechnical parameters, including elastic modulus, Poisson’s ratio, soil gravity, cohesion, angle of internal friction, and static earth pressure coefficient, are obtained and shown in Table 7. According to the parameters in Table 7, twenty numerical simulation experiments were performed to obtain the maximum jacking force values corresponding to the twenty groups of parameters (see Table 7). It should be noted that the water pressure is not considered for this paper. According to Figure 8, the calculated jacking force was compared with the measured jacking force of the No. 6 shaft in Section 2. The results show that among the twenty groups of parameters shown in Table 7, the calculated jacking forces of No. 12 and No. 17 coincide with the measured ones.

Then, multivariate ANOVA was performed on the seven parameters, as shown in Table 7, and the companion probabilities of the six influencing parameters relative to the jacking force were calculated (see Table 7). The criterion for determining the impact significance is as follows: if the companion probability of a parameter is less than the significance level α = 0.05, the parameter has a significant impact on the jacking force. Furthermore, the smaller the companion probability, the more significant the impact. If the companion probability of a parameter is greater than the significance level α = 0.05, the effect of the parameter on the jacking force is not significant.

The companion probability corresponding to the elastic modulus of the soil is , indicating that the elastic modulus has a very significant effect on the jacking force. For the angle of internal friction, the companion probability is , indicating that the angle of internal friction significantly influences the jacking force. The values corresponding to the other factors are all greater than 0.05. That is, the influence of the other factors on the jacking force is insignificant. The companion probabilities corresponding to the following six parameters are in descending order: the elastic modulus of soil, angle of internal friction, static earth pressure coefficient, unit weight of soil, Poisson’s ratio, and cohesion. For the jacking force of the shaft in the vertical tunneling process, the elastic modulus and the angle of internal friction of soil are the main influencing factors. The secondary influencing factors are the static earth pressure coefficient, unit weight of soil, Poisson’s ratio, and cohesion.

It can be seen from Figure 8 that the calculated jacking forces corresponding to No. 12 and No. 17 are the closest to the measured jacking force. According to the geomechanical parameters of No. 12 and No. 17, as shown in Table 7, the values of the secondary influencing factors are determined directly and the inversion analysis of the two main influencing factors (E and φ) is carried out. The optimum parameters can be inverted by ignoring the influence of the secondary influencing factors and only considering the influence of the main influencing factors. Finally, the optimum inversion geomechanical parameters of the poorly graded sand can be obtained as follows: cohesion c = 5.29 kPa, Poisson’s ratio μ = 0.24, soil gravity γ = 19.5 kN/m3, static earth pressure coefficient K0 = 0.32, angle of internal friction φ = 33°, and elastic modulus of soil E = 25 MPa.

5. Results and Discussion

5.1. Validation via Comparison with Field Measurement Data

The comparison between the calculated jacking force of the numerical model and the measured jacking force is shown in Figure 9. As shown in Figure 9, the abscissa is the jacking distance of the shaft and the ordinate is the jacking force. The changing trend of the numerical model and the measured jacking force is relatively consistent, and the maximum jacking force is almost the same.

It can be seen from Figure 9 that the measured jacking force first increases and then decreases with the increasing jacking distance. This is related to the compressing and shearing process of the soil above the shaft. The soil above the shaft would be compressed first. After that, an ultimate state may be achieved and a shear band may appear; then, the soil would slip along the shear band [6]. Due to this, the jacking force would increase first before the soil reaches the ultimate state and then decrease because of the soil’s slip.

There are some differences between the measured and calculated jacking force due to the ceiling segment being simplified in the numerical model. The ceiling segment on the top of the shaft (see Figure 3) is lifted upward together with the shaft in the whole process, and it will be removed by the diver after shaft construction is completed. In the numerical model in this paper, the ceiling segment is simplified to a pile shoe-like shape (see Figure 7). As a result, it is easier for the soil above the shaft to reach its ultimate state in the numerical model than in the field project. Thus, the maximum jacking force could be reached earlier in the numerical model than the measured one.

5.2. Effect of Elastic Modulus

Figure 10 shows the jacking force with respect to the jacking distance when the elastic modulus of soil E ranges from Es to 5Es (Es = 5.42 MPa). It can be seen from Figure 10 that the elastic modulus has a great influence on the jacking force. In the range 0 m∼10.5 m of the jacking distance, i.e., 0∼0.84 L0, where L0 is the total length of the shaft, the increase in the elastic modulus causes an increase in the jacking force, while in the range 10.5 m∼12.5 m, i.e., 0.84 L0L0, the increase in the elastic modulus of soil has little effect on the jacking force.

5.3. Effect of the Angle of Internal Friction

Figure 11 shows the effect of the angle of internal friction on the jacking force. It can be seen from Figure 11 that the angle of internal friction greatly influences the simulation results. When the jacking distance is 0∼11 m, that is, 0∼0.88 L0, the jacking force increases with an increase in the angle of internal friction; when the jacking distance is 11∼12.5 m, that is 0.88 L0L0, with an increase in the angle of internal friction, the change in the jacking force is slight.

5.4. Effect of the Unit Weight of Soil

Figure 12 shows the soil unit weight’s effect on the jacking force. It can be seen from Figure 12 that when the value is 15 kN/m3∼25 kN/m3, the unit weight of soil has little effect on the jacking force.

5.5. Effect of Cohesion

Figure 13 shows the effect of cohesion on the jacking force. It can be seen from Figure 13 that when the cohesion value ranges from 2 kPa to 10 kPa, cohesion has little effect on the jacking force.

5.6. Effect of Poisson’s Ratio and the Static Earth Pressure Coefficient

Figure 14 shows the effect of Poisson’s ratio and the static earth pressure coefficient on the jacking force. As shown in Figure 14, there is little effect of Poisson’s ratio and the static earth pressure coefficient on the jacking force when Poisson’s ratio ranges from 0.15 to 0.25 and the static earth pressure coefficient ranges from 0.175 to 0.25.

According to Figures 913, the elastic modulus and the angle of internal friction of the soil have the greatest influence on the jacking force, followed by Poisson’s ratio and the static earth pressure coefficient, and finally, the unit weight and cohesion. It is consistent with the conclusion of the model parameter inversion analysis as described in Section 3.3.

6. Conclusions

(1)The Latin hypercube sampling method and the multivariate ANOVA method were applied to study the influence of the soil geomechanical parameters on the jacking force. The result shows that the elastic modulus and the angle of internal friction of soil were the main influencing factors and that Poisson’s ratio, static Earth pressure coefficient, unit weight, and cohesion were the secondary factors. The optimum inversion geomechanical parameters of the poorly graded sand can be obtained as follows: cohesion c = 5.29 kPa, Poisson’s ratio μ = 0.24, soil gravity γ = 19.5 kN/m3, static earth pressure coefficient K0 = 0.32, angle of internal friction φ = 33°, and elastic modulus of soil E = 25 MPa.(2)Based on the numerical model, the effects of soil elastic modulus, angle of internal friction, unit weight, cohesion, and other factors on the jacking force were analyzed. The research results show that in the range of the jacking distance from 0 to 0.84 L0, the increase in the elastic modulus causes an increase in the jacking force, while in the range 0.84 L0L0, the increase in the elastic modulus of soil has little effect on the jacking force. When the jacking distance is 0∼0.88 L0, the jacking force increases with an increase in the angle of internal friction; when the jacking distance is 0.88 L0,∼L0, with an increase in the angle of internal friction, the change in the jacking force is slight. The unit weight of soil, cohesion, Poisson’s ratio, and static earth pressure coefficient have little effect on the jacking force.

Data Availability

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

Conflicts of Interest

The authors declare no conflicts of interest.

Acknowledgments

This work was financially supported by the Natural Science Foundation of Zhejiang Province (No. LHZ23E080001).