Casting and Solidification

Acceleration of Macrosegregation Simulation Based on Lattice Boltzmann Method

2018 Volume 58 Issue 1 Pages 114-122

Details

Abstract

Lattice Boltzmann Method (LBM), newly developing technique of computational fluid dynamics, is coupled with solute and energy conservation equations to develop a macrosegregation simulation model (LBM-coupled model) with high computational efficiency. LBM does not require time-consuming calculations for correction of velocity and pressure of fluid in contrast to methods directly solving Navier-Stokes (NS) equation and, therefore, LBM is computationally efficient. In this study, the accuracy of the LBM-coupled model is investigated by calculating the steady state flows and by comparing the results with those of analytical solutions and a conventional model based on the NS equation. The results between them are almost identical with each other and it indicates that the accuracy of the LBM-coupled model is sufficiently high. Moreover, a macrosegregation simulation is carried out for a simple case where the macrosegregation emerges only by natural convection, by means of the LBM-coupled and conventional models. The LBM-coupled model yields almost the same result with the one of NS-based model. Importantly, however, the simulation of LBM-coupled model is about five times faster than the one of NS-based model.

1. Introduction

Prediction and elimination of macrosegregation have been important issues in continuous casting and ingot casting processes.^{1,2,3,4)} Macrosegregation is non-uniform distribution of alloying elements over spatial scales ranging to ingot scales. It is caused by long-range advections of alloying elements due to flows of segregated liquid and motion of solid. It is not straightforward to understand and control macrosegregation behavior in detail because several different mechanisms are involved in occurrence of the advection of alloying elements such as thermosolutal convection, a flow due to solidification shrinkage, forced flows due to pouring, rotation, bending of cast and so on. Computational approaches have played an important role in this field as well as experimental and theoretical approaches.^{1,2,3,4)}

Several models have been developed for numerical simulations of formation processes of macrosegregation, for instance, so-called single domain models based on a mixture theory^{5)} or volume averaging,^{6,7,8)} two-phase model which accounts for relative movement of solid and liquid, and multi-phase models^{9,10,11)} where mesoscopic details such as the interdendritic liquid, the columnar and equiaxed dendrites are separately modelled. The simulation models have been gradually sophisticated and have been widely applied to a variety of complicated and realistic processes. However, there is a longstanding problem in practical use of these models, *i.e.*, high computational cost.^{4)} It restricts the computational system to a small size and/or it forces one to use coarse mesh size, often hampering accurate description of macrosegregation behavior.^{4)} Hence, it is important to develop a simulation model with high computational efficiency. The simulation models for macrosegregation consist of a set of momentum, mass, energy and solute conservation equations^{2)} and the momentum conservation is described by the Navier-Stokes (hereafter abbreviated as NS) equation for fluid flows. Although several numerical methods were developed for solving the NS equation such as MAC^{12)} and SIMPLE methods,^{13)} a time-consuming calculation is required for correction of the velocity and pressure in these methods. The computational burden for solving NS equation is generally high and it occupies a large part of the computational cost of macrosegregation simulations. Hence, a key to acceleration of the macrosegregation simulations lies in the calculation of fluid flow.

Lattice Boltzmann Method (LBM) has been attracting a great deal of attention as a promising method for efficiently computing the fluid dynamics.^{14,15,16,17,18,19,20)} It is a newly developing technique of computational fluid dynamics. LBM describes the time evolution of particle distribution function, from which one can calculate the macroscopic quantities such as the density, velocity and pressure of fluid. Importantly, an asymptotic analysis called Chapman-Enskog analysis shows that LBM recovers the NS equation at macroscopic scale by appropriately modeling the equilibrium distribution function. One of the advantages of LBM is that it does not require the time-consuming calculation for correction of velocity and pressure and, therefore, LBM is in general computationally more efficient than the methods directly solving the NS equation. In addition, LBM is suitable for parallel computing, which is a notable advantage in the light of recent rapid development of parallel computing technique. LBM has been increasingly utilized for simulations in a lot of important systems such as hydrodynamic systems, magnetohydrodynamic systems and also multiphase and multicomponent fluids.^{15,16,18)} Furthermore, it was recently coupled with a phase-field method to simulate solidification microstructures under fluid flows.^{21)}

