3.1. Experimental methodology

Here, we describe our experimental set-up and our procedure for tracking particles and reconstruct the concentration field surrounding the particle. Our experimental set-up is illustrated in figure 2(a). It consists of three main components: a randomly stirred mixing tank to generate homogeneous turbulence, dye-laden particles and a combined shadowgraph and PLIF imaging set-up to track particles to obtain particle-scale measurements of dye release. We now describe each.

Figure 2.

(a) Photograph of experimental apparatus, illustrating: impeller-stirred mixing tank; high-speed shadowgraph imaging arrangement including cameras, prisms and collimated backlight sources; and PLIF laser, optics and sheet. (b) Visualisation of the concentration boundary layer and wake around fine particles in turbulence, showing PLIF image (green channel) and shadowgraph image of particle (purple channel). Supplementary movies 1 and 2 available at https://doi.org/10.1017/jfm.2022.998 show sample PLIF/shadowgraph video sequences.

3.1.2. Dye-laden particles

To create spherical dye-laden particles, we load sodium fluorescein dye (Sigma Aldrich) onto commercially available anion-exchange resin beads (Dowex $1{\text{X}}2$, Fischer Scientific). The beads, which are supplied in $\textrm {Cl}^-$ form, undergo a reversible reaction as

(3.1)\begin{equation} \textrm{R}\textrm{Cl} + \textrm{Na}\textrm{Fl} \rightleftharpoons \textrm{R}\textrm{Fl} + \textrm{Na}\textrm{Cl}. \end{equation}

In a typical procedure, we prepare $100\,\textrm {ml}$ of a $0.2\,\textrm {M}$ solution of sodium fluorescein to which $15\,\textrm {g}$ of ion-exchange beads is added under stirring. This loading of dye is close to the capacity of the ion-exchange beads and the dye solution is visibly depleted afterwards. We have observed that sodium fluorescein is not absorbed well onto resin beads in $\textrm {OH}^-$ form.

Relevant physical properties of the dye-laden particles are presented in table 1. The mean particle radius $\left \langle a \right \rangle$ and standard deviation $\sigma _a$ were obtained by processing microscope images taken at a resolution of $400\,\textrm {px}\,\textrm {mm}^{-1}$ using a Hough-transform circle-detection routine. Statistics were based upon samples of 381 (50–100 mesh), 1350 (100–200 mesh) and 5152 (200–400 mesh) particles. While the larger 50–100 mesh particles exhibited little variation in size, the smallest 200–400 mesh particles had a relatively large coefficient of variation $\sigma _a/\left \langle a \right \rangle$ of around $24\,\%$. Particle density was measured in triplicate using a $5\,\textrm {ml}$ pycnometric flask and an analytical balance with $10\,\textrm {mg}$ resolution. The quiescent settling velocity in water at $20\,^\circ \textrm {C}$ was calculated from

(3.2)\begin{equation} w_0 = \sqrt{\frac{g \left\langle a \right\rangle}{3 C_D } \frac{\rho_p - \rho_f}{\rho_f}} \end{equation}

based on the empirical drag coefficient relation

(3.3)\begin{equation} C_D = \frac{24}{{\it Re} _{w}} \left( 1 + 0.1315 {\it Re} _{w}^{0.82 - 0.05\log_{10} {\it Re} _{w}} \right) , \end{equation}

which is applicable in the Reynolds number range ${\it Re} _{w} = 0.01$–$20$ (Clift, Grace & Weber Reference Clift, Grace and Weber1978).

Table 1.

Properties of fluorescent dye-laden particles. Uncertainty in particle density is ${\pm }2\,\textrm {kg}\,\textrm {m}^{-3}$. The quiescent settling velocity $w_0$ in water at $20\,^\circ \textrm {C}$ was calculated based on the particle density $\rho _p$ and mean radius $\left \langle a \right \rangle$ using (3.2).

3.2. Particle tracking and wake reconstruction

Table 2 summarises the experimental conditions tested. From our previous PIV characterisation (Lawson & Ganapathisubramani Reference Lawson and Ganapathisubramani2022) and the particle properties in table 1, we compute the particle Stokes number ${St}$, turbulent Péclet number ${Pe}$, sinking ratio ${Sr}$ and size relative to the Kolmogorov scale. Although the turbulent Péclet number is not controlled for in our experiments, the values are well within the convection-dominated regime (${Pe} \gg 1$). It can be seen that the particles in our experiments are comparable to or smaller than the smallest dynamically active motions of the turbulent flow. For all cases, the characteristic turbulent acceleration $a_\eta = (\epsilon ^3/\nu )^{1/4}$ is at least one order of magnitude smaller than the strength of gravitational acceleration, so slip due to turbulent acceleration is negligible. Furthermore, for the two smallest classes of particles, ${St} \ll 1$. By varying the turbulence intensity, we varied the sinking ratio ${Sr}$ by an order of magnitude over the range in which the transition in mass transfer mechanism is expected to occur. For each condition, we obtained time series of 2738 (5000 at 2.4 kHz frame rate) shadowgraph and PLIF images, corresponding to between 1.9 and 6.8 s of turbulent motion, equivalent to approximately 1.3–3.2 integral time scales.

Table 2.

Summary of experimental conditions.

