Numerical Computing Research on Tunnel Structure Cracking Risk under the Influence of Multiple Factors in Urban Deep Aquifer Zones

.

zones, nearly 90% of subway stations have leakage problems, and nearly 30% of the leakage occurs at deformation joints [3,4]. Once water seepage occurs, it will not only reduce the strength of the concrete structure but also bring safety hazards to later operations. To ensure the safety of tunneling operations, a new numerical analysis is required, by means of which the influence of multiple factors in urban deep aquifer zones and deformation joints on structure cracking risk can be reasonably researched.
At present, various methods, namely empirical methods, analytical methods, numerical methods, and model test methods have made great accomplishments in the study of tunnel cracking risk. Liu et al. [5] analyzed 109 tunnels with lining cracks and found that the main influencing factors for tunnel lining cracking are bias pressure, concrete shrinkage, back cavities, landslides, and relaxation pressure. Li et al. [6] found that temperature and concrete shrinkage were important reasons for concrete cracking in multi-purpose utility tunnels. These studies obtained the influencing factors of cracking based on empirical statistics, but cannot quantitatively analyze these factors. Based on double K theory [7][8][9], damage mechanics theory [10], and linear elastic fracture mechanics [11], some researchers have studied the propagation mechanism and theoretical model of concrete cracks from the concrete material itself. The model test method can truly reflect the cracking situation of the structure, but the testing process is relatively complex [12,13]. Liu et al. [14,15] carried out a series of model tests on structure cracking. Numerical simulation can consider tunnel cross-section shapes, geological properties, complex construction processes, and the coupling interaction between surrounding rock and tunnel structures [16][17][18]. It is currently a commonly used method for studying structure cracking risk. Based on the load-structure method, the effect of strata on the structure is simplified as the load acting on the structure [19,20]. Huang et al. [21] established a numerical model to analyze the influence mechanism of geological factors and temperature factors, and determined that concrete cracking depends on the maximum tensile stress of the structure. He et al. [22], Liu et al. [23], and Li et al. [24] also investigated the impact of different influencing factors on structure cracking based on the load-structure method. For the abovementioned studies, nevertheless, existing studies have generally focused on the effects of a single influencing factor, such as geological conditions [25], temperature [26,27], etc. on structure cracking risk, and lack of comprehensive analysis under the influence of multiple factors. At the same time, the load-structure method analyzes the concrete structure as a plate [28,29], ignoring the interaction between the concrete structure and the surrounding soil, which is not consistent with the actual engineering situation and results in inaccurate calculation results.
Based on the above shortcomings, a theoretical calculation model and a three-dimensional concrete-soil interaction thermo-mechanical coupling model was established and the maximum tensile stress of the tunnel structure was analyzed in this paper. Furthermore, the response mechanism of structural stress and deformation under the influence of the grade of rock and soil mass, overburden thickness, temperature difference, structure's length-height ratio, structure's thickness, and structure's elastic modulus was investigated, and the stress and deformation response characteristics of the structure with deformation joints were explored. The research results can provide important technical guidance for the cracking risk of tunnel structure under the influence of multiple factors in urban deep aquifer zones and for reducing the risk and cost of engineering construction.

Theoretical Model
The deformation of a concrete structure is constrained by the foundation. Based on this, a concrete wall model on an elastic foundation is established, as shown in Figure 1. It is assumed that the concrete wall and the foundation are deformable elastic structures, and the constraints between them can be reflected by the relationship between stress and deformation [30,31]. Assume that the horizontal shear stress is linear with the horizontal displacement: where τ is the shear stress on the interface between the concrete wall and foundation, u is the horizontal displacement of the foundation, and C x expresses the horizontal resistance coefficient of the foundation. Taking a unit at x point of the structure for stress analysis, the balance equation can be expressed as where N is the internal force of the unit, Q represents the shear stress of the foundation to the unit, σ x is the internal stress of the unit, and H and t are the height and width of the unit, respectively. The displacement at any point of the structure is composed of constrained displacement and free displacement: where u σ is the constrained displacement, is the linear expansion coefficient of concrete, and T is the temperature difference.
According to the generalized Hooke's law [32], the stress-strain relationship of the structure can be expressed as: where E is the elastic modulus of concrete. Substituting Equations (1) and (7) The general solution of Equation (9) is =+ cosh sinh uA βxB βx (10) where A and B are obtained by boundary conditions. At the 0 point, u = 0. At the L/2 point, σ x = 0. L is the length of the concrete wall. Then Equation (10) can be given as Substituting Equation (11) into Equation (10) obtains Then the horizontal stress of the concrete wall can be obtained: EαT L β (13) According to Equation (13) and the symmetry of the structure, the maximum horizontal stress of the concrete wall is

