Nomenclature

$\alpha $

Angle-of-attack

c

body length

Cp

pressure coefficient $\dfrac{2}{{\gamma M_\infty ^2}}\left\{ {\dfrac{p}{{{p_\infty }}} - 1} \right\}$

$C{m_{q\;}}$

$\dfrac{{2{V_\infty }}}{c}\dfrac{{\partial Cm}}{{\partial q}}$

$C{m_{\alpha \;}}$

$\;\dfrac{{{V_\infty }}}{c}\dfrac{{\partial Cm}}{{\partial \alpha }}$

h(x,t)

instantaneous height of a point on the body

l

Moment arm

$\rho $

Fluid density

$p$

Pressure

${P_{b,s}}$

Pressure with blowing or on solid surface

$p_{b,s}^{\prime}$

Pressure perturbation due to unsteadiness

${P_{ref}}$

Pressure at the reference location

${P_\infty }$

Free stream pressure

R _B

Base radius

R _N

Nose radius

S

Reference area

$T$

Temperature

X

Axial distance

V

Velocity magnitude

R(x)

Local radius

Re

Reynolds number

$M$

Mach number

M _x

Local Mach number

q

pitch rate

t

time for particle to move from apex to station x

$\Delta {\varepsilon _{p\;}}\left( x \right)\;$

change in angle of inclination of tangent to (projected) streamline (path line) of a particular fluid element between vertex and point (x, R(x), $\varphi$ ) on body surface, measured in osculating plane of reference

θ _b (x)

local angle of inclination of body surface to axis of symmetry

θ₀

slope of the body at the designated apex

ϑ

pitch angle

$\varphi$

azimuthal direction (360-θ) measured clockwise from the starboard radius on the horizontal diameter plane through the vertex

2.0 Computational approach to steady flow field with and without frontal blowing

The conservative form of the compressible Reynolds-Averaged Navier-Stokes equations as used in the algorithm is as follows:

(1) \begin{equation}\frac{{\partial \rho }}{{\partial t}} + \frac{{\partial\! \left[ {\rho {u_j}} \right]}}{{\partial {x_j}}} = 0 \end{equation}

(2) \begin{equation}\frac{\partial }{{\partial t}}\!\left( {\rho {u_i}} \right) + \frac{\partial }{{\partial {x_j}}}\left[ {\rho {u_i}{u_j} + p{\delta _{ij}} - \left( {1 + \frac{{{\mu _t}}}{\mu }} \right)\!{\tau _{ij}}} \right] = 0 \end{equation}

(3) \begin{equation}\frac{\partial }{{\partial t}}\!\left( {\rho E} \right) + \frac{\partial }{{\partial {x_j}}}\left[ {\rho {u_j}E + {u_j}p - \left( {1 + \frac{{{\mu _t}}}{\mu }} \right){u_i}{\tau _{ij}} - \left( {{c_p}\frac{\mu }{{\Pr }} + {c_p}\frac{{{\mu _t}}}{{{{\Pr }_t}}}} \right)\!\frac{{\partial T}}{{\partial {x_j}}}} \right] = 0 \end{equation}

ρ in these equations refers to the density, p the pressure and T is the temperature. E is the total energy and u is the velocity vector in tensor format. τ _ij is the viscous stress tensor which may be written as:

(4) \begin{equation}{\tau _{ij}} = \mu \left( {\frac{{\partial {u_i}}}{{\partial {x_j}}} + \frac{{\partial {u_j}}}{{\partial {x_i}}} - \frac{2}{3}\frac{{\partial {u_k}}}{{\partial {x_k}}}{\delta _{ij}}} \right) \end{equation}

The above Navier Stokes equations, formulated in strong conservation form to be used in the structured meshes of the present study, are first cast into curvilinear coordinate scheme following Vinokur [Reference Vinokur11] transformation method. The numerical algorithm used to solve the Navier Stokes equations is an implicit factorisation finite centre difference scheme about a regular rectangular prism, which is described at length by Beam and Warming [Reference Beam and Warming12]. Local time linearisation is applied to the non-linear terms following Pulliam and Chausee [Reference Pulliam and Chaussee13]. Approximate factorisation scheme is applied to the resulting matrices, which factorises the operator itself resulting in efficient matrix equations with narrow bandwidth. This results in block tridiagonal matrices, which are easy to solve. The spatial derivatives can thus be approximated using second order central differences. Explicit and implicit artificial dissipation terms are added to achieve nonlinear stability. A spatially variable time step is used to accelerate convergence to steady state solutions. The algorithm can be used to solve time accurate or steady state flow problems.

