# Analytical Solution of Three-Dimensional Heaving Displacement of Ground Surface during Tunnel Freezing Construction Based on Stochastic Medium Theory

Authors: Hai-bing Cai; Meng-kai Li; Xiao-fang Li; Sheng Li; Rong-bao Hong
Journal: Advances in Civil Engineering
 (2021)

---

## Abstract

Many studies indicate that the heaving displacement during tunnel freezing construction can be regarded as a plane strain problem. However, due to the fact that many of the frozen pipes are placed obliquely in practical engineering, it should be a three-dimensional problem. Therefore, considering stochastic medium theory, the analytical solution of three-dimensional heaving displacement of ground surface during tunnel freezing construction is established, and the analytical solution is applied to the project case in Beijing. The calculation results show that, considering the inclination angle between the freezing pipe and the center axis of the tunnel, the maximum heaving displacement of the ground surface occurs at the tunnel center. When the inclination angle is equal to 0°, the analytical solution of three-dimensional ground surface heaving displacement is in good agreement with that of the two-dimensional prediction model, which further verified the reliability of the established analytical solution of the three-dimensional frost heave prediction model.


---

## Body

## 1. Introduction

On the basis of the sandbox model test, the Polish scholar Litwiniszyn [1] proposed a stochastic medium theory to investigate the displacement of rock strata and ground surface during underground mining in the 1950s; in this theory, rock formation is regarded as a kind of “stochastic medium”, and the deformation of rock formation and surface caused by underground mining is regarded as a stochastic process and the prediction of rock formation and surface deformation can be realized by using probability and statistics method. After the development and improvement of Chinese scholar Liu et al. [2, 3], the application field of stochastic medium theory has developed from the analysis of surface deformation caused by underground mining to the analysis of surface deformation caused by open-pit mining, near surface excavation, and stratum drainage. Yang et al. [4] and Yang et al. [5] analyzed and discussed in detail the surface movement and deformation caused by urban tunnel construction by using stochastic medium theory, further expanding the application field of stochastic medium theory. Due to the deformation of the surrounding rock caused by shield construction of subway tunnel in soft stratum environment, surface subsidence and shield segment rupture problems occur frequently [6]. Shi et al. [7, 8] considered the time-space development process of ground movement and deformation in tunnel excavation and, based on the stochastic medium theory, systematically researched the issues related to ground movement and deformation caused by mining method and shield method construction of urban tunnel. At present, the theory has become one of the practical methods to predict the surface deformation in urban metro tunnel engineering in China.If ground frost heave induced by soil freezing is considered the inverse process of ground settlement [9] caused by soil excavation, then ground frost heave during tunnel freezing construction can also be predicted via the stochastic medium theory. Based on the stochastic medium theory and shear displacement method, Zeng and Huang [10] established a calculation model of soil deformation caused by double-O-tube shield tunnel construction. Yang et al. [4], Li [11], Tao et al. [12], and Zhou et al. [13] used stochastic medium theory to analyze the surface frost heaving deformation caused by the freezing method of subway tunnel construction and obtained some meaningful conclusions. However, the common deficiency is that freezing process is not considered in the analysis process. In the study of Li et al [14], considering the frost heaving force and frost heaving deformation of noncircular tunnel with lining support, the analytical solution of frost heaving force is obtained by the theory of complex function and continuity condition. In the study of Huang et al. [15], the analytical solution of circular tunnel frost heaving force is derived on the basis of the nonuniform frost heaving of tunnel surrounding rock. In the study of Cai [16–18] based on single-pipe freezing theory, flat-panel freezing theory, and calculation equation of frost heave ratio of soil, considering the formation process of the frozen wall and using stochastic medium theory, a two-dimensional analytical solution of ground frost heave in horizontal freezing period of tunnel is established and the analytical solution is improved in this study by considering the freezing point of soil and assuming the constant surface temperature of the freezing pipe, the analytical solution methods for the radius of the freezing front and the outer radius of the expanded area are determined and applied to practical engineering cases, and then, the reliability of the analytical solution was verified by comparing with the field measurement results [19].In the actual construction of the subway tunnel freezing method, the freezing pipe group is generally designed to be inclined and radial, which is not strictly horizontal freezing [20–23]. Therefore, the research only on the three-dimensional frost heave deformation of the ground can reflect the actual situation of the project based on the analysis of the three-dimensional freezing temperature field of the ground.In this paper, an analytical solution method of three-dimensional frost heave deformation is established considering the inclined radial layout characteristics of the freezing pipe group in the actual tunnel construction. It is an improvement on the analytical solution of two-dimensional with the stochastic medium theory that the temperature field before the intersection of the frozen wall is calculated by single-pipe freezing theory and the temperature field after the intersection of the frozen wall is calculated by flat-panel freezing theory.


## 2. Analytical Solution for Ground Frost Heave

### 2.1. Ground Surface Settlement

Figure1 shows that the length, width, and thickness of the microelement aredε, dξ, anddη, respectively. The distance between the center of the microelement and the surface is η. According to the stochastic medium theory, the vertical displacement of the soil at any point above the tunnel caused by the excavation of a single microelement can be obtained as(1)dWex,y,z=1r2zexp−πr2zx2+y2dεdξdη,where rz is the convergence radius.Figure 1
Unit volume diagram.The global coordinate systemX,Y,Z is transformed into cylindrical coordinate systemR,θ,z. The equation for converting the two coordinate systems can be expressed as follows:(2)Z=zX=RcosθY=Rsinθ.Therefore, in cylindrical coordinates, the vertical displacement of the upper soil caused by a single microelement is written as(3)WeR,θ,z=1r2zexp−πr2zR2.Based on the above assumptions, the volume of the microelement remains unchanged, that is, the volume of the soil in the whole excavation area remains unchanged, then(4)εeR=∂Ue∂Rεeθ=UeRεeZ=∂We∂Z,(5)εeR+εeθ+εeZ=0,where εeR, εeθ, and εeZ is determined by the volumetric strain of the microelement in the cylindrical coordinate system, respectively.Substituting equations (3) and (4) into equation (5), the horizontal displacement Ue can be obtained by the boundary conditions that R=0and R⟶±∞, Ue=0.(6)Ue=Rr3zdrzdZexp−πr2zR2.In the Cartesian coordinate system, the horizontal displacement of the overlying stratum after the excavation of microelements is written as follows:(7)Uex,y,z=x2+y2r3zdrzdZexp−πx2+y2r2zdεdξdη.Supposingrz is a proportional function of depth, the convergence radius can be obtained as(8)r0=R=ηtanβ,(9)rz=R=η−ztanβ,where β is the main influence angle of the overlying soil layer and rz is the horizontal convergence radius in depth.In actual engineering, excavation of the soil will make the soil around the tunnel contract inward. If the tunnel is excavated in the areaΩ and the area after the excavation surface is contracted and becomes Ψ, then the size of the contraction deformation is Δ=Ω−Ψ.In accordance with the superposition principle, the vertical displacement of the ground surface at any point will be the cumulative displacement of the upper soil sinking displacement, which is caused by excavation of each microunit body. It can be written as follows:(10)Wx,y,z=∭Δtan2βη−z2exp−πtan2βη−z2x−ε2+y−ξ2dεdξdη.The horizontal displacement of the upper soil settlement is expressed as follows:(11)Uxx,y,z=∭Δ−x−εtan2βη−z2exp−πtan2βη−z2x−ε2+y−ξ2dεdξdηUyx,y,z=∭Δ−y−εtan2βη−z2exp−πtan2βη−z2x−ε2+y−ξ2dεdξdη.In accordance with equation (10), the three-dimensional model of the stochastic medium theory is utilized to predict the settlement of the upper soil after tunnel excavation. Two parameters need to be determined, namely, the main influence angle β of rock and soil layer and the deformation area Δ after tunnel excavation. Δ mainly affects the influence range of the soil layer, while β affects the shape of the settlement range.From the above equations for calculating stratum settlement, there are many triple integral calculations involved, and the integral function of each triple integral is not integral. Therefore, when solving these triple integrals, a numerical integration method should be adopted, which can be solved by using Mathematical and Maple software.


