# Temperature Field of Artificially Frozen Gravel Formation and Optimization of Freezing Pipe Layout Parameters under Seepage Flow

Authors: Tian-liang Wang; Hong-fang Song; Yu Shu; Fei Zhang; Ya-meng He
Journal: Geofluids
 (2022)

---

## 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.


---

## Body

## 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 ([4–6], [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.Table 1
Physical and thermal properties of gravel soils.
PropertiesValueUnitDry density1.883g/cm3Porosity0.3−Specific gravity2.69−Permeability coefficient60m/dFreezing temperature-0.04°CThermal conductivity coefficient1.07J/(kg·°C)Specific volume heat1.30m2/(s·°C)Particles larger than 2 mm80%Particles smaller than 0.075 mm0%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.Table 2
Merged time and thickness of frozen wall in the physical modeling.
Frozen wallSeepage flow rate (m/d)01.252.503.755.00Merged time (h)9.1712.6716.8321.1729.33Upstream thickness (mm)52362513−Downstream thickness (mm)5157616672Total thickness (mm)10393867972


### 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].
(1)Ceq∂T∂t+ρlL∂θl∂T+∇ρlCluT−λeq∇T=Qt,where Ceq denotes the equivalent volume heat capacity (m2/(s·°C)), T is the temperature (°C), ρl denotes the density of flowing water (kg/m3), L represents the latent heat of ice-water phase change (J/kg), θl indicates the volume content of flowing water, λeq denotes the equivalent heat conductivity coefficient (J/(kg·°C)), Cl corresponds to the specific heat capacity of flowing water (m2/(s·°C)), u indicates the velocity vector of flowing water (m/d), ∇ signifies the Hamiltonian operator, and Qt denotes the heat source.According to the basic properties of the porous medium, the parameters were identified as follows:(2)dzρCPeq∂T∂t+dzρCPu∇T+∇⟶q=dzQ+q0+dzQvd,⟶q=−dzλeq∇T,ρCPeq=θPρPCP,p+1−θPρCP,λeq=θPλP+1−θPλ+λdisp,where dz denotes the thickness of the gravel formation (m); CPandCP,p denote the specific capacity of ice-water mixture and soil particles, respectively (m2/(s·°C)); u is the velocity vector (m/d); ρandρP signify the density of ice-water mixture and soil particles, respectively (kg/m3); θp represents the volume fraction of soil particles in the total volume; λandλP denote the heat conductivity coefficient of ice-water mixture and soil particles, respectively (J/(kg·°C)); λdisp denotes the tensor of dispersion heat conductivity coefficient (J/(kg·°C)); and Q,Qvd,andq0 correspond to heat sources.


#### 2.2.3. Governing Equation of Seepage Flow

According to the assumptions in Section2.2.1, the seepage differential equation is identified as follows [1, 3]:
(3)∂∂tεPρ+∇ρu=Qm,∂∂tεPρ=ρεPχf+1−εPχP∂p∂t,u=−Kμd∇p,where εP denotes the porosity of the gravel formation; Qm represents the heat source; p signifies the pore water pressure (Pa); χfandχP denote the compressibility of water and soil particles, respectively (1/Pa); K indicates the permeability coefficient of gravel formation (m/d); and μd 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].
(4)ρ=θ1ρ1+θ2ρ2,CP=1ρθ1ρ1CP,1+θ2ρ2CP,2+L1⟶2∂αm∂T,αm=12θ2ρ2−θ1ρ1θ1ρ1+θ2ρ2,λ=θ1λ1+θ2λ2,θ1+θ2=1,where ρ1andρ2 denote the density of ice and water during phase changing (kg/m3), θ1andθ2 denotes the volume content of ice and water, CP,1andCP,2 represent the specific capacities of ice and water, respectively (m2/(s·°C)), L1⟶2 indicates the latent heat of ice-water phase change (J/kg), am denotes the equivalent heat conductivity coefficient of ice-water mixture (J/(kg·°C)), and λ1andλ2 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 Section2.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.Figure 1
Temperature field at closure time of frozen wall.
(a)
0.25 m/d (1.25 m/d)(b)
0.50 m/d (2.50 m/d)(c)
0.75 m/d (3.75 m/d)Table 3
Differences between results from the numerical and physical model.
Frozen wallSeepage flow rate0 m/d0.25 (1.25) m/d0.50 (2.50) m/d0.75 (3.75) m/d1.00 (5.00) m/dMerged time (%)0.812.503.247.8416.75Upstream thickness (%)0.329.688.9611.2015.21Downstream thickness (%)0.275.764.1223.1226.58Total thickness (%)0.307.696.6218.6921.36


