**3. Axisymmetric fluid dynamics equations**

The one of the serious problems in transonic regime of the flight of a bulbous payload shroud is wall pressure fluctuations caused by shock wave-turbulent boundary layer interaction. A terminal shock wave of sufficient strength interacting with a boundary layer may cause flow separation and boundary layer may become unstable. The strength of the terminal shock and the mechanism of its interaction with the boundary layer are linked to a specific configuration of heat shield of a satellite launch vehicle. The shock wave turbulent boundary layer interaction unsteadiness may produce large amplitude fluctuations of the loads acting on the heat shield. The frequency band of the acoustic loads is typically in the range of several hundred Hz to several kHz. The experimental results obtained from the wind-tunnel at zero angle of incidence depict that the flow pattern remains the identical with reference to the wind-tunnel configu‐ ration even when the model is rotated. The measurements are made at two diametrically opposite locations indicate that the flow is axisymmetric. Therefore, a numerical solution of the time-dependent, compressible, turbulent, axisymmetric Reynolds-averaged Navier-Stokes equation is attempted to analyze the flow at transonic speeds over the hemisphere-cylinder and the bulbous heat shield of a typical launch vehicle. Now, Equation (1) can be written as

$$\frac{\partial \mathbf{U}}{\partial t} + \frac{\partial \mathbf{F}}{\partial \mathbf{x}} + \frac{1}{r} \frac{\partial \left( r \mathbf{G} \right)}{\partial r} + \frac{\mathbf{H}}{r} = \mathbf{0} \tag{10}$$

where

where *τij*

88 Computational and Numerical Simulations

dilatation

*Sij* <sup>=</sup> <sup>1</sup>

<sup>2</sup> ( <sup>∂</sup>*ui* ∂ *xj* + ∂*uj* ∂ *xi*

by equation of state

fluctuation, *u*'

The heat flux component is

is the stress tensor, which is proportional to the rate of strain tensor and the bulk

*j i k*

æ ö ¶ ¶ ¶ æ ö

ld

= ++ ç ÷ ç ÷ ç ÷ ç ÷ ¶¶ ¶ è ø è ø (4)

¶ = - ¶ (6)

è ø (7)

= = ò ò (8)

(5)

*ji k u u u xx x*

*ij ij*

*j*

where *k* is the coefficient of heat conduction. The pressure is related to the density and energy

Turbulent flows can be simulated by the Reynolds equations, in which statistical average are taken of rapidly fluctuating Reynolds stress terms which cannot be determined from the mean values of the velocity and density. Represent *u(t)* at a particular location *(x, y, z)*. Then the time

> 0 0 0 0 1 1 '2 '2 , *t T t T t t u udt u u dt T T*

+ +

where the integration interval *T* is chosen to be large than any significant period of the

. The integrals in the above equations are independent of starting time *t*0. The

ç ÷ = +W é ù ç ÷ ë û ¶è ø

*ij ij*

t m

) , *Ωij* <sup>=</sup> <sup>1</sup>

in which γ is the ratio of specific heats.

average and its mean square of *u* is defined as

<sup>2</sup> ( <sup>∂</sup>*ui* ∂ *xj* − ∂*uj* ∂ *xi* )

Usually *λ = -2μ/3*, and velocity gradient tensor can be represented as

*i*

*j <sup>u</sup> <sup>S</sup> x* æ ö ¶

*j*

*<sup>T</sup> q k x*

( ) <sup>1</sup> <sup>1</sup> <sup>2</sup> *i i p E uu*

æ ö =- - ç ÷

 r

g

where strain rate tensor *S*ij and rotation rate tensor *Ω*ij can be written as

$$\mathbf{U} = \begin{bmatrix} \rho \\ \rho u \\ \rho v \\ \rho E \end{bmatrix} \mathbf{F} = \begin{bmatrix} \rho u \\ \rho u^2 + p - \sigma\_{\infty} \\ \rho uv - \tau\_{\infty} \\ (\rho E + p)u - u\sigma\_{\infty} - v\tau\_{\infty} + q\_{\text{x}} \end{bmatrix} \mathbf{G} = \begin{bmatrix} \rho v \\ \rho uv - \tau\_{\text{x}r} \\ \rho v^2 + p - \sigma\_{\text{rr}} \\ (\rho E + p)v - u\tau\_{\text{xx}} - v\sigma\_{\text{rr}} + q\_{\text{r}} \end{bmatrix} \mathbf{H} = \begin{bmatrix} 0 \\ 0 \\ \sigma\_{\text{\ast}} \\ 0 \end{bmatrix}$$

