Abstract

An artificial ground freezing method is often applied to highly permeable gravel formations. Seepage flow increases energy consumption and engineering accidents under this condition. Based on the physical modeling tests, the numerical simulations were conducted. The physical modeling test examined the development of the temperature and frozen wall. Numerical modeling integrated with the ACO algorithm was established to optimize the layout parameters of the freezing pipes. The results indicate that the flowing water prolongs the closure time of the frozen wall. Meanwhile, the total thickness of the frozen wall is also reduced by the flowing water. There are significant differences in the development rates of the frozen wall in different zones. The thickness of the entire frozen wall is nonuniform owing to the seepage flow. Following optimization using the proposed algorithm, the closure time was shortened from 82.4 d to 56.9 d for the frozen wall. Moreover, the freezing efficiency increased by 30.95% after optimization, and consequently, the entire frozen wall was more uniform with a nonweak zone. A case study showed that this optimization system is an effective method for artificial ground freezing operations in geotechnical engineering.

1. Introduction

Artificial ground freezing (AGF) is a common measure in geotechnical interventions [1, 2], such as mining shafts, subway tunnels, foundation pits, and slope stabilization. Besides good environmental friendliness, AGF has the advantages of water blocking, high strength and stiffness of soil, and good formation adaptability. The AGF technique is often used in tunnel building with complex formation conditions: for example, in water-rich and low-strength formations [3]. Moreover, AGF has been extensively used as a technique to reliably reduce the risks of damage during the building of tunnels ([46], [7]).

In water-rich formations, extensive research has been conducted on the application of AGF in soft soil formation [8, 9]. In contrast to soft soil formation, sandy soil, pebble soil, and gravel soil formations have the problem of high permeability. When high permeability formation is encountered during the AGF process, the coupling effect of seepage flow and freezing pipes should not be neglected. In engineering practices, such formation often presents difficult merging of frozen wall and even disorganizes AGF working plan [3]. Combined with physical modeling, numerical simulation, and theoretical analysis, Pimentel et al. [2] found that the seepage flow significantly prolonged the closure time of the frozen wall in sand formation.

Previous studies have shown that seepage flow has great influences on the closure time, shape, and thickness of the frozen wall. Frozen walls require more freezing time, smaller spacing, or multirow freezing pipes to reach the designed thickness [10, 11]. However, these measures are not optimal solutions because of the significantly higher construction costs. Besides, most designers rely on empirical formulas to design AGF, which has great constraints in the face of complicated engineering geological and hydrological formation conditions. Thus, the application of an optimization algorithm provides an intelligent and efficient method. Several specific algorithms have been developed within the family metaheuristic optimization methods, such as the neural network method, genetic algorithm, simulated annealing, and ant colony optimization (ACO) [12].

Among these methods, ACO is aimed at obtaining the optimal path in a graph by simulating the behavior of ants seeking a path for food ([13, 14], and [15]). Moreover, ACO has been used to solve different engineering and scientific problems, such as geotechnical engineering, optimization of pipe layout, and commercial operation ([15, 16], and [17]). In geotechnical engineering, Tian and Zhou [18] implemented ACO to invert the soil mechanical parameters. Gao and Feng [17] suggested a new immunized continuous ant colony algorithm to solve the optimization problem by improving the calculation efficiency in the back analysis of geotechnical engineering. Shi et al. [19] integrated ACO with a genetic algorithm to seek the slip surface of a geotechnical slope. In conclusion, ACO is an effective and simple algorithm for optimizing the layout parameters of freezing pipes in this study.

Above all, this study is aimed at determining the optimal layout parameters of freezing pipes in a tunnel-connected aisle under seepage flow, such that the minimum freezing duration required to form a fully frozen wall is obtained. To this end, the numerical simulation model was integrated with the ACO algorithm. The study also briefly presents the laboratory physical modeling test, the development of numerical modeling, and the parameterization of freezing pipes in a tunnel-connected aisle. Lastly, issues regarding the optimization evaluation and improvement of the freezing efficiency are discussed.

2. Physical and Numerical Modelling of Horizontal Single-Row-Pipe Freezing

2.1. Physical Modelling Test

Considering the seepage flow, several large-scale physical modeling tests were performed. The result analysis of the physical modeling tests can be found in our published paper [20]. The physical and thermal properties of the gravel soil are summarized in Table 1 for the numerical modelling.

