Lattice Boltzmann Simulation for the Forming Process of Artificial Frozen Soil Wall

A lattice Boltzmann model is proposed to simulate the forming process of artificial frozen soil wall. The enthalpy method is applied to deal with the latent-heat source term, and the adjustable thermal diffusivity is utilized to handle the change of thermophysical parameters. The model is tested by the heat conduction with solid–liquid phase change in semi-infinite space, which shows a good consistence between the numerical and analytical solutions, and the mesh resolution has little effect on the numerical results. Lastly, the development of frozen soil wall is discussed when the freezing pipes are arranged in a square. The results show that the evolution of temperature field with time is closely related to the distance from the freezing pipe. For the soil near freezing pipe, the temperature gradient is larger, the soil temperature drops rapidly and freezes in a short time. The time history curve of temperature is relatively smooth. For the soil far away from freezing pipe, the temperature evolution curve has obvious multistage, which can be divided into four stages: cooling, phase change, partly frozen and completely frozen. The spacing of freezing pipes has a significant influence on the overlapping time of artificial frozen soil wall, and there is a power function relationship between them.


Introduction
As an effective temporary ground improvement technique, artificial ground freezing (AGF) has been widely adopted in geotechnical engineering [1,2], including departure and reception of shield tunnel, tunnels connecting passage in Metro, mine shaft sinking and municipal engineering, etc.It has advantages of strong stratum adaptability, good sealing performance, high strength and little influence on the surrounding environment.In ground improvement engineering, the development of the frozen soil wall determines its average temperature, thickness, physical and mechanical properties.These indexes reflect the strength and stability of the frozen soil wall, which is directly related to the scheduling management of projects.In the process of artificial freezing, the temperature evolution of soil is a transient heat conduction problem including ice-water phase change, latent heat release, internal heat source, moving boundary and irregular geometric boundary.The soil temperature distribution is also affected by the interaction of freezing pipes.Therefore, the forming process of frozen soil wall is very complicated, and a better understanding of heat transfer mechanisms is essentially important, which can provide necessary technical guarantees for the implementation of artificial freezing projects.
Based on the steady heat conduction theory, a large number of scholars have proposed the analytical methods of temperature evolution with a single pipe, single or double row pipes [3][4][5].
Energies 2019, 12, 46 2 of 15 However, these methods are limited to simplified boundary conditions and idealized initial conditions.Besides, the theoretical formula is too complicated for engineering application.In recent years, numerical methods have been widely applied in heat conduction problems with phase change.Singh [6] studied the flow and heat transfer characteristics of a phase transition, melting problem.Santos [7] applied the finite element method to predict freezing times of mushrooms.Farrokhpanah [8] introduced a new smoothed particle hydrodynamics (SPH) method to model the heat transfer with phase change considering the latent heat released (absorbed) during solidification (melting).Furenes [9] used the event location algorithm in the finite difference method for phase-change problems.
When compared with the traditional numerical methods, the lattice Boltzmann method (LBM) enjoys advantages of both macroscopic and microscopic approaches [10,11].It has clear physical conception, easy programming, high computational efficiency and is easy to apply for complex domains [12][13][14].So, LBM has been explored to deal with heat conduction problems with phase change.Miller [15] proposed a simple model for the liquid-solid phase change based on the lattice Boltzmann method with enhanced collisions.Jiaung [16] firstly developed an enthalpy-based lattice Boltzmann model for simulating solid/liquid phase change problem governed by the heat conduction equation.Huber et al. [17][18][19] improved this model, and used it to couple thermal convection and phase change of single-component systems.Eshraghi [20] developed a new variation to solve the heat conduction with phase change by treating implicitly the latent heat source term.Huang [21] proposed a new lattice Boltzmann model to treat the latent heat source term by modifying the equilibrium distribution function.Sadeghi [22,23] proposed a three-dimensional Boltzmann model to study the film-boiling phenomenon.Chatterjee [24] extended the lattice Boltzmann formulation to simulate three-dimensional heat diffusion coupled with solid-liquid phase change.Li [25] presented a three-dimensional multiple-relaxation-time lattice Boltzmann model for the solid-liquid phase change based on the enthalpy conservation equation.The above studies mainly neglected the change of thermal diffusivity for simplifying the calculation.However, the thermal diffusivity of liquid water is only 1/9 of ice, and the water content of artificial frozen soil is generally high, so it is necessary to consider the change of thermal diffusivity during the forming process of artificial frozen soil wall.
In this paper, the enthalpy approach is applied to treat the latent heat source term in the energy equation, the adjustable thermal diffusivity is utilized to simulate the change of thermophysical parameters, and a thermal lattice Boltzmann model is proposed to simulate the forming process of artificial frozen soil wall.The model is subsequently tested with the solid-liquid phase change of pure substance in semi-infinite space.Finally, the forming process of an artificial frozen wall is simulated when the four freezing pipes are arranged in a square, the evolution of the freezing front and the temperature distribution are analyzed during the artificial freezing process, which provides the theoretical basis for the design and construction of practical engineering.

