## Abstract

Engine durability tests are used by manufacturers to demonstrate engine life and minimum performance when subjected to doses of test dusts, often Arizona Road Dust. Grain size distributions are chosen to replicate what enters the engine; less attention is paid to other properties such as composition and shape. We demonstrate here the differences in the probability of interaction of a particle of a given particle Reynolds number on to a vane if particle shape, vane geometry, and flow Reynolds number are varied and discuss why the traditional definition of Stokes number is inadequate for predicting the likelihood of interaction in these flows. We develop a new generalized Stokes number for nozzle guide vanes and demonstrate its use through application to 2D sections of the General Electric E^{3} nozzle guide vane. The new Stokes number is used to develop a reduced-order probability curve to predict the interaction efficiency of spherical and nonspherical particles, independent of flow conditions and vane geometry. We show that assuming spherical particles instead of more realistic sphericity of 0.75 can lead to as much as 25% difference in the probability of interaction at Stokes numbers of around unity. Finally, we use a hypothetical size distribution to demonstrate the application of the model to predict the total mass fraction of dust interaction with a nozzle guide vane at design point conditions and highlight the potential difference in the accumulation factor between spherical and nonspherical particles.

## 1 Introduction

Predicting the damage to gas turbine engines operating in harsh environments is made difficult by the multiple variables that characterize the environmental dust. The dust undergoes significant changes in size and mineralogical distribution through the engine due to processes such as sorting by the particle separator or fan, fragmentation in the compressor, agglomeration in the combustor, phase change, and reactions along the gas path. Furthermore, the proportions of minerals vary across the range of particle sizes; the original ingested sediment may contain a high proportion of quartz in coarse grain fractions, but a much greater proportion of lower melting point clays in the smaller grain fractions.

This is significant because the likelihood of particle-wall interaction is a function of particle size, shape, and density. The smaller the particle, the greater the drag per unit volume thus tendency to follow fluid streamlines. However, the smaller the particle, the faster the time to soften or melt under heating. Hence, the probability of deposition is a function of the physical properties that influence the trajectory and the chemical properties that influence the phase change. To predict the risk of operating in a given environment, therefore, requires knowledge of the composition of the material at different points throughout the engine. Since the proportion of each mineral varies across the size range, it is likely that the composition of the material that interacts with the nozzle guide vane will be markedly different to the bulk dust first entering the engine. This has implications for engine subcomponent dust durability testing, where a suitable test dust must be chosen. Our objective in the current work is to develop a function that can be used to ascertain interaction probability for any given mineral, of any size in any flow, in order that the composition of the dust impacting the nozzle guide vane may be better predicted.

## 2 Background

There have been a number of incidents in the past of aircraft ingesting quantities of environmental particulate (EP) and suffering engine performance loss. Clarkson et al. [1] and Clarkson and Simpson [2] used these events in combination with several other, less severe encounters to develop a *Duration of Exposure* versus *Ash Concentration* chart to identify the potentially hazardous *dose* events. This approach correlates the severity level of exposure events with a numerical measure of the total mass of particulate matter ingested *per unit of engine volume flow rate* and is becoming more accepted among civil operators and manufacturers as a metric to determine operational limits.

If the dose of an event and the engine mass flow at the time are known, one can estimate the mass of dust ingested. In the case of a turbo-fan powered fixed-wing aircraft, this may be reduced by centrifuging action of the fan (see Ref. [3]) or in the case of a military helicopter by a particle separator [4]. In each case, the reduction is a function of the properties of the particulate and flow. This process is mirrored in other parts of the engine as particles deposit at different locations along the gas path.

### 2.1 Stokes Number.

If a particle reaches the turbine stages via the core flow, it may impact adhere to the blade surface, causing a build-up of material, as shown in Fig. 1, in this case on the leading edge [5]. In other designs, the deposit builds up on the vane midchord region, thought to be due to the position of the combustor flame relative to the blade. The presence of a surface deposit causes an increase in surface roughness and a narrowing of the throat area. The effect of the former is to thicken the boundary layer that increases the profile drag thereby reducing the efficiency of the stage. A consequence of the latter is to create a back pressure that narrows the surge margin. The deposit may also act as a tacking agent, thereby preventing further bounce. These effects are nonlinear and depend upon the rate of delivery and the way in which the particulate interacts with the surface.

The probability of interaction depends upon the approach position relative to the blade (i.e., a particle on a leading edge stagnation streamline is more likely to hit than one initially headed between two adjacent blades) and the Stokes number of the particle-flow system. A Stokes number greater than unity (*Stk* ≫ 1) implies a large deviation from the flow streamlines and a high probability of impact; a Stokes number much less than unity (*Stk* ≪ 1) will closely follow the fluid streamlines and avoid interaction. A Stokes number of around unity will have a probability of impact between zero and unity.

*generalized*Stokes number developed by Israel and Rosner [6] does and is defined as

*μ*

_{f}are the density and viscosity of the particle and flow, respectively,

*d*

_{p}is the particle diameter, $\tau pv$ is the particle velocity response time in Stokes flow, and $\tau fv$ is the flow velocity response time, defined by a characteristic velocity

*U*and length

*L*that are the characteristics of the fluid-obstacle interaction. The non-Stokes drag correction factor $\psi $ moderates the original Stokes number to reflect the inertial drag effects. It can be approximated by the following function:

_{p0}is the initial freestream particle Reynolds number, Re

_{p}is the instantaneous relative particle Reynolds number, and

*C*

_{D}is the drag coefficient. $\varphi $ decreases with increasing Re

_{p}, which means the difference between the generalized and original Stokes numbers gets larger with increasing Re

_{p}. It allows one to calculate an

*effective*Stokes number, i.e., what the Stokes number would have to be to reproduce the same particle-obstacle interaction outcome. If the correction factor were not applied, the drag per unit mass at a high particle Reynolds number would be much greater than that predicted by Stokes law, which would reduce the likelihood of interaction (due to a greater acceleration being set up in the presence of a fluid velocity gradient, all other parameters remaining equal). The traditional form of the Stokes number however is much simpler to calculate and is more widely used, e.g., Ref. [7].

### 2.2 Particle Reynolds Number.

*C*

_{D}is the particle drag coefficient,

*d*

_{p}is the particle diameter,

*μ*

_{f}is the fluid viscosity, and (

**u**−

**v**) is the relative velocity between the particle and fluid. The particle Reynolds number Re

_{p}is defined as

**u**and