Based on the physical modeling test presented in Wang et al. [20], the approximate closure times and thicknesses of the frozen walls are summarized in Table 2. Evidently, the closure time of the frozen wall increases due to the seepage flow, whereas the total frozen wall thickness decreases. Moreover, the development of the upstream thickness is restricted owing to the ‘scouring effect’ of the flowing water. However, the downstream thickness increases owing to the “water barrier effect” of upstream frozen wall. More specifically, the seepage flow is favorable to the downstream frozen wall but unfavorable to the upstream frozen wall.

2.2. Model Equations and Parameters
2.2.1. Model Assumptions

The application of AGF under seepage flow presents a complex thermohydraulic problem involving temperature and seepage boundary, phase change, and internal heat source. The numerical modeling was developed based on the following assumptions. (1)Gravel soil formation is considered a continuous, uniform, and isotropic saturated porous medium with interconnected pores. The thermohydraulic problem in gravel formation is assumed to be a two-dimensional (2D) model of seepage-heat transfer in a porous medium(2)The thermohydraulic problem is defined as seepage-heat transfer in saturated gravel formation, and the heat transfer mode is heat convection and heat conduction. The ice-water phase change occurs only in the temperature range of 0 to -0.04°C. Gravel soil density and mass migration with temperature were ignored during the AGF process(3)The seepage flow in the formation conforms to Darcy’s law. The model focuses only on the macroscopic water flow. Additionally, the temperatures of the flowing water and soil particles are identical

2.2.2. Governing Equation of Energy Conservation and Heat Conduction

Based on the heat conduction differential equation, the corresponding energy conservation equation is obtained by considering the ice-water phase change [3]. where denotes the equivalent volume heat capacity (m2/(s·°C)), is the temperature (°C), denotes the density of flowing water (kg/m3), represents the latent heat of ice-water phase change (J/kg), indicates the volume content of flowing water, denotes the equivalent heat conductivity coefficient (J/(kg·°C)), corresponds to the specific heat capacity of flowing water (m2/(s·°C)), indicates the velocity vector of flowing water (m/d), signifies the Hamiltonian operator, and denotes the heat source.

According to the basic properties of the porous medium, the parameters were identified as follows: where denotes the thickness of the gravel formation (m); denote the specific capacity of ice-water mixture and soil particles, respectively (m2/(s·°C)); is the velocity vector (m/d); signify the density of ice-water mixture and soil particles, respectively (kg/m3); represents the volume fraction of soil particles in the total volume; denote the heat conductivity coefficient of ice-water mixture and soil particles, respectively (J/(kg·°C)); denotes the tensor of dispersion heat conductivity coefficient (J/(kg·°C)); and correspond to heat sources.

2.2.3. Governing Equation of Seepage Flow

According to the assumptions in Section 2.2.1, the seepage differential equation is identified as follows [1, 3]: where denotes the porosity of the gravel formation; represents the heat source; signifies the pore water pressure (Pa); denote the compressibility of water and soil particles, respectively (1/Pa); indicates the permeability coefficient of gravel formation (m/d); and denotes the dynamic viscosity of flowing water (N·s/m2).

2.2.4. Thermohydraulic Coupling Equation

In the application of AGF in gravel formation under seepage flow, the ice-water phase change changes the gravel formation permeability and affects the seepage field. Meanwhile, the seepage flow directly affects the gravel formation temperature, and this thermohydraulic coupling is achieved by the governing equation of the ice-water phase change [1]. where denote the density of ice and water during phase changing (kg/m3), denotes the volume content of ice and water, represent the specific capacities of ice and water, respectively (m2/(s·°C)), indicates the latent heat of ice-water phase change (J/kg), denotes the equivalent heat conductivity coefficient of ice-water mixture (J/(kg·°C)), and are the heat conductivity coefficients of water and ice, respectively (J/(kg·°C)).

2.3. Model Validation

