Document Type : Full Lenght Research Article
Authors
^{1} MalekAshtar University of Technology, Tehran, Iran.
^{2} MalekAshtar University of Technology, Tehran, Iran
Abstract
Keywords
1.
History: Received 12 December 2013 Received in revised form 23 February 2014 Accepted 9 June 2014
Keywords: Numerical simulation convergingdiverging nozzle Divergence angle Thrust coefficient Discharge coefficient Static temperature 
A B S T R A C T Compressible gas flow inside a convergentdivergent nozzle and its exhaust plume at different nozzle pressure ratios (NPR) have been numerically studied with several turbulence models. The numerical results reveal that, the SST k–ω model give the best results compared with other models in time and accuracy. The effect of changes in value of divergence halfangle ( ) on the nozzle performance, thrust coefficient ( ) and discharge coefficient ( ) has been investigated numerically. The predicted results show that for a given divergence angle, the thrust coefficient ( ) increases by increasing nozzle pressure ratio. Also, for a given nozzle pressure ratio, the thrust coefficient increases as the nozzle divergence angle decreases. 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.
© 2014 Published by Semnan University Press. All rights reserved. 
Numerical simulation of turbulent compressible flows in a CD nozzle with different divergence angles
Mohammad Hadi Hamedi Estakhrsar, Mehdi Jahromi* ^{ } MalekAshtar University of Technology, Tehran, Iran



Journal of Heat and Mass Transfer Research
Journal homepage: http://jhmtr.journals.semnan.ac.ir 


Introduction
Flow of gases through a convergingdiverging 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 CD 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 nonaxisymmetric twodimensional convergentdivergent exhaust nozzles. The nozzle internal performance prediction code shows excellent agreement with experimental data in predicting the peak gross thrust coefficients for basic 2DCD 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 overexpanded nozzle by solving the ReynoldsAveraged NavierStokes (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 twoequation kω turbulence model in conservation law form and general curvilinear coordinate to predict the surface pressure distribution and internal thrust coefficient of a 2DCD nozzle. Hamed and Vogiatzis[5] by using the SpalartAllmaras 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 overexpanded 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 dualstream jets using RANS with linear twoequation and nonlinear twoequation 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 secondorder ReynoldsStress turbulence model proposed by Zhang et al. [9] as well as the kε model. Detailed comparison with experimental data show that the Reynoldsstress 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 twoequation 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 inletboundary 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 twodimensional 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 halfangle on thrust coefficient and discharge coefficient in CD nozzle is less studied.
In the current study the supersonic flow inside a 2Dconvergingdiverging 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 halfangle 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 Reynoldsaverage NavierStokes (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) NavierStokes equations are decomposed into the mean (ensembleaveraged or timeaveraged) and fluctuating components. For the velocity components:
(1) 
Where and are the mean and fluctuating velocity components (i = 1; 2; 3).
Likewise, for pressure and other scalar quantities:
(2) 
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 overbar on the mean velocity, ) yields the ensembleaveraged momentum equations.
They can be written in Cartesian tensor form as:
(3) 

