1 Introduction
Numerical simulations of multiphase systems in which small particles react with the embedding turbulent fluid are very complex. However, such systems are frequently encountered, especially in industry. Two excellent examples of this are pulverised coal combustion and fluidised bed combustion, which are two of the most popular technologies for thermal power generation. Given their popularity, it is not surprising that a lot of research is devoted to optimise and improve these processes. Representative examples of recent numerical studies are the works of Choi & Kim (Reference Choi and Kim2009) aimed at reduction of NO
$_{x}$
emission from pulverised coal combustion, Al-Abbas, Naser & Dodds (Reference Al-Abbas, Naser and Dodds2012) who focused on oxy-fuel combustion of low-rank coal, Gubba et al. (Reference Gubba, Ingham, Larsen, Ma, Pourkashanian, Tan, Williams and Zhou2012) who investigated pulverised coal and biomass co-firing and Adamczyk et al. who examined different combustion technologies in fluidised bed (Adamczyk et al.
Reference Adamczyk, Kozołub, Klimanek, Białecki, Andrzejczyk and Klajny2015) and pulverised coal (Adamczyk et al.
Reference Adamczyk, Bialecki, Ditaranto, Gladysz, Haugen, Katelbach-Wozniak, Klimanek, Sladek, Szlek and Wecel2017) boilers. In all of the aforementioned publications, the focus is on industrial-scale simulations, which would not be feasible without the aid of modelling. Heat and momentum exchange between phases as well as all the stages of coal conversion, such as drying, devolatilisation and char burnout, are usually taken into account. All of these processes can be strongly influenced by turbulence and these interactions have to be modelled. Thus, it is apparent that the accuracy of such simulations depends on the completeness of the models and their capability of reflecting real-life phenomena.
On the other hand, the amount of modelling can be greatly reduced by using an extremely fine computational grid such that the entire range of turbulent scales, together with the associated effects on mass, momentum and heat transfer, are resolved. Such an approach, called a direct numerical simulation (DNS), has lately been applied in studies of pulverised coal combustion by Luo et al. (Reference Luo, Wang, Fan and Yi2012, Reference Luo, Bai, Jin, Qiu and Fan2017), Brosh & Chakraborty (Reference Brosh and Chakraborty2014), Brosh et al. (Reference Brosh, Patel, Wacks and Chakraborty2015), Hara et al. (Reference Hara, Muto, Kitano, Kurose and Komori2015) and Muto, Yuasa & Kurose (Reference Muto, Yuasa and Kurose2017), allowing for very detailed examination of coal ignition, flame stabilisation or interactions between vortices and coal particles. However, the DNS technique is restricted in its use to small-scale simulations (or very low Reynolds numbers) since its application to industrial-scale problems is prohibitively expensive in terms of computational power. It should also be mentioned that the smallest fluid scale that is resolved in the simulations presented in the literature cited above is the Kolmogorov scale. The Kolmogorov scale is defined as the scale of the smallest turbulent eddies, i.e. the scale where the kinetic energy is dissipated into heat. The boundary layer of the particles, which are smaller than the Kolmogorov size, is, however, not resolved. Hence, the boundary layer effects have to be modelled.
Over the years, extensive investigations of turbulent flows have been conducted, resulting in an abundance of models accounting for the effect of turbulence on flow transport properties. Amongst the well-established ones are the
$k{-}\unicode[STIX]{x1D716}$
(Launder & Spalding Reference Launder and Spalding1974) and
$k{-}\unicode[STIX]{x1D714}$
(Wilcox Reference Wilcox1988) models, which relate the Reynolds stresses to the mean flow stresses through the eddy viscosity. There are also a number of models suitable for homogeneously reacting turbulent flows, such as the eddy dissipation model (Magnussen & Hjertager Reference Magnussen and Hjertager1977), probabilistic descriptions relating instantaneous scalar fluctuations to their mean values (Pope Reference Pope1985) or models for computing the turbulent flame speed (Zimont et al.
Reference Zimont, Polifke, Bettelini and Weisenstein1998). In addition to this, there is much more modelling required when dealing with reacting multiphase flows; for example, the particle dispersion caused by the turbulence in particle-laden flows can be described using stochastic transport models (Baxter & Smith Reference Baxter and Smith1993; Graham Reference Graham1996).
The above models are mature and well-tested, and most of them have numerous extensions. However, until very recently, no models existed that represent the impact of turbulence on heterogeneous reactions, such as char conversion. Therefore, in all industrial-scale simulations of pulverised coal combustion this effect is not captured. This knowledge gap was addressed by Krüger et al. (Reference Krüger, Haugen, Mitra and Løvås2017a
) who proposed the first model to account for the mass transfer rate modifications due to turbulence for isothermal reactions. This work was later extended to non-isothermal reactions (Krüger, Haugen & Løvås Reference Krüger, Haugen and Løvås2017b
). The numerical studies of Krüger et al. showed that under some circumstances, turbulence may cause the particles to form clusters, inside which the reactant is quickly consumed. As a result, there exist regions that are rich in reactants but lack particles, and regions where the particle density is very high but the reactant concentration is low. This leads to a considerable decrease in the overall mass transfer rate, which should not be neglected when performing simulations of turbulent flows containing reacting particles. The work of Krüger et al. (Reference Krüger, Haugen, Mitra and Løvås2017a
) was further extended by Haugen et al. (Reference Haugen, Krüger, Mitra and Løvås2018) who showed that turbulence may also enhance the mass transfer since it causes an increase in the relative velocity between phases. They proposed a model that accounts for both effects of turbulence on the reactant consumption rate. It was shown that which effect dominates depends on the Damköhler number, which is given by
$\mathit{Da}=\unicode[STIX]{x1D70F}_{L}/\unicode[STIX]{x1D70F}_{c}$
, where
$\unicode[STIX]{x1D70F}_{L}$
and
$\unicode[STIX]{x1D70F}_{c}$
are typical fluid (turbulence) and chemical time scales, respectively. For the remainder of this paper,
$\unicode[STIX]{x1D70F}_{L}$
is set to be equal to the turnover time of the integral-scale eddies. In these studies, monodisperse particle systems were studied, but in real situations, such as in pulverised coal-fired boilers, particles with different sizes are present. The current work aims to verify the previous findings of Krüger et al. (Reference Krüger, Haugen, Mitra and Løvås2017a
) and Haugen et al. (Reference Haugen, Krüger, Mitra and Løvås2018) for polydisperse systems and to further investigate how the mass transfer is affected by the turbulence under different conditions.
Clustering of particles due to turbulence in an embedding fluid has been studied in a large number of papers (Eaton & Fessler Reference Eaton and Fessler1994; Bec et al.
Reference Bec, Biferale, Cencini, Lanotte, Musacchio and Toschi2007; Calzavarini et al.
Reference Calzavarini, Kerscher, Lohse and Toschi2008; Toschi & Bodenschatz Reference Toschi and Bodenschatz2009; Baker et al.
Reference Baker, Frankel, Mani and Coletti2017; Haugen et al.
Reference Haugen, Krüger, Mitra and Løvås2018). It has been shown that for turbulence with a sufficient scale separation, clustering will occur for a range of Stokes numbers (Baker et al.
Reference Baker, Frankel, Mani and Coletti2017; Haugen et al.
Reference Haugen, Krüger, Mitra and Løvås2018). Conceptually, one can argue that particle clusters should be found at least for
$\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}\leqslant \unicode[STIX]{x1D70F}_{p}\leqslant \unicode[STIX]{x1D70F}_{L}$
, where
$\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$
is the Kolmogorov time scale and
$\unicode[STIX]{x1D70F}_{p}$
is the response time of the particles. Or, in other words, particle clusters should be found when
$\mathit{St}_{\unicode[STIX]{x1D702}}=\unicode[STIX]{x1D70F}_{p}/\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}>1$
and
$\mathit{St}_{L}=\unicode[STIX]{x1D70F}_{p}/\unicode[STIX]{x1D70F}_{L}<1$
. With this being said, it is clear that the sharpest clusters are found for
$\mathit{St}_{\unicode[STIX]{x1D702}}\approx 1$
(Bec et al.
Reference Bec, Biferale, Cencini, Lanotte, Musacchio and Toschi2007; Calzavarini et al.
Reference Calzavarini, Kerscher, Lohse and Toschi2008; Toschi & Bodenschatz Reference Toschi and Bodenschatz2009). Even though the sharpest clusters are found at scales around the Kolmogorov scale, this does not necessarily mean that these clusters have the most influence on the mass transfer coefficients. The reason for this is that in order for the clustering to slow down the mass transfer, the life time of the clusters should be longer than the time it takes for the particles to consume most of the reactants within the cluster. This is typically not the case for these smallest clusters since their life times are too short.
2 Governing equations and numerical methods
As mentioned before, there is a large number of applications where heavy inertial particles are embedded in a turbulent flow. For many of these applications, there is also mass transfer between the particles and the fluid. The mass transfer may yield a net molar production (e.g. char oxidation to carbon monoxide or the reduction of a metal oxide by natural gas to produce steam and carbon dioxide), a net molar reduction (e.g. oxidation of a metal oxide) or it may be molar neutral (e.g. oxidation of char to carbon dioxide). On top of this, the reactions may be exothermic (e.g. oxidation of char or a metal oxide) or endothermic (e.g. gasification of char). The mass transfer rate may also be dependent on temperature, through the kinetic reaction rates at the particle surface. An example of the complexity can be seen in the detailed conversion model for single point particle char as described in, for example, Haugen, Tilghman & Mitchell (Reference Haugen, Tilghman and Mitchell2014) and Haugen, Mitchell & Tilghman (Reference Haugen, Mitchell and Tilghman2015). In order to make the results obtained in the following as general as possible, and to be able to isolate the effects of the turbulence alone, we will here use a simplified description of the chemical reactions. The main simplifications made in this work are that (1) heterogeneous reaction kinetics is assumed to be infinitely fast, (2) the reaction occurs only on the external surface of the spherical particles, (3) the reaction is of one step, unimolar and isothermal and (4) the evolution of particle size and density is not accounted for, i.e. the particle acts as a catalyst.
2.1 Fluid-phase equations
Fluid motion is governed by the continuity

