Numerical Study of Velocity and Mixture Fraction Fields in a Turbulent Non-reacting Propane Jet Flow Issuing into Parallel Co-Flowing Air in Isothermal Condition through OpenFOAM

1. Introduction

Turbulent jets emerging from axisymmetric nozzles are pivotal flows classified as free shear layers. The dominating phenomenon is the turbulence due to the velocity gradient and several instabilities between the jet and the surrounding fluid [1]. Besides, re-acting flows are regarded as beneficial branches of turbulent jets, which have vast industrial usages in combustion, military industries engine and so on. Tackling these kinds of jets are generally intricate, because of interactions between turbulent mixing and heat release (due to chemical reactions) [2,3].

Although numerical modeling of turbulent reacting flows is substantially founded on models expanded for non-reacting constant-density flows [3], the use of such models may not be sufficiently authentic [4,5,6]. Having simplified such problem by a turbulent, non-reacting variable-density jet, it is simultaneously conceivable to preserve the intricacy of variable density, and eliminate the elaboration of coupling between turbulent mixing and chemical heat release [3]. In this regard, numerous experimental studies have been carried out on constant-density kinds of turbulent jets [7,8,9,10,11]. Moreover, extensive experiments have been conducted in the scope of strongly variable-density jets through axisymmetric nozzles [12,13,14,15,16,17,18,19,20]. Furthermore, by the advent of progressive laser instruments, it has been feasible to deeply scrutiny the re-acting type of turbulent jets [21,22].

In recent years, by the development of potent computers, the studies have been inclined toward the numerical simulations, due to be cheaper, faster and easier. In this scope, utilization of open source codes has been rapidly more popular, especially the OpenFOAM which has been written in C++ language and allows users to upgrade and develop the system further [1,23].

It should be noticed that the behavior of variable-density turbulent jets is intensively sensitive to the inlet conditions, particularly the injection ratio [24], the direction of the co-flowing and the geometry of the nozzle [2].

In this work, a numerical study has been carried out on a turbulent axisymmetric jet using Realizable k-ε eddy viscosity turbulence model through manipulating and regulating ReactingFoam solver in OpenFOAM 5x, with the aim of investigation of the capability of this procedure in simulating the experimental data which have been presented here: http://www.sandia.gov/TNF/DataArch/ProJet.html. It should be emphasized that ReactingFoam is one of the standard solvers in OpenFOAM which has been particularly extended to resolve the mixing problems of compressible flows involving combustion and reactions.

To simulate this problem, initially, sandiaD_LTS tutorial (in combustion>reactingFoam>RAS) has been regarded as the basement and then has been manipulated through deactivating combustion and reaction from corresponding files due to lack of these phenomena in the experiment.

2. Numerical basement

2.1. Mixture fraction theory

In contrast with premixed systems, in non-premixed combustion, fuel and oxidizer enter the reaction zone in distinct streams [25]. Under certain assumptions, the thermochemistry can be reduced to a single parameter: the mixture fraction, f, the mass fraction that originated from the fuel stream, which is a conserved scalar quantity, and therefore its governing transport equation does not have a source term. Consequently, combustion is simplified to a mixing problem. After mixing, the chemistry can be modeled through Equilibrium model, Steady Diffusion Flamelet or Unsteady Diffusion Flamelet, for being in chemical equilibrium, near that or being meaningfully far from that, respectively [25].

Although the assumption of equal diffusivities of species is problematic for laminar flows, it is generally acceptable for turbulent flows where turbulent convection overwhelms molecular diffusion. Consequently, the species equations can be reduced to a single transport equation for the mixture fraction, f. Non-premixed modeling involves the solution of this equation. Afterward, species concentration can be elicited from the predicted mixture fraction field. Moreover, under a particular set of simplifying assumptions, the instantaneous thermochemical state of the fluid is related to mixture fraction, [25].

The power of the mixture fraction modeling approach is that the chemistry is reduced to one or two conserved mixture fractions. Under the assumption of chemical equilibrium, all thermochemical scalars (species fractions, density, and temperature) are uniquely related to the mixture fraction(s) [25].

2.2. The transport equation for the mixture fraction

