A Study of RANS Turbulence Models in Fully Turbulent Jets: A Perspective for CFD-DEM Simulations

: This paper presents an analysis of linear viscous stress Favre averaged turbulence models for computational ﬂuid dynamics (CFD) of fully turbulent round jets with a long straight tube geometry in the near ﬁeld. Although similar work has been performed in the past with very relevant solutions, considerations were not given for the issues and limitations involved with coupling between an Eulerian and Lagrangian phase, such as in fully two-way coupled CFD-DEM. These issues include limitations on solution domain, mesh cell size, wall modelling, and momentum coupling between the two phases in relation to the particles size. Therefore, within these considerations, solutions are provided to the Navier–Stokes equations with various turbulence models using a three-dimensional wedge long straight tube geometry for fully developed turbulence ﬂow. Simulations are performed with a Reynolds number of 13,000 and 51,000 using two different tube diameters. It is found that a modiﬁed k - ε turbulence model achieved the most agreeable results for both the velocity and turbulent ﬂow ﬁelds between these two ﬂow regimes, while a modiﬁed k - ω SST/BSL also provided suitable results. numerical on turbulent jets. Studies include fundamental work and modiﬁcations that were performed to help mitigate the problems associated with spreading rates.


Introduction
Round turbulent jets have a wide range of applications that include combustion, jet cutting, and surface coating and they have been studied extensively, numerically and using experiments. The turbulent round jets can be categorized into three main groups: (i) particle-laden jets, (ii) impinging jets, and (iii) free round jets, each of which has distinct geometries. A long straight tube nozzle is commonly found in particle-laden jets applications with the particle mixture injected into the inlet. The choice of the particle-laden flow geometry is driven by the need to achieve fully developed flow before the nozzle exit to help quantify spreading rates with variations of particle diameters and Stokes number with a higher accuracy. This is demonstrated in the papers given by [1][2][3][4], in which 93, 90, and 163 nozzle diameters are used for the lengths of the nozzles, and the flow fully develops by the nozzle exit. Similar long nozzles are also shown in many other studies, including [5][6][7][8][9]. Using test cases that have fully developed flow provide consistency between the boundary conditions of the numerical simulations and input parameters of the experimental data.
The geometries for the application of impinging jets range from a smooth contraction to straight nozzles, and they may or may not have fully developed flow [10]. The free jets typically have contraction or sharp edged slots without fully developed flow [11].
There have been a plethora of numerical studies over the past years on jet turbulence modelling as shown in Table 1, but the geometries and operating parameters are not common across these studies. This difference is because the focus is on pure, single-phase CFD simulations without the considerations and issues indicative of CFD-DEM coupling.
The most common way to model a jet in numerical simulations is using a fixed boundary condition taken from DNS or experimental data for pipe flow. This circumvents Table 1. Relevant RANS numerical studies on turbulent jets. Studies include fundamental work and modifications that were performed to help mitigate the problems associated with spreading rates.

Reynolds Number Geometry Nozzle Type Turbulence Models Key Features
Pope [13] NA Axi-Symmetric Round jet Modified k-epsilon Suggested changing C ε1 = 1.6 and proposed vortex stretching term Givi and Ramos [15] 373,000 Axi-Symmetric Round jet Modified standard k-ε and k-ω turbulence models Tested k-epsilon constants and C ε1 = 1.52 with C ε2 = 1.90 provided most agreeable spreading rates Kim and Chung [16] 11,600 Axi-Symmetric  Comparable results between k-ω Wilcox and modified k-epsilon Quinn and Militzer [11] 208,000 Axi-Symmetric

Sharp Edged and Smoothly Contoured
Modified k-epsilon (C ε1 = 1. 48) Comparable results between experiments for velocity but not turbulence quantities. Researchers have implemented and tested many methods to circumvent this issue, such as adding vortex stretching, cross diffusion terms, and, also, modification of the empirical coefficients, and this has provided accurate results with a constant edged boundary condition. However, very few studies have been performed for long straight tube geometries that arise in fully coupled CFD-DEM simulations and optimization problems. Therefore, this present paper's contribution is modification and then testing of RANS turbulence models for the application of long straight tube jet flow that are required for CFD-DEM simulations to ensure that the results are agreeable to the experimental data.