Reynolds stresses and turbulent heat fluxes in the mean flow equations are modeled by introducing an isotropic eddy viscosity, *μ*<sup>t</sup> . Thus the viscous terms in the above equations become

in the outer region

**a.** lmaxFmax

[13, 32, 33].

**4. Numerical scheme**

**4.1. Spatial discretization**

the integral form over a finite volume as

**b.** 0.25lmax [(u2

+v2

The effective viscosity is given by

)]0.5/Fmax

( ) <sup>0</sup> <sup>0</sup> 0.0168(1.6) *w KIF*

The scale length *l*max is the maximum value of *l* when the function *F( = l/D*1*|ω|)* attains its

max

*l*


( ) *t i* min , <sup>0</sup>

 mm

This algebraic model, which utilizes the vorticity distribution to determine the scale lengths, has been extensively used in conjunction with the Reynolds-averaged Navier-Stokes equations

To facilitate the spatial discretization in the numerical scheme, Equation (10) can be written in

where Ω is the computational domain. Γ is the boundary domain. The contour integration

Figure 3 depicts a typical stencil of the computing cell which has four edges *(1− 4)*, four vortices *(A − D)* and a cell-centre grid point *P*. The spatial and temporal terms are decoupled using the

U FG H *d dr dx d* ( ) *<sup>t</sup>* WG W ¶ W+ - = W

around the boundary of the cell is taken in the anticlockwise sense.

<sup>1</sup> <sup>6</sup>

The coefficient of *F*w is calculated as the minimum of the following two values

maximum *F*max. The Klebanoff intermittency correction factor is given by

m

1 5.5 0.3 *KIF <sup>l</sup> <sup>F</sup>*

= *F F* (13)

Unsteady Flowfield Characteristics Over Blunt Bodies at High Speed

http://dx.doi.org/10.5772/57050

= (15)

¶ òò ò (16)

(14)

91

m

$$\begin{aligned} \sigma\_{xx} &= -\frac{2}{3} (\mu + \mu\_i) \nabla \cdot \phi + 2 (\mu + \mu\_i) \frac{\partial u}{\partial x} \\\\ \sigma\_{rr} &= -\frac{2}{3} (\mu + \mu\_i) \nabla \cdot \phi + 2 (\mu + \mu\_i) \frac{\partial v}{\partial x} \\\\ \tau\_{rx} &= \tau\_{xx} = (\mu + \mu\_i) \Big( \frac{\partial u}{\partial r} + \frac{\partial v}{\partial x} \Big) \\\\ \sigma\_{\ast} &= -p - \frac{2}{3} (\mu + \mu\_i) \nabla \cdot \phi + 2 (\mu + \mu\_i) \frac{v}{r} \\\\ \nabla \cdot \phi &= \frac{\partial u}{\partial x} + \frac{\partial v}{\partial r} + \frac{v}{r} \\\\ q\_x &= -cp \Big( \frac{\mu}{\text{Pr}} + \frac{\mu\_t}{\text{Pr}\_t} \Big) \frac{\partial T}{\partial x} \\\\ q\_r &= -cp \Big( \frac{\mu}{\text{Pr}} + \frac{\mu\_t}{\text{Pr}\_t} \Big) \frac{\partial T}{\partial r} \end{aligned}$$

Temperature is related to pressure and density by the perfect gas Eq. (7). The coefficient of molecular viscosity is calculated according to Sutherland's law. The value of the turbulent Prandtl number Prt is assumed to take a constant value of 0.90. The closure of the system of equations is achieved by introducing following the algebraic turbulence model of the Baldwin-Lomax [31]

$$\left(\mu\_0\right)\_i = \left(0.4D\_1l\right)\rho\left|\rho\right|\tag{11}$$

where *ω* is the vorticity function, *l* the normal distance to the model wall, and *D*<sup>1</sup> is the Van Driest's damping factor

$$D\_1 = 1 - \exp\left[-\left(\frac{\rho\_w \left|\rho\right|\_w}{\mu\_w}\right)^{0.5} \frac{l}{26}\right] \tag{12}$$

in the outer region

$$\left(\mu\_0\right)\_0 = 0.0168(1.6) F\_w F\_{KIF} \tag{13}$$

The coefficient of *F*w is calculated as the minimum of the following two values

**a.** lmaxFmax

U=

become

*<sup>σ</sup>xx* <sup>=</sup> <sup>−</sup> <sup>2</sup>

*<sup>σ</sup>rr* <sup>=</sup> <sup>−</sup> <sup>2</sup>

*<sup>σ</sup>*<sup>+</sup> <sup>=</sup> <sup>−</sup> *<sup>p</sup>* <sup>−</sup> <sup>2</sup>

