Temperature Field and Stability Analysis of the Frozen Wall Based on the Actual Position of Freezing Holes

Abstract

Taking the Qingdong Mine as the research object, combined with field measurement data, numerical simulation and theoretical analysis are used to examine the temperature field and stability of the frozen wall in the mine, respectively. The results show that during the active freezing period, under the same freezing time, the average temperature of the effective frozen wall of the fine sand layer is 0.2–1.0 and 0.5–2.5 °C lower than that of the sandy clay layer and clay layer, respectively. The effective frozen wall thickness of the fine sand layer is 0.04–0.17 and 0.17–0.33 m larger than that of the sandy clay layer and clay layer, respectively. The soil cooling between the two circles of freezing holes is the fastest. Due to the deflection of the freezing holes, the interface temperature field is asymmetrical. For deep clay with a depth of 200–250 m, it is most economical and reasonable for the brine temperature in the active freezing period to be −25 and −30 °C. At the designed brine temperature for cooling, during the excavation of the control layer of the topsoil layer (−216 m sandy clay), the side-wall temperature, average temperature, and thickness of the frozen wall meet the design requirements. The ultimate bearing capacity of the frozen wall is 3.20 MPa. When the well is empty for 30 h after excavation, the maximum radial displacement is 26.85 mm, so the frozen wall strength and stability are in a safe state. Overall, the findings of this study can serve as a useful reference for similar freezing projects.

1. Introduction

Artificial ground freezing (AGF) is a soil stabilization technique, which is used to freeze the water, rock, and soil in a formation by circulating a low-temperature refrigerant (usually low-temperature brine or liquid nitrogen) in the freezing tubes embedded in the rock and soil layer. The frozen wall with certain strength and low permeability can resist the earth pressure of groundwater, isolate the groundwater, and ensure the safe implementation of civil engineering projects [1]. This method is environmentally friendly and can be adapted for complex engineering conditions [2]. In the 1880s, this method was proposed for soft soil tunneling, mining, and excavation of underground regions. Since the 1950s, it has been widely used in the United States, Australia, China, and other countries [3].

In recent years, owing to its unique advantages, the AGF method has been increasingly used in the excavation of tunnels in regions rich in water and soft soil [4,5,6,7,8], deep foundation pit excavation [9,10], and other projects. Because coal mine shafts generally have a large depth, the topsoil is thick, and the geological and hydrogeological conditions of the topsoil are complex, so the freezing method is also widely used in coal mine shafts and has received considerable research attention. Hu et al. [11,12,13,14] derived the analytical solution of the steady-state temperature field under different arrangements of freezing pipes, which greatly promoted the theoretical application of the freezing method. Jiao et al. [15] analyzed the development law of the temperature field of frozen wall in the deep alluvial multi-circle holes and proposed a new numerical model for the freezing temperature field. The feasibility of the numerical model was verified through field measurement. Based on the actual engineering conditions of the Banji coal mine, Wang et al. [16] studied the application of the differential control freezing technology in the repair of damaged wellbore in deep alluvium by using numerical simulation from the perspective of temperature field development and distribution law. To understand the temporal and spatial evolution law of the frozen wall temperature field in the deep alluvial layer, it is essential to examine the characteristic surface temperature, average temperature, and thickness of the frozen wall of the multi-turn tubes in the deep alluvial layer [17,18,19]. In the construction of coal mine shafts using AGF method, if the groundwater velocity is greater than 5 m/d, the influence of the groundwater velocity on the freezing temperature field should be considered in the construction design of the freezing method and the analysis of the temperature field [20,21,22].

Wang et al. [23] studied the development law of the freezing temperature field of an expansive clay layer in the Yangcun Mine by using a large-scale physical model which can guide engineering applications. Based on the study of the freezing temperature field, the frozen wall is generally regarded as a temperature-dependent material, and its mechanical properties have been examined. Some researchers considered the temperature field of the frozen wall to be homogeneous and analyzed the elastoplasticity of the frozen wall to obtain the distribution of the displacement field and stress field under the action of external loads [24,25,26]. However, some scholars regarded the frozen wall as a non-homogeneous functionally graded material along the radial direction, simplified the radial temperature field of the frozen wall, and conducted an elastoplastic analysis of the frozen wall considering the excavation unloading effect to obtain the distribution law of stress field and displacement [27,28,29,30].

When the AGF method is used to excavate the shaft, due to the complexity of the deep soil geological environment and the diversity of soil properties, it is necessary to study the freezing temperature field of the deep soil. Therefore, the author takes Qingdong Mine as the research object, based on the actual location of the frozen holes in the field, combined with the soil thermophysical test parameters measured in the laboratory test, the numerical calculation software was used to study the deep soil freezing temperature field in the coal mine shaft. Based on the analysis of the temperature field of the frozen wall, the ultimate bearing capacity and creep deformation of the frozen wall at the topsoil control layer are studied, which provides a safety guarantee for the shaft excavation. The results provide an important reference for the design and construction of deep soil frozen walls.

2. Project Background

2.1. Shaft Geology

