INTERNATIONAL JOURNAL OF A publication of IIETA

HEAT AND TECHNOLOGY

ISSN: 0392-8764
Vol. 35, No. 3, September 2017, pp. 673-677
DOI: 10.18280/ijht.350328
Licensed under CC BY-NC 4.0 http://www.iieta.org/Journals/IJHT

Numerical simulation of the compressible flow in convergent-divergent

nozzle

Souici Mohamed1*, Aouissi Mokhtar1, T. Ben Chatti2

1
Laboratoire de Mécanique, Département de Mécanique, Universitéde Laghouat, Laghouat
03000, Algeria
²Ecole normale supérieure de Laghouat, Laghouat 03000, Algeria

Email: [email protected]

ABSTRACT

This work concerns the numerical simulation of a non-reactive turbulent flow of a compressible fluid air
(assumed as a perfect gas), using the MacCormack scheme for the discretization of the terms in the various
extras equations of the model, and this in a convergent divergent nozzle, in order to study the pressure
influence on the flow characteristics such as speed pressure and temperature ... etc. And particularly
delamination phenomenon. These results are represented by field’s contours and variable flow profiles of the
axis of the nozzle and near the wall the system of equations governing the flow is solved using the finite
volume method. The turbulence model used is the two transport equations of type (𝑘, 𝜔) SST (Shear-Stress-
Transport) of Menter.

Keywords: Converging-diverging Nozzle, Turbulence, Shock Wave, Supersonic, Compressible Flow, Finite
Volume.

1. INTRODUCTION 2. TURBULENCE MODELLING

The problems of fluid dynamics are often difficult to solve The k-ω SST shear-stress-transport model [8] combines
because the system of equations governing the phenomenon several desirable elements of existing two-equation models.
which is strongly nonlinear system. It is difficult to find exact The two major features of this model are a zone weighting of
solutions. However, digital technology especially the model coefficients and a limitation on the growth of the eddy
calculation applied to fluid dynamics change was successful, viscosity in rapidly strained flows. The zone modelling uses
these successes are due to the close interaction between Wilcox’s k-ω model near solid walls and the standard k-ε
theory, numerical simulation and experimentation in fluid model near boundary layer edges and in free-shear layers.
dynamics. On the one hand, experience is essential to test the This switching is achieved with a blending function of the
hypotheses and the results that emerge from the theory [1], model coefficients. The shear stress transport modelling also
on the other hand, the theory is needed to explain the results. modifies the eddy viscosity by forcing the turbulent shear
Emerge from the theory, on the other hand, the theory is stress to be bounded by a constant time the turbulent kinetic
needed to explain the results. Numerical simulation is energy inside boundary layers. This modification improves
independent of experience, it is necessary for the validation the prediction of flows with strong adverse pressure gradients
of experimental results. The results are presented for nozzle and separation.
flow subsonic - supersonic. Various digital tests presented in The k-ω SST model uses two additional transport
this study relate to the influence of the variation of the equations to describe the turbulence as summarized below:
geometry of the nozzle such that the angle of divergence as
well as the effect of the variation of the input of the D k U i       k  t  k 
characteristic quantities of the pressure flow.   ij   * k    (1)
Dt x j xk  xk 
General equations the balance equations of mass,
momentum and energy for asymmetrical laminar flows of
viscous compressible fluids can be written in the following
vector conservative form.

673
 D  
  ij
U i
   2 


 t
    k  k 

uv  C k (7)
Dt vt x j  xk   xk 
  (2) These two remarks led to introduction of weakly non-
2  1  F 1   2 k  linear turbulence models in which the 𝐶𝜇 factor is allowed to

ù x j x j vary according to:

The left-hand side of the Esq. (1) and (2) is the Lagrangian  
 
D  1
derivative  U i 
 xi
and the turbulent stress tensor C   min  0.09,  (8)
Dt t  1 
ô𝑖𝑗 is given by: 
  
A0  A s S 2  A   2 2 


 U i U j 2 U k  2 With:
 ij   t     ij    k  ij (3)
 x j  xi 3  x k  3
 1 2
S  2 S ij S ij  S 2 ,
The function 𝐹1 is designed to blend the model coefficients   *
3 kk
of the original k-ω model in boundary layer zones with the   
transformed k-ω model in free shear layer and free stream   1  U i  U j 
 S ij
zones. This function is expressed in term of local variables as:  2   x j  xi 
 
 (9)
   1
   2 k 500v  4   2k    
2 ij  ij ,
F 1  tanh  min  max   ,  (4)  *
  0.09 y  y 2  CD k y 2    
      1  U i U j 
 ij    
 2   x j  xi 