Background
In experiments, anisotropic turbulence behavior is observed in high-speed jet flow due to the strong boundary layer separation at the exit of the nozzle, which makes this type of flow difficult to model. This anisotropic turbulence has also been demonstrated in numerical simulations when using LES Bonelli et al. [22]. RANS turbulence models studied in this paper use the very common linear viscous Boussinesq approximation that assumes that the momentum transfer caused by turbulence can be modelled by an eddy viscosity. In using this model, we assume that the eddy viscosity is independent of direction and, hence, the rate of turbulence generation is uniform in all directions. This is not to say there is no anisotropic turbulence when using these types of models, because the mean strain rate tensor of the velocity field in the Navier-Stokes equations is not necessarily isotropic. Suffice to say, there is an expected error in using RANS approaches because of anisotropy, but if we can capture an average quantification, then we can achieve results that are sufficient for understanding of the flow behavior while achieving adequate single-phase results before the addition of particles in unresolved fully coupled CFD-DEM simulations.
Turbulent round jets have been difficult to model numerically since the beginning of RANS turbulence modelling by producing spreading rates that are not agreeable. This is shown in many papers, including Pope [13], Quinn and Militzer [11], Launder and Spalding [23], Morgans et al. [17], Givi and Ramos [15], and also stated in David C. Wilcox's famous book Turbulence modelling for CFD [14]. RANS turbulence models have empirical coefficients that have been adjusted for generalized flow conditions. In an ideal world, these coefficients would be set from principle values, such as in kinetic theory of gases. Instead the developed RANS turbulence models are derived from dimensional analysis and some intuition on theory and, therefore, the coefficients are merely empirical providing agreeable results for a wide range of generalized flow conditions. Free shear flows, such as jet flows, are one of these generalized flow conditions.
In 1978, Pope suggested that changing the constant, C ε1 = 1.6, in the k-ε model will solve the issue with round/plane jets, but states that this goes against any kind of generality [13,17]. That being said, researchers have adopted the method of changing C ε1 and it is perfectly valid because, as stated before, these coefficients are merely empirical relations to fit to a generalized flow regime [14,15,17,20]. Instead, Pope proposed a modification to the k-ε model by adding a vortex stretching variable to the destruction of turbulence term. The ε equation in the standard k-ε model is modified by where, χ p is the non-dimension measure of vortex stretching given by and Ω ij and S ij are the mean-rotation and mean-strain rate tensors defined as In the OpenFOAM tensor library, Equations (3) and (4) can be obtained by decomposing to the skew and symmetric parts of the gradient of velocity, U, respectively. That is: Ω ij = skew(∇U) and S ij = symm(∇U). The theory is that the primary mechanism of turbulent energy dissipation is through stretching of the vortex rings in three-dimensions, but not in two-dimensions. Rubel [24] tested this model and determined that while this correction helps rectify the problems with round jet modelling, plane jet modelling was still experiencing problems. Because, exclusively, three-dimensional wedge simulations are considered in this work, the issues with plane jet modelling are ignored.
Contrasting from the k-ε model, the Wilcox k-ω turbulence model produces spreading rates for the round nozzle that are comparable, but larger, than predicted spreading rates for the plane jets [25]. To solve this issue for the plane jets while conserving the correct predicted values for round jets, Wilcox realized that, when the constant β is reduced to 0.06 from 0.0708, the predicted spreading rates for both the round and plane jets were close to the experimental values. Rather then changing the value of the β constant, Wilcox proposed a more "general" solution that reformulates the k-ω model by where β 0 = 0.0708, Although this may solve the issue with the round/plane jet anomaly, the k-ω turbulence model is sensitive to free-stream boundary conditions; this could present an issue in determining a solution that is invariant with different initial free-stream conditions [25,26] and, to that affect, this turbulence model is not considered in this present work.
In this present work, we compare the velocity and kinetic energy profiles for various wall function RANS approaches to turbulence modelling. Two distinct test cases are considered with a Reynolds number of Re dh = 13,000 for the test case given by Hardalupas et al. [1] and Re dh = 51,000 for the test case given by Bogusławski and Popiel [27].

