2.1. Effects of boundary condition simplification

As mentioned in the previous section, most of the simulations in this paper consider the channel intersection as the computational geometry with Dirichlet boundary conditions at its inlets and outlets given by the Boussinesq flow in a straight, rectangular channel. Although this configuration is more physically accurate than the external flow used in previous theoretical works for the Stokes trap (e.g. Lin et al. Reference Lin, Kumar, Richter, Wang, Schroeder and Narsimhan2021), these boundary conditions should be imposed on a section of a channel branch away from the intersection to more accurately model typical experimental set-ups, either by simulating the droplet in the full microfluidic domain or by using a moving-frame construction like the one described in Roure et al. (Reference Roure, Zinchenko and Davis2023) (see figure 1c), in which the frame changes with time as the drop moves and/or deforms. When doing the latter, it is still possible to implement the methods proposed in the present paper, including the control algorithm, by using a transient mode superposition of pre-calculated potential densities. Here, we briefly assess the main quantitative differences that may appear when considering the alternative inlet/outlet boundary conditions on sections of the branches away from the intersection.

We start by looking at the effect of branch length (i.e. the distance from the intersection to where we apply the boundary conditions within the rectangular inlet/outlet branches) on the undisturbed flow (i.e. without the droplet). Figure 2(a) shows a comparison between the undisturbed flow field for the hexagonal prism domain shown in figure 1(b) (in black) and the flow field of the full geometry shown in figure 1(a) (in blue), for a channel with dimensionless branch lengths $L = 2.0$. As seen by the results in figure 2, although the flow fields are qualitatively similar, they present notable quantitative differences that arise between the direct application of the Boussinesq boundary conditions at the hexagonal computational cell versus at a dimensionless distance $L = 2.0$ from the hexagonal intersection. As expected, and seen in figure 2(b), increasing the length of the channel branches results in much smaller discrepancies between the velocity fields. For $L = 1.5$, these differences become very small, with discrepancies of less than $0.5\,\%$ for the bulk of the channel (except for the origin, where the velocity is zero), indicating that the limit of very long inlet/outlet channels is well reached for $L = 2.0$. As the results in this paper are more focused on the effects of flow-mode superposition on the drop deformation and mixing dynamics than the modelling of a more physically accurate channel flow, we use primarily the simplified hexagonal construction shown in figure 1(b). However, as quantitative differences in drop shape may arise from considering the full channel geometry, in § 3 we address some of these differences between the simplified external flow field and simulations of a droplet inside a full-channel representation of a Stokes trap with branch size $L = 2.0$ performed using the moving-frame approach from Roure et al. (Reference Roure, Zinchenko and Davis2023). As in our previous work, the computational frame expands as the drop stretches, and it moves if conditions are such that the drop drifts from its initial location. Overall, the results for the full-channel simulations are qualitatively similar to those for the simplified geometry, but with less-deformable droplets. As we show in the next section, this discrepancy is characterized by a quantitative difference in the components of the spherical-harmonic decomposition of the drop shape. These discrepancies are usually small for smaller droplets, but can become up to ${\sim }40\,\%$ in some extreme cases for larger, more-deformable droplets.

Figure 2. Effect of channel branch length on background flow in a Stokes trap. (a) Qualitative comparison between the vector fields for $L = 0$ (black) and $L = 2.0$ (blue) at the same points. (b) Colour maps show the scalar discrepancy between the vector fields for different channel branch lengths.

3.1. Characterizing drop deformation via spherical harmonics

Often in the literature, drop deformation is characterized by parameters such as the Taylor deformation, which measures the deviation of a drop from its spherical equilibrium shape. From figure 3, it is clear that the different drop shapes induced by the distinct flow modes show substantial variation. Hence to properly describe drop deformation, we need a more precise way to characterize the deformed drop shape. One way to do this, valid for star-shaped drop geometries, is to decompose the shape of the droplet in spherical harmonics. Namely, the shape of the droplet is described by the function $r(\theta,\varphi )$, where $r$ is the spherical radial coordinate from the drop centre, and $\varphi$ and $\theta$ are the azimuthal and polar angles, respectively. This function can be decomposed as

(3.9)\begin{equation} r(\theta,\varphi) = \sum_{\ell,m} c_{\ell m}\,Y_{\ell}^{m}(\theta,\varphi), \end{equation}

where $Y_{\ell }^m(\theta,\varphi )$ are the normalized spherical harmonics, and the coefficients $c_{\ell m}$ are

(3.10)\begin{equation} c_{\ell m} = \int_{S^2} r(\theta,\varphi)\,\bar{Y}_{\ell}^{m}(\theta,\varphi) \,{\rm d}\varOmega, \end{equation}