Lagrangian particle tracking was performed using an in-house code based on Lawson et al. (Reference Lawson, Bodenschatz, Lalescu and Wilczek2018). This implements a low-density, predictor–corrector particle tracking algorithm following Ouellette, Xu & Bodenschatz (Reference Ouellette, Xu and Bodenschatz2006). The main difference is in particle image registration. Shadowgraph images are first levelised against a white image obtained from the $95$th brightness percentile on a per-pixel basis. Subsequently, particle centres are identified using a Hough circle transform. Existing particle tracks are extrapolated linearly based on the last two frames and matched to the nearest corresponding particle centres (within a search radius of $20\,\textrm {px}$) and satisfying a triangulation tolerance ($5\,\textrm {px}$). New tracks are initiated by triangulating unmatched particles and are extrapolated using a nearest-neighbour interpolation of the particle velocity field, with a larger search tolerance. Tracks are terminated if a suitable match is not found or if matching is ambiguous.

We apply Taylor's frozen flow hypothesis to reconstruct the scalar field in the vicinity of particles as they cross the laser sheet. In order to do this, we identify when particles cross the laser sheet based on their position and a calibrated model of the laser sheet plane. This plane is identified using the procedure illustrated in figure 3. When particles cross the laser sheet, we observe a sharp increase in their brightness due to their fluorescence as well as specular reflections, as shown in figure 3(a). We identify a set of crossing events based on a threshold, which define a set of points known to lie upon the laser sheet plane, as illustrated in figure 3(b). We use a robust fitting procedure to fit a plane to these points.

Figure 3.

Identification of the laser sheet plane from particle trajectories. (a) An example of time-series measurements of particle brightness, showing the identification of particles which cross the laser sheet. (b) A 3-D reconstruction of the laser sheet plane (green), based on the position of identified crossings (black circles) of particle trajectories. Particle trajectories have been coloured according to their brightness. The white dotted line indicates the particle radius.

To reconstruct the 3-D scalar field, we identify the time $t_0$ when a particle trajectory $\boldsymbol {X}(t)$ crosses the laser sheet plane $\boldsymbol {X}(t_0)\boldsymbol {\cdot}\boldsymbol {e}_n + c = 0$. The PLIF images $I_{2D}(\boldsymbol {S}+\boldsymbol {s}, t)$ are then sampled in the vicinity of the projected position of the particle $\boldsymbol {S} = \boldsymbol {P}(\boldsymbol {X})$, where $\boldsymbol {s}$ is the two-dimensional (2-D) image coordinate relative to the particle centre and $\boldsymbol {P}(\boldsymbol {x})$ represents the projection operation from 3-D object space to the 2-D image space of the camera. Images are sampled over $\pm 5$ frames around each crossing event. Combined with the constraint that the fluorescence is emitted in the laser sheet plane, this allows us to reconstruct 11 planar slices of the scalar field $I_j(\boldsymbol {r})$ in the co-moving particle frame (2.5) for each crossing event $j$.

4.1. Numerical methodology

The numerical procedure to implement the model is similar to that of Lawson & Ganapathisubramani (Reference Lawson and Ganapathisubramani2021). We summarise the details here only briefly. The trajectories of particles sinking in homogeneous, isotropic turbulence at $R_\lambda = 433$ are obtained from the Johns Hopkins Turbulence Database (Li et al. Reference Li, Perlman, Wan, Yang, Meneveau, Burns, Chen, Szalay and Eyink2008). Ensembles of $1000$ Lagrangian particles are uniformly seeded throughout the simulation domain at time $t=0$ and their trajectories are integrated as

(4.1)\begin{equation} \dot{\boldsymbol{X}} = \boldsymbol{u}(\boldsymbol{X}(t), t) + \boldsymbol{w}_0, \end{equation}

where the settling velocity is given by (2.3). Equation (4.1) is integrated over the duration of the simulation $T_{sim} = 10.056$ using a second-order Runge–Kutta scheme with a time step ${{\rm \Delta} t} = 0.002 \approx 0.047\tau _\eta$, corresponding to the time resolution of the database. The velocity and velocity gradient are interpolated from the database at the tracer position using a fourth-order Lagrange polynomial interpolation. This yields a time series of the local velocity gradient $\boldsymbol{\mathsf{G}}(\boldsymbol {X},t_i)$ at $i = 0 \ldots 5028$ times $t_i = i{\rm \Delta} t$. The orientation of the body frame $\boldsymbol{\mathsf{R}}(t)$ (2.6) is integrated from an initially random (isotropic) orientation using an adaptive time step, third-order Runge–Kutta scheme using a piecewise, linear interpolation to interpolate $\boldsymbol{\mathsf{G}}(\boldsymbol {X},t)$. This provides time series of the particle–fluid relative velocity $\boldsymbol {v}_0(t)$ and gradient $\boldsymbol{\mathsf{E}}(t)$ in the co-rotating frame.

We track the motion (4.1) of seven separate ensembles of spherical particles with the same initial position but different sinking rates $w_0/u_\eta = (a/\eta ){Sr} = 10^{-3}\unicode{x2013}10^{0}$. Although $w_0/u_\eta$ is small, the difference in sinking rate can build up into a large difference in trajectory over time. From these we obtain the Lagrangian relative velocity history (2.11) of two sizes of small spheres $a/\eta = 10^{-2}$ and $10^{-1}$ at differing sinking ratio $\log _{10}{Sr} = -1, -0.5, 0, 0.5, 1$.

