A singularly altered streamline topology allows faster transport from deformed drops

Abstract

We analyse the effect of drop-deformation-induced change in streamline topology on the scalar transport rate (the Nusselt number $Nu$) in an ambient planar linear flow. The drop-phase resistance is assumed dominant, and the drop deformation is characterised by the capillary number ($Ca$). For a spherical drop ($Ca = 0$) in an ambient planar extension, closed streamlines lead to $Nu$ increasing with the Péclet number ($Pe$), from $Nu_0$, corresponding to purely diffusive transport, to $4.1Nu_0$, corresponding to a large-$Pe$ diffusion-limited plateau. For non-zero $Ca$, we show that the flow field consists of spiralling streamlines densely wound around nested tori foliating the deformed drop interior. Now $Nu$ increases beyond the aforementioned primary plateau, saturating in a secondary one that approaches $22.3Nu_0$ for $Ca \rightarrow 0$, $Pe\,Ca \rightarrow \infty$. The enhancement appears independent of the drop-to-medium viscosity ratio. We further show that this singular dependence, of the transport rate on drop deformation, is generic across planar linear flows; chaotically wandering streamlines in some of these cases may even lead to a tertiary enhancement regime.

1. Introduction

Transport of heat and mass in multiphase scenarios is relevant to diverse natural phenomena: growth of cloud condensation nuclei (Kinzer & Gunn 1951; Beard & Pruppacher 1971; Duguid & Stampfer 1971), nutrient uptake by microswimmers (Magar, Goto & Pedley 2003; Guasto, Rusconi & Stocker 2012; Stocker 2012) and industrial processes (fuel atomisation in internal combustion (IC) engines (Law 1982), spray drying (Patel, Patel & Suthar 2009), liquid–liquid extraction (Wegener, Paul & Kraume 2014) and suspension polymerisation (Vivaldo-Lima et al. 1997; Brooks 2010)). The underlying aim is to calculate the transport rate from a drop or a particle immersed in an ambient shearing flow. In non-dimensional terms, this amounts to determining the Nusselt ($Nu$) or Sherwood ($Sh$) number, the rate of heat and mass transport measured in diffusion units, respectively, as a function of the Péclet number ($Pe$). The latter is defined as $Pe = \dot {\gamma } a^2/{\mathcal {D}}$, and compares the relative magnitudes of convection and diffusion; $\dot {\gamma }$ here is a characteristic shear rate, $a$ the particle or (undeformed) drop radius and ${\mathcal {D}}$ the thermal or mass diffusivity. Mass transport in liquids usually occurs in the convection dominant limit $(Pe \gg 1):Pe\sim O(10^3)$ for small molecules in sheared aqueous environments, but can be several orders of magnitude larger for macromolecular solutes in viscous solvents. Herein, we show, for a drop in an ambient planar linear flow, that deformation-induced alteration of the interior streamline topology has a singular effect on the transport rate in the convection dominant limit. For sufficiently large $Pe$, transport from even a weakly deformed drop is enhanced more than five-fold, relative to a spherical one, for an ambient planar extension, and by nearly an order of magnitude for an ambient simple shear flow.

For rigid particles, transport for large $Pe$ occurs across a thin boundary layer, with $Nu$ scaling as the inverse boundary layer thickness. A balance of convection and diffusion time scales yields a boundary layer thickness of $O(a\, Pe^{-{1}/{3}})$, so $Nu \propto Pe^{{1}/{3}}$ (Leal 2007). For drops, when the ambient phase resistance to transport is dominant, the so-called exterior problem, the boundary layer thickness is $O(a\, Pe^{-{1}/{2}})$ due to the interfacial slip and $Nu \propto Pe^{{1}/{2}}$, being larger than that for a rigid particle at the same $Pe$. Transport for large $Pe$ also depends sensitively on the streamline topology (Subramanian & Koch 2006a,c; Krishnamurthy & Subramanian 2018b), and boundary-layer-driven enhancement above occurs only when the exterior flow has an open-streamline (Acrivos & Goddard 1965; Gupalo & Riazantsev 1972; Gupalo, Riazantsev & Ulin 1975; Polyanin 1984; Leal 2007; Krishnamurthy & Subramanian 2018a,b) or open-pathline (Banerjee & Subramanian 2021) topology. Owing to the unbounded domain, such a topology is a common occurrence for the exterior problem, either in the Stokes limit itself (Krishnamurthy & Subramanian 2018a; Banerjee & Subramanian 2021) or due to deviations from Stokesian hydrodynamics (Subramanian & Koch 2006b,c, 2007; Krishnamurthy & Subramanian 2018b). In contrast, a closed-streamline topology is more common for a confined domain such as the spherical drop interior. While a generic ambient linear flow leads to chaotically wandering streamlines in the Stokes limit (Stone, Nadim & Strogatz 1991; Sabarish 2021), many canonical flows including uniform flow, linear extensional flows and the family of planar linear flows lead to closed interior streamlines (Sabarish & Subramanian 2023). Earlier efforts analysing transport from a spherical drop in these flows, in the limit of a dominant drop-phase resistance (the interior problem), have shown $Nu$ to plateau at an order unity value for $Pe \rightarrow \infty$ (Newman 1931; Kronig & Brink 1951; Johns & Beckmann 1966; Watada, Hamielec & Johnson 1970; Brignell 1975; Prakash & Sirignano 1978), implying diffusion-limited transport due to the closed-streamline topology.

Herein, for the first time, we consider the interior problem for a deformed drop in an ambient linear flow. While there has been extensive research on drop deformation and breakup in linear flows over the last several decades (Taylor 1934; Rumscheidt & Mason 1961; Torza et al. 1971; Hakimi & Schowalter 1980; Rallison 1984; Bentley & Leal 1986b; Stone 1994), substantially less attention has been paid to the correlation between drop deformation and streamline topology. The importance of drop deformation is characterised by the capillary number, $Ca = \mu \dot {\gamma }a/T$, $\mu$ being the ambient fluid viscosity and $T$ the interfacial tension coefficient. For the interior problem in particular, one has closed streamlines and diffusion-limited transport for $Ca = 0$ (spherical drop), on account of Stokesian reversibility constraints (Newman 1931; Kronig & Brink 1951; Johns & Beckmann 1966; Watada et al. 1970; Brignell 1975; Prakash & Sirignano 1978), but this must no longer be true for finite $Ca$. However, earlier experiments (Torza et al. 1971) and computations (Kennedy, Pozrikidis & Skalak 1994; Komrakova et al. 2014) have suggested the persistence of closed streamlines within a deformed drop in simple shear flow. In contrast, we show, using the analytical velocity field for small $Ca$, and supporting finite-$Ca$ boundary integral computations, that drop deformation opens up closed streamlines. The aforementioned sensitive dependence of transport on streamline topology implies then that the $Nu$–$Pe$ relationship is qualitatively affected by the altered streamline topology. For an ambient planar extension, Langevin simulations show that the drop-deformation-induced spiralling-streamline topology leads to singularly enhanced transport, for $Pe \rightarrow \infty$, for any non-zero $Ca$. The singular enhancement persists for other ambient planar linear flows, being robust to the addition of vorticity.

This article is organised as follows. In § 2.1, we write down the equations governing the motion within and outside a neutrally buoyant drop immersed in an ambient shearing flow. Solving these equations for a weakly deformed drop, using a domain perturbation technique, leads to analytical expressions for the velocity and pressure fields which are given in Appendix A. The analytical interior velocity field is used to solve for the scalar (temperature or concentration) field in order to determine the rate of scalar transport. This is accomplished using Langevin simulations, details of which are given in § 2.2. Langevin simulations using the analytical velocity field, valid for weak deformations, are complemented by those that use the velocity field for finitely deformed drops, determined using boundary integration computations; details of these computations are given in § 2.3 and in Appendix B. Section 3.1 discusses the interior streamline topology within, and scalar transport from, a spherical drop in an ambient planar extensional flow. The deformation-induced alteration of the streamline topology, and the accompanying transport enhancement, are examined in § 3.2. Sections 4.1 and 4.2 examine the interior streamline topology and scalar transport rate for spherical and deformed drops, respectively, in other ambient planar linear flows. In § 5, we summarise our findings, while also emphasising the significance of the transport enhancement, found here for drops, in the context of other deformable microstructures.