with $S^2$ the unit sphere, and overbar denoting complex conjugation. Numerical implementation of (3.10) consists of first projecting the unstructured drop mesh on the unit sphere and using a linear quadrature (e.g. Zinchenko et al. Reference Zinchenko, Rother and Davis1997).

As an example of harmonic decomposition, we examine the case of a droplet undergoing a tri-axial extensional flow shown previously in figure 3. Figure 5 shows numerical results for (a) the harmonic decomposition of drop shape in a tri-axial extensional flow mode, and (b) reconstruction of drop shape by using the first few harmonics. As the tri-extensional flow is locally quadratic, we expect the main excited harmonic to be $Y_{33}$, as in the regime of small deformations. Indeed, as shown in figure 5(a), besides $Y_{00}$, the largest harmonic in the spectrum is $Y_{33}$, followed by $Y_{66}$. All other harmonics are substantially smaller, which is supported by the fact that we can reconstruct the drop shape with good accuracy by using only three harmonics, as shown in figure 5(b).

Figure 5. Harmonic decomposition of the shape of a droplet in a Stokes trap undergoing a tri-axial extensional flow for $Ca = 0.1$, $\lambda = 1$, $H = 1$ and $a = 0.5$. The results show (a) the evolution of the $Y_{33}$ and $Y_{66}$ harmonics with time as the drop extends, and (b) the reconstruction of drop shape from the three main harmonics for $t = 0.25$. The dashed curves in (a) are the same harmonics for a full-channel simulation with branch length $L = 2.0$, whereas the solid curves are for the simplified hexagonal geometry. The dashed shapes in (b) are the numerical drop shape, whereas the solid lines are the harmonic approximations. The meshed geometry in (b) is a three-dimensional visualization of the harmonic reconstruction using the main three modes. (c) A comparison between the simulations in the simplified hexagonal channel (solid contours) and full channel (dashed contours) from (a).

The results shown in figure 5 indicate that a good parameter to measure drop deformation in a tri-axial extensional flow is the imaginary part of $c_{33}$. The results for the full-channel simulations display a similar increasing behaviour, although with smaller drop deformation, which is characterized by lower absolute values of the harmonic coefficients $c_{33}$ and $c_{66}$. This discrepancy in deformation also leads to a difference in long-time behaviour, as the droplet in the full-channel simulation will eventually reach a stationary state, whereas the droplet in the simplified geometry escapes the hexagonal region in three different directions, suggesting that it will eventually break up. Hence we can perform numerical experiments with our trapped droplet. One possible example of such a numerical experiment is an oscillatory flow of the form

(3.11)\begin{equation} \boldsymbol{Q}_0(t) = \boldsymbol{Q}_{{tri}}\cos(\omega t), \end{equation}

where $\omega$ is the inlet frequency. Figure 6(a,b) show the response of $c_{33}$, normalized by the drop radius $a$, to the oscillatory tri-extensional flow for $Ca = 0.1$, $\lambda = 1$, $H = 1$ and $\omega = 3$ for different values of drop radius $a$. The vertical dashed lines in figure 6 indicate the points in time where $\boldsymbol {Q} = \boldsymbol {0}$. As expected, larger droplets are more deformable, which results in a larger values of $\text {Im}(c_{33})/a$. Moreover, for droplets with radii $a=0.4$ or less (figure 6a), besides an initial transient regime, the harmonic response displays a sinusoidal behaviour slightly out of phase with the flux, indicating a small lag in the drop response, which is expected given the elastic character of the interfacial tension. We also note that drop size slightly affects the phase of drop oscillation. For the full channel simulations, shown as the dashed curves, we again observe a smaller deformation, but with the same qualitative trends. For $a = 0.5$, for the full-channel simulations (dashed curves), we observe a similar trend, with a harmonic deformation. However, for the hexagonal channel simulation, where deformations are larger, we observe non-sinusoidal oscillations, as shown in figure 6(b), indicating a ‘nonlinear’ response. Note that the oscillations in figure 6(b), besides displaying non-harmonic behaviour, also have slightly increasing peaks, indicating that either a larger time is required for the droplet oscillation to reach periodic behaviour or these oscillations are unstable (in the sense that they may potentially lead to drop breakup).

Figure 6. Numerical results for the imaginary part of the harmonic amplitude $c_{33}$, normalized by the drop radius, versus time for a droplet undergoing an oscillatory tri-axial extensional flow $\boldsymbol {Q}_0 = \boldsymbol {Q}_{{tri}} \cos (\omega t)$. The results consider $Ca = 0.1$, $\lambda = 1$, $H = 1$, $\omega = 3$, and (a) $a = 0.25, 0.3, 0.4$, and (b) $a = 0.5$. (c) The harmonic response for $a = 0.4$, and $\omega = 0.0, 0.5, 1.0, 1.5, 2.0, 2.5$. (d) The frequency response for the drop sizes in (a). The solid curves are the results for the simplified geometry, whereas the dashed curves in (a,b,d) were obtained by a full-channel simulation with branch length $L = 2.0$. The dotted curves in (a,b) indicate the amplitude of the flux $Q_1$, whereas the vertical lines indicate the times where $\boldsymbol {Q} = \boldsymbol {0}$.