The Favre mean (density-averaged) mixture fraction equation is [25]:

$$\frac{\partial (\rho \overline{f})}{\partial t}+\nabla .(\rho \overrightarrow{v}\overline{f})=\nabla .(\rho \overrightarrow{v}\overline{f})+\frac{\partial }{\partial x_j}(\frac{{\mu }_l+{\mu }_t}{{\sigma }_t}\nabla \overline{f})$$

3. Description of computational domain and assumptions

By comparing the inner and outer dimensions of the tube with the aspects of the horizontal cross section of the experimental domain, and regarding the immense differentiation between the bulk velocity of Propane and co-flowing air, the assumption of an axisymmetric problem is not far from the reality. By this assumption, it is conceivable to massively alleviate the needed computational resources. Thus, the computational domain has been confined to a single 2D plane, i.e., a 5-degree sector of the biggest cylinder which can be embedded into the real domain.

There are two strategies for preparing the computational geometry. At first one, a single sector can be defined by neglecting the length of the nozzle. Thus, the flow within the nozzle will not be simulated. In this condition, the inlets boundaries of the Propane and the co-flowing air will be placed at the level of the nozzle exhaust.

In the other approach (as this study), a nozzle (for example, with 10 cm length) can be placed at the upstream to better establish of fully developed flow within the nozzle (before the injection into the principal computational domain), and to better form the boundary layer between the co-flowing air and the outer edge of the nozzle, a phenomenon which has been observed in the experiment. It is evident that the first approach is cheaper and faster because of reduced cell numbers, and presumably less accurate due to attenuate the simulation of the real boundary conditions around the nozzle's body.

Contrary to the Standard k-ε model, in the Realizable k-ε model, C_µ, (the coefficient which is imposed in the equation of turbulent viscosity, ν_t, is not a constant, but a function. It is worthy to notice that although the transport equation of turbulence kinetic energy, k, is the same as that in the Standard and the RNG k-ε models, except for the model constants, the transport equation of its Dissipation Rate , ε, has been modified by extracting the aforementioned equation from the exact equation of the transport of the mean-square vorticity fluctuations. Hence, the Realizable k-ε more accurately predicts the spreading rate of both planar and round jets and presents the best proficiency in comparison with all versions of the k-ε models [26].

By regarding the plotted figures of experimental data and with the aim of more minifying the domain, the length of the computational field has been confined to 100 cm (x/D=190), because no cognizable variations have been observed after x/D>80, in experimental data.

The volume fraction of N2 and O2 in the co-flowing air has been assumed as 0.7632 and 0.2368, respectively. Furthermore, the flow assumed to be compressible and isothermal (294 K), and the injecting fluids are considered to be ideal gasses.

To define the computational domain as a sector and its boundaries, the blockMeshDict, the changeDictionaryDict and the extrudeMeshDict files in the system folder have been utilized. Moreover, inletfuel (inlet boundary of Propane at upstream on the inner diameter of the nozzle), inletair (inlet boundary of co-flowing air at upstream, from outer ridge of the nozzle to the lateral border, leftside), outlet (the outlet boundary at downstream), axis (the axis of symmetry which the sector rotates 5 degrees around that), leftside (the lateral edge of domain), burnerwall (part of upstream boundaries which represents the body of the nozzle), front and back (to simulate a 2D problem) have been defined and adjusted as the boundaries through the blockMeshDict file and the corresponding files of each variable in the 0 folder.

Figure 1 shows a simplified scheme of the computational domain:

Eventually, in order to reduce the computational expenditures, and moreover, to obtain a proper aspect ratio, the grids have manipulated through predefined lengthGrading for both X and Y directions of each block, in blockMeshDict file.

To obtain turbulence kinetic energy, k, we have:

$$\begin{array}{l}{U}_{RMS}=\sqrt{\overline{{u^{\prime }}^2}},V_{RMS}=\sqrt{\overline{{v^{\prime }}^2}},W_{RMS}=\sqrt{\overline{{w^{\prime }}^2}}\\ \Rightarrow U_{RMS}^2=\overline{{u^{\prime }}^2},V_{RMS}^2=\overline{{v^{\prime }}^2},W_{RMS}^2=\overline{{w^{\prime }}^2}\\ k=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.(\overline{{u^{\prime }}^2}+\overline{{v^{\prime }}^2}+\overline{{w^{\prime }}^2})=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.(U_{RMS}^2+V_{RMS}^2+W_{RMS}^2)\end{array}$$

But, due to the 2D assumption which can be observed in experimental data, it is can be reduced to:

$$k=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.(U_{RMS}^2+V_{RMS}^2+W_{RMS}^2)=\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2$}\right.(U_{RMS}^2+2V_{RMS}^2)$$

