Lattice Boltzmann simulation of liquid–vapor system by incorporating a surface tension term*
Song Bao-Wei, Ren Feng, Hu Hai-Bao, Huang Qiao-Gao
School of Marine Science and Technology, Northwestern Polytechnical University, Xi’an 710072, China

Corresponding author. E-mail: songbaowei@nwpu.edu.cn

Project supported by the National Nature Science Foundation of China (Grant No. 51109178) and the Specialized Research Fund for the Doctoral Program of Higher Education of China (Grant No. 20116102120009).

Abstract

In this study, we investigate the pseudopotential multiphase model of lattice Boltzmann method (LBM) and incorporate a surface tension term to implement the particle interaction force. By using the Carnahan–Starling (CS) equation of state (EOS) with a proper critical pressure–density ratio, a density ratio over 160000 is obtained with satisfactory numerical stability. The added surface tension term offers a flexible choice to adjust the surface tension strength. Numerical tests of the Laplace rule are conducted, proving that smaller spurious velocity and better numerical stability can be acquired as the surface tension becomes stronger. Moreover, by wall adhesion and heterogeneous cavitation tests, the surface tension term shows its practical application in dealing with problems in which the surface tension plays an important role.

Keyword: 47.61.Jd; 68.03.Cd; 47.55.D–; 47.50.–d; lattice Boltzmann method; surface tension; pseudopotential model; numerical stability
1. Introduction

Multiphase problems involving liquid– vapor systems have grown to be a hot topic in scientific research and engineering application fields. Due to its good performance in dealing with mesoscopic and multiphase problems, the lattice Boltzmann method (LBM) emerges as an attractive method.[13] Among the multiphase models, three well-accepted models are: the color-gradient model, [4] the phase-field-based model, [5] and the pseudopotential model.[1] Despite the progress made by these models, they have revealed problems like small density ratio, large interface spurious velocity, small reduced temperature, etc, [6, 7] which restrict their practical applications.

Peng et al.[6] conducted simulations to compare several different EOSs in dealing with the problems mentioned above. A density ratio over 1000 was obtained. His work showed that the modification of interaction potential based on different EOSs would be helpful to simulate multiphase problems with a density ratio close to the real liquid– gas system. Kupershtokh also made a comparison between some EOSs.[810] In his method, a relatively large lattice speed was used to acquire good numerical stability of lattice Boltzmann (LB) algorithm. Besides, the scheme was validated with strict mathematical derivation. Liu et al.[11] analyzed the pseudopotential model and showed that a relatively small ∂ P/∂ ρ ratio of EOS would result in a better numerical stability.

Usually, when dealing with multiphase problems, the surface tension effect is very important, especially in microscopic problems. Using the mean field approximation for intermolecular attraction and Enskog’ s treatment of the exclusion-volume effect, He et al.[12] proposed the surface tension term. Later, Peng modified the term and conducted the Laplace rule test.[7] However, it should be noted that no theoretical derivation or validation by methods like molecular dynamics (MD) was given in Yuan’ s paper. Hence, further work considering the surface tension term needs to be done and compared as well as the numerical scheme.

In this paper, we concentrate on the pseudopotential model, incorporate the surface tension term, and deal with typical multiphase problems, namely, the Laplace rule, wall adhesion, and heterogeneous cavitation.

2. Methodology
2.1. Framework of LBM

Lattice Boltzmann equations (LBEs) can be derived from the BGK approximation of the Boltzmann equation. In the standard LBM without external forces, the evolution of the particles is described as

where τ is the dimensionless relaxation time. We adopt the D2Q9 model and the equilibrium distribution function is expressed as

where c is the lattice speed, and ω i is the weight factor equal to ω 0 = 4/9, ω 1− 4 = 1/9, and ω 5− 8 = 1/36 for the D2Q9 model. The equilibrium velocity ueq is derived from the interacting force, including cohesion force, wall adhesive force, and external force. Using the Shan– Doolen model to incorporate the forcing term, [13, 14] we have

The macroscopic density and velocity are expressed as

2.2. Multiphase pseudopotential model

The original pseudopotential model was proposed by Shan and Chen.[1] They approximated the interparticle force by using the nearest neighbor interaction potential, and the force was expressed as

where g is the interparticle interaction strength, and the effective mass ψ (ρ ) is formulated as

