An Analysis of CFD-DEM with Coarse Graining for Turbulent Particle-Laden Jet Flows

1. Introduction

Two-phase dense particle-laden jets are found in many modern industrial applications, including jet mills, sandblasting, cold spray additive manufacturing, and a new innovative method for particle comminution called High-Pressure Slurry Ablation^©. Such systems can be costly, both in equipment and time, to monitor and evaluate experimentally. Traditional experimental techniques may involve many design iterations and consider several operating parameters, where valuable information for optimization is collected at the end of the experiment or test, and midstream phenomena are difficult to measure and quantify. Monitoring particle collision statistics and the fluid flow fields inside industrial units is difficult, further decreasing the amount of information available for optimization.

For simple, diluted two-phase air-particle and water-particle flows, there are a few methods that can be used to quantify the particle and fluid flow fields, such as Laser Doppler velocimetry (LDV), Particle Image Velocimetry (PIV), and Planar Nephelometry (PN) [1,2,3,4,5,6]. All of these methods are limited to simplified flow problems where it is required to have a predictable path or sampling location for the particles and transparent containers and carrier fluids. Furthermore, once an additional phase is added to a flow, such as in air-liquid-particle flows, light refraction occurs at the free surface making these methods ineffective [7]. In summary, the experimental methods presented here can be used to collect results to validate this current work; however, they cannot be used when opaque three-phase systems and industrial flow optimization problems are considered.

Fully coupled CFD-DEM solvers are well developed and validated across many different flow regimes [8,9,10]. Although many CFD-DEM studies exist for two-phase jet flows, most are for applications in wear in which they implemented a dilute flow regime that is uncoupled [11,12,13,14]. In a dilute flow regime, with reasonable accuracy, it is possible to perform an uncoupled numerical simulation where only fluid-particle interactions are considered. The uncoupled approach drastically increases the simplicity of the flow problem, primarily because of fewer limitations of mesh cell size and the absence of forces, and therefore models, acting on the fluid and particles that can bring forth both instabilities and more complexity to the flow problem. Furthermore, one can track particles as a post-carrier-phase solution process with an uncoupled simulation. Other two-phase CFD-DEM numerical solutions found in the literature are for combustion applications, where the particles are assumed to be droplets [15].

In our preliminary study [16], a long straight tube was considered with a fully developed particle-laden air jet at a Reynolds number of $15,000$ and mass loading up to $0.86$ that was simulated using fully coupled un-resolved CFD-DEM approach. Samples of particle flux at axial and radial locations along the jet were compared to the experimental data of Hardalupas et al. [1]. The primary consideration for this work was the effect of the particle volume fraction on the fluid flow field because of the concern for accurate momentum coupling between the two phases. It was concluded that a CFD cell minimum side length of $1.86$ particle diameters could be used without changing results. However, it was also realized that the round/plane jet anomaly caused significant issues for the single-phase results in context to the two-phase problem that would not otherwise exist if this were a pure single-phase flow problem. As explained at length in Wilcox’s well-known book Turbulence Modelling for CFD [17], using the k-$\epsilon$ turbulence model for standard and plane jets produces agreeable results. Still, a spreading rate as high as 40% over the measured experimental values can be seen for round jets. Many solutions have already been proposed to solve this issue, such as adding a vortex stretching term by Pope [18] and changing the empirical coefficients in the works of Faghani et al. [19], Givi and Ramos [20], Morgans et al. [21], but no nozzle domain was included in the flow profiles presented in these works.

Because this work is used as a foundation for the study of industrial flow systems with fully coupled two-phase CFD-DEM simulations, it is required that the nozzle profile shape be included in the solution domain. This is partly due to the required location and initial parameters for particle insertion - particles could be inserted at the nozzle exit. However, this would require experimental data for each change in nozzle shape, which would be impractical for an industrial flow problem. Instead, work was performed to slightly modify the empirical constants in the standard Launder and Spalding [22] k-$\epsilon$ turbulence model from $C_{\epsilon 1}=1.44$ and $C_{\epsilon 2}=1.92$ to $C_{\epsilon 1}=1.52$ and $C_{\epsilon 2}=1.90$. This approach circumvents the issue with spreading rates and has been used in the works by Faghani et al. [19] and also validated extensively in our previous work [23] by comparing the numerical results to the well-cited works of Bogusławski and Popiel [24] and Hardalupas et al. [1]. Therefore, in this current work, we use the modified k-$\epsilon$ model for our Hardalupas et al. [1] test cases and also perform an analysis comparing unmodified to modified k-$\epsilon$ for the Lau and Nathan [7] test cases for additional validation.

For CFD-DEM to be deemed an effective method, it is important to first ensure the accuracy of the models by comparing the numerical results to known experimental results. After validation is achieved for the simplified flows, more solver development into three-phase flows can be performed with the end goal of optimization of an industrial unit. Therefore, this work presents a validation study for CFD-DEM in the application of turbulent high-speed jets with a look into coarse-graining to help with the computational expense.

2. Numerical Methodology

The computational fluid dynamics (CFD) portion of the simulations is performed by discretization of the fluid domain into computational cells, where the incompressible Navier-Stokes equations are solved given the boundary conditions. The Navier-Stokes equations are highly non-linear partial differential equations and, therefore, require numerical models to be solved. This Eulerian framework for flow calculation is coupled with Newton’s equation of motion, using the velocity Verlet algorithm, for particle tracking in a Lagrangian framework of the discrete element method (DEM). The two frameworks exchange momentum as the simulation runs using averaging and interpolation methods.

2.1. CFD

The carrier phase is described by the multi-phase fully coupled Navier-Stokes Equations given by

$$\nabla \left({\epsilon }_f{\rho }_f{u}_f\right)=0$$

(1)

$$\frac{\partial \left({\epsilon }_f{u}_f\right)}{\partial t}+\nabla ·\left({\epsilon }_f{u}_f{u}_f\right)=-switch·\nabla p-K_{sl}\left(u_f-u_s\right)/\phantom{K_{sl}\left(u_f-u_s\right){\rho }_f}{\rho }_f+\nabla ·\left(switch·\tau \right)+f$$