Having plotted some non-dimensional parameters, assessment of axial and radial profiles of mean velocities, turbulence energy, mean mixture fraction, and the mixture fraction half radius, and subsequently, the comparison have been fulfilled between corresponding numerical and experimental data.

Table 1 comprises a concise list of computational domain's dimensions and some inlet data:

4. Initial and boundary conditions of variables

Regarding the provided meshes, it is evident that the last computational node on the leftside will be placed out of the boundary layer. In other words, Yplus can't be confined between 30 and 100, which is essential for the correctness of assumption of utilizing the wall functions. Consequently, it is not allowed to use any wall functions. On the other side, due to numerical costs' considerations, at least, it is not logical to refine the mesh toward the outer walls, i.e., leftside.

Thus, in the first approach, by assuming the boundary condition as fixedValue(=9.2) for the velocity, U, and zeroGradient for other variables, it is tried to simulate slip condition on the leftside boundary, which is not far from the reality. Besides, it is assumed that the meshes are refine enough on the inner and outer sides of the nozzle's body, i.e., burnerwall. Afterwards, wall function has been imposed for alphat ,epsilon ,k and nut, and noSlip condition for U.

Moreover, to interrogate the authenticity of such assumptions on nozzle's body, i.e., burnerwall, in some other test cases (the second approach), the mentioned boundary has been partitioned to its three sides: burnerwall_jet, burnerwall_air and burnerwall_upper, and each one of them have been assigned its own boundary conditions for each of the variables. Although the corresponding numerical results have not been presented in this study, but it could not be perceived any significant differentiations, which may approve the assumptions of first's approach.

Table 2 and Table 3 comprise a brief list of initial and boundary conditions for diverse variables, for the first approach. But, similar information for the second one has not been presented here for abbreviation:

It should be noticed that all of the variables have been introduced as empty to the axis boundary because it is an edge of the domain, not a surface. Additionally, wedge has been assigned to front and back boundaries, simultaneously, because they have been defined as two repetitive sides of the sector.

In this study, at inletfuel, by assuming ${\mathrm{Re}}_{D,fuel}=68000$, $C_{\mu }=0.09$and $D_{h,fuel}=5.26mm$for fully turbulent Propane jet, the turbulent intensity (turbulence level) , $I$, the boundary condition for turbulence kinetic energy, k, and the boundary condition for its dissipation rate, ε, have been successively calculated as [25]:

$$I_{inlet,fuel}=0.16{\mathrm{Re}}_{D,fuel}^{\raisebox{1ex}{$-1$}\!\left/ \!\raisebox{-1ex}{$8$}\right.}\simeq 0.0398$$

$$k_{inlet,fuel}=\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.{(U.I)}^2\simeq 6.679$$

$$l_{inlet,fuel}=0.07D_{h,fuel}=0.000368m$$

$${\epsilon }_{inlet}=C_{\mu }^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$4$}\right.}{k}_{inlet,fuel}^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{l}_{inlet,fuel}^{-1}\simeq 5630$$

Moreover, at inletair, by assuming $I_{\mathit{inlet},air}=0.004$, $C_{\mu }=0.09$and $D_{h,air}=30cm$for turbulent co-flowing air, we have these boundary conditions:

$$k_{inlet,air}=\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.{(U.I)}^2\simeq 0.002$$

$$l_{inlet,air}=0.07D_{h,air}=0.021m$$

$${\epsilon }_{inlet,air}=C_{\mu }^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$4$}\right.}{k}_{inlet,air}^{\raisebox{1ex}{$3$}\!\left/ \!\raisebox{-1ex}{$2$}\right.}{l}_{inlet,air}^{-1}\simeq 0.0006998$$