and the momentum equation

where
$\text{D}/\text{D}t=\unicode[STIX]{x2202}/\unicode[STIX]{x2202}t+\boldsymbol{u}\boldsymbol{\cdot }\unicode[STIX]{x1D735}$
is the advective derivative,
$\unicode[STIX]{x1D70C}$
is the density of the fluid,
$\boldsymbol{u}$
is the velocity vector and
$\text{}\underline{\text{}\underline{\unicode[STIX]{x1D70E}}}$
is the total stress tensor, which is defined as

when
$p$
and
$\unicode[STIX]{x1D707}$
are pressure and dynamic viscosity, respectively. The rate of strain tensor is given by

The remaining terms in (2.2) are
$\boldsymbol{f}$
and
$\boldsymbol{F}$
. The former is a forcing function, which is responsible for injecting kinetic energy into the domain (see Haugen et al. (Reference Haugen, Kleeorin, Rogachevskii and Brandenburg2012) for details), while the latter represents momentum exchange between the fluid and particles and will be introduced in the next section.
Thermodynamic properties are coupled through the ideal gas equation:

in which
$c_{s}$
is the isothermal speed of sound. Additionally, the equation governing the reactant is given as

where
$X$
stands for the reactant mole fraction,
$D$
is its diffusivity and
$R$
represents the rate of reaction occurring on the particle surface.
2.2 Dispersed-phase equations
Particles are represented as point sources and are tracked in a Lagrangian reference frame. For the point particle approach to be applicable, the particle size cannot be greater than the size of a fluid grid cell. Furthermore, a two-way coupling is applied, i.e. there are mutual interactions between particles and the fluid phase. Since the particle distribution is assumed to be dilute, particle collisions are neglected. The
$i$
th particle obeys the following equations of motion:


where
$m_{i}$
denotes its mass while
$\boldsymbol{x}_{i}$
and
$\boldsymbol{v}_{i}$
are its position and velocity vectors, respectively. Assuming that no forces other than drag act on the particle, its acceleration can be expressed as

where the particle response time is given by