**v**are the velocities of the carrier and particle phases, respectively. In numerical simulations of particle transport, the drag coefficient is computed based on a model such as that of Haider and Levenspiel shown in Fig. 2 [8]. In this model, the inertial effects begin at an increasingly smaller Reynolds number as the particle becomes less spherical due to a larger surface area per unit volume. This translates to greater drag and an increased ability to follow fluid streamlines at the same Reynolds number. If the nonsphericity of the test dust is not considered or modeled, there may be an over-prediction of deposition rate.

### 2.3 Influence of Mineralogy.

Stokes number is often used to scale the boundary conditions in gas turbine component and subcomponent rig tests, such as in Ref. [7]. In this work, the authors used the original Stokes number (Eq. (1) excluding $\psi $) to tune the particle boundary conditions to simulate the full-scale conditions. The deposition rate was found to be most adequately simulated by matching flow and particle velocity and temperature. However, the same study also highlighted the issue of inertial effects. The absence of a non-Stokes drag correction neglects the effect of pressure, which is implicit in the fluid density in Eq. (4). If the pressure is not matched in the rig test, the density of the gas in will be lower than in the real engine, leading to lower drag and a reduced tendency to follow fluid streamlines. The result is a higher deposition rate, despite being at the Stokes number.

Particle Reynolds number is proportional to fluid density, particle diameter, and relative velocity and inversely proportional to fluid viscosity. The viscosity and fluid density are additionally dependent on the fluid temperature. The velocity and temperature field in a gas turbine engine can vary considerably in all directions, while the particle grain size on entry to the engine can be expected to be up to 100 *μ*m. The particle response time, which dictates the likelihood of interaction, is additionally influenced by the particle density and particle shape. These particle properties can also vary widely, for example, across geographical locations due to the underlying geology and localized historical variations in the wind- or water-borne processes. The mineralogy can also vary as the particles pass through the engine due to phase change, chemical reactions, or high impact fracturing, as revealed by Smialek et al. [9]. In this work, the chemical composition of deposits in the cooling holes and on the vane surfaces was found to differ to that of the raw dune dust. Clearly, the dust mineralogy by which Stokes and Reynolds number are influenced affects where contaminant particles are distributed throughout the engine.

### 2.4 Accumulation Factor.

*c*

_{vp}is the particulate concentration expressed as a volume per unit mass,

*W*

_{1}is the air mass flow rate through the engine core, and $\zeta NGV$ is what is termed the

*accumulation factor*. To complete the estimation of the accumulated mass, Clarkson and Simpson [2] include an efflux term to account for the shedding of mass through aerodynamic forcing or erosion of the deposited material. However, this approach does not currently contain enough information to differentiate the severity of the threat between EP type and engine. The accumulation factor and shedding rate are assumed to be constant for a given contaminant type, engine design, and engine operating point; this assumption has been used to back calculate doses of historical events to verify the regions of safety margin deterioration for the dose approach. To date, these parameters have been based on the empirical data, and the rig-level efforts are ongoing to refine the parameters to include more physics. Giehl et al. [10] associate the adhesion and shedding ability to the proportion of crystalline and quartz material in the dust. They demonstrate for a number of ash types; there is an increase in the accumulation rate due to surface wetting impact with increasing glass content. This behavior is easier to predict for amorphous ash with predictable variation in melt viscosity versus temperature. However, mineral dusts exhibit more complex changes in melt viscosity under heating; therefore, their deposition behavior is more difficult to predict.

The uncertainty in the accumulation factor is reflected in the wide range of values that are quoted in the literature, some of which are summarized in Table 1 [11–15]. The ranges cover all data points published by the authors. The variability arises due to the many dependencies: tests at different gas temperature, substrate temperature, substrate roughness, concentration, dust type, and Stokes number all give rise to variations in the deposition rate. There is currently no probabilistic model to correlate the accumulation factor with these properties; more experimental research work is required to understand the primary and secondary drivers. The use of reduced-order models in this effort may help to reveal other primary drivers, more closely linked to the particle mineralogy such as crystallographic structure or melting behavior. In the current work, the probability of particle accumulation is treated as a product of two probabilistic processes: *interaction* and *retention*; if the probability of interaction can be reduced to a function in one or two parameters, then the factors affecting the likelihood of a particle sticking to a surface may be more readily investigated.

## 3 Theory

### 3.1 Capture Efficiency.

Consider an NGV in a cross flow laden with liquid droplets and assume that any droplet that impacts the surface is considered *captured*. As the fluid streamlines pass around the NGV, the likelihood that a droplet will be captured will depend upon the following:

Generalized Stokes number, which accounts for particle and fluid freestream properties;

Initial starting position relative to the cylinder on an ordinate perpendicular to the freestream.

The Stokes number accounts for the inertia of the fluid–particle system. But without knowledge of the starting position, two particles with identical Stokes number can have entirely different fates: one may interact while the other may bypass the cylinder, as depicted in Fig. 3. The uncertainty introduced by starting position is only relevant to particles with a Stokes number of around unity: particles with *Stk* ≫ 1 with a starting position within the projected area shadow will impact the cylinder, while particles with *Stk* ≪ 1 will follow the fluid streamlines. This gives rise to a probability density curve for the *likelihood of interaction* that depends upon Stokes number alone. Israel and Rosner demonstrated this relationship for a cylinder in a cross flow of liquid aerosol particles if the Stokes number takes the form given in Eq. (1) [6]. In their study, five different *freestream* particle Reynolds numbers were selected over a range of generalized Stokes numbers to generate the data. The authors fitted a single-parameter curve to the data for future use as a single-parameter prediction.

In the present work, we aim to create a similar reduced-order curve fit for NGVs. The post-impact fate of the particle is not considered; instead, the focus is on the likelihood of a particle making contact with the NGV. In contrast to liquid aerosols, the fate of a dust particle after impact with an NGV surface is not straightforward to predict; this is the subject of ongoing work. Since naturally occurring minerals appear in different proportions across wide distributions of size and shape, it is likely that the composition of the dust interacting with the NGV compared with the dust entering the stage will be different. Since the mechanics of retention depends so greatly on the physicochemical properties of the particle (e.g., yield strength, Poisson’s ratio, Young’s modulus, and density—see Ref. [16]), this may permit a more accurate prediction of the overall mass that will accumulate on a nozzle guide vane.

### 3.2 Particle-Fluid Interaction.

**u −**)

*v*_{0}is the initial velocity difference and

