Next Article in Journal
Nanoparticle Delivery in Prostate Tumors Implanted in Mice Facilitated by Either Local or Whole-Body Heating
Next Article in Special Issue
Predictions of Vortex Flow in a Diesel Multi-Hole Injector Using the RANS Modelling Approach
Previous Article in Journal
Comparison of Semi-Empirical Single Point Wall Pressure Spectrum Models with Experimental Data
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

by
Dustin Steven Weaver
* and
Sanja Mišković
*
Norman B. Keevil Institute of Mining Engineering, Faculty of Applied Science, University of British Columbia, 517-6350 Stores Road, Vancouver, BC V6T 1Z4, Canada
*
Authors to whom correspondence should be addressed.
Submission received: 19 June 2021 / Revised: 23 July 2021 / Accepted: 28 July 2021 / Published: 31 July 2021
(This article belongs to the Collection Advances in Turbulence)

Abstract

:
This paper presents an analysis of linear viscous stress Favre averaged turbulence models for computational fluid dynamics (CFD) of fully turbulent round jets with a long straight tube geometry in the near field. 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 flow. Simulations are performed with a Reynolds number of 13,000 and 51,000 using two different tube diameters. It is found that a modified k- ε turbulence model achieved the most agreeable results for both the velocity and turbulent flow fields between these two flow regimes, while a modified k- ω SST/BSL also provided suitable results.

1. 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 the use of a nozzle in the solution domain and wall functions or eliminates the requirement for resolving the boundary layer. Although this is the best solution for fundamental single-phase flow problems, determining an optimum nozzle shape for particle-laden jets using CFD-DEM simulations requires that the nozzle profile shape be included in the solution domain. This is due in part to the location and initial operating parameters for particle insertion. Particles could be inserted at the nozzle exit, but this would require experimental data for each change in nozzle shape, which would be impractical for any type of industrial flow problem. Therefore, to tune these RANS turbulence models, a length of nozzle is included in the solution domain to ensure that the shear layer along the nozzle wall and the resulting boundary layer separation are being modelled correctly.
It is important to realize that with the current developments in computational resources, LES turbulence modelling is the ideal candidate for many different flow problems because of its ability to model anisotropic flow behaviors. However, for some applications, such as in fully coupled Eulerian–Lagrangian flow modelling, LES is impractical or impossible because of limitations related to mesh size for sufficient resolution of turbulent eddies and wall modelling. Therefore, for such applications, RANS turbulence models are a viable option. It should be noted that many studies have been performed using CFD-DEM and LES for jet flow, but a one-way coupled regime have been used. Although one-way coupled simulations may be adequate for dilute flows, particle–particle and particle–fluid interactions make a significant contribution to the physics of jets for dense flows [12].
RANS turbulence models are derived purely from dimensional analysis and some intuition in theory. Coefficients are then used to provide accurate results to a variety of general flow regimes. It would be ideal if these coefficients were set from principle values, such as kinetic theory of gases, but instead, they are empirically set for a variety of general flow regimes. This has been proven to be accurate in many cases but has shown inaccuracy in capturing the spreading rates for round jets [13,14].
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.

2. Numerical Methodology

