Generalized lattice Boltzmann method with multirange pseudopotential

The physical behavior of a class of mesoscopic models for multiphase flows is analyzed in details near interfaces. In particular, an extended pseudopotential method is developed, which permits to tune the equation of state and surface tension independently of each other. The spurious velocity contributions of this extended model are shown to vanish in the limit of high grid refinement and/or high order isotropy. Higher order schemes to implement self-consistent forcings are rigorously computed for $2d$ and $3d$ models. The extended scenario developed in this work clarifies the theoretical foundations of the Shan-Chen methodology for the lattice Boltzmann method and enhances its applicability and flexibility to the simulation of multiphase flows to density ratios up to $O\left(100\right)$.

Figure 1

Smoothing of interface properties. We show the static profiles for a flat interface with an improved SC model (14, 7). The parameters are chosen to produce a macroscopic pressure tensor (15) keeping fixed $A_1=-5.0$ (same density ratio) and at varying $A_2$: $A_2=-5.0$ $(◻)$ and $A_2=-30.0$ $(엯)$. The results of numerical simulations are compared with the analytical estimates (solid lines) obtained solving the mechanical stability equation (11) applied to (15). Notice the smoothing in the interface for a fixed density ratio due to a change in the surface tension. Inset: Laplace test for the same density ratio obtained from the numerical simulation of the improved SC model (14, 7). The parameters are chosen to reproduce a macroscopic pressure tensor (15) with $A_1=-5.0$ and different $A_2$: $A_2=-5.0$ $(◻)$, $A_2=-15.0$ $(엯)$, $A_2=-30.0$ $(▵)$. We plot the pressure drop $\mathrm{\Delta }P$ as a function of the inverse radius of the static drop. We also compare the results with the theoretical predictions (solid lines) given by the continuum analysis that leads to $\sigma (A_2=-5)=0.0398$, $\sigma (A_2=-15)=0.0716$, and $\sigma (A_2=-30)=0.100$ in lattice units. In all numerical cases the lattice Boltzmann equation (1) has been integrated in time in a fully periodic domain $L_x\times L_y=100\mathrm{\Delta }x\times 100\mathrm{\Delta }x$ and $\tau =0.7\mathrm{\Delta }t$.

Figure 2

Smoothing of interface properties. We show the static profiles for a flat interface with an improved SC model (14, 20). The parameters are chosen such as to produce a macroscopic pressure tensor (21) with $A_1=-5.0$, $A_2=-5.0∕{\rho }_0^2$, and $\lambda ={\rho }_0$ for different values of ${\rho }_0$: ${\rho }_0=1.0$ $(◻)$, ${\rho }_0=0.5$ $(엯)$. Results are also compared with the analytical estimate resulting from our continuum interpretation (solid line). The two profiles have been plotted by rescaling the lattice grid by a factor $1∕{\rho }_0$. Notice the increased interface grid resolution obtained at fixed density ratio. Inset: the results of the Laplace tests made on spherical droplets with the same density ratio and varying grid resolution is also plotted. We get the same surface tension with different interface resolutions, ${\rho }_0=1.0$ $(◻)$, ${\rho }_0=0.7$ $(엯)$, ${\rho }_0=0.5$ $(▵)$. We also compare the results with the theoretical prediction (solid line) given by the continuum analysis that leads to $\sigma (A_2=-5)=0.0398$. In all numerical cases the lattice Boltzmann equation has been integrated in time in a fully periodic domain $L_x\times L_y=100\mathrm{\Delta }x\times 100\mathrm{\Delta }x$ and $\tau =0.7\mathrm{\Delta }t$.

Figure 3

Spurious contributions around a static $2d$ drop for different density ratios with the standard SC model (6, 7). Top: we show the local Mach number defined as $\sqrt{u_x^2+u_y^2}∕c_s$ with $c_s$ the lattice sound speed. A drop of radius $15\mathrm{\Delta }x$ and density ratio $\mathrm{\Delta }\rho ∕{\rho }_g\sim 35$, obtained with ${\mathcal{G}}_b=-6.0$ in (6), is considered. Bottom: the same as in the case of the top figure with a higher density ratio $\mathrm{\Delta }\rho ∕{\rho }_g\sim 60$ obtained with ${\mathcal{G}}_b=-7.0$ in (6). Note in both plots the angular dependency due to lack of perfect radial symmetry in the forcing terms. In all numerical cases the lattice Boltzmann equation (1) has been integrated in time in a fully periodic $2d$ domain $L_x\times L_y=100\mathrm{\Delta }x\times 100\mathrm{\Delta }x$ with the drop put in the middle of the system. The relaxation time is $\tau =1.0\mathrm{\Delta }t$.

