Underexpanded jet

Underlying Flow Regime 1-01 © copyright ERCOFTAC 2004

Test Case

Brief description of the study test case

For the Moderately Under-expanded Jet, Seiner & Norum (1979,1980):

Jet flow is air, exhausting to atmosphere

Flow is from a convergent-divergent nozzle

Exit Mach number, M₁=2

Stagnation pressure, P₁/P_∞ = 1.45

Exit diameter, d=5.08 cm

Results available show the radial variation of Mach number at downstream distances of z/d =0.25, 0.75, 1.25, 1.75, 2.25, 2.75, 3.25, 3.75, 4.25, 4.75, 5.25, 5.75 (the four entries in bold type are compared with CFD by Cumber et al, 1994)

Results also show the magnitude of the fluctuating streamwise velocity and pressure with downstream distance, both of which are compared with CFD.

For the Highly Under-expanded Jet, Donaldson & Snedecker (1971):

(for illustration, see general arrangement in Figure 1, which was taken from this paper).

Jet flow is air, exhausting to atmosphere

Flow is from a convergent nozzle

Exit Mach number, M₁=1

Stagnation pressure, P₀/P_∞ = 6.56

Exit pressure, P₁/P_∞ = 3.57

Exit diameter, d=0.511 inch

Results available show the radial velocity profiles at downstream distances of z/d =1.96, z/d =3.92, z/d =7.32 and z/d =11.7 (all of which are compared with CFD by Cumber et al, 1995).

The first Mach disk is at z/d =1.58, a second normal shock is present at z/d =3.9, and the flow on the axis is just subsonic at z/d =23.5

Results are also available for the axial variation of velocity, pressure, and jet half-width.

Test Case Experiments

For the Moderately Under-expanded Jet: Seiner & Norum (1979,1980):

The experiments were performed in the unheated blowdown apparatus at NASA Langley Research center. This apparatus was able to provide, continuously, 4.1kg/s of dry air, with pressure control accurate to 0.3%. Static pressure measurements were made with a supersonic, conical, static tube; and this probe was fitted to a drive mechanism, with positioning accurate to 0.025mm. The pressures (static or total) were recored onto a computer via a digital multimeter, with time averaging chosen to reduce data scatter to the (0.2%) accuracy of the pressure transducers.

The nozzle used was a convergent-divergent nozzle. Geometrical details of the shape are not referenced. However the authors state that they used design guidelines (unreferenced) and a method of characteristics to ensure that the flow was parallel and supersonic at the stated Mach number at the nozzle exit. There is no reason to suspect that the nozzle was not well designed and gave the stated performance. This being the case, CFD simulations can be undertaken without geometrical details of the shape of this nozzle: The domain boundary can be coincident with the nozzle exit. Because the flow is choked at the nozzle exit, no information from the downstream computation can propagate upstream into the nozzle. This nozzle had an exit diameter of 5.08cm, and the exit Mach number was measured to be M=1.99.

Seiner and Norum (1979) present their data of the variations of both total and static pressure along the jet centreline, along with radial Mach number profiles at intervals of Δ(z/d)=0.50 from z/d=0.25 to z/d=5.75 (shown later, Figs 3 & 5).

The second paper from these authors, Seiner and Norum (1980) presents further data from their apparatus, obtained using hot-film anemometry, showing the fluctuating streamwise velocity along the jet axis for 0 < z/d < 10 (shown later, Fig 4).

This flow appears to have been essentially as per the target conditions. No checks were made on globally conserved quantities, however, in practice, the only relevant parameter that could have been considered was axial momentum.

For the Highly Under-expanded Jet: Donaldson & Snedeker (1971) :

Air was supplied from a storage tank, via an automatic regulator and settling chamber, to a convergent nozzle of exit diameter d=0.511 inch, before venting to atmosphere. No geometrical details are given for the convergent part of the nozzle. However there is no reason to suspect that the nozzle was not well designed and gave the performance stated below. CFD simulations can again be undertaken without geometrical details of the shape of this nozzle.

The local velocities were computed from Pitot and static pressure measurements. These probes were made from a 0.032inch diameter stainless-steel tube. Experiments were performed for all three types of jet (as illustrated in Figure 1). Here, the highly under-expanded jet will be studied, with a pressure ratio of P₀/P_∞ = 6.56, a nozzle pressure of P₁/P_∞ = 3.57, and the exit Mach number was M=1.00.