**g**is the acceleration due to gravity. For gas-particle flows, where the ratio of the carrier flow density to particle material density is very small (≃10

^{−3}), the BBO can justifiably be simplified to

^{−3}), the suspension can be considered dilute and particle–particle interactions are negligible, as discussed in Sec. 2.4. The remaining effects can be modeled in the drag coefficient, which is contained within the particle velocity response time $\tau pv$.

### 3.3 Response Time.

*μ*m particle of density 2551 kg m

^{−3}in a gas of viscosity 5.3 E−05 kg m

^{−1}s

^{−1}is 1.6 cm s

^{−1}. The transit time of a particle traveling at 100 ms

^{−1}from combustor to NGV (10 cm) is of the order of milliseconds; in this time, a particle traveling perpendicular to the gravitational field force would fall approximately 0.0016 cm if only subjected to Stokes drag. Hence, the effect of gravity is generally negligible for the deposition of small particles in turbomachinery flows (see also Ref. [7]). Equation (7) therefore reduces to the steady-state equation of motion for an isolated particle in a gas. Re-writing Eq. (7) in terms of drag coefficient and Reynolds number, we have

*C*

_{D}Re

_{p}/24 approaches unity. For particle Re

_{p}> 1, it is known as the

*drag factor*

*f*and is the ratio of drag coefficient to Stokes drag. The solution of Eq. (8) for constant

**u**and an initial particle velocity of zero is

*e*− 1)/

*e*) of the free stream velocity. In Stokes flow, this is simply the

*response time*of the particle. However, if the relative Reynolds number is greater than unity, the value of

*f*grows considerably and has the effect of shortening the time taken to achieve 63% and indeed 99% of the freestream velocity. The product $\tau vp/f$ could be termed the

*non-Stokesian*response time. There are several correlations available in the literature for

*f*as a function of Reynolds number. A simple yet relatively good correlation for Re

_{p}is given by Schiller and Nauman as reported in Ref. [18]:

*C*

_{D}-Re

_{p}curve between the initial particle Reynolds number Re

_{p0}and 0. This is the basis for the non-Stokes drag correction factor $\psi $ given by Israel and Rosner in Eq. (2). Performing the integration using the Schiller–Naumann correlation for the

*C*

_{D}(Re

_{p}) term yields a function for the non-Stokes drag correction factor as a function of Reynolds number alone:

_{p}= 800. This is demonstrated in Fig. 4, which shows the particle velocity response

*v*to a step change in carrier velocity

*u*for a particle of density 2550 kg m

^{−3}and diameter 10

*μ*m in a flow of viscosity 5.3 × 10

^{−05}kg m

^{−1}s

^{−1}and an initial particle Reynolds number of 100. To generate the data, a numerical integration was performed for the acceleration of a one-dimensional particle in a carrier flow of constant velocity with gravity and unsteady effects neglected. A time step of 1 × 10

^{−08}s was used, and the iteration stopped when the particle reached 99.9% of the freestream velocity. The drag coefficient was held constant for the Israel–Rosner case and can be compared with the exact solution using the Schiller–Naumann correlation for the drag coefficient. The Stokesian drag coefficient is plotted for comparison to show the implication of ignoring the Reynolds effects in cases where the particle will access the non-Stokesian range of the drag coefficient curve: the particle reaches 63% of the freestream velocity in a third of the time. The implication holds true if Stokesian response time alone is used in the computing of Stokes number; a particle of given Reynolds number will respond much quicker to changes in the flow momentum (direction or magnitude), and as such may avoid obstacles it was otherwise expected to hit. The implication for dynamic scaling in rig testing is elaborated by Sacco et al. [7].

*s*is the surface area of a sphere having the same volume as the particle and

*S*is the actual surface area of the particle. This is a far more complex expression and cannot be used to solve directly the integral in Eq. (1). For a spherical particle, there is a little difference between the Schiller–Naumann and Haider–Levenspiel drag correlations on the particle velocity history—see Fig. 4. However, when computing the time history of a nonspherical particle, there is a noticeable difference in the

*non-Stokesian*response time. Figure 5 shows the non-Stokesian response time ($\tau vp$) as a function of

*initial*particle Reynolds number, normalized with response time of a spherical particle of the same volume-equivalent diameter. A numerical integration was performed to produce each curve, with the iteration stopped once the particle had reached 63% of the freestream velocity.

Different composition minerals have different sphericities due to their crystallographic structure. The sphericity of a particle is observed to decrease with particle size for the same mineral. However, in this study, the sphericity values are fixed across the size ranges but chosen to compare the behavior of a pure sphere ($\varphi =1$) to more realistic mineral sphericities of $\varphi =0.75$ and $\varphi =0.54$, which for example are the unit cell structures of quartz and gypsum, respectively [19]. For the purposes of comparison, we also fix the density although in nature quartz (2.65 g cm^{−}^{3}) and gypsum (2.31 g cm^{−}^{3}) do not have identical density. What is interesting to note is that response time for a gypsum particle at an initial Reynolds number of 100 is half that of the spherical particle—i.e., twice as quick to respond to changes in the flow momentum. This reflects the differences in the drag coefficient curves presented in Fig. 2.

## 4 Methodology

The preceding discussion has demonstrated some of the shortcomings of using the Stokesian response time that can be overcome through the use of a correction factor. Using this to compute a generalized Stokes number, it has been shown that the prediction of cylinder capture efficiency can be reduced to a two-parameter equation. However, the particles are limited to spherical droplets. In the present study, we aim to

Develop a generalized Stokes number for nozzle guide vanes

Propose a reduced-order function to predict the interaction probability of a mineral dust with a nozzle guide vane

Investigate the effect of particle nonsphericity on interaction probability

We conduct a numerical investigation on the widely investigated GE-E^{3} turbine vane geometry using the flow solver and particle tracking algorithm in starccm+. Our simulations are limited to two dimensions in order to draw a comparison with the cylinder results and isolate the circumferential divergence effect of the obstacle, but analyze three 2D sections along the span of the nozzle guide vane: near-hub, midspan, and near-endwall. Three-dimensional effects are thus not modeled here. The comprehensive data set in Ref. [20] is used to validate the numerical simulations. A photograph of the nozzle guide vane row is shown in Fig. 6.

### 4.1 GE-E^{3} High-Pressure Turbine.