Besides drop size and physical parameters such as $Ca$ and $\lambda$, an important quantity in oscillatory flows that influences the extension of drop deformation is the oscillation frequency $\omega$. The results in figures 6(c) and 6(d) show, respectively, the influence of different frequencies on the harmonic $c_{33}$, and a frequency ramp for oscillation amplitudes for different drop sizes. As expected from classical results for droplets under simple oscillatory extensional flows, smaller oscillation frequencies result in larger drop deformations, which are characterized by larger values of $\text {Im}(c_{33})$.

We can also use our method to simulate the stretching and relaxation behaviour of a droplet under step-strain conditions, similar to the vesicle experiments performed in Kumar et al. (Reference Kumar, Richter and Schroeder2020b). Besides, the decomposition in spherical harmonics allows us to quantify drop deformation for both regular extensional ($\boldsymbol {Q}_{{ext}}$) and tri-extensional ($\boldsymbol {Q}_{{tri}}$) flow modes. To illustrate this type of numerical experiment, figure 7 shows numerical results for a droplet undergoing a step-strain, tri-axial extension for $Ca = 0.1$, $\lambda = 1$, $H = 1$ and $a = 0.5$. The flow configuration is given by $\boldsymbol {Q}_{{tri}}$ for $t \le 0.2$, and $\boldsymbol {0}$ for $t >0.2$. For $t \le 0.2$, the drop experiences a tri-axial extension that is characterized by an increase in the $Y_{33}$ harmonic, as shown in figure 5. As the external flow suddenly stops, the droplet shape relaxes towards its initial spherical configuration. The results for the full-channel simulation, represented by the dashed curve, display a similar qualitative behaviour, with amplitudes approximately 30 % lower.

Figure 7. Numerical results for the $Y_{33}$ harmonic response of a droplet undergoing a step tri-axial strain with $Ca = 0.1$, $\lambda = 1$, $a = 0.5$, $H = 1$, and $\boldsymbol {Q}_0 = \boldsymbol {Q}_{{tri}}$ for $t \le 0.2$, and $\boldsymbol {Q} = \boldsymbol {0}$ for $t > 0.2$. The result represented by the solid curve is for the simplified hexagonal domain, whereas the dashed curve is the result for the full-channel geometry with branch length $L = 2.0$.

3.2. Mode combination and shape manipulation

As shown in the previous subsection, subjecting the droplet to different flow modes in the Stokes trap may excite specific harmonics on the drop interface – a feature that can be used for shape manipulation. These modes can also be combined together to produce more complex drop shapes. To illustrate this observation, figure 8 shows numerical simulations for droplets undergoing a combination of (a) shear + tri-axial extension and (b) extension + tri-axial extension for $Ca = 0.05$, $\lambda = 1$ and different drop radii $a$.

Figure 8. Numerical results for drop shapes resulting from combinations (a) tri-extensional + shear and (b) tri-extensional + extensional flow modes, for $Ca = 0.05$, $\lambda = 1$, $H = 1$, and different drop radii. The solid contours represent steady shapes, whereas the long-dashed contours, corresponding to (a) $t = 0.925$ and (b) $t = 0.35$, represent larger drops that eventually escape the intersection, possibly leading to breakup. The short-dashed contours (overlapping the solid contours) were obtained by the harmonic superposition of the half modes from table 1, and essentially coincide with the full simulations.

The results shown in figure 8 show that the combination of different flow modes can result in asymmetric drop shapes such as the ones shown in figure 8(a), while the droplet is kept at the centre of the channel by the controller. This asymmetry is even more pronounced for larger droplets. However, above a certain radius threshold, the droplet becomes too elongated, escaping the intersection, possibly leading to breakup. For small radii $a$, for which the droplet undergoes small deformations, we observe a linear deformation regime, where the harmonics superpose linearly. This linear superposition for small droplets is shown in table 1, which gives numerical results for the main spherical harmonics present in the steady shapes shown in figure 8, with additional results from numerical simulation using the isolated ‘half modes’ $\boldsymbol {Q}_{{tri}}/2$, $\boldsymbol {Q}_{{ext}}/2$ and $\boldsymbol {Q}_{{sh}}/2$. The numerical results show that for small droplets (e.g. $a = 0.3$), the harmonic spectrum of the final shape due to a combined mode seems to be given by a linear superposition of the isolated modes, as expected from the theory of small deformations (Leal Reference Leal2007). In contrast, as in the case for the oscillatory flow, we start to observe a nonlinear behaviour for large droplets, which is characterized by the lack of linear harmonic superposition and the presence of extra harmonics. As an example, we can start to see some small discrepancies for $a= 0.4$ between the final harmonics and the ones from the isolated modes. This linear superposition of harmonic modes for moderate deformations is further supported by the reconstruction of the final shape by a linear combination of the half modes, which is shown by the short-dashed contours in figure 8. For $a = 0.5$, as in the harmonic response to the oscillatory flow presented in § 3.1, linearity is no longer observed.