Measurements are available showing the radial velocity profiles at z/d=1.96, 3.92, 7.32 and 11.7, along with the axial decay of velocity, pressure and the velocity half-width.

Photographic observations were made, which showed the presence of a normal shock (Mach disk) at a downstream location of z/d=1.58. The first radial velocity measurements were made just after this point, hence exhibiting a subsonic velocity on the jet axis (c.f. Figure2). A second normal shock was also observed from the pressure trace at a downstream distance of z/d=3.9m, again explaining the small subsonic region on the jet axis measured just beyond this point at z/d=3.92m. Beyond this point, the jets axis was found to remain supersonic until z/d=23.5 (though with the peak velocity occurring off-centre).

No quantified information is given about the numerical accuracy of the measurement technique, although a verbal discussion is provided. The authors note that the jet exhibited a flapping instability, which peaked at certain pressure ratios; and also that, of course, the measured profiles of the jet were further influenced by the effects of turbulence. Hence the authors comment that, firstly, their results are probably at their most accurate in regions where the turbulence is at its lowest, and secondly, that it is difficult to define a suitable averaging period — a long average would show the mean effects of the flapping instability, but this would be expected to be different from any snapshot image of the jet.

This flow also appears to have been essentially as per the target conditions, with the exception of the flapping instability described above. As with the moderately under-expanded jet above, no checks were made on globally conserved quantities.

CFD Methods

5.1 Description of Cumber et als (1994 and 1995) CFD work.

In the following sub-sections section, the CFD methodology used by Cumber et al. (1994 and 1995) will be described.

5.1.1 Numerical Methods (Cumber et als CFD)

This work required the solution of a series of six partial differential equations, which were implemented into a modified version of the CFD code ‘COBRA’ from Mantis Numerics Ltd. The underlying equations were applied in a cylindrical co-ordinate frame, and assuming axi-symmetry. Convergence was achieved by solving the time-dependant form of the equations, via a time-marching scheme, until the steady state was reached, although the point at which convergence was declared is not quantified.

A second-order accurate finite volume scheme was used to integrate the equations, on a co-located grid. For the diffusion and source-terms, the central differencing scheme was applied; whereas for the advective and pressure terms, a second-order accurate version of Godunovs method was used. In order to increase the convergence rate, a ‘diagonalised implicit’ scheme was adopted. This calculates a second-order accurate, explicit, rate of change; then applies these diagonalised implicit operators to make the scheme implicit for sound waves and diffusion. Hence, allowing the solution to proceed using much larger time steps than a purely explicit scheme would require.

5.1.2 Turbulence Modelling (Cumber et als CFD)

Correct computation of this flow is clearly dependent on correct treatment of variable-density effects. Here, the density-weighted (Favre) averaged approach (Jones & Whitelaw, 1982) is applied, along with the classic eddy-viscosity model (Jones & Launder, 1972), and solution of the energy equation. The constants applied in the model were the standard ones, shown to be applicable to a wide range of flows (Jones & Whitelaw, 1982), and for the energy equation, the turbulent Prandtl number (σ_E) was set at 0.9. Air was assumed to be an ideal gas.

The standard eddy viscosity model defines the turbulent viscosity (μ_t) from the solution of equations for the turbulence energy (k), and its dissipation rate (ε), thus :

\mu _{t}=C_{\mu }\rho {\frac {k^{2}}{\varepsilon }}

\left(1\right)

Whereas this model has been (and still is) extensively applied for incompressible flows, there are known problems for compressible flows — the turbulence energy (k) will be overpredicted, leading to increased μ_t and excessive turbulent mixing. It has been shown from a DNS study by Lee et al. (1991) that instantaneous ‘eddy shocklets’ exist where high and low speed flows mix, and this led to a suggested modification from Sarkar et al. (1991) to the eddy viscosity model. Sarkar suggests that an additional sink term (ε_c) is included to account for the turbulence dissipation from these eddy shocklets :

\varepsilon _{c}=\alpha M_{t}^{2}\varepsilon

\left(2\right)

where α is a constant, set to unity, and M_t is the turbulent Mach number :

M_{t}={\frac {\sqrt {2k}}{a}}