Assumptions
To develop the heat conduction model with phase change for soil, the following assumptions are made: (1) Ignoring the influence of pore structure, soil is regard as the continuous, homogeneous and isotropic medium.
(2) According to the physical state of water, the freezing soil is divided into two parts: solid and liquid zones, and the thermophysical parameters are constant in each zone.
(3) The freezing temperature of soil is constant, and the liquid phase fraction is applied to trace the solid-liquid interface.
(4) The coupling of temperature, stress and moisture is neglected during the artificial freezing process.

Mathematical Model
Compared with the three-dimensional model, the two-dimensional one could reflect the general evolution laws and save a lot of computing resources, so it is selected in the present paper.According to energy conservation, the mathematical model of soil freezing can be expressed as [26]: Liquid phase The solid-liquid interface where T is the temperature, t is the time.ρ,k and α are the density, thermal conductivity and thermal diffusivity, respectively, and all these parameters can be determined by the volume fraction of each component in soil [27,28].L a is the latent heat, S is the position of solid-liquid interface.r is the normal direction of frozen front.The subscript l represents that the water in the soil is the liquid phase, and s the solid phase.

Phase Change Treatment
The difficulties of temperature prediction lie in the treatment of latent heat and the movement of solid-liquid interface during the soil freezing process.In this paper, the enthalpy model proposed by Shamsundar [29] is adopt to develop the unified energy equation in the whole region (including liquid zone, solid zone and solid-liquid interface), and the solid-liquid interface is determined by solving the enthalpy parameter.This model does not need to separate the solid and liquid phases, and trace the solid-liquid interface.Mathematically, it has been proved to be equivalent to the heat conduction equation with phase change [30].So, the Equations ( 1)-(3) can be unified as where H is the total enthalpy, which can be divided into sensible and latent enthalpy components as: where C p is the specific heat capacity, ϕ is the liquid phase fraction, which is 0 for solid zone, 1 for liquid zone.
When the freezing temperature T f is constant, the relationship between the liquid phase fraction ϕ and the total enthalpy H can be expressed as the follows: Substituting Equation (5) into Equation (4) yields: Energies 2019, 12, 46 4 of 15 Equation ( 7) is applied to describe the evolution of temperature field during the soil freezing process.In terms of the liquid phase fraction ϕ, the thermal diffusivity α can be expressed as:

Dimensionless Treatment
To facilitate transformation between physical and lattice units, the following dimensionless parameters are introduced where X, θ are the dimensionless coordinate and temperature.Ste, F o are the Stefan number and Fourier number, respectively.L is the reference length, T 0 is the temperature of cold source, T i is the initial temperature of soil.

Lattice Boltzmann Equation
The d-dimensional m-speed (DdQm) model proposed by Qian et al. [31], is the basic model of LBM.Compared with D2Q9 model, the D2Q4 model (Figure 1) has not only comparable accuracy but also better computational efficiency for temperature evolution, so it is employed to solve the heat conduction equation, which is governed by Equation (7).And the enthalpy model is adopted to treat the latent heat source term.The discrete form of the lattice Boltzmann equation can be written as: where g i (r, t) is the temperature distribution function of the ith direction at the lattice site r and time t, g eq i (r, t) represents the equilibrium distribution function, τ is the dimensionless relaxation time, whose value should insure to be within (0.5,2) [32], e i is the discrete velocity in the lattice, which is composed of the velocity vectors: where c is the lattice speed, and c = δ x /δ t , δ x , δ t are the lattice space and time step, respectively, and they are usually taken as δ x = δ t = 1 for simplifying the calculation.
Energies 2018, 11, x FOR PEER REVIEW 4 of 16 Equation ( 7) is applied to describe the evolution of temperature field during the soil freezing process.In terms of the liquid phase fraction φ, the thermal diffusivity α can be expressed as:

Dimensionless Treatment
To facilitate transformation between physical and lattice units, the following dimensionless parameters are introduced ( ) where X , θ are the dimensionless coordinate and temperature.Ste , o F are the Stefan number and Fourier number, respectively.L is the reference length, T 0 is the temperature of cold source, i T is the initial temperature of soil.The d-dimensional m-speed (DdQm) model proposed by Qian et al. [31], is the basic model of LBM.Compared with D2Q9 model, the D2Q4 model (Figure 1) has not only comparable accuracy but also better computational efficiency for temperature evolution, so it is employed to solve the heat conduction equation, which is governed by Equation (7).And the enthalpy model is adopted to treat the latent heat source term.The discrete form of the lattice Boltzmann equation can be written as:

Lattice Boltzmann Equation
where ( ) r is the temperature distribution function of the ith direction at the lattice site r and g eq i (r, t) can be described as where ω i is the weight factor in the ith direction, and Sr is the heat source term, according to Equation ( 7), it can be expressed as: Based on the first-order expansion, the heat source term Sr is approximately discretized into: By Chapman-Enskog expansion, Equation ( 10) can be recovered to the macroscopic heat conduction Equation ( 7), the thermal diffusivity α is given by where c s is lattice sound speed, for the D2Q4 model, c 2 s = c 2 /2.The macroscopic temperature can be calculated as:

Boundary Conditions
The adiabatic boundary involved in the present study, is handled by the non-equilibrium extrapolation scheme, proposed by Guo et al. [33] in 2002, which has second order accuracy.The main idea of this approach is to divide the temperature distribution functions at the boundary node N B into its equilibrium and non-equilibrium parts, The equilibrium part g eq i (N B , t) can be got by Equation (12), and the non-equilibrium part g neq i (N B , t) is approximated by extrapolating from the neighboring node N O ,

Unit Conversion
For numerical methods, it is necessary to achieve the unit conversion between physical and lattice units.In the present study, the dimensionless treatment is carried out for all the parameters, which ensures the consistency of heat transfer criterion.The non-dimensional numbers such as Stefan number and Fourier number, are used as a bridge to realize the conversion between two unit systems.
For the thermophysical parameters (latent heat L a , specific heat C p ), the unit conversion can be achieved based on the Stefan number Ste.
where subscript p represents the physical unit, and L the lattice unit.
According to an artificial freezing project, the sandy silt [34] is taken as the research object.The unit conversion of thermophysical parameters is handled, and the comparison between physical and lattice units is shown in Table 1.For the physical time t p and lattice time steps N, the relationship between them can be established according to the Fourier number F o .
where t p and t L are the time in physical and lattice unit, t L = Nδ t .L p and L L are the reference length in physical and lattice units, if L p is the length of calculation domain, then L L = nδ x , n is the corresponding number of lattices.
If the physical model of 4 m × 4 m is divided into a lattice of 1000 × 1000 grid cells, and the values of α p and α L are assigned according to Table 1.The relationship between t p and t L can be deduced as It can be known that if δ t is selected as 1, the physical time 1 s is corresponding to 3.35-time steps and 1 day to 25,790-time steps.

Flowchart of Program Realization
Considering the effects of heat transfer, latent heat and movement of phase-change interface during the soil freezing process, the lattice Boltzmann model is proposed based on the enthalpy method, the corresponding flowchart is shown in Figure 2.

Verification
To verify accuracy of proposed model, LBM is applied to simulate the heat conduction problem with solid-liquid phase change in semi-infinite space, as shown in Figure 3.At the time t = 0, the uniform initial temperature is T i , and the substance is in the liquid state.The cold source temperature T 0 is set at the position of x = 0, and keeps constant at t > 0. The analytical solution of the solid-liquid interface S(t) and temperature T(x, t) are as follows [26]   Solid-liquid interface Solid-liquid interface

