Numerical simulation of turbulent compressible flows in a C-D nozzle with different divergence angles

1.

Introduction

Flow of gases through a converging-diverging nozzle is one of the important problems used for modeling compressible flow using computational fluid dynamics (CFD).The converging – diverging nozzles used to accelerate the fluid to supersonic speeds passed the throat of such a nozzle [1]. The efficiency of an engine nozzle, usually represented by the nozzle thrust coefficient (C_f), is defined as the ratio of the actual nozzle gross thrust to the ideal gross thrust.

Nowadays, numerical simulation of the flow field in C-D nozzles is one of the main steps in design process of engines. Geatz [2] describes the creation of a computational code that can be used to predict the thrust performance characteristics of non-axisymmetric two-dimensional convergent-divergent exhaust nozzles. The nozzle internal performance prediction code shows excellent agreement with experimental data in predicting the peak gross thrust coefficients for basic 2D-CD nozzle geometries.

Turbulence modeling plays a major role in simulating such a complicated flow. The previous numerical studies have assessed the accuracy of the turbulence model for predicting the flow filed and the thrust performance accurately. Q. Xiao et al. [3] investigated the compressible jet plume from a planar over-expanded nozzle by solving the Reynolds-Averaged Navier-Stokes (RANS) equations with several turbulence models. By comparing the simulation results with available experimental data, they show that, two- equation Shear Stress Model (SST) gives the best results and the simulations are able to predict the velocity profiles, total pressure decay and axial jet thickness distribution in the jet plume reasonably well.

Hamed and Vogiatzis[4] used the two-equation k-ω turbulence model in conservation law form and general curvilinear coordinate to predict the surface pressure distribution and internal thrust coefficient of a 2D-CD nozzle. Hamed and Vogiatzis[5] by using the Spalart-Allmaras turbulence model confirmed the importance of the turbulence model in producing realistically or unrealistically numerical results.

Balabelet al. [6], assessed several turbulence models in terms of their effects on the agreement between the experimental centerline pressure distribution and the 2D computational results at over-expanded conditions. Their results indicate that both the prediction of the shock location and pressure level behind the shock strongly depend on the applied turbulence model.

The accuracy of turbulent jet plume prediction significantly depends on the numerical method and turbulence model. Several authors have investigated the effect of different turbulence models on the prediction of jet mixing. Georgiadis and Papamoschou [7] studied single and coaxial dual-stream jets using RANS with linear two-equation and nonlinear two-equation explicit algebraic stress turbulence modelling. Their comparison of computed mean flow field development with experiments shows that the standard SST model provides the overall best agreement with experimental data.

Chenault and Beran [8] conducted a numerical investigation of supersonic injection using second-order Reynolds-Stress turbulence model proposed by Zhang et al. [9] as well as the k-ε model. Detailed comparison with experimental data show that the Reynolds-stress model simulation results are physically consistent. However, the simulations with the k-ε model resulted in nonphysical and inconsistent turbulence prediction.

Dembowski and Georgiadis [10] conducted a numerical study for supersonic axisymmetric jet flow using two-equation SST and k-ε model with and without compressibility correction. Their results indicate that these models do not predict supersonic nozzle flows accurately.

Nozzle performance was considered in many experimental and numerical studies, especially from the point of flow and heat transfer characteristics with various inlet-boundary conditions and flow geometries. Park et al. [11] investigated sonic nozzles that were applied to gas flow rate measurements and determined that the critical pressure ratio was highly dependent on the Reynolds number rather than area ratio especially in the cases with low flow velocity.

Paik et al. [12], who determined higher discharge coefficients with the increasing of mass flow rate, reported a variation of discharge coefficients for sonic nozzles with flow geometry and Reynolds number. In addition, Ahmad [13] correlated the variation of nozzle discharge coefficient for choked nozzle flows.

Kim et al.[14] considered several kinds of gases and turbulence models with a wide range of Reynolds numbers on different sonic nozzle geometries. Recently, Balabel et al. [15] applied the turbulent gas flow dynamics in a two-dimensional convergent– divergent rocket nozzle and numerically predicted the associated physical phenomena for various operating conditions.