2.1. 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
C ε 2 ε 2 k C ε 2 C ε 3 χ p ε 2 k
where, χ p is the non-dimension measure of vortex stretching given by
χ p = Ω i j Ω i k S k i ε / k 3
and Ω i j and S i j are the mean-rotation and mean-strain rate tensors defined as
Ω i j = 1 2 u i x j u j x i
S i j = 1 2 u i x j + u j x i
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: Ω i j = s k e w ( U ) and S i j = s y m m ( 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
β = β 0 f β
where
β 0 = 0.0708 , f β = 1 + 85 χ ω 1 + 100 χ ω , χ ω = Ω i j Ω j k S k i β ω 3
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 d h = 13,000 for the test case given by Hardalupas et al. [1] and Re d h = 51,000 for the test case given by Bogusławski and Popiel [27].

2.2. Numerical Model

The fluid flow in the present research is described by the Navier–Stokes equations given by
u j x j = 0
u i t + x j u j u i = x i p ρ 0 + x j ν m u i x j + u j x i 2 3 ν m u k x k δ i j
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
2 3 u k x k δ i j = 0
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 ¯ + u and p = p ¯ + p , and with some lengthy reductions the Navier–Stokes equations become the RANS equations given by
u ¯ i t + x j u ¯ j u ¯ i = x i p ¯ ρ 0 + x j ν m u ¯ i x j + u ¯ j x i 2 3 u ¯ k x k δ i j u i u j ¯
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
u i u j ¯ = ν t u i x j + u j x i 2 3 u k x k δ i j + 2 3 k δ i j
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.
u i t + x j u j u i = x i p ˜ ρ 0 + x j ν m + ν t u i x j + u j x i 2 3 u k x k δ i j
p ˜ = p ¯ + 2 3 ρ 0 k δ i j
where ν t is the turbulent viscosity and p ˜ 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
α l ρ l u l = 0
· ρ l α l u l u l = α l p + · α f τ l K s l u l u s + α / ρ l g + f
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 s l 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.

2.2.1. k ε Turbulence Models

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
k t + u j k x j = τ i j u j x j ε + x j ν + ν t σ k k x j + L k
ε t + u i ε x j = C ε 1 f 1 ε k τ i j u i x j C ε 2 f 2 ε 2 k + x j ν + ν t σ ε ε x j + L ε
where eddy viscosity is given by ν t = C μ f μ k 2 k 2 ε ε and the Boussinesq assumption is used to obtain
τ i j = ν t 2 S i j 2 3 u k x k δ i j 2 3 k δ i j
S i j = 1 2 u i x j + u j x i = s y m m ( U )
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 32 2 k 3 32 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
y = C y 1 Re y + C y 2 Re y 2 Re y = k y k y ν ν , C y 1 = 2.4 , C y 2 = 0.003
It should be noted that f μ 1 at y + = 300 . For the RNG turbulence model η = S k k ε ε , η 0 = 4.38 , S = 2 S i j S i j 0.5 .
For completeness, we will state the previously mentioned Pope correction to the standard Launder and Spalding [23] k- ε turbulence model given by
C ε 2 ε 2 k C ε 2 C ε 3 χ p ε 2 k , χ p = Ω i j Ω i k S k i ε / k 3
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].

2.2.2. 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 two-equation 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
D k D t = P β k ω + x j ν + σ k ν t k x j
D ω D t = γ ν t P β ω 2 + + x j ν + σ ω ν t ω x + 2 1 F 1 σ ω 2 1 ω k x j ω x j
F 1 = tanh arg 1 4 arg 1 = min max k 0.09 ω y , 500 ν y 2 ω , 4 σ ω 2 k C D k ω y 2 C D k ω = max 2 σ ω 2 1 ω k x j ω x j , 10 20
where, γ , β , σ k , and σ ω are constants, ϕ , that are calculated using the equation ϕ = F 1 ϕ 1 + 1 F 1 ϕ 2 . The k- ω model, ϕ 1 , constants are taken directly from Wilcox [38] to be σ k 1 = 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 σ k 2 = 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
D k D t = P ε + x j ν + ν t σ k k x j D k D t = P C ϕ k ω + x j ν + ν t σ k k x j
and for rate of dissipation
D ε D t = C ε 1 f 1 ε k P C ε 2 f 2 ε 2 k + x j ν + ν t σ ε ε x j D ω D t = C ε 1 f 1 1 ω k P C ϕ C ε 2 f 2 1 ω 2 + x j ν + ν t σ ε ω x 2 k ν + ν t σ ε k x j ω x j + ω k x j 1 σ ε 1 σ k ν t k x j
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
2 ν k k x j ω x j
and terms for turbulent diffusion given by
ω k x j 1 σ ε 1 σ k ν t k x j
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 σ ε σ ε = 1 1 σ k σ 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 σ ω 2 d e r i v e d = 1 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 m o d 1 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
σ k 1 = 0.85 , σ ω 1 = 0.5 , β 1 = 0.0750 , a 1 = 0.31 β = 0.09 , κ = 0.41 , γ 1 = β 1 β 1 β β σ ω 1 κ σ ω 1 κ β β 0.55
and turbulent viscosity is given by
ν t = a 1 k max a 1 ω , Ω F 2 F 2 = tanh arg 2 2 arg 2 = max 2 k 0.09 ω y , 500 ν y 2 ω
where, Ω is the vorticity magnitude calculated as Ω = 2 Ω i j : Ω i j , 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 k ω ω ), the production term in the dissipation equation for k- ω SST model is equal to γ ω ω k k P instead of γ γ ν t ν 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.
ν t = a 1 k max a 1 ω , S F 2
where, S = 2 S i j S i j = 2 S i j : S i j . A turbulent production limiter is also used to avoid turbulence production in the stagnation regions for both k and ω , as given by
P = min P , 10 β ρ ω k
C D k ω is changed by the use of 10 10 instead of 10 20 ,
C D k ω = max 2 σ ω 2 1 ω k x j ω x j , 10 10
and some constants are changed, specifically, γ 1 = 5 59 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.