<sup>∇</sup> <sup>⋅</sup>*<sup>ϕ</sup>* <sup>=</sup> <sup>∂</sup>*<sup>u</sup>* ∂ *x* + ∂*v* ∂*r* + *v r*

*qx* <sup>=</sup> <sup>−</sup>*cp*( *<sup>μ</sup>*

*qr* <sup>=</sup> <sup>−</sup>*cp*( *<sup>μ</sup>*

Lomax [31]

Pr <sup>+</sup>

Pr <sup>+</sup>

Prandtl number Prt

Driest's damping factor

*μt* Pr*<sup>t</sup>* ) <sup>∂</sup>*<sup>T</sup>* ∂ *x*

*μt* Pr*<sup>t</sup>* ) <sup>∂</sup>*<sup>T</sup>* ∂*r*

*<sup>τ</sup>rx* <sup>=</sup>*τxr* =(*<sup>μ</sup>* <sup>+</sup> *μt*)( <sup>∂</sup>*<sup>u</sup>*

*ρ ρu ρv ρE*

, F=

90 Computational and Numerical Simulations

*ρu <sup>ρ</sup><sup>u</sup>* <sup>2</sup> <sup>+</sup> *<sup>p</sup>* <sup>−</sup>*σxx ρuv* −*τrx* (*ρE* + *p*)*u* −*uσxx* −*vτrx* + *qx*

introducing an isotropic eddy viscosity, *μ*<sup>t</sup>

<sup>3</sup> (*<sup>μ</sup>* <sup>+</sup> *μt*)<sup>∇</sup> <sup>⋅</sup>*<sup>ϕ</sup>* <sup>+</sup> 2(*<sup>μ</sup>* <sup>+</sup> *μt*)

<sup>3</sup> (*<sup>μ</sup>* <sup>+</sup> *μt*)<sup>∇</sup> <sup>⋅</sup>*<sup>ϕ</sup>* <sup>+</sup> 2(*<sup>μ</sup>* <sup>+</sup> *μt*)

∂*r* + ∂*v* <sup>∂</sup> *<sup>x</sup>* )

<sup>3</sup> (*<sup>μ</sup>* <sup>+</sup> *μt*)<sup>∇</sup> <sup>⋅</sup>*<sup>ϕ</sup>* <sup>+</sup> 2(*<sup>μ</sup>* <sup>+</sup> *μt*)

, G=

∂*u* ∂ *x*

∂*v* ∂*r*

> *v r*

Reynolds stresses and turbulent heat fluxes in the mean flow equations are modeled by

Temperature is related to pressure and density by the perfect gas Eq. (7). The coefficient of molecular viscosity is calculated according to Sutherland's law. The value of the turbulent

equations is achieved by introducing following the algebraic turbulence model of the Baldwin-

where *ω* is the vorticity function, *l* the normal distance to the model wall, and *D*<sup>1</sup> is the Van

*w w w*

m

<sup>1</sup> 1 exp <sup>26</sup>

é ù æ ö ê ú =- -ç ÷ ç ÷ è ø ë û

*<sup>l</sup> <sup>D</sup>* r w

r w

0.5

( ) ( ) 0 1 0.4 *<sup>i</sup>*

= *D l*

m

is assumed to take a constant value of 0.90. The closure of the system of

(11)

(12)

*ρv ρuv* −*τxr <sup>ρ</sup><sup>v</sup>* <sup>2</sup> <sup>+</sup> *<sup>p</sup>* <sup>−</sup>*σrr* (*ρE* + *p*)*v* −*uτxr* −*vσrr* + *qr*

, H=

. Thus the viscous terms in the above equations

> **b.** 0.25lmax [(u2 +v2 )]0.5/Fmax

The scale length *l*max is the maximum value of *l* when the function *F( = l/D*1*|ω|)* attains its maximum *F*max. The Klebanoff intermittency correction factor is given by

$$F\_{\rm KIF} = \left[ 1 + 5.5 \left( 0.3 \frac{l}{l\_{\rm max}} \right)^{6} \right]^{-1} \tag{14}$$

The effective viscosity is given by

$$
\mu\_t = \min\left(\mu\_i, \mu\_0\right) \tag{15}
$$

This algebraic model, which utilizes the vorticity distribution to determine the scale lengths, has been extensively used in conjunction with the Reynolds-averaged Navier-Stokes equations [13, 32, 33].

### **4. Numerical scheme**

#### **4.1. Spatial discretization**

To facilitate the spatial discretization in the numerical scheme, Equation (10) can be written in the integral form over a finite volume as