and the correction factor to account for the Reynolds number effect, as given by the empirical correlation of Schiller & Naumann (Reference Schiller and Naumann1933), is
$f_{i}=0.15Re_{i}^{0.687}$
. Here,
$Re_{i}=\unicode[STIX]{x1D70C}|\boldsymbol{u}(\boldsymbol{x}_{i})-\boldsymbol{v}_{i}|d_{p,i}/\unicode[STIX]{x1D707}$
is the particle Reynolds number. The Stokesian particle response time is given by

where
$d_{p,i}$
is the particle diameter and the particle material density is
$\unicode[STIX]{x1D70C}_{p,i}$
. Momentum exchange is incorporated in the fluid-phase equations through the last term on the right-hand side of (2.2):

in which
$V_{c}$
is the volume of the relevant grid cell and the summation is over all particles located inside the same grid cell. Similarly, the term representing the rate of reaction in (2.6) is given as

where
$X$
is the reactant mole fraction in the cell in which the particle is located,
$A_{i}=\unicode[STIX]{x03C0}d_{p,i}^{2}$
is the particle surface area and

stands for the mass transfer coefficient, while

is the Sherwood number, as given by the empirical correlation of Ranz & Marshall (Reference Ranz and Marshall1952), and
$\mathit{Sc}=\unicode[STIX]{x1D707}/(D\unicode[STIX]{x1D70C})$
is the Schmidt number. In obtaining (2.13) it has been assumed that the particle reactions are diffusion controlled, i.e. that a reactant will react instantly when it reaches the particle surface.
2.3 Characteristic scales and non-dimensional numbers
Before presenting the results, it is worthwhile introducing some relevant dimensionless numbers. There are two characteristic time scales for turbulent flows: the turnover time of the integral scale

and the time scale of the Kolmogorov scale eddies at which energy is dissipated

Here,
$u_{rms}$
is the root-mean-square velocity and the integral length scale of turbulence is given by (Pope Reference Pope2000)

where
$\unicode[STIX]{x1D708}=\unicode[STIX]{x1D707}/\unicode[STIX]{x1D70C}$
represents kinematic viscosity and
$\unicode[STIX]{x1D716}$
is the dissipation rate of turbulent kinetic energy. Note that the above definition of the integral scale is different from the one employed by Krüger et al. (Reference Krüger, Haugen, Mitra and Løvås2017a
) and Haugen et al. (Reference Haugen, Krüger, Mitra and Løvås2018). In their work, a turbulent forcing scale
$L_{f}=L_{x}/(2\unicode[STIX]{x03C0}k_{f})$
, where
$L_{x}$
is the domain size and
$k_{f}$
is the forcing wavenumber, was used as the integral scale. The two scales,
$L$
and
$L_{f}$
, are not the same. Indeed,
$L\approx 2\unicode[STIX]{x03C0}L_{f}$
, which means that also several non-dimensional numbers, such as the Damköhler number and the Stokes number based on the integral scale, are roughly a factor
$2\unicode[STIX]{x03C0}$
different from the ones used by Krüger et al. (Reference Krüger, Haugen, Mitra and Løvås2017a
) and Haugen et al. (Reference Haugen, Krüger, Mitra and Løvås2018). The reason for defining the integral scale differently from what was done in Krüger et al. (Reference Krüger, Haugen, Mitra and Løvås2017a
) and Haugen et al. (Reference Haugen, Krüger, Mitra and Løvås2018) was that
$L_{f}$
was not really an integral scale but merely a scale that was used to define a suitable Reynolds number. It can be seen from comparing the integral wavenumber,
$k_{L}=2\unicode[STIX]{x03C0}/L$
, with the energy spectrum, which is what is done in figure 1, that the current integral scale does indeed represent a size that is reasonably close to the real integral scale. In the figure, also the Taylor microscale and the Kolmogorov scales are shown, for both
$\mathit{Re}_{\unicode[STIX]{x1D706}}=150$
(red) and
$\mathit{Re}_{\unicode[STIX]{x1D706}}=350$
(black).

Figure 1. Kinetic energy spectrum; black:
$Re_{\unicode[STIX]{x1D706}}\approx 350$
, red:
$Re_{\unicode[STIX]{x1D706}}\approx 150$
(cases A and M in table 1). The wavenumbers
$k_{L}$
,
$k_{\unicode[STIX]{x1D706}}$
and
$k_{\unicode[STIX]{x1D702}}$
correspond to the integral, Taylor and Kolmogorov scales, respectively.
Table 1. Simulation parameters. The total number of grid points varies between
$64^{3}$
and
$512^{3}$
, depending on
$Re_{\unicode[STIX]{x1D706}}$
and
$Da$
. In all cases, the Schmidt number is equal to 0.2. The parameter
$Da_{2}$
quantifies the effect of particle clustering and is explained in a later sub-section.

Under the assumption of Stokes flow, the particle time scale (
$\unicode[STIX]{x1D70F}_{p}$
) is given by (2.11). Yet another scale, characteristic for reacting flows, is the chemical time scale. The reactant consumption rate for homogeneous reactant and particle distributions can be defined as

where overline means volume averaging and (2.13) and (2.14) have been employed. If we assume that the relative velocity between the particles and the fluid is negligible, such that
$\mathit{Sh}=2$
(see (2.15)), the reactant consumption rate then becomes

Here,
$n_{p}$
is the particle number density and
$\overline{d_{p}}$
is the average particle diameter. The chemical time scale is now defined as

By combining the above time scales, one can define several dimensionless numbers, such as the Damköhler number

the Stokes number based on the integral scale

and the Kolmogorov-based Stokes number

Naturally, in a polydisperse system, a range of Stokes numbers can be identified since the Stokes number is a function of the particle diameter. However, to make the analysis clear, the mean Stokes number will in the following be calculated based on a mean particle diameter, such that

when

