Wetting and Spreading Behavior of Axisymmetric Compound Droplets on Curved Solid Walls Using Conservative Phase Field Lattice Boltzmann Method

Abstract

Compound droplets have received increasing attention due to their applications in many several areas, including medicine and materials. Previous works mostly focused on compound droplets on planar surfaces and, as such, the effects of curved walls have not been studied thoroughly. In this paper, the influence of the properties of curved solid wall (including the shape, curvature, and contact angle) on the wetting behavior of compound droplets is explored. The axisymmetric lattice Boltzmann method, based on the conservative phase field formulation for ternary fluids, was used to numerically study the wetting and spreading of a compound droplet of the Janus type on various curved solid walls at large density ratios, focusing on whether the separation of compound droplets occurs. Several types of wall geometries were considered, including a planar wall, a concave wall with constant curvature, and a convex wall with fixed or variable curvature (specifically, a prolate or oblate spheroid). The effects of surface wettability, interfacial angles, and the density ratio (of droplet to ambient fluid) on the wetting process were also explored. In general, it was found that, under otherwise identical conditions, droplet separation tends to happen more likely on more hydrophilic walls, under larger interfacial angles (measured inside the droplet), and at larger density ratios. On convex walls, a larger radius of curvature of the surface near the droplet was found to be helpful to split the Janus droplet. On concave walls, as the radius of curvature increases from a small value, the possibility to observe droplet separation first increases and then decreases. Several phase diagrams on whether droplet separation occurs during the spreading process were produced for different kinds of walls to illustrate the influences of various factors.

1. Introduction

The movement of compound droplets composed of different components occurs commonly in nature and engineering applications. Due to the potential applications in manufacturing of inkjet printing [1], drug delivery [2], microfluidic preparation [3], and food and cosmetics [4], the research on compound droplets continues to draw significant attention. Stone et al. [5] used the boundary integral method with axisymmetric assumption to numerically study the core-shell compound droplet with large deformation in linear flow; they also studied its rupture in the flow field. Hua et al. [6] used a numerical model based on the immersed boundary method to study the effects of droplet radius, surface tension ratio, and internal droplet position on the deformation of compound droplets in shear flows in both two and three dimensions and calculated the inclinations of internal and external droplets. Xu et al. [7] studied the coalescence kinetics of immiscible droplets through experiments, measured the evolution of the liquid bridge, compared it with miscible droplets, and proposed a theoretical model to analyze the influence of immiscibility in droplet coalescence. Yang et al. [8] used the phase field model to numerically study the dynamics of axisymmetric compound liquid threads and found that a larger inner liquid radius would significantly delay the evolution of compound liquid threads. They also found that when the inner liquid line radius was larger, the middle compound liquid thread would not shrink into droplets, while larger viscosity and surface tension ratios could delay the evolution of compound liquid threads. Wöhrwag et al. [9] proposed a model to simulate ternary fluid systems by incorporating a new free energy formula into the entropic lattice Boltzmann method (LBM). This model allows us to simulate problems with large density ratios and large surface tension ratios covering partial and complete wetting states.

All the above studies focused on the motion of compound droplets without a solid boundary. The fluid-solid interaction during the spreading of compound droplets on solid walls is more complex, and the related numerical work is more challenging. As the solid wall may have different shapes, curvatures, and wettabilities, it is worth conducting a lot of research in this area. Zhang et al. [10] proposed a geometric wetting condition for simulating ternary fluid flow by using the weighted contact angle model; they verified this model by simulating compound droplets on a substrate. Zhang et al. [11] used a diffuse interface method to simulate the morphology transformation of compound droplets on a straight wall, determined the irreversible and reversible configuration transformations, and verified the final equilibrium morphology of the droplets. Shi and Wang [12] proposed a phase field model for the dynamics of three-component immiscible fluids on a solid wall, used it to simulate the compound drop dynamics on a plane, and showed how the stability and deformation of the compound drop depend on the viscosity ratio of different fluids. In addition, they proposed an effective adaptive mesh refinement technique to improve the computation speed. Bhopalam et al. [13] combined a phase-field model using the ternary Navier–Stokes–Cahn–Hilliard equations with a neo-Hookean model for the solid, and studied the static wetting of a Janus compound droplet on a soft solid. They also simulated the interesting capillary origami of different compound droplets and the configuration shift of compound droplets on a soft solid under different physical parameters. Yang et al. [14] developed a novel diffuse-interface model to describe compound droplets in contact with solid. The wetting behavior of compound droplet on flat, inclined, spherical, and rough substrates was studied to verify the accuracy of their model. Huang [15] developed a hybrid lattice-Boltzmann finite-difference method for the simulation of ternary fluids near various immersed solid objects. Li et al. [16] proposed a numerical method for simulating three-phase flow with moving contact lines on complex surfaces within the framework of the color-gradient LBM. They validated their method through the simulations of a Janus droplet resting on a flat surface, a Janus droplet deposited on a cylinder, and the capillary intrusion of ternary fluids for various viscosity ratios. It was then used to study the dynamics of a compound droplet passing through an array of cylinders subject to a uniform incoming flow. Zhang and Huang [17] used the numerical method proposed in [15] to explore the wetting and spreading behavior of compound droplets on a two-dimensional wedge. Chowdhury et al. [18] studied the dynamics of capturing bubbles at a liquid-liquid interface under axisymmetric conditions by using the ternary Navier–Stokes–Cahn–Hilliard equations. They used a cone with varying wettability and a certain cone angle to remove the bubble, and then performed detailed numerical analyses of bubble detachment for a wide range of flow configurations.

There are many studies on three-phase flows, but there are still some insufficiencies in the existing works. For example, in Refs. [10,11,12], the solid walls involved are all planar, and the wetting of compound droplets on curved solid walls was not studied. Ref. [13] considered the challenging fluid-solid interactions, but it only studied two-dimensional problems. The test cases in Ref. [14] were mostly two-dimensional (only one was three-dimensional). Ref. [15] adopted the phase field theory based on the Cahn–Hilliard equation, and used the finite difference method to solve the phase field equations; it could not handle problems with large density ratios. Refs. [16,17] only studied problems with small to intermediate density ratios in two dimensions. Unlike most previous works, the numerical method in this paper is an axisymmetric LBM based on the conservative Allen–Cahn equations (CACEs) and it can simulate more realistic ternary fluid problems with large density ratios under axisymmetric conditions. It greatly saves computational resources when compared with full three-dimensional simulations. In addition, one major focus of this study is on the influence of curved solid wall properties (including the shape, curvature, and contact angle) on the wetting behavior of compound droplets, which has been rarely investigated in depth in the literature. This study was not only driven by curiosity, but also motivated by the following considerations. First, the knowledge on such wetting behavior may be helpful to understand the interaction between compound droplets and micro solid particles of various shapes. Second, the topological changes of compound droplets during such processes might be useful in some industrial processes aiming to separate different fluid components. Third, the diverse equilibrium morphologies of compound droplets on different curved surfaces could be employed to fabricate certain small components with special shapes (like microlens) in future.

In the remainder of this paper, Section 2 gives the numerical method used in the simulations, Section 3 presents some verifications of the method and the analyses on the results of compound droplets spreading on various curved walls, and Section 4 concludes this paper.

2. Numerical Method

The numerical method includes two parts: (1) the lattice Boltzmann equations (LBEs) for the axisymmetric CACEs used for the interface dynamics of a ternary fluid system, and (2) the LBEs for the axisymmetric incompressible NSEs to simulate the hydrodynamics with interfacial tension effects. The two parts are coupled with each other through the velocity, the physical properties of the fluids (i.e., density and viscosity), and the interfacial tension force. The two components are briefly introduced as follows.