$$\frac{\partial}{\partial t} \int\_{\Omega} \mathbf{U}d\Omega + \int\_{\Gamma} \left(\mathbf{F}dr - \mathbf{G}d\mathbf{x}\right) = \int\_{\Omega} \mathbf{H}d\Omega \tag{16}$$

where Ω is the computational domain. Γ is the boundary domain. The contour integration around the boundary of the cell is taken in the anticlockwise sense.

Figure 3 depicts a typical stencil of the computing cell which has four edges *(1− 4)*, four vortices *(A − D)* and a cell-centre grid point *P*. The spatial and temporal terms are decoupled using the method of lines. The flux vector is divided into the inviscid and viscous components. A cellcentered scheme is used to store the flow variables [34] – [35]. The discretization of inviscid fluxes is performed using the cell average scheme. When the integral governing Eq. (16) is applied separately to each cell in the computational domain, we obtain a set of coupled differential equations of the form

$$A\_{\hat{i}\_{\hat{i}\_{\hat{i}}}\hat{j}} \frac{\partial \mathbf{U}\_{\hat{i}\_{\hat{i}}\hat{j}}}{\partial t} + Q\left(\mathbf{U}\_{\hat{i}\_{\hat{i}}\hat{j}}\right) - V\left(\mathbf{U}\_{\hat{i}\_{\hat{i}}\hat{j}}\right) + D\left(\mathbf{U}\_{\hat{i}\_{\hat{i}}\hat{j}}\right) + A\_{\hat{i}\_{\hat{i}}\hat{j}}\left(\mathbf{U}\_{\hat{i}\_{\hat{i}}\hat{j}}\right) = 0 \tag{17}$$

involving second differences is then 'switch-on' to damp oscillations in the vicinity of shock waves. This switching is achieved by means of a shock sensor based on the local second differences of pressure. Since the computational domain is having structured grids, the cell centers are defined by two indices *(i,j)*in these coordinate directions. The dissipation term are

( ) , <sup>D</sup>

tion terms are composed of first and second differences of dependent variables, e.g.

D -+-

*i,j t*

,

where *Δt*i,j is the local cell-centre time step. The cell-edge components of the artificial dissipa‐

<sup>2</sup> 1, , 1,

where *κ* (2) and *κ*(4) are constants, taken equal to *1/4* and *1/256* respectively. The scaling quantity *(ΔA/Δt)*i,j in Eq.(18) confirms the inclusion of the cell volume in the dependent variable of Eq. (16). The blend of second and fourth differences provides third-order background dissipation

The spatial discretization can be summarized here which is employed in numerical simula‐ tions. The convective terms are nonlinear, hyperbolic and grid dependent. A structured nonoverlapping quadrilateral cell is used in the numerical simulations. The diffusive terms are quasi-linear, elliptic, grid independent, cell centered use of dual control volume for evaluating the gradients at a given location. Thus, the discretized solution to the governing equations results in a set of volume-averaged state variables of mass, momentum, and energy which are balance with their area-averaged fluxes (inviscid and viscous) across the cell faces [34]. The

*p pp i j ij i j*

, <sup>2</sup> 1, , 1,

in smooth regions of the flow and first-order dissipation at shock waves.

*i j p pp i j ij i j*


*i j*

<sup>=</sup> <sup>D</sup> (18)

Unsteady Flowfield Characteristics Over Blunt Bodies at High Speed

http://dx.doi.org/10.5772/57050

93

(19)

*Ad d d d i j CD BC AB DA*

written in terms of differences of cell-edge values as

*dAB* =*dAB* (2) −*dAB* (4)

(Ui+1,j−Ui,j)

The adaptive coefficients

=max(0, *κ* (4)

max(*υi*+1, *υi*, *<sup>j</sup>*

−*ε*<sup>2</sup> (2) )

(Ui+2,j−3Ui+1,j + 3Ui,j−Ui-1,j)

)

are switched on or off by use of the shock wave sensor ν, with

u

with

*dAB* (2) =*ε*<sup>2</sup> (2)

*dAB* (4) =*ε*<sup>2</sup> (4)

*ε*2 (2) =*κ* (2)

*ε*2 (4)

Where *Ai,j* is the area of the computational cell, *Q(Ui,j)* and *V(Ui, j)* are inviscid and viscous fluxes respectively, and *D(Ui, j)* is the artificial dissipation flux added for numerical stability.

**Figure 3.** Stencil of the computational cell.

#### **4.2. Artificial dissipation**