### 2.2. Analytical Solution for Three-Dimensional Ground Frost Heave

Figure2 shows that the angle α is located between the nth freezing pipe and the center axis of the tunnel. After the elapse of time t, as the cooling capacity is continuously transferred to the freezing pipe, the negative temperature in the freezing pipe will exchange heat with the surrounding soil layer, and the freezing will be diverged along the periphery of the freezing pipe, and the freezing radius of the freezing pipes gradually increases. The inner and outer freezing fronts of the frozen soil wall tend to rapidly smoothen. The frozen soil wall is assumed to expand uniformly from Rt to RΔt due to the frost heave characteristics of soil and the expanded area is Δ.Figure 2
3D layout of freezing pipe.In accordance with the three-dimensional model of stochastic medium theory, when in the Cartesian coordinate systemεξη, the vertical heaving displacement of the ground surface during the freezing time of tunnel can be expressed as follows:(12)Wx,y,z=−∭Δtan2βη−z2exp−πtan2βη−z2x−ε2+y−ξ2dεdξdη.Since the integral region is a hollow cylinder, the space rectangular coordinate system is transformed four times in order to solve the integral equation. Figure3(a) shows that the Cartesian coordinate system εξηxyz is converted into the Cartesian coordinate system ε1ξ1η1. The conversion equation of the two coordinate systems is(13)ε=ε1+X1ξ=ξ1η=η1+Z1.Figure 3
Coordinate transformation diagram: (a) Coordinate systemε1ξ1η1. (b) Coordinate system ε2ξ2η2.
(a)(b)Figure3(b) shows that the coordinate system ε1ξ1η1 is converted into the coordinate system ε2ξ2η2. And the center of the freezing pipe is connected with the center of the tunnel as an axis η2, the tangent of the freezing circle is taken as an axis ε2, and the axis ζ2 coincides with the direction of ζ1. The angle between the axis η and the axis η2 is set as γ. The conversion equation of the two coordinate systems is expressed as follows:(14)ε1,ξ1,η1,1=ε2,ξ2,η2,1cosγ0sinγ00100−sinγ0cosγ00001.Figure4(a) shows that the coordinate system ε3ξ3η3 is obtained by rotating the coordinate system ε2ξ2η2 around the axis ε2 with angle α. And the conversion equation of the two coordinate systems can be written by(15)ε2,ξ2,η2,1=ε3,ξ3,η3,110000cosαsinα00−sinαcosα00001.Figure 4
Coordinate transformation diagram: (a) Coordinate systemε3ξ3η3. (b) Coordinate system rφξ.
(a)(b)Figure4(b) shows that the integral region is a hollow cylinder, φ is the angle from oε3 to om in the polar coordinate roφ, and the coordinate system ε3ξ3η3 is converted into a cylindrical coordinate system rφζ for convenient calculation. The conversion equation of the two coordinate systems is as follows:(16)ε3=rcosφη3=−rsinφξ3=−ξ.From equations (12) to (16), the conversion equation from the Cartesian coordinate system εξη into the cylindrical coordinate system rϕζ can be expressed as follows:(17)ε=rcosφcosγ+ζsinα+rsinφcosαsinγ+X1,ξ=rsinφsinα−ζcosα,η=rcosφsinγ−ζsinα+rsinφcosαcosγ+Z1.In accordance with the substitution equation of the triple integral, equation (12) may be written as follows:(18)Wx,y,z=−∭Δttan2βrcosϕsinγ−ζsinα+rsinϕcosαcosγ+Z1−z2exp−πtan2βrcosφsinγ−ζsinα+rsinφcosαcosγ+Z1−z2X−rcosϕcosγ−ζsinα+rsinϕcosαsinγ−X12+y−rsinφsinα+ζcosα2rdrdφdζ.The integral area is a hollow cylinder, the radius is fromRt to RΔt, the length of the freezing pipe is l, and the vertical frost heave displacement of each point on the surface of the tunnel during the freezing period in the cylindrical coordinate system rφζ can be written as follows:(19)Wx,y,z=−∫RtRΔt∫02π∫0ltan2βrcosϕsinγ−ζsinα+rsinϕcosαcosγ+Z12exp−πtan2βrcosϕsinγ−ζsinα+rsinϕcosαcosγ+Z1−z2X−rcosφcosφ−ζsinα+rsinφcosαsinγ−X12+y−rsinφsinα+ζcosα2rdrdφdζ,where xyz are the coordinates of the points on the earth’s surface in the Cartesian coordinate system (xyz). And X1Z1 are the abscissa and vertical coordinates of the freezing pipe center in the Cartesian system (xyz), respectively. α is the angle between the freezing pipe and the central line of the tunnel, which is positive in elevation and negative in depression. β is the primary influence angle of ground heaving. rφζ are polar radius, polar angle, and vertical coordinates in cylindrical coordinates, respectively. And γ is the angle between the connecting line and the axis η between the center of the freezing pipe and the center of the tunnel. l is the length of the freezing pipe.


### 2.3. Analytical Solution before the Closure of the Frozen Wall

Suppose that a circle of freezing pipes is arranged around the tunnel, and the freezing pipes are arranged uniformly. Freezing pipes are arranged horizontally and the soil properties are the same. If the number of frozen pipes isn, in accordance with the superposition principle, the displacement of ground frost heave caused by tunnel freezing before the closure of the frozen wall is equal to the sum of ground frost heave displacement caused by n single-pipe freezing. Then the freezing radius of all freezing pipes at t time is Rt, and the frost heave radius is RΔt.The serial number of freezing pipes is recorded as 1, 2, ...,n. In the Cartesian coordinate system εξηxyz, the coordinate of the center of the nth freezing pipe is defined asxi,zi,0, which can be obtained as(20)Xi=Rcosi2πn,Zi=h−Rdsini2πn,where h is the depth between the tunnel center to the ground surface and Rd is the radius of the freezing pipe layout circle.Before calculating the frost heave of the ground surface caused by theith freezing pipe, the Cartesian coordinate system needs to be converted to the cylindrical coordinate system rφζ. The conversion equation is(21)ε=rcosϕcosγ+ζsinα+rsinϕcosαsinγ+Xi,(22)ξ=rsinϕsinα−ζcosα,(23)η=rcosϕsinγ−ζsinα+rsinϕcosαcosγ+Zi.Substituting equation (20) into equations (21)–(23), respectively, then(24)ε=rcosϕcosγ+ζsinα+rsinϕcosαsinγ+Xi,(25)ζ=rsinφsinα−ζcosα,(26)η=rcosϕsinγ−ζsinα+rsinϕcosαsinγ+Zi.Hence,(27)γ=i2πn.Substituting equation (27) into equations (23)–(25), respectively, then(28)ε=rcosϕcosi2πn+ζsinα+rsinϕcosαsini2πn+Rdcosi2πn,(29)ξ=rsinϕsinα−ζcosα,(30)η=rcosϕsini2πn−ζsinα+rsinϕcosαcosi2πn+Rdcosi2πn.From equation (19), under the condition of space problem, the vertical frost heave displacement of each point on the ground surface caused by the ith freezing pipe can be expressed as(31)Wx,y,z=−∫RtRΔt∫02π∫0ltan2βrcosϕsini2π/n−ζisinα+risinϕicosαcosi2π/n+Zi−z2exp−πtan2βricosφsini2π/n−ζsinα+risinφicosαcosi2π/n+Zi−z2X−ricosϕcosi2π/n−ζsinα+risinϕicosαsini2π/n−Xi2+y−risinφisinα+ζicosα2rdrdφdζ.Substituting equations (28)–(30) into equation (31), respectively, then(32)Wx,y,z=−∫RtRΔt∫02π∫0ltan2βrcosφsini2π/n−ζisinα+risinφicosαcosi2π/n+h−Rdsini2π/n−z2exp−πtan2βricosφsini2π/n−ζsinα+risinφicosαcosi2π/n+h−Rdsini2π/n−z2X−ricosφcosi2π/n−ζsinα+risinφicosαsini2π/n−Rdcosi2πn2+y−risinϕisinα+ζicosα2rdrdϕdζ.In accordance with the superposition principle, the displacement of ground frost heave before the closure of the frozen wall during tunnel freezing is equal to the sum of heaving displacement caused byn single pipes, which can be obtained by(33)Wx,y,z=−∑i=1n∫RtRΔt∫02π∫0ltan2βrcosφsini2π/n−ζisinα+risinφicosαcosi2π/n+h−Rdsini2π/n−z2exp−πtan2βricosφsini2π/n−ζsinα+risinφicosαcosi2π/n+h−Rdsini2π/n−z2X−ricosφcosi2πn−ζsinα+risinφicosαsini2πn−Rdcosi2πn2+y−risinφisinα+ζicosα2rdrdφdζ.Similarly, the horizontal frost heave displacement of the ground surface before the closure of frozen wall during tunnel freezing construction can be expressed by(34)Ux,y,z=−∑i=1n∫02π∫RtRΔt∫0l−x−ricosφicosi2πn−ζisinα+risinφicosαsini2πn−Rdcosi2πntan2βricosφisini2π/n−ζisinα+risinφicosαcosi2π/n+Rdsini2π/n−z3exp−πtan2βricosϕisini2πn−ζsinα+risinϕicosαcosi2πn+h−Rdsini2πn−z2x−ricosφicosi2πn−ζsinα+risinφicosαsini2πn−Rdcosi2πn2+y−risinϕisinα+ζicosα2rdrdϕdζ.(35)Ux,y,z=−∑i=1n∫02π∫RtRΔt∫0l−y−risimϕsinα+ζicosαtan2βricosϕisini2π/n−ζisinα+risinϕicosαcosi2π/n+Rdsini2π/n−z3exp−πtan2βricosϕisini2πn−ζsinα+risinφicosαcosi2πn+h−Rdsini2πn−z2x−ricosϕcosϕicosi2πn−ζsinα+risinϕicosαsini2πn−Rdcosi2πn2+y−risinφisinα+ζicosα2rdrdφdζ.