While it is generally accepted that the particle clustering is strongest when
$St_{\unicode[STIX]{x1D702}}\approx 1$
, we chose to present the results in terms of
$St_{L}$
. This is, as already discussed in the introduction, because we expect these larger-scale clusters (i.e. clusters formed by the larger eddies) to have a greater influence on the mass transfer rate due to their longer life times.
2.4 Reactant decay rate
It has been confirmed by numerous numerical and experimental studies that particles embedded in a turbulent flow have a tendency to cluster (Eaton & Fessler Reference Eaton and Fessler1994; Yoshimoto & Goto Reference Yoshimoto and Goto2007; Monchaux, Bourgoin & Cartellier Reference Monchaux, Bourgoin and Cartellier2012). The formation of such clusters is primarily due to the turbulent eddies that have time scales that are comparable to the particle response time. From this, it follows that also the cluster size depends on the Stokes number, as shown by Haugen et al. (Reference Haugen, Krüger, Mitra and Løvås2018). For very low
$St$
, mainly the smallest eddies lead to clustering, which results in many small clusters present in the flow. When the Stokes number is close to or above unity, clustering is caused by the integral scale eddies, and thus the clusters are larger.
As Krüger et al. (Reference Krüger, Haugen, Mitra and Løvås2017a
) demonstrated, the extent to which particle clustering affects the rate of mass transfer depends on the Damköhler number. If
$\unicode[STIX]{x1D70F}_{c}>\unicode[STIX]{x1D70F}_{L}$
, i.e.
$Da$
is low, the amount of reactant inside the particle clusters is similar to what it is outside of the clusters; thus, the effect of the clustering is negligible and the reactant consumption rate is similar to what is found for homogeneous particle and reactant distributions, and as such, can be computed from (2.19). However, if
$\unicode[STIX]{x1D70F}_{c}<\unicode[STIX]{x1D70F}_{L}$
, which corresponds to high
$Da$
, a significant fraction of the reactant inside the clusters is consumed during the life time of the cluster. In this situation, clustering plays an important role as it leads to a decrease of the rate at which the reactant is consumed because the average particle sees a reactant concentration that is significantly lower than the average reactant concentration in the fluid. Krüger et al. (Reference Krüger, Haugen, Mitra and Løvås2017a
) proposed to incorporate this effect in a formulation of the reactant decay rate, such that the real decay rate is given by

where
$\unicode[STIX]{x1D6FC}_{c}$
may be interpreted as a cluster characteristic decay rate. This parameter is generally dependent on the cluster properties, such as size, shape and number density.
While clustering of particles will decrease the mass transfer rate between particles and fluid, there is another turbulence-induced effect that will enhance the mass transfer rate. This effect is due to the fact that particles will not be significantly accelerated by turbulent eddies that have time scales that are shorter than the response time of the particles. This means that such fast turbulent eddies will yield a relative velocity difference between the particles and the fluid. Haugen et al. (Reference Haugen, Krüger, Mitra and Løvås2018) demonstrated the effect of this and described how it can be modelled (see their equation (3.14)).
Finally, Haugen et al. (Reference Haugen, Krüger, Mitra and Løvås2018) showed that the combined effect of turbulence on the mass transfer rate between particles and fluid can be expressed solely through a modified Sherwood number, which is given by

Using the above while employing (2.19)–(2.22) together with (2.27) yields

In the following, the predictions of (2.29) will be compared with the time-averaged reactant decay rate obtained from the DNS

where
$\langle \cdot \rangle$
represents time averaging,

is the volume-averaged decay rate and the integration starts at
$t_{0}$
, which is when the statistically stationary state is reached and
$\overline{X}$
is initialised to unity.
2.5 Computational methods and simulation set-up
A cubic computational domain with sides of length
$2\unicode[STIX]{x03C0}$
and periodic conditions prescribed for all boundaries is considered. Initially, particles are randomly distributed in a turbulent flow field. This is achieved by performing simulations only of the fluid phase until the turbulence is statistically stationary before inserting the particles. Also, the turbulent field is statistically invariant under translations and rotations, i.e. it is homogeneous and isotropic.
The modelling approach employed in the current work is almost identical to the one used by Krüger et al. (Reference Krüger, Haugen, Mitra and Løvås2017a ) and Haugen et al. (Reference Haugen, Krüger, Mitra and Løvås2018) with the difference that now a polydisperse particle system is considered. The key points of the approach are as follows. The DNS technique is applied, which means that the flow is resolved on the numerical grid down to the smallest turbulent scales. Therefore, since all the scales of the turbulence are resolved on the numerical mesh, only the interactions between phases (such as mass and momentum exchange between fluid and particles) require modelling due to the point-particle approximation. The accuracy of results is ensured by using high-order numerical methods, i.e. a third-order Runge–Kutta scheme for time advancement and a sixth-order finite-difference scheme for spatial derivatives, as implemented in the Pencil Code (Pencil Code), which is the software used for all simulations.
Various particle size distributions are studied, where the upper and lower particle diameter cut-offs are given by
$d_{max}$
and
$d_{min}$
, respectively. The first distribution is referred to as the ‘uniform distribution’ and corresponds to the situation when the probability of picking a random particle with a given diameter is equal for all diameters between
$d_{max}$
and
$d_{min}$
. Since the Damköhler number associated with a given particle size is proportional to the particle diameter, it follows that for this particle size distribution, the contribution to the total Damköhler number is larger for the larger particles. To see the effect when the contribution to the total Damköhler number is the same for all particle sizes, the particles are distributed such that the particle number density is given by
$n_{p}(d_{p})\sim 1/d_{p}$
, which is referred to as the ‘compensated distribution’. Finally, a Gaussian distribution (with a standard deviation
$\unicode[STIX]{x1D70E}=1/12(d_{max}-d_{min})$
) is considered. The above-mentioned particle size distributions are shown in figure 2. A Dirac distribution is also included, which corresponds to a single particle size, as used by Haugen et al. (Reference Haugen, Krüger, Mitra and Løvås2018) for simulations of monodisperse systems.