In cell-centered spatial discretization schemes, such as above which is non-dissipative, therefore, artificial are added to Eq. (17). The approach of Jameson et al. [36] is adopted to construct the dissipative function *D*i,j consisting of a blend of second and fourth differences of the vector conserved variables *Ui,j* . Fourth differences are added everywhere in the flow domain where the solution is smooth, but are 'switched-off' in the region of shock waves. A term involving second differences is then 'switch-on' to damp oscillations in the vicinity of shock waves. This switching is achieved by means of a shock sensor based on the local second differences of pressure. Since the computational domain is having structured grids, the cell centers are defined by two indices *(i,j)*in these coordinate directions. The dissipation term are written in terms of differences of cell-edge values as

$$\mathbf{D}\_{\dot{i},\dot{j}} = \frac{\Delta A\_{\dot{i},\dot{j}} \left(d\_{AB} - d\_{\text{CD}} + d\_{\text{BC}} - d\_{\text{DA}}\right)}{\Delta t\_{\dot{i},\dot{j}}} \tag{18}$$

where *Δt*i,j is the local cell-centre time step. The cell-edge components of the artificial dissipa‐ tion terms are composed of first and second differences of dependent variables, e.g.

*dAB* =*dAB* (2) −*dAB* (4)

with

(17)

method of lines. The flux vector is divided into the inviscid and viscous components. A cellcentered scheme is used to store the flow variables [34] – [35]. The discretization of inviscid fluxes is performed using the cell average scheme. When the integral governing Eq. (16) is applied separately to each cell in the computational domain, we obtain a set of coupled

> ( ) ( ) ( ) ( ) <sup>U</sup> U U U U0 , *i,j A QV DA i j i,j i,j i,j i,j i,j <sup>t</sup>*

respectively, and *D(Ui, j)* is the artificial dissipation flux added for numerical stability.

+-++ =

Where *Ai,j* is the area of the computational cell, *Q(Ui,j)* and *V(Ui, j)* are inviscid and viscous fluxes

In cell-centered spatial discretization schemes, such as above which is non-dissipative, therefore, artificial are added to Eq. (17). The approach of Jameson et al. [36] is adopted to construct the dissipative function *D*i,j consisting of a blend of second and fourth differences of

where the solution is smooth, but are 'switched-off' in the region of shock waves. A term

. Fourth differences are added everywhere in the flow domain

differential equations of the form

92 Computational and Numerical Simulations

**Figure 3.** Stencil of the computational cell.

the vector conserved variables *Ui,j*

**4.2. Artificial dissipation**

¶

¶

$$\begin{aligned} d\,\,^{(2)}\_{AB} &= \varepsilon\_2^{(2)} \{ \mathbf{U}\_{\mathbf{i}+1,\mathbf{j}} - \mathbf{U}\_{\mathbf{i},\mathbf{j}} \} \\ d\,\,^{(4)}\_{AB} &= \varepsilon\_2^{(4)} \{ \mathbf{U}\_{\mathbf{i}+2,\mathbf{j}} - \mathbf{3} \mathbf{U}\_{\mathbf{i}+1,\mathbf{j}} + \mathbf{3} \mathbf{U}\_{\mathbf{i},\mathbf{j}} - \mathbf{U}\_{\mathbf{i}-1,\mathbf{j}} \} \end{aligned}$$

The adaptive coefficients

$$
\varepsilon\_2^{(2)} = \kappa^{(2)} \max(\upsilon\_{i+1}, \upsilon\_{i,j}),
$$

$$
\varepsilon\_2^{(4)} = \max\{0, \kappa^{(4)} - \varepsilon\_2^{(2)}\}
$$

are switched on or off by use of the shock wave sensor ν, with

$$\nu\_{i,j} = \left| \frac{p\_{i+1,j} - 2p\_{i,j} + p\_{i-1,j}}{p\_{i+1,j} + 2p\_{i,j} + p\_{i-1,j}} \right| \tag{19}$$

where *κ* (2) and *κ*(4) are constants, taken equal to *1/4* and *1/256* respectively. The scaling quantity *(ΔA/Δt)*i,j in Eq.(18) confirms the inclusion of the cell volume in the dependent variable of Eq. (16). The blend of second and fourth differences provides third-order background dissipation in smooth regions of the flow and first-order dissipation at shock waves.