### 2.4. Analytical Solution after the Closure of the Frozen Wall

As the frozen soil columns gradually enlarge with prolonged freezing time, the adjacent frozen soil columns interconnect to form a continuous frozen soil wall around the tunnel. The inner and outer freezing fronts of the frozen soil wall tend to rapidly smoothen. Figure5 shows that the frozen soil wall is assumed to expand uniformly from Rt to RΔt due to the frost heave characteristics of soil and the expanded area is Δ. Ground frost heave is induced by the expansion effect of the entire frozen wall after the closure of the frozen wall.Figure 5
Coordinate system transformation after frozen wall intersection.The Cartesian coordinate systemεξη is converted to the cylindrical coordinate system rθζ and θ is the angle of polar coordinate ro’θ. The axis ξ is perpendicular to the outside of the paper, and the axis ζ is in the same direction as the axis ξ. Then the conversion equation is(36)ε=rcosθ,(37)η=h−rsinθ,(38)ξ=ζ.In accordance with the substitution equation of the triple integral, equation (12) may be written in the cylindrical coordinate system rφζ as(39)Wx,y,z=−∭Δtan2βh−rsinθ−z2exp−πtan2βh−rsinθ−z2x−rcosθ2+y−ζ2rdrdθdζ.When the angle between the freezing pipe and the horizontal direction isα, the integral region is approximately a hollow frustum, and the vertical heaving displacement of the ground surface is written by(40)Wx,y,z=∫RtRΔt∫02π∫−1/2Lcosα1/2Lcosαfr,θ,ζrdrdθdζ.(41)fr,θ,ζ=−tan2βh−rsinθ−z2exp−πtan2βh−rsinθ−z2x−rcosθ2+y−ζ2.


### 2.5. Solution for the Radius of the Freezing Front

According to the paper written by Cai and Jiang et al. [17, 24, 25], it can be known that the single-pipe freezing theory can be applied to two-dimensional heat conduction problems, while the flat-panel freezing theory is applied to one-dimensional semi-infinite heat conduction problems.Therefore, the radius of the freezing front is the radius of a single frozen soil column before the closure of the frozen wall, which can be solved using the single-pipe freezing theory.After the closure of the frozen wall, the radii of the freezing front are the inner and outer radii of the frozen wall, which can be solved using the flat-panel freezing theory.


#### 2.5.1. Single-Pipe Freezing Theory

Figure6 shows that the temperature field can be divided into frozen and unfrozen areas according to the location of the freezing front. The temperature of the soil in the frozen area is Tf and that in the unfrozen area is Tu, both of which are functions of time t and radial coordinate r. The radius of the frozen soil column is rt, and soil temperature at the freezing front is freezing point Td. The initial temperature of soil is T0 and the outer radius of the freezing pipe is r0. The surface temperature of the freezing pipe is assumed constant and its value is Tc.Figure 6
Schematic diagram of single tube freezing.For this two-dimensional heat conduction problem, the differential equation of heat transmission in frozen and unfrozen areas can be written as follows:(42)∂Tf∂t=αf∂2Tf∂r2+1r∂Tf∂rr0≤r≤rt,(43)∂Tu∂t=αu∂2Tu∂r2+1r∂Tu∂rrt≤r<∞,where αf and αu are the thermal diffusivity of frozen and unfrozen soils, respectively:(44)αf=kfcfρf,αu=kucuρu,where kf and ku are the thermal conductivity of frozen and unfrozen soils, respectively. cf and cu are the specific heat of frozen and unfrozen soils. ρf and ρu are the density of frozen and unfrozen soils.The initial conditions for the differential equations are as follows:(45)t=0,Tu=T0.The boundary conditions of the differential equations are as follows:(46)r=r0,Tf=Tc,(47)r=rt,Tf=Tu=Td,(48)r=∞,Tu=T0.The surface temperature of the freezing pipe is easily measured using a temperature sensor. Therefore, the surface temperature of the freezing pipe is assumed constant in this study. Equation (47) is selected as the boundary condition of the differential equation.The heat balance equation at the phase transition boundary is given by(49)kf∂Tf∂r−ku∂Tu∂rr=rt=Ldrtdt,where L is the volumetric latent heat of soil, which is given by(50)L=Lwρdw0−wu,where Lw is the latent heat of water, ρd is the dry density of the soil, w0 is the initial water content in the soil, and wu is the unfrozen water content in the frozen soil.Combining the above initial and boundary conditions to solve the partial differential equations, the temperature field distributions in frozen and unfrozen areas can be obtained as follows:(51)Tf=Tc+Td−TcEir02/4αft−Eir2/4αftEir02/4αft−Eirt2/4αftr0≤r≤rt,(52)Tu=T0+Td−T0Eir2/4αutEirt2/4αutrt≤r<∞,where Eix is the exponential integral function, which is given by(53)Eix=∫x∞e−ηηdη.Substituting equations (51) and (52) into equation (49) yields(54)kfTd−Tce−A2/4αfEir02/4αft−EiA2/4αf+kuTd−T0e−A2/4αuEiA2/4αu=A24L,where A is an undetermined constant that varies with time, which can be(55)rt=At.In accordance with single-tube freezing theory, the radius of the freezing front (i.e., the radius of the frozen soil column) may be solved by using equations (54) and (55) before the closure of the frozen wall when the surface temperature of the freezing pipe is constant.The adjacent frozen soil columns begin to interconnect, which indicates that the frozen wall has just closed, and the closure time of the frozen wall can be calculated by the following equation:(56)ti=l2A2,where tiis the closure time of the frozen wall and l is the space between adjacent freezing pipes.