The parameter g in Eqs. (5) and (6) is adjusted to ensure that the whole term inside the square root is positive. Similarly, the wall adhesive force takes the form of

where gads is the fluid– solid interaction strength. The solid wall is regarded as a phase with constant density s, which is 1 for a solid node and 0 for a fluid node.

To implement Eq. (6), we apply a proper EOS, which describes the relationship between pressure, density (or volume), and temperature for a given substance. In our work, we adopt the CS EOS to acquire larger density ratio and much smaller spurious interface velocity, [6, 8, 15] which is formulated as

To satisfy the critical point conditions, we have

By setting Tc = 1, ρ c = 1, and Pc = 0.01, the parameters are calculated as a = 0.03853462257, b = 0.1304438842, and R = 0.02785855166. In simulation cases of this paper, the reduced temperature is Tr = T/Tc = 0.59. In addition, we introduce a dimensionless coefficient proposed first by Kupershtokh[8]

Here K represents a crucial factor in the numerical stability of multiphase cases, where smaller K generally results in better numerical stability. Instead of decreasing K by decreasing the δ t/δ x ratio, namely increasing the lattice speed c, we keep a uniform lattice velocity, which is a common practice in most LB simulations and can avoid too many time steps in unsteady problems.

To theoretically predict the density of both liquid phase and gas phase, we solve the Maxwell construction using the integral technique.[16] The density ratio of up to approximately 160000 (Tc = 0.55) is achieved successfully, and further adjustment of K will lead to an even larger density ratio in our tests. The LBM simulated results fit well with the Maxwell construction in Fig. 1.

Fig. 1. Coexistence curve of theoretical and calculated density of the CS EOS.

2.3. Surface tension term

Physically, the surface tension is caused by the attraction between molecules of fluid, and it is an effect within the surface layer of a fluid that causes the layer to behave as an elastic sheet.[11] The original pseudopotential model does not explicitly contain a surface tension term. Shan and Chen showed that the surface tension term could be exclusively formulated as[12]

where κ determines the strength of surface tension. This scheme offers a flexible choice to adjust the surface tension strength in situations where the surface tension is important.

Numerically, we refer to the finite difference method (FDM) and use the five-point scheme for the calculation of 2ρ

The lattice discretized scheme for the gradient term ∇ ∇ 2ρ (assuming ϕ = 2ρ ) is

By expanding the right term, it can be proved that equation (13) is identical to the six-point method, [6] while remaining an easy form to calculate in the lattice evolution loop. By incorporating the surface tension term in the interparticle interaction force, equation (5) is implemented as follows:

3. Results and discussion
3.1. Laplace rule

The typical Laplace rule (Fig. 2) indicates that the pressure difference across the drop/bubble interface has an inverse proportional relation to the radius, namely,

Fig. 2. Droplet at T = 0.59Tc. (a) Density distribution; (b) spurious velocity in the x direction.

Figure 2 shows the density and velocity distribution, and we should note that the spurious velocity is unavoidable due to the modeling of multiphase interaction. Besides, the maximum spurious velocity is less than the maximum spurious velocity (namely, 0.092571) using the CS EOS when κ = 0 in Ref. [6]. In Fig. 3, drops/bubbles of various sizes are simulated, and each line satisfies the linear relation well with the slope equal to the surface tension coefficient σ . By changing κ in Eq. (11), we can increase σ by up to 55.4% and decrease σ by about 2% without modifying the initial density ratio.

Fig. 3. Plot of curvature (1/R) versus pressure difference.

Another form of surface tension term proposed by Peng is described as Fs = κ ψ (ρ )∇ ∇ 2ψ (ρ ).[7] However, our simulation results show that the surface tension coefficient increases by 19.2% at most and decreases by 15.0% if the initial density ratio is not modified.

It should be noted here that the interface becomes thicker once the surface tension is strengthened, as shown in Table 1. Since the pseudopotential model is a diffuse-interface method, it is usual that the interface occupies several lattice units. Huang et al. constructed an analytical expression of density distribution of Laplace rule[14]

Table 1. Properties of droplets (R ≈ 50) for different κ . w denotes the interface thickness.

where (x0, y0) is the center position of domain, and the radius of the droplet is implicitly located, namely, ρ = (ρ l + ρ v)/2. Figure 4 shows the coexistence curve in the x direction, and the curve in the range of 0 ≤ x ≤ 100 is just symmetrical to that of 100 ≤ x ≤ 200.