The spatial discretization can be summarized here which is employed in numerical simula‐ tions. The convective terms are nonlinear, hyperbolic and grid dependent. A structured nonoverlapping quadrilateral cell is used in the numerical simulations. The diffusive terms are quasi-linear, elliptic, grid independent, cell centered use of dual control volume for evaluating the gradients at a given location. Thus, the discretized solution to the governing equations results in a set of volume-averaged state variables of mass, momentum, and energy which are balance with their area-averaged fluxes (inviscid and viscous) across the cell faces [34]. The finite volume code constructed in this manner reduces to a central difference scheme and is second-order accurate provided that the mesh is smooth enough [34]. The cell-centered spatial discretization scheme is non-dissipative; therefore, artificial dissipation terms are included as a blend of a Laplacian and biharmonic operator in a manner analogous to the second and fourth difference. The artificial dissipation term [36] was added explicitly to prevent numerical oscillations near the shock waves to damp high-frequency modes.

### **5. Multi-stage time-stepping scheme**

The spatial discretiztion described above reduces the governing flow equations to semidiscrete ordinary differential equations. The integration is performed employing an efficient multi-stage scheme [36]. The following three-stage, time-stepping scheme is used for the numerical simulation (for clarity, the subscripts *i* and *j* are neglected here)

$$\begin{aligned} \mathbf{U}^{(0)} &= \mathbf{U}^{\mathrm{II}} \\ \mathbf{U}^{(1)} = \mathbf{U}^{\mathrm{v}} - 0.6\Delta t \Big( \mathbf{R}^{(0)} - \mathbf{D}^{(0)} \Big) \\ \mathbf{U}^{(2)} = \mathbf{U}^{\mathrm{n}} - 0.6\Delta t \Big( \mathbf{R}^{(1)} - \mathbf{D}^{(0)} \Big) \\ \mathbf{U}^{(3)} = \mathbf{U}^{\mathrm{n}} - 1.0\Delta t \Big( \mathbf{R}^{(2)} - \mathbf{D}^{(0)} \Big) \\ \mathbf{U}^{\mathrm{n}+1} = \mathbf{U}^{(3)} \end{aligned} \tag{20}$$

**5.1. Initial and boundary conditions**

**5.2. Geometrical details of the Model**

of conservative vector, *U*.

*Hemisphere-cylinder model*

The boat tail angle is 150

**Figure 4.** Geometry of the bulbous heat shield.

The L/D ratios of the spike are 0.5, 1.0 and 2.0.

*Spike geometry*

*Heat shield geometry*

is *20*<sup>0</sup>

Conditions corresponding to a freestream Mach number were given as initial conditions. On the surface, no slip condition is considered together with an adiabatic wall condition. The symmetric conditions were applied on the centerline. For the transonic flow simulations, nonreflecting far field boundary conditions are applied at the outer boundary of the computational cell. For supersonic flow, all the flow variables are extrapolated at the outflow from the vector

The dimension of the hemisphere-cylinder model is taken as *2.54×10*-2 m and *25.4×10*-2 m length.

The maximum payload shroud diameter *D* of the model is *0.04* m and the booster diameter *d* is *0.035* m. The symbols *D* and *d* are depicted in Fig. 4. The inclination at the fore body junction

The dimensions of the spiked-blunt body are depicted in Fig. 5. The model is axisymmetric, the main body has a hemispherical-cylinder nose, and the diameter D is 7.62×10-2 m. The aerospike consists of a conical part and a cylindrical part. The angle of the spike's cone is 100 and the diameter of the spike is 0.1D. The aerospike model has a simple stick configuration.

flow direction. The values of *l*1*/D* and *l*2*/D* are 0.96 and 1.4, respectively.

and the total length of the shroud from the stagnation point to the boat tail is *0.083* m.

measured clock-wise from the axis with reference to the on-coming

Unsteady Flowfield Characteristics Over Blunt Bodies at High Speed

http://dx.doi.org/10.5772/57050

95

where *n* is the current time level, *n + 1* is the new time level, and residual *R* is the sum of the inviscid and viscous fluxes. The multi-stages time-stepping scheme has been proved to be second-order accurate in time for a linear system of one-dimensional equation [35]. The artificial dissipation is evaluated only at the first stage. The permissible time step of an explicit scheme is limited by the Courant-Friedrichs-Lewy condition, which states that a difference scheme cannot be convergent and stable approximation unless its domain of dependence contains the domain of dependence of the corresponding differential equation. A conservative choice of the Courant number is made in the simulation to achieve a stable numerical solution. A global time-step is used rather than the grid-varying time-step to numerically simulate a time-accurate solution and is computed using following expression [37]

$$\left(\Delta t\right)\_{i,j} = \min \left[ \frac{|\mu|}{\Delta x} + \frac{|v|}{\Delta r} + c \sqrt{\frac{1}{\left(\Delta x\right)^2 + \left(\Delta r\right)^2}} \right]^{-1} \tag{21}$$

where *i,j* are grid point as shown in Fig. 3.

### **5.1. Initial and boundary conditions**