The main objective of this study is to develop a macrosegregation model coupled with LBM with a view to accelerating the macrosegregation simulation. In this study, the LBM-coupled model is constructed by combining LBM for fluid flow with energy and solute conservation equations. We investigate the accuracy of the LBM-coupled model by comparing the results of the steady state flows with the results of analytical solutions and a conventional model based on the NS equation. Moreover, we carry out the macrosegregation simulation for a simple case, where the macrosegregation emerges only by natural convection, by means of the LBM-coupled and conventional models to compare their computational speeds. This paper is organized as follows. The conventional macrosegregation model based on the NS equation is first explained in the next section. LBM and the LBM-coupled model are described in section 3, followed by the results and discussion in section 4. The conclusions are given in the last section.

2. Conventional Simulation Model Based on the Navier-Stokes Equation

The focus of this study is acceleration of the macrosegregation simulation. It is tackled by coupling LBM with a solidification model, more precisely, energy and solute conservation equations. The present modeling rests on a solidification model in a NS-based macrosegregation model developed by Sawada *et al.*^{22,23)} which was chosen in the light of usability and ease of numerical implementation in this study. This model is hereafter called the NS-based model and it is briefly explained below.

In the NS-based model, the incompressible fluid flow is calculated by the NS equation within the Boussinesq approximation with addition of the Darcy’s flow term, which is given by,

(1) |

(2) |

(3) |

(4) |

Several numerical algorithms can be employed for solving this model especially Eqs. (1) and (2). For instance, the simulation of channel segregation in a Sn–Bi alloy^{23)} was carried out based on the SOLA (HSMAC) method with a semi-implicit algorithm developed for stabilization of calculation of Darcy’ term.^{22)} Also, the SIMPLE method was employed in the simulation of center-line segregation in steel.^{24)} Note that the time-consuming calculation for correction of **u** and *p* is necessary in these methods. The acceleration of the simulation requires improvement in calculation efficiency for **u**, which is tackled by coupling LBM with this model as detailed in the next section.

In this study, the numerical simulation of above-mentioned NS-based model was carried out to investigate the accuracy of the LBM-coupled model. All calculations in this study were conducted in two-dimensional systems. According to the previous work,^{22)} the fluid flow was calculated based on the SOLA method with the semi-implicit algorithm which improves the numerical stability associated with the Darcy’s term. Equations (3) and (4) were discretized based on second-order finite difference formulas with a square grid spacing of *δx*. The time evolutions of *T* and *c _{l}* were solved using a simple first-order Euler scheme. The temperature recovery method was employed for Eq. (3). The time change of

Note that the accuracy and the computational burden are determined by the criterion of convergence in the iterative calculation for correction of velocity and pressure in the SOLA method. The criterion is given by a parameter *α* as *α*<∇·**u** and it was set to *α* = 1.0 × 10^{−6} in all the simulations of the NS-based model.

3. Lattice Boltzmann Method (LBM)-coupled Model

LBM can be viewed as a discrete version of the Boltzmann equation. It is a direct extension of the lattice gas cellular automaton (Boolean type model), a special kind of cellular automaton for fluid simulation.^{15)} LBM is a relatively new technique of computational fluid dynamics and it has attracted considerable attention in a variety of fields. Moreover, LBM has been extended to describe a reaction diffusion equation with advection term such as Eqs. (3) and (4).^{26)} Therefore, the NS-based model explained in section 2 can be fully recast into the LBM framework. However, our particular attention in this study is directed at the acceleration of the computation of fluid flow. Therefore, LBM for incompressible fluid flow is coupled with the solidification model described by Eqs. (3) and (4) in this study.

In LBM, the fluid consists of fictive and microscopic particles moving at specified (discrete) velocities in specified (discrete) directions on a lattice and their collective behavior determines the macroscopic quantities of fluids such as the density, velocity and pressure. More specifically, a particle distribution function for the discrete velocity vector is calculated by the following lattice Boltzmann equation,

(5) |

(6) |

(7) |

(8) |

(9) |

(10) |

Accordingly, *τ* is given by

(11) |

The above-mentioned model is one of the standard LBMs employed for analyses on fluid dynamics. To construct the LBM-coupled model for macrosegregation, influences of solid fraction, Darcy’s flow and thermosolutal convection must be taken into account. For this, we utilize a model proposed for isothermal incompressible flow in porous media^{27)} and slightly modify this model to make it suitable for coupling with the solidification model in section 2. In the present LBM, the force term is included in the lattice Boltzmann equation and Eq. (5) is accordingly redefined as follows,

(12) |

(13) |

(14) |

(15) |

(16) |