(2)

where ${\epsilon }_f$ is the liquid phase volume fraction, defined as the ratio of the volume of fluid to the volume of the cell, $u_l$ is the liquid fluid velocity, $u_s$ is the mean solid particle velocity, $\tau$ is the liquid phase stress tensor, and f is an explicit force term that can be used to exchange momentum between the solid and carrier phase directly. The variable, $switch$, becomes either the void fraction, $switch={\epsilon }_f$, or one, $switch=1$, depending on whether Set I(Model Bfull or B) or Set II (Model A) is used for coupling [25,26] under the CFDEM${}^{\circledR }$ code, which will be explained further in Section 2.3[10]. The implicit momentum coupling term, $K_{sl}$, is given by

$$K_{sl}=\frac{{\epsilon }_f{\displaystyle \sum _i}{\overline{F}}_d}{V_c·\left|u_f-u_p\right|}$$

(3)

where $V_c$ is the volume of the cell and ${\displaystyle \sum _i}{\overline{F}}_d$ is the summation of all the forces acting on the fluid by the particles, which are defined in Section 2.3.

The liquid phase volume fraction, ${\epsilon }_f$, is calculated by determining the volume of particles inside of the cell. The volume of particles inside the cell can be calculated exactly, but this would be extremely computationally expensive. Therefore, the particles’ volume is calculated using a model that divides the particles into 29 non-overlapping regions of equal volume. The centroid of the particle is then used to determine the respective cell in which the new region exists. The approach allows accurate modeling while also providing decreased computational expense. Finally, the volumes of all regions inside a cell are added to obtain the total volume of particles, from which the liquid phase volume fraction can be calculated given the cell volume. There are errors associated with using a model for volume fraction calculation in an unresolved coupled CFD-DEM simulation. However, they are kept at a minimum with a sufficiently large cell size compared to particle size. Calculation of volume fraction is critical for accurate momentum coupling, so through our previous work, where similar simulations are performed for jet flow, it was realized that a minimum cell side length of $1.87$ particle diameters provided consistent results [16].

Turbulence is modeled using the "standard" k-$\epsilon$ model of Launder and Spalding [22] with slight modifications to the empirical coefficients. The turbulence model is given by

$$\frac{\partial \left(k\right)}{\partial t}+\frac{\partial \left(u_{j}k\right)}{\partial x_j}={\tau }_{ij}\frac{\partial u_j}{\partial x_j}-\epsilon +\frac{\partial }{\partial x_j}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_k}\right)\frac{\partial k}{\partial x_j}\right]+L_k$$

(4)

$$\frac{\partial \epsilon }{\partial t}+\frac{\partial \left(u_i\epsilon \right)}{\partial x_j}=C_{\epsilon 1}\frac{\epsilon }{k}{\tau }_{ij}\frac{\partial u_i}{\partial x_j}-C_{\epsilon 2}\frac{{\epsilon }^2}{k}+\frac{\partial }{\partial x_j}\left[\left(\nu +\frac{{\nu }_t}{{\sigma }_{\epsilon }}\right)\frac{\partial \epsilon }{\partial x_j}\right]+L_{\epsilon }$$

(5)

where eddy viscosity is given by ${\nu }_t=C_{\mu }{f}_{\mu }{k}^2/\phantom{k^2\epsilon }\epsilon$ and the Boussinesq assumption is used to obtain

$${\tau }_{ij}={\nu }_t\left(2S_{ij}-\frac{2}{3}\frac{\partial u_k}{\partial x_k}{\delta }_{ij}\right)-\frac{2}{3}k{\delta }_{ij}$$

(6)

$$S_{ij}=\frac{1}{2}\left(\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right)=symm(\nabla U)$$

(7)

The turbulent kinetic energy is given by k, the variable, $\epsilon$, is the turbulent dissipation rate, u is the fluid velocity, f is a damping function, and the turbulence length scale is given by $l=C_{\mu }{k}^{3/2}/\phantom{k^{3/2}\epsilon }\epsilon$[23]. The first term in Equation 6 is the deviatoric part and the second term is the isotropic part of the fluid stress tensor. The last term, k, is kinetic energy, which in OpenFOAM and most other CFD codes gets absorbed into pressure. Since we are exclusively considering incompressible flow, the following term becomes

$$-\frac{2}{3}\frac{\partial u_k}{\partial x_k}{\delta }_{ij}=0$$

(8)

but is still included on the CFD side for solution boundedness. It is worth noting that Equation 8 is not used to calculate the viscous particle force, as shown in Section 2.2. Instead, the deviatoric part of the fluid stress tensor is used for this calculation.

The most widely accepted empirical coefficients for the Launder and Spalding [22] k-$\epsilon$ turbulence model are given by $C_{\epsilon 1}=1.44$ and $C_{\epsilon 2}=1.92$, which are changed to $C_{\epsilon 1}=1.52$ and $C_{\epsilon 2}=1.90$ to provide agreeable results for a long-straight tube nozzle single-phase flow regime [20,23].

2.2. DEM

Lagrangian particle translations through time are governed by Newton’s Second Law.

$$m_p\frac{d{V}_p}{dt}=m_{p}g+f_{\nabla ·\tau }+f_{\nabla p}+f_l+\sum _{N_p}{f}_{p-p}+\sum _{N_p}{f}_{p-w}+f_d$$

(9)

$$I_p\frac{d{\omega }_p}{dt}=T_p$$

(10)

Velocity Verlet integration is then used on these equations to calculate the new particle positions and velocities [27,28].

The pressure force is given by

$$f_{\nabla p}=-V_s\nabla p$$

(11)

where $V_s$ is provided by the divided void fraction model to be the total volume of solid particles (both whole and partial particle volumes) inside the cell.

The viscous force models the inter-molecular forces between the particle and fluid and is given by

$$f_{\nabla ·\tau }=-\left(\nabla ·\tau \right)V_s$$

(12)

