Isotropy and spurious currents in pseudo-potential multiphase lattice Boltzmann models

TL;DR

This paper investigates the causes of spurious currents in multiphase lattice Boltzmann models, highlighting limitations of increasing isotropy and proposing a new implementation to better suppress these artifacts.

Contribution

It introduces a consistent implementation that more effectively reduces spurious currents and clarifies the limitations of isotropy enhancement in eliminating them.

Findings

Higher isotropy alone cannot eliminate certain low-level spurious currents.

A new implementation improves suppression of spurious currents.

Theoretical analysis shows the importance of local hydrostatic balance.

Abstract

The spurious currents observed in multiphase flow simulations with pseudo-potential lattice Boltzmann (LB) models are usually understood to be the result of the lack of isotropy of the model-generated interaction force between phases. Remedies have been proposed to utilize larger stencils to compute the interaction force with higher orders of isotropy. In this short communication, we point out the incompleteness in the current understanding and propose a new consistent implementation to more effectively suppress the spurious currents. We also demonstrate theoretically that certain low-level spurious currents cannot be eliminated by increasing isotropy if the local hydrostatic balance inside the diffuse interface is not established in the LB models.

Tables

Table 1. Table 1 : Parameters in the D2Q9, D2Q13 and D2Q25 models

	$c_s^2$	$w\left(0\right)$	$w\left(1\right)$	$w\left(2\right)$	$w\left(3\right)$	$w\left(4\right)$	$w\left(5\right)$	$w\left(6\right)$	$w\left(7\right)$	$w\left(8\right)$
D2Q9	1/3	4/9	1/9	1/36
D2Q13	2/5	2/5	8/75	1/25		1/300
D2Q25	4/7	72/245	16/147	16/315		1/105	8/2205			1/8820

Equations

$$\mathbf{F}\left(\mathbf{x},t\right)=-G\psi \left(\mathbf{x},t\right)\sum _{\alpha }{w}_{\alpha }\left(|{\mathbf{e}}_{\alpha }|\right)\psi \left(\mathbf{x}+{\mathbf{e}}_{\alpha }\delta t,t\right){\mathbf{e}}_{\alpha },$$

$$\begin{array}{rl}& F_i=-G\psi \sum _{\alpha }{w}_{\alpha }{e}_{\alpha i}[\psi +\delta t{e}_{\alpha j}{\partial }_j\psi \\ & +\frac{1}{2}\delta t^2{e}_{\alpha j}{e}_{\alpha k}{\partial }_j{\partial }_k\psi +\frac{1}{6}\delta t^3{e}_{\alpha j}{e}_{\alpha k}{e}_{\alpha l}{\partial }_j{\partial }_k{\partial }_l\psi +\cdots ].\end{array}$$

$$\begin{array}{rl}& f_{\beta }\left(\mathbf{x},t+\delta t\right)=f_{\beta }\left(\mathbf{x}-\mathbf{c}\delta t,t\right)\\ & -\frac{1}{\tau }\left[f_{\beta }\left(\mathbf{x}-\mathbf{c}\delta t,t\right)-f_{\beta }^{(eq)}\left(\mathbf{x}-\mathbf{c}\delta t,t\right)\right]+{\varphi }_{\beta }\left(\mathbf{x}-\mathbf{c}\delta t,t\right)\end{array}$$

$$\begin{array}{rl}& f_{\beta }^{(eq)}=\rho w_{\beta }\left[1+\frac{\mathbf{c}\cdot \mathbf{u}}{c_s^2}+\frac{{\left(\mathbf{c}\cdot \mathbf{u}\right)}^2}{2c_s^4}-\frac{\mathbf{u}\cdot \mathbf{u}}{2c_s^2}\right],\\ & {\varphi }_{\beta }=\left(1-\frac{1}{2\tau }\right)w_{\beta }\left[\frac{\mathbf{c}-\mathbf{u}}{c_s^2}+\frac{\left(\mathbf{c}\cdot \mathbf{u}\right)\mathbf{c}}{c_s^4}\right]\cdot \mathbf{F}\delta t,\end{array}$$