2.1. Axisymmetric LBEs for the CACEs

In the conservative phase-field formulation for ternary fluids, the volume fraction (order parameter) $c_i(i=\mathrm{1,2},3)$ is used to identify fluid $\mathrm{i}$ and $c_i$ and satisfies $0\le c_i\le 1$ and the relation ${\sum }_{i=1}^3{c}_i=1$. At a given point, $c_i=1$ means that it is completely occupied by fluid i and $c_i=0$ means that fluid i is totally absent. There exist transition regions with $0<c_i<1$ that contain partially fluid i and also one or two of the other fluids. The interface is represented by the contour at $c_i=0.5$. Because ${\sum }_{i=1}^3{c}_i=1$, only two volume fractions are independent. Without loss of generality, we choose $c_1$ and $c_2$ here.

The CACEs for $c_i$ (i = 1, 2) with convection included are given by [19],

$${\partial }_t{c}_i+\overrightarrow{u}·\nabla c_i=\nabla ·\left[{\mathrm{m}}_0\right(\nabla c_i-({\mathsf{\lambda }}_{\mathrm{i}}{\overrightarrow{n}}_i-{\overrightarrow{\beta }}_i)\left)\right],$$

(1)

where $\overrightarrow{u}$ is the fluid velocity, $m_0$ is the (constant) mobility, ${\lambda }_i$ is defined as $4c_i(1-c_i)/W$ with W being the interface thickness, ${\overrightarrow{n}}_i=\nabla c_i/(\left|\nabla c_i\right|+{10}^{-12})$ (${10}^{-12}$ is added to avoid division-by-zero) is the unit normal vector of the interface of fluid i, and ${\overrightarrow{\beta }}_i$ is a Lagrange multiplier to enforce the condition ${\sum }_{i=1}^3{c}_i=1$, ${\overrightarrow{\beta }}_i=c_i{\sum }_{j=1}^3{\lambda }_j{\overrightarrow{n}}_j$. As for incompressible fluids where the velocity field is divergence free, $\nabla ·\overrightarrow{u}=0$, one can rewrite the convective term as $\overrightarrow{u}·\nabla c_i=\nabla {·(\overrightarrow{u}c}_i)$. In Cartesian coordinates, the CACE for $c_i$ may be written as

$${\partial }_t{c}_i+{\partial }_{\alpha }\left(u_{\alpha }{c}_i\right)={\partial }_{\alpha }\left(m_0{\partial }_{\alpha }{c}_i\right)-m_0{\partial }_{\alpha }({\lambda }_i{n}_{i,\alpha }-{\beta }_{i,\alpha }),$$

(2)

and in cylindrical coordinates for axisymmetric problems, it reads [19],

$${\partial }_t{c}_i+{\partial }_{\alpha }\left(u_{\alpha }{c}_i\right)+\frac{1}{r}{c}_i{u}_r={\partial }_{\alpha }\left(m_0{\partial }_{\alpha }{c}_i\right)+\frac{m_0}{r}[{\partial }_r{c}_i-({\lambda }_i{n}_{i,r}-{\beta }_{i,r}\left)\right]-m_0{\partial }_{\alpha }({\lambda }_i{n}_{i,\alpha }-{\beta }_{i,\alpha }),$$

(3)

where $\mathrm{r}$ is the coordinate in the radial direction (the other coordinate z is in the axial direction). Note that Equation (3) can be rewritten in conservative form as,

$${\partial }_t\left(r{c}_i\right)+{\partial }_{\alpha }(r{u}_{\alpha }{c}_i+m_0{c}_i{\delta }_{\alpha r})={\partial }_{\alpha }\left[{m_0\partial }_r\right(r{c}_i)-m_{0}r({\lambda }_i{n}_{i,\alpha }-{\beta }_{i,\alpha }\left)\right].$$

(4)

We first present the axisymmetric LBEs using the multiple-relaxation-time (MRT) collision model [20] to solve the axisymmetric CACEs. Two sets of distribution functions (DFs), $g_l$ and $h_l$, are used for the volume fractions $c_1$ and $c_2$, respectively. Here we only give the details on $g_l$ for $c_1$ (the formulations on $h_l$ for $c_2$ are similar and omitted for brevity). Inspired by [19,20], we adopt the following MRT version of the LBEs on $g_l$ for the axisymmetric CACE for $c_1$,

$$g_l(\overrightarrow{x}+{\overrightarrow{e}}_l{\delta }_t,t+{\delta }_t)-g_l(\overrightarrow{x},t)={-\left(M^{-1}{S}_{g}M\right)}_{lm}\left[g_m\right(\overrightarrow{x},t)-g_m^{eq}(\overrightarrow{x},t\left)\right]+{\delta }_t{G}_l(\overrightarrow{x},t),$$

(5)

where ${\overrightarrow{e}}_l$ is the lattice velocity along the direction $l$ ($l=\mathrm{0,1},\cdots ,8$ for the D2Q9 velocity model adopted here), ${\delta }_t$ is the time step, and $M$ is a matrix that transforms the vector of DFs into a vector of moments [21]. The volume fraction $c_1$ is obtained from the DFs $g_l$ as

$$c_1=\frac{1}{r}{\sum }_l{g}_l$$

(6)

For brevity, the details of the MRT collision model, the boundary conditions near the solid wall, and the calculation of the spatial derivatives are given in Appendix A.

2.2. Axisymmetric LBEs for the NSEs

The axisymmetric LBEs for the NSEs are quite similar to those in [22] except for a few changes related to the extension from binary to ternary fluids. For completeness, it is briefly introduced as follows. In LBM, the incompressible NSEs for ternary fluid flows with variable fluid properties are usually approximated by [23],

$$\frac{\partial p}{\partial t}+\rho (c_1,c_2)c_s^2\nabla ·\overrightarrow{u}=0,$$

(7)

$$\rho (c_1,c_2)(\frac{\partial \overrightarrow{u}}{\partial t}+\overrightarrow{u}·\nabla \overrightarrow{u})=-\nabla p+{\overrightarrow{F}}_{ST}+\nabla ·\left[\eta \right(c_1,c_2\left)\right(\nabla \overrightarrow{u}+{(\nabla \overrightarrow{u})}^T\left)\right],$$

(8)

where $\mathrm{p}$ is the hydrodynamic pressure, $\rho (c_1,c_2)$ is the density linearly interpolated from the volume fractions $c_1$ and $c_2$ as $\rho \left(c_1,c_2\right)={\rho }_1{c}_1+{\rho }_2{c}_2+{\rho }_3\left(1-c_1-c_2\right)$(${\rho }_i$ is the density of the i-th fluid), and $\eta (c_1,c_2)=\rho (c_1,c_2)\upsilon (c_1,c_2)$ is the dynamic viscosity with $\mathsf{\upsilon }({\mathrm{c}}_1,{\mathrm{c}}_2)$ being the kinematic viscosity found from those of fluid $\mathrm{i}(\mathrm{i}=1,2,3)$ as ${\left[\upsilon \right(c_1,c_2\left)\right]}^{-1}=c_1{\upsilon }_1^{-1}+c_2{\upsilon }_2^{-1}+(1-c_1-c_2){\upsilon }_3^{-1}$ [23]. The dynamic viscosity of fluid $\mathrm{i}$ is ${\eta }_i={\rho }_i{\upsilon }_i$, and ${\overrightarrow{F}}_{ST}$ is the interfacial tension force; its expression is given in Appendix A.

In Cartesian coordinates, the LBEs using the single-relaxation-time (SRT) collision model to recover Equations (7) and (8) read [24],