where $\tau ={\rho }_f\nu \left(\nabla u+{\left(\nabla u\right)}^T\right)$ and, therefore,

$$f_{\nabla ·\tau }=-\nabla ·\left(\nabla u+{\left(\nabla u\right)}^T\right){\rho }_f\nu =-{\rho }_f\nu \left(\Delta u\right)-{\rho }_f\nu \left(\nabla ·{\left(\nabla u\right)}^T\right)$$

(13)

as indicated in the CFDEM${}^{\circledR }$ code [10].

The lift force is calculated from Archimedes’ principle by stating that a solid body (particle) in a fluid will experience an opposing force to gravity because of the density disparity between the particle and fluid. It is given by

$$f_l={\rho }_{f}g{V}_p$$

(14)

The particle-particle interaction force, $f_{p-p}$ , has contribution of normal and tangential forces,

$$\sum _{N_p}{f}_{p-p}=\sum _{N_p}\left(F_{c,i}^n+F_{c,i}^t\right)$$

(15)

CFDEM${}^{\circledR }$ software uses LIGGGHTS${}^{\circledR }$ for particles, and the particle-particle contact model used is a non-linear dashpot type known as the Hertz-Mindlin model. Details about this contact model can be found in Di Renzo and Di Maio [29], Di Renzo and Di Maio [30]. The normal and tangential forces are contributions of spring and damping forces given, respectively, by

$$F_{c,i}^n=\left(K_n{\delta }_n-{\gamma }_n{v}_{rn}\right)$$

(16)

$$F_{c,i}^t=\left(K_t{\delta }_t-{\gamma }_t{v}_{rt}\right)$$

(17)

where $K_{n,t}$ is the stiffness coefficient, $\delta$ is the overlap distance, $v_r$ is the relative velocity between the two particles, and $\gamma$ is the damping coefficient. Further information on these coefficients is found in Appendix A.

The contact model used for particle-wall normal and tangential interactions, ${\displaystyle \sum _{N_p}}{f}_{p-w}$, is the same as the particle-particle contact model given by Equations 15-17, which are further defined in the coarse graining Section 2.4.

The force due to drag is given by

$$\begin{array}{c}{f}_d=\frac{{\beta }_p}{1-{\epsilon }_f}{V}_s{u}_r\\ u_r=u_f-u_s\end{array}$$

(18)

$$\begin{array}{c}\frac{{\beta }_p}{1-{\epsilon }_f}=\frac{0.75{\rho }_f{\epsilon }_f{C}_d\left|u_r\right|}{\left(d_g/f_{cg}\right){\epsilon }_f^{2.65}},{\epsilon }_f>0.8\end{array}$$

(19)

$$\begin{array}{c}{C}_d\left|u_r\right|=\frac{24{\nu }_f}{\left(d_g/f_{cg}\right){\epsilon }_f}\left(1+0.15{Re}_p^{0.687}\right)\end{array}$$

(20)

$$\begin{array}{c}\frac{{\beta }_p}{1-{\epsilon }_f}=\frac{150{\nu }_f{\rho }_f\left(1-{\epsilon }_f\right)}{{\epsilon }_f{\left({\varphi }_p{d}_g/f_{cg}\right)}^2}+\frac{1.75{\rho }_f\left|u_r\right|}{{\varphi }_p{d}_g/f_{cg}},{\epsilon }_f\le 0.8\end{array}$$

(21)

where $V_p$ is the volume of the particle, $u_r$ is the relative velocity, $u_f$ is either cell or interpolated fluid velocity (depending on if point forcing or volume averaging is chosen for fluid-particle coupling force), and $u_s$ is the solid velocity. It should be noted that the solid velocity, $u_s$, is the individual particle velocity when applying the force to the particles but is the average solid velocity of all particles inside a cell when applying the force to the fluid.

2.3. CFD-DEM Coupling

CFD-DEM coupling can be classified into unresolved and resolved methods. Fully resolved methods use an immersed boundary approach, where the particles are seen as a boundary condition for the CFD simulation. This method is restricted to larger particles that cover at least ten computational cells [10]. An unresolved approach uses the volume fraction of solids inside a cell and averaging methods to exchange momentum between the two phases. An unresolved method requires multiple particles to fit inside a computational cell. For our type of turbulent jet flow, a minimum cell side length of 1.86 particle diameters was found to be sufficient [16].

Coupling using an unresolved approach between the carrier (CFD) and particulate (DEM) phases can be performed in many ways. In CFDEM${}^{\circledR }$ software, the coupling can be achieved using three similar methods, with the difference being in the methodology used for calculating viscous and pressure forces by setting the $switch$ variable in Equation 2 equal to 1 or ${\epsilon }_f$ for Set I and Set II, respectively. The full description of the models is given in detail by Zhou et al. [26]. We tested Set I and Set II in these jet flow simulations and determined there were only minute differences in results. Set II goes against Newton’s third law and therefore Set I was chosen for the work. The summation of all forces for the coupling term in Set I are given by

$$\begin{array}{c}{\displaystyle \sum _i}{\overline{F}}_d=f_{\nabla ·\tau }+f_{\nabla p}+f_{drag}\end{array}$$

(22)

where the viscous and pressure terms in the momentum equation are left unchanged as indicated by the $switch=1$ variable in Equation 2.

2.4. Coarse Graining

Coarse-graining is a technique that allows a single particle (grain) to represent a cluster of multiple individual particles by scaling the size and contact parameters. Coarse-graining is very similar to the multiphase particle-in-cell (MPPIC) method. The difference is that MPPIC is more of a combination of a two-fluid model and uses indirect stress models, whereas coarse graining uses explicit DEM models for the stresses [31,32]. Even with the availability of modern computing resources, it is still costly to run simulations that include tens of millions of particles. With the coarse-graining approach, it is possible to reduce the tracking of, say, 20 million particles to tracking of $20/\phantom{20f_{cg}^3}{f}_{cg}^3$ million particles, where $f_{cg}$ is a coarse-graining factor that is typically equal to two or three. Therefore, it would be a great benefit to realize that coarse-graining is an effective method.

