Buoyancy-driven attraction of active droplets

TL;DR

This paper investigates how buoyancy effects influence active oil droplet behavior, revealing that buoyancy can cause attraction and collisions contrary to the repulsive Marangoni flow, with implications for understanding droplet clustering.

Contribution

It provides a numerical analysis of buoyancy-driven attraction in active droplets, incorporating parameters like Peclet, Galileo, and Rayleigh numbers, and derives relationships governing droplet interactions.

Findings

Buoyancy induces attraction, leading to droplet collisions.

Marangoni flow causes droulsion, counteracted by buoyancy effects.

Transition between regimes is characterized by Pe ~ Ra^{0.63}.

Abstract

Active oil droplets in a liquid are believed to repel due to the Marangoni effect, while buoyancy effects caused by the density difference between the droplets, diffusing product, and ambient fluid are usually overlooked. Recent experiments have observed active droplet clustering phenomena due to buoyancy-driven convection (Kruger et al. Eur. Phys. J. E, vol. 39, 2016, pp.1-9). In this study, we numerically analyze the buoyancy effect in addition to Marangoni flow, characterized by Peclet number $Pe$. The buoyancy effects originate from (i) the density difference between the droplet and the ambient liquid, which is characterized by Galileo number $Ga$, and (ii) the density difference between the diffusing product (i.e. filled micelles) and the ambient liquid, characterized by a solutal Rayleigh number $Ra$. We analyze how the attracting and repulsing behavior depends on the control parameters $Pe$, $Ga$, and $Ra$. We find that while Marangoni flow causes repulsion, the buoyancy effect leads to attraction, and even collisions can take place at high Ra. We also observe a delayed collision as $Ga$ increases. Moreover, we derive that the attracting velocity, characterized by a Reynolds number $Re_d$, is proportional to $Ra^{1/4}/(l/R)$, where $l/R$ is the normalized distance by radius between neighboring droplets. Finally, we obtain repulsive velocity, characterized by $Re_{rep}$, as proportional to $PeRa^{-0.38}$. The balance of attractive and repulsive effects results in $Pe \sim Ra^{0.63}$, which agrees with the transition curve between regimes with and without collision.

Equations

$$D\frac{\partial \hat{c}}{\partial \hat{n}}=-\alpha ,$$

$${\hat{u}}_s=M{\nabla }_s\hat{c},$$

$$\hat{\rho }(\hat{c})={\hat{\rho }}_0(1+\beta \hat{c}),$$

$$\frac{\partial c}{\partial t}+\mathbf{u}\cdot \nabla c=\frac{1}{Pe}{\nabla }^{2}c,$$

$$\frac{d{\mathbf{U}}_d}{dt}=\frac{3Sc}{4\pi Pe}\int (\mathbf{\tau }\cdot \mathbf{n})dS-\frac{G{a}^{2}S{c}^2}{P{e}^2}{\mathbf{e}}_z,$$

$$Ra=\frac{c_0\beta R^{3}g}{\nu D}\text{,}$$

$$Pe=\frac{\alpha MR}{D^2},$$

$$Sc=\frac{\nu }{D},$$

$$Ga=\frac{\sqrt{|{\hat{\rho }}_d/{\hat{\rho }}_0-1|g{R}^3}}{\nu }.$$

$$\frac{\partial c}{\partial n}=-1.$$

$$u_s={\nabla }_{s}c.$$

$$\frac{\partial c}{\partial z}=0,$$

$$\mathbf{u}=0.$$

$${\ell }_{\infty }=\underset{t\to \infty }{\lim }\ell (t)=\underset{t\to \infty }{\lim }(\text{min}({\ell }_1(t),{\ell }_2(t))).$$

$$v_z\sim \frac{D}{R}R{a}^{1/2},$$

$$h/R\sim R{a}^{-1/4}.$$

$$v_y(r)\times 2\pi r{h}_2\sim v_z\times \pi h_1^2/4,$$

$$R{e}_y(r)=\frac{v_y(r)R}{\nu }\sim \frac{1}{Sc}\frac{R{a}^{1/4}}{r/R}.$$

$$\frac{R{e}_y(r)}{R{a}^{1/4}}\sim \frac{R}{r}-\frac{R}{L_y-r}.$$

$$U_1^0=U_2^0=0,$$

$$u_2^0\sim \frac{DR{a}^{1/4}}{\ell }$$

$$U_{2}^{1}=\Big{(}1+\frac{R^{2}}{6}\nabla^{2}\Big{)}u_{2}^{0}.$$

$$u_{2}^{1}=\Big{(}\frac{3R}{2\ell}-\frac{R^{3}}{2\ell^{3}}\Big{)}U_{1}^{1}.$$

$$U_{2}^{2}=\Big{(}1+\frac{R^{2}}{6}\nabla^{2}\Big{)}u_{2}^{1}.$$

$$U_{2}^{n}\sim O\Big{(}\Big{(}\frac{R}{\ell}\Big{)}^{n-1}\Big{)}.$$