$$f_l\left(\overrightarrow{x}+{\overrightarrow{e}}_l{\delta }_l,t+{\delta }_l\right)-f_l(\overrightarrow{x},t)=-\frac{1}{{\tau }_f}(f_l-f_l^{eq})+(1-\frac{1}{2{\tau }_f})({\overrightarrow{e}}_l-\overrightarrow{u})\cdot [\nabla \rho c_s^2({\Gamma }_l-{\Gamma }_l\left(0\right))+{\overrightarrow{F}}_{ST}{\Gamma }_l]{\delta }_t,$$

(9)

where $f_l$ and $f_l^{eq}$ are the DFs and equilibrium DFs for the hydrodynamic variables. $f_l^{eq}$ is given by,

$$f_l^{eq}=w_l[p+\rho c_s^2(\frac{1}{c_s^2}{e}_{l\alpha }{u}_{\alpha }+\frac{1}{2c_s^4}(e_{l\alpha }{e}_{l\beta }-c_s^2{\delta }_{\alpha \beta })u_{\alpha }{u}_{\beta }\left)\right],$$

(10)

and the relaxation parameter ${\tau }_f$ is related to the kinematic viscosity as $\upsilon =c_s^2({\tau }_f-0.5){\delta }_t$. ${\Gamma }_l$ is the dimensionless DFs given by: ${\Gamma }_l=w_l(1+\frac{1}{c_s^2}{e}_{l\alpha }{u}_{\alpha }+\frac{1}{2c_s^4}{e}_{l\alpha }{u}_{\alpha }{e}_{l\beta }{u}_{\beta }-\frac{1}{2c_s^2}{u}_{\alpha }{u}_{\alpha })$. The fluid pressure and momentum are calculated as,

$$p={\sum }_l{f}_l+\frac{1}{2}{\delta }_t(\overrightarrow{u}·\nabla \rho c_s^2),$$

(11)

$$\rho \overrightarrow{u}=\frac{1}{c_s^2}{\sum }_l{f}_l{\overrightarrow{e}}_l+\frac{1}{2}{\delta }_t{\overrightarrow{F}}_{ST},$$

(12)

To simulate axisymmetric three-phase flows, the equations must be properly changed to account for axisymmetric effects. Specifically, the axisymmetric formulation used here is modified from that in [22]. When azimuthal flows are absent, the target axisymmetric NSEs to be recovered by the axisymmetric LBEs are

$$\frac{1}{c_s^2}\frac{\partial p}{\partial t}+\frac{\partial }{{\partial }_{xa}}\left(\rho u_a\right)=-\frac{\rho u_r}{r},$$

(13)

$$\rho (\frac{\partial u_{\alpha }}{\partial t}+u_{\beta }\frac{\partial u_{\alpha }}{\partial x_{\beta }})=-\frac{\partial p}{\partial x_{\alpha }}+F_{ST,,\alpha }+\frac{\partial }{\partial x_{\beta }}\left[\eta \right(\frac{\partial u_{\alpha }}{\partial x_{\beta }}+\frac{\partial u_{\beta }}{\partial x_{\alpha }}\left)\right]+F_{axisym,\alpha }-\frac{1}{r}\rho u_r{u}_{\alpha },$$

(14)

where the additional force due to axisymmetric effects $F_{axisym,\alpha }$ is given by,

$$F_{axisym,\alpha }=\frac{\eta }{r}(\frac{\partial u_{\alpha }}{\partial r}+\frac{\partial u_r}{\partial x_{\alpha }}),$$

(15)

As compared with the expression in [25], Equation (14) is more simplified because a different form of interfacial tension force is used here.

The axisymmetric LBEs for Equations (13) and (14) read [25],

$$f_l\left(\overrightarrow{x}+{\overrightarrow{e}}_l{\delta }_t,t+{\delta }_t\right)-f_l(\overrightarrow{x},t)=-\frac{1}{{\tau }_f}(f_l-f_l^{eq})+(1-\frac{1}{2{\tau }_f})[S_l+S_l^{\prime }+S_l^{″}]{\delta }_t,$$

(16)

where $S_l=(\overrightarrow{e_l}-\overrightarrow{u})·[\nabla \rho c_s^2({\Gamma }_l-{\Gamma }_l\left(0\right))+{\overrightarrow{F}}_{ST}{\Gamma }_l]$ is a source term having the same form as in Equation (9), $S_l^{\prime }$ and $S_l^{″}$ are another two source terms added due to the axisymmetric effects [25],

$$S_l^{\prime }=w_l{c}_s^2(-\frac{\rho u_r}{r}),S_l^{″}=({\overrightarrow{e}}_l-\overrightarrow{u})·({\overrightarrow{F}}_{axisym}-\frac{\rho u_r\overrightarrow{u}}{r}){\Gamma }_l,$$

(17)

The fluid pressure and momentum in the axisymmetric formulation are computed from

$$p={\sum }_l{f}_l+\frac{1}{2}{\delta }_t(\overrightarrow{u}·\nabla \rho c_s^2-\rho c_s^2\frac{u_r}{r}),$$

(18)

$$\rho \overrightarrow{u}=\frac{1}{c_s^2}{\sum }_l{f}_l{\overrightarrow{e}}_l+\frac{1}{2}{\delta }_t({\overrightarrow{F}}_{ST}+{\overrightarrow{F}}_{axisym}-\frac{\rho u_r\overrightarrow{u}}{r}),$$

(19)

It is noted that the new pressure and velocity are computed in an iterative manner as the above two equations are coupled with each other. The velocity gradient terms in Equation (15) are evaluated by the 2nd-order central difference schemes. For conciseness, the MRT model [21] and the Lax–Wendroff streaming scheme used to suppress the chequerboard effect are described in detail in Appendix A.

3. Results and Discussions

This work mainly studies the wetting and spreading of compound droplets on different types of solid walls and focuses on whether droplet separation occurs during the wetting process. Commonly, compound droplets may take two different configurations in the absence of solid walls, namely, the core-shell type and the Janus type. In this work, we concentrate exclusively on Janus droplets. For small droplets with relatively large interfacial tension, the effect of gravity can be ignored. In this section, we first present the calculation of some important quantities and the physical model, and then verify our numerical method by checking the shape of compound droplets on different walls in equilibrium state. After that, we will discuss the influence of several factors on the droplet separation phenomenon, including the wall type, curvature, contact angle, the interfacial angles, and the density ratio.

In phase field simulations of ternary fluid flows, the free energy functional $\mathcal{F}$ is defined by [26]

$$\mathcal{F}(c_1,c_2,c_3,\nabla c_1,\nabla c_2,\nabla c_3)={\int }_V\left(\frac{12}{W}{F}_0\right(c_1,c_2,c_3)+\frac{3}{4}W{\gamma }_1{\left|\nabla c_1\right|}^2+\frac{3}{4}W{\gamma }_2{\left|\nabla c_2\right|}^2+\frac{3}{4}W{\gamma }_3{\left|\nabla c_3\right|}^2)dV,$$

(20)

where the bulk free energy density $F_0$ is given by

$$F_0(c_1,c_2,c_3)={\gamma }_1{c}_1^2{(1-c_1)}^2+{\gamma }_2{c_2^2(1-c_2)}^2+{\gamma }_3{c}_3^2{(1-c_3)}^2,$$

(21)

and the latter three terms in Equation (20) represent the interfacial energies associated with $c_i(i=1,2,3)$. Note that an additional term in [26] for total spreading cases is omitted because here we only consider partial spreading cases. With Equation (21), one may separate $\mathcal{F}$ into three parts as $\mathcal{F}=\sum _{i=1}^3{\mathcal{F}}_i$ where ${\mathcal{F}}_i$ represents the contribution related to $c_i$,