Figure 2. Particle size distributions used in the current work (for
$Re_{\unicode[STIX]{x1D706}}=150$
,
$\overline{\mathit{St}}_{L}=0.1$
).
Although we characterise the polydisperse system using the Stokes number defined based on the average-sized particle, it is important to remember that a broad range of Stokes numbers is contained in such a system. Figure 2 shows the normalised particle number density as a function of the Stokes numbers computed based on the integral (
$St_{L}$
, middle axis) and the Kolmogorov (
$St_{\unicode[STIX]{x1D702}}$
, top axis) time scales. (The values of Stokes numbers on the axes are relevant to the case in which the compensated distribution was used and the average Stokes number was approximately equal to 0.1.) It can be seen that for the smallest particles
$St_{\unicode[STIX]{x1D702}}$
is close to unity, while for large particles
$St_{L}$
is around 0.3. In order for particles to cluster, the time scale of the particles must be similar to some of the time scales in the flow (i.e.
$\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}\leqslant \unicode[STIX]{x1D70F}_{p}\leqslant \unicode[STIX]{x1D70F}_{L}$
). This means that in the case of these polydisperse systems, a wide range of eddies can contribute to clustering. Furthermore, since the size of the cluster depends on the Stokes number, clusters of various sizes may be present.
3 Results
In this section, the influence of turbulence on the mass transfer in polydisperse particle systems is studied by numerical simulations and the results are discussed and compared with previous results for monodisperse systems (Krüger et al.
Reference Krüger, Haugen, Mitra and Løvås2017a
; Haugen et al.
Reference Haugen, Krüger, Mitra and Løvås2018). Simulation parameters are presented in table 1. For each case, a series of simulations with different Damköhler numbers were performed. The Damköhler number was altered by changing the particle number density. To control the Stokes number, the particle material density was varied. In most cases, the Taylor microscale Reynolds number, defined as
$Re_{\unicode[STIX]{x1D706}}=u^{\prime }\unicode[STIX]{x1D706}/\unicode[STIX]{x1D708}$
, where the Taylor microscale is given by
$\unicode[STIX]{x1D706}=u^{\prime }\sqrt{15\unicode[STIX]{x1D708}/\unicode[STIX]{x1D716}}$
and
$u^{\prime }=u_{rms}/\sqrt{3}$
, is approximately 150. Additionally, for selected cases the domain size and the root-mean-square velocity (
$u_{rms}$
) were increased, resulting in a higher Reynolds number (
$Re_{\unicode[STIX]{x1D706}}\approx 350$
). This was done in order to obtain a larger-scale separation in the flow in order to better understand the effect of a polydisperse particle size distribution. Kinetic energy spectra for both Reynolds numbers are shown in figure 1. It can be seen that for the higher Reynolds number the spectra extends to lower wavenumbers, while the wavenumbers of the Kolmogorov scale (
$k_{\unicode[STIX]{x1D702}}$
) are similar for both
$Re$
.
3.1 Shape of size distribution
The normalised reactant decay rate (
$\tilde{\unicode[STIX]{x1D6FC}}=\langle \unicode[STIX]{x1D6FC}\rangle /\unicode[STIX]{x1D6FC}_{qsc}$
) obtained from the DNS for different particle size distributions is shown in figure 3 as a function of Damköhler number. Since the normalising factor is the decay rate computed as if the flow was quiescent,
$\tilde{\unicode[STIX]{x1D6FC}}$
shows the range of
$Da$
for which the mass transfer is enhanced (
$\tilde{\unicode[STIX]{x1D6FC}}>1$
) or decreased (
$\tilde{\unicode[STIX]{x1D6FC}}<1$
) due to turbulence. For the purpose of comparison, the equivalent results obtained by Haugen et al. (Reference Haugen, Krüger, Mitra and Løvås2018) with monodisperse particles are included as well (‘Dirac distr.’). The error bars in figure 3, as well as in all of the following figures, represent the standard deviation of the results. In the model (2.29), the reactant consumption rate for high
$Da$
is controlled by the cluster decay rate (
$\unicode[STIX]{x1D6FC}_{c}$
). This parameter is generally unknown; therefore, it was used as a fit parameter to make the model predictions (2.29) fit the simulation results (2.30). Its values are listed in table 1.
From figure 3, it can be seen that there is no significant difference between the normalised decay rates computed for
$\overline{\mathit{St}}_{L}=0.1$
with Gaussian, uniform and compensated distributions (cases K, J and D, respectively). The effect of the Dirac distribution will be discussed in a later sub-section.

Figure 3. Normalised decay rate as a function of Damköhler number obtained for different particle size distributions for
$\overline{\mathit{St}}_{L}=0.1$
(cases D, J, K).
3.2 Particle clustering and the effect of momentum back-reactions
It is not possible to change the Damköhler number while maintaining
$Re_{\unicode[STIX]{x1D706}}$
and
$St$
constant, since variations in particle mass loading will influence the root-mean-square velocity. Indeed, as was shown by Krüger et al. (Reference Krüger, Haugen, Mitra and Løvås2017a
),
$u_{rms}$
, and consequently also the Reynolds and Stokes numbers, decrease when the mass loading is high. In order to verify if this impacts the results, a series of simulations without back-reaction from the particles to the fluid, i.e. where
$\boldsymbol{F}=0$
in (2.2), were performed. Without particle feedback, the results are less representative in terms of physics but
$u_{rms}$
remains unchanged irrespective of
$Da$
. Figure 4 compares the obtained normalised reactant consumption rate as a function of
$Da$
. It can be seen that as long as the Stokes number is low, the difference between the results with and without back-reactions from particles to fluid is almost negligible. However, this is not the case for higher
$\overline{\mathit{St}}_{L}$
, for which the difference becomes significant at high Damköhler numbers. This is because the mass loading is proportional to the Stokes number for a given Damköhler number. An expression for the mass loading as a function of non-dimensional numbers is given in appendix A.

Figure 4. Comparison of the normalised decay rate obtained with and without particle back-reaction for simulations with
$\overline{\mathit{St}}_{L}=0.05$
(a), 0.35 (b) and 1.2 (c) (cases C, E, F).
The assumption has been that particles cluster purely due to turbulent eddies that are associated with a time scale that is similar to the particle response time. Then, if the location of turbulent eddies of different sizes were not correlated, particles of different sizes should cluster in uncorrelated positions. This is however not the case, as can be seen from the particle number density plots in figure 5. Here, the three upper (lower) panels present the results of a single simulation in which particle back-reactions were (were not) included, and the reported Stokes numbers correspond to different particle sizes. The strength of the dependency between particle locations can be demonstrated in a quantitative manner using correlation numbers. The correlation number of the number density of particles with different sizes is given as

