2.1 CFD model and boundary condition

Figure 2 shows the numerical model setup for the current research interest, including a ground plane and an isolated NACA-1-85-100-contraction1.25 nacelle inlet [Reference Re and Abeyounis28]. The inlet model has a thickness of 30.5 mm, and the outer diameter at the AIP (D) is 457.2 mm. The sizes of the computational domain are given in Fig. 2, and the descriptions of the symbols in the figure are given in Table 1. The origin of the coordinate system is located at the centre of the cross section at the lip of the nacelle inlet.

Figure 2.

Illustration of the CFD model of the nacelle inlet near the ground.

The boundary surface of the domain intersecting the cylindrical portion of the model and the outlet of the model were set as the ‘Pressure-Outlet’ boundary condition. The ground plane and all the inlet surfaces were set as the ‘No-Slip-Wall’ boundary condition. All the other surfaces of the domain were set as ‘Pressure-Far-Field’ boundary condition. The ambient pressure was fixed at 1 atm (101,325 Pa at sea level). The pressure at the outlet of the model was set to a negative pressure to form the internal flow in the inlet.

2.2 Fluid governing equations

The flow field of the fluid (air) phase needs to be solved first, of which the results will serve as the input of the subsequent particle tracking calculation. The steady compressible RANS equation was solved by the pressure-based solver in the commercial package ANSYS Fluent 2022 R1 [Reference Fluent29] to obtain the airflow field of the inlet near the ground. The governing equations of the fluid phase are expressed as follows.

Continuity equation:

(1) \begin{align}\nabla \cdot ( {\rho \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over u} }) = 0\end{align}

Momentum equation:

(2) \begin{align}\nabla \cdot ({\rho \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over u} \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over u} }) = - \nabla p + \mu \nabla \cdot \left[ {\left( {\nabla \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over u} + \nabla {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over u} }^T}} \right) - \frac{2}{3}\nabla \cdot \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over u} I} \right]\end{align}

where $\rho $ is the fluid density, $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over u} $ is the fluid velocity, $\mu $ is the fluid dynamic viscosity and $I$ is the unit vector.

The shear stress transport (SST) $k - \omega $ model [Reference Menter30] was chosen to simulate the turbulence of the air flow field. The SIMPLEC scheme was applied to solve the pressure-velocity coupling equation. The gradient term was discretised by the Least-Square-Cell-based scheme, the pressure term by a second-order interpolation scheme and all the other spatial terms of the RANS and turbulence equations by the second-order upwind scheme.

Table 1.

Definitions and dimensions of the model geometry of the nacelle inlet

Table 1. Long description

A table describing the dimensions and positions of various parts of a nacelle inlet model. The table has 7 rows and 2 columns. The first column is labeled Description and the second column is labeled Symbol. The rows are as follows: Row 1: Outer diameter of the AIP, D = 457.2 (mm). Row 2: Length of the cowl from the lip (highlight) to the AIP, L = 1D. Row 3: Length of the cylindrical portion, L_cy = 2.5D. Row 4: Position of the lip, a. Row 5: Position of the throat, b. Row 6: Position of the AIP (i.e. start of the cylindrical portion), c. Row 7: Position of the exit of the external cowl surface, d. Row 8: Position of the exit of the cylindrical portion, e.

2.3 Mesh and independence examination

A structured-mesh strategy was adopted to discretise the computational domain in ANSYS ICEM [31]. The layout of the mesh for the nacelle inlet in ground effect is shown in Fig. 3. To confirm the above results are independent of the mesh configuration, three sets of meshes with coarse, fine and dense resolutions were generated, as shown in Table 2, and then sent to calculations at the condition of freestream Mach number of ${M_\infty } = 0.79$ , mass flow ratio of $MFR = 0.67$ and angle-of-attack $\alpha = 2.1^\circ $ [Reference Re and Abeyounis28]. It can be seen from Fig. 4 that the results of the wall static pressure coefficient calculated by the three meshes agree well with each other, except for a slight difference at the top lip of the cowl. By careful inspection of the flow field, it is found that the coarse mesh produces a significant pressure difference from both the other two meshes on the top and bottom suction surfaces of the inlet, as circled in Fig. 5a. Overall, the difference between the results calculated by the fine and dense meshes is small and negligible. Hence, the fine mesh was finally chosen for all the subsequent calculations, which yields the maximum value of about 18 for the wall y⁺, meeting the near-wall turbulence modeling requirement of y ⁺ < 30.