The high-pressure turbine for the General Electric Energy Efficient Engine is a two-stage through-flow design of moderate loading. A full-scale warm air turbine test rig with cooling flows was run to evaluate the aerodynamic performance of the HPT and is reported in Ref. [20]. To achieve this, the design point efficiency was determined, and the turbine was mapped over a large range of operating points including the sub-idle and several off-design Reynolds numbers. This gives us a suitable range of operating conditions with which to generate a range of Stokes numbers and particle Reynolds numbers to generate the interaction probability curves. The vane and stage design properties are given in Table 2.

Property | Unit | Hub | Mid | Tip |
---|---|---|---|---|

Number of vanes | – | 46 | 46 | 46 |

Radial position | cm | 32.576 | 34.576 | 36.576 |

Pitch (spacing) | cm | 4.450 | 4.723 | 4.996 |

Axial width | cm | 3.376 | 3.378 | 3.383 |

Chord length | cm | – | 6.201 | – |

Trailing edge height | cm | – | 4.00 | – |

Trailing edge thickness | cm | – | 0.0965 | – |

Uncovered turning | deg | 9.2 | 8.4 | 8.7 |

Trailing edge wedge angle | deg | 10.2 | 9.2 | 9.0 |

Property | Unit | Hub | Mid | Tip |
---|---|---|---|---|

Number of vanes | – | 46 | 46 | 46 |

Radial position | cm | 32.576 | 34.576 | 36.576 |

Pitch (spacing) | cm | 4.450 | 4.723 | 4.996 |

Axial width | cm | 3.376 | 3.378 | 3.383 |

Chord length | cm | – | 6.201 | – |

Trailing edge height | cm | – | 4.00 | – |

Trailing edge thickness | cm | – | 0.0965 | – |

Uncovered turning | deg | 9.2 | 8.4 | 8.7 |

Trailing edge wedge angle | deg | 10.2 | 9.2 | 9.0 |

### 4.2 Computational Mesh.

Three 2D cross-sections from the GE-E^{3} HPT geometry coordinates given in Ref. [20] were selected as test cases. The planes are depicted in Fig. 7. A computational domain was constructed around each 2D geometry. Each domain extends one vane chord length upstream and half a vane chord length downstream of the row. The confined space for the wake to develop reflects the position of the next blade row. The inlet and outlet were designated pressure boundaries, while the north and south boundaries were designated as periodic to reflect the repeating nature of the annular row and allowing for a reduction in computational expense. The domain height reflects the spacing between each vane in the row (vane pitch). The vane surface was designated as a perfectly smooth, no-slip wall.

Grid generation for the three domains was carried out using the in-built meshing features of starccm+. The grids were all unstructured and composed of tetrahedral elements. A structured 40-layer prism mesh was employed to resolve the flow in the near-wall boundary layer. The growth of the prism layers was tailored to give a wall *y* + of less than unity for all converged solutions. This ensured that the turbulent boundary layer was always fully resolved, hence eliminating the need to apply wall functions. For each of the three geometry cases, a grid independence check was conducted by monitoring the skin friction coefficient profile and the average skin friction coefficient at the NGV wall. Grid independence was achieved with meshes of 108,050, 118,272, and 121,888 cells for the hub, mid, and tip sections, respectively. The continuous phase took 0.5–1.6 CPU hours to reach convergence. The meshed domain of the midspan section is shown in Fig. 8.

### 4.3 Continuous Phase Modeling.

The conservation equations for mass, momentum, and energy for a continuous, compressible, and turbulent flow were solved to give the pressure and velocity distributions around the NGV at each desired operating point. The Navier–Stokes equations were solved in a segregated manner using a second-order upwind discretization. Pressure-velocity coupling was handled using the hybrid Rhie-Chow predictor-corrector approach with the SIMPLE algorithm contained within the starccm+ code. Temperature variations within the flow domain were determined using a segregated fluid temperature approach. This model solved the total energy equation with temperature as the solved variable. Enthalpy was then computed from this using the equation of state which in this case was the ideal gas law.

Turbulent fluctuations in the velocity field were computed using the shear-stress transport k–*ω* form of the Reynolds-averaged Navier–Stokes equations. A second-order scheme using this formulation was selected due to its favorability for situations containing nonzero pressure gradients [21]. The turbulence quantities of turbulent intensity and length scale were specified based upon those determined experimentally in Ref. [16]. The turbulence intensity was specified at the inlet and outlet boundaries to be 5%, and the turbulent length scale based upon the integral length scale of 2.216 cm determined in the Ohio State TuRFR deposition rig [16].

For each case, the continuous phase was considered converged when the residuals decreased over five orders of magnitude and when the average skin friction, lift, and drag coefficients on the vane showed less than 0.001% change over ten solution iterations. The converged solution for the continuous phase at a vane Reynolds number of 176,480 has been validated for all three geometries. The NGV surface Mach number distributions at the hub, tip, and midspan from the experimental test of the GE-E^{3} HPT (Ref. [20]) have been compared with the computational solution, as shown in Fig. 9. The surface Mach number distribution over the midspan vane geometry is compared with the experimental data and original GE-E^{3} design intent for the *solid vane*. The numerical simulations show good agreement with the experimental data in the leading 65% portion of the suction surface and the majority of the pressure surface. However, an additional supersonic region is predicted close to the trailing edge suction side of the numerical solution, as shown in Fig. 10. This is thought to be due to neglecting the surface roughness introduced by the cooling holes and frictional losses due to the endwalls. Since the interacting particles impact almost exclusively on the pressure surface (unless bouncing is modeled), we do not expect this difference to impact the results of the present study.

### 4.4 Dispersed Phase.

The dispersed phase in the simulations was modeled using the *Lagrangian Multiphase* model in starccm+. The maximum mass loading of ≃10^{−10} is small enough that the momentum exchange with the fluid can be considered one way (see Sec. 2.4). Once the flow field was solved, particles were injected into the domain inlet with uniform spacing. Particles of constant diameter and density were injected at a rate of one per iteration from 2001 discrete, equally spaced locations along the inlet boundary. They were then tracked through the domain until they either impacted the NGV or exited the domain. Each particle was injected from rest, in both the same direction and in thermal equilibrium with the continuous phase. The particle trajectories were then computed by advancing the particles through the domain over a series of time steps. The size of this step is adjusted dynamically by the code according to the maximum and minimum Courant numbers for the given grid and flow velocity. Additional limits are placed upon the particle time step size by the maxima and minima of the momentum relaxation and eddy interaction times for the particle drag and turbulent dispersion models.

