Mass transfer from small spheroids suspended in a turbulent fluid

Abstract

By coupling direct numerical simulation of homogeneous isotropic turbulence with a localised solution of the convection–diffusion equation, we model the rate of transfer of a solute (mass transfer) from the surface of small, neutrally buoyant, axisymmetric, ellipsoidal particles (spheroids) in dilute suspension within a turbulent fluid at large Péclet number, $\textit {Pe}$. We observe that, at $\textit {Pe} = O(10)$, the average transfer rate for prolate spheroids is larger than that of spheres with equivalent surface area, whereas oblate spheroids experience a lower average transfer rate. However, as the Péclet number is increased, oblate spheroids can experience an enhancement in mass transfer relative to spheres near an optimal aspect ratio $\lambda \approx 1/4$. Furthermore, we observe that, for spherical particles, the Sherwood number $\textit {Sh}$ scales approximately as $\textit {Pe}^{0.26}$ over $\textit {Pe} = 1.4\times 10^{1}$ to $1.4\times 10^{4}$, which is below the $\textit {Pe}^{1/3}$ scaling observed for inertial particles but consistent with available experimental data for tracer-like particles. The discrepancy is attributed to the diffusion-limited temporal response of the concentration boundary layer to turbulent strain fluctuations. A simple model, the quasi-steady flux model, captures both of these phenomena and shows good quantitative agreement with our numerical simulations.

1. Introduction

When rigid particles are immersed in a turbulent fluid, they may slip, spin, tumble and reorient themselves relative to their chaotic surroundings (Voth & Soldati 2017). Simultaneously, soluble material may be transferred away from the surface by convection and diffusion. Practical examples include industrial processes such as crystallisation (Myerson 2002) or the dissolution of fine solids (Armenante & Kirwan 1989), the transfer of nutrients to planktonic osmotrophs (Karp-Boss, Boss & Jumars 1996), the leaching of pollutants from microplastics (Law 2017; Seidensticker et al. 2017) and the encounter rates of bacteria with marine viruses (Guasto, Rusconi & Stocker 2012). Such particles are rarely ever spherical. Often, the particle is sufficiently small so that its motion approximates that of a tracer embedded in a linear shear. Also, the material may diffuse slowly, so that convection is the dominant mechanism of mass transfer away from the particle. For this case, where the particle Reynolds number is small, we seek to examine two questions: What is the rate of mass transfer (e.g. the solute flux) from the surface? And how does it depend upon particle shape?

There are multiple approaches to this and related problems in the literature (see e.g. the works reviewed by Clift, Grace & Weber (1978) and Crowe et al. (2012)). All seek to relate the Sherwood number $\textit {Sh}$, a dimensionless measure of the mass transfer rate, to the fluid, flow and particle properties. The relevant independent dimensionless groups are as follows. The turbulence is characterised by the Taylor-microscale Reynolds number $\textit {Re}_\lambda$ (Pope 2000). The particle Stokes number ${{\textit {St}}} \equiv \tau _p/\tau _\eta$ parametrises the response time of the particle $\tau _p$ with respect to the Kolmogorov time scale $\tau _\eta$ of the turbulence, characteristic of the average shear rate experienced by the particle (Toschi & Bodenschatz 2009). When ${{\textit {St}}} \ll 1$, the particle motion is tracer-like. The relative time scale of diffusion at the particle length scale $r$ is described by the Péclet number $\textit {Pe} \equiv r^{2}/\kappa \tau _\eta$, where $\kappa$ is the mass diffusion coefficient of the solute. When $\textit {Pe} \gg 1$, convection dominates the mass transfer away from the particle. Sometimes, these ratios are alternatively described by a particle Reynolds number $\textit {Re} = r^{4/3}\left \langle {\epsilon } \right \rangle ^{1/3}/\nu \propto {{\textit {St}}}^{2/3}$ and Schmidt number $\textit {Sc} = \nu /\kappa$, where $\left \langle {\epsilon } \right \rangle$ is the mean kinetic energy dissipation rate and $\nu$ is the kinematic viscosity of the fluid. Jointly, the assumption of large $\textit {Pe}$ and small ${{\textit {St}}}$ implies that the Schmidt number is large; this is in fact often the case for solutes in liquids, e.g. $\textit {Sc} \approx 10^{3}$ for NaCl in water (Lide 2013).

At an aggregate scale, several experimental curve fits (‘correlations’) are available which describe the mass transfer rate as a function of the bulk properties of the particle, fluid and turbulence (Harriott 1962; Levins & Glastonbury 1972; Sano, Yamaguchi & Adachi 1974; Asai, Konishi & Sasaki 1988; Armenante & Kirwan 1989). However, there has been a historical paucity of data describing the average mass transfer rate for small, tracer-like particles in turbulent flows with ${{\textit {St}}} \ll 1$. In a comprehensive review of the literature available to them, out of thousands of data points published, Armenante & Kirwan (1989) identified fewer than 26 data points corresponding to particles with $\textit {Re} < 0.1$. They therefore conducted a comprehensive set of experiments to determine the mass transfer rate of small, spherical particles with ${{\textit {St}}} \ll 1$ and large Péclet number. They demonstrated that, in the tracer regime, the available correlations for inertial particles considerably underestimated the mass transfer rate of tracers, which was attributed to a difference in the convective transfer mechanism. Contemporaneously, Asai et al. (1988) conducted similar experiments and came to the same conclusion: because the relative translation between a tracer particle and its surrounding is negligible, convective mass transfer is provided by the linear shear in which it is embedded. These data provide a valuable reference at an aggregate scale, but do not provide information at the scale of individual particles, which is necessary for the closure of point-particle direct numerical simulation models (Subramaniam 2013; Deen et al. 2014). Moreover, as recently discussed by Oehmke & Variano (2021), there are very few data concerning the mass transfer from non-spherical particles embedded in turbulent flows. Their experimental data suggest that large (Taylor-microscale-sized), neutrally buoyant, rod-like particles dissolve more rapidly than disc-like particles of equivalent volume in isotropic turbulence. However, a systematic investigation of the effects of particle shape and size upon the average transfer rate remains unexplored.

At the scale of the particle, this problem is less well studied. From a theoretical perspective, the flux to small particles in steady flows is well understood (Leal 2012) and is known to scale as $\textit {Sh} = \alpha \textit {Pe}^{1/3} + O(1)$ when $\textit {Pe} \gg 1$ provided the streamlines surrounding the particle are open. Here, the Péclet number is suitably redefined in terms of some characteristic strain rate $E$ of the steady flow problem. Depending upon the particle geometry and relative flow field, various asymptotic expressions are available for the coefficient $\alpha$ and subsequent corrections (Acrivos & Taylor 1962; Acrivos & Goddard 1965; Sehlin 1969; Gupalo, Polianin & Riazantsev 1976; Poe & Acrivos 1976; Batchelor 1979; Acrivos 1980; Myerson 2002; Dehdashti & Masoud 2020; Lawson 2021). This is supplemented by a wealth of experimental and numerical data at finite Reynolds and Péclet numbers (Clift et al. 1978; Sparrow, Abraham & Tong 2004; Kishore & Gu 2011; Ke et al. 2018; Ma & Zhao 2020).