Verification
Solid phase Liquid phase where λ is the unknown parameter, which can be obtained by the following transcendental equation.
In this case, the entire domain has a size of L × H = 2 m × 0.16 m, and it is discretized using 125 × 10, 250 × 20 and 500 × 40 grid cells with lattice resolution of 16 mm, 8 mm and 4 mm, respectively.The thermophysical parameters are set according to Table 1.The temperature is fixed at T 0 on the left side, and the right boundary is adiabatic.The lattice spaces are equal in horizontal and vertical direction, δ x = δ y = 1.0 and the time step is δ t = 1.0.
Figure 4 shows the evolution of solid-liquid interface when L × H = 125 × 10, 250 × 20 and 500 × 40, respectively.It can be seen that the results in the present study are in good agreement with the analytical ones, and the mesh resolution has little effect on the numerical results.So the proposed model can accurately simulate the movement of solid-liquid interface during the freezing process.The temperature and error distribution are shown in Figure 5 after freezing 10 days.The slope change of temperature distribution is observed at the solid-liquid interface, which is induced by the variety of thermal diffusivity.A good consistence can also be seen between the numerical results and analytical ones, which indicates the validity of the proposed model in handling the heat conduction problem with phase change.The errors are defined as the numerical solutions of temperature minus the analytical ones.As shown in Figure 5b, the errors fluctuate near the solid-liquid interface, the finer grid resolution has the smaller error, and the farther the distance from the interface, the smaller the error.The maximum errors at the interface are only −0.2 • C, 0.05 • C and 0.025 • C for L × H = 125 × 10, 250 × 20 and 500 × 40, respectively, which are acceptable for engineering application.And the lattice resolution of 4 mm is selected in the following sections.of thermal diffusivity.A good consistence can also be seen between the numerical results and analytical ones, which indicates the validity of the proposed model in handling the heat conduction problem with phase change.The errors are defined as the numerical solutions of temperature minus the analytical ones.As shown in Figure 5b, the errors fluctuate near the solid-liquid interface, the finer grid resolution has the smaller error, and the farther the distance from the interface, the smaller the error.The maximum errors at the interface are only −0.2 °C, 0.05 °C and 0.025 °C for L × H = 125 × 10, 250 × 20 and 500 × 40, respectively, which are acceptable for engineering application.And the lattice resolution of 4 mm is selected in the following sections.

Results and Discussion
In practical engineering, the freezing pipes are usually arranged in a rectangular (or diamond) shape.In this paper, the four freezing pipes arranged in a square are selected as an example, which is shown in Figure 6.The development of frozen soil wall and temperature distribution are studied during the freezing process.The dimension of physical model is 4.0 m × 4.0 m.The spacing of freezing pipes is 1.2 m, and the outer diameter is 0.12 m.To ensure the mesh accuracy of the freezing pipes, the entire domain is divided into a lattice of 1000 × 1000 grid cells.The temperature of freezing pipes T is kept at −30 °C, the freezing temperature f T is 0 °C, and the initial temperature of soil T is

