Multicomponent flow on curved surfaces: A vielbein lattice Boltzmann approach

We develop and implement a finite difference lattice Boltzmann scheme to study multicomponent flows on curved surfaces, coupling the continuity and Navier-Stokes equations with the Cahn-Hilliard equation to track the evolution of the binary fluid interfaces. The standard lattice Boltzmann method relies on regular Cartesian grids, which makes it generally unsuitable to study flow problems on curved surfaces. To alleviate this limitation, we use a vielbein formalism to write the Boltzmann equation on an arbitrary geometry, and solve the evolution of the fluid distribution functions using a finite difference method. Focusing on the torus geometry as an example of a curved surface, we demonstrate drift motions of fluid droplets and stripes embedded on the surface of the torus. Interestingly, they migrate in opposite directions: fluid droplets to the outer side while fluid stripes to the inner side of the torus. For the latter we demonstrate that the global minimum configuration is unique for small stripe widths, but it becomes bistable for large stripe widths. Our simulations are also in agreement with analytical predictions for the Laplace pressure of the fluid stripes, and their damped oscillatory motion as they approach equilibrium configurations, capturing the corresponding decay timescale and oscillation frequency. Finally, we simulate the coarsening dynamics of phase separating binary fluids in the hydrodynamics and diffusive regimes for tori of various shapes, and compare the results against those for a flat two-dimensional surface. Our finite difference lattice Boltzmann scheme can be extended to other surfaces and coupled to other dynamical equations, opening up a vast range of applications involving complex flows on curved geometries.

Figure 1

(a) Time evolution of stripe center ${\theta }_c$ for stripes initialized at ${\theta }_0=\pi /2$ on the torus with $r=0.8$ and $R=2$ ($a=0.4$), with initial widths of $\mathrm{\Delta }\theta =0.95\pi$ (${\theta }^{\mathrm{eq}}<\pi$), $\mathrm{\Delta }\theta =0.65\pi$ (${\theta }^{\mathrm{eq}}>\pi$), $\mathrm{\Delta }\theta =0.6\pi$ (${\theta }^{\mathrm{eq}}<\pi$, with one oscillation), and $\mathrm{\Delta }\theta =0.3\pi$ (${\theta }^{\mathrm{eq}}=\pi$). (b) Diagram indicating the location of the equilibrium position ${\theta }_c^{\mathrm{eq}}$ as a function of the stripe width $\mathrm{\Delta }{\theta }_0$ and the radii ratio $a=r/R$, for stripes initialized at ${\theta }_0=\pi /2$. (c)–(e) Examples of stripes equilibrated at (c) ${\theta }_c^{\mathrm{eq}}>\pi$ ($\mathrm{\Delta }{\theta }_0=0.65\pi$), (d) ${\theta }_c^{\mathrm{eq}}=\pi$ ($\mathrm{\Delta }{\theta }_0=0.3\pi$), and (e) ${\theta }_c^{\mathrm{eq}}<\pi$ ($\mathrm{\Delta }{\theta }_0=0.6\pi$). (f)–(h) Interface length ${\ell }_{\mathrm{total}}$ as a function of the stripe center position (solid line) for the stripe parameters considered in (c)–(e). The symbols highlight the interface lengths at maximum oscillation amplitude at initialization (0) and after each half period (1, 2, etc.).

Figure 2

Diagram indicating the location of the equilibrium position ${\theta }_c^{\mathrm{eq}}$ as a function of the stripe width $\mathrm{\Delta }{\theta }_0$ and the radii ratio $a=r/R$, for stripes initialized at ${\theta }_0=\pi /2$.

Figure 3

Comparison of the Laplace pressure obtained numerically (dashed lines and circles) against the analytic formula (3.12) for $\kappa =5\times {10}^{-4}$ (dotted lines and filled circles) and $2.5\times {10}^{-4}$ (dotted lines and empty circles). The analytic predictions (solid black lines) are almost everywhere overlapped with the numerical results.

Figure 4

Time evolution of the stripe center ${\theta }_c$ for stripes initialized at (a) ${\theta }_0=0.95\pi$ with $\mathrm{\Delta }{\theta }_0=0.280\pi$ (equilibrating at ${\theta }_c^{\mathrm{eq}}=\pi$) and (b) ${\theta }_0=0.7\pi$ with $\mathrm{\Delta }{\theta }_0=0.796\pi$ (equilibrating at ${\theta }_c^{\mathrm{eq}}=3\pi /4$). The numerical results are shown using dotted lines and symbols, while the analytic solutions are shown using solid lines.

Figure 5

The damping coefficient $\alpha$ obtained by fitting Eq. (3.13) to the simulation results (points), for stripes initialized at (a) ${\theta }_0=0.95\pi$ with ${\theta }_c^{\mathrm{eq}}=\pi$ and (b) ${\theta }_0=0.7\pi$ with ${\theta }_c^{\mathrm{eq}}=0.75\pi$. The dashed lines represent the viscous damping coefficient ${\alpha }_{\nu }$, given in Eq. (3.14).

Figure 6

The angular frequency ${\omega }_0$, obtained by fitting Eq. (3.13) to the simulation results (points). The black dashed-dotted curves correspond to the analytic expressions, as given by Eq. (3.15) for (a), when ${\theta }_c^{\mathrm{eq}}=\pi$ and Eq. (3.16) for (b), when ${\theta }_c^{\mathrm{eq}}=3\pi /4$.

Figure 7