However, the relative flow experienced by a tracer in a turbulent flow is not steady and the theoretical description of unsteady mass transfer is much less complete. Analytical expressions for the transfer rate are available only in the low-Péclet-number regime (Pozrikidis 1997; Feng & Michaelides 1998; Michaelides 2003). At high Péclet number, Batchelor (1979) and later Lawson (2021) showed that, for time-periodic motion of small particles in a steady linear shear $\boldsymbol {u}(\boldsymbol {x}) = \boldsymbol{\mathsf{E}}^{0} \boldsymbol {x} + \tfrac 12\boldsymbol {\omega }^{0} \times \boldsymbol {x}$, unsteady concentration fluctuations arise due to the periodic motion of the particle under the action of the strain $\boldsymbol{\mathsf{E}}^{0}$ and vorticity $\boldsymbol {\omega }^{0}$. Here, in a frame of reference co-rotating with the particle, the particle surface appears stationary but is subject to an unsteady velocity field $\boldsymbol {v}(\,\boldsymbol {y},t)$, periodic in time $t$, where $\boldsymbol {y}$ is the spatial coordinate in the co-rotating body frame. Remarkably, the concentration fluctuations scale as $\textit {Pe}^{-1/3}$ in magnitude relative to the mean, provided that the period $T$ of the motion is much smaller than $|\boldsymbol{\mathsf{E}}^{0}|^{-1} \textit {Pe}^{1/3}$. As such, the mean flow field in the body frame of the particle $\tilde {\boldsymbol {v}}(\boldsymbol {t}) = ({1}/{T})\int _{0}^{T}\boldsymbol {v} \,\mathrm {d}t$ determines the mean transfer rate.

Classical asymptotic methods for steady flows can then be applied to determine the mean transfer rate (Acrivos & Goddard 1965; Leal 2012; Lawson 2021). A phenomenon known as the partial suppression of convection by rotation then emerges, whereby the preferential alignment and rotation of the particle with respect to its surroundings alters the mean flow field it perceives. For small spherical and spheroidal particles embedded in rotation-dominated linear shear $|\boldsymbol {\omega }^{0}| \gg |\boldsymbol{\mathsf{E}}^{0}|$, the transfer rate is shown to scale with a Péclet number based on the rate of vortex stretching $\textit {Pe}_{\omega } = E_\omega r^{2}/\kappa$, where $E_\omega = \omega _i^{0} E_{ij}^{0}\omega _j^{0}/ |\boldsymbol {\omega }^{0}|^{2}$ (Batchelor 1979; Lawson 2021), since particles spin aligned with the vorticity. More generally, for axisymmetric particles spinning along a fixed direction $\boldsymbol {p}$, the transfer rate scales with the strain rate $E_p = E_{ij}^{0} p_i p_j$ perceived in that direction (Lawson 2021). Therefore, orientation dynamics plays an important role in determining the transfer rate from small particles.

In analogy to time-periodic flows, Batchelor (1980) argued that, for small spheres in isotropic turbulence, the mean relative flow field experienced by the sphere is an axisymmetric strain whose magnitude is proportional to the mean rate of vortex stretching $E_\omega$. This argument relies upon averaging the flow in the body frame over a long time with respect to $\tau _\eta$, which is permissible in the limit of $\textit {Pe} \rightarrow \infty$, since the concentration boundary layer should be insensitive to velocity fluctuations faster than $\tau _\eta \textit {Pe}^{1/3}$. Since the mean vortex stretching scales with $1/\tau _\eta$, this argument predicts that the transfer rate scales as $\textit {Sh} \approx 0.55 \textit {Pe}^{1/3} + O(1)$. However, at large but finite Péclet number, the concentration boundary layer may remain sensitive to velocity fluctuations occurring on dynamically active time scales. The characteristic strain rate of the flow field to which the concentration boundary layer is sensitive may therefore depend upon $\textit {Pe}$, so an additional $\textit {Pe}$ dependence may be introduced before the asymptotic behaviour is reached. This brings the application of the asymptotic scaling $\textit {Sh} \sim \textit {Pe}^{1/3}$ to large but finite $\textit {Pe}$ into some doubt. On the other hand, the experimental data of Armenante & Kirwan (1989) and Asai et al. (1988) suggest that the scaling exponent is close to  $1/3$.

An additional complication arises when the particle is not spherical. Small aspherical particles (e.g. spheroids) will spin, tumble and preferentially align themselves with respect to their turbulent environment differently, depending upon their shape (Zhang et al. 2001; Mortensen et al. 2008; Pumir & Wilkinson 2011; Parsa et al. 2012; Byron et al. 2015; Zhao et al. 2015; Voth & Soldati 2017). These preferential alignments appear due to the coupling between the Stokesian dynamics of small particles and the turbulence under Jeffery's equation (Jeffery 1922). In isotropic turbulence, elongated rod-like particles tend to orient themselves parallel to the vorticity vector and weakly align themselves with the intermediate eigenvector of the rate-of-strain tensor (Pumir & Wilkinson 2011). In contrast, flattened, disc-like particles align their symmetry axis parallel to the most compressive direction and orthogonal to the vorticity (Voth & Soldati 2017). Therefore, the average straining flow perceived by a particle has a dependence upon shape in turbulent flows. On this basis, one expects that aerodynamic effects due to shape may play an important role in determining the transfer rate from small particles in turbulence.

The above discussion shows that convective mass transfer from small particles with ${{\textit {St}}} \ll 1$ and $\textit {Pe} \gg 1$ embedded in turbulence occurs through linear shear, that the relative flow field experienced by the particle depends upon its shape, and that particle-scale models are required to capture the role of particle shape in the mass transfer process. In this paper, we address this problem by numerically simulating the transfer of a passive scalar (the solute) from a dilute, one-way-coupled suspension of small spheroidal particles in turbulence. This model couples the numerical solution of the unsteady scalar transport from a spheroidal particle embedded in a time-varying linear shear, which is obtained from the Lagrangian time history of spheroidal tracer particles in homogeneous, isotropic turbulence. We validate this approach through comparisons of the average mass transfer rate against experimental data and examine the statistics of the transfer rate as a function of particle shape. These observations reveal the shape dependence of the mean transfer rate and attribute deviations from the classic $\textit {Pe}^{1/3}$ scaling to the time-limited diffusive response of the concentration boundary layer to turbulent velocity fluctuations. Finally, we introduce a particle-scale, analytical model for the transfer rate from spheroidal particles, which we validate against the results of our numerical simulations.

The paper is structured as follows. The numerical simulation procedure is described in § 2. We present a validation of this procedure for spherical particles in § 3.1, then examine the role of particle shape in § 3.2 and the particle-scale dependence of the mass transfer rate in § 3.3. We introduce and validate the quasi-steady flux model in § 4. We present conclusions in § 5.

2. Numerical simulation

2.1. Dilute suspension model