Although there is a large amount of studies concerning the flow through the CD nozzles over a wide range of nozzle pressure ratio (NPRs), the detailed investigation on the effects of the divergence half-angle on thrust coefficient and discharge coefficient in CD nozzle is less studied.

In the current study the supersonic flow inside a 2Dconverging-diverging nozzle and its exhaust plumes simulated numerically with three different turbulence models and the results were compared with the numerical results of Krishnamurty and Shy [16] and the experimental data of Mason et al. [17].According to the applied analysis, the proper turbulence model, which gives the best prediction of the flow filed, is selected. Then the parametric study is carried out in which the divergence half-angle and nozzle pressure ratio (NPR) are varied. Attention is also paid to the effect of nozzle angle on nozzle performance and changes in the value of thrust coefficient ( ) and discharge coefficient ( ).

2. Governing Equation

The FLUENT software is used to solve the Reynolds-average Navier-Stokes (RANS) equation with turbulence models. The governing equations, which include the conservation equations for mass, momentum and energy, along with the equation of state, are treated in generalized coordinates and in conservative form.

In Reynolds averaging, the solution variables in the instantaneous (exact) Navier-Stokes equations are decomposed into the mean (ensemble-averaged or time-averaged) and fluctuating components. For the velocity components:

Where and  are the mean and fluctuating velocity components (i = 1; 2; 3).

Likewise, for pressure and other scalar quantities:

Where  denotes a scalar such as pressure, temperature, or species concentration.

Substituting expressions of this form for the flow variables into the instantaneous continuity and momentum equations and taking a time (or ensemble) average (and dropping the over-bar on the mean velocity, ) yields the ensemble-averaged momentum equations.

They can be written in Cartesian tensor form as:

Equations (3) and (4) are called Reynolds-averaged Navier-Stokes (RANS) equations. Additional terms now appear that represent the effects of turbulence.

These Reynolds stresses, , must be modeled in order to close Equation (4) [18].Several available turbulence models are employed to predict the flow behavior in the considered physical domain.

Three turbulence models are evaluated for the present case, namely, the Reynolds stress model (RSM) [19], the standard k-ε (STD) model [20] and the shear-stress transport (SST) k-ω model [21].

The Reynolds stress model (RSM) [19] is the most elaborate turbulence model that FLUENT provides. Abandoning the isotropic eddy-viscosity hypothesis, the RSM closes the Reynolds-averaged Navier-Stokes equations by solving transport equations for the Reynolds stresses, together with an equation for the dissipation rate. This means that five additional transport equations are required in 2D flows and seven additional transport equations must be solved in 3D.This modeling approachwasdeveloped by launder et al. [19].

The standard k-ɛmodel [20] is a semi-empirical model based on model transport equations for the turbulent kinetic energy (k) and its dissipation rate (ɛ). The model transport equation for k is derived from the exact equation, while the model transport equation for ɛ is obtained using physical reasoning and bears little resemblance to its mathematically exact counterpart.

The shear-stress transport (SST) k-ω model was developed by Menter [21] to effectively blend the robust and accurate formulation of the k-ω model in the near-wall region with the free-stream independence of the k-ε model in the far field. In this model, the definition of the turbulent viscosity is modified to account for the transport of the turbulent shear stress. Details of the closure equations for mentioned turbulence models can be found in Ref. [18].

3. Grid independence Study

Mason et al. [17] conducted an experiment to determine the effect of throat contouring on the nozzle internal performance. They tested five non-axisymmetric converging-diverging nozzles in the static test facility of the Langley 16-foot transonic tunnel and recorded internal performance data at nozzle pressure ratios up to 9.0.

In the current study, the geometry of CD nozzle, which Mason et al. [17] had used for their experiment, is used as baseline nozzle geometry. The schematic of nozzle geometry is shown in Figure 1.The design parameters for configuration of baseline nozzle geometry

Table1. The design parameters of baseline nozzle geometry (all dimensions are in cm)

are provided in Table 1 as well.

The baseline geometry has the constant divergence half-angle. In the current research, the parametric study of the divergence half- angle effect ranging from 5º to 20º is presented

In this parametric study, the nozzle pressure ratio (NPR) has been changed in the range of 6-10. The expansion ratio (ratio of exit area to throat area) of the baseline CD nozzle is 1.8, so by using the gas dynamics relations, the value of the design NPR will be 8.81.