We then solve the convection diffusion problem (2.8) numerically using a modified version of OpenFOAM scalarTransportFoam solver used in Lawson (Reference Lawson2021) and Lawson & Ganapathisubramani (Reference Lawson and Ganapathisubramani2021). This implements a second-order finite-volume solver for $\theta$, with a standard linear upwind Gaussian discretisation of the convective term in (2.8) and linear Gaussian scheme for the diffusive term with explicit non-orthogonal correction to maintain second-order accuracy. Convective forcing (2.11) is implemented by specifying the relative velocity field as the superposition of Stokes flow solutions, obtained from analytical expressions given in Kim & Karrila (Reference Kim and Karrila1991). Equation (2.8) is discretised onto a structured grid in spherical coordinates $(r,\psi,\phi )$ with $150\times 64\times 64$ cells, where $r$ is the radial coordinate, $\psi$ is the polar coordinate and $\phi$ the azimuthal. Cells are spaced uniformly in the azimuthal direction $\phi _i = 0\dots 2{\rm \pi}$ and polar direction $\cos \psi _i = -1 \dots 1$. To adequately resolve the thin concentration boundary layer, we employ a mesh refinement ${\rm \Delta} r_{i+1} = \alpha r_i$, where ${\rm \Delta} r_i = r_{i+1}-r_i$ is the spacing between cells in the radial direction. We choose $\alpha = 1.0724$. Based on an estimated boundary layer thickness $\delta = {Pe}_w^{-1/3}$ at the largest Péclet number (${Pe}_w = 10^5, {Pe} = 10^4$) tested, there are at least $32$ cells within a distance $\delta$ from the surface. We impose the Dirichlet boundary condition $\theta = 1$ at the particle surface $r/a=1$ and the Neumann boundary condition on the outer surface $r/a = 100$ to approximate the zero-concentration boundary far from the particle. Time stepping is performed using an implicit Euler scheme with a time step ${{\rm \Delta} t}^{\star } = 0.0236$ from the initial condition $\theta ({\boldsymbol {y}}^{\star }, 0) = 0$.

As in Lawson & Ganapathisubramani (Reference Lawson and Ganapathisubramani2021), there is an initial transient in the ensemble-average mass transfer rate $\left \langle {Sh} \right \rangle$ obtained from the simulation, which corresponds to the diffusive growth of the concentration boundary layer. An example of this transient is provided in Appendix A, where a stationary state can be observed by ${t}^{\star } > 75$. Averaging over the time interval $75 < {t}^{\star } < {T}^{\star }_{sim}$, the deviation of the ensemble average $\left \langle Sh \right \rangle$ from the time average $\left \langle \overline {Sh} \right \rangle$ does not exceed $2\,\%$ and the statistical uncertainty in $\left \langle \overline {Sh} \right \rangle$ is below ${\pm }0.5\,\%$ (based on a 95 % confidence interval).

4.2. Average mass transfer rate

We first examine the average mass transfer rate predicted by our numerical simulations. Figure 7(a) shows the ensemble-average mass transfer rate to sinking particles of different sizes as a function of the turbulent Péclet number ${Pe}$ and sinking ratio ${Sr}$. We observe that the average mass transfer rate is insensitive to the particle size relative to the Kolmogorov scale $(a/\eta = 10^{-2}, 10^{-1})$, provided the sinking ratio ${Sr}$ is held constant. The square markers in figure 7(a) show experimental measurements of the mass transfer rate to small spheres in turbulence with ${St} < 0.1$ and ${Sr} < 0.12$ obtained from data tabulated in Armenante (Reference Armenante1983) and later reported in Armenante & Kirwan (Reference Armenante and Kirwan1989). At small ${Sr} = 0.1$, our numerical simulations are in good agreement with the experimental data. For fixed ${Pe}$, a significant enhancement in convective mass transfer is seen as ${Sr}$ is increased beyond ${Sr} \approx 0.33$.

Figure 7.

Ensemble-average Sherwood number as a function of (a) turbulent Péclet number ${Pe}$ and (b) quiescent settling Péclet number ${Pe}_w$. The markers correspond to: simulations with $a/\eta = 10^{-2}\ (\circ )$ and $10^{-1}\ (\times )$; experimental data of small spheres $(\square )$ with ${St} < 0.1$ and ${Sr} < 0.12$ from Armenante (Reference Armenante1983). The dashed line corresponds to the asymptotic expression (4.2) for small particles sinking in quiescent flow (Acrivos & Goddard Reference Acrivos and Goddard1965).

Figure 7(b) shows the same simulation data as a function of the quiescent settling Péclet number ${Pe}_w$. When convection due to sinking is strong in comparison with the turbulence (i.e. ${Sr} \gg 1$), the data approach the high-Péclet-number asymptotic expression

(4.2)\begin{equation} {Sh} = 0.6245 {Pe}_w^{1/3} + 0.461 + O({Pe}_w^{{-}1/3}) \end{equation}