3. 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.

3.1. 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 d h = 13,000.

3.2. 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 50 D 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 d h = 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 0 m = 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.16 Re d h 1 / 8 0.04 with Re d h = 50,000.

3.3. 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. OpenFOAM 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.16 Re d h 1 / 8 0.05 with Re d h = 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 k l 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 ω = ε ε k C μ k C μ . 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 .

3.4. 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.

4. Results and Discussion

4.1. 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, G C I f i n e , is then calculated for each point and averaged along the sampling line. The G C I f i n e 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 G C I f i n e 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.

4.2. 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 x D 20 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 Figure 3b and Figure 4b. It was observed that the flow profiles for the Hardalupas et al. [1] test cases follow the 1/7th 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
U U U = 1 2 r 2 r D D U = 1 2 r 2 r D D 1 7
Typically the 1/7th power law velocity is normalized with the maximum centreline velocity at the nozzle exit, U 0 m , 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 0 m , because the plots never go above unity. Therefore, going forward, we will normalize the numerical data with the maximum centreline velocity, U 0 m . 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.

4.3. 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
u i u j ¯ = ν t u i x j + u j x i 2 3 k δ i j
where k is the turbulent kinetic energy defined by the one-half of trace of the Reynolds stress tensor
k = 1 2 t r a c e 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
k = 1 2 u i u i ¯ = 1 2 u x 2 ¯ + u y 2 ¯ + u z 2 ¯
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 root-mean-square (RMS) components of velocity given in the Hardalupas et al. [1] paper’s notation as u = R M S u = u 2 , which is the square root of the variance equivalent to the present author’s paper notation as u a x i a l 2 ¯ and u r a d i a l 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 a x i a l 2 + u r a d i a l 2
where, to obtain Equation (39), we have a third component of turbulent kinetic energy distribution given by u θ 2 = u a x i a l 2 + u r a d i a l 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.

5. 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.

Author Contributions

Conceptualization, D.S.W. and S.M.; methodology, D.S.W.; software, D.S.W.; validation, D.S.W.; formal analysis, D.S.W.; investigation, D.S.W.; resources, S.M.; data curation, D.S.W.; writing—original draft preparation, D.S.W.; writing—review and editing S.M. and D.S.W.; visualization, D.S.W. and S.M.; supervision, S.M.; project administration, S.M.; funding acquisition, S.M. All authors have read and agreed to the published version of the manuscript.

Funding

The research is supported and funded by the MITACS Accelerate program and Disa Technologies (contract No. IT23227). We gratefully acknowledge the computing support of Compute Canada via their WestGrid program and ARC group at the University of British Columbia for providing the access to the Sockeye cluster.

Data Availability Statement

The cases, models, and data solutions are available at https://0-doi-org.brum.beds.ac.uk/10.6084/m9.figshare.15039798.v1, accessed on 19 June 2021.

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.