$$f_{\beta }^{(eq)}=\rho w_{\beta },{\varphi }_{\beta }=\left(1-\frac{1}{2\tau }\right)w_{\beta }\frac{\mathbf{c}\cdot \mathbf{F}}{c_s^2}\delta t,$$

$$\begin{array}{rl}& \rho \left(\mathbf{x},t+\delta t\right)=\sum _{\beta }{w}_{\beta }\left[\rho \left(\mathbf{x}-\mathbf{c}\delta t,t\right)+\frac{1}{2}\frac{c_{\beta i}{F}_i\left(\mathbf{x}-\mathbf{c}\delta t,t\right)}{c_s^2}\delta t\right]\\ & =\sum _{\beta }{w}_{\beta }\left[\rho +\frac{1}{2}\delta t^2{c}_{\beta j}{c}_{\beta k}{\partial }_j{\partial }_k\rho +\frac{1}{24}\delta t^4{c}_{\beta j}{c}_{\beta k}{c}_{\beta l}{c}_{\beta m}{\partial }_j{\partial }_k{\partial }_l{\partial }_m\rho +\cdots \right]\\ & -\frac{1}{2}\sum _{\beta }{w}_{\beta }\left[\frac{\delta t^2}{c_s^2}{\partial }_j{F}_i{c}_{\beta i}{c}_{\beta j}+\frac{\delta t^4}{6c_s^2}{\partial }_j{\partial }_k{\partial }_l{F}_i{c}_{\beta i}{c}_{\beta j}{c}_{\beta k}{c}_{\beta l}+\cdots \right],\end{array}$$

$$\begin{array}{rl}& \rho u_i\left(\mathbf{x},t+\delta t\right)\\ & =\sum _{\beta }{w}_{\beta }\left[\rho \left(\mathbf{x}-\mathbf{c}\delta t,t\right)+\frac{1}{2}\frac{c_{\beta j}{F}_j\left(\mathbf{x}-\mathbf{c}\delta t,t\right)}{c_s^2}\delta t\right]c_{\beta i}+\frac{1}{2}{F}_i\left(\mathbf{x},t\right)\delta t\\ & =\sum _{\beta }\left(-w_{\beta }{c}_{\beta i}{c}_{\beta j}\delta t{\partial }_j\rho -\frac{1}{6}{c}_{\beta i}{c}_{\beta j}{c}_{\beta k}{c}_{\beta l}\delta t^3{\partial }_j{\partial }_k{\partial }_l\rho +...\right)\\ & +\frac{1}{2}\sum _{\beta }\left(\frac{\delta t{F}_j}{c_s^2}{w}_{\beta }{c}_{\beta i}{c}_{\beta j}+\frac{\delta t^3{\partial }_k{\partial }_l{F}_j}{2c_s^2}{c}_{\beta i}{c}_{\beta j}{c}_{\beta k}{c}_{\beta l}+\cdots \right)+\frac{1}{2}{F}_i\delta t\\ & =\delta t\left(F_i-c_s^2{\partial }_i\rho \right)+{\delta }_t^3\frac{3c_s^2}{4}{\partial }_k{\partial }_k\left(F_i-\frac{2c_s^2}{3}{\partial }_i\rho \right)+...\end{array}$$

Full text

††journal: Computers & Fluids

1 Introduction