#### 2.5.2. Flat-Panel Freezing Theory

Figure7 shows that the flat panel is assumed to generate heat exchange only with the soil on the right side, the temperature of the flat panel is TC′, and the temperature of the soil in the freezing and unfrozen areas is Tf and Tu′, respectively, both of which are functions of time t and coordinate x. The distance between the freezing front and the flat panel is Xt, and the temperature of soil at the freezing front is the freezing point Td. The initial temperature of soil is T0.Figure 7
Plate freezing schematic.For the flat-panel freezing problem, when0≤x≤Xt, the differential equations of the frozen area can be written as(57)∂Tf∂t=αf∂2Tf∂x2.Similarly, whenXt≤x<∞., the differential equation is(58)∂Tu∂t=αu∂2Tu∂x2.The initial conditions for the differential equations are as follows:(59)t=0,Tu=T0.The boundary conditions of the differential equations are as follows:(60)x=0,Tf=Th′,(61)x=Xt,Tf=Tu=Td,(62)x=∞,Tu=T0.In addition to the above boundary conditions, the heat balance equation at the phase transition boundaryXt can be expressed as(63)kf∂Tf∂x−ku∂Tu∂xx=Xt=LdXtdt.Combining the above initial and boundary conditions to solve the partial differential equations, the temperature field distributions in the frozen area can be obtained as follows0≤x≤Xt:(64)Tf=Tc′+Td−Tc′Φx/2αftΦXt/2αf.The temperature field distributions in the unfrozen area can be expressed as followsXt≤x<∞:(65)Tu=T0+Td−T01−Φx/2αut1−ΦXt/2αu,where Φy is Gaussian integral function, which is given by(66)Φy=2π∫0ye−η2dη.Substituting equations (64) and (65) into equation (63) yields(67)kfTd−Tc′e−B2/4αfαfΦB/2αf+kuTd−T0e−B2/4αuαu1−ΦB/2αu=π2BL,where B is a constant to be determined, and the following condition exists:(68)Xt=Bt.Therefore, the flat panel is regarded as the freezing pipe circle after the closure of the frozen wall. The radius of the freezing fronts is calculated as follows:(69)Rt=Rd+Bt,(70)Rct=Rd−Bt,where Rd is the radius of the freezing pipe circle, Rt is the radius of the outer freezing front, and Rct is the radius of the inner freezing front.The thickness of the frozen wallEt is obtained by(71)Et=Rt−Rct=2Bt.In accordance with the flat-panel freezing theory, the radii of the freezing fronts after the closure of the frozen wall can be solved by using equations (67) to (70) when the thermophysical parameters of the frozen and unfrozen soils are determined. It is generally considered that the flat-pane cooling temperature Tc′ is equal to the freezing tube wall temperature Tc.


### 2.6. Solution for the Outer Radius of the Expanded Area

The frost heave rate in engineering is generally expressed as(72)εf=Δhfh×100%,where εf is the frost heave ratio of the soil, Δhf is the longitudinal frost heave of the sample, and h is the original height of the sample.Liu [26] thinks that the frost heave at a certain point during the freezing process of the soil is the integral of the frost heave ratio of the frozen soil layer along the thickness of the frozen layer, which is not directly related to the frost heave force. According to its principle, the frozen wall expands evenly outward, then(73)RΔt=R1t+∫R1tRt1+εfdr=R1t+Rt−R1t1+εf.Many experimental studies on frost heave of soils show that the frost heave of soils is restrained by loads, and the frost heave ratio of soils decreases with the increase of overlying loads in exponential function.(74)εf=εf0exp−bP,where εf0 is the frost heave ratio of soil without load, P is the load, kPa, and b is empirical constant, kPa−1.Substituting equation (74) into equation (73) yields(75)RΔt=R1t+Rt−R1t1+εf0exp−bP.


## 2.1. Ground Surface Settlement

Figure1 shows that the length, width, and thickness of the microelement aredε, dξ, anddη, respectively. The distance between the center of the microelement and the surface is η. According to the stochastic medium theory, the vertical displacement of the soil at any point above the tunnel caused by the excavation of a single microelement can be obtained as(1)dWex,y,z=1r2zexp−πr2zx2+y2dεdξdη,where rz is the convergence radius.Figure 1
Unit volume diagram.The global coordinate systemX,Y,Z is transformed into cylindrical coordinate systemR,θ,z. The equation for converting the two coordinate systems can be expressed as follows:(2)Z=zX=RcosθY=Rsinθ.Therefore, in cylindrical coordinates, the vertical displacement of the upper soil caused by a single microelement is written as(3)WeR,θ,z=1r2zexp−πr2zR2.Based on the above assumptions, the volume of the microelement remains unchanged, that is, the volume of the soil in the whole excavation area remains unchanged, then(4)εeR=∂Ue∂Rεeθ=UeRεeZ=∂We∂Z,(5)εeR+εeθ+εeZ=0,where εeR, εeθ, and εeZ is determined by the volumetric strain of the microelement in the cylindrical coordinate system, respectively.Substituting equations (3) and (4) into equation (5), the horizontal displacement Ue can be obtained by the boundary conditions that R=0and R⟶±∞, Ue=0.(6)Ue=Rr3zdrzdZexp−πr2zR2.In the Cartesian coordinate system, the horizontal displacement of the overlying stratum after the excavation of microelements is written as follows:(7)Uex,y,z=x2+y2r3zdrzdZexp−πx2+y2r2zdεdξdη.Supposingrz is a proportional function of depth, the convergence radius can be obtained as(8)r0=R=ηtanβ,(9)rz=R=η−ztanβ,where β is the main influence angle of the overlying soil layer and rz is the horizontal convergence radius in depth.In actual engineering, excavation of the soil will make the soil around the tunnel contract inward. If the tunnel is excavated in the areaΩ and the area after the excavation surface is contracted and becomes Ψ, then the size of the contraction deformation is Δ=Ω−Ψ.In accordance with the superposition principle, the vertical displacement of the ground surface at any point will be the cumulative displacement of the upper soil sinking displacement, which is caused by excavation of each microunit body. It can be written as follows:(10)Wx,y,z=∭Δtan2βη−z2exp−πtan2βη−z2x−ε2+y−ξ2dεdξdη.The horizontal displacement of the upper soil settlement is expressed as follows:(11)Uxx,y,z=∭Δ−x−εtan2βη−z2exp−πtan2βη−z2x−ε2+y−ξ2dεdξdηUyx,y,z=∭Δ−y−εtan2βη−z2exp−πtan2βη−z2x−ε2+y−ξ2dεdξdη.In accordance with equation (10), the three-dimensional model of the stochastic medium theory is utilized to predict the settlement of the upper soil after tunnel excavation. Two parameters need to be determined, namely, the main influence angle β of rock and soil layer and the deformation area Δ after tunnel excavation. Δ mainly affects the influence range of the soil layer, while β affects the shape of the settlement range.From the above equations for calculating stratum settlement, there are many triple integral calculations involved, and the integral function of each triple integral is not integral. Therefore, when solving these triple integrals, a numerical integration method should be adopted, which can be solved by using Mathematical and Maple software.


## 2.2. Analytical Solution for Three-Dimensional Ground Frost Heave