## 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.Table 1
Physical and thermal properties of gravel soils.
PropertiesValueUnitDry density1.883g/cm3Porosity0.3−Specific gravity2.69−Permeability coefficient60m/dFreezing temperature-0.04°CThermal conductivity coefficient1.07J/(kg·°C)Specific volume heat1.30m2/(s·°C)Particles larger than 2 mm80%Particles smaller than 0.075 mm0%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.Table 2
Merged time and thickness of frozen wall in the physical modeling.
Frozen wallSeepage flow rate (m/d)01.252.503.755.00Merged time (h)9.1712.6716.8321.1729.33Upstream thickness (mm)52362513−Downstream thickness (mm)5157616672Total thickness (mm)10393867972


## 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].
(1)Ceq∂T∂t+ρlL∂θl∂T+∇ρlCluT−λeq∇T=Qt,where Ceq denotes the equivalent volume heat capacity (m2/(s·°C)), T is the temperature (°C), ρl denotes the density of flowing water (kg/m3), L represents the latent heat of ice-water phase change (J/kg), θl indicates the volume content of flowing water, λeq denotes the equivalent heat conductivity coefficient (J/(kg·°C)), Cl corresponds to the specific heat capacity of flowing water (m2/(s·°C)), u indicates the velocity vector of flowing water (m/d), ∇ signifies the Hamiltonian operator, and Qt denotes the heat source.According to the basic properties of the porous medium, the parameters were identified as follows:(2)dzρCPeq∂T∂t+dzρCPu∇T+∇⟶q=dzQ+q0+dzQvd,⟶q=−dzλeq∇T,ρCPeq=θPρPCP,p+1−θPρCP,λeq=θPλP+1−θPλ+λdisp,where dz denotes the thickness of the gravel formation (m); CPandCP,p denote the specific capacity of ice-water mixture and soil particles, respectively (m2/(s·°C)); u is the velocity vector (m/d); ρandρP signify the density of ice-water mixture and soil particles, respectively (kg/m3); θp represents the volume fraction of soil particles in the total volume; λandλP denote the heat conductivity coefficient of ice-water mixture and soil particles, respectively (J/(kg·°C)); λdisp denotes the tensor of dispersion heat conductivity coefficient (J/(kg·°C)); and Q,Qvd,andq0 correspond to heat sources.


### 2.2.3. Governing Equation of Seepage Flow

According to the assumptions in Section2.2.1, the seepage differential equation is identified as follows [1, 3]:
(3)∂∂tεPρ+∇ρu=Qm,∂∂tεPρ=ρεPχf+1−εPχP∂p∂t,u=−Kμd∇p,where εP denotes the porosity of the gravel formation; Qm represents the heat source; p signifies the pore water pressure (Pa); χfandχP denote the compressibility of water and soil particles, respectively (1/Pa); K indicates the permeability coefficient of gravel formation (m/d); and μd 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].
(4)ρ=θ1ρ1+θ2ρ2,CP=1ρθ1ρ1CP,1+θ2ρ2CP,2+L1⟶2∂αm∂T,αm=12θ2ρ2−θ1ρ1θ1ρ1+θ2ρ2,λ=θ1λ1+θ2λ2,θ1+θ2=1,where ρ1andρ2 denote the density of ice and water during phase changing (kg/m3), θ1andθ2 denotes the volume content of ice and water, CP,1andCP,2 represent the specific capacities of ice and water, respectively (m2/(s·°C)), L1⟶2 indicates the latent heat of ice-water phase change (J/kg), am denotes the equivalent heat conductivity coefficient of ice-water mixture (J/(kg·°C)), and λ1andλ2 are the heat conductivity coefficients of water and ice, respectively (J/(kg·°C)).


