# Lattice Boltzmann kinetic modeling and simulation of thermal liquid-vapor system

###### Abstract

We present a highly efficient lattice Boltzmann (LB) kinetic model for thermal liquid-vapor system. Three key components are as beow: (i) a discrete velocity model by Kataoka et al. [Phys. Rev. E 69, 035701(R)(2004)]; (ii) a forcing term aiming to describe the interfacial stress and recover the van der Waals equation of state by Gonnella et al. [Phys. Rev. E 76, 036703 (2007)]; and (iii) a Windowed Fast Fourier Transform (WFFT) scheme and its inverse by our group [Phys. Rev. E 84, 046715 (2011)] for solving the spatial derivatives, together with a second-order Runge-Kutta (RK) finite difference scheme for solving the temporal derivative in the LB equation. The model is verified and validated by well-known benchmark tests. The results recovered from the present model are well consistent with previous ones[Phys. Rev. E 84, 046715 (2011)] or theoretical analysis. The usage of less discrete velocities, high-order RK algorithm and WFFT scheme with 16th-order in precision makes the model more efficient by about times and more accurate than the original one.

Keywords: lattice Boltzmann; liquid-vapor system; computational efficiency

PACS Nos.: 47.11.-j, 47.20.Hw, 05.70.Ln

## 1 Introduction

In the recent two decades the lattice Boltzmann (LB) method has been becoming a promising simulation tool for various complex systems, ranging from the low Mach number nearly incompressible flows to high speed compressible flows under shocks, from simple fluids to multiphase and/or multi-component fluids, from phase transition kinetics to hydrodynamic instabilities, etc. The LB methods in literature can be roughly classified into two categories, new solvers of various partial differential equations and kinetic models of various complex systems. In a recent mini-review it was pointed out that the LB kinetic model can be used to investigate the macroscopic behavior of the system due to its deviating from local thermodynamic equilibrium, including cases where the Navier-Stokes models have already be there. This idea was further specified in following works and considerable new physical insights for various systems were obtained from a more fundamental level.

Most early LB models for multiphase flows were for isothermal system and have been successfully applied to a wide variety of flow problems, such as drop collisions, wetting, phase separation and phase ordering, etc. In recent years, some thermal LB multiphase models have appeared. Those models can be roughly classified into three categories, i.e., the hybrid approach, the double-distribution-function approach, and the multispeed-extra-force (MEF) approach. In this paper we focus on the MEF approach which has a rigorous theoretical background. In the MEF approach, the one was developed by Gonnella, Lamura and Sofonea (GLS) is typical since it considers all physical factors, e.g., viscous dissipation, interfacial tension, and compression work done by the pressure, in both the momentum and energy equations.

A GLS-like model is formulated in two steps: firstly compose a thermal LB (TLB) model for ideal gas system, then an extra force term is added into the LB equation, where is responsible for the interfacial stress and recovering the nonideal equation of state (EOS). The GLS model was developed from the two-dimension 33-discrete-velocity (D2V33) TLB model by Watari and Tsutahara, and utilized to reproduce, in the continuum limit, the full set of thermohydrodynamic equations with the stress terms developed by Onuki. In a recent work, to make better the energy conservation and damp the spurious velocities to be negligible small in practical simulations, the Windowed Fast Fourier Transform (WFFT) and its inverse were introduced to calculate the spatial derivatives in the LB equation. For convenience of description, we refer to this model as the D2V33-FFT-TLB model.

The D2V33-FFT-TLB model is still subject to at least the following two constraints that greatly hamper its wider applications: (i) limited density ratio and temperature range; (ii) low computational efficiency. In the present study, we address mainly the later restriction from both the physical and computational point of views. In fact, computational efficiency is of essential importance, especially for models proposed to mimic the multiphase system. In a multiphase system without “fast factors” like shear, shock or strong convection, the evolution of the system depends mainly on some “slow mechanisms”, such as interface tension, viscous force, thermal diffusion, etc. Generally, long lasting simulations are needed to establish the growth properties.

## 2 The model