for small particles in uniform flow (Acrivos & Goddard Reference Acrivos and Goddard1965). By comparing data at similar ${Pe}_w$ but different ${Sr}$, we note that the general trend is for strong turbulence (i.e. ${Sr} < 1$) to increase the mass transfer rate compared to sinking alone. However, at the very largest Péclet numbers ${Pe}_w \ge 10^4$, we observe a small but statistically significant decrease ($0.20 \pm 0.035$, or around $1.5\,\%$) in $\left\langle \overline{Sh} \right\rangle$ at $Sr = 1$ compared with $Sr = 10$. This is consistent with the assertion of Batchelor (Reference Batchelor1980) that, for fixed ${Pe}_w \gg 1$, the effect of turbulent stirring should initially decrease the average rate of mass transfer. Figures 7(a) and 7(b) therefore demonstrate that ${Sr}$ effectively interpolates between two regimes of mass transfer: strain-dominated, where ${Sr} \ll 1$ and the mass transfer rate primarily depends upon the turbulent strain rate; and slip-dominated, where ${Sr} \gg 1$ and convective mass transfer is determined by the slip velocity.

This observation raises the possibility of finding an empirical function which smoothly interpolates between these two limiting regimes. Motivated by the success of previous empirical expressions of the form (1.1) with $\beta =1/2$ and $\gamma =1/3$ (Frössling Reference Frössling1938; Ranz & Marshall Reference Ranz and Marshall1952), we plot the normalised mass transfer rate $(\left \langle {Sh} \right \rangle -1)/{Pe}^{1/3}$ as a function of the sinking ratio in figure 8. Although this normalisation does not incorporate the anomalous Péclet number scaling we have observed due to the convection-suppression effect (Lawson & Ganapathisubramani Reference Lawson and Ganapathisubramani2021), we find it provides a satisfactory collapse of our data over $10^2 \le {Pe} \le 10^4$. A fit to the data in this range of the form

(4.3)\begin{equation} {Sh} = 1 + {Pe}^{1/3}\frac{\left( \alpha' {Sr} + \beta' \right)^{4/3}}{{Sr} + \gamma'} \end{equation}

with $\alpha ' = 0.6245^{3/4}$ provides the asymptotic scaling (4.2) at large ${Sr}$ accurate to leading order in ${Pe}_w$. The coefficients $\beta '=1.293$ and $\gamma '=3.269$ are chosen by minimising the sum of squares error in $\left \langle \overline {{Sh}} \right \rangle$. The fit to our numerical data in this range is good, with the largest residual being $5.2\,\%$.

Figure 8.

Normalised mass transfer rate $\left \langle {Sh}-1 \right \rangle /{Pe}^{1/3}$ shown as a function of the sinking ratio ${Sr}$. Filled markers correspond to numerical simulation data at fixed ${Pe}$, whereas open markers correspond to experimental data at fixed ${Pe}_w$ digitally extracted from Ohashi et al. (Reference Ohashi, Sugawara and Kikuchi1981). The solid black line shows the best fit of our simulation data to (4.3) over the range $10^2 \le {Pe} \le 10^4$.

Compensating for Péclet number in this manner, we are also able to directly compare our numerical simulations with the experimental data of Ohashi et al. (Reference Ohashi, Sugawara and Kikuchi1981), who measured the mass transfer rate from freely suspended spherical particles sinking in vertical turbulent pipe flow. In these experiments, three sizes of spherical ion-exchange particles were injected into a vertical turbulent pipe flow containing a sodium hydroxide electrolyte. The mass transfer was determined by the diffusion-limited exchange of sodium/potassium cations, measured by atomic absorption analysis. The mass transfer coefficient $k = {Sh} D / a$ was determined at fixed ${Sc} = 368$ for varying bulk upflow velocity $U_b$, and therefore varying turbulent dissipation rate $\epsilon = 2f U_b^3/D_p$, where $f$ is the Fanning friction factor and $D_p$ is the pipe diameter. The quiescent sinking Reynolds and Péclet numbers span ${\it Re} _w = 5.5, 20, 42$ and ${Pe}_w = 1.0\times 10^3, 3.6\times 10^3, 7.7\times 10^3$. We digitally extracted the mass transfer coefficient $k$ and bulk flow velocity $U_b$ from their figure 1 and computed the quiescent sinking velocity of their ion-exchange particles using (3.2). Following Ohashi et al. (Reference Ohashi, Sugawara and Kikuchi1981), we obtained the friction factor from the Blasius correlation $f = 0.079 (U_b D_p/\nu )^{-1/4}$. The quantitative agreement at the smallest ${\it Re} _w = 5.5$ is good, despite this being beyond the Stokes regime of our simulation. The quantitative agreement at larger Reynolds numbers is fair, but shows the same trends. This is expected, given the reported Reynolds-number dependence of mass transfer in uniform flow and the well-established appearance of a wake recirculation zone at ${\it Re} _w \approx 20$ (Clift et al. Reference Clift, Grace and Weber1978). These observations provide further confidence in our numerical simulation approach and confirm the existence of a transition in the mass transfer mechanism parametrised by ${Sr}$.

4.3. Mechanism of mass transfer

Figure 9 presents a visualisation of the concentration boundary layer and wake around a simulated particle as the sinking ratio is increased. When ${Sr} = 0.1$, fluid is swept approximately along streamlines into two, nearly symmetric wakes. This closely resembles the topology of the concentration field expected under steady, extensional strain. As the sinking ratio is increased, the topology of the wake changes towards a single wake in a near-uniform flow field. This closely resembles the topology of the concentration field expected under steady, uniform flow driven by slip. We therefore qualitatively confirm the existence of a transition in the dominant mass transfer mechanism between the small and large ${Sr}$ regimes.

Figure 9.

Visualisation of the concentration boundary layer and wake around a spherical particle at ${Pe} = 10^{4}$ and $a/\eta = 10^{-2}$ with (a) ${Sr} = 0.1$, (b) ${Sr} = 1$ and (c) ${Sr} = 10$. The visualisation shows the concentration field around a particle at ${t}^{\star }=100$, which was seeded at the same position at ${t}^{\star }=0$ but subject to different sinking rates.

To examine the transition in mass transfer mechanism quantitatively, we consider a budget of the diffusive, convective and unsteady solute fluxes in a control volume surrounding the particle. The flow field relative to the particle (2.7) can be decomposed into linear contributions $\boldsymbol {v} = \boldsymbol {v}_\epsilon + \boldsymbol {v}_w$ from the relative slip due to gravity $\boldsymbol {v}_w$ and strain $\boldsymbol {v}_\epsilon$ imposed by turbulent eddies. For a spherical shell control volume $\boldsymbol{\mathcal{V}} = \{\boldsymbol {y} : a \le |\boldsymbol {y}| \le r \}$ with variable outer radius $r$, we obtain the temporally averaged balance of convective and diffusive solute fluxes from (2.8) as

(4.4)\begin{equation} {A}^{{\star}} + {Q}^{{\star}}_\epsilon + {Q}^{{\star}}_{w} + {Q}^{{\star}}_{D} - 4{\rm \pi}\overline{\left\langle {Sh} \right\rangle} = 0, \end{equation}

where we have defined

(4.5)\begin{equation} \left.\begin{gathered} {A}^{{\star}}({r}^{{\star}}) = \int_{\boldsymbol{\mathcal{V}}} \overline{\left\langle \frac{\partial \theta}{\partial {t}^{{\star}}} \right\rangle} \mathrm{d}\boldsymbol{y}, \\ {Q}^{{\star}}_{\epsilon}({r}^{{\star}}) = \int_{|{\boldsymbol{y}}^{{\star}}|={r}^{{\star}}} \overline{\left\langle {\boldsymbol{v}_\epsilon}^{{\star}}\theta \right\rangle}\boldsymbol{\cdot} \mathrm{d}\boldsymbol{S},\quad {Q}^{{\star}}_{w}({r}^{{\star}}) = \int_{|{\boldsymbol{y}}^{{\star}}|={r}^{{\star}}} \overline{\left\langle {\boldsymbol{v}_w}^{{\star}}\theta \right\rangle}\boldsymbol{\cdot} \mathrm{d}\boldsymbol{S}, \\ {Q}^{{\star}}_D({r}^{{\star}}) ={-}\frac{1}{{Pe}} \int_{|{\boldsymbol{y}}^{{\star}}|={r}^{{\star}}} \boldsymbol{\nabla}_\star\overline{\left\langle \theta \right\rangle} \boldsymbol{\cdot} \mathrm{d}\boldsymbol{S}, \end{gathered}\right\} \end{equation}