Figure 3.

Computational mesh of the nacelle inlet close to the ground: (a) global view, (b) axial view at the inlet lip and (c) lateral view at the inlet throat.

Table 2.

Parametric settings of the three leveled meshes for the independence examination

Table 2. Long description

A table comparing the parameters of three different mesh resolutions used in a computational domain. The table has 8 rows and 3 columns. The columns are labeled Coarse mesh, Fine mesh, and Dense mesh. The rows are labeled with different variables: Number of nodes in the boundary layer, Number of nodes in the inlet circumference, Number of nodes from the lip to the AIP, Number of nodes from the AIP to the outlet, Height of the cell attached to the inlet wall (mm), Total number of nodes, and Total number of elements. The values in the table are as follows: Row 1: Number of nodes in the boundary layer, Coarse mesh: 35, Fine mesh: 50, Dense mesh: 50. Row 2: Number of nodes in the inlet circumference, Coarse mesh: 84, Fine mesh: 120, Dense mesh: 176. Row 3: Number of nodes from the lip to the AIP, Coarse mesh: 50, Fine mesh: 68, Dense mesh: 80. Row 4: Number of nodes from the AIP to the outlet, Coarse mesh: 80, Fine mesh: 103, Dense mesh: 110. Row 5: Height of the cell attached to the inlet wall (mm), Coarse mesh: 0.05, Fine mesh: 0.02, Dense mesh: 0.02. Row 6: Total number of nodes, Coarse mesh: 2,198,406, Fine mesh: 4,082,484, Dense mesh: 6,702,399. Row 7: Total number of elements, Coarse mesh: 2,253,040, Fine mesh: 4,169,214, Dense mesh: 6,829,276.

Figure 4.

Comparison of the inlet wall static pressure coefficient results calculated by the three resolved meshes: (a) Φ = 0°, (b) Φ = 180°.

Figure 5.

Results of pressure distribution calculated by (a) coarse mesh, (b) fine mesh and (c) dense mesh.

2.4 Fluid solver validation

To validate the accuracy of the fluid solver, a testing condition of freestream Mach number of ${M_\infty } = 0.79$ , mass flow ratio of $MFR = 0.67$ and angles of attack of $\alpha = 0^\circ $ and 2.1° was chosen to calculate the flow field of the inlet. The calculated results of the wall static pressure distribution are compared with the experimental measurements by Re and Abeyounis [Reference Re and Abeyounis28], as shown in Fig. 6. There is a slight difference between both kinds of results on the suction surfaces of the cowl lip at both angles of attack, which may come from the interference of the three struts used in the experiment to instrument the pressure probes [Reference Re and Abeyounis28]. In general, the numerical results are in good agreement with the experimental measurements at both angles of attack.

Figure 6.

Comparison of the inlet wall static pressure coefficient between the experimental measurements and current results at (a) α = 0°, (b) α = 2.1°.

Figure 6. Long description