The fact that the Navier Stokes equations are formulated in a strong conservation form as well as the implicit treatment of the artificial dissipation and the inclusion of a second order dissipation terms helps in the shock capturing at high Mach numbers.

We selected k- $\varepsilon $ and the Wilcox k- $\omega $ turbulence models for this investigation. For the two equations turbulence models, the convection and diffusion terms of their transport equations are negligible in the inertial sublayer so that local equilibrium prevails, which implies that the production of the turbulent kinetic energy k is equal to the dissipation rate. The local equilibrium condition leads to two simple relations that can be used as boundary conditions for k and the dissipation terms for both incompressible and compressible flows. The compressible wall functions have been successfully applied to both attached and separated flows under Mach number ranging from 0.1 to 10.

While the purpose in this study was not to scrutinise the capabilities of the turbulence models, nevertheless for these simple bodies, very little difference in results was observed from the two turbulence models. The Sutherland law is implemented to compute the laminar viscosity coefficient ( $\mu $ ). Since the flows computed involve hypersonic Mach number the choice of the Prandtl number is selected to suite the Mach number whereas the use of specific heat and the specific heat ratio c _p /c _v (γ) in the algorithm for one species gas must also respect the high Mach number flow conditions subject to the relations by Grundmann [Reference Grundmann14] and Gordon and McBride [Reference Gordon and McBride15]:

\begin{equation*}{c_v}\!\left( T \right) = R\left( {\dfrac{5}{2} + \;\dfrac{{\left( {{{\dfrac{{hv}}{{kT}}}^2}} \right){e^{{\raise0.7ex\hbox{${hv}$} \!\mathord{\big/ }\!\lower0.7ex\hbox{${kT}$}}}}}}{{{{\left( {{e^{\dfrac{{hv}}{{kT}}}} - 1} \right)}^2}}}} \right),\;\;\;\;\;\;{c_p}\!\left( T \right) - {c_v}\!\left( T \right) = R\end{equation*}

(5) \begin{equation}\frac{{{c_p}\!\left( T \right)}}{R} = \;\mathop \sum \limits_{n = 0}^{n = 4} {a_n}{T^n}\;\end{equation}

Where T is temperature, R is the real gas constant, h is Planck’s constant, v is the fundamental vibrational frequency and k is Boltzman’s constant, and where ${a_n}$ are curve-fit constants. As explained in Khalid and Juhany [Reference Khalid and Juhany16] the input flow conditions (ρ, ρu _i , E) for the outflow surface are treated like a new ‘slip’ condition, which remains fixed during the flow applied run. The algorithm iterates between the successive grid planes near the exit control boundary until the desired flow conditions are matched at the exit boundary as the convergence criterion is reached. For the corresponding solid surface run this surface would understandably carry a no-slip boundary condition. Since the present applications of the code deals with axi-symmetric flow fields produced when a simple two-dimensional mesh is rotated about the body axis through a 45° azimuthal angle. The flow has to be carefully treated near the polar stations. In such cases a computational coordinate surface may represent a physical 3D line. For example, with an “O” grid the circles (2D) or cylinders (3D) collapse to points and lines, respectively. The flow variables at the singular boundary are the average of the points surrounding the singular point.

The flow field as recovered from using above CFD method on a solid blunt cone of bluntness 0.34 and the same body with frontal blowing rate of $\dfrac{{\rho {u_{\;local}}}}{{{\rho _\infty }{u_\infty }}} = 0.576\;$ are shown in Fig. 1. The outflow at the leading edge was only implemented around the front blunt nose. Following the hypersonic heat alleviation experiment [Reference Nowak6], the efflux rate here corresponds to mass rate of 2.09 kg/s. A similar study was conducted on a blunt cone of bluntness 0.12.

Figure 1.

Mach number (solid surface) and velocity vectors (with frontal blowing) past a blunt cone of bluntness 0.32.