Numerical Model
The fluid flow in the present research is described by the Navier-Stokes equations given by where, u is the fluid velocity, p is the pressure, ρ 0 is the density of the fluid, ν m is the viscosity of the fluid, and x is a unit in space. It should be noted that the assumption of an incompressible flow is made in this work since the considered velocities are far from the sonic velocity and, therefore, the following term becomes but is still included in the formulation for solution boundedness. Turbulent flows can be calculated for each point in the flow field by solving the Navier-Stokes equations (7) and (8), directly; this approach is called direct numerical simulations (DNS). For this to be accurate, all the spatial and temporal turbulent scales need to be resolved, requiring an extremely fine mesh with a number of mesh points satisfying N 3 ≥ Re 9/4 , which is extremely computationally expensive. Less costly alternatives are to solve for the mean pressure and velocity fields using a Reynolds averaged approach (RANS), or to solve for the large scales directly while modelling the subgrid scale turbulence by the use of a large eddy simulation approach (LES).
Assuming that the velocity and pressure can be described by the mean and fluctuation components, u =ū + u and p =p + p , and with some lengthy reductions the Navier-Stokes equations become the RANS equations given by where, u i u j RANS Reynolds stress term. By examining the difference between the Navier-Stokes (Equations (7) and (8)) and the RANS it can be realized that the difference between them is a Reynolds stress term, and the use of mean pressure and mean velocities instead of instantaneous. This is an important realization during code development because we can merely specify an additional term when RANS or LES filtering is desired for turbulence modelling. It should be noted that many CFD codes use this method, including OpenFOAM. The Reynolds stress term, u i u j , now needs to be modelled to close the RANS equations. We can calculate the Reynolds Stress term directly with the addition of six transport equations, but this approach is expensive. There are single transport equation models, such as the mixing-length and Sparlat-Allmaras, but their accuracy is lacking. The most common method, with a good balance in required computational resources and accuracy, is the use of an eddy viscosity turbulence model with two transport equations, such as the k-ε, k-ω, and their variants. The Boussinesq eddy viscosity assumption is now used to model the Reynolds stress tensor for both subgrid LES and RANS turbulence models. The Boussinesq eddy viscosity assumption for the Reynolds Stress with incompressibility is given by where the first term on the right hand side is the deviatoric part, the second is the isotropic part, and k is specific kinetic energy. We then plug Equation (11) into Equation (10) to obtain the Navier-Stokes momentum equation in the OpenFOAM code.
where ν t is the turbulent viscosity andp is a modified mean pressure that has absorbed the isotropic part of the Reynolds stress term.
It is now prudent to discuss the Navier-Stokes in relation to fully coupled unresolved CFD-DEM simulations to gain a better understanding on why a wall function approach with a coarser mesh is required. Coupling between the carrier phase and solid phase involves the use of volume fraction of solids with the Navier-Stokes equations given by where α l is the liquid phase volume fraction, u l is the liquid fluid velocity, u s is the mean solid particle velocity, τ l is the liquid phase stress tensor, K sl is the implicit momentum source term that is a volume average of all interacting forces acting on the solid phase due to the motion of the fluid in each cell, and f is any other explicit forces acting on the fluid. As can be realized by the Equations (14) and (15), the momentum of the fluid is highly dependent on the volume fraction of solids, α l , and hence, the calculation of this term is critical to obtain accurate physics. The most accurate solution to this would be to calculate the volume of the particle ratios directly in the carrier phase cells, but this is impractical because it would be extremely inefficient. Instead, a volume fraction model is used that allows efficient computation while providing reasonable results for accurate flow physics.
To ensure that this model is correct, the particle size to cell size ratio must have minimal particle overlap [28,29]. It is, therefore, unfeasible to fully resolve the boundary layer and have an overall fine mesh in these types of simulations unless the desired particle size is extremely small The most common two-equation eddy viscosity turbulence models can be categorized into k-ε and k-ω turbulence models with variants. Here, an equation for turbulent kinetic energy, k , and for a turbulence length scale, typically using ε or ω, are used. The k-ε model is the more popular of the two, with the generalized formulation given by where eddy viscosity is given by ν t = C µ f µ k 2 ε and the Boussinesq assumption is used to obtain and k is the turbulent kinetic energy, ε is the turbulent dissipation rate, u is fluid velocity, f is a damping function, and the turbulence length scale is given by l = C µ k 3/2 ε. Table 2 shows the values or formulations of the generalized variables for the standard Launder and Spalding [23], variants, and modified k-ε turbulence models.
For the Chien model we make an estimation for the non-dimensional wall distance, y + . Wallin and Johansson [34] demonstrated that the non-dimensional wall distance as a function of kinetic energy is approximately equal to the non-dimensional wall distance as a function of fluid velocity for y + ≤ 100; that is, y * ≈ y + for y + ≤ 100. With this in mind, Wallin and Johansson [34] suggested using an alternative scaling of viscosity given by It should be noted that f µ ≈ 1 at y + = 300. For the RNG turbulence model η = Sk ε, For completeness, we will state the previously mentioned Pope correction to the standard Launder and Spalding [23] k-ε turbulence model given by where the non-dimensional measure of vortex stretching term, χ p , is computed by first performing the inner product of the two mean-rotation rate tensors and then by calculating the double inner product of the resultant and the mean-strain rate tensor. The coefficients for the Pope correction are equal to that of the standard k-ε model with the additional constant C ε3 = 0.79 providing the most agreeable results for the respective geometry and boundary condition tested in the original paper. The coefficients provided for these models are empirical relations that have been fit for generalized flow problems. Realizing this, along with the issues with the jet flow modelling, Tam and Ganesan [35] performed many jet flow experiments to determine a new set of coefficients that will allow accurate modelling of the jet flow for a variety of Reynolds numbers and geometries. This set of coefficients are given by C µ = 0.0874, C ε1 = 1.40, C ε2 = 2.02, σ k = 0.324, σ ε = 0.377, C ε3 = 0.822. It is noted that the diffusion coefficients σ k and σ ε are significantly smaller than the standard coefficients and can be characterized by an increase in diffusion of k and ε because of their location in the denominators of Equations (16) and (17). These coefficients will be used in the present work with the naming as Tam-Thies KE. Slight modifications to the empirical constants in the Launder and Spalding [23] are also used with "mod1" having the change of C ε1 = 1.50 similarly performed by Quinn and Militzer [11] and Faghani et al. [20] and "mod2" having the change of C ε1 = 1.52 and C ε2 = 1.90 similarly performed by Givi and Ramos [15].

