In order to numerically solve the heat conduction equation while considering the phase change, the following basic assumptions were made: 1) The soil was considered a homogeneous isotropic continuous medium, which was divided into solid and liquid phase zones according to the physical state of the internal moisture. The thermophysical parameters of each phase zone were fixed at constant values and were made independent of the temperature; 2) The effects of moisture migration, soil frost heaving, and stress field change were ignored. Whereas the heat conduction effect during the freezing process was only considered; 3) The heat conduction and the movement of the phase change interface during the soil’s freezing process are along a single direction; 4) The effect of the solid-liquid phase change interface’s thickness in the soil body is not considered; 5) The soil’s freezing temperature is constant. Besides, the freezing process was completed instantaneously, ignoring its physical evolution process.

The mathematical equation for the temperature field evolution can be simplified according to the temperature-induced change in the physical state of water in the soil, which is divided into the solid-phase zone and liquid-phase zones (Ozisik, 1993) as follows:

Solid-phase zone: ∂ T s x , t ∂ t = α s ∂ 2 T s x , t ∂ x 2   0 < x < s t (1)

Liquid-phase zone: ∂ T l x , t ∂ t = α l ∂ 2 T l x , t ∂ x 2   s t < x < L (2) where T s and T l are the soil’s temperatures in the solid-phase zone and liquid-phase zone, respectively; α s and α l are the thermal diffusion coefficients of the soil in the solid-phase zone and liquid-phase zone, respectively, in which α = k / ρ C p , where ρ is the density, C p is the specific heat capacity, and k is the thermal conductivity that can be obtained from the volume fraction of water (or ice), soil particles, and air in the soil and their respective thermophysical parameters (Radoslwa et al., 2006); t is the freezing time.

In order to study the movement process for the phase change interface during soil freezing, the enthalpy parameter was introduced into the heat conduction equation according to the enthalpy model proposed by Shamsundar and Sparrow (Shamsundar and Sparrow, 1975) as follows: ρ ∂ H ∂ t = k ∂ 2 T ∂ x 2 (3) where H is the enthalpy; C p is the soil’s specific heat capacity, which is assumed to not change with temperature; H can be expressed as: H = C p T + φ L a (4) where φ is the liquid-phase fraction, in which φ = 1 in the liquid phase and φ = 0 in the solid phase.

Substituting Eq. 4 into Eq. 5 yields the following: ∂ T ∂ t = α ∂ 2 T ∂ x 2 − L a C p ∂ φ ∂ t (5)

Comparing Eq. 5 to the heat conduction in Eqs. 1, 2, it can be seen that Eq. 5 adds a source term related to the latent heat of the phase change L a in the heat conduction equation.

When the soil’s freezing temperature is a constant value T f , the relationship between the liquid-phase fraction and the enthalpy can be defined as follows: φ = 0   H < C p T f H − C p T f L a   C p T f ≤ H ≤ C p T f + L a 1   H > C p T f + L a (6)

The temperature of the soil’s cold source can be defined as follows: T 0 , t = T 0 (7)

The soil body has a constant initial temperature that is defined as follows: T x , 0 = T i (8)

At the solid-liquid phase change interface y = s t : T s x , t = T l x , t = T f (9) k s ∂ T s x , t ∂ x − k l ∂ T l x , t ∂ x = ρ s L a d s t d t (10)

In order to simplify the calculation, the soil’s heat conduction equation and the corresponding boundary conditions are normalized as follows: X = x L S = s L F 0 = α s t L 2 θ = T − T 0 T i − T 0 (11) where X is the dimensionless length; S is the dimensionless solid-liquid phase change position; F 0 is the dimensionless time; θ is the dimensionless temperature.

Substituting Eq. 11 into Eqs. 1, 2 yields the following:

Solid-phase zone: ∂ θ s X , F 0 ∂ F 0 = ∂ 2 θ s X , F 0 ∂ X 2   0 < X < S F 0 (12)

Liquid-phase zone: ∂ θ l X , F 0 ∂ F 0 = α s α l ∂ 2 θ l X , F 0 ∂ X 2   S F 0 < X < 1 (13)

The corresponding boundary condition transformations are cold source temperature θ s 0 , F 0 = 0 and initial temperature θ l X , 0 = 1 .

The following equations are valid at the solid-liquid phase change interface X = S F 0 : θ s X , F 0 = θ l X , F 0 = T f − T i T 0 − T i (14) ∂ θ s X , F 0 ∂ X − k l k s ∂ θ l X , F 0 ∂ X = T f − T i T 0 − T i S t e d S F 0 d F 0 (15) where S t e is the Stephen number, a dimensionless quantity that is related to the phase change process, S t e = C p T f − T 0 L a .