Conditions corresponding to a freestream Mach number were given as initial conditions. On the surface, no slip condition is considered together with an adiabatic wall condition. The symmetric conditions were applied on the centerline. For the transonic flow simulations, nonreflecting far field boundary conditions are applied at the outer boundary of the computational cell. For supersonic flow, all the flow variables are extrapolated at the outflow from the vector of conservative vector, *U*.

### **5.2. Geometrical details of the Model**

### *Hemisphere-cylinder model*

The dimension of the hemisphere-cylinder model is taken as *2.54×10*-2 m and *25.4×10*-2 m length.

### *Heat shield geometry*

finite volume code constructed in this manner reduces to a central difference scheme and is second-order accurate provided that the mesh is smooth enough [34]. The cell-centered spatial discretization scheme is non-dissipative; therefore, artificial dissipation terms are included as a blend of a Laplacian and biharmonic operator in a manner analogous to the second and fourth difference. The artificial dissipation term [36] was added explicitly to prevent numerical

The spatial discretiztion described above reduces the governing flow equations to semidiscrete ordinary differential equations. The integration is performed employing an efficient multi-stage scheme [36]. The following three-stage, time-stepping scheme is used for the

1 (0) 0

*() n ( )*

*t*

*t*

*t n ()*

where *n* is the current time level, *n + 1* is the new time level, and residual *R* is the sum of the inviscid and viscous fluxes. The multi-stages time-stepping scheme has been proved to be second-order accurate in time for a linear system of one-dimensional equation [35]. The artificial dissipation is evaluated only at the first stage. The permissible time step of an explicit scheme is limited by the Courant-Friedrichs-Lewy condition, which states that a difference scheme cannot be convergent and stable approximation unless its domain of dependence contains the domain of dependence of the corresponding differential equation. A conservative choice of the Courant number is made in the simulation to achieve a stable numerical solution. A global time-step is used rather than the grid-varying time-step to numerically simulate a

*() n*

= = -D -

<sup>0</sup> U U U U 0.6 R D

2 n (1) 0

*( ) ( )*

= -D -

U U 0.6 R D

U U 1.0 R D 1 3 U U

time-accurate solution and is computed using following expression [37]

, 2 2 <sup>1</sup> min *i j u v t c*

ê ú D = ++

*x r x r*

D D D +D ë û


( )

where *i,j* are grid point as shown in Fig. 3.

3 (2) 0

*() n ( )*

= -D - <sup>+</sup> <sup>=</sup>

( )

( )

(20)

(21)

( )

( ) ( )

1

oscillations near the shock waves to damp high-frequency modes.

numerical simulation (for clarity, the subscripts *i* and *j* are neglected here)

**5. Multi-stage time-stepping scheme**

94 Computational and Numerical Simulations

The maximum payload shroud diameter *D* of the model is *0.04* m and the booster diameter *d* is *0.035* m. The symbols *D* and *d* are depicted in Fig. 4. The inclination at the fore body junction is *20*<sup>0</sup> and the total length of the shroud from the stagnation point to the boat tail is *0.083* m. The boat tail angle is 150 measured clock-wise from the axis with reference to the on-coming flow direction. The values of *l*1*/D* and *l*2*/D* are 0.96 and 1.4, respectively.

**Figure 4.** Geometry of the bulbous heat shield.

### *Spike geometry*

The dimensions of the spiked-blunt body are depicted in Fig. 5. The model is axisymmetric, the main body has a hemispherical-cylinder nose, and the diameter D is 7.62×10-2 m. The aerospike consists of a conical part and a cylindrical part. The angle of the spike's cone is 100 and the diameter of the spike is 0.1D. The aerospike model has a simple stick configuration. The L/D ratios of the spike are 0.5, 1.0 and 2.0.

and the implementation of boundary and initial conditions become important factor for the selection of the computational region. The known physically acceptable far-field boundary conditions usually limit the flow variables to asymptotic values at large distances from the payload shroud. For the supersonic speeds, the computational domain is kept 4 to 6 times the maximum diameter *D*. Figure 6 depicts the computational grid in the physical domain of the hemisphere-cylinder, the heat shield and the conical spiked attached to the blunt body. The

Unsteady Flowfield Characteristics Over Blunt Bodies at High Speed

http://dx.doi.org/10.5772/57050

97

grid independent is carried out using the above mentioned numerical algorithm.

(a) Hemisphere-cylinder (b) Bulbous heat shield

(c) Spiked-blunt body (d) Close-up view of grid over spiked-blunt body

**Figure 6.** Computational grids.

**Figure 5.** Dimensions of the spiked-blunt body.