Here we present a highly efficient TLB model for van der Waals (VDW) fluids, including the following three parts: (i) two-dimension 16-discrete-velocity (D2V16) LB model for ideal gas system, (ii) an appropriate interparticle force accounting for the nonideal gas effects, and (iii) a higher-order windowed FFT (WFFT) scheme and a 2nd Runge-Kutta finite difference scheme. For convenience of description, we refer to this model as the D2V16-FFT-TLB model.

### 2.1 D2V16 TLB model for ideal gas system

Following the multispeed approach, Kataoka and Tsutahara proposed a TLB model for compressible flows. The model utilizes a discrete-velocity model (DVM) that has only 16 discrete velocities. The discrete equilibrium distribution function (DEDF) is calculated by a polynomial of the flow velocity up to the fourth order. Hydrodynamic quantities, such as density , velocity , and temperature are calculated by and , where is the discrete distribution function, represents the total degrees of freedom, is a free parameter introduced to describe the extra degrees of freedom corresponding to molecular rotation and/or vibration. Specifically, for and for .

To the best of our knowledge, the D2V16 DVM is the one with the least number of discrete velocities that can correctly recover the hydrodynamic equations at the NS level. Physically, it is much faster than other TLB models.

### 2.2 The forcing term for multiphase flow system

Although originally works only for ideal gas system, the multispeed model can be extended to multiphase system using the extra force method. For instance, GLS improved the D2V33 LB model by introducing an appropriate force term, , to describe the nonideal gas effects. Then the improved equation reads,

(1) |

where takes the following form:

(2) |

with , , and .

In the continuum limit Eq. (1) can recover the following NSEs for VDW fluids:

(3) | |||||

where is the reversible part of stress, comprising the VDW EOS and the nonideal gas interaction term , with , is the surface tension coefficient and a constant, the unit tensor. is the usual viscous stress tensor with the shear and bulk viscosities and . is the total energy density.

By applying Chapman-Enskog procedure, it is interesting to find that the incorporation of the force term in the LB equation makes no additional constrains on the requirements or the formulation process of the DEDF. In other words, if the DEDF satisfies the seven kinetic moments, i.e., Eqs. (1)-(7) as listed in Ref. [16], then the “original LB model + extra forcing term” is suitable for describing the multiphase flow system, without considering the specific formulation of the DEDF and the configuration of DVM. Consequently, the forcing term can be directly applied to other LB models or DVMs under the same framework, e.g., the aforementioned D2V16 model or the D2V61 model so as to reduce the computational load or improve the isotropy. Here we choose the D2V16 TLB model to replace the D2V33 TLB model so as to increase the computational efficiency.

### 2.3 Numerical schemes for spatial and temporal derivatives

From Ref. [8], it is well known that decides the accuracy of the FFT scheme. Here we further develop the WFFT scheme through applying the with 16th order in accuracy, i.e., takes the following form:

(4) |

where , and can be obtained in a similar way. Some simple derivations demonstrate that the FFT approach with has a 16th-order accuracy in space. The time evolution of is derived by the 2nd finite difference scheme and , where the superscripts , indicate the consecutive two iteration steps.

## 3 Simulation results and physical analysis

In this section, several numerical investigations are conducted to validate the physical properties of the D2V16-FFT-TLB model.

### 3.1 Liquid-vapor coexistence curve

To check if the D2V16-FFT-TLB model can reproduce the correct equilibrium thermodynamics of the multiphase system, simulations about the liquid-vapor coexistence curve at various temperatures are performed. Simulations are carried out over a domain with Periodic Boundary Conditions (PBCs) in both the - and -directions. The initial conditions are set as if , if ; and if , where , and . Parameters are , , , , , , , . Here the initial temperature is . It will be dropped by a small value when the equilibrium state has been achieved. Figure 1(a) shows the liquid-vapor coexistence curve obtained from the LB simulations at various temperatures, compared to the theoretical predictions from Maxwell construction. The two sets of results have a satisfying agreement. Moreover, the lowest temperature that the D2V16-FFT-TLB model can simulate is about , while for the D2V33-FFT-TLB model, it is about . So the former is stuitable for a wider temperature range and a wider range of density ratio.