2. Problem formulation

2.1. Governing equations and the small $Ca$ expansion

Consider a neutrally buoyant Newtonian drop of viscosity $\lambda \mu$ in an ambient Newtonian fluid of viscosity $\mu$ undergoing a planar linear flow, $\boldsymbol {u}^\infty = \boldsymbol {\varGamma } \boldsymbol{\cdot } \boldsymbol {x}$, with $\boldsymbol {\varGamma }=2\dot {\gamma }\left [\begin{smallmatrix} 0 & 1 & 0 \\ \beta & 0 & 0\\ 0 & 0 & 0 \end{smallmatrix}\right ]$ being the (transpose of the) velocity gradient tensor. The associated rate of strain and vorticity tensors are given by $\boldsymbol {E}=\dot {\gamma }(1+\beta )\left [\begin{smallmatrix} 0 & 1 & 0 \\ 1 & 0 & 0\\ 0 & 0 & 0 \end{smallmatrix}\right ]$ and $\boldsymbol {\varOmega }=\dot {\gamma }(1-\beta )\left [\begin{smallmatrix} 0 & 1 & 0 \\ -1 & 0 & 0\\ 0 & 0 & 0 \end{smallmatrix}\right ]$, respectively. Planar linear flows are a one-parameter family with $\beta \in [-1,1]$; $\beta = -1$, $0$ and $1$ correspond to solid-body rotation, simple shear and planar extension, respectively. The interval $-1 \leq \beta < 0$ corresponds to elliptic linear flows with closed streamlines, and $0 < \beta \leq 1$ corresponds to hyperbolic planar linear flows with open streamlines; see figure 1. The schematic in figure 2 shows a drop, slightly deformed from its native spherical form (of radius $a$), in an ambient hyperbolic planar linear flow. While our focus in this work is on planar extensional flow ($\beta = 1$), we also discuss the streamline topology for non-unity $\beta$ in § 4.2, along with the implications for transport.

Figure 1. Members of the one-parameter planar linear flow family: (a) solid-body rotation ($\beta =-1$), (b) an elliptic linear flow with $\beta =-0.5$, (c) a hyperbolic planar linear flow with $\beta =0.5$ and (d) planar extension ($\beta =1$).

Figure 2. Schematic of a slightly deformed drop (of undeformed radius $a$) suspended in an ambient hyperbolic planar linear flow.

Assuming viscous effects to be dominant, the governing equations, in non-dimensional form, are

(2.1)$$\begin{gather} \lambda \nabla^{2}\hat{\boldsymbol{u}}-\boldsymbol{\nabla} \hat{p}=0, \quad \boldsymbol{\nabla} \boldsymbol{\cdot } \hat{\boldsymbol{u}}=0, \end{gather}$$

(2.2)$$\begin{gather}\nabla^{2}\boldsymbol{u}-\boldsymbol{\nabla} p=0, \quad \boldsymbol{\nabla} \boldsymbol{\cdot } \boldsymbol{u}=0, \end{gather}$$

with $\hat {\boldsymbol {u}} (\boldsymbol {u})$ and $\hat {p} (\,p)$ being the interior (exterior) velocity and pressure fields. Here, $a$, $\dot {\gamma }a$ and $\mu \dot {\gamma }$ have been used as the characteristic length, velocity and stress scales, respectively. Equations (2.1)–(2.2) are subject to the following boundary conditions:

(2.3)\begin{equation} \left.\begin{array}{l} \boldsymbol{u} \rightarrow \boldsymbol{u}^{\infty} \quad \text{for} \ \boldsymbol{x} \rightarrow \infty, \\ \boldsymbol{u}=\hat{\boldsymbol{u}}\ (\text{Velocity continuity at the interface}),\\ \dfrac{1}{| \boldsymbol{\nabla} F |} \dfrac{\partial F}{\partial t }+ Ca (\boldsymbol{u} \boldsymbol{\cdot } \boldsymbol{n})=0\ (\text{Kinematic boundary condition at the interface}),\\ (\boldsymbol{\sigma}-\hat{\boldsymbol{\sigma}})\boldsymbol{\cdot }\boldsymbol{n} =\dfrac{1}{Ca}(\boldsymbol{\nabla} \boldsymbol{\cdot } \boldsymbol{n})\boldsymbol{n}\ (\text{Interfacial stress balance}).\\ \end{array}\right\} \end{equation}

Here, $\boldsymbol {n}$ is the outer unit normal to the drop surface, $\hat {\boldsymbol {\sigma }}=-\hat {p} \boldsymbol {I}+ \lambda (\boldsymbol {\nabla } \hat {\boldsymbol {u}}+\boldsymbol {\nabla } \hat {\boldsymbol {u}}^{\textrm {T}})$ and $\boldsymbol {\sigma }=-p \boldsymbol {I}+ (\boldsymbol {\nabla } \boldsymbol {u}+\boldsymbol {\nabla } \boldsymbol {u}^{\textrm {T}})$ are the interior and exterior stress tensors, respectively, and $t$ in the kinematic condition is non-dimensionalised by $\dot {\gamma }^{-1}$. Under the action of the imposed flow, the drop deforms, attaining a steady shape for $Ca$ less than a critical value $Ca_{c}$ ($Ca_{c} \simeq 0.06$ for $\lambda =1$ and an ambient planar extension, and is higher for other members of the planar linear flow family; see Rallison 1981; Bentley & Leal 1986a). For small $Ca$, one may describe the deformed drop shape using a scalar function, $F(\boldsymbol {x},t)=r-(1+Ca\,f (\boldsymbol {n},t ))=0$, to $O(Ca)$. It is then straightforward, albeit tedious, to calculate the disturbance fields and perturbed drop shape for $Ca$, $\lambda Ca \ll 1$, using a domain perturbation technique (Leal 2007), as has been done across several earlier efforts (Cox 1969; Barthes-Biesel & Acrivos 1973; Rallison 1980; Leal 2007; Vlahovska, Bławzdziewicz & Loewenberg 2009; Ramachandran & Leal 2012). The results, to $O(Ca)$, have been tabulated in Appendix A.

2.2. Langevin simulations for the interior problem

The scalar transport for the interior problem of interest is governed by the convection–diffusion equation for the scalar concentration field, with convection being driven by the Stokesian interior velocity field:

(2.4)\begin{equation} \frac{\partial c}{\partial \tau}+ Pe\,\hat{\boldsymbol{u}}(\boldsymbol{x}) \boldsymbol{\cdot } \boldsymbol{\nabla}c= \nabla^{2} c. \end{equation}

In (2.4), velocity and time are scaled by $\dot {\gamma }a$ and $a^2 /\mathcal {D}$, respectively, $\mathcal {D}$ being the diffusivity coefficient; the Péclet number is defined as $Pe=\dot {\gamma } a^2 / \mathcal {D}$; $Pe\gg 1$ denotes the convection dominant limit. In our simulations, $\hat {\boldsymbol {u}}(\boldsymbol {x})$ in (2.4) is either the small-$Ca$ analytical form, $\hat {\boldsymbol {u}}^{(0)}(\boldsymbol {x}) +Ca\,\hat {\boldsymbol {u}}^{(1)}(\boldsymbol {x})$ discussed in § 2.1 above, with expressions for $\hat {\boldsymbol {u}}^{(0)}(\boldsymbol {x})$ and $\hat {\boldsymbol {u}}^{(1)}(\boldsymbol {x})$ given in Appendix A; or obtained numerically using the boundary integral method described in § 2.3. With $\hat {\boldsymbol {u}}(\boldsymbol{ x})$ obtained in this manner, (2.4) is solved starting from an initially uniform distribution, with the drop surface acting as a perfectly absorbing boundary. The suitably normalised concentration satisfies $c=1$ within the drop at $\tau =0$, and $c=0\ \forall \tau \geq 0$ at the surface of the deformed drop which, for small $Ca$, is as defined in § 2.1. The latter absorption boundary condition corresponds to the complete neglect of any transport resistance in the exterior fluid, which in turn arises from assuming the scalar diffusivity within the drop to be much smaller than that in the ambient fluid. In this limit, one may neglect the asymptotically small change in concentration (from the surface to the ambient value) in the drop exterior, and this, in turn, is equivalent to setting $c$ to be zero at the drop surface, as mentioned previously.

