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 Reference Voth and Soldati2017). Simultaneously, soluble material may be transferred away from the surface by convection and diffusion. Practical examples include industrial processes such as crystallisation (Myerson Reference Myerson2002) or the dissolution of fine solids (Armenante & Kirwan Reference Armenante and Kirwan1989), the transfer of nutrients to planktonic osmotrophs (Karp-Boss, Boss & Jumars Reference Karp-Boss, Boss and Jumars1996), the leaching of pollutants from microplastics (Law Reference Law2017; Seidensticker et al. Reference Seidensticker, Zarfl, Cirpka, Fellenberg and Grathwohl2017) and the encounter rates of bacteria with marine viruses (Guasto, Rusconi & Stocker Reference Guasto, Rusconi and Stocker2012). 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 (Reference Clift, Grace and Weber1978) and Crowe et al. (Reference Crowe, Schwarzkopf, Sommerfeld and Tsuji2012)). 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 Reference Pope2000). 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 Reference Toschi and Bodenschatz2009). 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 Reference Lide2013).
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 Reference Harriott1962; Levins & Glastonbury Reference Levins and Glastonbury1972; Sano, Yamaguchi & Adachi Reference Sano, Yamaguchi and Adachi1974; Asai, Konishi & Sasaki Reference Asai, Konishi and Sasaki1988; Armenante & Kirwan Reference Armenante and Kirwan1989). 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 (Reference Armenante and Kirwan1989) 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. (Reference Asai, Konishi and Sasaki1988) 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 Reference Subramaniam2013; Deen et al. Reference Deen, Peters, Padding and Kuipers2014). Moreover, as recently discussed by Oehmke & Variano (Reference Oehmke and Variano2021), 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 Reference Leal2012) 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 Reference Acrivos and Taylor1962; Acrivos & Goddard Reference Acrivos and Goddard1965; Sehlin Reference Sehlin1969; Gupalo, Polianin & Riazantsev Reference Gupalo, Polianin and Riazantsev1976; Poe & Acrivos Reference Poe and Acrivos1976; Batchelor Reference Batchelor1979; Acrivos Reference Acrivos1980; Myerson Reference Myerson2002; Dehdashti & Masoud Reference Dehdashti and Masoud2020; Lawson Reference Lawson2021). This is supplemented by a wealth of experimental and numerical data at finite Reynolds and Péclet numbers (Clift et al. Reference Clift, Grace and Weber1978; Sparrow, Abraham & Tong Reference Sparrow, Abraham and Tong2004; Kishore & Gu Reference Kishore and Gu2011; Ke et al. Reference Ke, Shu, Zhang, Yuan and Yang2018; Ma & Zhao Reference Ma and Zhao2020).
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 Reference Pozrikidis1997; Feng & Michaelides Reference Feng and Michaelides1998; Michaelides Reference Michaelides2003). At high Péclet number, Batchelor (Reference Batchelor1979) and later Lawson (Reference Lawson2021) 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 Reference Acrivos and Goddard1965; Leal Reference Leal2012; Lawson Reference Lawson2021). 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 Reference Batchelor1979; Lawson Reference Lawson2021), 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 Reference Lawson2021). Therefore, orientation dynamics plays an important role in determining the transfer rate from small particles.
In analogy to time-periodic flows, Batchelor (Reference Batchelor1980) 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 (Reference Armenante and Kirwan1989) and Asai et al. (Reference Asai, Konishi and Sasaki1988) 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. Reference Zhang, Ahmadi, Fan and McLaughlin2001; Mortensen et al. Reference Mortensen, Andersson, Gillissen and Boersma2008; Pumir & Wilkinson Reference Pumir and Wilkinson2011; Parsa et al. Reference Parsa, Calzavarini, Toschi and Voth2012; Byron et al. Reference Byron, Einarsson, Gustavsson, Voth, Mehlig and Variano2015; Zhao et al. Reference Zhao, Challabotla, Andersson and Variano2015; Voth & Soldati Reference Voth and Soldati2017). These preferential alignments appear due to the coupling between the Stokesian dynamics of small particles and the turbulence under Jeffery's equation (Jeffery Reference Jeffery1922). 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 Reference Pumir and Wilkinson2011). In contrast, flattened, disc-like particles align their symmetry axis parallel to the most compressive direction and orthogonal to the vorticity (Voth & Soldati Reference Voth and Soldati2017). 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 Reference Feng and Michaelides2009; Deen et al. Reference Deen, Peters, Padding and Kuipers2014; Derksen Reference Derksen2014) 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 (Reference Harriott1962) 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. Reference Li, Perlman, Wan, Yang, Meneveau, Burns, Chen, Szalay and Eyink2008). 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 Reference Jeffery1922). 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
\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
\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 Reference Lawson2021)
\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
\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 (Reference Kim and Karrila1991). This co-rotating frame and surrounding Stokes flow correspond to the configuration studied by Batchelor (Reference Batchelor1980) 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 Reference Leal2012)
\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
\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
\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 (Reference Lawson2021) 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 (Reference Armenante and Kirwan1989) (
$\triangle$) and Asai et al. (Reference Asai, Konishi and Sasaki1988) (
$\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 (Reference Armenante and Kirwan1989). 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. (Reference Asai, Konishi and Sasaki1988). 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 (Reference Armenante and Kirwan1989), 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. (Reference Asai, Konishi and Sasaki1988) is less good; we note that these data suggest systematically smaller transfer rates than in Armenante & Kirwan (Reference Armenante and Kirwan1989) 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. (Reference Clift, Grace and Weber1978), 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 Reference Lawson2021).
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 Reference Batchelor1979; Lawson Reference Lawson2021), 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
\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 Reference Batchelor1979) 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 Reference Batchelor1979; Lawson Reference Lawson2021).) 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. Reference Yeung, Pope, Lamorgese and Donzis2006). 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
\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
\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 (Reference Batchelor1980) 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
\begin{equation} \widetilde{{\boldsymbol{v}}}^{\star}\, \!\boldsymbol{\cdot} \boldsymbol{\nabla}\theta = \frac{1}{\textit{Pe}} \nabla^{2} \theta, \end{equation}where
\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 Reference Lawson2021) using well-known asymptotic methods (Leal Reference Leal2012). At large Péclet number, the transfer rate may be written as
\begin{equation} \textit{Sh} = \alpha(\widetilde{\boldsymbol{\mathsf{A}}}/\widetilde {E}; \lambda) \textit{Pe}_{\widetilde {E}}^{1/3} + O(1) , \end{equation}where
\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 (Reference Lawson2021) 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. Reference Clift, Grace and Weber1978)
\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
\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. Reference Clift, Grace and Weber1978).
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
\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 (Reference Batchelor1979). 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. Reference Karp-Boss, Boss and Jumars1996). 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.
Acknowledgements
The authors acknowledge the use of the IRIDIS High Performance Computing Facility, and associated support services, at the University of Southampton, in the completion of this work. Pertinent data for this paper are available at doi:10.5258/SOTON/D1976.
Funding
This project has received funding from the European Union's Horizon 2020 research and innovation programme under the Marie Skłodowska-Curie grant agreement no. 846648.
Declaration of interests
The authors report no conflict of interest.
Appendix A. Asymptotic approximation of the quasi-steady mass transfer rate for a spheroid
In this section, we outline the procedure detailed by Lawson (Reference Lawson2021) to obtain the asymptotic solution for the quasi-steady mass transfer rate from a small, torque-free spheroid embedded in a linear velocity gradient at large 
$\textit {Pe}$. This poses the quasi-steady problem (4.1) as
\begin{equation} \boldsymbol{w} \boldsymbol{\cdot} \boldsymbol{\nabla}\theta = \frac{1}{\textit{Pe}_{\widetilde {E}}} \nabla^{2} \theta, \end{equation}
where 
$\boldsymbol {w} = (\tau _\eta \widetilde {E})\widetilde{{\boldsymbol {v}}}^{\star}$, i.e. the time scale for non-dimensionalisation is 
$1/\widetilde {E}$ rather than 
$\tau _\eta$.
The general solution is obtained using a general curvilinear coordinate system 
$(\xi,\eta,\zeta )$ as illustrated in figure 9. This coordinate system is chosen so that the 
$\xi$ coordinate measures the distance normal to the surface, 
$\eta$ coordinate pathlines follow ‘surface streamlines’ tangent to the surface shear 
$\boldsymbol {\tau } = \partial \boldsymbol {w}/\partial \xi |_{\xi =0}$ and 
$\zeta$ coordinate pathlines coincide with contours of constant surface pressure. Thus, on the surface of the spheroid 
$\xi = 0$, lines of constant 
$\zeta$ are surface streamlines and 
$\eta$ is a coordinate along that streamline. It is useful to describe this coordinate system in terms of the covariant coordinate vectors
\begin{equation} \boldsymbol{h}_\xi = \frac{\partial {\boldsymbol{y}}^{{\star}}}{\partial \xi},\quad \boldsymbol{h}_\eta = \frac{\partial {\boldsymbol{y}}^{{\star}}}{\partial \eta},\quad \boldsymbol{h}_\zeta = \frac{\partial {\boldsymbol{y}}^{{\star}}}{\partial \zeta}, \end{equation}
whose inner products define the metric tensor 
$g_{ij}$. Using this coordinate system, the general expression for the transfer rate is
\begin{equation} \textit{Sh} = \frac{0.808 \textit{Pe}_{\widetilde {E}}^{1/3}}{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 + O(1), \end{equation}
where 
$\rho = \sqrt {\det {g}}$ and 
$F$ is the covariant component of the shear 
$\boldsymbol {\tau } \equiv F \boldsymbol {h}_\eta$ along surface streamlines. Physically speaking, 
$\rho$ represents the density of surface area 
$\delta A = \rho \delta \eta \delta \zeta$ of the surface element measuring 
$\delta \eta \times \delta \zeta$. The inner integral in (A3) is carried out along surface streamlines, beginning and terminating at critical points 
$(0,\eta _0,\zeta )$ and 
$(0,\eta _1,\zeta )$ where the surface shear vanishes. The outer integral is a summation across all surface streamlines covering the body.
Figure 9. Illustration of the curvilinear coordinate system 
$(\xi,\eta,\zeta )$ defined on the surface (
$\xi =0$) of a spheroid in an arbitrary velocity gradient. Thick black lines show 
$\eta$-coordinate pathlines (
$\xi =0,\zeta =\textrm {const.}$), which are tangent to the local shear stress on the surface. Thin grey lines are 
$\zeta$-coordinate pathlines (
$\xi = 0, \eta =\textrm {const.}$), which correspond to contours of constant surface pressure. Three fixed points of the surface streamlines are shown: red (source), green (saddle) and blue (sink).
To proceed, we adopt a numerical approach to discretise the surface into an 
$n_\eta \times n_\zeta$ mesh of points 
${\boldsymbol {y}}^{\star }_{i,j} = {\boldsymbol {y}}^{\star }(0, \eta _i, \zeta _j)$. To do this, we numerically integrate a set of 
$j = 1, \ldots, n_\zeta$ streamlines covering the body, label each streamline with a unique 
$\zeta _j$ and evaluate the position at 
$i = 1, \ldots, n_\eta$ different locations 
$\eta _i$ along the streamline. This discretisation allows us to numerically approximate the metric 
$\rho$, which can then be used to numerically integrate (A3).
The surface coordinate system for a spheroid embedded in a linear velocity gradient is defined as follows. The shear at the particle surface is
\begin{equation} \boldsymbol{\tau} = \left.\frac{\partial \boldsymbol{w}}{\partial \xi} \right|_{\xi = 0} = {\varPhi}\boldsymbol{h}_\xi - \boldsymbol{h}_\xi (\boldsymbol{h}_\xi \boldsymbol{\cdot} (\sfbfPhi \boldsymbol{h}_\xi)),\quad \sfbfPhi = \frac{3 r^{3}}{4{\rm \pi} ac^{2}} \frac{\boldsymbol{\mathsf{S}}}{\widetilde {E}} , \end{equation}
where 
$\boldsymbol{\mathsf{S}}$ is the stresslet induced by the particle, whose elements are a linear combination of the components of the velocity gradient 
$\widetilde {\boldsymbol{\mathsf{A}}}$ with geometry-dependent coefficients. Complete expressions for 
$\boldsymbol{\mathsf{S}}$ can be found in Kim & Karrila (Reference Kim and Karrila1991). The surface normal and streamwise covariant coordinate vectors are
\begin{equation} \boldsymbol{h}_\xi = \frac{\boldsymbol{\mathsf{D}}{\boldsymbol{y}}^{{\star}}}{|\boldsymbol{\mathsf{D}}{\boldsymbol{y}}^{{\star}}|},\quad \boldsymbol{h}_\eta = \frac{\partial {\boldsymbol{y}}^{{\star}}}{\partial \eta} = \frac{\boldsymbol{\tau}}{\boldsymbol{\tau} \boldsymbol{\cdot} \boldsymbol{\nabla} \eta}, \end{equation}
where (in the body system) 
$\boldsymbol{\mathsf{D}}$ is a diagonal matrix whose entries are 
$a^{-2}$, 
$c^{-2}$ and 
$c^{-2}$. A suitable definition for the streamwise coordinate is
\begin{equation} \eta(\,{\boldsymbol{y}}^{{\star}}) = \boldsymbol{h}_\xi \boldsymbol{\cdot} (\sfbfPhi \boldsymbol{h}_\xi), \end{equation}
which happens to correspond to the surface pressure perturbation. A fundamental property of this definition is that 
$\eta$ is monotonic, increasing along streamlines from 
$\eta _0 = \phi _1$ to 
$\eta _1 = \phi _3$, where 
$\phi _1 \leqslant \phi _2 \leqslant \phi _3$ are the eigenvalues of 
$\sfbfPhi$ associated with eigenvectors 
$\boldsymbol {q}_i$. Provided these eigenvalues are distinct, 
$\boldsymbol {h}_\xi \rightarrow \pm \boldsymbol {q}_1$ as 
$\eta \rightarrow \phi _1$ and 
$\boldsymbol {h}_\xi \rightarrow \pm \boldsymbol {q}_3$ as 
$\eta \rightarrow \phi _3$, which correspond to critical points upon the surface in figure 9 where the surface shear stress vanishes. This divides the surface into four quadrants, depending upon the basin of attraction of each critical point. Since 
$\zeta$ is constant along streamlines and every streamline passes through 
$\eta = \phi _2$, 
$\zeta$ can be thought of as a coordinate along the curve 
${\boldsymbol {y}}^{\star }_0(\zeta ) = {\boldsymbol {y}}^{\star }(0, \phi _2, \zeta )$. The coordinate 
$\zeta = \boldsymbol {h}_\xi \boldsymbol {\cdot } \boldsymbol {q}_2$ is then sufficient to define a point 
$(\eta, \zeta ) \in [\phi _1,\phi _3] \times [-1,1]$ on each quadrant.
A technical point remains that 
$\eta _i$ and 
$\zeta _j$ should be chosen so that the surface is densely covered in mesh points 
${\boldsymbol {y}}^{\star }_{i,j}$. This is achieved by generating a uniform sampling of seed points on the surface, which are integrated numerically (A2) towards 
$\eta = \phi _2$ to construct 
$\zeta _j$. The 
$\eta _i$ are chosen as
\begin{equation} \eta_i = \begin{cases} \phi_1 \cos^{2}\psi_i + \phi_2 \sin^{2}\psi_i, & \text{ if } \psi_i < 0, \\ \phi_2 \cos^{2}\psi_i + \phi_3 \sin^{2}\psi_i, & \text{ if } \psi_i > 0, \end{cases} \end{equation}
for 
$\psi _i$ evenly distributed on the interval 
$[-{\rm \pi} /2,{\rm \pi} /2]$. We find 
$n_\eta = 513$ and 
$n_\zeta = 564$ to be sufficient to approximate (A3).
Appendix B. Parametrisation of the local shear
For a torque-free particle in a steady linear shear, the relative velocity gradient 
$\boldsymbol{\mathsf{A}}$ can be specified by the rate-of-strain tensor 
$\boldsymbol{\mathsf{E}}$ in the body frame. This can be seen by noting that, in this co-rotating frame, the relative vorticity vector must satisfy
\begin{equation} 0 = \frac{1}{2}\boldsymbol{\omega} + \frac{\lambda^{2}-1}{\lambda^{2}+1} \boldsymbol{i} \times \boldsymbol{\mathsf{E}}\boldsymbol{i}, \end{equation}
where 
$\boldsymbol {i}$ denotes the orientation of the particle's symmetry axis in the body frame. As a result, the relative velocity field can be described completely in terms of the eigenvalues 
$\sigma _i$ of the rate-of-strain tensor 
$\boldsymbol{\mathsf{E}}$ and the orientation of the spheroid with respect to its eigenvectors 
$\boldsymbol {e}_i$. The local Péclet number 
$\textit {Pe}_{E} = r^{2} |\boldsymbol{\mathsf{E}}|/\kappa$ can therefore be defined in terms of the magnitude of the local strain rate 
$|\boldsymbol{\mathsf{E}}|^{2} = {\mathsf{E}}_{ij}{\mathsf{E}}_{ij}$, which is an invariant across reference frames. Continuity requires 
$\sum \sigma _i = 0$, so that, up to a choice of scale described by 
$\textit {Pe}_{E}$, the eigenvalues of 
$\boldsymbol{\mathsf{E}}$ can be parametrised by (Lund & Rogers Reference Lund and Rogers1994)
\begin{equation} s ={-}\frac{3\sqrt{6} \, \sigma_1 \sigma_2 \sigma_3}{(\sigma_1^{2} + \sigma_2^{2} + \sigma_3^{2})^{3/2}}, \end{equation}
where 
$-1 \leqslant s \leqslant 1$.
The orientation of the spheroid with respect to the strain eigenvectors can be parametrised in terms of spherical polar coordinates. Identifying the eigenvectors by the ordering 
$\sigma _1 \geqslant \sigma _2 \geqslant \sigma _3$, we may choose to write
\begin{equation} \boldsymbol{i} = \sin\psi\cos\phi \sin\phi \,\boldsymbol{e}_1 + \sin\psi\sin\phi ~\boldsymbol{e}_2 + \cos\psi \,\boldsymbol{e}_3, \end{equation}
so that 
$\psi$ is the angle between the symmetry axis and the most compressive strain direction 
$\boldsymbol {e}_3$ and 
$\phi$ is the azimuth measured from the most extensional strain direction 
$\boldsymbol {e}_1$. Owing to axisymmetry, the orientation of the equatorial axes does not affect the mass transfer rate, which is specified by 
$s, \psi, \phi, \lambda$ and 
$\textit {Pe}_{{E}}$. Because of fore–aft symmetry, it is only necessary to consider orientations 
$\psi \in [0, {\rm \pi}/2]$ and 
$\phi \in [0, {\rm \pi}]$ in the computation of 
$\alpha (s,\psi,\phi ;\lambda )$ in (4.3).
Open access
