Improved locality of the phase-field lattice-Boltzmann model for immiscible fluids at high density ratios

Based on phase-field theory, we introduce a robust lattice-Boltzmann equation for modeling immiscible multiphase flows at large density and viscosity contrasts. Our approach is built by modifying the method proposed by Zu and He [Phys. Rev. E 87, 043301 (2013)] in such a way as to improve efficiency and numerical stability. In particular, we employ a different interface-tracking equation based on the so-called conservative phase-field model, a simplified equilibrium distribution that decouples pressure and velocity calculations, and a local scheme based on the hydrodynamic distribution functions for calculation of the stress tensor. In addition to two distribution functions for interface tracking and recovery of hydrodynamic properties, the only nonlocal variable in the proposed model is the phase field. Moreover, within our framework there is no need to use biased or mixed difference stencils for numerical stability and accuracy at high density ratios. This not only simplifies the implementation and efficiency of the model, but also leads to a model that is better suited to parallel implementation on distributed-memory machines. Several benchmark cases are considered to assess the efficacy of the proposed model, including the layered Poiseuille flow in a rectangular channel, Rayleigh-Taylor instability, and the rise of a Taylor bubble in a duct. The numerical results are in good agreement with available numerical and experimental data.

Figure 1

Effect of using different interpolation schemes for updating the relaxation time in the layered Poiseuille flow at ${\mu }^{*}=100$ (${\tau }_{\mathrm{H}}=0.5$ lu) with (a) ${\rho }^{*}=1$ and (b) ${\rho }^{*}=10$. The FD solution is shown by the solid back line, the black cross symbols represent the use of the harmonic interpolation in Eq. (22), the blue symbols ($\mathsf{I}$) represent the use of the linear interpolation in Eq. (23), and the red circles represent the use of the dynamic viscosity to update the relaxation time according to Eq. (25).

Figure 2

The behavior of the relaxation time using different interpolation schemes (${\rho }^{*}=10$ and ${\mu }^{*}=100$).

Figure 3

Normalized velocity profile for the two-phase Poiseuille flow at ${\mu }^{*}=100$ (${\tau }_{\mathrm{L}}=0.5$ lu) and density ratio of (a) 10, (b) 100, and (c) 1000. The FD solutions are shown by solid back lines, the LB results of Ref. [22] are labeled “central” and shown by green squares, the LB results of Ref. [40] are labeled “mixed” and shown by blue triangles, and the current results are labeled “current” and shown by red circles.

Figure 4

Convergence study for the layered Poiseuille flow at ${\rho }^{*}=1000$ and ${\mu }^{*}=100$ (${\tau }_{\mathrm{L}}=0.5$ lu).

Figure 5

The evolution of a single mode Rayleigh-Taylor instability at $\text{At}=0.500$ (${\rho }^{*}=3$), $\text{Re}=3000$, ${\mu }^{*}=1$, $\text{Ca}=0.26$, and $\text{Pe}=1000$.

Figure 6

Time evolution of the Rayleigh-Taylor instability at $\text{At}=0.500$ (${\rho }^{*}=3$), $\text{Re}=3000$, ${\mu }^{*}=1$, $\text{Ca}=0.26$, and $\text{Pe}=1000$ for (a) bubble front position and the (b) liquid front position. Comparative results were extracted from Refs. [18, 36, 43]. The value of $y_i$ defines the interface position at (a) $x=0$ and (b) $x=L/2$ during the simulation.

Figure 7

The evolution of a single mode Rayleigh-Taylor instability at $\text{At}=0.998$ (${\rho }^{*}=1000$), $\text{Re}=3000$, ${\mu }^{*}=100$, $\text{Ca}=0.44$, and $\text{Pe}=1000$.

Figure 8

Domain schematic of the slug flow tests for the Taylor bubble rise. The fluid domain size is $10L\times L$, and the initial bubble size is $3L\times 4L/5$.

Figure 9

The time evolution of the planar Taylor bubble with snapshots taken at $t^{*}=0,4,8,12,16,20$. The fluid properties are defined by ${\rho }^{*}=1000$ and ${\mu }^{*}=100$, while the flow condition is specified through ${\text{Re}}_r=200$ and $\text{Eo}=100$.

Figure 10

Contours of the phase field for a Taylor bubble at $t^{*}=20$ with $\text{Eo}=100$ and ${\text{Re}}_r=200$. The results from Ref. [50] were supplied by Dr. J. Fabre allowing for the current LBM outputs to be compared with the profile obtained using the boundary element method. The values $x_i$ and $y_i$ are used to define the interface location with respect to the bubble nose located at (3,0).

Figure 11

Rise of the Taylor bubble front ($\varphi ={\varphi }_0=0.5$) vs time, where $\text{Eo}=100$ and ${\text{Re}}_r=200$.

Figure 12

The strong scalability of the models on the TCLB solver on CPU.

Figure 13

The strong scalability of the models on the TCLB solver on GPU.