Three different grids are generated to simulate the gas flow inside the baseline geometry of CD nozzle and its exhaust plume. A 2D computational domain with the assigned boundary conditions is shown in Figure 2.

The results are presented for a course mesh, with an interval size of 0.96 mm, a medium mesh, with an interval size of 0.6 mm and a fine mesh, with an interval size of 0.4 mm inside the convergent–divergent nozzle which correspond to a total computational grid (inside and outside the nozzle) of 7 000 cells (Grid A), 23 000 cells (Grid B) and 56 000 cells (Grid C), respectively. Figure 3 shows the details of two-dimensional mesh near the nozzle (Grid B) of the computational domain.

The comparison (Figure 4) shows that the results

obtained using Grid B and Grid C, are very close. Therefore, the computations are performed on Grid B in a compromise between computational time and accuracy.

Figure 4 shows the predicted Mach number on the centerline of computational domain (symmetry line) using three different meshes at design nozzle pressure ratio (NPR=8.81).

4. Boundary Conditions

Boundary conditions are specified for the inlet, outlets, and solid nozzle walls (Figure 2). At the nozzle inlet, total pressure and total temperature are specified. At the outlet, the boundary condition of pressure far-filed with very small Mach number ( ) is performed and the ambient pressure is set to the atmospheric pressure value.

The simulation is two-dimensional and all of the walls are considered adiabatic. At walls, no-slip boundary conditions are imposed and symmetry boundary condition is applied at the nozzle centerline.

5. Validation of Numerical Model

Numerical simulations of the supersonic gas flow inside the converging-diverging nozzle have been done using FLUENT 6.3.26. In this study, three different turbulence models (STD k-ɛ, SST k-ω and RSM) are used to simulate the flow field inside the CD nozzle and its exhaust plume.

Airflows with viscosity of and with the range of nozzle pressure ratio 6 to 10 are simulated. The ambient temperature of 300 K and the ambient pressure of 101325 Pa are considered. Turbulence intensity of 10% and hydraulic diameter of are assumed for calculation of turbulence quantities at the inlet. The compressible, steady-state, Reynolds-averaged, Navier-Stokes equations are solved using a density-based algorithm with a second order accuracy in space. The residual of 10^-5 is selected as convergence criterion for termination of the solution procedure of the equations.

Figure 5 compares the predictions of the three turbulence models for the static pressure on the upper nozzle sidewall with the experimental measurements [17] and numerical results [16] at design nozzle pressure ratio (NPR=8.81). The NPR is defined as the ratio of the total (or stagnation) pressure at the nozzle inlet (P_o) to the ambient pressure (P_a).

The static wall pressure is normalized using the total pressure P_o and is plotted against the dimensionless axial location x/L, where L is the length of the nozzle.

The results confirmed that, the SST k-ω turbulence model showed the better predictions than STD k-ɛ model and RSM. Therefore, this model is used to investigate the effect of nozzle half-angles on the nozzle performance.

For better comparison of different turbulence models performance, the numerical results of Mach number of gas flow on outlet section of C-D nozzle are plotted in Fig.6.

6. Result and Discussion

The nozzle thrust coefficient, ,is defined as the ratio of the actual nozzle thrust to the computed ideal nozzle thrust ( ). In addition, the discharge coefficient, , is defined as the ratio of the measured mass-flow rate to ideal mass-flow rate ( ).For the 2D nozzle, they are defined by

Where, ρ, u, and p are local density, axial velocity, and pressure, respectively. The subscript, “ideal”, refers to

quantities calculated for an ideal nozzle based on one-dimensional isentropic flow assumption.

The isentropic thrust ( ), and mass flow rate, ( ), are calculated by the following equations [2, 17].

Where is the specific heat ratio ( ), R is the gas constant ( ), and are stagnation conditions.

Figure 7 provides a comparative view of the thrust coefficient prediction to the computational results of Geatz [2] and the experimental data of Mason et al. [17] for baseline nozzle configuration.  One can see that the nozzle performance prediction provides the results that are well in agreement with the experimental and numerical data over an NPR range of 6-10. It is seen that the present numerical results over-predict the experimental results by up to 1%.