References

  1. Hardalupas, Y.; Taylor, A.; Whitelaw, J. Velocity and particle-flux characteristics of turbulent particle-laden jets. Proc. R. Soc. Lond. A Math. Phys. Sci. 1989, 426, 31–78. [Google Scholar] [CrossRef]
  2. Lau, T.C.; Nathan, G.J. Influence of Stokes number on the velocity and concentration distributions in particle-laden jets. J. Fluid Mech. 2014, 757, 432–457. [Google Scholar] [CrossRef] [Green Version]
  3. Modarress, D.; Wuerer, J.; Elghobashi, S. An Experimental Study of a Turbulence Round Two-Phase Jet. Chem. Eng. Commun. 1984, 28, 341–354. [Google Scholar] [CrossRef]
  4. Lau, T.C.W.; Nathan, G.J. The effect of Stokes number on particle velocity and concentration distributions in a well-characterised, turbulent, co-flowing two-phase jet. J. Fluid Mech. 2016, 809, 72–110. [Google Scholar] [CrossRef] [Green Version]
  5. Fan, J.; Zhang, X.; Cen, K. Laser Diffraction Method Measurements of Particle-Gas Dispersion Effects in a Coaxial Jet. Aerosol Sci. Technol. 1997, 26, 447–458. [Google Scholar] [CrossRef] [Green Version]
  6. Prevost, F.; Boree, J.; Nuglisch, H.J.; Charnay, G. Measurements of fluid/particle correlated motion in the far field of an axisymmetric jet. Int. J. Multiph. Flow 1996, 22, 685–701. [Google Scholar] [CrossRef]
  7. Sheen, H.; Jou, B.; Lee, Y. Effect of particle size on a two-phase turbulent jet. Exp. Therm. Fluid Sci. 1994, 8, 315–327. [Google Scholar] [CrossRef]
  8. Shuen, J.S. A Theoretical and Experimental Investigation of Dilute Particle-Laden Turbulent Gas Jets (Two-Phase Flow). Ph.D. Thesis, The Pennsylvania State University, State College, PA, USA, 1984. [Google Scholar]
  9. Tsuji, Y.; Morikawa, Y.; Tanaka, T.; Karimine, K.; Nishida, S. Measurement of an axisymmetric jet laden with coarse particles. Int. J. Multiph. Flow 1988, 14, 565–574. [Google Scholar] [CrossRef]
  10. Ashforth-Frost, S.; Jambunathan, K.; Whitney, C. Velocity and turbulence characteristics of a semiconfined orthogonally impinging slot jet. Exp. Therm. Fluid Sci. 1997, 14, 60–67. [Google Scholar] [CrossRef]
  11. Quinn, W.; Militzer, J. Effects of nonparallel exit flow on round turbulent free jets. Int. J. Heat Fluid Flow 1989, 10, 139–145. [Google Scholar] [CrossRef]
  12. Faeth, G. Mixing, transport and combustion in sprays. Prog. Energy Combust. Sci. 1987, 13, 293–345. [Google Scholar] [CrossRef] [Green Version]
  13. Pope, S.B. An explanation of the turbulent round-jet/plane-jet anomaly. AIAA J. 1978, 16, 279–281. [Google Scholar] [CrossRef]
  14. Wilcox, D.C. Turbulence Modeling for CFD, 3rd ed.; Number Book, Whole; DCW Industries: La Cãnada, CA, USA, 2006. [Google Scholar]
  15. Givi, P.; Ramos, J. On the calculation of heat and momentum transport in a round jet. Int. Commun. Heat Mass Transf. 1984, 11, 173–182. [Google Scholar] [CrossRef]
  16. Kim, K.Y.; Chung, M.K. Calculation of a strongly swirling turbulent round jet with recirculation by an algebraic stress model. Int. J. Heat Fluid Flow 1988, 9, 62–68. [Google Scholar] [CrossRef]
  17. Morgans, R.C.; Dally, B.B.; Nathan, G.J.; Lanspeary, P.V.; Fletcher, D.F. Applications of the Revised Wilcox 1998 K-Omega Turbulence Model to a Jet in Co-Flow. In Proceedings of the 2nd International Conference on CFD in the Minerals and Process Industries, Melbourne, Australia, 6–8 December 1999; pp. 479–484. [Google Scholar]
  18. Post, S.; Iyer, V.; Abraham, J. A Study of Near-Field Entrainment in Gas Jets and Sprays Under Diesel Conditions. J. Fluids Eng. 2000, 122, 385–395. [Google Scholar] [CrossRef]
  19. Georgiadis, N.J.; Yoder, D.A.; Engblom, W.A. Evaluation of Modified Two-Equation Turbulence Models for Jet Flow Predictions. AIAA J. 2006, 44, 3107–3114. [Google Scholar] [CrossRef] [Green Version]
  20. Faghani, E.; Saemi, S.D.; Maddahian, R.; Farhanieh, B. On the effect of inflow conditions in simulation of a turbulent round jet. Arch. Appl. Mech. 2011, 81, 1439–1453. [Google Scholar] [CrossRef]
  21. Wang, C.; Barghi, S.; Zhu, J. Hydrodynamics and reactor performance evaluation of a high flux gas-solids circulating fluidized bed downer: Experimental study. AIChE J. 2014, 60, 3412–3423. [Google Scholar] [CrossRef]
  22. Bonelli, F.; Viggiano, A.; Magi, V. High-speed turbulent gas jets: An LES investigation of Mach and Reynolds number effects on the velocity decay and spreading rate. Flow Turbul. Combust. 2021. [Google Scholar] [CrossRef]
  23. Launder, B.E.; Spalding, D.B. The numerical computation of turbulent flows. Comput. Methods Appl. Mech. Eng. 1974, 3, 269–289. [Google Scholar] [CrossRef]
  24. Rubel, A. On the vortex stretching modification of the k-epsilon turbulence model—Radial jets. AIAA J. 1985, 23, 1129–1130. [Google Scholar] [CrossRef]
  25. Wilcox, D.C. Formulation of the k-w Turbulence Model Revisited. AIAA J. 2008, 46, 2823–2838. [Google Scholar] [CrossRef] [Green Version]
  26. Menter, F.R. Influence of freestream values on k-omega turbulence model predictions. AIAA J. 1992, 30, 1657–1659. [Google Scholar] [CrossRef]
  27. Bogusławski, L.; Popiel, C.O. Flow structure of the free round turbulent jet in the initial region. J. Fluid Mech. 1979, 90, 531–539. [Google Scholar] [CrossRef]
  28. Weaver, D.; Miskovic, S. Analysis of Coupled CFD-DEM Simulations in Dense Particle-Laden Turbulent Jet Flow. In Proceedings of the ASME 2020 Fluids Engineering Division Summer Meeting, 13–15 July 2020. Virtual, Online. [Google Scholar] [CrossRef]
  29. Peng, Z.; Doroodchi, E.; Luo, C.; Moghtaderi, B. Influence of void fraction calculation on fidelity of CFD-DEM simulation of gas-solid bubbling fluidized beds. AIChE J. 2014, 60, 2000–2018. [Google Scholar] [CrossRef]
  30. Chien, K.Y. Predictions of Channel and Boundary-Layer Flows with a Low-Reynolds-Number Turbulence Model. AIAA J. 1982, 20, 33–38. [Google Scholar] [CrossRef]
  31. Arbiter, N. Development and scale-up of large flotation cells. Min. Eng. 2000, 52, 28. [Google Scholar]
  32. Launder, B.E.; Sharma, B.I. Application of the energy-dissipation model of turbulence to the calculation of flow near a spinning disc. Lett. Heat Mass Transf. 1974, 1, 131–137. [Google Scholar] [CrossRef]
  33. Thies, A.T.; Tam, C.K.W. Computation of turbulent axisymmetric and nonaxisymmetric jet flows using the K-epsilon model. AIAA J. 1996, 34, 309–316. [Google Scholar] [CrossRef]
  34. Wallin, S.; Johansson, A.V. An explicit algebraic Reynolds stress model for incompressible and compressible turbulent flows. J. Fluid Mech. 2000, 403, 89–132. [Google Scholar] [CrossRef]
  35. Tam, C.K.W.; Ganesan, A. Modified kappa-epsilon Turbulence Model for Calculating Hot Jet Mean Flows and Noise. AIAA J. 2004, 42, 26–34. [Google Scholar] [CrossRef]
  36. Kolmogorov, A.N. Equations of turbulent motion in an incompressible fluid. Dokl. Akad. Nauk SSSR 1941, 30, 299–303. [Google Scholar]
  37. Saffman, P.G. A Model for Inhomogeneous Turbulent Flow. Proc. R. Soc. Lond. Ser. A Math. Phys. Sci. 1970, 317, 417–433. [Google Scholar] [CrossRef]
  38. Wilcox, D.C. Reassessment of the scale-determining equation for advanced turbulence models. AIAA J. 1988, 26, 1299–1310. [Google Scholar] [CrossRef]
  39. Menter, F.R. Two-equation eddy-viscosity turbulence models for engineering applications. AIAA J. 1994, 32, 1598–1605. [Google Scholar] [CrossRef] [Green Version]
  40. Menter, F.R.; Kuntz, M.; Langtry, R. Ten Years of Industrial Experience with the SST Turbulence Model. Heat Mass Transf. 2003, 4, 625–632. [Google Scholar]
  41. Weller, H.G.; Tabor, G.; Jasak, H.; Fureby, C. A tensorial approach to computational continuum mechanics using object-oriented techniques. Comput. Phys. 1998, 12, 620. [Google Scholar] [CrossRef]
  42. Issa, R. Solution of the implicitly discretised fluid flow equations by operator-splitting. J. Comput. Phys. 1986, 62, 40–65. [Google Scholar] [CrossRef]
  43. Spalding, D.B. A Single Formula for the “Law of the Wall”. J. Appl. Mech. 1961, 28, 455–458. [Google Scholar] [CrossRef]
  44. Kalitzin, G.; Medic, G.; Iaccarino, G.; Durbin, P. Near-wall behavior of RANS turbulence models and implications for wall functions. J. Comput. Phys. 2005, 204, 265–291. [Google Scholar] [CrossRef]
  45. Popovac, M.; Hanjalic, K. Compound Wall Treatment for RANS Computation of Complex Turbulent Flows and Heat Transfer. Flow Turbul. Combust. 2007, 78, 177–202. [Google Scholar] [CrossRef] [Green Version]
  46. Celik, I.B. Procedure for Estimation and Reporting of Uncertainty Due to Discretization in CFD Applications. J. Fluids Eng. 2008, 130, 078001. [Google Scholar] [CrossRef] [Green Version]
  47. Rajaratnam, N. Chapter 2 The Circular Turbulent Jet. In Turbulent Jets; Developments in Water Science; Elsevier: Amsterdam, The Netherlands, 1976; Volume 5, pp. 27–49. ISSN 0167-5648. [Google Scholar] [CrossRef]
  48. Rodriguez, G.; Weheliye, W.; Anderlei, T.; Micheletti, M.; Yianneskis, M.; Ducci, A. Mixing time and kinetic energy measurements in a shaken cylindrical bioreactor. Chem. Eng. Res. Des. 2013, 91, 2084–2097. [Google Scholar] [CrossRef]