Figure2 shows that the angle α is located between the nth freezing pipe and the center axis of the tunnel. After the elapse of time t, as the cooling capacity is continuously transferred to the freezing pipe, the negative temperature in the freezing pipe will exchange heat with the surrounding soil layer, and the freezing will be diverged along the periphery of the freezing pipe, and the freezing radius of the freezing pipes gradually increases. The inner and outer freezing fronts of the frozen soil wall tend to rapidly smoothen. The frozen soil wall is assumed to expand uniformly from Rt to RΔt due to the frost heave characteristics of soil and the expanded area is Δ.Figure 2
3D layout of freezing pipe.In accordance with the three-dimensional model of stochastic medium theory, when in the Cartesian coordinate systemεξη, the vertical heaving displacement of the ground surface during the freezing time of tunnel can be expressed as follows:(12)Wx,y,z=−∭Δtan2βη−z2exp−πtan2βη−z2x−ε2+y−ξ2dεdξdη.Since the integral region is a hollow cylinder, the space rectangular coordinate system is transformed four times in order to solve the integral equation. Figure3(a) shows that the Cartesian coordinate system εξηxyz is converted into the Cartesian coordinate system ε1ξ1η1. The conversion equation of the two coordinate systems is(13)ε=ε1+X1ξ=ξ1η=η1+Z1.Figure 3
Coordinate transformation diagram: (a) Coordinate systemε1ξ1η1. (b) Coordinate system ε2ξ2η2.
(a)(b)Figure3(b) shows that the coordinate system ε1ξ1η1 is converted into the coordinate system ε2ξ2η2. And the center of the freezing pipe is connected with the center of the tunnel as an axis η2, the tangent of the freezing circle is taken as an axis ε2, and the axis ζ2 coincides with the direction of ζ1. The angle between the axis η and the axis η2 is set as γ. The conversion equation of the two coordinate systems is expressed as follows:(14)ε1,ξ1,η1,1=ε2,ξ2,η2,1cosγ0sinγ00100−sinγ0cosγ00001.Figure4(a) shows that the coordinate system ε3ξ3η3 is obtained by rotating the coordinate system ε2ξ2η2 around the axis ε2 with angle α. And the conversion equation of the two coordinate systems can be written by(15)ε2,ξ2,η2,1=ε3,ξ3,η3,110000cosαsinα00−sinαcosα00001.Figure 4
Coordinate transformation diagram: (a) Coordinate systemε3ξ3η3. (b) Coordinate system rφξ.
(a)(b)Figure4(b) shows that the integral region is a hollow cylinder, φ is the angle from oε3 to om in the polar coordinate roφ, and the coordinate system ε3ξ3η3 is converted into a cylindrical coordinate system rφζ for convenient calculation. The conversion equation of the two coordinate systems is as follows:(16)ε3=rcosφη3=−rsinφξ3=−ξ.From equations (12) to (16), the conversion equation from the Cartesian coordinate system εξη into the cylindrical coordinate system rϕζ can be expressed as follows:(17)ε=rcosφcosγ+ζsinα+rsinφcosαsinγ+X1,ξ=rsinφsinα−ζcosα,η=rcosφsinγ−ζsinα+rsinφcosαcosγ+Z1.In accordance with the substitution equation of the triple integral, equation (12) may be written as follows:(18)Wx,y,z=−∭Δttan2βrcosϕsinγ−ζsinα+rsinϕcosαcosγ+Z1−z2exp−πtan2βrcosφsinγ−ζsinα+rsinφcosαcosγ+Z1−z2X−rcosϕcosγ−ζsinα+rsinϕcosαsinγ−X12+y−rsinφsinα+ζcosα2rdrdφdζ.The integral area is a hollow cylinder, the radius is fromRt to RΔt, the length of the freezing pipe is l, and the vertical frost heave displacement of each point on the surface of the tunnel during the freezing period in the cylindrical coordinate system rφζ can be written as follows:(19)Wx,y,z=−∫RtRΔt∫02π∫0ltan2βrcosϕsinγ−ζsinα+rsinϕcosαcosγ+Z12exp−πtan2βrcosϕsinγ−ζsinα+rsinϕcosαcosγ+Z1−z2X−rcosφcosφ−ζsinα+rsinφcosαsinγ−X12+y−rsinφsinα+ζcosα2rdrdφdζ,where xyz are the coordinates of the points on the earth’s surface in the Cartesian coordinate system (xyz). And X1Z1 are the abscissa and vertical coordinates of the freezing pipe center in the Cartesian system (xyz), respectively. α is the angle between the freezing pipe and the central line of the tunnel, which is positive in elevation and negative in depression. β is the primary influence angle of ground heaving. rφζ are polar radius, polar angle, and vertical coordinates in cylindrical coordinates, respectively. And γ is the angle between the connecting line and the axis η between the center of the freezing pipe and the center of the tunnel. l is the length of the freezing pipe.


## 2.3. Analytical Solution before the Closure of the Frozen Wall

Suppose that a circle of freezing pipes is arranged around the tunnel, and the freezing pipes are arranged uniformly. Freezing pipes are arranged horizontally and the soil properties are the same. If the number of frozen pipes isn, in accordance with the superposition principle, the displacement of ground frost heave caused by tunnel freezing before the closure of the frozen wall is equal to the sum of ground frost heave displacement caused by n single-pipe freezing. Then the freezing radius of all freezing pipes at t time is Rt, and the frost heave radius is RΔt.The serial number of freezing pipes is recorded as 1, 2, ...,n. In the Cartesian coordinate system εξηxyz, the coordinate of the center of the nth freezing pipe is defined asxi,zi,0, which can be obtained as(20)Xi=Rcosi2πn,Zi=h−Rdsini2πn,where h is the depth between the tunnel center to the ground surface and Rd is the radius of the freezing pipe layout circle.Before calculating the frost heave of the ground surface caused by theith freezing pipe, the Cartesian coordinate system needs to be converted to the cylindrical coordinate system rφζ. The conversion equation is(21)ε=rcosϕcosγ+ζsinα+rsinϕcosαsinγ+Xi,(22)ξ=rsinϕsinα−ζcosα,(23)η=rcosϕsinγ−ζsinα+rsinϕcosαcosγ+Zi.Substituting equation (20) into equations (21)–(23), respectively, then(24)ε=rcosϕcosγ+ζsinα+rsinϕcosαsinγ+Xi,(25)ζ=rsinφsinα−ζcosα,(26)η=rcosϕsinγ−ζsinα+rsinϕcosαsinγ+Zi.Hence,(27)γ=i2πn.Substituting equation (27) into equations (23)–(25), respectively, then(28)ε=rcosϕcosi2πn+ζsinα+rsinϕcosαsini2πn+Rdcosi2πn,(29)ξ=rsinϕsinα−ζcosα,(30)η=rcosϕsini2πn−ζsinα+rsinϕcosαcosi2πn+Rdcosi2πn.From equation (19), under the condition of space problem, the vertical frost heave displacement of each point on the ground surface caused by the ith freezing pipe can be expressed as(31)Wx,y,z=−∫RtRΔt∫02π∫0ltan2βrcosϕsini2π/n−ζisinα+risinϕicosαcosi2π/n+Zi−z2exp−πtan2βricosφsini2π/n−ζsinα+risinφicosαcosi2π/n+Zi−z2X−ricosϕcosi2π/n−ζsinα+risinϕicosαsini2π/n−Xi2+y−risinφisinα+ζicosα2rdrdφdζ.Substituting equations (28)–(30) into equation (31), respectively, then(32)Wx,y,z=−∫RtRΔt∫02π∫0ltan2βrcosφsini2π/n−ζisinα+risinφicosαcosi2π/n+h−Rdsini2π/n−z2exp−πtan2βricosφsini2π/n−ζsinα+risinφicosαcosi2π/n+h−Rdsini2π/n−z2X−ricosφcosi2π/n−ζsinα+risinφicosαsini2π/n−Rdcosi2πn2+y−risinϕisinα+ζicosα2rdrdϕdζ.In accordance with the superposition principle, the displacement of ground frost heave before the closure of the frozen wall during tunnel freezing is equal to the sum of heaving displacement caused byn single pipes, which can be obtained by(33)Wx,y,z=−∑i=1n∫RtRΔt∫02π∫0ltan2βrcosφsini2π/n−ζisinα+risinφicosαcosi2π/n+h−Rdsini2π/n−z2exp−πtan2βricosφsini2π/n−ζsinα+risinφicosαcosi2π/n+h−Rdsini2π/n−z2X−ricosφcosi2πn−ζsinα+risinφicosαsini2πn−Rdcosi2πn2+y−risinφisinα+ζicosα2rdrdφdζ.Similarly, the horizontal frost heave displacement of the ground surface before the closure of frozen wall during tunnel freezing construction can be expressed by(34)Ux,y,z=−∑i=1n∫02π∫RtRΔt∫0l−x−ricosφicosi2πn−ζisinα+risinφicosαsini2πn−Rdcosi2πntan2βricosφisini2π/n−ζisinα+risinφicosαcosi2π/n+Rdsini2π/n−z3exp−πtan2βricosϕisini2πn−ζsinα+risinϕicosαcosi2πn+h−Rdsini2πn−z2x−ricosφicosi2πn−ζsinα+risinφicosαsini2πn−Rdcosi2πn2+y−risinϕisinα+ζicosα2rdrdϕdζ.(35)Ux,y,z=−∑i=1n∫02π∫RtRΔt∫0l−y−risimϕsinα+ζicosαtan2βricosϕisini2π/n−ζisinα+risinϕicosαcosi2π/n+Rdsini2π/n−z3exp−πtan2βricosϕisini2πn−ζsinα+risinφicosαcosi2πn+h−Rdsini2πn−z2x−ricosϕcosϕicosi2πn−ζsinα+risinϕicosαsini2πn−Rdcosi2πn2+y−risinφisinα+ζicosα2rdrdφdζ.