The image contains two line graphs labeled (a) and (b), comparing experimental measurements and computational fluid dynamics (CFD) results of the inlet wall static pressure coefficient at different angles. Panel A: The graph shows the static pressure coefficient (C_p) on the y-axis ranging from -2.0 to 1.5 and the normalized distance along the inlet wall (X/L) on the x-axis ranging from 0 to 100. The data points for experimental measurements at angles 0 degrees and 180 degrees are represented by black circles and red triangles, respectively. The solid black line represents CFD results at 0 degrees, and the dashed red line represents CFD results at 180 degrees. The graph indicates that the static pressure coefficient varies significantly along the inlet wall, with notable differences between experimental and CFD results at different angles. Panel B: This graph also shows the static pressure coefficient (C_p) on the y-axis ranging from -2.0 to 1.5 and the normalized distance along the inlet wall (X/L) on the x-axis ranging from 0 to 100. Similar to Panel A, the data points for experimental measurements at angles 0 degrees and 180 degrees are represented by black circles and red triangles, respectively. The solid black line represents CFD results at 0 degrees, and the dashed red line represents CFD results at 180 degrees. The graph shows a similar trend to Panel A, with variations in the static pressure coefficient along the inlet wall and differences between experimental and CFD results at different angles.

3.1 Sampled particle size composition

As for the later sand ingestion study, the parameters related to particle shape distribution, drag model and particle-wall rebound model in the current particle tracking method were all based on the sand sample meeting the requirements of China GJB 1171-91 standard. The requirements of composition and content, as well as the particle size distribution, are listed in Tables 3 and 4, respectively.

3.2 Particle shape distribution probability

For a real sediment bed, there are various shapes of the sand particles. In this study, we abstracted out the cross-sectional shapes of the sampled particles to be sphere, triangular prism, quadrangular prism and pentagonal prism by scanning electron microscope. Note that each of the shapes may have different aspect ratios and circularities. Therefore, there exists a distribution probability of each particle shape, which is calculated by statistics of the respective proportion of each shape from a number of random subsamples. In DPMTrace, a probability method was used to determine the distribution probability of the particle shapes, of which the principle is depicted as follows:

1. (1) Assume that an event has an occurrence probability of 30%. A random integer a is taken for the event and divided by 100 to find the remainder b. If 0 ≤ b ≤ 29, the event occurs; otherwise, it does not.

2. (2) Similarly, assume that Event A has an occurrence probability of 20%, as well as Event B of 30% and Event C of 50%. The parameters of a and b are obtained in the same way as above. If 0 ≤ b ≤ 19, Event A occurs. If 20 ≤ b ≤ 49, Event B occurs. If 50 ≤ b ≤ 99, Event C occurs.

3. (3) Assume that the probability of an event occurring is 30.5%. A random integer a is taken for the event and divided by 1000 to find the remainder b. If 0 ≤ b ≤ 304, the event occurs. In other words, if the probability of the event occurring in the percent sign is not an integer, but with n decimal places, the remainder to simulate the probability is obtained by dividing the random number by 10²⁺ⁿ .

Based on the above principle, the distribution probabilities of the particle shapes were calculated, and the simulated values are fairly close to the original values of the sampled sand particles, with the maximum relative error of about 0.56%. It is noted that the particle shape distribution probabilities used in DPMTrace were obtained from sand samples collected from a representative region, as detailed in Fig. 7a.

Table 3.

Particle composition and content

Table 3. Long description

A table showing the composition and mass percentage of various oxides in a sand sample. The table has two columns and six rows. The first column is labeled Composition and lists the oxides SiO2, Al2O3, CaO, Fe2O3, MgO, Na2O, and K2O. The second column is labeled Mass percentage and lists the corresponding mass percentages for each oxide. The values are as follows: SiO2: 68-70 percent, Al2O3: 10-11 percent, CaO: 5-6 percent, Fe2O3: 3-4 percent, MgO: 2-3 percent, Na2O: 2-3 percent, K2O: 1-2 percent.

Table 4.

Particle size distribution

Table 4. Long description

A table with two columns and six rows. The first column lists particle size ranges in micrometers (μm): 0-75, 75-125, 125-200, 200-400, 400-600, 600-900, and 900-1000. The second column lists the corresponding mass percentages with uncertainties: 12 ± 5 percent, 20 ± 3 percent, 26 ± 3 percent, 25 ± 3 percent, 12 ± 3 percent, 4 ± 1 percent, and 1 percent. The table shows the distribution of mass percentages across different particle size ranges.