$${\mathcal{F}}_i={\int }_V(\frac{12}{W}{\gamma }_i{c}_i^2{(1-c_i)}^2+\frac{3}{4}W{\gamma }_i{\left|\nabla c_i\right|}^2)dV=\frac{3{\gamma }_i}{4W}{\int }_V(16c_i^2{(1-c_i)}^2+W^2{\left|\nabla c_i\right|}^2)dV,$$

(22)

In the conservative phase field framework, the free energies mostly come from the interfacial regions because they almost vanish in the bulk regions where $c_i\approx 0$ or $c_i\approx 1$, $\left|\nabla {\mathrm{c}}_{\mathrm{i}}\right|\approx 0$. In actual simulations,${\mathcal{F}}_i$ is easily obtained from the fields for the volume fraction $c_i$ and its gradient. However, for specific problems, one is usually more interested in the interfacial energies associated with an interface, e.g., $E_{s,ij}$ between fluids $\mathrm{i}$ and $\mathrm{j}$ (than in ${\mathcal{F}}_i$). Fortunately, $E_{ij}$ can be found from ${\mathcal{F}}_i$ as

$$E_{s,ij}=({\gamma }_i+{\gamma }_j)(\frac{{\mathcal{F}}_i}{{\gamma }_i}+\frac{{\mathcal{F}}_j}{{\gamma }_j}-\frac{{\mathcal{F}}_k}{{\gamma }_k}),$$

(23)

where $i,j$, and $k$ are different. The total kinetic energy and the kinetic energy of fluid i may be calculated as

$$E_{k,t}={\int }_V\frac{1}{2}\rho (c_1,c_2){\left|\overrightarrow{u}\right|}^{2}dV,$$

(24)

$$E_{k,i}={\int }_{V}N\left(c_i\right)\frac{1}{2}{\rho }_i{\left|\overrightarrow{u}\right|}^{2}dV,$$

(25)

where $N\left(c_i\right)=1$ for $c_i>0.5$ and $N\left(c_i\right)=0$ otherwise. Note that in cylindrical coordinates for axisymmetric problems, ${\int }_V\left(\right)dV={\int }_V\left(\right)\left(2\pi r\right)dzdr$.

3.1. Physical Model and Parameters

Several kinds of solid walls with different geometries will be considered. To illustrate the basic physical settings, we choose the wetting of compound droplets on an ellipsoidal wall (corresponding to a prolate spheroid in three dimensions) as an example. As shown in Figure 1, the initial compound droplet of the Janus type is composed of two droplets of different components (fluid 1 on the left and fluid 2 on the right) that have been fused together to be in static equilibrium. The volumes of the two droplets are both equal to 4.189 (about the volume of a sphere with a radius of $R_0=1$). In other words, we chose the reference length as $L_r=R_0$. Note that the initial state may be also realized by placing two spherical droplets of radius 1, one of fluid 1 and the other of fluid 2, in contact and letting the system evolve freely under the action of interfacial tension (with dissipation) to equilibrium. To circumvent this step, we employed the analytical solution for the final equilibrium configuration under given interfacial angles. The surrounding environment is fluid 3. The angles ${\phi }_i(i=1,2,3)$ at the triple junction measured in fluid $\mathrm{i}$ in Figure 1 are called interfacial angles and, for cases other than total spreading, they satisfy

$${\sum }_{i=1}^3{\phi }_i=2\pi ,sin{\phi }_3/{\sigma }_{12}=sin{\phi }_2/{\sigma }_{13}=sin{\phi }_1/{\sigma }_{23}.$$

(26)

The domain used in our simulations is a rectangle with a size of $L_x\times L_y.$ The semi-axis lengths in the x-direction, y-direction, and z-direction of the prolate spheroid are $l_x=1.310$, $l_y=l_z=0.8736$ (in order to meet the axisymmetric requirement, it is necessary to have $l_y=l_z$), and its volume of the ellipsoid is $V=4\pi l_x{l}_y{l}_z/3=4.189$ (the same as the volume of a unit sphere). For conciseness, in the following only $l_y$ is given for a prolate or oblate spheroid and $l_z=l_y$ will be omitted. Note that the x- and y-axes in Figure 1 correspond to the z- and r-axes in the cylindrical coordinate system, respectively. Associated with the three kinds of interfaces in a ternary fluid system, there are three different contact angles (CAs) on a wall and these CAs are not independent. From the force balance of a compound droplet on a wall in equilibrium, they must satisfy the following condition [15]

$${sin\phi }_{2}cos{\theta }_{13}-{sin\phi }_{3}cos{\theta }_{12}-{sin\phi }_{1}cos{\theta }_{23}=0.$$

(27)

where ${\theta }_{ij}$ represents the contact angle for the interface between fluids $\mathrm{i}$ and $\mathrm{j}$ on the wall measured in fluid $\mathrm{i}$ (${\theta }_{ij}+{\theta }_{ji}=\pi$). Periodic boundary conditions are used on the left and right sides of the domain. Solid wall boundary conditions are used on the top side, and symmetrical boundary conditions are used on the bottom side. For simplicity, the density ratio and the kinematic viscosity ratio of fluid 1 to fluid 2 are fixed at $r_{\rho 12}=1$ and $r_{\upsilon 12}=1$ (note that the present method can handle other values of ${\mathrm{r}}_{\mathsf{\rho }12}$ and ${\mathrm{r}}_{\mathsf{\upsilon }12}$). The density ratio and the kinematic viscosity ratio of fluid 1 to fluid 3 are $r_{\rho 13}$ and $r_{\upsilon 13}$, respectively.

It is noted that the present setting in Figure 1 is similar to “Configuration L (lens)” in [11], which also considered axisymmetric problems. However, they mainly focused on the equilibrium states of compound droplets on a planar wall at some fixed interfacial angles. In contrast, we will study the dynamic wetting processes on different curved walls at a variety of conditions, during which the compound droplet may undergo significant topological changes. In addition, [13,16] also considered a Janus droplet near a wall. However, the arrangement of the droplets with respect to the wall is different from the present one. In their studies, both constituent droplets wetted the wall. The problem studied by them was two dimensional and does not actually exist in reality.

The reference velocity is chosen as $U_r={\sigma }_{12}/{\eta }_1$, which is derived from the interfacial tension ${\sigma }_{12}$ and the dynamic viscosity of fluid 1 ${\eta }_1(={\rho }_1{\nu }_1)$, and the reference time is derived from $L_r$ and $U_r$ as $T_r=\frac{L_r}{U_r}=\frac{L_r{\eta }_1}{{\sigma }_{12}}$. All length and time quantities are measured in $L_r$ and $T_r$, respectively. From these, one can calculate the Reynolds number as $Re=\frac{{\rho }_r{U}_r{L}_r}{{\eta }_r}=\frac{{\rho }_1{\sigma }_{12}{L}_r}{{\eta }_1^2}$ and the Weber number as $We=\frac{{\rho }_r{U}_r^2{L}_r}{{\sigma }_{12}}=Re$. The capillary number is found to be $Ca=\frac{{\eta }_r{U}_r}{{\sigma }_{12}}=\frac{We}{Re}=1$ and the Ohnesorge number is $Oh=\frac{{\eta }_1}{\sqrt{{\rho }_1{\sigma }_{12}{L}_r}}=\frac{1}{\sqrt{Re}}$. In the simulations, $L_r$ was discretized by $N_L$ grids and $T_r$ was discretized by $N_t$ time steps. Then, the grid size and time step are obtained as: ${\delta }_x=L_r/N_L$, ${\delta }_t=T_r/N_t$. In phase-field simulations, the interface thickness and mobility must be properly set. Regarding the former, $W/{\delta }_x=$ 5 was used as in many phase-field simulations in the literature and $N_L$ was set to be larger enough to make the Cahn number $C_n=W/L_r=(W/{\delta }_x)/N_L$ small. For the latter, we used the mobility $m_0=0.1$ (in lattice units) empirically (which provides reasonable results under most circumstances). In order to ensure the accuracy of the simulation, tests on the grid density and domain size were carried out for a typical case (see Appendix B and Appendix C for details). Based on the grid density tests, to balance the requirements of accuracy and computational overhead, $N_L=50$ was used in the following.