Results and Discussion
In practical engineering, the freezing pipes are usually arranged in a rectangular (or diamond) shape.In this paper, the four freezing pipes arranged in a square are selected as an example, which is shown in Figure 6.The development of frozen soil wall and temperature distribution are studied during the freezing process.The dimension of physical model is 4.0 m × 4.0 m.The spacing of freezing pipes is 1.2 m, and the outer diameter is 0.12 m.To ensure the mesh accuracy of the freezing pipes, the entire domain is divided into a lattice of 1000 × 1000 grid cells.The temperature of freezing pipes T 0 is kept at −30 • C, the freezing temperature T f is 0 • C, and the initial temperature of soil T i is 10 • C. The thermophysical parameters of soil are shown in Table 1.The freezing pipes are set as the constant temperature and the four side boundaries of the model are thermally insulated.The temporal evolutions of frozen zone are presented in Figure 7, and the time-history curves of temperature at the points A, B, C, D and O are shown in Figure 8.It can be seen that the temperature at the point A, which is closer to the freezing pipe, has larger thermal gradient, the soil freezes quickly, the latent heat has little influence on it, and the time-history curve is smooth.The temperature at the points C, D and O, which is farther from the freezing pipe, has a similar temperature evolution trend.The time-history curves show strong multistage, and for point C it can be divided into four stages as shown in Figure 8: (1) Cooling: the temperature drops rapidly at the early stage of artificial freezing, and reaches 0 °C in about 10 days.(2) Phase change: when the temperature drops to 0 °C, soil begins to freeze and releases the latent heat.The further away from the freezing pipe, the slower the energy transfers, and the longer the persistent time of phase change stage.(3) Partly frozen: the temperature descends faster in this stage, because there is larger temperature gradient and the thermal diffusivity of frozen soil is greater than that of unfrozen soil.(4) Completely frozen: the temperature evolution is mainly affected by the thermal diffusivity of frozen soil in this stage, the overall trend is relatively stable.For the point B, the distance from freezing pipe is moderate, the temperature is somewhere in between, and shows insignificant multistage.3) Partly frozen: the temperature descends faster in this stage, because there is larger temperature gradient and the thermal diffusivity of frozen soil is greater than that of unfrozen soil.(4) Completely frozen: the temperature evolution is mainly affected by the thermal diffusivity of frozen soil in this stage, the overall trend is relatively stable.For the point B, the distance from freezing pipe is moderate, the temperature is somewhere in between, and shows insignificant multistage.
and reaches 0 °C in about 10 days.(2) Phase change: when the temperature drops to 0 °C, soil begins to freeze and releases the latent heat.The further away from the freezing pipe, the slower the energy transfers, and the longer the persistent time of phase change stage.(3) Partly frozen: the temperature descends faster in this stage, because there is larger temperature gradient and the thermal diffusivity of frozen soil is greater than that of unfrozen soil.(4) Completely frozen: the temperature evolution is mainly affected by the thermal diffusivity of frozen soil in this stage, the overall trend is relatively stable.For the point B, the distance from freezing pipe is moderate, the temperature is somewhere in between, and shows insignificant multistage.Figure 9 shows the temperature distribution in the main section.Under the action of freezing pipes, the soil temperature decreases rapidly in 10 days, and there is funnel-shaped distribution around the freezing pipes.In about 20 days, the frozen soil wall overlaps in the main section, after that the temperature drops rapidly until soil completely frozen between double rows of freezing pipes in about 40 days, and then the temperature decrease rate slows down gradually.Figure 9 shows the temperature distribution in the main section.Under the action of freezing pipes, the soil temperature decreases rapidly in 10 days, and there is funnel-shaped distribution around the freezing pipes.In about 20 days, the frozen soil wall overlaps in the main section, after that the temperature drops rapidly until soil completely frozen between double rows of freezing pipes in about 40 days, and then the temperature decrease rate slows down gradually.

Cooling Phase change
Partly frozen Completely frozen Figure 9 shows the temperature distribution in the main section.Under the action of freezing pipes, the soil temperature decreases rapidly in 10 days, and there is funnel-shaped distribution around the freezing pipes.In about 20 days, the frozen soil wall overlaps in the main section, after that the temperature drops rapidly until soil completely frozen between double rows of freezing pipes in about 40 days, and then the temperature decrease rate slows down gradually.Figure 10 shows the temperature distribution in the intersection.The distance from the freezing pipes is relatively farther, so the temperature development in the intersection is obviously slower than that in the main section.In about 30 days, most of soil has been frozen in the intersection, then the temperature decreases rapidly, and the stable frozen soil wall forms in about 40 days.The time-history curves of freezing front between two freezing pipes are present in Figure 11, which show that the closer the spacing of freezing pipe is, the faster the freezing front develops, but the whole difference is not significant.Before the frozen soil wall overlapped, the spacing of freezing pipes has little effect on the evolution of freezing front.Figure 10 shows the temperature distribution in the intersection.The distance from the freezing pipes is relatively farther, so the temperature development in the intersection is obviously slower than that in the main section.In about 30 days, most of soil has been frozen in the intersection, then the temperature decreases rapidly, and the stable frozen soil wall forms in about 40 days.Figure 10 shows the temperature distribution in the intersection.The distance from the freezing pipes is relatively farther, so the temperature development in the intersection is obviously slower than that in the main section.In about 30 days, most of soil has been frozen in the intersection, then the temperature decreases rapidly, and the stable frozen soil wall forms in about 40 days.The time-history curves of freezing front between two freezing pipes are present in Figure 11, which show that the closer the spacing of freezing pipe is, the faster the freezing front develops, but the whole difference is not significant.Before the frozen soil wall overlapped, the spacing of freezing pipes has little effect on the evolution of freezing front.The time-history curves of freezing front between two freezing pipes are present in Figure 11, which show that the closer the spacing of freezing pipe is, the faster the freezing front develops, but the whole difference is not significant.Before the frozen soil wall overlapped, the spacing of freezing pipes has little effect on the evolution of freezing front.Figure 13 shows the relationship between the overlapping time and the spacing of freezing pipes.The overlapping time of frozen soil wall is greatly influenced by the spacing.The freezing times at the point C and point O increase with the increasing of the spacing, and there is a power function relationship between them.Therefore, at the design stage of artificial freezing engineering, the spacing of freezing pipes should be decided according to the overlapping time of artificial frozen soil wall.Figure 12 shows the thickness evolution of frozen soil wall at the point C along the x direction.The thickness develops faster at the early stage of frozen soil wall overlapped, and the developing speed gradually slows down as time goes on.In general, the thickness evolution has a similar tendency, but for the different freezing pipe spacing, there is a significant difference in overlapping time at the point O.When the spacing is 1.0 m, 1.2 m, 1.4 m, respectively, the required time for soil freezing is 8 day, 11 day and 16 day from point C to point O. Figure 13 shows the relationship between the overlapping time and the spacing of freezing pipes.The overlapping time of frozen soil wall is greatly influenced by the spacing.The freezing times at the point C and point O increase with the increasing of the spacing, and there is a power function relationship between them.Therefore, at the design stage of artificial freezing engineering, the spacing of freezing pipes should be decided according to the overlapping time of artificial frozen soil wall.Figure 13 shows the relationship between the overlapping time and the spacing of freezing pipes.The overlapping time of frozen soil wall is greatly influenced by the spacing.The freezing times at the point C and point O increase with the increasing of the spacing, and there is a power function relationship between them.Therefore, at the design stage of artificial freezing engineering, the spacing of freezing pipes should be decided according to the overlapping time of artificial frozen soil wall.