(17) |

Note that it is not necessary to carry out the time-consuming calculation for correction of **u** and *p* in solving LBM, in contrast to the conventional numerical schemes for the NS equation. The lattice Boltzmann Eq. (12) exhibits the discretized form and therefore it is directly solved by the numerical calculation in a simple explicit manner. In this regards, the numerical implementation of this model is basically not complicated. Furthermore, although it is not the focus of this study, this model is suitable for parallel computing.

Let us explain some points about numerical simulations of the LBM-coupled model. First as mentioned above, the present model is based on D2Q9 model in two-dimensional square lattice with the lattice spacing *δx*. In the numerical calculation of the LBM-coupled model, all the quantities, *i.e.*, *f _{i}* and thereby

(18) |

(19) |

4. Results and Discussion

4.1. Accuracy of LBM-coupled Model for Steady State Flow
We investigate the accuracy of the LBM-coupled model, focusing on the contributions to the steady state flow one by one. We first focus on the steady state flow between parallel planes (Poiseuille flow) in a fluid system where the calculation accuracy of viscous flow can be examined. Then, we calculate the steady state flow between parallel planes in a solid and liquid two-phase system (equivalent to porous media) where the Darcy’s term plays an important role. In addition, we examine natural convention in a fluid system driven by the temperature difference, where the contribution of buoyancy term becomes important. These are discussed in this subsection. The accuracy and computational efficiency for the macrosegregation simulations are described in the next subsection. As described in section 2, all the simulations of NS-based model were carried out based on the SOLA method.

The analytical solution of the NS equation can be obtained for the steady state flow between parallel planes in a fluid system. The fluid velocity in *x*-direction *u _{x}* between two planes parallel in

(20) |

The computational system for this analysis is schematically shown in **Fig. 1**. It consists of 128 and 16 lattice points in *x*- and *y*-directions, respectively, which is sandwiched between two planes parallel in *x*-direction. A non-slip boundary condition was applied to both planes (*i.e.*, at *y* = 0 and *y* = 16*δx*) based on the bounce-back boundary scheme.^{28)} For the sake of convenience, all quantities are normalized by *c _{s}* (or

Fig. 1.

Schematic illustration of two-dimensional system for steady-state flow in a fluid between two parallel planes.

Fig. 2.

Calculated results of steady-state velocity in a fluid between two parallel planes. The results of LBM, NS and Eq. (20) are represented by filled circle, open triangle and solid line, respectively.

Next, we focus on the steady state flow between parallel planes in a solid and liquid two-phase system. In this case, the Darcy’s term is considered in the NS equation. The analytical solution of this problem is given as^{22)}

(21) |

Fig. 3.

Calculated results of steady-state velocity in a fluid and solid two-phase (*f _{s}* = 0.5) between two parallel planes. The results of LBM, NS and Eq. (21) are represented by filled circle, open triangle and solid line, respectively.

Also, we investigate the steady state flow in a fluid system under a temperature gradient where the accuracy associated with the buoyancy term can be discussed. *δρ* in Eq. (14) (also in Eq. (1)) was given by *δρ*/*ρ* = 1 −*β _{T}*(

Fig. 4.

Schematic illustration of two-dimensional system for steady-state flow under temperature gradient.

Fig. 5.

Results of temperature distribution and the fluid velocity during steady-state in a fluid system under the temperature gradient, calculated by (a) the LBM-coupled model and (b) NS-based model. The temperature distribution is represented by the gray scale and the contour line at 0.002 *T*/*T*_{0} interval. The velocity vector is indicated by the arrow, the length of which represents the magnitude of the velocity.

In this subsection, we investigated the accuracy of the LBM-coupled model, focusing on the contributions to the steady state flow one by one. This investigation shows that the accuracy of the LBM-coupled model is sufficiently high. Since the main concern of this subsection was the accuracy of the present LBM, the computational speed was not discussed. Actually, the calculation of LBM-coupled model was always faster than the one of NS-based model in all three cases discussed above. However, the calculation times for the steady state flow were quite short and they are not suitable for accurate assessment of the computational efficiency of the present LBM. In the next subsection, we focus on the macrosegregation simulations and discuss the computational speed as well as the accuracy.