While researchers are developing new coarse-graining models every year, this current work uses the technique currently implemented in CFDEM${}^{\circledR }$, which is proposed by Radl et al. [33]. The basic premise of this method is to equate the density, energy densities, translational velocity, and total rotational kinetic energy of the original and coarse-grained particles and merely scale the radius of the particles. [34]. This is performed by keeping the coefficient of restitution, coefficient of friction, density, and Young’s modulus constant while scaling the radius of the particles by a coefficient $f_{cg}$. Dimensional analysis and a significant amount of reduction is then performed to determine that particle contact forces using the Hertz-Minlin non-linear contact model can be scaled directly using Equation A4 and no other modification is needed. The full derivation is shown in Appendix A.

The scaling of particle coupling forces is performed slightly differently than for the contact forces. The drag force acting on the particles was previously stated in Equations 18-21. The drag coefficient is not scaled linearly with coarse-graining as shown through the $d_g/\phantom{d_g{f}_{cg}}{f}_{cg}$ term in Equations 19 and 21. Instead, it is scaled when calculating for total drag force using the total solid volume, $V_s$, in Equation 18.

3. Simulation Setup

Seven simulations are performed for validation and are compared to the experimental work by Hardalupas et al. [1] and Lau and Nathan [7]. Table 1 summarizes particle, operating, and flow conditions for the selected validation cases.

3.1. Generalized CFD Setup

Due to the nature of CFD-DEM, a complete 3-D circular domain is needed to accurately capture the physics of the particles inside of the fluid domain. The geometries of domains considered in this work are described in Figure 1.

Discretization is performed by producing an all-hexahedral mesh using Cubit^© meshing software. A mesh sensitivity analysis was performed in our previous work, where it was concluded that an average cell size of $0.001$ should be used. This produced a minimum side length for the smallest cells of 357 $\mu m$, which is well above the requirement of the minimum cell side length of 1.86 particle diameters as noted from our previous work [16]. This mesh resolution ensured that achieve mesh-independent results was achieved while also being large enough to allow the effective coupling to the Lagrangian phase [16]. Max orthogonality of all meshes is kept below 30, as defined as the deviation angle between the line connecting two cell centers and the face normal of one [35]. Skewness is kept below $0.6$ for all meshes as defined as $\left|\overline{d}\right|/\left|\overline{p}\right|$, where $\overline{d}$ is the distance to the common face center and $\overline{p}$ is the distance between the two cell centers. Further mesh quality information can be found in Table 2.

Temporal discretization is performed using the first-order implicit Euler scheme. Spatial discretization is performed using a second-order linear upwind scheme for all velocity and turbulent fields. This native OpenFOAM${}^{\circledR }$ scheme for spatial discretization uses upwind interpolation weights depending on the gradient of the fluid velocity [36]. The linear second-order scheme is used for all interpolations between cell centered to face center values. A final solver tolerance of $1·{10}^{-8}$ is used for all scalar and vector fields.

Boundary conditions are shown in Table 3. For the velocity boundary condition, the mean velocity found in Table 1 is used to calculate a turbulence power-law profile. First, the maximum velocity is calculated by

$$\begin{array}{c}{u}_{max}=\frac{u_{mean}}{2}\left(\frac{1}{n}-1\right)\left(\frac{1}{n}+2\right)=u_{mean}\frac{\left(n+1\right)\left(2n+1\right)}{2n^2}\end{array}$$

(23)

and then the turbulence power law is used with this maximum velocity for the profile. A great description of the power law is given by Jaroslav [37]. The turbulence intensity is set to $i=0.05$, turbulence mixing length is set to $l=0.07D$, and $C_{\omega }=0.09$ for all cases.

A wall function approach is used near the wall to allow larger cell sizes and achieve effective momentum coupling between the fluid and particulate phases. It is a continuous wall function modeled from Spalding’s equation that fits a curve of the relationship between $u+$ and $y+$[23,38]. Typically, wall functions constructed using the Launder and Spalding methods require a first cell height that produces a $y+$ in the logarithmic layer $\left(30\le y+\le 300\right)$ but OpenFOAM${}^{\circledR }$ uses adaptive wall functions proposed by Kalitzin et al. [38], which obtains agreeable solutions from $11\le y+\le 111$. A $y+$ between 24 and 76 for all simulations is obtained as indicated in Table 2. Therefore, we are confident that the wall boundary layer is being effectively modeled, as was also validated in our previous publications [16,23].

3.2. Generalized DEM Setup

Typically, glass beads have a very high Young’s Modulus of $6.89·{10}^{10}$ Pascals (Pa). To fully resolve the particle collisions with this Young’s Modulus would require an extremely small time step, which is impractical and could be over ${10}^3$ times more computationally expensive depending on the leeway given for the critical time step criterion. Chen et al. [39] tested a wide range of Young’s Moduli, reducing them to $0.00001E$, and concluded that although Young’s Modulus has significant effects on single particle collisions, the bulk behavior of particles in the system (a rotating drum) was conserved. Lommen et al. [40] performed a study testing the accuracy of reducing Young’s Modulus with the intent of speeding up DEM simulations. It was determined that, at a minimum, if Young’s Modulus is kept greater than or equal to ${10}^7$ Pa, there is no change in the simulation results. There are very small relative velocities with the simulations in this current work and, therefore, it is possible to reduce this Young’s Modulus even further. A comparison of simulations with a Young’s Modulus of $5·{10}^6$ and ${10}^7$ is first performed to determine if there is any effect of Young’s Modulus on simulation results. An average relative error of 0.06% for particle velocities and 1.49% for particle fluxes was realized, with the majority of an error occurring at the outer regions of the jet where particle flux approaches zero. Therefore, it was decided that a Young’s Modulus of $5·{10}^6$ would suffice to reduce the computational cost.