Conclusion
(1) Based on the enthalpy method, a lattice Boltzmann model is proposed to simulate the heat conduction problem with phase change.The model is applied to test the solid-liquid phase change of pure substance, and the results show that the evolution of both temperature distribution and solidliquid interface are in good agreement with the analytical solutions, and the mesh resolution has little effect on the numerical results.
(2) The temperature evolution of soil is associated with the distance from freezing pipe.When it is closer to the freezing pipe, the time-history curve of temperature is smoother, which is less affected by the latent heat.While it is farther, the time-history curve shows strong multistage, which can be divided into four stages: cooling, phase change, partly frozen, and completely frozen.
(3) Due to the effect of freezing pipes, the soil temperature in the main section decreases rapidly, there is funnel-shaped distribution around the freezing pipes, and the frozen soil wall is overlapped in about 20 days.The temperature development in the intersection is obviously slower compared with that in the main section.In about 30 days, most of soil has been frozen in the intersection, then the temperature decreases rapidly, and the stable frozen soil wall forms in about 40 days.
(4) The spacing of the freezing pipes has a significant influence on the overlapping time of artificial frozen soil wall, and there is a power function relationship between them.But it has little effect on the evolution of freezing front and the thickness of frozen soil wall.

Conclusion
(1) Based on the enthalpy method, a lattice Boltzmann model is proposed to simulate the heat conduction problem with phase change.The model is applied to test the solid-liquid phase change of pure substance, and the results show that the evolution of both temperature distribution and solid-liquid interface are in good agreement with the analytical solutions, and the mesh resolution has little effect on the numerical results.
(2) The temperature evolution of soil is associated with the distance from freezing pipe.When it is closer to the freezing pipe, the time-history curve of temperature is smoother, which is less affected by the latent heat.While it is farther, the time-history curve shows strong multistage, which can be divided into four stages: cooling, phase change, partly frozen, and completely frozen.
(3) Due to the effect of freezing pipes, the soil temperature in the main section decreases rapidly, there is funnel-shaped distribution around the freezing pipes, and the frozen soil wall is overlapped in about 20 days.The temperature development in the intersection is obviously slower compared with that in the main section.In about 30 days, most of soil has been frozen in the intersection, then the temperature decreases rapidly, and the stable frozen soil wall forms in about 40 days.
(4) The spacing of the freezing pipes has a significant influence on the overlapping time of artificial frozen soil wall, and there is a power function relationship between them.But it has little effect on the evolution of freezing front and the thickness of frozen soil wall.

16 Figure 2 .
Figure 2. Flowchart of program realization for soil freezing.

Figure 2 .
Figure 2. Flowchart of program realization for soil freezing.

Figure 2 .
Figure 2. Flowchart of program realization for soil freezing.