Although it is clear that the steady-state results depend on the boundary conditions, not the initial values, it is tried to assign the logical initial values for various variables.

5. Other governing equations [25]

5.1. Mass Conservation Equation

$$\frac{\partial \rho }{\partial t}+\nabla .(\rho \overrightarrow{U})=0$$

5.2. Energy Conservation Equation

$$\frac{\partial }{\partial t}(\rho E)+\nabla .(\overrightarrow{U}(\rho E+p))=-\nabla {\displaystyle {\sum }_j{h}_j}{J}_j+S_h$$

5.3. Momentum Conservation Equations

$$\frac{\partial (\rho \overrightarrow{U})}{\partial t}+\nabla .(\rho \overrightarrow{U}\overrightarrow{U})=-\nabla p+\rho \overrightarrow{g}+\overrightarrow{F}$$

(1)

5.4. Transport equations for turbulence kinetic energy, k, and its dissipation rate, ε, [26,27,28]

6. Numerical discretizing schemes

Table 4 comprises a detailed list of assigned discretizing schemes for diverse terms:

7. Solution procedures

Due to lack of combustion and reaction in the experiment, the files of the constant folder have been modified to eliminate these phenome in the numerical test cases. The set of governing equations have been resolved through PIMPLE algorithm by setting nOuterCorrectors, nCorrectors, nNonOrthogonalCorrections, and relaxationFactors to 1, 2, 0 and 1, respectively.

In an entirely transient simulation of such problem, at first, slower fluid's requisite time of passing the entire of the domain should be obtained. Typically, 15 to 20 times of that is sufficient as the total required time of simulating. Moreover, the corresponding time step can be elicited by limiting the Courant number to 3 or 4 (it depends on the nature of the problem and trial and error). But in this study, because just the steady-state needs to be regarded, to accelerate a solution to such condition, local-time stepping method has been used (i.e., pseudo-transient simulation), by adjusting the discretizing scheme of temporal terms to localEuler. Therefore, both endTime and deltaT which are introduced to controlDict file have no physical meaning and represent the number of iterations.

Moreover, in the fvSolution file, tolerance (the convergence level) has been adjusted to 10e-08 and 10e-06 for Yi (species) and all of the other variables, respectively.

In the sampleDict file of the system folder, required data have been introduced on x/D=0, 4, 15, 30 and 50, and also adjacent to axis boundary (y/D=0), each on 101 points. Subsequently, although these data have been plotted in the createGraphs folder of each test case by the gnuplot, the corresponding data have been obtained from the results which have been recorded by OpenFOAM, for the last iteration, in postProcessing>sample>30000 folder of each test cases. Eventually, the mentioned data have been transferred to Excel and plotted alongside the experimental data to compare with each other, as have been shown in Figures 6 to 28 at the end of this study.

Initially, as it has been demonstrated in Table 5, four meshes with different cells and aspect ratios have been created:

Afterward, all the tests cases have been performed for 30,000 iterations. Subsequently, the massGraph file has been plotted as massGraph.eps for each test case to investigate needed iterations for converging (Figures 2, 3, 4 and 5). Finally, the steady-state (converged) results of all test cases have been compared to each other to pick out the independent grid.

It should be noticed that all the phases of the performance for each test case have been carried out through Allrun script file, in serial mode (on single CPU), on an ordinary laptop.

To assess the accuracy of the numerical measurements, it has been observed that for each test case and on 4 reference levels (x/D=0, 4 , 15 and 30), massGraph.eps (convergence diagram of mass flux) and momentumGraph.eps (convergence diagram of total momentum) plots have revealed the convergence and conversation of Propane and conservation of total momentum.

It can be perceived that after passing sufficient requisite iterations in each test case, total mass flux of numerical domain (5 degrees sector) at all 4 reference levels, has been converged to 2.952e-02 g/s (=2.125 gr/s for total reference assumed cylinder) which is 7.6% less than measured mass flux through flow meter (2.8% less than measured mass flux by data extracted through Rayleigh Scattering System and Laser Doppler Velocimetry).

8. Results and discussions