where CDkω is a cross diffusion term added in Eq. (2)  
According to Bradshaw’s assumption the eddy viscosity is
defined in the following way: And 𝐴0 = 0, 𝐴𝑠 = 2√3 , 𝐴𝛺 = 0 in this case the Bradshaw
coefficient (0.31) is substituted by Cμ1/2 in the formulation of
a1k the eddy viscosity.
μt  (5)
max  a1ω, Ω  F 2

2. NUMERICAL METHOD
where F2 is a function that is one for the boundary-layer
flows and zero for the free shear layers?
With:

  2ij ij

1  U i U j 
ij  


2   x j  xi 


 , F 2  tanh ar g 2
2 
 2 k 500v 
ar g 2  max  
 0.09 y  y 2 
Figure 1. Nozzle geometry
2.1 Reliability condition in turbulence models
The Navier-Stokes equations are solved on a
The two-equation turbulence models are based on the computational domain of variables ζ and η (transformed
Boussinesq assumptions where the Reynolds stresses is coordinates of the physical domain), by the use of finite
expressed as a linear function of the mean strain tensor: volumes predictor-corrector. The new system of equations is
solved by using MacCormack's explicit-implicit scheme [4].
2C   k 2  1  2 This algorithm is second-order accurate in space and time.
  u iu j   S ij  S ii    k  ij (6) The basic discretization for the convective fluxes is modified
  3  3 to account for the physical properties of information
propagation, as done initially by Steger & Warming [5]. The
where Cμ As shown by Moore & Moore [3] these equations flux splitting is made second order accurate, but in shock
can give negative values of the normal stress if Sll is too large. regions where it is lowered to first order. The viscous terms
Bradshaw has noticed that in two-dimensional boundary are cantered and the axisymmetric source terms are
layers submitted to a strong pressure gradient the shear stress integrated at the Centre of each control volume in both the ζ
was approximately proportional to the turbulent kinetic and η directional sweeps. To reach a steady-state solution
energy with: with a minimum number of iterations, the explicit

674
discretization is complemented with an implicit numerical [7] path that begins just after the col of the nozzle the density
approximation which is free from stability conditions. thus undergoes a small sharp increase then continues its
decrease until the exit of the nozzle .

4. DISCUSION VALIDATION OF RESULTS

The validation of the results of our work is very important

because it allows us to verify and compared it our results
obtained with numerical simulation by other experimental
results, theoretical and numerical. Our work has been
completed with the numerical results of Back and all [4].
They obtained profiles of ratio of the pressure in a supersonic
nozzle, as shown in Figures 2.a and 2.b they acceptable in
good agreement with our results.

Figure 4.a Contours of Static Pressure (pa)

Figure 2.a Ratio of the pressure comparison between our

works And that of Back and all [4]

Figure 4.b Evolution of Static Pressure center and cell wall

The pressure drop inside the nozzle schematically in Figures
4.a and 4.b. qualitatively this profile follows almost the same
shape of density profile

Figure 2.b Ratio of the pressure comparison between our

works and that of Back and all [4]

Figure 5.a Contours of Temperature (K)

The temperature distribution and its evolution in shown

respectively in Figures 5.a and 5.b , it was seen that the
temperature within the nozzle undergoes three main phases
of decline, the first phase is stable in the form of a straight
line, it is at the level of converging up to the col of the nozzle,
Figure 3. Contours of Density (kg / m3) the second phase of decline is a sharp decrease in
temperature at the col and the last phase is in the divergent, it
Figure 3 notes in this Figure that the profile of the density continues to decrease until the exit of the nozzle with small
(colors change) takes two different paths, the first path is sharp increase in the input of the divergent and the other in
from the inlet of the nozzle up to its col, the density in this the middle of the, on the other hand the proportion steady
portion is almost constant at its maximum value. The second

675
parietal distribution along the nozzle because friction observations for the evolution of Mach number except
between the flow and the walls. parietal disturbance are observed in the vicinity of the col due
to friction between the fluid and the wall, and the convergent-
divergent profile of the nozzle accelerates the gas to a
subsonic velocity to supersonic velocity.

5. CONCLUSIONS

In recent years, remarkable progress has been made in the