Figure 3 .
Figure 3. Diagram of heat conduction problem with phase change in semi-infinite space.To verify accuracy of proposed model, LBM is applied to simulate the heat conduction problem with solid-liquid phase change in semi-infinite space, as shown in Figure3.At the time 0 t = , the uniform initial temperature is i T , and the substance is in the liquid state.The cold source temperature 0 T is set at the position of 0 x= , and keeps constant at 0 t > .The analytical solution of the solid-liquid interface ( ) S t and temperature ( , ) T x t are as follows[26]

Figure 3 .
Figure 3. Diagram of heat conduction problem with phase change in semi-infinite space.

Figure 4 .
Figure 4. Comparisons between lattice Boltzmann method (LBM) and analytical solutions for the solid-liquid interface.

Figure 4 .Figure 5 .
Figure 4. Comparisons between lattice Boltzmann method (LBM) and analytical solutions for the solid-liquid interface.Energies 2018, 11, x FOR PEER REVIEW 9 of 16

Figure 6 .
Figure 6.Schematic diagram of the calculation model.The temporal evolutions of frozen zone are presented in Figure7, and the time-history curves of temperature at the points A, B, C, D and O are shown in Figure8.It can be seen that the temperature at the point A, which is closer to the freezing pipe, has larger thermal gradient, the soil freezes quickly, the latent heat has little influence on it, and the time-history curve is smooth.The temperature at the points C, D and O, which is farther from the freezing pipe, has a similar temperature evolution trend.The time-history curves show strong multistage, and for point C it can be divided into four stages as shown in Figure8:(1) Cooling: the temperature drops rapidly at the early stage of artificial freezing, and reaches 0 • C in about 10 days.(2) Phase change: when the temperature drops to 0 • C, soil begins to freeze and releases the latent heat.The further away from the freezing pipe, the slower the energy transfers, and the longer the persistent time of phase change stage.(3) Partly frozen: the temperature descends faster in this stage, because there is larger temperature gradient and the thermal diffusivity of frozen soil is greater than that of unfrozen soil.(4) Completely frozen: the temperature evolution is mainly affected by the thermal diffusivity of frozen soil in this stage, the overall trend is relatively stable.For the point B, the distance from freezing pipe is moderate, the temperature is somewhere in between, and shows insignificant multistage.

Figure 8 .
Figure 8.The time-history curves of temperature at the points A, B, C, D and O.

Figure 8 .
Figure 8.The time-history curves of temperature at the points A, B, C, D and O.

Figure 8 .
Figure 8.The time-history curves of temperature at the points A, B, C, D and O.

Figure 9 .
Figure 9. Temperature distribution in the main section.

Figure 10 .
Figure 10.Temperature distribution in the intersection.

Figure 9 .
Figure 9. Temperature distribution in the main section.

Figure 10 .
Figure 10.Temperature distribution in the intersection.

Figure 10 .
Figure 10.Temperature distribution in the intersection.

Figure 11 .
Figure 11.The time-history curve of freezing front.

Figure 12
Figure12shows the thickness evolution of frozen soil wall at the point C along the x direction.The thickness develops faster at the early stage of frozen soil wall overlapped, and the developing speed gradually slows down as time goes on.In general, the thickness evolution has a similar tendency, but for the different freezing pipe spacing, there is a significant difference in overlapping time at the point O.When the spacing is 1.0 m, 1.2 m, 1.4 m, respectively, the required time for soil freezing is 8 day, 11 day and 16 day from point C to point O.

Figure 12 .
Figure 12.The thickness evolution of frozen soil wall.

Figure 11 .
Figure 11.The time-history curve of freezing front.

Figure 12 16 Figure 11 .
Figure12shows the thickness evolution of frozen soil wall at the point C along the x direction.The thickness develops faster at the early stage of frozen soil wall overlapped, and the developing speed gradually slows down as time goes on.In general, the thickness evolution has a similar tendency, but for the different freezing pipe spacing, there is a significant difference in overlapping time at the point O.When the spacing is 1.0 m, 1.2 m, 1.4 m, respectively, the required time for soil freezing is 8 day, 11 day and 16 day from point C to point O.

Figure 12 .
Figure 12.The thickness evolution of frozen soil wall.

Figure 12 .
Figure 12.The thickness evolution of frozen soil wall.

Figure 13 .
Figure 13.Relationship between the overlapping time and the spacing of freezing pipes.

Figure 13 .
Figure 13.Relationship between the overlapping time and the spacing of freezing pipes.

Table 1 .
Comparisons of thermophysical parameters between physical and lattice units.