The particle trajectory is computed by summing the forces on the particle that arise from its position in the solved continuous phase and its velocity relative to the carrier fluid, and calculating the resultant acceleration, from which a change in velocity can be found. The new particle velocity and position are calculated, and the process continued until the particle either collides with the vane or leaves the domain via the outlet. If the particle passes through the north periodic boundary, its position is communicated to the south periodic boundary and vice versa.

The drag force is the most dominant in the transport of particles $dp>1\mu m$ in this context (see Sec. 2.4). Turbulent dispersion is the dominant transport mechanism for submicron particles. Both of these forces have been included in the analysis along with the pressure gradient force given the nonzero pressure gradient in the continuous phase. Additional forces such as thermophoresis and gravity were neglected as these dominate for particles with diameters less than 0.1 *μ*m [22] and greater than 25 *μ*m, which are not relevant in the current context. The effects of the virtual mass force were neglected due to the large ratio between the dispersed and continuous phase densities.

Turbulent dispersion of particles due to local velocity fluctuations in the boundary layer was modeled using the stochastic discrete random walk model. Isotropic turbulence levels are assumed throughout the domain, and the instantaneous particle velocity is computed from the sum of the local Reynolds-averaged velocity and the turbulence-induced fluctuation. Due to the random nature of the turbulent velocity fluctuations, two identical particles injected at the same location may follow different paths. Therefore, an increased number of particles need to be injected and tracked from each injection point for the results to be statistically meaningful. To cater for this, 25 identical particles are released from each injection point and used to compute a mean Lagrangian trajectory.

### 4.5 Interaction Probability.

To record the interaction probability or analogously the capture efficiency of the NGV, a condition was put on to the particle trajectory computation to discontinue if the particle position coincided with the vane surface. These are effectively *captured* by the NGV. Particles that do not intercept the vane depart the domain. Particles are not permitted to bounce off vane surfaces, which may lead to an underprediction of deposition, especially on the suction surface which is otherwise shielded from dust exposure. However, since the outcome of the interaction is so heavily dependent on the particle properties and impact properties, it would be futile to include bounce without the inclusion of an accurate rebound model for all expected particles and their states. Such models are still in development, although Singh and Tafti [23] have developed and demonstrated good agreement viscosity-based sticking models for noncrystalline materials such as coal ash. In the present work, it should be remembered that the purpose of the interaction probability is to predict more accurately the starting composition reaching the vane and understand how it varies from the composition of the material entering the stage.

### 4.6 Capture Efficiency.

*capture fate*of each particle injected into the domain will be a function of the initial

*y*-coordinate of the particle and its velocity response time relative to that of the flow. Since the flow response time is only a

*characteristic*of the bulk flow, if each particle injected has the same response time, each particle will also have the same generalized Stokes number. Yet, we know that some of these particles will not interact with the NGV; two particles of the same Stokes number or two particles of the same injection position but different Stokes number can have entirely different fates, as illustrated in Fig. 3. Thus, for a given Stokes number, the NGV has a capture

*efficiency*:

*n*

_{captured}is the number of particles that come into contact with the NGV, and

*n*

_{injected}is the number of particles injected into the domain. Other authors have adopted the same definition, e.g., Refs. [23,21]. The capture efficiency is a property of the NGV and depends on the Stokes number. Hence by generating a range of Stokes numbers, we can investigate the capture efficiency for a range of particle and flow properties.

### 4.7 Test Matrix.

We utilized the three test points from the Reynolds number experiment in Ref. [20] to obtain the boundary conditions for the continuous phase in the numerical simulations. The Reynolds numbers in that experiment were based on characteristic velocity at the throat, assumed to be the minimum throat gap width. The domain boundary conditions for the three Reynolds number chosen (the design point plus two off-design points) are presented in Table 3. The main difference between cases is the flow Reynolds number, which is varied in the experiment by modulating the inlet pressure at design point values of pressure ratio and speed to change flow density. As can be noted from Table 3, the throat velocity does not change greatly. However, the change in flow density would affect the particle trajectory owing to differences in the drag. It is worth pointing out that we have borrowed test points in which the turbine efficiency is independent of Reynolds number, see Ref. [20].

Re_{th} | – | 123,520 | 141,426 | 176,480 |
---|---|---|---|---|

P_{S,outl} | Pa | 144,749 | 165,371 | 205,858 |

T_{T,outl} | K | 711.1 | 711.0 | 710.7 |

W_{1,outl} | kg s^{−1} | 7.59 | 8.66 | 10.84 |

Re_{th} | – | 123,520 | 141,426 | 176,480 |
---|---|---|---|---|

P_{S,outl} | Pa | 144,749 | 165,371 | 205,858 |

T_{T,outl} | K | 711.1 | 711.0 | 710.7 |

W_{1,outl} | kg s^{−1} | 7.59 | 8.66 | 10.84 |

*particle*Reynolds numbers to inject into each continuous phase simulation test point. We selected 15 “target” particle Reynolds numbers based on throat velocity between Re

_{p}= 10–1000 and calculated the particle diameter required to achieve this at each

*throat*Reynolds number Re

_{th}as fixed by the continuous phase test matrix. Particle diameter is then calculated from

_{th}refers to a value at the minimum throat width, drawn perpendicular to the suction surface of the vane. We performed four discrete phase tests to investigate

the effect of throat Reynolds number Re

_{th}the effect of section radial position

the effect of particle sphericity $\varphi $

the choice of flow velocity response time $\tau fv$

on the suitability of a reduced-order function. All of the diameters were injected within the applicability limits for neglecting thermophoresis and gravity forces. A check on Knudsen number for each case confirms that the no-slip assumption holds for each case. The baseline case was the “Mid” section geometry with an inlet Reynolds number of 176,500, using spherical particles. This inlet Reynolds number is the GE-E^{3} design point. The test matrix is shown in Fig. 11 for the inlet Reynolds number test.

It is worth noting that there is not a great deal of difference in the computed particle sizes between each flow Reynolds number. The difference appears to increase with increasing flow Reynolds number. This is reflected in the test matrices of the hub and tip, which points to the fact that each case has a similar throat velocity and throat width. In other vane geometries, this may not be the case.

In each of the above cases, all particles were prescribed the same density, equal to 2551 kg m^{−3}. The complete set of 15 runs at the midsection was repeated for particle sphericities of 1.00 (spherical), 0.75, and 0.54, the two latter of which loosely resemble typical shapes of quartz and gypsum, respectively [19]. The effect of sphericity on the particle drag coefficient was modeled by implementing Eq. (12) into the Lagrangian particle tracking algorithm. This allowed us to test the application of the generalized Stokes number for nonspherical particles.

