Multicomponent Flow on Curved Surfaces: A Vielbein Lattice Boltzmann Approach

TL;DR

This paper introduces a novel lattice Boltzmann method using vielbein formalism to simulate multicomponent fluid flows on curved surfaces like tori, capturing complex interface dynamics and phase separation behaviors.

Contribution

The authors develop a new lattice Boltzmann scheme for curved geometries using vielbein formalism, enabling accurate simulation of multicomponent flows on arbitrary surfaces.

Findings

Fluid droplets and stripes migrate in opposite directions on a torus.

The global minimum configuration for stripes is unique at small widths but bistable at larger widths.

Simulations agree with analytical predictions for Laplace pressure and oscillatory motion.

Abstract

We develop and implement a novel 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. 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 down the Boltzmann equation on an arbitrary geometry, and solve the evolution of the fluid distribution functions using a finite difference method. Focussing 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 a 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 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.

Equations

$$\mathbf{u}(\mathbf{x})=u^a(q^b){\partial }_{\mathbf{a}},$$

$${\mathbf{u}}^2=g_{ab}{u}^a{u}^b.$$

$$⟨{\mathbf{\omega }}^{\hat{\mathbf{a}}},{\mathbf{e}}_{\hat{\mathbf{b}}}⟩\equiv {\omega }_a^{\hat{a}}{e}_{\hat{b}}^a={\delta }^{\hat{a}}{}_{\hat{b}},{\omega }_a^{\hat{a}}{e}_{\hat{a}}^b={\delta }^b{}_a,g_{ab}{e}_{\hat{a}}^a{e}_{\hat{b}}^b={\delta }_{\hat{a}\hat{b}}.$$

$$\mathbf{u}=u^{\hat{a}}{\mathbf{e}}_{\hat{\mathbf{a}}},$$

$$u^{\hat{a}}={\omega }_a^{\hat{a}}{u}^a,u^a=e_{\hat{a}}^a{u}^{\hat{a}},$$

$${\mathbf{u}}^2={\delta }_{\hat{a}\hat{b}}{u}^{\hat{a}}{u}^{\hat{b}}.$$

$$c_{\hat{a}\hat{b}}{}^{\hat{c}}=⟨{\mathbf{\omega }}^{\hat{\mathbf{c}}},[{\mathbf{e}}_{\hat{\mathbf{a}}},{\mathbf{e}}_{\hat{\mathbf{b}}}]⟩={\omega }_a^{\hat{c}}([e_{\hat{a}},e_{\hat{b}}]{)}^a,$$

$${\Gamma }^{\hat{d}}{}_{\hat{b}\hat{c}}={\delta }^{\hat{d}\hat{a}}{\Gamma }_{\hat{a}\hat{b}\hat{c}},{\Gamma }_{\hat{a}\hat{b}\hat{c}}=\frac{1}{2}(c_{\hat{a}\hat{b}\hat{c}}+c_{\hat{a}\hat{c}\hat{b}}-c_{\hat{b}\hat{c}\hat{a}}).$$

$$\Psi ={\int }_V\left[\frac{A}{4}(1-{\varphi }^2{)}^2+\frac{\kappa }{2}(\nabla \varphi {)}^2\right]dV,$$

$${\xi }_0=\sqrt{\frac{\kappa }{A}},\gamma =\sqrt{\frac{8\kappa A}{9}}.$$

$$\mu (\mathbf{x})=\frac{\delta \Psi }{\delta \varphi (\mathbf{x})}=-A\varphi (1-{\varphi }^2)-\kappa \Delta \varphi .$$

$${\partial }_t\varphi +{\nabla }_{\hat{a}}(u^{\hat{a}}\varphi )={\nabla }_{\hat{a}}(M{\nabla }^{\hat{a}}\mu ),$$

$${\partial }_{t}n+{\nabla }_{\hat{a}}(u^{\hat{a}}n)=0,nm\frac{D{u}^{\hat{a}}}{Dt}=-{\nabla }_{\hat{b}}{T}^{\hat{a}\hat{b}}+n{F}^{\hat{a}},$$

$$n{F}^{\hat{a}}=-\varphi {\nabla }^{\hat{a}}\mu =-{\nabla }^{\hat{a}}{p}_{\mathrm{binary}}+\kappa \varphi {\nabla }^{\hat{a}}\Delta \varphi ,$$

$$p_{\mathrm{binary}}=A\left(-\frac{1}{2}{\varphi }^2+\frac{3}{4}{\varphi }^4\right).$$

$$\frac{\partial f}{\partial t}+\frac{1}{\sqrt{g}}\frac{\partial }{\partial q^b}\left(v^{\hat{a}}{e}_{\hat{a}}^{b}f\sqrt{g}\right)+\frac{\partial }{\partial v^{\hat{a}}}\left[\left(\frac{F^{\hat{a}}}{m}-{\Gamma }^{\hat{a}}{}_{\hat{b}\hat{c}}{v}^{\hat{b}}{v}^{\hat{c}}\right)f\right]=J[f],$$

$$\begin{gathered} \frac{\partial f}{\partial t}+\frac{v^{\hat{\phi }}}{R+r\cos \theta }\frac{\partial f}{\partial \phi }+\frac{v^{\hat{\theta }}}{r(1+a\cos \theta )}\frac{\partial [f(1+a\cos \theta )]}{\partial \theta }+\frac{F^{\hat{\phi }}}{m}\frac{\partial f}{\partial v^{\hat{\phi }}}+\frac{F^{\hat{\theta }}}{m}\frac{\partial f}{\partial v^{\hat{\theta }}} \\ -\frac{\sin \theta }{R+r\cos \theta }\left[v^{\hat{\phi }}\frac{\partial (f{v}^{\hat{\phi }})}{\partial v^{\hat{\theta }}}-v^{\hat{\theta }}\frac{\partial (f{v}^{\hat{\phi }})}{\partial v^{\hat{\phi }}}\right]=-\frac{1}{\tau }[f-f^{(\mathrm{eq})}]. \end{gathered}$$

$$J[f]=-\frac{1}{\tau }[f-f^{\mathrm{eq}}].$$

$$\nu =\frac{\eta }{nm}=\frac{{\eta }_v}{nm}=\frac{\tau k_{B}T}{m},$$

$$\left(\begin{array}{c}{v}_1\\ v_2\\ v_3\\ v_4\end{array}\right)=\left(\begin{array}{c}-\sqrt{3+\sqrt{6}}\\ -\sqrt{3-\sqrt{6}}\\ \sqrt{3-\sqrt{6}}\\ \sqrt{3+\sqrt{6}}\end{array}\right).$$

$$n=\sum _{\mathbf{k}}{f}_{\mathbf{k}},n\mathbf{u}=\sum _{\mathbf{k}}{f}_{\mathbf{k}}{\mathbf{v}}_{\mathbf{k}}.$$

$$f_{\mathbf{k}}^{\mathrm{eq}}=n{w}_{k_{\theta }}{w}_{k_{\phi }}\left\{1+{\mathbf{v}}_{\mathbf{k}}\cdot \mathbf{u}+\frac{1}{2}[({\mathbf{v}}_{\mathbf{k}}\cdot \mathbf{u}{)}^2-{\mathbf{u}}^2]+\frac{1}{6}{\mathbf{v}}_{\mathbf{k}}\cdot \mathbf{u}[({\mathbf{v}}_{\mathbf{k}}\cdot \mathbf{u}{)}^2-3{\mathbf{u}}^2]\right\}.$$

$$w_k=\frac{Q!}{H_{Q+1}^2(v_k)},$$

$$w_1=w_4=\frac{3-\sqrt{6}}{12},w_2=w_3=\frac{3+\sqrt{6}}{12}.$$