The computed Mach number for different Nozzle Pressure Ratios is plotted in Fig.8. This figure describes briefly the influence of nozzle pressure ratio on gas flow behavior inside a symmetric CD nozzle and the exhaust plume. As shown, by increasing the NPR, the amplitude of oscillations will decrease.

Various divergence half-angles of CD nozzle are considered in simulations. The contour of Mach number with various diverging half-angles at nozzle pressure ratio of 10 is shown in Fig.9. By changing the nozzle half-angles, the other nozzle parameters are held constant.  The CD nozzle with (ɛ=10.85º) is the baseline geometry which is studied by Mason et al [17] and

Krishnamurty et al [16] (see Fig.9b). As shown, the divergence nozzle half- angle has a great effect on the nature of the oblique shock waves.

Thermal field behavior variation due to divergence angle and shock wave locations shown by the contour of static temperature for various nozzle geometries in Fig.10.

In Figure 11, the thrust coefficient for the nozzles with various divergence angles are plotted against the nozzle pressure ratio. For each configuration, the value of increases from a minimum at the lowest nozzle pressure ratios to a peak level near the design nozzle pressure ratio of 8.81.

This figure shows that, for a given divergence nozzle angle, the thrust coefficient increases with increasing nozzle pressure ratio, and, for a given nozzle pressure ratio, the thrust coefficient increases as the nozzle half- angle decreases.

The discharge coefficient data in Figure 12 shows some variation with nozzle geometry. However, as

should be expected, is almost independent of nozzle pressure ratio since the nozzles are choked for all cases in the numerical simulations. As can be seen for a given nozzle pressure ratio, the discharge coefficient increases as the nozzle angle increases.

Figure 13 shows that by increasing the divergence angle, the first oblique shock wave, moves closer to the outlet section of nozzle. This phenomenon is more evident for nozzle half angle of 20°.  For this reason, the value of thrust coefficient for divergence nozzle angle of 20° has a difference (about 2%) from other nozzles.

Comparison of internal static temperature distributions along the upper nozzle sidewall is presented in Figure 14 at nozzle pressure ratio of 10. The data are presented as local static temperature normalized by inlet total temperature, , and are plotted as a function of x normalized by length of the nozzle. This figure shows

that by increasing the value of divergence angle, the shock waves will be weaker and also the gas flow behavior will be closer to the ideal conditions.

7. Conclusions

In this paper, the effects of divergence half-angles on the nozzle performance at different nozzle pressure ratios (NPR) in the range of 6-10 are investigated numerically. Three different two-equation turbulence models are used for simulating the compressible gas flow inside the two-dimensional converging-diverging nozzle and its exhaust plume. For validation of the numerical method, the predicted results of thrust coefficient are compared with available experimental data of Mason et al. [17] and numerical results of Geatz [2]. The numerical results revealed that, the SST k–ω model gave the best results compared with other models.

The effect of changes in value of divergence half-angle ( ) on the nozzle performance, thrust coefficient ( ) and discharge coefficient ( ) has been investigated numerically.

Some important findings from the gas flow simulations in CD nozzle are as the following:

• The predicted results show that for a given divergence nozzle angle, the thrust coefficient increases with increasing nozzle pressure ratio, and, for a given nozzle pressure ratio, the thrust coefficient increases as the nozzle angle decreases.
• For each configuration, the value of increases from a minimum at the lowest nozzle pressure ratio to a peak level near the design nozzle pressure ratio of 8.81.
• When the CD nozzle is chocking, the value of discharge coefficient is independent of nozzle pressure ratio and also for a given nozzle pressure ratio, the discharge coefficient increases as the divergence nozzle angle ( ) increases.
• Static temperature distribution along the upper nozzle sidewall shows that by increasing the value of divergence angle, the shock waves will be weaker and the gas flow behavior will be closer to the ideal conditions.

References

[1].    K. K. Quintao, Design Optimization of Nozzle Shapes for Maximum Uniformity of Exit Flow,Florida International University,M.S. Thesis, (2012).

[2].    A. M. Geatz, A Prediction Code for the Thrust Performance of Two Dimensional, Non-Axisymmetric, Converging Diverging Nozzles,Air Force Institute of Technology,PH.D Thesis, (2005).