The independent grid's converged results have been presented through relevant figures at the end of this section to abbreviate the discussion.

Table 6 comprises a detailed list of assigned discretizing schemes for various terms, which have been placed as the foundation of the discussion:

8.1. Umean

As it has been apparently demonstrated in Figure 18, due to imposing a mean uniform profile for axial velocity to the numerical model at inletfuel, the jet has been authorized to establish a fully turbulent profile while passing through the nozzle.

Consequently, at the end of the 10 cm nozzle (i.e., at the inject point) fully-developed pipe flow theory is instituted and Umean's profile has the maximum value of 64.028 m/s (7.2% less than the maximum of experimental Umean).

The numerical Umean remains nearly constant for approximately 6 jet diameters downstream of the nozzle outlet, before decreasing rapidly to approach the outer co-flowing air velocity of 9.2 m/s farther downstream.

It should be emphasized that the propellant of the mixing between injected Propane and the co-flowing air is nothing but turbulence which is due to vorticities.

The transport equation of the aforementioned vorticities has a source term which involves the gradient of the mean velocities.

Accordingly, these vorticities (and consequently, the turbulence and the mixing from there), will be intensified where the gradients of velocities are more augmented.

It can be perceivable from Figures 18, 20, 21, 22 and 23 that at first, the radial profile of numerical Umean has a superb adaptation with corresponding experimental data in the vicinity of the jet (i.e., in x/D<10).

But subsequently, having marched toward the outlet, this conformity will be gradually diminished (i.e., in 10<x/D<30).

Eventually, this coincidence will be revitalized and be radiantly improved (i.e., in x/D>30).

It should be kept in mind that each radial distance, y, is normalized by the inner diameter of the nozzle, D.

In numerical mean radial profiles of Umean (Figures 20, 21, 22 and 23), it is observable that at increasing radial distances, although the mean velocities rapidly approach free stream values of 9.2 m/s, the rate of this acceding will be pacified.

It should be accentuated that this behavior of velocity is similar than demeanor of Mixed Fraction of Propane in L_f's diagram (Figure 6), i.e., having receded from the nozzle, the width of the jet will be gently expanded.

Finally, in radial profile of Umean at x/D=0, (Figure 19), it can be noticed that the numerical results are in an eligible conform with corresponding experimental data, especially from the outer diameter of the nozzle (i.e., y/D=1.7) toward the leftside boundary.

8.2. Vmean

It can be perceivable from Figures 25, 26, 27 and 28, the numerical Vmean has a conceptual adaptation with both AIR seeded experimental results, and the corresponding JET seeded data.

However, having approached toward the outlet, the scale of this conformity will be leniently assuaged.

Unfortunately, in the axial profile of Vmean at y/D=0 (Figure 24), it is not conceivable to attain a logical coincidence between numerical results and experimental ones.

It may be caused by debilities of RANS models in such problem, or by excessive simplification of the problem by assumptions as mentioned earlier in the present numerical model.

It should be emphasized that capturing radial velocity (Vmean) seems to be frequently so demanding, even with more sophisticated models, for example, LES [29].

8.3. Turbulence Kinetic Energy (k)

It can be perceivable from Figures 12, 14, 15, 16 and 17 that the radial profiles of turbulence kinetic energy (k) have an appropriate overall adaptation with corresponding experimental data.

Having approached toward the outlet, this conformity will be gradually improved, due to alleviated fluctuations of velocity. (approximately after x/D<50).

In numerical mean radial profiles of k (Figures 14, 15, 16 and 17), it is observable that at increasing radial distances, the mean turbulence kinetic energy, k, will be gradually alleviated.

Again, as like as Umean, having receded from the nozzle, the width of the jet will be steadily expanded, which is conform with the manner of Mixed Fraction of Propane in L_f's diagram(Figure 6).

Eventually, in radial profile of turbulence kinetic energy, k, at x/D=0, (Figure 13), as like as Umean it can be noticed that the numerical results are in a good adaptation with corresponding experimental data, especially from the outer diameter of the nozzle (i.e., y/D=1.7) toward the leftside boundary.

8.4. Mixture Fraction of Propane