3.2. Numerical Validation

First, the equilibrium shapes of compound droplets on several types of solid walls are calculated by numerical simulations and compared with the corresponding analytical solutions based on the equilibrium conditions (see Appendix D and Appendix E for the derivation of the analytical solutions of the equilibrium state of the compound droplet on different walls). The total simulation time in each case is $t=600$. At this time, the total kinetic energy of the compound droplet has approached zero. The comparisons between the numerical and analytical results are shown in Figure 2, where the blue dashed lines are the interfaces obtained by our simulations and the red solid lines are from the corresponding analytical solutions. The common parameters used in the simulations are: ${\theta }_{13}=60°$, $r_{\rho 12}=1$, $r_{\rho 13}=825$, $r_{\upsilon 12}=1$, $r_{\upsilon 13}=0.06$, $Re=400$ ($Oh=0.05$). Note that the other two CAs in Equation (27) are not important here because only the left droplet wets the wall. The interfacial angles of Figure 2a,d are:${\phi }_1=150°$, ${\phi }_2=130°$, ${\phi }_3=80°$. Interfacial angles of Figure 2b,c are:${\phi }_1=155°$, ${\phi }_2=135°$, ${\phi }_3=70°$. The solid wall in Figure 2a is concave (i.e., having negative curvature) and its radius of curvature is $R_c=4$. In other words, the droplet wets the inner surface of a sphere. The solid walls in Figure 2b–d are all convex (i.e., the droplet wets the outer surface of an ellipsoid or a sphere). In Figure 2b the prolate spheroid has semi-axis lengths in the x-direction and y-direction $l_x=1.310$, $l_y=0.8736$. In Figure 2c, the oblate spheroid has semi-axis lengths $l_x=0.6$, $l_y=1.291$. In Figure 2d, the sphere has a radius of $R=1$. From Figure 2, it is observed that the equilibrium interface positions by our simulation are close to the theoretical solutions for all cases. Thus, our numerical method should be reliable to predict the wetting characteristics of Janus droplets on curved walls. It can also be seen from Figure 2 that as the curvature of solid wall varies, the shape of droplet 2 changes slightly whereas the shape of droplet 1 changes more significantly. When the wall curvature decreases with all other conditions fixed, the contact area between droplet 1 and the solid wall increases as found through the comparison between Figure 2a,d (or through the comparison between Figure 2b,c).

3.3. Wetting of Compound Droplets on Three Types of Solid Walls

In this section, three types of solid walls are considered: a concave spherical surface with a radius of curvature $R_c=4$, a plane, and a prolate spheroid with the semi-axis lengths $l_x=1.310$, $l_y=0.8736.$ The purpose is to explore the effect of wall shape on the dynamic behavior of a Janus droplet during its spreading on the wall. Figure 3 shows the wetting processes on the three kinds of walls. The common parameters used in the simulations are:${\phi }_1=150°$, ${\phi }_2=130°$, ${\phi }_3=80°$, ${\theta }_{13}=60°$, $r_{\rho 12}=1$, $r_{\rho 13}=825$, $r_{\upsilon 12}=1$, $r_{\upsilon 13}=0.06$, $Re=400$ ($Oh=0.05$). It is noted that the Reynolds number and the physical properties of the fluids (including the density and viscosity ratios) in this section are also used in all subsequent sections (unless specified otherwise).

It can be seen from Figure 3 that the motions of the compound droplet for all cases follow a similar sequence. First, the left droplet spreads on the wall. This initiates a capillary wave propagating from the contact line along the interface between fluids 1 and 3. After some time, the wave reaches the three-phase point (where the three fluids meet) and continues in two directions: one along the interface between the two drops (i.e., fluids 1 and 2) and the other along that between fluids 2 and 3. During its propagation, the capillary wave is damped to some extent due to viscous dissipation. It also takes some time for it to arrive at the three-phase point and subsequently affect the right droplet. As the right droplet has no direct contact with the wall, it is only affected by the interfacial tension of the interface between the two droplets (for convenience, called “fusion interface” below). Thus, its deformation and motion are lagging the left. When the intrinsic contact angle is the same, the concave wall “bends” the interface between fluid 1 and 3 more heavily than the planar and convex walls because the initial shape of the left droplet deviates from its (imagined) equilibrium configuration on the concave wall most severely (in other words, for the concave wall case the system initially has the largest interfacial energy potential). This causes the left droplet to spread the fastest and deform the most violently on the concave wall. During this process, the left droplet disconnects with the right droplet the earliest (see the snapshot at $t=50$ in Figure 3a). For the planar wall case, the initial potential to drive the left droplet is not as large as the concave wall case but is still enough to split the Janus droplet at a later time (see the snapshot at $t=90$ in Figure 3b). For the convex wall, the driving potential is the smallest and is insufficient to split the compound droplet (as see in Figure 3c). From Figure 3, it is also observed that, in the end, the two droplets remain in the form of a compound droplet for all three cases. On both the concave and planar walls, after the separation occurs, the left droplet slows down whereas the right droplet keeps moving towards the left because it was accelerated to obtain a certain momentum by the interfacial tension force from the left droplet. After some time, the right droplet touches the left one again and they merge to form a different Janus droplet on the wall that eventually reaches static equilibrium. Thus, the observed splitting of Janus droplet is only transient. Making use of this splitting stage to permanently separate the two fluids is beyond the scope of this work. Here we simply focus on the events of topological changes during the wetting process, which are of interest on their own.

Here it is helpful to examine the changes in the energy of the whole system (including the wall) between the initial state and the final equilibrium state for the three cases. In these two states, the velocity is zero everywhere, and therefore the system energy only consists of the interfacial energies. Denote the total area of the solid surface as $A_s$ and the area of the surface wetted by the left droplet (of fluid 1) as $A_w$. The area of the interface between fluids $\mathrm{i}$ and $\mathrm{j}$ is $A_{ij}$. The interfacial tension between the solid wall and fluid $\mathrm{i}$ is ${\sigma }_{is}$. In the initial state (labelled with a superscript ${ }^0$), the system energy may be expressed as $E_t^0={\sigma }_{12}{A}_{12}^0+{\sigma }_{13}{A}_{13}^0+{\sigma }_{23}{A}_{23}^0+{\sigma }_{3s}{A}_s$. In the final state (labelled with a superscript ${ }^{\mathrm{f}}$), it becomes $E_t^f={\sigma }_{12}{A}_{12}^f+{\sigma }_{13}{A}_{13}^f+{\sigma }_{23}{A}_{23}^f+{\sigma }_{1s}{A}_w+{\sigma }_{3s}{(A}_s-A_w)$. Thus, the energy change from the initial state to the final state is $\mathsf{\Delta }{E}_t^{0\to f}=E_t^0-E_t^f={\sigma }_{12}\left(A_{12}^0-A_{12}^f\right)+{\sigma }_{13}{(A}_{13}^0-A_{13}^f)+{\sigma }_{23}{(A}_{23}^0-A_{23}^f)+({\sigma }_{3s}-{\sigma }_{1s})A_w$. From Young’s equation on the contact angle, one has ${\sigma }_{3s}={\sigma }_{1s}+{\sigma }_{13}cos{\theta }_{13}$. Then, one has $∆E_t^{0\to f}={\sigma }_{12}∆A_{12}+{\sigma }_{13}∆A_{13}+{\sigma }_{23}∆A_{23}+{\sigma }_{13}cos{\theta }_{13}{A}_w$ where $∆A_{ij}=A_{ij}^0-A_{ij}^f$. Without loss of generality, we consider the energy change scaled by the interfacial tension between fluids 1 and 2, $\mathsf{\Delta }{E}_t^{0\to f}/{\sigma }_{12}$. The initial areas $A_{ij}^0$ are all given. The final areas $A_{ij}^f$ and the wetted area $A_w$ are obtained from the analytical solutions. In the end, the energy changes are found to be $∆E_t^{0\to f}/{\sigma }_{12}=2.887$ for the case on the concave wall, $∆E_t^{0\to f}/{\sigma }_{12}=2.864$ for the planar wall, and $∆E_t^{0\to f}/{\sigma }_{12}=1.792$ for the ellipsoidal wall. These results provide quantitative evidence for the above discussions on the driving potential.