Engineering Background
Hangzhou Road multi-purpose utility tunnel is located in Hangzhou Road, Hightech Industrial Development Zone, Liaocheng City ( Figure 2). The utility tunnel starts from Guangyue Road in the west and ends at Lushan Road in the east. The Hangzhou Road multi-purpose utility tunnel has a total length of 3162.7 m, with a single-cabin rectangular section and cast-in-site reinforced concrete structure. About 30 m of the tunnel is divided into one construction section, and deformation joints are set in two adjacent sections. The main structure concrete of the utility tunnel adopts self-waterproof C35 concrete. The Quaternary geomorphic units along the proposed tunnel belong to the central part of the Yellow River alluvial plain, with relatively flat terrain and low terrain. The strata are mainly composed of silty clay, silt, and sand.
Liaocheng's river systems mainly belong to the Yellow River and the Haihe River, each with numerous tributaries. The types of groundwater in Liaocheng City are relatively complex. According to the occurrence conditions of groundwater, the water-bearing rock groups are divided into rock groups bearing loose rock pore water, and rock groups bearing carbonate-rock-like fissure karst water. The loose rock pore water aquifers can be divided into shallow unconfined aquifers, intermediate confined water aquifers, and deep confined water aquifers.

Numerical Model
To further study the stress and displacement characteristics of the tunnel, the simulation software of FLAC3D is applied in this study for further numerical work. FLAC3D is a commercial finite difference numerical simulation software that includes 18 elastic and plastic constitutive models, five calculation modes, and various modes that can be coupled with each other. Compared with other software, FLAC3D has the advantages of fast solving speed and accurate calculation results and can simulate complex geotechnical problems, and has been widely used in the field of geotechnical engineering [33][34][35][36]. Wang [37] et al. carried out a FLAC3D numerical simulation study of the deformation of a coal mine's deep buried track, and the calculation results are consistent with the field measurement results, with a small error (within 5%). Unlu and Gercek [38] employed FLAC3D software to calculate the radial boundary displacement along the longitudinal direction of a circular tunnel. Through the built-in Fish programming language, batch modeling was carried out and structured hexahedral-only volume meshes were created in this study.

Model Validation
Model validation is an important step to ensure appropriate and acceptable results of computational simulation, which is usually carried out by comparing model simulation data with independent experimental or field datasets [39][40][41]. In Huang [42], the vertical displacement of a utility tunnel project in Chongqing was measured. Based on this reference, the same numerical model was established to simulate the tunnel displacements. The length of the tunnel is 8.9 m, the width is 3.9 m, and the buried depth is 4.2 m. The thickness of the tunnel is 0.4 m, the total length is 30 m, the entire model size is 64 m × 30 m × 20 m. Eight-node brick elements are used to simulate the soil and concrete. Figure 3 shows the simulation and measured results of the maximum vertical displacement of the top and bottom of the structure. The error of the maximum vertical displacement between the numerical result and the actual measured value is small, which is less than 10%. Therefore, the numerical model is reliable.

Three-Dimensional Numerical Model
The schematic diagram of the numerical model is shown in Figure 4. Based on the project design drawings, the cross-section of the concrete structure is a rectangle, L = 10 m, h = 4 m, and w = 30 m. The thickness of the structure is D = 0.5 m, buried depth is 5 m. According to Saint Venant's principle [43], the whole model size is 70 m × 21 m × 30 m, which is composed of 32,860 grid points and 29,640 zones. Eight-node brick elements are used to simulate the soil and concrete. All tunnel meshes established in the study adopt the same size. The mesh is divided into 1 m in the longitudinal direction of the tunnel, 0.5 m in the length and height directions, and 0.2 m in the thickness direction. The deformation joints are mainly connected by waterstops, and the liner structural element is used to simulate the deformation joints. The deformation joints are set in the middle of the structure, as shown in Figure 5, with a width of 30 mm and an elastic modulus of 2.6 MPa [44,45]. The mechanical behavior of concrete is described through the elastic constitutive model, and the rock and soil mass's behavior is simulated with the Mohr-Coulomb constitutive model. The structure adopts C35 concrete, and the surrounding soil adopts Grade IV rock and soil mass for calculation. Related physical and mechanical parameters are listed in Table 1 [46,47]. The surface between the structure and the soil is simulated by the interface element, and the friction angle is 15° [48]. Assume that the temperature of the structure decreases gradually from outside to inside, and the temperature difference is taken as 5 °C [49].  For the whole model, the soil is homogeneous, continuous, and isotropic, regardless of groundwater. The top surface of the model is set as a free boundary, the bottom of the model is a fixed boundary, and normal constraint is applied to the rest of the surfaces. The isotropic conduction model is used in the thermo-mechanical coupling calculation, the internal and external temperatures of the structure are fixed, and adiabatic boundaries are prescribed around the model surfaces.

