Effect of Mesh Size on the Viscous Flow Parameters of an Axisymmetric Nozzle
 Author: Haoui Rabah
 Organization: Haoui Rabah
 Publish: International Journal Aeronautical and Space Sciences Volume 12, Issue2, p149~155, 30 June 2011

ABSTRACT
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 NavierStokes equations were resolved using the finite volume method in order to determine the supersonic flow parameters at the exit of the convergingdiverging 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.

KEYWORD
Supersonic flow , Viscous flow , Finite volume , Nozzle

1. Introduction
The calculation of viscous flow in an axisymmetric convergingdiverging 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 nonlinear 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 threedimensional 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 Navierstokes equations in a fluxvector 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 andc_{v} is the specific heat at constant volume3. 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 andH 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 VanLeer
In this study, the VanLeer 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
W_{E} (Euler variable) is writtenW_{E} ^{R} at the new reference mark:where V^{→}_{n} obtained from V^{→}, via the rotation
, through:R where:
The overall transformation
is written asR 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 flowf (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} (W^{R} ) andF^{?}_{VL} (W^{R} ), where WR is defined like the transform ofW by rotationR , are:where
,
u_{n} andv_{n} 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
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 noslip 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 ydirection 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 att (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 onedimensional 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
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 105 to 106. Calculations were stopped when the residue was equal to 105.
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.05m , 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
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.
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
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 nonstationary 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 105. 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.

[Fig. 1.] Boundary conditions.

[Fig. 2.] Grid of solution domain.

[Fig. 3.] Grid with refinement.

[Fig. 4.] Velocity profile with various residues.

[Fig. 5.] Velocity profile with various meshes size.

[Fig. 6.] Velocity profile with various refinements near the wall.

[Fig. 7.] Profile of shear stress with various refinements near the wall.

[Fig. 8.] Viscous stress at the wall with refinement.

[Fig. 9.] Temperature profile with various refinements.

[Fig. 10.] Temperature contours.

[Fig. 11.] The comparison of Mach number contours.