$$U(\ell)=U_{2}=U_{2}^{0}+U_{2}^{1}+U_{2}^{2}+O\Big{(}\frac{R}{\ell}\Big{)}=u_{2}^{0}+O\Big{(}\frac{R}{\ell}\Big{)}\sim\frac{DRa^{1/4}}{\ell},$$

$$R{e}_d(\ell )=\frac{UR}{\nu }\sim \frac{D}{\nu }\frac{R{a}^{1/4}}{\ell /R}=\frac{1}{Sc}\frac{R{a}^{1/4}}{\ell /R}.$$

$$\frac{R{e}_d}{R{a}^{1/4}}\sim \frac{R}{\ell }-\frac{R}{L_y-\ell }.$$

$$\frac{R{e}_d}{R{a}^{1/4}}\sim \frac{R}{\ell }-\frac{R}{L_y-2\ell }.$$

$$U_{\text{rep}}=\frac{1}{4\pi R^2}{\int }_S{\text{u}}_{s}dS.$$

Full text

1 Introduction

The fundamental principles of microorganisms propulsion have gained attention across disciplines over the past few decades (Brennen & Winet, 1977; Stone & Samuel, 1996; Lauga & Powers, 2009; Marchetti et al., 2013; Li & Ardekani, 2016; Blackiston et al., 2021). Given the abundance of such microorganisms such as bacteria and plankton in our ecosystem (Hays et al., 2005), studying their individual and collective motion is critical for understanding the dynamics of the entire ecosystem (Guasto et al., 2012). The interactions between microorganisms can be purely physical, i.e. based on hydrodynamics (Ramia et al., 1993; Ishikawa et al., 2006), or biological, i.e. based on visual signals (Trushin, 2004) or by chemical signals (Adler, 1975). Disentangling these effects makes it difficult to analyze the various interactions in real microorganism colonies. To reduce the complexity, in recent years artificial microswimmers as a simplified model have been investigated in order to understand the interactions between living microorganisms (Pedley, 2016; Maass et al., 2016; Datt & Elfring, 2019; Hokmabad et al., 2019; Gompper et al., 2020; Chen et al., 2021; Li, 2022). Such artificial microswimmers are designed to propel themselves by converting free energy from the environment into kinetic energy (Ogrin et al., 2008). Similar interactions as those between living microorganisms are observed, such as chemotaxis, collective entrainment, and cluster formation (Maass et al., 2016; Lozano et al., 2016; Jin et al., 2017; Lohse & Zhang, 2020; Jin et al., 2021).

One extensively studied type of artificial microswimmer is a dissolving active oil droplet floating in water (Maass et al., 2016). The driving mechanism behind the propulsion of such active droplets is the Marangoni effect. The basic feature is that whenever there is an inhomogeneity of surfactant concentration at the surface of the droplet, the consequent surface tension difference generates a tangential Marangoni flow adjacent to the surface, which leads to the self-propulsion of the droplet (Herminghaus et al., 2014; Morozov & Michelin, 2019a, b; Michelin, 2022). This effect can also be generalised to other coupled systems such as the particles with catalytic surfaces. The resulting flows are then referred to as diffusio-phoretic flow (Anderson, 1989). With the Marangoni effect or diffusiophoresis be present, such active droplets become repulsive. This simply happens because the Marangoni flow or the diffusiophoretic flow will propel the active droplet or particles towards higher surfactant concentration direction (where the surface tension is lower), and the concentration of surfactant molecules is lower between two close-by droplets or particles than that in the periphery.

Repulsive interactions induced by Marangoni effects between active droplets have been well studied in numerous experimental and theoretical works. A clear experimental observation of the repulsion was conducted by Moerman et al. (2017), who quantitatively measured the repulsive velocity for a pair of active droplets and analyzed the relations between the repulsive force and their distance. Later on, Lippera et al. (2020) theoretically analyzed the repulsive interaction between a pair of droplets and identified different motion modes. In a further study, Lippera et al. (2021) investigated the repulsive interactions for a pair of obliquely-colliding droplets, and identified whether the droplets interact directly or through their chemical wake. Besides the direct interaction between active droplets, Jin et al. (2017) found the active droplets also show trail avoidance behavior. They reported that the active droplet emitted filled micelles in the wake play a role as chemical repellents and cause trajectory avoidance. Based on such observation, Daftari & Newhall (2022) developed a mathematical model to mimic the trail-avoidance with multiple active droplets. They observed that the active droplets could trap themselves due to the trail-avoidance, a feature which has been called transit self-caging behavior.

However, besides the Marangoni effect caused by the difference in surface tension, also buoyancy effect caused by the inhomogeneity in the density field can affect the flow, namely by natural convection. One interesting example of a droplet for which the interplay between Marangoni and buoyancy forces leads to rich dynamics is the phenomenon of the ”jumping droplet” (Li et al., 2019, 2021, 2022). In that case the droplet repeatedly jumps up due to the Marangoni effect, and then slowly sinks due to buoyancy (Li et al., 2019, 2022). Another example is the evaporation of a binary micro-droplet, where buoyancy competes with Marangoni forces and can drive the convection inside the droplet, playing a crucial role in the evaporation process (Edwards et al., 2018; Li et al., 2018, 2020; Diddens et al., 2021). Finally for a pair of fixed droplets, Lopez de la Cruz et al. (2022) reported an oscillatory flow near the droplet, which again is triggered by the competition between Marangoni and buoyancy effects.