k-ω Turbulence Models
The standard k-ω model is given by a production term of turbulent kinetic energy, k, and frequency of dissipation, ω. Kolmogorov, in 1941, was the first to suggest the twoequation turbulence model for kinetic energy and dissipation per unit turbulence kinetic energy [36]. Unaware of Kolmorgorov's work, Saffman [37] produced a k-ω model that had better agreement than Kolmorgorov's, and has served as the basis of "modern" k-ω turbulence models. One of the most influential contributors to the k-ω turbulence model was the work performed by David C. Wilcox. A sequential timeline of Wilcox's work will not be stated, and since no modifications to this model have been performed in the present work, the reader is encouraged to review the paper given in Wilcox [25] for details on the model.
The goal of the k-ω SST and BSL turbulence models is to use the k-ω model due to its robustness near the wall, and then transition to the k-ε model for its effectiveness in the far field. This is achieved using a blending function, F 1 , that is a function of wall distance with k-ω and "transformed" k-ε models. When the value of F 1 is equal to one, then the k-ω proposed by Wilcox [38] is used. When F 1 is equal to zero, then a transformed version of the standard k-ε model of Launder and Spalding [23] is used. It should be noted here that the "transformed" k-ε model is typically stated to be an exact transformation of the k-ε model, but this is erroneous to state because terms are dropped in the derivation. This will be discussed further in this paper.
To better understand this model and determine its ability to adapt to turbulent jet flow conditions, the k − ω BSL/SST turbulence models [39] with the updated modifications [40] will first be stated, followed by an in depth study on the origins of the model to determine if the coefficients are likely candidates for modification. The k − ω BSL model is given by where, γ, β, σ k , and σ ω are constants, φ, that are calculated using the equation φ = The k-ω model, φ 1 , constants are taken directly from Wilcox [38] to be σ k1 = 0.5, σ ω1 = 0.5, β 1 = 0.0750 and γ 1 = β 1 /β * − σ ω1 κ/ β * ∼ = 0.56 by ensuring that κ = 0.41. The constants for the "transformed" k-ε model, φ 2 , are given by σ k2 = 1.0, σ ω2 = 0.856, β 2 =0.0828, and γ 2 = β 1 /β * − σ ω2 κ/ β * ∼ = 0.44 with ν t = k/ω . These are used to mimic the constants given by the standard k-ε model of Launder and Spalding [23]. As previously stated, the k-ω BSL/SST turbulence model uses the standard 1988 Wilcox k-ω turbulence model in the near field and transitions to a "transformed" k-ε in the far field. The "transformed" k-ε equation is derived by a change in dependent variables from ε to ω using the relation defined by ε = C φ kω with the standard k-ε Launder and Spalding [23] model. Specifically, for kinetic energy the model becomes and for rate of dissipation If we then compare this fully derived transformed k-ε formulation with the model given from Menter [26], it is observed that the equation for turbulent kinetic energy is equivalent, but that some terms are dropped in the rate of dissipation equation. Specifically, fluid viscosity in the cross diffusion term given by and terms for turbulent diffusion given by In his book "Turbulence Modelling for CFD", Wilcox [14] states that "focusing on free shear flow, we can ignore molecular viscosity". Perhaps the contribution of molecular viscosity has an insignificant contribution with the application to free sheer flows, but this does not explain why the model can be applied to applications other than free shear flows. Wilcox also states that if we assume 1/σ ε = 1/σ k for simplicity, then the terms for turbulent diffusion given by Equation (28) can be dropped. It can, therefore, be argued that this is not an "exact" transformation of the standard k-ε model. That being said, this is not a problem, since the k-ε is not any more fundamental than the k-ω model; with the only difference in transformation method. Therefore, the modifications that have been performed for the standard k-ε model in jet flow, both in this work and by other researchers, do not necessarily carry over to the "transformed" k-ε model.
We can dive into this analysis further by relating the empirical coefficients between the standard and "transformed" k-ε turbulence models. It is realized that the coefficients chosen for the "transformed" k-ε model directly mimic that of the standard k-ε model for λ 2 and β 2 , where σ ω2 = 0.856 for the "transformed" and σ ω2derived = 1/σ ε = 0.769 for the derived version. It is still unclear whether this adjustment makes a significant contribution to a solution, but it should be noted nevertheless.
All this being said, if the standard k-ε model produces spreading rates that are not agreeable, then it makes sense that the k-ω BSL model will produce spreading rates similar to the standard k-ε and, therefore, some modification of the empirical coefficients need to be performed to provide "more" agreeable results. To directly mimic the modification of the the k-ε model, we set C ε1 = 1.5 and relate the constants from the k-ε to the transformed k-ε model by C ε1 − 1 = 0.5 = γ 2 . This will be aptly named mod1 for the k-ω SST/BSL models in the solution results. It should be noted that although other researchers have modified the empirical constants for the k-ε in the past [11,15,17,20], to the present authors' knowledge there has been no one that has modified the constants in k-ω SST/BSL models other than the developer of the turbulence model themselves. Therefore, results for the modified k-ω SST/BSL models should not be indicative as a general solution for a variety of jet flows unless further testing is performed.
The k − ω SST model is the most commonly used model proposed by Menter [39] that is differentiated from the BSL model through a change in coefficients, turbulence viscosity, and blending functions. The coefficients for the k − ω SST turbulence model are given by and turbulent viscosity is given by where, Ω is the vorticity magnitude calculated as Ω = √ 2 Ω ij : Ω ij , and y is the distance from the centre of the cell to the wall. It should be noted that because of the change in eddy viscosity relative to the BSL model (ν t = k ω), the production term in the dissipation equation for k-ω SST model is equal to γω k P instead of γ ν t P.
There have since been updates to the original [39] model given in Menter et al. [40] that provide more agreeable solutions in industrial flow applications and is considered the "standard" model, so we will name this model SST std. to avoid confusion. The most significant change is that of eddy viscosity, from the use of the magnitude of vorticity to the magnitude of the strain invariant.
where, S = 2S ij S ij = √ 2 S ij : S ij . A turbulent production limiter is also used to avoid turbulence production in the stagnation regions for both k and ω, as given by CD kω is changed by the use of 10 −10 instead of 10 −20 , and some constants are changed, specifically, γ 1 = 5/9 and γ 2 = 0.44. The default k-ω SST that is implemented in OpenFOAM is a modification of the original k − ω SST model. The origins of these modifications are unclear because of lack of a reference, but it is very similar to the original SST model and will be labelled SST def. in all the results.