The coefficient values in the present study (except Young’s Modulus) are taken directly from Lorenz et al. [41] for glass beads in which experiments were performed colliding two nearly spherical glass particles. Poisson’s ratios for glass beads vary depending on the paper, but a value of $0.3$ is common throughout [42,43]. Tang et al. [43] performed a wonderful study on testing for these parameters between a glass bead and glass plate and found a Poisson’s ratio, $\nu$, equal to $0.245$, a coefficient of restitution, ${\epsilon }_n$, equal to $0.926$, a coefficient of friction, $c_f$, equal to $0.18$, and a coefficient of rolling friction, $c_{rf}$, equal to $0.01$. Whether this transforms over to particle-particle contact is uncertain. Another study found a ${\epsilon }_n$ equal to about $0.8$ for varying sizes of glass beads in the range from 2 to 4 mm with the value decreasing with a reduction in glass beads size. These coefficients are not highly critical because the volume fraction of particles is sufficiently low, where the drag force will dominate the flow. This is another example of a need for research that has detailed and defined inputs for the material properties of the particles.

The experimental data at the nozzle exit is used to obtain particle insertion velocities and fluctuations. DEM particles are inserted on a surface and the insertion velocities and fluctuations do not vary along that surface. Therefore, the mean of the experimental velocities in the radial direction of the nozzle is used with the inputs shown in Table 4.

3.3. Hardalupas Et al. DEM Setup

The solids inserted into the system are "mono-sized" 40 and 80 µm spherical glass beads with the velocities given in Table 4. It should be noted that it is stated to be "mono-sized," but it is extremely difficult to obtain true mono-sized particles because of the expense and quantity needed for these types of experiments. Therefore, a size range is used, with the mean being the mono-sized value and a range small enough to be considered negligible to overall results. The 40 µm particles of a nominal range of 37-44 µm and the 80 µm have a nominal range of 60-95 µm. The information on the particle size distribution can be found in Hardalupas et al. [1].

3.4. Lau and Nathan DEM Setup

Polymer particles of density $1200kg/\phantom{kg{m}^3}{m}^3$ are used for this study with the inlet velocities given in Table 7 . These particles have an extremely low nominal range as compared to the particle size distribution (PSD) of Hardalupas et al. [1][7]. Other than the PSD of the particles, the paper does not go into details on the specific spheres used, but in a previous paper by Lau and Nathan, with what appears as the same setup, they state that they are Microbeads Spheromers [44]. The Microbead Spheromers are made with polymethylmethacrylate (acrylic). Polymethylmethacrylate has a density of 1.18 and Flexural Modulus $2.9·{10}^9$ Pascals. Similar to the Hardalupas et al. [1] test cases, we reduce the Young’s Modulus to $5·{10}^6$ for help with the computational expense. Poisson’s ratio of acrylic is about $0.37$, and this is what is used in the DEM simulations. The coefficient of restitution is given by ${\epsilon }_n=0.934$ and coefficient of friction by $c_f=0.096$, which are taken directly from Lorenz et al. [41] in which experiments are performed on the collision of two $4mm$ acrylic spheres. It is assumed that the coefficient of rolling friction is extremely low for these acrylic spheres, and therefore Coefficient of rolling friction, $c_{rf}$, is set to the low value of $0.01$

3.5. Coupling Setup

Coupling is performed implicitly between the two phases. The soft sphere non-linear spring-dashpot contact model of Hertz-Mindlin is used. To fully model and ensure particles are observing one another, a DEM time-step lower than 20% Rayleigh and Hertzian critical timesteps is used [28]. This produces a very small timestep; therefore, a coupling interval to the CFD side is set to 10 to facilitate less computational expense while keeping the fluid Courant number less than $0.5$. DEM time-steps and specifics for percentages of Rayleigh and Hertzian time-step criterion are shown in Table 4. The Gidaspow drag model is used, which combines the drag models of Ergun [45] and Wen and Yu [46]. The full model description can be found in Zhu et al. [47]. Pressure, viscous, and Archimedes lift forces are used and calculated from Equations 11, 13, and 14, respectively. Second-order linear interpolation is used to obtain face-centered values from cell centers for velocity, void fraction, and pressure in all calculations of forces on the particles.

4. Results and Discussion

Mesh sensitivity analysis was first performed to determine if the numerical solution is independent of domain discretization (mesh resolution). Three systematically refined meshes of $2.88$, $3.80$, and $6.45$ million elements were used with the same operating parameters. Richardson extrapolation was going to be used, but it was found that the solution differed so minutely that Richardson extrapolation would hold no quantitative value for the solution of mesh convergence. Instead, it was deemed that the medium mesh size of $3.80$ million elements would be used. Simulations were performed in our previous work that demonstrates a particle diameter to minimum side cell size ratio of up to $1.86$ particles diameters can be used while still conserving the accurate momentum exchange between the two phases.

4.1. Particle Full Development

A fully developed flow is critical in these simulations. If we did not have a fully developed flow before the outlet of a simulation, then particle flux and velocities would entirely depend on both the insertion location and initial velocities. It is difficult to conclude that the experimental work performed by Hardalupas et al. [1] is a fully developed flow. The pipe used in the study is $93D$ in length with a Reynolds number of 13000. For single-phase flows, this Reynolds number is well into the turbulent zone; therefore, no laminar boundary layer forms and there will be a turbulent velocity profile. To further demonstrate this, an idealized entry length is calculated until fully developed flow by using an assumed velocity profile at the inlet and the classical solution developed by Latzko [48] given as

$$\begin{array}{c}x/\phantom{xD}D=0.636{Re}^{0.25}\end{array}$$

(24)

where $Re=\left(dG\right)/\mu$ for the circular tube. Calculating this, an entry length of $x/d\approx 6.8$ is obtained for an idealized single-phase flow. The pipe used for Hardalupas et al. [1] is well above this entry length, but this does not indicate for certain that the particles are fully developed because of particle slippage during particle insertion and mixing. Indeed, this is briefly discussed in the paper when the results for the highest particle size ($200\mu m$) showed signs of a lack of fully developed flow. In the work performed by Mena and Curtis [49], it is postulated and concluded that 51 pipe diameters is sufficient for fully developed flow and they measure this by the use of pressure transducers placed periodically along the pipe. In the study, they use water as the carrier fluid; therefore, viscous effects are significantly higher than air and particle slip velocity would be less of a concern.