In this study, the temperature in the macroscopic heat conduction equation is taken as a scalar, and the D2Q4 model proposed by Qian et al. (Qian et al., 1992) is used for the soil’s temperature field evolution process during freezing (Figure 2). The time evolution of the particle temperature distribution function g i r , t can be expressed by the discrete lattice Boltzmann equation in the form of BGK (Bhatnagar, Gross, and Krook) as: g i r + e i δ t , t + δ t − g i r , t = − g i r , t − g i e q r , t τ + S r i δ t   i = 0 , 1 , 2 , 3 (16) where g i r , t is the particle temperature distribution function along direction i at lattice point r and moment t ; e i is the discrete velocity, which consists of a set of velocity vectors in four directions defined as follows: e = c 1,0 , − 1,0 , 0,1 , 0 , − 1 (17) where c is the lattice velocity, c = δ x δ t ; δ x and δ t are the discrete lattice step and time step, respectively; τ is the dimensionless relaxation time, in which a stable solution is usually obtained for τ = 1.0 (Sukop and Thorne, 2006); g i e q r , t is the equilibrium state distribution function. Since the convection effect was not considered in the soil freezing process, the D2Q4 model g i e q r , t can be expressed by the following equation: g i e q r , t = ω i T r , t (18) where ω i is the weight factor, ω 0 = ω 1 = ω 2 = ω 3 = 1 4 .

The discrete source term S r i in Eq. 16 is defined as: S r i = ω i S r (19) where S r can be obtained from Eq. 5 as: S r = L a C p ∂ φ ∂ t (20)

Substituting Eqs. 18, 19, and 20 into Eq. 16, the final discrete lattice Boltzmann equation for the freezing process of the soil is obtained as follows: g i r + e i δ t , t + δ t − g i r , t = − g i r , t − g i e q r , t τ + ω i L a C p φ t + δ t − φ t   i = 0 , 1 , 2 , 3 (21)

In this study, Eq. 21 was implemented in two steps for the convenience of the programming calculations. These steps are as follows:

(1) Collision

Based on the Chapman-Enskog multiscale expansion method, the lattice Boltzmann Eq. 21 can be reduced to Eq. 5 that contains the phase change by considering the relationship between the macroscopic temperature T and the particle temperature distribution function (i.e., T = ∑ i = 0 3 g i r , t ). The expression for the thermal diffusion coefficient α was obtained (Mohamad, 2011) as follows: α = c s 2 τ − 1 2 δ t (24) where c s is the lattice speed of sound. For the D2Q4 model, c s 2 = 1 2 c 2 .

Since the soil’s thermal diffusion coefficients in the liquid and solid phase zones during freezing soil differ, the relaxation times also vary. In this paper, the variation of the variable relaxation time with the liquid-phase fraction was considered and set according to the thermal diffusion coefficients of different phases as follows: τ = φ α l + 1 − φ α s δ t c s 2 + 1 2 (25)

The soil’s thermal diffusion coefficient is found according to the percentage of the liquid phase fraction among different phases. Consequently, a unified equation can solve the thermal diffusion coefficient of different phases according to the changes in liquid phase fraction. During the numerical calculation, the relaxation time is obtained from Eq. 25, which is then brought into Eq. 16 for the evolution of the temperature field.

Accurate simulation of the boundary conditions is an important part of the numerical calculation that significantly influences accuracy, efficiency, and stability. The following determines the corresponding particle temperature distribution functions according to the macroscopic boundary conditions:

(1) Initial conditions

Assuming that the temperature of the soil body at the initial moment is constant T i , and the particle temperature distribution function at each grid point is in equilibrium, then: g i x , 0 = ω i T i   i = 0,1,2,3 (26)

(2) Boundary conditions

The non-equilibrium state extrapolation format proposed by Guo et al. (Guo et al., 2002) was used to handle the seepage field boundary of the soil. For the temperature field boundary, there is no essential difference between the two in terms of application (Yong et al., 2009; Gao et al., 2014). This method produces the distribution function with overall second-order accuracy. The basic idea is that the particle temperature distribution function g i B , t on the boundary lattice point B can be decomposed into two parts: the equilibrium state g i e q B , t and the non-equilibrium state g i n e q B , t , that is: g i B , t = g i e q B , t + g i n e q B , t (27)

For the equilibrium part g i e q B , t , the cold source constant temperature boundary (temperature is always T 0 ) can be obtained from Eq. 18. The adiabatic boundary can be taken as g i e q O , t with its neighbouring grid point O . The solution of the non-equilibrium part g i n e q B , t is complicated. However, the calculation can be simplified by approximating g i n e q O , t of the adjacent lattice points in the temperature field as follows: g i B , t = g i e q B , t + g i O , t − g i e q O , t (28)