$${\left(\frac{\partial f}{\partial v^{\hat{\theta }}}\right)}_{k_{\theta }{k}_{\phi }}=\sum _{k_{\theta }^{\prime }=1}^Q{\mathcal{K}}_{k_{\theta },k_{\theta }^{\prime }}^H{f}_{k_{\theta }^{\prime },k_{\phi }},{\left(\frac{\partial f}{\partial v^{\hat{\phi }}}\right)}_{k_{\theta }{k}_{\phi }}=\sum _{k_{\phi }^{\prime }=1}^Q{\mathcal{K}}_{k_{\phi },k_{\phi }^{\prime }}^H{f}_{k_{\theta },k_{\phi }^{\prime }},$$

$${\left(\frac{\partial (f{v}^{\hat{\phi }})}{\partial v^{\hat{\phi }}}\right)}_{k_{\theta }{k}_{\phi }}=\sum _{k_{\phi }^{\prime }=1}^Q{\tilde{\mathcal{K}}}_{k_{\phi },k_{\phi }^{\prime }}^H{f}_{k_{\theta },k_{\phi }^{\prime }},$$

$${\mathcal{K}}_{i,m}^H=$$

$${\tilde{\mathcal{K}}}_{i,m}^H=$$

$${\mathcal{K}}_{i,m}^H=$$

$${\tilde{\mathcal{K}}}_{i,m}^H=$$

$${\theta }_{-}={\theta }_c-\Delta \theta /2,{\theta }_{+}={\theta }_c+\Delta \theta /2.$$

$$\Delta A=2\pi rR{\int }_{{\theta }_{-}}^{{\theta }_{+}}d\theta (1+a\cos \theta )=2\pi rR[\Delta \theta +2a\sin (\Delta \theta /2)\cos {\theta }_c],$$

Full text

I Introduction

Hydrodynamics on curved manifolds is relevant for a wide range of physical phenomena. Examples range from the motion of electrons in graphene at the micro-scale (Giordanelli et al., 2018), through thin liquid films (Schwartz & Weidner, 1995; Howell, 2003), confined active matter (Keber et al., 2014; Henkes et al., 2018; Janssen et al., 2017) and bio-membranes (Henle & Levine, 2010; Arroyo & Desimone, 2009) at the meso-scale, to relativistic flows in astrophysics (Martí et al., 2015) and at the cosmological scale (Ellis et al., 2012). However, despite its importance, the study of flows on curved space has received much less attention when compared to corresponding investigations on two- and three-dimensional flat space. Suitable numerical approaches to study these problems are also still limited, especially when the flow phenomena of interest involve several fluid components.

Here our focus is on multicomponent flow on curved two-dimensional surfaces. An important motivation to study such problem arises from biological membranes and their synthetic counterparts. Experimentally it has been observed that self-assembled lipid and polymer membranes can adopt an astonishing range of shapes and morphologies (Seddon & Templer, 1995), from single bilayers to stacks and convoluted periodic structures. Moreover, these membranes are usually comprised of several species, which can mix or demix depending on the thermodynamic conditions under which they are prepared (Baumgart et al., 2003; Bacia et al., 2005; Aufderhorst-Roberts et al., 2017). The interplay between curvature and composition is a ubiquitous structural feature for bio-membranes, and they are key to biological functions and synthetic membrane-based applications (McMahon & Gallop, 2005; Parthasarathy et al., 2006; Pontani et al., 2013).

There is much interest to understand this interplay between membrane curvature and composition. However, to date continuum modelling of membranes with several lipid components have largely focussed on their equilibrium configurations (Jülicher & Lipowsky, 1996; Hu et al., 2011; Paillusson et al., 2016; Fonda et al., 2018). Several dynamic studies of phase separation on curved surfaces have been carried out in the literature. However, apart from a few exceptions (Nitschke et al., 2012), they usually involve diffusive dynamics and ignore the importance of hydrodynamics (Marenduzzo & Orlandini, 2013; Jeong & Kim, 2015; Gera & Salac, 2017). The aim of this paper is to develop a flexible lattice Boltzmann framework to simulate multicomponent flow on arbitrary curved surfaces. For simplicity, here we will assume the two-dimensional flow is Newtonian. For lipid membranes, this assumption is supported by both Molecular Dynamics simulations and experimental observations (Dimova et al., 2006; den Otter & Shkulipa, 2007; Cicuta et al., 2007; Camley & Bown, 2011).

Our approach is based on the lattice Boltzmann method (LBM) (Krüger et al., 2017; Succi, 2001), which has recently become increasingly popular to study multicomponent flow phenomena, with good agreement against experiments and other simulation methods, including for drop dynamics, liquid phase separation, microfluidics and porous media (Liu et al., 2016; Sadullah et al., 2018; Varagnolo et al., 2013; Liu et al., 2015). Within the lattice Boltzmann literature, there are several models for multicomponent flow, including the so-called free energy (Swift et al., 1996), pseudo-potential (Shan & Chen, 1993) and color (Gunstensen et al., 1991) models. In this work, we have chosen to employ the free energy model, though our framework can be adapted to account for the pseudo-potential and color models. Our approach can also be extended to account for more fluid components (Semprebon et al., 2016; Wörhwag et al., 2018; Abadi et al., 2018; Liang et al., 2016; Ridl & Wagner, 2018), as well as coupled to other dynamical equations, including those for liquid crystals (Spencer & Care, 2006; Denniston et al., 2001) and viscoelastic fluids (Malaspinas et al., 2010; Gupta et al., 2015).

Standard lattice Boltzmann method is based on regular Cartesian grids, and this makes it unsuitable for use on curved surfaces. Thus, an important contribution of this work is the use of vielbein formalism (Busuioc & Ambru\cbs, 2019) to solve the hydrodynamics equations of motion on curved surfaces, which we combine with the free energy binary fluid model. To our knowledge, this is the first time this has been done. Previous lattice Boltzmann simulations on curved surfaces have been carried out for single-component flows (Mendoza et al., 2013; Busuioc & Ambru\cbs, 2019; Hejranfar et al., 2017). Since it is not always possible to have a regular grid on curved surfaces and for the discrete velocity sets to coincide with the lattice discretizations, here we solve the discrete Boltzmann equation using a finite difference approach, rather than a collision-propagation scheme. The latter is usually the case in standard lattice Boltzmann implementation.

The capabilities of our new method are demonstrated using several problems. Firstly, we study drift motion of fluid droplets and stripes when placed on the surface of a torus. This drift is due to non-uniform curvature, and as such, is not present on flat space, or for surfaces with uniform curvature (e.g. a sphere). For the stripes, analytical results are available for their equilibrium configuration, Laplace pressure, and relaxation dynamics (Busuioc et al., 2019b), thus providing an excellent platform to systematically examine the accuracy of our method. We demonstrate that these predictions are accurately captured in our simulations. Secondly, we simulate binary phase separation on the surface of a torus for equal and unequal compositions, both in diffusive and hydrodynamic regimes. We compare and contrast the results for tori of various shapes against those for flat two-dimensional surface (Bray, 2002; Kendon et al., 2001; Wagner & Yeomans, 1997).

II Computational Model and Method

In this section we develop a framework that allows simulations of multicomponent flow on arbitrary curved surfaces. Our vielbein lattice Boltzmann approach has three key features. Firstly, similar to standard lattice Boltzmann method, we exploit the Boltzmann equation to solve the continuum equations of motion, and we use a discrete and finite set of fluid distribution functions. Secondly, unlike standard lattice Boltzmann method, the discrete velocity sets do not coincide with the neighbouring lattice points. Thus, rather than solving the Boltzmann equation using a sequence of collision and propagation steps, we take advantage of a finite difference method. Thirdly, to describe the curved surface, we employ a vielbein field, which decouples the velocity space from the coordinate space (Cardall et al., 2013; Busuioc & Ambru\cbs, 2019). This simplifies the formulation and computation of the governing Boltzmann equation.