## 2.4. Analytical Solution after the Closure of the Frozen Wall

As the frozen soil columns gradually enlarge with prolonged freezing time, the adjacent frozen soil columns interconnect to form a continuous frozen soil wall around the tunnel. The inner and outer freezing fronts of the frozen soil wall tend to rapidly smoothen. Figure5 shows that the frozen soil wall is assumed to expand uniformly from Rt to RΔt due to the frost heave characteristics of soil and the expanded area is Δ. Ground frost heave is induced by the expansion effect of the entire frozen wall after the closure of the frozen wall.Figure 5
Coordinate system transformation after frozen wall intersection.The Cartesian coordinate systemεξη is converted to the cylindrical coordinate system rθζ and θ is the angle of polar coordinate ro’θ. The axis ξ is perpendicular to the outside of the paper, and the axis ζ is in the same direction as the axis ξ. Then the conversion equation is(36)ε=rcosθ,(37)η=h−rsinθ,(38)ξ=ζ.In accordance with the substitution equation of the triple integral, equation (12) may be written in the cylindrical coordinate system rφζ as(39)Wx,y,z=−∭Δtan2βh−rsinθ−z2exp−πtan2βh−rsinθ−z2x−rcosθ2+y−ζ2rdrdθdζ.When the angle between the freezing pipe and the horizontal direction isα, the integral region is approximately a hollow frustum, and the vertical heaving displacement of the ground surface is written by(40)Wx,y,z=∫RtRΔt∫02π∫−1/2Lcosα1/2Lcosαfr,θ,ζrdrdθdζ.(41)fr,θ,ζ=−tan2βh−rsinθ−z2exp−πtan2βh−rsinθ−z2x−rcosθ2+y−ζ2.


## 2.5. Solution for the Radius of the Freezing Front

According to the paper written by Cai and Jiang et al. [17, 24, 25], it can be known that the single-pipe freezing theory can be applied to two-dimensional heat conduction problems, while the flat-panel freezing theory is applied to one-dimensional semi-infinite heat conduction problems.Therefore, the radius of the freezing front is the radius of a single frozen soil column before the closure of the frozen wall, which can be solved using the single-pipe freezing theory.After the closure of the frozen wall, the radii of the freezing front are the inner and outer radii of the frozen wall, which can be solved using the flat-panel freezing theory.


### 2.5.1. Single-Pipe Freezing Theory

Figure6 shows that the temperature field can be divided into frozen and unfrozen areas according to the location of the freezing front. The temperature of the soil in the frozen area is Tf and that in the unfrozen area is Tu, both of which are functions of time t and radial coordinate r. The radius of the frozen soil column is rt, and soil temperature at the freezing front is freezing point Td. The initial temperature of soil is T0 and the outer radius of the freezing pipe is r0. The surface temperature of the freezing pipe is assumed constant and its value is Tc.Figure 6
Schematic diagram of single tube freezing.For this two-dimensional heat conduction problem, the differential equation of heat transmission in frozen and unfrozen areas can be written as follows:(42)∂Tf∂t=αf∂2Tf∂r2+1r∂Tf∂rr0≤r≤rt,(43)∂Tu∂t=αu∂2Tu∂r2+1r∂Tu∂rrt≤r<∞,where αf and αu are the thermal diffusivity of frozen and unfrozen soils, respectively:(44)αf=kfcfρf,αu=kucuρu,where kf and ku are the thermal conductivity of frozen and unfrozen soils, respectively. cf and cu are the specific heat of frozen and unfrozen soils. ρf and ρu are the density of frozen and unfrozen soils.The initial conditions for the differential equations are as follows:(45)t=0,Tu=T0.The boundary conditions of the differential equations are as follows:(46)r=r0,Tf=Tc,(47)r=rt,Tf=Tu=Td,(48)r=∞,Tu=T0.The surface temperature of the freezing pipe is easily measured using a temperature sensor. Therefore, the surface temperature of the freezing pipe is assumed constant in this study. Equation (47) is selected as the boundary condition of the differential equation.The heat balance equation at the phase transition boundary is given by(49)kf∂Tf∂r−ku∂Tu∂rr=rt=Ldrtdt,where L is the volumetric latent heat of soil, which is given by(50)L=Lwρdw0−wu,where Lw is the latent heat of water, ρd is the dry density of the soil, w0 is the initial water content in the soil, and wu is the unfrozen water content in the frozen soil.Combining the above initial and boundary conditions to solve the partial differential equations, the temperature field distributions in frozen and unfrozen areas can be obtained as follows:(51)Tf=Tc+Td−TcEir02/4αft−Eir2/4αftEir02/4αft−Eirt2/4αftr0≤r≤rt,(52)Tu=T0+Td−T0Eir2/4αutEirt2/4αutrt≤r<∞,where Eix is the exponential integral function, which is given by(53)Eix=∫x∞e−ηηdη.Substituting equations (51) and (52) into equation (49) yields(54)kfTd−Tce−A2/4αfEir02/4αft−EiA2/4αf+kuTd−T0e−A2/4αuEiA2/4αu=A24L,where A is an undetermined constant that varies with time, which can be(55)rt=At.In accordance with single-tube freezing theory, the radius of the freezing front (i.e., the radius of the frozen soil column) may be solved by using equations (54) and (55) before the closure of the frozen wall when the surface temperature of the freezing pipe is constant.The adjacent frozen soil columns begin to interconnect, which indicates that the frozen wall has just closed, and the closure time of the frozen wall can be calculated by the following equation:(56)ti=l2A2,where tiis the closure time of the frozen wall and l is the space between adjacent freezing pipes.


### 2.5.2. Flat-Panel Freezing Theory