To the current authors’ knowledge, Lau and Nathan [44] and their future work of Lau and Nathan [7] are the only studies to test for fully developed flow. They perform particle experiments with $L/d\approx 163.8$ and $L/d\approx 301.3$ with identical Stokes numbers to test for particle slippage. They realized that the particle concentration and velocities converged to similar values and therefore concluded that $L/d\approx 163.8$ was a sufficiently long enough nozzle for full flow development at Stokes numbers of $0.3-11.2$ and Reynolds numbers of 10000 and 20000. In the study performed by Hardalupas et al. [1], a Stokes number of $St=50$ and $St=10.27$ with a Reynolds number of 13000 is used. There is a possibility that with the larger Stokes number of 50 and the fact that Lau and Nathan [44] has a pipe length almost $1.75$ times that of Hardalupas et al. [1] that we may not have fully developed flow. Therefore, no conclusion can be realized.

To test for fully developed flow for the numerical Hardalupas et al. [1] test case, simulations were performed using merely the nozzle or pipe flow with the same operating parameters as previously outlined. Still, instead of particles inserted at $45D$, they are inserted at $60D$ with sampling at locations at $0D$ and negative $5D$, $10D$, $15D$, $20D$, $25D$, $30D$, $35D$, $40D$, $45D$, $50D$, and $55D$, as referenced from the nozzle exit with the positive in the flow direction. In addition, two simulations were performed to test for the location of full development and determine whether inlet velocity makes a significant difference in flow development.

It was realized that most of the flow development occurs in the first 40 nozzle diameters, where the last $20D$ experiences very minute changes in all parameters (Figure 2). The inside of the nozzle for the Lau and Nathan (Lau3, Table 1 and Table 4) test case are also examined. It was realized that full development occurs in about the first 30 nozzle diameters with only minute differences observed in the last 10 nozzle diameters Figure 3.

4.2. Hardalupas Et al. Results

It is realized that the results at the nozzle exit for all the Hardalupas et al. [1] test cases provide a similar trend for the mean and root mean square velocity near the central regions of the jet stream. Still, towards the outer regions near the nozzle wall, it is observed that the particles move more towards a laminar shape profile rather than a more turbulent shape as demonstrated by the experiments (Figure 4a and Figure 5a). This is believed to be caused by the interpolation method used for velocity in particle-fluid interactions.

To achieve effective momentum coupling, we require a cell size with a minimum length side of at least $1.86$ times the particle diameter for these types of flow problems [16]. To model the boundary layer at the wall of a CFD simulation, we either need to fully integrate into the wall or use an algebraic function. To fully integrate to the wall, we require an extremely fine mesh with, at the very least, 8-10 cells located below the $y+=11.5$ [23]. Achieving this, while also achieving effective momentum coupling, would require extremely small particles, far smaller than those used in many cases including the experimental studies considered in this current work. Therefore, a wall approach near the wall is used in this current work to accurately model the strong viscous effects on fluid velocity while achieving effective momentum coupling. That said, there isn’t a perfect solution here because there is also concern with fluid velocity and void fraction interpolation when calculating the drag particle-fluid coupling force. A wall function is used in this current work to model the fluid boundary layer in the larger cell. Still, for interpolation, that same wall function is not used for particle coupling parameter interpolation, but instead, a linear interpolation method.

In CFD-DEM, the interpolation method used is a linear cell point method. This method breaks each face into triangles to define tetrahedra and the cell center point; then, it cycles through the tetrahedra to determine where the point or particle cell point center is located. An inverse distance linear interpolation is then used from known velocity points of the cell centers. This produces a different distribution of coupling velocities than is modeled by our wall function. Optimally, we would want to use an interpolation method that follows the wall function but is out of this work’s scope.

Although there are discrepancies at the nozzle exit it is observed that further in the flowfield there are very similar trends between the numerical and experimental results for the lower (0.23) and higher (0.86) mass fractions, respectively.

In general, for the $80\mu m$ particle simulations (Hardalupas1 and Hardalupas2), discrepancies are observed at the nozzle exit associated with the interpolation method for both mean (rms) velocities and particle number fluxes (Figure 4a). As previously outlined, this can be attributed to the issues with the interpolation method near the wall. That being said, relative accuracy is achieved for the velocities where the general change in results between the two mass fractions is similar to that of the experiments. At location $10D$, similar slopes are achieved between all numerical and experimental results. The numerical results show a slight decrease in slope relative to the experimental results near the outer regions of the jet. At $20D$, higher mean velocities are observed for the numerical results with root mean square velocities closely following and smaller spread with the particle number fluxes. It is concluded that with $80\mu m$ particles, up to $10D$ should see viable results for the case of optimization of an industrial unit, but anything further downstream would need more analysis. Work needs to be performed for pipe flows with wall functions in CFD-DEM to solve the interpolation issues, but it is observed that the particles “correct” themselves after they exit the nozzle. This could be attributed to having relatively low Stokes numbers, and therefore the particles closely follow the highly accurate and validated fluid flow where wall effects are not present.

For the $40\mu m$ particles, the same errors at the nozzle exit are observed because of interpolation for all fields but with the same relative accuracy or change in trends between the two mass fluxes (Figure 5ab). An interesting note for the nozzle exit is that with smaller mass loading (Hardalupas3), a particle number flux with a more flat profile than that of the larger mass loading (Hardalupas4) is observed. It is postulated that if mass loading continued to increase, a more triangular profile for particle number flux would be achieved. Further downstream, similar slopes in the results for mean velocity and particle number flux are observed but a significantly lower root mean square velocity for the numerical results. It is curious to note that these smaller particle trends line up quite well up to $20D$, but with the larger particles results can only be trusted up to $10D$. It can be concluded that as Stokes number increases, where the particles tend to "go their own path," the numerical results start deviating from the experimental results.

4.3. Lau and Nathan Single-Phase Results