### 3.2 Phase separation after quench

Here the thermal phase separation after quench is employed to validate if the new model can recover the accurate physical scenario. Simulations are performed on a lattice nodes. The initial conditions are set as , where is a random density noise with an amplitude of . Parameters are set to be , for D2V33-FFT-TLB model and for D2V16-FFT-TLB model, , . The FTT scheme with 8th order in precision and the 1st Euler scheme are used for the D2V33-FFT-TLB model, while the FFT scheme with 16th order in precision and the 2nd RK scheme are for the D2V16-FFT-TLB model.

Typical snapshots of configurations are plotted in Fig. 2. After quench the fluid begins to separate spontaneously into small regions with higher and lower densities. The density difference increases with time. Under the actions of surface tension and viscous force, small domains coalesce to each other and larger ones appear at about . From patterns at and , the liquid and vapor domains evolve in an equal way, leading to an interwoven bicontinuous pattern. The growth of domains continues at , and eventually at , the system reaches a completely separated state with least surface.

Compared to the D2V33-FFT-TLB model, the D2V16-FFT-TLB model recovers nearly the same results, particularly at former times. With the accumulation of physical and numerical differences between the two models, the patterns present larger differences. For example, at , the vapor phase at the right boundary calculated from the D2V16-FFT-TLB model is much larger than that calculated from the D2V33-FFT-TLB model. And at , the liquid band simulated by the D2V16-FFT-TLB model is located lower than the one mimicked by the later. However, if we focus on the evolution mechanism of the system, rather than the morphologies at specific times, the differences can be neglected. Furthermore, physically, due to the utilization of the DVM with the least discrete velocities and, numerically, due to application of higher-order RK scheme for temporal derivative, the D2V16-FFT-TLB model is about 10 times faster than the D2V33-FFT-TLB model. This is of essential importance for studying a phase-separating system, because the enhancement in computational efficiency means the decrease in simulation time we needed, the increase in lattice nodes we employed and the feasibility of the simulation on a personal computer.

Figure 1(b) depicts variations of the total energy for the process shown in Fig. 2. To be seen is that the D2V16-FFT-TLB model not only owns higher efficiency but also higher accuracy and better energy conservation in simulations. When the new model is used, remarkably reduces from to .

### 3.3 Large deformation of a rectangular droplet

To further assess if the new model is appropriate for simulating large deformation problem and examine its ability for refraining spurious velocity, a large deformation problem under surface tension is studied. Initially a rectangular droplet with length and height is placed at the center of the whole domain, where the lattice unit is used. The liquid and vapor densities are and , respectively. Throughout the simulation, the temperature is fixed to be . Parameter are , , and . For both models, the FFT scheme with 16 order in accuracy and the 2nd RK are employed. Figure 3 illustrates snapshots of the droplet shapes at various times. When the simulation starts, the droplet will firstly transform its configuration at the left and right sides, where the surface tension reaches its maximum due to the largest curvature. As time evolves, the droplet gradually shrinks to a dumbbell shape for minimizing the net interfacial energy. At , it becomes approximately an ellipse. After that, it oscillates through changing the length of the major axis and minor axis. Due to the effects of surface tension and viscosity, the deformed droplet relaxes to a circular equilibrium shape at about . The simulation results obtained from the two models are consistent with each other, and are qualitatively agreement with Lamb’s theoretical analysis and other researchers’ results.

Evolution of the maximum spurious velocity during the large deformation procedure is illustrated in Fig. 1(c). When , the system approximately achieve its equilibrium and almost remains unchanged. With the identical numerical algorithms, the new model has the similar ability to refrain .

## 4 Conclusions

In the paper, a highly efficient LB model for thermal liquid-vapor system is proposed. Due to the applications of DVM with less discrete velocities and higher-order schemes, the new model is more efficient, has more potential to keep conservation of the total energy and refrain spurious velocity near the liquid-vapor interface. Therefore, it is more suitable for studying thermohydrodynamics of multiphase flows. In future studies, by inserting more realistic EOS, we will increase the density ratio that the new model can simualte, and will investigate the hydrodynamic instabilities which are highly interested in science and engineering and where various nonequilibrium effects are pronounced.