II.1 Brief Introduction to Vielbein Fields

Let us begin by considering a two-dimensional curved surface embedded in three dimensions. Vector fields, such as the velocity field $\bm{u}(\bm{x})$, on the two-dimensional surface can be expressed in the curvilinear coordinate system using

[TABLE]

where $u^{a}(q^{b})$ represent the components of the velocity field on a manifold parametrised using the coordinates $q^{b}$ ($1\leq a,b\leq 2$ for two-dimensional manifolds). Furthermore, the squared norm of the velocity field $\bm{u}$ can be computed as

[TABLE]

$g_{ab}$ is called the metric tensor. This description of vector fields in curvilinear coordinates can become inconvenient for practical computations. This is because the elements of the metric tensor $g_{ab}$ may become singular at various points due to the choice of surface parametrisation. In such instances, the contravariant components $u^{a}$ of the velocity must diverge in order for squared norm $\bm{u}^{2}$ to remain finite.

The difficulty described above can be alleviated by introducing, as an interface between the coordinate space and the velocity space, the vielbein vector fields (frame) $\bm{e_{{\hat{a}}}}=e_{{\hat{a}}}^{a}\bm{\partial_{a}}$. Dual to the vielbein vector fields are the vielbein one-forms (co-frame) $\bm{\omega^{\hat{a}}}=\omega^{\hat{a}}_{a}\bm{dq^{a}}$. We reserve the hatted indices to denote the vielbein framework. The vielbein frame and co-frame have to satisfy the following relations

[TABLE]

With the above vielbein frame and co-frame, the vector field $\bm{u}$ can be written as

[TABLE]

where the vector field components are

[TABLE]

and the squared norm

[TABLE]

In the vielbein framework, the information on the metric tensor is effectively absorbed in the components of the vector field, which makes the formulation and derivation of the lattice Boltzmann approach significantly less cumbersome.

In the lattice Boltzmann implementation used in this paper, we need to introduce two more geometrical objects. First, the Cartan coefficients $c_{{\hat{a}}{\hat{b}}}{}^{\hat{c}}$ are defined as

[TABLE]

with the commutator $([\bm{e_{\hat{a}}},\bm{e_{\hat{b}}}])^{a}=e_{\hat{a}}^{b}\partial_{b}e_{\hat{b}}^{a}-e_{\hat{b}}^{b}\partial_{b}e_{\hat{a}}^{a}$. Second, $\Gamma^{\hat{a}}{}_{{\hat{b}}{\hat{c}}}$ and $\Gamma_{{\hat{a}}{\hat{b}}{\hat{c}}}$ represent the connection coefficients, which are defined as

[TABLE]

In Appendix A, we detail the application of the vielbein formalism for a torus. It is worth noting that our approach is general and other curved geometries can be handled in a similar way.

II.2 Binary Fluid Model and Equations of Motion

We consider a binary mixture of fluids $A$ and $B$, characterised by an order parameter $\phi$, such that $\phi=1$ corresponds to a bulk $A$ fluid and $\phi=-1$ to a bulk $B$ fluid. A simple free energy model that allows the coexistence of these two bulk fluids is given by the following Landau free energy (Briant & Yeomans, 2004; Krüger et al., 2017)

[TABLE]

where $A$ and $\kappa$ are free parameters, which are related to the interface width $\xi_{0}$ and surface tension $\gamma$ through

[TABLE]

The chemical potential can be derived by taking the functional derivative of the free energy with respect to the order parameter, giving

[TABLE]

The evolution of the order parameter $\phi$ is specified by the Cahn-Hilliard equation. In covariant form it is given by

[TABLE]

where the hatted indices are taken with respect to the orthonormal vielbein basis. Equivalently, indices with respect to the coordinate basis can be used, e.g. $\nabla_{{\hat{a}}}(u^{\hat{a}}\phi)=\nabla_{a}(u^{a}\phi)$. In the above, $M$ is the mobility parameter, $\mu$ is the chemical potential, and the fluid velocity $\bm{u}$ is a solution of the continuity and Navier-Stokes equations

[TABLE]

where $D/Dt=\partial_{t}+u^{\hat{b}}\nabla_{\hat{b}}$ is the material (convective) derivative, $m$ is the particle mass, $n$ is the number density and $T^{{\hat{a}}{\hat{b}}}=p_{\rm i}\delta^{{\hat{a}}{\hat{b}}}+\sigma^{{\hat{a}}{\hat{b}}}$ is the ideal gas stress tensor. $T^{{\hat{a}}{\hat{b}}}$ consists of the ideal gas pressure $p_{\rm i}=nk_{B}T$ and the viscous stress for the Newtonian fluid $\sigma^{{\hat{a}}{\hat{b}}}=-\eta(\nabla^{\hat{a}}u^{\hat{b}}+\nabla^{\hat{b}}u^{\hat{a}}-\delta^{{\hat{a}}{\hat{b}}}\nabla_{\hat{c}}u^{\hat{c}})-\eta_{v}\delta^{{\hat{a}}{\hat{b}}}\nabla_{\hat{c}}u^{\hat{c}}$. The latter is written in terms of the dynamic (shear) and volumetric (bulk) viscosities $\eta$ and $\eta_{v}$. The thermodynamic force term $F^{\hat{a}}$ takes the following form

[TABLE]

A summary on how the differential operators must be applied for the cases of the Cartesian and torus geometries is provided in the Supplementary Information (Ambru\cbs et al., 2019).

II.3 The Vielbein Lattice Boltzmann Approach

In this paper, we employ the lattice Boltzmann approach to solve the hydrodynamics equations [Eq. (13)], while the Cahn-Hilliard equation [Eq. (12)] is solved directly using a finite difference method. The details of the numerical implementation are discussed in the Supplementary Information (Ambru\cbs et al., 2019). It is possible to solve the Cahn-Hilliard equation using a lattice Boltzmann scheme, and on flat manifolds, it has been suggested that extension to more fluid components is more straightforward in this approach (Li & Wagner, 2007; Ridl & Wagner, 2018). However, for our purpose here, it is more expensive and require us to use a higher order quadrature.

We use a discretised form of the Boltzmann equation that reproduces the fluid equations of motion in the continuum limit. In covariant form, the Boltzmann equation on an arbitrary geometry is given by (Busuioc & Ambru\cbs, 2019):

[TABLE]

where $\sqrt{g}$ is the square root of the determinant of the metric tensor, and $J[f]$ is the collision operator.

For the specific case of a torus, the Boltzmann equation reads

[TABLE]

The steps needed to derive Eq. (16) from Eq. (15) are summarised in Appendix A. Here $r$ and $R$ represent the inner (small) and outer (large) radii, $a=r/R$ is the radii ratio, while the angle $\theta$ goes round the inner circle and $\varphi$ covers the large circle. The range for both $\theta$ and $\varphi$ is $[0,2\pi)$ and the system is periodic with respect to both these angles. The last term on the left hand side of Eq. (16) corresponds to inertial and reaction forces that arise when we have flow on curved surfaces, since fluid motion is constrained on the surface.

As commonly the case in the lattice Boltzmann literature, we employ the BGK approximation for the collision operator,

[TABLE]

The relaxation time $\tau$ is related to the fluid kinematic viscosity $\nu$, dynamic viscosity $\eta$ and volumetric viscosity $\eta_{v}$ by (Krüger et al., 2017; Dellar, 2001)