4.2. Macrosegregation SimulationWe have carried out the macrosegregation simulations using the LBM-coupled and NS-based models to assess the computational speed and the accuracy of the former model. In these calculations, the permeability was given by the Kozeny-Carman model as described in the previous section. The density change in the buoyancy term was given as (*ρ* + *δρ*)/*ρ* ≈ 1−*β _{c}*(

For the sake of convenience, we focus on a model alloy system consisting of A and B atoms. The parameters employed in the simulation are shown in **Table 1**. Typical values in metallic alloy systems were assigned to most of parameters and some values such as Δ*H*, *β _{T}* and

Table 1. Parameters employed for macrosegregation simulations.

Parameters | Value |
---|---|

Lattice spacing, δx [m] | 1.0×10^{−4} |

Time step, δt [s] | 1.0×10^{−4} |

Reference density, ρ [kg/m^{3}] | 8200 |

Alloy composition, c_{0} [mass fraction] | 0.1 |

Liquidus temperature, T_{0} [K] | 480 |

Latent heat, ΔH [J/kg] | 60000 |

Heat capacity [J/(kg K)], C_{p} | 1000 |

Equilibrium partition coefficient, k [–] | 0.2 |

Kinematic viscosity, v [m^{2}/s] | 5.0×10^{−6} |

Thermal diffusivity, a [m_{T}^{2}/s] | 1.2×10^{−6} |

Coefficient of volume expansion due to temperature difference, β [1/K]_{T} | 1.0×10^{−4} |

Coefficient of volume expansion due to concentration difference, β [–]_{c} | −0.244 |

Dendrite arm spacing, λ [m]_{s} | 1.7×10^{−4} |

Heat transfer coefficient, h [W/(m^{2} K)] | 10000 |

The result of the LBM-coupled model is shown in **Fig. 6**. These snapshots represent the distributions of temperature, solid fraction and fluid velocity at *t* = 150, 160 and 170 s. The temperature is represented by the gray scale and the solid fraction is shown by the contour line at intervals of 0.01. The fluid velocity is indicated by the arrow, the length of which represents its magnitude. Since the partition coefficient is less than 1 in this alloy, the solute atom (B atom) is gradually enriched in the liquid as the solidification proceeds. Then, the downward flow of B-rich liquid occurs, because the density of this alloy increases with increase in the concentration of B. The anti-clockwise flow is observed. It washes away the B-rich liquid from the upper region and the solidification in upper region accordingly proceeds faster than the lower region. **Figure 7** shows the result of NS-based model. Although very slight differences appear between Figs. 6 and 7, the result of LBM-coupled model agrees well with the one of NS-based model.

Fig. 6.

Result of macrosegregation simulation by the LBM-coupled model. The snapshots (a), (b) and (c) correspond to the results at *t* =150, 160 and 170 s, respectively. In each figure, the temperature is represented by gray scale and the solid fraction is shown by the contour line at intervals of Δ*f _{s}* = 0.01. The maximum and minimum values for the contour lines are indicated. The fluid velocity is shown by the arrow, the length of which indicates the magnitude of the velocity.

Fig. 7.

Result of macrosegregation simulation by the NS-based model. The snapshots (a), (b) and (c) correspond to the results at *t* = 150, 160 and 170 s, respectively. In each figure, the temperature is represented by gray scale and the solid fraction is shown by the contour line at intervals of Δ*f _{s}* = 0.01. The maximum and minimum values for the contour lines are indicated. The fluid velocity is shown by the arrow, the length of which indicates the magnitude of the velocity.

**Figure 8** shows the results of the distribution of segregation ratio in the end of each simulation. In this study, the segregation ratio is defined as the local concentration divided by the average concentration. The contour line for the segregation ratio is shown at intervals of 0.01. The results of the LBM-coupled and NS-based models yield almost identical results where the negative and the positive segregations form in the upper and lower regions, respectively. From a closer look at Fig. 8, one can find slight differences in contour lines between both models. The minimum and maximum values of segregation ratio are 0.969 and 1.071, respectively, in the LBM-coupled model, while those are 0.962 and 1.093 in the NS-based model. Note that the time evolutions of *T* and *c _{l}* were calculated in the same manner in both models. Hence, the slight differences in the segregation ratio in Fig. 8 originate only from the accuracy of numerical calculation of fluid flow in both approaches. In principle, complete agreement between both models cannot be expected because the numerical scheme of NS equation as well as LBM basically provides not a rigorous solution but an approximate solution, each of which suffers from different kinds of errors in numerical calculations. Although such small difference exists between these results, it can be concluded that the accuracy of the LBM-coupled model is as high as the one of NS-based model.

Fig. 8.

Distribution of segregation ratio in the end of simulations by (a) LBM-coupled model and (b) NS-based model. The contour line for segregation ratio is drawn at intervals of 0.01. The maximum and minimum values for the contour line are indicated.

The main purpose of the present study is the acceleration of macrosegregation simulation as described in the introduction. In this study, all the calculations were carried out with CPU of Xenon W5590 3.33 GHz. The computational times for the macrosegregation simulation shown in Figs. 6, 7, 8 were 48 and 245 min in the LBM-coupled and NS-based models, respectively. Hence, the LBM-coupled model is about five times faster than the NS-based model and the acceleration is successfully achieved. Note that the computational cost of the NS-based model mainly originates from the iterative calculation for correction of **u** and *p*. Hence, the computational cost of NS-based model can be lowered by relaxing the convergence criterion, *i.e.*, increasing the value of *α*. However, it causes low accuracy of the NS-based model. Also, in general, once the solidification starts in local region of an ingot, the computational burden for the iterative calculation in the NS-based model increases due to the contribution of the Darcy’s term. Hence, when both models are applied to the simulations of solidification for a longer time in larger ingot than the present case, the difference in the computational time between these models should be remarkably increased. In addition, LBM is suitable for parallel computing and it is basically straightforward to achieve good performance of parallel scaling.^{29)} Therefore, the advantage of the LBM-coupled model over the NS-based model becomes more salient in parallel computation.

Before closing this section, let us point out an important issue about numerical stability of LBM. As mentioned in section 3, there is a restriction in the Prandtl number in the LBM-coupled model as given by Eq. (19). Hence, a high value was assigned to *ν* in this study. This problem arises because the explicit scheme was utilized for solving Eq. (3) and it can be avoided using an implicit scheme. However, the implicit scheme is time-consuming in some cases. Moreover, the explicit scheme is rather preferable when the parallel computing is considered. Also, there is another problem in the present LBM-coupled model. Except for inclusion of Darcy’s term and buoyancy term in a two-phase system (porous media), LBM that the present modeling relies upon corresponds to the standard LBM based on the BGK approximation. It is well-known that the BGK-based model exhibits the numerical instability at high Reynolds number *Re*. This problem causes a limitation in the application range of LBM.^{30)} In fact, we employed small system size to avoid this problem in the present macrosegregation simulation. In this regard, new LBMs that are stable even at high *Re* have recently been developed such as the multi-relaxation time (MRT) collision model^{31)} and the modified BGK collision model.^{30)} In these models, moreover, the numerical simulation can be made stable even when *τ* is very close to 0.5, which enables the calculation for low value of *ν*. Therefore, when one of these models is utilized, it is possible to remove the limitations of the present LBM-coupled model (*i.e.*, high *v* and low *Re*). Although computational burdens of these models are generally higher than the one of the present LBM, they do not require the time-consuming calculation for correction of velocity and pressure. Therefore, one can expect that a LBM-coupled model based on one of these new LBMs should be computationally more efficient than the conventional method based on the NS equation. This improvement remains to be tackled in a future work.

5. Conclusions

Computational approaches offer effective way to analyze and predict the macrosegregation behavior in casting. One of the common drawbacks in simulation models for macrosegregation is high computational cost associated with the calculation of fluid flow. Instead of directly solving the Navier-Stokes (NS) equation, one can calculate the fluid flow by solving the simple finite difference equation for distribution function of fictive particles in the lattice Boltzmann method (LBM). In order to develop a macrosegregation simulation model with high computational efficiency, in this study, we coupled LBM including the Darcy’s flow and buoyancy flow in a solid-liquid two-phase system with the energy and solute conservation equations. The accuracy of this model (LBM-coupled model) was investigated in detail by performing the simulations for three types of steady state flows and by comparing the results with those of analytical solutions and the conventional model (NS-based model). The results between them are almost identical with each other, demonstrating that the accuracy of the LBM-coupled model is sufficiently high.

The LBM-coupled model was applied to a simulation of macrosegregation in ingot of a binary alloy. The result is consistent with the one of NS-based model and this fact also supports the accuracy of the LBM-coupled model. Importantly, the simulation of LBM-coupled model is about five times faster than the one of NS-based model and the acceleration of macrosegregation simulation was successfully achieved. This advantage of the LBM-coupled over NS-based models should become more significant when a more realistic case, *i.e.*, solidification for a longer time in a larger ingot is considered and when parallel computational technique is utilized. However, the numerical stability of LBM restricts the application range of the present LBM-coupled model. More specifically, this model cannot be applied to flows at high Reynolds number in a fluid with low kinematic viscosity. In this regard, the improvement of the model based on recently developed LBM remain as an important work.

Acknowledgements

This research was supported by 25th ISIJ research promotion grant. The authors wish to thank Dr. Tomoki Sawada for providing them with detailed information about the NS-based model. Also, the authors are thankful to Prof. Kiyotaka Matsuura for stimulating discussions and for providing valuable comments on the manuscript.

References

- 1) M. C. Flemings:
*ISIJ Int.*,**40**(2000), 833. - 2) C. Beckermann:
*Int. Mater. Rev.*,**47**(2002), 243. - 3) G. Lesout:
*Mater. Sci. Eng. A*,**413–414**(2005), 19. - 4) E. D. Pickering:
*ISIJ Int.*,**53**(2013), 935. - 5) W. D. Bennon and F. P. Incropera:
*Int. J. Heat Mass Transf.*,**30**(1987), 2161. - 6) C. Beckermann and R. Viskanta:
*Physicochem. Hydrodyn.*,**10**(1988), 195. - 7) C. Beckermann and R. Viskanta:
*Appl. Mech. Rev.*,**46**(1993), 1. - 8) M. C. Schneider and C. Beckermann:
*ISIJ Int.*,**35**(1995), 665. - 9) J. Ni and C. Beckermann:
*Metall. Trans. B*,**22B**(1991), 349. - 10) M. Wu and A. Ludwig:
*Metall. Mater. Trans. A*,**37A**(2006), 1613. - 11) M. Ahmadein, M. Wu and A. Ludwig:
*J. Cryst. Growth*,**417**(2015), 65. - 12) A. A. Amsden and F. H. Harlow:
*J. Comput. Phys.*,**6**(1970), 322. - 13) S. V. Patankar: Numerical Heat Transfer and Fluid Flow, Hemisphere, Washington, DC, (1980), 126.
- 14) B. Benzi, S. Succi and M. Vergassola:
*Phys. Rep.*,**222**(1992), 145. - 15) S. Chen and G. D. Doolen:
*Annu. Rev. Fluid Mech.*,**30**(1998), 329. - 16) C. K. Aidun and J. R. Clausen:
*Annu. Rev. Fluid Mech.*,**42**(2010), 439. - 17) X. He and L. S. Luo:
*J. Stat. Phys.*,**88**(1997), 927. - 18) X. He and L. S. Luo:
*Phys. Rev. E*,**56**(1997), 6811. - 19) T. Inamuro, M. Yoshino and F. Ogino:
*Phys. Fluids*,**9**(1997), 3535. - 20) Z. Guo, B. Shi and N. Wang:
*J. Comput. Phys.*,**165**(2000), 288. - 21) R. Rojas, T. Takaki and M. Ohno:
*J. Comput. Phys.*,**298**(2015), 29. - 22) T. Sawada, K. Oikawa, K. Anzai, F. Takahashi, K. Kajikawa and H. Yamada:
*J. JFS*,**81**(2009), 283. - 23) T. Sawada, K. Oikawa, K. Anzai, F. Takahashi, K. Kajikawa and H. Yamada:
*J. JFS*,**81**(2009), 177. - 24) T. Murao, T. Kajitani, H. Yamamura, K. Anzai, K. Oikawa and T. Sawada:
*ISIJ Int.*,**54**(2014), 359. - 25) T. Sawada:
*Jpn. Steel Works Tech. Rep.*,**64**(2013), 9. - 26) Z. Chai and T. S. Zhao:
*Phys. Rev. E*,**87**(2013), 063309. - 27) Z. Guo and T. S. Zhao:
*Phys. Rev. E*,**66**(2002), 036304. - 28) X. He, Q. Zou, L. S. Luo and M. Dembo:
*J. Stat. Phys.*,**87**(1997), 115. - 29) M. D. Mazzeo and P. V. Coveney:
*Comput. Phys. Commun.*,**178**(2008), 894. - 30) L. Wang, J. Mi and Z. Guo:
*Int. J. Heat Mass Transf.*,**94**(2016), 269. - 31) Q. Liu, Y.-L. He, Q. Li and W.-Q. Tao:
*Int. J. Heat Mass Transf.*,**73**(2014), 761.

© 2018 by The Iron and Steel Institute of Japan