The single-phase results for the Hardalupas et al. [1] test cases were validated in our previous work of Weaver and Mišković [23], but Lau and Nathan [7] results were not. Therefore, three different simulation setups were used to determine that single-phase results were accurate before adding the particulate phase. If this is not done, then there is a risk of obtaining a model that may have the perfect combination of parameters for the two-phase flow but not necessarily because each of the phases separately is accurate. This would make it difficult to conclude with any confidence that this model can be generally applied to other similar flows and accurately represents the physics.

All cases start with the Lau3 test case ($48m/\phantom{ms}s$) but without the inclusion of particles. The default Launder and Spalding [22] k-$\epsilon$ is used for one case. For the second case, a change of the empirical coefficients is made to the turbulence model by $C_{\epsilon 1}=1.52$ and $C_{\epsilon 2}=1.90$. This modification is researched and validated extensively in our previous work of Weaver and Mišković [23]. The third simulation uses the modified k-$\epsilon$ turbulence model but without the use of co-flow. This is to test the effect that co-flow has on these types of numerical simulations.

There is very little difference between all simulations at the nozzle exit for both mean and root mean square velocities as shown in Figure 6ad. At location $10D$, the unmodified k-$\epsilon$ turbulence model produces spreading that more closely follows the experiment for mean velocity, but the modified k-$\epsilon$ produces a drop in axial velocity that more closely matches the experiment. Therefore, it is up to the researcher to weigh the pros and cons of using the modified k-$\epsilon$; whether it is important for accurate spreading or accurate axial drop in mean velocities. That being said, in the current authors’ previous work, it was found that both spreading rates and axial drop in mean velocity were both improved with the slight modification of the empirical coefficients [23] for two different sets of data provided by Bogusławski and Popiel [24] and Hardalupas et al. [1]. For root mean square velocities, the modified k-$\epsilon$ produces spreading rates that more closely match experiments at all sampling locations. There is no appreciable difference in all sampling locations in adding a co-flow to simulations. Considering all this, it was decided to use the modified k-$\epsilon$ with a co-flow for all subsequent simulations for the Lau and Nathan [7] test cases.

4.4. Lau and Nathan Results

As previously noted, there will be errors at the nozzle exit associated with the interpolation method and this also extends to the Lau and Nathan [7] results as indicated in Figure 7ab. A more laminar flow profile is realized for the mean particle velocity with a higher mean velocity in the central region while decreasing rapidly towards the wall relative to the experimental results. That being said, it is observed that the overall trend when changing the Stokes number agrees with the experimental results. For example, the lowest Stokes number has the highest mean particle velocity in the central region. Conversely, the highest Stokes number has the lowest mean particle velocity for both the experiments and numerical simulations (Figure 7a). At the location $10D$ we observe a slightly steeper slope for the mean particle velocities and a slightly higher root mean square velocity with the same slope. It is interesting to note that mean particle velocity changes at $10D$ depending on the Stokes number. In contrast, in the experiments, the normalized mean particle velocity minimally changes between Stokes numbers. This could be attributed to the errors at the nozzle exit, but this is difficult to postulate with confidence because we are not experiencing this behavior with the Hardalupas et al. [1] test cases.

4.5. Coarse Graining Results

Some DEM simulations use a large number of particles, upwards of 6 million. This would require a significant amount of compute cores to track and calculate all collisions between particles. Also, developing codes to have the ability to scale to a very large amount of compute cores can be difficult because of the issues associated with the slowed communication between compute nodes. To circumvent this, many researchers have turned to coarse-graining methods. The essence of coarse-graining is that multiple particles are represented as a single grain and then the grain is tracked through time. This significantly decreases the number of trajectories tracked throughout the system and, depending on the type of coarse-graining model, can reduce the number of tracking points by $n_g=n_p/\phantom{n_p{f}_{cg}^3}{f}_{cg}^3$. We can see that with a coarse-graining factor of 2, $f_{cg}=2$, the example of 6 million particles reduces down to 750 thousand particles. This reduces overall computational cost significantly. Also, additional computational savings can be realized by increasing the integration time-step and the duration of the soft sphere collisions will generally be smaller by orders of magnitude [32].

As previously mentioned in Section 2.4, the type of coarse-graining used for this analysis was first proposed by Radl et al. [33] and further expanded to include the Hertz nonlinear contact model by Nasato et al. [34]. Radl et al. [33] conclude that the major hurdle in using a coarse-graining approach is to correctly compute the collision rate and inter-parcel or inter-grain stress, with there being no easy way to scale the interaction parameters in an inertial flow regime [33]. This conclusion is significant because it will be difficult to relate this current work collision frequencies, forces, and stresses from an unscaled system (no coarse-graining) to a scaled (coarse-grained) one. That being said, if overall trends are followed, we may find the usefulness of coarse-graining in an industrial application.

To optimize an industrial scale unit, a statistical analysis needs to be performed on parameters, including collision frequencies and forces acting on particles. Therefore, at first glance, it would seem that coarse-graining shouldn’t be used in this type of setting because of the disparity between inter-grain stresses from the un-scaled to scaled system. That being said, the goal of numerical simulations, and any experimental study for that matter, is to find the optimum setup through an analysis of changes in geometric and operating parameters. Since this is the ultimate practical goal, and it may not be necessary to obtain the exact answers of the physical space, we can, for example, compare two different geometries of an unscaled system to a scaled system using those same geometries.

To investigate the accuracy of coarse-graining, two sets of simulations with the same operating parameters as the Hardalupas1 case are performed. The first set of simulations is unscaled, while the second set is scaled. In each set, a plate is inserted at two different axial locations: $5D$ and $15D$ relative to the nozzle exit, resulting in four simulations in total. The collision frequencies and force statistics applied to particles when they hit the plate for both the $5D$ and $15D$ cases are then outputed. The ratio of results obtained from the $5D$ and $15D$ simulations is calculated and compared the ratios between the scaled and unscaled systems to determine the accuracy of coarse-graining. We then analyze the collision frequencies and mean/variances of the normal/tangential forces for all collisions (Figure 8a), and also analyzed the results grouped by location with intervals of $0.5D$ on the plate. These results and then used to determine whether the coarse-graining method accurately predicted the results that follow the simulation trend of an unscaled system.

