Numerical Simulation Study on Frost Heave during the Freezing Phase of Shallow-Buried and Undercut Tunnel Using the Freeze-Sealing Pipe Roof Method

Abstract

This article begins with the engineering geological conditions and freezing design scheme of the Gongbei Tunnel’s underground excavation section, then applies the mathematical model theory of the horizontal freezing-tunnel-formation freezing temperature field and frost heave displacement field, and builds a coupled two-dimensional finite element calculation model. The development law of the frozen soil curtain and the variation law of frost heave displacement during the freezing phase were studied by comparing on-site observation data. According to the findings of this study, the construction of the artificial frozen curtain is mostly based on two types of freezing tubes that freeze the soil between jacked pipes and seal the water. At 90 days, the thickness of the frozen soil curtain ranges from 2.32 m to 2.58 m, guaranteeing that its strength fulfills water-sealing safety criteria. The distribution and variation of frost heave displacement are highly related to engineering geological circumstances, the freezing scheme, the frozen soil curtain development process, and the pipe curtain structure. The maximum vertical frost heave displacement value at any time is located at the centerline, which is 155.67 mm at 90 d. The numerical simulation findings are acceptable and can potentially be utilized for predicting frost heave in subsequent projects. More research is required to effectively represent complicated working conditions and to develop more exact large-scale numerical models for tunnel excavation, support structure building, and other situations.

1. Introduction

With the steady rise of China’s economic aggregate over the last decade, the urbanization rate of permanent residents has exceeded 65% [1,2]. The continuous growth of urban populations and the increasing density of construction facilities create many pressures such as land resource scarcity, crowded living space, traffic congestion, and environmental deterioration, which must be alleviated by the utilization and development of urban underground space as well as the acceleration of urban underground transportation construction [3].

Construction of urban underground tunnels, as one of the important solutions, currently faces major technical challenges, such as the safety risks of underpass operating tunnels or existing buildings, tunnel construction in soft soil layers with high water content, and predicting and regulating settlement deformation [4,5]. Tunnel excavation frequently results in the deterioration of rock and soil mass and surface deformation of the project and its surrounding environment, resulting in underground pipeline rupture, serious ground settlement, or even collapse, cracking, and damage to adjacent structures [6,7].

As an underground excavation construction technology, the pipe roof method (PRM) is commonly utilized in the building of large cross-section tunnels in soft soil locations [8]. Since its first use in the building of subterranean tunnels in Hong Kong, China, in 1984, this approach has been used in a variety of cities, including Taipei [9], Shanghai [10], Beijing [11], Shenyang [12], Zhengzhou [13], Xiamen [14], and Taiyuan [15]. In comparison to other construction methods, the PRM can effectively control the stratum deformation caused by tunnel excavation by forming a temporary support structure with high rigidity, reduce the damage caused by surface settlement to adjacent buildings, and be applied flexibly in response to complex and changing engineering geological conditions [16].

However, when this method is used for soft soil layers with high compressibility, high water content, and low strength in China’s coastal areas such as Tianjin, Shanghai, and Guangdong, there are relatively large construction risks, which primarily include: uncontrolled jacking direction of steel pipe, sharp increase in jacking force, excessive surface settlement, pipe curtain damage, and poor water sealing performance [17,18]. Surface settling and water sealing performance, in particular, frequently determine project success or failure.

Nowadays, Chinese engineering experts have suggested a novel combination of artificial ground freezing (AGF) technology with PRM, resulting in the formation of a new construction technique known as the freeze-sealing pipe roof method (FSPR). The freezing tube is cleverly installed inside the jacked pipe in FSPR, and the soil between the adjacent jacked pipes is artificially frozen, which can not only replace the traditional PRM’s water stop lock catch but also form a high-strength composite supporting structure of “steel pipe curtain–frozen soil curtain”.

The practicality of FSPR was essentially determined by Hu et al. [19], who outlined the fundamental idea of “freezing, anti-weakening, and limiting frost heave” and performed a preliminary study on several potential freezing systems using the finite element program ANSYS. Following the scale model test, three unique structures—the “main circular freezing-tube”, the “strengthening profiled freezing-tube”, and the “hot brine limiting-tube”—were created for the first time. Additionally, the impact of freezing the water between the pipes was proven [20]. Based on the results of large-scale model testing, Ren et al. [21] concentrated on analyzing the strengthening freezing impact of profiled freezing tubes and proposed opening them in advance during construction. Kang et al. [22] studied the changes in the frozen soil curtain throughout the freezing, tunnel excavation, and thawing phases using a thermal coupling analysis approach. They emphasized that the temperature field is an important aspect in forecasting the thickness of the frozen soil curtain. Hu et al. [23] investigated the bearing capacity of “steel pipe-frozen soil” composite constructions at four different temperatures and presented limit-state criteria for water sealing. Zhang et al. [24] analyzed the measured data of surface displacement, tunnel arch displacement, and horizontal convergence during the construction of the Gongbei Tunnel and proposed that when using the step method to excavate the upper part of a large cross-section tunnel, attention should be paid to the arch displacement caused by excavation construction and the frequency of on-site real-time monitoring be increased. The author’s team also conducted theoretical and experimental studies on the freezing temperature field of the FSPR method [25,26,27], the improved design of special-shaped freezing pipes, and the mechanical properties of jacked pipe-frozen soil composite structural mechanics, and made corresponding guiding suggestions for the specific construction process.

According to the research findings, the technological obstacles encountered by the FSPR approach in its use are mostly exhibited in three areas: First and foremost, it has a completely novel freezing idea. Second, the influence of frost heave on the strata, surface, and surrounding buildings (structures) must be considered. Third, there are few examples of this approach in use both at home and abroad, and there are still certain practical issues that need to be addressed and improved in building. As a result, extensive research on its freezing temperature field and frost heave displacement field is required.