Coming back to dissolving active droplet, also here, despite the repulsive interaction by the Marangoni effect, Krueger et al. (2016) observed opposing collective behaviours when buoyancy is significant in the collective droplet system. They reported that the active droplets attract each other and form clusters hovering in the fluid and ascribed this finding to buoyancy effect because the collective behavior only occurs when the density difference between the solvent and the droplet is above a certain threshold. In a further study, Hokmabad et al. (2022) investigated the spontaneous rotation of the cluster formed by the attracted droplets. However, an explanation of the detailed mechanism of the attraction is still missing, especially on how the buoyancy effect drives the attraction and overtakes the repulsion driven by the Marangoni effect.

Inspired by the above-mentioned studies, we focus on the collective behavior of active droplets with buoyancy effects. Note that there are two different types of buoyancy effects, either by the density difference between the droplet and ambient fluid, or by the solutal density difference between the dissolving product (i.e. filled micelles) and the ambient fluid. Hereafter, we call them the ”droplet buoyancy effect” and the ”product buoyancy effect”, respectively. In this work, we will quantitatively analyze the interplay between the Marangoni effect, the droplet buoyancy effect, and the product buoyancy effect. We first simulate the interaction between a pair of active droplets. The numerical simulations allow us to capture the flow field around the droplet and how it induces the droplet interaction. Then we develop a model to predict the attracting velocity based on the method of reflections and on Faxen’s law. We test our model with systems of two and three interacting droplets. Then we analyze the repulsive velocity (based on the Marangoni flow) from simulations of a pair of fixed droplets. Finally, by comparing these results with those from the attractive velocity model due to buoyancy effect, we obtain a good prediction for the regime transition of droplet collision.

The paper is organized as follows: We first describe the problem setup in Section 2. The numerical method and validations of the numerical scheme are provided in Section 3. We first qualitatively analyze the role of the diffusiophoretic effect (characterized by $Pe$) and the product buoyancy effect in Section 4 and then we analyze the role of the droplet buoyancy effect in Section 5. Next, we develop a model to explain the attraction and calculate the attracting velocity and how it scales in section 6.1 and 6.2 . The model is then further tested with the system of a pair of droplets and three droplets in section 6.3. Then the repulsive effect is analyzed with cases of a pair of fixed droplets in section 7. Finally, concluding remarks are given in section 8.

2 Setup and control parameters

We start with a pair of active droplets in the surfactant solution sketched in figure 1. The gradual solubilization of the oil into the surfactant micelles causes a repulsive interaction via the Marangoni effect (Jin et al., 2017). Simultaneously, the oil-filled micelles are generated near the droplet surface. Besides that, the droplet buoyancy effect and product buoyancy effect are taken into consideration in the simulations.

Considering similarities between diffusiophoresis and the Marangoni effect (Desai & Michelin, 2021), for simplicity, we focus on the phoretic effect induced by the concentration gradient of the filled micelles and will use the corresponding terminology. The physical variables to describe the system are the solutal concentration $\hat{c}$ and the velocity $\hat{\bf{u}}$. Note that all dimensional physical fields are marked with a hat (e.g. $\hat{c}$, $\hat{\bf{u}}$), while the dimensionless ones are without a hat (e.g. $c$, $\bf{u}$).

The droplets emit a solute (filled micelles) at a rate $\alpha>0$. The concentration boundary condition at the droplet surface is given by

[TABLE]

where $D$ is the diffusion coefficient of the dissolution product, $\alpha$ the dissolution rate at the surface, and $\partial\hat{c}/\partial\hat{n}$ the concentration gradient normal to the surface.

The tangential concentration gradient at the surface induces a slip velocity, which is the so-called diffusiophoretic flow. The magnitude of the slip velocity $u_{s}$ is proportional to the local tangential concentration gradient, given by

[TABLE]

where $M$ is the mobility and $\nabla_{s}$ represents the tangential gradient. Since the filled micelles are emitted as a chemical repellent, the case of $M>0$ is considered.

We define $\hat{\rho}_{0}$ as the density of the surrounding fluid without any dissolved product and $\hat{\rho}_{d}$ as the density of the droplet itself. Note that in the experiments by Krueger et al. (2016), the density difference among $\hat{\rho}_{0}$, $\hat{\rho}_{d}$ and the density of the dissolving product (filled micelles) $\hat{\rho}$ is lower than $3\%$. Therefore, we consider the density difference within the Boussinesq approximation, i.e. the density of the fluid $\hat{\rho}$ is assumed to be linearly proportional to the filled micelle concentration