The numerical modelling of the single-row pipe freezing was built, and its performance was investigated. The model parameters are the same as those of the physical model, as presented in Section 2.1. Figure 1 illustrates the spatial distribution of the gravel formation temperature at the closure time of the frozen wall. In Figure 1, the dimensions and seepage flow rate are calculated using the similarity constant, as listed in [20]. Namely that, the seepage flow rate in is the rate in physical modeling test. The evolution of the temperature distribution shows good agreement with the physical modeling test data. The numerical results replicated the physical model test results well. The modeling results reveal that the seepage flow delays the closure time of the frozen wall and causes a decrease in the frozen wall thickness. In Table 3, the result differences between the numerical and physical models are summarized. The differences are obvious under higher seepage flow.

3. Numerical Simulation

3.1. Model Setup

The numerical model of the tunnel connected aisle was established by referring to a tunnel-connected aisle with a single-row-pipe freezing method in Jinan, China. The dimensions of the gravel formation have a width of 30 m and a thickness of 40 m, and the related gravel soil properties are listed in Table 1. The designed frozen wall is shown in Figure 2. Evidently, the frozen wall in vertical direction is 6 m in height, the radius of the arch frozen wall at the top arch zone is 3 m, and the frozen wall of the bottom zone is 6 m in width. As a design requirement, the thickness of the frozen wall in each zone should be larger than 1.0 m.

In the numerical simulation model, the groundwater flows horizontally from left to right at a flow rate of 0.50 m/d. The spacing between the freezing pipes was 100 mm with equidistant distribution, the diameter of each freezing pipe was 89 mm, and the temperature of the freezing pipes was −20°C. COMSOL Multiphysics 5.4 was adopted to develop the model. The numerical simulation model adopted an adiabatic boundary with a free triangular mesh division, as shown in Figure 3. Gravel formation around the freezing pipes was established as boundary layer meshing to fit the heat diffusion mode of the freezing pipes.

3.2. Frozen Wall Merging of Tunnel Connected Aisle

Figure 4 illustrates the simulation results of the gravel formation temperature for different freezing durations. Since the flowing water provides constant heat, the heat diffused by the freezing pipes is carried to the downstream zone, which acts on a zone other than the designed frozen wall.

Initially, the flowing water almost homogeneously flowed passed the whole cross-section (see Figure 4(a)). As the frozen columns develop, the seepage flow rate increases greatly within the spaces between the freezing pipes, inhibiting the merging of the frozen wall between adjacent freezing pipes. Meanwhile, the development of the frozen wall presents inhomogeneity at different zones owing to the different spatial layouts of the freezing pipes. The development of the frozen wall at both the top arch and bottom zones is faster than that at the vertical upstream and downstream zones (Figures 4(b) and 4(c)). Once the total frozen wall is connected, there is no more flowing water within the interior of the frozen wall, and then, the effect of flowing water on the temperature is considerably reduced (Figure 4(d)). Thereafter, the frozen wall in the upstream zone exhibited a much faster growth towards the inward direction than outward direction, while the frozen wall in the downstream zone grew much faster towards the outward direction than inward. Finally, the designed thickness was reached with a nonuniform frozen wall thickness.

4. Optimization Process

Based on the previous sections, it can be concluded that seepage flow considerably delays or even prevents the merging of the designed frozen wall. An equidistant distribution of the freezing pipes is not an optimal method under the effect of seepage flow. Therefore, the spatial layout of the freezing pipes should be optimized to improve the freezing efficiency and obtain a uniform frozen wall. The ACO algorithm method was employed to search for the optimal location of freezing pipes in a tunnel-connected aisle.

In this study, the location of each freezing pipe is considered a node to a path sought by artificial ants, and a type of freezing pipe layout is obtained after all the locations of each freezing pipe are confirmed. Subsequently, this type of freezing pipe layout is calculated using the numerical simulation model to evaluate the freezing efficiency, and the evaluation results are the corresponding pheromone levels. According to the pheromone level, artificial ants choose the path with many pheromones as the current optimal solution in each iteration. And then, repeated iteration calculation can optimize the layout parameters of freezing pipes.

The numerical simulation model is integrated within the ACO algorithm to determine the optimal location of the freezing pipes. In doing so, the layout optimization system of the freezing pipes can be established. For the ACO algorithm developed in this study, the parameterization of the freezing pipe layout, pheromone update, and behavior of ants seeking a path are established in the following sections.

4.1. Parameterization of Freezing Pipe Layout

For the layout of freezing pipes, the ACO algorithm depends on an accurate description of the solution domain. The domain is the behavior environment of artificial ants. Therefore, it is necessary to parameterize the domain to be optimized.