\left(3\right)

(k being the turbulence energy, and a the local speed of sound).

This extra dissipation term, (ε_c), is applied into the classic (Jones & Launder, 1972) eddy-viscosity model, and manifests itself in two places. Firstly, obviously as an extra sink term for the turbulence energy equation :

S_{k}=-\rho M_{t}^{2}\varepsilon

\left(4\right)

and secondly, in the calculation of the eddy viscosity, from eq. (1) (effectively adding ε_c to ε):

\mu _{t}=C_{\mu }\rho {\frac {k^{2}}{(1+M_{t}^{2})\varepsilon }}

\left(5\right)

The results presented in Section 6 show a comparison between the classic eddy viscosity model, and the above (Sarkar er al. 1991) modified eddy viscosity model.

5.1.3 Boundary Conditions (Cumber et als CFD)

The CFD simulations were performed on a two-dimensional axi-symmetric grid

The jet centreline boundary (z axis at r=0) was modelled as a symmetry boundary.

The free boundary parallel to the z-axis (r=r_max) was modelled as a fixed pressure boundary (which in practice accurately describes the real boundary condition).

The downstream boundary, at z=z_max, was also modelled as a fixed pressure boundary (which in practice accurately describes the real boundary condition)

The jet inlet boundary, at z=0, was modelled as a hole (for the jet), in a flat plate. The jet inlet profiles were all defined as ‘flat’ across the jet (which is not unreasonable), and a turbulence energy was assumed taking 5% of the jet velocity, and its dissipation rate from a length scale of 0.05d. These assumptions for turbulence levels were subsequently shown to be un-influential except before the first shock.

The location of the two pressure boundaries was tested to ensure that they were not too close to the jet to affect the computed results.

5.1.4 Grid Used (Cumber et als CFD)

Cumber et al. tested the grid resolution needed for grid-independent results, and concluded that, if using a regular mesh, a grid spacing of Δr =Δz = d/64 was needed, in order to capture the shock structure. This resulted in 1.23 x 10⁵ grid nodes for the Seiner & Norum Jets.

However, this number of grid cells was reduced to 5.2x10⁴ grid nodes by the use of an adaptive mesh. This algorithm overlayed successively finer meshes (each containing cells of half the spacing of its predecessor), up to a maximum of five levels. The authors state that this algorithm is described in more detail in the paper of Falle & Giddings (1992).

5.1.5 Discussion (Cumber et als CFD)

The CFD methodology employed by Cumber et al. (1994 and 1995) appears to have been performed to a high standard. The modelling, and numerical techniques, used were recorded in their paper, and the key points noted above.

A plot of the grid used (especially the resulting adapted mesh) is not shown, however this is likely to have been for practical reasons, since the fine mesh required would not have reproduced well in print. However, the paper does present results to show the effects of different mesh densities, with these results converging towards the chosen mesh size. Indeed it is recorded that for the Donaldson & Snedeker jet, there was only a 2% difference in peak velocity between a mesh spacing of d/32, and the chosen d/64.

The overall size of the domain was tested to ensure this did not affect the results, and where assumptions have been made (for example in the level of turbulence energy), the choices made have been recorded.

5.2 Description of Bartosiewicz et als (2002) CFD work.

In the following sub-sections section, the CFD methodology used by Bartosiewicz et al. (2002) will be described.

5.2.1 Numerical Methods (Bartosiewicz et als CFD)

The governing equations were solved using a second-order accurate finite volume scheme, on a co-located computational mesh. A coupled approach was used for the conservation equations (continuity, momentum and energy), whereas the turbulence terms were solved separately. Diffusion was calculated using the second-order central difference scheme.

Convergence was achieved using a time-marching approach until the steady-state was reached, and convergence was declared when the normalised global mass imbalance was less than 10^-5, and the energy imbalance less than 10^-4.

The flow was first computed using a standard k-ε scheme, and the resulting solution used for initialisation of the Reynolds stress model.

5.2.2 Turbulence Modelling (Bartosiewicz et als CFD)

Bartosiewicz et al. apply a second moment closure methodology to solve the RANS equations. Taking the transport equation for the turbulent stresses $({\overline {u_{i}u_{j}}})$ :