[TABLE]

such that $\sigma^{{\hat{a}}{\hat{b}}}=-\eta(\nabla^{{\hat{a}}}u^{\hat{b}}+\nabla^{{\hat{b}}}u^{\hat{a}})$.

Rather than considering fluid distribution functions $f(\bm{v})$ with continuous velocity space $\bm{v}=(v^{{\hat{\theta}}},v^{{\hat{\varphi}}})$, we discretise the velocity space using $\bm{v}_{\bm{k}}=(v_{k_{\theta}},v_{k_{\varphi}})$. Due to the inertial and reaction force terms in Eq. (15), we need at least a fourth order quadrature ($Q=4$) when non-Cartesian coordinates are employed. We obtain inaccurate simulation results when third order quadrature (or lower) is used. The 16 velocity directions, corresponding to $Q=4$, are illustrated in Fig. 1 (Sofonea et al., 2018). The possible values of $v_{k_{\theta}}$ and $v_{k_{\varphi}}$ ($1\leq k_{\theta},k_{\varphi}\leq 4$) are given as the roots of the fourth order Hermite polynomial,

[TABLE]

We use Hermite polynomials orthogonal with respect to the weight function $e^{-v^{2}/2}/\sqrt{2\pi}$, which are described in detail, e.g. in the Appendix of Ref. (Sofonea et al., 2018). It is worth noting that, unlike standard lattice Boltzmann algorithm, in general the velocity directions do not coincide with the neighbouring lattice points. For simplicity, we have set $k_{B}T=1$ and $m=1$, such that the reference scale for the velocity $v_{\rm ref}=(k_{B}T/m)^{1/2}=1$, which is also the sound speed in an isothermal fluid.

The particle number density $n$ and velocity ${\bm{u}}$ can be computed as zeroth and first order moments of the distribution functions

[TABLE]

With the discretisation of the velocity space, we also replace the Maxwell-Boltzmann equilibrium distribution with a set of distribution functions $f^{\rm eq}_{\bm{k}}$ corresponding to the discrete velocity vectors $\bm{v}_{\bm{k}}$. Due to the use of the vielbein formalism, the expression for $f^{\rm eq}_{\bm{k}}$ coincides with the one employed on the flat Cartesian geometry (Sofonea et al., 2018)

[TABLE]

The quadrature weights $w_{k}$ can be computed using the formula

[TABLE]

where $Q$ is the order of the quadrature and $H_{n}(x)$ is the Hermite polynomial of order $n$. For $Q=4$, the weights have the following values (Sofonea et al., 2018)

[TABLE]

To compute the force terms in the Boltzmann equation, Eq. (16), we consider a unidimensional expansion of the distribution with respect to the velocity space degrees of freedom (Busuioc & Ambru\cbs, 2019; Busuioc et al., 2019a). In particular, the velocity derivatives appearing in Eq. (16) can be computed as

[TABLE]

where the kernels $\mathcal{K}^{H}_{i,m}$ and $\widetilde{\mathcal{K}}^{H}_{i,m}$ can be written in terms of Hermite polynomials (Busuioc & Ambru\cbs, 2019)

[TABLE]

We list below the components of the above matrices for the case of $Q=4$:

[TABLE]

III Drift Dynamics of Fluid Stripes and Droplets

In this section we begin by studying the behaviour of fluid stripes on the torus geometry. By minimising the interface length subject to area conservation, we find there is a second order phase transition in the location of the equilibrium position as we vary the stripe area. In particular we observe bistability when the stripe area exceeds a critical value. We validate the ability of our method to capture this effect in Subsec. III.1. We then consider the Laplace pressure test in Subsec. III.2. The Laplace pressure takes a different form on a curved torus geometry compared to that on a flat geometry (Busuioc et al., 2019b). Furthermore, the approach to equilibrium configuration through a damped harmonic motion is investigated in Subsec. III.3. We show that we recover the damping coefficient and the angular frequency as derived in Busuioc et al. (2019b). Finally, we contrast the drift dynamics of fluid stripes with droplets on the torus in section III.4. While the former drift to the inside of the torus, the latter move to the outside of the torus.

III.1 Equilibrium positions of fluid stripes

The basic idea behind establishing the equilibrium position of fluid stripes is that the interface length must attain a minimum for a fixed stripe area. We consider a stripe of angular width $\Delta\theta$, centred on $\theta=\theta_{c}$, such that its interfaces are located at

[TABLE]

As a convention, here the stripe is identified with the minority, rather than the majority, fluid component. The area $\Delta A$ enclosed between the upper and lower interfaces can be obtained as follows

[TABLE]

where $a=r/R$. The preservation of the area allows the variation of the stripe width $\Delta\theta$ to be related to a variation of the stripe centre $\theta_{c}$. Setting $d\Delta A=0$,

[TABLE]

The total interface length $\ell_{\rm total}=\ell_{+}+\ell_{-}$ can be computed as

[TABLE]

Imposing $d\ell_{\rm total}=0$ yields an equation involving the stripe width $\Delta\theta_{eq}$ and stripe centre $\theta_{c}^{eq}$ at equilibrium

[TABLE]

The above equation has different solutions depending on the stripe width. For narrow stripes, the equilibrium position is located at $\theta_{c}^{eq}=\pi$. There is a critical point corresponding to stripe width $\Delta\theta_{eq}=\Delta\theta_{\rm crit}=2\arccos(a)$, or alternatively stripe area

[TABLE]

For stripes with areas larger than this critical value, two equilibrium positions are possible, namely

[TABLE]

We now reproduce the above phenomenon using our lattice Boltzmann approach. Unless stated otherwise, in section III, we use a torus with $r=0.8$ and $R=2$ ($a=r/R=0.4$). We set the parameters in our free energy model, Eq.(9), to $\kappa=5\times 10^{-4}$ and $A=0.5$, and set the kinematic viscosity $\nu=2.5\times 10^{-3}$ and mobility parameter in the Cahn-Hilliard equation $M=2.5\times 10^{-3}$. Due to its homogeneity with respect to $\varphi$, the system is essentially one dimensional, such that a single node is used on the $\varphi$ direction (i.e., $N_{\varphi}=1$). The discretisation along the $\theta$ direction is performed using $N_{\theta}=320$ nodes. Throughout this paper we ensure that our discretization is such that the spacing is always smaller than the interface width $\xi_{0}$, as given in Eq.(10). The time step is set to $\delta t=5\times 10^{-4}$.

We initialise the fluid stripes using a hyperbolic tangent profile

[TABLE]

where $\phi_{0}$ is an offset due to the Laplace pressure (see next subsection)

[TABLE]

We consider stripes having the same initial position centred at $\theta_{0}=\pi/2$, but initialised with different initial widths $\Delta\theta_{0}$. The area of these stripes is given by

[TABLE]

The equilibrium positions $\theta_{c}^{eq}$ for four different stripes are shown in Fig. 2(a). The first case corresponds to a very large stripe ($\Delta\theta_{0}=0.95\pi$, $\Delta A\simeq 1.88\Delta A_{\rm crit}$), for which the possible equilibria $\theta^{eq}_{c}$ are close to $\pi/2$ and $3\pi/2$. Due to the initial condition, the stripe is attracted by the equilibrium point on the upper side of the torus, where it will eventually stabilise. As the stripe size decreases, its kinetic energy as it slides towards the equilibrium point will be sufficiently large for it to go over the “barrier” at $\theta_{c}=\pi$ to the lower side of the torus. Because of energy loss due to viscous dissipation, its kinetic energy may be insufficient to overcome this barrier again, so the stripe remains trapped on the lower side. This is the case for the second stripe having $\Delta\theta_{0}=0.65\pi$ ($\Delta A\simeq 1.29\Delta A_{\rm crit}$). Further decreasing the stripe size causes the peak at $\theta_{c}=\pi$ to also decrease, allowing the stripe to overcome it a second time as it migrates back towards the upper side. The third stripe, initialised with $\Delta\theta_{0}=0.6\pi$ ($\Delta A\simeq 1.19\Delta A_{\rm crit}$), stabilises on the upper side of the torus. Finally, the fourth stripe is initialised with $\Delta\theta_{0}=0.3\pi$, such that its area $\Delta A\simeq 0.59\Delta A_{\rm crit}$ is below the critical value. Thus, the fourth stripe will perform oscillations around the equilibrium at $\theta_{c}=\pi$, where it will eventually stabilise.