As mentioned in the previous sections, the freezing duration should last for 82.4 d to achieve the designed frozen wall thickness for the tunnel connected aisle. The development of the frozen wall is helpful in defining the domain to be optimized for the layout of the freezing pipes. The frozen wall edges at freezing durations of 40 d and 50 d are shown in Figure 5. Evidently, the upstream and downstream frozen wall thicknesses are significantly different. The center line of the entire frozen wall is toward the downstream side due to seepage flow. Furthermore, the development rate of the frozen wall at different locations also exhibited differences, leading to evident weak zones and longer freezing durations. Thus, the center line offset and uneven thickness of the frozen wall are the solution domains for the ACO algorithm, and the solution domain can be divided into upstream and downstream zones, top arch zones, and bottom zones, as shown in Figure 5. Two parameters, , are defined as the offset of the freezing pipes to determine the solution domain, and the spatial distribution density of the freezing pipes is parameterized using the mean () and standard deviation () of the Gaussian distribution.

4.1.1. Layout Parameterization of Freezing Pipes at the Vertical Frozen Wall of Upstream and Downstream Zones

The horizontal coordinate is determined by where the symbols represent the upstream and downstream zones, respectively, in the following parameters: and correspond to the horizontal coordinates of the freezing pipes, denote the horizontal coordinates of the frozen wall center, and represent the horizontal coordinate offsets of the freezing pipes, and denote the labels of the freezing pipes, and correspond to the numbers of the freezing pipes, and is defined as where denote the offsets of the freezing pipes, the ranges of which depend on the designed thickness of the frozen wall .

The vertical coordinate is determined by where denote the vertical coordinates of the freezing pipes, denote the vertical coordinate offsets of the freezing pipes, represents the vertical coordinate of the bottom frozen wall center, and is defined as where and denote the mean of the freezing pipes and and denote the standard deviations of the freezing pipes with representing the Gaussian distribution.

4.1.2. Layout Parameterization of Freezing Pipes at the Frozen Wall of Top Arch Zone

The horizontal coordinate is determined by

The vertical coordinate is determined by where denotes the label of the freezing pipes, represents the vertical coordinate of the top arch frozen wall center, denotes the number of freezing pipes, and the radial coordinate offsets of the freezing pipes is defined as where denote the offsets of the freezing pipes, the ranges of which depend on the designed thickness of the frozen wall . The angular spacing function of the freezing pipes is defined as where denotes the mean of the freezing pipes and corresponds to the standard deviation of the freezing pipes with representing the Gaussian distribution in Equation (9).

4.1.3. Layout Parameterization of Freezing Pipes at the Frozen Wall of Bottom Zone

The horizontal coordinate is determined by where denotes the horizontal coordinate of the freezing pipes, indicates the mean of the freezing pipes, represents the standard deviation of the freezing pipes with representing the Gaussian distribution in Equation (9), and denotes the number of freezing pipes.

The vertical coordinate is determined by where corresponds to the vertical coordinate of the freezing pipes.

4.2. Pheromone Update and Behavior of Ants Seeking Path
4.2.1. Pheromone Update Mechanism

Following the evaluation of the previous solution of the freezing pipe layout by the artificial ants, the artificial ants deposit pheromones along their paths in the solution domain. The artificial ant performs tasks and seeks a path by identifying the pheromones in each path. Accordingly, the parameter is defined as the pheromone on node to a path at iteration . Initially, the amount of pheromone on each path is assumed to be equal in the solution domain, namely, . In the iteration, the artificial ants seek paths by identifying the pheromones on a node to a path and deposit new pheromones with the evaporation of the original pheromones. Thus, the pheromone update scheme in iteration is defined as follows: where denotes the pheromone evaporation ratio of pheromone in the range (0, 1) and represents the pheromone increment deposited by the artificial ant on node .

According to the pheromone update strategy, the ACO algorithm can be classified as an Ant-Cycle model [14], Ant-Density model, and Ant-Quantity model [13]. In these models, both Ant-Density and Ant-Quality use local information, whose search is not governed by any measures of the final results. However, the Ant-Cycle uses global information; that is, its ants lay an amount of trail proportional to the quality of the solution produced. Thus, both the Ant-Density and Ant-Quality models yield worse results than those obtained with the Ant-Cycle model [14]. The Ant-Cycle model was applied in this study and is described as follows: where denotes the pheromone strength with a constant value and represents the calculation result in iterations conducted by the artificial ant .