{\frac {\partial }{\partial t}}\left(\rho {\overline {u_{i}u_{j}}}\right)+{\frac {\partial }{\partial x_{k}}}\left(u_{k}\rho {\overline {u_{i}u_{j}}}\right)=P_{ij}+D_{ij}^{L}+D_{ij}^{T}+\phi _{ij}+\varepsilon _{ij}

\left(6\right)

The first two terms on the right hand side of the equation represent the production and laminar diffusion terms, and can be implemented exactly into the solver thus :

P_{ij}=-\rho \left({\overline {u_{i}u_{k}}}{\frac {\partial u_{j}}{\partial x_{k}}}+{\overline {u_{j}u_{k}}}{\frac {\partial u_{i}}{\partial x_{k}}}\right)

\left(7\right)

D_{ij}^{L}={\frac {\partial }{\partial x_{k}}}\left[\mu {\frac {\partial }{\partial x_{k}}}{\Bigl (}{\overline {u_{i}u_{j}}}{\Bigr )}\right]

\left(8\right)

For the turbulent diffusion term $(D_{ij}^{T})$ Bartosiewicz declares that they have used the same model as described by Lien and Leschziner (1994). On inspection of this sub-reference, the modeling used appears to be the standard GGDH (Generalised Gradient Diffusion Hypothesis) model of Daly & Harlow (1970) :

D_{ij}^{T}=C_{s}{\frac {\partial }{\partial x_{k}}}\left({\frac {k}{\varepsilon }}{\overline {u_{k}u_{j}}}{\frac {\partial }{\partial x_{j}}}{\overline {u_{i}u_{j}}}\right)

\left(9\right)

, C_s=0.22 [9]

Likewise for the Pressure-Strain term (), Bartosiewicz refers the reader to Lien and Leschziner (1994), in which, the reader finds the classic ‘Basic’ second moment closure, containing destruction-of-production and redistribution components :

C₁=1.8, C₂=0.6 [10]

NB: Lien and Leschziner also include wall reflection terms in their pressure-strain model, however, since this is a free flow, away from any wall boundaries, one presumes Bartosiewicz omitted these terms.

Finally, the dissipation term (ε_ij) was assumed isotropic, though including the extra contribution to account for dissipation in the eddy shocklets (as Sarkar, 1991) :

[11]

with M_t as equation (3)

5.2.3 Boundary Conditions (Bartosiewicz et als CFD)

The computational domain used was 2D axi-symmetric, extending 20 diameters downstream and 4 diameters from the jet axis.

The jet centerline boundary (z axis at r=0) was modeled as a symmetry boundary.

The free boundary parallel to the z-axis (r=r_max) was modeled with a parallel co-flow.

The inlet face (z=0), beyond the inlet jet itself (r>r_jet), was also modeled with a parallel co-flow.

The downstream boundary, at z=z_max, was modeled as an outlet. In subsonic regions, this was set to ambient pressure, in supersonic regions, the properties were ‘extrapolated from the interior of the domain’.

For the jet inlet itself, the turbulent intensity was set to 5%, and the turbulent viscosity ratio was set to 2. However, Cumbers observations (section 5.1.3) were that the overall prediction of the jet was not sensitive to the definition of inlet turbulence.

5.2.4 Grid Used (Bartosiewicz et als CFD)

Bartosiewicz et al. adopt an adaptive gridding strategy on a triangular mesh, which refined the mesh in areas of strong Mach number and density gradients. The resulting mesh contained 30,000 cells. The quality of the mesh was assessed by comparing the axial and radial Mach number profiles between successive meshing steps. When the difference in profiles became less than 5%, the resulting solution was deemed grid-independent.

5.2.5 Discussion (Bartosiewicz et als CFD).

Bartosiewiczs CFD methodology appears to be robust, although not going into the same detail as Cumber (for example the software solver used is not named, nor is any reassurance given that the location of outer domain boundary does not influence the final results). The modelling used is described, although leaving the reader to check sub-references in order to try and ascertain the equations actually solved. The discretisation, gridding, and convergence strategies are declared.

Bartosiewiczs work was only in fact published after this QNET-CFD document was originally written. However, the inclusion by Bartosiewicz of more advanced turbulence models, which are used to compute the same experimental flows as already described herein, warranted a partial re-write of this QNET document.

Contributors: Chris Lea - Health & Safety Laboratory