where
${\tilde{n}}_{i}=n_{i}-\overline{n_{i}}$
and
$n_{i}$
is the particle number density of particles with size
$i$
. The corresponding values are presented in table 2 for simulations with and without back-reactions (i.e. one-way and two-way coupling). The Stokes numbers of the three different particle sizes are 0.04, 0.2 and 1, the same as in figure 5. It is clear from the table that the clusters of particles with different particle sizes are more correlated when the fluid can feel the presence of the particles (two-way coupling). The exception is the correlation between the two largest particle sizes, which is independent of back-reactions.

Figure 5. Particle number density for different Stokes numbers (case H). (a–c) Simulations with two-way coupling between particles and fluid; (d–f) only one-way coupling (no back-reactions from particles on fluid) is accounted for. White corresponds to zero particle number density, while black represents a particle number density that is 23 times larger than the average.

Figure 6. Particle number density for different Stokes numbers (case Q), one-way coupling,
$Re_{\unicode[STIX]{x1D706}}\approx 350$
.
It is also clear, however, that there is still some correlation between the position of clusters of different sizes also for the case with only one-way coupling. This correlation is due to the fact that the locations of turbulent eddies of different sizes are (weakly) correlated, even for isotropic turbulence that does not feel the presence of the particles.
Table 2. The correlation number between different Stokes numbers (
$St=0.04$
, 0.2 and 1), for simulations with and without back-reactions.

By inspecting the ‘strength’ or ‘sharpness’ of the clusters in figure 5 it can be observed that for the one-way coupling the sharpest structures are found for the smallest Stokes number. This is as expected since this case has a Kolmogorov-based Stokes number of
$\mathit{St}_{\unicode[STIX]{x1D702}}\approx 1$
, which is known to yield the sharpest clustering (Bec et al.
Reference Bec, Biferale, Cencini, Lanotte, Musacchio and Toschi2007). For the one-way coupling, the sharpness of the clusters is decreasing with increasing Stokes numbers. If we now consider the case with two-way coupling, we see that the trends are different. For the two-way coupling, the sharpest clusters are found for the intermediate Stokes number (
$\mathit{St}_{\unicode[STIX]{x1D702}}\approx 5$
), while also the largest Stokes number has sharper structures than is the case for the largest Stokes number with one-way coupling.
It is now interesting to see what happens if the Reynolds number is increased while
$\mathit{St}_{L}$
is kept unchanged. In figure 6 a contour plot of the particle number density is shown for
$\mathit{Re}_{\unicode[STIX]{x1D706}}\approx 350$
when particle back-reactions are turned off (one-way coupling). This is comparable to what is shown for
$\mathit{Re}_{\unicode[STIX]{x1D706}}\approx 150$
in figure 5(d–f). In figure 6, the integral-based Stokes numbers are the same as they are in figure 5, but the Kolmogorov-based Stokes numbers are different. From the figure, we see that for the case with
$\mathit{Re}_{\unicode[STIX]{x1D706}}\approx 350$
the clusters for the smallest Stokes number are not as sharp as for
$\mathit{Re}_{\unicode[STIX]{x1D706}}\approx 150$
. This difference is actually not due to the difference in Reynolds number, but rather due to the fact that the Kolmogorov-based Stokes number is 2 for the case with
$\mathit{Re}_{\unicode[STIX]{x1D706}}\approx 350$
while it is 1 for the smallest Stokes number case with
$\mathit{Re}_{\unicode[STIX]{x1D706}}\approx 150$
. This difference in sharpness when increasing the Kolmogorov-based Stokes number slightly beyond unity is consistent with previous findings by Bec et al. (Reference Bec, Biferale, Cencini, Lanotte, Musacchio and Toschi2007). For the larger Stokes numbers we see that, within error bars, the clustering is independent of
$\mathit{Re}_{\unicode[STIX]{x1D706}}$
. This is reasonable since the particles with
$\unicode[STIX]{x1D70F}_{p}\gg \unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$
are essentially independent of what happens at the Kolmogorov scale.
It can also be observed from figure 5 that cluster size and strength depend on the Stokes number and the way the fluid–particle coupling is handled. These features can be analysed using the auto-correlation function of the particle number density, which is given by

and is shown in figure 7. The faster the auto-correlation function decreases the more compact are the clusters. The distance at which
$C(r)$
crosses the
$x$
-axis indicates the cluster length scale. Thus, for the smallest particles the clusters are larger and more diffusive in the case of two-way coupling, whereas the opposite is true for larger particles. It is also interesting to note that for the larger particle sizes, the clustering is stronger when two-way coupling is invoked. The reason for the stronger clustering with two-way coupling is most likely that turbulence inside the heavy clusters will be suppressed by the presence of the particles, which means that the clusters will be relatively stable, while more particles can still be transported to the cluster. This transport of particles from volumes of high turbulence intensity to volumes of low turbulence is similar to the turbophoretic transport that is experienced for isotropic and non-homogeneous turbulence (Mitra, Haugen & Rogachevskii Reference Mitra, Haugen and Rogachevskii2018).

Figure 7. The auto-correlation function of the particle number density for simulations with and without back-reactions and for different particle sizes (case H).
In order to understand this effect better, we plot the probability density function (p.d.f.) of the particle number density for simulations with different mean Stokes numbers and for realisations with and without momentum back-reactions in figure 8. By inspecting the figure, it can be seen that the p.d.f. of the particle number density for simulations with
$\mathit{St}_{L}=1.2$
(figure 8
a,c) is significantly wider for the case with back-reactions (figure 8
a,b) than for the case without back-reactions (figure 8
c,d). In particular we see that for the smallest particle sizes, there is a much higher probability of finding sub-volumes where there are no particles. Likewise, the probability of finding sub-volumes with high particle number densities is also higher. This means that the particle clustering is stronger for cases with back-reactions. It can also be seen that back-reactions have a stronger effect for large
$\overline{\mathit{St}}_{L}$
(figure 8
a,c), while the effect is much less pronounced for smaller
$\overline{\mathit{St}}_{L}$
(figure 8
b,d). This is due to the fact that the mass loading scales with the Stokes number, which means that for small Stokes numbers the effect of back-reactions on the fluid is negligible.