Recalling the objectives given in Sec. 4, we propose a generalized Stokes number by which to reduce the number of dependent variables $\eta NGV$. We adopt the form presented in Eq. (1) and adopt the correction factor $\psi $ given in Eq. (2). The particle velocity response time $\tau pv$ is determined by the particle diameter, particle density, and fluid viscosity, which is referenced at the throat. To investigate the influence of flow velocity response time $\tau fv$, we produce three sets of characteristic velocity *U* and characteristic length *L*, which are stated in Table 4.

Re_{th} | – | 123,520 | 141,426 | 176,480 |
---|---|---|---|---|

$\tau fv,01$ | ||||

$(LthUin)hub$ | s | $0.12686.3$ | $0.12686.8$ | $0.12687.0$ |

$(LthUin)mid$ | s | $0.12580.2$ | $0.12580.7$ | $0.12581.0$ |

$(LthUin)tip$ | s | $0.12374.7$ | $0.12375.1$ | $0.12375.5$ |

$\tau fv,02$ | ||||

$(LthUth)hub$ | s | $0.126360.7$ | $0.126374.8$ | $0.126401.4$ |

$(LthUth)mid$ | s | $0.125361.6$ | $0.125375.8$ | $0.125401.1$ |

$(LthUth)tip$ | s | $0.123368.5$ | $0.123383.4$ | $0.123409.0$ |

$\tau fv,03$ | ||||

$(LLEUin)hub$ | s | $0.13686.3$ | $0.13686.8$ | $0.136401.4$ |

$(LLEUin)mid$ | s | $0.13680.2$ | $0.13680.7$ | $0.136401.1$ |

$(LLEUin)tip$ | s | $0.14274.7$ | $0.14275.1$ | $0.142409.0$ |

Re_{th} | – | 123,520 | 141,426 | 176,480 |
---|---|---|---|---|

$\tau fv,01$ | ||||

$(LthUin)hub$ | s | $0.12686.3$ | $0.12686.8$ | $0.12687.0$ |

$(LthUin)mid$ | s | $0.12580.2$ | $0.12580.7$ | $0.12581.0$ |

$(LthUin)tip$ | s | $0.12374.7$ | $0.12375.1$ | $0.12375.5$ |

$\tau fv,02$ | ||||

$(LthUth)hub$ | s | $0.126360.7$ | $0.126374.8$ | $0.126401.4$ |

$(LthUth)mid$ | s | $0.125361.6$ | $0.125375.8$ | $0.125401.1$ |

$(LthUth)tip$ | s | $0.123368.5$ | $0.123383.4$ | $0.123409.0$ |

$\tau fv,03$ | ||||

$(LLEUin)hub$ | s | $0.13686.3$ | $0.13686.8$ | $0.136401.4$ |

$(LLEUin)mid$ | s | $0.13680.2$ | $0.13680.7$ | $0.136401.1$ |

$(LLEUin)tip$ | s | $0.14274.7$ | $0.14275.1$ | $0.142409.0$ |

## 5 Results and Discussion

The capture efficiency was recorded for every test point, which is presented in Fig. 12. The range of particle Reynolds numbers, characterized at the throat, produce a good variation in capture efficiency from almost zero to one. The range of capture efficiencies produced at the same particle Reynolds number illustrates the influence of additional factors, such as, for example, particle sphericity for Re_{p} = 100, $\u223c70%$ of spherical particles have collided with the vane, whereas $\u223c48%$ of particles with a sphericity of 0.75 (e.g., Quartz) and just $\u223c30%$ of particles with $\varphi =0.54$ (e.g., Gypsum) have collided with the vane. This is due to the difference in the drag curves for each sphericity, as shown in Fig. 2. Looking at Fig. 5, the inertial effects appear at lower Reynolds number for nonspherical particles (Re_{p} > 30 for $\varphi =0.54$ and Re_{p} > 70 for $\varphi =0.54$) which explains the much reduced capture efficiency of these particles at Re_{p} = 100. This demonstrates that the reference dimensions in the definition of particle Reynolds number for this numerical experiment are at least appropriate.

### 5.1 Effect of Throat Reynolds Number, Radial Position, and Particle Shape.

To investigate the suitability of a generalized Stokes number for NGVs, the abscissae data shown in Fig. 12 are re-cast using Eq. (1), using the throat gap as a characteristic length and inlet velocity as a characteristic velocity. The results are plotted in Fig. 13. The reference particle diameter is given by Eq. (15). It can be seen that for spherical particles, the data points collapse on to a single trendline, which is independent of the throat Reynolds number and span-wise position. The effect of vane section radial position in Fig. 13 (circle, cross, star) results in a maximum difference of around 15% in capture efficiency at the hub compared with the tip at the same throat Reynolds number and particle Reynolds number, most likely due to the smaller inlet area rather than vane geometry. This demonstrates the Stokes number’s ability to “correct” span-wise effects and may also broaden its applicability for three-dimensional vanes.

Similarly, the effect of changes to operating mass flow rate (i.e., throat Reynolds number—see Table 3) is also catered for by the Stokes number (circle, plus, diamond symbols). For example, we see a large difference in capture efficiency by the midspan geometry for a particle with a Reynolds number of 75 across the three throat Reynolds numbers: 43%, 71%, and 93% for Re_{th} = 123,520, 141,426, and 176,480, respectively. The particle diameters used in each case were 7.6 *μ*m, 6.6 *μ*m, and 5.3 *μ*m, respectively. It is interesting to note that the capture efficiency is greater for the case with the lower throat Reynolds number, despite the throat velocity in this case being the smallest. Since particle Reynolds number is designed to be the same for each particle at the throat, we attribute this difference to the large difference in flow pressure in the rest of the domain (see Table 3). This would result in a larger flow density hence greater drag force at the same drag coefficient and relative velocity, which would decrease the particle response time, therefore leading to few particle-vane interactions. This highlights the importance of flow density in achieving dynamic similarity in rig scale tests and illustrates how generalized Stokes number can be used as a matching parameter to “correct” this effect.