The interior streamlines, for the planar linear flows under consideration, do not conform to any obvious symmetries, precluding the use of standard coordinate systems. Therefore, we use Langevin simulations to simulate the transport process that, on average, is described by (2.4). The Langevin equations (Gardiner 2004) for the individual tracers are

(2.5)\begin{equation} {\rm d}\kern0.7pt \boldsymbol{x}= Pe\,\hat{\boldsymbol{u}}(\boldsymbol{x})\,{\rm d}\tau+ \sqrt{2}\,{\rm d}\boldsymbol{W}(\tau). \end{equation}

Here, $\hat {\boldsymbol {u}}(\boldsymbol {x})$ convects the tracer particle at position $\boldsymbol {x}(\tau )$, and the final term, $d\boldsymbol {W}(\tau )$, corresponds to the standard Wiener process satisfying $\langle \textrm {d}\boldsymbol {W}(\tau )\rangle =0$, $\langle (\textrm {d}\boldsymbol {W})^2\rangle =\textrm {d} \tau$. We solve the Langevin equations using a fourth-order Runge–Kutta technique with a time step that satisfies $Pe\,\textrm {d}\tau \ll 1$, ensuring numerical convergence (Wilkie 2004). In the numerical implementation, $\textrm {d}\boldsymbol {W}(\tau ) = \sqrt {\textrm {d}\tau } N(0,1)$, $N(0,1)$ being a Gaussian random variable with zero mean and unit variance. For a sufficiently large initial number of tracers, the ensemble average of the tracer concentration converges to the solution of the convection–diffusion equation. In all our runs, we start from $10^5$ tracers uniformly distributed within the drop, with the absorbing boundary condition implemented at each time step by removing any tracers that end up crossing the drop surface. To address the decrease in the number of tracers with time, and the associated increase in the amplitude of the tracer number fluctuations, we implemented a replenishment protocol where, when the tracers reduce to half the initial number, an equal number of new tracers are introduced at the same locations as the original ones. This ensures that the tracer number never dips below 50 000 during the course of any simulation, in turn ensuring accurate long-time averages required for the $Nu$ calculation (see the following). Further, the manner of introduction also ensures that the addition of new tracers only affects the amplitude of the concentration profile, and not its shape.

Using the time-dependent volume-averaged concentration obtained from the simulations, we calculate the non-dimensional transport rate as $Nu=({2}/{\bar {c}}) F$ (Kronig & Brink 1951; Juncu 2010), where $\bar {c}(\tau )= ({1}/{V_d})\int _{V_d} c(\boldsymbol {x},\tau )\,\textrm {d}V$ is the volume-averaged concentration and $F=-({1}/{4{\rm \pi} }) \int _{S} \boldsymbol {n} \boldsymbol{\cdot } \boldsymbol {\nabla } c(\boldsymbol {x},\tau )\,\textrm {d}S$ is the surface averaged flux, with the area of the undeformed drop ($4{\rm \pi}$) used for normalisation; incompressibility implies that the drop volume $V_d= {4{\rm \pi} }/{3}$. We argue later that this definition, which is also that used in the earlier literature (see, for instance, Kronig & Brink 1951) yields the same $Pe$-dependence for the transport rate, as would be obtained from the usual definition, based on the normal derivative of the concentration at the drop surface, when a true steady state is made possible due to a compensating source field. The choice of the undeformed-drop-based normalisation above is because, on one hand, surface areas of the deformed drops remain close to the spherical value over the range of $Ca$ examined; even for the highest $Ca\ (=0.04)$, the deformed drop area exceeds the spherical value only by $1.3\,\%$, implying the transport enhancements found below owe their origin to the change in streamline topology, and not an increased surface area. On the other hand, a normalisation that accounts for the drop deformation would involve solving for the transport rate, for the fictitious problem of diffusion from a deformed drop.

Using the divergence theorem, the above expression for $F$ is written as $F=- ({1}/{4{\rm \pi} }) \int _{V_d} \nabla ^2 c(\boldsymbol {x},\tau )\,\textrm {d}V$, over the volume of the drop ($V_d$), and further use of the convection–diffusion equation reduces this to $F=- ({1}/{4{\rm \pi} }) \int _{V_d} ({\partial c}/{\partial \tau })(\boldsymbol {x},\tau )\,\textrm {d} V$. Since the drop attains a steady non-spherical shape with volume equal to the original undeformed sphere (${4{\rm \pi} }/{3}$), one may write

(2.6)\begin{equation} F=- \frac{1}{3} \frac{{\rm d}}{{\rm d}\tau} \left[ \frac{1}{V_d} \int_{V_d} c(\boldsymbol{x},\tau) \,{\rm d}V\right] =- \frac{1}{3} \frac{\partial \bar{c}}{\partial \tau}. \end{equation}

Hence, one finally obtains

(2.7)\begin{equation} Nu=- \frac{2}{3} \frac{1}{\bar{c}} \frac{{\rm d}\bar{c} }{{\rm d}\tau} =- \frac{2}{3} \frac{{\rm d}\log \bar{c}}{{\rm d}\tau}. \end{equation}

Since the absorption boundary condition renders the interior problem an inherently unsteady one, both $c(\boldsymbol {x},\tau )$ and $Nu$ are in general functions of time ($\tau$). It is important, however, to note that, for all cases examined, $c(\boldsymbol {x},\tau )$ approaches a simple exponential decay for sufficiently long times, and correspondingly, $Nu$ approaches a time-independent plateau value that equals two-thirds of the largest (least negative) eigenvalue of the convection–diffusion operator. In each of our simulations, therefore, $Nu$ exhibits an algebraic decrease proportional to $\tau ^{-{1}/{2}}$ for short times, independent of $Pe$, corresponding to an initial diffusive transport across a thin boundary layer of thickness $\tau ^{{1}/{2}} \approx ( \mathcal {D}t )^{{1}/{2}}$ (see the next paragraph), but eventually asymptotes to a $Pe$-dependent plateau with superposed small amplitude fluctuations that arise from the finite number of tracers simulated. All of the $Nu$-values used later in this paper, in the context of the $Nu$–$Pe$ relationships, correspond to time-independent estimates, obtained by taking a decade-long time average in this plateau regime.

To validate our simulations, we consider the analytically solvable pure diffusion problem corresponding to $Pe=0$. For this case, the time-dependent concentration field is given by $c_0(\boldsymbol {x},\tau )=({2}/{{\rm \pi} r}) \sum _{n=1}^{\infty } ({(-1)^{n+1}}/{n}) \sin (n{\rm \pi} r) \,\textrm {e}^{-n^{2}{\rm \pi} ^{2} \tau }$ (Carslaw & Jaeger 1948; Kronig & Brink 1951), with the volume-averaged concentration given by $\bar {c}_0(\tau )=6 \sum _{n=1}^{\infty } {\textrm {e}^{-n^2 {\rm \pi}^2 \tau }}/{n^2 {\rm \pi}^2}$; accordingly, $\lim _{\tau \to \infty } Nu_{0}=(2/3) {\rm \pi}^{2}$, the subscript ‘$0$’ denoting the diffusion limit. For sufficiently short times, the diffusive transport occurs across a thin $O(\tau ^{{1}/{2}})$ concentration boundary layer just beneath the drop surface that may now be approximated as an infinite plane; this yields $Nu_0(\tau ) \approx 2/\sqrt {{\rm \pi} \tau }$ for $\tau \ll 1$. Figure 3 shows $Nu_0(\tau )$ obtained from the Langevin simulations, alongside both the exact analytical solution and the short-time asymptote above, and highlights the excellent agreement between them. As mentioned previously, the short-time asymptote continues to be valid for all $Pe$, although for large $Pe$, it gives way to convection-induced oscillations for $\tau \gtrsim O(Pe^{-1})$ (see the inset in figure 4b).