Figure 8. The p.d.f. of
$n_{p}/\overline{n}_{p}$
for different particle sizes,
$Re_{\unicode[STIX]{x1D706}}\approx 150$
, compensated distr., (a,b) with back-reaction, (c,d) no back-reaction, (a,c)
$\overline{\mathit{St}}_{L}=1.2$
, (b,d)
$\overline{\mathit{St}}_{L}=0.05$
, in all cases
$Da\approx 35$
(cases C and F).
3.3 Effect of back-reactions on reactant concentrations
From figure 4 it is seen (through the fact that
$\tilde{\unicode[STIX]{x1D6FC}}$
is above unity for larger
$\mathit{Da}$
) that the effect of clustering is less when momentum back-reactions from the particles to the fluid are neglected, i.e. when
$\boldsymbol{F}$
in (2.2) is set to zero. The reason for this can be understood from the previous sub-section, which showed that particle clustering is stronger due to momentum back-reactions. More insight into this can be gained from figure 9, which shows the p.d.f. of the reactant mole fraction normalised by its mean value. Here, the curves referred to as ‘constrained’ are given by

and represent the p.d.f. obtained if data are collected at the position of the particles,
$\unicode[STIX]{x1D6FF}$
is the Dirac delta function,
$N_{part}$
is the total number of particles in the domain and
$X(\boldsymbol{r}_{j})$
is the reactant concentration at the position of particle
$j$
. The ‘not constrained’ curves are obtained based on the reactant mole fraction present in the entire domain and are given by

when
$N_{grid}$
is the total number of grid points and
$X_{i}$
is the reactant concentration in grid cell
$i$
. It becomes clear from figure 9 that the reactant distribution is significantly narrower when
$\boldsymbol{F}=0$
, which is particularly pronounced for
$\overline{\mathit{St}}_{L}=1.2$
(figure 9
a). A narrow distribution means that all particles have access to a similar amount of reactant. Hence, the clustering does not slow down the conversion rate so much.

Figure 9. The p.d.f. of
$X/\overline{X}$
, (a)
$\overline{\mathit{St}}_{L}=1.2$
, (b)
$\overline{\mathit{St}}_{L}=0.005$
, compensated distribution,
$Da\approx 35$
(cases C and F).
3.4 The effect of turbulence on the overall mass transfer rate
Having gained an understanding of the interactions between differently sized particles in a polydisperse system, we can now study how the mass transfer rate is influenced due to turbulence. We begin by analysing p.d.f.s of particle number density and reactant mole fraction for different particle sizes for case N, in which
$\overline{\mathit{St}}_{L}=0.04$
and
$Re_{\unicode[STIX]{x1D706}}\approx 350$
. These results, presented in figure 10, are the most representative for the considered polydisperse particle system since
$\unicode[STIX]{x1D70F}_{p,min}\approx \unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}$
and
$\unicode[STIX]{x1D70F}_{p,max}<\unicode[STIX]{x1D70F}_{L}$
, which means that for all particle time scales there exists a turbulent flow time scale of the same order. Hence, turbulence can potentially make all particle sizes cluster. Despite the fact that the range of flow scales was significantly narrower for cases studied in previous sections (with
$Re_{\unicode[STIX]{x1D706}}\approx 150$
), the p.d.f. of the particle number density presented in figure 10(a) verifies previous findings for the smaller Reynolds number. Here, again, the broadest p.d.f. is obtained for particles with
$\mathit{St}_{L}$
of the order of 0.1 (which corresponds to the largest particles in this simulation). As a consequence, we expect that the mass transfer rate will be mostly affected by clustering of these largest particles. The corresponding p.d.f.s of the reactant mole fraction (figure 10
b) show that there is indeed slightly less reactant available at the locations of the largest particles. This is not surprising since they are more clustered and the clusters are larger, yielding longer cluster life times.

Figure 10. The p.d.f.s of (a)
$n_{p}/\overline{n}_{p}$
and (b)
$X/\overline{X}$
for different particle sizes,
$Re_{\unicode[STIX]{x1D706}}\approx 350$
, compensated distr., with back-reaction,
$\overline{\mathit{St}}_{L}=0.04$
,
$Da=25$
(case N).
Normalised reactant decay rates for a compensated particle size distribution,
$Re_{\unicode[STIX]{x1D706}}\approx 150$
and different mean Stokes numbers (corresponding to cases A, B, D and F) are shown in figure 11(a). The tendency here is the same as previously observed by Haugen et al. (Reference Haugen, Krüger, Mitra and Løvås2018). As long as the Damköhler number is low (or, in other words, as long as the mass transfer is not affected by the particle clustering), a higher normalised reactant decay rate is obtained for higher
$St$
, which is a consequence of the fact that the mean Sherwood number increases with Stokes number (see table 1). For higher Damköhler numbers, the reactant consumption rate starts to be dependent on the cluster decay rate. As expected, in the case of
$\overline{\mathit{St}}_{L}=0.005$
,
$\tilde{\unicode[STIX]{x1D6FC}}$
begins to decrease only for relatively high
$Da$
, which can be interpreted as a weak effect of particle clustering due to the fact that for smaller Stokes numbers the corresponding clusters have shorter life times. As the Stokes number increases, the effect of clustering leads to a fast decrease of the normalised decay rate. This effect is again less when
$\overline{\mathit{St}}_{L}$
is further increased beyond unity since these heavy particles are less sensitive to the flow. The equivalent results for
$Re_{\unicode[STIX]{x1D706}}\approx 350$
(cases L, M, O and P) are presented in figure 11(b) from which it can be seen that the conclusions drawn for
$Re_{\unicode[STIX]{x1D706}}\approx 150$
are also true for the higher Reynolds number.