The generalized NGV Stokes number is less useful for correcting the effect of particle shape. Despite identical particle Reynolds number and flow conditions, the capture efficiency at a given Stokes number is shown to vary greatly—as much as 50% at a Stokes number of around 1 when the particle sphericity is set to 0.54. This is explained by the larger surface area, hence drag per unit mass of a nonspherical particle. Interestingly, the test points for the $\varphi =0.75$ particles show a similar shape to the spherical particles, while the $\varphi =0.54$ particles appear to exhibit a log-linear relationship of capture efficiency with Stokes number above *Stk*_{NGV} = 1. This can be explained by the lower transition Reynolds number of the latter case, as can be seen in Figs. 2 and 5, which has the effect of creating more drag per unit mass thanks to the elevated drag coefficient at the same particle Reynolds number. This effect is less visible in the $\varphi =0.75$ case, probably as a result of the particle Reynolds number along the trajectory not exceeding the transition Reynolds number of that shape. This may be important when scaling for dynamic similarity—for example, the *effective* Stokes number of a quartz particle will be lower than that calculated if the particle is assumed spherical. The consequence of this would be a larger than the expected mean particle diameter of dust interacting with the vane and a reduction in the overall vane dose.

### 5.2 Interaction Probability.

The importance of the drag coefficient on particle-vane interaction is mentioned in the work of Sacco et al. [7], in which the authors discuss the concept of the “effective” Stokes number to account for non-Stokesian drag. Nevertheless, the Stokes number is defined and used without correction. In the present work, we adopt the same definition of Stokes number and apply the correction factor in order to compare it with other potential Stokes number definitions. Rather encouragingly, there appears to be a very similar correlation between Stokes number and capture efficiency when using either throat gap or vane leading edge diameter in the flow response time, irrespective of mass flow rate or radial position. The curve fit shown in Fig. 14 is through the points defined by the throat gap; the points defined by leading edge diameter have a slightly worse fit to the curve. When throat velocity is used instead of inlet velocity, the curve fit retains its shape but is shifted to the right and exhibits a slightly poorer fit to that of the flow response time defined by throat gap.

*capture efficiency*is a term borrowed from previous work on aerosol capture by solid cylinders, it is not applicable in the present case, where the probability of a particle being retained by a vane after impact is determined by a whole new set of variables relating to its mineral phase, yield strength, Young’s modulus, and impact energy, in addition to properties of the surface itself such as temperature, roughness, material, and so on. Hence, we adopt the term

*interaction probability*$\eta NGV$ to describe the likelihood of a particle coming into contact with a vane surface. It is a number distribution expressed as a function of generalized Stokes number

*Stk*

_{NGV}generated from the curve fits shown in Fig. 14. We recommend the following form:

*a*,

*b*,

*c*, and

*d*are dependent on the flow response time chosen (subject to availability of flow data). Curve fit parameters for the three response times definitions are given in Table 5. We also suggest coefficients for a sphericity of $\varphi =0.75$ and 0.54, although this is to demonstrate the influence of shape on overall interaction efficiency of an NGV for a given distribution of particle size. The Stokes number is of the form presented in Eq. (2), with characteristic values given in Table 4 for the definition of flow response time and particle Reynolds number. The curve fits are applicable for Re

_{p}between 10 and 1000. The curve fits for $\varphi =1.00$ and $\varphi =0.75$ display a high goodness, but the curve for $\varphi =0.54$ is not good, most likely due to the transition point occurring at a lower Reynolds number (see Fig. 2). This may point to a limit of applicability of the other curve fits: should flow path accelerations lead to particle Reynolds numbers well in excess of the transition point, a new curve fit may be required.

a | b | c | d | |
---|---|---|---|---|

$\tau fv,1,\varphi =1.00$ | −0.961 | −2.7 | 1.6 | 0.039 |

$\tau fv,2,\varphi =1.00$ | 1.010 | 0.51724 | 0.6082 | −0.010 |

$\tau fv,3,\varphi =1.00$ | −0.9610 | −0.2 | 1.6 | 0.039 |

$\tau fv,1,\varphi =0.75$ | 0.9729 | 1.17775 | 1.2568 | 0.0278 |

$\tau fv,1,\varphi =0.54$ | 1.010 | 0.51724 | 0.6082 | −0.010 |

a | b | c | d | |
---|---|---|---|---|

$\tau fv,1,\varphi =1.00$ | −0.961 | −2.7 | 1.6 | 0.039 |

$\tau fv,2,\varphi =1.00$ | 1.010 | 0.51724 | 0.6082 | −0.010 |

$\tau fv,3,\varphi =1.00$ | −0.9610 | −0.2 | 1.6 | 0.039 |

$\tau fv,1,\varphi =0.75$ | 0.9729 | 1.17775 | 1.2568 | 0.0278 |

$\tau fv,1,\varphi =0.54$ | 1.010 | 0.51724 | 0.6082 | −0.010 |

### 5.3 Validation.

To the authors’ knowledge, there are no experiments currently in the literature that can be used to fully validate the interaction probability. The most appropriate existing experimental data that could help us to validate our results are those generated by researchers at Ohio State University working on representative nozzle guide vanes using the *turbine reacting flow rig* (TuRFR). The experiments that are appropriate to the present study report an accumulation factor which, as we describe above, is a product of interaction probability and retention probability. Since we do not know the retention probability, we cannot yet isolate the interaction probability. Nevertheless, we may examine these works to afford some validity to the interaction probability concept.

There are no experiments that subject the NGV row to a single particle size, so we are unable to compare across the range of Stokes numbers studied. However, if the flow Reynolds number and Stokes number *are* reported, they can be used to calculate the non-Stokes drag correction factor and the interaction probability according to Eqs. (11) and (16), respectively. Experimental data from a recent study by Whitaker et al. [24] are used for this purpose. Germane parameters are collected in Table 6. The particle sphericity is not reported, but in an earlier study by the same authors, the shape is reported to be irregular [25]. For this reason, we have calculated $\eta NGV$ as if the sphericity were 0.75 and compared with the capture efficiency quoted in the experiments. In the two cases conducted at a turbulence level similar to this study, the interaction probability exceeds the accumulation factor. This partly validates the new model and suggests that of the material that interacts with the vane, around 18% of the 6.5 *μ*m particles and 15% of the 4.6 *μ*m particles would deposit.

### 5.4 Application.

*μ*

_{geo}is the geometric mean and $\sigma geo$ is the geometric standard deviation. The subscript [.]

_{3}denotes that this is a frequency distribution

*by mass*; size distribution functions may also be expressed

*by number*. The two minerals have different bulk densities and weightings. In real world operations, the ingested dust or ash may contain a number of minerals, with a range of densities, shapes, textures as well as sizes; understanding all the processes that take place along the gas path—such as fragmentation, phase change, reactions, agglomeration—before reaching the NGV, is an open area of research.