[TABLE]

where $\beta$ is the proportionality constant between the density and the product concentration.

The velocity field outside the droplets is governed by the Navier-Stokes equations and the product concentration field by the advection-diffusion equation. The equations are non-dimensionalized by $R$ for lengths, $c_{0}=\alpha R/D$ for concentrations, and $\alpha M/D$ for velocities. Then the non-dimensional governing equations can be written as

[TABLE]

{subeqnarray}∂u∂t

+(u⋅∇)u=-∇p+ScPe∇^2 u-RaScPe2ce_z, ∇ ⋅u=0,

\returnthesubequationand the velocity of the droplet $\boldsymbol{U}_{d}$ satisfies

[TABLE]

where $\int(\tau\cdot n)dS$ is the force integrated over the surface of the droplet and $\boldsymbol{e}_{z}$ is the unit vector of the $z$ axis.

The dimensionless control parameters of these equations are the Rayleigh number $Ra$, which represents the strength of the product buoyancy effect,

[TABLE]

the Péclet number $Pe$, which indicates the strength of the diffusiophoretic effect,

[TABLE]

the Schmidt number $Sc$

[TABLE]

which is kept as a constant in our study, and the Galileo number $Ga$, which represents the strength of the droplet buoyancy effect,

[TABLE]

The concentration boundary condition (equation (1)) at the droplet surface reads in non-dimensional form

[TABLE]

The non-dimensional version of the velocity boundary condition at the droplet surface (equation (2)) is

[TABLE]

We apply periodic boundary conditions along the horizontal directions ($x$ and $y$) of the domain (figure 1 (b)), and solid wall boundary condition at the top and bottom of the domain. The concentration and velocity boundary conditions at the top and bottom are, respectively,

[TABLE]

and

[TABLE]

3 Numerical methods and validation

The Navier-Stokes equations and advection-diffusion equation are solved using direct numerical simulation (DNS) in Cartesian coordinates. We spatially discretize the equations using the central second-order finite difference scheme. Uniform staggered grids are used in all directions. The time integration is accomplished by using a fractional-step method. The non-linear terms are computed explicitly by a low-storage third-order Runge-Kutta scheme and the viscous and diffusion terms by a Crank-Nicolson scheme (Verzicco & Orlandi, 1996; van der Poel et al., 2015; Ostilla-Mónico et al., 2015; Spandan et al., 2018). The model for the droplet-droplet and droplet-wall collision is based on the spring-dashpot model by Costa et al. (2015).

Because the vicinity of the surface is adopted to satisfy the constant normal fluxes and slip velocity boundary condition, we cannot allow the gap between the droplets to reach zero. Therefore, when the gap width is below 2 grid spacings, we assume that the droplets are in contact.

The numerical setup is shown in figure 1. The radius $R$ of the droplet is the characteristic length of the system, and the domain size is $L_{x}\times L_{y}\times L_{z}=16\times 16\times 24$. Two droplets of unit radius are initially aligned along the $y$ axis at the center of the domain with an initial distance $L_{0}=4$. We use uniform grids $N_{x}\times N_{y}\times N_{z}=161\times 161\times 241$. Since there is the periodic boundary condition along the $y$ axis, the simulations with two droplets are actually a part of a series of droplets aligning along the $y$ axis. In $x$ direction, the periodic boundary conditions result in a balanced force. We use $\ell_{1}(t)$ and $\ell_{2}(t)$ to denote the distance between the droplet and its two neighboring droplets along the $y$ axis, and we further define $\ell(t)=\text{min}(\ell_{1}(t),\ell_{2}(t))$ as the droplet’s distance to its nearest neighbor as shown in figure 1 (b).

We will take a range of parameters based on the data in the experiment by Krueger et al. (2016): Péclet number $0.5\leq Pe\leq 10$, Rayleigh number over the range $0.1\leq Ra\leq 245$ and Galileo number $0\leq Ga\leq 0.19$. The Schmidt number in the experiments is at the order of $10^{4}$. However, the simulations at such high $Sc$ are challenging due to the very small diffusivity compared to viscosity. Therefore, in our study, we set the Schmidt number $Sc=100$.

Our code has been used to simulate diffusiophoretic particles. For the corresponding code validation, we refer the readers to our previous work (Chen et al., 2021). As a further validation, we test our code by simulating particle-laden flow with both droplet buoyancy effect and product buoyancy effect, and comparing with the existing results from the literature (Majlesara et al., 2020). These authors consider the cases about sedimenting cold/hot (fixed temperature) spherical particles in a long vertical fluid channel and study their terminal velocity. In that work, the product buoyancy effect is induced by the temperature variation, characterized by the Grashof number $Gr=Ra/Pr$. The particle buoyancy effect is characterized by the Reynolds number $Re$, which carries the same information as the Galileo number, $Re=\frac{2\sqrt{3}}{3}Ga$. We apply our code to simulate the same cases as Majlesara et al. (2020), and compare the normalized terminal velocity $U_{t}/U_{0}$ for various $Re$ and $Gr$, where $U_{0}$ is the characteristic buoyancy velocity. The numerical results are plotted in figure 2; they agree very well with those by Majlesara et al. (2020).