The Dongfeng well of Qingdong Mine is a newly built wellbore. The alluvial layer of the wellbore is buried at a depth of 234.55 m. The topsoil of the wellbore is Quaternary and Neogene from top to bottom. The Quaternary and Neogene clay layers are thick. They account for a large proportion and have a large amount of swelling, especially in the lower Neogene clay and calcareous clay layers, where semi-consolidated clays are the mainstay, with low water content and local disintegration when exposed to water. The freezing temperature of Neogene clay strata is low, between −1.82 and −2.17 °C, and the deep soil layer is more difficult to be frozen. The maximum thickness of the single layer of clay layer is 42.1 m, and the hole depth is 171.15–213.25 m. The freezing method is adopted for the excavation of the shaft.

2.2. Design of Freezing Scheme

According to the bidding documents, the freezing depth is determined to be 305 m. Based on the well wall structure diagram and the geological data of the well inspection holes provided in the bidding documents, after the optimization of the plan, the comprehensive freezing construction scheme is adopted in the Dongfeng Well, where the main freezing holes, auxiliary freezing holes, and prevent piece holes are operated together with the three circles of holes. The differential freezing scheme is employed in the main freezing holes, while the “plum blossom pile” arrangement is adopted for the auxiliary and prevent piece holes. To monitor the real-time development of the freezing temperature field, temperature measurement holes (C1–C3) and hydrological holes (S1–S2) are set in the freezing area. Among them, the S1 hole (hole depth: 40 m) is designed with four flower tube parts at a depth of 6.5–11.5, 25.5–30.5, 31.5–33, and 35.0–39.0 m. A partition at a vertical depth of 25 m is established, and the casing is lowered to obtain the change in the water level layer by layer (the casing size is Φ38 × 3.5 mm). For S2 hole (depth: 125 m), a flower tube at a position of 118.0–123.0 m is designed. The parameters of the freezing holes are listed in

Table 1. Parameters of the freezing holes.

Freezing Holes	Circle Diameter (m)	Hole Spacing (m)	Depth (m)	Freezing Method
Main freezing holes	13.2	1.296	305/255	Differential freezing
Auxiliary freezing holes	10.1	2.115	240	Full depth freezing
Prevent piece holes	9.5	1.990	150	Full depth freezing

, and their layout and cross-section are shown in Figure 1. Based on the depth of the wellbore alluvial layer and the control formation (−216 m sandy clay), combined with the relevant specification [26], the brine temperature during the active freezing period is determined to be −25 to −30 °C. The average design temperature of the effective frozen wall of the Neogene control horizon (−216.45 m sandy clay) is −12 °C, and the design thickness is calculated with Domke formula [31].

$$E=R\times \left[0.29\times \left(\frac{P}{R}\right)+2.3\times {\left(\frac{P}{K}\right)}^2\right]$$

(1)

where E is calculated frozen wall thickness, m; R is shaft radius, R = 4.053 m; P is calculated horizon pressure, MPa; P = 0.013H = 2.81 MPa; K is calculated strength of frozen soil, K = 4.93 MPa.

After calculation, E = 3.70 m. According to the calculation results and combined with the construction unit’s experience in freezing construction of vertical shafts with similar geological conditions, the thickness of the frozen wall E = 4.2 m is determined.

3. Analysis of Temperature Field of Frozen Wall

3.1. Mathematical Model of Freezing Temperature Field

The determination of freezing temperature field is a transient heat conduction problem with internal heat source and phase transition interface movement. The differential equations governing the two-dimensional (2D) freezing temperature field in the cartesian coordinate system are shown as follows:

In the unfrozen areas,

$$\rho C_u\frac{\partial T_u}{\partial t}=k_u\left(\frac{{\partial }^2{T}_u}{\partial x^2}+\frac{{\partial }^2{T}_u}{\partial y^2}\right)$$

(2)

In the frozen areas,

$$\rho C_f\frac{\partial T_f}{\partial t}=k_f\left(\frac{{\partial }^2{T}_f}{\partial x^2}+\frac{{\partial }^2{T}_f}{\partial y^2}\right)$$

(3)

At the phase transition interface,

$$\frac{k_f}{\cos \phi +\sin \phi }\left(\frac{\partial T_f}{\partial x}+\frac{\partial T_f}{\partial y}\right)-\frac{k_u}{\cos \phi +\sin \phi }\left(\frac{\partial T_u}{\partial x}+\frac{\partial T_u}{\partial y}\right)=\psi \frac{d{\xi }_N(x_N,y_N)}{dt}$$

(4)

$$T_f=T_u$$

(5)

where ρ is the density of soil, kg/m³; T_u and T_f are the temperatures of the unfrozen and frozen regions, respectively, °C; C_u and C_f are the specific heat capacities of unfrozen soil and frozen soil, respectively, J/(kg·K); k_u and k_f are the thermal conductivities of unfrozen soil and frozen soil, respectively, W/(m·K); x, y is plane position coordinates, m; $\phi$ is the circumferential coordinate with the freezing hole center as the coordinate origin, rad; $\left(x_N,y_N\right)$ is the coordinate of the frozen front, m; ${\xi }_N$ is the position of the frozen front, m; $\psi$ is the latent heat of phase transformation per unit volume of rock and soil mass, J/m³; t is the time, s.

The initial condition of soil temperature field is

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

(6)

where T₀ is the initial temperature of soil, °C.

The boundary conditions of soil temperature field are as follows:

$${\left.T\right|}_{\left(x_p,y_p\right)}=T_c\left(t\right)$$