Numerical Simulation Result
After the casting and backfilling of the open-cut cast-in situ concrete structure, the longitudinal stress and deformation of the structure are shown in Figure 6.
As seen in Figure 6a, the maximum longitudinal stress is concentrated at the top and bottom of the structure, and the top and bottom are in tension on the inside and compression on the outside. When the structure has no deformation joints, the maximum tensile stress of the upper surface is 2.37 MPa, which is greater than that of the lower surface of 1.87 MPa. With deformation joints, the tensile stress is distributed in a V-shaped curve along the Y direction of the structure, the minimum tensile stress appears at the deformation joint, and the maximum tensile stress appears at both ends of the model. This kind of distribution characteristic of tensile stress was also reported by Zhuo et al. [20]. The maximum tensile stress of the upper surface is 2.07 MPa, which is greater than that of the lower surface of 1.63 MPa. Compared with no deformation joints, the maximum tensile stress of the upper and lower surface is reduced by 12.7% and 12.8%, respectively, and the tensile stress at the deformation joints of the upper and lower surface is reduced by 58.6% and 64.7%, respectively. Taking the horizontal resistance coefficient C x as 1.0 N/mm 3 [50], calculating the maximum tensile stress of the lower surface using Equation (14) gives 1.50 MPa, which is slightly smaller than the numerical simulation result, with an error of 8.0%, and the theoretical calculation value is relatively safe.
As shown in Figure 6b, the maximum vertical displacement is concentrated at the top and bottom of the structure; the top of the structure subsides, and the bottom of the structure uplifts. When the structure has no deformation joints, the maximum vertical displacement of the upper surface is 9.75 mm, which is greater than that of the lower surface of 5.94 mm. With deformation joints, the vertical displacement changes greatly at the deformation joints, and the maximum vertical displacement of the upper surface is 10.44 mm, which is greater than that of the lower surface of 6.38 mm. Compared with no deformation joints, the maximum vertical displacement of the upper and lower surface increases by 7.1% and 7.4%, respectively.
Since the structural strength at the deformation joints is smaller than that of the concrete on both sides, the stress is released at the deformation joints, and the tensile stress of the structure is significantly reduced, while the deformation at the deformation joints increases. The increase in the deformation at the deformation joints will further cause stress release in the structure, leading to a decrease in stress at the deformation joints. The deformation joints can effectively reduce the longitudinal tensile stress of the structure. This conclusion is basically in agreement with previous research [21][22][23][24]. However, they did not consider the influence of multiple factors on the stress and deformation of the structure.

Rock and Soil Properties Effect
Based on the "Specifications for design of highway tunnels" (JTG 3370-2018) in China, the rock and soil mass can be divided into four grades (III, IV, V, VI) according to the strength from strong to weak [47]. The physical and mechanical parameters of different grades of rock and soil mass are shown in Table 2.  Figure 7 shows the longitudinal tensile stresses of the structure under different rock and soil mass grades. The density of the rock and soil mass directly affects the load of the overlying soil on the structure. At the same time, the better the strength of the rock and soil, the stronger the self-stabilization ability of the rock and soil, and the smaller the additional stress of the overlying soil on the structure. Therefore, as the rock and soil mass grade increases, the maximum longitudinal tensile stress of the structure gradually increases, and the risk of structural cracking increases. The properties of the rock and soil mass around the structure are directly related to the structural stress, which is consistent with the results of Weng et al. [12] and Huang et al. [21]. Without deformation joints, the maximum longitudinal tensile stress under grade Ⅵ rock and soil mass is 3.00 MPa, which is 36.4% higher than that under grade Ⅲ rock and soil mass of 2.20 MPa. With deformation joints, the maximum longitudinal tensile stress under grade Ⅲ rock and soil mass is reduced to 1.92 MPa, which is a reduction of 12.7%. The maximum longitudinal tensile stress under grade Ⅵ rock and soil mass is reduced to 2.63 MPa, which is a reduction of 12.3%. The settlement displacements under different rock and soil mass grades are shown in Figure 8. The same as the influence law of the maximum tensile stress of the structure, the rock and soil mass grades increase, and the maximum settlement displacement of the structure increases. Therefore, by improving the strength properties of the soil around the structure, the settlement displacements can be effectively reduced, and this result is consistent with that of Wang et al. [50] and Tan et al. [51]. Without deformation joints, the maximum settlement value under grade Ⅵ rock and soil mass is −13.50 mm, which is 89.3% higher than that under grade Ⅲ rock and soil mass of −7.13 mm. The deformation growth is relatively obvious. With deformation joints, the maximum settlement value increases by 10.1% under grade Ⅲ rock and soil mass. The growth magnitude of the maximum settlement value gradually decreases while the rock and soil mass grade increases.