The program’s parameters are typically taken in lattice units during the numerical calculation of the lattice Boltzmann method. In contrast, physical problems are usually taken in physical units. Hence, a conversion relationship between physical and lattice units must be established. For this purpose, all physical and lattice units involved in the calculation must be dimensionless separately to ensure that the heat conduction criterion before and after dimensionless processing remains the same. Accordingly, a dimensionless parameter was used as a connection bridge to realize the unit conversion between the physical and lattice units.

Based on the dimensionless time, the relationship between the physical time t and the lattice time step can be established as follows: F 0 = α p t p L p 2 = α L t L L L 2 (29) where α p and α L are the thermal diffusion coefficients of the physical and lattice units, respectively; t p and t L are the times of the physical and lattice units, respectively ( t L = N δ t , where N is the number of calculation time steps); L P and L L are the reference lengths of physical and lattice units, respectively.

If L P is the length of the whole calculation domain, then L L = n δ x , where n is the number of lattices in the calculation domain.

The correspondence between the physical unit phase change parameters (latent heat of phase change L α , specific heat capacity C p and the lattice unit can be established by combining the dimensionless temperature based on the dimensionless S t e number as follows: S t e = C p p T f p − T 0 p L a p = C p L T f L − T 0 L L a L (30) where C p p and C p L are the specific heat capacities of the physical and lattice units, respectively; L a p and L a L are the latent heats of phase change in the physical and lattice units, respectively; T f p and T f L are the phase change temperatures of the physical and lattice units, respectively; T 0 p and T 0 L are the cold source temperatures of the physical and lattice units, respectively.

In this study, a lattice Boltzmann model based on the enthalpy method is established, and the corresponding calculation procedure is prepared by considering the effects of heat conduction, latent heat release from solid-liquid phase change, and phase change interface movement in the soil’s freezing process. The adopted flowchart is shown in Figure 3.

Under a unilateral cold source, the temperature of the soil body decreases rapidly, leading to the freezing front’s formation by condensing internal water into ice. Figure 6 shows the evolution of the freezing front in the soil body with time. The figure shows that at the early stage of freezing, the freezing front in the soil body moves faster, and the evolution rate slows down with time. Besides, the frozen front moved 0.37, 0.52, 0.74, 0.89, and 1.02 m at 5, 10, 20, 30, and 40 days, respectively.

Figure 7 shows the development of the soil’s temperature field under the unilateral cold source. The soil temperature decreases rapidly at the early stage of freezing because the cold source temperature is − 30 ℃ , much lower than the initial temperature of 10 ℃ . Thereafter, over time, the temperature field decreases because the temperature difference in the soil gradually becomes smaller. Considering the freezing front as the boundary, the two sides present a significant difference in the trend of the temperature field due to the influence of heat release from the solid-liquid phase change and different thermal diffusion coefficients. The soil temperature in the solid-phase zone develops faster and is approximately linear; in the liquid-phase zone, the soil develops slowly and gradually converges to the initial temperature of the soil. According to the code (DG/T J08-902, 2016), when the soil excavation depth is 12–30 m, the design reference value of the average temperature in the freezing wall must be taken as − 8 ∼ − 10 ℃ . Therefore, when − 10 ℃ is taken as the limit for freezing wall temperature, the thickness of the freezing wall developed at 5, 10, 20, 30, and 40 days under the unilateral cold source is 0.24, 0.33, 0.47, 0.57, and 0.66 m, respectively.

In artificial freezing works, multiple tubes are usually used to form the soil’s freezing curtain. Therefore, the interaction between two or more cold sources must be studied. For this purpose, the temperature of the cold source was applied on both sides of the soil within 2.0 m width to study the variation of the temperature field in the soil. Figure 8 shows the trend in the variation of the liquid-phase fraction in the soil for 5, 10, 20, and 30 days. It can be seen that the position of the freezing front moves 0.37, 0.54, 0.78, and 0.98 m under the action of the bilateral cold sources for the time periods of 5, 10, 20, and 30 days. Compared to the unilateral cold source, the freezing front’s movement has moved faster with time.

Figure 9 shows that due to the bilateral cold sources, the temperature in the soil gradually decreases from both sides to the middle. The overall temperature in the soil drops below 0 ℃ in about 30 days, and all the water condenses into ice. The temperature in the soil decreases rapidly after that, and at 35 days and 45 days, the overall temperature in the soil drops below − 13.6 ℃ and − 26.4 ℃ , which are less than the average temperature of the frozen wall design’s reference value. Hence, the excavation work of underground works can be carried out.