(7)

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

(8)

where(x_p, y_p) is the freezing holes boundary coordinates, m; T_c(t) is the real-time temperature of the freezing holes, °C.

Equations (2)–(8) constitute the definite solution of the freezing temperature field in 2D cartesian coordinate system.

3.2. Numerical Calculation Model of Freezing Temperature Field

3.2.1. Basic Assumptions

• Unfrozen soil and frozen soil are homogeneous isotropic bodies.
• The influence of groundwater seepage on the temperature field is ignored.
• The vertical heat transfer of shaft formation and freezing holes is ignored, and the determination of freezing temperature field is simplified to a 2D plane problem.
• The measured brine temperature is directly applied to the boundary of the freezing holes as a temperature load.
• The latent heat of phase transition of the soil is released in a certain temperature range near the freezing temperature during freezing.

3.2.2. Numerical Calculation Model

A circular domain with a radius of 20 m is taken as the solution area of the model, and a numerical calculation model is established based on the actual positions of hole formation on site. The mesh division and geometric size of the model are shown in Figure 2.

3.2.3. Initial Value and Boundary Value Conditions of the Model

According to the in situ measured initial temperature of the formation before freezing, the initial temperature of −120.5 m fine sand, −190 m clay, and −216 m sandy clay are 21, 22.2, and 23.1 °C, respectively. The boundary temperature of freezing holes is the brine temperature measured on site, as shown in Figure 3.

3.2.4. Thermal Physical Parameters of Soil

Based on the measured data of on-site temperature measurement points, the parameters in the numerical calculation model are adjusted using the soil thermophysical parameters measured in the laboratory. The parameters of the numerical calculation model are shown in

Table 2. Thermophysical parameters of soil.

Depth (m)	Soil Property	Density (kg·m−3)	Thermal Conductivity (W·m−1·K−1)		Specific Heat Capacity (J·kg−1·K−1)		Latent Heat of Phase Transformation (kJ·kg−1)
			Frozen	Unfrozen	Frozen	Unfrozen
120.5	Fine sand	2003	1.735	1.412	1024	1170	28.6
190	Clay	1987	1.598	1.318	1199	1298	34.4
216	Sandy clay	1995	1.706	1.391	1080	1198	30.6

.

3.3. Numerical Results and Discussion

3.3.1. Model Feasibility Verification

According to the measured data in Figure 4a,b, the soil temperature at the C1 measurement point reaches the freezing point of soil when frozen for 54 d, while that at the C3 measurement point reaches the freezing point of soil when frozen for 30 d. This is because the measurement point inside the freezing hole ring diameter is affected by the superposition effect of cooling capacity during multi-hole freezing, so the soil temperature inside the diameter drops faster than that outside the diameter. To verify the feasibility of the numerical model, the measured temperature points of different layers and the measured side-wall temperature during excavation are compared with the corresponding numerical results. As shown in Figure 4 and Figure 5, the errors between the measured values in different formations and the calculated values are approximately 1 °C, which means that the calculated values are highly consistent with the measured values. This indicates that the numerical model is feasible and reliable for examining the distribution and development of freezing temperature field.

3.3.2. Average Temperature and Thickness of Effective Frozen Wall

The effective frozen wall represents the original frozen wall in each layer of this study minus part of frozen soil removed by shaft excavation and the retained frozen wall. Due to the random deflection of the freezing holes, the thickness of the frozen wall in different directions is not uniform, so the thickness in four different directions is taken in this study, and the average thickness of the four directions is taken as the research index of the frozen wall thickness. The average temperature of the frozen wall is the weighted average of the temperature in its area. By extracting the effective frozen wall thickness and average temperature data from the numerical model, the variation in the effective frozen wall thickness and average temperature with time is obtained, as shown in Figure 6.

(1)

Average temperature of frozen wall

It can be seen from Figure 6 that the average temperature of the frozen wall can be roughly divided into five stages. The first stage is from 40 to 60 d, during which the average temperature of the frozen wall decreases rapidly with the freezing time. The second stage is from 60 to 80 d; during this period, the average temperature of the frozen wall decreases slightly with the freezing time and generally remains stable. The average temperature of the frozen wall increases with the freezing time from 80 to 96 d, and it decreases with the freezing time from 96 to 109 d. It again increases with the freezing time from 109 to 120 d. The freezing process can be divided into active freezing period and maintenance freezing period. After the on-site freezing begins, the freezing unit pumps the continuously cooled brine through the circulation system to the freezing pipe embedded in the soil. By circulating the cold brine to cool the formation, the water in the formation is frozen, and a closed frozen wall is formed around the working face to resist the water and soil pressure from the formation. According to the brine temperature measured on site, after freezing for 20 to 30 d, the brine temperature in the main freezing holes and auxiliary holes reaches the designed low temperature level. Further, the brine temperature remains basically stable and then enters the active freezing stage (80 d before freezing).