Table 1. Numerical results for the steady-state harmonic decomposition of different simple and combined flow modes for $Ca = 0.05$, $\lambda = 1$ and $H = 1$. The results are rounded to three decimal places.

3.3. A hydrodynamic ‘three-phase rotor’

Drop rotation often plays an important role in improving the chaotic mixing inside the droplets (Stone, Nadim & Strogatz Reference Stone, Nadim and Strogatz1991), which is crucial for applications in drop-based microreactors. It is known from the literature that droplets in a simple shear flow display a tumbling motion for high viscosity ratios (Bilby & Kolbuszewski Reference Bilby and Kolbuszewski1977; Wetzel & Tucker Reference Wetzel and Tucker2001; Oliveira & Cunha Reference Oliveira and Cunha2015). In contrast, droplets with low viscosity ratios reach a steady-state deformation, where they present a ‘tank-treading’ motion (Kennedy, Pozrikidis & Skalak Reference Kennedy, Pozrikidis and Skalak1994). Whereas a four-roll mill can produce rotational flow, it is difficult, and possibly not feasible, to produce a flow that is locally rotational at the centre of a Stokes trap. For example, the model used by Shenoy et al. (Reference Shenoy, Kumar, Hilgenfeldt and Schroeder2019) represents the external flow inside a Stokes trap as a superposition of potential sources, which is irrotational. In fact, the mode shown in figure 3(b), which we labelled as ‘shear’, is locally an extensional flow.

However, the six-branch Stokes trap allows us to generate a rotating extensional external flow by combining the three different possible modes for the ‘shear’ mode, $\boldsymbol {Q}_{{sh}}$, $\mathcal {S}\boldsymbol {Q}_{{sh}}$ and $\mathcal {S}^2 \boldsymbol {Q}_{{sh}}$, off phase by $2 {\rm \pi}/3$, in the form

(3.12)\begin{equation} \boldsymbol{Q}_{{rotor}} = \tfrac{1}{3} \left[ \boldsymbol{Q}_{{sh}} \cos(\omega t) + \mathcal{S} \boldsymbol{Q}_{{sh}} \cos(\omega t + 2 {\rm \pi}/3) + \mathcal{S}^2 \boldsymbol{Q}_{{sh}} \cos(\omega t + 4 {\rm \pi}/3) \right]. \end{equation}

Note, for example, that for $t = {\rm \pi}/2$, $\boldsymbol {Q}_{{rotor}} \propto \boldsymbol {Q}_{{ext}}$. Numerical results for a single droplet in a Stokes trap undergoing the external flow produced by the flux configuration described in (3.12) with $\omega = 3$ are shown in figure 9. The droplet starts with a spherical shape and undergoes a transient extension regime until it reaches a periodic regime, with frequency $\omega$, similar to the wobbling motion observed in high-viscosity-ratio droplets undergoing a simple shear flow. An animated version of the drop motion can be found in the supplementary material.

Figure 9. Numerical results for a droplet undergoing a three-phase extensional flow for $Ca = 0.05$, $\lambda = 1$, $a = 0.4$, $H = 1$ and $\omega = 3$. The motion of the drop is comprised of a short transient regime, shown in (a), where the droplet transitions from a spherical to an ovoid shape, and a periodic wobbling regime, shown in (b). The timeline at the bottom displays the full motion of the droplet. The solid drop shape for each part represents the first drop configuration for that part (i.e. first and fourth panels), whereas the dashed shape corresponds to the last configuration displayed on the timeline for that part (i.e. third and seventh panels). (c) The effect of the frequency $\omega$ on the maximum amplitude of the $c_{22}$ harmonic.

From the frequency response curve shown in figure 9(c), we see that, as in the case of the oscillatory tri-axial flow in § 3.1, the increase in rotation frequency leads to a smaller drop deformation. One interesting feature of this type of flow is that at each time step, the external flow acts similarly to an extensional flow. As the internal flow inside a droplet undergoing an external extensional flow has four circulation regions, a continuously rotating extensional flow will continuously change these invariant mixing regions, making it possible to observe effective mixing inside the droplet. In the next section, we investigate how this rotating flow mode may be used to produce active mixing inside the droplet.