To model a dilute suspension of particles, we adopt a point-particle, Lagrangian approach coupled with a numerical solution of the convection–diffusion equation in a reference frame co-moving and co-rotating with the particle. The basic idea is to use the time history of the relative velocity field experienced by the particle as a forcing of the convection–diffusion equation in its immediate vicinity. The relative velocity field is prescribed by the Lagrangian time history of the relative velocity gradient, plus a Stokesian perturbation due to the presence of the particle. We then solve the convection–diffusion problem in the immediate vicinity of the particle. The boundary condition far from the particle effectively treats incoming fluid as uncontaminated, such that the turbulence provides a continuous supply of fresh fluid. In doing so, we retain the advantage of particle-resolved direct numerical simulation (Feng & Michaelides 2009; Deen et al. 2014; Derksen 2014) that the concentration boundary layer immediately surrounding the particle is resolved, at the expense of disregarding the larger-scale turbulent mixing and transport. We expect this approximation to have a negligible effect upon modelling the mass transfer rate of dilute suspensions of non-inertial particles. For example, Harriott (1962) note that there is no appreciable concentration effect upon mass transfer rates for volume fractions of particles up to $15\,\%$, which is well beyond the dilute suspension regime.

The suspension is modelled as follows. For our turbulent forcing, we use the Lagrangian time history of $4096$ tracer particles in isotropic turbulence at $R_\lambda \approx 433$ obtained from the Johns Hopkins University (JHU) Turbulence Database (Li et al. 2008). Lagrangian tracers are uniformly seeded at positions $\boldsymbol {x}_0$ throughout the simulation domain at time $t = 0$ and their trajectories $\boldsymbol {X}(t;\boldsymbol {x}_0)$ are integrated as tracers $\dot {\boldsymbol {X}}(t; \boldsymbol {x}_0) = \boldsymbol {u}(\boldsymbol {X}(t; \boldsymbol {x}_0), t)$ following the flow. We then obtain the time history of the local velocity gradient $\boldsymbol{\mathsf{G}}(\boldsymbol {X},t)$ experienced by each particle from the database.

This Lagrangian time history is then coupled to the Stokesian rotational dynamics of small spheroids using Jeffery's equation (Jeffery 1922). In the laboratory frame, the spheroid of aspect ratio $\lambda =a/c$ is centred at $\boldsymbol {X}$ and the orientation of its principal semi-axes of length $(a,c,c)$ are described by the rotation matrix $\boldsymbol{\mathsf{R}} = [\boldsymbol {p},\boldsymbol {q},\boldsymbol {r}]$. Jeffery showed that the solid-body rotation rate $\boldsymbol {\varOmega }$ of a torque-free spheroid embedded in an arbitrary linear shear is given by

(2.1)\begin{equation} \boldsymbol{\varOmega} = \frac{1}{2} \boldsymbol{\omega}^{0} + \frac{\lambda^{2} - 1}{\lambda^{2} + 1} \boldsymbol{p}\times \boldsymbol{\mathsf{E}}^{0} \boldsymbol{p}, \end{equation}

where $\omega _i^{0} = -\epsilon _{ijk}G_{jk}$ and $\boldsymbol{\mathsf{E}}^{0} =( {\boldsymbol{\mathsf{G}}}^{\textrm{T}}+\boldsymbol{\mathsf{G}})/2$ are the vorticity and rate-of-strain tensor perceived by the particle in the laboratory frame. The body frame therefore rotates as

(2.2)\begin{equation} \dot{\boldsymbol{\mathsf{R}}} = [\boldsymbol{\varOmega}]_\times \boldsymbol{\mathsf{R}} \end{equation}

so that the relative velocity gradient experienced by the particle in the co-rotating frame is (Lawson 2021)

(2.3)\begin{equation} \boldsymbol{\mathsf{A}} = {\boldsymbol{\mathsf{R}}}^{\textrm{T}}(\boldsymbol{\mathsf{G}} - [\boldsymbol{\varOmega}]_\times)\boldsymbol{\mathsf{R}}. \end{equation}

As before, we apply a similar decomposition of the relative velocity gradient tensor $\boldsymbol{\mathsf{A}}$ into a symmetric strain tensor $\boldsymbol{\mathsf{E}} = (\boldsymbol{\mathsf{A}}+ {\boldsymbol{\mathsf{A}}}^{\textrm{T}})/2$ and perceived vorticity vector $\boldsymbol {\omega }$ in the co-rotating frame. The relative velocity $\boldsymbol {v}$ experienced in the vicinity of the particle is therefore

(2.4)\begin{equation} \boldsymbol{v}(\,\boldsymbol{y},t) = \boldsymbol{\mathsf{E}}\boldsymbol{y} + \tfrac{1}{2} \boldsymbol{\omega} \times \boldsymbol{y} + \boldsymbol{v}'(\,\boldsymbol{y}; \boldsymbol{\mathsf{A}}), \end{equation}

where $\boldsymbol {y}$ is the spatial coordinate in the co-rotating body frame and $\boldsymbol {v}'$ is the Stokes perturbation to the velocity field due to the presence of the particle. The complete expression for $\boldsymbol {v}'$ is rather involved, but can be readily computed from equations given by Kim & Karrila (1991). This co-rotating frame and surrounding Stokes flow correspond to the configuration studied by Batchelor (1980) in his theoretical analysis of the transfer rate from spherical particles.

The time history of the relative velocity gradient provides the convective forcing term for the convection–diffusion problem in the vicinity of the particle. In the body frame, the transport equation for the concentration field $\theta (\,{\boldsymbol {y}}^{\star }, {t}^{\star })$ can be written in dimensionless form as (Leal 2012)

(2.5)\begin{equation} \frac{\partial \theta}{\partial {t}^{{\star}}} + {\boldsymbol{v}}^{{\star}}\boldsymbol{\cdot} \boldsymbol{\nabla} \theta = \frac{1}{\textit{Pe}} \nabla^{2} \theta, \end{equation}

where ${\boldsymbol {v}}^{\star } = \boldsymbol {v}\tau _\eta / r$ is the non-dimensionalised relative velocity, ${t}^{\star }=t/\tau _\eta$ and ${\boldsymbol {y}}^{\star } = \boldsymbol {y}/r$ are the dimensionless spatial and temporal coordinates, and $r$ is the linear dimension of the particle. The turbulent Péclet number $\textit {Pe} \equiv r^{2}/\kappa \tau _\eta = r^{2} \left \langle {\epsilon } \right \rangle ^{1/2}/\kappa \nu ^{1/2}$ therefore forms a control parameter for the convection–diffusion problem, where $\kappa$ is the diffusion coefficient of the solute and $\tau _\eta \equiv (\nu /\left \langle {\epsilon } \right \rangle )^{1/2}$ is the Kolmogorov time scale.

Because this frame is co-rotating, $\boldsymbol {v} = 0$ on the surface of the particle $S_p$, satisfying the no-slip and no-penetration conditions. We model a uniform concentration boundary condition on the surface of the particle