For this problem, we monitored the centroid velocities and positions of the droplets along the x-direction $u_i$ and $x_i$ ($\mathrm{i}=1$ for the left droplet and $\mathrm{i}=2$ for the right droplet). Take the left droplet as an example. Its centroid velocity $u_1$ was calculated by,

$$u_1=\frac{{\int }_{A\left|c_1>0.5\right.}yu(x,y)dxdy}{{\int }_{A\left|c_1>0.5\right.}ydxdy},$$

(28)

where $A\left|c_1>0.5\right.$ represents the region where $c_1>0.5$. Figure 4 shows the evolutions of the relative velocity and position between the left and right droplets’ centroids along the x-axis, $u_1-u_2$ and $x_1-x_2$, respectively, for the three cases in Figure 3. It can be seen from Figure 4a that in the early stage ($t<80$), the relative velocities are negative for all cases, and its magnitude first increases and then decreases. During this stage, the maximum magnitude of the relative velocity (abbreviated as MMRV for brevity below) for the third case (the ellipsoidal wall) is the largest, the MMRV for the planar wall is the smallest, and that for the concave wall is in between. Interestingly, although the MMRV is larger on the concave wall than that on the ellipsoidal wall, the maximum distance between the two droplets’ centroids is smaller on the ellipsoidal wall during the early stage, as seen from Figure 4b. This is because the change in distance is determined by the integration of the relatively velocity in time, not by the MMRV. It can also be seen in Figure 4b that the maximum distance is the largest for the flat wall. From Figure 3 and Figure 4, it seems that whether the Janus droplet splits is more related with the MMRV than with the maximum distance between the two droplets.

For convenience, in subsequent discussions we set the flat wall as a baseline. The concave wall extends towards the upper and right side (in regions not far away from the axis, predominantly towards the upper side for $R_c=4$) on the right of the baseline. This limits the leftwards motion of the left droplet to some extent. At the same time, the interfacial tension force near the contact line gives the most violent pull to the left droplet towards the upper and right side, somehow tearing the left droplet off after a certain time. As we consider axisymmetric problems here, the upper direction is the radial direction, and the contact line corresponds to a circle in three dimensions. When the contact line moves further away from the axis, the perimeter of the circle (along which the tugging force acts) increases (proportional to the radial coordinate). In contrast, the ellipsoidal wall extends towards the left and upper side. For the specific prolate spheroid considered here, the direction is predominantly towards the left. Therefore, the increase in the perimeter of the contact line during wetting is less significant than that in the first case. In addition, the tugging force on the contact line is smaller in the third case. The differences in these two factors partially explain why the Janus droplet splits in the first case but not in the third. Intuitively, one can imagine two extreme scenarios (corresponding to two types of mechanisms) for the separation to occur. The first is the “relative motion induced separation”, in which the left droplet is accelerated leftwards extremely fast and the right droplet almost stays in its original place due to inertia. The second is the “deformation induced separation”, in which the left droplet deforms (to become relatively flat) very quickly and breaks the connection with the right droplet. Of course, in reality, both mechanisms may play some role concurrently. For the above three cases, the second mechanism seems to be more effective than the first one; as seen in Figure 3, the left droplet experiences significant deformations on the flat and the concave walls and separation occurs in these two cases. The situation on the convex wall is the opposite.

Figure 5 shows the evolution of the interfacial energy of the fusion interface between the two droplets with time for the three cases on different walls. The expression of this interfacial energy is given by Equation (23) with $i=1$, $j=2$. Because, during the period of droplet separation, the fuse interface disappears and the interfacial energy $E_{s,12}$ becomes zero, it is straightforward to determine from Figure 5 when the compound droplet splits and at what time the left and right droplets contact each other again. As seen in Figure 5, initially the interfacial energy decreases the fastest on the concave wall and the slowest on the ellipsoid. The change of $E_{s,12}$ on the planar wall looks close to that on the concave wall in the initial stage, during which $E_{s,12}$ vanishes for a certain time. In contrast, $E_{s,12}$ always remains positive on the ellipsoid, indicating that droplet separation never occurs. These results are consistent with the previous findings. Through the comparison of the interfacial energy in the late stage when the system approaches equilibrium on the three kinds of solid walls, it can also be found that $E_{s,12}$ on the concave wall is the largest and that on the ellipsoid is the smallest.

The evolution of the kinetic energy of the two droplets was also monitored, as shown in Figure 6. From this figure, it is found that the maximum kinetic energy of the left droplet (observed during the initial stage when it wets the surface) on the ellipsoid is much smaller than those on the concave and planar surfaces. This is not only because the attraction of the ellipsoidal wall to the droplet is the weakest (due to its particular shape), but also because the interaction between the two droplets on the ellipsoid has always been relatively large (as no separation occurs), and the left droplet is retarded by the right droplet the most. On the other hand, the maximum kinetic energy of the right droplet (observed after the initial stage) on the ellipsoid is the largest and occurs the earliest among all three cases. The reason is as follows. Because the compound droplet on the ellipsoid does not split, the right droplet is always in an accelerated state in the early stage of the wetting process. In contrast, for the planar and concave walls, the right droplet accelerates first, then decelerates after separation, and then accelerates again after the two droplets reconnect. For all three walls, the kinetic energies of both droplets become very low after a long time (e.g., $\mathrm{t}>400$), as seen in Figure 6.

As seen from the above, it is more difficult to make the compound droplet separate on an ellipsoidal wall than on a concave wall or a planar wall. Next, we further explore the wetting of a Janus droplet on ellipsoidal walls with different curvatures. Three different ellipsoids with the same volume (4.189) were studied. The first is an oblate spheroid with $l_x=0.6$, $l_y=1.291$. The second is a unit sphere. The third is the prolate spheroid just studied previously. The interfacial angles are slightly different from those in the above: ${\phi }_1=155°$, ${\phi }_2=135°$, ${\phi }_3=70°$. The other parameters are the same as those in Figure 3c. The flow field evolutions for the three cases are shown in Figure 7.