(4) 
Equations (3) and (4) are called Reynoldsaveraged NavierStokes (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 shearstress transport (SST) kω model [21].
The Reynolds stress model (RSM) [19] is the most elaborate turbulence model that FLUENT provides. Abandoning the isotropic eddyviscosity hypothesis, the RSM closes the Reynoldsaveraged NavierStokes 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 semiempirical 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 shearstress transport (SST) kω model was developed by Menter [21] to effectively blend the robust and accurate formulation of the kω model in the nearwall region with the freestream 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 nonaxisymmetric convergingdiverging nozzles in the static test facility of the Langley 16foot 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
Fig.1 Schematic of nozzle geometry [17] 
Table1. The design parameters of baseline nozzle geometry (all dimensions are in cm)
Parameter 

1.8 
3.52 
1.37 
2.46 
0.68 
10.8 
20.8 
Value 
are provided in Table 1 as well.
The baseline geometry has the constant divergence halfangle. 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 610. 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 twodimensional mesh near the nozzle (Grid B) of the computational domain.
The comparison (Figure 4) shows that the results
Fig.2 2D computational domain with the assigned boundary conditions 
Fig.3 Mesh generated for flow computation 
Fig.4 Mach number distribution on the centerline of computational domain 
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 farfiled with very small Mach number ( ) is performed and the ambient pressure is set to the atmospheric pressure value.
The simulation is twodimensional and all of the walls are considered adiabatic. At walls, noslip 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 convergingdiverging 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, steadystate, Reynoldsaveraged, NavierStokes equations are solved using a densitybased 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 halfangles 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 CD 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 massflow rate to ideal massflow rate ( ).For the 2D nozzle, they are defined by
(10) 

And 

(11) 
Where, ρ, u, and p are local density, axial velocity, and pressure, respectively. The subscript, “ideal”, refers to
Fig.5 Internal staticpressure distributions along the nozzle sidewall for three different turbulence models 
quantities calculated for an ideal nozzle based on onedimensional isentropic flow assumption.
The isentropic thrust ( ), and mass flow rate, ( ), are calculated by the following equations [2, 17].
(12) 

And 

(13) 
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 610. It is seen that the present numerical results overpredict 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 halfangles of CD nozzle are considered in simulations. The contour of Mach number with various diverging halfangles at nozzle pressure ratio of 10 is shown in Fig.9. By changing the nozzle halfangles, 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
Fig.7 Comparison of the thrust coefficient prediction 
Fig.8 The variations of flow Mach number with nozzle pressure ratio 
Fig.9 The contour of Mach number with various divergent halfangles at nozzle pressure ratio of 10,(a) ɛ=5º, (b) ɛ=10.85º, (c) ɛ=15º and (d) ɛ=20º 
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
Fig.10 The contour of static temperature with various divergent halfangles at nozzle pressure ratio of 10, (a) ɛ=5º, (b) ɛ=10.85º, (c) ɛ=15º and (d) ɛ=20º 
Fig.11 The thrust coefficient ( ) for the nozzles with various divergence angles (ɛ) 
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
Fig.12 The discharge coefficient ( ) for the nozzles with various divergence angles (ɛ) 
Fig.13 Effect of divergence angle on Mach number 
Fig.14 Comparison of internal static temperature distributions along the nozzle sidewall for various divergence angles 
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 halfangles on the nozzle performance at different nozzle pressure ratios (NPR) in the range of 610 are investigated numerically. Three different twoequation turbulence models are used for simulating the compressible gas flow inside the twodimensional convergingdiverging 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 halfangle ( ) 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:
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, NonAxisymmetric, 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 Overexpanded Nozzle, 45th AIAA Aerospace Sciences Meeting and Exhibit, Nevada, AIAA, (2007).
[4]. A.Hamed,C.Vogiatzist, Overexpanded TwoDimensional ConvergentDivergent Nozzle Performance, Effects of ThreeDimensional Flow Interactions, Journal ofPropulsion and Power,14, 234240, (1998).
[5]. A.Hamed, C.Vogiatzist, ThreeDimensional Flow Computations and Thrust Predictions in 2DCD Overexpanded Nozzles, AIAA paper 970030, (1997).
[6]. A. Balabel, A.M. Hegab, S. Wilson, M. Nasr, S. ElBehery, 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 HighSpeed Dual Stream Jets, AIAA 20033311, (2003).
[8]. C. F.Chenault, P.Beran, Numerical Investigation of Supersonic Injection Using a ReynoldsStress Turbulence Model, AIAA Journal, 37, 12571269, (1999).
[9]. H.Zhang, R.So, T.Gatski, C. Speziale, a NearWall SecondOrder Closure for Compressible Turbulent Flows, NearWall Turbulent Flows, edited by R. So, C. Speziale, and B. Launder, Elsevier, New York, 209218,(1993).
[10]. M.A. Dembowski, N.J. Georgiadis, an Evaluation of Parameters Influencing Jet Mixing Using the WIND NavierStokes Codes, NASA/TM211727,(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, InterLaboratory 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. ElBehery, Assessment of Turbulence Modeling for Gas Flow in TwoDimensional 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 ConvergingDiverging Nozzles, Journal of Propulsion and Power, 13, 753762, (1997).
[17]. M.L. Mason, L.E. Putnam,J.R. Richard, The Effect of Throat Contouring on Two Dimensional ConvergingDiverging 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 aReynoldsStress Turbulence Closure, J. Fluid Mech., 68, 537566, (1975).
[20]. B. E. Launder, D. B. Spalding, Mathematical Models of Turbulence, AcademicLondon, 169189, (1972).
[21]. F. R.Menter, TwoEquation EddyViscosity Turbulence Models for Engineering Applications, AIAA journal, 32, 15981605, (1994).