4 Effect of Péclet number and of Rayleigh number

In this section, we first investigate the role of the Péclet number and of the Rayleigh number by simulating the interaction between a pair of droplets with $Sc=100$, $0.5\leq Pe\leq 10$, and $0.1\leq Ra\leq 245$. The droplet buoyancy effect is absent in this subsection ($Ga=0$), and will be analyzed in subsection 5.

To demonstrate different interaction modes, we first focus on the cases $Pe=5$ and $Ra=0.1$, $2$ and $245$ in figure 3, all for $Sc=100$ as throughout in this paper. For $Ra=0.1$, the diffusiophoretic effect is dominant. The mutual repulsion drives the droplets to the horizontally balanced positions ($1/4L_{y}$ and $3/4L_{y}$). The repulsion of the neighbouring droplets acts as restoring force to the balanced position, while the droplets also experience a damping force due to the viscous drag. Therefore the droplets perform a damped oscillation in horizontal direction around the balanced positions. In vertical direction, the droplets rebound from the walls because of the concentration accumulation in between.

As $Ra$ increases to $2$, the droplets approach their neighoring droplets. The reason is that the equi-distance balanced position becomes an unstable equilibrium due to the attraction between the droplets. In the end, the two droplets reach a new balanced point of finite distance $\ell<L_{y}/2$. Along the vertical direction ($z$), due to the stronger product buoyancy effect, the droplets first sediment to the bottom. Then the droplets gradually float up as the concentration between the wall and the droplets accumulate. As $Ra$ further increases, the terminal distance between the neighbouring droplets decreases. For $Ra=245$, the product buoyancy effect becomes even stronger. The droplets are more attractive to each other and collide in the end.

From the results, we find that different strengths of the product buoyancy effect lead to different terminal distances between the droplets along the horizontal axis. Therefore, we define the terminal distances $\ell_{\infty}$ to quantify that strength:

[TABLE]

The dependence of $\ell_{\infty}$ on the Rayleigh number $Ra$ is shown in figure 4 (a). We identify two different types of interaction according to $\ell_{\infty}$: (a) $Ra<50$: without collision, where the droplets remain at an equilibrium distance without colliding with each other. The distance $\ell_{\infty}$ between the droplets reduces as $Ra$ increases. (b) $Ra\geq 50$, with collision, where the droplets collide due to the sufficiently strong attraction driven by the product buoyancy effect.

We also simulate cases for different $Pe$ and $Ra$. The results can be classified into the two mentioned interaction modes, which are presented by different symbols in figure 4 (b). The results indicate that there is competition between repulsion by diffusiophoresis and attraction by the product buoyancy effect. Higher $Pe$ results in a higher $Ra$ threshold, above which the collision occurs. This complies with the experimental results by Krueger et al. (2016), who find that the surfactant concentration ($Pe$) is increased, higher density differences ($Ra$) are needed for collective behavior to occur.

In summary, we numerically observe the interaction between droplets. We find very similar features as in the experiments by Krueger et al. (2016). While the diffusiophoretic effect (characterised by $Pe$) results in repulsion between droplets, the product buoyancy effect (characterised by $Ra$) leads to their attraction. We identify two different interaction modes: when the diffusiophoretic effect is dominant, the droplets reach a finite distance without collision; when the product buoyancy effect is dominant, the droplets collide in the end. In the next section, we will further investigate the role of the droplet buoyancy effect.

5 Effect of increasing Galileo number

In this section, we analyze the influence of the droplet buoyancy effect, as quantified by the Galileo number $Ga$. We numerically investigate the interactions between a pair of droplets with $Pe=5$, $Ra=245$ and $0\leq Ga\leq 0.19$.

Snapshots at different times are plotted in figure 5. For all examined cases with different $Ga$, the droplets collide in the end. To further analyze the interaction, we will examine the temporal change of the horizontal distance $\ell$ and the vertical height of the droplets.

We first plot the horizontal distance $\ell$ versus time $t$ in figure 6 (a). The plot indicates that as $Ga$ increases, the waiting time for the collision to occur is longer. Next, we have a close inspection of the movement of droplets near the moment of collision through plotting $\ell-\ell_{c}$ as a function of $t_{c}-t$ in log-log scale in figure 6 (b), where $\ell_{c}$ is the collision distance and $t_{c}$ is the collision time. Remarkably, all curves collapse on each other near the collision point. It suggests that the attracting behavior of the droplets are mainly determined by $Ra$ and $Pe$, while the change of $Ga$ only leads to a delayed collision.

We also plot the height $h$ (right y axis) along $t_{c}-t$ in figure 6 (c). As $Ga$ increases, the droplets wait for longer time before the occurrence of the approaching stage, and the rising velocity is also smaller for larger $Ga$. This is because a longer time is needed to build up a sufficiently large vertical concentration gradient to lift up a heavier droplet.