Fig. 4. Coexistence curve of LBM simulated and theoretical density distribution in the x direction.

Since the cohesive force in Eq. (5) is the main force that determines the phase interaction, our analysis will concentrate on it. Equation (5) is a discretized form. In fact, there is another equivalent form

In the x direction, we have

In our work, the ratio of critical pressure and critical density is set to be 0.01, ensuring a small ∂ p/∂ ρ ratio. Moreover, a thicker interface means a smaller density gradient (∂ ρ /∂ x) around the interface, making the interaction force smaller. Usually, the density of vapor phase is small and the velocity of vapor phase is easy to vibrate even under small perturbation, which causes the spurious velocity and numerical instability. In our experience, when numerical instability occurs, density changes sharply near the interface, making the interface rather thin. The above analysis can also explain the situation that numerical stability becomes better when κ becomes larger, since the density of vapor phase is larger.

3.2. Wall adhesion

According to Young’ s equation, the contact angle of a droplet on a smooth plate board will satisfy

where σ l, v is the liquid– vapor interfacial tension, σ s, v and σ s, l are the solid– vapor and solid– liquid interaction strength, respectively. Usually, the contact angle is determined by the solid– liquid interaction strength gads only.[17] By changing κ , we can also acquire different wettability, or hydrophobic/hydrophilic property, on the solid wall, as shown in Fig. 5.

Fig. 5. Wall adhesion (the gray arrows denote the velocity vector, with the reference vector of unit magnitude) (a) gads = − 2.0, κ = 0, and θ = 77° ; (b) gads = − 2.0, κ = 0.02, and θ = 114.2.

In our work, although the curve in Fig. 5(b) does not fit linearly as Fig. 5(a), we can also obtain the corresponding relationship before simulating cases involving wall contact angle parameters. Besides, by changing κ , numerical stability is improved.

3.3. Heterogeneous cavitation

Cavitation is a catastrophic transition from liquid to vapor, and it is called heterogeneous cavitation when it is nucleated by pre-existing bubbles or other disruptions in the structure of the liquid.[18] Simulations on heterogeneous cavitation usually concentrated on the critical radius[19, 20] according to the energy theory by Ou and Truller, [20] who gave

where R* is the critical radius, and Δ P is the outer pressure difference acting on the bubble. If given the pressure difference, a bubble with a radius less than R* will condensate, while a bubble with a radius larger than R* will nucleate cavitation. Our simulations deal mainly with different phenomena caused by different surface tension strength. In all cases, the density of the liquid– vapor system is (3.116, 0.08), and the pressure difference is 5.162× 10− 3. To gain different surface tension strengths, the parameter κ in the surface tension term is made 0 and 0.06, individually.

From Figs. 7 and 8, we know that the radius of a bubble changes with time. It should be noted that the radius in every time step is obtained through tracking the interface where the density of the node is closest to 0.5(ρ l + ρ v). Thus, only integral data are captured, and the tendency is drawn clearer with polynomial fitting.

Fig. 6. Coexistence curve of calculated and fitted data of contact angles. (a) κ = 0 and gads varies from − 2.2 to − 1.6. (b) gads = − 2.0 and κ varies from − 0.002 to 0.06.

Fig. 7. Effect of surface tension on evolution. (a)– (d) κ = 0, bubble nucleates cavitation. (a) t = 200, (b) t = 2000, (c) t = 3000, (d) t = 3440. (e)– (h) κ = 0.06, bubble dissipates due to condensation. (e) t = 200, (f) t = 500, (g) t = 800, (h) t = 940.

Fig. 8. Radius of bubble against time. (a) κ = 0; (b) κ = 0.06.

When the surface tension is strong, more energy for the bubble is needed to expand and it has a tendency to condensate. On the contrary, bubbles will nucleate cavitation. Our simulation results are consistent with the theoretical prediction that the surface tension coefficient lies between 0.09213 and 0.14607 (see Table 1), which corresponds to κ . Thus equation (20) is validated from the prospect of surface tension that the original pseudopotential cannot do.

4. Conclusions

