Computer Modeling in Engineering & Sciences

DOI: 10.32604/cmes.2022.020388

ARTICLE

A Frost Heaving Prediction Approach for Ground Uplift Simulation Due to Freeze-Sealing Pipe Roof Method

Shengjun Deng1,2,3, Haolin Chen1, Xiaonan Gong2, Jiajin Zhou2, Xiangdong Hu4 and Gang Jiang1,*

1Institute of Geotechnical Engineering, Nanjing Tech University, Nanjing, 210009, China
2Research Center of Coastal and Urban Geotechnical Engineering, Zhejiang University, Hangzhou, 310058, China
3Department of Civil and Environmental Engineering, Princeton University, Princeton, NJ, 08544, USA
4Key Laboratory of Geotechnical and Engineering of Ministry of Education, Department of Geotechnical Engineering, Tongji University, Shanghai, 200092, China
*Corresponding Author: Gang Jiang. Email: g.jiang@njtech.edu.cn
Received: 20 November 2021 Accepted: 27 December 2021

1  Introduction

Freeze-sealing pipe roof (FSPR) method as an innovative pre-supporting method in tunnel engineering is applied in the Gongbei tunnel of Hongkong-Zhuhai-Macau Bridge, which is the first application in the world. Engineering applications using FSPR method has been appeared in recent years with the rapid development of underground space technology. FSPR method combines pipe-roofing method and artificial ground freezing method. The definition of FSPR method is that jacking pipes with large diameter are laid out in a circle along the cross section of tunnel in advance, then the artificial ground freezing method is adopted to freeze soil between jacking pipes to prevent water from entering into tunnel during excavation. The pipe jacking mainly plays a role of load bearing, while ground freezing mainly plays a role of water sealing [1–3]. However, ground surface uplift occurs during tunnel construction due to the ground freezing. The ground uplift is mainly caused by frost heave of frozen soil between steel pipes. Particular emphasis is placed on the ground surface deformation, therefore a rational method to predict the deformation of frost heave is required. Generally, the available frost heave model can be classified into three categories: the theoretical model, semi-empirical formula model and numerical model [4,5]. The theoretical model focuses on describing the physical changes in matter during the soil freezing process. The early proposal of this model was presented by Miller, who caught the attention of physicists and engineers by rigid ice model [6]. Subsequently, the model was further developed [7]. The semi-empirical formula model is not only used to explain the frost heaving phenomenon, but also to the frost heave calculation [8,9]. Typical models include segregation potential model proposed by Konrad [10–13], the rigid-ice model proposed by Miller and developed by O'Neill et al. [14], the coupled heat-fluid model proposed by Harlan [15], the thermal mechanical model proposed by Fromend et al. [16]. The numerical model is proposed based on combined thermal-hydraulic-mechanical field [17,18], elastic foundation beam theory [19] and other similar methods [20–22]. Also, numerous scholars have devoted their efforts to calculation methods for frost heave force, such as elastoplastic calculation model for surrounding rock proposed by Feng et al. [23], a mathematical mechanical model for tunnel in cold region proposed by Lai et al. [24], a coupled analytical approach for underground excavation problems proposed by Yang et al. [25]. A practical method was proposed by combining the advantages of laboratory testing or numerical simulation to simulate and predict the ground surface deformation during the AGF construction process [26,27]. From the perspective of engineering application, a major limitation in model analysis of frost heave is that model parameters have no practical physical meaning, and even if it is, it is difficult to measure. The accuracy of the model depends on the accuracy of the test parameters. There are some studies about the buried pipeline in frozen ground involving the test parameters [28,29], but lack a practical frost heaving calculation model [30]. Accordingly, a new model is proposed to predict the frost heave deformation. The model proposes the coefficient of cold expansion related with frost heaving rate from test as calculation parameter, which makes the model parameters correspond to the actual engineering. The ground surface uplift distribution of the Gongbei tunnel is analyzed based the proposed model. Meanwhile, the simulation result is verified by comparing with the field monitoring data.

2  The Proposed Frost Heaving Prediction Approach

ABAQUS as a three-dimensional finite element program is employed currently to carry out the numerical study. It is well known that temperature variation can lead to heat-induced expansion in frozen soil. Here we ignore the effect of temperature on strain ratio because it is small enough compared to frost heave. We propose the coefficient of cold expansion to simulate the expansion process, which is related to temperature. The basic framework of the model is as follows.