6 Attraction model with buoyancy

In this section, we further investigate the origin of attraction and develop a model to estimate the attracting velocity in the buoyancy-dominant cases. Employing the point heat source model, we first derive a scaling law for the horizontal velocity around a droplet, and then calculate the attractive velocity, using the methods of reflections and Faxen’s law. Since the droplet buoyancy effect only leads to delayed collision, we neglect it, i.e. we assume $Ga=0$ throughout this section.

6.1 The velocity field near a single droplet

We start by simulating a single droplet to investigate the flow around it. Figure 7 (a) shows the concentration and horizontal velocity ($v_{y}$) around a single droplet at $Ra=245$. From the fields, we observe a strong downwards plume, which leads to a higher concentration underneath the droplet. In the meantime, a horizontal flow is induced sidewards of the droplet. This inward flow drives the attraction between two nearby droplets.

We represent the buoyancy-driven flow near a single droplet by the schematics in figure 7 (b). Since the flow around a single droplet is axisymmetric, it is illustrated in cylindrical coordinates $(r,z)$. There is horizontal inward flow sidewards of the droplet, and vertical downwards flow under the droplet.

A similar case that has been well studied is the natural convection near a heat source or dissolutions source (Fujii, 1963; Moses et al., 1993; Dietrich et al., 2016). Fujii (1963) theoretically studied the buoyancy-driven convection near a fixed heat source in the fluid, and quantitatively obtained the buoyancy driven velocity. The theoretical results were later verified in experiments with a heating sphere in a fluid by Moses et al. (1993). Both the velocity and the width of the plume scale with the Rayleigh number $Ra$ (Fujii, 1963; Moses et al., 1993; Dietrich et al., 2016):

[TABLE]

[TABLE]

We define the width $h_{1},h_{2}$ of each flow branch as the distance between $10\%$ of the maximum vertical and horizontal velocity. Due to the limited domain size, the vertical and horizontal flow cannot attain the asymptotic velocity of 0, preventing us from using a smaller threshold for the $h_{1},h_{2}$ definition. The normalized values $h_{1}/h_{1}(Ra=245)$ and $h_{2}/h_{2}(Ra=245)$ versus $Ra$ obtained in simulations are plotted in figure 8. When $Ra$ is large enough ($Ra\geq 100$), the width of the channel well agrees with (16). For $Ra<100$, the numerical results deviate, because of the existence of a strong enough diffusiophoretic effect.

Given the width of both the horizontal and the vertical flow, by continuity, we can further derive the relationship between the strength of the two velocities ($v_{y}$ for horizontal and $v_{z}$ for vertical), namely

[TABLE]

where $r$ refers to the horizontal distance from the droplet center (along the red dashed curve in figure 7 (a)). Then we define the local Reynolds number $Re_{y}(r)$ using the horizontal velocity $v_{y}(r)$. With (15) and (17), we obtain:

[TABLE]

We verify (18) with the numerical simulations of a single droplet in the domain. Note that due to the periodic boundary condition, the horizontal velocity is also influenced by the neighboring droplets outside the domain,

[TABLE]

The results for different $Ra$ and $Pe$ are shown in figure 9. The numerical results agree well with the theory equation (19) for $r/R>4$. The results deviate near the droplet surface $r/R<4$, because the horizontal velocity reduces to zero approaching the droplet surface.

6.2 Droplet velocity using method of reflections and Faxen’s law

In this section, we apply Faxen’s law and the method of reflections to account for the interactions between multiple droplets. The principle of the method of reflections is to perform successive approximations for the interaction of droplets within the fluid (Guazzelli & Morris, 2011; Varma et al., 2018). The velocity of the droplet is calculated iteratively, and in each step, the velocity of the droplet is updated with the disturbance from other droplets using Faxen’s law (Guazzelli & Morris, 2011). Despite the far-field assumption of the method, even for close distance $\ell/R\sim O(1)$, it reaches a surprisingly accurate result (Ishikawa et al., 2006; Spagnolie & Lauga, 2012).

First, we consider a pair of active droplets (droplet 1 and 2) far apart. Since there is no external force and the droplet is isotropic, the droplets are stationary:

[TABLE]

where $U_{i}^{j}$ represents the velocity of droplet $i$ after the $j$th reflection process.

Then in first reflection, we suppose that the droplets are only moderately far apart, and each droplet makes a disturbance at the velocity of the other. From (18), droplet 1 causes a fluid velocity disturbance at droplet 2:

[TABLE]

where $u_{i}^{j}$ is the fluid velocity disturbance at the center of droplet $i$ caused by the other droplet after the $j$th reflection. According to the Faxen’s law, the velocity of droplet 2 due to the velocity disturbance caused by droplet 1 is (Guazzelli & Morris (2011), p.87):

[TABLE]

Since the two droplets are equivalent, the same velocity is obtained for droplet 1 after the first reflection.