### **5.3. Computational grid**

One of the controlling factors for the numerical simulation is the proper grid arrangement. The following procedure is adopted to generate grid in the computational domain of the model. The computational region is divided into a number of non-overlapping zones. The mesh points are generated in each zone using finite element methods [38] in conjunction with the homotopy scheme [39]. The above models are defined by a number of mesh points in the cylindrical coordinate system. Using these surface points as the reference node, the normal coordinate is then described by the exponentially stretched grid points extending towards up to an outer computational boundary.

$$\begin{aligned} \mathbf{x}\_{i,j} &= \mathbf{x}\_{i,o} \left[ \frac{e^{\left(j-1\right)\beta f\left(wr-1\right)} - 1}{e^{\beta} - 1} \right] + \mathbf{x}\_{i,w} \left[ 1 - \frac{e^{\left(j-1\right)\beta f\left(wr-1\right)} - 1}{e^{\beta} - 1} \right] \\\ r\_{i,j} &= r\_{i,o} \left[ \frac{e^{\left(j-1\right)\beta f\left(wr-1\right)} - 1}{e^{\beta} - 1} \right] + r\_{i,w} \left[ 1 - \frac{e^{\left(j-1\right)\beta f\left(wr-1\right)} - 1}{e^{\beta} - 1} \right] \\\ i &= 1,2,3...,nr\_{i} \qquad j = 1,2,3...,nr \end{aligned} \tag{22}$$

where subscripts *o* and *w* are wall and outer surface points, respectively, *β* is a stretching factor. *nx* and *nr* are total number of grid points in *x* and *r* directions, respectively.

The outer boundary of the computational domain is varied from 5 to 8 times the cylinder diameter, *D*, and the grid-stretching factor *β* in the radial direction varied from 1.5 to 5.0. At transonic freestream Mach number, the computational domain of dependence is unbounded, and the implementation of boundary and initial conditions become important factor for the selection of the computational region. The known physically acceptable far-field boundary conditions usually limit the flow variables to asymptotic values at large distances from the payload shroud. For the supersonic speeds, the computational domain is kept 4 to 6 times the maximum diameter *D*. Figure 6 depicts the computational grid in the physical domain of the hemisphere-cylinder, the heat shield and the conical spiked attached to the blunt body. The grid independent is carried out using the above mentioned numerical algorithm.

(c) Spiked-blunt body (d) Close-up view of grid over spiked-blunt body

**Figure 6.** Computational grids.

**Figure 5.** Dimensions of the spiked-blunt body.

One of the controlling factors for the numerical simulation is the proper grid arrangement. The following procedure is adopted to generate grid in the computational domain of the model. The computational region is divided into a number of non-overlapping zones. The mesh points are generated in each zone using finite element methods [38] in conjunction with the homotopy scheme [39]. The above models are defined by a number of mesh points in the cylindrical coordinate system. Using these surface points as the reference node, the normal coordinate is then described by the exponentially stretched grid points extending towards up to an outer

( ) ( ) ( ) ( )


*j nr j nr*

*e e*

*e e*

*nx* and *nr* are total number of grid points in *x* and *r* directions, respectively.

= =

1 1 1 1

1 1

1 1 1,2,3...., , 1,2,3....,

<sup>1</sup> <sup>1</sup> <sup>1</sup>

b

b

b

b

(22)

<sup>1</sup> <sup>1</sup> <sup>1</sup>

( ) ( ) ( ) ( )


*j nr j nr*

*e e i nx j nr*

<sup>é</sup> ù é <sup>ù</sup> - - <sup>=</sup> <sup>ê</sup> ú ê + - <sup>ú</sup> <sup>ê</sup> - ú ê - <sup>ú</sup> <sup>ë</sup> û ë <sup>û</sup>

*e e*

<sup>é</sup> ù é <sup>ù</sup> - - <sup>=</sup> <sup>ê</sup> ú ê + - <sup>ú</sup> <sup>ê</sup> - ú ê - <sup>ú</sup> <sup>ë</sup> û ë <sup>û</sup>

1 1 1 1

where subscripts *o* and *w* are wall and outer surface points, respectively, *β* is a stretching factor.

The outer boundary of the computational domain is varied from 5 to 8 times the cylinder diameter, *D*, and the grid-stretching factor *β* in the radial direction varied from 1.5 to 5.0. At transonic freestream Mach number, the computational domain of dependence is unbounded,

, , ,

*x x x*

b

*ij io i w*

b

, , ,

*r r r*

b

*ij io i w*

b

**5.3. Computational grid**

96 Computational and Numerical Simulations

computational boundary.