Figure 11. Normalised decay rate as a function of Damköhler number obtained for different
$\overline{\mathit{St}}_{L}$
using a compensated distribution for (a)
$Re_{\unicode[STIX]{x1D706}}\approx 150$
(cases A, B, D, F) and (b)
$Re_{\unicode[STIX]{x1D706}}\approx 350$
(cases L, M, O, P).
In order to assess if the effect of particle clustering is significant for a given set of parameters we employ the following quantitative description. It follows from (2.29) that when
$\mathit{Sh}_{mod}/\overline{Sh}=1/2$
, i.e. when the effect of particle clustering has reduced the mass transfer rate to half of what it would have been if particle clustering was neglected, the Damköhler number is given by

This is presented in figure 12, from which it is seen that in all of the investigated cases the effect of turbulence is at its strongest when
$\overline{\mathit{St}}_{L}\approx 0.1$
, which is where
$\mathit{Da}_{2}$
takes the lowest value (note that the lower the value of
$\mathit{Da}_{2}$
, the greater the influence of particle clustering). This is due to the fact that this yields a combination of relatively large and sharp clusters. For
$\overline{\mathit{St}}_{L}<0.1$
the clusters are smaller, while for
$\overline{\mathit{St}}_{L}\gg 0.1$
the clusters are weaker because
$\mathit{St}_{L}>1$
for the largest particles, meaning that the particle concentrations are more homogeneous. For
$Re_{\unicode[STIX]{x1D706}}\approx 350$
(figure 12
b), there is almost no difference between mono- and polydisperse particle systems. The same trend is also seen for
$Re_{\unicode[STIX]{x1D706}}\approx 150$
(figure 12
a), although slightly lower values of
$\mathit{Da}_{2}$
indicate a stronger tendency to clustering of the monodisperse system. It can also be observed that the effect of turbulence is slightly stronger for higher
$Re_{\unicode[STIX]{x1D706}}$
. This is probably due to the fact that in this case a broader range of flow scales contribute to particle clustering, i.e. there are more particles for which
$\unicode[STIX]{x1D70F}_{\unicode[STIX]{x1D702}}<\unicode[STIX]{x1D70F}_{p}<\unicode[STIX]{x1D70F}_{L}$
. While for
$\overline{\mathit{St}}_{L}<0.1$
,
$\mathit{Da}_{2}$
is not much affected when the back-reaction from particles is neglected, there is a significant difference in the results when
$\overline{\mathit{St}}_{L}>0.1$
, which is in line with the conclusions already drawn in previous sections.

Figure 12. Damköhler number
$Da_{2}$
as a function of Stokes number, all cases with compensated and Dirac distributions.
We have shown that the effect of turbulence on reactive particles is due to a modified mass transfer coefficient, i.e. a modified Sherwood number. Chemical kinetics, however, is not directly affected by the turbulence, only indirectly through the reactant concentration. For purely kinetics-controlled reactions, this means that turbulence would not affect the conversion rate at all. The results presented in this work are formally applicable only to cases with diffusion-controlled reactions (fast kinetics). In reality, however, using the modified Sherwood number as described in this work will be correct also when the reactions are not diffusion-controlled. That is, for all kinetic rates, the correct rate of the reactions is obtained by using the traditional approach, which combines both kinetics and diffusion (see e.g. Haugen et al. Reference Haugen, Mitchell and Tilghman2015), as long as the diffusion is described with the use of the modified Sherwood number.
4 Conclusions
The effect of turbulence on the mass transfer rate in a dilute, polydisperse particle system was analysed over a range of conditions. We show that for polydisperse systems, the reaction rate is affected by particle clustering in the same way as for monodisperse systems. Even though particles of various sizes differ in the way they are distributed in the domain, this effect can be as strong as in a monodisperse system, provided that the scale separation in the flow is sufficiently large, or in other words, the Reynolds number is sufficiently high. When
$\overline{\mathit{St}}_{L}$
is of the order of
$10^{-1}$
, the rate of mass transfer can be reduced by 50 % even for
$Da$
lower than 10. It is therefore clear that, when studying real systems, the effect of turbulence on the overall mass transfer rate should be accounted for. The model given by (2.28) and (2.29) allows the incorporation of this effect directly in Reynolds-averaged Navier–Stokes-based codes so that the model can be applied in simulations of large-scale reactors and boilers. It is found that the mass transfer rate is not very dependent on the shape of the particle size distribution, it is the width that matters.
We also observed that, despite having different sizes and hence different Stokes numbers, particles cluster in correlated positions. This correlation is greater when momentum transfer between fluid and particles is two-way. Both the distribution of reactant and particle number density are found to be broader when two-way coupling is applied. A direct consequence of this is that the mass transfer rate between the particles and the fluid is reduced due to the back-reaction of particle momentum to the fluid when the mass loading is significant.
Finally, it is worth pointing out that in reality the particle diameter (and/or material density) decreases as the surface reaction progresses. This will cause the Stokes number and the Damköhler number to decrease since
$Da\sim d_{p}$
and
$St\sim d_{p}^{2}$
. As a result, the effect of turbulence may be different at initial and final stages of particle conversion. If the particle shrinkage model is used in a simulation, the effect it has on the particle conversion rate will be automatically included by recomputing the reactant decay rate every time step, such that
$\tilde{\unicode[STIX]{x1D6FC}}=\tilde{\unicode[STIX]{x1D6FC}}(d_{p}(t))$
.
Acknowledgements
The research leading to these results has received funding from the research project ‘Gaspro’, financed by the Research Council of Norway (267916). N.E.L.H. also acknowledges the Research Council of Norway under the FRINATEK grant 231444. This research was supported in part with computational resources provided by NOTUR (grant NN9405K). The work of A.K. was supported by the statutory research fund of the Silesian University of Technology, Faculty of Energy and Environmental Engineering, Gliwice, Poland.
Appendix A. Mass loading
The mass loading is given by

when

and
$V_{p,i}=(\unicode[STIX]{x03C0}/6)d_{p,i}^{3}$
,
$n_{p,i}$
and
$d_{p,i}$
are the volume, radius and particle number density of particles with size
$i$
. The chemical time scale is now given by

such that

when

By this, the average Stokes number becomes

Combining the above equations yields

since
$\overline{\mathit{St}}_{L}=\overline{\mathit{St}_{\unicode[STIX]{x1D702}}}/\sqrt{\mathit{Re}}$
. One can also avoid all reference to any turbulence property by introducing the non-dimensional number

such that
