Effect of Mesh Size on the Viscous Flow Parameters of an Axisymmetric Nozzle

  • cc icon
  • 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 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.


  • KEYWORD

    Supersonic flow , Viscous flow , Finite volume , Nozzle

  • 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

    image

    (Schlichting, 1979).

    The Navier-stokes equations in a flux-vector formulation in the Cartesian coordinate system is given by

    image

    Where the vectors and are given by

    image
    image
    image
    image

    The heat flux vector q has three components, qx qy and qz, given by Fourier’s law of heat conduction, which relates the heat flux to the local temperature gradient.

    image

    Where k denotes the coefficient of thermal conductivity and is a function of the Prandtl number Pr, viscosity μ and specific heat at constant pressure cp

    image

    The energy per unit of mass e is defined as the sum of internal energy and kinetic energy:

    image

    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:

    image

    where mes(Ci, j ) is the measurement (in m3) of an infinitely small volume at a center (i, j). aire(Ci, 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(Ci, j ) and mes(Ci, 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:

    image
    image
    image
    image

    4. Discretization in Time

    The present numerical method is based on an explicit approach in time and space. The time step Δt is:

    image

    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:

    image

    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:

    image

    where Vn obtained from V, via the rotation R, through:

    image

    where:

    image
    image

    The overall transformation R is written as

    image
    image

    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:

    image

    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:

    image
    image

    where

    image

    , 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

    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:

    image

    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 rdiv = 0.05m for a throat radius r* = 0.01m. The diverging length was thus ldiv = 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 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.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 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.

  • 1. Ferziger J. H, Peric M 2002 Computational Methods for Fluid Dynamics P.217-259 google
  • 2. Goudjo J, Desideri A 1989 A Finite Volume Schemes to Resolution an Axisymmetric Euler Equation(Research report INRIA 1005) google
  • 3. Haoui R 2009 Application of the finite volume method for the supersonic flow around the axisymmetric cone body placed in free stream. [The 14th International Conference on Computational Methods and Experimental Measurements] P.379-388 google
  • 4. Haoui R 2010 Finite volumes analysis of a supersonic non-equilibrium flow around the axisymmetric blunt body. [International Journal of Aeronautical and Space Sciences] Vol.11 P.59-68 google doi
  • 5. Haoui R, Gahmousse A, Zeitoun D 2001 Chemical and vibrational nonequilibrium flow in a hypersonic axisymmetric nozzle. [International Journal of Thermal Sciences] Vol.40 P.787-795 google doi
  • 6. Haoui R, Gahmousse A, Zeitoun D 2003 Condition of convergence applied to an axisymmetric reactive flow. [The 16th French Congress of Mechanic] google
  • 7. Hoffmann K. A 1995 Computational Fluid Dynamics for Engineers Vol. II. Chap. 14. Engineering Education System P.202-235 google
  • 8. Schlichting H 1979 Boundary-Layer Theory google
  • 9. Van Leer B 1982 Flux-vector splitting for the Euler equations. In E. Krause ed. [Eighth International Conference on Numerical Methods in Fluid Dynamics. Lecture Notes in Physics Vol. 170.] P.507-512 google doi
  • [Fig. 1.] Boundary conditions.
    Boundary conditions.
  • [Fig. 2.] Grid of solution domain.
    Grid of solution domain.
  • [Fig. 3.] Grid with refinement.
    Grid with refinement.
  • [Fig. 4.] Velocity profile with various residues.
    Velocity profile with various residues.
  • [Fig. 5.] Velocity profile with various meshes size.
    Velocity profile with various meshes size.
  • [Fig. 6.] Velocity profile with various refinements near the wall.
    Velocity profile with various refinements near the wall.
  • [Fig. 7.] Profile of shear stress with various refinements near the wall.
    Profile of shear stress with various refinements near the wall.
  • [Fig. 8.] Viscous stress at the wall with refinement.
    Viscous stress at the wall with refinement.
  • [Fig. 9.] Temperature profile with various refinements.
    Temperature profile with various refinements.
  • [Fig. 10.] Temperature contours.
    Temperature contours.
  • [Fig. 11.] The comparison of Mach number contours.
    The comparison of Mach number contours.