According to Figure 3, after freezing for 80 d, the brine temperature circulating in the auxiliary holes is in a state of constant fluctuation, so 80–120 d is the maintenance freezing stage. In this stage, the average temperature of the frozen wall rises and falls with the temperature of brine circulating in the frozen hole and shows the same variation trend. Therefore, the average temperature of the frozen wall rises (drops) as the temperature of brine circulating in the freezing holes increases (decreases). According to the above analysis, in the active freezing stage, the average temperature of the frozen wall decreases rapidly after overlap and then remains stable after reaching a certain level. In the maintenance freezing stage, the average temperature of the frozen wall changes with the temperature of the circulating brine and the two temperatures are synchronized. For the same freezing time, the average temperature of fine sand is the lowest, followed by that of sandy clay, and the average temperature of clay is the highest. In the active freezing period, the average temperature of the fine sand frozen wall is 0.2–1.0 and 0.5–2.5 °C lower than that of the sandy clay layer and clay layer, respectively. The difference decreases with the freezing time. During the maintenance freezing period, the average temperature of the fine sand frozen wall remains basically the same as that of the sandy clay layer, but it is 0.3–0.5 °C lower than that of clay frozen wall. Therefore, under the same freezing time and freezing scheme, the freezing effect of fine sand is the best, while that of sandy clay layer is the worst. For the control layer, i.e., the −216 m sandy clay layer, the design requirement of the average temperature of the frozen wall is −12 °C, so the average frozen wall temperature in the control layer meets the requirement.

(2)

Frozen wall thickness

It can be seen from Figure 6 that the effective frozen wall thickness is positively correlated with the time; the longer the freezing time, the greater the thickness of the frozen wall. The thickness of the frozen wall increases rapidly during the initial stage of the frozen wall overlap. This is mainly because the frozen wall expands on both sides from outside and within the freezing hole ring diameter during the early stage of overlap. When the inner frozen wall extends to the excavation path, the thickness of the effective frozen wall is only determined by the expansion rate of frozen soil outside the freezing hole ring, and the lateral frozen soil directly contacts the soil in the normal temperature area, so the frozen wall expands slowly and the thickness of the frozen wall gradually increases with time. Under the same cooling condition, the thickness of the frozen wall in the fine sand layer is the largest, followed by that of the sandy clay layer, and the frozen wall in the clay layer has the smallest thickness. This is mainly determined by the intrinsic soil characteristics. Specifically, the thermal conductivity of fine sand > thermal conductivity of sandy clay > thermal conductivity of clay, and the higher the thermal conductivity, the faster the heat transfer, so the thickness of the frozen wall increases faster. Under the same freezing time, the effective frozen wall thickness of the fine sand layer is 0.04–0.17 and 0.17–0.33 m larger than that of the sandy clay layer and clay layer, respectively.

3.3.3. Spatial-Temporal Variation of Temperature Field at the Main Surface and Interface

The temperature data on the main surface and interface at different times in the temperature field were extracted, and the obtained radial temperature distribution of the frozen wall at the main surface and interface in different soil layers is shown in Figure 7. According to the freezing holes, the temperature field of frozen wall can be divided into three zones: zone I, zone II, and zone III. Under the same freezing time, the cooling in zone II is the fastest, followed by that in zone I and zone III. This is because during the freezing process, there is no external heat supply in zone II, and this area is also affected by the superposition of cooling capacity of auxiliary holes and main freezing holes. Although zone I is less affected by the superposition of cooling capacity of double-loop pipes during the freezing process, there is no external source to supplement heat in this zone. The freezing process of zone III is not affected by the superposition of the cooling capacity of the double-ring freezing holes, and the surrounding soil with constant ground temperature can be used as an external source to supplement heat for the soil in the zone. During the freezing process, the closer the region to the freezing hole, the faster the cooling. At the same layer and under the same freezing time, the temperature at the main surface of the frozen wall is obviously lower than that at the interface. When frozen for 70 d, the minimum temperature of the main surface is approximately 7 °C lower than that of the interface. In the same layer, during the later stage of freezing, the lowest temperature on the main surface is always observed at the frozen hole, while on the interface, due to the random deflection of the freezing holes, the lowest temperature is located at a certain position between the two circles of freezing holes. For different soil layers, regardless of the main surface or interface, under the same freezing time, the cooling of fine sand is the fastest, followed by that of sandy clay, and clay layer has the slowest cooling. When frozen for 70 d, the minimum temperature of the fine sand layer is 0.5 and 1 °C lower than that of the sandy clay layer and clay layer, respectively. In the case of deflection of the freezing holes, the temperature field distribution of the main surface is symmetric; especially, the soil temperature field in the II zone has high symmetry. However, the temperature field at the interface is not symmetrical.

3.3.4. Influence of Brine Temperature on Temperature Field of Frozen Wall