3.3 Particle drag force calculation

The DPMTrace code is based on a one-way coupled Eulerian–Lagrangian particle tracking method, i.e. the effect of the particle existence on the continuous phase is not considered. The particle’s equation of motion in DPMTrace is based on Newton’s second law of motion and is written as follows:

(3) \begin{align}\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over a} = \frac{{d{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over u} }_p}}}{{dt}} = \frac{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over u} - {{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over u} }_p}}}{{{\tau _r}}} + \frac{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over g} ({{\rho _p} - \rho })}}{{{\rho _p}}} + \frac{\rho }{{{\rho _p}}}\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over u} \nabla \mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over u} \end{align}

where $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over a} $ is the particle acceleration, ${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over u} _p}$ is the particle velocity, ${\rho _p}$ is the particle density, $t$ is time and $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over g} $ is the gravity vector. The left-hand-side term is the acceleration, the first term on the right-hand side accounts for the acceleration due to the aerodynamic drag, the second term on the right-hand side is the resultant acceleration due to gravity and buoyancy, and the last term on the right-hand side is the pressure gradient force.

${\tau _r}$ is particle relaxation time calculated by Ref. [Reference Gosman and Loannides32]:

(4) \begin{align}{\tau _r} = \frac{{{\rho _p}{d_p}^2}}{{18\mu }}\frac{{24}}{{{C_d}{\mathop{\rm Re}\nolimits} }}\end{align}

where $\mu $ is the molecular viscosity of the fluid, ${d_p}$ is the particle diameter, and ${\mathop{\rm Re}\nolimits} $ is the relative Reynolds number defined as:

(5) \begin{align}{\mathop{\rm Re}\nolimits} \equiv \frac{{\rho {d_p}\left| {\vec u - {{\vec u}_p}} \right|}}{\mu }\end{align}

If the shape of the particle is spherical [Reference Liu, Mather and Reitz33]:

(6) \begin{align}{C_d} = \left\{ \begin{array}{l@{\quad}l} \dfrac{24}{\rm Re} \bigg(1 + \dfrac{1}{6}{\mathrm{Re}^{2/3}}\bigg), & {\mathrm{Re}} \le {\rm{1000}}\\[12pt]0.424, & \mathrm{Re} \gt 1000\end{array} \right.\end{align}

For non-spherical objects, the drag coefficient ${C_d}$ is a product of two components:

(7) \begin{align}{C_d} = {C_{d1}} \cdot {C_{d2}}\end{align}

where ${C_{d1}} = f\left( {{\mathop{\rm Re}\nolimits} ,{\theta _r}} \right)$ is the component of the drag coefficient as a function of the relative Reynolds number ${\mathop{\rm Re}\nolimits} $ and the relative angle between the fluid’s and particle’s velocity vectors ${\theta _r}$ . ${C_{d2}} = f\left( {{\mathop{\rm Re}\nolimits} ,h/d} \right)$ is the component of the drag coefficient as a function of the relative Reynolds number ${\mathop{\rm Re}\nolimits} $ and the stretching ratio $h/d$ , where $h$ is the three-dimensional stretching height and $d$ is the equivalent diameter of a round particle whose area is equal to that of the real particle’s cross-section. For a given sand sample, the fitted formulae of both ${C_{d1}}$ and ${C_{d2}}$ are derived by pursuing the minimum relative error between each of them and the total drag coefficient ${C_d}$ .

Finally, the drag coefficient of non-spherical objects is written as

(8) \begin{align}{{C}_{d}}={{C}_{d1}}\cdot {{C}_{d2}}=\frac{{{\text{p}}_{1}}+{{\text{p}}_{2}}{{\theta }_{r}}+{{\text{p}}_{3}}{{\theta }_{r}}^{2}+{{\text{p}}_{4}}\operatorname{Re}+{{\text{p}}_{5}}{{\operatorname{Re}}^{2}}}{\text{1+}{{\text{p}}_{6}}{{\theta }_{r}}+{{\text{p}}_{7}}{{\theta }_{r}}^{2}+{{\text{p}}_{8}}\operatorname{Re}+{{\text{p}}_{9}}{{\operatorname{Re}}^{2}}} \centerdot \frac{{{\text{p}}_{10}}+{{\text{p}}_{11}}(h/d)+{{\text{p}}_{12}}\operatorname{Re}}{\text{1+}{{\text{p}}_{13}}(h/d)+{{\text{p}}_{14}}\operatorname{Re}+{{\text{p}}_{15}}{{\operatorname{Re}}^{2}}}\end{align}

For the current sample, $h/d$ is taken randomly in the interval [0.2, 1], ${{\rm{p}}_1}$ to ${{\rm{p}}_{15}}$ are the 15 coefficients fitted by the experimental measurements of the drag coefficients for modelled particles of sphere, triangular prism, quadrangular prism and pentagonal prism. As shown in Fig. 7b–d, the irregular particle models were enlarged and tested in a water–glycerol mixture, from which the drag coefficients were measured and incorporated into the drag model.

3.4 Particle-wall rebound probability

The rebound characteristics of sand particles upon walls are dependent on many factors; thus, there is some uncertainty for a particle to achieve a certain rebound characteristic, i.e. rebound angle and rebound velocity. Figure 8 shows the schematic diagram of a particle impact and rebound on the wall. $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over v} $ is the incident velocity, ${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over v} _{\rm{r}}}$ the rebound velocity, $\left| {{{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over v} }_{\rm{r}}}/\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over v} } \right|$ the velocity restitution coefficient, ${\alpha _i}$ is the incident angle and ${\alpha _r}$ the rebound angle. The experimental procedure for the particle-wall rebound model used in DPMTrace is shown in Fig. 9.