Simulation Setup
Two distinct cases are considered that are separated by geometry and flow conditions. Namely, cases that follow the work performed by Hardalupas et al. [1] and Bogusławski and Popiel [27]. The operating parameters chosen for the simulations are taken directly from the inputs in the respective papers and if an assumption was needed then an explanation and justification will be shown.

Hardalupas et al. Test Case
Hardalupas et al. [1] use a phase-Doppler anemometer to measure the velocity and flux distributions in both clean jet (with seeding powders) and particle-laden high-speed air jets. The flow is fully developed in a 15 mm diameter nozzle of 93 nozzle diameters in length. An average velocity of 13 m/s is used that correlates to a Reynolds number of Re dh = 13,000.

Bogusławki and Popiel Test Case
The experiments performed in the study of Bogusławski and Popiel [27] use an air jet issued from a long straight tube of 38.5 mm in diameter and 50D in length. The velocity and turbulence characteristics are measured using a constant temperature thermo-anemometer. This paper defines u inf = 20 m/s that corresponds to a Reynolds number, Re dh = 51,000, but then uses this same u inf as the maximum centreline velocity at the nozzle exit. This is somewhat of an ambiguity because the mean free-stream velocity in a nozzle will not be the maximum centreline velocity and, therefore, to avoid confusion the free-stream velocity, in this case will be considered the maximum centreline velocity, u inf = u 0m = 20 m/s. The fixed value inlet boundary condition for the numerical simulation similarly follow the Hardalupas et al. [1] test case in that the velocity profile varies radially to the power law for turbulent jet flow with an exponent, n = 6.6. The boundary conditions were set to that of Table 3, where I = 0.16Re −1/8 dh ≈ 0.04 with Re dh = 50,000.

Model Setup
Two domains are created by taking a 5 degree wedge slice of a circular 3D domain. Figure 1 shows the geometry for both the Hardalupas et al. [1] and Bogusławski and Popiel [27] test cases. A nozzle length of 100D is used with a free-stream domain half width of 4D. A nozzle wall thickness is include in the domain with thickness D/6 and extended out towards the nozzle inlet by D/2. This is to potentially model any affects caused by boundary layer attachment at the nozzle wall. Discretization is performed using hexahedral cells generated in Ansys ® ICEM CFD. Uniformly space cells of size 0.00075, 0.0006, and 0.0005 are used with a single cell being used for the width in the wedge direction. The implicit second order backward scheme is used for temporal discretization. Open-FOAM second order linearUpwind scheme is used for all velocity and turbulent fields. This scheme uses upwind interpolation weights depending on the gradient of the fluid velocity [41]. The widely used pressure-based PISO (pressure-implicit split-operator) algorithm proposed by Issa [42] is used for the fluid solver. The OpenFOAM PCG solver with a GAMG pre-conditioner is used for pressure and PBiCG with DILU pre-conditioner used for all other fields [41]. Final solver tolerance is set to 1 × 10 −8 for all scalar and vector fields. Table 3 shows the boundary conditions for the Hardalupas et al. [1] test case. The eddy viscosity near the wall is calculated from a velocity based continuous wall function that follows the Spalding's equation [43]. The inlet velocity is set to the power law for fully developed turbulent pipe flow with the exponent, n, set equal to 6.6 and a maximum velocity of u max = 14.7. A "calculated" entry is used to signify that the turbulent viscosity is calculated from k, , and ω in the turbulence model. In Table 3, the turbulence intensity is equal to I = 0.16Re −1/8 dh ≈ 0.05 with Re dh = 1.3 × 10 4 . Entrainment boundary conditions are used for pressure and velocity at the domain outlet. Pressure is set to zero for outflow and set to −0.5 U 2 for inflow. Velocity is set to zero gradient at all times except the tangential component is set to zero. This is used to ensure that the inflow velocity can "find" the normal component, facilitating a proper entrainment boundary.
There are many ways to calculate an initial ω, such as ω = C −1/4 µ √ k l ∼ = 17.9, but the most certain way is to a run a k-ε simulation and use the free-stream results to calculate the initial free-stream ω for the k-ω with ω = ε kC µ . It was determined that an average value of 1597 will suffice for the initial free-stream values. It should be noted that this value is difficult to discern because of the boundary layer separation at the exit of the nozzle that may cause the k-ω Wilcox turbulence model to be inaccurate from sensitivity to free-stream input conditions [26]. The initial free-stream values were set to 0.1 for kinetic energy, k, that corresponds to a turbulent intensity of 0.05.