## Acknowledgements

This work is partly supported by National Natural Science Foundation of China (under Grant Nos. 11071024, 11202003 and 11203001), Foundation of State Key Laboratory of Explosion Science and Technology (under Grant No. KFJJ14-1M), Natural Science Foundation of Hebei Province (under Grant No. A2013409003), Foundation for Outstanding Young Scholars of Hebei Province, Excellent Youth Foundation of Hebei Educational committee (under Grant No. YQ2013013) and Technology Support Program of LangFang (under Grant Nos. 2011011023, 2012011021/26/30).

## References

- 1. S. Succi, The Lattice Boltzmann Equation for Fluid Dynamics and Beyond (Oxford University Press, New York, 2001).
- 2. A. Xu, G. Zhang, Y. Li and H. Li, Prog. Phys. (2014, in press)
- 3. J. M. V. A. Koelman, EPL 15, 603 (1991).
- 4. A. Xu, G. Zhang, Y. Gan, F. Chen and X. Yu, Front. Phys. 7, 582 (2012).
- 5. Y. Gan, A. Xu, G. Zhang and Y. Yang, EPL 103, 24003 (2013).
- 6. http://www-thphys.physics.ox.ac.uk/people/JuliaYeomans/.
- 7. V. Sofonea, A. Lamura, G. Gonnella and A. Cristea, Phys. Rev. E 70, 046702 (2004).
- 8. A. Xu, G. Gonnella and A. Lamura, Phys. Rev. E 67, 056105 (2003); Phys. Rev. E 74, 011505 (2006); Physica A 331, 10 (2004).
- 9. G. Gonnella, A. Lamura and V. Sofonea, Phys. Rev. E 76, 036703 (2007).
- 10. Y. Gan, A. Xu, G. Zhang, Y. Li and H. Li Phys. Rev. E 84, 046715 (2011); Y. Gan, A. Xu, G. Zhang, Y. Li Commun. Theor. Phys. 57, 681 (2012); Front. Phys. 7, 481 (2012).
- 11. Y. Gan, A. Xu, G. Zhang and Y. Li, Phys. Rev. E 83, 056704 (2011).
- 12. B. Yan, A. Xu, G. Zhang, Y. Ying, H. Li, Front. Phys. 8, 94 (2013); C. Lin, A. Xu, G. Zhang G, Y. Li and S. Succi, Phys. Rev. E 89, 013307 (2014); F. Chen, A. Xu, G. Zhang, Y. Wang, Front. Phys. 9, 246 (2014).
- 13. K. N. Premnath and J. Abraham, Phys. Rev. E 71, 056706 (2005).
- 14. J. Hyväuoma, P. Raiskinmäi, A. J äberg, A. Koponen, M. Kataja, and J. Timonen, Phys. Rev. E 73, (2006) 036705.
- 15. R. Zhang and H. Chen, Phys. Rev. E 67, 066711 (2003).
- 16. G. Gonnella, A. Lamura, A. Piscitelli and A. Tiribocchi, Phys. Rev. E 82, 046302 (2010).
- 17. A. Márkus and G. Házi, Phys. Rev. E 83, 046705 (2011).
- 18. M. Watari and M. Tsutahara, Phys. Rev. E 67, 036306 (2003).
- 19. A. Onuki, Phys. Rev. Lett. 94, 054501 (2005) ; Phys. Rev. E 75, 036304 (2007).
- 20. T. Kataoka and M. Tsutahara, Phys. Rev. E 69, 035701(R) (2004).
- 21. A. Xu, EPL 69, 214 (2005).
- 22. W. H. Press, S. A. Teukolsky, W. T. Vetterling and B. P. Flannery Numerical Recipes (Cambridge University Press, Cambridge, 2007).
- 23. H. Lamb, Hydrodynamics (Dover, New York, 1932).
- 24. X. Yang, J. J. Feng, C. Liu and J. Shen, J. Comput. Phys. 218, 417 (2006).