4.2.2. Behavior of Artificial Ants Seeking a Path

The artificial ants obtain the pheromone on a node to a path in the iteration and seek the corresponding node to the path after the probability calculation. The behavior can be described as follows: where denotes the probability that the artificial ant seeks the path node in the iteration, represents all the path nodes in step , indicates the path node sought by the artificial ant of step in iteration , and and denote the random numbers of (0, 1).

The path-seeking adjustment mechanism and the analysis of historical path information are introduced to further optimize the path-seeking behavior and to increase the diversity of path seeking and avoid the occurrence of the local optimal solution in the iteration process. The adjustment mechanism is described as follows: where rep represents the number of historical paths in the paths sought by artificial ants in the iteration, represents the path sought by the artificial ant in the iteration, denotes the collection of historical paths, denotes the path sought by m artificial ants in the t iteration, and rand is a random factor.

5. Optimization Result Analysis and Evaluation

The optimization parameters are . Table 4 lists the ranges of values for the layout parameters of the freezing pipes. The number of artificial ants selected was 10 as referred by Marwan et al. [1] and trial calculation. The termination criterion for the ACO algorithm is the maximum number of iterations, which was predefined according to the seepage flow rate. The convergence history is presented in Figure 6 to assess the performance of the ACO algorithm. The time required to obtain a completely frozen wall decreased significantly with an increase in iterations, and the iterations converged to an optimum solution after approximately 150 iterations. The ACO algorithm is relatively efficient in seeking the optimum solution for the layout of freezing pipes layout [1].

Figures 79 illustrate the correlation between the frozen wall and seepage flow at different freezing duration. It is evident from the figures that the closure time of frozen wall is shortened from 82.4 d to 56.9 d after layout optimization of the freezing pipes, and the freezing efficiency increases by 30.95%. Therefore, a layout of the freezing pipes determined by the ACO algorithm significantly reduced the freezing duration.

In Figure 7, the frozen columns are approximately formed in the designed boundary of the frozen wall and exhibit a more even distribution. The freezing duration is amended by the rapid development of the local frozen wall and the phenomenon of heat loss. When the freezing duration reached 40 d, the top arch and the bottom frozen wall were almost merged under the equidistant distribution of freezing pipes with weak zones, as shown in Figure 8. However, the frozen wall evenly developed with nonweak zones under the optimized layout. When the freezing duration reached 56.9 d, four weak zones still exist with nonuniform frozen walls under equidistant layout, requiring more freezing time for a fully frozen wall. Compared to the condition of equidistant layout, the completely frozen wall is uniformly merged with non-weak zones under the optimized layout, as shown in Figure 9. The deviation of the frozen wall from the designed position was significantly reduced. Meanwhile, the seepage flow rate decreased steadily in the upstream and downstream directions with the development of the frozen wall.

6. Conclusions

Seepage flow considerably influences the development of frozen walls. Considering the tunnel-connected aisle as an example, the numerical modelling is integrated within the ACO algorithm to determine the optimal layout of the freezing pipes. The conclusions drawn from the results are as follows. (1)In the presence of seepage flow, the closure time increased linearly, while the thickness decreased with an increase in the seepage flow(2)In the tunnel connected aisle, the development of the frozen wall presents inhomogeneity in different zones. Consequently, the thickness of the entire frozen wall was non-uniform in different zones(3)Four parameters were defined to parameterize the spatial distribution density of the freezing pipes. The coupled thermohydraulic model, together with the ACO algorithm, was applied to the optimization(4)The closure time of the frozen wall was shortened from 82.4 d to 56.9 d, and the freezing efficiency increased by 30.95% after optimization. The entire frozen wall was uniformly developed with a nonweak zone

Data Availability

The data that support the findings of this study are available from the corresponding author upon reasonable request, as the data also forms part of an ongoing study.

Conflicts of Interest

The authors declare that they have no conflicts of interest.

Acknowledgments

This research was funded by the National Natural Science Foundation of China (NSFC) (under grant Nos. 51978426 and 51708369) and Natural Science Foundation of Hebei Province (under grant Nos. E2021210007 and E2021210128).