The temperature of circulating cooling agent (brine) governs the formation of frozen wall. Here, −216 m sandy clay layer control is taken as the research object, and based on the actual position of the freezing holes, the active freezing period under different brine temperatures as well as the variations in the characteristic parameters of the temperature field of the frozen wall (effective frozen wall average temperature, effective thickness of frozen wall) with time are examined. The brine temperature in the active freezing period is −20, −25, −30, and −35 °C. It is assumed that the on-site circulating brine temperature for the first 19 d is the same as that in the actual project, and the brine temperature of the freezing holes remains unchanged after the 20th day of brine circulation. As shown in Figure 8 and Figure 9, the variation in the brine temperature during the active freezing period affects the average temperature and thickness. Under the same freezing time, the lower the brine temperature, the lower the effective frozen wall temperature and the greater the thickness. When frozen for 130 d, the brine temperature decreases from −20 and −35 °C, and the average temperature of effective frozen wall decreases from −12.16 to −21.09 °C, corresponding to a decrease of 73.4%. Under freezing for 130 d, the effective frozen wall thickness increases from 4.05 m to 4.65 m, which corresponds to an increase of 14.8%. For the sandy clay layer at the control horizon of topsoil (−216 m), during excavation (frozen for 130 d), the designed average temperature of effective frozen wall is −12 °C, and the designed thickness is 4.2 m. If the brine temperature in the active freezing period is −20 °C, the average temperature of the effective frozen wall just reaches −12 °C during the excavation of this layer, but the thickness cannot reach the design requirement of 4.2 m. Therefore, the salt water temperature in the active freezing period is −20 °C, which cannot meet the design requirements. If the brine temperature during the active freezing period is set to −25 °C or below, as shown in Figure 8 and Figure 9, the average temperature and thickness of effective frozen wall during the excavation of the control layer meet the design requirements. Therefore, the brine temperature during the active freezing period is set to −25 °C or below. It can be seen from Figure 8 and Figure 9 that when the brine temperature is −25 °C and below, the average temperature of effective frozen wall reaches −15 to −21 °C and the thickness is 4.29–4.65 m when frozen for 130 d; under the premise of ensuring safety, this can avoid increasing the difficulty in excavation due to the entry of excessive frozen soil in the excavation section as well as save the freezing costs. Therefore, after entering the maintenance freezing period, the brine temperature in the freezing hole should be adjusted in real-time while maintaining the average temperature and thickness of the effective frozen wall. According to Figure 3, the brine temperature during the active freezing period is −25 to −27 °C. The brine temperature of the main freezing holes remains unchanged after freezing for 80 days, and the brine temperature of auxiliary holes is controlled, which ensures the thickness of the effective frozen wall, reduces the expansion of frozen soil into the excavation section, decreases the difficulty in excavation, and saves freezing cost. Therefore, the theoretical analysis is consistent with the engineering practice.

4. Analysis of Viscoelastic-Plastic Stability of Frozen Wall

After shaft excavation, the frozen wall, as a temporary support structure, resists the soil-water pressure in the stratum, and its safety is necessary for the personnel security and working face construction progress. The frozen wall after excavation should have sufficient bearing capacity to resist the external load of the frozen wall (which can be considered as permanent horizontal ground pressure acting on the outer frozen wall). At the same time, as the frozen wall is a rheological material, to ensure its stability during the period without support after excavation, its creep deformation should meet certain requirements. The stress–strain relationship and deformation of the frozen wall are related to the physical and mechanical parameters of the material, and these parameters are directly determined by the negative temperature level of the frozen wall. Therefore, in this study, the security of the frozen wall is evaluated in terms of the temperature field and stability of the frozen wall. The frozen wall of the control horizon of the topsoil of Qingdong Mine −216 m sandy clay horizon after excavation is selected for analysis. When frozen for 131 d, the soil of this layer is excavated, the excavation radius a = 4.053 m, the frozen wall radius b = 8.403 m, and the horizontal ground pressure borne by the frozen wall is P, MPa. The horizontal ground pressure P = 0.013 H [31], where H is the buried depth, m. The horizontal pressure P = 2.81 MPa in the control horizon. The mechanical model of the frozen wall is shown in Figure 10.

4.1. Evaluation of the Temperature Field of Frozen Wall

According to the relevant specification [31,32] and the freezing scheme design, the effective frozen wall thickness of the control layer of the wellbore surface soil layer is 4.2 m, the average temperature is −12 °C, and the inner edge temperature of driving section in the −150 to −250 m clay layer is −4 to −8 °C. Based on the linear interpolation method, the edge temperature of the clay horizon with control horizon of −216 m is −6.6 °C. According to the numerical results, after freezing time of 131 d, the excavation begins, and the effective frozen wall thickness is 4.35 m. The average temperature of effective frozen wall is −14.14 °C. The wellbore temperature is below −7 °C, as shown in Figure 11. To summarize, all the characteristic indexes of the frozen wall temperature field meet the design requirements of the freezing scheme.

4.2. Elastic-Plastic Analysis of Frozen Wall

4.2.1. Elastic Stress State of Frozen Wall

In the elastic state, the stress field of the frozen wall is expressed as follows [33]:

$$\left\{\begin{array}{cc}{\sigma }_r^e=\frac{b^{2}P}{b^2-a^2}\left(1-\frac{a^2}{r^2}\right)& \\ {\sigma }_{\theta }^e=\frac{b^{2}P}{b^2-a^2}\left(1+\frac{a^2}{r^2}\right)& \end{array}\right.$$

(9)

4.2.2. Elastic-Plastic Stress Analysis of Frozen Wall

When the external load of the frozen wall exceeds its ultimate elastic load P_e, the frozen wall enters the elastic-plastic state. The frozen wall begins to yield from the inner wall. With the increase in the external load of the frozen wall, the plastic zone continues to expand toward the outer wall, and the plastic zone radius of the frozen wall is assumed to be r_c. Further, it is assumed that the frozen wall satisfies the Mohr–Coulomb criterion after entering the plastic state, i.e.,

$${\sigma }_1=\frac{1+\sin \phi }{1-\sin \phi }{\sigma }_3+\frac{2c\cos \phi }{1-\sin \phi }$$

(10)

where c is cohesive force of frozen soil, MPa; φ is the angle of internal friction in the frozen soil, °. For plane strain state, ${\sigma }_z=\frac{1}{2}\left({\sigma }_r+{\sigma }_{\theta }\right)$, and it is assumed that the compression is positive and tension is negative. Therefore, ${\sigma }_{\theta }>{\sigma }_z>{\sigma }_r$, i.e.,${\sigma }_1={\sigma }_{\theta }$, ${\sigma }_2={\sigma }_z$, and ${\sigma }_3={\sigma }_r$. In the elastic-plastic state, ${\sigma }_1={\sigma }_{\theta }^p$ and ${\sigma }_3={\sigma }_r^p$. Substituting these expressions in Equation (10), we obtain

$${\sigma }_{\theta }^p-{\sigma }_r^p=\frac{2\sin \phi }{1-\sin \phi }{\sigma }_r^p+\frac{2c\cos \phi }{1-\sin \phi }$$

(11)

For axisymmetric plane strain structure, the equilibrium equation is expressed as follows:

$$\frac{d{\sigma }_r^p}{dr}+\frac{{\sigma }_r^p-{\sigma }_{\theta }^p}{r}=0$$

(12)

Substituting Equation (11) into Equation (12), we get

$$\frac{d{\sigma }_r^p}{dr}-\frac{1}{r}\frac{2\sin \phi }{1-\sin \phi }=\frac{1}{r}\frac{2c\cos \phi }{1-\sin \phi }$$

(13)

Using the boundary condition: ${\sigma }_r^p\left|{}_{r=a}\right.=0$, it can be obtained that

$${\sigma }_r^p=\frac{c\cos \phi }{\sin \phi }\left[a^{-\frac{2\sin \phi }{1-\sin \phi }}{r}^{\frac{2\sin \phi }{1-\sin \phi }}-1\right]$$

(14)

Substituting Equation (14) into Equation (11), we get

$$\begin{array}{cc}\left\{\begin{array}{c}{\sigma }_r^p=\frac{c\cos \phi }{\sin \phi }\left[a^{-\frac{2\sin \phi }{1-\sin \phi }}{r}^{\frac{2\sin \phi }{1-\sin \phi }}-1\right]\\ {\sigma }_{\theta }^p=\frac{c\left(1+\sin \phi \right)\cos \phi }{\sin \phi \left(1-\sin \phi \right)}\left[a^{-\frac{2\sin \phi }{1-\sin \phi }}{r}^{\frac{2\sin \phi }{1-\sin \phi }}-1\right]+\frac{2c\cos \phi }{1-\sin \phi }\end{array}\right.& r\in \left[a,r_c\right]\hfill \end{array}$$

(15)

When $r\in (r_c,b]$, the frozen wall can be regarded as a thick-walled elastic cylinder, whose inner and outer walls are subjected to ${\sigma }_{r_c}$ and P, respectively, and the analytical expression of its stress is [24]

$$\begin{array}{cc}\left\{\begin{array}{c}{\sigma }_r^e=\frac{b^{2}P-r_c^2{\sigma }_{r_c}}{b^2-r_c^2}-\frac{b^2{r}_c^2\left(P-{\sigma }_{r_c}\right)}{r^2\left(b^2-r_c^2\right)}\\ {\sigma }_{\theta }^e=\frac{b^{2}P-r_c^2{\sigma }_{r_c}}{b^2-r_c^2}+\frac{b^2{r}_c^2\left(P-{\sigma }_{r_c}\right)}{r^2\left(b^2-r_c^2\right)}\end{array}\right.& \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em} r\in (r_c,b]\hfill \end{array}$$

(16)

At the interface of the elastic-plastic zone of the frozen wall, when r = r_c,

$$\left\{\begin{array}{cc}{\sigma }_{r_c}^e={\sigma }_{r_c}^{}& \\ {\sigma }_{\theta }^e=\frac{2b^{2}P-\left(b^2+r_c^2\right){\sigma }_{r_c}^{}}{b^2-r_c^2}& \end{array}\right.$$

(17)

Equation (17) satisfies Equation (10), and ${\sigma }_1={\sigma }_{\theta }^e$, ${\sigma }_3={\sigma }_{r_c}^e$. Substituting Equation (17) into Equation (10), we get

$$\frac{b^2(P-{\sigma }_{r_c})}{b^2-r_c^2}=\frac{b^{2}P-r_c^2{\sigma }_{r_c}}{b^2-r_c^2}\sin \phi +c\cos \phi$$

(18)

According to the stress continuity at the interface of elastic-plastic zone of frozen wall,

$${\sigma }_{r_c}={\sigma }_{r_c}^e={\sigma }_{r_c}^p=\frac{c\cos \phi }{\sin \phi }\left[a^{-\frac{2\sin \phi }{1-\sin \phi }}{r}_c{}^{\frac{2\sin \phi }{1-\sin \phi }}-1\right]$$

(19)

Using Equations (18) and (19), the plastic zone radius r_c of the frozen wall under the action of external load P can be obtained.

When the inner wall of the frozen wall just enters the plastic zone, r_c = a. At this time, the frozen wall is in the elastic limit state, and the external load borne by the frozen wall is denoted as the elastic ultimate bearing capacity P_e. When r_c = a, ${\sigma }_{r_c}=0$. According to Equation (18), the elastic ultimate bearing capacity of frozen wall can be obtained as follows [24]:

$$P_e=\frac{c\cos \phi \left(b^2-a^2\right)}{b^2(1-\sin \phi )}$$

(20)

When the frozen wall is completely plastic, r_c = b. In this case, the external load borne by the frozen wall is called the plastic ultimate bearing capacity P_p. According to the first formula in Equation (15), the plastic ultimate bearing capacity of the frozen wall can be written as follows [24]:

$$P_p=\frac{c\cos \phi }{\sin \phi }\left[a^{\frac{-2\sin \phi }{1-\sin \phi }}{b}_{}^{\frac{2\sin \phi }{1-\sin \phi }}-1\right]$$

(21)

The excavation radius of wellbore a = 4.053 m, and the outer radius of frozen wall b = 8.403 m. According to the indoor triaxial shear test of frozen soil, when the frozen soil temperature is −14.14 °C, its cohesion c =1.39 MPa, and the internal friction angle $\phi$ = 13.18°. Substituting the above parameters into Equations (20) and (21), the elastic ultimate bearing capacity can be obtained as Pe = 1.35 MPa, and the plastic ultimate bearing capacity Pp = 3.20 MPa. The external load of the frozen wall is P = 2.81 MPa, which is greater than its elastic ultimate bearing capacity but less than its plastic ultimate bearing capacity. Therefore, the frozen wall is in the safe state of bearing capacity.

4.3. Viscoelastic Analysis of Frozen Wall

It is assumed that the frozen soil is incompressible [34], i.e., ${\epsilon }_r+{\epsilon }_{\theta }+{\epsilon }_z=0$, and the in-plane strain ${\epsilon }_z=0$. In this case, the geometric equation can be obtained as follows:

$$\frac{du}{dr}+\frac{u}{r}=0$$

(22)

The solution of the above equation is

$$u=\frac{k}{r}$$

(23)

where k is the integral constant.

The creep of frozen soil can be expressed by a unified rheological equation as follows [35]:

$$\epsilon =A{\sigma }^B{t}^C$$

(24)

In practice, the frozen wall is in a triaxial stress state. Therefore, Equation (24) can be extended to a three-dimensional (3D) state through equivalent stress and equivalent strain. For the 3D state, the constitutive equation of frozen soil creep is

$${\epsilon }_i=A{\sigma }_i^B{t}^C$$

(25)

where ${\epsilon }_i=\frac{\sqrt{2}}{3}\sqrt{{\left({\epsilon }_{\theta }-{\epsilon }_r\right)}^2+{\left({\epsilon }_r-{\epsilon }_z\right)}^2+{\left({\epsilon }_z-{\epsilon }_{\theta }\right)}^2}$, ${\sigma }_i=\frac{\sqrt{2}}{2}\sqrt{{\left({\sigma }_{\theta }^{ne}-{\sigma }_r^{ne}\right)}^2+{\left({\sigma }_r^{ne}-{\sigma }_z^{ne}\right)}^2+{\left({\sigma }_z^{ne}-{\sigma }_{\theta }^{ne}\right)}^2}$

Using the geometric equation and ${\sigma }_z=\frac{1}{2}\left({\sigma }_r^{ne}+{\sigma }_{\theta }^{ne}\right)$ [35], it can be obtained that ${\epsilon }_i=\frac{2}{\sqrt{3}}\frac{k}{r^2}$; ${\sigma }_i=\frac{\sqrt{3}}{2}\left({\sigma }_{\theta }^{ne}-{\sigma }_r^{ne}\right)$.

Substituting the above expressions of equivalent strain ${\epsilon }_i$ and equivalent stress ${\sigma }_i$ into Equation (25) and simplifying it, we get

$${\sigma }_{\theta }^{ne}-{\sigma }_r^{ne}={\left(\frac{2}{\sqrt{3}}\right)}^{\frac{B+1}{B}}{\left(\frac{k}{A}\right)}^{\frac{1}{B}}{r}^{-\frac{2}{B}}{t}^{-\frac{C}{B}}$$

(26)

Substituting Equation (26) into Equation (12), we obtain

$$\frac{d{\sigma }_r^{ne}}{dr}={\left(\frac{2}{\sqrt{3}}\right)}^{\frac{B+1}{B}}{\left(\frac{k}{A}\right)}^{\frac{1}{B}}{r}^{-1-\frac{2}{B}}{t}^{-\frac{C}{B}}$$

(27)

The general solution of the above equation can be obtained as follows:

$${\sigma }_r^{ne}=-\frac{B}{2}{\left(\frac{2}{\sqrt{3}}\right)}^{\frac{B+1}{B}}{\left(\frac{k}{A}\right)}^{\frac{1}{B}}{r}^{-\frac{2}{B}}{t}^{-\frac{C}{B}}+d$$

(28)

where d is the integral constant, which can be obtained by using boundary conditions. In other words, when r = a, ${\sigma }_r^{ne}=0$, then

$$d=\frac{B}{2}{\left(\frac{2}{\sqrt{3}}\right)}^{\frac{B+1}{B}}{\left(\frac{k}{A}\right)}^{\frac{1}{B}}{a}^{-\frac{2}{B}}{t}^{-\frac{C}{B}}$$

(29)

Substituting Equation (29) into Equation (28), we obtain

$${\sigma }_r^{ne}=-\frac{B}{2}{\left(\frac{2}{\sqrt{3}}\right)}^{\frac{B+1}{B}}{\left(\frac{k}{A}\right)}^{\frac{1}{\mathrm{B}}}{t}^{-\frac{C}{B}}\left(r^{-\frac{2}{B}}-a^{-\frac{2}{B}}\right)$$

(30)

where the constant K can be used as the boundary condition: when r = b, ${\sigma }_r=P$.

$$k=\frac{\sqrt{3}}{2}A{t}^C{\left[\frac{-\sqrt{3}P}{B(b^{-\frac{2}{B}}-a^{-\frac{2}{B}})}\right]}^B$$

(31)

The radial creep displacement of the frozen wall can be obtained using Equations (23) and (31) as follows:

$$u=\frac{\sqrt{3}}{2}\frac{1}{r}A{t}^C{\left[\frac{-\sqrt{3}P}{B(b^{-\frac{2}{B}}-a^{-\frac{2}{B}})}\right]}^B$$

(32)

The creep test parameters can be determined according to the creep test of frozen soil under uniaxial compression: A = 0.00195, B = 1.591, and C = 0.317. Substituting r = a = 4.053 m, b = 8.403 m, P = 2.81 MPa, and the values of A, B, and C into Equation (30), the variation in the side-wall displacement with time can be obtained, as shown in Figure 12. It is evident that when the shaft wall is exposed for 30 h, its displacement is 26.85 mm. According to the Coal Mine Technical specification for vertical shaft Sinking by Freezing Method: MT/T1124-2011, in the thick clay layer with buried depth greater than 200 m, the exposure time of the shaft wall is limited to 30 h, and the maximum radial displacement of the shaft wall is not more than 50 mm. Therefore, the displacement of the frozen wall meets the requirements of the code and it is in a safe stable state.

5. Conclusions

Based on the actual position of freezing holes in the Qingdong Mine, combined with the real-time monitoring of the brine temperature of freezing holes, a numerical model was established for calculating the transient thermal temperature field considering latent heat of phase transformation. The calculated values of temperature measurement points and side-wall temperature were verified by the field measured data. The distribution and development of freezing temperature field of fine sand, sandy clay, and clay layers were studied by the proposed numerical method. At the same time, the stability of the frozen wall during the excavation of the control layer was analyzed theoretically. The main results of this study are summarized as follows:

• Under the same freezing time, the average effective freezing temperature of the fine sand layer was 0.2–1 and 0.5–2.5 °C lower than that of the sandy clay layer and clay layer, respectively. The effective frozen wall thickness of the fine sand layer was 0.04–0.17 and 0.17–0.33 m larger than that of the sandy clay layer and clay layer, respectively. The temperature field at the main surface and interface of the fine sand layer had the fastest decrease, while that of the sandy clay layer had the slowest decrease. Therefore, the fine sand layer had the best freezing effect, while the sandy clay layer had the worst freezing effect.
• Based on the main freezing holes and auxiliary holes, the freezing temperature field could be divided into three characteristic areas in the radial direction. Due to the superposition effect of the cooling capacity of double-ring freezing holes, the cooling of soil in the frozen zone II was the fastest, followed by that in the frozen zone I without external heat source, and the frozen zone III without superposition effect of the cooling capacity of freezing holes with external heat source had the slowest cooling. The deflection of the freezing holes had a minor influence on the symmetry of the temperature field II zone of the main surface, but it had a significant influence on the symmetry of the temperature field II zone of the interface. Therefore, the drilling process should be strictly controlled in the freezing construction to minimize the deflection of freezing holes.
• The brine temperature had a considerable influence on the temperature field of the frozen wall. At a high brine temperature, the average temperature and thickness of the effective frozen wall could not meet the design requirements. A lower brine temperature ensured the safety of shaft excavation, but greatly increased the cost of freezing construction. According to the theoretical analysis, the temperature of brine during the active freezing period of the project should be −25 to −30 °C, and the temperature of brine during the active freezing period was finally determined to be −25 to −27 °C based on the actual on site situation in the middle and late period of freezing during the maintenance freezing period. Consequently, to ensure an appropriate effective thickness and average temperature of frozen wall, the brine temperature of auxiliary holes should be controlled in real time to reduce the expansion of frozen soil into the excavation section, decrease the difficulty in excavation, and save the freezing cost.
• According to the actual freezing situation at the site, the safety of the frozen wall in the −216 m sandy clay control layer during excavation was evaluated. The results showed that during the excavation of this formation, the side-wall temperature was lower than −7 °C, which met the design requirement of −6.6 °C. The average effective frozen wall temperature was −14.14 °C, which also met the design requirement of −12 °C. The effective frozen wall thickness was 4.35 m, which met the design requirement of 4.2 m. The external load of frozen wall was 2.81 MPa, which was lower than its plastic ultimate bearing capacity (3.20 MPa). The maximum radial displacement of the frozen side wall at 30 h was 26.85 mm, which was less than the maximum limit of 50 mm stipulated in the code. Therefore, in this kind of freezing condition, the frozen wall of control layer is expected to be in a safe state during excavation.