in conjunction with the boundary condition $\boldsymbol {v}=\boldsymbol {v}_w=\boldsymbol {v}_\epsilon = 0$ on $r=a$. Numerically, this is obtained through volume integrals of the convection–diffusion terms in (2.8), which is solved using a conservative, finite-volume method. Thus, mass conservation is implicitly satisfied. The terms in (4.4) can be understood as follows. The diffusive flux of material into $\boldsymbol {\mathcal {V}}$ from the particle surface is $4{\rm \pi} \overline {\left \langle {Sh} \right \rangle }$. This must be balanced by the convective flux due to slip ${Q}^{\star }_w$ and strain ${Q}^{\star }_\epsilon$, the diffusive flux ${Q}^{\star }_D$ at $|\boldsymbol {y}|=r$ or the time-average unsteady accumulation ${A}^{\star }$ of material within $\boldsymbol {\mathcal {V}}$. Because $\boldsymbol {v}_w(\boldsymbol {y},t)=\boldsymbol {v}_w(-\boldsymbol {y},t)$ is even in $\boldsymbol {y}$, only the odd, parity-antisymmetric component of $\theta (\boldsymbol {y},t)$ contributes to the convective flux due to slip (4.5). Likewise, since the strain component $\boldsymbol {v}_\epsilon (\boldsymbol {y},t)=-\boldsymbol {v}_\epsilon (-\boldsymbol {y},t)$ is odd in $\boldsymbol {y}$, only the parity-symmetric, even component of $\theta (\boldsymbol {y},t)$ contributes to the convective flux due to strain.

Figure 10(a–c) shows the balance of these terms (4.4) in the mass transport away from a spherical particle at ${Pe} = 10^4$ and ${Sr} = 0.32, 1.0$ and $3.2$. The sum of the convective terms ${Q}^{\star }_C = {Q}^{\star }_w + {Q}^{\star }_\epsilon$ is also shown. Near the particle, diffusion dominates mass transfer; far from the particle, the diffusive flux ${Q}^{\star }_D$ vanishes by construction at the boundary condition $r/a=100$. We therefore define a diffusion thickness, $\delta _{D}$, such that ${Q}^{\star }_{D}(a+\delta _{D}) = 0.01 (4{\rm \pi} \left \langle \overline {{Sh}} \right \rangle )$. This defines a measure of the thickness of the concentration boundary layer, at the edge of which the diffusive flux of solute away from the particle is negligible. We note that this is analogous to defining a $99\,\%$ thickness in terms of the scalar gradient $\boldsymbol {\nabla }_\star \overline {\left \langle \theta \right \rangle }$. Beyond this scale, convection dominates the mass transfer out of $\boldsymbol {\mathcal {V}}$ and the control volume is in equilibrium, i.e. the unsteady accumulation ${A}^{\star }$ is negligible. At very large $r/a \approx 58$, the equilibrium approximation breaks down and the unsteady accumulation of material accounts for over $10\,\%$ of the solute flux from the particle. This is expected given the finite run-time of our simulation. Nonetheless, we are able to observe a wide range of scales over which convection ${Q}^{\star }_C \approx 4{\rm \pi} \left \langle \overline {{Sh}} \right \rangle$ balances the diffusive flux from the particle surface.