Judging by the number of times that the stripe centre $\theta_{c}$ crosses the barrier at $\theta_{c}=\pi$, two types of stripes having $\Delta A>\Delta A_{\rm crit}$ can be distinguished: (i) the ones that cross the $\theta_{c}=\pi$ line an even number of times stabilise on the upper half of the torus, while (ii) the ones that cross it an odd number of times stabilise on the lower half of the torus. This is presented in Fig. 2(b), where the equilibrium position $\theta_{c}^{eq}$ for stripes initialised at $\theta_{0}=\pi/2$ is represented as a function of $\Delta\theta_{0}$ in comparison with the analytical predictions in Eq. (33).

Panels (c-e) in Fig. 2 illustrate the three scenarios where the stripes are equilibrated at $\theta_{c}^{eq}>\pi$, $\theta_{c}^{eq}=\pi$, and $\theta_{c}^{eq}<\pi$ respectively. The total interface lengths $\ell_{\rm total}$ ($\sim\Psi$) for the stripes shown in (c-e) are represented in panels (f-h) of Fig. 2. The interface lengths at the equilibrium positions corresponding to the initial state as well as to the turning points corresponding to half-periods are also shown using symbols, numbered sequentially in the legend ([math] corresponds to the initial state). It can be seen that $\ell_{\rm total}$ measured at these turning points decreases monotonically. When $\ell_{\rm total}$ decreases below its value at $\theta_{c}=\pi$, the stripe centre can no longer cross the $\theta_{c}=\pi$ line and becomes trapped in one of the minima.

Fig. 3 further summarises the location of the equilibrium stripe position as a function of the stripe width $\Delta\theta_{0}$ and the radii ratio $a=r/R$. Our simulations are performed by keeping $R=2$ constant, such that the various values of $a$ are obtained by changing $r$. As before, the stripe is initialised at $\theta_{0}=\pi/2$. Moving from the top right corner of the diagram towards the bottom left corner, the subsequent regions distinguish between whether the stripes stabilise on the top half ($<\pi$) or on the bottom half ($>\pi$), depending on the number of times that $\theta_{c}$ crosses $\pi$. In the bottom left corner, the stripes stabilise at $\theta_{c}^{eq}=\pi$. The black region between the purple band and the lower left region corresponds to stripes that cross $\pi$ more than $3$ times but stabilise away from $\pi$ ($\theta_{c}^{eq}\neq\pi$). Due to the diffuse nature of the interface, the stripes evaporate when $r\Delta\theta\lesssim 5\xi_{0}$ ($\xi_{0}=\sqrt{\kappa/A}\simeq 0.031$). These regions correspond to the top left and bottom right corners of the diagram and are shown in red.

III.2 Laplace pressure test

Since the stripe interfaces have a non-vanishing curvature, it can be expected that there will be a pressure difference across this interface. This pressure difference is often termed the Laplace pressure. This pressure difference was recently derived analytically on a torus and the result is (Busuioc et al., 2019b)

[TABLE]

This expression can be simplified for the two types of minima highlighted in the previous subsection

[TABLE]

We remind the readers that, on the first branch, $\theta_{c}^{eq}=\pi$. On the second branch, the equilibrium position is determined via $a\cos\theta_{c}^{eq}+\cos(\Delta\theta_{eq}/2)=0$.

In order to validate our numerical scheme against the Laplace pressure test on the torus, we perform numerical simulations for two values of $\kappa$ in our free energy model, $\kappa=2.5\times 10^{-4}$ and $5\times 10^{-4}$. These effectively change the surface tension and interface width in our simulations, see Eq. (10). All of the other simulation parameters are kept the same as in the previous subsection: $R=2$, $r=0.8$, $A=0.5$, $\nu=2.5\times 10^{-3}$ and $M=2.5\times 10^{-3}$. We consider stripes of various areas $\Delta A$ in Fig. 4. After the stationary state is reached, we measure the total pressure $p=p_{\rm i}+p_{\rm binary}=nk_{B}T+A(-\frac{1}{2}\phi^{2}+\frac{3}{4}\phi^{4})$ in the interior and exterior of the stripe, and compute the difference $\Delta p$ between these two values. The simulation results are shown using dashed lines and symbols in Fig. 4. We observe an excellent agreement with the analytic results, Eq. (38), which are shown using the solid lines.

III.3 Approach to equilibrium

For stripes close to their equilibrium position, the time evolution of the departure $\delta\theta=\theta_{c}-\theta_{c}^{eq}$ can be described as a damped harmonic oscillation:

[TABLE]

where the damping coefficient $\alpha=\alpha_{\nu}+\alpha_{\mu}$ receives contributions from the viscous damping due to the fluid (Busuioc et al., 2019b)

[TABLE]

as well as from the diffusion due to the mobility of the order parameter, $\alpha_{\mu}$ (Busuioc et al., 2019b). In the applications considered in this paper, $\alpha_{\mu}\ll\alpha_{\nu}$, such that we will only consider the approximation $\alpha\simeq\alpha_{\nu}$. For the case of subcritical stripes ($\Delta A<\Delta A_{\rm crit}$), which equilibrate at $\theta_{c}^{eq}=\pi$, the oscillation frequency is (Busuioc et al., 2019b)

[TABLE]

For the supercritical stripes, ($\Delta A>\Delta A_{\rm crit}$), when the equilibrium position is at $a\cos\theta_{c}^{eq}+\cos(\Delta\theta_{eq}/2)=0$, $\omega_{0}^{2}$ is given by

[TABLE]

We will now demonstrate that our lattice Boltzmann implementation captures the dynamical approach to equilibrium as described by the analytical results. First, we consider a torus with $r=0.8$ and $R=2$ ($a=0.4$), and set $\kappa=5\times 10^{-4}$, $A=0.5$ and $\tau=M=2.5\times 10^{-3}$. The number of nodes is $N_{\theta}=320$, and the order parameter $\phi$ is initialised according to Eq. (34), where the stripe centre is located at an angular distance $\delta\theta_{0}=\theta_{c}-\theta_{c}^{eq}=-\pi/20$ away from the expected equilibrium position. Fig. 5 shows a comparison between the numerical and analytical results for the time evolution of $(\theta_{c}^{eq}-\theta_{c})/\pi$ for the cases (a) $\theta_{c}^{eq}=\pi$ with initial stripe width $\Delta\theta_{0}=0.28\pi$, and (b) $\theta_{c}^{eq}=3\pi/4$ with $\Delta\theta_{0}=0.786\pi$. For the analytical solution, the angular velocity $\omega_{0}$ is computed using Eqs. (41) and (42) for cases (a) and (b) respectively, and the damping factor $\alpha\simeq\alpha_{\nu}$ is computed using Eq. (40). We have also set the offset to $\varsigma=0$. It can be seen that the analytic expression provides an excellent match to the simulation results for the stripe that goes to $\theta_{c}^{eq}=\pi$. For the stripe equilibrating to $3\pi/4$, we observe a small discrepancy, especially during the first oscillation period. However, the overall agreement is still very good.