Figure 3. Comparison of $Nu_0(\tau )/Nu_{0}(\tau \rightarrow \infty )$ obtained from the Langevin simulations for $Pe=0$ with the exact analytical solution, $(2/3 ) (\sum _{n=1}^{\infty } \textrm {e}^{-n^2 {\rm \pi}^2 \tau } )/( \sum _{n=1}^{\infty } {\textrm {e}^{-n^2 {\rm \pi}^2 \tau }}/{n^2 {\rm \pi}^2})$ (Kronig & Brink 1951), and the corresponding short-time asymptote $3/({\rm \pi} ^2\sqrt {{\rm \pi} \tau })$.

Figure 4. (a) Closed streamlines in a spherical drop octant ($Ca=0$, $\lambda =1$), in an ambient planar extension, along with the corresponding in-plane streamlines (black); the black dotted fixed-point circle corresponds to $E=-0.1628$, $D\in [0,0.5706]$. Inset (i) shows the inplane streamline pattern; inset (ii) shows all streamlines as ordered pairs on the $D$–$E$ plane (same colour code as the main plot). (b) Plot of $\lim _{\tau \rightarrow \infty } Nu/Nu_0$ vs $Pe$ for planar extension ($Ca = 0$, $\lambda =1$); the $Nu(\tau )/Nu_0$-curves appear in an inset.

2.3. Boundary element method

The boundary integral formulation for the Stokes equations expresses the velocity field at any point in terms of distributions of point forces (the single-layer potential) and point-force dipoles/sources (the double-layer potential) on the surface of the drop. For an interior point $\boldsymbol {x}_{0}$, it may be written in non-dimensional form as (Pozrikidis 1992, 2002; Kennedy et al. 1994)

(2.8)\begin{align} u_{j}(\boldsymbol{x}_{\boldsymbol{0}})&=\frac{Ca}{\lambda} u_{j}^{\infty}(\boldsymbol{x}_{\boldsymbol{0}}) -\frac{1}{8 {\rm \pi}\lambda} \int_{S} \Delta f_{i}(\boldsymbol{x}) G_{ij}(\boldsymbol{x},\boldsymbol{x}_{\boldsymbol{0}}) \,{\rm d}S(\boldsymbol{x})\nonumber\\ &\quad + \frac{1-\lambda}{8 {\rm \pi}\lambda} \int_{S} u_{i}(\boldsymbol{x}) T_{ijk}(\boldsymbol{x}, \boldsymbol{x}_{\boldsymbol{0}}) \boldsymbol{n}_{k}(\boldsymbol{x}) \,{\rm d}S(\boldsymbol{x}). \end{align}

In the above relation, $S$ is the surface of the deformed drop and $\boldsymbol {n}$ is the unit normal pointing into the ambient fluid. The first term on the right-hand side corresponds to the imposed flow $\boldsymbol {u}^{\infty }(\boldsymbol {x})= \boldsymbol {\varGamma } \boldsymbol{\cdot } \boldsymbol {x}$ with $\boldsymbol {\varGamma }$ defined in § 2.1. In the second term, $\Delta f_{i}(\boldsymbol {x})=(\sigma _{ij}- \hat {\sigma }_{ij}) n_{j}$ is the jump in the interfacial force density, while the kernels of the single- and double-layer potentials, $\boldsymbol {G}$ and $\boldsymbol {T}$, are the Stokeslet and the associated stress tensor, given by $G_{ij}(\boldsymbol {x},\boldsymbol {x}_{\boldsymbol {0}})={\delta _{ij}}/{|\boldsymbol {\hat {x}}|}+{\hat {x}_{i} \hat {x}_{j}}/{|\boldsymbol {\hat {x}}|^{3}}$ and $T_{ijk}(\boldsymbol {x},\boldsymbol {x}_{\boldsymbol {0}})=-6({\hat {x}_{i} \hat {x}_{j} \hat {x}_{k}}/{|\hat {\boldsymbol {x}}|^{5}})$, where $\hat {\boldsymbol {x}}=\boldsymbol {x}-\boldsymbol {x}_{\boldsymbol {0}}$. Combining the boundary integral representations for an interior and exterior point $\boldsymbol {x}_0$ (the latter being given by (2.8) but without the $\lambda$) and, then, suitably combining the limiting forms of these representations for $\boldsymbol {x}_{0}$ approaching the surface of the drop, one obtains the following Fredholm integral equation of the second kind governing the interfacial velocity field:

(2.9)\begin{align} u_{j}(\boldsymbol{x}_{\boldsymbol{0}})&=\frac{2 Ca }{1+\lambda} u_{j}^{\infty}(\boldsymbol{x}_{\boldsymbol{0}}) -\frac{1}{4 {\rm \pi}(1+\lambda)} \int_{S} \Delta f_{i}(\boldsymbol{x}) G_{ij}(\boldsymbol{x},\boldsymbol{x}_{\boldsymbol{0}}) \,{\rm d}S(\boldsymbol{x})\nonumber\\ &\quad + \frac{1-\lambda}{4 {\rm \pi}(1+\lambda)} \int_{S}^{PV} u_{i}(\boldsymbol{x}) T_{ijk}(\boldsymbol{x}, \boldsymbol{x}_{\boldsymbol{0}}) \boldsymbol{n}_{k}(\boldsymbol{x}) \,{\rm d}S(\boldsymbol{x}), \end{align}

where $PV$ signifies the improper integral that results from the double-layer potential for $\boldsymbol {x}_{0}$ on $S$ (Pozrikidis 1992). For $\lambda =1$, (2.9) reduces to a much simpler explicit relation for the velocity field:

(2.10)\begin{equation} u_{j}(\boldsymbol{x}_{\boldsymbol{0}})=Ca\, u_{j}^{\infty}(\boldsymbol{x}_{\boldsymbol{0}}) -\frac{1}{8 {\rm \pi}} \int_{S} \Delta f_{i}(\boldsymbol{x}) G_{ij}(\boldsymbol{x},\boldsymbol{x}_{\boldsymbol{0}}) \,{\rm d}S(\boldsymbol{x}). \end{equation}

For finite $Ca$, we obtain the interior velocity field by numerically solving either (2.9) or (2.10), and then use it as an input for the Langevin simulations to obtain $Nu$.

The boundary element method (BEM) simulations start from an initially spherical drop, with this surface discretised into a large number of triangular elements. Each element is assigned an array of six marker points (nodes), leading to a system of linear equations for the nodal velocity fields which has to be solved at every time step, for $\lambda \neq 1$. The nodal velocity field values are then used in evaluating the surface integrals in (2.9) or (2.10) and, thence, the interfacial velocity field. One then uses the latter to evolve the surface, the time step chosen for this evolution being small enough for numerical convergence. This eventually leads to a steadily deformed surface for $Ca$ values less than the threshold for breakup. Note that the aforementioned procedure is along the same lines as in earlier efforts; see, for instance, Pozrikidis (1992, 2002) and Kennedy et al. (1994). Once steady state is achieved, the velocity field at a given interior point may be obtained by using the original boundary integral representation, (2.8); the contribution of each triangular element to the surface integrals in (2.8) is accounted for via an interpolation process. It is this (appropriately interpolated) interior velocity field, associated with the steadily deformed drop, that is then used to convect the tracers in the Langevin simulations, with the absorbing boundary condition being implemented at the deformed drop boundary. It is worth mentioning that the streamline topology for small but finite $Ca$ is particularly sensitive to the number of surface elements, requiring sufficiently high resolution for accurate results. The lack of resolution is very likely the reason the deformation-induced change in streamline topology went unnoticed in earlier efforts (Rallison 1981; Kennedy et al. 1994). We choose a much larger number of surface elements compared to the said efforts, and for each such choice, two sets of interpolation points ($13$ and $19$). Additional information regarding the BEM computations is given in Appendix B.

3. Scalar transport in an ambient planar extensional flow

3.1. Transport from a spherical drop

