- Email: [email protected]
A lattice Boltzmann model for heat and mass transfer phenomena with phase transformations in unsaturated soil during freezing process
International Journal of Heat and Mass Transfer 94 (2016) 29–38
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
A lattice Boltzmann model for heat and mass transfer phenomena with phase transformations in unsaturated soil during freezing process Wenyu Song, Yaning Zhang, Bingxi Li ⇑, Xinmeng Fan School of Energy Science and Engineering, Harbin Institute of Technology, 92 West Dazhi Street, Nan Gang District, Harbin 150001, China
a r t i c l e
i n f o
Article history: Received 18 May 2015 Received in revised form 29 October 2015 Accepted 3 November 2015
Keywords: Lattice Boltzmann method Frozen soil Stochastic structure Water content distribution
a b s t r a c t In this paper, a mesoscopic numerical model for simulating the heat and moisture transport phenomena in frozen soil during freezing process is presented. The model includes a stochastic generation-growth method for reconstructing the soil structure with given macroscopic geometric parameters, a lattice Boltzmann method (LBM) model for solving multiphase fluid flow and a LBM model for solving heat conduction with phase change in porous media. By using the present model, freezing process in sandy loam soil is demonstrated. Two types of boundary condition, prescribed temperature and adiabatic boundary conditions are employed in the given cases. All the numerical results show good agreement with experimental results. Considering the mesoscopic nature of LBM, the proposed model is a potential alternative to traditional continuum models. Ó 2015 Elsevier Ltd. All rights reserved.
1. Introduction Frozen ground region is generally distributed on earth. Permafrost region covers approximately 23% of the global land surface; seasonally frozen ground region can be found all over the regions above 24° latitude [1]. The freezing–thawing process of frozen soil has great effects on regional climate and hydrology [2,3]. Moreover, the abundant greenhouse gases in frozen soil are also important to global warming [4,5]. However, the coexistence of heat and mass transfer, phase change and multiphase flow during the freezing–thawing process makes it considerably difficult to give a simple and complete description of water content and temperature distribution in frozen soil. Therefore, a favorable heat and moisture transfer model shows its great importance to frozen soil research. The coupled heat and mass transport process was introduced into frozen soil model to solve the heat and moisture balance problems in the early 1970s [6–8]. In the earlier studies, certain processes, for instance, vapor flux and its corresponding phase change were omitted for simplicity [6–9]. During the past decades, more evidence revealed the complexity of heat and mass transport in freezing–thawing process of frozen soil [10–12]. Thus far, various numerical models with different levels of complexity have been developed. In term of governing equations, the most complex as well as the most comprehensive model includes ten variables determined by four prognostic equations and five diagnostic ⇑ Corresponding author. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2015.11.008 0017-9310/Ó 2015 Elsevier Ltd. All rights reserved.
equations with one assumption [11–13]; in the relatively simple model which is the most commonly used, only two prognostic equations are used for solving temperature and liquid water content. In term of physical processes, some models take into account the vapor movement and its corresponding phase change in the water and heat balances [10,14] while others neglect those processes [15–17]. Furthermore, even with the most complex model, the estimation of ice content and position of the ice lens remains a challenging task. One reason is the highly nonlinear relationship between temperature, ice content and liquid water content [13]. On the other hand, fundamentally, soil is granular porous material. Considerable research shows that the heat and mass transport in porous material is strongly affected by its pore structure [18–20]. The underlying problems raise the need for alternative approaches to understand the connection between the macroscopic phenomena and the underlying microscopic dynamic at a more fundamental level [21]. The conventional numerical methods such as finite difference method (FDM), finite element method (FEM) and finite volume method (FVM) are based on the discretization of macroscopic continuum equations. This scheme makes conventional numerical methods face great challenges while solving flow with complex interface or flow with complex boundary. Since the late 1980s, a new mesoscopic method base on the discrete kinetic theory, named the lattice Boltzmann method (LBM), has become a promising tool for the investigations of a wide range of complex flows, including multiphase flows, porous flows, thermal flows and reactive transport [22]. The fundamental idea of LBM is quite
30
W. Song et al. / International Journal of Heat and Mass Transfer 94 (2016) 29–38
Nomenclature c cs C e f f1 F g g g0 G H L s t T u u0 DV x, y x y0
a v
Dt Dx e e0
m q
lattice speed speed of sound specific heat discrete velocity density distribution function of fluid lattice liquid fraction force density distribution function of thermal lattice gravitational acceleration sum of unknown distribution functions adjustment coefficient enthalpy latent heat tag variable time temperature velocity common velocity of the fluid mixture difference of volume axial coordinate lattice site axial coordinate of boundary thermal diffusivity mass fraction time step lattice spacing porosity given porosity kinematic viscosity density
straightforward: the macroscopic dynamics of fluid are the result of the collective behavior of many microscopic particles in the system and the macroscopic dynamics are not sensitive to the underlying details in microscopic physics [23]. Thus, LBM can be considered as a bridge sitting in the intermediate area between microscopic molecular dynamics and macroscopic fluid dynamics. The major advantage of LBM is that it behaves like a solver for the conservation equations in the bulk flow, meanwhile the kinetic equation offers many advantages of molecular dynamics, including clear physical pictures, easy implementation of complex boundary conditions and fully parallel algorithms [22]. In this work, an advanced mesoscopic numerical model based on LBM for simulating the heat and mass transfer phenomena with phase transformation in frozen soil during freezing process is presented. The fundamental idea of the proposed model is analogous to LBM: the complex phenomena in frozen soil during freezing process are caused by the interactions between several common physical phenomena and the stochastic structure of soil. The model is validated by comparing the numerical results with existed experimental results. This model should not only provide an advanced tool for the frozen soil research but also for the investigation of the similar process in generalized porous media as well. 2. Numerical models 2.1. Reconstruction of soil structure In nature, soils are heterogeneous assemblies of materials, forming porous media [24]. Various studies have shown the importance of the stochastic approaches for flow, heat transfer and solute transport in porous media [24–26]. In this paper, a
q0 s
ud ug
w w x
given density relaxation time core distribution probability growth probability effect mass water potential weight factor
Superscripts ads adhesive force b body force coh cohesive force eq equilibrium state k number of iteration tmp temperature Subscripts A air pressure i index of speed direction of lattice l liquid state P pressure s solid state S osmotic Tot total tmp temperature Z position a thermal conductivity r first phase r second phase m kinematic viscosity
stochastic generation-growth method is adopted for soil structure reconstruction based on the given porosity. The reconstruction process is conducted as follows: (i) Randomly locate the growing cores based on a given core distribution probability ud. The value of ud represents the quantity of particles in soil and further determines the size of particles at a given porosity e0. (ii) Define the growing probability of growing cores in each direction ug,i. For a D2Q9 lattice Boltzmann model, there are 8 directions on the growing cores according to the velocity directions of a D2Q9 lattice. (iii) Grow particles from the cores along the given directions; the newly grown particles are used as growing cores in the next cycle; calculate porosity e, stop growing once e meets the given porosity e0, otherwise go back to step (i). (iv) For some critical circumstances, e.g. for the permeable structure with low porosity, an extra step is executed to identify connected regions. If the model meets the requirement, stop growing, otherwise go back to step (i). Thus, for an area of a given size, the quantity of particles is determined by ud. The shape of a particle is defined by ug, isotropic distribution gives the particles round-like shape; anisotropic distribution gives the particles square-like shape, triangular-like shape or fiber-like shape as illustrated in Fig. 2. Therefore, by using a combination of different ug, a soil model with more realistic structure can be generated. Substantially, for the D2Q9 model, the generated soil structure is a D(4, 82) bond percolation network. The first row of Fig. 1 illustrates the soil structure generated in 300 300 grids while the second row shows the connectivity of
31
W. Song et al. / International Journal of Heat and Mass Transfer 94 (2016) 29–38
(a)
(b)
(c)
Fig. 1. Soil structure generated by stochastic generation-growth method and connectivity of pores. (a) e = 0.55, ud = 0.5; (b) e = 0.40, ud = 0.3; (c) e = 0.40, ud = 0.35.
Fig. 2. Structures generated by anisotropic growing probability ug.
pores with an assumption of fluid inlet from the lower boundary. In the first row, black represents solid elements and white represents fluid elements; in the second row, black represents impermeable region and white represents permeable region. It indicates that porosity e, core distribution probability ud and growing probability ug all have strong influences on fluid flow in porous media. For porous media with a low porosity, ud plays a more important role on permeability.
2.2. The Shan and Chen pseudopotential LBM model The pseudopotential model was proposed by Shan and Chen in 1993 [27–30]. Due to its simplicity and versatility, it is the most widely used multiphase LBM model [21] to the best knowledge
of the author. We adopt a single relaxation time collision operator (BGK) LB equation, which can be expressed as follows:
f r;i ðx; t þ DtÞ ¼ f r;i ðx; tÞ
1 eq f r;i ðx; tÞ f r;i ðx; tÞ
sr;m
f r;i ðx þ ei Dt; t þ DtÞ ¼ f r;i ðx; t þ DtÞ
ð1Þ ð2Þ
where f r;i ðx; tÞ is the density distribution function of the ith direction at the lattice site x and time t (see Fig. 3); feq represents the equilibrium distribution function; c = Dx/Dt is the lattice speed, in which Dx = Dt = 1 represents the lattice spacing and time step respectively; sm is the dimensionless relaxation time. In the computer algorithm, the collision step is represented by Eq. (1), which means the system reaches the local Maxwellian equilibrium on a time scale s. Eq. (2) represents the streaming step in which
32
W. Song et al. / International Journal of Heat and Mass Transfer 94 (2016) 29–38
fluid particles are streamed from one lattice to a neighboring lattice according to the discrete velocity ei in the lattice. For a D2Q9 lattice,
8 ð0;0Þ i¼0 > > > < i ¼ 1;2;3;4 ei ¼ ðcos ½ða 1Þp=2;sin ½ða 1Þp=2Þ > > > : pffiffiffi 2ðcos ½ða 5Þp=2þ p=4;sin ½ða 5Þp=2þ p=4Þ i ¼ 5;6;7;8
ð3Þ The equilibrium distribution function feq is given by
3 9 3 eq 2 2 f r;i ¼ xi qr 1 þ 2 ðei ueq Þ þ 4 ðei ueq Þ 2 ðueq Þ c 2c 2c
ð4Þ
mr ¼
sr;m 0:5ÞDt
ð5Þ
ð6Þ
X f r;i
ð7Þ
i
qr ur ¼
X f r;i ei
ð8Þ
i
the mass fraction v of phase r is calculated by
vr ¼
qr
In order to introduce the force term (e.g. fluid–fluid interaction force, fluid–solid interaction force, gravitational force and other bulk forces), a velocity shift scheme [27] is adopted. In this scheme, the parameter ueq r in Eq. (4) is determined by the relation
ueq ¼ u0 þ
sr;m Fr qr
ð10Þ
where u0 is the common velocity of the fluid mixture, which is given by
P u0 ¼
qr ur r s P qrr;m r sr;m
ð11Þ
In Eq. (10), Fr is the sum of all interactive forces between fluid ads b tmp coh particles. In this model, Fr ¼ Fcoh r þ Fr þ Fr þ Fr , Fr represents the cohesive force between particles from the same phase; Fads r stands for the adhesive force between fluid particles and solid boundary; Fbr is the sum of bulk forces such as gravitational force and buoyancy force; Ftmp r represents the bulk force caused by temperature potential change. In the current model, vapor in soil is neglected. Therefore the fluid system can be treated as a binary mixture of water and air. Only repulsion between different components is taken into account, the interactive force applied on the particles of one component is carried out by coh Fcoh r ðxÞ ¼ wr ðx; tÞG
X
xa wr ðx þ ei Þei
where si works as a tag, the value of si is
si ðx þ ei Þ ¼
0
node x þ ei is fluid lattice node
ð14Þ
1 node x þ ei is solid boundary
Similar to Eq. (12), Gads r adjusts the interaction force between the particles of component r and solid boundary. For wetting fluid, ads Gads r is chosen positive while for non-wetting fluid, Gr is chosen negative. The effect mass is
q w ¼ q0 1 exp
ð15Þ
Fr ¼ qr g b
ð16Þ
In frozen soil, the moisture moves along the temperature gradient from high temperature to low temperature during both freezing and thawing processes. Therefore, an accurate numerical model to simulate the temperature distribution in frozen soil is of great importance. In LB model, the heat conduction process is performed by applying another layer of thermal lattice. Similar to Section 2.2, a D2Q9 lattice is adopted; the LB equation for thermal diffusion can be expressed as follows:
g i ðx þ ei Dt; t þ DtÞ ¼ g i ðx; tÞ ð9Þ
qr þ qr
ð13Þ
2.3. LBM model for heat conduction problem with phase change
in which cs represents the lattice sound speed. The macroscopic density and momentum at each lattice node can be calculated as
qr ¼
xi si ðx þ ei Þei
i
Then, the gravitational force of component r is given by
The viscosity in lattice unit is given by
c2s ð
X
q0
where in Eq. (4), the weight factor xi is given by
8 i¼0 > < 4=9 xi ¼ 1=9 i ¼ 1; 2; 3; 4 > : 1=36 i ¼ 5; 6; 7; 8
ads Fads r ðxÞ ¼ wr ðx; tÞG
ð12Þ
i
represents the other component of the fluid mixture; Gcoh where r r is a adjuster to tune the surface tension between two components. The cohesive force between fluid particles and solid boundary can be expressed analogously by
1
sa
g i ðx; tÞ g eq i ðx; tÞ
ð17Þ
The discrete velocity ei in the lattice is expressed in Eq. (3). The equilibrium distribution function g eq i is given by
g eq i ¼ xi T
ð18Þ
The weight factor xi is given in Eq. (5). In analogy with the kinematic viscosity m, the thermal diffusivity a is calculated as
a ¼ c2s ðsa 0:5ÞDt
ð19Þ
The macroscopic temperature can be calculated as the 0th order momentum of gi as
X gi ¼ T
ð20Þ
i
In order to solve the heat conduction problem with phase change, Jiaung [31] proposed a LB model with a source term in Eq. (17), which can be written as
g i ðx þ ei Dt; t þ DtÞ ¼ g i ðx; tÞ þ xi
1
sa
g i ðx; tÞ g eq i ðx; tÞ
L
f ðt þ DtÞ f l;i ðtÞ C l;i
ð21Þ
In order to solve the source term, a sub-loop is introduced into every time step, which is executed as follows: (i) The value of the particle distribution function at the new time level t + Dt in the kth iteration of sub-loop is calculated according to Eq. (21) as
g ki ðx þ ei Dt; t þ DtÞ ¼ g i ðx; tÞ þ xi
1
sa
g i ðx; tÞ g eq i ðx; tÞ
L k1 f l;i ðt þ DtÞ f l;i ðtÞ C
ð22Þ
where the superscripts k and k1 represent the current and previous iterations, respectively.
W. Song et al. / International Journal of Heat and Mass Transfer 94 (2016) 29–38
2
6
33
respectively. For low salinity soil under standard atmospheric pressure, WS and WA can be neglected, hence
5
WTot ¼ WP þ WZ þ Wtmp
ð28Þ
ads In the LBM model, Wp is carried out by Fcoh r and Fr (caused by
0
3
capillary action); WZ is represented by Fbr and Wtmp is represented by Ftmp r . When the soil temperature is above freezing temperature, the gradient of water potential caused by temperature change can be derived from Clapeyron’s equation, which can be written as
1
rW ¼
7
4
(ii) The macroscopic temperature in the kth iteration is calculated by Eq. (20). (iii) The total enthalpies in the kth iteration is calculated as k1
ð23Þ
L
(iv) The liquid fraction in the current iteration is calculated as k
fl ¼
8 > < k0
H Hs > Hl Hs
:
1
H < Hs Hs 6 H 6 Hl
ð24Þ
H > Hl
where Hs and Hl represent the enthalpies at the freezing and melting temperatures, respectively. (v) Repeat steps (i)–(iv) until convergence is reached, the convergence criterion is expressed as
! f k f k1 T k T k1 8 min k1 ; < 10 f Tk
ð25Þ
rW ¼
ð26Þ
in which the dimensionless scale T* = (T T0)/(T1 T0); t* = at/l2; l is the natural length of the system; T0 and T1 are the initial temperature and the temperature applied on the boundary, respectively; St = C(T0 T1)/L is the Stefan number. Thus the total time steps needed for the simulation can be estimated by using Eq. (26).
Ls rT TjDV s j
ð30Þ
where Ls and DVs represent the latent heat of fusion and specific volume change during freezing. Thus the movement of water in soil can be quantified by rW. It is noted that In Eq. (30), Ls and DVs are known as the property of water, T = 273 K is the freezing point of water in Kelvin. Therefore, Eq. (30) becomes
rW ¼ 13:48 106 rT
ð31Þ
The temperature water potential has the same value with the suction pressure caused by freezing process can be written as
Pr ¼
Ls q u2 ¼ r r TjDV s j 2
ð32Þ
in lattice unit the suction pressure is also carried out by
P¼
In addition, when the liquid fraction fl in a thermal lattice reaches 0, the corresponding fluid lattice is converted to ‘‘Bounce back” lattice, which ensures the conservation of ice content. It is noteworthy that the number of thermal diffusivity in real physical units is so small that if we introduce such small thermal diffusivity the code will blow up. However, by using the total time of the evolution tmax as the timescale, we can get a nondimensional version of the differential equation of Eq. (21) which proposed by Huber et al. [32] is given by
@T at max 2 1 @f l ¼ 2 r T St @t @t l
ð29Þ
where Ll and DVl represent the latent heat of condensation and specific volume change during vapor condensation, respectively. Analogously, when the soil temperature reaches the freezing temperature, the gradient of water potential caused by temperature can be carried out by
8
Fig. 3. The layout of velocity directions of a D2Q9 lattice.
Hk ¼ CT k þ f l
Ll rT T DV l
qu2
ð33Þ
2
hence, the suction pressure in lattice unit becomes
P¼
qu2 P qr u2r r
ð34Þ
we have known that
ur csr ¼ u cs
ð35Þ
assuming that in soil, the real speed of sound csr = 150 m/s, the real density of water qr = 1000 kg/m3, it can be calculated that P = 0.199. However, for the real soil, Eq. (31) needs to be corrected. The correction factor depends on the type of soil [33]. Here we propose a relation to implement rW in lattices tmp F tmp r ¼ G rWDx tmp
where G
ð36Þ
is determined by the type of soil.
3. Verification of LBM model, results and discussion 2.4. The implementation of soil water potential in LBM model 3.1. Introduction of laboratory experiment In soil, water moves from one position to another due to gravity, pressure, capillary action, osmosis and so on. The tendency of water movement is quantified by water potential W. The total water potential WTot is given by
WTot ¼ WP þ WZ þ WS þ WA þ Wtmp
ð27Þ
where WP, WZ, WS, WA and Wtmp stand for pressure, gravitational, solute (osmotic), air pressure and temperature potentials,
Mizoguchi performed a classic laboratory experiment in 1990 [34], the result has been adopted in various studies [35,36]. In his research, Kanagawa sandy loam soil was packed in four identical cylinders. The dimensions of all four cylinders were identical, which was 20 cm long with an internal diameter of 8 cm. The specimens were prepared with the same initial conditions, which involved a uniform temperature of 6.7 °C and close volumetric
34
W. Song et al. / International Journal of Heat and Mass Transfer 94 (2016) 29–38
2
6
A
3
5
7
4
8
6
2
5
y
B
1
0
0
D
x Fig. 4. Boundary conditions in Mizoguchi’s experiment and in numerical simulation.
3
1
7
C
8
4
Fig. 5. Illustration of unknown distribution functions at boundaries AB and DC.
water content of 0.33. Three specimens were used for performing freezing test while the fourth was used to measure the initial condition. The side faces and bottom faces of the cylinders were thermally isolated. The top faces were covered by metal slabs, which contained circulating fluid at a temperature of 6 °C. Thus the specimens were frozen from the top down. Thermal couples were installed in the cylinders at different depths to obtain the temperature distributions. Then the specimens were removed from cooling slabs after 12, 24 and 50 h, respectively. After removing the cylinder shells, the samples were divided into 1 cm slices for obtaining water content distributions.
Table 1 Thermal properties of the materials for numerical simulation. Thermal conductivity k/(W m1 K1) Sandy loam soil Water Ice Air
1.62 0.56 2.24 0.026
Specific heat C/(J kg1 K1)
Density (kg m3)
850
2215
4180 2050 1005
1000 917 1.29
Latent heat H/(kJ kg1)
336 336
3.2.1. Soil model generation The computational domain is divided by 200 500 lattices in order to maintain the same aspect ratio as the cylinder in Mizoguchi’s experiment as shown in Fig. 4. In addition, 200 500 lattices also provide the sufficient representativeness of the soil model generated by the algorithm. In Section 2.1 Mizoguchi didn’t perform porosity test in his research, according to the exist research, the porosity of sandy loam soil is between 0.4 and 0.5 [37]. Therefore 0.45 is adopted in the current research. The value of core distribution probability ud is given as 0.35; the growing probabilities of growing cores ug,i are random in each direction of each growing core. Thus, a chaotic structure with various particle sizes is generated. Substantially, the generated structure is a 200 500 binary matrix, in which 1 represents solid particle while 0 stands for pore. 3.2.2. Initial conditions and boundary conditions The computational domain is given a uniform temperature of T0 = 6.7 °C; the initial water content is 0.33, which is initialized by providing a density ratio of qr =qr ¼ 1=0:36 to the fluid lattices of each phase. The adjustment coefficient for temperature potential is given by Gtmp = 0.2. For the fluid lattices whose value of corresponding lattices in soil model is 1, are assigned bounce-back boundary condition. Fluid lattices on boundary AB and DC are also applied bounce-back boundary condition in order to maintain mass conservation. In addition, fluid lattices on boundary AD and BC whose value of corresponding lattices in soil model is 0, are applied periodic boundary condition in order to maintain the continuity of flow.
Frozen soil
Frozen fringe
Freezing direction
3.2. Soil model generation and parameters for numerical simulation
Unfrozen soil
Fig. 6. Illustration of three regions in freezing process.
For the thermal lattices, a Dirichlet boundary condition with a constant temperature of 6 °C is applied to boundary AB as shown in Fig. 4. At boundary AB, the distribution functions g4, g7 and g8 are unknown as shown in Fig. 5. The unknown distribution functions can be determined by using the known distribution functions and the given constant temperature at the boundary. According to Eq. (20), the prescribed temperature boundary condition at boundary AB can be written as
W. Song et al. / International Journal of Heat and Mass Transfer 94 (2016) 29–38
T AB ¼
8 X
g i ¼ g 0 þ g 1 þ g 2 þ g 3 þ g 5 þ g 6 þ g 0 ðx4 þ x7 þ x8 Þ
35
ð37Þ
i¼0
where g0 represents the sum of residual distribution functions. Then g0 can be calculated by
g0 ¼
T AB ðg 0 þ g 1 þ g 2 þ g 3 þ g 5 þ g 6 Þ x4 þ x7 þ x8
ð38Þ
It is assumed that the residual distribution functions satisfy the relation
g i ¼ xi g 0
ð39Þ
Finally, the unknown distribution functions can be calculated by
x4 ¼ 19 g 0 x7 ¼ 361 g 0 x8 ¼ 361 g 0
ð40Þ
As for boundary DC, an adiabatic boundary condition is applied as shown in Fig. 4, which can be written as Fig. 7. Comparison between experimental and numerical results of water content distribution after 12 h freezing.
@T ¼0 @y y¼y0
ð41Þ
Considering the given temperature gradient at boundary DC, the temperature at boundary DC can be calculated with a three-point finite difference scheme [38], which is carried out by
T DC ¼
2 @T 4Tðy0 þ DyÞ Tðy0 þ 2DyÞ þ 3 @y y¼y0 3
Fig. 8. Comparison between experimental and numerical results of water content distribution after 24 h freezing.
Fig. 9. Comparison between experimental and numerical results of water content distribution after 50 h freezing.
Fig. 10. Water content distribution after 12 h freezing.
ð42Þ
36
W. Song et al. / International Journal of Heat and Mass Transfer 94 (2016) 29–38
Vast research has revealed that from an engineering point of view the whole freezing system can be divided into three zones as shown in Fig. 6. As mentioned in Section 2.3, water flows in both
unfrozen soil and frozen fringe are driven by the temperature gradient. As a matter of fact, freezing and drying have the same effect in regard to removing water from soil. However, freezing provides a much quicker process than drying, extremely high water potential gradients leads to very rapid upward water flow. According to Eqs. (28) and (29), the driven force in frozen fringe is one order of magnitude higher than that in unfrozen soil, which is clearly evident in Figs. 7–9. The freezing front is clearly visible in Figs. 7–9, where the water accumulates rapidly. The comparisons of experimental results and simulation results at 12, 24 and 50 h are shown in Figs. 7–12 respectively. The simulation results are obtained from three independent simulations; the soil model is regenerated in each simulation. All the simulation results are in good agreement with the experimental results. The rapid decrease of water content in frozen fringe and slow recovery in deeper position are well predicted. Moreover, by the courtesy of mesoscopic nature of LBM, the present model is not only able to provide an accurate prediction of macroscopic water content, but also offers a detailed description of how the microscopic structure affects the water distribution during the freezing process, which is absent in the existing macroscopic models. It should be noted that, despite the accuracy of the present model in predicting the water content, the size of the frozen fringe in simulation is smaller than that in the experiment. In reality, water in frozen fringe does not freeze instantly, the mushy state actually provides an extra resistance to water flow in this zone. Therefore, the water flow in the frozen fringe in particular may need further study. As mentioned previously, water flows in both unfrozen soil and frozen fringe are driven by the temperature gradient.
Fig. 11. Water content distribution after 24 h freezing.
Fig. 12. Water content distribution after 50 h freezing.
Thus, the unknown distribution functions g3, g5 and g6 at boundary DC as showed in Fig. 5 can be calculated with the same scheme as used at boundary AB. 3.2.3. Thermal properties The thermal properties for simulation are given in Table 1. It is noteworthy that in fluid region where vr < 1, the thermal conductivity, specific heat, density and latent heat need to be corrected. The correction of specific heat can be carried out by
C ¼ C r vr þ C r vr
ð43Þ
Analogously, the correction of density and latent heat can be calculated by
q ¼ qr vr þ qr vr
ð44Þ
H ¼ H r vr
ð45Þ
The correction of thermal conductivity can be expressed by Maxwell’s equation, which can be written as
k ¼ kr
2kr þ kr 2vr ðkr kr Þ 2kr þ kr þ vr ðkr kr Þ
ð46Þ
where v is obtained from Eq. (9). 4. Results and discussion
W. Song et al. / International Journal of Heat and Mass Transfer 94 (2016) 29–38
37
LBM, the proposed model provides a great potential for simulating heat and mass transfer phenomena with phase transformations in porous media.
Conflict of interest None declared. References
Fig. 13. Comparison between experimental and numerical results of temperature distribution after 12, 24 and 45 h freezing.
Consequently, to obtain an accurate temperature distribution is of equal importance with the accurate water content distribution. Comparisons of temperature distribution in experimental results and simulation results at 12, 24 and 45 h are shown in Fig. 13. It can be observed that the simulation results at 24 and 45 h are in good agreement with the experimental results; on the other hand, at 12 h the temperature of simulation in frozen soil is lower than those from experiment. Since in the present model, considering that the speed of water flow is considerably low, the heat flux caused by water flux therefore is neglected. This assumption introduces inaccuracy during the first hours of freezing when the water flow flux is relatively high. 5. Conclusions While the traditional models for solving the heat and mass transfer problem in frozen soil take macroscopic approaches, the present model adopted a mesoscopic approach. This paper presented a stochastic generation-growth method for reconstructing the soil structure with given macroscopic geometric parameters, a lattice Boltzmann model for solving multiphase fluid flow and a lattice Boltzmann model for solving heat conduction with phase change in porous media. The present model was validated by a classic experiment performed by Mizoguchi. A good agreement between numerical and experimental results was observed when comparing water content, temperature distribution and the position of the freezing front. Moreover, the result proved that the complex phenomena in frozen soil during freezing process could be divided into the interaction between several common physical phenomena and the stochastic structure of soil. However, segregation between fluid lattice and thermal lattice introduced an inaccuracy when facing the relatively high flow speed. On the other hand, this shortcoming can be overcome by introducing a macroscopic total fluid speed into the equilibrium distribution function of thermal lattice. In addition, segregation between fluid lattice and thermal lattice provides extra flexibility: the collision and streaming processes are fully local for different lattices, which means new transport phenomena can be introduced easily without too much implementation. Therefore, fluid lattice and thermal lattice can adopt different structures. For example, in the present model, both fluid and thermal lattices use D2Q9 lattice. In fact, D2Q5 lattice for thermal lattices can provide a higher performance without compromising any accuracy. Considering the mesoscopic nature of
[1] X. Xiaozu, W. Jiadeng, Z. Lixin, Physics of Frozen Soils, second ed., Science Press, Beijing, 2010. [2] P. Viterbo, A. Beljaars, J.F. Mahfouf, J. Teixeira, The representation of soil moisture freezing and its impact on the stable boundary layer, Q. J. R. Meteorol. Soc. 125 (559) (1999) 2401–2426. [3] E. Poutou, G. Krinner, C. Genthon, N. de Noblet-Ducoudré, Role of soil freezing in future boreal climate change, Clim. Dyn. 23 (6) (2004) 621–639. [4] P.M. Cox, R.A. Betts, C.D. Jones, S.A. Spall, I.J. Totterdell, Acceleration of global warming due to carbon-cycle feedbacks in a coupled climate model, Nature 408 (6809) (2000) 184–187. [5] P. Friedlingstein, L. Bopp, P. Ciais, J.L. Dufresne, L. Fairhead, H. LeTreut, P. Monfray, J. Orr, Positive feedback between future climate change and the carbon cycle, Geophys. Res. Lett. 28 (8) (2001) 1543–1546. [6] J. Cary, H. Mayland, Salt and water movement in unsaturated frozen soil, Soil Sci. Soc. Am. J. 36 (4) (1972) 549–555. [7] R. Harlan, Analysis of coupled heat-fluid transport in partially frozen soil, Water Res. Res. 9 (5) (1973) 1314–1323. [8] G.L. Guymon, J.N. Luthin, A coupled heat and moisture transport model for arctic soils, Water Res. Res. 10 (5) (1974) 995–1001. [9] M. Fuchs, G. Campbell, R. Papendick, An analysis of sensible and latent heat flow in a partially frozen unsaturated soil, Soil Sci. Soc. Am. J. 42 (3) (1978) 379–385. [10] G. Flerchinger, K. Saxton, Simultaneous heat and water model of a freezing snow-residue-soil system I. Theory and development, Trans. ASAE 32 (2) (1989) 565–571. [11] R. Jordan, A one-dimensional temperature model for a snow cover: technical documentation for SNTHERM. 89, DTIC Document, 1991. [12] L. Zhao, D. Gray, A parametric expression for estimating infiltration into frozen soils, Hydrol. Process. 11 (13) (1997) 1761–1775. [13] Q. Li, S. Sun, Y. Xue, Analyses and development of a hierarchy of frozen soil models for cold region study, J. Geophys. Res. 115 (D3) (2010). [14] P.-E. Jansson, Coupled heat and mass transfer model for soil-plant-atmosphere systems, 2001. [15] H. Engelmark, U. Svensson, Numerical modelling of phase change in freezing and thawing unsaturated soil, Nord. Hydrol. 24 (2–3) (1993) 95–110. [16] A. Slater, A. Pitman, C. Desborough, Simulation of freeze-thaw cycles in a general circulation model land surface scheme, J. Geophys. Res. (1984–2012) 103 (D10) (1998) 11303–11312. [17] K.A. Cherkauer, D.P. Lettenmaier, Hydrologic effects of frozen soils in the upper Mississippi River basin, J. Geophys. Res. (1984–2012) 104 (D16) (1999) 19599–19610. [18] D.B. Ingham, I. Pop, Transport Phenomena in Porous Media III, Elsevier, 2005. [19] K. Vafai, Handbook of Porous Media, CRC Press, 2005. [20] P. Vadász, Emerging Topics in Heat and Mass Transfer in Porous Media: From Bioengineering and Microelectronics to Nanotechnology, Springer Science Business Media, 2008. [21] L. Chen, Q. Kang, Y. Mu, Y.-L. He, W.-Q. Tao, A critical review of the pseudopotential multiphase lattice Boltzmann model: methods and applications, Int. J. Heat Mass Transf. 76 (2014) 210–236. [22] S. Chen, G.D. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30 (1998) 329–364. [23] L.P. Kadanoff, On two levels, Phys. Today 39 (1986) 7. [24] B. Yaron, R. Calvet, R. Prost, The soil as a porous medium, in: Soil Pollution, Springer, Berlin Heidelberg, 1996, pp. 3–24. [25] M. Wang, J. He, J. Yu, N. Pan, Lattice Boltzmann modeling of the effective thermal conductivity for fibrous materials, Int. J. Therm. Sci. 46 (9) (2007) 848–855. [26] S. Schlüter, H.-J. Vogel, On the reconstruction of structural and functional properties in random heterogeneous media, Adv. Water Resour. 34 (2) (2011) 314–325. [27] X. Shan, H. Chen, Lattice Boltzmann model for simulating flows with multiple phases and components, Phys. Rev. E 47 (3) (1993) 1815–1819. [28] X. Shan, H. Chen, Simulation of nonideal gases and liquid–gas phase transitions by the lattice Boltzmann equation, Phys. Rev. E 49 (4) (1994) 2941–2948. [29] X. Shan, G. Doolen, Multicomponent lattice-Boltzmann model with interparticle interaction, J. Stat. Phys. 81 (1–2) (1995) 379–393. [30] X. Shan, G. Doolen, Diffusion in a multicomponent lattice Boltzmann equation model, Phys. Rev. E 54 (4) (1996) 3614–3620. [31] J.-R.H.C. Wen-Shu Jiaung, Lattice Boltzmann method for the heat conduction problem with phase change, Numer. Heat Transf. Part B 39 (2) (2001) 167– 187.
38
W. Song et al. / International Journal of Heat and Mass Transfer 94 (2016) 29–38
[32] C. Huber, A. Parmigiani, B. Chopard, M. Manga, O. Bachmann, Lattice Boltzmann model for melting with natural convection, Int. J. Heat Fluid Flow 29 (5) (2008) 1469–1480. [33] R.W. Style, S.S.L. Peppin, The kinetics of ice-lens growth in porous media, J. Fluid Mech. 692 (2012) 482–498. [34] M. Mizoguchi, Water Heat and Salt Transport in Freezing Soil (Ph.D. thesis), University of Tokyo, Tokyo, 1990. [35] K. Hansson, J. Šimu˚nek, M. Mizoguchi, L.-C. Lundin, M.T. van Genuchten, Water flow and heat transport in frozen soil, Vadose Zone J. 3 (2) (2004) 693–704.
[36] X. Tan, W. Chen, H. Tian, J. Cao, Water flow and heat transport including ice/ water phase change in porous media: numerical simulation and application, Cold Reg. Sci. Technol. 68 (1–2) (2011) 74–84. [37] S.A. Shahid, A.A. Qidwai, F. Anwar, I. Ullah, U. Rashid, Improvement in the water retention characteristics of sandy loam soil using a newly synthesized poly (acrylamide-co-acrylic acid)/AlZnFe2O4 superabsorbent hydrogel nanocomposite material, Molecules 17 (8) (2012) 9397–9412. [38] C.-H. Liu, K.-H. Lin, H.-C. Mai, C.-A. Lin, Thermal boundary conditions for thermal lattice Boltzmann simulations, Comput. Math. Appl. 59 (7) (2010) 2178–2193.
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
A lattice Boltzmann model for heat and mass transfer phenomena with phase transformations in unsaturated soil during freezing process Wenyu Song, Yaning Zhang, Bingxi Li ⇑, Xinmeng Fan School of Energy Science and Engineering, Harbin Institute of Technology, 92 West Dazhi Street, Nan Gang District, Harbin 150001, China
a r t i c l e
i n f o
Article history: Received 18 May 2015 Received in revised form 29 October 2015 Accepted 3 November 2015
Keywords: Lattice Boltzmann method Frozen soil Stochastic structure Water content distribution
a b s t r a c t In this paper, a mesoscopic numerical model for simulating the heat and moisture transport phenomena in frozen soil during freezing process is presented. The model includes a stochastic generation-growth method for reconstructing the soil structure with given macroscopic geometric parameters, a lattice Boltzmann method (LBM) model for solving multiphase fluid flow and a LBM model for solving heat conduction with phase change in porous media. By using the present model, freezing process in sandy loam soil is demonstrated. Two types of boundary condition, prescribed temperature and adiabatic boundary conditions are employed in the given cases. All the numerical results show good agreement with experimental results. Considering the mesoscopic nature of LBM, the proposed model is a potential alternative to traditional continuum models. Ó 2015 Elsevier Ltd. All rights reserved.
1. Introduction Frozen ground region is generally distributed on earth. Permafrost region covers approximately 23% of the global land surface; seasonally frozen ground region can be found all over the regions above 24° latitude [1]. The freezing–thawing process of frozen soil has great effects on regional climate and hydrology [2,3]. Moreover, the abundant greenhouse gases in frozen soil are also important to global warming [4,5]. However, the coexistence of heat and mass transfer, phase change and multiphase flow during the freezing–thawing process makes it considerably difficult to give a simple and complete description of water content and temperature distribution in frozen soil. Therefore, a favorable heat and moisture transfer model shows its great importance to frozen soil research. The coupled heat and mass transport process was introduced into frozen soil model to solve the heat and moisture balance problems in the early 1970s [6–8]. In the earlier studies, certain processes, for instance, vapor flux and its corresponding phase change were omitted for simplicity [6–9]. During the past decades, more evidence revealed the complexity of heat and mass transport in freezing–thawing process of frozen soil [10–12]. Thus far, various numerical models with different levels of complexity have been developed. In term of governing equations, the most complex as well as the most comprehensive model includes ten variables determined by four prognostic equations and five diagnostic ⇑ Corresponding author. http://dx.doi.org/10.1016/j.ijheatmasstransfer.2015.11.008 0017-9310/Ó 2015 Elsevier Ltd. All rights reserved.
equations with one assumption [11–13]; in the relatively simple model which is the most commonly used, only two prognostic equations are used for solving temperature and liquid water content. In term of physical processes, some models take into account the vapor movement and its corresponding phase change in the water and heat balances [10,14] while others neglect those processes [15–17]. Furthermore, even with the most complex model, the estimation of ice content and position of the ice lens remains a challenging task. One reason is the highly nonlinear relationship between temperature, ice content and liquid water content [13]. On the other hand, fundamentally, soil is granular porous material. Considerable research shows that the heat and mass transport in porous material is strongly affected by its pore structure [18–20]. The underlying problems raise the need for alternative approaches to understand the connection between the macroscopic phenomena and the underlying microscopic dynamic at a more fundamental level [21]. The conventional numerical methods such as finite difference method (FDM), finite element method (FEM) and finite volume method (FVM) are based on the discretization of macroscopic continuum equations. This scheme makes conventional numerical methods face great challenges while solving flow with complex interface or flow with complex boundary. Since the late 1980s, a new mesoscopic method base on the discrete kinetic theory, named the lattice Boltzmann method (LBM), has become a promising tool for the investigations of a wide range of complex flows, including multiphase flows, porous flows, thermal flows and reactive transport [22]. The fundamental idea of LBM is quite
30
W. Song et al. / International Journal of Heat and Mass Transfer 94 (2016) 29–38
Nomenclature c cs C e f f1 F g g g0 G H L s t T u u0 DV x, y x y0
a v
Dt Dx e e0
m q
lattice speed speed of sound specific heat discrete velocity density distribution function of fluid lattice liquid fraction force density distribution function of thermal lattice gravitational acceleration sum of unknown distribution functions adjustment coefficient enthalpy latent heat tag variable time temperature velocity common velocity of the fluid mixture difference of volume axial coordinate lattice site axial coordinate of boundary thermal diffusivity mass fraction time step lattice spacing porosity given porosity kinematic viscosity density
straightforward: the macroscopic dynamics of fluid are the result of the collective behavior of many microscopic particles in the system and the macroscopic dynamics are not sensitive to the underlying details in microscopic physics [23]. Thus, LBM can be considered as a bridge sitting in the intermediate area between microscopic molecular dynamics and macroscopic fluid dynamics. The major advantage of LBM is that it behaves like a solver for the conservation equations in the bulk flow, meanwhile the kinetic equation offers many advantages of molecular dynamics, including clear physical pictures, easy implementation of complex boundary conditions and fully parallel algorithms [22]. In this work, an advanced mesoscopic numerical model based on LBM for simulating the heat and mass transfer phenomena with phase transformation in frozen soil during freezing process is presented. The fundamental idea of the proposed model is analogous to LBM: the complex phenomena in frozen soil during freezing process are caused by the interactions between several common physical phenomena and the stochastic structure of soil. The model is validated by comparing the numerical results with existed experimental results. This model should not only provide an advanced tool for the frozen soil research but also for the investigation of the similar process in generalized porous media as well. 2. Numerical models 2.1. Reconstruction of soil structure In nature, soils are heterogeneous assemblies of materials, forming porous media [24]. Various studies have shown the importance of the stochastic approaches for flow, heat transfer and solute transport in porous media [24–26]. In this paper, a
q0 s
ud ug
w w x
given density relaxation time core distribution probability growth probability effect mass water potential weight factor
Superscripts ads adhesive force b body force coh cohesive force eq equilibrium state k number of iteration tmp temperature Subscripts A air pressure i index of speed direction of lattice l liquid state P pressure s solid state S osmotic Tot total tmp temperature Z position a thermal conductivity r first phase r second phase m kinematic viscosity
stochastic generation-growth method is adopted for soil structure reconstruction based on the given porosity. The reconstruction process is conducted as follows: (i) Randomly locate the growing cores based on a given core distribution probability ud. The value of ud represents the quantity of particles in soil and further determines the size of particles at a given porosity e0. (ii) Define the growing probability of growing cores in each direction ug,i. For a D2Q9 lattice Boltzmann model, there are 8 directions on the growing cores according to the velocity directions of a D2Q9 lattice. (iii) Grow particles from the cores along the given directions; the newly grown particles are used as growing cores in the next cycle; calculate porosity e, stop growing once e meets the given porosity e0, otherwise go back to step (i). (iv) For some critical circumstances, e.g. for the permeable structure with low porosity, an extra step is executed to identify connected regions. If the model meets the requirement, stop growing, otherwise go back to step (i). Thus, for an area of a given size, the quantity of particles is determined by ud. The shape of a particle is defined by ug, isotropic distribution gives the particles round-like shape; anisotropic distribution gives the particles square-like shape, triangular-like shape or fiber-like shape as illustrated in Fig. 2. Therefore, by using a combination of different ug, a soil model with more realistic structure can be generated. Substantially, for the D2Q9 model, the generated soil structure is a D(4, 82) bond percolation network. The first row of Fig. 1 illustrates the soil structure generated in 300 300 grids while the second row shows the connectivity of
31
W. Song et al. / International Journal of Heat and Mass Transfer 94 (2016) 29–38
(a)
(b)
(c)
Fig. 1. Soil structure generated by stochastic generation-growth method and connectivity of pores. (a) e = 0.55, ud = 0.5; (b) e = 0.40, ud = 0.3; (c) e = 0.40, ud = 0.35.
Fig. 2. Structures generated by anisotropic growing probability ug.
pores with an assumption of fluid inlet from the lower boundary. In the first row, black represents solid elements and white represents fluid elements; in the second row, black represents impermeable region and white represents permeable region. It indicates that porosity e, core distribution probability ud and growing probability ug all have strong influences on fluid flow in porous media. For porous media with a low porosity, ud plays a more important role on permeability.
2.2. The Shan and Chen pseudopotential LBM model The pseudopotential model was proposed by Shan and Chen in 1993 [27–30]. Due to its simplicity and versatility, it is the most widely used multiphase LBM model [21] to the best knowledge
of the author. We adopt a single relaxation time collision operator (BGK) LB equation, which can be expressed as follows:
f r;i ðx; t þ DtÞ ¼ f r;i ðx; tÞ
1 eq f r;i ðx; tÞ f r;i ðx; tÞ
sr;m
f r;i ðx þ ei Dt; t þ DtÞ ¼ f r;i ðx; t þ DtÞ
ð1Þ ð2Þ
where f r;i ðx; tÞ is the density distribution function of the ith direction at the lattice site x and time t (see Fig. 3); feq represents the equilibrium distribution function; c = Dx/Dt is the lattice speed, in which Dx = Dt = 1 represents the lattice spacing and time step respectively; sm is the dimensionless relaxation time. In the computer algorithm, the collision step is represented by Eq. (1), which means the system reaches the local Maxwellian equilibrium on a time scale s. Eq. (2) represents the streaming step in which
32
W. Song et al. / International Journal of Heat and Mass Transfer 94 (2016) 29–38
fluid particles are streamed from one lattice to a neighboring lattice according to the discrete velocity ei in the lattice. For a D2Q9 lattice,
8 ð0;0Þ i¼0 > > > < i ¼ 1;2;3;4 ei ¼ ðcos ½ða 1Þp=2;sin ½ða 1Þp=2Þ > > > : pffiffiffi 2ðcos ½ða 5Þp=2þ p=4;sin ½ða 5Þp=2þ p=4Þ i ¼ 5;6;7;8
ð3Þ The equilibrium distribution function feq is given by
3 9 3 eq 2 2 f r;i ¼ xi qr 1 þ 2 ðei ueq Þ þ 4 ðei ueq Þ 2 ðueq Þ c 2c 2c
ð4Þ
mr ¼
sr;m 0:5ÞDt
ð5Þ
ð6Þ
X f r;i
ð7Þ
i
qr ur ¼
X f r;i ei
ð8Þ
i
the mass fraction v of phase r is calculated by
vr ¼
qr
In order to introduce the force term (e.g. fluid–fluid interaction force, fluid–solid interaction force, gravitational force and other bulk forces), a velocity shift scheme [27] is adopted. In this scheme, the parameter ueq r in Eq. (4) is determined by the relation
ueq ¼ u0 þ
sr;m Fr qr
ð10Þ
where u0 is the common velocity of the fluid mixture, which is given by
P u0 ¼
qr ur r s P qrr;m r sr;m
ð11Þ
In Eq. (10), Fr is the sum of all interactive forces between fluid ads b tmp coh particles. In this model, Fr ¼ Fcoh r þ Fr þ Fr þ Fr , Fr represents the cohesive force between particles from the same phase; Fads r stands for the adhesive force between fluid particles and solid boundary; Fbr is the sum of bulk forces such as gravitational force and buoyancy force; Ftmp r represents the bulk force caused by temperature potential change. In the current model, vapor in soil is neglected. Therefore the fluid system can be treated as a binary mixture of water and air. Only repulsion between different components is taken into account, the interactive force applied on the particles of one component is carried out by coh Fcoh r ðxÞ ¼ wr ðx; tÞG
X
xa wr ðx þ ei Þei
where si works as a tag, the value of si is
si ðx þ ei Þ ¼
0
node x þ ei is fluid lattice node
ð14Þ
1 node x þ ei is solid boundary
Similar to Eq. (12), Gads r adjusts the interaction force between the particles of component r and solid boundary. For wetting fluid, ads Gads r is chosen positive while for non-wetting fluid, Gr is chosen negative. The effect mass is
q w ¼ q0 1 exp
ð15Þ
Fr ¼ qr g b
ð16Þ
In frozen soil, the moisture moves along the temperature gradient from high temperature to low temperature during both freezing and thawing processes. Therefore, an accurate numerical model to simulate the temperature distribution in frozen soil is of great importance. In LB model, the heat conduction process is performed by applying another layer of thermal lattice. Similar to Section 2.2, a D2Q9 lattice is adopted; the LB equation for thermal diffusion can be expressed as follows:
g i ðx þ ei Dt; t þ DtÞ ¼ g i ðx; tÞ ð9Þ
qr þ qr
ð13Þ
2.3. LBM model for heat conduction problem with phase change
in which cs represents the lattice sound speed. The macroscopic density and momentum at each lattice node can be calculated as
qr ¼
xi si ðx þ ei Þei
i
Then, the gravitational force of component r is given by
The viscosity in lattice unit is given by
c2s ð
X
q0
where in Eq. (4), the weight factor xi is given by
8 i¼0 > < 4=9 xi ¼ 1=9 i ¼ 1; 2; 3; 4 > : 1=36 i ¼ 5; 6; 7; 8
ads Fads r ðxÞ ¼ wr ðx; tÞG
ð12Þ
i
represents the other component of the fluid mixture; Gcoh where r r is a adjuster to tune the surface tension between two components. The cohesive force between fluid particles and solid boundary can be expressed analogously by
1
sa
g i ðx; tÞ g eq i ðx; tÞ
ð17Þ
The discrete velocity ei in the lattice is expressed in Eq. (3). The equilibrium distribution function g eq i is given by
g eq i ¼ xi T
ð18Þ
The weight factor xi is given in Eq. (5). In analogy with the kinematic viscosity m, the thermal diffusivity a is calculated as
a ¼ c2s ðsa 0:5ÞDt
ð19Þ
The macroscopic temperature can be calculated as the 0th order momentum of gi as
X gi ¼ T
ð20Þ
i
In order to solve the heat conduction problem with phase change, Jiaung [31] proposed a LB model with a source term in Eq. (17), which can be written as
g i ðx þ ei Dt; t þ DtÞ ¼ g i ðx; tÞ þ xi
1
sa
g i ðx; tÞ g eq i ðx; tÞ
L
f ðt þ DtÞ f l;i ðtÞ C l;i
ð21Þ
In order to solve the source term, a sub-loop is introduced into every time step, which is executed as follows: (i) The value of the particle distribution function at the new time level t + Dt in the kth iteration of sub-loop is calculated according to Eq. (21) as
g ki ðx þ ei Dt; t þ DtÞ ¼ g i ðx; tÞ þ xi
1
sa
g i ðx; tÞ g eq i ðx; tÞ
L k1 f l;i ðt þ DtÞ f l;i ðtÞ C
ð22Þ
where the superscripts k and k1 represent the current and previous iterations, respectively.
W. Song et al. / International Journal of Heat and Mass Transfer 94 (2016) 29–38
2
6
33
respectively. For low salinity soil under standard atmospheric pressure, WS and WA can be neglected, hence
5
WTot ¼ WP þ WZ þ Wtmp
ð28Þ
ads In the LBM model, Wp is carried out by Fcoh r and Fr (caused by
0
3
capillary action); WZ is represented by Fbr and Wtmp is represented by Ftmp r . When the soil temperature is above freezing temperature, the gradient of water potential caused by temperature change can be derived from Clapeyron’s equation, which can be written as
1
rW ¼
7
4
(ii) The macroscopic temperature in the kth iteration is calculated by Eq. (20). (iii) The total enthalpies in the kth iteration is calculated as k1
ð23Þ
L
(iv) The liquid fraction in the current iteration is calculated as k
fl ¼
8 > < k0
H Hs > Hl Hs
:
1
H < Hs Hs 6 H 6 Hl
ð24Þ
H > Hl
where Hs and Hl represent the enthalpies at the freezing and melting temperatures, respectively. (v) Repeat steps (i)–(iv) until convergence is reached, the convergence criterion is expressed as
! f k f k1 T k T k1 8 min k1 ; < 10 f Tk
ð25Þ
rW ¼
ð26Þ
in which the dimensionless scale T* = (T T0)/(T1 T0); t* = at/l2; l is the natural length of the system; T0 and T1 are the initial temperature and the temperature applied on the boundary, respectively; St = C(T0 T1)/L is the Stefan number. Thus the total time steps needed for the simulation can be estimated by using Eq. (26).
Ls rT TjDV s j
ð30Þ
where Ls and DVs represent the latent heat of fusion and specific volume change during freezing. Thus the movement of water in soil can be quantified by rW. It is noted that In Eq. (30), Ls and DVs are known as the property of water, T = 273 K is the freezing point of water in Kelvin. Therefore, Eq. (30) becomes
rW ¼ 13:48 106 rT
ð31Þ
The temperature water potential has the same value with the suction pressure caused by freezing process can be written as
Pr ¼
Ls q u2 ¼ r r TjDV s j 2
ð32Þ
in lattice unit the suction pressure is also carried out by
P¼
In addition, when the liquid fraction fl in a thermal lattice reaches 0, the corresponding fluid lattice is converted to ‘‘Bounce back” lattice, which ensures the conservation of ice content. It is noteworthy that the number of thermal diffusivity in real physical units is so small that if we introduce such small thermal diffusivity the code will blow up. However, by using the total time of the evolution tmax as the timescale, we can get a nondimensional version of the differential equation of Eq. (21) which proposed by Huber et al. [32] is given by
@T at max 2 1 @f l ¼ 2 r T St @t @t l
ð29Þ
where Ll and DVl represent the latent heat of condensation and specific volume change during vapor condensation, respectively. Analogously, when the soil temperature reaches the freezing temperature, the gradient of water potential caused by temperature can be carried out by
8
Fig. 3. The layout of velocity directions of a D2Q9 lattice.
Hk ¼ CT k þ f l
Ll rT T DV l
qu2
ð33Þ
2
hence, the suction pressure in lattice unit becomes
P¼
qu2 P qr u2r r
ð34Þ
we have known that
ur csr ¼ u cs
ð35Þ
assuming that in soil, the real speed of sound csr = 150 m/s, the real density of water qr = 1000 kg/m3, it can be calculated that P = 0.199. However, for the real soil, Eq. (31) needs to be corrected. The correction factor depends on the type of soil [33]. Here we propose a relation to implement rW in lattices tmp F tmp r ¼ G rWDx tmp
where G
ð36Þ
is determined by the type of soil.
3. Verification of LBM model, results and discussion 2.4. The implementation of soil water potential in LBM model 3.1. Introduction of laboratory experiment In soil, water moves from one position to another due to gravity, pressure, capillary action, osmosis and so on. The tendency of water movement is quantified by water potential W. The total water potential WTot is given by
WTot ¼ WP þ WZ þ WS þ WA þ Wtmp
ð27Þ
where WP, WZ, WS, WA and Wtmp stand for pressure, gravitational, solute (osmotic), air pressure and temperature potentials,
Mizoguchi performed a classic laboratory experiment in 1990 [34], the result has been adopted in various studies [35,36]. In his research, Kanagawa sandy loam soil was packed in four identical cylinders. The dimensions of all four cylinders were identical, which was 20 cm long with an internal diameter of 8 cm. The specimens were prepared with the same initial conditions, which involved a uniform temperature of 6.7 °C and close volumetric
34
W. Song et al. / International Journal of Heat and Mass Transfer 94 (2016) 29–38
2
6
A
3
5
7
4
8
6
2
5
y
B
1
0
0
D
x Fig. 4. Boundary conditions in Mizoguchi’s experiment and in numerical simulation.
3
1
7
C
8
4
Fig. 5. Illustration of unknown distribution functions at boundaries AB and DC.
water content of 0.33. Three specimens were used for performing freezing test while the fourth was used to measure the initial condition. The side faces and bottom faces of the cylinders were thermally isolated. The top faces were covered by metal slabs, which contained circulating fluid at a temperature of 6 °C. Thus the specimens were frozen from the top down. Thermal couples were installed in the cylinders at different depths to obtain the temperature distributions. Then the specimens were removed from cooling slabs after 12, 24 and 50 h, respectively. After removing the cylinder shells, the samples were divided into 1 cm slices for obtaining water content distributions.
Table 1 Thermal properties of the materials for numerical simulation. Thermal conductivity k/(W m1 K1) Sandy loam soil Water Ice Air
1.62 0.56 2.24 0.026
Specific heat C/(J kg1 K1)
Density (kg m3)
850
2215
4180 2050 1005
1000 917 1.29
Latent heat H/(kJ kg1)
336 336
3.2.1. Soil model generation The computational domain is divided by 200 500 lattices in order to maintain the same aspect ratio as the cylinder in Mizoguchi’s experiment as shown in Fig. 4. In addition, 200 500 lattices also provide the sufficient representativeness of the soil model generated by the algorithm. In Section 2.1 Mizoguchi didn’t perform porosity test in his research, according to the exist research, the porosity of sandy loam soil is between 0.4 and 0.5 [37]. Therefore 0.45 is adopted in the current research. The value of core distribution probability ud is given as 0.35; the growing probabilities of growing cores ug,i are random in each direction of each growing core. Thus, a chaotic structure with various particle sizes is generated. Substantially, the generated structure is a 200 500 binary matrix, in which 1 represents solid particle while 0 stands for pore. 3.2.2. Initial conditions and boundary conditions The computational domain is given a uniform temperature of T0 = 6.7 °C; the initial water content is 0.33, which is initialized by providing a density ratio of qr =qr ¼ 1=0:36 to the fluid lattices of each phase. The adjustment coefficient for temperature potential is given by Gtmp = 0.2. For the fluid lattices whose value of corresponding lattices in soil model is 1, are assigned bounce-back boundary condition. Fluid lattices on boundary AB and DC are also applied bounce-back boundary condition in order to maintain mass conservation. In addition, fluid lattices on boundary AD and BC whose value of corresponding lattices in soil model is 0, are applied periodic boundary condition in order to maintain the continuity of flow.
Frozen soil
Frozen fringe
Freezing direction
3.2. Soil model generation and parameters for numerical simulation
Unfrozen soil
Fig. 6. Illustration of three regions in freezing process.
For the thermal lattices, a Dirichlet boundary condition with a constant temperature of 6 °C is applied to boundary AB as shown in Fig. 4. At boundary AB, the distribution functions g4, g7 and g8 are unknown as shown in Fig. 5. The unknown distribution functions can be determined by using the known distribution functions and the given constant temperature at the boundary. According to Eq. (20), the prescribed temperature boundary condition at boundary AB can be written as
W. Song et al. / International Journal of Heat and Mass Transfer 94 (2016) 29–38
T AB ¼
8 X
g i ¼ g 0 þ g 1 þ g 2 þ g 3 þ g 5 þ g 6 þ g 0 ðx4 þ x7 þ x8 Þ
35
ð37Þ
i¼0
where g0 represents the sum of residual distribution functions. Then g0 can be calculated by
g0 ¼
T AB ðg 0 þ g 1 þ g 2 þ g 3 þ g 5 þ g 6 Þ x4 þ x7 þ x8
ð38Þ
It is assumed that the residual distribution functions satisfy the relation
g i ¼ xi g 0
ð39Þ
Finally, the unknown distribution functions can be calculated by
x4 ¼ 19 g 0 x7 ¼ 361 g 0 x8 ¼ 361 g 0
ð40Þ
As for boundary DC, an adiabatic boundary condition is applied as shown in Fig. 4, which can be written as Fig. 7. Comparison between experimental and numerical results of water content distribution after 12 h freezing.
@T ¼0 @y y¼y0
ð41Þ
Considering the given temperature gradient at boundary DC, the temperature at boundary DC can be calculated with a three-point finite difference scheme [38], which is carried out by
T DC ¼
2 @T 4Tðy0 þ DyÞ Tðy0 þ 2DyÞ þ 3 @y y¼y0 3
Fig. 8. Comparison between experimental and numerical results of water content distribution after 24 h freezing.
Fig. 9. Comparison between experimental and numerical results of water content distribution after 50 h freezing.
Fig. 10. Water content distribution after 12 h freezing.
ð42Þ
36
W. Song et al. / International Journal of Heat and Mass Transfer 94 (2016) 29–38
Vast research has revealed that from an engineering point of view the whole freezing system can be divided into three zones as shown in Fig. 6. As mentioned in Section 2.3, water flows in both
unfrozen soil and frozen fringe are driven by the temperature gradient. As a matter of fact, freezing and drying have the same effect in regard to removing water from soil. However, freezing provides a much quicker process than drying, extremely high water potential gradients leads to very rapid upward water flow. According to Eqs. (28) and (29), the driven force in frozen fringe is one order of magnitude higher than that in unfrozen soil, which is clearly evident in Figs. 7–9. The freezing front is clearly visible in Figs. 7–9, where the water accumulates rapidly. The comparisons of experimental results and simulation results at 12, 24 and 50 h are shown in Figs. 7–12 respectively. The simulation results are obtained from three independent simulations; the soil model is regenerated in each simulation. All the simulation results are in good agreement with the experimental results. The rapid decrease of water content in frozen fringe and slow recovery in deeper position are well predicted. Moreover, by the courtesy of mesoscopic nature of LBM, the present model is not only able to provide an accurate prediction of macroscopic water content, but also offers a detailed description of how the microscopic structure affects the water distribution during the freezing process, which is absent in the existing macroscopic models. It should be noted that, despite the accuracy of the present model in predicting the water content, the size of the frozen fringe in simulation is smaller than that in the experiment. In reality, water in frozen fringe does not freeze instantly, the mushy state actually provides an extra resistance to water flow in this zone. Therefore, the water flow in the frozen fringe in particular may need further study. As mentioned previously, water flows in both unfrozen soil and frozen fringe are driven by the temperature gradient.
Fig. 11. Water content distribution after 24 h freezing.
Fig. 12. Water content distribution after 50 h freezing.
Thus, the unknown distribution functions g3, g5 and g6 at boundary DC as showed in Fig. 5 can be calculated with the same scheme as used at boundary AB. 3.2.3. Thermal properties The thermal properties for simulation are given in Table 1. It is noteworthy that in fluid region where vr < 1, the thermal conductivity, specific heat, density and latent heat need to be corrected. The correction of specific heat can be carried out by
C ¼ C r vr þ C r vr
ð43Þ
Analogously, the correction of density and latent heat can be calculated by
q ¼ qr vr þ qr vr
ð44Þ
H ¼ H r vr
ð45Þ
The correction of thermal conductivity can be expressed by Maxwell’s equation, which can be written as
k ¼ kr
2kr þ kr 2vr ðkr kr Þ 2kr þ kr þ vr ðkr kr Þ
ð46Þ
where v is obtained from Eq. (9). 4. Results and discussion
W. Song et al. / International Journal of Heat and Mass Transfer 94 (2016) 29–38
37
LBM, the proposed model provides a great potential for simulating heat and mass transfer phenomena with phase transformations in porous media.
Conflict of interest None declared. References
Fig. 13. Comparison between experimental and numerical results of temperature distribution after 12, 24 and 45 h freezing.
Consequently, to obtain an accurate temperature distribution is of equal importance with the accurate water content distribution. Comparisons of temperature distribution in experimental results and simulation results at 12, 24 and 45 h are shown in Fig. 13. It can be observed that the simulation results at 24 and 45 h are in good agreement with the experimental results; on the other hand, at 12 h the temperature of simulation in frozen soil is lower than those from experiment. Since in the present model, considering that the speed of water flow is considerably low, the heat flux caused by water flux therefore is neglected. This assumption introduces inaccuracy during the first hours of freezing when the water flow flux is relatively high. 5. Conclusions While the traditional models for solving the heat and mass transfer problem in frozen soil take macroscopic approaches, the present model adopted a mesoscopic approach. This paper presented a stochastic generation-growth method for reconstructing the soil structure with given macroscopic geometric parameters, a lattice Boltzmann model for solving multiphase fluid flow and a lattice Boltzmann model for solving heat conduction with phase change in porous media. The present model was validated by a classic experiment performed by Mizoguchi. A good agreement between numerical and experimental results was observed when comparing water content, temperature distribution and the position of the freezing front. Moreover, the result proved that the complex phenomena in frozen soil during freezing process could be divided into the interaction between several common physical phenomena and the stochastic structure of soil. However, segregation between fluid lattice and thermal lattice introduced an inaccuracy when facing the relatively high flow speed. On the other hand, this shortcoming can be overcome by introducing a macroscopic total fluid speed into the equilibrium distribution function of thermal lattice. In addition, segregation between fluid lattice and thermal lattice provides extra flexibility: the collision and streaming processes are fully local for different lattices, which means new transport phenomena can be introduced easily without too much implementation. Therefore, fluid lattice and thermal lattice can adopt different structures. For example, in the present model, both fluid and thermal lattices use D2Q9 lattice. In fact, D2Q5 lattice for thermal lattices can provide a higher performance without compromising any accuracy. Considering the mesoscopic nature of
[1] X. Xiaozu, W. Jiadeng, Z. Lixin, Physics of Frozen Soils, second ed., Science Press, Beijing, 2010. [2] P. Viterbo, A. Beljaars, J.F. Mahfouf, J. Teixeira, The representation of soil moisture freezing and its impact on the stable boundary layer, Q. J. R. Meteorol. Soc. 125 (559) (1999) 2401–2426. [3] E. Poutou, G. Krinner, C. Genthon, N. de Noblet-Ducoudré, Role of soil freezing in future boreal climate change, Clim. Dyn. 23 (6) (2004) 621–639. [4] P.M. Cox, R.A. Betts, C.D. Jones, S.A. Spall, I.J. Totterdell, Acceleration of global warming due to carbon-cycle feedbacks in a coupled climate model, Nature 408 (6809) (2000) 184–187. [5] P. Friedlingstein, L. Bopp, P. Ciais, J.L. Dufresne, L. Fairhead, H. LeTreut, P. Monfray, J. Orr, Positive feedback between future climate change and the carbon cycle, Geophys. Res. Lett. 28 (8) (2001) 1543–1546. [6] J. Cary, H. Mayland, Salt and water movement in unsaturated frozen soil, Soil Sci. Soc. Am. J. 36 (4) (1972) 549–555. [7] R. Harlan, Analysis of coupled heat-fluid transport in partially frozen soil, Water Res. Res. 9 (5) (1973) 1314–1323. [8] G.L. Guymon, J.N. Luthin, A coupled heat and moisture transport model for arctic soils, Water Res. Res. 10 (5) (1974) 995–1001. [9] M. Fuchs, G. Campbell, R. Papendick, An analysis of sensible and latent heat flow in a partially frozen unsaturated soil, Soil Sci. Soc. Am. J. 42 (3) (1978) 379–385. [10] G. Flerchinger, K. Saxton, Simultaneous heat and water model of a freezing snow-residue-soil system I. Theory and development, Trans. ASAE 32 (2) (1989) 565–571. [11] R. Jordan, A one-dimensional temperature model for a snow cover: technical documentation for SNTHERM. 89, DTIC Document, 1991. [12] L. Zhao, D. Gray, A parametric expression for estimating infiltration into frozen soils, Hydrol. Process. 11 (13) (1997) 1761–1775. [13] Q. Li, S. Sun, Y. Xue, Analyses and development of a hierarchy of frozen soil models for cold region study, J. Geophys. Res. 115 (D3) (2010). [14] P.-E. Jansson, Coupled heat and mass transfer model for soil-plant-atmosphere systems, 2001. [15] H. Engelmark, U. Svensson, Numerical modelling of phase change in freezing and thawing unsaturated soil, Nord. Hydrol. 24 (2–3) (1993) 95–110. [16] A. Slater, A. Pitman, C. Desborough, Simulation of freeze-thaw cycles in a general circulation model land surface scheme, J. Geophys. Res. (1984–2012) 103 (D10) (1998) 11303–11312. [17] K.A. Cherkauer, D.P. Lettenmaier, Hydrologic effects of frozen soils in the upper Mississippi River basin, J. Geophys. Res. (1984–2012) 104 (D16) (1999) 19599–19610. [18] D.B. Ingham, I. Pop, Transport Phenomena in Porous Media III, Elsevier, 2005. [19] K. Vafai, Handbook of Porous Media, CRC Press, 2005. [20] P. Vadász, Emerging Topics in Heat and Mass Transfer in Porous Media: From Bioengineering and Microelectronics to Nanotechnology, Springer Science Business Media, 2008. [21] L. Chen, Q. Kang, Y. Mu, Y.-L. He, W.-Q. Tao, A critical review of the pseudopotential multiphase lattice Boltzmann model: methods and applications, Int. J. Heat Mass Transf. 76 (2014) 210–236. [22] S. Chen, G.D. Doolen, Lattice Boltzmann method for fluid flows, Annu. Rev. Fluid Mech. 30 (1998) 329–364. [23] L.P. Kadanoff, On two levels, Phys. Today 39 (1986) 7. [24] B. Yaron, R. Calvet, R. Prost, The soil as a porous medium, in: Soil Pollution, Springer, Berlin Heidelberg, 1996, pp. 3–24. [25] M. Wang, J. He, J. Yu, N. Pan, Lattice Boltzmann modeling of the effective thermal conductivity for fibrous materials, Int. J. Therm. Sci. 46 (9) (2007) 848–855. [26] S. Schlüter, H.-J. Vogel, On the reconstruction of structural and functional properties in random heterogeneous media, Adv. Water Resour. 34 (2) (2011) 314–325. [27] X. Shan, H. Chen, Lattice Boltzmann model for simulating flows with multiple phases and components, Phys. Rev. E 47 (3) (1993) 1815–1819. [28] X. Shan, H. Chen, Simulation of nonideal gases and liquid–gas phase transitions by the lattice Boltzmann equation, Phys. Rev. E 49 (4) (1994) 2941–2948. [29] X. Shan, G. Doolen, Multicomponent lattice-Boltzmann model with interparticle interaction, J. Stat. Phys. 81 (1–2) (1995) 379–393. [30] X. Shan, G. Doolen, Diffusion in a multicomponent lattice Boltzmann equation model, Phys. Rev. E 54 (4) (1996) 3614–3620. [31] J.-R.H.C. Wen-Shu Jiaung, Lattice Boltzmann method for the heat conduction problem with phase change, Numer. Heat Transf. Part B 39 (2) (2001) 167– 187.
38
W. Song et al. / International Journal of Heat and Mass Transfer 94 (2016) 29–38
[32] C. Huber, A. Parmigiani, B. Chopard, M. Manga, O. Bachmann, Lattice Boltzmann model for melting with natural convection, Int. J. Heat Fluid Flow 29 (5) (2008) 1469–1480. [33] R.W. Style, S.S.L. Peppin, The kinetics of ice-lens growth in porous media, J. Fluid Mech. 692 (2012) 482–498. [34] M. Mizoguchi, Water Heat and Salt Transport in Freezing Soil (Ph.D. thesis), University of Tokyo, Tokyo, 1990. [35] K. Hansson, J. Šimu˚nek, M. Mizoguchi, L.-C. Lundin, M.T. van Genuchten, Water flow and heat transport in frozen soil, Vadose Zone J. 3 (2) (2004) 693–704.
[36] X. Tan, W. Chen, H. Tian, J. Cao, Water flow and heat transport including ice/ water phase change in porous media: numerical simulation and application, Cold Reg. Sci. Technol. 68 (1–2) (2011) 74–84. [37] S.A. Shahid, A.A. Qidwai, F. Anwar, I. Ullah, U. Rashid, Improvement in the water retention characteristics of sandy loam soil using a newly synthesized poly (acrylamide-co-acrylic acid)/AlZnFe2O4 superabsorbent hydrogel nanocomposite material, Molecules 17 (8) (2012) 9397–9412. [38] C.-H. Liu, K.-H. Lin, H.-C. Mai, C.-A. Lin, Thermal boundary conditions for thermal lattice Boltzmann simulations, Comput. Math. Appl. 59 (7) (2010) 2178–2193.