Density gradient calculation in a class of multiphase lattice Boltzmann models

The multiphase lattice Boltzmann (LB) models based on pairwise interactions show great potential for simulating multiphase flows due to the conceptual and computational simplicity. Although the dynamics of multiphase flows are reproduced by the pairwise interaction force, the gradient of density (or effective density, i.e., pseudopotential) is implicitly involved in these models via the specialized forcing scheme or the consistent scheme for ${\varepsilon }^3\text{-order}$ term. This work focuses on the calculation of density gradient in this class of multiphase LB models. Theoretical analyses are first carried out to reveal the involvement and calculation of density gradient. On the basis of a low Mach number approximation, an improved scheme is then proposed to calculate the density gradient for the recent LB model with self-tuning equation of state. Analytical and numerical calculations show that the improved scheme is more accurate and can help to reduce the numerical error when the reduced temperature is relatively low.

Figure 1

Numerical results for the simulations of one-dimensional flat interface by the pseudopotential LB model with consistent scheme for ${\varepsilon }^3\text{-order}$ term. (a) The coexistence curves obtained by simulation and the analytical ones given by mechanical stability condition. (b) The relative error of gas density $E_{\rho ,\text{gas}}=|{\rho }_{g,\text{numerical}}^{}-{\rho }_{g,\text{analytical}}^{}|/{\rho }_{g,\text{analytical}}^{}$ versus the analytical density ratio ${\rho }_{l,\text{analytical}}^{}/{\rho }_{g,\text{analytical}}^{}$.

Figure 2

Comparisons between the denominator $G^2{\delta }_x^2\rho +\frac{2}{3}\frac{d{p}_{\text{EOS}}^{}}{d\rho }$ of Eq. (28) and the denominator $G^2{\delta }_x^2\rho$ of Eq. (7). The scaling factor $K_{\text{INT}}$ is determined by prescribing $W=5{\delta }_x$, and the scaling factor $K_{\text{EOS}}$ is fixed at 1. In the $\rho -T_r$ plane, the coexistence curve by Maxwell construction is shown by the solid line and the region with $\frac{d{p}_{\text{EOS}}^{}}{d\rho }<0$ is depicted by the gray area.

Figure 3

The profiles of normalized density ${\rho }^{*}$ and density gradient $d{\rho }^{*}/d{x}^{*}$ across phase interface. The solid lines are the analytical profiles, the dot denotes the density ${\rho }^{*}$ at the grid point set to its analytical value, and the cross and circle denote the density gradient $d{\rho }^{*}/d{x}^{*}$ at the grid point calculated by Eqs. (7) and (28), respectively. The reduced temperature is chosen as (a) $T_r=0.9$, (b) $T_r=0.8$, (c) $T_r=0.7$, (d) $T_r=0.6$, (e) $T_r=0.5$, and (f) $T_r=0.4$, respectively.

Figure 4

Coexistence curves obtained by the simulations of one-dimensional flat interface with different schemes for $\mathbf{\nabla }\rho$ and given by the Maxwell construction in thermodynamic theory. The surface tension is fixed at $\sigma =0.01$ and the interface thickness is set to (a) $W=5{\delta }_x$, (b) $W=10{\delta }_x$, (c) $W=20{\delta }_x$, and (d) $W=40{\delta }_x$, respectively.

Figure 5

Relative errors of gas density $E_{\rho ,\text{gas}}=\left({\rho }_{g,\text{numerical}}^{}-{\rho }_{g,\text{Maxwell}}^{}\right)/{\rho }_{g,\text{Maxwell}}^{}$ obtained by the simulations of one-dimensional flat interface with different schemes for $\mathbf{\nabla }\rho$. The surface tension is fixed at $\sigma =0.01$ and the interface thickness is set to (a) $W=5{\delta }_x$, (b) $W=10{\delta }_x$, (c) $W=20{\delta }_x$, and (d) $W=40{\delta }_x$, respectively.

Figure 6

Local linearization of the nonlinear EOS for multiphase fluids around the gas and liquid points. The density range with $\frac{d{p}_{\text{EOS}}^{}}{d\rho }<0$ is depicted by the gray area.

Figure 7

Variations of the gas density ${\rho }_g^{}$ with the interface curvature $1/r$ obtained by the simulations of two-dimensional circular interface with different schemes for $\mathbf{\nabla }\rho$ and analytically determined by Eq. (41). The interface thickness and surface tension are fixed at $W=10{\delta }_x$ and $\sigma =0.01$, respectively, and the reduced temperature is set to (a) $T_r=0.7$ and (b) $T_r=0.6$, respectively.

Figure 8

Evolutions of the total kinetic energy $E_k$ with the time-step number $t/{\delta }_t$ in the simulations of a moving droplet with different schemes for $\mathbf{\nabla }\rho$ and different moving velocities. The initial radius, interface thickness, and surface tension are fixed at $r_0=128{\delta }_x, W=10{\delta }_x$, and $\sigma =0.01$, respectively, and the reduced temperature is set to (a) $T_r=0.7$ and (b) $T_r=0.6$, respectively.