2.1 Equation of Thermal Equilibrium

For the isotropic homogeneous material, equation of thermal equilibrium is expressed as Eq. (1).

ρc∂T∂t=λ×∇2T+qv(1)

where ρ is mass density, c is specific heat, λ is thermal conductivity, t is the time, T is the temperature, qv is internal heat source.

2.2 A New Approach for Frost Heaving Prediction

(1)   One-dimensional model

For transient heat conduction without internal heat source in semi-infinite domain, the equation can be expressed as Eq. (2).

∂T∂t=a⋅∇2T=a⋅∂2T∂x2(2)

where a=λ/ρc

Fig. 1 is the temperature distribution at a special time in semi-infinite domain. Assume that T0 is the initial temperature for the whole domain, Tf(Tf<T0) is the given boundary temperature at x=0 when freezing begins, the thermophysical parameters are constant.

Figure 1: Temperature distribution at a special time t in semi-infinite domain

The equation of initial condition can be represented as Eq. (3).

T=T0,0≤x≤+∞,t=0(3)

The equation of boundary conditions can be represented as Eq. (4).

{T=Tf,x=0,t>0T=T0,x=+∞,t>0(4)

On the basis of θ=T−Tf, the following Eq. (5) can be obtained by the Laplace transformed and Eqs. (2)–(4).

θθ0=T−TfT0−Tf=erf(x2at)=erf(φ)(5)

where Gaussian error function erf(φ) can be expressed as Eq. (6).

erf(φ)=2π∫0φe−z2dz(6)

When x2at=2, T−TfT0−Tf≈1, T=T0, that is x=δ=4atδ, the temperature at this point has not changed exactly at the time tδ, which is the critical point of temperature change. Accordingly, the semi-infinite domain can be divided into two zones by the x=δ, shown in Fig. 1. One is the left zone between x=0 and x=δ, where temperature changes. The other one is the right zone start from x=δ to x=∞, For the temperature distribution at this time, it can be described by Eq. (7).

∇2T=∂2T∂x2=0(7)

The boundary conditions are as shown in Eq. (8).

{T=Tf,x=0T=T0,x=δ(8)

The Eq. (8) is substituted into Eq. (7), then the formula of temperature distribution is as follows:

T=Tf+xδ(T0−Tf)(0≤x≤δ)(9)

δ=4atδ is substituted into Eq. (9), which is equivalent to

T=Tf+x4atδ(T0−Tf)(0≤x≤δ)(10)

Hereafter, the formula of temperature change can be expressed as

ΔT=T−T0=(T0−Tf)⋅(x4atδ−1)(0≤x≤δ)(11)

According to the following statements, a new frost heaving prediction model is established. (1) The pore water in the soil becomes ice, causing the volume of the soil to expand. The moisture migration causing the frost heaving phenomenon is the main factor. The findings above are discovered by Taber and Beskow through the soil frost heaving test [31,32]. (2) The moisture migration is subject to the temperature gradient [33,34]. In this model, temperature gradient appears in the left zone between x=0 and x=δ of Fig. 1 where we assume that the frost heave occurs.

Therefore, the linear frost heaving rate can be calculated as follows:

α=η⋅∫0+∞|ΔT|dx=η⋅∫0δ|ΔT|dx(12)

In which, α is linear frost heaving rate, η represents the relative change of frost heave generated per unit temperature reduction in the x direction. The Eq. (13) can be further derived through substituting Eq. (11) into Eq. (12).

α=η⋅δ2(T0−Tf)=η⋅2atδ(T0−Tf)(13)

Based on the assumptions above, the maximum frost heaving rate can be expressed as Eq. (14).

αf=η⋅ΔTf=η⋅(T0−Tf)(14)

where αf is the maximum linear frost heaving rate.

Based on the frost heaving test of saturated soil under hydration conditions, αf of different soil samples can be measured. Moreover, the coefficient η of different soil samples corresponding to different soil layers can be calculated at a special initial temperature T0 and given boundary temperature Tf. According to the Eq. (14), the coefficient of cold expansion η corresponding to different soil layers can be expressed as follows:

η=αf/(T0−Tf)(15)

where we name the coefficient η(η>0) as the coefficient of cold expansion.

(2)   Two-dimensional model

For two-dimensional model, assume that (Tf−T0) is the temperature change at the surface of pipe (r=r0) when freezing begins. Position r=δ is the temperature change boundary, where the temperature change is zero. Accordingly, the two-dimensional temperature change distribution can be divided into two zones by the annulus r=δ, shown in Fig. 2. One is the internal zone between r=r0 and r=δ, where temperature changes. The other one is the external zone start from r=δ to ∞, where there is no change in temperature.

Figure 2: Temperature change distribution based on single pipe in two dimensions

According to the above, the equation of boundary conditions can be represented as

{ΔT=Tf−T0,r=r0,t>0ΔT=0,r=δ,t>0(16)

In Eq. (16), ΔT can be expressed as

ΔT=(Tf−T0)⋅ln⁡δrln⁡δr0(17)

In which, δ is the diameter of the temperature change boundary annulus, r is the diameter of the annulus, r0 is the diameter of pipe.

Therefore, the surface frost heaving rate can be calculated as follows:

β=η⋅∬|ΔT|dxdy=η⋅∫02πdθ∫r0δ|ΔT|rdr(18)

where β represents the surface frost heaving rate, θ represents the angle in polar coordinates, r represents the diameter in polar coordinates. The formula can be further derived through substituting Eq. (17) into Eq. (18).

β=η⋅(T0−Tf)⋅2πln⁡δr0[12r02ln⁡r0δ+14(δ2−r02)](19)

where the coefficient of cold expansion η as the basic characteristics of material can be obtained according to Eq. (15).

(3)   Three-dimensional model

The three-dimensional model is composed of multiple two-dimensional models along the longitudinal direction. The numerical simulation of frost heave is based on the model above.

(4)   Finite element calculation

Since the elastic modulus varies with temperature, the frost heave strain includes the elastic strain and the temperature strain.

Δεn=Δεne+ΔεnT(20)

where Δεn represents the frost heave strain, Δεne represents the elastic strain increment, ΔεnT represents the temperature strain increment.

The elastic strain increment can be calculated by Hooke's law as follows:

Δεne=1EeQΔσn(21)

Q=[1−μs−μs000−μs1−μs000−μs−μs10000002(1+μs)0000002(1+μs)0000002(1+μs)](22)

where Ee is the modulus, μs is the poisson's ratio.

Based on the obtained temperature field and the concept of linear expansion coefficient, temperature strain increment can be calculated. The results are as follows:

ΔεnT={αΔTn,αΔTn,αΔTn,0,0,0}T(23)

where α is the coefficient of cold expansion that is the expansion in ABAQUS software. ΔTn is the temperature change of the node.

At any period Δt, based on the basic assumption of elastic theory, the physical equation (stress-strain relationship) can be expressed as the following equation:

Δσn=(D¯n)(Δεn−ΔεnT)(24)

D¯n=EnQ−1(25)

The node force of the element

ΔFe=∭BTΔσdxdydz(26)

The stiffness matrix of the element

ke=∭BTD¯ndxdydz(27)

Thus, the whole equilibrium equation is as follows:

KΔδn=ΔPnL+ΔPnT(28)

In which, K is the global stiffness matrix, ΔPnL is the incremental node load caused by external load, ΔPnT is the incremental node load caused by temperature.

The displacement increment Δδn of each node can be obtained by the Eq. (28), and then substituted into the Eq. (24), the stress increment ΔQn of each node can be obtained.

3  Validation of the Frost Heaving Prediction Model: Gongbei Tunnel Case

3.1 Engineering Description

To verify the validity and feasibility of the proposed frost heaving prediction model, numerical results of the ground uplift due to freeze-sealing pipe roof method is compared with the field measured data. The freeze-sealing pipe roof (FSPR) method combining pipe roofing method and artificial ground freezing method is firstly applied in the Gongbei tunnel with the length of 255 m, which is a critical link of Hongkong-Zhuhai-Macau Bridge. Fig. 3 shows the FSPR structure, which consists of frozen soil and 36 pieces of 1620 mm diameter roofing pipes around the circumference of excavation. The spacing between the pipe edges ranges between 355 to 358 mm. In the normal direction of the tunnel circumference, the odd number roof pipes are shifted by 30 cm to the tunnel direction with respect to the even diameter pipes. The excavation cross-sectional area is approximately 345 m2 with 18.8 m in width and 20.6 m in height, which is the largest single tunnel excavation in China. The tunnel is buried within soft sandy and silty clay about 4.5 m below the ground surface. The purpose of ground freezing is preventing groundwater from entering the tunnel. A freezing process will last for 50 days before excavation. However, the ground surface uplift occurs due to the frost heave. And the amount of frost heave is simulated by the frost heaving prediction model. The correspondingly outcomes are employed to validate the proposed model.

Figure 3: The cross section of FSPR structure

3.2 Geometric Modeling

3.2.1 Geometric Model and Material

The ground deformation due to frost heave is simulated using the finite element software ABAQUS. It is assumed that the soil layers are horizontally distributed, isotropic and homogeneous, while roofing pipes and concrete are isotropic homogeneous elastic materials. Refer to the detailed investigation survey data of the Gongbei tunnel, the adjacent soil layers are combined into one layer for calculation. There are five soil layers in the model. The geometry of the model (255 m × 150 m × 60 m in length × width × height) and its cross section are respectively shown in Figs. 4 and 5. The upper surface of the model is ground surface, which is freely constrained. The bottom surface is fully constrained, and the two sides parallel to the tunnel axis are bounded by the normal direction in the x direction. The two sides perpendicular to the tunnel axis are the working well edges constrained by both y and z directions. Tie constraint is adopted between the roofing pipes and surrounding soil. According to the temperature of Zhuhai City in June 2016, the initial temperature of the model is set as 28°C. Since the positive temperature does not cause frost heaving, the value of the initial temperature has little effect on the frost heave calculation. During the freezing phase, the outer surface temperature of roofing pipe is set as −15°C. The soil along the coastal areas expresses elastoplastic properties and can be considered as homogeneous materials. According to the Subsection 2.2 Eq. (15), the coefficient of cold expansion η is obtained by the maximum frost heaving rate αf of the freezing test for each soil layer under the hydration condition when the initial temperature T0 = 28°C and given boundary temperature Tf = −15°C. Thus, the mechanical and thermophysical parameters of the model can be respectively found in Tables 1 and 2.

Figure 4: Geometry of the Gongbei tunnel model

Figure 5: Layout of the roofing pipes in layered soil

In which, the maximum frost heaving rate is obtained from a model calibration based on soil frost heave test under hydration conditions for each soil layer in laboratory. And the coefficient of cold expansion is calculated according to Eq. (15).

3.2.2 Feature Sections

In order to facilitate the analysis of ground surface uplift law, the feature sections are selected, including lateral direction and longitudinal direction. All sections are symmetrically distributed. There are 24 lateral monitoring sections (YK + 390∼YK + 630) along the longitudinal direction, where YK + 400 is the first monitoring section, 10 m from the boundary. Meanwhile, the monitoring data on the corresponding section (YK + 400, YK + 430, YK + 440) is obtained separately, shown in Fig. 6.

Figure 6: Feature sections

3.3 Finite Element Model

The finite element model consists of soil part, pipe-roofing part, concrete part, and lining part, with a total of approximately 489016 elements, including C3D8T elements and C3D6T elements, shown in Fig. 7. The soil in the coastal area exhibits elastoplastic properties and is considered as a homogeneous material. The Mohr-Coulomb model is used for soil. The pipe roofing, concrete and lining are simulated by an isotropic elastic model. Considering the influence of temperature and frost heave, the coupled temperature-displacement element is chosen and the implementation of frost heave depends on the coefficient of thermal expansion in the ABAQUS software.

Figure 7: Finite element model

There are 5 steps throughout the whole construction process, containing initial ground stress balance, pipe jacking, backfilling pipe with concrete, soil freezing, tunnel excavation. The ground stress is first released 40% during excavation, and then fully released after support. The main calculation process in Fig. 8 is as follows:

Figure 8: Simulation process

3.4 Model Validation

3.4.1 Temperature Field

Temperature is one of the factors affecting the frost heave effect, which indirectly affects the soil deformation, therefore the temperature field is the basis for calculating the soil deformation. The initial temperature of the model is 28°C, and the positive temperature does not cause frost heaving. Thus, the value of the initial temperature has little effect on the frost heave calculation. After freezing process of 50 days, the thickness of the frozen soil reaches the required width of 2 m. The temperature distribution is shown in Fig. 9.

Figure 9: Temperature distribution when freezing for 50 days

3.4.2 Ground Surface Uplift

The ground equilibrium is achieved by ABAQUS when the surface is horizontal. Then the frost heaving model is applied to predict the vertical uplift of ground surface due to frost heave. Fig. 10 shows the deformation of ground surface uplift after 50 days of freezing. The maximum uplift appears at the center of the model with the value of 285.6 mm, gradually decreasing along the longitudinal and lateral directions. The results along the longitudinal direction are mainly affected by the constraints of the two sides, while the one along the lateral direction is caused by the temperature field distribution.

Figure 10: Deformation drawing of ground surface uplift (Zoom in 20 times)

Fig. 11 illustrates the comparison between the results obtained by the model prediction and the field monitoring data at cross section YK2 + 400, YK2 + 430, YK2 + 440 after 50 days of freezing. Both the model simulation results and field monitoring data show that the maximum value of ground surface uplift appears at the symmetry axis of the tunnel section, and the amount of frost heaving on both sides gradually becomes smaller. The curve is shaped like normal distribution. The maximum values of the prediction model curve and the field measured curve at section YK2 + 400 are 75.0 and 72.7 mm, the values of YK2 + 430 are 180.0 and 168.2 mm, and the values at YK2 + 440 are 195.8 and 178.5 mm, respectively. Therefore, the two curves match well in terms of data and shape. Due to the different freezing scheme during different construction period, the curve of the field measurement and model simulation is asymmetric.

Figure 11: Comparison of model prediction and field monitoring along the cross section

The vertical uplift of ground surface after 50 days of freezing due to frost heave along the longitudinal direction of the tunnel is also analyzed, shown in Fig. 12. The shape of the curve is similar to the parabola with the opening downward. The vertical uplift in the middle span of the tunnel is larger, and the amount of ground uplift near the working wells at both ends is gradually reduced because of the restraining effect from the two sides. Due to the different excavation progress, when the excavation speed is slow, the freezing intensity is relatively high, and the frost heave is relatively large, and vice versa. Therefore, there is a big difference between the field measurement and numerical simulation at 100 and 200 m in the longitudinal direction of the tunnel.

Figure 12: Comparison of model prediction and field monitoring along the longitudinal

A good agreement is discovered between the model prediction and field measured data, except for some misfits that are due to inhomogeneous soil and local grouting that affecting heat transfer and groundwater seepage. In addition, the weight and dynamic load of the surface buildings make the monitoring data slightly smaller than the model result in some local location. The favorable agreement indicates that the proposed frost heaving prediction model can well predict the ground surface uplift due to frost heave.

4  Conclusions

The ground deformation due to FSPR method involves a complicated hydrothermal coupling process, which is difficult to accurately simulate. According to that, a frost heaving prediction approach is proposed to simulate ground deformation of the Gongbei tunnel. The proposed model deals with the frost heave calculation by combining the coefficient of cold expansion based on the calculation equation and the frost heaving test under hydration conditions. From the result analysis of prediction model and field monitoring, some conclusions are drawn:

(1)   The coefficient of cold expansion in model and the frost heaving rate from the test under the hydration condition can achieve a good correspondence that makes the calculation result closer to the actual engineering.

(2)   The numerical simulation and monitoring data of the Gongbei tunnel are employed to validate the proposed model, and a good agreement between the frost heaving prediction model and the field measured data verifies the rationality and applicability of the proposed model to predict the ground deformation.

(3)   The maximum uplift of the Gongbei tunnel appears at the center of the model, gradually decreasing along the lateral and longitudinal directions. The result in the lateral direction represents a normal distribution due to the influence of the constraint of two sides, while the one along the lateral direction shapes like a parabola with the opening downward due to the temperature field distribution.

It should be noted that the application of the proposed frost heaving prediction model to the Gongbei tunnel will be a reference for similar cases.

Acknowledgement: The author wishes to thank the high-performance computing cluster of Computer Science Department at Princeton University and the editor Dr. Lynn Wang for the careful format review.

Funding Statement: This work is supported by the financial support from National Natural Science Foundation of China (No. 51478340), Natural Science Foundation of Jiangsu Province (No. BK20200707), The Natural Science Foundation of the Jiangsu Higher Education Institutions of China (No. 20KJB560029), China Postdoctoral Science Foundation (No. 2020M671670), Key Laboratory of Soft Soils and Geoenvironmental Engineering (Zhejiang University), Ministry of Education (No. 2020P04), the support above is gratefully acknowledged.

Conflicts of Interest: The authors declare that they have no conflicts of interest to report regarding the present study.

This work is licensed under a Creative Commons Attribution 4.0 International License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.