## 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].
(1)Ceq∂T∂t+ρlL∂θl∂T+∇ρlCluT−λeq∇T=Qt,where Ceq denotes the equivalent volume heat capacity (m2/(s·°C)), T is the temperature (°C), ρl denotes the density of flowing water (kg/m3), L represents the latent heat of ice-water phase change (J/kg), θl indicates the volume content of flowing water, λeq denotes the equivalent heat conductivity coefficient (J/(kg·°C)), Cl corresponds to the specific heat capacity of flowing water (m2/(s·°C)), u indicates the velocity vector of flowing water (m/d), ∇ signifies the Hamiltonian operator, and Qt denotes the heat source.According to the basic properties of the porous medium, the parameters were identified as follows:(2)dzρCPeq∂T∂t+dzρCPu∇T+∇⟶q=dzQ+q0+dzQvd,⟶q=−dzλeq∇T,ρCPeq=θPρPCP,p+1−θPρCP,λeq=θPλP+1−θPλ+λdisp,where dz denotes the thickness of the gravel formation (m); CPandCP,p denote the specific capacity of ice-water mixture and soil particles, respectively (m2/(s·°C)); u is the velocity vector (m/d); ρandρP signify the density of ice-water mixture and soil particles, respectively (kg/m3); θp represents the volume fraction of soil particles in the total volume; λandλP denote the heat conductivity coefficient of ice-water mixture and soil particles, respectively (J/(kg·°C)); λdisp denotes the tensor of dispersion heat conductivity coefficient (J/(kg·°C)); and Q,Qvd,andq0 correspond to heat sources.


## 2.2.3. Governing Equation of Seepage Flow

According to the assumptions in Section2.2.1, the seepage differential equation is identified as follows [1, 3]:
(3)∂∂tεPρ+∇ρu=Qm,∂∂tεPρ=ρεPχf+1−εPχP∂p∂t,u=−Kμd∇p,where εP denotes the porosity of the gravel formation; Qm represents the heat source; p signifies the pore water pressure (Pa); χfandχP denote the compressibility of water and soil particles, respectively (1/Pa); K indicates the permeability coefficient of gravel formation (m/d); and μd 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].
(4)ρ=θ1ρ1+θ2ρ2,CP=1ρθ1ρ1CP,1+θ2ρ2CP,2+L1⟶2∂αm∂T,αm=12θ2ρ2−θ1ρ1θ1ρ1+θ2ρ2,λ=θ1λ1+θ2λ2,θ1+θ2=1,where ρ1andρ2 denote the density of ice and water during phase changing (kg/m3), θ1andθ2 denotes the volume content of ice and water, CP,1andCP,2 represent the specific capacities of ice and water, respectively (m2/(s·°C)), L1⟶2 indicates the latent heat of ice-water phase change (J/kg), am denotes the equivalent heat conductivity coefficient of ice-water mixture (J/(kg·°C)), and λ1andλ2 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 Section2.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.Figure 1
Temperature field at closure time of frozen wall.
(a)
0.25 m/d (1.25 m/d)(b)
0.50 m/d (2.50 m/d)(c)
0.75 m/d (3.75 m/d)Table 3
Differences between results from the numerical and physical model.
Frozen wallSeepage flow rate0 m/d0.25 (1.25) m/d0.50 (2.50) m/d0.75 (3.75) m/d1.00 (5.00) m/dMerged time (%)0.812.503.247.8416.75Upstream thickness (%)0.329.688.9611.2015.21Downstream thickness (%)0.275.764.1223.1226.58Total thickness (%)0.307.696.6218.6921.36


