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.

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 c_v is the specific heat at constant volume

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³) 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:

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.

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 W_E (Euler variable) is written W_E ^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(W^R) and F^?_VL(W^R), where WR is defined like the transform of W by rotation R, are:

where

, u_n and v_n are the velocity at the reference mark of the interface.

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.

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.

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).

The conditions at the symmetry boundary are: (i) no flow across the boundary and (ii) no scalar flux across the boundary.

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.

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

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.

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.