Wall Modelling Concerns and Mesh Resolution
It is the present authors' position that this work will serve as a basic for future work in multiphase particle-laden flows. Therefore, the meshing thought process utilized in this research has the consideration of limitations in mesh size relative to a particle diameter and the primary focus of this paper is the use of wall function approaches. In eddy viscosity RANS turbulence models, it is required to either fully integrate to the wall, or use wall function to estimate the near wall properties where viscous affects cannot be ignored. Fully integrating to the wall requires the first cell height to produce a y+ of about 1 or less, and the first 8-10 cells located below a y+ of about 11.5. This can cause issues with the addition of particles because of the momentum coupling between the two phases. Therefore, a good option is to use a wall function approach that will allow a larger cell near the wall with effective momentum transfer between the two phases.
The first wall functions were constructed using the Launder and Spalding methods in which it is required that the first cell height produce a y+ in the logarithmic layer (30 < y+ < 300) [23]. Newer adaptive wall functions have since been proposed, which allow a wide range of y+ values that include the more illusive buffer region. OpenFOAM uses adaptive wall function approaches proposed by Kalitzin et al. [44]. These wall functions were tested in the same work and it was realized that solutions were agreeable through a spread of y+ values from 11 to 111, with the highest error found at 5 < y+ < 11. This is the wall function used in the present work.
A zero gradient boundary condition is used for kinetic energy at the wall, with wall functions for both ε and ω. In OpenFOAM, the epsilon wall function provides a constraint at the wall for both turbulent kinetic energy, k, and dissipation, [45]. The Omega wall function uses a blending of the viscous and logarithmic layer with a blending function and, theoretically, can be accurate in the buffer region. That being said, much is still unknown about the buffer region that is between the viscous and log-layer and, therefore, it is the best practice to stay out of this region as much as possible. A continuous wall function for turbulent viscosity is used, which is modelled from Spalding's equation that fits a curve of the relationship between u+ and y+ [43] for the full range of the boundary layer. It should be noted that this is a function that is based on velocity or shear stress at the wall and not kinetic energy.

Grid Independence
A grid sensitivity analysis is performed using the Hardalupas et al. [1] test case to determine if there are errors associated with the discretization of the fluid domain into a mesh and that the solution is independent of the mesh cell size. Three systematically refined uniformly sized meshes of 26, 85, and 100 thousand hexahedral cells are generated in ICEM © CFD meshing software with a resulting y+ based on velocity of 27, 15, and 13, respectively.
The default k-ε turbulence model available in OpenFOAM is used for the grid independence with the simulation setup as described in Section 3 of this paper. It should be noted that the default k-ε turbulence model is the model proposed by Launder and Spalding [23]. Richardson extrapolation, as recommended by the Journal of Fluids Engineering described in Celik [46], was used to determine an extrapolated solution along a sampling line for both fluid velocity and turbulent kinetic energy ( Figure 2). The fine-grid convergence index, GCI f ine , is then calculated for each point and averaged along the sampling line. The GCI f ine for velocity has a value of 0.2% for the average and 0.8% for the maximum at 10D and 0.03% for the average and 0.2% for the maximum at 20D. The GCI f ine for kinetic energy has a value of 0.4% for the average and 4.4% for the maximum at 10D and 0.05% for the average and 0.1% for the maximum at 20D.
The mesh of 100, 000 was chosen for all proceeding simulations because it was the finest mesh available that allows a reasonable y+ value for accurate wall function modelling.

