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.
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
Figure 1 shows that the length, width, and thickness of the microelement are, , and, 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 aswhere is the convergence radius.
The global coordinate system is transformed into cylindrical coordinate system. The equation for converting the two coordinate systems can be expressed as follows:
Therefore, in cylindrical coordinates, the vertical displacement of the upper soil caused by a single microelement is written as
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, thenwhere , , and 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 can be obtained by the boundary conditions that and , .
In the Cartesian coordinate system, the horizontal displacement of the overlying stratum after the excavation of microelements is written as follows:
Supposing is a proportional function of depth, the convergence radius can be obtained aswhere is the main influence angle of the overlying soil layer and 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:
The horizontal displacement of the upper soil settlement is expressed as follows:
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
Figure 2 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 to due to the frost heave characteristics of soil and the expanded area is .
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:
Since the integral region is a hollow cylinder, the space rectangular coordinate system is transformed four times in order to solve the integral equation. Figure 3(a) shows that the Cartesian coordinate system is converted into the Cartesian coordinate system . The conversion equation of the two coordinate systems is
(a)
(b)
Figure 3(b) shows that the coordinate system is converted into the coordinate system . And the center of the freezing pipe is connected with the center of the tunnel as an axis , the tangent of the freezing circle is taken as an axis , and the axis coincides with the direction of . The angle between the axis and the axis is set as . The conversion equation of the two coordinate systems is expressed as follows:
Figure 4(a) shows that the coordinate system is obtained by rotating the coordinate system around the axis with angle . And the conversion equation of the two coordinate systems can be written by
(a)
(b)
Figure 4(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 is converted into a cylindrical coordinate system for convenient calculation. The conversion equation of the two coordinate systems is as follows:
From equations (12) to (16), the conversion equation from the Cartesian coordinate system into the cylindrical coordinate system can be expressed as follows:
In accordance with the substitution equation of the triple integral, equation (12) may be written as follows:
The integral area is a hollow cylinder, the radius is from to , 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 can be written as follows:where are the coordinates of the points on the earth’s surface in the Cartesian coordinate system (xyz). And 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. 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 is n, 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 , and the frost heave radius is .
The serial number of freezing pipes is recorded as 1, 2, ..., n. In the Cartesian coordinate system , the coordinate of the center of the nth freezing pipe is defined as, which can be obtained aswhere is the depth between the tunnel center to the ground surface and is the radius of the freezing pipe layout circle.
Before calculating the frost heave of the ground surface caused by the ith freezing pipe, the Cartesian coordinate system needs to be converted to the cylindrical coordinate system . The conversion equation is
Substituting equation (20) into equations (21)–(23), respectively, then
Hence,
Substituting equation (27) into equations (23)–(25), respectively, then
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
Substituting equations (28)–(30) into equation (31), respectively, then
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 by n single pipes, which can be obtained by
Similarly, the horizontal frost heave displacement of the ground surface before the closure of frozen wall during tunnel freezing construction can be expressed by
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. Figure 5 shows that the frozen soil wall is assumed to expand uniformly from to 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.
The Cartesian coordinate system is converted to the cylindrical coordinate system 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
In accordance with the substitution equation of the triple integral, equation (12) may be written in the cylindrical coordinate system as
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
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
Figure 6 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 and that in the unfrozen area is , both of which are functions of time t and radial coordinate r. The radius of the frozen soil column is , and soil temperature at the freezing front is freezing point . The initial temperature of soil is and the outer radius of the freezing pipe is . The surface temperature of the freezing pipe is assumed constant and its value is .
For this two-dimensional heat conduction problem, the differential equation of heat transmission in frozen and unfrozen areas can be written as follows:where and are the thermal diffusivity of frozen and unfrozen soils, respectively:where and are the thermal conductivity of frozen and unfrozen soils, respectively. and are the specific heat of frozen and unfrozen soils. and are the density of frozen and unfrozen soils.
The initial conditions for the differential equations are as follows:
The boundary conditions of the differential equations are as follows:
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 bywhere is the volumetric latent heat of soil, which is given bywhere is the latent heat of water, is the dry density of the soil, is the initial water content in the soil, and 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:where is the exponential integral function, which is given by
Substituting equations (51) and (52) into equation (49) yieldswhere is an undetermined constant that varies with time, which can be
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:where is the closure time of the frozen wall and l is the space between adjacent freezing pipes.
2.5.2. Flat-Panel Freezing Theory
Figure 7 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 , and the temperature of the soil in the freezing and unfrozen areas is and , respectively, both of which are functions of time t and coordinate x. The distance between the freezing front and the flat panel is , and the temperature of soil at the freezing front is the freezing point . The initial temperature of soil is .
For the flat-panel freezing problem, when , the differential equations of the frozen area can be written as
Similarly, when , the differential equation is
The initial conditions for the differential equations are as follows:
The boundary conditions of the differential equations are as follows:
In addition to the above boundary conditions, the heat balance equation at the phase transition boundary can be expressed as
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 follows :
The temperature field distributions in the unfrozen area can be expressed as follows :where is Gaussian integral function, which is given by
Substituting equations (64) and (65) into equation (63) yieldswhere is a constant to be determined, and the following condition exists:
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:where is the radius of the freezing pipe circle, is the radius of the outer freezing front, and is the radius of the inner freezing front.
The thickness of the frozen wall is obtained by
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 is equal to the freezing tube wall temperature .
2.6. Solution for the Outer Radius of the Expanded Area
The frost heave rate in engineering is generally expressed aswhere is the frost heave ratio of the soil, is the longitudinal frost heave of the sample, and 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
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.where is the frost heave ratio of soil without load, is the load, kPa, and is empirical constant, kPa−1.
Substituting equation (74) into equation (73) yields
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. Figure 8 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].
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 Figure 8.
For the convenience of calculation, the simplified thermal and physical parameters of the tunnel are shown in Tables 1 and 2.
The tangent value of the primary influence angle of the tunnel is, and the convergence radius can be obtained by
Therefore, and are taken as the calculation area. In accordance with equation (54), the thermal diffusivity of frozen and unfrozen soils may be obtained by
In accordance with equation (50), the volumetric latent heat of soil is
The axis x 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.
Figure 9 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.
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 10–14, respectively.
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 , 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.
As shown in Figures 10–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 Figures 20 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.
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.Nomenclature
| Ue: | The horizontal displacement |
| We: | The vertical displacement |
| dε: | The length of the microelement |
| dξ: | The width of the microelement |
| dη: | The thickness of the microelement |
| r(z): | Convergence radius |
| ε(eR, eθ, ez): | The volumetric strain of the microelement |
| XYZ: | Global coordinate system |
| RθZ: | Cylindrical coordinate system |
| εξη: | Cartesian coordinate system |
| rφζ: | Cylindrical coordinate system |
| β: | The main influence angle |
| Ψ: | The contract area |
| Ω: | The excavation area |
| Δ: | The deformation area |
| α: | Rotation angle |
| γ: | Rotation angle |
| l: | The length of the freezing pipe |
| h: | The depth between the tunnel center to the ground surface |
| r: | Radial coordinate |
| R(t): | The freezing radius |
| RΔ(t): | The frost heave radius |
| Rd: | The radius of the freezing pipe layout circle |
| Tf: | The temperature of the soil in the frozen area |
| Tu: | The temperature of the soil in the unfrozen area |
| T0: | The initial temperature of soil |
| Td: | The temperature of the freezing front |
| Tc: | The temperature of the freezing pipe wall |
| T′c: | The flat-pane cooling temperature |
| αf: | The thermal diffusivity of frozen soils |
| αu: | The thermal diffusivity of unfrozen soils |
| kf: | The thermal conductivity of frozen soils |
| ku: | The thermal conductivity of unfrozen soils |
| cf: | The specific heat of frozen soils |
| cu: | The specific heat of unfrozen soils |
| ρf: | The density of frozen soils |
| ρu: | The density of unfrozen soils |
| ρd: | The dry density of the soil |
| : | The initial water content in the soil |
| : | The unfrozen water content in the frozen soil |
| : | The latent heat of water |
| L: | The volumetric latent heat of soil |
| A: | An undetermined constant that varies with time |
| B: | A constant to be determined |
| P: | The load (kPa) |
| b: | Empirical constant (kPa−1) |
| εf: | Frost heave ratio of soil |
| εf0: | Frost heave ratio of soil without load |
| h: | The original height of the sample |
| Δhf: | The longitudinal frost heave of the sample. |
Data Availability
The datasets generated and analyzed during the current study are available from the corresponding author upon reasonable request.
Conflicts of Interest
The authors declare that there are no conflicts of interest regarding the publication of this paper.
Acknowledgments
This research was supported by the National Natural Science Foundation of China (Grant no.51778004), Research Activities Fund Project for Reserve Candidate of Academic and Technical Leaders of Anhui Province, China (Grant no. 2018H170), and Academic Funding for Top-notch Talents in University Disciplines (Majors) of Anhui Province, China (Grant no. gxbjZD10).