Then we start with the second reflection, the velocity of the droplet obtained in the first reflection will cause disturbance to the other one. The fluid velocity caused by droplet 1 at the center of droplet 2 is (Lamb (1932), p. 599):

[TABLE]

Again with Faxen’s law, the velocity disturbance of the droplet 2 after the second reflection is given by:

[TABLE]

For higher-order reflection, it is found that the velocity disturbance after reflection

[TABLE]

Therefore, we neglect the higher-order small terms, and the velocity of the droplet is approximated as

[TABLE]

We define the Reynolds number of the droplet $Re_{d}$ by the droplet velocity $U$:

[TABLE]

$Re_{d}(\ell)$ is different from $Re_{y}(r)$, where $Re_{d}(\ell)$ expresses the velocity of a droplet influenced by the other droplet at distance $\ell$, while $Re_{y}(r)$ corresponds to the fluid velocity at distance $r$ away from a single droplet.

Equation (27) considers the influence from only one neighboring droplet. Note that the lateral boundaries are periodic. We consider the influence from the two neighboring droplets and obtain

[TABLE]

6.3 Model validation

The $Re_{d}$ of the droplets at different distances $\ell$ of different $Ra$ and $Pe$ obtained from simulations are shown in figure 10. The numerical results collapse for sufficiently large droplet separation $\ell$. For large distances $\ell/R$, the numerical results can be described by equation (28) which excellently agrees with the data for the velocity of the droplet especially for high $Ra$. For low $Ra$ ($Ra\leq 50$), there is a deviation between the numerical results and the relationship (28) near the droplet, which can be explained by the influence of diffusiophoretic flow near the droplet.

To further test our theory, we simulate the case of three droplets initially located at the center of the domain with $Ra=245$, $Sc=100$, and $Pe=5$, where the snapshots are given in figure 11(a). The horizontal velocity of the middle droplet remains at zero due to the symmetry about the middle axis. With our model of subsection 6.2, $Re_{d}$ follows:

[TABLE]

Indeed, in figure 11 (b), for large $\ell/R$ again an excellent agreement is seen between the numerical results and the prediction of equation (29).

7 Repulsive effect by diffusiophoresis

Given the good agreement between the attraction model with numerical results, now we further study the repulsive velocity from simulations.

We simulate a pair of droplets fixed at the center of the domain (figure 1 (a)) with a horizontal distance $\ell=3$ at $y$ direction. Since the repulsive diffusiophoretic motion mainly comes from the slip velocity induced by concentration field, we estimate the repulsive velocity by the integral of slip velocity at the surface (Stone & Samuel, 1996):

[TABLE]

We define $Re_{\text{rep}}=U_{\text{rep}}R/\nu$ to represent the repulsive interaction. We simulate the cases of different $Pe$ and $Ra$, and the resulting concentration field is shown in figure 12. The relationship between $Re_{\text{rep}}$ and $Pe$, $Ra$ is shown in figure 13. Figure 13 (a) shows the relationship between the normalized $Re_{\text{rep}}$ and $Pe$. The results indicate that the diffusiophoretic effect leads to repulsive motion, which agrees with our conclusions in Section 4, and $Re_{\text{rep}}$ is proportional to $Pe$. Figure 13 (b) shows $Re_{\text{rep}}$ for different $Ra$, and we find that a stronger buoyancy effect reduces the repulsive velocity between droplets. We fit the results with a power law ansatz and get

[TABLE]

To better understand the decrease of $Re_{\text{rep}}$ with increasing $Ra$, we first study the influence of $Ra$ on its surrounding concentration field. We plot the concentration gradient near a single droplet of different $Ra$ at $Pe=5$ from simulations in Figure 14 (a). It indicates that the concentration gradient has a significant drop as $Ra$ increases. This can be rationalized as follows: As $Ra$ increases, buoyancy-driven convection reduces the thickness of the concentration boundary layer (Fujii, 1963; Dietrich et al., 2016). As the surface concentration gradient remains constant (equation (10)), a reduction in the boundary layer thickness leads to a lower local concentration gradient near the droplet.

Moreover, we find that buoyancy also influences the position of the plume at the droplet surface. Through the concentration field in figure 12, as $Ra$ increases, the plume moves closer to the bottom of the droplet. To evaluate the effect, we define $\theta$ as the angle between the maximum concentration point and the droplet bottom point to represent the position of the plume as indicated in figure 12. Figure 14 (b) shows the change of the plume position for different $Ra$ at $Pe=5$. This finding thus suggests that a stronger buoyancy effect (higher $Ra$) pulls the plume towards the bottom point and this can reduce the horizontal component of the repulsive diffusiophoretic velocity.

We acknowledge that the arguments above are handwaving and qualitative. The complex system dynamics resulting from the coupling between convection and the concentration field makes a theoretical derivation of the relationship between $Re_{\text{rep}}$ and $Ra$ too challenging.

However, if we combine the equations for the attractive (27) and repulsive velocities (31), we obtain

[TABLE]