The pseudo-potential multiphase lattice Boltzmann (LB) models, (also known as the Shan-Chen models [1, 2], have been widely applied to study a wide range of multiphase flow problems. Despite their successes, the existence of non-physical flux around a steady and static two-phase interface, known as spurious currents, still plagues most multiphase applications and remains an unresolved problem.

Many efforts were made to explore the origin of these spurious currents and to suppress them. Those efforts are comprehensively reviewed in the literature [3, 4]. Wagner pointed out that spurious currents were introduced by an incompatible discretization of the interaction force in the two-phase LB model [5]. Shan [6] and Li and Fischer [7] both realized that when the discretization schemes of the interaction force in the multiphase LB models lack isotropy, spurious currents would emerge. Yuan and Schaefer reported that certain equation of state (EOS) could potentially reduce the level of spurious currents [8]. Yu and Fan found that, compared to the single-relaxation-time (SRT) LB models, the multiple-relaxation-time (MRT) LB models could be used to suppress the spurious currents by tuning the relaxation parameters irrelevant to the Navier-Stokes equation [9]. Guo et al. concluded that spurious currents were inevitable due to intrinsic imbalance of interaction force and the density gradient in the pseudo-potential multiphase LB models [10]. Mattila et al. suspected that the existence of spurious currents could be associated with the second-order accuracy of LB models due to the trapezoidal time-integration scheme, they therefore proposed to use higher-order LB models to suppress the spurious currents.

2 Shan’s improvement and its incompleteness

Among all these explanations, a rather well-known explanation of the origin of spurious currents in the psuedo-potential LB models was given by Shan, who realized that the high-order terms in the Taylor series of the interaction force, i.e., the interaction force exerted on one phase due to the existence of another phase around, lack the required isotropy [6]. The interaction force ${\bf F}$ in the pseudo-potential multiphase LB models is computed as

[TABLE]

where $G$ is a parameter measuring the intensity of the interaction, $\psi$ is the field potential that is a function of local fluid density $\rho$, ${\bf e}_{\alpha}$ is the vector stencil employed in the computation of interaction force, which is not necessarily the same as the discrete velocity set in LB models, $w_{\alpha}$ is the corresponding weighting factor and $\delta t$ is the time step size. Eq. (1) can be expanded in terms of a Taylor series at ${\bf x}$ and $t$, using tensor notations, as

[TABLE]

Except the first term, each term in the above Taylor series contains a part $\sum_{\alpha}w_{\alpha}\underbrace{e_{\alpha i}e_{\alpha j}e_{\alpha k}e_{\alpha l}\cdots}_{n~{}e_{\alpha}}$, which is a $n$th order tensor. Shan pointed out that spurious currents were originated from the lack of complete isotropy of these high-order tensors, when the number of ${\bf e}_{\alpha}$ is finite. As a remedy, Shan employed larger stencils to compute ${F}_{i}$, which allowed additional tensors to be isotropic and increased the order of isotropy in the computed interaction force $F_{i}$. For example, in two space dimensions, the highest order of the isotropy realizable can be increased from fourth with the velocity stencil shown in Fig. 1a to eighth with the stencil in Fig. 1c.

In this short communication, we would like to point out the incompleteness in Shan’s recommended remedy and its implementation. In fact, there is a second aspect in terms of isotropy that has usually been ignored but plays an important role in inducing spurious currents. To explain this, let us recall the algorithm of the LB method

[TABLE]

where ${\bf c}_{\beta}$ is the discrete velocity set in LB model that may be different from ${\bf e}_{\alpha}$, $\tau$ is the relaxation time. The equilibrium distribution $f_{\beta}^{(eq)}$ and the forcing function $\phi_{\beta}$ are defined as

[TABLE]

where $c_{s}$ is the speed of sound, which is an input parameter in a specific LB model based on numerical quadrature requirements. The forcing function in Eq. (4) is the one proposed by Guo [11], which ensures a second-order accurate body force term in the reproduced Navier-Stoke equation. We further assume a zero velocity field ${\bf u}(t)=0$ is reached at the current time $t$, then Eq. (4) is simplified as

[TABLE]

Finally, for demonstration purposes, a special case $\tau=1$ is assumed, which allows great simplification for mathematics. The density at $t+\delta t$ is calculated as

[TABLE]

where all the odd-order terms in the Taylor series are zero due to the symmetry of the discrete velocity set. Physically, we should then require the $n$th order tensors $\sum_{\beta}w_{\beta}\underbrace{c_{\beta i}c_{\beta j}c_{\beta k}c_{\beta l}\cdots}_{n~{}c_{\beta}}$ to be isotropic, as otherwise the density field $\rho(t+\delta t)$ would have already contained errors due to deviations from isotropy. Having established that, it is straightforward to realize the incompleteness of Shan’s remedy: without concurrently enforcing the isotropy in Eq. (6) to a similar order at the same time, the isotropy of a multiphase LB simulation would not be effectively improved. This is probably the reason why spurious currents were not reduced as significantly as expected with Shan’s remedy, only by a factor of 3 when the isotropy of the interaction force was increased from 4th-order to 8th-order [6].

3 A complete remedy for the lack of isotropy

Essentially, this second aspect of isotropy concerns the distribution of the interaction force back to the lattice nodes, while the first aspect in Shan’s analysis concerns the calculation of the interaction force. This second isotropy is determined by $w_{\beta}$ and ${\bf c}_{\beta}$, which could not be remedied without expanding the discrete velocity sets in the LB model. It has been proven that the nearest-neighboring LB models, such as D2Q9, D3Q15, D3Q19, and D3Q27 can only achieve a fourth-order isotropy with their discrete velocity sets [6]. However, with the stencils in Fig. 1b and Fig. 1c, we can easily construct a D2Q13 model and a D2Q25 model to increase the second isotropy to 6th- and 8th-order, respectively. These two models can have the same form of equilibrium distribution function and forcing function as the regular LB models, except that the weighting factor $w_{\beta}$ and the speed of sound $c_{s}$ have to be redefined. These parameters are given in Table 1. Other models allowing even higher-orders of isotropy can be formulated using the same philosophy. The velocity stencils and the corresponding weighting factors given in literature [6, 12] can be used as references, the remaining job is to design the equilibrium distribution function and the forcing function i.e., $f_{\beta}^{(eq)}$ and $\phi_{\beta}$ to satisfy the constraints that lead to the reproduction of the Navier-Stokes equation. It is worth mentioning that Mattila et al. also mentioned that LB models with larger sets of discrete velocities could be employed to reduce the level of spurious currents [13]. However, their models were designed to incorporate higher-order time-integration schemes to replace the trapezoidal rule in standard LB models rather than to maximize the order of isotropy. Therefore, their models with the largest number of discrete velocities did not show the most significant reduction of spurious currents. The relationship between the lack of isotropy in the distribution of interaction force and the appearance of spurious currents was not explicitly stated.

We employ a simple test case of a droplet suspended in a 2D periodic domain in vapor. The grid resolution of the test is $60\times 60$, and a droplet with an initial radius of $r_{0}=15$ is placed at the center of the domain $(x_{c},y_{c})$. The initial density distribution is defined as $\rho_{0}(x,y)=\frac{\rho_{l}+\rho_{v}}{2}-\frac{\rho_{l}-\rho_{v}}{2}\tanh{\left[\frac{2\sqrt{\left(x-x_{c}\right)^{2}+\left(y-y_{c}\right)^{2}-r_{0}}}{W}\right]}$, where $\rho_{l}$ and $\rho_{v}$ are the liquid and vapor density, respectively, at a given temperature $T$ below the critical temperature $T_{c}$. The equation of state (EOS) used in the simulations is the Peng-Robinson (P-R) EOS [14], $p=\frac{\rho RT}{1-b\rho}-\frac{a\alpha(T)\rho^{2}}{1+2b\rho-b^{2}\rho^{2}}$, where $a=2/49$, $b=2/21$, $R=1$, as defined in Yuan and Schaefer [8], $\alpha(T)=[1+(0.37464+1.54226\omega-0.26992\omega^{2})(1-\sqrt{T/T_{c}})]^{2}$, $\omega$ is chosen to be 0.344 for water. The pseudo-potential function $\psi$ is calculated as $\psi=\sqrt{\frac{2(p-c_{s}^{2}\rho)}{Gc_{s}^{\prime 2}\delta t}}$, $c_{s}^{2}=\sum_{\beta}w_{\beta}c_{\beta i}c_{\beta i}$, $c_{s}^{\prime 2}=\sum_{\alpha}w_{\alpha}e_{\alpha i}e_{\alpha i}$. The exact difference method (EDM) [15] is adopted as the forcing scheme to distributed the interaction force in LB models. At a reduced temperature $T/T_{c}=0.725$, the steady state density contours and velocity fields generated with the standard D2Q9 model, the modified D2Q9 model under Shan’s improvement with the eighth-order isotropy to compute the interaction force, and D2Q25 model are shown in Fig. 2. Clearly, although Shan’s improvement was able to reduce the magnitude of the spurious velocity to certain extent, it is not so effective compared to the D2Q25 model. More importantly, the azimuthal-dependent flow patterns still exist with Shan’s best improvement, which indicates that there are remaining errors due to the lack of isotropy. The same flow patterns are no longer exist with the D2Q25 model, the remaining spurious currents are almost perpendicular to the interface.

To quantify the reduction of spurious currents by the proposed models, we calculated the maximum spurious velocity and the field-averaged spurious velocity for a larger range of the reduced temperature. These results are shown in Fig. 3. For D2Q13 and D2Q25 models, two different relaxation times are used. One is chosen identical to the relaxation time in the D2Q9 models $\tau=1.0$, which are labeled as D2Q13T and D2Q25T. The other is designed to result in the same viscosity as in the D2Q9 models, i.e., $\tau=11/12$ for the D2Q13 model and $\tau=19/24$ for the D2Q25 model. The latter two simulations are labeled as D2Q13V and D2Q25V. For the examined reduced temperature range, $0.6\leq T/Tc\leq 0.95$, which covers a range of liquid-to-vapor density ratio (from 4 to 800), D2Q13 and D2Q25 models always reduce the magnitude of the spurious velocity by another order of magnitude compared to Shan’s improvement. It is worth noticing that Shan’s improvement becomes ineffective at small reduced temperatures. As we shall see shortly, this is probably because Shan’s improvement increases a thermodynamic inconsistency that offsets the benefit of increasing the isotropy in calculating the interaction force.

The thermodynamic inconsistency is another critical issue in pseudo-potential models. The lack of thermodynamic inconsistency can usually be seen from the deviations of the numerically obtained liquid and vapor densities from the corresponding values obtained by the Maxwell equal area rule, for a given EOS of a pure substance. In the literature, there are many attempts to quantify the magnitude of such derivations in the suspending droplet case shown above. However, the droplet case should not be used to evaluate such deviations from the results of the Maxwell equal area rule, as the curved liquid-vapor interface in this case leads to different pressures in the liquid and vapor bulk phases while the Maxwell equal area rule is based on a same pressure in the two phases. A more appropriate case to measure the thermodynamic inconsistency in the pseudo-potential models is a 1D flat interface case. In this case, we have measured the numerically obtained liquid and vapor densities at different reduced temperatures with P-R EOS. The new models (D2Q13, D2Q25) in general do not alter the coexisting liquid-vapor densities significantly. In fact, they slightly improve the liquid-vapor densities at smaller reduced temperatures compared to Shan’s improvements (D2Q9 ISO6, D2Q9 ISO8). The larger deviations from the thermodynamic consistency with Shan’s improvement may explain its failure to reduce the spurious velocity at small reduced temperatures observed in Fig. 3. It is worth mentioning that there are also many available ways to improve the thermodynamic consistency in the literature, such as using a coupled form to construct the intermolecular force [15], modifying EOS [16], and introducing correction terms in the forcing schemes of LBM [17]. These available improvements can be used to improve the thermodynamic consistency in our models. Therefore, we do not address the issue of thermodynamic inconsistency in the present study.

The major side effect of using D2Q13 and D2Q25 modes is the increased interface thickness. In the 1D flat interface case, the interface thickness, defined as the region with $1.05\rho_{v}\leq\rho\leq 0.95\rho_{l}$, where $\rho_{v}$ and $\rho_{l}$ are the numerically obtained vapor and liquid densities in the two-phase bulk regions with each model, are measured and shown in Fig. 5. The increased interface thickness comes from two aspects. First, by using large stencils to compute the interaction force in the pseudo-potential models, the local force depends on the potential $\psi$ from more neighboring grid points, which makes the interface broader. This aspect also impacts the interface thickness when using Shan’s improvement, but it has not been emphasized in the literature. Second, the use of large lattice models allows a local grid point to directly communicate with not just the nearest neighboring nodes, but also the next layer of neighboring nodes, which adversely affects the locality of LBM. This aspect also makes the interface thicker. If the thickening interface has to be avoided in a certain application, the D2Q13 model could be a better choice compared to the D2Q25 model.

Finally, we would like to briefly comment on the no-slip boundary treatment for the proposed models. In general, the no-slip boundary can still be treated following the bounce-back schemes. However, the uses of certain bounce-back schemes, such as the half-way bounce-back may not be straightforward, as it is difficult to place a solid wall precisely half-way for all links. On the other hand, schemes such as the standard bounce-back [18] and modified bounce-back [19] are not affected. Additional attention to be paid is that the unknown boundary distribution functions have to be constructed on the first two layers of interior grid points. As a demonstration, a test case of a droplet contacting with a flat wall has been added. The P-R EOS is still used and the temperature is set to $T=0.9T_{c}$. The standard bounce-back scheme is adopted to enforce the no-slip condition. The varying wettability in the two new models can still be achieved by tuning the virtual density of the wall $\rho_{w}$ as in Ref. [20].

4 Conclusion and discussion

In this work, we point out the incompleteness of the previous understanding of the spurious currents in the multiphase flow simulation with the pseudo-potential LB models. There are two types of isotropy requirements: the first concerns the macroscopic force calculation, and the second is the mesoscopic redistribution. The two should be considered together in order to more effectively reduce spurious currents. We proposed two LB models with more discrete velocities that can reduce spurious currents by one order of magnitude.

The remaining spurious velocities are associated with the discrete nature in LB models. Following the same analysis in Eq. (6), the momentum at $({\bf x},t+\delta t)$ can be computed as

[TABLE]

Regardless of whether the truncated high-order tensors are isotropic or not, $\rho u_{i}({\bf x},t+\delta t)$ appears always non-zero since the first two terms cannot be zero at the same time. Even if we only keep the leading-order term in Eq. (7), the coefficient $\delta tF_{i}/c_{s}^{2}-\delta t\partial_{j}\rho\delta_{ij}$ may not be precisely zero. Ideally, on the N-S equation level this coefficient should vanish when a hydrostatic balance is established, but once discretized, the precise balance is usually violated. The spurious currents due to this imbalance are aligned along with the direction of pressure gradient, which corresponds to what we observed in Fig. 2c. On a 1D flat interface, this spurious current further reduces to what Guo et al. reported as inevitable artificial velocities in LB simulations [10]. The pseudo-potential LB models calculate the interaction force in a discretized form and again distribute this force to the discrete distribution functions, both inducing errors that could result in spurious currents. A consistent consideration of the two processes together may help further suppress or even remove the spurious current, which will be pursued in the future.

Acknowledgements: Authors would like to thank Prof. Zhaoli Guo at Huazhong University of Science and Technology, China, and Prof. Haibo Huang at University of Science and Technology of China for the inspiring discussions. Funding support from Energi Simulation and the William A. Fustos Family Professorship in Energy and Mineral Engineering at Penn State University are gratefully acknowledged.