The interior velocity field $\hat{\boldsymbol{u}}(\boldsymbol{x})$, given in Appendix A, may be used to characterise the interior streamline topology. For $Ca = 0$, almost all interior streamlines are closed, and may be obtained as the curves of intersection of two one-parameter families of invariant surfaces (Cox, Zia & Mason 1968; Torza et al. 1971; Powell 1983; Krishnamurthy & Subramanian 2018a,b). Denoting the parameters as $D$ and $E$, the invariant surfaces are $x_{3}=D(1-r^2)^{-{1}/{3}}$ and

(3.1)\begin{equation} x_{2}={\pm} r\left( \frac{1}{2}-\frac{(1-\beta)(1+\lambda)}{2 r^2(1+\beta)}+\frac{E}{r^2(1-r^2)^{2/3}}\right)^{{1}/{2}}, \end{equation}

with $-0.1628 \leq E \leq 0.1628$ and $-0.5706 \leq D \leq 0.5706$; for $\beta = 1$, both families are independent of $\lambda$. Each ordered pair $(D,E)$ corresponds to a closed streamline (Krishnamurthy & Subramanian 2018a), the exceptions being fixed (stagnation) points, and streamlines connecting them. Figure 4(a) depicts the closed streamlines in an octant of a spherical drop, along with the quarter-circle of fixed points that ‘threads’ these streamlines, for an ambient planar extensional flow. The in-plane streamlines ($D=0$) in the corresponding quadrant are also shown (in black). The streamlines in other octants may be obtained by symmetry; the quarter-circle, when extended thus, leads to two orthogonally oriented fixed-point circles of identical radii that intersect at $(0,0,\pm \sqrt {{3}/{5}})$.

Figure 4(b) shows $Nu/Nu_0$ in the long-time limit as a function of $Pe$, for $Ca = 0$ and $\lambda = 1$, for an ambient planar extension; $Nu_0 = (2/3){\rm \pi} ^2$ being the long-time diffusive rate of transport ($Pe = 0$). The normalised transport rates for arbitrary times are shown in an inset. For small $Pe$, the time-dependent curves in the inset decay monotonically to a long-time plateau, reflecting the self-adjointness of the diffusion operator (real and negative eigenvalues). For large $Pe$, the decay has an oscillatory character, reflecting transient adjustment of the iso-scalar contours and streamlines on time scales of $O(\dot {\gamma }^{-1}Pe^{{1}/{3}})$, driven by shear-enhanced diffusion (Christov & Homsy 2009; Juncu 2010). The main figure shows $\lim _{\tau \to \infty } Nu/Nu_0$ increasing from unity to a plateau value of $4.1$ for $Pe \gtrsim O(100)$. The plateau arises because tracers are rapidly convected around closed streamlines for larger $Pe$, with the transport controlled by the much slower rate at which they diffuse across these streamlines. Thus, the plateau value above corresponds to the rate of diffusion across closed streamlines that are coincident with iso-scalar contours (Frankel & Acrivos 1968; Acrivos 1971). It is worth mentioning that one may, in principle, derive a streamline-averaged diffusion equation in the coordinates $D$ and $E$, that governs the scalar concentration for $Pe \rightarrow \infty$, and whose solution should yield the 2$D$-diffusion-limited plateau; however, the streamline-averaged diffusivity tensor has to be determined numerically even for the simplest scenario of an ambient uniform flow (Kronig & Brink 1951).

The $\lambda$-independence of the streamline topology for planar extension implies that the effect of $\lambda$ may be simply accounted for by replacing $Pe$ by $Pe(1+\lambda )^{-1}$, as in figure 4(b). While the interior streamline topology exhibits a non-trivial dependence on $\lambda$ for other planar linear flows (Powell 1983; Krishnamurthy & Subramanian 2018b), precluding the above use of a rescaled $Pe$, almost all interior streamlines are still closed for $Ca = 0$, and as a result, $Nu$ exhibits a $Pe$-dependence analogous to planar extension (see figure 8). The large-$Pe$ diffusion-limited plateau equals $1.08$ for $\lambda = 1$ for simple shear; as expected, $Nu/Nu_0=1$ for solid-body rotation, regardless of $Pe$ or $\lambda$.

3.2. Transport from a deformed drop

Figure 5(a) depicts a pair of streamlines within a weakly deformed drop ($Ca$ small but finite) in an ambient planar extension. Each streamline has a tightly spiralling character, winding densely around an invariant torus for sufficiently long times, with these tori foliating the interior of the drop octant. The outer limiting torus comprises the surfaces of the deformed octant (that include portions of the $x_1\pm x_2=0$ planes), and a perturbed form of the quarter-circle of fixed points in figure 4(a); the inner limiting torus is a closed curve (shown in green in figure 5a) that loops around the perturbed quarter-circle. Here $D$ and $E$ are no longer constant along a finite-$Ca$ streamline, but instead vary on a time scale of $O(Ca^{-1}\dot {\gamma }^{-1})$ that, for small $Ca$, is much longer than the $O(\dot {\gamma }^{-1})$ period of circulation around the zero-$Ca$ closed streamlines. This two-time-scale structure is illustrated in figure 5(b), where small-amplitude fast oscillations are seen to be superposed on a slow $O(Ca)$ drift in the $D$–$E$ plane. The amplitude of these oscillations decreases with decreasing $Ca$, as illustrated via the pair of insets in figure 5(b). For $Ca \rightarrow 0$, but for times much greater than $O(Ca^{-1}\dot {\gamma }^{-1})$, the streamlines must asymptote to smooth closed contours in the $D$–$E$ plane. These contours may, in principle, be obtained from solving an autonomous system of ordinary differential equations (ODEs) in $(D,E)$, derived using the method of averaging, implying the existence of an adiabatic invariant for small but finite $Ca$. The latter is a function of $D$ and $E$ that remains constant over times of $O(Ca^{-1}\dot {\gamma }^{-1})$ and, therefore, serves to parameterise the said contours: the method of averaging has been applied in Sabarish & Subramanian (2023) to the exterior finite-$Ca$ streamlines. A closed-form expression for the analogue of the above adiabatic invariant has been derived for the simpler problem of a translating drop (Bajer & Moffatt 1990, 1992).

Figure 5. (a) A pair of spiralling streamlines, in a deformed-drop octant, for an ambient planar extension, for $Ca=0.005$, $\lambda =1$: the blue and magenta curves correspond to initial points $(0.7,0,0.05)$ and $(0.5,0,0.3)$, respectively, with black arrows indicating the direction of spiralling; the dark green curve is the inner limiting torus. (b) Same streamlines as curves in the $D$–$E$ plane; insets highlight the two-time-scale structure of the $D$–$E$ plane curves, and the decreasing amplitude of the fast wiggles when $Ca$ reduces from $0.005$ to $0.001$. (c) $\lim _{\tau \rightarrow \infty } Nu/Nu_{0}$ vs $Pe\,Ca(1+\lambda )^{-1}$ for various $Ca$ for an ambient planar extension; $\triangledown$ and $\bigcirc$ correspond to $\lambda =0$ and $1$, respectively. Inset shows the corresponding curves from BEM-cum-Langevin simulations for $Ca = 0.005$ (compared with the $O(Ca)$ theory to illustrate good agreement) and $Ca= 0.04$; the double-headed arrows indicate the lengths of the respective primary plateaus.

Figure 5(c) shows the long-time limiting values of $\lim _{\tau \to \infty } Nu/Nu_0$ as a function of $Pe\,Ca(1+\lambda )^{-1}$ for $Pe \lesssim 10^8$, and for $Ca$ ranging from $5\times 10^{-4}$ to $0.01$. The latter value is about six times less than the critical $Ca$ for breakup (Bentley & Leal 1986a), ensuring the validity of the weak deformation assumption and, thence, of the $O(Ca)$ velocity field. Interestingly, for a deformed drop in an ambient planar extension, $Nu/Nu_0$ is no longer bounded from above by $4.1$ which now only serves as a primary plateau attained when $Pe \gg 1$, $Pe\,Ca \ll 1$. As $Pe\,Ca$ increases beyond $O(100)$, $Nu/Nu_0$ rises above $4.1$, eventually saturating in a higher secondary plateau for $Pe\,Ca\rightarrow \infty$. While the primary plateau is determined by diffusion along the two coordinates orthogonal to the one along the zero-$Ca$ closed streamlines, the secondary plateau is determined by one-dimensional diffusion across invariant tori traced out by finite-$Ca$ spiralling streamlines in figure 5(a). The secondary plateau value is $Ca$-dependent in general, but approaches $22.3$ for $Ca \rightarrow 0$. Thus, remarkably, for any finite $Ca$ however small, one obtains a nearly six-fold enhancement over the primary plateau ($\approx$4.1) for sufficiently large $Pe (\gg Ca^{-1})$.

