Chromodynamic multirelaxation-time lattice Boltzmann scheme for fluids with density difference

We develop, after Dellar [Phys. Rev. E. 65, 036309 (2002); J. Comput. Phys. 190, 351 (2003)], a multiple-relaxation-time (MRT), chromodynamic, multicomponent lattice Boltzmann equation (MCLBE) scheme for simulation of isothermal, immiscible fluid flow with a density contrast. It is based on Lishchuk's method [Brackbill, Kothe, and Zemach, J. Comp. Phys. 100, 335 (1992); Lishchuk, Care, and Halliday, Phys. Rev. E. 67, 036701, (2003)] and the segregation of d'Ortona et al. [Phys. Rev. E. 51, 3718, (1995)]. We focus on fundamental model verifiability but do relate some of our data to that from previous approaches, due to Ba et al. [Phys. Rev. E 94, 023310 (2016)] and earlier Liu et al. [Phys. Rev. E 85, 046309 (2012)], who pioneered large density difference chromodynamic MCLBE and showed the practical benefits of an MRT collision model. Specifically, we test the extent to which chromodynamic MCLBE MRT schemes comply with the kinematic condition of mutual impenetrability and the continuous traction condition by developing analytical benchmarking flows. We conclude that our data, taken with those of Ba et al., verify the utility of MRT chromodynamic MCLBE.

Figure 1

Schematic. Square $D2Q9$ lattice with our indexing convention. Odd values of $i$ identify the longer links.

Figure 2

Velocity variation for the test illustrated in Fig. 9. Simulation data are represented by crosses and semianalytic theory (which accounts for the transverse variation of the density) by the continuous line. For these data, $L_x=200, {\alpha }_B=0.2, {\alpha }_R=0.9, \mathrm{\Lambda }=8, {\nu }_R={\nu }_B=0.333, \beta =0.4$.

Figure 3

Flow at a flat, sheared interface; comparison of simulation data (crosses) and a semianalytical solution (see Appendix pp3), for a large range of density contrasts, $\mathrm{\Lambda }$. For these data, shear viscosity $\eta =0.166$, while kinematic viscosity $\nu =\frac{\eta }{\rho }$ changes. The interface centers on $x=500$ lattice units, with the fluid on the right the denser. The fluid on the left (right) is moving south (north) initially. For panels (a)–(d), $\mathrm{\Lambda }=1,20,31.25,50$, respectively. These data confirm continuity of velocity and correct transmission of stress across a flat, sheared interface.

Figure 4

Normalized microcurrent flow excerpt for $\mathrm{\Lambda }=10$ in the vicinity of a drop, radius $R=20$ for (a) constant curvature $K=\frac{1}{R}$ and (b) numerically calculated curvature $K={\mathbf{\nabla }}_s{\rho }^N$. See Tables 3 and 4 to scale these velocity fields. Note, the circulation in the case of fixed curvature [panel (a)] is more localized.

Figure 5

Low Re internal flow past a tethered, cylindrical drop for $\mathrm{\Lambda }=5$. Flow outside the drop has been suppressed in these data. Panel (a) shows the total flow, in which the velocity field clearly has a nonphysical component perpendicular to the interface. Panel (b) shows the microcurrent error, measured from the frozen phase field from panel (a), without external flow, and panel (c) shows the physical flow exposed by subtracting the microcurrent in panel (b). The solid black line represents the center of the interface between the fluids (${\rho }^N=0$ contour). The internal flows are clearly parallel to the interface, as is the external flow (not shown).

Figure 6

Flow at a curved, sheared interface; comparison of simulation data (crosses) and semianalytical solution (see Appendix pp3) for a range of density contrasts, $\mathrm{\Lambda }$. For these data, shear viscosity $\eta =$ constant, while kinematic viscosity $\nu =\frac{\eta }{\rho }$ changes. In all panels the interface centers on $r=120$ lattice units, with fluid on the left the denser. The fluid on the left (right) is in uniform rotation (at rest), initially. For panels (a)–(d), $\mathrm{\Lambda }=1,2,3,5$, respectively, which is smaller than the range of $\mathrm{\Lambda }$ shown in Fig. 3. These data confirm continuity of velocity and correct transmission of stress across a curved interface.

Figure 7

Schematic. Geometry and initial conditions. Fluid is in unidirectional flow, $u\left(x\right){\widehat{e}}_y$. There is translational invariance in the $y$ direction and density stratification, with the red (blue) fluid assumed to have density ${\rho }_1$ (${\rho }_2$). Initially, the red (blue) fluid moves in the $y$ ($-y$) direction.

Figure 8

Schematic. Geometry and initial conditions. A weakly compressible fluid is in rotational flow $u_{\varphi }\left(r\right){\widehat{e}}_{\varphi }$. $r$ represents a “transverse” coordinate. A no-slip boundary is located at $r=\left[R\right]$. There is transverse density stratification with the red (blue) fluid having density ${\rho }_1$ (${\rho }_2$). Initially, the red (blue) fluid moves in the ${\widehat{e}}_{\varphi }$ (is at rest).

Figure 9

Schematic representation of the geometry and transverse density stratification used in the pressure gradient (white arrow) driven flow tests.