Optimal surface-tension isotropy in the Rothman-Keller color-gradient lattice Boltzmann method for multiphase flow

The Rothman-Keller color-gradient (CG) lattice Boltzmann method is a popular method to simulate two-phase flow because of its ability to deal with fluids with large viscosity contrasts and a wide range of interfacial tensions. Two fluids are labeled red and blue, and the gradient in the color difference is used to compute the effect of interfacial tension. It is well known that finite-difference errors in the color-gradient calculation lead to anisotropy of interfacial tension and errors such as spurious currents. Here, we investigate the accuracy of the CG calculation for interfaces between fluids with several radii of curvature and find that the standard CG calculations lead to significant inaccuracy. Specifically, we observe significant anisotropy of the color gradient of order $7\%$ for high curvature of an interface such as when a pinchout occurs. We derive a second order accurate color gradient and find that the diagonal nearest neighbors can be weighted differently than in the usual color-gradient calculation such that anisotropy is minimized to a fraction of a percent. The optimal weights that minimize anisotropy for the smallest radius of curvature interface are found to be $w=(0.298,0.284,0.275)$ for diagonal nearest neighbors for the cases of the interface smoothing parameter $\beta =(0.5,0.7,0.99)$, somewhat higher than the $w=0.25$ value derived by Leclaire et al. [Leclaire, Reggio, and Trepanier, Computers and Fluids 48, 98 (2011)] based on obtaining isotropic errors to second order. We find that use of these optimal $w$ values yields over a factor of 10 decrease in anisotropy and over a factor of 30 decrease in mean anisotropy relative to using the standard $w=1$ value. And we find a factor of about 2 decrease in the anisotropic error and up to factor 15 decrease in mean anisotropic error relative to the choice of $w=0.25$ for small radius of curvature interfaces. The improved CG calculations will allow the method to be more reliably applied to studies of phenomenology and pore scale processes such as viscous and capillary fingering, and droplet formation where surface-tension isotropy of narrow fingers and small droplets plays a crucial role in correctly capturing phenomenology. We present an example illustrating how different phenomena can be captured using the improved color-gradient method. Namely, we present simulations of a wetting fluid invading a fluid filled pipe where the viscosity ratio of fluids is unity in which droplets form at the transition to fingering using the improved CG calculations that are not captured using the standard CG calculations. We present an explanation of why this is so which relates to anisotropy of the surface tension, which inhibits the pinchouts needed to form droplets.

Figure 1

The $D2Q9$ lattice.

Figure 2

Snapshots of color difference $D=({\rho }_r-{\rho }_b)$ for small droplets of red fluid with three different radii of 2, 4, and 8 from left to right.

Figure 3

Plots of color difference $D(x,z_0)=({\rho }_r-{\rho }_b)(x,z_0)$ and cubic spline fits of the droplets with three different radii of $r_0=2,4$, and 8 from left to right.

Figure 4

Plots of the normalized color gradient calculated using Eqs. (37, 38, 39) for three different radii of droplets $r_0$ and three different values of $w$. From left to right: Color gradient (CG) using $w=1$ (original CG calculation), calculation using $w=0.25$ (diagonals weighted down by a factor of $w=0.25$), and calculation using $w=0$ (CG via finite differences along coordinate axes). From top to bottom: Smallest radius ($r_0=2$), middle radius ($r_0=4$), and largest radius ($r_0=8$).

Figure 5

Plots of the anisotropic error $E\left(w\right)$ in the space domain calculated using different values of $w$ calculated using Eq. (42). From left to right: Error for $w=1$ (original CG calculation), calculation using $w=0.25$ (diagonals weighted down by a factor of $w=0.25$), and calculation using $w=0$ (CG via Cartesian finite differences). From top to bottom: Smallest radius ($r_0=2$), middle radius ($r_0=4$), and largest radius ($r_0=8$).

Figure 6

Plots of the maximum anisotropy vs $w$ for three different radii $r_0$ of a droplet.