In summary, the pseudopotential multiphase model has been investigated and modified based on the CS EOS. With the parameters proposed in our work, large liquid– vapor density ratio can be achieved, which can be up to 160000 when the reduced temperature is 0.55. Our modification deals mainly with the surface tension term, which offers a flexible way to adjust surface tension strength in liquid– vapor systems. The spurious velocity around the liquid– vapor interface is reduced and the numerical stability is improved if the surface tension strength increases. Furthermore, simulations on wall adhesion and heterogeneous cavitation are conducted, showing the practical application of the surface tension term in dealing with liquid– vapor problems.

Reference
1 Shan X and Chen H 1993 Phys. Rev. E 47 1815 DOI:10.1103/PhysRevE.47.1815 [Cited within:3] [JCR: 2.313]
2 Michael R S, Osborn W R and Yeomans J M 1995 Phys. Rev. Lett. 75 830 DOI:10.1103/PhysRevLett.75.830 [Cited within:1] [JCR: 7.943]
3 Haihu L, Albert J V, Qinjun K and Charles W 2013 Transp Porous Med. 99 555 [Cited within:1]
4 Haihu L and Yonghao Z 2010 J. Comput. Phys. 229 9166 DOI:10.1016/j.jcp.2010.08.031 [Cited within:1] [JCR: 2.138]
5 Haihu L and Albert J V 2012 Phys. Rev. E 85 046309 DOI:10.1103/PhysRevE.85.046309 [Cited within:1] [JCR: 2.313]
6 Peng Y and Schaefer L 2006 Phys. Fluids 18 042101 DOI:10.1063/1.2187070 [Cited within:]
7 Peng Y 2005 Thermal Lattice Boltzmann Two-Phase Flow Model for Fluid Dynamics(Ph. D. Dissertation) Pennsylvania University of Pittsburgh [Cited within:3]
8 Kupershtokh A L, Medvedev D A and Karpov D I 2009 Comput. Math. Appl. 58 965 DOI:10.1016/j.camwa.2009.02.024 [Cited within:3] [JCR: 0.75]
9 Kupershtokh A L 2010 Comput. Math. Appl. 59 2236 DOI:10.1016/j.camwa.2009.08.058 [Cited within:1] [JCR: 0.75]
10 Kupershtokh A L 2011 Comput. Math. Appl. 61 3537 DOI:10.1016/j.camwa.2010.06.032 [Cited within:1] [JCR: 0.75]
11 Liu M L, Yu Z, Wang T F, Wang J F and Fan L S 2010 Chemical Engineering Science 65 5615 DOI:10.1016/j.ces.2010.08.014 [Cited within:2] [JCR: 2.386]
12 He X, Chen S and Zhang R 1999 J Comput Phys 152 642 DOI:10.1006/jcph.1999.6257 [Cited within:2] [JCR: 2.138]
13 Shan X and Doolen G D 1995 J. Stat. Phys. 81 379 DOI:10.1007/BF02179985 [Cited within:1] [JCR: 1.4]
14 Huang H, Krafczyk M and Lu X 2011 Phys. Rev. E 84 046710 DOI:10.1103/PhysRevE.84.046710 [Cited within:2] [JCR: 2.313]
15 Chen X, Zhong C and Yuan X 2011 Comput. Math. Appl. 61 3577 DOI:10.1016/j.camwa.2010.07.018 [Cited within:1] [JCR: 0.75]
16 Elhassan A E, Craven R J B and de Reuck K M 1997 Fluid Phase Equilibria 130 167 DOI:10.1016/S0378-3812(96)03222-0 [Cited within:1] [JCR: 2.379]
17 Benzi R, Biferale L, Sbragaglia M, Succi S and Toschi F 2006 Math. Comput. Sim. 72 84 DOI:10.1016/j.matcom.2006.05.034 [Cited within:1]
18 Sukop M and Or D 2005 Phys. Rev. E 71 046703 DOI:10.1103/PhysRevE.71.046703 [Cited within:1] [JCR: 2.313]
19 Zhong M, Zhong C W and Bai C R 2012 ACSA 1 73 [Cited within:1]
20 Zhang M, Zhou C Y, Shams I and Liu J Q 2009 Acta Phys. Sin. 58 8406(in Chinese) [Cited within:2] [JCR: 1.016] [CJCR: 1.691]
21 Or D and Tuller M 2002 Water Resour. Res. 38 1601 [Cited within:1] [JCR: 3.149]