While the zero-$Ca$ velocity field for planar extension exhibits a simple $(1+\lambda )^{-1}$-scaling, the $O(Ca)$ correction has a more complicated dependence on $\lambda$ (see Appendix A). Nevertheless, the secondary plateaus for $\lambda = 0$ and $1$, in figure 5(c), nearly coincide for small $Ca$, suggesting that the underlying adiabatic invariant parameterising the small-$Ca$ tori is only weakly dependent on $\lambda$. Towards validating our results based on the $O(Ca)$ velocity field, and extending them to larger $Ca$, we have used the boundary integral method (or BEM; see § 2.3 and Appendix B) to compute the finite-$Ca$ velocity field that is then used in the Langevin simulations. The results of the BEM-cum-Langevin simulations are shown as an inset in figure 5(c), and confirm the increase in $Nu/Nu_0$ beyond the primary plateau even for finite $Ca$. The results for the higher $Ca\ (=0.04)$ show the gradual disappearance of the separation between the primary and secondary plateaus with increasing $Ca$, and an increase in the secondary plateau value.

To understand the mechanism underlying the emergence of the secondary plateaus in figure 5(c), we first consider a spherical drop ($Ca = 0$). Figure 6(a) shows the (normalised) angle-averaged interior concentration profiles for pure diffusion and an ambient planar extension, as a function of the radial distance ($r$), for $Ca = 0$. These correspond to the long-time limit when the concentration profile in a given simulation decreases in amplitude, while preserving its shape; as already seen in Figure 4(b) (inset), this long-time regime manifests as a plateau in the $Nu$ vs time plots. In figure 6(a), the concentration decreases monotonically from the drop center for pure diffusion, but is a non-monotonic function of $r$ for planar extensional flow, with the peak concentration occurring along the fixed point circles of radius $\sqrt {{3}/{5}}$ identified in figure 4(a); the concentration peak at $r = \sqrt {{3}/{5}}$ persists despite the angle-averaging. In light of this, one may interpret the long-time transport, on the zero-$Ca$ large-$Pe$ primary plateau, as occurring due to diffusion from the fixed point circles to the drop surface. The angle-averaged diffusion length ($1-\sqrt {{3}/{5}}$) being smaller than the drop radius (unity), leads to the enhanced transport associated with the large-$Pe$ plateau in figure 4(b): $\lim _{\tau \to \infty } Nu/Nu_0 = 4.1$ as opposed to $Nu/Nu_0 = 1$ for pure diffusion.

Figure 6. (a) Angle-averaged concentration profiles corresponding to the slowest decaying eigenmodes, plotted as a function of $r$, for pure diffusion ($Pe = 0$; $c(r) \propto {\sin ({\rm \pi} r)}/{r}$) and for an ambient ambient planar extension; $Pe=10^{6}$, $\lambda =1$ in the latter case. (b) The angle-averaged concentration profiles, as a function of radial distance $r$, for the slowest decaying eigenmodes for $Ca=0$ ($\lambda = 1$, $Pe = 3000$) and $0.001$ ($\lambda = 1$, $Pe = 3000$, $10^8$), for an ambient planar extension.

Next, in figure 6(b), we compare the angle-averaged concentration profiles, again for an ambient planar extension, for $Ca = 0$ and $0.001$. The concentration profiles are plotted for two $Pe$ values chosen such that the smaller $Pe (=3000)$ corresponds to the primary plateau ($Pe \gg 1,\ Pe\, Ca \ll 1$), and the larger $Pe\ (=10^8)$ corresponds to the secondary plateau ($Pe\,Ca \gg 1$). The zero-$Ca$ and finite-$Ca$ profiles are virtually coincident at the smaller $Pe$, implying that the primary plateau is insensitive to drop deformation. In contrast, at $Pe = 10^8$, the finite-$Ca$ profile differs markedly from the zero-$Ca$ one, this difference being a reflection of the qualitatively different finite-$Ca$ streamline topology. The concentration now peaks along the innermost torus, the singular green curve shown in figure 5(a), and this leads to the angle-averaged concentration profile exhibiting two peaks, corresponding to the inner and outer halves of the aforementioned singular curve. The outer peak is much closer to the drop surface than the original zero-$Ca$ peak, resulting in an $O(1)$ reduction in the effective diffusion length for any non-zero $Ca$ and, thus, a substantial transport enhancement associated with the secondary plateau.

The concentration profiles, for $Pe$ pertaining to the primary and secondary plateau regimes, are illustrated in figure 7 via contour plots; figure 7(a) corresponds to $Ca = 0$, $Pe = 5\times 10^6$ and figure 7(b) corresponds to $Ca = 0.001$, $Pe = 10^8$. Since there is no difference between the large-$Pe$ concentration profiles for $Ca = 0$, and that on the primary plateau for non-zero $Ca$, figure 7(a) may be regarded as representative of the primary plateau. Note that both figures are for a fixed non-zero $x_3\ (=0.5)$ because the innermost torus (singular curve) in figure 5(a), that must correspond to the concentration maximum for $Pe$ on the secondary plateau, does not intersect the $x_1$–$x_2$ plane. Expectedly, the concentration maxima in figure 7(a) correlate with the points of intersection of the two fixed point circles with the plane $x_3 = 0.5$ (shown as pink dots). In contrast, figure 7(b) shows the bifurcation of the original zero-$Ca$ peak value into a pair of (slightly smeared out) maxima for small but finite $Ca$, in each of the four quadrants. The bifurcated maxima correlate well with the points of intersection of the singular curve (shown as red dots), and are consistent with the twin-peaked concentration profile in figure 6(b).

Figure 7. Concentration contours for an ambient planar extension, at the cross-section $x_{3}=0.5$, for (a) $Ca=0$, $Pe=5 \times 10^{6}$, $\lambda =1$ and (b) $Ca=0.001$, $Pe=10^{8}$, $\lambda =1$; colour bars on the right-hand side denote concentration magnitudes. Magenta dots denote the intersections of the zero-$Ca$ fixed-point circles with the plane $x_3 = 0.5$, and correlate with the concentration maxima in (a). Red dots, corresponding to the intersections of the innermost torus (singular curve), for $Ca =0.001$, correlate with the maxima in (b); all interior streamlines are depicted as black curves.

4. Scalar transport in other ambient planar linear flows

4.1. Transport from a spherical drop

Figure 8(a) depicts the inplane streamline pattern for $Ca = 0$, while figure 8(b,c) characterise the geometry of the off-plane zero-$Ca$ interior streamlines, for a drop immersed in an ambient hyperbolic planar linear flow (with $\beta = 0.5$). The absence of the four-fold rotational symmetry characteristic of planar extension implies that the interior streamline pattern no longer exhibits an eight-fold symmetry as in figure 4(a). The analogues of the two fixed-point circles in figure 4(a) now have differing radii, and organise the off-plane streamlines into three distinct groups, and this is indicated by the numerical labels in figure 8(a). Each of the invariant regions in figure 8(a) is characterised by a different $\phi$-interval, $\phi$ being the azimuthal angle in the $x_1$–$x_2$ plane. Region $1$ corresponds to $\phi _{1} \leq \phi \leq {\rm \pi}-\phi _{1}$ and ${\rm \pi} +\phi _{1} \leq \phi \leq 2 {\rm \pi}-\phi _{1}$; region $2$ corresponds to $0 \leq \phi \leq 2 {\rm \pi}$; while region $3$ corresponds to $2 {\rm \pi}-\phi _{2} \leq \phi \leq \phi _{2}$ and ${\rm \pi} -\phi _{2} \leq \phi \leq {\rm \pi}+ \phi _{2}$, where $\phi _1$ and $\phi _2$ are functions of $\beta$ and $\lambda$; a fixed-point analysis of the system ${\textrm {d}\boldsymbol {x}}/{\textrm {d}t} = \hat {\boldsymbol {u}}^{(0)}(\boldsymbol {x})$ yields $\phi _1=\tan ^{-1} \sqrt {{(5+\beta +2 \lambda (1- \beta ))}/{(5 \beta +1-2\lambda (1-\beta ))}}$ and $\phi _2= \tan ^{-1}\sqrt {{(2 \beta -\lambda (1-\beta ))}/{(2+\lambda (1-\beta ))}}$. Off-plane streamlines belonging to two of the aforementioned regions (regions $1$ and $3$) are shown in figure 8(b), while those belonging to the third group (region $2$) appear in figure 8(c). Taken together, these three groups of streamlines fill up a drop hemisphere, with the streamlines in the other hemisphere being obtained via a reflection-symmetry.

