Abstract
Taking into account moisture migration and heat change during the soil freezing process, as well as the influence of absolute porosity reduction on seepage during the freezing process, we construct a numerical model of hydrothermal coupling using laws of conservation of energy and mass. The model is verified by the results of large-scale laboratory tests. By applying the numerical calculation model to the formation of artificial shaft freezing temperature fields under the action of large-flow groundwater, we conclude that groundwater with flow rates of less than 5 m/d will not have a significant impact on the artificial freezing temperature field. The maximum flow rates that can be handled by single-row freezing pipes and double-row freezing pipes are 10 m/d and 20 m/d, respectively, during the process of freezing shaft sinking. By analyzing the variation of groundwater flow rate during freezing process, we find that the groundwater flow velocity can reach 5–7 times the initial flow velocity near the closure moment of the frozen wall. Finally, in light of the action characteristics of groundwater on the freezing temperature field, we make suggestions for optimal pipe and row spacing in freezing pipe arrangement.
1. Introduction
Artificial ground freezing is currently the main method for the construction of shafts in water-rich soft soil layers. This method involves the construction of one or more loops of freezing pipes around a shaft wellbore; through continuous heat exchange between the low-temperature refrigerant in the freezing pipe and the frozen soil, the circulation of fluids results in a frozen wall with a certain strength and good sealing performance, providing a stable environment for the construction of the wellbore [1–4]. However, the increase of the groundwater flow velocity in the soil layer when the traditional freezing pipe arrangement scheme is adopted can lead to increased closure times or even failure to close. This leads to large economic losses when soil excavation and the construction of the wellbore cannot be carried out as scheduled [5–7]. As a result, it is of great interest to study the formation law of the artificial shaft freezing temperature field under the action of groundwater with a large flow rate.
Our work builds on past research. Zhou et al. [8] studied the influence of seepage and pipe spacing on the closure time of saturated sand in the conventional brine freezing process and obtained the development law of temperature field in different regions using the double-pipe freezing orthogonal model test. Huang et al. [9] carried out a model test on a single freezing pipe under seepage conditions and found that the action of water flow resulted in different degrees of reduction in the frozen area upstream and downstream of the freezing pipe compared with nonseepage conditions. Lao et al. [10] obtained a mathematical model of a three-pipe freezing temperature field under the action of groundwater with different flow rates through model testing. Wang et al. [11] studied the effect of groundwater flow rates on the freezing temperature field under liquid nitrogen (−80°C) freezing conditions. The results showed that flowing water with a flow rate of 10 m/d or more had a significant influence on the liquid nitrogen freezing effect. Vitel et al. [12] designed a mathematical model of hydrothermal coupling that was completely consistent with thermodynamic principles. A three-dimensional freezing test under high-permeability flow conditions provided a key parameter basis for the study of hydrothermal coupling in phase transitions of porous media. Anagnostou et al. [13] summarized previous scholars’ models describing the formation of freezing temperature fields under the action of groundwater and designed new experimental devices on the basis of previous solutions to problems of freezing-cold volume loss. Large-scale model tests of artificial ground freezing under seepage conditions with flow rates of 0, 1, 1.5, 2.0, and 2.1 m/d were performed. The test results provide an important basis for the numerical calculation and engineering application of the freezing method.
The large similar model test has limitations. Most notably, it can only partially restore the development law of the freezing temperature field. However, numerical calculation can solve this problem. As a result, the numerical calculation has always been an important technical means to study the hydrothermal coupling problem.
Along these lines, Yang and Pi [14] established a mathematical model for the development of the frozen peak surface of a single frozen pipe under the action of groundwater using the theory of heat transport in porous media and Darcy’s law and analyzed variations in the temperature field and groundwater flow field during the freezing process. Gao et al. [15] and Liu et al. [16] used the finite element method to study the evolution of temperature fields in vertical and horizontal tunnel freezing under the action of groundwater. Vitel et al. [12, 17, 18] constructed a hydrothermal numerical model consistent with thermodynamics to simulate the artificial ground freezing of saturated nondeformable porous media under seepage conditions. This numerical model has been well verified by the results of three-dimensional ground freezing experiments under high seepage velocity conditions. Alzoubi et al. [19, 20] conceived and developed a controlled laboratory scale AGF experimental rig and developed a three-dimensional conjugate mathematical and numerical model of the bayonet freeze pipes and porous ground structure using enthalpy-porosity method. The model was further validated against global heat balance and local temperature distributions from experiments at various operating conditions. This study provides an important reference for design and optimization of industrial AGF system. Lin et al. [21] addressed the difficult closure of a frozen wall in a coal Mine shaft due to excessive seepage velocity in an aquifer when the aquifer is penetrated via the artificial freezing method. Huang et al. [22–24] developed a hydrothermal coupling model to simulate the influence of water flow on the freezing process by considering the water/ice phase transition and combined the model with the Nelder–Mead simplex method based on the COMSOL Multiphysics platform. The position of the freezing pipes around the circular tunnel was optimized. Ahmed et al. [25] used the “ant colony algorithm” to optimize the layout of the freezing pipes under the action of low-flow groundwater, thereby shortening the closure time of the frozen wall and making the strength of the frozen wall more uniform. This method led to new ideas for the optimal design of the layout of freezing pipes.
By studying moisture migration and heat changes in the soil during the freezing process and considering the influence of the absolute porosity reduction on seepage during the freezing process, we construct a numerical calculation model for hydrothermal coupling based on the laws of conservation of energy and mass. Then, we verify the model using large-scale laboratory tests. Finally, we calculate and analyze the formation law of the artificial frozen wall temperature field under the action of large-flow groundwater and propose an optimal freezing pipe arrangement.
2. The Coupled Hydrothermal Model
2.1. Temperature Field Equation
In the porous medium, a control body is taken and simultaneously subjected to a seepage field and a temperature field. According to the conservation of energy,where and are the heat supplied to the control body by an external heat source and the heat generated by the heat source inside the control body, respectively; , , and are the heat that needs to be consumed to change the temperature of the solid skeleton, ice, and water, respectively; and is the heat generated during the water-ice phase transition.
The heat from the external heat source flowing into the control unit per unit time iswhere is the heat flux density vector, which satisfies the Fourier heat transfer law
Regardless of the source influence inside the porous medium,
The amounts of heat required to change the temperature of the solid skeleton and ice per unit time are expressed as
Considering the convective heat transfer of water, the amount of heat required to change the temperature of water of per unit time is expressed aswhere is the porosity of the porous media; is the volumetric content of unfrozen water in the frozen soil; and is the relative velocity vector of the fluid; and are density and specific heat, respectively; and subscripts s, i, and l represent solid skeleton, ice, and water, respectively.
During the freezing process, when the temperature of the control body drops to the freezing point, the water in the control body is gradually converted into ice. During this process, the temperature remains just above freezing point, and no significant temperature change occurs. This is the water-ice phase transformation. Water will release latent heat continuously in the stage.where is the content of liquid water and is the latent heat of the water-ice phase transition.
Substituting equations (2)–(8) in equation (1), we obtain the differential equation of heat conduction for low-temperature porous media:
Since the control body selection is arbitrary, convert equation (1) towhere is the equivalent heat capacity.
During the phase-change process, the proportion of water to ice in the control body changes continuously. We assume that the content of unfrozen water is only a function of temperature T. The water-ice phase transformation process only occurs in the interval [0, −1]°C. To reasonably describe this process, we introduce the Heaviside function:where is the temperature at a certain location, is the phase transition temperature (and generally ), and is the radius of the phase-change interval.
Tests show that when the geotechnical temperature drops below the freezing temperature, a small amount of unfrozen water remains inside [7]. The variation of the unfrozen water content in the space of the control body with temperature can be expressed as
The water content of the whole control body can be expressed as
Furthermore, the ice content of the whole control body can be obtained:
The content of solid particles is as follows:
The distribution of water, ice, and solid particles in the control body during the freezing process is shown in Figure 1.
According to the volume-averaging method, the equivalent heat capacity of the soil can be obtained as follows:
The thermal conductivity K refers to the heat passing through the unit area soil per unit time under the unit gradient. The thermal conductivity inside the porous media by the volume weighted average method [7, 22, 23] is
Again, subscripts s, i, and l represent the solid skeleton, ice, and water, respectively.
2.2. Water Flow Field Equation
The surface element is taken on the outer surface of the control body , and the fluid mass passing through the surface of per unit time is expressed aswhere is the flow velocity through the section and is the outer normal direction of the section unit .
In the control body , for any volume unit selected, a change in the ratio of water to ice causes the density of the medium to change with time, causing the mass of the entire control body to change. The process can be described as
Regardless of the internal source of the porous medium, using conservation of mass, we can obtain
Water seepage in porous media is assumed to satisfy Darcy’s law. And the rate of water seepage in porous media is [21–23]where is the intrinsic permeability matrix, in which is the permeability of saturated porous media before freezing; is the water pressure; and is the gravity acceleration vector.
is the viscosity of liquid water which could be expressed as a function of freezing temperature [22, 23]: is the relative permeability which could be expressed as a function of water saturation [12, 22, 23]:where is the material constant. The water saturation is equal to unfrozen water content during the freezing process of porous media [12, 22]; then, we can derive the expression of relative permeability combining with equation (12):
When the temperature of the frozen soil drops to the freezing temperature, the liquid water in the soil turns into ice, and the relative porosity of the soil approaches zero; from equation (21) we can find that the seepage velocity becomes zero at this time too. This method that reflected the change mechanism of the seepage velocity in the freezing process by relative permeability is called apparent-heat capacity method. The existing research results are based on the apparent-heat capacity method for mathematical modeling when using Comsol Multiphysics for hydrothermal coupling calculation [12, 15, 17, 21–24]. The main difference between them is the expression of unfrozen water content; however, Vitel et al. discovered through research that the influence of the function shape is very limited and can be neglected, and reasonable results can be obtained by using the Heaviside function to smoothly approximate the curve of unfrozen water content [17], and this processing method can effectively improve the calculation efficiency. Therefore, we chose this method to express the unfrozen water content in this study, as shown in equation (12).
Substituting equation (21) in equation (20), the continuous equation of the seepage field is obtained:
It should be noted that when performing hydrothermal coupling numerical calculations, there is another commonly modeling methodology called the enthalpy-porosity method [19, 20], where the water-ice phase-change interface is modeled as a mushy zone. The transformation from water to ice in this zone is considered as a porous medium, where a modified Darcy source term is used to simulate motion in the phase-change region. The main difference between enthalpy-porosity method and apparent-heat capacity method which is used in this paper is the mechanism that is going to force the velocity to a value of zero in the freezing zone. The mechanism of the seepage velocity change of enthalpy-porosity method is described as follows.
The conservation equation of momentum is [19, 20]where , , , are the Darcy, inertial, buoyancy, and mushy source terms, respectively. They are given as
The Darcy term and inertial term make up the total resistance to the flow. At no-seepage, or low seepage velocity, the inertial term is so low that can be safely neglected [19, 20]. is the Ergun coefficient. The buoyancy source term is used to induce the natural convection within the voids, where is the gravitational acceleration, and is the thermal expansion coefficient. The mushy source term is a modified Darcy source term that is introduced to control the freezing process within the mushy zone [19, 20]. As the temperature of the soil gradually decreases, the liquid fraction within the pore fluid gradually reduces to 0, and mushy source term dominates all other terms in the momentum equation (26) and forces the superficial velocity close to zero [19, 20].
In this paper, all numerical calculations are implemented by Comsol Multiphysics, so the apparent-heat capacity method that depends on the relative permeability approach was chosen in mathematical modeling.
2.3. Numerical Implementation
In this study, the hydrothermal coupling numerical calculation model was constructed by the PDE module in COMSOL Multiphysics, and this model was solved by the fully coupled solver. The solution process is shown in Figure 2.
In Figure 2, is the relative permeability of the frozen soil. Since changes with the content of unfrozen soil during the freezing process, the change of the relative permeability affects the flow field, so acts as an intermediate variable to couple the seepage field to the temperature field.
3. Model Validation
3.1. Design of Large Similar Model Test
The test device consists of a test chamber, a freeze simulation system, a groundwater simulation system, and a data test system. The test chamber has a length of 2,400 mm, a width of 1,500 mm, and a height of 800 mm. Among them, the length of the intermediate sand layer test area is 2000 mm. The test device system is shown in Figure 3.
The geometric similarity ratio selected in the experiment is . According to the similar model test theory, the main parameters in the test are shown in Table 1.
The arrangements of the freezing pipes and the temperature measuring points are shown in Figure 4. The two pipes on the axis ML2 are used to simulate conventional freezing pipes, and another pipe located on the axis ML1 is used to simulate an auxiliary freezing pipe in an optimized freezing pipe arrangement. This paper mainly studies the development law of the freezing temperature field under the action of groundwater, so only the two freezing pipes located on ML2 are functioning.
Medium-coarse sand (particle size range: 0.25 mm–0.5 mm) was used to simulate the porous medium. The tank was filled with a 50 mm thick sand layer each time, which was vibrated and compacted to ensure uniformity. When the sand layer was filled to an intermediate position of the height of the tank, the temperature measuring points were arranged. To facilitate the accurate positioning of the measuring points and prevent the position of the measuring points from changing due to the flow of water during the test, the thermocouple strings were tied to fine bamboo sticks, and the sand layer of the specified plane was leveled. The thermocouple string was laid out as shown in Figure 5. The minimum distance from the measuring points (nos. 23 and 25) to the boundary of the box reaches 510 mm, which is 4.25 times the distance between the freezing pipe and the measuring points (nos. 23 and 25), and the test is carried out in winter, the indoor temperature being stable at 6∼8°C, so the interference of the external environment on the test results can be ignored.
The test system after assembly is shown in Figure 6.
3.2. Establishment of Numerical Calculation Model
The geometric model in the test is three-dimensional. If we neglect freezing pipe deflection during the vertical freezing process, the axially uneven distribution of the freezing pipe wall temperature, and the heterogeneity of the soil, the model becomes two-dimensional. Therefore, the temperature measurement surface is selected for this research.
The transient analysis capabilities of the porous media and subsurface flow module and the hear transfer in porous media module of the COMSOL Multiphysics package are used in this study. The parametric variables are brought into the equation to realize the transient solution of the hydrothermal coupling mathematical model during the freezing process.
The geometric model of the numerical calculation is shown in Figure 7.
We obtain the physical parameters experimentally as follows: the latent heat of the water/ice phase transition is 334 kJ/kg, and the freezing temperature of the sand is −1°C. The relevant physical quantity values involved in the calculation are shown in Table 2.
3.3. Comparison of Numerical Calculation Results and Model Tests
To study the formation law of freezing temperature field under the action of groundwater under different flow rates, the two freezing pipes in the downstream direction are opened in both the test and numerical calculation. The accuracy of the hydrothermal coupling numerical calculation model constructed in this paper is judged by two comparison methods: typical position and frozen range.
3.3.1. Comparison of Typical Position Measurement Data
Measuring point no. 13 is located at the midpoint of the area where the two freezing pipes are located. Changes in the temperature measured at this measuring point with time reflect the formation of the frozen wall formed by the adjacent freezing pipes.
It can be seen from Figure 8 that the temperature variations at point no. 13 calculated by the model presented in this paper closely match the experimentally observed variations for a groundwater flow velocity of 0 m/d. More specifically, (1) the temperature changes are consistent—that is, both experience the first rapid temperature drop phase, the phase transition phase, the second rapid temperature drop phase, and the stable freezing phase—and (2) the time corresponding to each stage is basically the same—the first rapid temperature drop phase lasts for about 100 minutes, the closure of the frozen column formed by adjacent frozen pipes occurs at 250 minutes, and the end of the subrapid temperature drop phase occurs at 600 minutes. The test data showed slight fluctuations in the later stage, but the overall change law remained consistent with the simulation results. The comparison between the numerical calculation results and the experimental data of other measuring points is shown in Figure 9. It should be noted that the position of the measuring points around the two freezing pipes is symmetrically distributed along the axis ML1. Therefore, only the data of the measuring points around one of the freezing pipes are analyzed in this part, and the overall distribution law of the freezing temperature field will be presented in the following section by the form of an isotherm. By comparison, it can be found that the numerical calculation results have a high consistency with the overall development law of the test data of most of the measuring points under the condition that the seepage velocity is zero, and only the difference between the numerical calculation result and the experimental result of the no. 15 measuring point is obvious, the reason for the difference may be that the detection data of the measuring point are inaccurate due to some uncertain factors. Since the sand layer is not completely uniform during the test, there are some slight differences in the test data of the measuring points symmetrical on both sides of the axis ML2.
(a)
(b)
(c)
(d)
It can be seen from Figure 10 that the overall development trend of the temperature at measuring point no. 13 is consistent between the numerical calculation and experimental tests. However, the calculated value is generally slightly higher than the test result, and the calculated end moment of the phase change is delayed by about 20 min. There are several reasons for this. First, in the test, the water flow velocity is controlled by the pressure of the inlet section, and the sand layer is assumed to be uniform in each cross section. However, due to the large size of the test box, the uniformity of the sand layer cannot be guaranteed during the process of sand filling and compaction. Therefore, there is a discrepancy between the actual flow rate of the water flowing between the two adjacent freezing pipes and the rate used in calculations. When the control flow rate is slightly smaller than the design value, the experimentally obtained closure time will be earlier than the numerical calculation result. The comparison between the numerical calculation results and the experimental data of other measuring points is shown in Figure 11. By comparison, it can be found that the numerical calculation results of most measuring points have a high consistency with the overall development law of the test data under the condition that the seepage velocity is 5 m/d except the measuring point no. 15, the reason for the obvious difference between the numerical calculation result and the experimental result of the only one measuring point may be that the detection data of the measuring point are inaccurate due to some uncertain factors. Due to the action of the water flow, the temperature of the measuring point on the right side of the axis ML2 is significantly lower than that of the left side.
(a)
(b)
(c)
(d)
3.3.2. Comparison of Frozen Range
Figure 12 shows that the frozen wall formed by the two freezing pipes has been closed when the freezing time reaches 350 min. In addition, the −1°C isotherm at the upstream position has a pronounced concave shape in the middle of the action zone of the two frozen pipes. Due to the action of the water flow, the frozen area on the right side of the axis ML2 is significantly larger than the frozen area on the left side. As the freezing time continues to increase, the frozen area gradually expands. When the freezing time reaches 500 min, the −1°C isotherm at the upstream position is smoother, but the middle of the action zone of the two frozen pipes still has a concave shape along the direction of water flow. After 650 min, the frozen range was further expanded, but the overall shape of the frozen area did not change significantly.
(a)
(b)
(c)
In order to further verify the accuracy of the numerical calculation model, we compare the frozen range obtained by numerical calculation after 350, 500, and 650 minutes under the action of groundwater moving at 5 m/d with the frozen range obtained by the test (Figures 13 and 14).
(a)
(b)
(c)
When the freezing time was 350 min, 500 min, and 650 min, the differences between the frozen range calculated by the numerical calculation and the frozen range obtained by the test were only 2.4%, 4.2%, and 6.2%, respectively, and the shape of the frozen area obtained by the two was more consistent. This shows that the mathematical model and the selected solution method meet the actual needs of the project. The numerical calculation model can provide strong support for the study of the formation of an artificial freezing temperature field under the action of groundwater.
4. Engineering Applications
4.1. Study of the Freezing Rule for Single-Row Pipes under the Action of Groundwater
The central wind shaft of Panyi Mine in Huainan Mining Area was constructed by the freezing method. The wellbore passes mainly through Cenozoic alluvium and Permian Shihezi coal-bearing strata. The thickness of the Cenozoic alluvium is 205.8 m, and it is composed mainly of clay and sand. The sand layer is rich in groundwater and is the key control layer of the whole freezing project. The freezing parameters of the central wind shaft in Panyi Mine are shown in Table 3.
We now apply the numerical calculation model established in the previous section to the arrangement of pipes found in the central wind shaft of the Panyi Mine. And the sand-permeable stratum at a depth of 200 m is selected for analysis. The numerical calculation model of stratum freezing under the action of groundwater is constructed as shown in Figure 15.
In this model, the calculated control range is 28 m × 28 m. The initial temperature of both the groundwater and the formation is 20°C. The axis L1 is parallel to the direction of the water flow, the axis L2 is perpendicular to the direction of the water flow, and the intersection of L1 and L2 coincides with the center of the freezing pipe arrangement circle. For the convenience of analysis, the left and right intersections of the L1 and the freezing pipe arrangement circle are, respectively, defined as the “upstream position” and “downstream position,” and the intersection point of L2 and the freezing pipe arrangement circle is defined as the “side position.”
According to test results of freezing temperature of saturated sand, an effective frozen wall is formed when the temperature of the frozen sand layer reaches −1°C; this means that the −1°C isotherm is the contour of the frozen wall. The distributions of the freezing temperature field under the action of different flow rates at 30 d, 60 d, and 90 d are plotted in Figure 16.
(a)
(b)
(c)
(d)
(e)
Using the data in Figure 16, we can analyze variations in the artificial freezing temperature field under the action of groundwater with different flow rates. We find that, due to the convective heat transfer of the water flow, part of the cooling capacity of the upstream freezing pipes is carried by the water flow to the downstream region. The end result is that the closure time of the downstream position of the frozen wall is earlier than that of the upstream position. For t = 30 d, under the action of groundwater less than 10 m/d, the downstream area of the frozen wall has completed the intersection, while the upstream position of the frozen wall has not been closed. For t = 60 d, the downstream area of the frozen wall has been closed, and the upstream position can only be closed under the condition that the flow rate is less than 7 m/d. For t = 90 d, the frozen range of the downstream area has formed a distinct convex shape along the direction of the water flow, and the frozen range on both sides has expanded to the intended excavation position, while the upstream position can only be closed if the groundwater flow rate is less than 9 m/d.
Numerical calculation shows that groundwater flowing at less than 5 m/d has little effect on the freezing temperature of single-row pipes. As a result, in this paper, only the influence of groundwater with a flow rate greater than 5 m/d is studied. Figure 17 shows that the closure time of the frozen wall formed by the single-row freezing pipes increases with the flow rate. The closure time under the action of 5 m/d water flow is 30 d, while the closure time under the action of 11 m/d water flow is 171 d.
Whether the downstream area can be closed plays a decisive role in the formation of the entire frozen wall. After the frozen wall of the downstream zone is closed, there is no longer movement of water throughout the frozen zone inside the circle of the freezing pipe arrangement; as a result, the effect of water flow on the freezing temperature field is weakened. Therefore, as long as the frozen wall in the downstream position can complete the intersection, the entire frozen wall will be closed. Through numerical calculation, we find that the closure time of the frozen wall at the downstream location shows a parabolic relationship with the water flow velocity. When the water flow velocity is less than 7 m/d, the closure time of the frozen wall in the downstream region decreases with increasing flow velocity; when the flow velocity is greater than 7 m/d, the closure time increases with increasing flow velocity. The reasons for this phenomenon are as follows: due to convective heat transfer, part of the cooling released by the freezing pipes is used to reduce the temperature of the water flow. When this part of the low-temperature water flows through a higher temperature region, part of the cooling capacity is transferred to this area, and the water flow becomes the flow medium of the heat. This realizes the process of transferring the cooling capacity of the upstream freezing pipes to the downstream area. Therefore, when the water flow velocity is within a certain range, the increase of the flow velocity promotes the process of energy transfer. However, as the flow rate is further increased, the replenishing cold from upstream cannot compensate for the energy loss caused by the convective heat transfer of the water stream to the downstream region. The effect of this part of the heat from the upstream on the freezing process in the downstream area is weakened. When the groundwater flow rate is further increased, the closure time of the frozen wall in the downstream area increases.
In fact, when the flow rate of groundwater is large, the unfrozen area in the upstream position is still subjected to the impact of a large water flow after the downstream frozen wall is closed; this position will not be closed until the internal extent of the freeze pipe arrangement circle is frozen. This process requires continuous cooling for a long time, resulting in large energy losses. During actual construction, the soil inside the pipe arrangement is not allowed to completely frozen, in order to prevent adverse effects on the excavation of the soil and the construction of the wellbore. Therefore, it is necessary to comprehensively consider the above factors in determining the water flow rate allowed by the freezing method. In the calculations herein, we find that for single-row freezing pipes, the maximum flow rate is 10 m/d.
Figure 18 shows the development law of the thickness of the frozen wall under different flow rates at a freezing time of 90 d. Through comparative analysis, we find that the groundwater flow rate has the most significant effect on the thickness of the frozen wall at the upstream position; the thickness of the frozen wall of this position decreases sharply with the increase of the water flow velocity. When the water flow velocity is greater than 9 m/d, the thickness of the frozen wall at this position is 0, indicating that the frozen wall in the upstream position cannot be closed within 90 d. The frozen wall on both sides is less affected by the flow rate. When the groundwater flow rate is increased from 5 m/d to 12 m/d, the thickness of the frozen wall in this area decreases with the increase of the flow velocity, but the decrease is not obvious. Due to the superposition effect of the cold volume from the upstream and downstream regions, the thickness of the frozen wall in the downstream region is always greater than the thickness of the frozen wall at other locations when the freezing time reaches 90 d.
It can be seen from Figure 19 that when the initial flow rate of groundwater is large, the flow rate will undergo the following changes as the freezing time increases: (1) In the first stabilization phase, the frozen area and the resistance to water flow are small. The groundwater flow rate will stabilize near the initial flow rate. Since the placement of the freezing pipes has a certain influence on the permeability of the soil, the initial flow rate in the middle position of the adjacent freezing pipes of the downstream region is slightly larger than the design rate. (2) The second stage involves a rapid increase in flow rate. Since the permeability of the frozen area is small, the water flow can only pass through the unfrozen area in the middle of the two adjacent freezing pipes. Therefore, under a constant pressure head, the smaller the seepage area, the greater the water flow rate. (3) The next stage is a slow growth phase. When the frozen wall formed by the adjacent freezing pipes is about to be closed, the water flow rate at the midpoint of the line connecting the two freezing pipes increases with the increase of the freezing time. As can be seen from Figure 20, once the water flow speed increases to a certain value, it decreases sharply to 0, and this turning point is the closure time of the frozen wall. The reason for this change in the flow rate is that as the volume of the frozen column formed by the adjacent two freezing pipes increases, the permeable region of the cross section decreases continuously, but the water pressure from the upstream does not change, resulting in a gradually increasing water flow. Through the calculations in this paper, we find that groundwater with an initial flow rate of 10 m/d can reach 7 times the initial flow rate at the moment of closure.
We now analyze the closure time of the frozen wall for different pipe spacings (Figure 21). Through comparative analysis, we find that the spacing of the freezing pipes has a significant influence on the final closure time of the frozen wall. Under the action of groundwater with a flow rate of 5 m/d, when the spacing of the freezing pipes is reduced from 1.15 m to 1.05 m, the closure time of the entire frozen wall is decreased by 11 d, and when the spacing of the freezing pipes is increased from 1.15 m to 1.30 m, the closure time of the entire frozen wall is increased by 20 d. With the increase in the flow rate, the influence of the change of the freezing pipe spacing on the freezing efficiency of the frozen wall becomes more apparent. Under the action of groundwater with a flow rate of 9 m/d, when the spacing of the freezing pipes is reduced from 1.15 m to 1.05 m, the closure time of the entire frozen wall is decreased by 52 d, and when the spacing of the freezing pipes is increased from 1.15 m to 1.30 m, the closure time of the entire frozen wall is increased by 34 d. When the flow rate of groundwater reaches 10 m/d, the frozen wall formed by the freezing scheme with the freezing pipe spacing of 1.30 m cannot been closed within 180 d. Therefore, when there is a large flow rate of groundwater in the local layer, the spacing of the freezing pipes should be reduced in the design of the freezing scheme to improve the freezing efficiency of the entire region. We conclude that when the single-row pipe arrangement is adopted, the spacing of the freezing pipes should not be greater than 1.3 m. However, considering the possible existence of freezing pipe deflection, the spacing of the freezing pipes should not be less than 1 m. Therefore, the spacing of the freezing pipes can be set to 1–1.3 m.
4.2. Study on the Freezing Law of Double-Row Pipe under the Action of Groundwater with Different Flow Rates
The wind well in the north area of Zhangji Mine in Huainan Mining Area was constructed by the freezing method. The wind well intersects Quaternary sediments to a depth of 271.25 m and Permian strata from 271.25∼1067.80 m; the total thickness of these Permian sediments is 796.55 m. The bedrock weathering zone has a bottom boundary of 317 m and a thickness of 45.75 m, and the strong weathering zone has a depth of 293.25 m and a thickness of 22 m. At the depth range of 380 m–430 m is a water-rich sand layer, which is the key control layer of the whole freezing project. The freezing scheme design parameters are shown in Table 4.
Based on the design parameters of the freezing scheme, the sand-permeable stratum at a depth of 400 m is selected for analysis, and the numerical calculation model of the double-row pipes stratum freezing under the action of groundwater is constructed as shown in Figure 22.
In this calculation model, the calculated control range is 34 m × 34 m. The initial temperature of the groundwater and the formation is 20°C. The temperature of the freezing pipe is −32°C. The axis L1 is parallel to the direction of the water flow, the axis L2 is perpendicular to the direction of the water flow, and the intersection of L1 and L2 coincides with the center of the freezing pipe arrangement circle. For the convenience of analysis, the left and right intersections of the L1 and the main freezing pipe arrangement circle are, respectively, called the “upstream outer ring position” and the “downstream outer ring position,” and the left and right intersections of L1 and the auxiliary freezing pipe arrangement circle are, respectively, the “upstream inner position” and “downstream inner ring position.”
Using Figures 23 and 24, the variation law of temperature field formed by double-row freezing pipes under the action of groundwater with different flow rates can be analyzed.
(a)
(b)
(c)
(d)
(e)
In the case that there is no significant difference in the spacing of the main freezing pipes, since a circle of auxiliary freezing pipes is arranged inside an arrangement ring of the main freezing pipes, the resistance of the entire frozen area to groundwater is significantly improved. Groundwater with a flow rate below 10 m/d does not have a significant effect on the freezing time of the frozen wall formed by double-row freezing pipes, and when the groundwater flow rate reaches 20 m/d, a frozen wall can still form successfully. When the freezing time reaches 30 d, the outer ring area of the frozen wall downstream under the action of groundwater with a flow rate less than 20 m/d has already been closed, and the frozen wall formed by the side inner-ring freezing pipes has been connected to the frozen wall formed by the outer-ring freezing pipes. Increases in the flow rate have a significant effect on the frozen wall in the downstream area, but the frozen walls on both sides do not change significantly.
When the freezing time reaches 60 d, the frozen area under the action of groundwater with a flow rate of less than 14 m/d has been closed, and the extent of the frozen wall is close to the intended excavation boundary. When the flow rate is greater than 14 m/d, the frozen wall of the upstream region is not closed, and as the flow rate increases, the unfrozen region increases. When the freezing time reaches 90 d, the frozen area under the action of groundwater with a flow rate of less than 16 m/d has been closed, and the extent of the frozen wall has exceeded the control boundary of the proposed excavation area. Under the action of groundwater with a flow rate greater than 16 m/d, the frozen area of the upstream position is not closed, but the frozen range on both sides has developed to the inside of the boundary to be excavated.
Figure 24 shows the closure times of different positions of the frozen wall. We find that the closure time of the downstream outer ring position is earlier than that of the downstream inner ring and the upstream position. When the water flows through this position, it has passed the precooling action and the blocking action of the three rows of freezing pipes, and the convective heat transfer effect of the water flow with the frozen area is weakened, so the closure time of the position is less affected by the flow rate. When the water flow rate is less than 14 m/d, the closure time at this position remains unchanged at 17 d. When the groundwater flow rate is greater than 14 m/d, the closure time increases with increasing flow rate. When the groundwater flow rate reaches 20 m/d, the closure time of the downstream outer ring position is 34 d. The closure time of the position of the outer ring and the inner ring upstream maintains the same development trend throughout the freezing process. When the flow rate is low, the influence of the water flow on the freezing efficiency is small. Since the spacing of the auxiliary freezing pipes is larger than that of the main freezing pipes, the closure time of the upstream outer ring position is earlier than that of the inner ring position. With increasing flow rate, the convective heat transfer of the water flow causes the cooling capacity of the freezing pipes at the outer ring position upstream to be largely lost, while the cooling loss of the frozen region at the upstream inner ring position is small. Therefore, the closure time of the position of the upstream inner ring is earlier than that of the outer ring area.
However, since the row spacing between the freezing pipes is small, after the frozen wall formed by the inner ring pipes is closed, the water flow bypasses the frozen zone and flows to the downstream zone, and the influence of the water flow on the action area of the freezing pipes of upstream outer ring is weakened. Therefore, the frozen wall at this position is closed in a short time after the closure of the frozen wall of the upstream inner ring. When the groundwater flow rate is greater than 10 m/d, the intersection time of the downstream inner ring position is earlier than that of the upstream position. After the frozen wall is closed, the unfrozen area inside the frozen wall is no longer affected by flowing groundwater; therefore, the time when the frozen area reaches the intended excavation boundary will not change significantly under the flow velocity within 15 m/d. When the flow rate is further increased, the interval between the closure time of the frozen wall of the upstream position and downstream area is longer. Therefore, before the frozen wall in the upstream position is closed, the frozen areas on both sides and downstream have developed to the position to be excavated. In this case, the closure time of the frozen wall of the upstream outer ring is basically the same as the time when the entire frozen wall extends to the boundary to be excavated.
To study the influence of the row spacing of the freezing pipes on the development law of the freezing temperature field, the row spacing is adjusted to 2.5 m and 1.5 m while keeping the number of auxiliary freezing pipes unchanged. The closure time and the time when the frozen wall develops to the boundary to be excavated corresponding to the two cases are compared with those of the original freezing scheme (Figure 25; the curve formed by the solid figure represents the variation of the time required for the frozen wall to reach the excavation boundary with the flow rate, and the curve formed by the hollow figure represents the variation law of the closure time of the frozen wall with the flow rate).
According to Figure 25, the influence of the row spacing on the closure time of the frozen wall can be found. Under the action of groundwater with a small flow rate, when the row spacing increases, the time to complete freezing becomes longer. The reason is that when the groundwater flow rate is small, the time when the frozen wall formed by the main freezing pipes and the freezing wall formed by the auxiliary freezing pipes are integrated is the main factor affecting the final closure time of the entire frozen wall. When the spacing of the rows is increased, intersection of the two takes longer, which leads to a longer closure time of the entire frozen wall.
Under the action of groundwater with a large flow rate, when the row spacing increases, the closure time of the entire frozen wall becomes shorter. The reason for this is that as the groundwater flow rate increases, the intersecting time of the frozen wall formed by the upstream main freezing pipes becomes the main factor affecting the closure time of the entire frozen wall. The increase of the row spacing leads to the reduction of the pipe spacing of the auxiliary freezing pipes, and the intersection time of the frozen wall formed by the auxiliary freezing pipes becomes shorter. After that, the impact of the water flow on the frozen zone of the main freezing pipes is weakened, thereby increasing the freezing efficiency of the main freezing pipe and shortening the closure time of the entire frozen wall.
It can be seen from Figure 25 that the larger the row spacing, the shorter the time required for the frozen wall to develop to the excavation boundary. The reason is that after the row spacing is increased, the distance between the auxiliary freezing pipe arrangement circle and the excavation boundary is shortened, so the time required for the frozen wall to expand to the excavation boundary is reduced.
Further analysis shows that when the row spacing of the freezing pipes is large, it takes a long time for the frozen wall formed by the auxiliary freezing pipes to intersect with the frozen wall formed by the main freezing pipes, resulting in the excessive expansion of the frozen wall to the interior of the excavation boundary in the late stage of freezing. At the same time, due to the large spacing between the rows, the temperature of the frozen soil between the main freezing pipe arrangement ring and the auxiliary freezing pipe arrangement ring is relatively high, and the overall strength of the frozen wall cannot meet the design requirements, which adversely affects the later excavation. When the row spacing of the freezing pipe is small, the cooling capacity will be concentrated near the freezing pipe arrangement circle, and the time required for the freezing wall to reach the design thickness becomes long, which causes waste of energy and delays the construction period of the excavation.
In order to achieve a reasonable distribution of cooling capacity and ensure that the strength and thickness of the frozen wall meet the design requirements, the row spacing of the freezing pipes should be controlled within a reasonable range. Through calculation, the reasonable range of the row spacing of the double-row freezing pipe design scheme is 2∼2.5 m.
5. Conclusions and Recommendations
We constructed a numerical calculation model of hydrothermal coupling based on the law of conservation of energy and mass. This model considers groundwater flow, heat change during the freezing process, and the influence of the absolute porosity reduction on seepage during the freezing process. The accuracy of the model was verified by large-scale laboratory tests. Based on the actual layout scheme of the frozen pipes, numerical calculation models of the freezing temperature field of the single-row pipes and double-row pipes under the action of groundwater were constructed. The conclusions drawn from this research are as follows:(1)Groundwater with a flow rate of less than 5 m/d does not have a significant effect on the artificial freezing temperature field. The layout scheme of single-row freezing pipes is suitable for the formation of freezing under the action of the groundwater with a flow rate of less than 10 m/d, while the layout scheme of double-row freezing pipes is suitable for formation freezing under the action of the groundwater with a flow rate of less than 20 m/d.(2)During the development of the freezing temperature field of the single-row pipes, the groundwater flow rate increases as the range of the unfrozen area decreases because groundwater flows through the middle of the frozen area between the adjacent two freezing pipes. When the frozen wall is about to be closed, the groundwater flow rate at this position can reach 5 to 7 times the initial flow rate, resulting in slower expansion of the freezing range just before closure.(3)During the freezing process, the frozen range upstream is most affected by the water flow; increasing water flow will cause the frozen area upstream to decrease sharply, while the frozen wall areas on both sides are less affected by the water flow. Therefore, in the design of the arrangement of the freezing pipes, the installation density of the freezing pipes in the upstream regions should be appropriately increased, thereby improving the freezing efficiency of the entire freezing pipe action area.(4)The pipe spacing and the spacing of the rows have significant effects on the formation of the entire frozen wall. In the design of the single-row pipes, a reasonable range of pipe spacing is 1–1.3 m. In the design of the double-row pipes, a reasonable range of the row spacing is 2–2.5 m.(5)A “group-pipe effect” will occur, when the three freezing pipes are arranged in a triangle. In the future research, we will study the formation law of this “group-pipe effect” by changing the row spacing and the pipe spacing under the action of groundwater with different flow rates.
Data Availability
The data used to support the findings of this study are available from the corresponding author upon request.
Conflicts of Interest
The authors declare that they have no conflicts of interest.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (grant nos. 51878005 and 51374010) and National Science and Technology Major Project (no. 2016YFC0600902-01).