Due to technical problems of corresponding experiments, it is not conceivable to elicit experimental data of Propane's mean mixture fraction at the inlet (i.e., x/d=0).

Thus, corresponding axial and all the radial comparisons have been carried out at x/D=4, 15, 30, 50 (Figures 6, 7, 8, 9, 10 and 11).

As it is manifestly observable in these figures, a magnificent conformation between numerical results and experimental data have been revealed in the sphere of mean mixture fraction and its half radius (L_f), either for axial profile or radial ones.

It should be noticed that in numerical results at x/D=4, the mean mixture fraction values of unity (corresponding to pure Propane) announces the existence of a potential core, which extends approximately 0.2 diameters from the centerline.

The numerical mean mixture fraction remains approximately constant over this potential core region, which extends nearly 6 jet diameters downstream of the nozzle's outlet, before drastically decreasing due to co-flowing air is entrained because of Kelvin-Helmholtz instabilities (due to drastic difference between Jet's and air's velocities).

8.5. Mixture Fraction Half Radius (L_f)

L_f, ascertains the jet spreading rate, would be defined as the radial location at which the mixture fraction is equal to half its value at the centerline. Figure 11, shows the variations of mixture fraction half radius, L_f, (normalized by D) with axial distance. It is explicit that it reveals a phenomenal conformity with the corresponding experimental data.

9. Conclusion

Having manipulated reactingFoam (one of the standard solvers in OpenFOAM) through eliminating reaction and combustion, some numerical simulation of velocity and mixture fraction fields have been carried out in a turbulent non-reaction Propane jet flow issuing into parallel co-flowing air, in isothermal condition, and compared with experimental data which had been previously described here: http://www.sandia.gov/TNF/DataArch/ProJet.html.

In this regard, a two-equation Realizable k-ε eddy viscosity turbulence model has been used to simulate the turbulent flow field on a 2D plane (i.e., on a 5-degree sector of the experimental domain) to attain a better understanding of the flow structure and mixing process in a situation which there is no chemical interaction and heat-transfer.

Having compared the elicited numerical axial and radial profiles of mean velocities, turbulence energy, mean mixture fraction, and the mixture fraction half radius (L_f) with corresponding experimental data, superb conformations have been revealed for mixture fraction fields and mixture fraction half radius (L_f). Moreover, in contrast with Vmean, it has been relatively observed appropriate adaptations for Umean, especially when the investigated areas had been inclined toward the nozzle exit.

Altogether, regarding the extreme shortened numerical expenditures (due to the significant simplifying assumption in both test cases' diverse adjustments and turbulence model) in comparison with more sophisticated turbulence models (for example LRR, SSG, DES, LES, and so on), the elicited numerical results have eminently demonstrated relevant precisions and conformity with the experimental data in aforesaid spheres.

Figure 2. Convergence Diagram of Mass Flux for Mesh A.

Figure 2. Convergence Diagram of Mass Flux for Mesh A.

Figure 3. Convergence Diagram of Mass Flux for Mesh B.

Figure 3. Convergence Diagram of Mass Flux for Mesh B.

Figure 4. Convergence Diagram of Mass Flux for Mesh C.

Figure 4. Convergence Diagram of Mass Flux for Mesh C.

Figure 5. Convergence Diagram of Mass Flux for Mesh D.

Figure 5. Convergence Diagram of Mass Flux for Mesh D.

Figure 6.

Figure 6.

Figure 7.

Figure 7.

Figure 8.

Figure 8.

Figure 9.

Figure 9.

Figure 10.

Figure 10.

Figure 11.

Figure 11.

Figure 12.

Figure 12.

Figure 13.

Figure 13.

Figure 14.

Figure 14.

Figure 15.

Figure 15.

Figure 16.

Figure 16.

Figure 17.

Figure 17.

Figure 18.

Figure 18.

Figure 19.

Figure 19.

Figure 20.

Figure 20.

Figure 21.

Figure 21.

Figure 22.

Figure 22.

Figure 23.

Figure 23.

Figure 24.

Figure 24.

Figure 25.

Figure 25.

Figure 26.

Figure 26.

Figure 27.

Figure 27.

Figure 28.

Figure 28.