CC BY-NC-ND
0THEORETICAL & APPLIED MECHANICS LETTERS 2, 022003 (2011)
Migration-collision scheme of Lattice Boltzmann method for heat conduction problems involving solidification
Xia Li,1 Yang Liu,2 Yousheng Xu,1, a) MeiyingYe,1 and Fengmin Wu1
Department of Physics, Zhejiang Normal University, Jinhua 321004, China ^Department of Mechanical Engineering, The Hong Kong Polytechnic University, Kowloon, Hong Kong, China
(Received 01 November 2010; accepted 26 January 2011; published online 10 March 2011)
Abstract The phase-change problem is solved by the migration-collision scheme of lattice Boltzmann method. After formula derivation, we find that this method can give a rigorous numerical value for the phase-change temperature, which is of crucial importance. One-dimensional solidification in half-space and two-dimensional solidification in a corner are simulated. The phase change temperature and the liquid-solid interface are both obtained, and the results conform to the analytical solution. ( 2011 The Chinese Society of Theoretical and Applied Mechanics. [doi:10.1063/2.1102203]
Keywords Lattice Boltzmann method, migration-collision scheme, collision-migration scheme
Matter phase transition is closely related to industrial production and daily life. It is a nonlinear problem1-4 and thus should be settled by a high-efficiency and reliable numerical method. Lattice Boltzmann method (LBM) is a relatively new approach, which is applied widely to the numerical simulation of micro fluidic systems.5-10 This paper is dedicated to the simulation of pure substance phase transition via LBM.
Wolf-Gladrow11 firstly derived the LB equation for diffusion to solve heat conduction problem. Because the melting and solidification were analogous to exothermic and endothermic chemical reactions, De Fabritiis et al.12 proposed a generalized mesoscopic LB model for describing flows with solid-liquid phase transitions, which was represented by a chemical term. They used two types of quasi-particles for liquid and solid phases, respectively. Miller et al.13 simplified and extended their model subsequently. In 2001, Jiaung et al.14 firstly extended Lattice Boltzmann equation with an enthalpy formulation for treating the phase change region. Later, the proposed model was utilized to simulate the melting process of a generic laser surface in the two-dimensional and three-dimensional models by Chatterjee et al..15'16 And Huber et al.17 developed the lattice Boltzmann method coupled with natural convection. However, the extended lattice Boltzmann method developed by Jiaung14 had two problems. Firstly, it was not able to give the precise phase change temperature while handling the phase transition region. Secondly, the temperature of the point adjacent to the phase change region is sometimes higher than the temperature of surrounding points. The present work is dedicated to improving the numerical computation method for the phase change model. The present method can remedy the defects and make the numerical results preferably close to the truth.
Generally, the thermal diffusion equation of total enthalpy for phase change problem can be expressed as
= V- (kVT).
a)Corresponding author. Email: ysyou@zjnu.cn.
The total enthalpy H has two parts, namely, sensible enthalpy CT and latent enthalpy L. So H can be written as
H = CT + fi - L,
where fl is the liquid fraction and C is the constant pressure specific heat.
It is known to all that the liquid fraction fi has different values in different states of matter. fi is zero for the solid state and unity for the liquid state, respectively, and the value of fl is between 0 and 1 in mushy state (phase change region). So the total enthalpy in energy equation can be resolved into
CT + L, CTm + fl CT,
T=T T < T
where Tm is the phase-change temperature. Substitution Eq. (3) into Eq. (1) yields
dT/dt = a • V2T, T > Tm or T < Tm,
(L/C)(dfi/dt) = a •V2T, T = Tm. (4)
Equation (4) demonstrates that the heat source term (L/C)(dfi/dt) takes effect only in phase change region. Here a is the thermal diffusivity and a = k/pC. And the liquid fraction is defined as14
H - H s Hi — H s '
where Hs and Hl represent the enthalpy values at the beginning and ending of solidification respectively.
The discrete LBE was firstly written to describe the dynamics of the distribution function in the lattice gas automaton (LGA). Then a Bhatnagar-Gross-Krook (BGK) collision model was later adopted to address the complicated collision term. Recently, it has been argued that the LBE method can be derived from the continuum Boltzmann equation with BGK collision model.18
The Boltzmann BGK equation is
f + e -Vf = - I(f - f(0)).
In order to describe the evolution of microscopic particle distribution function fi(x,t), the BGK equation is discredited.
The equilibrium state particle distribution function is defined as
f(eq)(x,t) = uT (x,t),
where are the weight factors.
So the discrete BGK equation is
dfi (x,t) dt
- -[fi(x,t) - f(eq)(x,t)].
+ Ci ■ fi(x,t)
f(°\x + Ci At,t + At)].
And this scheme is equivalent to the collision-migration scheme, because Eq. (11) can be written as
fi(x') = fi(x>,t')-T[fi(x>,t)-f(0\x',t')], (12)
where x' represents (x + ci ■ At), t' represents t + At. After second migration
f i(x' + Ci ■ At, t' + At) f (x' ,t' ) - 1[fi (x' ,t' ) - fV ,t' )].
By using x, t, f to substitute x', t', f', Eq. (13) is changed to
fi(x + Ci ■ At, t + At)
=fi(x,t) --[fi(x,t) - f(0)(x,t)}.
Regarding phase change, the heat source term
L dfi L f - fi
= = —41 = Ui—Jl - Jl i i i C dt i C At
is introduced to the LB equation.14 The LB equation of (14) is extended as
The temperature of the phase change is defined as
T (x,t) = J2 fi(x,t).
The next step is the time discretion of LBM. Collision and migration are two relatively independent processes of particles motion, and LBM has thus an evolutionary process, like ... collision-migrationcollision-migration ..., in all time discrete sequence. So the collision-migration scheme and migration-collision scheme are equivalent. After time discretion, the migration-collision scheme is selected. The evolutionary process is divided into two steps
fi (x + a ■ At, t + At) = fi(x,t), (10)
the particles move along the graticules with velocity
ci = Ax/At
fi(x + a ■ At, t + At) =fi(x + Ci ■ At,t + At) - -f (x + Ci ■ At,t + At)-
fi(x + a ■ At, t + At) =fi(x,t) - A [fi(x,t) - ffW)] - At^i. (15)
The LB equation of (15) was proven correctly with Chapman-Enskog Expansion by Jiaung.14
Migration-collision scheme and collision-migration scheme have the same physical significance, but there are some differences in their specific methods for the computation of phase-change region. The particles of the whole system without external temperature were at dynamic balance in mesoscope. In this state, migration and collision happened simultaneously. However, fixing a temperature on the wall of the system, the particles on the boundary will migrate and transfer temperature into the interior. Because of the temperature change, each particle starts to collide. Hence the migrationcollision scheme fits the fact. The essential difference between the two schemes is how to process the latent heat.
In the migration-collision scheme, after migration, the temperature T inside of the model become
fi(x,t + At).
If T' < Tm in the x grid, it has started to change the phase, so
H = CT' + ft ■ L (17)
and if Hs < H <Hi
H - H s Hi — H s
H - H s L
CT ' + fi ■ L - CTm L
Then after collision f' (x,t + At)
fi(x,t + At)--[fi(x,t + At)-
f(0)(x,t + At)] - Atui Ct-L.
Then the summation over all directions of every term of Eq. (19), the xth grid temperature is obtained as
T =£ f',(x,t + At)
T' - 0 - C
--T' - T' + Tm
CT ' + fi ■ L - CTm
-m •
Here the temperature of the xth grid becomes the phase change temperature exactly. Because the temperature of T - Tm has been compensated to the latent heat by colliding in Eq. (19). On the other hand, the collisionmigration scheme cannot give the exact value of phase
change temperature. Furthermore, since the first step is collision, the latent heat will migrate to the neighbour girds in the change phase region, the neighbour girds will sometimes have a higher temperature than surrounding.
The boundary condition of the phase change is fixed by particles flowing through grids of the boundary in LBM. The boundary condition of Jiaung14 is employed in this work. Take the one-dimensional lattice for ex-ample(shown in Fig. 1).
Boundary
J Ax/2 « ]
Fig. 1. The distribution function fj (T) is fixed if the temperature on the boundary is kept at Tj .
If the distribution function fj (T) at the boundary lattice j is known, f2j-+1 can be identified as
f2,j+1 = f1,j,
flj = fj (T ) + ^ fj (T ) - f2
(21) (22)
where f2j- denotes the particle migrating from lattice j + 1. In the two-dimensional lattice, subscript 1 and 2 represent two opposite directions.
To test and verify the accuracy of the modified scheme, two models (solidification in a half -space and solidification from a corner in a quarter-space) are imitated. It is found that the migration-collision scheme describes the models well. And the computational results of the migration-collision scheme is in better agreement with the analytical solutions than Jiaung's.14
Model 1 is solidification in a half -space. Initially, the semi-infinite system is filled with the liquid at temperature T^. T^ is set to be 0.3, which is higher than the phase-change temperature Tm = 0. This is a one-dimensional model. The temperature on the left wall is kept at T0 = -1. The diagram of the model is shown in Fig. 2.
As the solidification begins, the low temperature of the left boundary is going to be transferred into the system. A liquid-solid interface is generated. The temperatures in solid and liquid regions are both variable and different from Tm, whose analytical solutions are obtained1
[Ts(x,t) - To]/[Tm - To] = [erf (x/2^ast)]/[erf (A)], [Ti (x,t) Tœ] / [Tm Tœ ] = [erfc(x/2yOlt )]/
[erfc(A^as/ai )].
The position of the liquid-solid interface is defined as
M (t) = 2\Vatt, (24)
where A is a parameter of the transcendental equation.14 The Stefan number is the ratio of sensible heat over latent heat, written as St = pCs(T0 — Tm)/L.
The relaxation time of every state is gained from
T = Ts + (t - Ts )fl.
Fig. 2. Schematic of solidification in a half-space.
Figures 3 and 4 are diagrams of the temperature distribution and the interfacial position at different times for as/ai = 1/1, 1/2 and 1/4. The other parameters are set as: Ax = At = 0.005, as = 0.005, Cs = 1, Cs/Ci = 1, St = 0.2. Numerical simulations are performed with 200 lattices in D1Q2, with temperature distributions at different times (t = 0.5, 1.0 and 1.5) shown in Fig. 3. The numerical solutions and analytical solutions have good coherence by comparison.
Because the three solid diffusivities as are all equal to 0.005, the temperature diffusions of solid phase (T G [—1, 0]) (Fig. 3) and the liquid-solid interfacial positions (Fig. 4) exhibit small changes. The temperature of liquid phase (T G [0,0.3]) is changed significantly in Fig. 3. The bigger the ai, the more rapidly the temperature diffuses. The longer the evolution time, the lower the inside temperature is. These conclusions are in good accordance with the actual situations. Furthermore, the phase change temperature Tm = 0 can be definitely found in Fig. 3.
Model 2 is solidification in a corner. This model19 is a two-dimensional problem. The model does not have top and right boundaries as a corner. After initialization, the temperatures on the left and bottom boundaries are both equal to T0. The values of the other temperatures are the same as used in Fig. 5.
To investigate the effect of a, the calculation is firstly carried out at as1 = 0.001 25, as2 = 0.002 5, as3 = 0.003 75, as4 = 0.005. The numerical values of St = 4.0, as/ai = 0.25, Cs = 1, Cs/Ci = 1, and Ax = At are set for computations, as shown in Fig. 6 and 400 x 400 lattices in D2Q9 is used.
From observation, the figures agree with the analytical solutions. The bigger the as, the faster the low
(a)The temperature distribution for as/al = 1/1, al = 0.005.
(b)The temperature distribution for as/al = 1/2, al = 0.01.
(c)The temperature distribution for as/al = 1/4, al = 0.02.
Fig. 3. The temperature distribution
Fig. 4. The panel showing the comparison of the liquidsolid interfacial positions between present computational solutions and analytical solutions for as/ai = 1/1, 1/2, 1/4.
Too =0.3
Liquid
Tm=0 I T=Ts(x, y, t)
M (x, y, t)
Interface
T=Ts(x, y, t)
x T0 = -1 x
Fig. 5. Schematic of solidification in a half -space.
temperature diffuses, and the more quickly the temperature inside the system drops down. There are 4 unambiguous lines of phase change temperature Tm shown in Fig. 6, where Tm are all equal to zero.
The solidification models manifest two predictable trends. The first one is that the internal temperature gets lower as time goes on. The other is that bigger a makes temperature diffuse more rapidly. These two points agree with the actual situations. The conclusion of these models manifests that the migration-collision scheme of LBM is not only available, but also more rigorous than Jiaung's.14 It can offer explicit liquid-solid interface and phase-change temperature Tm. Migrationcollision scheme is simple but powerful in dealing with the phase change problem of pure substance. The natural convection for phase change will be studied in the next work. The purpose of the present research is to offer scientific judgement for natural gas hydrate exploitation.
This work was supported by the National Natural Science Foundation of China (10932010,11072220 ) and the Natural Science Foundation of Zhejiang Province (Y607425, Z6090556 )
1. L. Jia and Z. Fang, Advanced Heat Transfer, (Higher Education Press, Beijing, 2008).
2. M. N. Ozisik, Heat Conduction, (Wiley-Interscience, 1993).
3. T. J. Lu, J. H. Du, S. Y. Lei, and B. X. Wang, Heat Mass Transfer. 37, 237 (2001).
4. V. R. Voller, C. R. Swaminathan, and B. G. Thomas, Int. J. Numer. Meth. Eng. 30, 875 (1990).
5. Y. L. He, Y. Wang, and Q. Li, The Theory and Applications of Lattice Boltzmann, (Science Press, Beijing, 2008).
6. Z. L. Guo and C. G. Zhen, The Principles and Applications of Lattice Boltzmann, (Science Press, Beijing, 2009).
7. X. W. Shan and H. D. Chen, Phys. Rev. E 47, 1815 (1993).
8. W. W. Yan, Y. Liu, Z. L. Guo and Y. S. Xu, Int. J. Mod. Phys. C 17, 771 (2006).
9. S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond, (Oxford University Press, USA, 2001).
10. D. Wolf-Gladrow, Lattice-Gas Cellular Automata and Lattice Boltzmann Models: An Introduction, (Springer-Verlag, 2000).
(a)Diagram showing temperature distribution (b)Diagram showing temperature distribution for aa = 0.00125. for a.s = 0.0025.
(c)Diagram showing temperature distribution (d)Diagram showing temperature distribution for as = 0.00375. for as = 0.005.
Fig. 6. Diagrams showing temperature distribution
11. D. Wolf-Gladrow, J. Stat. Phys. 79, 1023 (1995).
12. G. De Fabritiis, A. Mancini, D. Mansutti, and S. Succi, Int. J. Modern Phys. C 9, 1405 (1998).
13. W. Miller, S. Succi, and D. Mansutti, Phys. Rev. Lett. 86, 3578 (2001).
14. W. S. Jiaung, J. R. Ho, and C. P. Kuo, Numer. Heat Transfer B 39, 167 (2001).
15. D. Chatterjee and S. Chakraborty, Phys. Lett. A 86 320 (2005).
16. D. Chatterjee and S. Chakraborty, Int. J. Therm. Sci. 47, 552 (2008).
17. C. Huber, A. Parmigiani, B. Chopard, M. Manga, and O. Bachmann, Int. J. Heat Fluid Flow 29, 1469 (2008).
18. X. Shan and X. Y. He, Phys. Rev. Lett. 80, 65 (1998).
19. J. R. King, D. S. Riley, and A. M. Wallman, Proc. R. Soc. Lond. A 445, 3449 (1999).