4.1. General remarks

Besides influencing drop shape and deformation, the different flow modes investigated in the prior portion of this paper also affect the internal circulation. In this section, we investigate how these induced internal flows influence mixing inside the droplets.

Due to the large number of microfluidic applications, chaotic mixing inside droplets has been investigated extensively in the literature, both theoretically and experimentally (Sarrazin et al. Reference Sarrazin, Loubiere, Prat, Gourdon, Bonometti and Magnaudet2006; Blanchette Reference Blanchette2010; Zhao et al. Reference Zhao, Riaud, Luo, Jin and Cheng2015; Chen et al. Reference Chen, Zhao, Wang, Zhu, Tian, Xu, Wang and Huang2018). For example, chaotic mixing inside isolated spherical droplets induced by linear external flows was investigated by Stone et al. (Reference Stone, Nadim and Strogatz1991), inducing chaos by applying external, non-aligned extensional and rotation flows. For quadratic flows, even earlier results by Bajer & Moffatt (Reference Bajer and Moffatt1990) show a stretch-twist-fold chaotic dynamics. Later, Stone & Stone (Reference Stone and Stone2005) investigated a transient combination of shear and uniform flows to emulate the conditions in a typical microfluidic serpentine mixer. This work was also extended to direct numerical simulations of deformable two-dimensional droplets in a serpentine channel (Muradoglu & Stone Reference Muradoglu and Stone2005) using a finite-volume/front-tracking scheme. The characterization of mixing inside droplets in different microfluidic channels is still being explored (Fu et al. Reference Fu, Wang, Zhang, Bai, Jin and Cheng2019; Belousov et al. Reference Belousov, Filatov, Kukhtevich, Kantsler, Evstrapov and Bukatin2021; Cao et al. Reference Cao, Zhou, Yu and Liu2021). Recently, Gissinger, Zinchenko & Davis (Reference Gissinger, Zinchenko and Davis2021) used boundary-integral simulations to investigate mixing inside a droplet trapped in constrictions formed by rigid particles and fibres. Beyond the context of microchannels, the work by Watanabe, Hasegawa & Abe (Reference Watanabe, Hasegawa and Abe2018) also explored active mixing inside droplets in an acoustic trap.

Here, similarly to the previous work by Gissinger et al. (Reference Gissinger, Zinchenko and Davis2021), we use boundary-integral simulations to study the mixing dynamics inside the droplet in a Stokes trap. One of the advantages of boundary-integral methods for mixing simulations is that the velocity of the fluid at any point inside the droplet can be determined by the numerical evaluation of (2.4), without requiring interpolation. In contrast to the problem considered in Gissinger et al. (Reference Gissinger, Zinchenko and Davis2021), our droplet undergoes continuous deformation, not being confined to a steady state. To illustrate the velocity calculation inside the droplets for deformable droplets, figure 10(a) shows the short-time evolution of streamlines inside an initially spherical droplet. It reveals the formation of six circulation regions at the midplane $z = 0$ inside a droplet undergoing a tri-axial extensional flow. At $t = 0$, the flow inside the droplet resembles the undisturbed external flow shown in figure 3(a). As the droplet deforms, we start to see six saddle-like fixed points moving from the drop centre towards its boundary, giving rise to the six circulation regions. For droplets with small deformations (e.g. $Ca \to 0$), this formation happens almost instantly. These circulation regions are very similar to those observed inside a spherical droplet subjected to an external quadratic flow.

Figure 10. Flow inside an initially spherical droplet subject to an external tri-axial extensional flow with ${Ca = 0.1}$, $\lambda = 1$, $a = 0.4$ and $H = 1$. (a) The transient formation of six circulation regions inside the droplet. (b) The details of the mixing simulations, including the regions $V_{{dye}}$ (in black) and $V_{{clear}}$ (in white) used in the calculation of the mixing number. The final configuration is calculated by backtracing the centres of cells in a Cartesian grid to their initial positions.