Another blunt Gasjet configuration with a flat frontal surface of the same circular cross section as the 0.12 blunt cone was computed at a Mach number of M = 6.8. The Mach number distribution and velocity vectors for the truncated blunt frontal nose are shown in Fig. 2.

Figure 2.

Mach number distribution past a solid blunt cone of (left) and with blowing (right).

A feature of particular bearing on flow stability is the inherent unsteady flow field (Fig. 3, right), which evolves as a result of the high Mach number flow field being confronted by the outflow immediately in front of the flat front surface. As discussed by Khalid and Juhany [Reference Khalid and Juhany17], the flow seems to be numerically instable as the solution seems to converge followed by a breakdown of the flow field uniformity and deterioration of convergence. This behaviour seems to repeat near the convergence. This type of flow behaviour as discussed in reference [Reference Garicano-Mena, Lani and Deconinck18, Reference Ismail and Roe19] is due to the carbuncle instability phenomenon.

Figure 3.

Gasjet experiment by Nowak (Ref. Reference Panaras and Drikakis8) M = 6.8, with a frontal blowing rate of, $\dfrac{{\rho {u_{\;local}}}}{{{\rho _\infty }{u_\infty }}} = 0.15$ (left). Note the unsteady flow near the leading edge (right).

It was observed in both computations and the shadowgraph that the outflow at the nose causes the stagnation region to be blown outwards in a spherical bubble shape pushing against the frontal shock. Note the presence of the supersonic bubble inside of the main shock. The flow region in between the main jet interface and the nozzle shock contains the slower expanded flow in both simulation and the experiment.

The final set of results deal with a 0.32 bluntness cone with shoulder peripheral blowing just past the spherical nose. The flow field with and without the blowing from the shoulder jet of the Slatjet is shown in Fig. 4. It is clear from the viewgraph with shoulder blowing that the boundary layer is blown off from the surface with the strong peripheral jet from the shoulder.

Figure 4.

Slatjet solid model (top) and with shoulder blowing (bottom), blowing rate = $\dfrac{{\rho {u_{\;local}}}}{{{\rho _\infty }{u_\infty }}} = 0.15$ .

The most important aspect of the present study was to be able to recover the steady pressure fields from above cases which are then in turn used in the following analysis to recover dynamic stability characteristics of blunt bodies. The pressure fields for the above blunt-nosed case bodies are shown in the following Fig. 5. While grid sensitivity study and the numerical error has been discussed in Khalid and Juhany [Reference Khalid and Juhany16, Reference Khalid and Juhany17], it can be seen that the prediction especially on the straight portion of the model lies within 1–2% of the experiment.

Figure 5.

Pressure distributions on blunt-nosed bodies: R3 (a), R1 (b), Gasjet (c) and Slat-jet (d).

The pressure distributions as computed on the blunt models showed satisfactory agreement with the measurements on the straight potion away from the blunt nose. The comparison for the blunt region was difficult, as no measurements had been recorded for this region.

3.0 Analytic approach to dynamic stability from steady flow field

The analytic approach relies on being able to obtain of a reliable pressure flow field and then using the classic shock expansion approach to perturb this solution. The angular displacement of a typical streamline near the surface of the body can be evaluated in terms of the pitch rate, the body half-angle and the flow velocity. The steady pressure distribution in all cases as discussed in previous discussion was used from existing results on heat alleviation study [Reference Khalid and Juhany16, Reference Khalid and Juhany17] on blunt conical bodies in hypersonic flow. The unsteady pressure at small angles of attack on an oscillating blunt body with or without blowing may be expressed in terms of the steady and unsteady component:

(6) \begin{equation}{P_{b,s}}\!\left( {x,\alpha ,t} \right) = \;{P_{b,s}}\!\left( {x,\alpha } \right) + \;p_{b,s}^{\prime}\!\left( {x,\alpha } \right)\end{equation}

Where, the subscripts b and s refer to either solid or frontal blowing condition and $p_{b,s}^{\prime}$ refer to the instantaneous perturbation value as results of the body oscillations.

The unsteady pressure perturbation $p_{b,s}^{\prime}$ may be evaluated from the shock expansion expression as developed by Eggers [Reference Eggers20]. The shock expansion as used here will be drawn from Refs [Reference Khalid and East21, Reference Khalid and East22] without entering too much detail. Hence:

(7) \begin{equation}p_{b,s}^{\prime}\!\left( {x,\alpha } \right) = \;\frac{{\rho V{{\left( x \right)}^2}}}{{\sqrt {(M_x^2} - 1)}}\Delta {\varepsilon _p}\end{equation}

which gives,

(8) \begin{equation}\frac{{{P_{b,s}}\!\left( {x,\alpha ,t} \right)}}{{{P_{b,s}}\!\left( {x,\alpha } \right)}} = 1 + \;\;\frac{{\gamma M_x^2}}{{\sqrt {(M_x^2} - 1)}}\Delta {\varepsilon _p}\!\left( x \right)\end{equation}

Where $\Delta {\varepsilon _p}$ is the angular change of a particle from the body vertex to a station x along the body subject to a motion schedule involving plunging and or pitching motion. This angular change may be written in terms of the instantaneous height of the particle h(x,t) and the expansion of the streamline along which it traverses. Thus following Scott [Reference Scott23],

(9) \begin{equation}\Delta {\varepsilon _{p\;}}\left( x \right) = \;\left| {ta{n^{ - 1\;}}} \right.\left\{ {\frac{{dh\!\left( {x,\tau\! \left( x \right)} \right)}}{{dx}}.\;sin\varphi + \;\frac{{dR\!\left( x \right)}}{{dx}}\;\;} \right\}\end{equation}

This inclination angle $\Delta {\varepsilon _{p\;}}\left( x \right)$ of a path line projected on the surface of the body in purely pitching motion can be shown [Reference Khalid and Juhany16] as:

\begin{equation*}\Delta {\varepsilon _{p\;}} = \;\left( {{\theta _b}\!\left( x \right) - \;{\theta _0}} \right) + 2\;\vartheta \left( {{\theta _b}\!\left( x \right) - \;{\theta _0}} \right)\sin {\theta _0}\cos {\theta _0}\sin \varphi - \;\frac{q}{{{{\bar V}_\infty }}}\cos {\theta _0}\sin \varphi .\end{equation*}

(10) \begin{align}&\Bigg\{ 2x + \;\tan {\theta _{0\;}}.\left[ {\mathop \smallint \nolimits_0^x \left( {{\theta _b}\!\left( x \right) - \;{\theta _0}} \right)dx - \left( {x - \;{x_0}} \right).\left( {{\theta _b}\!\left( x \right) - {\theta _0}} \right)} \right] \nonumber\\ & + \frac{1}{{\sqrt {M_\infty ^2 - 1} }}.\;\;\left[ {\mathop \smallint \nolimits_0^x \left( {{\theta _b}\!\left( x \right) - \;{\theta _0}} \right)dx - \left( {x - \;{x_0}} \right).\left( {{\theta _b}\!\left( x \right) - {\theta _0}} \right)} \right] \Bigg\}\end{align}

where, ${\theta _0}$ is the angle at the vertex and ${\theta _b}\!\left( x \right)$ is the inclination of the surface at the stations x. $\;q$ is the body pitch rate, is the magnitude of the free stream velocity at the nose and ${M_\infty }$ is the Mach number, which corresponds to the free stream velocity at the nose. For a conical body with constant surface slope ( ${\theta _b}\!\left( x \right) = {\theta _0}),$ the above expression simplifies to:

(11) \begin{equation}\Delta {\varepsilon _{p\;}} = \; - 2x\frac{{\cos {\theta _0}\sin \varphi }}{{{V_\infty }}}\end{equation}

and the effective angle-of-attack of the body pitching at q and instantaneous pitch angle $\vartheta $ is given by:

(12) \begin{equation}{\alpha _e} = \;\vartheta - \frac{{q\!\left( {x\sec {\theta _0} + {\theta _0}} \right)}}{{{V_\infty }}}\end{equation}

The unsteady pressure can thus be written in terms of the instantaneous pressure perturbation and the effects from pitch rate and the effective angle-of-attack.