(2.6)\begin{equation} \left\{\begin{array}{@{}ll} \theta(\,{\boldsymbol{y}}^{{\star}}, {t}^{{\star}}) = 1, & {\boldsymbol{y}}^{{\star}} \in S_p, \\ \theta(\,{\boldsymbol{y}}^{{\star}}, {t}^{{\star}}) \rightarrow 0, & |{\boldsymbol{y}}^{{\star}}| \rightarrow \infty, \end{array} \right. \end{equation}

which corresponds to non-dimensionalising the concentration field as $C(\,\boldsymbol {y},t)$ as $\theta = (C-C_0)/(C_1-C_0)$, where $C_1$ and $C_0$ are the concentration of solute at the surface and infinity in physical units. Far from the particle, the concentration vanishes; this models the turbulent flow as providing a continuous supply of fresh fluid to the surface of the particle.

The non-dimensional measure of the mass transfer rate (solute flux) $Q$ is the local Sherwood number

(2.7)\begin{equation} \textit{Sh} = \frac{Q}{4{\rm \pi} r\kappa(C_1 - C_0)} ={-}\frac{1}{4{\rm \pi}} \iint_{S_p} \boldsymbol{\nabla} \theta \boldsymbol{\cdot} \mathrm{d}\boldsymbol{S}. \end{equation}

Conventions differ on the choice of characteristic length scale $r$; here, we adopt the definition $r \equiv \sqrt {A/4{\rm \pi} }$, where $A$ is the surface area of the particle, so that $\textit {Sh} = 1$ for purely diffusive flux of material from a spherical particle.

2.2. Numerical solution

The numerical procedure to implement the model is as follows. To obtain the trajectories of tracer particles from the JHU database, we seeded Lagrangian tracers uniformly across the simulation domain at positions $\boldsymbol {x}_0$ at $t = 0$ and their trajectories $\boldsymbol {X}(t;\boldsymbol {x}_0)$ were integrated using a second-order Runge–Kutta method over the duration of the simulation $T_{sim} = 10.056$, approximately $5.05$ integral time scales. The integration time step is $\Delta t = 0.002 \approx 0.047 \tau _\eta$, which corresponds to the available time steps stored in the database. The velocity and velocity gradient field are interpolated from the database at the tracer position using a fourth-order Lagrange polynomial interpolation. This yields a time series of the velocity gradient $\boldsymbol{\mathsf{G}}(\boldsymbol {X}, t_i)$ at $i = 0,\ldots,5028$ times $t_i = {i}\Delta t$. The orientation of the body frame $\boldsymbol{\mathsf{R}}(t)$ (2.2) is integrated using an adaptive time-step, third-order Runge–Kutta scheme using a piecewise linear interpolation to interpolate $\boldsymbol{\mathsf{G}}(\boldsymbol {X},t)$. The initial condition $\boldsymbol{\mathsf{R}}_0 = \boldsymbol{\mathsf{R}}(0)$ is an isotropically distributed random orientation.

For each particle time history, we solve (2.5) numerically using a second-order finite volume method using a modified version of the scalarTransportFoam solver of OpenFOAM. The mesh and solver are the same as those used in Lawson (2021) and we shall review only the salient details here. Equation (2.5) is discretised on a structured grid in prolate, spherical or oblate spheroidal coordinates $(\mu,\psi,\phi )$, depending upon the particle aspect ratio. In the spheroidal system, the $\mu$ direction is analogous to a radial coordinate, the $\psi$ direction parametrises the polar angle from the symmetry axis, and $\phi$ is an azimuthal angle. The inner boundary at $\mu = \mu _0$ corresponds to the surface of the spheroid, whereas the outer boundary at $\mu = \mu _N$ corresponds to fluid far from the particle. The dimensions of the spheroid are chosen such that the surface area is $4{\rm \pi}$, equivalent in surface area to that of a sphere with unit radius. The outer boundary is chosen such that its largest dimension is $100$ and is very slightly oblate or prolate, having an aspect ratio between $0.999$ and $1.001$.

The mesh is discretised into $150\times 64\times 64$ cells in the $(\mu,\psi,\phi )$ directions, respectively, with uniform spacing in the $\psi$ and $\phi$ directions. To adequately resolve the thin concentration boundary layer, which is of thickness $\delta$, we employ a mesh refinement in the $\mu$ direction such that the grid spacing is $\Delta \mu _{i+1} = R \Delta \mu _i$, where $\Delta \mu _{i+1} = \mu _{i+1}-\mu _i$ is the spacing between adjacent cells in the $\mu$ direction. Owing to the nature of the spheroidal coordinate system chosen, the resolution varies across the surface of the particle. The initial spacing $\Delta \mu _1$ is chosen such that the thickness of the largest cell near the particle surface is at most $2\times 10^{-4} \min (a,c)$ and the mesh refinement factor $R$ is chosen accordingly. Based on an estimated boundary layer thickness $\delta = \textit {Pe}^{-1/3}$ at the largest Péclet number tested, for the most extreme aspect ratios $\lambda = 16$ ($\lambda = 1/16$), this yields between $26~(45)$ and $68~(82)$ cells within a distance $\delta \approx 0.0414$ from the surface.

We impose the Dirichlet boundary condition $\theta = 1$ (2.6) on $S_p$ and the von Neumann boundary condition on the outer boundary of the simulated domain to approximate the zero-concentration boundary far from the particle. Time stepping is performed using an implicit Euler scheme and a time step of $\Delta t = 0.001$ from the initial condition ${\theta (\,{\boldsymbol {y}}^{\star }, 0) = 0}$. Statistics of the transfer rate are gathered in the interval $75\tau _\eta \leqslant t \leqslant T_{sim}$, which is sufficiently long to permit the statistics of the transfer rate to reach a steady state.

3. Results

3.1. Validation

We first examine the approach of the simulation to a statistically stationary state. In figure 1(a), we show the time dependence of the ensemble-average Sherwood number for spherical particles. After an initial transient, which corresponds to the diffusive growth of the concentration boundary layer, the ensemble-average Sherwood number reaches a steady state after ${\sim }75\tau _\eta$. Averaging over the interval $75\tau _\eta \leqslant t \leqslant T_{sim}$, the deviation of the ensemble-average Sherwood number $\left \langle {\textit {Sh}} \right \rangle$ from the time average of the ensemble $\left \langle {\overline {\textit {Sh}}} \right \rangle$ does not exceed $2.3\,\%$. The variation is comparable for other aspect ratios. By applying a bootstrap resampling to our data, we estimate the statistical uncertainty in $\left \langle {\overline {\textit {Sh}}} \right \rangle$ to be below $\pm 0.20\,\%$ (based on a 95 % confidence interval).

Figure 1. Sherwood number for small spheres in a turbulent flow. (a) Time dependence of the ensemble-average Sherwood number over the duration of the simulation, reaching a statistically stationary state around ${\sim }75\tau _\eta$. Solid lines show ensemble average at varying $\textit {Pe}$; dashed lines show average over the statistically stationary period. (b)  Comparison of the $\textit {Pe}$ scaling of the average Sherwood number obtained from the model against literature data of spheres with ${{\textit {St}}} < 0.1$. Open markers show experimental data of Armenante & Kirwan (1989) ($\triangle$) and Asai et al. (1988) ($\times$), solid black marker shows numerical data, dashed line shows Batchelor's asymptotic scaling ($\textit {Sh} \approx 0.55 \textit {Pe}^{1/3} + O(1$)) and dotted line shows the best power-law fit.

As a validation of our numerical model, we compare the average Sherwood number predicted by the model against experimental data of small spheres with ${{\textit {St}}} < 0.1$ at high Péclet number in turbulent flow. In figure 1(b), we compare the results of our numerical simulations for spherical particles against data digitally extracted from figure 5 of Armenante & Kirwan (1989). These data correspond to experiments performed in a baffled, stirred tank using glycerol–water solutions of varying kinematic viscosity (and therefore Schmidt number) as the working fluid. Additionally, we present similar data extracted from figure 3 of Asai et al. (1988). Markers of similar colour correspond to fixed Schmidt number $\textit {Sc}$, and a comparison across datasets at fixed $\textit {Pe} = (r/\eta )^{2} \textit {Sc}$, where $\eta$ is the Kolmogorov length scale, allows particle size effects to be inferred. To within the scatter of the data from Armenante & Kirwan (1989), the agreement between the experimental and numerical data is good, which serves to validate the modelling approach. The agreement with the data of Asai et al. (1988) is less good; we note that these data suggest systematically smaller transfer rates than in Armenante & Kirwan (1989) and in some cases indicate (unphysically) that $\left \langle {\textit {Sh}} \right \rangle < 1$. The dotted line in figure 1(b) shows the best power-law fit $\left \langle {\textit {Sh}} \right \rangle = 0.99 \textit {Pe}^{0.26}$ to our numerical data for spheres. This fit is observed to be in reasonable agreement with the experimental data over the range shown (of similar success to Batchelor's $\textit {Pe}^{1/3}$ result) and shows a weaker dependence upon $\textit {Pe}$ than expected. We shall examine the physical origin of this anomalous scaling in § 3.3.

3.2. Shape dependence

We now examine the dependence of the average mass transfer rate upon particle shape, as parametrised by the spheroid aspect ratio. Figure 2(a) shows the scaling of the ensemble-average Sherwood number with the turbulent Péclet number for prolate and oblate spheroids with aspect ratios $\lambda = 2^{n}$, $n \in \{ -4,-2,0,2,4 \}$. We observe that, with the possible exception of very strongly flattened spheroids ($\lambda = 1/16$), the slope of the power-law scaling $\left \langle {\textit {Sh}} \right \rangle \sim \textit {Pe}^{\gamma }$ at large $\textit {Pe}$ appears to be less than $1/3$, consistent with the anomalous scaling behaviour seen for spheres.

Figure 2. Mass transfer coefficient for small spheroids in a turbulent flow. (a) The scaling of ensemble-average Sherwood number with turbulent Péclet number for five different aspect ratios $\lambda$, compared to the expected $\left \langle {\textit {Sh}} \right \rangle \sim \textit {Pe}^{1/3}$ scaling for steady flows. (b) The dependence of the ensemble-average mass transfer rate $\left \langle {\overline {\textit {Sh}}} \right \rangle$ upon the aspect ratio $\lambda$, compensating for the expected $\textit {Pe}^{1/3}$ scaling.

To better examine the shape dependence, figure 2(b) shows the variation of the average mass transfer rate as a function of particle aspect ratio. We have compensated the abscissa with the expected $\textit {Pe}^{1/3}$ scaling for steady and time-periodic flows. If this scaling holds, one expects that the curves should approach a constant asymptotic form as $\textit {Pe}$ is increased. However, there is a continuous evolution in the shape dependence with increasing $\textit {Pe}$. We observe that, at fixed Péclet number, there is always a tendency for prolate spheroids to experience a greater average mass transfer rate than spherical particles of equivalent surface area. In contrast, oblate particles do not consistently show the same trend across all $\textit {Pe}$. At moderate Péclet number, oblate particles tend to experience lower mass transfer rate compared to an equivalent sphere. This trend is found to be consistent with the diffusion-limited behaviour at $\textit {Pe} = 0$ (Clift et al. (1978), not shown). However, at $\textit {Pe} = O(10^{3})$, this trend is reversed for spheroids near an optimal aspect ratio $\lambda \approx 1/4$. Furthermore, the relative mass transfer advantage of adopting a strongly prolate shape over a spherical one is increased. For example, a spheroid with $\lambda =16$ experiences roughly ${\sim }31\,\%$ greater flux than an equivalent sphere at $\textit {Pe}=1.4\times 10^{1}$, but $42\,\%$ greater flux at $\textit {Pe} = 1.4\times 10^{4}$. The relative advantage is smaller for less extreme differences in aspect ratio and $\textit {Pe}$. We note that the same qualitative behaviour is observed for spheroids in steady, rotation-dominated ($|\boldsymbol {\omega }^{0}| \gg |\boldsymbol{\mathsf{E}}^{0}|$) flow, where oblate particles also exhibit maximal mass transfer near $\lambda \approx 0.31$ under vortex compression at large $\textit {Pe}$ (Lawson 2021).

3.3. Local dependence of the mass transfer rate

Thus far, we have only considered statistics of the mean mass transfer rate. We now examine the fluctuations in the mass transfer rate. Figure 3(a) shows the standard deviation of mass transfer fluctuations $\textit {Sh}' = \textit {Sh} - \left \langle {\textit {Sh}} \right \rangle$ relative to the mean transfer rate for different aspect ratios and Péclet numbers. Across the range of $\textit {Pe}$ tested, the fluctuation magnitude is between 12 % and 19 % of the mean mass transfer rate, with a peak around $\textit {Pe} = O(10^{3})$. Thus, turbulent fluctuations in the mass transfer rate are relatively small in comparison to the mean and decay at large  $\textit {Pe}$.

Figure 3. Second-order statistics of transfer rate fluctuations $\textit {Sh}'$. (a) The standard deviation of transfer rate fluctuations with respect to the mean, as a function of Péclet number. (b) The power spectral density of the transfer rate fluctuations for spherical particles at varying  $\textit {Pe}$.

Figure 3(b) shows the power spectral density $E(f)$ of the transfer rate fluctuations for a spherical particle, where $f$ is the frequency. At low $\textit {Pe}$, the spectrum rolls off around $f\tau _\eta \approx 2$, corresponding to the most rapid velocity fluctuations in the turbulent forcing. As the Péclet number is increased, the spectral content at high frequencies is attenuated, indicating that the concentration boundary layer becomes unresponsive to rapid velocity fluctuations beyond a certain time scale. This phenomenon can be demonstrated analytically in time-periodic flows (Batchelor 1979; Lawson 2021), where velocity fluctuations occurring on a dimensionless time scale $\mathcal {T} \ll \textit {Pe}^{1/3}$ result in a negligible contribution to the transfer rate at large  $\textit {Pe}$.

To examine the temporal response of the concentration boundary layer in more detail, we consider the cross-correlation $R_{SP}(\tau )$ between the instantaneous transfer rate $\textit {Sh}(t+\tau )$ at time $t+\tau$ and the local Péclet number $\textit {Pe}_{E}(t) \equiv r^{2} |\boldsymbol{\mathsf{E}}| / \kappa$ at time $t$. Here, $|\boldsymbol{\mathsf{E}}|(t) = (E_{ij}E_{ij})^{1/2}$ is the magnitude of the local strain rate experienced by the particle at time $t$ and is an invariant across reference frames ($|\boldsymbol{\mathsf{E}}^{0}| = |\boldsymbol{\mathsf{E}}|$). Figure 4(a) shows the behaviour of this cross-correlation as a function of the time lag $\tau$ for spheres at varying turbulent $\textit {Pe}$. We observe that, as the average Péclet number of the ensemble increases, the local transfer rate correlates less strongly with the local strain rate and over a longer time scale. Furthermore, there is a time delay between fluctuations in the local strain rate and the local transfer rate, evidenced by the location of the correlation peak in figure 4. In other words, the concentration boundary layer takes time to respond to changes in the local strain rate and the time scale of the response increases with $\textit {Pe}$. Figure 4(b) demonstrates that this behaviour is consistent across different spheroid aspect ratios, although there are some quantitative differences in the response time scale.

Figure 4. Correlation between the local transfer rate $\textit {Sh}(t+\tau )$ and local Péclet number $\textit {Pe}_{E}(t)$ as a function of time lag $\tau$ for (a) spheres at varying turbulent Péclet number and (b) spheroids at $\textit {Pe} = 1.4\times 10^{4}$.

To identify this response time scale, we consider the correlation between the transfer rate and the local strain rate filtered at time scale $\tau _f$. We therefore introduce the filtered, relative velocity gradient

(3.1)\begin{equation} \widetilde{\boldsymbol{\mathsf{A}}}(t) = \frac{1}{\tau_f} \int_{-\tau_f}^{0} \boldsymbol{\mathsf{A}}(t+\tau) \,\mathrm{d}\tau. \end{equation}

Thus, $\widetilde {\boldsymbol{\mathsf{A}}}$ is the average velocity gradient recently experienced by the particle over the time scale $\tau _f$. Figure 5(a) shows the peak correlation coefficient between the mass transfer rate of spherical particles and the local Péclet number $\textit {Pe}_{\widetilde {E}} \equiv r^{2} \widetilde {E}/\kappa$, which is based on the magnitude of the filtered strain rate $\widetilde {E}^{2} = \widetilde {{\mathsf{E}}}_{ij}\widetilde {{\mathsf{E}}}_{ij}$. The data for spheroids with other aspect ratios are qualitatively similar. We note that, under this filtering operation, $|\widetilde {\boldsymbol{\mathsf{E}}^{0}}| \ne |\widetilde {\boldsymbol{\mathsf{E}}}|$, emphasising that the local Péclet number represents the average rate of strain recently perceived by the particle. We observe that, for any given $\textit {Pe}$, there is an optimal filter time scale $\tau _{d}$ where the transfer rate correlates most strongly with the filtered velocity gradient. Furthermore, this correlation is strong: for spheres, fluctuations in $\widetilde {E}$ are sufficient to explain roughly $80\,\%$ of the variance in $\textit {Sh}'$. This shows that the mass transfer rate is strongly deterministic and that temporally local information (the average velocity gradient recently experienced by the particle) is sufficient to predict the transfer rate. We therefore identify $\tau _d$ as the response time scale of the concentration boundary layer to velocity gradient fluctuations.

Figure 5. Identification of the concentration boundary layer response time scale. (a) Correlation coefficient between local Péclet number $\textit {Pe}_{\widetilde {E}}$ and mass transfer rate $\textit {Sh}$ for spherical particles. (b) Scaling of the response time scale $\tau _{d}$ with Péclet number for $\lambda = 1/16, \ldots, 16$, showing the best power-law fit for spheres $\tau _d/\tau _\eta = \beta \textit {Pe}^{\gamma }$ with $\beta = 0.84$ and $\gamma = 0.35$, and the best fit $\beta = 0.96$ with fixed $\gamma = 1/3$.

Figure 5(b) shows the dependence of the response time scale $\tau _d$ upon $\textit {Pe}$ for all aspect ratios tested. In all cases the response time scale is well approximated by a power-law fit $\tau _d/\tau _\eta = \beta \textit {Pe}^{\gamma }$, which is illustrated for spheres. In analogy to time-periodic flows, we might expect $\tau _d \sim \textit {Pe}^{1/3}$. The best-fit exponents for $\lambda = 1/16, \ldots, 16$ are $\gamma = 0.32, 0.31, 0.35, 0.32$ and $0.24$. The discrepancy between the best-fit exponent and $1/3$ is small and fixing $\gamma = 1/3$ also provides an adequate fit to the data, with $\beta =0.96$ for spheres. This suggests that, despite the fact that the relative velocity field is not time-periodic, the concentration boundary layer remains insensitive to velocity fluctuations occurring on time scales faster than $O(\textit {Pe}^{1/3}\tau _\eta )$. Alternatively, this can be interpreted as an increase in sensitivity of the concentration boundary layer to lower-frequency motions at larger $\textit {Pe}$. There is some additional dependence of the response time scale upon the aspect ratio: mass transfer from prolate spheroids appears to be sensitive to strain fluctuations on shorter time scales than from oblate spheroids at large $\textit {Pe}$. However, this dependence is weaker than the $\textit {Pe}$ dependence.

To examine this local dependence more closely, we present the conditional average of the local transfer rate to spherical particles given the local Péclet number $\textit{Pe}_{\widetilde{E}} = \textit {Pe} \widetilde {E}\tau _\eta$ in figure 6(a), where the local filter time scale is chosen at the correlation optimum $\tau _f = \tau _d$. We observe that, when the local Péclet number is large relative to the mean, the conditional transfer rate $\left \langle {\textit {Sh}} | {\textit {Pe}_{\widetilde {E}}} \right \rangle$ collapses onto a single curve independently of $\textit {Pe}$. We note that the background flow field experienced by a torque-free spherical tracer is always a pure straining flow. The conditional transfer rate coincides with that obtained for a sphere in steady axisymmetric strain $\widetilde {{\mathsf{A}}}_{11}/2 = -\widetilde {{\mathsf{A}}}_{22} = -\widetilde {{\mathsf{A}}}_{33} = \widetilde {E}/\sqrt {6}$ and approaches the asymptotic scaling $\textit {Sh} \approx 0.905 \textit {Pe}_{\widetilde {E}}^{1/3}$ (Batchelor 1979) when the local Péclet number is made very large. (For spheres embedded in other strain topologies (e.g. planar strain) at equivalent $\textit {Pe}_{\widetilde {E}}$, the transfer rate differs by less than $1\,\%$ in comparison to axisymmetric strain (Batchelor 1979; Lawson 2021).) A similar $\textit {Pe}_{\widetilde {E}}^{1/3}$ scaling of the conditional transfer rate can also be observed for prolate and oblate spheroids (not shown).

Figure 6. Conditional statistics of the local mass transfer rate for spherical tracers. (a) Coloured curves show the conditional average of the mass transfer rate, given the local Péclet number $\textit {Pe}_{\widetilde {E}}$ with filter time scale $\tau _f = \tau _d$. The black dashed line shows the transfer rate for a sphere in steady axisymmetric strain, whilst the dotted line shows the asymptotic scaling. (b)  Probability distribution of the local shear rate $\widetilde {E}^{\star } = \textit {Pe}_{\widetilde {E}} / \textit {Pe}$ at different $\textit {Pe}$ (coloured curves). The black dashed parabola shows a log-normal distribution fit to the data at $\textit {Pe} = 1.4 \times 10^{1}$ for comparison.

The distribution of the local shear rate $\widetilde {E}^{\star } = \widetilde {E}\tau _\eta$ is shown in figure 6(b) for spheres. As one might expect, the shear rate is close to log-normally distributed but skewed towards smaller values of ${\widetilde {E}}^{\star }$, which is consistent with unfiltered measurements of the local kinetic energy dissipation rate $\epsilon = 2\nu E^{2}$ (Yeung et al. 2006). As $\textit {Pe}$ and therefore $\tau _d$ increases, the effective shear rate ${\widetilde {E}}^{\star }$ to which the concentration boundary layer is found to be responsive is reduced in magnitude. Jointly, these observations explain the anomalous mass transfer rate scaling in the range of $\textit {Pe}$ observed. At large $\textit {Pe}_{\widetilde {E}}$, the conditional transfer rate scales as $\textit {Pe}_{\widetilde {E}}^{1/3}$. However, the convection suppression effect attenuates the magnitude of the turbulent shear fluctuations to which the particle is responsive as $\textit {Pe}$ increases. Therefore, the average transfer rate grows less strongly than $\textit {Sh} \sim \textit {Pe}^{1/3}$, because the effective shear rate the boundary layer perceives decreases at increasing $\textit {Pe}$. This can be seen more formally through the statistical identity

(3.2) \begin{equation} \left\langle {\textit{Sh}} \right\rangle = \int\left\langle {\textit{Sh}} | {\textit{Pe}_{\widetilde {E}}} \right\rangle P({\widetilde {E}}^{{\star}})\,\mathrm{d}\widetilde{{E}^{{\star}}} \sim \textit{Pe}^{1/3} \left\langle { ({\widetilde {E}}^{{\star}})^{1/3} } \right\rangle, \end{equation}

where we have substituted the scaling

(3.3)\begin{equation} \left\langle {\textit{Sh}} | {\textit{Pe}_{\widetilde {E}}} \right\rangle \sim \textit{Pe}_{\widetilde {E}}^{1/3}. \end{equation}

For spheres, we find $\langle {({\widetilde {E}}^{\star })^{1/3}}\rangle \sim \textit {Pe}^{-0.053}$ over the range of $\textit {Pe}$ tested, which is in good agreement with the observed departure from the expected $1/3$ scaling law seen in § 3.1. The remaining discrepancy in the exponent ($0.28$ compared to $0.26$) may be attributable to finite Péclet-number effects. We caution that to extrapolate this result to the asymptotic limit requires verification of the local scaling (3.3) and a scaling law for $\langle {({\widetilde {E}}^{\star })^{1/3}}\rangle$ in the limit $\textit {Pe} \rightarrow \infty$, which is beyond the scope of this study. The result of Batchelor (1980) is compatible with (3.2) provided $\langle {({\widetilde {E}}^{\star })^{1/3}}\rangle = O(1)$ at $\textit {Pe} \rightarrow \infty$, as Batchelor's analysis of the average relative flow field predicts.

4. Quasi-steady flux model

In this section, we introduce a simple model for the average mass transfer rate from a particle in a turbulent flow. It is based upon the observation that the local transfer rate correlates most strongly with the filtered relative velocity gradient $\widetilde {\boldsymbol{\mathsf{A}}}$, rather than the instantaneous relative velocity gradient $\boldsymbol{\mathsf{A}}$. We therefore hypothesise that the local solution to (2.5) is approximated by the steady-state solution to

(4.1)\begin{equation} \widetilde{{\boldsymbol{v}}}^{\star}\, \!\boldsymbol{\cdot} \boldsymbol{\nabla}\theta = \frac{1}{\textit{Pe}} \nabla^{2} \theta, \end{equation}

where

(4.2)\begin{equation} \widetilde{{\boldsymbol{v}}}^{\star}\, \!(\,\boldsymbol{y}^{{\star}}) = \frac{\tau_\eta}{r} [\widetilde{\boldsymbol{\mathsf{A}}} \boldsymbol{y} + \boldsymbol{v}'(\,\boldsymbol{y}; \widetilde{\boldsymbol{\mathsf{A}}})] \end{equation}

is the (dimensionless) recent time-average relative velocity field experienced by the particle. This is the quasi-steady assumption: it expresses the concept that the concentration boundary layer is in local equilibrium with velocity fluctuations occurring on time scales longer than a diffusive response time $\tau _d$. The filter time scale $\tau _f$ should therefore be matched to the diffusive time scale $\tau _d$. The orientation dynamics of the particle is accounted for implicitly, since the transfer rate is expressed in terms of the recent time history of the relative velocity gradient $\widetilde {\boldsymbol{\mathsf{A}}}$.

We therefore require a solution for the mass transfer rate from a spheroid embedded in an arbitrary linear shear at large Péclet number. A solution to this problem was recently provided by the author (Lawson 2021) using well-known asymptotic methods (Leal 2012). At large Péclet number, the transfer rate may be written as

(4.3)\begin{equation} \textit{Sh} = \alpha(\widetilde{\boldsymbol{\mathsf{A}}}/\widetilde {E}; \lambda) \textit{Pe}_{\widetilde {E}}^{1/3} + O(1) , \end{equation}

where

(4.4)\begin{equation} \alpha = \frac{0.808}{4{\rm \pi}} \int \left( \int_{\eta_0}^{\eta_1} \rho^{3/2} F^{1/2} \,\mathrm{d}\eta \right)^{2/3}\,\mathrm{d}\zeta. \end{equation}

The expression for the prefactor $\alpha$ involves an integration of the dimensionless surface shear rate $F$ over the particle's surface using a general curvilinear coordinate system $(\xi,\eta,\zeta )$ whose metric tensor $g_{ij}$ defines $\rho = \sqrt {\det {g}}$. For an arbitrary particle geometry and surface shear distribution, (4.4) can be evaluated numerically using the procedure detailed in Lawson (2021) and briefly outlined here in Appendix A.

At first glance, the parameter space of $\alpha$ is large; $\widetilde {\boldsymbol{\mathsf{A}}}$ has nine degrees of freedom in addition to the aspect ratio $\lambda$. Firstly, we note that the rescaling by the local shear rate $\widetilde {E}$ in (4.3) eliminates the dependence of $\alpha$ upon the magnitude of $\widetilde {\boldsymbol{\mathsf{A}}}$. Continuity, axisymmetry and the torque-free condition (2.1) further constrain five degrees of freedom such that the prefactor $\alpha = \alpha (s, \psi, \phi ; \lambda )$ can be written in terms of four parameters which describe the geometry of the flow. These are: the strain topology $s$, a pair of angles $\psi$ and $\phi$ describing the orientation of the spheroid relative to the strain field, and the aspect ratio $\lambda$. We provide a complete definition of these geometric parameters in Appendix B. This parameter space is sufficiently small to allow $\alpha$ to be tabulated for all local flow topologies.

In contrast, in the limit of purely diffusive flux, the transfer rate is prescribed by (Clift et al. 1978)

(4.5)\begin{equation} \textit{Sh}_0 = \frac{c}{r} \frac{(\lambda^{2} - 1)^{1/2}}{\cosh^{{-}1}(\lambda)}. \end{equation}

Interpolating between these two limits, we approximate

(4.6)\begin{equation} \textit{Sh}(t) = \frac{\textit{Sh}_0}{2} + \left( \left(\frac{\textit{Sh}_0}{2}\right)^{3} + \alpha^{3} \textit{Pe}_{\widetilde {E}} \right)^{1/3}, \end{equation}

which provides a close approximation (within 10 %) of the transfer rate for prolate and oblate spheroids in uniform, creeping flow with $\textit {Pe} < 70$ (Clift et al. 1978).

The model is therefore specified by (4.6), up to the dependence of the filter time scale $\tau _f$ upon $\textit {Pe}$. On the basis of the arguments presented in the previous section, we choose

(4.7)\begin{equation} \tau_f = \tau_d \approx \tau_\eta \beta\,\textit{Pe}^{1/3}, \end{equation}

where $\beta =0.96$ is the best-fit coefficient to the boundary layer response time scale in § 3.3. In principle, we might choose to make $\beta$ a function of $\lambda$ to capture the aspect-ratio dependence seen in figure 5(b). However, we find that the model prediction is relatively insensitive to the choice of $\beta$ and adopt this formulation in favour of simplicity.

Figure 7 shows the results of evaluating the quasi-steady flux model. In figure 7(a), we see that the model reproduces the characteristic trends identified in § 3.2 well. The quantitative agreement between the two is excellent. This is remarkable since no free parameters were used to fit the data. The largest error (${\sim }14\,\%$) is observed at the lowest Péclet numbers, where higher-order terms are required to accurately approximate the transfer rate in (4.3) and the error may be introduced by the interpolation (4.6). By construction, the model reproduces the anomalous scaling of the mass transfer rate. Furthermore, the model reproduces the appearance of a local maximum in the transfer rate near $\lambda \approx 1/4$ at large Péclet number, as well as capturing the relative increase in the transfer rate experienced by prolate spheroids in comparison to oblate spheroids. The reversal in the trend for oblate spheroids can therefore be seen as a distinction between two limiting behaviours.

Figure 7. Mass transfer rate predicted by the quasi-steady flux model for small spheroids in a turbulent flow. (a)  Dependence of the transfer rate upon particle aspect ratio at varying turbulent Péclet number. Open markers show numerical simulation results from § 3.2, and lines show the prediction of the quasi-steady flux model (4.6). (b)  Comparison of the relative magnitude of $\textit {Sh}$ fluctuations predicted by the model and numerical simulation.

Figure 7(b) shows the relative magnitude of transfer rate fluctuations $\sigma _{\textit {Sh}} = \langle {\textit {Sh}'}^{2} \rangle ^{1/2}$ predicted by the model (4.6) and that obtained from numerical simulation. Across all aspect ratios and Péclet numbers tested, there is a good agreement (within 40 %) of the magnitude of predicted fluctuations, although there is a tendency to overestimate the magnitude of the fluctuation. This agreement is further supported by measurements of the correlation coefficient $R_{SM}$ between the instantaneous transfer rate predicted by the model and simulation, shown in figure 8. We observe that there is a strong correlation between the two signals, which becomes marginally weaker at more extreme aspect ratio and Péclet number. This demonstrates that the quasi-steady approximation (4.1) holds remarkably well and that the model can predict also the instantaneous mass transfer rate to a good approximation, although the magnitude of $\textit {Sh}'$ fluctuations is slightly overestimated.

Figure 8. Correlation coefficient between transfer rate obtained from quasi-steady flux model and numerical simulation.

5. Conclusions

In this paper, we have presented one-way coupled simulations of the mass transfer rate (solute flux) from neutrally buoyant, spheroidal tracer particles in isotropic turbulence. This approximates the mass transfer process of small particles in dilute turbulent suspension. The simulation is based upon the local solution of the convection–diffusion equation on a conformal grid around a particle, forced by the Lagrangian time history of the flow field perceived by a spheroidal tracer. The numerical simulations are shown to be in good agreement with experimental data from small, spherical particles in highly turbulent flows.

The simulation approach allows us to examine the role of particle shape in convection-dominated mass transfer at large Péclet numbers ($\textit {Pe} = 1.4 \times 10^{1}$ to $1.4 \times 10^{4}$) and small Stokes number. For spherical particles, we observe an anomalous scaling of the ensemble-average mass transfer rate $\left \langle {\textit {Sh}} \right \rangle \sim \textit {Pe}^{0.26}$, which is consistent with the experimental data for spheres. For other aspect ratios, the best-fit scaling exponent differs, but nonetheless grows less strongly than $\left \langle {\textit {Sh}} \right \rangle \sim \textit {Pe}^{1/3}$. Prolate spheroidal particles are shown to experience a greater mass transfer per unit area than spherical particles of equivalent surface area at all Péclet numbers tested. In contrast, oblate spheroidal particles experience a lower mass transfer per unit area at $\textit {Pe} = O(10)$, but exhibit an optimal aspect ratio $\lambda \approx 1/4$ at large $\textit {Pe}$ where the mass transfer rate is locally maximised.

Statistics of the instantaneous transfer rate show that local fluctuations in $\textit {Sh}$ are attenuated at increasing $\textit {Pe}$. The transfer rate is shown to be most receptive to velocity fluctuations occurring on a time scale $\sim \tau _\eta \textit {Pe}^{1/3}$. This is a demonstration of the convection suppression effect first identified by Batchelor (1979). Defining a local Péclet number $\textit {Pe}_{\widetilde {E}}$ based on a recent time history of the local shear, the conditional transfer rate $\left \langle {\textit {Sh}} | {\textit {Pe}_{\widetilde {E}}} \right \rangle$ is shown to approach the asymptotic scaling $\textit {Pe}_{\widetilde {E}}^{1/3}$ at large $\textit {Pe}_{\widetilde {E}}$. However, the effect of increasing $\textit {Pe}$ is to reduce the magnitude of ${\widetilde {E}}^{\star }$. Jointly, these observations explain the anomalous $\textit {Sh}$ scaling observed in the range of $\textit {Pe}$ tested, because the characteristic shear rate to which the transfer rate is responsive decreases in magnitude with increasing  $\textit {Pe}$.

Based on these observations, we introduce the quasi-steady flux model. This is a particle-scale model of the mass transfer to small tracer particles ${{\textit {St}}} \ll 1$ embedded in unsteady flows at large Péclet numbers $\textit {Pe} \gg 1$. The model approximates the local, unsteady solution of the scalar field surrounding the particle with the steady solution of the convection–diffusion problem based on the recent time history of the relative velocity field. The steady solution is obtained from an interpolation between the purely diffusive solution at $\textit {Pe} = 0$ and the asymptotic solution at $\textit {Pe} \rightarrow \infty$. This simple model is shown to be in good, quantitative agreement with our numerical simulations. The largest discrepancy occurs at low Péclet numbers, where the error introduced by the interpolation procedure is the largest. By construction, the model reproduces the anomalous scaling of the mass transfer rate and reproduces the same characteristic trends in shape dependence for oblate and prolate particles. Furthermore, comparisons of the instantaneous transfer rate predicted by model and simulation are strongly correlated and in good quantitative agreement.

We anticipate that our results may be directly applied to modelling the mass transfer from spheroidal particles in numerical simulations of turbulent flows using point-particle methods. In addition, we envisage several extensions to the present study. For example, our numerical simulation approach and quasi-steady flux model could be extended to model the evolution in shape of a dissolving particle, or the behaviour of particles in anisotropic and wall-bounded turbulent flows where non-spherical particles exhibit different alignment statistics to isotropic turbulence. One application of particular interest is nutrient acquisition by planktonic osmotrophs in turbulent environments, where optimisation of nutrient flux by convection has long been conjectured to drive the adaptation of cell shape and chain formation (Karp-Boss et al. 1996). Here, by incorporating a relative slip velocity into the particle flow field, our approach may quantify the combined effects of particle inertia, shape and sedimentation upon nutrient acquisition. We hope to explore these and related problems in future work.