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.

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 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 Re λ (Pope 2000). The particle Stokes number St ≡ τ p /τ η parametrises the response time of the particle τ p with respect to the Kolmogorov time scale τ η of the turbulence, characteristic of the average shear rate experienced by the particle (Toschi & Bodenschatz 2009). When St 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 Pe ≡ r 2 /κτ η , where κ is the mass diffusion coefficient of the solute. When Pe 1, convection dominates the mass transfer away from the particle. Sometimes, these ratios are alternatively described by a particle Reynolds number Re = r 4/3 1/3 /ν ∝ St 2/3 and Schmidt number Sc = ν/κ, where is the mean kinetic energy dissipation rate and ν is the kinematic viscosity of the fluid. Jointly, the assumption of large Pe and small St implies that the Schmidt number is large; this is in fact often the case for solutes in liquids, e.g. Sc ≈ 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 St 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 Re < 0.1. They therefore conducted a comprehensive set of experiments to determine the mass transfer rate of small, spherical particles with St 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 Sh = αPe 1/3 + O(1) when Pe 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 α 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 u(x) = E 0 x + 1 2 ω 0 × x, unsteady concentration fluctuations arise due to the periodic motion of the particle under the action of the strain E 0 and vorticity ω 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 v( y, t), periodic in time t, where y is the spatial coordinate in the co-rotating body frame. Remarkably, the concentration fluctuations scale as Pe −1/3 in magnitude relative to the mean, provided that the period T of the motion is much smaller than |E 0 | −1 Pe 1/3 . As such, the mean flow field in the body frame of the particleṽ(t) = (1/T) T 0 v dt 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 |ω 0 | |E 0 |, the transfer rate is shown to scale with a Péclet number based on the rate of vortex stretching (Batchelor 1979;Lawson 2021), since particles spin aligned with the vorticity. More generally, for axisymmetric particles spinning along a fixed direction p, the transfer rate scales with the strain rate E p = E 0 ij 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 ω . This argument relies upon averaging the flow in the body frame over a long time with respect to τ η , which is permissible in the limit of Pe → ∞, since the concentration boundary layer should be insensitive to velocity fluctuations faster than τ η Pe 1/3 . Since the mean vortex stretching scales with 1/τ η , this argument predicts that the transfer rate scales as Sh ≈ 0.55Pe 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 929 A19-3 boundary layer is sensitive may therefore depend upon Pe, so an additional Pe dependence may be introduced before the asymptotic behaviour is reached. This brings the application of the asymptotic scaling Sh ∼ Pe 1/3 to large but finite 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 St 1 and Pe 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 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.

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 λ ≈ 433 obtained from the Johns Hopkins University (JHU) Turbulence Database (Li et al. 2008). Lagrangian tracers are uniformly seeded at positions x 0 throughout the simulation domain at time t = 0 and their trajectories X (t; x 0 ) are integrated as tracersẊ (t; x 0 ) = u(X (t; x 0 ), t) following the flow. We then obtain the time history of the local velocity gradient G(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 λ = a/c is centred at X and the orientation of its principal semi-axes of length (a, c, c) are described by the rotation matrix R = [p, q, r]. Jeffery showed that the solid-body rotation rate Ω of a torque-free spheroid embedded in an arbitrary linear shear is given by where ω 0 i = − ijk G jk and E 0 = (G T + G)/2 are the vorticity and rate-of-strain tensor perceived by the particle in the laboratory frame. The body frame therefore rotates aṡ so that the relative velocity gradient experienced by the particle in the co-rotating frame is (Lawson 2021 As before, we apply a similar decomposition of the relative velocity gradient tensor A into a symmetric strain tensor E = (A + A T )/2 and perceived vorticity vector ω in the co-rotating frame. The relative velocity v experienced in the vicinity of the particle is therefore where y is the spatial coordinate in the co-rotating body frame and v is the Stokes perturbation to the velocity field due to the presence of the particle. The complete expression for 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 θ( y , t ) can be written in dimensionless form as (Leal 2012) ∂θ ∂t where v = vτ η /r is the non-dimensionalised relative velocity, t = t/τ η and y = y/r are the dimensionless spatial and temporal coordinates, and r is the linear dimension of the particle. The turbulent Péclet number Pe ≡ r 2 /κτ η = r 2 1/2 /κν 1/2 therefore forms a control parameter for the convection-diffusion problem, where κ is the diffusion coefficient of the solute and τ η ≡ (ν/ ) 1/2 is the Kolmogorov time scale. Because this frame is co-rotating, 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 which corresponds to non-dimensionalising the concentration field as C( 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 Conventions differ on the choice of characteristic length scale r; here, we adopt the definition r ≡ √ A/4π, where A is the surface area of the particle, so that Sh = 1 for purely diffusive flux of material from a spherical particle.

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 x 0 at t = 0 and their trajectories X (t; 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 t = 0.002 ≈ 0.047τ η , 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 G(X , t i ) at i = 0, . . . , 5028 times t i = i t. The orientation of the body frame R(t) (2.2) is integrated using an adaptive time-step, third-order Runge-Kutta scheme using a piecewise linear interpolation to interpolate G(X , t). The initial condition R 0 = 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 (μ, ψ, φ), depending upon the particle aspect ratio. In the spheroidal system, the μ direction is analogous to a radial coordinate, the ψ direction parametrises the polar angle from the symmetry axis, and φ is an azimuthal angle. The inner boundary at μ = μ 0 corresponds to the surface of the spheroid, whereas the outer boundary at μ = μ N corresponds to fluid far from the particle. The dimensions of the spheroid are chosen such that the surface area is 4π, 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 × 64 × 64 cells in the (μ, ψ, φ) directions, respectively, with uniform spacing in the ψ and φ directions. To adequately resolve the thin concentration boundary layer, which is of thickness δ, we employ a mesh refinement in the μ direction such that the grid spacing is is the spacing between adjacent cells in the μ direction. Owing to the nature of the spheroidal coordinate system chosen, the resolution varies across the surface of the particle. The initial spacing μ 1 is chosen such that the thickness of the largest cell near the particle surface is at most 2 × 10 −4 min(a, c) and the mesh refinement factor R is chosen accordingly. Based on an estimated boundary layer thickness δ = Pe −1/3 at the largest Péclet number tested, for the most extreme aspect ratios λ = 16 (λ = 1/16), this yields between 26 (45) and 68 (82) cells within a distance δ ≈ 0.0414 from the surface.
We impose the Dirichlet boundary condition θ = 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 t = 0.001 from the initial condition θ( y , 0) = 0. Statistics of the transfer rate are gathered in the interval 75τ η t T sim , which is sufficiently long to permit the statistics of the transfer rate to reach a steady state.

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 ∼75τ η . Averaging over the interval 75τ η t T sim , the deviation of the ensemble-average Sherwood number Sh from the time average of the ensemble Sh 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 Sh to be below ±0.20 % (based on a 95 % confidence interval).
As a validation of our numerical model, we compare the average Sherwood number predicted by the model against experimental data of small spheres with 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 Sc, and a comparison across datasets at fixed Pe = (r/η) 2 Sc, where η 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 Sh < 1. The dotted line in figure 1(b) shows the best power-law fit Sh = 0.99Pe 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 Pe 1/3 result) and shows a weaker dependence upon Pe than expected. We shall examine the physical origin of this anomalous scaling in § 3.3.

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 λ = 2 n , n ∈ {−4, −2, 0, 2, 4}. We observe that, with the possible exception of very strongly flattened spheroids (λ = 1/16), the slope of the power-law scaling Sh ∼ Pe γ at large Pe appears to be less than 1/3, consistent with the anomalous scaling behaviour seen for spheres.
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 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 Pe is increased. However, there is a continuous evolution in the shape dependence with increasing 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 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 Pe = 0 (Clift et al. (1978), not shown). However, at Pe = O(10 3 ), this trend is reversed for spheroids near an optimal aspect ratio λ ≈ 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 λ = 16 experiences roughly ∼31 % greater flux than an equivalent sphere at Pe = 1.4 × 10 1 , but 42 % greater flux at Pe = 1.4 × 10 4 . The relative advantage is smaller for less extreme differences in aspect ratio and Pe. We note that the same qualitative behaviour is observed for spheroids in steady, rotation-dominated (|ω 0 | |E 0 |) flow, where oblate particles also exhibit maximal mass transfer near λ ≈ 0.31 under vortex compression at large Pe (Lawson 2021).

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 Sh = Sh − Sh relative to the mean transfer rate for different aspect ratios and Péclet numbers. Across the range of Pe tested, the fluctuation magnitude is between 12 % and 19 % of the mean mass transfer rate, with a peak around Pe = O(10 3 ). Thus, turbulent fluctuations in the mass transfer rate are relatively small in comparison to the mean and decay at large 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 Pe, the spectrum rolls off around f τ η ≈ 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 T Pe 1/3 result in a negligible contribution to the transfer rate at large Pe.
To examine the temporal response of the concentration boundary layer in more detail, we consider the cross-correlation R SP (τ ) between the instantaneous transfer rate Sh(t + τ ) at time t + τ and the local Péclet number Pe E (t) ≡ r 2 |E|/κ at time t. Here, |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 (|E 0 | = |E|). Figure 4(a) shows the behaviour of this cross-correlation as a function of the time lag τ for spheres at varying turbulent 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 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.
To identify this response time scale, we consider the correlation between the transfer rate and the local strain rate filtered at time scale τ f . We therefore introduce the filtered, relative velocity gradient Thus, A is the average velocity gradient recently experienced by the particle over the time scale τ f . Figure 5(a) shows the peak correlation coefficient between the mass transfer rate of spherical particles and the local Péclet number Pe E ≡ r 2 E/κ, which is based on the magnitude of the filtered strain rate E 2 = E ij E ij . The data for spheroids with other aspect ratios are qualitatively similar. We note that, under this filtering operation, | E 0 | / = | 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 Pe, there is an optimal filter time scale τ d where the transfer rate correlates most strongly with the filtered velocity gradient. Furthermore, this correlation is strong: for spheres, fluctuations in E are sufficient to explain roughly 80 % of the variance in 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 τ d as the response time scale of the concentration boundary layer to velocity gradient fluctuations. Figure 5(b) shows the dependence of the response time scale τ d upon Pe for all aspect ratios tested. In all cases the response time scale is well approximated by a power-law fit τ d /τ η = βPe γ , which is illustrated for spheres. In analogy to time-periodic flows, we might expect τ d ∼ Pe 1/3 . The best-fit exponents for λ = 1/16, . . . , 16 are γ = 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 γ = 1/3 also provides an adequate fit to the data, with β = 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(Pe 1/3 τ η ). Alternatively, this can be interpreted as an increase in sensitivity of the concentration boundary layer to lower-frequency motions at larger 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 Pe. However, this dependence is weaker than the 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 Pe E = Pe Eτ η in figure 6(a), where the local filter time scale is chosen at the correlation optimum τ f = τ d . We observe that, when the local Péclet number is large relative to the mean, the conditional transfer rate Sh|Pe E collapses onto a single curve independently of 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 A 11 /2 = − A 22 = − A 33 = E/ √ 6 and approaches the asymptotic scaling Sh ≈ 0.905Pe 1/3 E (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 Pe E , the transfer rate differs by less than 1 % in comparison to axisymmetric strain (Batchelor 1979;Lawson 2021).) A similar Pe 1/3 E scaling of the conditional transfer rate can also be observed for prolate and oblate spheroids (not shown).
The distribution of the local shear rate E = Eτ η 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 E , which is consistent with unfiltered measurements of the local kinetic energy dissipation rate = 2νE 2 (Yeung et al. 2006). As Pe and therefore τ d increases, the effective shear rate E 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 Pe observed. At large Pe E , the conditional transfer rate scales as Pe 1/3 E . However, the convection suppression effect attenuates the magnitude of the turbulent shear fluctuations to which the particle is responsive as Pe increases. Therefore, the average transfer rate grows less strongly than Sh ∼ Pe 1/3 , because the effective shear rate the boundary layer perceives decreases at increasing Pe. This can be seen more formally through the statistical identity where we have substituted the scaling For spheres, we find ( E ) 1/3 ∼ Pe −0.053 over the range of 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 ( E ) 1/3 in the limit Pe → ∞, which is beyond the scope of this study. The result of Batchelor (1980) is compatible with (3.2) provided ( E ) 1/3 = O(1) at Pe → ∞, as Batchelor's analysis of the average relative flow field predicts.

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 A, rather than the instantaneous relative velocity gradient A. We therefore hypothesise that the local solution to (2.5) is approximated by the steady-state solution to 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 τ d . The filter time scale τ f should therefore be matched to the diffusive time scale τ 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 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 where α = 0.808 4π The expression for the prefactor α involves an integration of the dimensionless surface shear rate F over the particle's surface using a general curvilinear coordinate system (ξ, η, ζ ) whose metric tensor g ij defines ρ = √ 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 α is large; A has nine degrees of freedom in addition to the aspect ratio λ. Firstly, we note that the rescaling by the local shear rate E in (4.3) eliminates the dependence of α upon the magnitude of A. Continuity, axisymmetry and the torque-free condition (2.1) further constrain five degrees of freedom such that the prefactor α = α(s, ψ, φ; λ) 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 ψ and φ describing the orientation of the spheroid relative to the strain field, and the aspect ratio λ. We provide a complete definition of these geometric parameters in Appendix B. This parameter space is sufficiently small to allow α 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) Interpolating between these two limits, we approximate which provides a close approximation (within 10 %) of the transfer rate for prolate and oblate spheroids in uniform, creeping flow with Pe < 70 (Clift et al. 1978).
The model is therefore specified by (4.6), up to the dependence of the filter time scale τ f upon Pe. On the basis of the arguments presented in the previous section, we choose where β = 0.96 is the best-fit coefficient to the boundary layer response time scale in § 3.3. In principle, we might choose to make β a function of λ 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 β 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 (∼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 λ ≈ 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(b) shows the relative magnitude of transfer rate fluctuations σ Sh = Sh 2 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 Sh fluctuations is slightly overestimated.

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 (Pe = 1.4 × 10 1 to 1.4 × 10 4 ) and small Stokes number. For spherical particles, we observe an anomalous scaling of the ensemble-average mass transfer rate Sh ∼ 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 Sh ∼ 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 Pe = O(10), but exhibit an optimal aspect ratio λ ≈ 1/4 at large Pe where the mass transfer rate is locally maximised.
Statistics of the instantaneous transfer rate show that local fluctuations in Sh are attenuated at increasing Pe. The transfer rate is shown to be most receptive to velocity fluctuations occurring on a time scale ∼ τ η Pe 1/3 . This is a demonstration of the convection suppression effect first identified by Batchelor (1979). Defining a local Péclet number Pe E based on a recent time history of the local shear, the conditional transfer rate Sh|Pe E is shown to approach the asymptotic scaling Pe 1/3 E at large Pe E . However, the effect of increasing Pe is to reduce the magnitude of E . Jointly, these observations explain the anomalous Sh scaling observed in the range of Pe tested, because the characteristic shear rate to which the transfer rate is responsive decreases in magnitude with increasing 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 St 1 embedded in unsteady flows at large Péclet numbers Pe 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 Pe = 0 and the asymptotic solution at Pe → ∞. 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.
Appendix A. Asymptotic approximation of the quasi-steady mass transfer rate for a spheroid In this section, we outline the procedure detailed by Lawson (2021) 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 Pe. This poses the quasi-steady problem (4.1) as where w = (τ η E) v , i.e. the time scale for non-dimensionalisation is 1/ E rather than τ η . The general solution is obtained using a general curvilinear coordinate system (ξ, η, ζ ) as illustrated in figure 9. This coordinate system is chosen so that the ξ coordinate measures the distance normal to the surface, η coordinate pathlines follow 'surface streamlines' tangent to the surface shear τ = ∂w/∂ξ | ξ =0 and ζ coordinate pathlines coincide with contours of constant surface pressure. Thus, on the surface of the spheroid ξ = 0, lines of constant ζ are surface streamlines and η is a coordinate along that streamline. It is useful to describe this coordinate system in terms of the covariant coordinate vectors h ξ = ∂y ∂ξ , h η = ∂y ∂η , h ζ = ∂y ∂ζ , (A2a-c) whose inner products define the metric tensor g ij . Using this coordinate system, the general expression for the transfer rate is where ρ = √ det g and F is the covariant component of the shear τ ≡ Fh η along surface streamlines. Physically speaking, ρ represents the density of surface area δA = ρδηδζ of the surface element measuring δη × δζ . The inner integral in (A3) is carried out along surface streamlines, beginning and terminating at critical points (0, η 0 , ζ ) and (0, η 1 , ζ ) where the surface shear vanishes. The outer integral is a summation across all surface streamlines covering the body.
To proceed, we adopt a numerical approach to discretise the surface into an n η × n ζ mesh of points y i,j = y (0, η i , ζ j ). To do this, we numerically integrate a set of j = 1, . . . , n ζ streamlines covering the body, label each streamline with a unique ζ j and evaluate the position at i = 1, . . . , n η different locations η i along the streamline. This discretisation allows us to numerically approximate the metric ρ, 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 where S is the stresslet induced by the particle, whose elements are a linear combination of the components of the velocity gradient A with geometry-dependent coefficients.
Complete expressions for S can be found in Kim & Karrila (1991). The surface normal and streamwise covariant coordinate vectors are where (in the body system) D is a diagonal matrix whose entries are a −2 , c −2 and c −2 . A suitable definition for the streamwise coordinate is which happens to correspond to the surface pressure perturbation. A fundamental property of this definition is that η is monotonic, increasing along streamlines from η 0 = φ 1 to η 1 = φ 3 , where φ 1 φ 2 φ 3 are the eigenvalues of associated with eigenvectors q i . Provided these eigenvalues are distinct, h ξ → ±q 1 as η → φ 1 and h ξ → ±q 3 as η → φ 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 ζ is constant along streamlines and every streamline passes through η = φ 2 , ζ can be thought of as a coordinate along the curve y 0 (ζ ) = y (0, φ 2 , ζ ). The coordinate ζ = h ξ · q 2 is then sufficient to define a point (η, ζ ) ∈ [φ 1 , φ 3 ] × [−1, 1] on each quadrant.
A technical point remains that η i and ζ j should be chosen so that the surface is densely covered in mesh points y i,j . This is achieved by generating a uniform sampling of seed points on the surface, which are integrated numerically (A2) towards η = φ 2 to construct ζ j . The η i are chosen as for ψ i evenly distributed on the interval [−π/2, π/2]. We find n η = 513 and n ζ = 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 A can be specified by the rate-of-strain tensor E in the body frame. This can be seen by noting that,