As seen in Figure 7, droplet separation only occurs near the oblate spheroid. During its spreading on the unit sphere and on the prolate spheroid, the left droplet is always connected with the right one. On the unit sphere, the interfacial area between the two droplets shrinks to a quite small value, from $\mathrm{t}=70$ to $\mathrm{t}=90$, approaching the critical state of droplet separation. However, at that critical moment, the distance between the two droplets reaches its maximum and the deformations of the interfaces also reach the largest degree. Thus, there is no further bending or distortion of the interface between the two droplets to make them disconnect. The situation on the prolate spheroid is similar to that in the above. These observations further confirm the importance of the second mechanism (the deformation induced separation). The oblate spheroid extends away from the axis the most among all three solid objects of the same volume and provides the largest area for the left droplet to spread upwards, leading to sufficiently large distortion of the left droplet in a short enough time. These are key to breaking up the Janus droplet.

3.4. Effects of the Radius of Curvature and Contact Angle of the Wall on Droplet Separation

This section explores the effects of the radius of curvature of the wall together with its wettability on the separation behavior of compound droplets over a wider parameter range. Both an ellipsoidal wall with variable positive curvature and a concave wall with a constant negative curvature were considered. The ellipsoid volume is fixed at 4.189. Several lengths of the semi-axis in the x-direction $l_x=0.7$, $0.8$, $0.9$, $1$, $1.1$, $1.2$, $1.3$ were tested. The corresponding values of $l_y$ are $l_y=1.195$, $1.118$, $1.054$, $1$, $0.9535$, $0.9129$, and $0.8771$, respectively. The contact angle ${\theta }_{13}$ varies from $50°$ to $75°$. Other parameters are the same as those used in Figure 7 in Section 3.3. Figure 8 shows the phase diagram of the separation state on ellipsoids with different $l_x$ and ${\theta }_{13}$. In this figure, the cross symbol means no separation occurs and the filled circle means the contrary. With the decrease of $l_x$ (concurrently the increase of $l_y$), the radius of curvature of the ellipsoid near the axis ($r_c$) increases. The corresponding values of $r_c$ are ${\mathrm{r}}_{\mathrm{c}}=2.040$, $1.562$, $1.234$, $1$, $0.8265$, $0.6945$, and $0.5918$, respectively. From Figure 8 it can be seen that smaller contact angles are conducive to droplet separation. This is because more hydrophilic walls have a stronger attraction force on the left droplet, making its acceleration and deformation during the wetting process larger. In addition, at each contact angle, to decrease $l_x$ facilitates the occurrence of droplet separation. As the contact angle increases, the critical $l_x$ for droplet separation decreases. It means that to observe droplet separation during the spreading on a less hydrophilic ellipsoid, one must increase the radius of curvature. All these agree with previous findings and analyses.

In addition to ellipsoidal walls, we also varied the radius of curvature for a concave wall having a constant curvature. Based on previous results, droplet separation is more likely to occur on a concave wall. To encompass a broader range of parameters, the interfacial angles in this part are chosen as ${\phi }_1=150°$, ${\phi }_2=130°$, and ${\phi }_3=80°$. The contact angle varies from $50°$ to $75°$. The radius of curvature changes from 1.32 to 64. Other parameters are not altered. The phase diagram of the separation state on a concave wall under different combinations of $R_c$ and ${\theta }_{13}$ are shown in Figure 9. Note that the increments in $R_c$ are nonuniform and the horizontal coordinates are not to scale. As found from this figure, when the wall is not so hydrophilic (e.g., ${\theta }_{13}=65°$ and $70°$) and the radius of curvature is large (e.g., $R_c>2$), one can see a transition from “no separation” to “separation” by properly reducing the radius of curvature. However, this effect of reducing ${\mathrm{R}}_{\mathrm{c}}$ fails to work when $R_c\le 1.4$. At relatively small $R_c$, the volume enclosed by the concave wall is also small, and there is not enough space to accommodate the left droplet spreading outwards (i.e., away from the axis). When the left droplet does not have enough deformation, its connection with the right droplet is hard to break. Therefore, if one wants to change the radius of curvature of a concave wall to promote the separation of Janus droplets, a moderate value for $R_c$ should be taken.

3.5. Effects of the Interfacial Angles and Wall Contact Angle on Droplet Separation

In addition to the geometrical factors and wettability of the wall, the interfacial angles may also influence the droplet separation behavior. Figure 10a,b shows two typical wetting and spreading processes on a prolate spheroid with ${\phi }_1=155°$, ${\phi }_2=135°$, ${\phi }_3=70°$ and ${\phi }_1=135°$, ${\phi }_2=115°$, ${\phi }_3=110°$. The contact angle ${\theta }_{13}$ is fixed at ${\theta }_{13}=50°$ for these two cases. In Figure 10c, the contact angle ${\theta }_{13}$ is ${\theta }_{13}=60°$, and the interfacial angles are the same as in Figure 10a. All other parameters are the same as in Figure 3c. As found in Figure 10a, when ${\phi }_1=155°$, ${\phi }_2=135°$, ${\phi }_3=70°$, the two droplets separated for a period of time. In Figure 10b, when ${\phi }_1=135°$, ${\phi }_2=115°$, ${\phi }_3=110°$, the two droplets are always tightly connected throughout the wetting process. The significance of the interfacial angles may be demonstrated through the energy difference between two static states: one is the initial state and the other is an imaginary state (labelled with a superscript ${ }^{\mathrm{s}\mathrm{e}\mathrm{p}}$) in which the left and right droplets are two separate free spheres. For the time being, we assume that no wall is present. Then, the system energy for the former state is $E_t^0={\sigma }_{12}{A}_{12}^0+{\sigma }_{13}{A}_{13}^0+{\sigma }_{23}{A}_{23}^0$ and that for the latter is $E_t^{sep}={\sigma }_{13}{A}_{13}^{sep}+{\sigma }_{23}{A}_{23}^{sep}$ with $A_{13}^{sep}=A_{23}^{sep}=4\pi R_0^2$. The scaled energy difference is found to be $\mathsf{\Delta }{E}_t^{0\to sep}/{\sigma }_{12}=-0.1091$ for ${\phi }_1=155°$, ${\phi }_2=135°$, ${\phi }_3=70°$ and $\mathsf{\Delta }{E}_t^{0\to sep}/{\sigma }_{12}=-0.8599$ for ${\phi }_1=135°$, ${\phi }_2=115°$, ${\phi }_3=110°$. From the energy perspective, the minimum energy required to split the Janus droplet for ${\phi }_1=155°$, ${\phi }_2=135°$, ${\phi }_3=70°$ is much smaller than that for ${\phi }_1=135°$, ${\phi }_2=115°$, ${\phi }_3=110°$. This explains the differences between Figure 10a,b. By comparing Figure 10a,c, it can be found that the wetting behavior at ${\mathsf{\theta }}_{13}=50°$ is obviously different from that under the same interfacial angles at ${\mathsf{\theta }}_{13}=60°$. This is due to the more hydrophilic wall’s stronger attraction, as stated before.

To know the effects of the interfacial angles and wall wettability on droplet separation in a broader parameter regime, a number of simulations on compound droplet wetting on the above prolate spheroid under different interfacial angles and contact angles were carried out. The interfacial angle ${\phi }_1$ varies from $150°$ to $164°$, and the contact angle ${\theta }_{13}$ varies from $45°$ to $70°$. Note that the difference between the first two interfacial angles (${\phi }_1-{\phi }_2$) is fixed at $20°$. Other parameters are the same as those given before. Based on the simulation outcome, a phase diagram on whether droplet separation occurs during the wetting process is established in Figure 11. As shown in Figure 11, in general when the interfacial angle ${\phi }_1$ increases, the critical contact angle for separation increases. For instance, at ${\phi }_1=154°$, the critical contact angle is between $50°$ and $55°$ and at ${\phi }_1=158°$ it is between $60°$ and $65°$. These results indicate that, at a larger interfacial angle ${\phi }_1$, the bond between the left and right droplets weakens and may be broken by the wetting on a less hydrophilic wall. Besides, the different variation ranges of the interfacial angle and contact angle in Figure 11 suggest that the interfacial angle has greater impact on the droplet separation phenomenon than the contact angle. While the separation phenomenon depends on both the intrinsic properties of both a Janus droplet, like the interfacial angles, and external factors, like the wall wettability, it seems that the former plays a more influential role than the latter.