Figure 7

Plots of the normalized color gradient as a function of direction $\theta$ at the radius of maximal color gradient $r_F={max}_r{F}_{r,0^{\circ }}\left(w\right)$ for the original color-gradient calculation ($w=1$) and the proposed optimal value of $w=0.298$ for the smallest droplet with $r_0=2$.

Figure 8

Plots of the normalized color gradient as a function of direction $\theta$ at the radius of maximal anisotropy $r^{*}={max}_r{E}_{r,{45}^{\circ }}-E_{r,0^{\circ }}$ for the original color-gradient calculation ($w=1$) and the proposed optimal value of $w=0.298$ for the smallest droplet with $r_0=2$.

Figure 9

Plots of the anisotropic error $E\left(w\right)$ in the space domain calculated using the proposed optimal value of $w=0.298$. From left to right: Error for $r_0=2$, 4, and 8.

Figure 10

Plots of the mean anisotropic error as a function of direction $\theta$ at the radius of maximal anisotropy $r^{*}={max}_r(E_{r,{45}^{\circ }}-E_{r,0^{\circ }})$ for the original color-gradient calculation ($w=1$) and the proposed optimal value of $w=0.298$ for the smallest droplet with $r_0=2$.

Figure 11

Plot of the anisotropic error at ${45}^{\circ }$ together with the color gradient for the two cases of $w=1$ (original CG calculation) and $w=0.298$ (optimal $w$ that minimizes anisotropy) for the smallest droplet with $r_0=2$. The color gradient is scaled by factors of $s_1$ and $s_2$ to enable the shape of the color-gradient magnitude $\left|\mathbf{F}\right(w\left)\right|$ to be visualized on the same plot as the error.

Figure 12

Plots for the case of $w=0.25$ for the smallest radius of droplet with $r_0=2$ for comparison with Figs. 7, 8, 10, and 11 which were for the optimal case of $w=0.298$. Upper left: The normalized color gradient as a function of direction $\theta$ at the radius of maximal color gradient $r_F={max}_r{F}_{r,0^{\circ }}\left(w\right)$ for the original color-gradient calculation ($w=1$) and the value of $w=0.25$. Upper right: The normalized color gradient as a function of direction $\theta$ at the radius of maximal anisotropy $r^{*}={max}_r{E}_{r,{45}^{\circ }}-E_{r,0^{\circ }}$ for the original color-gradient calculation ($w=1$) and the value of $w=0.25$. Lower left: The mean anisotropic error as a function of direction $\theta$ for the original color-gradient calculation ($w=1$) and the value of $w=0.25$. Lower right: The anisotropic error at ${45}^{\circ }$ together with the color gradient for the two cases of $w=1$ (original CG calculation) and $w=0.25$. The color gradient is scaled by factors of $s_1$ and $s_2$ to enable the shape of the color-gradient magnitude $\left|\mathbf{F}\right(w\left)\right|$ to be visualized on the same plot as the error.

Figure 13

Plots of the anisotropic error $E\left(w\right)$ in the space domain calculated using the value of $w=0.25$ for comparison with Fig. 9 calculated with the optimal value of $w=0.298$. From left to right: Error for $r_0=2$, 4, and 8.

Figure 14

Snapshots of flow in a pipe for wetting angle of $\theta =0$ and various capillary numbers using the optimal value of $w=0.298$ to calculate the color gradient.

Figure 15

Snapshots of flow in a pipe for wetting angle of $\theta =0$ and various capillary numbers using the standard color-gradient calculations corresponding to a value of $w=1$.

Figure 16

Mechanism explaining how the anisotropy in the standard color-gradient calculation with $w=1$ resists a finger head from breaking off and forming a droplet.

Figure 17

Snapshots of flow in a pipe for wetting angle of $\theta =0$ and various capillary numbers using Leclaire's value of $w=0.25$ to calculate the color gradient.