Although there are multiple ways to characterize mixing in fluids, one that is particularly interesting in our case is the mixing number, introduced by Stone & Stone (Reference Stone and Stone2005). This particular choice comes from the visual intuitiveness of such a parameter, and its connection to the convective transport of a fluid region. This parameter was also shown to be correlated (or inversely correlated) to other quantities such as entropy, intensity of segregation, and other mixing quantifiers (Muradoglu & Stone Reference Muradoglu and Stone2005; Hoeijmakers et al. Reference Hoeijmakers, Araujo, van Heijst, Nijmeijer and Trieling2010; Gissinger et al. Reference Gissinger, Zinchenko and Davis2021; Roshchin & Patlazhan Reference Roshchin and Patlazhan2023). To illustrate the definition of this metric, we focus on the classical example of an initially spherical droplet with passive dye completely filling one of its hemispheres. This problem is illustrated in figure 10(b). As time passes, the dye is advected with the same velocity as the flow inside the droplet. As the dye is assumed to be non-diffusing, the sets $V_{{dye}}$, consisting of the dyed points, and $V_{{clear}}$, consisting of the clear points, are disjoint. The mixing number is a measure of closeness between the two disjoint sets. In our specific case, we define it in the following grid-independent form:

(4.1)\begin{equation} m(t) = \frac{1}{a^2 V_d} \left[\int_{V_{{dye}}} d^2(\boldsymbol{x},V_{{clear}}) \,{\rm d} V + \int_{V_{{clear}}} d^2(\boldsymbol{x},V_{{dye}}) \,{\rm d} V \right], \end{equation}

where $d^2(\boldsymbol {x},A) = \inf _{\boldsymbol {y}\in A} d^2(\boldsymbol {x},\boldsymbol {y})$ is the square of the distance between the point $\boldsymbol {x}$ and the set $A$. The normalization factor $a^2$ is used to make the mixing number non-dimensional and to avoid an extra drop-size dependency in $m(t)$. If the system is well mixed, then the two sets become strongly intertwined and the mixing number approaches zero, hence providing an inverse measure of mixing between the two sets.

In practical applications, as the nonlinear dynamics of the system must be solved numerically, the implementation of (4.1) is performed in a coarser domain given by a Cartesian grid formed by cubic cells of the same volume. Under these circumstances, the mixing number is given by

(4.2)\begin{equation} m(t) \approx \frac{1}{a^2 N_{g}} \sum_{k=1}^{N_{g}} d^2(\boldsymbol{x}_k,\text{Opp}(\boldsymbol{x}_k)), \end{equation}

which is the same expression as in Stone & Stone (Reference Stone and Stone2005). Here, the summation is over the grid cells, $N_g$ is the total number of cells, $\boldsymbol {x}_k$ is the midpoint of a cell $k$, and