Comparison between the values of ${\omega }_0$ obtained by fitting Eq. (3.13) to the numerical results, shown with points, and the analytic expressions (3.15) for $\mathrm{\Delta }A<\mathrm{\Delta }{A}_{\mathrm{crit}}$ and (3.16) for $\mathrm{\Delta }A>\mathrm{\Delta }{A}_{\mathrm{crit}}$, shown with solid black lines.

Figure 8

(a) Time evolution of the position of the center ${\theta }_c/\pi$ for drops initialized according to Eq. (3.17) with $({\theta }_0,R_0)\in \{(5\pi /10,0.938),(7\pi /10,0.924),(9\pi /10,0.910)\}$. (b)–(d) Snapshots of the evolution of the drop corresponding to ${\theta }_0=9\pi /10$ for $t=0$, 650, and 1775.

Figure 9

Growth of the fluid domain size $L_d\left(t\right)$ for an even mixture in two dimensions in (a) the inertial-hydrodynamics and (b) the diffusive regimes. For the diffusive regime, we remove the convective term in the Cahn-Hilliard equation. (c)–(f) Snapshots of the typical fluid configurations at $t=40$, 100, 250, and 3000 corresponding to the case indicated in (a). (g)–(j) Snapshots of the fluid configurations corresponding to the case indicated in (b), at times $t=300$, 2700, $15000$, and $168000$. These are selected such that $L_d\left(t\right)$ matches the values corresponding to (c)–(f).

Figure 10

Growth of the fluid domain size $L_d\left(t\right)$ for an even mixture on the surface of a torus with $R=2.5$ and $r=1$ in (a) the inertial-hydrodynamics and (b) the diffusive regimes. For the diffusive regime, the convective term in the Cahn-Hilliard equation is removed. (c)–(e) Snapshots of the typical fluid configurations at $t=40$, 100, and 250 corresponding to the case indicated in (a). (f)–(h) Snapshots of the typical fluid configurations at $t=350$, 4250, and $25500$ corresponding to the case indicated in (b). The times are chosen such that the values of $L_d$ match the ones corresponding to (c)–(e).

Figure 11

Growth of the fluid domain size $L_d\left(t\right)$ for an even mixture on (a) a thick torus ($R=2$, $r=1.25$) and (b) a thin torus ($R=5$, $r=0.5$). (c)–(e) and (f)–(h) Snapshots of the typical fluid configurations at $t=40$, 100, and 250 corresponding to the cases indicated in (a) and (b), respectively.

Figure 12

(a) Growth of the fluid domain $L_d\left(t\right)$ for an uneven mixture ($\overline{\varphi }=-0.3$) in two dimensions with and without hydrodynamics. In both cases, an exponent of $\frac{1}{3}$ characteristic of the diffusive regime is observed at late times. (b)–(e) Snapshots of the typical fluid configurations at times $t=50$, 250, 500, and 4000, corresponding to the case with hydrodynamics. (g)–(j) Snapshots of the fluid configurations corresponding to the case without hydrodynamics, at times $t=150$, 1050, 1500, and 6500. These are selected such that the values of $L_d\left(t\right)$ correspond to those in (b)–(e).

Figure 13

(a) Growth of the fluid domain size $L_d\left(t\right)$ for an uneven mixture on a torus with $R=2.5$ and $r=1$ with and without hydrodynamics. (b)–(d) Snapshots of the typical fluid configurations at $t=40$, 250, and 1500 corresponding to the case with hydrodynamics. (e)–(g) Snapshots of the typical fluid configurations at $t=100$, 1250, and 3500 corresponding to the case without hydrodynamics. The times are chosen such that the values of $L_d$ match those in (b)–(d).

Figure 14

(a) The average distribution of the component $(\langle \varphi \rangle +1)/2$ as a function of $\theta$ at various times. (b)–(d) Snapshots of the fluid configurations at $t=2500$, 5000, and $18000$.

Figure 15

The discrete velocity set employed by the lattice Boltzmann models based on the third (a) and fourth (b) order Gauss-Hermite quadratures. The filled black circle in the center of the figure corresponds to a lattice point in space. In our off-lattice implementation, the velocity directions do not coincide with neighboring lattice points.

Figure 16

Time evolution of the stripe and drop center ${\theta }_c$, for a stripe initialized at ${\theta }_0=0.1\pi$ and a drop initialized at ${\theta }_0=0.9\pi$, using both $Q=3$ and 4 implementations. The results overlap for all times.

Figure 17

Comparison of the growth of the fluid domain size $L_d\left(t\right)$ in the inertial-hydrodynamics regime in two dimensions for (a) an even mixture and (b) an uneven mixture, using the $Q=3$ and 4 implementations. Minor discrepancies can be observed at late times.

Figure 18

Comparison of the growth of the fluid domain size $L_d\left(t\right)$ in the inertial-hydrodynamics regime on a torus for even (a), (c), (d) and uneven (b) mixtures on the tori with $R=2.5$ and $r=1$ (a), (b), $R=2$ and $r=1.25$ (c), as well as $R=5$ and $r=0.5$ (d). The results corresponding to $Q=4$ and 3 are shown with purple $+$ signs and cyan $\times$ signs, respectively. The solid line shows the best fit curve corresponding to $a{t}^{\alpha }$, with $\alpha =\frac{2}{3}$ for (a), (c), (d) and $\frac{1}{3}$ for (b), to the $Q=4$ data on the interval $30\le t\le 200$ for the even mixture and $200\le t\le 2000$ for the uneven mixture.

Figure 19

Snapshots of the typical fluid configurations at $t=40$, 100, and 250, using the $Q=3$ (top line) and $Q=4$ (bottom line) implementations.