Fluid Velocity
The potential core of the jet is the region in which there is a small difference in the velocity of the bulk and jet regions. The core is characterized by a conical shape narrowing from the nozzle exit to the end of the developing region. After the potential core region, the flow develops into what is known as the fully established flow field [47]. Figure 3a illustrates the velocity as a function of axial non-dimensional distance from the nozzle exit for the Hardalupas et al. [1] test case. Although it is difficult to discern exactly where the potential core ends for the experimental data because of the limited number of sampling points in the axial direction, if an extrapolation line is fitted to follow the available points, we can comment that the potential core has a length that is the most similar to the standard k-ε, Pope Corr. k-ε, and k-ω BSL turbulence model. The k-ε mod1 model with C ε1 = 1.52 and C ε1 = 1.90 is the most similar throughout the entire range of velocities, but has a slightly longer potential core, as shown by Givi and Ramos [15]. The conclusion, therefore, is that if merely the modelling of the potential core is desired and not the fully established flow field, then a variety of turbulence models can be used with sufficient accuracy. Nonetheless, if the full range of the near field, from the potential core through to the fully established flow (0 ≥ x/D ≥ 20), is modelled, then the k-ε mod1 turbulence model is the most desired with the k-ω SST mod1 and k-ω BSL mod1 coming in a close second. This is shown for the axial and radial velocity profiles by observing Figure 3.
The 1/7th power law has been shown to be an accurate and simple model for profiles of fully developed turbulence pipe flow for a variety of Reynolds numbers and is plotted in orange for Figures 3b and 4b. It was observed that the flow profiles for the Hardalupas et al. [1] test cases follow the 1/7 th power law well for the inner region of the flow at the nozzle exit. However, the profile has a slight "U" shape, which is characteristic of a more laminar flow without a boundary layer. By observing Figure 3b it is realized that the numerical results agree more with the Hardalupas et al. [1] test case then the power law. The Bogusławski and Popiel [27] test case exhibits different behavior in that the velocity at the nozzle exit follows the 1/7th power law for both the experimental and numerical results.  It was originally thought that k-ω SST turbulence model, with its variants, would be the best option for these types of flow because of its ability to model the near field effectively using k-ω and the outer field using k-ε. This was found to be untrue because of the issues with empirical coefficients associated with these models. In the derivation of the k-ω SST/BSL turbulence models some terms are dropped and seem to be absorbed into the cross diffusion empirical coefficient of the transformed k-ε given by σ ω2 = 0.856 that are provided in the k-ω BSL/SST turbulence models.
It was found that the default k-ω SST model compared very well to the standard k-ε model in the far field. It is stated by Menter [26] that the blending function goes to zero at the log region. We observed that the y+ value for the mesh and Reynolds numbers was about 13 for the Hardalupas et al. [1] test case and 30 for the Bogusławski and Popiel [27] test case. With this y+, the wall functions are fully modelling the flow field up to, or very close to, the log region for the Bogusławski and Popiel [27] test case and into the buffer region for the Hardalupas et al. [1]. Hence, for the Hardalupas et al. [1], the flow field is calculated slightly with the k-ω turbulence model and the rest with the "transformed" k-ε, while the flow field is calculated almost exclusively from the "transformed" k-ε for the Bogusławski and Popiel [27] test case.
Hence, the flow field is calculated from the transformed k-ε turbulence model for the k-ω SST turbulence model. The derivation of the k-ω SST model shows that two terms are dropped, and this does affect the results slightly with the difference between the modified and k-ε equation. It is thought that the constants can be further modified to account for this slight change. No previous research has investigated and reported if modifying the constants of the k-ω SST model would gain any benefit in the predicted flow field, and an extensive study would be required to fully understand the behavior with this type of modification and this particular aspect was not researched in this study. That being said, sufficiently accurate results using the modified k-ε turbulence model were achieved, which can be used for future work in round turbulence dense CFD-DEM turbulence jet flows.
It should be noted that the stated inlet velocity for the experiments performed by Bogusławski and Popiel [27] is U ∞ = 20 m/s with the profile at the nozzle exit following the 1/7th power law given by Typically the 1/7th power law velocity is normalized with the maximum centreline velocity at the nozzle exit, U 0m , and not the infinite, or mean, free-stream velocity, U ∞ , and, therefore, the maximum velocity is slightly greater in magnitude than the normalized velocity. The data given from the experiments demonstrate that the definition of U ∞ is not the free-stream velocity but, actually, the maximum centreline velocity at the nozzle exit, that is U ∞ = U 0m , because the plots never go above unity. Therefore, going forward, we will normalize the numerical data with the maximum centreline velocity, U 0m . This data treatment allows consistent normalization across each of the turbulence models because the maximum centreline velocity changes slightly between turbulence models with a constant free-stream velocity value at the inlet.