Through observation of the physically tested probability curves of velocity restitution coefficient for incident angles ${\alpha _i} = 10^\circ \sim 90^\circ $ it is found that the velocity restitution coefficient assumes a normal distribution at any incident angle, as shown by Fig. 10a for an example of ${\alpha _i} = 40^\circ $ . Under the same condition, the mean value of the velocity restitution coefficients after repeated times of rebound tests of sampled particles will approach the mean value of the normal distribution. This characteristic can also explain the phenomenon that the mean velocity restitution coefficient of sand particles fluctuates basically steadily with the increase of incident velocity at any incident angle. Since the velocity restitution coefficient appears randomly between 0 and 1 at any incident angle, this is a phenomenon that satisfies the typical continuous probability distribution condition.

Figure 10.

(a) Probability of the velocity restitution coefficient and (b) probability of the rebound angle of the sampled particles at $\alpha = 40^\circ $ .

To fit the probability density function of the velocity restitution coefficient, the following probability density function for normal distributions was adopted:

(9) \begin{align}y = {y_0} + A{e^{ - \frac{{{{\left( {x - {x_c}} \right)}^2}}}{{2{w^2}}}}}\end{align}

where ${y_0}$ is the offset, $A$ is the amplitude, $w$ is the width and ${x_c}$ is the expectation value. These parameters were obtained by fitting the above equation with our experimental measurements of particle rebound characteristics using the GaussAmp mathematic function in the Origin software [Reference May and Stevenson34].

Through observation of the physically tested probability curves of rebound angle for incident angles ${\alpha _i} = 10^\circ \sim 90^\circ $ it is found that the rebound angle also assumes a normal distribution at any incident angle, as shown by Fig. 10b for example of ${\alpha _i} = 40^\circ $ . Since the rebound angle appears randomly between 0° and 180° at any incident angle, this is also a phenomenon that satisfies the typical continuous probability distribution condition. The probability density function of the rebound angle adopts the same form of probability density function as that of the velocity restitution coefficient.

The process of simulating the rebound probability densities in DPMTrace is basically depicted as follows:

1. (1) Assume the probability density function of a continuous variable $x$ is $f(x)$ , then the probability that $x$ appears in the interval $[{a,b}]$ is $\int_a^b {f(x)} dx$ .

2. (2) For particle-wall rebound, its rebound angle and velocity restitution coefficient are fixed values; that is, it needs to calculate the probability of occurrence of a certain rebound angle or restitution coefficient. DPMTrace uses the probability of an interval to replace the probability of occurrence of a certain value. For example, for the rebound angle, the interval [−1, 0] is used to replace the occurrence probability of a rebound angle of 0°, the interval [0, 1] is used to replace the occurrence probability of a rebound angle of 1°, and so on.

3. (3) For the probability density function $f(x)$ , $\int_{ - \infty }^{ + \infty } {f(x)} dx = 1$ . Because the velocity restitution coefficient and rebound angle both belong to a finite interval, i.e. [0, 1] and [0°, 180°], respectively, both integrations of probability density functions of the velocity restitution coefficient and rebound angle are smaller than 1.

4. (4) In order to simulate an event with the sum of probabilities of 1, the practice of normalisation processing is adopted. For example, assume that ${P_i}$ is the probability of a rebound angle $i^\circ $ , then the sum of probabilities of the rebound angle in the interval [0°, 180°] will be ${P_{sum}} = {P_0} + {P_i} + \cdots + {P_{180}}$ . Finally, ${P_i}/{P_{sum}}$ is used to represent the probability of the rebound angle of $i^\circ $ . Likewise, the probability of the velocity restitution coefficient is also normalised in the same manner.

5. (5) The probability of the velocity restitution coefficient and rebound angle of a certain value is simulated using the method described in Subsection 3.2.

Figure 10 shows a comparison of the velocity restitution coefficient and rebound angle probabilities between the fitted results of Equation (9) by GaussAmp in Origin and the simulation results by the currently developed code DPMTrace. It can be clearly seen that our code successfully simulates both rebound probabilities well, indicating good simulation accuracy of the rebound probability by DPMTrace.

3.5 Particle injection types

In DPMTrace, there are four particle injection types available to determine the initial distribution of particles, namely, monitoring plane, group, conical surface and cone, as shown in Fig. 11.

Figure 11.

Four types of particle injection models available in the present particle tracking method: (a) monitoring plane, (b) group, (c) conical surface and (d) cone. Images (e) and (f) are respectively the uniform and random particle distribution models determining the initial particle position.

The monitoring surface injection type directly employs the planes that have been generated in the phase of meshing and yields a random distribution of the initial position of the particles at one monitoring surface (Fig. 11a). In the schematic of the group injection type (Fig. 11b), ${{\rm{P}}_{\rm{1}}}$ is the starting position, ${{\rm{P}}_{\rm{2}}}$ is the ending position and ${V_{d0}}$ is the initial injection velocity of particles. The particles are uniformly distributed between ${{\rm{P}}_{\rm{1}}}$ and ${{\rm{P}}_{\rm{2}}}$ at the starting moment. In the schematic of the conical surface injection type (Fig. 11c), $O$ is the origin of the cone. $\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over n} $ is the axial direction of the conical surface. $R$ is the diameter of the conical surface. $\theta $ is the half-cone angle. The particles are uniformly distributed on the circumference of the bottom of the conical surface. In the schematic of the cone injection type (Fig. 11d), ${R_1}$ and ${R_2}$ are the outer and inner diameters of the cone, respectively. Besides, the particles can be initially distributed uniformly or randomly in the volume of the cone (Fig. 11e and f).

3.6 Particle tracking strategy

DPMTrace tracks the particle trajectory on a cell-by-cell level, as illustrated for a particle in Fig. 12. The tracking steps are described as follows:

1. (1) Assume that the particle with the initial velocity of ${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over v} _0}$ is initially located at the point ${\rm{P0}}$ in the body cell of Cell1.