Figure 10.

Decomposition of the scalar transport from a spherical particle as a function of control volume radius $r/a$ at $a/\eta = 10^{-2}$. (a–c) Full breakdown of unsteady, diffusive and convective terms in (4.4) for ${Pe} = 10^4$ and (a) ${Sr} = 0.32$, (b) ${Sr} = 1.0$ and (c) ${Sr} = 3.2$. The black dotted line indicates $r=a+\delta _D$, where ${Q}^{\star }_D(r)/4{\rm \pi} \left \langle \overline {Sh} \right \rangle = 0.01$. (d) Breakdown of convective terms only for ${Pe} = 10^1\unicode{x2013}10^4$ and ${Sr} = 1.0$. Markers indicate the turbulent Péclet number, whilst colour indicates terms as in (a–c).

By comparing figure 10(a–c), we discern the effect of the sinking ratio upon the mechanism of convective mass transfer away from the particle. At large turbulent Péclet number, we identify a large range of scales over which convection dominates the mass transfer away from the particle. The lower end of this range corresponds to the thickness of the concentration boundary layer, $\delta _D$, whilst the upper end of this range is limited by the time needed for our simulation to reach equilibrium. Across this range, the balance between the convective flux due to slip ${Q}^{\star }_w$ and strain ${Q}^{\star }_\epsilon$ is not uniform. In all cases, convection by strain dominates the solute flux away from the particle at large scales. However, there exists a subrange of scales over which the relative contributions of convection due to strain and slip are approximately constant. The lower extent of this range corresponds to the diffusive thickness of the concentration boundary layer, $\delta _D$, whereas the upper extent of this range increases as the sinking rate ${Sr}$ is increased.

The balance of convection and diffusion across scales also depends upon the turbulent Péclet number, as shown in figure 10(d). As the turbulent Péclet number is increased, the concentration boundary layer thins and we observe an increasing range of scales over which convection dominates the mass transfer away from the particle. The relative contributions of the convective fluxes exhibit a remarkable collapse beyond $r/a \approx 3$ for ${Sr} = 1$. Figure 11(a) shows the thickness of the concentration boundary layer in terms of the conventional 99 % thickness $\delta _{99}$ and the diffusion layer thickness $\delta _D$. We observe that the two definitions provide comparable measures of the boundary layer thickness and converge towards similar values at large ${Pe}$. The data are also in agreement with the expected $\delta _{99} \sim {Pe}^{-1/3}$ scaling at large ${Pe}$, although the range over which this scaling is applicable is arguably only valid for ${Pe} > 10^3$ such that $\delta _{99}/a \ll 1$.

Figure 11.

(a) The variation of the concentration boundary layer thickness based on the diffusive flux thickness $\delta _D$ and the $99\,\%$ thickness $\delta _{99}$ with turbulent Péclet number and sinking ratio. (b) The breakdown of fluxes due to strain (${Q}^{\star }_\epsilon$, dashed line), slip (${Q}^{\star }_w$, solid line) and unsteady accumulation (${A}^{\star }$, dotted line) at the edge of the concentration boundary layer $r = a + \delta _D$, as a function of the sinking ratio ${Sr}$. Colours correspond to varying turbulent Péclet number.

From our discussion so far, we have seen that the mechanism of convective solute transport away from the particle depends upon scale. We have also seen that the diffusion thickness $\delta _D$ serves as a measure of the thickness of the concentration boundary layer, which lies within a region where convection is dominant and the relative contributions due to slip and strain are approximately independent of $r$. It is therefore instructive to compare the relative contribution of these convective fluxes within this range. Figure 11(b) shows this balance at $r = a+\delta _D$. By definition, diffusion accounts for 1 % of the normalised flux at this location. As the turbulent Péclet number is increased, the relative contributions of each of these mechanisms approaches an asymptotic value. As the sinking ratio is increased, there exists a transition in the mechanism of convective mass transfer at the edge of the concentration boundary layer: at ${Sr} = 0.1$, virtually all mass transfer is attributable to the turbulence; at ${Sr} = 10$, virtually all mass transfer occurs due to slip. Around ${Sr} \approx 1$, the contributions are comparable. This quantitatively demonstrates a change in the mechanism of mass transfer at the scale of the particle from strain-dominated to slip-dominated as the sinking ratio ${Sr}$ is increased.

4.4. Fluctuations in the mass transfer rate