[3].    Q. Xiao, H.M. Tsai , D. Papamoschou,  Numerical Study of Jet Plume Instability from an Over-expanded Nozzle, 45th AIAA Aerospace Sciences Meeting and Exhibit, Nevada, AIAA, (2007).

[4].    A.Hamed,C.Vogiatzist, Over-expanded Two-Dimensional Convergent-Divergent Nozzle Performance, Effects of Three-Dimensional Flow Interactions, Journal ofPropulsion and Power,14, 234-240, (1998).

[5].    A.Hamed, C.Vogiatzist, Three-Dimensional Flow Computations and Thrust Predictions in 2DCD Over-expanded Nozzles, AIAA paper 97-0030, (1997).

[6].    A. Balabel, A.M. Hegab, S. Wilson, M. Nasr, S. El-Behery, Numerical Simulation of Turbulent Gas Flow in a Solid Rocket Motor Nozzle, 13th International Conference on Aerospace Science and Aviation Technology, Egypt, (2009).

[7].    N. J. Georgiadis, D. Papamoschou, Computational Investigation of High-Speed Dual Stream Jets, AIAA 2003-3311, (2003).

[8].    C. F.Chenault, P.Beran, Numerical Investigation of Supersonic Injection Using a Reynolds-Stress Turbulence Model, AIAA Journal, 37, 1257-1269, (1999).

[9].    H.Zhang, R.So, T.Gatski, C. Speziale, a Near-Wall Second-Order Closure for Compressible Turbulent Flows, Near-Wall Turbulent Flows, edited by R. So, C. Speziale, and B. Launder, Elsevier, New York, 209-218,(1993).

[10]. M.A. Dembowski, N.J. Georgiadis, an Evaluation of Parameters Influencing Jet Mixing Using the WIND Navier-Stokes Codes, NASA/TM-211727,(2002).

[11]. K. A. Park, Y. M. Choi, H. M. Choi, T. S. Cha, B.H. Yoon,The Evaluation of Critical Pressure Ratios of Sonic Nozzles at Low Reynolds Numbers. Flow Measure, Instrum,12, 37–41, (2001).

[12]. J. S. Paik, K. A. Park, J. T. Park, Inter-Laboratory Comparison of Sonic Nozzles atKRISS,Flow Measure, Instrum, 11, 339–344, (2000).

[13]. R. A. Ahmad, Discharge Coefficients and Heat Transfer for Axisymmetric Supersonic Nozzles. Heat Transfer Eng, 22, 40–61, (2001).

[14]. H. D. Kim, J. H. Kim, K. A. Park, T.Setoguchi, S. Matsuo, Computational Study of the Gas Flow Through a Critical Nozzle. Proc. Instn. Mech. Engrs, 217, 1179–1189, (2003).

[15]. A. Balabel, A.M. Hegab, M. Nasr, Samy M. El-Behery, Assessment of Turbulence Modeling for Gas Flow in Two-Dimensional Convergent–Divergent Rocket Nozzle, Applied Mathematical Modelling, 35, 3408–3422, (2011).

[16]. V.S,Krishnamurty, W.Shyy, Effect of Wall Roughness on the Flow through Converging-Diverging Nozzles, Journal of Propulsion and Power, 13, 753-762, (1997).

[17]. M.L. Mason, L.E. Putnam,J.R. Richard, The Effect of Throat Contouring on Two Dimensional Converging-Diverging Nozzles at Static Condition, NASA Technical Paper 1704, (1980).

[18]. Fluent, User’s Guide Fluent 6.3.26, Fluent Incorporated, Lebanon, NH, (2006).

[19]. B. E. Launder, G. J. Reece, and W. Rodi, Progress in the Development of aReynolds-Stress Turbulence Closure, J. Fluid Mech., 68, 537-566, (1975).

[20]. B. E. Launder, D. B. Spalding, Mathematical Models of Turbulence, AcademicLondon, 169-189, (1972).

[21]. F. R.Menter, Two-Equation Eddy-Viscosity Turbulence Models for Engineering Applications, AIAA journal, 32, 1598-1605, (1994).