Turbulence Properties
As outlined before, when Favre averaging the Navier-Stokes equations, an additional Reynolds stress term is added that accounts for the fluctuations in the flow field. The Reynolds stress is given as a tensor by u i u j =   u x u x u x u y u x u z u y u x u y u y u y u z u z u x u z u y u z u z   This is a symmetric tensor and so we have gained six additional scalar transport equations in addition to the Navier-Stokes. Instead of modelling each of these six transport equations, we choose a more heuristic approach by modelling the Reynolds stress directly with the Boussinesq approximation with a linear function of the mean-strain-rate tensor using an eddy viscosity that was previously stated by Equation (11) but will also be stated here for completeness where k is the turbulent kinetic energy defined by the one-half of trace of the Reynolds stress tensor k = 1 2 trace u i u j Therefore, the kinetic energy can be expressed by one-half of the sum of the variances, or one-half the mean square magnitude of the velocity fluctuations given as It should be noted that the eddy viscosity is isotropic but not necessarily homogeneous and, therefore, the turbulence itself is not necessarily isotropic because the mean-strain-rate tensor of the Navier-Stokes equations is not isotropic. This is an important distinction because eddy viscosity is merely independent of direction.
To calculate turbulence kinetic energy for the experiments we use the provided rootmean-square (RMS) components of velocity given in the Hardalupas et al. [1] paper's notation as u = RMS(u) = √ u 2 , which is the square root of the variance equivalent to the present author's paper notation as (u axial ) 2 and (u radial ) 2 . Only the radial and axial velocity components are given and, therefore, kinetic energy is in "cylindrical" coordinates and transformation can be challenging. One estimation is the so-called two-dimension approximation that is proposed by Rodriguez et al. [48] k = 3 4 u axial 2 + u radial 2 (39) where, to obtain Equation (39), we have a third component of turbulent kinetic energy distribution given by (u θ ) 2 = (u axial ) 2 + (u radial ) 2 0.5 . This is assuming that the turbulent kinetic energy values are the same magnitude in the θ direction, which is not necessarily the case. It is observed by Figure 3e-h that this approximation provides results that are agreeable for the k-ε mod1 turbulence model but it would be unwise to make a definite conclusion with the two-dimensional assumption and, therefore, no definite conclusion for kinetic energy can be realized purely from this test case.
The turbulent kinetic energy for both cases drastically changes from a "U" shape to a "W" shape between −0.1D to 0.1D as the flow leaves the nozzle exit. It is interesting to note that a "W" shape is achieved at the nozzle exit within the range of a single nozzle diameter, 1D, in the numerical results. However, the experimental results demonstrate a "W" shape, but outside the range of 1D. This is a very important observation because it appears that the location of the hump in the W shape is not significantly different between the numerical and experimental results further in the flow field.
From Figure 4e-h it is observed that the turbulent kinetic energy in the experimental results of Bogusławski and Popiel [27] follows the overall trend well but the magnitude is significantly lower than the numerical results for all of the turbulence models tested. This was also observed in the work performed in Faghani et al. [20] and Quinn and Militzer [11]. Contrasting to this, the use of the two-dimensional approximation causes a higher magnitude in results for the Hardalupas et al. [1] test case. It is thought that this approximation is more for swirling flows causing more turbulent kinetic energy generation in the θ direction. The turbulence intensity in the axial direction lined up well with the experimental data, but the intensity in the other directions was about half of those achieved with the numerical results. This can be expected because of the isotropic assumption of the linear viscous RANS turbulence models ensures that the turbulence modulation is equal in all directions. Therefore, it is up to the researcher to determine if the benefits of using this type of model outweigh the potential errors. This notion is true for the use of any linear viscous RANS turbulence model because our goal is to capture the average rather then the instantaneous flow behavior.

Conclusions
Obtaining accurate single-phase spreading rates for jets is critical before the inclusion of particles in a CFD-DEM simulations because, otherwise, it would be difficult to attribute an agreeable answer with confidence in the CFD-DEM simulation to the correct use of models and not a combination of models that merely work for one case but not necessarily another. With this in mind, single-phase simulations have been performed in this current work to obtain agreeable single-phase results within limitations of CFD-DEM before the inclusion of particles. It is the authors intention to use this work as a base and perspective for future work in fully coupled dense two-phase CFD-DEM jet flow simulations by addressing the limitations on geometry, mesh cell size, wall modelling, and momentum coupling that occur. It was realized that a slight modification of the empirical constant C ε1 = 1.5 in the standard Launder and Spalding [23] k-ε model provided the most agreeable results for both test cases for fluid velocity. A change in coefficients of the k-ω SST v1/BSL turbulence models that mimic the C ε1 = 1.5 turbulence model by setting γ = 0.5 provided suitable results as well. It is unclear whether RANS models in general will accurately model the turbulent kinetic energy in all directions for jet flow but, as a general conclusion, the kinetic energy from the numerical results will be higher in magnitude then the experiments because this is observed for the Bogusławski and Popiel [27] test case and the work performed previously by Faghani et al. [20] and Quinn and Militzer [11]. It is, therefore, up to the researcher to weigh the costs and benefits of using these models with considerations to help quantify the errors in kinetic energy along with whether merely the velocity field is desired.

Conflicts of Interest:
The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results