Thus far, we have examined only the average rate of mass transfer from the particle. We now turn our attention to the fluctuations in the mass transfer rate and their correlation with the convective forcing. Figure 12 shows the magnitude of the mass transfer fluctuations with respect to the mean as a function of the turbulent Péclet number and sinking ratio. At small sinking ratio ${Sr} \ll 1$, we observe first an increase in the magnitude of the mass transfer fluctuations with increasing ${Pe}$, with a shallow peak around ${Pe} \approx 10^2\unicode{x2013}10^3$, followed by a decrease in fluctuation magnitude with increasing ${Pe}$. This is consistent with the behaviour we have observed for neutrally buoyant spherical and spheroidal particles, which arises from the convection-suppression effect (Lawson & Ganapathisubramani Reference Lawson and Ganapathisubramani2021). The relative magnitude of ${Sh}$ fluctuations is slightly larger in the weakly sinking case ${Sr} = 0.1$ compared with our previous data at ${Sr} = 0$; the difference is larger than the statistical uncertainty in the data. As the sinking ratio is further increased, the relative magnitude of the mass transfer fluctuations is decreased. Furthermore, the ${Pe}$ dependence is also changed. At the largest sinking rates ${Sr} = 3.2, 10$ we observe a monotonic increase in the fluctuation magnitude with increasing ${Pe}$. However, the magnitude of these fluctuations remains quite small (2 %–5 % of the mean) at the largest turbulent Péclet numbers tested.

Figure 12.

Standard deviation $\sigma _{{Sh}}$ of fluctuations in the local mass transfer rate ${Sh}$, relative to the mean, as a function of the turbulent Péclet number and sinking ratio ${Sr}$. The open markers correspond to simulations with $a/\eta = 10^{-2}\ (\circ )$ and $10^{-1}\ (\times )$; the black filled markers show data for spherical particles at ${Sr} = 0$ from Lawson & Ganapathisubramani (Reference Lawson and Ganapathisubramani2021).

For sinking particles, turbulence introduces fluctuations in the solution of (2.8) through two mechanisms. Firstly, fluctuations in strain $\boldsymbol{\mathsf{E}}$ may directly cause fluctuations in the mass transfer rate. Secondly, although the particle sinks at a constant rate $\boldsymbol {w}_0$ in the laboratory frame, turbulence causes the particle to rotate with angular velocity ${\boldsymbol {\varOmega } = \boldsymbol {\omega }/2}$. Therefore, the slip velocity in the particle frame $\boldsymbol {v}_0$ also fluctuates in time, since $\dot {\boldsymbol {v}_0} = -w_0 \boldsymbol{\mathsf{R}}^{{\rm T}} (\boldsymbol {\varOmega } \times \boldsymbol {e}_z)$. The magnitude of these fluctuations is therefore determined by the rate of rotation $\varOmega _{\perp } = |\boldsymbol {\varOmega } \times \boldsymbol {e}_z|$ perpendicular to the direction of gravity. The two mechanisms may have significant correlation with each other, since the strain rate $E$ and vorticity magnitude $|\boldsymbol{\omega} |$ are strongly correlated in homogeneous turbulence (Zeff et al. Reference Zeff, Lanterman, McAllister, Roy, Kostelich and Lathrop2003). Indeed, in our data, we find that $E$ and $\varOmega _\perp$ are interdependent with product-moment correlation coefficient $R_{E\varOmega } \approx 0.54$.

In figure 13, we present the cross-correlation functions $R_{E}(\tau )$ and $R_{\varOmega }(\tau )$ between the instantaneous mass transfer rate ${Sh}(t+\tau )$ evaluated at a time lag $\tau$ relative to the strain $E(t) = |\boldsymbol{\mathsf{E}}|$ (figure 13a–c) and rotation rate $\varOmega _{\perp }(t)$ (figure 13d–f). For ${Sr} = 0.1$, where the average mass transfer rate is dominated by strain, we observe that the strain $E$ is strongly, positively correlated with the instantaneous transfer rate. There is also a weaker, positive correlation with the rotation rate $\varOmega _\perp$. As the turbulent Péclet number is increased, the local transfer rate correlates less strongly with the local strain rate (and rotation rate) and over a longer time scale. A similar behaviour can also be observed at ${Sr} = 1$. This is an example of the convection-suppression effect, where the concentration boundary layer becomes unresponsive to velocity fluctuations occurring on time scales faster than $O(\tau _\eta {Pe}^{1/3})$, in agreement with simulations of neutrally buoyant spheroids (Lawson & Ganapathisubramani Reference Lawson and Ganapathisubramani2021).

Figure 13.

Cross-correlation function of the instantaneous mass transfer rate ${Sh}(t+\tau )$ against (a–c) local strain magnitude $E(t)$ and (d–f) local rotation rate $\varOmega _\perp (t)$ at time $t$. Line colours correspond to different turbulent Péclet numbers. The circular markers indicate the time lag $\tau$ at which the cross-correlation is maximised.

On the other hand at ${Sr} = 10$, where the average mass transfer rate is dominated by slip, fluctuations in the mass transfer rate are strongly, negatively correlated with the rotation rate $\varOmega _\perp$ and to a lesser extent the strain rate $E$. Furthermore, increasing the turbulent Péclet number increases the strength of this effect and also increases the time scale over which the mass transfer rate correlates with the rotation rate. We therefore conclude that particle rotation by turbulence acts to decrease the mass transfer rate by slip. This is consistent with Batchelor's argument that, in the limit of high turbulent Péclet number, ‘the rotation of the particle and the ambient fluid has the effect of suppressing the convective transport due to … streaming motions’ (Batchelor Reference Batchelor1980). However, as we observe in figures 7(b) and 12, the net effect upon the mean and fluctuations of the mass transfer rate is small when ${Sr} \gg 1$.