Figure7 shows that the flat panel is assumed to generate heat exchange only with the soil on the right side, the temperature of the flat panel is TC′, and the temperature of the soil in the freezing and unfrozen areas is Tf and Tu′, respectively, both of which are functions of time t and coordinate x. The distance between the freezing front and the flat panel is Xt, and the temperature of soil at the freezing front is the freezing point Td. The initial temperature of soil is T0.Figure 7
Plate freezing schematic.For the flat-panel freezing problem, when0≤x≤Xt, the differential equations of the frozen area can be written as(57)∂Tf∂t=αf∂2Tf∂x2.Similarly, whenXt≤x<∞., the differential equation is(58)∂Tu∂t=αu∂2Tu∂x2.The initial conditions for the differential equations are as follows:(59)t=0,Tu=T0.The boundary conditions of the differential equations are as follows:(60)x=0,Tf=Th′,(61)x=Xt,Tf=Tu=Td,(62)x=∞,Tu=T0.In addition to the above boundary conditions, the heat balance equation at the phase transition boundaryXt can be expressed as(63)kf∂Tf∂x−ku∂Tu∂xx=Xt=LdXtdt.Combining the above initial and boundary conditions to solve the partial differential equations, the temperature field distributions in the frozen area can be obtained as follows0≤x≤Xt:(64)Tf=Tc′+Td−Tc′Φx/2αftΦXt/2αf.The temperature field distributions in the unfrozen area can be expressed as followsXt≤x<∞:(65)Tu=T0+Td−T01−Φx/2αut1−ΦXt/2αu,where Φy is Gaussian integral function, which is given by(66)Φy=2π∫0ye−η2dη.Substituting equations (64) and (65) into equation (63) yields(67)kfTd−Tc′e−B2/4αfαfΦB/2αf+kuTd−T0e−B2/4αuαu1−ΦB/2αu=π2BL,where B is a constant to be determined, and the following condition exists:(68)Xt=Bt.Therefore, the flat panel is regarded as the freezing pipe circle after the closure of the frozen wall. The radius of the freezing fronts is calculated as follows:(69)Rt=Rd+Bt,(70)Rct=Rd−Bt,where Rd is the radius of the freezing pipe circle, Rt is the radius of the outer freezing front, and Rct is the radius of the inner freezing front.The thickness of the frozen wallEt is obtained by(71)Et=Rt−Rct=2Bt.In accordance with the flat-panel freezing theory, the radii of the freezing fronts after the closure of the frozen wall can be solved by using equations (67) to (70) when the thermophysical parameters of the frozen and unfrozen soils are determined. It is generally considered that the flat-pane cooling temperature Tc′ is equal to the freezing tube wall temperature Tc.


## 2.5.1. Single-Pipe Freezing Theory

Figure6 shows that the temperature field can be divided into frozen and unfrozen areas according to the location of the freezing front. The temperature of the soil in the frozen area is Tf and that in the unfrozen area is Tu, both of which are functions of time t and radial coordinate r. The radius of the frozen soil column is rt, and soil temperature at the freezing front is freezing point Td. The initial temperature of soil is T0 and the outer radius of the freezing pipe is r0. The surface temperature of the freezing pipe is assumed constant and its value is Tc.Figure 6
Schematic diagram of single tube freezing.For this two-dimensional heat conduction problem, the differential equation of heat transmission in frozen and unfrozen areas can be written as follows:(42)∂Tf∂t=αf∂2Tf∂r2+1r∂Tf∂rr0≤r≤rt,(43)∂Tu∂t=αu∂2Tu∂r2+1r∂Tu∂rrt≤r<∞,where αf and αu are the thermal diffusivity of frozen and unfrozen soils, respectively:(44)αf=kfcfρf,αu=kucuρu,where kf and ku are the thermal conductivity of frozen and unfrozen soils, respectively. cf and cu are the specific heat of frozen and unfrozen soils. ρf and ρu are the density of frozen and unfrozen soils.The initial conditions for the differential equations are as follows:(45)t=0,Tu=T0.The boundary conditions of the differential equations are as follows:(46)r=r0,Tf=Tc,(47)r=rt,Tf=Tu=Td,(48)r=∞,Tu=T0.The surface temperature of the freezing pipe is easily measured using a temperature sensor. Therefore, the surface temperature of the freezing pipe is assumed constant in this study. Equation (47) is selected as the boundary condition of the differential equation.The heat balance equation at the phase transition boundary is given by(49)kf∂Tf∂r−ku∂Tu∂rr=rt=Ldrtdt,where L is the volumetric latent heat of soil, which is given by(50)L=Lwρdw0−wu,where Lw is the latent heat of water, ρd is the dry density of the soil, w0 is the initial water content in the soil, and wu is the unfrozen water content in the frozen soil.Combining the above initial and boundary conditions to solve the partial differential equations, the temperature field distributions in frozen and unfrozen areas can be obtained as follows:(51)Tf=Tc+Td−TcEir02/4αft−Eir2/4αftEir02/4αft−Eirt2/4αftr0≤r≤rt,(52)Tu=T0+Td−T0Eir2/4αutEirt2/4αutrt≤r<∞,where Eix is the exponential integral function, which is given by(53)Eix=∫x∞e−ηηdη.Substituting equations (51) and (52) into equation (49) yields(54)kfTd−Tce−A2/4αfEir02/4αft−EiA2/4αf+kuTd−T0e−A2/4αuEiA2/4αu=A24L,where A is an undetermined constant that varies with time, which can be(55)rt=At.In accordance with single-tube freezing theory, the radius of the freezing front (i.e., the radius of the frozen soil column) may be solved by using equations (54) and (55) before the closure of the frozen wall when the surface temperature of the freezing pipe is constant.The adjacent frozen soil columns begin to interconnect, which indicates that the frozen wall has just closed, and the closure time of the frozen wall can be calculated by the following equation:(56)ti=l2A2,where tiis the closure time of the frozen wall and l is the space between adjacent freezing pipes.


## 2.5.2. Flat-Panel Freezing Theory

Figure7 shows that the flat panel is assumed to generate heat exchange only with the soil on the right side, the temperature of the flat panel is TC′, and the temperature of the soil in the freezing and unfrozen areas is Tf and Tu′, respectively, both of which are functions of time t and coordinate x. The distance between the freezing front and the flat panel is Xt, and the temperature of soil at the freezing front is the freezing point Td. The initial temperature of soil is T0.Figure 7
Plate freezing schematic.For the flat-panel freezing problem, when0≤x≤Xt, the differential equations of the frozen area can be written as(57)∂Tf∂t=αf∂2Tf∂x2.Similarly, whenXt≤x<∞., the differential equation is(58)∂Tu∂t=αu∂2Tu∂x2.The initial conditions for the differential equations are as follows:(59)t=0,Tu=T0.The boundary conditions of the differential equations are as follows:(60)x=0,Tf=Th′,(61)x=Xt,Tf=Tu=Td,(62)x=∞,Tu=T0.In addition to the above boundary conditions, the heat balance equation at the phase transition boundaryXt can be expressed as(63)kf∂Tf∂x−ku∂Tu∂xx=Xt=LdXtdt.Combining the above initial and boundary conditions to solve the partial differential equations, the temperature field distributions in the frozen area can be obtained as follows0≤x≤Xt:(64)Tf=Tc′+Td−Tc′Φx/2αftΦXt/2αf.The temperature field distributions in the unfrozen area can be expressed as followsXt≤x<∞:(65)Tu=T0+Td−T01−Φx/2αut1−ΦXt/2αu,where Φy is Gaussian integral function, which is given by(66)Φy=2π∫0ye−η2dη.Substituting equations (64) and (65) into equation (63) yields(67)kfTd−Tc′e−B2/4αfαfΦB/2αf+kuTd−T0e−B2/4αuαu1−ΦB/2αu=π2BL,where B is a constant to be determined, and the following condition exists:(68)Xt=Bt.Therefore, the flat panel is regarded as the freezing pipe circle after the closure of the frozen wall. The radius of the freezing fronts is calculated as follows:(69)Rt=Rd+Bt,(70)Rct=Rd−Bt,where Rd is the radius of the freezing pipe circle, Rt is the radius of the outer freezing front, and Rct is the radius of the inner freezing front.The thickness of the frozen wallEt is obtained by(71)Et=Rt−Rct=2Bt.In accordance with the flat-panel freezing theory, the radii of the freezing fronts after the closure of the frozen wall can be solved by using equations (67) to (70) when the thermophysical parameters of the frozen and unfrozen soils are determined. It is generally considered that the flat-pane cooling temperature Tc′ is equal to the freezing tube wall temperature Tc.


## 2.6. Solution for the Outer Radius of the Expanded Area