Next we consider three tori having radii ratio $a=r/R=0.4$, with $r=0.8$, $1$ and $1.2$, and perform two sets of simulations. In the first set of simulations, the initial configuration corresponds to a stripe centred on $\theta_{0}=0.95\pi$, with initial width $\Delta\theta_{0}=0.28\pi$. These stripes relax towards $\theta_{c}^{eq}=\pi$. In the second set of simulations, the stripes are initially centred at $\theta_{0}=0.7\pi$, and they equilibrate at $\theta_{c}^{eq}=3\pi/4$, with initial width $\Delta\theta_{0}=0.786\pi$. The simulations are performed using $N_{\theta}=320$, $400$ and $480$ nodes for $r=0.8$, $1$ and $1.2$, respectively. The best-fit values of $\alpha$ and $\omega_{0}$ for the three torus geometries are shown in Fig. 6 and Fig. 7 respectively as function of the kinematic viscosity $\nu$ (varying between $2.5\times 10^{-3}$ and $7.5\times 10^{-3}$) at $\kappa=5\times 10^{-4}$; and the surface tension parameter $\kappa$ (varying between $2.5\times 10^{-4}$ and $6.25\times 10^{-4}$) at $\nu=2.5\times 10^{-3}$. For each simulation, Eq. (39) is fitted to the numerical data for the time evolution of the stripe centre as it relaxes towards equilibrium, using $\alpha$ and $\omega$ as free parameters, while $\varsigma=0$. For simplicity, we have used $M=\nu$ and $A=0.5$ in Fig. 6 and Fig. 7. Panels (a) in Fig. 6 and Fig. 7 correspond to stripes equilibrating at $\theta_{c}^{eq}=\pi$, while panels (b) in Fig. 6 and Fig. 7 are for $\theta_{c}^{eq}=3\pi/4$. It can be seen that the analytic expressions are in good agreement with the numerical data in all instances simulated.

Finally we investigate the applicability of Eqs. (41) and (42) with respect to various values of the stripe area, $\Delta A$. The simulations are now performed on the torus with $r=0.8$ and $R=2$, using $\kappa=5\times 10^{-4}$, $A=0.5$, $\tau=M=2.5\times 10^{-3}$. Fig. 8 shows the values of $\omega_{0}$ obtained by fitting Eq. (39) to the numerical data (points) and the analytic expressions (solid lines). As before, for the fitting, we set $\varsigma=0$, and use $\alpha$ and $\omega_{0}$ as free parameters. An excellent agreement can be seen, even for the nearly critical stripe, for which $\omega_{0}$ is greatly decreased.

III.4 Droplets on Tori

We will now show that, when placed on a torus, a fluid droplet will also exhibit a drift motion. However, in contrast to stripes, the drops will move towards the outer rather than the inner side of the torus. To study this phenomenon quantitatively, we initialise drops on a torus using the following equation

[TABLE]

where $r=\sqrt{(x-x_{c})+(y-y_{c})+(z-z_{c})}$ is the Euclidean distance between the point with coordinates $(x,y,z)$ and the centre of the drop $(x_{c},y_{c},z_{c})$, corresponding to $(\theta,\varphi)$ and $(\theta_{0},0)$ in polar coordinates respectively. The relation between the Cartesian and polar coordinates are given in Eq. (49) of Appendix A. The parameter $\theta_{0}$ represents the center of the drop, while $R_{0}$ is a measure of its radius. $\xi_{0}$ is the interface width derived for the Cartesian case. In principle the interfacial profile will be different on a torus, but currently we are not aware of a closed analytical formula. We also do not introduce in Eq. (43) the offset $\phi_{0}$ responsible for the Laplace pressure difference, since the analysis of this quantity is less straightforward than for the azimuthally-symmetric stripe domains discussed in the previous subsections.

In order for the drops to have approximately the same areas, for a given value of $\theta_{0}$, $R_{0}$ is obtained as a solution of

[TABLE]

where the first term in the parenthesis corresponds to the configuration when the droplet is centred on the outer equator and has $R_{0}=30\xi_{0}$. The drift phenomenon we report here is robust with respect to the drop size, but we choose a relatively large drop size because small drops are known to evaporate in diffuse interface models. The simulation parameters are the same as in Sec. III.3, namely $r=0.8$, $R=2$, $\kappa=5\times 10^{-4}$, $A=0.5$ and $\tau=M=2.5\times 10^{-3}$.

As shown in Fig. 9(a), similar to the stripe configuration in the previous sub-section, we observe a damped oscillatory motion. Here the three drops are initialised at different positions on the torus. Moreover, as is commonly the case for an underdamped harmonic motion, the drops initially overshoot the stable equilibrium position, but they eventually relax to the minimum energy configuration. For the drops we find that all drops eventually drift to $\theta=0$ (the outer side of the torus). Typical drop configurations during the oscillatory motion are shown in Fig. 9(b-d). Compared to the oscillatory dynamics for the stripe configurations, we also observe that the oscillation dies out quicker for the drops.

IV Phase Separation

In this section we investigate binary phase separation on the torus and compare the results against those on flat surfaces. We consider hydrodynamics and diffusive regimes for both even (section IV.1) and uneven (section IV.2) mixtures.

The fluid order parameter at lattice point $(s,q)$ is initialised as

[TABLE]

where $\overline{\phi}$ is a constant and $(\delta\phi)_{s,q}$ is randomly distributed between $(-0.1,0.1)$. We characterise the coarsening dynamics using the instantaneous domain length scale $L_{d}(t)$, computed using the following function:

[TABLE]

where $A_{\rm total}$ is the total area of the simulation domain. The total interface length at time $t$, $L_{I}(t)$, is computed by visiting each cell $(s,q)$ exactly once, starting from the bottom left corner, where $s=q=1$, and progressing towards the top right corner, where $s=N_{1}$ and $q=N_{2}$. $N_{1}=N_{x}$ and $N_{2}=N_{y}$ for the Cartesian domains and $N_{1}=N_{\varphi}$ and $N_{2}=N_{\theta}$ for the torus domains. For each cell where $\phi_{s,q}\times\phi_{s+1,q}<0$, the length of the vertical interface between the $(s,q)$ and $(s+1,q)$ is added to $L_{I}$. In the case of the Cartesian geometry, this length is $\delta y$, while for the torus, the length is given by $r\delta\theta$. Similarly, if $\phi_{s,q}\times\phi_{s,q+1}<0$, the length of the horizontal interface ($\delta x$ for the Cartesian case and $(R+r\cos\theta_{q+1/2})\delta\varphi$ for the torus case, where $\theta_{q+1/2}=\theta_{q}+\delta\theta/2$ is the coordinate of the cell interface) is added to $L_{I}$. The periodic boundary conditions allow the cells with $(N_{1}+1,q)$ and $(s,N_{1}+1)$ to be identified with the cells $(1,q)$ and $(s,1)$, respectively.

Unless specified otherwise, we use the following parameters in this phase separation section: $M=\tau=2.5\times 10^{-3}$, $\delta t=5\times 10^{-4}$, $A=0.5$ and $\kappa=5\times 10^{-4}$. In the initial state, the distributions for the LB solver are initialised using Eq. (21) with a constant density $n_{0}=20$ and vanishing velocity.

IV.1 Even mixtures

IV.1.1 Cartesian Geometry