2. (2) Taking ${\rm{P0}}$ as the base point and the vector ${\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over v} _0}$ as the direction, make a ray to intersect with the face element of Cell1 to get the point ${\rm{P1}}$ as the next position of the particle. The velocity of the particle at ${\rm{P1}}$ is:

   (10) \begin{align}{\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over v} _1} = {\mathord{\buildrel{\lower3pt\hbox{$\scriptscriptstyle\rightharpoonup$}} \over v} _0} + {a_1} \cdot {t_1}\end{align}
   where the acceleration ${a_1}$ is obtained from Equation (1) with the aerodynamic parameters taken from Cell1. Time ${t_1} = {{\rm{s}}_{\rm{0}}}/{v_0}$ , where ${{\rm{s}}_{\rm{0}}}$ is the distance between point ${\rm{P0}}$ and point ${\rm{P1}}$ .

3. (3) Similarly, all the subsequent positions, ${\rm{P2}}$ and ${\rm{P3}}$ , of the particle can be determined. Here, the connecting line of ${{P0\!-\!P1\!-\!P2\!-\!P3}}$ is formed and represents the trajectory of the particle.

Figure 12.

Tracking method of the trajectory of a single particle.

Figure 13.

Flowchart of the current particle calculation method.

If the particle hits the boundaries of the cells during the movement, different treatments need to be taken, depending on the different types of boundaries: (1) if the particle hits an inlet or outlet of the computational domain, then it stops moving; (2) if the particle hits a monitoring plane, it continues to move and the particle information (position, velocity, etc.) is recorded when the particle crosses the monitoring plane; (3) if the particle hits a wall, the velocity and angle of the particle after rebound are calculated according to the wall rebound model introduced in Subsection 2.4.

3.7 Particle calculation flowchart

The in-house code ‘DPMTrace’ is purely a particle tracking code, namely the air flow field simulation and particle motion simulation are completely separated. The flowchart of particle calculation in this study is shown in Fig. 13. Generally, there are three major steps corresponding to the three calculation modules:

1. (1) Mesh generation. This is the first step to input the mesh generated by the commercial meshing software ANSYS ICEM [31]. Note that the current Eulerian–Lagrangian sediment transport method is based on fully structured meshes.

2. (2) Fluid solution. This is the second step to input the fluid solution data calculated by ANSYS Fluent [Reference Fluent29]. The necessary fluid solution data include the three components of the three-dimensional flow field, the density and laminar viscosity of the fluid.

3. (3) Particle tracking. This is the core step to calculate the particle characteristics by the currently developed code DPMTrace. Herein, there are three sub-steps to follow sequentially, corresponding to the three modules in DPMTrace. The first module is file input, including the inputs of the mesh and fluid solution data, as well as the setting of the particle tracking results output path. Then it moves to the second module, i.e. parametric settings, including the settings of the monitoring planes, the control parameters, the particle parameters, as well as the drag and rebound parameters. The monitoring planes can be the boundary surfaces, like the outlet of the nacelle inlet or the interior planes that have been intentionally built for monitoring. The control parameters include the gravitational acceleration, the number of threads for parallel calculation and the interval of output of particle trajectory. The particle parameters include particle density, number of particles injected into the flow field, particle diameter distribution conforming to both the uniform and Rosin–Rammler (R–R) distributions, particle shape distribution probability and injection parameters. The drag and rebound parameters are included in the configuration data file under the installation folder. After setting all the above necessary parameters, it comes to the final module of calculation, i.e. calculation of the particles’ three-dimensional characteristics including the trajectory and concentration distribution, as well as the calculation of the separation efficiency for bifurcated flow channels, such as an inertial particle separator (IPS) [Reference Barone, Loth and Snyder35] for turboshaft engines and a turboprop engine inlet with bypass duct [Reference Zheng, Zhenlong, Huijun, Zhang and Kun36]. A local view of the DPMTrace interface is given in Fig. 14. It is noted that the developed DPMTrace software adopts a one-way coupled Eulerian–Lagrangian framework, in which the feedback effect of particles on the gas-phase flow field is neglected. This assumption is reasonable for the low particle-loading conditions considered in the present study, where particle-induced momentum exchange is negligible, and the gas-phase flow dominates particle transport. Therefore, the present method is mainly applicable to low- and moderate-concentration sand-ingestion conditions. Under extremely high particle-loading conditions, the momentum exchange between the particle and gas phases may become significant, and a two-way coupled framework would be more appropriate.