(4.3) \begin{equation} \text{Opp}(\boldsymbol{x}) = \left\{\begin{array}{@{}ll@{}} V_{{dye}}, & \text{if} \ \boldsymbol{x} \in V_{{clear}},\\ V_{{clear}}, & \text{if} \ \boldsymbol{x} \in V_{{dye}}, \end{array}\right. \end{equation}

where $V_{{dye}}$ and $V_{{clear}}$ are the coarse-grained versions of their continuous counterparts. For a two-dimensional system, the definition of the mixing number (4.1) would be similar but with areas instead of volumes. Its numerical counterpart (4.2), however, would remain unaltered, with the exception that the Cartesian grid would now be two-dimensional.

The regions $V_{{dye}}$ and $V_{{clear}}$ are determined by tracing the centre point of each cell to its starting position and using the initial condition for dye concentration (see figure 10b). This method is referred to as the backward Poincaré cell method (Wang et al. Reference Wang, Fan and Chen2001) and has been used in previous works to obtain graphical representations for the chaotic mixing inside droplets. For incompressible flows, such a method yields more accurate results for the mixing number when compared to forward propagation, as it consists in a direct discretization of the final concentration profile obtained by the method of characteristics for the advective transport equation (Roure & Davis Reference Roure and Davis2021). Namely, if the dye concentration $c(\boldsymbol {x},t)$ undergoes a purely advective transport, with $\partial c/\partial t + \boldsymbol {u} \boldsymbol {\cdot } \boldsymbol {\nabla } c = 0$ and initial condition $c_0(\boldsymbol {x})$, then the concentration at a time $t$ is given by $c_0(\varPsi _t^{-1} (\boldsymbol {x}))$, where $\varPsi _t$ is the time evolution of the dynamical system from a starting position at $t = 0$. Note that forward propagation is still necessary to calculate quantities such as the mixing entropy (Muradoglu & Stone Reference Muradoglu and Stone2005), which we do not explore in this work.

For calculation of both regular and backward trajectories, we use a second-order Runge–Kutta scheme. The fluid velocity field inside the droplet at each time step is calculated by the numerical evaluation of (2.4). To this end, the drop shapes and potential densities for the relevant time steps are pre-calculated and stored by solving the boundary-integral problem. To keep track of the points inside the droplet, we use an indicator function

(4.4)\begin{equation} I(\kern0.7pt\boldsymbol{y}) = \frac{1}{4 {\rm \pi}} \int_{S_d} \frac{\boldsymbol{n}(\boldsymbol{x})\boldsymbol{\cdot}(\boldsymbol{x} - \boldsymbol{y})}{\| \boldsymbol{x} - \boldsymbol{y} \|^3} \,{\rm d}S_{\boldsymbol{x}}, \end{equation}

which is 1 for $\boldsymbol {y}\in V_d$, and 0 for $\boldsymbol {y} \notin V_d$. For drop-surface discretization into flat mesh triangles $\triangle$, each triangle contribution to (4.4) is simply the signed area of the spherical triangle obtained by projecting $\triangle$ onto the unit sphere centred at $\boldsymbol {y}$, and this area is found analytically from spherical trigonometry (as in Zinchenko & Davis Reference Zinchenko and Davis2013).

4.2. Mixing in deformable droplets

We now analyse how different flow modes may influence the mixing inside deformable droplets. As we are considering the droplet to be neutrally buoyant and centred at $z = 0$, the plane $z = 0$ is a two-dimensional invariant manifold for all the flow modes considered in this paper. As the mixing in this two-dimensional submanifold often correlates with overall three-dimensional mixing inside the droplet (Stone & Stone Reference Stone and Stone2005), we focus our mixing analysis on the cross-section $z = 0$.

One of the main consequences of drop deformation is the breaking of the kinematic reversibility usually present in Stokes flows, which may improve the mixing inside the droplet for specific cases. To illustrate the effects of irreversibility, we return to the oscillatory flow problem discussed in § 3.1. For a perfectly spherical droplet (e.g. $Ca = 0$), the linearity of Stokes equations implies that an oscillatory flow mode would not produce any type of effective mixing inside the droplet. Instead, all points would return to their initial position after one period. In contrast, for a deformable droplet, this kinematic symmetry is broken by the drop deformation, meaning that a given material point inside the droplet will display non-periodic dynamics, as illustrated in figure 11(a). This symmetry breaking is similar to that observed in the relative trajectories of particles leading to hydrodynamic dispersion (Cunha & Hinch Reference Cunha and Hinch1996; Loewenberg & Hinch Reference Loewenberg and Hinch1997; Roure & Cunha Reference Roure and Cunha2018).

Figure 11. Symmetry breaking of kinematic reversibility caused by drop deformation. (a) A droplet with $a = 0.4$, $\lambda = 1$ and $Ca = 0.1$ undergoing a periodic deformation caused by an external oscillatory tri-axial extension flow. After one period, the material point presents a displacement from its initial position. (b) A Poincaré section at $z=0$ for three initial positions (A, B, C).

To better visualize the global effects of this symmetry breaking for the oscillatory tri-axial extension flow mode, figure 11(b) shows a Poincaré section for three different starting points inside the droplet at the $z=0$ plane and $t = 1.625$, where the droplet shape is approximately spherical and the drop displays a periodic motion. Each point in the discrete trajectories shown in figure 11 corresponds to the material particle position after one period of oscillation. The results in figure 11(b) show the existence of non-periodic orbits, which are caused by drop deformation. From figure 11, we see that near the centre of the droplet (i.e. away from the surface), the symmetry breaking is small, as indicated by the points very close to each other. In contrast, near the surface, we observe a more noticeable deviation from periodicity, as indicated by the presence of two attractor-like structures. Due to the nature of our numerical method, it is hard to tell precisely if the Poincaré map is spiralling down to an attractor or if the structure consists of quasi-periodic orbits. In both cases, the breaking of periodicity is clear.

Although the breaking of the kinematic reversibility due to drop deformation may potentially improve mixing locally, it alone does not guarantee a full mixing inside the droplet, especially in the plane $z = 0$. For example, for the tri-axial extension mode, as in the problem of a spherical droplet under a quadratic flow, the internal dynamics is constrained to the six symmetry quadrants inside the droplet. One way to overcome this issue and to induce a more effective mixing even in the midplane $z = 0$ is to use a time-dependent combination of flow modes, which is, in fact, the main strategy used in traditional microfluidic mixers. One such alternative would be the previously discussed three-phase mode $\boldsymbol {Q}_{{rotor}}$ discussed in § 3.3. Figure 12 shows numerical simulations of mixing inside a droplet undergoing a three-phase extensional external flow mode for $a = 0.4$, $Ca = 0.1$ and $H = 1$ for different times and different values of viscosity ratio and frequencies. The number below each droplet is the mixing number $m(t)$, calculated using (4.2).

Figure 12. Numerical simulations of mixing inside a droplet undergoing a three-phase extensional flow for $a = 0.4$, $Ca = 0.1$, $H = 1$ at different times for distinct values of viscosity ratio and frequency $\omega$. The results are for the midplane $z = 0$. The number below each droplet is the mixing number $m(t)$, calculated using (4.2). Droplets with a lower viscosity ratio present a better mixing, which is indicated by a smaller mixing number.

As in the example shown in figure 10(b), the droplet starts with black points in the lower region and white points in the upper region. The final configuration of the points is calculated by using the backward Poincaré cell method described in this section. One immediate result, expected from the results found in Stone & Stone (Reference Stone and Stone2005), is that mixing is more effective for less-viscous droplets. This result is indicated by the very small mixing numbers and happens because the lower viscosity of the droplet results in a faster internal advection.

Another important factor in the mixing induced by the three-phase extensional flow is the frequency of the flow. Namely, for high frequencies such as $\omega = 6$, we observe very little overall mixing. However, for low frequencies (e.g. $\omega = 1.5$), we observe a more effective mixing even for $\lambda = 1$. The reason behind the better mixing effectiveness is similar to the increase in effectiveness caused by lowering the viscosity ratio, namely the interplay between internal circulation and drop rotation. For high values of $\omega$, the droplet rotates much faster than the time it takes for the internal flow to advect the passive dye. For lower values of $\omega$, internal advection happens faster than rotation, allowing for a more effective mixing in a shorter time. Of course, for $\omega = 0$, the inner flow becomes steady, meaning that an effective mixing in the $z=0$ plane is impossible. Hence there should be an optimal value of $\omega$ to promote mixing.

As mentioned in the beginning of this subsection, drop deformation often plays an important role in mixing. Although our results from figure 11 indicate that drop deformation can potentially aid mixing inside the droplet by breaking the kinematic reversibility of Stokes flow, earlier results by Muradoglu & Stone (Reference Muradoglu and Stone2005) show an opposite trend. In fact, in our system, we also observe situations where drop deformation slows down mixing. As an example, figure 13 shows the results for numerical mixing simulations of a droplet subject to a three-phase extensional flow with $\omega = 3$, $\lambda = 1$, $H = 1$, $a = 0.4$ and $Ca = 0.05$. Comparing the mixing numbers with the result shown in the second row of figure 12, we see that, like the results in Muradoglu & Stone (Reference Muradoglu and Stone2005), a smaller $Ca$ results in a better mixing – although the difference between the two cases is less pronounced in our system. This result is characterized by the lower mixing numbers at most time steps. Similarly to the effect of frequency on mixing, this improvement on mixing for less-deformable droplets can be explained physically by an interplay between surface deformation velocity and inner advection. For higher values of $Ca$, drop deformation happens faster than inner advection, resulting in a less effective mixing at larger scales. Hence, although drop deformation breaks the kinematic reversibility of the Stokes flow inside the droplet, larger deformations can decrease mixing, as observed previously for passive mixers in Muradoglu & Stone (Reference Muradoglu and Stone2005).

Figure 13. Numerical simulations of mixing inside a droplet undergoing a three-phase extensional flow for $a = 0.4$, $Ca = 0.05$, $H = 1$, $\omega = 3$ and $\lambda = 1$ at different times. The results are for the midplane $z = 0$. The number below each droplet is the mixing number $m(t)$, calculated using (4.2).

Another way to combine the different modes to enhance active mixing is to alternate between different modes, as is usually done in passive mixing (e.g. serpentine channels) and in the investigation of Stone & Stone (Reference Stone and Stone2005), where a spherical droplet was subjected to alternating uniform and shear flows. In our system, one possibility is to alternate between three-phase extension mode and the tri-axial extension. To illustrate this improvement, figure 14 shows numerical results for the mixing inside a droplet for $Ca = 0.1$, $\lambda = 1$, $\omega = 3$, and external flow given by

(4.5) \begin{equation} \boldsymbol{Q}_0(t) = \left\{\begin{array}{@{}ll@{}} \boldsymbol{Q}_{{rotor}}, & \text{for}\ t \le 2 {\rm \pi}/\omega \pmod{4 {\rm \pi}/\omega}, \\ \boldsymbol{Q}_{{tri}}, & \text{for}\ t > 2 {\rm \pi}/\omega \pmod{4 {\rm \pi}/\omega}. \end{array} \right. \end{equation}

The results shown in figure 14 indicate that even for viscosity ratios of ${O}(1)$, it is possible to get a more effective mixing by alternating between equal periods of the two flow modes. Besides the mixing simulations for the calculation of the mixing number, we included animations of trajectories of multiple tracer particles for different modes that can be found in the supplementary material.

Figure 14. Numerical simulation of mixing inside a droplet for an external flow alternating between three-phase extension and tri-axial extension modes for $Ca = 0.1$, $\lambda = 1$, $H = 1$ and $\omega = 3$ for the midplane $z = 0$. The number below each droplet is the mixing number $m(t)$, calculated using (4.2).