Overburden Thickness Effect
In previous numerical studies [19][20][21][22][23][24], there have been few studies on the impact of overburden thicknesses on the structural cracking risk. Figure 9 shows the longitudinal tensile stresses of the structure under different overburden thicknesses when the grade of rock and soil mass is Ⅳ. The overburden thickness directly affects the load of the upper soil on the structure, so as the overburden thickness increases, the longitudinal tensile stress of the structure increases, and the risk of structural cracking increases. Without deformation joints, the maximum longitudinal tensile stress increases by 13.1% as the overburden thickness increases from 5 m to 20 m, which indicates that overburden thicknesses have a relatively small influence on structural cracking. With deformation joints, the maximum longitudinal tensile stress is reduced by 12.7% under an overburden thickness of 5 m. The decreased magnitude of longitudinal tensile stress gradually decreases as the overburden thickness increases.  Figure 10 shows the settlement displacements of the structure under different overburden thicknesses. Without deformation joints, the maximum settlement displacement of the structure increases by 60.6% as the overburden thickness increases from 5 m to 20 m. With deformation joints, the maximum settlement displacement increases. As the overburden thickness increases, the growing magnitude of the maximum settlement displacement decreases.  Figure 11 shows the longitudinal tensile stresses of the structure under different temperature differences when the grade of rock and soil mass is Ⅳ and the overburden thickness is 5 m. The longitudinal tensile stress of the structure increases with the increase in temperature difference. Without deformation joints, the maximum longitudinal tensile stress is 0.82 MPa under a temperature difference of 0 °C, and the maximum longitudinal tensile stress is 3.14 MPa under a temperature difference of 7.5 °C. The maximum longitudinal tensile stress of the structure increases by 2.8 times. The temperature difference has a significant impact on the tensile stress of the structure, so the change in temperature difference is the main factor leading to structural cracking. The observed influence law of temperature differences in this study is similar to that reported by Liu et al. [5] and Li et al. [6], but the results obtained through numerical calculations in this article are more intuitive and obvious. Based on the works of Huang et al. [21], Wang [30], and De Schutter [31], the phenomenon above can be attributed to the expansion deformation of concrete structures and the increase in temperature stress. Meanwhile, it can be seen from Equation (14) that the maximum tensile stress of the structure is linearly related to the temperature difference, and the numerical simulation results are consistent with the theoretical calculation results.