Lawson & Ganapathisubramani (Reference Lawson and Ganapathisubramani2021) showed that, for neutrally buoyant spheroids, the instantaneous mass transfer rate is most strongly correlated with the velocity field $\tilde{\boldsymbol{v}}$ (2.13) filtered on a finite time scale $\tau _f \sim \tau _\eta {Pe}^{1/3}$ due to the convection-suppression effect. To test this hypothesis for sinking particles, we introduce the filtered, Lagrangian forcing terms

(4.6a,b)\begin{equation} \tilde{\boldsymbol{\mathsf{E}}}(t) = \frac{1}{\tau_f} \int_{-\tau_f}^{0} \boldsymbol{\mathsf{E}}(t+\tau)\,\mathrm{d}\tau \quad\textrm{and}\quad \widetilde{\boldsymbol{v}_0}(t) = \frac{1}{\tau_f} \int_{-\tau_f}^{0} \boldsymbol{v}_0(t+\tau)\,\mathrm{d}\tau , \end{equation}

which specify the filtered relative velocity field in (2.13). We note that $\boldsymbol {v}_0$ only fluctuates in the particle-aligned frame and is a constant in the laboratory frame, so fluctuations in $\boldsymbol {v}_0$ are induced by turbulent rotation of the particle. We then measure the multiple correlation coefficient

(4.7)\begin{equation} R^2 = \frac{r^2_{{Sh},\tilde{E}} + r^2_{{Sh},\tilde{v}_0} - 2 r_{{Sh},\tilde{E}} r_{{Sh},\tilde{v}_0} r_{\tilde{E},\tilde{v}_0}}{1 - r^2_{\tilde{E},\tilde{v}_0}} \end{equation}

between the instantaneous mass transfer rate ${Sh}$ and the magnitudes of the filtered strain $\tilde {E} = |\tilde {\boldsymbol{\mathsf{E}}}|$ and relative slip velocity $\tilde {v}_0 = |\tilde {\boldsymbol {v}}_0|$ as a function of the filter time scale. In (4.7), $r_{a,b}$ are the ordinary product-moment correlation coefficients of the variables indicated in the subscripts. Essentially, this measures the extent to which the linear model $\left \langle Sh \right \rangle + a_1 (\tilde {v}_0-\left \langle \tilde {v}_0 \right \rangle ) + a_2 (\tilde {E}-\langle {\tilde {E}}\rangle )$ predicts variations in ${Sh}$ due to variations in $\tilde {E}$ and $\tilde {v}_0$. The result is shown in figure 14. For ${Sr} = 0.10$, we observe a strong correlation $r_{Sh,\tilde {E}}$ between the filtered strain $\tilde {E}$ and the instantaneous mass transfer rate, similar to neutrally buoyant particles. As the turbulent Péclet number is increased, the time scale $\tau _f$ which maximises the strength of this correlation also increases. The additional information from $\tilde {v}_0$ only modestly improves the ability to predict fluctuations in ${Sh}$. However, as the sinking ratio is increased, the strength of the correlation coefficient $r_{Sh,\tilde {E}}$ decreases and fluctuations in $\tilde {E}$ become anti-correlated with ${Sh}$ at ${Sr} = 10$, which is consistent with our observations in figure 13 that turbulent rotation suppresses convection. In contrast, a much stronger joint correlation is observed when both the filtered strain $\tilde {E}$ and relative velocity $\tilde {v}_0$ are used to predict ${Sh}$. The general trend is for the peak correlation to become weaker as the sinking ratio is increased and as the turbulent Péclet number is increased. However, an exception to this trend is seen at ${Sr} = 10$, where the correlation strength increases with ${Pe}$.

Figure 14.

Joint cross-correlation coefficient of the filtered velocity field with the instantaneous transfer rate at sinking ratios of (a) ${Sr} = 0.1$, (b) ${Sr} = 1$ and (c) ${Sr} = 10$. The solid lines show the multiple correlation coefficient $R$ (4.7) for varying filter time scale $\tau _f$. The dashed lines show the ordinary correlation coefficient $r_{{Sh},\tilde {E}}$ between the instantaneous transfer rate and the filtered strain magnitude $\tilde {E}$.

In figure 15, we identify the time scale $\tau _d = \arg \max _{\tau _f} R$ which maximises the multiple correlation coefficient $R$. We interpret this time scale as the response time of the concentration boundary layer to velocity fluctuations. In figure 15(a), we observe that the response time increases with increasing ${Pe}$ approximately as ${Pe}^{1/3}$ and decreases with increasing ${Sr}$. This scaling is confirmed in figure 15(b), where we observe that the data collapse well when $\tau _d$ is normalised by $\tau _\eta {Pe}^{1/3}$. Therefore, our data are consistent with Batchelor's prediction of the convection-suppression effect in turbulent flows, provided that the time scale for averaging the relative velocity field is chosen in proportion to $\tau _\eta {Pe}^{1/3}$

Figure 15.

Identification of the response time scale $\tau _d$ of the concentration boundary layer to the fluctuating velocity field. (a) Response time scale $\tau _d$ shown as a function of ${Pe}$ for varying ${Sr}$. The black dashed line shows the scaling $\tau _d \sim {Pe}^{1/3}$. (b) Response time scale $\tau _d$ shown as a function of sinking ratio ${Sr}$, compensated for the expected $\tau _d \sim {Pe}^{1/3}$ scaling.