3.8 Particle tracking validation

To validate the DPMTrace code, the IPS case studied comprehensively by Prof. Loth’s research group [Reference Barone, Loth and Snyder37] is chosen and calculated, of which the CFD model is shown in Fig. 15. The definitions and dimensions of the model are given in Table 5. The leftmost boundary plane is the entrance of the model, and there are two ducts, i.e. the core duct and scavenge duct, downstream of the IPS. The red cross-section is a virtual plane for generating the mesh of the particle injection plane. As for the boundary conditions, the entrance of the model was set as ‘Pressure-Inlet’, as well as both the exits of the core and scavenge ducts set as ‘Pressure-Outlet’, and all the rest of the boundary planes set as ‘No-Slip-Wall’.

Figure 14.

Interface of the DPMTrace software.

Figure 14. Long description

A diagram of the DPMTrace software interface. The interface includes several tabs and windows for setting parameters and running simulations. Panel A shows the main interface with tabs for File, Parameters, and Calculate. The File tab includes options to read mesh, read flowfield result, and set output path. The Parameters tab includes options for monitoring plane setting, controlling parameter setting, particle parameter setting, drag rebound modification, and saving parameter settings. The Calculate tab includes options to start the calculation and check separation efficiency. Panel B shows the particle parameter setting window, which includes fields for density, number, size distribution, shape distribution probability, injection distribution, and monitoring plane settings. Panel C shows a simulation of particle flow with a 3D model and particle trajectories. The simulation includes labels for various dimensions and parameters such as height, length, diameter, and velocity. The diagram also includes text indicating the progress of reading the flowfield result file, saving parameter settings, starting and finishing calculations, and monitoring plane information.

Figure 15.

Illustration of the CFD model of Loth’s IPS.

A fully structured mesh was generated to discretise the computational domain of the IPS model, as shown in Fig. 16. It has a total of 2,262,429 nodes and 2,373,448 elements. The height of the cell attached to the walls was set to 0.02 mm, yielding the maximum value of about 12 for the wall y ⁺. The steady compressible RANS equation was solved by the pressure-based solver in the commercial package ANSYS Fluent 2019 [Reference Fluent29] to obtain the air flow field of the IPS. The ST $k - \omega $ model [Reference Menter30] was chosen to simulate the turbulence of the air flow field. The spatial discretisation schemes were set the same as at the nacelle inlet, as described in Subsection 2.2. For the sand particle simulation, the properties of particles were set to match the experimental condition [Reference Barone, Loth and Snyder35].

Figure 17 shows qualitative comparisons of the particle velocity fields of the IPS between Loth’s experimental measurements [Reference Barone, Loth and Snyder38] and our numerical results. Both the magnitude and direction of the fluid and particle velocities calculated by the current method are in good agreement with the referenced experimental results. To further quantitatively validate the DPMTrace code, the results of the particle separation efficiency $\eta $ of the IPS were calculated with both the Discrete Phase Model DPM in ANSYS Fluent [Reference Fluent29] and the currently developed DPMTrace code shown in Fig. 18, together with Prof. Loth’s experimental measurements [Reference Barone, Loth and Snyder35]. Although the trend of the results calculated by Fluent is generally consistent with that of the experimental measurements, the magnitude is overestimated at all the simulation points, with a maximum of about 11.7%. In contrast, both the magnitude and trend of the results calculated by DPMTrace assume a fairly good agreement with the experimental measurements. Overall, it can be concluded that the currently developed DPMTrace scheme has a good performance in particle dynamics predictions.