domain of propulsion in general and aerospace propulsion in
particularly, but all this progress and the studies of these
physical phenomena at high velocity and very high
temperature in nozzles often misunderstood, this applies
especially several theoretical, numerical and even
experimental studies to understand the maximum of these
flow phenomena in nozzles used. These progress; numerical
Figure 5.b Evolution of Temperature center and cell wall (K)
study of flow inside the nozzles that can address the
problems of physical phenomena without losing a lot of
money and make a lot of the time, so we can say that the
numerical simulation of turbulent compressible flows is a
better means to better understand the physical phenomena.
This computational result examined the effects of several
parameters on the dynamic and thermal characteristics of
flow through a cooled convergent-divergent nozzle. The
computational results indicated the following:
For an initial pressure Po <Pc (critical pressure) It is noted
that the separation of the position is affected by this change
and moves the nozzle throat; by cons if one increases the
initial pressure as Po> Pc, then the release moves the lip of
the nozzle until it disappear which is desirable. Therefore,
delamination which is an undesirable phenomenon appears in
Figure 6.a Contours of Mach number some flows, particularly when the initial pressure is below
the critical pressure, this could be the cause of degradation of
aerodynamic performance of the nozzle, it yields to losses
and efficiency and causes noise and vibration as structures
where it is advantageous to avoid this by increasing the initial
pressure and make a wise choice of geometry.

REFERENCES

[1] Back L.H., Massier P.F., Gier H.L. (1964).

Convective heat transfer in a convergent-divergent
nozzle, Int. J. Heat Mass Transfer, Vol. 7, pp. 549-
568.
[2] Lutum E., Wolfersdorf J., Semmler K., Dittmar J.,
Figure 6.b Evolution of Mach number center and cell wall Weigand B. (2001). An experimental investigation of
film cooling on a convex surface subjected to
The evolution of the Mach number in (central and parietal) favourable pressure gradient flow, Int. J. Heat Mass
of the nozzle shown in Figures 6.a and 6.b respectively, we Transfer, pp. 44939-951.
observe that the regime at the inlet of the nozzle remains [3] Lebedev V.P., Lemanov V.V., Terekhov V.I. (2006).
almost steady or invariable up to the first contact of the Film-cooling efficiency in a laval nozzle under
abscissa of point of tangency at the level of convergent and conditions of high free stream turbulence, Journal of
we can say that the regime is subsonic because the value of Heat Transfer, Vol. 128, pp.571-579.
the Mach number _ is strictly less than unity, then it follows [4] Back L.H., Massier P.F., Gier H.L. (1965).
by a sharp increase in the vicinity of col of the nozzle exactly Comparison of measured and predicted flows through
between the first point of tangency at the level of converging conical supersonic nozzles, with emphasis on the
and the first point of tangency in the divergent in this zone transonic region, AIAA Journal, Vol. 3, No. 9, pp.
the regime becomes transonic saw his Mach number, as there 1606-1614.
was a small sharp decrease in the center nozzle, at the inlet to [5] Cuffel R.F., Back L.H., Massie P.F. (1968). Transonic
the zone of the divergent ,the value of Mach number on the flow field in a supersonic nozzle with small throat
inside of the nozzle continues to increase up to the outlet of radius of curvature, AIAA Paper, Vo1. 7, No. 7, pp.
the nozzle, in this zone the regime is hypersonic.. The same 1364-1366.

676
[6] Sauer R. (1947). General characteristics of the flow 𝜌 density
through nozzles at near critical speeds, NACA TM Ω scalar measure of the vortices tensor
1147. Ω𝑖𝑗 vortices tensor
[7] Lorenzini G., Saro O. (2016). Analysis of water 𝑌 specific heat ratio
droplet evaporation through a theoretical numerical 𝑎1 brads how constant
model, International Journal of Heat and Technology, 𝑟 radius, radial coordinate, recovery factor
Vol. 34, Special Issue 2, pp. S189 - S198. DOI: 𝑥 axial coordinate
10.18280/ijht34S201 𝑇 temperature
[8] Wang X., Wang X.R., Zhang J.X., Wang J.G. (2015). 𝑀 Mach number
Nozzle diameter influence on spray characteristics in a 𝑃 pressure
constant volume combustion chamber, International 𝑃𝑟𝑡 turbulent Brandtl number
Journal of Heat and Technology, Vol. 2, No. 3, pp. 9- 𝑈𝑖 mean velocities
12. DOI: 10.18280/mmep.020303 𝑢𝑖 fluctuating velocities
ℎ heat transfer coefficient
𝑃𝑟 Brandtl number
NOMENCLATURE 𝐹1 , 𝐹2 auxiliary functions in turbulence model
𝑡 time
𝐶𝑃 specific heat at constant pressure Subscripts
𝐶𝑣 specific heat at constant volume 0 nozzle entrance condition
𝑅 gas constant 𝑡ℎ throat position
𝐷 diameter 𝑠𝑒𝑝 separation position
𝑘 turbulent kinetic energy 𝑤 parameters on the wall surface
𝜔 specific turbulent dissipation rate 𝑎𝑤 adiabatic wall
𝜇 dynamic viscosity 𝑒 free stream condition
𝜇t turbulent viscosity

677