Figure 8. Closed streamlines within a spherical drop ($Ca = 0$), in an ambient planar linear flow with $\beta =0.5$; $\lambda =1$. (a) The inplane streamline pattern serves to demarcate three invariant regions (regions $1$, $2$ and $3$), separated by pairs of homoclinic (regions $1$ and $2$) and heteroclinic (regions $2$ and $3$) orbits which are shown as solid black curves; the expressions for the threshold angles $\phi _1$ and $\phi _2$ indicated, in terms of $\beta$ and $\lambda$, are given in the text; $\phi _{1}=58.19^{\circ }$ and $\phi _{2}=24.09^{\circ }$ in the figure. (b) and (c) Depictions of the three-dimensional closed streamlines in the three different regions; each of these subfigures includes a pair of black dotted curves corresponding to the fixed-point circles. The streamlines (in red and blue) that encircle the fixed-point circles appear in (b); those (in pink) whose projections onto the $x_1$–$x_2$ plane encircle the drop centre appear in (c). (d) Plot of $\lim _{\tau \rightarrow \infty } Nu/Nu_{0}$ vs $Pe$, for $\lambda = 1$, for a drop immersed in different members of the planar linear flow family, ranging from simple shear flow ($\beta =0$) to planar extension ($\beta =1$), for $Ca =0$; the primary plateau values have been indicated.

As already mentioned, for ambient planar linear flows other than planar extension, the interior streamline pattern exhibits a non-trivial $\lambda$-dependence. Nevertheless, and despite the additional geometrical complexity, almost all interior streamlines remain closed curves both for $\beta = 0.5$, as evident from figure 8(a–c), and for other $\beta$ values (not shown). As a result, scalar transport remains diffusion-limited in the limit $Pe \rightarrow \infty$. Thus, as shown in figure 8(d) for $\beta = 0$, 0.25, 0.5, 0.75 (and $1$), with $\lambda = 1$, the $\lim _{\tau \rightarrow \infty } Nu/Nu_0-Pe$ curves, for $Ca = 0$, exhibit a $Pe$-dependence analogous to that already seen in figure 4(b) for $\beta =1$. All of the $Nu$-curves start from unity for $Pe \rightarrow 0$, and eventually saturate in a $(\lambda,\beta )$-dependent large-$Pe$ diffusion-limited plateau. For all $\lambda$, the plateau value is a monotonically decreasing function of $\beta$. For $\lambda = 1$ in particular, the plateau reduces from $4.1$ for planar extension, to $1.08$ for simple shear flow and, finally, to unity (no convective enhancement) for solid-body rotation (not shown).

4.2. Transport from a deformed drop

Figure 9(a,b) depict successive portions of a single finite-$Ca$ interior streamline, that starts from $(x_1,x_2,x_3) \equiv (0.01,0.05,0.1)$, again for an ambient hyperbolic planar linear flow with $\beta = 0.5$. Similar to the planar extensional case in figure 5(a), the finite-$Ca$ streamline is no longer closed, and has a tightly spiralling character instead. A single turn of the spiral is completed in a time of $O(\dot {\gamma }^{-1})$, while the spiralling streamline drifts slowly across a volume comparable to any of the invariant regions (with projections as shown in figure 8a) in an asymptotically longer time of $O(Ca^{-1}\dot {\gamma }^{-1})$; similar to figure 8(a), the direction of the drift is indicated by black arrows. In addition, however, the finite-$Ca$ streamline is also seen to cross over from one invariant region to another on time scales even longer than $O(Ca^{-1}\dot {\gamma }^{-1})$, eventually accessing the entire (deformed) hemisphere. Such crossings are marked by a change in colour of the spiralling streamline in figure 9(a,b). Thus, in figure 9(a), the streamline starts off (in blue) from the point $(0.01,0.05,0.1)$ in region $1$, crossing over from region $1$ to region $2$ (marked by a change from blue to pink), and then crossing over to region $3$ (change from pink to red) and reaching the point $(-0.85,-0.27,0.04)$. In figure 9(b), the streamline continues on from $(-0.85,-0.27,0.04)$, is seen to return to region $2$, and then crosses over to region 3 (change from pink to red), finally ending up at $(0.88,-0.37,0.08)$.

Figure 9. (a) A finite-$Ca$ spiralling streamline, within a deformed-drop hemisphere, in an ambient planar linear flow with $\beta =0.5$; $Ca=0.01$, $\lambda =1$. The streamline starts from $(x_1,x_2,x_3) \equiv (0.01,0.05,0.1)$, with black arrows indicating the direction of spiralling in different regions of the drop hemisphere; the black dashed curve denotes the deformed drop profile in the $x_1$–$x_2$ plane, with the pair of black dotted curves denoting fixed-point circles. (b) Continuation of the spiralling streamline in (a) for longer times. Note that the colour of the spiralling streamline in both (a) and (b) changes in accordance with the region that it samples; the colour-coding follows figure 8(a–c). (c) The apparently chaotic region crossings, of a single finite-$Ca$ streamline, starting from $(0.01, 0.05, 0.2)$, illustrated via a $\phi$ vs $\tau$ plot for $\beta =0.5$, $\lambda =1$ and $Ca=0.005$; $\tau =Ca\,t$. Numerical labels denote the particular region of the zero-$Ca$ streamline pattern sampled (see figure 8a); accompanying insets highlight the multiple time scale structure of the finite-$Ca$ streamline.

The above region crossings are evidently absent for $Ca = 0$, and were absent for planar extension even for finite $Ca$, being precluded by the greater (four-fold) symmetry of the ambient flow. Thus, the crossings arise, for $\beta \neq 1$ owing to the drop-deformation-induced perturbation breaking open the homoclinic and heteroclinic trajectories present in the originally spherical drop; recall that these trajectories were shown as solid black curves in figure 8(a). The times corresponding to the aforementioned region crossings appear to be random, this randomness being indicative of the emergence of Lagrangian chaos for small but finite $Ca$ and $\beta \neq 1$. The onset of chaotic dynamics may be understood from the fact that broken homoclinic/heteroclinic connections must lead to transversely intersecting stable and unstable manifolds within the deformed drop (Rom-Kedar, Leonard & Wiggins 1990). A short stretch of a numerically generated sequence of region crossings, for $Ca = 0.005$, appears in figure 9(c), where the crossings have been characterised via a plot of $\phi$ as a function of time. This depiction exploits the fact that, independent of $x_3$, each of the three originally invariant regions is characterised by a different $\phi$-interval. The threshold angles $\phi _1$ and $\phi _2$ demarcating the different $\phi$-intervals were defined earlier, and for the particular choice of parameters in figure 9(c), are $\phi _{1}=58.19^{\circ }$ and $\phi _{2}=24.09^{\circ }$.