Figure 1. Slice of the 5 degree wedge simulation domain. D is the diameter of the nozzle, T is the length of the nozzle, and L is equal to 20D for Hardalupas et al. [1] test case and 12D for Bogusławski and Popiel [27] test case.
Figure 1. Slice of the 5 degree wedge simulation domain. D is the diameter of the nozzle, T is the length of the nozzle, and L is equal to 20D for Hardalupas et al. [1] test case and 12D for Bogusławski and Popiel [27] test case.
Fluids 06 00271 g001
Figure 2. Solutions for fine, medium, and coarse grids and Richardson extrapolation at (a,b) 10D and (c,d) 20D normal planes.
Figure 2. Solutions for fine, medium, and coarse grids and Richardson extrapolation at (a,b) 10D and (c,d) 20D normal planes.
Fluids 06 00271 g002
Figure 3. Test case Hardalupas et al. [1] normalized velocity profiles at (a) axial and (b) 0.1D, (c) 10D, and (d) 20D normal planes, and normalized turbulent kinetic energies at (e) axial, and (f) 0.1D, (g) 10D, and (h) 20D normal planes.
Figure 3. Test case Hardalupas et al. [1] normalized velocity profiles at (a) axial and (b) 0.1D, (c) 10D, and (d) 20D normal planes, and normalized turbulent kinetic energies at (e) axial, and (f) 0.1D, (g) 10D, and (h) 20D normal planes.
Fluids 06 00271 g003
Figure 4. Test case Bogusławski and Popiel [27] normalized velocity profiles at (a) axial and (b) 0D, (c) 6D, and (d) 12D normal planes, and normalized turbulent kinetic energies at (e) axial, and (f) 0D, (g) 6D, and (h) 12D normal planes.
Figure 4. Test case Bogusławski and Popiel [27] normalized velocity profiles at (a) axial and (b) 0D, (c) 6D, and (d) 12D normal planes, and normalized turbulent kinetic energies at (e) axial, and (f) 0D, (g) 6D, and (h) 12D normal planes.
Fluids 06 00271 g004
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.
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.
StudyReynolds NumberGeometryNozzle TypeTurbulence ModelsKey Features
Pope [13]NAAxi-SymmetricRound jetModified k-epsilonSuggested changing
C ε 1 = 1.6 and proposed vortex stretching term
Givi and Ramos [15]373,000Axi-SymmetricRound jetModified standard k- ε and k- ω turbulence modelsTested k-epsilon constants and C ε 1 = 1.52 with
C ε 2 = 1.90 provided most agreeable spreading rates
Kim and Chung [16]11,600Axi-SymmetricEdge fixed value set from experimentsStandard k-epsilon, algebraic stress modelHighly agreement results for swirling jets
Morgans et al. [17]50,000Axi-SymmetricFully developed straight tube round jet nozzleModified k-epsilon
( C ε 1 = 1.6), k- ω Wilcox
Comparable results between k- ω Wilcox and modified k-epsilon
Quinn and Militzer [11]208,000Axi-SymmetricSharp Edged and Smoothly ContouredModified k-epsilon
( C ε 1 = 1.48)
Comparable results between experiments for velocity but not turbulence quantities.
Post et al. [18]10,000Axi-SymmetricEdge fixed valueLaunder and Spalding, WFAcheived results of lower entertainment from experiment
Georgiadis et al. [19]610,950Axi-SymmetricConvergent contractionMany two-equation RANS turbulence modelsTam-Thies k-epsilon model demonstrated greatest agreement
Faghani et al. [20]20,800Axi-SymmetricSharp edge, round, and contoured nozzles with edge fixed value set to experimental data nozzle exit conditionsModified k-epsilon with C ε 1 = 1.5Significant differences in velocity, turbulence intensity, and turbulence kinetic energy with a change in nozzle shape
Wang et al. [21]20,000Impingement Three-DimensionalRound jet with fully developed flow inletMany two-equation RANS turbulence modelsk- ω SST demonstrated greatest agreement
Table 2. k- ε turbulence models that are used in the present study. Re t = k 2 k 2 ν ε ν ε .
Table 2. k- ε turbulence models that are used in the present study. Re t = k 2 k 2 ν ε ν ε .
TermLaunder and Spalding [23], mod1, mod2Chien [30]RNG Arbiter [31]Launder and Sharma [32]Thies and Tam [33]
C ε 1 1.44, 1.50, 1.521.351.421.441.40
C ε 2 1.92, 1.92, 1.901.801.681.922.02
f μ 1.0 1 e 0.0115 d + 1.0 exp 3.4 1 + R e t R t 50 50 2 1.0
f 1 1.01.01.01.01.0
f 2 1.0 1 0.4 1.8 e Re t Re T 36 36 C ε 2 + C μ η 3 1 η η η 0 η 0 1 + β η 3 1 0.3 exp Re t 2 1.0
σ k 1.01.00.71941.00.324
σ ε 1.31.31.30.71940.377
L k 0.0 2 ν k d 2 0.0 2 ν k y 2 0.0
L ε 0.0 2 ν ε d 2 e d + d + 2 2 0.0 2 ν ν t 2 u y 2 2 0.0
Table 3. Boundary conditions for Hardalupas et al. [1].
Table 3. Boundary conditions for Hardalupas et al. [1].
LocationVelocityPressureEddy ViscosityKinetic EnergyEpsilonOmega
WallNo SlipZero GradientSpalding Wall FunctionZero GradientEpsilon Wall FunctionOmega Wall Function
Inlet u max 1 r R 1 / n Zero GradientCalculated 1.5 I u 2 C ω 3 / 4 k 3 / 2 k 3 / 2 l l k 0.5 k 0.5 C μ 0.25 l C μ 0.25 l
OutletEntrainment VelocityTotal PressureCalculatedZero GradientZero GradientZero Gradient
Initial free-streamUniform 0.0Uniform 0.00 k = 1.5 I u 2 C ω 3 / 4 k 3 / 2 k 3 / 2 l l 1597
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Weaver, D.S.; Mišković, S. A Study of RANS Turbulence Models in Fully Turbulent Jets: A Perspective for CFD-DEM Simulations. Fluids 2021, 6, 271. https://0-doi-org.brum.beds.ac.uk/10.3390/fluids6080271

AMA Style

Weaver DS, Mišković S. A Study of RANS Turbulence Models in Fully Turbulent Jets: A Perspective for CFD-DEM Simulations. Fluids. 2021; 6(8):271. https://0-doi-org.brum.beds.ac.uk/10.3390/fluids6080271

Chicago/Turabian Style

Weaver, Dustin Steven, and Sanja Mišković. 2021. "A Study of RANS Turbulence Models in Fully Turbulent Jets: A Perspective for CFD-DEM Simulations" Fluids 6, no. 8: 271. https://0-doi-org.brum.beds.ac.uk/10.3390/fluids6080271

Article Metrics

Back to TopTop