(13) \begin{equation}\frac{{{P_{b,s}}\!\left( {x,\alpha ,t} \right)}}{{{P_\infty }}} = \frac{{{P_{b,s}}\!\left( {x,\alpha ,t} \right)}}{{{P_{b,s}}\!\left( {x,\alpha } \right)}}\frac{{{P_{b,s}}\!\left( {x,\alpha } \right)}}{{{P_{x,\alpha = 0}}}}\;\frac{{{P_{x,\alpha = 0}}}}{{{P_{x,ref}}}}\frac{{{P_{x,ref}}}}{{{P_\infty }}}\end{equation}

The pressure distribution of a body at angle–of-attack with respect to the pressure at zero angle-of-attack is given by [Reference Krasnov24]

(14) \begin{equation}{P_{,s}}\!\left( {x,\alpha } \right) = \;{P_{x,\alpha = 0}} - \;\eta \;{\alpha _e}\;Sin\;\varphi + \;\alpha _e^2\left( {D - Cos2\varphi } \right)\end{equation}

and where $\dfrac{{{P_{x,ref}}}}{{{P_\infty }}}$ represents the location on the spherical nose where the supersonic condition are first reached. This leads us to write the unsteady pressure coefficient as:

(15) \begin{equation}{C_{p\;}} = \;\frac{2}{{\gamma M_\infty ^2}}\left\{ {\frac{{{P_{b,s}}\!\left( {x,\alpha ,t} \right)}}{{{P_\infty }}}\frac{{{P_{b,s}}\!\left( {x,\alpha ,t} \right)}}{{{P_{b,s}}\!\left( {x,\alpha } \right)}}\frac{{{P_{b,s}}\!\left( {x,\alpha } \right)}}{{{P_{x,\alpha = 0}}}}\;\frac{{{P_{x,\alpha = 0}}}}{{{P_{x,ref}}}}\frac{{{P_{x,ref}}}}{{{P_\infty }}} - 1} \right\}\end{equation}

Static and pitching moment derivatives can now be obtained by obtaining $\dfrac{{\partial {C_{p\left( {b,s} \right)}}}}{{\partial \alpha }}$ and $\dfrac{{\partial {C_{p\left( {b,s} \right)}}}}{{\partial q}}$ and integrating the expression over the complete surface of the blunt body under solid and frontal blowing conditions, i.e.

(16) \begin{equation}C{m_{\alpha \;}} = \;\;\frac{1}{{S\;l}}\mathop \smallint \nolimits_0^l \mathop \smallint \nolimits_0^{2\pi } \frac{{\partial Cp\!\left( {b,s} \right)}}{{\partial \alpha }}\;R\!\left( x \right)\;l\!\left( x \right)\sin \varphi d\varphi \;dx\;\end{equation}

(17) \begin{equation}_{q\;} = \;\;\frac{{2{V_\infty }}}{{S\;l}}\mathop \smallint \nolimits_0^l \mathop \smallint \nolimits_0^{2\pi } \frac{{\partial Cp\!\left( {b,s} \right)}}{{\partial q}}\;R\!\left( x \right)\;l\!\left( x \right)\sin \varphi d\varphi \;dx\end{equation}

$l\!\left( x \right)$ is the moment ars $l = x\;Se{c^2}{\theta _0} - \;{x_0}$ and $R\!\left( x \right)$ is the local radius at station x. Pressure distribution for the blunt body under solid and blowing condition were successively introduced in the expression for the static and dynamic stability derivatives to obtain the stability characteristics. The above expression for $C{m_{q\;}}$ assumes that the $\dot \alpha $ effects in the expression $C{m_{\dot \theta }}$ are very small and, in fact, $C{m_{\dot \theta }}$ = $C{m_{q\;}}$ .

5.0 Conclusion

The present analytic modeling shows that the static stability derivative is not too adversely by blowing as the expression for $C{m_{\alpha \;}}$ is driven largely by the steady pressure distribution and not by the pitching conditions q or $\;\dot \vartheta $ . The dynamic stability appears to depreciate with both bluntness and blowing rate. Cones of smaller bluntness cannot be equipped to force excessive outflow from small blunt regions and thus their stability is not affected by blowing.

Dynamic stability of cones with sharp square-edged frontal bluntness equipped with blowing suffer more than the cones with tangential shoulder blowing. This is because the sharp edge in front traps a larger region of slower air bubble than simple rounded nose, which smoothly diverts the air downstream. The incoming flow is disturbed with larger bluntness and blowing leading to depreciation of $\;C{m_q}$ .