## 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 Table1. 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.Figure 2
Dimensions of the gravel formation and tunnel connected aisle.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 Figure3. Gravel formation around the freezing pipes was established as boundary layer meshing to fit the heat diffusion mode of the freezing pipes.Figure 3
Numerical simulation model and its meshing.


### 3.2. Frozen Wall Merging of Tunnel Connected Aisle

Figure4 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.Figure 4
Gravel formation temperature field at different freezing duration.
(a)
20 d(b)
40 d(c)
60 d(d)
80 dInitially, the flowing water almost homogeneously flowed passed the whole cross-section (see Figure4(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.


## 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 Table1. 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.Figure 2
Dimensions of the gravel formation and tunnel connected aisle.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 Figure3. Gravel formation around the freezing pipes was established as boundary layer meshing to fit the heat diffusion mode of the freezing pipes.Figure 3
Numerical simulation model and its meshing.


## 3.2. Frozen Wall Merging of Tunnel Connected Aisle

Figure4 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.Figure 4
Gravel formation temperature field at different freezing duration.
(a)
20 d(b)
40 d(c)
60 d(d)
80 dInitially, the flowing water almost homogeneously flowed passed the whole cross-section (see Figure4(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 Figure5. 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, d1andd2, 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.Figure 5
Frozen wall edges at different freezing duration.
(a)
40 d(b)
50 d


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

The horizontal coordinate is determined by(5)xi_u/d=xu/d−δxi_u/d,i_u/d=1,2,⋯,N_u/d,where the symbols u/d represent the upstream and downstream zones, respectively, in the following parameters: xi_u and xi_d correspond to the horizontal coordinates of the freezing pipes, xuandxd denote the horizontal coordinates of the frozen wall center, δxi_u and δxi_d represent the horizontal coordinate offsets of the freezing pipes, i_u and i_d denote the labels of the freezing pipes, N_u and N_d correspond to the numbers of the freezing pipes, and δxi_u/d is defined as
(6)δxi_u/d=d1_u/d,ifd1_u/d<d1_u/d_max,F_th+d2_u/d,else, where d1_u,d1_d,d2_u,andd2_d denote the offsets of the freezing pipes, the ranges of which depend on the designed thickness of the frozen wall F_th.The vertical coordinate is determined by(7)yi_u/dyb,ifi_u/d=1,yi_u/d−1+δyi_u/d−1,else,where yi_uandyi_d denote the vertical coordinates of the freezing pipes, δyi_uandδyi_d denote the vertical coordinate offsets of the freezing pipes, yb represents the vertical coordinate of the bottom frozen wall center, and δyi_u/d−1 is defined as
(8)δyi_u/d−1=yt−ybGi_u/d−1,u_u/d,σ_u/d∑j=1N−1j,u_u/d,σ_u/d,where μ_u and μ_d denote the mean of the freezing pipes and σ_u and σ_d denote the standard deviations of the freezing pipes with GX,μ,σ representing the Gaussian distribution.
(9)GX,μ,σ=1σ2πexp−X−μ/2σ2.


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

The horizontal coordinate is determined by(10)xi_t=xd+δRθi_t×cosθi_t,i_t=1,2,⋯,N_t.The vertical coordinate is determined by(11)yi_t=yt+xd+δRθi_t×sinθi_t,i_t=1,2,⋯,N_t,where i_t denotes the label of the freezing pipes, yt represents the vertical coordinate of the top arch frozen wall center, N_t denotes the number of freezing pipes, and the radial coordinate offsets of the freezing pipes δRθi_t is defined as
(12)δRθi_t=d1_t+d2_t2θi_tπ−d1_t,ifθi_t<π2,d2_t,ifθi_t³π2,where d1_tandd2_t denote the offsets of the freezing pipes, the ranges of which depend on the designed thickness of the frozen wall F_th. The angular spacing function of the freezing pipes θi_t is defined as
(13)θi_t=0,ifi_t=1,θi_t−1+δθi_t−1,ifi_t>1,(14)δθi_t−1=πGi_t−1,μ_t,σ_t∑j=1N−1Gj,μ_t,σ_t,where μ_t denotes the mean of the freezing pipes and σ_t corresponds to the standard deviation of the freezing pipes with GX,μ,σ 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(15)xi_b=xd,ifi_b=1,xi_b−1−δxi_b−1,if1<i_b≤N_b,(16)δxi_b−1=xu−xdGi_b−1,μ_b,σ_b∑j=1N−1Gj,μ_b,σ_b,where xi_b denotes the horizontal coordinate of the freezing pipes, μ_b indicates the mean of the freezing pipes, σ_b represents the standard deviation of the freezing pipes with GX,μ,σ representing the Gaussian distribution in Equation (9), and N_b denotes the number of freezing pipes.The vertical coordinate is determined by(17)yi_b=yb,where yi_b 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 antkk=1,2,⋯,m performs tasks and seeks a path by identifying the pheromones in each path. Accordingly, the parameter τijt is defined as the pheromone on node Cij to a path at iteration t. Initially, the amount of pheromone on each path is assumed to be equal in the solution domain, namely, τij0=const. In the t iteration, the m artificial ants seek paths by identifying the pheromones τijt on a node Cij to a path and deposit new pheromones Δτijt with the evaporation of the original pheromones. Thus, the pheromone update scheme in iteration t+1 is defined as follows:
(18)τijt+1=1−γ∙τijt+∆τijt,∆τijt=∑k=1m∆τijkt,where γ denotes the pheromone evaporation ratio of pheromone in the range (0, 1) and ∆τijkt represents the pheromone increment deposited by the artificial ant k on node Cij.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:
(19)Δτijkt=QEtk,iftheartificialantkpassesbythenodeCijiniterationt,0,else,where Q denotes the pheromone strength with a constant value and Etk represents the calculation result in t iterations conducted by the artificial ant k.


#### 4.2.2. Behavior of Artificial Ants Seeking a Path

The artificial ants obtain the pheromone on a node to a path in thet iteration and seek the corresponding node to the path after the probability calculation. The behavior can be described as follows:
(20)pijkt=τijt∑1sizeMiτijt,S_Cikt=arglargepi1kt:pijkt,1,ifr<β,arglargepi1kt:pijkt,2,else,where pijkt denotes the probability that the artificial ant k seeks the path node Cij in the t iteration, sizeMi represents all the path nodes in step i, S_Cikt indicates the path node sought by the artificial ant k of step i in iteration t, and r 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:(21)rep=sumtabukt∈tabu_h,tabut=rand,ifrep≥m−1m×m,tabut,else,where rep represents the number of historical paths in the paths sought by m artificial ants in the t iteration, tabukt represents the path sought by the artificial ant k in the t iteration, tabu_h denotes the collection of historical paths, tabut denotes the path sought by m artificial ants in the t iteration, and rand is a random factor.


## 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 Figure5. 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, d1andd2, 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.Figure 5
Frozen wall edges at different freezing duration.
(a)
40 d(b)
50 d


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

The horizontal coordinate is determined by(5)xi_u/d=xu/d−δxi_u/d,i_u/d=1,2,⋯,N_u/d,where the symbols u/d represent the upstream and downstream zones, respectively, in the following parameters: xi_u and xi_d correspond to the horizontal coordinates of the freezing pipes, xuandxd denote the horizontal coordinates of the frozen wall center, δxi_u and δxi_d represent the horizontal coordinate offsets of the freezing pipes, i_u and i_d denote the labels of the freezing pipes, N_u and N_d correspond to the numbers of the freezing pipes, and δxi_u/d is defined as
(6)δxi_u/d=d1_u/d,ifd1_u/d<d1_u/d_max,F_th+d2_u/d,else, where d1_u,d1_d,d2_u,andd2_d denote the offsets of the freezing pipes, the ranges of which depend on the designed thickness of the frozen wall F_th.The vertical coordinate is determined by(7)yi_u/dyb,ifi_u/d=1,yi_u/d−1+δyi_u/d−1,else,where yi_uandyi_d denote the vertical coordinates of the freezing pipes, δyi_uandδyi_d denote the vertical coordinate offsets of the freezing pipes, yb represents the vertical coordinate of the bottom frozen wall center, and δyi_u/d−1 is defined as
(8)δyi_u/d−1=yt−ybGi_u/d−1,u_u/d,σ_u/d∑j=1N−1j,u_u/d,σ_u/d,where μ_u and μ_d denote the mean of the freezing pipes and σ_u and σ_d denote the standard deviations of the freezing pipes with GX,μ,σ representing the Gaussian distribution.
(9)GX,μ,σ=1σ2πexp−X−μ/2σ2.


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

The horizontal coordinate is determined by(10)xi_t=xd+δRθi_t×cosθi_t,i_t=1,2,⋯,N_t.The vertical coordinate is determined by(11)yi_t=yt+xd+δRθi_t×sinθi_t,i_t=1,2,⋯,N_t,where i_t denotes the label of the freezing pipes, yt represents the vertical coordinate of the top arch frozen wall center, N_t denotes the number of freezing pipes, and the radial coordinate offsets of the freezing pipes δRθi_t is defined as
(12)δRθi_t=d1_t+d2_t2θi_tπ−d1_t,ifθi_t<π2,d2_t,ifθi_t³π2,where d1_tandd2_t denote the offsets of the freezing pipes, the ranges of which depend on the designed thickness of the frozen wall F_th. The angular spacing function of the freezing pipes θi_t is defined as
(13)θi_t=0,ifi_t=1,θi_t−1+δθi_t−1,ifi_t>1,(14)δθi_t−1=πGi_t−1,μ_t,σ_t∑j=1N−1Gj,μ_t,σ_t,where μ_t denotes the mean of the freezing pipes and σ_t corresponds to the standard deviation of the freezing pipes with GX,μ,σ 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(15)xi_b=xd,ifi_b=1,xi_b−1−δxi_b−1,if1<i_b≤N_b,(16)δxi_b−1=xu−xdGi_b−1,μ_b,σ_b∑j=1N−1Gj,μ_b,σ_b,where xi_b denotes the horizontal coordinate of the freezing pipes, μ_b indicates the mean of the freezing pipes, σ_b represents the standard deviation of the freezing pipes with GX,μ,σ representing the Gaussian distribution in Equation (9), and N_b denotes the number of freezing pipes.The vertical coordinate is determined by(17)yi_b=yb,where yi_b corresponds to the vertical coordinate of the freezing pipes.


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

The horizontal coordinate is determined by(5)xi_u/d=xu/d−δxi_u/d,i_u/d=1,2,⋯,N_u/d,where the symbols u/d represent the upstream and downstream zones, respectively, in the following parameters: xi_u and xi_d correspond to the horizontal coordinates of the freezing pipes, xuandxd denote the horizontal coordinates of the frozen wall center, δxi_u and δxi_d represent the horizontal coordinate offsets of the freezing pipes, i_u and i_d denote the labels of the freezing pipes, N_u and N_d correspond to the numbers of the freezing pipes, and δxi_u/d is defined as
(6)δxi_u/d=d1_u/d,ifd1_u/d<d1_u/d_max,F_th+d2_u/d,else, where d1_u,d1_d,d2_u,andd2_d denote the offsets of the freezing pipes, the ranges of which depend on the designed thickness of the frozen wall F_th.The vertical coordinate is determined by(7)yi_u/dyb,ifi_u/d=1,yi_u/d−1+δyi_u/d−1,else,where yi_uandyi_d denote the vertical coordinates of the freezing pipes, δyi_uandδyi_d denote the vertical coordinate offsets of the freezing pipes, yb represents the vertical coordinate of the bottom frozen wall center, and δyi_u/d−1 is defined as
(8)δyi_u/d−1=yt−ybGi_u/d−1,u_u/d,σ_u/d∑j=1N−1j,u_u/d,σ_u/d,where μ_u and μ_d denote the mean of the freezing pipes and σ_u and σ_d denote the standard deviations of the freezing pipes with GX,μ,σ representing the Gaussian distribution.
(9)GX,μ,σ=1σ2πexp−X−μ/2σ2.


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

The horizontal coordinate is determined by(10)xi_t=xd+δRθi_t×cosθi_t,i_t=1,2,⋯,N_t.The vertical coordinate is determined by(11)yi_t=yt+xd+δRθi_t×sinθi_t,i_t=1,2,⋯,N_t,where i_t denotes the label of the freezing pipes, yt represents the vertical coordinate of the top arch frozen wall center, N_t denotes the number of freezing pipes, and the radial coordinate offsets of the freezing pipes δRθi_t is defined as
(12)δRθi_t=d1_t+d2_t2θi_tπ−d1_t,ifθi_t<π2,d2_t,ifθi_t³π2,where d1_tandd2_t denote the offsets of the freezing pipes, the ranges of which depend on the designed thickness of the frozen wall F_th. The angular spacing function of the freezing pipes θi_t is defined as
(13)θi_t=0,ifi_t=1,θi_t−1+δθi_t−1,ifi_t>1,(14)δθi_t−1=πGi_t−1,μ_t,σ_t∑j=1N−1Gj,μ_t,σ_t,where μ_t denotes the mean of the freezing pipes and σ_t corresponds to the standard deviation of the freezing pipes with GX,μ,σ 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(15)xi_b=xd,ifi_b=1,xi_b−1−δxi_b−1,if1<i_b≤N_b,(16)δxi_b−1=xu−xdGi_b−1,μ_b,σ_b∑j=1N−1Gj,μ_b,σ_b,where xi_b denotes the horizontal coordinate of the freezing pipes, μ_b indicates the mean of the freezing pipes, σ_b represents the standard deviation of the freezing pipes with GX,μ,σ representing the Gaussian distribution in Equation (9), and N_b denotes the number of freezing pipes.The vertical coordinate is determined by(17)yi_b=yb,where yi_b 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 antkk=1,2,⋯,m performs tasks and seeks a path by identifying the pheromones in each path. Accordingly, the parameter τijt is defined as the pheromone on node Cij to a path at iteration t. Initially, the amount of pheromone on each path is assumed to be equal in the solution domain, namely, τij0=const. In the t iteration, the m artificial ants seek paths by identifying the pheromones τijt on a node Cij to a path and deposit new pheromones Δτijt with the evaporation of the original pheromones. Thus, the pheromone update scheme in iteration t+1 is defined as follows:
(18)τijt+1=1−γ∙τijt+∆τijt,∆τijt=∑k=1m∆τijkt,where γ denotes the pheromone evaporation ratio of pheromone in the range (0, 1) and ∆τijkt represents the pheromone increment deposited by the artificial ant k on node Cij.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:
(19)Δτijkt=QEtk,iftheartificialantkpassesbythenodeCijiniterationt,0,else,where Q denotes the pheromone strength with a constant value and Etk represents the calculation result in t iterations conducted by the artificial ant k.


### 4.2.2. Behavior of Artificial Ants Seeking a Path

The artificial ants obtain the pheromone on a node to a path in thet iteration and seek the corresponding node to the path after the probability calculation. The behavior can be described as follows:
(20)pijkt=τijt∑1sizeMiτijt,S_Cikt=arglargepi1kt:pijkt,1,ifr<β,arglargepi1kt:pijkt,2,else,where pijkt denotes the probability that the artificial ant k seeks the path node Cij in the t iteration, sizeMi represents all the path nodes in step i, S_Cikt indicates the path node sought by the artificial ant k of step i in iteration t, and r 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:(21)rep=sumtabukt∈tabu_h,tabut=rand,ifrep≥m−1m×m,tabut,else,where rep represents the number of historical paths in the paths sought by m artificial ants in the t iteration, tabukt represents the path sought by the artificial ant k in the t iteration, tabu_h denotes the collection of historical paths, tabut denotes the path sought by m artificial ants in the t iteration, and rand is a random factor.


## 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 antkk=1,2,⋯,m performs tasks and seeks a path by identifying the pheromones in each path. Accordingly, the parameter τijt is defined as the pheromone on node Cij to a path at iteration t. Initially, the amount of pheromone on each path is assumed to be equal in the solution domain, namely, τij0=const. In the t iteration, the m artificial ants seek paths by identifying the pheromones τijt on a node Cij to a path and deposit new pheromones Δτijt with the evaporation of the original pheromones. Thus, the pheromone update scheme in iteration t+1 is defined as follows:
(18)τijt+1=1−γ∙τijt+∆τijt,∆τijt=∑k=1m∆τijkt,where γ denotes the pheromone evaporation ratio of pheromone in the range (0, 1) and ∆τijkt represents the pheromone increment deposited by the artificial ant k on node Cij.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:
(19)Δτijkt=QEtk,iftheartificialantkpassesbythenodeCijiniterationt,0,else,where Q denotes the pheromone strength with a constant value and Etk represents the calculation result in t iterations conducted by the artificial ant k.


## 4.2.2. Behavior of Artificial Ants Seeking a Path

The artificial ants obtain the pheromone on a node to a path in thet iteration and seek the corresponding node to the path after the probability calculation. The behavior can be described as follows:
(20)pijkt=τijt∑1sizeMiτijt,S_Cikt=arglargepi1kt:pijkt,1,ifr<β,arglargepi1kt:pijkt,2,else,where pijkt denotes the probability that the artificial ant k seeks the path node Cij in the t iteration, sizeMi represents all the path nodes in step i, S_Cikt indicates the path node sought by the artificial ant k of step i in iteration t, and r 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:(21)rep=sumtabukt∈tabu_h,tabut=rand,ifrep≥m−1m×m,tabut,else,where rep represents the number of historical paths in the paths sought by m artificial ants in the t iteration, tabukt represents the path sought by the artificial ant k in the t iteration, tabu_h denotes the collection of historical paths, tabut 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 ared1,d2,μ,andσ. 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].Table 4
Rang of values for freezing pipe layout parameters.
ParametersThe zone of the freezing pipes located in the tunnel connected aisleUpstream zone (_u)Downstream zone (_d)Top arch zone (_t)Bottom zone (_b)d10,d_p/2,F_th/2−d_p/20,d_p/2,F_th/2−d_p/20,d_p/2,F_th/2−d_p/2—d20,d_p/2,F_th/40,d_p/2,F_th/40,d_p/2,F_th/2−d_p/2—μ1,0.25,N_u1,0.25,N_d1,0.25,N_t+21,0.25,N_b+2σ1,0.25,N_u1,0.25,N_d1,0.25,N_t+21,0.25,N_b+2Note: whered_p denotes the diameter of the freezing pipe.Figure 6
Iteration process of AOC algorithm.Figures7–9 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.Figure 7
Frozen wall and seepage flow at freezing duration of 20 d.
(a)
Equidistant distribution(b)
Optimized layoutFigure 8
Frozen wall and seepage flow at freezing duration of 40 d.
(a)
Equidistant distribution(b)
Optimized layoutFigure 9
Frozen wall and seepage flow at freezing duration of 56.9 d.
(a)
Equidistant distribution(b)
Optimized layoutIn Figure7, 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



---

*Source: 1597645-2022-03-24.xml*