which perfectly describes the transition between the attracting and the repelling regimes, see figure 4 (b). This plot nicely reflects that the mechanism behind the interaction between droplets is the competition between the attractive force by buoyancy and the repulsive force by diffusiophoresis.

8 Summary & Conclusions

We have studied the interaction between droplets with diffusiophoretic effect, droplet buoyancy effect and product buoyancy effect. The corresponding parameters are Péclet number ($Pe$), Galileo number ($Ga$) and Rayleigh number ($Ra$). We have simulated the cases over a range of $Pe$, $Ra$, and $Ga$, with $Sc$ being fixed at $100$.

For a pair of droplets, using numerical simulations, we have found that the product buoyancy effect leads to the attractive motion between droplets, while the Marangoni/diffusiophoretic effect results in repulsion. A larger $Pe$ results in a larger $Ra$ threshold, above which droplet collision occurs. If the Rayleigh number is sufficiently small, the distance between droplets reaches an equi-distance equilibrium, and as $Ra$ increases, the closest balanced distance decreases, which indicates that the product buoyancy weakens the repulsion caused by the Marangoni/diffusiophoretic effect. For sufficiently high Rayleigh numbers ($Ra\geq 50$), the droplets collide with each other. Then we investigated the influence of droplet buoyancy effect and found that the attracting behavior is similar for different $Ga$, and the change of $Ga$ only leads to a delayed collision.

With the simulation of a single droplet, we have found that the attraction originates from convective flow induced by the density difference between the dissolving product and ambient fluid. Based on this, we have created a simple model which well describes the horizontal velocity near the droplet. The local Reynolds number is inversely proportional to the distance from the droplet as shown in equation (18).

With the above model as a starting point, we have obtained the equation for the attracting velocity of the droplet at high $Ra$ with Faxen’s law and the method of reflections. The attracting velocity is proportional to $Ra^{1/4}$ and inversely proportional to the distance between the droplets. The results have been verified by the simulation results for cases with two and three droplets.

Then we have investigated the repulsive effect by simulating the case of a pair of fixed droplets and the repulsive velocity was approximated by the integral of the slip velocity (30). We have found that $Re_{\text{rep}}$, which represents the repulsive velocity, is proportional to $PeRa^{-0.38}$. The linear dependence of $Re_{\text{rep}}$ on $Pe$ is simply due to a larger diffusiophoretic repulsive force for larger $Pe$. In contrast, the $Ra$-dependence of $Re_{\text{rep}}$ is more complicated. It reflects that an increasing $Ra$ leads to a smaller horizontal concentration gradient and favours the plume to be closer to the bottom point of droplet, which reduces the repulsive velocity.

Combining the scaling relations of the attractive and repulsive velocity, we obtain $Pe\sim Ra^{0.63}$, which perfectly describes the transition curve between the attractive and repulsive regime in figure 4 (b). This indicates that the mechanism behind the interaction between droplets are the competition between attractive buoyancy force and repulsive diffusiophoretic force.

The present work contributes to the understanding of the interaction between active droplets, and specifically reveals the significant role played by the dissolving product buoyancy. It shows that product buoyancy can lead to attractive motion between active particles, which helps us understand the attraction of active droplet in the experiments of Krueger et al. (2016). We have proposed a simple model to predict the velocity of the interacting active droplets. The results provide a framework to understand the droplet attraction induced by the convective flow. Moreover, the present work reveals a possible way to change the collective behaviors by tuning the buoyancy.

In our simulation, the propulsion of the active droplet is simply modelled as diffusiophoreis. Alternatively, we could have taken Marangoni flow. Until now, buoyancy-driven attractive motions are only observed in the cases of active droplets but scarcely in phoretic particles, possibly due to the difficulties to generate large enough density difference between the product and ambient fluid by phoretic particles.

Many questions remain open. For example, how to determine the cluster size for multiple droplets? How does the flow field change if droplets are near a fluid-air interface? How to quantitatively determine the threshold Rayleigh number above which the droplets collide with each other and show collective behaviours? With the obtained insights into the attraction here, we hope it is seen as worthwhile to further explore the formation and motion of a cluster of active particles.

Acknowledgements

We greatly appreciate valuable discussions with Martin Assen, Utkarsh Jain, and Corinna Maass. We acknowledge the support from the Netherlands Center for Multiscale Catalytic Energy Conversion (MCEC), an NWO Gravitation program funded by the Ministry of Education and support from the ERC-Advanced Grant ”DDD” under the project number 740479. The simulations in this work were carried out on the national e-infrastructure of SURFsara, a subsidiary of SURF cooperation, the collaborative ICT organization for Dutch education and research, MareNostrum 4 which is based in Spain at the Barcelona Computing Center (BSC) under PRACE projects 2018194742, 2020225335 and 2020235589, on Irene at Trés Grand Centre de calcul du CEA (TGCC) under PRACE project 2019215098 and on Marconi successor at CINECA, Italy under PRACE project 2019204979. K. L. Chong is supported by the Natural Science Foundation of China under Grant No. 92052201.

Declaration of interests

The authors report no conflict of interest.