In addition to a prolate spheroid, we also performed similar investigations for a concave wall with a radius of curvature ${\mathrm{R}}_{\mathrm{c}}=4$ and generated a phase diagram in the ${\phi }_1-{\theta }_{13}$ plane under the same physical conditions, as shown in Figure 12. Like Figure 11, the range of variation for the contact angle is also from $45°$ to $70°$. As found from Figure 11 and Figure 12, when other parameters are equal, the interfacial angle ${\phi }_1$ required for the separation of compound droplets on the concave wall is smaller than that on the prolate spheroid. As mentioned above, a smaller interfacial angle ${\phi }_1$ corresponds to a stronger bond between the two droplets. This indicates that under the same wall wettability, the concave surface drags and deforms the left droplet more severely than the convex surface (such as the ellipsoid). This agrees with previous analyses in Section 3.3. By comparing the variation ranges of the interfacial angle in Figure 11 and Figure 12, it can be seen that the influence of the interfacial angle on droplet separation is even more significant on the concave surface than on the ellipsoid. For example, on the concave surface at ${\phi }_1=146°$, the critical contact angle for droplet separation is between $50°$ and $55°$ and at ${\phi }_1=148°$ it is between $60°$ and $65°$. Roughly speaking, to have a $~10°$ change in the critical contact angle, it takes only a ~$2°$ change in ${\phi }_1$ on the concave surface, whereas it requires a ~$4°$ change in ${\phi }_1$ on the ellipsoid. These findings indicate that the influence of the interfacial angle ${\phi }_1$ is amplified to some extent by the concave surface. In other words, the geometry of the wall not only affects whether droplet separation happens under the same wettability, but also affects the sensitivity of this phenomenon to the interfacial angle ${\phi }_1$. The stronger sensitivity to the interfacial angle ${\phi }_1$ on the concave wall can be understood from the following perspective. Because in the early state the left droplet experiences more extensive spreading and larger deformation on the concave wall, the interfacial area between the two droplets shrinks quickly to a low value (e.g., see Figure 3a in Section 3.3). Under such circumstances, even a seemingly minor change in the interfacial angle ${\phi }_1$ can determine whether it will further shrink to zero or not.

The above findings from Figure 12 may also be understood from the other side. That is, given the same change in the interfacial angle, the influence of the contact angle is smaller on the concave surface than on the ellipsoidal wall. This is because the wetting effect of the concave solid wall itself on the droplets has been relatively good, and the influence of continuing to enhance its wetting effect on the separation of droplets will not be too great. On the other hand, it can be seen from Figure 3a in Section 3.3 that when the droplets on the concave surface are separated, droplet 1 has been wetted on the wall to a considerable extent. The possibility of separation mainly depends on whether the liquid bridge between two droplets will break when droplet 1 spreads to the maximum on the solid wall; this condition is closely related to the interfacial angle, so the interfacial angle is the most influential factor in the process of droplet separation on the concave surface whose ${\mathrm{R}}_{\mathrm{c}}=4$.

3.6. Effect of the Density Ratio of Droplet to Ambient Fluid on Droplet Separation

In addition to the wall properties and interfacial angles, the density ratio may also have certain impacts on the wetting behavior of a compound droplet. This section explores the influence of the density ratio of droplet to ambient fluid (${\mathrm{r}}_{\mathsf{\rho }13}$) on the separation behavior of Janus droplets. The interfacial angles are fixed at ${\phi }_1=155°$, ${\phi }_2=135°$, ${\phi }_3=70°$. The same prolate spheroid as above is used. In addition to ${\mathrm{r}}_{\mathsf{\rho }13}$, the contact angle ${\theta }_{13}$ is also varied. Other parameters are the same as those given in Section 3.3. Note that the kinematic viscosity ratio ${\mathrm{r}}_{\mathsf{\upsilon }13}$ is fixed at $0.06$. Thus, as the density ratio ${\mathrm{r}}_{\mathsf{\rho }13}$ changes, the dynamic viscosity of the ambient fluid varies. Based on the simulation results using different combinations of ${\mathrm{r}}_{\mathsf{\rho }13}$ and ${\theta }_{13}$, another phase diagram on droplet separation on an ellipsoidal wall is established in Figure 13. From this figure, one can see that, in general, a larger density ratio ${\mathrm{r}}_{\mathsf{\rho }13}$ is favorable for droplet separation at a given contact angle ${\theta }_{13}$. A larger density ratio ${\mathrm{r}}_{\mathsf{\rho }13}$ corresponds to a lighter ambient fluid. The dynamic viscosity of the ambient fluid is also smaller at a higher ${\mathrm{r}}_{\mathsf{\rho }13}$ (with ${\mathrm{r}}_{\mathsf{\upsilon }13}$ being the same). These two factors reduce the obstruction of the ambient fluid to the motion of the compound droplet, making droplet separation easier. From the results at ${\mathsf{\theta }}_{13}=55°$, it is observed that when the contact angle is not small enough, droplet separation will not occur even if the density ratio is very large.

4. Conclusions

To conclude, the wetting and spreading behavior of a Janus droplet on solid walls of various shapes and wettabilities was investigated by axisymmetric LBM simulations, mostly at high density ratios. The evolutions of the interfaces and interfacial energies were examined to determine whether the two constituent droplets separate during the wetting process. The main factors found to promote the occurrence of droplet separation include (a) a concave wall with moderate radius of curvature, (b) a convex wall with large radius of curvature, (c) a small contact angle for the droplet near the wall, (d) large interfacial angles in the droplets, and (e) a large density ratio of the droplet to ambient fluid.

Comparisons between the wetting processes on three types of walls suggest that the moderately outward extending concave wall has better apparent wettability and enhances the spreading speed of the wetting droplet, leading to its large deformation and large velocity difference between the two constituent droplets. In the study of solid walls of various spheroid shapes, it was found that the larger the radius of curvature near the wetting droplet is, the faster its contact line moves on the wall, and the two droplets are more likely to separate. The influence of curvature for concave walls is non-monotonic. Only when the curvature is small, droplet separation is more likely to occur with its increase. When the curvature increases to a certain extent, the small radius of curvature limits the spreading of the wetting droplet and prevents droplet separation. When the contact angle of the wetting droplet is smaller, it experiences stronger attraction forces towards the wall and the distance between the two droplets increases more rapidly. As the interfacial angles in the constituent droplets increase, less energy is needed to split the compound droplet. When the density of the ambient fluid decreases, the compound droplet is less hindered, and the wetting droplet is allowed to move faster during its spreading. The five phase diagrams on droplet separation state, based on numerous simulations at different conditions, clearly show the influence of various factors for different walls. These results and findings may be helpful to future studies of Janus droplets near solid walls and could also advance the understanding of the interaction between compound droplets and solid microparticles. In addition, the study of the equilibrium morphology of compound droplets on curved surfaces might find applications in microfabrication of certain small components with special shapes (e.g., microlens).

Finally, the present work has some limitations. First, this study only discussed the wetting and spreading of compound droplets on curved surfaces under axisymmetric settings. To capture the behavior of compound droplets in many real situations, full three-dimensional simulations may be necessary. Second, the droplet separation reported in this work is only transient. It requires further investigations on how to sustain the separation of two constituent fluids.