We begin by considering the coarsening dynamics of a phase separating binary fluid with even mixtures on a flat two-dimensional surface. We use a simulation domain of $N_{x}\times N_{y}=512\times 512$ with a grid spacing of $\delta x=\delta y=0.02$. The linear size of the simulation domain is $L=512\times 0.02=10.24$ and its total area is $A_{\rm total}=L^{2}$.

As shown in Fig.  10(a), we observe that the fluid domain grows with an exponent of $2/3$. This exponent is often associated with the so-called inertial-hydrodynamics regime for binary fluid phase separation in three dimensions (Bray, 2002; Kendon et al., 2001). However, in two dimensions, it has been argued that self-similar growth in the inertial-hydrodynamics regime may be absent (Wagner & Yeomans, 1997). The apparent exponent of $2/3$ is really due to a mixture of viscous exponent of $1$ for the growth of the connected domains and an exponent of $1/3$ for the diffusive dissolution of circular droplets.

Classical morphologies typical of a spinodal decomposition phenomenon are shown in Fig. 10(c-f). The deviation from this apparent scaling law is observed at early times when domains of fluid components A and B are formed from the initial perturbation, and at late times when the domains become comparable in size to the simulation box. For the latter, there are very few domains left (see e.g. Fig. 10(f)) and coarsening slows down because of the lack of coalescence events between the fluid domains.

To access the diffusive regime, in this work we remove the advection term in the Cahn-Hilliard equation and decouple it from the Navier-Stokes equation. In this case, coarsening can only occur via diffusive dynamics, and indeed we do observe a growth exponent of $1/3$, as shown in Fig. 10(b), as expected for diffusive dynamics (Bray, 2002; Kendon et al., 2001). Representative configurations from the coarsening evolution are shown in Fig. 10(g-j). These snapshots look somewhat similar to those shown in Fig. 10(c-f) for the apparent $2/3$ scaling regime. The key difference between the morphologies is that more small droplets are accumulated during coarsening when hydrodynamics is on. It is also worth noting that the coarsening dynamics are much slower in the diffusive regime. At late times we see a deviation from the diffusive scaling exponent, where $L_{d}(t)$ appears to grow faster than $1/3$ exponent. In this limit, as illustrated in Fig. 10(j), the increase in $L_{d}(t)$ is primarily driven by finite size effects.

IV.1.2 Torus Geometry

We now consider the coarsening dynamics of a phase separating binary fluid on the surface of a torus. Initially we simulate a torus domain with $R=2.5$ and $r=1$ ($a=r/R=0.4$). These parameters are chosen such that the total area, $A_{\rm total}=4\pi^{2}rR$, is close to the one employed in the Cartesian case. The $\varphi$ direction is discretised using $N_{\varphi}=800$ nodes, while the $\theta$ direction is discretised using $N_{\theta}=400$ nodes. The fluid order parameter at lattice point $(s,q)$ is initialised according to Eq. (45) with $\overline{\phi}=0$.

Our simulation results are shown in Figs. 11(a) and 11(b) respectively for cases with and without coupling to hydrodynamics. Qualitatively we find a similar behaviour to the results obtained in the Cartesian case, Fig. 10. In panel (a), it can be seen that $L_{d}(t)$ grows with an apparent exponent of $2/3$ when hydrodynamics is on. Turning off the hydrodynamics, the $1/3$ diffusive exponent emerges, as demonstrated in panel (b). The coarsening dynamics is also much faster with hydrodynamics on the torus. Snapshots of the order parameter configuration at various times for the case of the even mixture with and without hydrodynamics are shown in panels (c-e) and (f-h) respectively.

Quantitatively, we observe that finite size effects occur earlier (smaller $L_{d}$) for the torus considered in Fig. 11 compared to the Cartesian case. This is expected since the effective length scale in the poloidal direction, $2\pi r$, is smaller than the width of the simulation box in the Cartesian case, even though the total surface areas are comparable. Indeed, we can observe that the departure from the $2/3$ (panel a) and $1/3$ (panel b) exponents occur when the fluid domains start to wrap around the circle in the poloidal direction.

In Fig. 12 we further show simulation results for a thicker ($R=2$ and $r=1.25$; $a=0.625$) and a thinner ($R=5$ and $r=0.5$; $a=0.1$) torus, having a total area equal to the one considered at the beginning of this section. The simulation parameters are kept the same as before, except that for the thicker torus, the time step must be decreased down to $\delta t=5\times 10^{-5}$ since the minimum spacing along the $\varphi$ direction occurring on the inner equator is $2\pi(R-r)/N_{\varphi}\sim 0.00589$. Comparing Figs. 11(a), 12(a) and 12(b), we can further conclude that finite size effects appear sooner for the thinner torus and later for the thicker one. This further strengthens the argument that the determining lengthscale for the finite size effects is the circumference in the poloidal direction, rather than the circumference on the inner side of the torus (at $\theta=\pi$), $2\pi(R-r)$. Otherwise, the thicker torus should display finite size effects the earliest among the three geometries simulated.

Given the fluid stripes are generally formed in the poloidal rather than the toroidal direction during phase separation, the drift phenomenon reported in sub-section III.3 for stripe configurations cannot be clearly visualised. However, domain drifts for drops, as reported in section III.4, can be seen in Figs. 11 and 12 during the late stages of the coarsening phenomenon. This drift phenomenon can be observed even clearer when we study uneven mixtures, as discussed in the next sub-section.

IV.2 Uneven Mixtures

IV.2.1 Cartesian Geometry

The simulation results for a mixture with asymmetric composition are shown in Fig. 13. We use the same simulation parameters as in Fig 10, except that $\overline{\phi}=-0.3$. Fig. 13(a) shows how the typical domain size scales with time both when hydrodynamics is turned on and off. Interestingly, in both cases we observe an exponent of $1/3$, albeit with different prefactors. This is in contrast to our results for the even mixtures, when an apparent exponent of $2/3$ is obtained with hydrodynamics. It has been suggested in the literature that the effect of hydrodynamics decreases as a function of the asymmetry of the mixture, though we do not yet know of a convincing systematic study of this effect. For example, Wagner & Cates (2001) showed that at high concentrations droplets with hydrodynamics exhibit the viscous hydrodynamic coarsening regime, but as droplet coalescence is reduced at lower volume fractions the effect of hydrodynamics diminishes. Here we observe the limit where the scaling is typical of that for diffusive dynamics.

The fluid configurations at various times in the simulation are shown in Fig. 13, panels (b-e), when hydrodynamics is taken into account. These can be compared to Fig. 13, panels (f-i), when the advection term is switched off in the Cahn-Hilliard equation. The differences are mainly that the morphologies with hydrodynamics are coarsening faster. In the non-hydrodynamic simulation there are more coalescence events visible because the restoration of a round shape takes more time. Thus, while the scaling exponent is the same with and without hydrodynamics, hydrodynamics still plays an important role in that it allows coalescing droplets return to a round shape more quickly.

IV.2.2 Torus Geometry

Here we consider a torus geometry with $R=2.5$ and $r=1$ (same geometry and simulation parameters as in Fig. 11), and the order parameter is initialised according to Eq. (45) with $\overline{\phi}=-0.3$. The simulation results for the uneven mixture are shown in Fig.  14. Quantitatively, we find a similar behaviour as for the Cartesian case. Both when hydrodynamics is turned on and off, we observe a $1/3$ exponent in our simulations. Similar to the even mixture shown in Fig. 11, we also find that finite size effects occur earlier (smaller $L_{d}$) for the torus compared to the Cartesian geometry. As discussed in the case of even mixtures, this occurs when the fluid domains start to wrap around the circle in the poloidal direction. Snapshots of the fluid configurations during phase separation are shown in panels (b-d) and (e-g) respectively for simulations with and without hydrodynamics.