Mineral: | A | B | A + B | |
---|---|---|---|---|

μ_{geo} | μm | 1.0 | 2.0 | – |

$\sigma geo$ | μm | 0.4 | 0.3 | – |

$\rho p,bulk$ | kg m^{−3} | 2500 | 3000 | – |

Weighting | (–) | 40% | 60% | 100% |

$\chi \varphi =1.00$ | (–) | 6.0% | 48.0% | 54.0% |

$\chi \varphi =0.75$ | (–) | 4.7% | 36.1% | 40.8% |

Mineral: | A | B | A + B | |
---|---|---|---|---|

μ_{geo} | μm | 1.0 | 2.0 | – |

$\sigma geo$ | μm | 0.4 | 0.3 | – |

$\rho p,bulk$ | kg m^{−3} | 2500 | 3000 | – |

Weighting | (–) | 40% | 60% | 100% |

$\chi \varphi =1.00$ | (–) | 6.0% | 48.0% | 54.0% |

$\chi \varphi =0.75$ | (–) | 4.7% | 36.1% | 40.8% |

*re-cast*as the generalized NGV Stokes number

*Stk*

_{gen}for the given operating conditions. To facilitate this,

*Stk*

_{gen}can be expressed in terms of known engine parameters:

*N*is the number of vanes,

*h*is the vane height, and

*W*

_{1}is the engine mass flow rate. The particle Reynolds number is given by rearrangement of Eq. (4), and the non-Stokes drag correction factor $\psi $ is a function of Re

_{p}as given by Eq. (2). For a given size distribution, the corresponding Stokes number distribution is then calculated, and subsequently, the interaction probability using Eq. (16). The product of

*n*

_{3}and $\eta NGV$ is the frequency distribution by the mass of the interacting particulate, which is illustrated in Fig. 15 for the blend of two minerals given in Table 7. The darker shaded area represents the portion of the original distribution that interacts with the NGV. Integrating this function between the upper and lower bounds of a given dust mass distribution gives the overall

*interaction factor by mass*$\chi 3,NGV$ for the NGV:

## 6 Conclusions

The purpose of the current work was to investigate the factors affecting the particulate trajectories in a nozzle guide vane array, in order to reduce the number of variables upon which the accumulation rate on a gas turbine nozzle guide vane depend.

The probability of spherical particle interaction with a nozzle guide vane can be predicted as a function of Stokes number alone. The *generalized NGV Stokes number**Stk*_{NGV} is defined by the product of a non-Stokesian drag correction factor and the ratio of characteristic particle to flow response times. The use of this number reduces the dependence of the interaction probability to one variable: Stokes number, thereby allowing for a quick estimation of the amount and composition of dust that will interact and potentially adhere to the vane and ultimately a better prediction of accumulation factor $\zeta $.

Particle shape should not be ignored when scaling for dynamic similarity. The effect of nonsphericity is to increase the drag per unit mass, which leads to a greater tendency of particles to follow fluid streamlines, thereby avoiding interaction. For a particle sphericity of 0.75, typical of a quartz particle, the reduction in interaction probability can be as much as 25%; for a gypsum particle of sphericity 0.54, the reduction is up to 50% at the same generalized NGV Stokes number. The consequence of assuming spherical particles in the calculation of Stokes number would be a larger than the expected mean particle diameter of dust interacting with the vane and a reduction in the overall vane dose. This would probably affect the rate of deposition.

Interaction probability can be predicted via a reduced-order function. Using the generalized Stokes number as a single-dependent variable, the interaction probability can be estimated for any given operating condition and any NGV geometry using a curve fit function. Three separate curve fits have been generated for sphericities of 1, 0.75, and 0.54. The characteristic flow response time used in the calculation of Stokes number can be chosen to suit the available data; if NGV throat velocity is used in the absence of inlet velocity, the interaction probability curve retains its shape but shifts to the right.

Interaction factor is the mass proportion of dust interacting with the vane and can be calculated for a given dust. A multimineral dust arriving at the NGV stage will possess a range of particle sizes. Each mineral’s size distribution can be represented by a probability density function such as a log-normal distribution. Using this size data, the generalized Stokes number can be calculated for an engine given operating point, and the overall *interaction factor*$\chi NGV$ estimated over a particle size range. This could be used to better control or predict the dust dose that leads to deposition of mineral dusts or volcanic ash on a gas turbine nozzle guide vane.

## Nomenclature

*c*=concentration

*d*=diameter

*f*=drag factor

*m*=mass

*n*=frequency distribution

*L*=characteristic length

*U*=characteristic velocity

*C*_{D}=drag coefficient

*P*_{S}=static pressure

*T*_{T}=total temperature

*W*_{1}=engine core air mass flow rate

- [.]
_{0}= initial

- [.]
_{3}= by mass

- [.]
_{bulk}= of the bulk dust

- [.]
_{exp.}= experimental

- [.]
_{f}= of the fluid

- [.]
_{fv}= of the flow velocity

- [.]
_{gen}= generalized

- [.]
_{hub}= of the hub section

- [.]
_{in}= of the domain inlet

- [.]
_{LE}= of the NGV leading edge

- [.]
_{mid}= of the midsection

- [.]
_{NGV}= of the nozzle guide vane

- [.]
_{out}= of the domain outlet

- [.]
_{p}= of the particle

- [.]
_{pred.}= predicted

- [.]
_{pv}= of the particle velocity

- [.]
_{th}= of the throat

- [.]
_{tip}= of the tip section

- [.]
_{vp}= of the particle volume

- Re =
Reynolds number

*Kn*=Knudsen number

*EP*=environmental particulate

*NGV*=nozzle guide vane

*Stk*=Stokes number

*Tu*=turbulence intensity

*ζ*=accumulation factor

*η*=capture efficiency or interaction probability

*μ*=viscosity

*μ*_{g}=geometric mean

*ρ*=density

*ϕ*=particle sphericity

*χ*=interaction factor

*ψ*=non-Stokes drag correction factor

## Acknowledgment

Part of this work was joint funded by the UK Engineering and Physical Sciences Research Council (Funder ID: 10.13039/501100000266) (EPSRC) and the Defence Science and Technology Laboratory (Funder ID: 10.13039/100010418) (DSTL), who are funding the doctoral studies of a co-author, Matthew Ellis, who performed the numerical simulations and post-processing for this work as part of EPSRC grant EP/P510579/1. The dataset generated during the current study is available in the Mendeley repository, http://dx.doi.org/10.17632/bzm2hwwt3s.