Temperature Difference Effect
With deformation joints, the maximum longitudinal tensile stress of the structure increases by 2.35 times as the temperature difference increases from 0 °C to 7.5 °C. The deformation joints can effectively reduce the influence of temperature differences on the structural tensile stress, which is consistent with the results of Dong et al. [1]. Compared with no deformation joints, the maximum longitudinal tensile stress is reduced by 2.4% under a temperature difference of 0 °C, and 14.6% under a temperature difference of 7.5 °C. The greater the temperature difference, the greater the reduction magnitude of the longitudinal tensile stress, and the more significant the effect of deformation joints.  Figure 12 shows the settlement displacements of the structure under different temperature differences. Without deformation joints, as the temperature difference increases from 0 °C to 7.5 °C, the maximum settlement displacement of the structure increases by 1.7%, and temperature difference has little effect on structural settlement displacement, and this result is consistent with that of Huang et al. [21]. With deformation joints, the maximum settlement value of the structure increases by 10.4% as the temperature difference increases from 0 °C to 7.5 °C. Compared with no deformation joints, the larger the temperature difference, the greater the growth magnitude of structural deformation.  Figure 13 shows the longitudinal tensile stresses of the structure under different length-height ratios when the height of the structure is 4 m. With the increase in the length-height ratio, the longitudinal tensile stress of the structure increases, and the risk of structural cracking increases, which is consistent with the influence law presented in He et al. [22]. However, in his research, there was no quantitative analysis of the influence degree of structural length-height ratios. Without deformation joints, the maximum longitudinal tensile stress of the structure increases by 43.9% as the length-height ratio increases from 1.5 to 3. With deformation joints, the maximum longitudinal tensile stress is reduced by 16.0% under a length-height ratio of 3. As the length-height ratio increases, the decreased magnitude of longitudinal tensile stress increases gradually. In Figure 14, with the increase in the length-height ratio of the structure, the settlement displacement increases from −0.80 mm to −19.97 mm, and the increment is very significant. With deformation joints, the maximum settlement displacement increases. The larger the length-height ratio, the smaller the deformation growth rate.  Figure 15 shows the longitudinal tensile stresses under different structure thicknesses when the length-height ratio of the structure is 2.5. The greater the structure thickness, the smaller the longitudinal tensile stress of the structure, which is consistent with the influence law presented by Equation (14) and Liu et al. [23]. Without deformation joints, the maximum longitudinal tensile stress of the structure decreases by 14.9% as the structure thickness increases from 0.3 m to 0.6 m. With deformation joints, the maximum longitudinal tensile stress is reduced by 14.4% under a structure thickness of 0.6 m. The greater the thickness of the structure, the greater the reduction magnitude of the longitudinal tensile stress. In Figure 16, when without deformation joints, the ability of the structure to resist deformation increases as the structure thickness increases [22,23]. The maximum settlement displacement of the structure decreases by 66.0% as the structure thickness increases from 0.3 m to 0.6 m, with a significant change. With deformation joints, the maximum settlement displacement increases by 9.9% under a structure thickness of 0.6 m. The change magnitude of the maximum settlement displacement increases with the increase in structure thickness.

Structure Material Effect
The main difference between different concrete materials is the elastic modulus, and the elastic moduli of different concretes are shown in Table 3.  Figure 17 shows the the longitudinal tensile stresses under different concretes when the length-height ratio of the structure is 2.5 and the structure thickness is 0.5 m. With the increase in the elastic modulus, the longitudinal tensile stress of the structure increases, mainly because the elastic modulus is positively correlated with the temperature stress of concrete [26,27]. Without deformation joints, the maximum longitudinal tensile stress of the structure increases by 7.1% as the elastic modulus increases from 28 GPa to 32.5 GPa. With deformation joints, the maximum longitudinal tensile stress is reduced by 12.8% under an elastic modulus of 28 GPa, and 12.4% under an elastic modulus of 32.5 GPa. In Figure 18, when without deformation joints, the elastic modulus of the structure increases, and the ability of the structure to resist deformation increases. The maximum settlement displacement of the structure decreases by 5.1% as the elastic modulus increases from 28 GPa to 32.5 GPa. With deformation joints, the maximum settlement displacement increases by 7.2% under an elastic modulus of 32.5 GPa. The growth magnitude of the maximum settlement displacement increases with the increase in the elastic modulus.

Conclusions
(1) After the casting and backfilling of the open-cut cast-in situ concrete structure, the maximum longitudinal stress and vertical displacement are concentrated at the top and bottom of the structure. The maximum longitudinal tensile stress and vertical displacement of the upper surface are greater than those of the lower surface.
(2) The maximum longitudinal tensile stress, settlement displacement, and the cracking risk of the structure increase with the increase in the grade of rock and soil mass, overburden thickness, temperature difference, and the structure's length-height ratio, while decreasing with the increase in the structure's thickness. With the increase in the concrete's elastic modulus, the maximum longitudinal tensile stress increases, and the maximum settlement displacement decreases. The effect of temperature difference on longitudinal tensile stress is the most significant, with the maximum tensile stress of the structure increasing by 2.8 times.
(3) When with deformation joints, the tensile stress is distributed in a V-shaped curve along the Y-direction of the structure, with the minimum tensile stress occurring at the deformation joint and the maximum tensile stress occurring at both ends of the model. The deformation joints can effectively reduce the longitudinal tensile stress of the structure, and the reduction magnitude of the tensile stress is the largest at the deformation joints, which is 64.7%. The vertical displacement of the structure changes abruptly at the deformation joints, and the maximum vertical displacements at the top and bottom of the structure increase to a certain extent compared with a structure with no deformation joints.