After analyzing the engineering geological conditions and freezing design scheme of the Gongbei Tunnel’s underground section, a thermodynamic coupling two-dimensional finite element calculation model based on the mathematical model theory of the stratum freezing temperature field and frost heave displacement field of the horizontal freezing tunnel was built in this article. The construction law of frozen soil curtain and the modification of frost-heave displacement during the freezing period were researched using the calculation platform of finite element software COMSOL Multiphysics 5.6 and integrated with the project’s field-observed data.

2. Project Background

The Gongbei Tunnel stretches from the waters of Gongbei Bay to the Maoshengwei Management Area of the Guangdong Provincial Public Security Border Protection Fifth Brigade in Xiangzhou District, Zhuhai City, near Macau. As illustrated in Figure 1, the concealed excavation portion crossing through the Gongbei Port is a two-way, six-lane overlapping tunnel built with the FSPR method. It is a typical shallow-buried and concealed excavation tunnel with a burial depth of 4~5 m, a length of 255 m, and a cross-section of 345 m². The plane line is of the “transition curve + circular curve” type, with a curvature radius of R = 885.852 m~906.298 m and a longitudinal slope of 0.35%.

2.1. Geological Conditions

The tunnel buried excavation portion, as illustrated in Figure 2, is in the transition zone between the continental edge’s shallow sea and the coastal alluvial plain. The sedimentary environment is relatively complicated, and there are geological changes between the tunnel’s left and right sides. The soil types in the formation, according to the exploration report, mostly comprise artificial backfill soil, silt and muddy soil, clay and silty clay, sandy gravel, weathered granite, etc., and their characteristics are shown in

Table 1. Soil layer characteristics in concealed excavation section.

Number	Soil Layer Name	Soil Layer Characteristics
Ⅰ	Artificial backfill	Composed primarily of granite residual soil layer, crushed stone, and fine sand, it is mildly to moderately compacted and has high permeability.
Ⅲ-1	Muck and mucky soil	Mostly silt and mucky sand, with consistent development, uniform thickness, flow plastic shape, high water content, high compressibility, and poor strength.
Ⅲ-2	Clay and silty clay	Mostly silty clay and silt, with iron–manganese nodules and a lot of sand, loose to slightly dense, with a thickness gradually increasing from west to east.
Ⅲ-3	Gravelly sand	Primarily made up of gravel sand, loose to slightly dense, well-permeable, and high in water content.
Ⅳ-3	Mucky soil	Mostly mucky soil and humus silt, flow plastic, steady development, high water content, high compressibility, low strength, with soft soil characteristics.
Ⅴ-1	Silty clay	Mostly silty clay, unevenly developed, with a maximum thickness of 4.5 m.
Ⅴ-2	Coarse and gravel sand	Mostly medium sand, coarse sand, and gravel sand, widely distributed, mainly in a medium-dense state, slightly sticky.
Ⅴ-3	Mucky soil	Mostly muddy soil with high water content, high compressibility, flow plastic, under consolidation, and low strength.
Ⅵ-2	Coarse and gravel sand	Mostly coarse and gravel sand, with uneven development and significant interface undulation, with medium density.
Ⅶ-1	Residual soil	Mostly weathered residual soil of granite, unevenly developed, moist, moderately dense.
Ⅷ-1	Weathered granite	Mostly made up of Biotite porphyry granite that has fully weathered. After absorbing water, the weathered rock mass has become a loose structure that is simple to soften.

.

The surface layer is primarily composed of soft soil layers such as silt and muddy soil, which have the characteristics of wide distribution, large thickness, and significant local changes; high compressibility and sensitivity; and a humid climate with long and abundant rainfall, which forms a large supply of groundwater. The mucky soil in layers III-1, IV-3, and V-3, as well as the silt layer with greater humus, have engineering features such as high water content, high compressibility, and poor strength, which are detrimental to subterranean engineering excavation and stability.

2.2. Freezing Program for the FSPR Method

A horizontal zoning and vertical segmentation freezing method is used to avoid an excessive volume of frozen soil, according to the freezing construction organization design papers and pertinent requirements of the Gongbei Tunnel. As indicated in Figure 3, horizontal zoning freezing is mostly used to collaborate with underground excavation construction by separating the pipe curtain section into five zones: A, B1, B2, B3, and C. Vertical zoning freezing entails establishing two freezing stations on the east and west sides of the tunnel near the working shafts at both longitudinal ends, each with two systems (a refrigeration system and a limit system). To manage the opening and shutting of freezing tubes and limiting tubes, the pipe curtain is separated into three zones: freezing zone 1 (64 m), freezing zone 2 (128 m), and freezing zone 3 (64 m) from east to west.

Excessive frost heave can easily result in project risks since the tunnel passes through the customs port complex, near a significant number of subterranean pipelines and building piles, and is primarily surrounded by water-rich and soft soil layers. There are stringent criteria for mitigating the influence of ground frost heave, and the thickness of the frozen soil curtain should not be excessive. To achieve the above two requirements, careful real-time monitoring of the freezing temperature field is required during the building process:

(1)

Water sealing between jacked pipes should be safe with a minimum design value of frozen soil curtain thickness. The minimum thickness of the frozen soil curtain is planned to be 2 m, considering design characteristics such as pipe diameter, average separation between neighboring pipes, and misalignment angle;

(2)