The frost heave rate in engineering is generally expressed as(72)εf=Δhfh×100%,where εf is the frost heave ratio of the soil, Δhf is the longitudinal frost heave of the sample, and h is the original height of the sample.Liu [26] thinks that the frost heave at a certain point during the freezing process of the soil is the integral of the frost heave ratio of the frozen soil layer along the thickness of the frozen layer, which is not directly related to the frost heave force. According to its principle, the frozen wall expands evenly outward, then(73)RΔt=R1t+∫R1tRt1+εfdr=R1t+Rt−R1t1+εf.Many experimental studies on frost heave of soils show that the frost heave of soils is restrained by loads, and the frost heave ratio of soils decreases with the increase of overlying loads in exponential function.(74)εf=εf0exp−bP,where εf0 is the frost heave ratio of soil without load, P is the load, kPa, and b is empirical constant, kPa−1.Substituting equation (74) into equation (73) yields(75)RΔt=R1t+Rt−R1t1+εf0exp−bP.


## 3. Case Analysis

The return line tunnel of a passenger transport station in Beijing, China, with a length of 40.0 m, a width of 6.0 m, and a buried depth of 13.0 m, is a shallow-buried tunnel with a large cross section. Given that tunnel will pass beneath existing pipelines, artificial ground freezing and undermining methods, instead of the open cut method, are selected for tunnel construction. Figure8 shows that the sand layer is frozen and reinforced, and the lower clay layer is closed, and the frozen wall is semicircular arched along the tunnel [27].Figure 8
Layout of the freezing pipes in the tunnel.After the active freezing period of the tunnel for 35 days, a nearly semicircular arched horizontal frozen wall with a thickness of 1.2–1.6 m and an arch span of 6.0 m was finally formed. The frozen wall was a nonclosed structure, and the average temperature of the frozen wall was −10°C. The layout of the freezing pipes in the tunnel is shown in Figure8.For the convenience of calculation, the simplified thermal and physical parameters of the tunnel are shown in Tables1 and 2.Table 1
Parameters of the tunnel.
ParametersValueTunnel radiusR2.25 mThe depth from tunnel center to ground surfaceh13.0 mRadius of freeze piper00.054 mRadius of freeze pipes circleRd3.25 mSpacing between freeze pipesl0.680 mNumber of freeze pipesn8Table 2
Parameters of soil.
ParametersValueDensityρ1958 kg/m3Dry densityρd1489 kg/m3Thermal conductivity of frozen soilkf31.5 kcal/(m·d·°C)Thermal conductivity of unfrozen soilku25.7 kcal/(m·d·°C)Specific heat of frozen soilcf0.30 kcal/(kg·°C)Specific heat of unfrozen soilcu0.45 kcal/(kg·°C)Water contentw023.0%Unfrozen water contentwu1.0%Freezing pointTd0°CSurface temperature of freezing pipeTc−25°CInitial temperature of soilT020°CFrost heave ratioεf0.56%The tangent value of the primary influence angle of the tunnel istanβ=0.8, and the convergence radius can be obtained by(76)rη=ηcotβ=18.75m.Therefore,−20m≤x≤20m and −10m≤y≤10m are taken as the calculation area. In accordance with equation (54), the thermal diffusivity of frozen and unfrozen soils may be obtained by(77)αf=kfcfρf=0.534m2/d,αu=kucuρu=0.0292m2/d.In accordance with equation (50), the volumetric latent heat of soil is(78)L=79.6×1489×0.23−0.01=26082Kcal/m3.The axisx is perpendicular to the freezing direction and the axis y is the direction along the freezing wall. And the frost heave radius and related parameters of different freezing time were introduced into the three-dimensional frost heave prediction model.Figure9 shows that the inclination angle α exists between the freezing pipe and the center axis of the tunnel. Because of the horizontal freezing, the inclination angle is equal to 0°, and the length of the frozen area is 20 m. In accordance with the Legendre–Gauss integration method, the numerical programs for the displacement of the ground surface are compiled using Maple software, and the results are presented in MATLAB mathematical software.Figure 9
Single-tube frozen longitudinal profile.The heaving displacement of the ground surface due to freezing for 30 days, 45 days, 60 days, 75 days, and 90 days is shown in Figures10–14, respectively.Figure 10
Heaving displacement of ground frost heave caused by freezing for 30 days (mm,α = 0°).Figure 11
Heaving displacement of ground frost heave caused by freezing for 45 days (mm,α = 0°).Figure 12
Heaving displacement of the ground surface caused by freezing for 60 days (mm,α = 0°).Figure 13
Heaving displacement of the ground surface caused by freezing for 75 days (mm,α = 0°).Figure 14
Heaving displacement of the ground surface caused by freezing for 90 days (mm,α = 0°).However, if the inclination angleα is not 0°, three-dimensional frost heave prediction model will be more complex in the same case. Since the inclination angle given in reference [28] is 4∘∼14∘, for the convenience of calculation, the inclination angle α is set to be 0°. Thus, the heaving displacement of the ground surface due to freezing for 30 days, 45 days, 60 days, 75 days, and 90 days is shown in Figures 15–19, respectively.Figure 15
Heaving displacement of ground frost heave caused by freezing for 30 days (mm,α = 10°).Figure 16
Heaving displacement of ground frost heave caused by freezing for 45 days (mm,α = 10°).Figure 17
Heaving displacement of ground frost heave caused by freezing for 60 days (mm,α = 10°).Figure 18
Heaving displacement of ground frost heave caused by freezing for 75 days (mm,α = 10°).Figure 19
Heaving displacement of ground frost heave caused by freezing for 90 days (mm,α = 10°).As shown in Figures10–14, when the inclination angle α between the freezing pipe and the center axis of the tunnel is 0°, the heaving displacement of the ground surface along the x-axis is in agreement with the results of the two-dimensional heaving displacement, which gradually decreases with an increase in distance from the tunnel center. Along the y-axis direction, that is, the tunnel axis direction, the heaving displacement of the ground surface in the middle part rarely changes, basically keeps horizontal. And the heaving displacement at both ends is obviously reduced. Since the length of the frozen pipe is 20 meters, the corresponding heaving displacement is the same, which can be regarded as the result of the translation of the same section. However, in two ends of the freezing pipe, the temperature of the freezing pipe is affected by the surrounding soil, the frost heave radius changes and the heaving displacement decreases.While the inclination angle is equal to 10°, the maximum heaving displacement of the ground surface does not occur at the central axis of the tunnel. Due to the influence of the inclination angle of the freezing pipe, one end of the freezing pipe is relatively close to the ground surface, the frost heave phenomenon will be more obvious, and the heaving displacement of the ground surface is larger than that of the two-dimensional prediction model.As shown in Figures20 and 21, the heaving displacement of the ground surface by freezing for 90 days is symmetrically distributed on the x-axis and y-axis at the inclination angle of 0°. However, in the angle of 10°, the heaving displacement of the ground surface is only symmetrical about the x-axis, and the center of the maximum displacement is shifted to the y-axis.Figure 20
Contour figure of ground frost heave caused by freezing for 90 days (mm,α = 0°).Figure 21
Contour figure of ground frost heave caused by freezing for 90 days (mm. 0 537α = 10°).


## 4. Conclusion

(1)
A three-dimensional frost heave prediction model for the tunnel freezing period is established based on the stochastic medium theory.(2)
The proposed three-dimensional frost heave prediction model is applied to an actual tunnel freezing project. Then, the three-dimensional heaving displacement of the ground surface is calculated. The results show that, considering the inclination angleα between the freezing pipe and the center axis of the tunnel, the maximum heaving displacement of the ground surface does not occur at the central axis of the tunnel. When the inclination angle is equal to 0°, the analytical solution of three-dimensional ground surface heaving displacement is in good agreement with that of the two-dimensional prediction model, which further verified the reliability of the established analytical solution of the three-dimensional frost heave prediction model.



---

*Source: 6668483-2021-01-30.xml*