The viscous flow in an axisymmetric nozzle was analyzed while accounting for the mesh sizes in both in the free stream and the boundary layer. The Navier-Stokes equations were resolved using the finite volume method in order to determine the supersonic flow parameters at the exit of the converging-diverging nozzle. The numerical technique in the aforementioned method uses the flux vector splitting of Van Leer. An adequate time stepping parameter, along with the Courant, Friedrich, Lewis coefficient and mesh size level, was selected to ensure numerical convergence. The boundary layer thickness significantly affected the viscous flow parameters at the exit of the nozzle. The best solution was obtained using a very fine grid, especially near the wall at which a strong variation of velocity, temperature and shear stress was observed. This study confirmed that the boundary layer thickness can be obtained only if the size of the mesh is lower than a certain value.
The nozzles are used at the exit of the shock tube in order to obtain supersonic flows for various tests. They also used in propulsion to obtain the thrust necessary to the displacement of the vehicles.
1. Introduction
The calculation of viscous flow in an axisymmetric converging-diverging nozzle is presented. A numerical technique that simulates the viscous supersonic flow and the boundary layer thickness in the nozzle, especially at the exit, was employed. Such as working fluid was air in standard state composed of 21% O2 and 79% N2, which is regarded as a perfect gas. The stagnation parameters were 2,000 K and 500 bars, thus the vibration and dissociation of molecules were neglected. These parameters were selected so that the flow at the exit of the nozzle was 300 K, 1 bar and yielded a Mach number of 5.
The non-linear partial derivative systems of equations,which governs viscous flow, were solved by Goudjo’s explicit unsteady numerical scheme Goudjo and Desideri (1989)and the finite volume method presented in previous works by Haoui (2001, 2009, 2010). Van Leer’s flux vector splitting method was also used (Van Leer, 1982). An adequate time stepping parameter and the mesh size level were selected to ensure numerical convergence (Haoui et al., 2003).The stationary solution obtained depended on the size of the mesh used in the numerical discretization (Haoui et al.,2003). Inviscid flow convergence was tested using a refined grid, which enabled us to obtain an exact solution. Then,we refined the grid again near the wall to determine the boundary layer thickness and other parameters.
2. Governing Equations
In a Newtonian fluid the viscous stresses are proportional to the rates of deformation. The three-dimensional form of Newton’s law of viscosity for compressible flows involves two constants of proportionality: dynamic viscosity μ, which relates stress to linear deformation, and viscosity λ, which relates stress to volumetric deformation.
For gases, a good working approximation for viscosity is
(Schlichting, 1979).
The Navier-stokes equations in a flux-vector formulation in the Cartesian coordinate system is given by
Where the vectors and are given by
The heat flux vector q has three components, q
x
q
y
and q
z
, given by Fourier’s law of heat conduction, which relates the heat flux to the local temperature gradient.
Where
k
denotes the coefficient of thermal conductivity and is a function of the Prandtl number Pr, viscosity
μ
and specific heat at constant pressure c
p
The energy per unit of mass
e
is defined as the sum of internal energy and kinetic energy:
T
is the temperature and
cv
is the specific heat at constant volume
3. Axisymmetric Formulation
The Institute National de Recherche en Informatique et en Automatique (INRIA) Sinus project developed a method that translates 3d axisymmetric flow to 2D axisymmetric flow using a technique involving domain disturbance.
The system of Eq. (1) can be written as:
where mes(C
i
,
j
) is the measurement (in m
3
) of an infinitely small volume at a center (i, j). aire(C
i
,
j
) is the surface of the symmetry plane passing through the center of the elementary volume. η
a
is the integrated normal. Goudjo presented a detailed calculation of, η
a
, aire(C
i
,
j
) and mes(C
i
,
j
) (Goudjo and Desideri, 1989). The third term of the equation expresses the axisymmetric flow condition. Flows,
W, E, F
and
H
are given by:
4. Discretization in Time
The present numerical method is based on an explicit approach in time and space. The time step
Δt
is:
Courant, Friedrich, Lewis (
CFL
) is a stability factor(Hoffmann, 1995).
V
is the velocity of the flow and a is the speed of sound.
Δx
is an incremental length of the mesh at the point (
i
,
j
). At each time step and for each point (
i
,
j
), the system of Eq. (15) can be written as:
Grid choice plays an important role in determining in the convergence of calculations. Therefore, sufficiently refined meshes are necessary at the places where the gradients of the flow parameters are significantly large.
5. Decomposition of Van-Leer
In this study, the Van-Leer decomposition (Van Leer, 1982)was selected, namely a decomposition of flows into two parts F
?
VL
and F
+
VL
. The flow at the interface between two cells must be calculated in order to utilize Van Leer’s decomposition method in 2d axisymmetric viscous flow through a nozzle.Moreover, we must know the normal direction of the interface; thus, the normal direction of the interface is indicated by a reference mark and an intermediary rotation R (Haoui, 2010) .
The vector
WE
(Euler variable) is written
WE R
at the new reference mark:
where V
→
n
obtained from V
→
, via the rotation
R
, through:
where:
The overall transformation
R
is written as
Moreover, at each interface (
i
+ 1/2), two neighbor states(
i
) and (
i
+ 1) are known. Thus, one can calculate the onedimensional flow F through the interface, total flow
f
(
W
, η
¯
) being deduced from F by applying the opposite rotation, as:
Thus, one component of flow f (F for example) can be used to define the decomposition of flow in two dimensions. The expressions of
F+VL
(
WR
) and
F?VL
(
WR
), where WR is defined like the transform of
W
by rotation
R
, are:
where
,
un
and
vn
are the velocity at the reference mark of the interface.
6. Boundary Conditions
CFD problems are defined by initial conditions and boundary conditions. In transient problems, the initial values of the flow variables must be specified at all solution
Boundary conditions.
points in the flow domain. The present work describes the implementation of the following most common boundary conditions in the discretised equations of the finite volume method: inlet, outlet, wall and symmetry,
Fig. 1
.
- 6.1 Inlet boundary conditions
The pressure and temperature were fixed at the inlet;however, the velocity module was extrapolated from the interior of the solution domain. Thus, the flow rate can be adjusted.
- 6.2 Body surface
The no-slip condition was applied to the body surface. The temperature gradient at the wall was zero, in accordance with the Fourier equation of heat conduction in the y-direction together with the assumption of zero heat flux at the wall. In this study, we assumed that the temperature at the wall was equal the stagnation temperature of the free stream. The wall shear stress was calculated by:
Here, we assumed that the coordinate of the unit vector
t
was in the direction of the shear force at the wall and the unit vector n was normal at
t
(Ferziger and Peric, 2002).
- 6. 3 Axis of symmetry
The conditions at the symmetry boundary are: (i) no flow across the boundary and (ii) no scalar flux across the boundary.
- 6. 4 Outlet boundary conditions
At the exit of the computational domain, the flow was supersonic and the values of the flow parameters were extrapolated from the interior values, including in the boundary layer.
7. Results and Interpretations
The nozzle comprised a convergent conicity of α
conv
= 45
°
followed by an arc of radius r=2r* and a divergent conicity of α
div
= 10
°
(
Figs. 2
and
3
). The stagnation pressure and temperature were 100 bars and 2,000 K, respectively. The Mach number desired at the exit of the nozzle was 5. The simple laws of a one-dimensional isentropic flow provided us a radius at the exit of the nozzle r
div
= 0.05m for a throat radius r* = 0.01m. The diverging length was thus l
div
= 0.228m.In our calculations we used several grid sizes beginning with
Grid of solution domain.
Grid with refinement.
Velocity profile with various residues.
116 × 10 (
Fig. 2
), in which 116 units were measured along the axis and 10 units along the radius. The mesh was then refined to sizes 223 × 20, 350 × 30, 466 × 40, 583 × 50 and 700 × 60 in order to observe how the refinement affected the results.
The residue values from which the results remain unchanged must be fixed. For this purpose, we used the 116× 10 grid since it was the least refined. The main parameter is the velocity profile which contained the boundary layer thickness, as shown in
Fig. 4
. We observed that the velocity profile was almost the same when the order of the residue was 10-5 to 10-6. Calculations were stopped when the residue was equal to 10-5.
The size of the grid of the calculation field from which the results remain unchanged must also be fixed, without refinement in the boundary layer. Thus, six grid sizes were tested for the same residue (
Fig. 5
). The velocity profile became flat near the wall,
r
= 0.05
m
, when the grid became more refined. The 350 × 30 grid was selected since it yielded good results and required less time computing.
The grid near the wall must also be refined in order to adequately simulate the flow parameters in the boundary
Velocity profile with various meshes size.
Velocity profile with various refinements near the wall.
layer, especially the velocity profile. Forty percent of the mesh near the wall must be multiplied by the constants 1,2, 3, 4, and 5 until the results show no changes. This process must be completed for the 12 meshes. The refinement thickness must be larger than the boundary layer thickness.We observed that the velocity profile became flat while approaching the wall, and the velocity on the axis almost remained unchanged. We selected the 350 × 78 grid, ni = 350 according to x and nj = 78according to y, for the final results.The boundary layer thickness at the exit of the nozzle can be deduced when the velocity reached 95% of the velocity on the axis. The boundary layer thickness was 10.7 mm, i.e. 21%of the radius. If the boundary layer was 99% of velocity on the axis, its thickness would become 18.3 mm and would consume 36% of the radius.
Another significant parameter in viscous flow must be presented, it is the stress τ
xy
.
Figure 7
shows the variations of shear stress along the radius according to the refinement of the grids in the boundary layer. This profile itself converged to the exact solution for the 350 × 78 grid. The intensity of the stress increased quickly while approaching the wall.The viscous stress at the wall can be calculated from all the stresses at the same point.
Figure 8
shows the variation of the stress τ
wall
along the wall of the nozzle with and without refinement in the boundary layer.
Profile of shear stress with various refinements near the wall.
Viscous stress at the wall with refinement.
Temperature profile with various refinements.
In regards to temperature profile, the solution also converged to the exact solution using the 350 × 78 grid (
Fig.9
), the wall of the nozzle was adiabatic and the profile of the temperature was thus perpendicular to the wall.
The flow in the nozzle was then compared to inviscid flow.
Figure 10
shows the temperature distribution in the nozzle.The thermal boundary layer thickness was visible near the wall.
A comparison with the inviscid flow is represented in
Temperature contours.
The comparison of Mach number contours.
Fig. 11
. The upper region of the chart represents the Mach contours of the viscous flow and the lower region represents the Mach contours of the inviscid flow. The boundary layer clearly influenced the flow parameters in the nozzle. For example, the Mach number at the exit was lower for the flow in the nozzle in comparison to that of the inviscid flow.
8. Conclusions
In conclusion, the results obtained in viscous flow strongly depended on the mesh size during the numerical calculations. The program converged due to the size of the meshes used; but, an exact solution was obtained only if the grid, especially near the wall, was very refined. The approximation by the infinite volumes method with the non-stationary scheme yielded good results. Our code was stable, consistent and the solution converged to the exact solution when the grid was very small. The exactitude of our code was carried out by using a 350 × 78 mesh with a residue of 10-5. We saw that the mesh size significantly influenced the flow parameters in the nozzle, and even on the wall stress. The computer codes which do not take into account the refinement of the grid near the wall, their results are not precise. The refinement of the grid in viscous flow has an effect on the results obtained.
Ferziger J. H
,
Peric M
2002
Computational Methods for Fluid Dynamics
3rd rev. ed.
Springer
Berlin
217 -
259
Goudjo J
,
Desideri A
1989
A Finite Volume Schemes to Resolution an Axisymmetric Euler Equation(Research report INRIA 1005)
National Institute of Research in Informatics and Automatic (NIRIA)
Haoui R
2009
Application of the finite volume method for the supersonic flow around the axisymmetric cone body placed in free stream.
Algarve Portugal
The 14th International Conference on Computational Methods and Experimental Measurements
379 -
388
Haoui R
2010
Finite volumes analysis of a supersonic non-equilibrium flow around the axisymmetric blunt body.
International Journal of Aeronautical and Space Sciences
11
59 -
68
DOI : 10.5139/IJASS.2010.11.2.059
Haoui R
,
Gahmousse A
,
Zeitoun D
2001
Chemical and vibrational nonequilibrium flow in a hypersonic axisymmetric nozzle.
International Journal of Thermal Sciences
40
787 -
795
DOI : 10.1016/S1290-0729(01)01265-0
Haoui R
,
Gahmousse A
,
Zeitoun D
2003
Condition of convergence applied to an axisymmetric reactive flow.
Nice France
The 16th French Congress of Mechanic
Hoffmann K. A
1995
Computational Fluid Dynamics for Engineers Vol. II. Chap. 14. Engineering Education System
Engineering Education System
Wichita KS
202 -
235
Schlichting H
1979
Boundary-Layer Theory
7th ed.
McGraw-Hill
New York
Van Leer B
1982
Flux-vector splitting for the Euler equations. In E. Krause ed.
Springer Berlin
Heidelberg
Eighth International Conference on Numerical Methods in Fluid Dynamics. Lecture Notes in Physics Vol. 170.
507 -
512
DOI : 10.1007/3-540-11948-5_66