The maximum design value of frozen soil curtain thickness cannot induce severe frost heave surface deformation. Because of the proportional relationship between the volume of frozen soil and the amount of surface frost heave, if the thickness of the frozen soil curtain is too thick, it will cause damage or even destruction to adjacent buildings (structures) and surrounding pipelines, affecting Gongbei Customs’ daily operations. The design thickness of the frozen soil curtain in the upper part of the tunnel (Zones A and B1) does not exceed 2.3 m, and the design thickness of the frozen soil curtain in the middle and lower parts of the tunnel (Zone B2C) does not exceed 2.6 m, according to the “Preliminary Research Report on Frost Heave Deformation during the Construction of the Gongbei Tunnel with the FSPR Method” and combined with relevant freezing engineering experience.

During the freezing phase, pick the freezing tube opening sequence and keep the temperature between −25 °C and −30 °C. Determine if the thickness of the frozen soil curtain matches the design criteria based on real-time data monitoring findings of temperature measurement locations in the soil layer. When the thickness of the frozen soil curtain approaches the minimum design value, the refrigeration system’s operating parameters should be modified promptly to prevent excessive development of the frozen soil curtain and to guarantee that its maximum thickness is between 2.3 m and 2.6 m.

3. Coupling Mathematical Model for Tunnel Freezing–Frost Heave

3.1. Mathematical Model for the Freezing Temperature Field

The transient heat conduction issue with phase transition includes the fluctuation of the freezing temperature field in the horizontal freezing tunnel [28,29]. Based on heat transfer and permafrost theories, and assuming that the temperature gradient of the low-temperature refrigerant in the freezing tube along the longitudinal direction (length direction) of the tunnel is approximately zero during the freezing process [30,31], the two-dimensional heat balance control differential equation of the freezing temperature field can be expressed as follows:

$$\left\{\begin{array}{l}{C}_f\frac{\partial T_f}{\partial t}=\frac{\partial }{\partial x}\left({\lambda }_f\frac{\partial T_f}{\partial x}\right)+\frac{\partial }{\partial y}\left({\lambda }_f\frac{\partial T_f}{\partial y}\right), i{n}^{} {\Omega }_f\\ C_u\frac{\partial T_u}{\partial t}=\frac{\partial }{\partial x}\left({\lambda }_u\frac{\partial T_u}{\partial x}\right)+\frac{\partial }{\partial y}\left({\lambda }_u\frac{\partial T_u}{\partial y}\right), in {\Omega }_u\\ C_L\frac{\partial T}{\partial t}=\frac{\partial }{\partial x}\left({\lambda }_L\frac{\partial T}{\partial x}\right)+\frac{\partial }{\partial y}\left({\lambda }_L\frac{\partial T}{\partial y}\right)+L_w{\rho }_i\frac{\partial {\theta }_i}{\partial t}, i{n}^{} {\Omega }_L\end{array}\right.$$

(1)

where

• ${\Omega }_f,{\Omega }_u,{\Omega }_L$—${\mathrm{Frozen}}^{}\mathrm{zone},$ ${\mathrm{Unfrozen} \mathrm{zone}, }^{}\mathrm{Phase} \mathrm{transition} \mathrm{zone}$
• $C_f,C_u,C_L—\mathrm{Specific} \mathrm{heat} \mathrm{of} \mathrm{soil} \mathrm{volume} \mathrm{in} {\Omega }_f,{\Omega }_u,{\Omega }_L,{\mathrm{kJ}/(\mathrm{m}}^3\cdot °\mathrm{C})$
• $T_f,T_u—\mathrm{Soil} \mathrm{temperature} \mathrm{in} {\Omega }_f,{\Omega }_u,°\mathrm{C}$
• ${\lambda }_f,{\lambda }_u,{\lambda }_L—\mathrm{Soil} \mathrm{thermal} \mathrm{conductivity} \mathrm{in} {\Omega }_f,{\Omega }_u,{\Omega }_L,{\mathrm{W}/(\mathrm{m}}^3\cdot °\mathrm{C})$
• $L_w—\mathrm{Latent} \mathrm{heat} \mathrm{of} \mathrm{phase} \mathrm{transition} \mathrm{of} \mathrm{water},333\mathrm{kJ}/\mathrm{kg}$
• ${\rho }_i—\mathrm{The} \mathrm{density} \mathrm{of} \mathrm{ice},{\mathrm{kg}/\mathrm{m}}^3$
• ${\theta }_i—\mathrm{Volume} \mathrm{ice} \mathrm{content} \mathrm{of} \mathrm{soil},\%$
• $t—\mathrm{Freezing} \mathrm{time},\mathrm{s}$

Assuming that the soil’s freezing and melting temperatures are T_d and Tr, respectively, and that no phase transition occurs when the soil temperature T < T_d or T > T_r, the phase transition latent heat of water in Ω_L = [T_d, T_r] in Equation (1) can be derived and transformed into the phase transition latent heat of water-bearing soil, namely:

$$L_w{\rho }_i\frac{\partial {\theta }_i}{\partial t}=-L_w{\rho }_w\frac{\partial {\theta }_w}{\partial T}\frac{\partial T}{\partial t}=-L_w{\rho }_w\frac{{\rho }_d(w_0-w_u)}{{\rho }_w(T_r-T_d)}\frac{\partial T}{\partial t}=-\frac{L}{T_r-T_d}\frac{\partial T}{\partial t}$$

(2)

where

• ${\rho }_w,{\rho }_d—{\mathrm{Density} \mathrm{of} \mathrm{water}, \mathrm{Dry} \mathrm{density} \mathrm{of} \mathrm{soil}, \mathrm{kg}/\mathrm{m}}^3$
• ${\theta }_w—\mathrm{Volume} \mathrm{moisture} \mathrm{content} \mathrm{of} \mathrm{soil}, \%$
• $w_0,w_u—\mathrm{Initial} \mathrm{moisture} \mathrm{content} \mathrm{of} \mathrm{soil}, \mathrm{Unfrozen} \mathrm{water} \mathrm{content} \mathrm{in} \mathrm{frozen} \mathrm{soil}, \%$
• $L—\mathrm{Latent} \mathrm{heat} \mathrm{of} \mathrm{phase} \mathrm{transition} \mathrm{in} \mathrm{soil},{\mathrm{kJ}/\mathrm{m}}^3$

The differential equation for managing the heat balance of permafrost in the phase transition zone Ω_L = [T_d, T_r] can be found by substituting the derivation of Equation (2) into Equation (1):

$$\left(C_L+\frac{L}{T_r-T_d}\right)\frac{\partial T}{\partial t}=\frac{\partial }{\partial x}\left({\lambda }_L\frac{\partial T}{\partial x}\right)+\frac{\partial }{\partial y}\left({\lambda }_L\frac{\partial T}{\partial y}\right)$$

(3)

By linear interpolating according to the temperature interval, the mass volume-specific heat C_L and thermal conductivity λ_L in Equation (3) can be obtained:

$$C^{*}=\left\{\begin{array}{ccc}{C}_f& (T<T_d)& \\ \frac{C_f+C_u}{2}+\frac{L}{T_r-T_d}& (T_d\le T\le T_r)& \\ C_u& (T>T_d)& \end{array}\right.$$

(4)

$${\lambda }^{*}=\left\{\begin{array}{ccc}{\lambda }_f& (T<T_d)& \\ {\lambda }_f+\frac{{\lambda }_u-{\lambda }_f}{T_r-T_d}(T-T_d)& (T_d\le T\le T_r)& \\ {\lambda }_u& (T>T_d)& \end{array}\right.$$

(5)

where

• $C^{*}—\mathrm{Equivalent} \mathrm{volumetric} \mathrm{specific} \mathrm{heat} \mathrm{of} \mathrm{soil}$
• ${\lambda }^{*}—\mathrm{Equivalent} \mathrm{thermal} \mathrm{conductivity} \mathrm{of} \mathrm{soil}$

The two-dimensional heat balance control differential equation of the freezing temperature field may be consistently stated by synthesizing Equations (1)–(5):

$$C^{*}\frac{\partial T}{\partial t}=\frac{\partial }{\partial x}\left({\lambda }^{*}\frac{\partial T}{\partial x}\right)+\frac{\partial }{\partial y}\left({\lambda }^{*}\frac{\partial T}{\partial y}\right)$$

(6)

Equation (6)’s initial temperature requirement:

$${\left.T\right|}_{t=0}=T_0$$

(7)

where T₀ is the initial temperature of the soil, and a uniform temperature distribution is assumed. The temperature field in Equations (3)–(8) must meet continuity and energy conservation constraints at the phase transition boundary, i.e., at the frozen front S(t), namely:

$$T_f[S(t),t]=T_d$$

(8)

$${\lambda }_f\frac{\partial T_f}{\partial n}-{\lambda }_u\frac{\partial T_u}{\partial n}=L\frac{dS(t)}{dt}$$

(9)

where n is the phase transition boundary’s normal vector.

Because the freezing tube is significantly smaller than the tunnel and stratum, it may be treated as a single low-temperature point source, and the boundary conditions at the freezing tube are as follows:

$${\left.T\right|}_{(x_d,y_d)}=T_b(t)$$

(10)

where

• $x_d,y_d—\mathrm{Center} \mathrm{coordinate} \mathrm{of} \mathrm{the} \mathrm{freezing} \mathrm{tube},\mathrm{m}$
• $T_b—\mathrm{Temperature} \mathrm{of} \mathrm{refrigerant} \mathrm{inside} \mathrm{the} \mathrm{freezing} \mathrm{tube},°\mathrm{C}$

The boundary conditions for soil at an infinite distance from the tunnel section:

$${\left.T\right|}_{(x=\infty ,y=\infty )}=T_0$$

(11)

When considering the heat exchange between the surface and the surrounding atmosphere, the surface may be seen as a convective heat transfer boundary, namely:

$${\lambda }^{*}\frac{\partial T}{\partial n}=\alpha (T_a-T)$$

(12)

where

• $T_a—\mathrm{Atmospheric} \mathrm{temperature},°\mathrm{C}$
• $\alpha —\mathrm{Convective} \mathrm{heat} \mathrm{transfer} \mathrm{coefficient} \mathrm{between} \mathrm{atmosphere} \mathrm{and} \mathrm{surface},$ ${\mathrm{kJ}/(\mathrm{m}}^2\cdot \mathrm{s}\cdot °\mathrm{C})$

The preceding two-dimensional heat balance control differential equation, as well as its initial and boundary conditions, provide a mathematical description of the stratum’s transient freezing temperature field.

3.2. Coupling Mathematical Model for Freezing–Frost Heave

Thermal elastic–plastic analysis is performed with the features of the frozen soil curtain in mind [32,33]. Assuming frozen soil is an isotropic linear elastic material, its total strain {ε} may be thought of as a combination of elastic strain increment {ε}_e induced by load and instantaneous bulk strain increment {ε}_v caused by temperature change [34,35], namely:

$$\left\{\epsilon \right\}={\left\{\epsilon \right\}}_e+{\left\{\epsilon \right\}}_v$$

(13)

The elastic strain increment {ε}_e can be expressed:

$${\left\{\epsilon \right\}}_e=[D{]}^{-1}\left\{\sigma \right\}$$

(14)

where [D] is the elastic matrix, and its expression is as follows:

$$[D]=\frac{E(1-\mu )}{(1+u)(1-2\mu )}\left\{\begin{array}{ccc}1& \mu /(1-\mu )& 0\\ \mu /(1-\mu )& 1& 0\\ 0& 0& (1-2\mu )/[2(1-\mu )]\end{array}\right\}$$

(15)

The thermal strain component produced by temperature change may be described in the local coordinate system as follows:

$$\left\{\begin{array}{cc}{\epsilon }_{11}^T& {\epsilon }_{12}^T\\ {\epsilon }_{21}^T& {\epsilon }_{22}^T\end{array}\right\}=\left[\begin{array}{cc}1/3& 0\\ 0& 1/3\end{array}\right]{\left\{\epsilon \right\}}_v$$

(16)

Hooke’s elastic theorem states that the total strain component of frozen soil in the local coordinate system is expressed as follows:

$$\left\{\begin{array}{l}{\epsilon }_{11}=\frac{1}{E}\left[{\sigma }_{11}-\mu ({\sigma }_{22}+{\sigma }_{33})\right]+\frac{1}{3}{\left\{\epsilon \right\}}_v\\ {\epsilon }_{22}=\frac{1}{E}\left[{\sigma }_{22}-\mu ({\sigma }_{11}+{\sigma }_{33})\right]+\frac{1}{3}{\left\{\epsilon \right\}}_v\\ {\epsilon }_{12}={\epsilon }_{21}=\frac{2(1+\mu )}{E}{\tau }_{12}\end{array}\right.$$

(17)

For the planar strain issue, with ε₁₃ = ε₂₃ = ε₃₃ = 0, it can be obtained that:

$${\sigma }_{33}=\mu ({\sigma }_{11}+{\sigma }_{22})-\frac{E}{3}{\left\{\epsilon \right\}}_v$$

(18)

The formula for the total strain of frozen soil in the local coordinate system is found by inserting Equation (18) into Equation (17):

$$\left\{\begin{array}{l}{\epsilon }_{11}=\frac{1-{\mu }^2}{E}\left({\sigma }_{11}-\frac{\mu }{1-\mu }{\sigma }_{22}\right)+\frac{1+\mu }{3}{\left\{\epsilon \right\}}_v\\ {\epsilon }_{22}=\frac{1-{\mu }^2}{E}\left({\sigma }_{22}-\frac{\mu }{1-\mu }{\sigma }_{11}\right)+\frac{1+\mu }{3}{\left\{\epsilon \right\}}_v\\ {\epsilon }_{12}={\epsilon }_{21}=\frac{2(1+\mu )}{E}{\tau }_{12}\end{array}\right.$$

(19)

The stress field equilibrium equation in the global coordinate system for the planar strain issue is as follows:

$$\left\{\begin{array}{l}\frac{\partial {\sigma }_x}{\partial x}+\frac{\partial {\tau }_{xy}}{\partial y}=0\\ \frac{\partial {\tau }_{xy}}{\partial x}+\frac{\partial {\sigma }_y}{\partial y}=-\gamma \end{array}\right.$$

(20)

where γ is the soil gravity, and the geometric equation is

$$\left\{\begin{array}{l}{\epsilon }_x=\frac{\partial u}{\partial x}\\ {\epsilon }_y=\frac{\partial v}{\partial y}\\ {\epsilon }_{xy}=\frac{\partial u}{\partial x}+\frac{\partial v}{\partial y}\end{array}\right.$$

(21)

The displacement boundary conditions:

$$\left\{\begin{array}{l}u=\overline{u}(x,y);\mathrm{Horizontal} \mathrm{displacement}\\ v=\overline{v}(x,y);\mathrm{Vertical} \mathrm{displacement}\end{array}\right.$$

(22)

Given that mechanical properties of frozen soil—such as elastic modulus E and Poisson’s ratio μ—fluctuate with temperature, the thermodynamic coupling constitutive equation can be found by rewriting Equation (19):

$$\left\{\begin{array}{l}{\epsilon }_x=\frac{1-\mu {(T)}^2}{E(T)}\left[{\sigma }_x-\frac{\mu (T)}{1-\mu (T)}{\sigma }_y\right]+\frac{1+\mu (T)}{3}{\left\{\epsilon \right\}}_v\\ {\epsilon }_y=\frac{1-\mu {(T)}^2}{E(T)}\left[{\sigma }_y-\frac{\mu (T)}{1-\mu (T)}{\sigma }_x\right]+\frac{1+\mu (T)}{3}{\left\{\epsilon \right\}}_v\\ {\epsilon }_{xy}=\frac{2\left[1+\mu (T)\right]}{E(T)}{\tau }_{xy}\end{array}\right.$$

(23)

According to the incremental theory [36,37,38], the stress–strain relationship for the assumed isotropic elastic–plastic unfrozen soil is as follows [39]:

$$\left\{\begin{array}{l}d\left\{\sigma \right\}={[D]}_{ep}d\left\{\epsilon \right\}\\ {[D]}_{ep}=[D]-{[D]}_p\end{array}\right.$$

(24)

where [D]_ep is the elastic–plastic matrix and [D]_p is the plastic matrix. And, there is

$${[D]}_p=\frac{[D]\left\{\frac{\partial P}{\partial \sigma }\right\}{\left\{\frac{\partial F}{\partial \sigma }\right\}}^T[D]}{{\left\{\frac{\partial F}{\partial \sigma }\right\}}^T[D]\left\{\frac{\partial P}{\partial \sigma }\right\}+A}$$

(25)

where F represents the yield function, P represents the plastic potential function, and A represents the hardening function. The coupling mathematical model of frost heave comprises Equations (6), (20) and (24).

4. Establishment of Finite Element Model

4.1. Numerical Calculation Model

Due to the tunnel’s complicated construction and geological circumstances, one of the cross-sections and its surrounding strata were chosen as the numerical computation object, as illustrated in Figure 4. The model has an overall dimension of 120 m × 60 m, and the geological condition has been simplified and separated into seven levels from the surface down. To assure the convergence and correctness of the calculation results, a planar triangular element was utilized for mesh division, and mesh refinement was performed in the pipe curtain area. The computation model has 215,091 elements, with an average element quality of 0.8264. For coupling, the calculation module employs the Solid Mechanics module and the Porous Medium Heat Transfer module [40,41].

4.2. Selection of Model Materials and Parameters

The yield criteria are based on the Drucker–Prager criterion, which assumes that all soil layers are ideal elastic–plastic constitutive models. Each soil layer’s material properties are established using the “Test Report of Physical and Mechanical Parameters of Artificial Frozen Soil at Gongbei Port on the Zhuhai side of the Hong Kong-Zhuhai-Macao Bridge”. The thermophysical parameters are listed in

Table 2. Thermophysical parameters of soil layer materials.

Soil Layer		Density (kg·m−3)	Moisture Content (%)	Thermal Conductivity	Specific Heat Capacity	Latent Heat of Phase Transition	Freezing Point
				(W·m−1·K−1)	(J·kg−1·K−1)	(kJ·m−3)	(°C)
Artificial backfill	Unfrozen	1660	16.05	1.398	1420	36.32	−0.5
	Frozen			1.690	1360
Gravelly sand	Unfrozen	2000	13.54	1.217	1410	31.48	−0.4
	Frozen			1.758	1340
Mucky soil	Unfrozen	1820	47.6	1.485	1510	41.26	−1.1
	Frozen			1.772	1450
Silty clay	Unfrozen	2010	26.37	1.537	1560	56.24	−1.6
	Frozen			1.824	1470
Medium sand	Unfrozen	2020	17.92	1.439	1430	32.87	−1.5
	Frozen			1.756	1360
Gravel clay	Unfrozen	1890	31.98	1.442	1580	48.25	−1.5
	Frozen			1.790	1470
Weathered granite	Unfrozen	2040	17.31	1.317	1470	33.06	−1.7
	Frozen			1.887	1380
Concrete		2450	—	1.74	920	—
Jacked pipe		7800	—	44.5	318	—

.

Table 3. Mechanical parameters of materials.

Soil Layer		Elastic Modulus (MPa)	Poisson’s Ratio	Cohesion (kPa)	Internal Friction Angle (°)
Artificial backfill	Unfrozen	20	0.3	12	15
	Frozen	10.67 + 4.3|T|	0.26	1730	4.98
Gravelly sand	Unfrozen	25	0.25	0	24
	Frozen	19.3 + 4.7|T|	0.29	1860	5.69
Mucky soil	Unfrozen	4	0.35	16	14
	Frozen	10.67 + 7.8|T|	0.28	990	4.49
Silty clay	Unfrozen	28	0.32	16	10
	Frozen	4.67 + 6.6|T|	0.31	1770	7.05
Medium sand	Unfrozen	25	0.25	0	24
	Frozen	16 + 5.2|T|	0.3	1910	4.49
Gravel clay	Unfrozen	30	0.25	12.5	10
	Frozen	34 + 1.2|T|	0.3	1450	4.49
Weathered granite	Unfrozen	29	0.26	80	40
	Frozen	29 + 1.2|T|	0.3	900	6.83
Concrete		25 × 103	0.2	—	—
Jacked pipe		210 × 103	0.27	—	—

shows the mechanical properties of frozen soil, where the elastic modulus E of frozen soil is calculated using the interpolation procedure in the table, the freezing temperature T (°C) is negative, and the jacked pipe and concrete are considered elastic materials.

4.3. Boundary Conditions and Calculation Assumptions

(1)

As shown in Figure 5, set the initial temperature of all soil layers to 20 °C, and the cooling plans for the circular freezing tubes in the concrete pipe and the profiled freezing tubes in the hollow pipe shall be taken based on their corresponding freezing system temperatures. Based on the engineering ambient temperature, the surface boundary (boundary 1) is set as a convective heat transfer boundary, while the borders on both sides (boundaries 2 and 3) and the bottom (boundary 4) are set as adiabatic, as shown in Figure 4a;

(2)

Without considering the surface load, the surface boundary (boundary 1) is regarded as a free boundary. Roller support is installed on both sides (boundaries 2 and 3), and only horizontal displacement is restricted. To restrict horizontal and vertical movement, fixed constraints are established at the bottom (boundary 4);

(3)

Assuming that the thermal resistance between different soil layers is zero, and ignoring the effect of soil moisture migration on the freezing process. The freezing point of all soil layers is uniformly set at −1.5 °C for the convenience of studying the size and shape of the frozen soil curtain.

To begin, an initial geostress field equilibrium analysis is performed to eliminate the initial vertical displacement of soil caused by gravity stress. Second, before the frost heave analysis, the model is subjected to volume force (soil gravity) and corresponding displacement boundary conditions to analyze the internal stress field of the soil at this time, and the data are derived and applied to the subsequent frost heave analysis model as the initial stress field of the soil. The soil maintains its equilibrium state without initial displacement due to the combined action of the initial stress field and the external volume force (gravity field), i.e., the initial vertical displacement of the soil is eliminated.

5. Analysis and Discussion of Calculation Results

5.1. Formation Law of Frozen Soil Curtain

Figure 6 depicts the cloud map of the temperature field distribution throughout the freezing phase. To study the changing features of the frozen soil curtain, the soil layer region around the 6# hollow pipe and 7# concrete pipe were intercepted. The frozen soil boundary (−1.5 °C) is denoted by the black line in each illustration.

According to the cooling plan, only the circular freezing tubes in the concrete pipe were opened in the first 32 days, as shown in Figure 6a. As the freezing time progresses, frozen soil forms near the circular freezing tubes at the two waist positions of each concrete pipe, and its range gradually expands from the pipe to the inside and outside of the pipe curtain. On day 32, the soil temperature in the pipe curtain area significantly decreased, particularly around the concrete pipes. The frozen soil curtain region has begun to “wrap” the concrete pipe and “overlap” with the neighboring hollow pipe but its size and temperature have not yet met the design requirements.

The soil temperature in the pipe curtain region rapidly decreased after the opening of the profiled freezing tubes in the hollow pipe from day 32 to day 45, and the distribution of the frozen soil curtain exhibited substantial changes compared to previously, as shown in Figure 6b. A frozen soil curtain of varying thickness has formed around the concrete pipes, with an average thickness range of 1.80 m to 1.87 m; the part of the hollow pipe located on the inner side of the pipe curtain has also been “wrapped” by the frozen soil; while a small portion on the outside has not been “wrapped”, the inner frozen soil curtain is developing faster than the outer side. The frozen soil between neighboring pipes also develops quite quickly, with an average thickness range of 1.82 m to 1.91 m and a temperature range of −0.2 °C to −17.3 °C.

After the frozen soil curtain has completely “wrapped” the hollow pipe, the frozen soil range on the inner and outer sides of the pipe curtain continues to expand but at a slower rate, and the temperature of the frozen soil between the adjacent pipes decreases further. When frozen for 90 days, the development of the frozen soil curtain begins to stabilize, as shown in Figure 6c. The average thickness increased somewhat, ranging from 2.32 m to 2.58 m, and the temperature range of the frozen soil between pipes is −7.2 °C to −22.1 °C, meeting the freezing design requirements.

The interface and axial plane between 6# hollow pipe and 7# concrete pipe were chosen as research subjects to explore the spatiotemporal impacts of temperature variations in frozen soil between adjacent pipes, as illustrated in Figure 7a. The distance between interface points was 0.4 m, while the distance between axial plane points was 0.09 m.

Figure 7b depicts the temperature–time curves for each point on the interface. Because of the varying distances from point 0#, each curve shows “layering” and “symmetry” features from 0 d to 32 d. The average cooling rate at point 0# is −0.83 °C/d; the average cooling rates of Z3 and Z4, Z2 and Z5, and Z1 and Z6 are −0.74 °C/d, −0.57 °C/d, and −0.41 °C/d, respectively. When the profiled freezing tube in the hollow pipe was opened for operation after 32 days, all curves revealed a dramatic change in curvature, indicating that the cooling rate at each location had been increased once again. Simultaneously, because the position of the profiled freezing tube is closer to the inner side of the pipe curtain, the curve of Z3 is substantially lower than that of Z4. The average cooling rates of 0#, Z3, and Z4 were −0.76 °C/d, −0.85 °C/d, and −0.6 °C/d, respectively, during 32 to 45 days, indicating that the profiled freezing tube has a strong freezing strengthening effect. The total temperature range of the interface was −7.3 °C~−19.7 °C after 90 days, and the temperature on the inner side of the pipe curtain was slightly lower than that on the outer side. The overall distribution was relatively uniform.

Figure 7c depicts the temperature–time curve for each point on the axial plane. The difference from the interface temperature distribution with the average cooling rate in 32 d~45 d is that the closer the point is to the 6# hollow pipe, the higher the average cooling rate. The cooling rates of 1#, 2#, 0#, 3#, and 4# are −1.17 °C/d, −0.97 °C/d, −0.76 °C/d, −0.56 °C/d, and −0.34 °C/d, respectively. The temperature curve variations at each location eventually tend to be horizontal and overlap after 70 days, showing that the temperature differential between adjacent jacked pipes steadily narrows. By 90 days, the temperature at each location is nearly identical, with an overall value of −19.9 °C.

According to the above analysis, the process of forming a frozen soil curtain and the distribution features of the freezing temperature field are directly connected to the layout and opening time sequence of the two types of freezing tubes in hollow and concrete pipes. The results of the calculations demonstrate that the thickness and temperature distribution of the frozen soil curtain can match the design requirements, and that it has high water sealing and safety performance, producing ideal support conditions for future tunnel section excavation.

5.2. Distribution Law of Stratum Frost Heave Displacement

The size of the frozen soil curtain grows larger with time, causing vertical and horizontal displacement of the stratum due to frost heave. Figure 8 depicts the distribution features of the vertical frost heave displacement field. The highest vertical frost heave displacement at 32 d is 26.08 mm. The formation of the frozen soil curtain is relatively slow during the first 32 days because only the circular freezing tube inside the concrete pipe is opened for freezing; after 32 days, the formation rate accelerates due to the opening of the profiled freezing tube to strengthen the freezing. The maximum displacement is 75.5 mm after 45 days and 133.69 mm after 70 days. Following that, the size of the frozen soil curtain changes insignificantly, with a maximum displacement of 167.49 mm at 90 days.

Because of the model’s symmetry, the vertical frost heave displacement field is horizontally symmetrical in general, and its value increases with the time of freezing. The maximum positive value of vertical displacement always appears near the top of the pipe curtain, related to the tunnel’s shallow buried depth, the thin covering soil layer, and the comparatively large volume of frozen soil. The vertical frost heave displacement decreases as soil depth increases, indicating that the larger effective stress of deep soil partially offsets the frost heave force, which has a limiting effect on frost heave deformation until the vertical displacement at the bottom of the pipe curtain shows a small negative value. The vertical frost heave displacement in the tunnel excavation segment is quite minor due to the overall stiffness of the pipe curtain.

As shown in Figure 9, the horizontal frost heave displacement field, like the vertical frost heave displacement field, changes slightly within 0–32 days, with a maximum value of 20.64 mm at 32 d. It increased rapidly from 32 d to 70 d, with maximum values of 39.9 mm and 56.8 mm at 45 d and 70 d, respectively. The development then slows in 70 to 90 days, reaching a maximum of 72.3 mm at 90 days.

Overall, the stratum’s horizontal frost heave displacement field has a horizontally antisymmetric distribution, and its absolute magnitude increases with the duration of continuous freezing. In addition to the ground surface, areas with high displacement values are distributed beyond the pipe curtain’s waist, where there is primarily a silty clay layer with high water content and a high frost heave rate, causing significant deformation as the thickness of the frozen soil curtain increases during the freezing process. The horizontal frost heave displacement is quite minor in the tunnel excavation segment.

5.3. Distribution Law of Ground Surface Frost Heave Displacement

Figure 10a depicts the ground surface’s vertical frost heave displacement distribution curve. The tunnel centerline (x = 0 m) serves as the symmetry axis for all the curves. The vertical displacement of each point rises with freezing duration, and the effect is particularly noticeable in the −30 m ≤ x ≤ 30 m. range. The highest vertical uplift value of the surface at any moment in time is found at the center line, and the farther the location is from the center line, the lower the value is and ultimately approaches zero.

Within 32 days to 70 days, there is a considerable rise, with maximum values of 68.96 mm and 122.4 mm at 45 and 70 days, respectively, and 10.71 mm and 19.38 mm at ±30 m from the centerline. From 70 d to 90 d, the range of variation was rather minimal, with a maximum value of 155.67 mm at 90 d and 25.36 mm at 30 ±m from the center line.

The ground surface’s horizontal frost heave displacement distribution curve is shown in Figure 10b, with the tunnel centerline (x = 0 m) as the axis of symmetry—indicating an anti-symmetric distribution—and the displacement at the centerline being zero at all times. The displacement values on both sides of the centerline first gradually increase with distance, reaching the extreme value at x = ±15.74 m, and then gradually decreasing to zero.

The horizontal displacement of the ground surface changes relatively little from 0 d to 32 d, with an average increase of 0.26 mm/d and an extreme value of 8.17 mm at 32 d. It quickly increases from 32 to 45 days, with an average rate of 1.25 mm/d and an extreme value of 24.47 mm at 45 days. It slows down from 45 d to 70 d, with an average rate of 0.88 mm/d and an extreme value of 46.44 mm at 70 d; from 70 d to 90 d, with an average rate of 0.66 mm/d and an extreme value of 59.63 mm at 90 d.

Comparing the numerical calculation results of vertical frost heave displacement at the surface midpoint of the selected model cross-section with the field-measured values, as shown in Figure 11, the variation patterns of the two are basically consistent, with some points having deviations less than 15 mm. The displacement value fluctuates relatively little from 0 d to 10 d, with an average increase rate of 0.07 mm/d. The displacement value varied relatively little from 0 to 10 d, with an average growth rate of 0.07 mm/d, then gradually increased to 0.93 mm/d from 10 d to 32 d. The displacement curve changes significantly from 32 d to 45 d, with an average increase rate of 3.47 mm/d, which then decreases to 2.14 mm/d from 45 d to 70 d, and further decreases to 1.54 mm/d from 70 d to 90 d. The numerical calculation result was 155.67 mm at 90 d, whereas the measured datum was 157.41 mm, showing that the calculation results are acceptable.

The findings of the foregoing study show that the distribution and change of frost heave displacement during the freezing period are directly connected to the engineering geological circumstances, freezing scheme, forming process of frozen soil curtain, and pipe curtain structure.

6. Conclusions

(1)

The formation of the frozen curtain during the freezing process is mostly dependent on two types of freezing tubes to freeze the soil between the jacked pipes and fulfill the objective of sealing water. At 90 days, the thickness of the frozen soil curtain varies from 2.32 m to 2.58 m, and the lower and more homogeneous frozen soil temperature between adjacent pipes ensures that its strength fulfills the water-sealing safety requirements;

(2)

The frost heave displacement distribution and fluctuation throughout the freezing phase are highly connected to geological conditions, freezing scheme, frozen soil curtain development, and pipe curtain structure. The ground surface’s vertical frost heave displacement curve is normally distributed with the tunnel centerline (x = 0 m) as the symmetry axis, and the change is more noticeable in the range of −30 m ≤ x ≤ 30 m; the maximum value at any time is located at the centerline, which is 155.67 mm at 90 d. The horizontal displacement at the ground surface’s center line is always zero, while the values on both sides of the centerline progressively rise with distance, reaching an extreme value at x = 15.74 m, 59.63 mm at 90 d;

(3)

In this study, the geological conditions, temperature gradient, and numerical model are appropriately simplified, and the calculation object is the two-dimensional cross-section of the tunnel during the freezing phase. Although the calculation results are in good agreement with the measured data and have good guiding value for engineering construction, it is still a topic worth further research on how to more accurately reflect complex actual working conditions and focus on establishing more precise large-scale numerical models considering subsequent construction phases, such as tunnel excavation, support structure construction, and ground thawing settlement.