The ratios we are comparing for collision frequencies are given by

$$\begin{array}{c}rati{o}_{cf,cg0}=\frac{C{F}_{cg0,5D}}{C{F}_{cg0,15D}},rati{o}_{cf,cg2}=\frac{C{F}_{cg2,5D}}{C{F}_{cg2,15D}}\end{array}$$

(25)

where $C{F}_{cg0,5D}$ and $C{F}_{cg0,15D}$ are the collision frequencies for the unscaled DEM method for simulations with 5D and 15D plate, respectively. $C{F}_{cg2,5D}$ and $C{F}_{cg2,15D}$ are the collision frequency for the scaled system with the coarse-graining factor of 2 for simulations with 5D and 15D plate, respectively. Similarly, for the force statistics we have

$$\begin{array}{c}rati{o}_{F_{n,CG0}}=\frac{F_{n,cg0,5D}}{F_{n,cg0,15D}},rati{o}_{F_{n,CG0}}=\frac{F_{n,cg2,5D}}{F_{n,cg2,15D}}\end{array}$$

(26)

$$\begin{array}{c}rati{o}_{F_{t,CG0}}=\frac{F_{t,cg0,5D}}{F_{t,cg0,15D}},rati{o}_{F_{t,CG0}}=\frac{F_{t,cg2,5D}}{F_{t,cg2,15D}}\end{array}$$

(27)

The collision frequencies and force statistics are then output for all collisions (Figure 8a) and concentric regions (Figure 8b), and calculated ratios and associated errors as shown in Table 5.

It should be noted, for completeness, that the percent error is calculated from the traditional formula given by

$$\begin{array}{c}\%Error=\frac{100·\left|rati{o}_{cg0}-rati{o}_{cg2}\right|}{rati{o}_{cg0}}\end{array}$$

(28)

When analyzing all the collisions acting on the plate, it is realized that there is an error of less than $4\%$ for all parameters as indicated by Table 5. These results demonstrate accuracy in comparing ratios of scaled to un-scaled systems when analyzing all collisions on the plate.

To further analyze the trends of particle collisions with the plate, the collisions are independently considered within specific concentric regions on the plate as shown in Figure 8b and perform the same analysis of collision frequency and force statistics.

In Figure 9 the relative errors are plotted between the scaled and unscaled system ratios for each region. A slightly different conclusion than when analyzing all collisions can be realized. The majority of mean values have an error less than $10\%$ with the only region above being at $r/\phantom{rd}d=3.25$ with an error of just above $20\%$. The highest errors in the force statistics are the variance errors which are in line with the collision frequency errors; most errors are observed in regions where the collision frequency is low. To investigate this further, (Figure 9b) shows individual collision frequencies for both unscaled and scaled systems at 5D and 15D plate positions with corresponding collision frequency ratio relative errors. It is realized that collision frequency is high near the center region with a very lower percentage error. Conversely, collision frequency is lowest in the outer regions with a relatively high percentage error. Considering all this, along with the very low error of less than $4\%$ for all collisions, it can be concluded that we have achieved relative accuracy that provides significant results for the practical application of an industrial flow problem.

A kernel density estimate plot (KDE) is then used to visualize the distribution of collisions on the plate (Figure 10). This plot uses a Gaussian smoothing algorithm to estimate the probability density function. This allows visualization of the distribution of particle impacts across the plate. It is realized that the probability of particle collisions between $CG=0$ and $CG=2$ are very comparable with the largest difference towards the outer regions as demonstrated when comparing at location $15D$ (Figure 10b and Figure 10d). This reassuringly confirms the use of coarse-graining for an industrial flow optimization problem.

Of special interest is comparing the flow profiles both inside of the nozzle and outside of the jet. This could be beneficial in future simulations because if flow development occurs sooner for coarse-graining, that further reduces the computational cost by allowing a shorter nozzle.

Very little difference in the length of full development for the particles (Figure 11 in all parameters is observed. Furthermore, it is realized that there is also very little difference between the scaled and unscaled systems in the jetstream for all fields, with the highest error being observed for root mean square velocity (Figure 12).

If the desire is to optimize an industrial unit, then the goal is to pinpoint a design that performs better than another. Therefore, from a practical standpoint, if relatively accurate answers are achieved that follow the same trends as changes in the design of a physical system, then the goal can be reached from purely numerical simulations alone. With that in mind, coarse-graining is deemed a very useful tool in the optimization of an industrial unit for high-speed jet flows because of the significant reduction in simulation cost. That being said, these results and conclusions shouldn’t be realized for other types of systems, in particular, systems with a large amount of particle-particle collisions in a semi-quasi steady state inertial system that significantly dictates the bulk flow behavior.

5. Conclusion

This paper has performed fully coupled CFD-DEM simulations and validated them against two different experimental studies, namely Hardalupas et al. [1] and Lau and Nathan [7]. Coupling parameters, models, and operating conditions were all considered. However, issues were encountered with the interpolation method used in fully coupled simulations, which produced a more laminar profile shape for the particles at the nozzle exit. Future work can investigate this issue further.

For the Hardalupas et al. [1] test cases, a Stokes number of 50 can be trusted up to 10D from the nozzle exit, and a Stokes number of 10.27 can be trusted up to 20D if only mean and particle number fluxes are desired. The modified k-epsilon models are consistent with the Lau and Nathan [7] single-phase results, providing further validation of previous work [23]. For the Lau and Nathan [7] test cases at 10D, there was a slightly steeper slope for mean particle velocity with a higher root mean square velocity. Overall, it is concluded that CFD-DEM accuracy generally increases with increasing Stokes number up to 50, based on both studies.

Finally, it is noted that coarse-graining is a useful method for practical applications such as optimization of industrial-scale units, where relative results between designs are desired. However, if the goal is to measure the total energy dissipation of a single experiment, coarse-graining falls short.