The altered finite-$Ca$ streamline topology described previously leads to a transport enhancement, for sufficiently large $Pe$, similar to the planar extensional case discussed in § 3.2. This is seen in figure 10 which shows the $\lim _{\tau \rightarrow \infty } Nu/Nu_0-Pe$ curves, for $Ca = 0.001$ and $0.005$ with $\lambda = 1$, possessing both primary and secondary plateaus; the primary plateaus are those associated with a spherical drop, and were already shown in figure 8(d). In fact, for $0 < \beta < 1$, the difference between the primary and secondary plateaus is seen to be greater than that for planar extension: for instance, for simple shear flow ($\beta = 0$), for $Ca=0.005$, the secondary plateau value ($9.2$) is almost nine times the primary one ($1.08$). An intriguing feature of the $Nu$–$Pe$ curve, for $\beta = 0.75$, $Ca = 0.005$, is the continued slow increase of $Nu$ even for the largest $Pe$ values, indicating a much more gradual approach towards a secondary plateau; unlike the planar extension, the angle-averaged concentration profile (not shown) continuously changes along the slowly rising curve. The zero-$Ca$ streamline pattern for $\beta = 0.75$ has homoclinic and heteroclinic connections similar to those shown in figure 8(a), but with the different invariant regions being comparable in size, which in turn leads to a greater frequency of (random) region crossings. The delayed plateauing of the $Nu/Nu_0$ curve, for $\beta = 0.75$, might therefore be related to a tertiary enhancement regime at still larger $Pe$, arising from Lagrangian chaos.

Figure 10. Plot of $\lim _{\tau \rightarrow \infty } Nu/Nu_{0}$ vs $Pe\,Ca$, for $\lambda = 1$, for a drop immersed in different members of the planar linear flow family; circles ($\circ$) and crosses ($\times$) correspond to $Ca=0.001$ and $0.005$, respectively.

5. Conclusions

By considering a drop in a range of ambient shearing flows, we have shown analytically and numerically, for the first time, the transition of closed interior streamlines, to spiralling ones, with the onset of drop deformation. Remarkably, and in contrast to earlier efforts, the altered streamline topology leads to a singular enhancement in transport rates: transport from even an infinitesimally deformed drop is manyfold greater than that from a spherical one. Importantly, the singular role of streamline topology, and the associated transport enhancement, extends beyond drops. Deformation-induced transport enhancement must occur for virtually all other suspended elastic microstructures, including capsules, vesicles and red blood cells in the tank-treading mode. There also exist implications for cytoplasmic streaming in large Eukaryotic cells, where closed streamlines lead to diffusion-limited transport (Goldstein, Tuval & van de Meent 2008).

The interior transport problem examined here is an inherently unsteady one, with the concentration field decaying to zero, for sufficiently long times, for all cases considered. While $Nu$ has been defined in a manner as to asymptote to a time-independent plateau, one may nevertheless wonder if there exists an alternate steady-state definition of scalar transport. A steady state may be achieved by including, for instance, a spatially uniform source within the drop, that, for sufficiently long times, compensates for the absorption at the surface. This is equivalent to a constant forcing on the right-hand side of (2.1). At steady state, the flux at the surface may then be computed in the usual manner. For pure diffusion, the steady-state concentration field may be obtained analytically, being proportional to $(1-r^2)$, and the associated flux may then be used for purposes of normalisation, when defining $Nu$. While an analytical expression is not available for a deformed drop in an ambient linear flow, the steady-state concentration field in this general case may nevertheless be shown to be the time integral of the transient one (one in the absence of the source forcing), and may therefore be obtained via a numerical integration (from $t = 0$ to $\infty$) of the transient field obtained via Langevin simulations (Sabarish 2021). In this way, one may evaluate the steady-state flux at finite $Pe$, with the ratio of this flux to the steady state diffusive one above yielding an alternate non-dimensional measure of the transport enhancement, one that only involves steady-state quantities. The two enhancement measures above are not the same. That based on the transient concentration field, and that has been used to generate the $Nu$ vs $Pe$ plots in this paper (figures 4b, 5c, 8d and 10), involves only the slowest decaying eigenmode, while the steady-state measure involves all the eigenmodes of the convection–diffusion operator. Nevertheless, the nature of the $Nu$–$Pe$ relationship is expected to remain insensitive to this change in definition, owing to the close resemblance of the slowest decaying eigenmode to the steady-state concentration profile, as illustrated in figure 11(a,b), for pure diffusion and for an ambient planar extension.

Figure 11. The concentration profiles corresponding to steady state (due to a spatially uniform source), and that corresponding to the slowest decaying eigenmode, as a function of radial distance $r$: (a) pure diffusion ($c(r)\propto {(1-r^2)}/{6}$ and $c(r)\propto {\sin ({\rm \pi} r)}/{r}$) and (b) an ambient planar extension ($Pe=10^{6}$, $\lambda =1$).

The current literature does not recognise the role of streamline topology in transport processes, implicitly assuming drop deformation to only have a perturbative effect on the scalar transport rate. In order to further emphasise the critical role of streamline topology, we therefore briefly examine a drop in an ambient biaxial axisymmetric extension. Axisymmetry of the ambient flow in this case implies that the finite-$Ca$ streamlines continue to be closed curves in a meridional plane, leading to the large-$Pe$ transport being diffusion-limited even for a deformed drop. Indeed, as shown in figure 12(a,b), the $\lim _{\tau \to \infty } Nu/Nu_0$ vs $Pe$ curves, for biaxial axisymmetric extension, only show a single large-$Pe$ plateau for both zero and finite $Ca$, with the plateau value only exhibiting a perturbative increase from $4.5$ for $Ca = 0$ to about $4.9$ for $Ca =0.01$. The insets in the said figures show that the nature of the concentration contours remain virtually unaltered for non-zero $Ca$. In fact, the lack of appreciation of the role of streamline topology in convective enhancement has led to the few efforts in literature, which analyse transport from deformed drops, to only consider an ambient (uniaxial or biaxial) axisymmetric extension; see, for instance, Favelukis & Lavrenteva (2013, 2014), Favelukis (2015, 2016) and Liu et al. (2018, 2019). As expected from the above discussion, all of these efforts find drop deformation to only have a perturbative effect. Note that, unlike $Ca = 0$, the finite-$Ca$ transport is not invariant to ambient flow reversal in the limit $Pe \rightarrow \infty$ (Brenner 1967). Nevertheless, since the effect of drop deformation remains perturbative for both uniaxial and biaxial extension, the respective large-$Pe$ plateaus are expected to differ only by $O(Ca)$.

Figure 12. (a) Plot of $\lim _{\tau \rightarrow \infty } Nu/Nu_0$ vs $Pe$ for an axisymmetric biaxial extension for $Ca = 0$; the inset shows the closed streamline pattern in a meridional plane, with black arrows denoting the direction of circulation, along with the concentration contours. (b) Plot of $\lim _{\tau \rightarrow \infty } Nu/Nu_{0}$ vs $Pe\,Ca$, for various $Ca$ and $\lambda = 1$, for a deformed drop in an axisymmetric biaxial extension; the inset emphasises the invariance of the closed-streamline pattern, and the concentration contours, to drop deformation. Both insets correspond to $Pe=10^{8}$, $\lambda =1$.

The significance of the disparate $Nu$–$Pe$ relationships for planar (figure 5c) and axisymmetric (figure 12b) extensional flow is best appreciated when one accounts for the virtually identical character of these two flows in all other scenarios. Both are prototypical strong flows that lead to material points separating exponentially in time. Polymer molecules in both flows undergo a coil-stretch transition, leading to an extension-thickening rheology and a dramatic departure from Newtonian behaviour. In addition, the critical $Ca$, for breakup, exhibits an identical $\lambda$-dependence for both flows, transitioning from a $\lambda ^{-{1}/{6}}$-scaling regime for small $\lambda$ to a large-$\lambda$ plateau (Acrivos & Lo 1978; Hinch & Acrivos 1980), pointing to the analogous nature of drop deformation and breakup. Indeed, the original asymptotic analysis for the critical $Ca$ in the limit $\lambda \rightarrow 0$, for planar extension, relied on treating this flow as an ambient axisymmetric extension plus a three-dimensional perturbation (Hinch & Acrivos 1980); the results for axisymmetric extension had been derived in an earlier effort (Acrivos & Lo 1978), and those for planar extension were obtained for finite amplitudes of the said perturbation. Despite the closely analogous behaviour in all of the aforementioned scenarios, it is shown here that the differing nature of the finite-$Ca$ streamlines leads to a profound difference in scalar transport characteristics.