At late times, the effect of the curvature on the domain dynamics becomes important. In Sec. III.4 we discussed how droplet domains migrate to the outer side of the torus. To quantify this effect during phase separation of uneven mixtures, we consider the average of $\phi$ with respect to the azimuthal angle $\varphi$:

[TABLE]

The discrete equivalent of the above relation is

[TABLE]

We plot $(\braket{\phi}+1)/2$ as a function of the poloidal angle $\theta$ at various times in Fig. 15(a). At late times, see e.g. Fig. 15(c), the typical configuration corresponds to the majority phase ($\phi=-1$) forming a continuum with several large droplets of the minority phase ($\phi=+1$) primarily in the outer side of the torus. At $t=18000$ [Fig. 15(d)], when the steady state is reached, the inner stripe spans $0.65\pi\lesssim\theta\lesssim 1.35\pi$. Our convention is to identify the stripe with the minority fluid component. The maximum of $(\braket{\phi}+1)/2$ is clearly reached at $\theta=0$, indicating that the outer side of the torus is populated by droplets centered on $\theta=0$.

V Conclusions

In this work we developed a vielbein lattice Boltzmann scheme to solve the hydrodynamics equations of motion of a binary fluid on an arbitrary curved surface. To illustrate the application of our vielbein lattice Boltzmann method to curved surfaces, here we focussed on the torus geometry and studied two classes of problems. First, due to the non-uniform curvature present on a torus, we showed drift motions of fluid droplets and stripes on a torus. Such dynamics are not present on a flat surface or on surfaces with uniform curvature. Interestingly the fluid droplets and stripes display preference to different regions of the torus. Fluid droplets migrate to the outer side of the torus, while fluid stripes move to the inner side of the torus. The exhibited dynamics are typical of a damped oscillatory motion. Moreover, for the fluid stripes, the corresponding dynamics can effectively be reduced to a one-dimensional problem by taking advantage of the symmetry with respect to the azimuthal angle. Our simulation results are in excellent agreement with the analytical predictions for the equilibrium position of the stripes, the Laplace pressure difference between the inside and outside of the stripes, and the relaxation dynamics of the stripes towards equilibrium.

We also studied phase separation dynamics on tori of various shapes. For even mixtures, $2/3$ and $1/3$ scaling exponents characteristic of hydrodynamics and diffusive regimes are observed. In contrast, for uneven mixtures, we only observe a $1/3$ scaling exponent both when hydrodynamics is turned on and off. Compared to Cartesian geometry, we saw that finite size effects kick in earlier for the torus geometry. By comparing the results for three torus aspect ratios, we conclude that the determining lengthscale for the finite size effects seems to be the perimeter in the poloidal direction, corresponding to fluid domains wrapping around the circle in the poloidal direction. That the stripes are observed to form in the poloidal rather than the toroidal direction prevents the observation of drift motion of fluid stripes towards the inner side of the torus during phase separation. However, the domain drifts for fluid drops to the outer side of the torus can be clearly observed at the late stage of phase separation.

While we focussed on the torus geometry, our approach can be applied to arbitrary curved geometry. Moreover, one interesting area for future work is to expand the method to account for unstructured mesh, where the geometrical objects needed for the Boltzmann equation must be evaluated numerically. A major challenge is to construct a numerical scheme which is accurate to second order or higher. Another important avenue for future investigations is to couple the hydrodynamics equations of motion with more complex dynamical equations, such as those for (active and passive) liquid crystals and viscoelastic fluids. We believe this work extends the applicability of the lattice Boltzmann approaches to a new class of problems, complex flows on curved manifolds, which are difficult to carry out using the standard lattice Boltzmann method.

Acknowledgements: We acknowledge funding from EPSRC (HK; EP/J017566/1 and EP/P007139/1), Romanian Ministry of Research and Innovation (VEA and SB; CCCDI-UEFISCDI, project number PN-III-P1-1.2-PCCDI-2017-0371/VMS, within PNCDI III), and the EU COST action MP1305 Flowing Matter (VEA and HK; Short Term Scientific Mission 38607). VEA gratefully acknowledges the support of NVIDIA Corporation with the donation of a Tesla K40 GPU used for this research. VEA and SB thank Professor Victor Sofonea (Romanian Academy, Timi\cbsoara Branch) for encouragement, as well as for sharing with us the GPU infrastructure available at the Timi\cbsoara Branch of the Romanian Academy.

Appendix A Application of the vielbein method to the torus geometry

The derivation of the Boltzmann equation, Eq. (15), written in conservative form with respect to vielbein vector fields is discussed in Busuioc & Ambru\cbs (2019). Using Eq. (15) as a starting point, here we present generic main steps required to write down the Boltzmann equation for any arbitrary curved surface. For concreteness, we focus on the torus geometry in this paper.

*Parametrising the surface. * As a two-dimensional manifold, a surface needs two coordinates $q^{1}$ and $q^{2}$ to be parametrised. In the case of a torus of inner radius $r$ and outer radius $R$, the parametrisation can be chosen in terms of the angles $\theta\in[0,2\pi)$ and $\varphi\in[0,2\pi)$ as follows:

[TABLE]

and the system is periodic with respect to both of these angles. 2. 2.

*Writing down the line element. * Differentiating the functions $x$, $y$ and $z$ with respect to $q^{1}$ and $q^{2}$ yields the formula

[TABLE]

where $\{i,j\}\in\{1,2,3\}$, $1\leq a,b\leq 2$ and $g_{ab}$ are the components of the metric tensor. In the case of Eq. (49), the line element becomes:

[TABLE]

leading to the metric tensor components

[TABLE] 3. 3.

Constructing the vielbein field. The vielbein vector frame consists of the vectors $\bm{e_{\hat{a}}}=e_{\hat{a}}^{a}\bm{\partial_{a}}$ which satisfy:

[TABLE]

Since Eq. (53) is invariant under the action of the orthogonal group with respect to the hatted indices, the vielbein is defined up to an arbitrary rotation. After fixing the vielbein, the vielbein one-form co-frame denoted via $\bm{\omega^{\hat{a}}}=\omega^{\hat{a}}_{a}\bm{dq^{a}}$ is uniquely fixed by Eq. (3).

For the torus geometry, the natural choice is to take

[TABLE] 4. 4.

Computing the Cartan coefficients. The commutator of two vector fields $\bm{u}$ and $\bm{v}$ is another vector field, denoted by $[\bm{u},\bm{v}]=(u^{a}\partial_{a}v^{b}-v^{a}\partial_{a}u^{b})\bm{\partial_{b}}$. The contraction between the co-frame one-form $\bm{\omega^{\hat{c}}}$ and the commutator of the tetrad frame vector fields $\bm{e_{\hat{a}}}$ and $\bm{e_{\hat{b}}}$ defines the Cartan coefficient $c_{{\hat{a}}{\hat{b}}}{}^{{\hat{c}}}$ (7), via the following relation:

[TABLE]

In the case of the torus, the commutator of the vielbein vectors $\bm{e_{\hat{\theta}}}$ and $\bm{e_{\hat{\varphi}}}$ is

[TABLE]

leading to the Cartan coefficients

[TABLE] 5. 5.

Computing the connection coefficients. On curved surfaces, the ordinary partial derivative operator $\partial_{a}$ must be replaced by a covariant derivative which ensures that the resulting vector or tensor is still contained in the tangent space. This is ensured starting from the covariant derivative of the basis vectors:

[TABLE]

In the case of the torus, the only non-vanishing connection coefficients are

[TABLE] 6. 6.

Writing the Boltzmann equation. Plugging Eq. (59) into Eq. (15) yields the Boltzmann equation for the torus geometry, Eq. (16).