Figure 4

Spurious contributions in lattice Boltzmann and their continuum interpretation. Top: we report the structure of the velocity field around a static drop of radius $15\mathrm{\Delta }x$ and density ratio $\mathrm{\Delta }\rho ∕{\rho }_g\sim 60$. The data are the same as the bottom panel of Fig. 3. Bottom: for a qualitative comparison we have solved the Laplace equation together with the continuity equations with a matching condition in a radial ring of a given width. The velocity field is obtained in an iterative way by first solving the Poisson equation ${\nabla }^2\varphi =\nabla \bullet \mathit{u}$ and then renewing the velocity filed as $\mathit{u}\to \mathit{u}-\nabla \varphi$. For both $\mathit{u}$ and $\varphi$ periodic boundary conditions are imposed in the horizontal and vertical directions. On the ring of range $17.5\mathrm{\Delta }x<r<20.5\mathrm{\Delta }x$ the matching condition, $u_r=\cos 4\theta$, $u_{\theta }=0$ is imposed being $u_r$ and $u_{\theta }$ the radial and azimuthal components of the velocity field, respectively.

Figure 5

(Color online) The grid points identifying the set of velocity fields for a $2d$ case. With reference to the weights reported in Table , different degrees of isotropy can be achieved: fourth order [up to $w\left(2\right)$], sixth order [up to $w\left(4\right)$], eighth order [up to $w\left(8\right)$], 10th order [up to $w\left(10\right)$], 12th order [up to $w\left(17\right)$], 14th order [up to $w\left(25\right)$], and 16th order [up to $w\left(32\right)$].

Figure 6

Reduction of spurious currents with higher order isotropic tensors in the forcing terms. The spurious contributions for a static drop are analyzed for a density ratio $\mathrm{\Delta }\rho ∕{\rho }_g\sim 60$ obtained for a standard SC model (6, 7). The parameter chosen is ${\mathcal{G}}_b=-7.0$. We show the vertical velocity, $u_y$, normalized with the lattice speed of sound, $c_s$, at $x=L_x∕2$ and as a function of $y∕\mathrm{\Delta }x$. The different plots correspond to different degrees of isotropy in the forcing term (6): fourth order $(◻)$, sixth order $(엯)$, eighth order $(▵)$, and tenth order $(\times )$. Notice the reduction of the spurious contributions in the proximity of the droplet surface, $y∕\mathrm{\Delta }x\sim 30,70$. Data are obtained from lattice Boltzmann equation (1) in a fully periodic domain $L_x\times L_y=100\mathrm{\Delta }x\times 100\mathrm{\Delta }x$ and $\tau =1.0\mathrm{\Delta }t$ in lattice units with a drop in the middle of the domain.

Figure 7

Reduction of the spurious currents with grid refinement at changing ${\rho }_0$. The spurious contributions for a static drop are analyzed for a density ratio $\mathrm{\Delta }\rho ∕{\rho }_g\sim 60$ obtained for the improved SC model (14, 20). The parameters are chosen in such a way to reproduce (15) with $A_1=-7.0$. We show the vertical velocity, $u_y$, normalized with the lattice speed of sound, $c_s$, at $x=L_x∕2$ and as a function of $y∕\mathrm{\Delta }x$. The different plots correspond to different degrees of refinement obtained with $A_2=-7.0∕{\rho }_0^2$ and $\lambda ={\rho }_0$ in (21): ${\rho }_0=1.0$ $(◻)$, ${\rho }_0=0.75$ $(엯)$, and ${\rho }_0=0.5$ $(▵)$. Details of the numerical simulations are the same as Fig. 6.

Figure 8

Reduction of the spurious currents with refinement in the forcing terms. The spurious contributions for a static drop are analyzed for various density ratios obtained in the improved SC model (14, 20). The parameters are chosen in such a way as to reproduce (15) with different $A_1$, thus changing the density ratio $\mathrm{\Delta }\rho ∕{\rho }_g$. We show the maximum velocity due to spurious contributions, normalized with the lattice speed of sound $\left(c_s\right)$ as a function of $\mathrm{\Delta }\rho ∕{\rho }_g$. The two different plots correspond to different degree of refinement obtained with $A_2=7.0∕{\rho }_0^2$ and $\lambda ={\rho }_0$ in (21) and different values of ${\rho }_0$: ${\rho }_0=1.0$ ($◻$, standard SC model), ${\rho }_0=0.5$ ($엯$, improved version). Notice the net reduction of the spurious contributions at fixed macroscopic physics (same density ratio and surface tension) obtained through the rescaling $A_2\sim 1∕{\rho }_0^2$. The numerical details are the same as described in Fig. 6.