## 1. Introduction

The characteristics of the unsteady flowfield over a hemisphere-cylinder model and a bulbous heat shield of a satellite launch vehicle at turbulent flow and at transonic speeds and a laminar flow over a conical spike attached to a forward facing blunt body at supersonic speed are numerically simulated by solving time-dependent compressible axisymmetric Navier-Stokes equations.

Hsieh [1] has conducted wind-tunnel tests of a hemisphere-cylinder model at zero angle of attack and freestream Mach number M_{∞} = 0.7 − 1.3 to investigate viscous-inviscid interaction. Hsieh [2] has solved full potential equation to analyze wind-tunnel results and found that the inviscid analysis is unable to predict the external flowfield satisfactory. This is due to the fact that the shock wave-turbulent boundary layer interaction causes a separated flow between M_{∞} = 0.80 − 0.90 on the hemisphere-cylinder model in a high speed wind-tunnel testing. The numerical simulations analyze the unsteady flow caused by shock wave-turbulent boundary layer at transonic Mach numbers.

A bulbous payload shroud is generally selected to accommodate an increased payload volume of the satellite in a launch vehicle. All launch vehicles require a heat shield to protect the satellite from aerodynamic loading, heating, aero-acoustic vibration, and other environmental conditions during the ascent phase of the flight and to provide aerodynamic forward surface. The wind-tunnel tests for Titan I − IV were conducted during 1955 to 1996 and summarized by Brower [3]. The estimation of flowfield characteristics around such a heat shield configuration is of great aerodynamic importance, as well as research interest. For the ascent flight, during the transonic speed range, their study is particular important because of such resulting phenomena as terminal shock wave movements, frequently coupled with substantial freestream dynamic pressure. Flow induced vibrations are important design requirement of aerospace launch vehicles. Awrejcewicz, and Krysko [4] have developed numerical simulation of a cylindrical panel within transonic ideal-gas flow stream and solved dynamics for all intervals of the frequency. These parameters directly depend on the intensity of the vorticity components of the turbulence, the strength of the shock wave, and the mechanism of their interaction, all of which are implicitly linked to the specific configuration of a bulbous heat shield of a satellite launch vehicle. The numerical flow simulation over a bulbous payload shroud at transonic Mach number range is very useful to decide the geometrical configuration for minimum buffeting load and minimum aerodynamic drag requirement. The terminal shock wave of sufficient strength to interact with the boundary layer can cause flow separation and flowfield may become unstable as observed in the high speed cinematography [5]. Therefore, it is desirable to determine the location of the terminal shock wave on the heat shield and the strength of the terminal shock wave as a function of transonic Mach numbers range. The strength of the terminal shock wave and the mechanism of their interaction are related to the specific configuration of the heat shield satellite launch vehicle. Fluctuations of pressure level in shock waves and in separated flow regions can cause flow instabilities and then leads to buffeting phenomenon [6] – [7]. The features of the transonic flowfield can be delineated through the wind-tunnel data such as schlieren photographs and oil flow patterns. It is characterized by a normal or a terminal shock wave, supersonic pocket on the cylindrical region of the heat shield, shock wave/turbulent boundary layer interaction, and a separation bubble may be caused by the shock wave/turbulent boundary layer interaction on the cylindrical section of the heat shield. The main features of transonic flowfield around a bulbous heat shield are illustrated in Fig. 1. The terminal shock waves are an essential feature of transonic flow [8]. As the freestream Mach number increases from subsonic values a shock wave system appears near the shoulder. The flow is called transonic if both subsonic (M < 1) and supersonic (M > 1) regions are present in the field.

The transonic range begins when the highest Mach number reaches unity (M = 1) on the heat shield. The general features of the flow are as present once the sonic velocity occurs at the shoulder and remains throughout the whole transonic range. There is a local supersonic region ahead of the main shoulder shock. The near normal shock wave grows and moves downstream as the freestream Mach number increases. The difficulties to analyze the flowfield are associated with the detail design and a quantitative theoretical prediction. For the former, a sufficient wind-tunnel test data is required; the latter is due to nonlinearity of the equation governing transonic flow requires Computational Fluid Dynamics approach. In the boat-tail region, a local separation results, due to sharp discontinuity in the longitudinal of the payload shroud. The regions of flow separation impose additional complexity to aerodynamic and structural design aspects [9] – [13]. The complex flowfield at the transonic speeds is also observed during the experimental investigation of the bulbous heat shield. Experimental studies [14] – [15] have been made to understand flow behavior at transonic Mach numbers. These experimental investigations were limited to the measurement of surface pressure distribution, oil flow patterns, shadowgraphs and schlieren pictures for various heat shield models at transonic Mach number range. Recently analyses of Ares launch vehicle are carried out in the transonic speed and reported in a series of papers by Pinier [16], Piatak et al. [17] and Sekula et al. [18]. The current work reveals the paramount importance of aerodynamics at transonic Mach range.

A high-speed flow over a blunt body generates a bow shock wave in front of it, which causes a rather high surface pressure, and as a result, high aerodynamic drag. The surface pressure on the blunt body can be reduced if a conical shock wave is generated by mounting a forward facing spike. The aerospike produces a region of re-circulating separated flow that shields the blunt-nosed body from the incoming flow as shown in Fig. 2. The literature review reveals that addresses the mechanism how the unstable flow is initiated and how it persists [19] – [20]. The combination of the numerical simulations with experimental investigations has found to be a powerful tool to analyze unsteady flow and first results of a renewed investigation of the aerospike problem. The aerospike has been known since the 1950’s to cause an unstable flow [21]. Wood [22] has distinguished five different types of flow regimes over spiked cones based on the semi-cone angle and flow characteristics which may correspond to the fluctuation and oscillation regimes.

Bogdonoff and Vas [23] were the first to identify the two modes of unstable axisymmetric separation, the fluctuation mode and the oscillation mode by varying flat faced and hemispherical blunt bodies. The flowfield problem associated with the blunt-nosed spike bodies can be distinguished by a conical blunt body with a total angles of the conical faces varied from 30^{0} to 180^{0} [22] or a hemisphere-cylinder body [24]. Kabelitz [25] has observed two distinct unsteady flow modes, namely, oscillation and pulsation [26] in the spike attached to the blunt-nosed (flat-faced) cylindrical body. Experimental studies have been focused on identifying the boundaries of the unsteady region. The flow just outside the separated shear layer approaching the body’s shoulder can be turned by an attached conical shock, and then the shock structure is stable because an equilibrium condition is reached between escaping and recirculating flows in the separated region. Kistler [27] was the first to make detailed fluctuating wall pressure measurements under the separated supersonic turbulent boundary layer upstream of a forward step.

For a range of spike lengths the flow can became unsteady with two modes of instability observed. The oscillation mode involves the motion of the fore-shock due to the spike tip. The pulsation mode features a large-scale motion of the bow shock associates with blunt body. Spike length to diameter (*L/D*) ratio of 0.9, half-cone angle of blunt body *70*^{0}*, M*_{∞} *= 2.21, Re*_{D} *= 0.12 × 10*^{6} for oscillation modes and *L/D* ratio of *1.0*, half-cone angle of blunt body *90*^{0}*, M*_{∞} *= 6, Re*_{D} *= 0.13 × 10*^{6} for the pulsation modes are numerically investigated by Badcock et al. [28]. Feszty et al. [29] have conducted a computational analysis of the pulsation mode using computational fluid dynamics.

Flowfield over a conical spike attached to a blunt body is analyzed to understand the periodic oscillations of flowfield. The laminar Navier-Stokes equations are solved using multi-stage Runge-Kutta time stepping method. If the turning angle of the flow is too large to be accomplished by an attached conical shock wave, a detached strong shock is generated, which pushes high-pressure flow from the reattachment zone at the body’s face into the recirculating region of the separated shear layer. This high-pressure flow that gets into the separated flow region inflates the separation bubble, and the shock structure is pumped upstream. This gives rise to self-excited shock oscillations during which the conical fore-shock wave and the accompanying shear layer oscillate laterally and their shape changes periodically from concave to convex. This type of flowfield is unsteady in nature. The separated shear layer with an inflection point in the velocity profile is inherently unstable [21], and when this hits the body at the reattachment point selective amplification of the disturbances takes place, and this would cause the surface pressure to fluctuate in the flow separation region. The point of reattachment could be shifting to and fro along the body surface because of these shock oscillations. Because of this unsteady oscillation of the separation bubble, pronounced variations in the locations of separation shock, conical shock wave ahead of the blunt cone, and the reattachment point on the blunt cone surface are observed in different models with identical freestream conditions. Panaras et al. [30] have numerically simulated unsteady flows at high speeds around spiked-blunt bodies. The experimental studies are also carried out to know the effect of variations to the spike diameter to blunt body diameter ratio.

The main aim of the present Chapter to analyze the unsteady flow characteristics and wall pressure fluctuations and oscillations over the hemisphere-cylinder, the bulbous payload shroud of a typical satellite launch vehicle and the conical spike attached to the forward facing blunt body. The numerical simulations present glimpse of the instantaneous flowfield features over various models at high speeds.

## 2. Governing fluid dynamics equations

The Navier-Stokes equations describe the motion of a viscous, heat conducting compressible fluid. In the Cartesian tensor notation, let *x*_{j}, be the coordinates, *p*, *ρ*, *T* and *E* the pressure, density, temperature, and total energy, and *u*, the velocity components. The governing fluid dynamics equations can be written as

where

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

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

The heat flux component is

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

in which γ is the ratio of specific heats.

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 average and its mean square of *u* is defined as

where the integration interval *T* is chosen to be large than any significant period of the fluctuation, *u*^{’}. The integrals in the above equations are independent of starting time *t*_{0}. The fluctuations are said to be statistically stationary. The root-mean-square value of *u*^{’} is defined as

The statistical theory needs the statistical properties of the fluctuations, such as frequency correlation. Estimates of the Reynolds stress terms must be provided by a turbulence model. The simplest turbulence models augment the molecular viscosity by an eddy viscosity, *μ*_{t} that approximately represents the effects of turbulent mixing, and is estimated with some characteristic length scale such as boundary layer thickness. Baldwin-Lomax [31] proposed algebraic or zero-equation turbulence for the outer law, eliminating boundary layer thickness and momentum thickness, in favor of a certain maximum function occurring in the boundary layer. A typical transonic flow patter over a bulbous heat shield of satellite launch vehicle is illustrated in Fig. 1. The effects of compressibility start to cause a radical change in the flow. This occurs when embedded pocket of supersonic flow appear, generally in the terminal shock wave.

## 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 configuration 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

where

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

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 Pr_{t} 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]

where *ω* is the vorticity function, *l* the normal distance to the model wall, and *D*_{1} is the Van Driest’s damping factor

in the outer region

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

l

_{max}F_{max}0.25l

_{max}[(u^{2}+v^{2})]^{0.5}/F_{max}

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

The effective viscosity is given by

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

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 cell-centered 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

Where *A**i,j* is the area of the computational cell, *Q(U**i,*_{j)} and *V(U**i, j)* are inviscid and viscous fluxes respectively, and *D(U*_{i, j)} is the artificial dissipation flux added for numerical stability.

### 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 *U*_{i,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

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

with

The adaptive coefficients

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

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 simulations. The convective terms are nonlinear, hyperbolic and grid dependent. A structured non-overlapping 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 semi-discrete 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)

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]

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, non-reflecting 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*

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*^{0} and the total length of the shroud from the stagnation point to the boat tail is *0.083* m. The boat tail angle is 15^{0} 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.

*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 10^{0} 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.

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

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.

## 6. Results and discussion

### 6.1. Hemisphere-cylinder model

A terminal shock wave of sufficient strength interacting with a boundary layer can cause flow separation and the process can become unsteady [40]. The numerical procedure described in the previous section is applied to compute the flowfield over the hemisphere-cylinder at *M*_{∞}* = 0.90* and Reynolds number *5.1×10*^{6}. Figure 7 depicts the close-up view of velocity field. The strong attached flow near the hemisphere-cylinder junction and an expansion due to the hemisphere geometry makes way to a separation following a terminal shock wave accompanying with supersonic pocket. The separation and reattachment points are indicted by symbols *S* and *R*, respectively. The separation is confined to a short distance and flow reattach at *x/D* = 1.07 for M_{∞} = 0.90. Figure 7 also shows the comparison between the numerical and the experimental results of Hsieh [1] – [2]. All the essential flowfield features of the transonic flow, such as supersonic pocket, terminal shock wave and expansion region are well captured and compare well with the shadowgraph picture.

Once the initial phase of the computation is over 16 axial location (*x/D* = *0.16 - 2.50*) along the hemisphere-cylinder measured from the stagnation point of the hemisphere are selected to study the sensitivity of the unsteadiness in the flow. Figure 8 depicts instantaneous vector plot at *M*_{∞}* = 0.90*. The calculated surface pressure data are analyzed for the time mean, root mean pressure and, fast Fourier transform FFT [41] – [42].

### 6.2. Spectral analysis

A spectral analysis is carried out on the computed pressure-time data for all possible modes of fluctuations employing fast Fourier transform [41], which converts the pressure history from time domain into frequency domain. Figure 9 shows the spectrum of sound pressure level SPL over the hemisphere-cylinder model. The pressure values have been converted from Pascal to decibel (dB) of surface pressure levels.

The surface pressure levels are computed in terms of the pressure reference at *20×10*^{-6} Pa. The frequencies for which assessment were carried out at the multiples of the fundamental frequency of 26 Hz for the hemisphere-cylinder model. A high value of SPL is found at 200 Hz at *x/D* = *1.7996*, *2.0660* and *2.2630*. A significant SPL of *168 dB* at about *26* Hz is found at *x/D = 2.611* which may be attributed to the separated flow associated with the terminal shock. The SPL value increases gradually in the local supersonic pocket.

The function is non-periodic, a period *T* is to be assumed. This is dictated by the lowest or fundamental frequency that is to be considered in the analysis

where *ω* is angular frequency. The period is divided into N equal intervals of *Δt* and the function is sampled at time *t*_{m}* = mΔt.* The Fast Fourier Transform FFT is a computer algorithm for calculating Discrete Fourier Transform DFT. The FFT computes in the frequency domain can be written as

where *N* is the number of data points, *Δt* is sampling interval and *T = NΔt* is record length. The FFT offers an enormous reduction of computer time as compared to DFT.

### 6.3. Bulbous heat shield of a satellite launch vehicle

Figure 10 depicts the density contour plots for the freestream transonic Mach number of 0.80 and 0.90 for the bulbous heat shield and compare the density plots with schlieren pictures. The strength of the terminal shock wave initially increases with Mach number and later decreases. It can be observed from the figures that all of the essential flowfield features of the transonic flow, such as supersonic pocket, normal shock, and expansion and compression regions are very well captured and compare well with the schlieren pictures. The density contour plots reveal that the supersonic pocket increases with increasing freestream Mach number, and as a result, the terminal shock moves downstream with increasing freestream Mach number. It is important to mention here that the density increases ahead of the stagnation region of the heat shield moves close to the heat shield with the increasing transonic Mach number. It also depends on the cone angle of the heat shield as observed in the density contours plots of the flowfield.

The general flowfield along the payload shroud is shown in Fig. 11 for freestream Mach number of *0.80 – 1.00*. Figure 12 depicts comparison between density contours and schlieren pictures. All of the essential flow characteristics of the transonic flow, such as the supersonic pocket, the terminal shock wave, and the expansion and compression regions, are well captured by present numerical simulations. It can be seen from the density contour plots of heat shield that the formation of the terminal shock and supersonic region over the payload shroud is a function of the geometrical parameters of the heat shield [43]. Envelope of Mach lines leads to formation of the terminal shock wave. The expansion regions of the axisymmteric supersonic flow are characterized by diverging Mach lines, whereas in the compression region they converge to form a shock wave. In the shock wave theory [8] in fact there is no “expansion shock”. The flow separation on the payload shroud is caused by the terminal shock wave. The shock-induced separated flow on the cone-cylinder is not found for the heat shield with a cone angle of *15* deg. As freestream Mach number increases, the terminal shock moves downstream and the local supersonic zone increases rapidly. The terminal shock becomes so strong that, as a result of a shock wave-boundary layer interaction, boundary layer separation occurs and it is the function of shock strength, geometrical parameter of the heat shield and freestream Mach number. The density contour plots reveals that a shear layer is formed which accommodates the recirculating flow for the transonic speeds. The downstream boundary layer is found to be thick, which is nearly the boat-tail height. It is worth to mention here that the main purpose to introduce the boat-tail is to increase payload volume of a satellite launch vehicle.

The shock wave separated boundary layer and flow separation caused by boat tail geometry of heat shield generated high and low frequency pressure fluctuations. Flow induced vibrations are important issues to be taken into account the design requirement of the satellite. Shock wave separated boundary layer and flow separation caused by boat tail geometry of the heat shield generated high and low frequency pressure fluctuations. Analyses of the time-dependent flowfield feature are essential to design the bulbous heat shield of a satellite launch vehicle. Fluctuations of pressure level in shock waves and in separation areas induce flow instabilities and then structural vibration leading to the buffeting phenomenon.

In the boat-tail region, a local flow separation occurs, due to sharp discontinuity in the longitudinal curvature. The flow reattachment length (*X*_{r}* – X*_{c}) is normalized by the boat-tail height *H*, where *X*_{c} is the location of the boat tail shoulder and *X*_{r} is the reattachment point. The prediction of reattachment point is the point where the axial component of the velocity along the downstream wall changes from negative to positive. The non-dimensional separation length (*X*_{r}* – X*_{c})/*H* is found *11.7* from the velocity vector plot as shown in Fig. 13(a). The location of the flow reattachment point is required to analyze the wall pressure fluctuations. The density contour plots of Fig. 11 reveal that the supersonic region increases with increasing freestream Mach number, and as a result, the terminal shock moves downstream with increasing freestream Mach number as shown in Fig. 13(b).

Once the initial phase of the computation was completed, some unsteadiness in the flow characteristics was observed. The sixteen locations (*x/D*) along the heat shield measured from the stagnation point are selected to study the unsteadiness of the flowfield. Figure 14 shows the instantaneous flow separation and pressure distribution at different location of the heat shield. Before the analysis of the amplification factor and sound pressure levels are initiated, a statistical approach was employed to ensure that the computed data are free from transitional phase, i.e., the pressure values are representative of the numerical results, if the computation had continued for a very long period of time.

A set of histograms of *Cp* is depicted in Fig. 15. The numerical data in the separation region and other stations exhibit a Gaussian distribution. A spectral analysis was carried out on the computed pressure data for all possible modes of oscillations using fast Fourier transform FFT of MATLAB [44]. A cyclic behaviour of the pressure coefficient is observed in the vicinity of the separation and reattachment points. From the spectrum analysis low frequency pressure fluctuations and sound pressure levels are found at freestream Mach number 0.95.

The wall pressure fluctuations may arise as an effect of unsteady pressure associated to the turbulent velocity field. The Mach variation on the flow physics is the change of the location of the intense shock wave which is originally generated on the fairing and moves towards the booster for Mach approaching unity. From the acoustic point of view, it is observed that the most critical situation correspond to *M*_{∞}* = 0.80*, where a significant unsteadiness of the shock wave is observed. These behaviors can be qualitatively seen in the schlieren picture. The visualizations suggest that a shock wave is generated in the fairing region and then it moves downstream for increasing Mach. The present chapter is focalized on the transonic range of flow conditions. The characterization of the pressure fluctuations is accomplished by statistical analysis, and flow visualizations are used to the physical interpretation of the results. The overall sound pressure level is mentioned in the density contour plots. The overall sound pressure level OSPL reaches the maximum amplitude at transonic conditions.

The numerical simulation is used to analysis the unsteady flowfield characteristics of the bulbous payload shroud.

### 6.4. Statistics analysis

A statistical approach was used in order to ensure that the data are free from transitional phase, i.e., the pressure values are representive of the data, if computational had continued for a long time. The computed surface pressure data along the shroud were analyzed for the time mean and standard deviation values using the following relations:

The total time period* nΔt* (*Δt = 0.8 × 10*^{-6}*s*) and is considered after 47 900 time step computation is over. Higher order moments of pressure fluctuations, the skewness coefficient, and the kurtosis coefficient are expressed as

The skewness coefficient expresses the asymmetric of fluctuations around the mean value and the kurtosis coefficient expresses the symmetric property around the mean value. Figure 16 shows variation of mean pressure, rms, skewness, and Kurtosis coefficient.

### 6.5. Flowfield over the spiked-blunt body

The flow is assumed to be laminar for the spiked-blunt body, which is consistent with the experimental study of Crawford [24] and Kenworthy [45] and the numerical simulation of Yamauchi et al. [46], Hankey and Shang [47] and Badcock et al. [28]. Therefore, the turbulent viscosity* µ*_{t} is neglected in Eq. (10). The time-dependent axisymmetric compressible laminar Navier-Stokes equations are employed to analyze unsteady low over the spiked-blunt body. The coefficient of molecular viscosity *µ* is calculated using Sutherland’s law.

Figure 17 shows the enlarged view of the computed density contour and velocity plots for semi-cone angle of spike* α* = 10, 15, 20 and 30 deg at M_{∞} = 6.0 and *L/D* = 0.5. Figure 18 shows the pressure variation (*p/pa*) along the surface of the spiked-blunt body for different semi-cone angle of the spike. The wall pressure is normalized by freestream pressure *pa*. The *x* = *0* location is the spike/nose-tip junction. The location of the maximum pressure on the surface of the spiked-blunt body is at a body angle of about 40 deg for all the semi-cone angle of the conical spike. This location corresponds to the reattachment point. A wavy pressure distribution is observed on the spike, which may be attributed the separated flowfield behavior. The maximum pressure level is occurred at the same location on the blunt body. The flowfield can be studied using the shock polar diagram in conjunction with the Computational Fluid Dynamics approach [48]. The computed conical shock wave angles are compared with Ref. [49] and found good agreement between them.

It can be observed from the figure that interaction between the conical oblique shock wave emanating from the tip of the spike and the reattachment shock wave on the blunt body is seen. The reflected reattachment shock wave and shear layer from the interaction are seen behind the reattachment shock wave. A large separated region is observed in front of the blunt body. Flow patterns are same as that for *L/D = 0.5*. When the spike is long, the angle of the conical shock wave emanating from the spike-tip decreases and flow separation occurs slightly downstream. Since the reattachment point moved backward and the spike is long, the length of the separated region extended.

### 6.6. Flow characteristics for the spiked-blunt body

Figure 19 show the enlarged view of the density contours and velocity vector plots for spike lengths of *L/D = 0.5* and the semi-cone angle of the spike 10 deg (Fig. 5) at *M*_{∞} = 2.01, 4.15 and 6.80. It can be seen from the figures that interaction between the conical shock emanating from the tip of the spike and the reattachment shock wave on the blunt body is observed. The reflected reattachment shock wave and shear layer from the interaction are observed behind the reattachment shock wave. A large separated flow region is visible in the vector plots. In the separation region, a number of vortices exist and the velocity magnitude is very low.

Flowfield was analyzed for *M*_{∞}* = 2.01, 4.15, 6.80* for *L/D = 0.5* and for the *Re* = *0.14 × 10*^{6} based on the diameter of the hemispherical body D. Computed density contour plots with schlieren pictures are compared in Fig. 20 for L/D = 0.5, 1.0, 2.0 for M∞ = 6.80. The computed flowfield shows agreement with the schlieren photographs taken in the experiment by Yamauchi et al. [46] and Crawford [24]. Figure 21 reveals the effects of ratio of *L/D* and M_{∞} over the flow over the spiked-blunt body.

Once the oscillatory motion is established in the flow, as can be visualized in the instantaneous velocity vector plots in Fig. 22, the periodic phenomenon is investigated by a spectral analysis to obtain information on the frequency and amplitude for various modes of oscillations. The time steps *Δt* were 5.0 × 10^{-7}*s* taken for M_{∞} = 6.8 and *L/D = 0.5*.

The periodic phenomenon is investigated by a spectral analysis to obtain information on the frequency and amplitude for various modes of oscillation. Figure 23 show the pressure coefficient [*Cp = 2{(p/p*_{∞}*) – 1}/γM*_{∞}^{2}] variation with respect to time on the spiked-blunt body at *M*_{∞}* = 4.15* and *L/D = 0.5*. The interaction of the conical shock wave emanating from the spike tip with the separation vortices governs the pressure oscillations. The pressure oscillations are found more cyclic in nature and function of the *L/D* ratio and *M*_{∞}. These pressure oscillations are low in amplitude and depend on the location on the spike. This unsteady behavior of the flowfield is caused by the separated region enclosed inside the reattachment.

A spectral analysis is carried out on the computed pressure history employing FFT of MATLAB [44]. The pressure amplitude versus frequency and phase plots for *L/D = 0.5* and *M*_{∞}* = 6.8* are computed using pressure oscillations data *p(t).* In the spectrum plots, there are pressure amplitude peaks of dominant frequency and multiples of the dominant frequency at various locations of the spike. The spectral analysis of the pressure reveals that the discrete frequencies of higher modes of oscillations are multiples of the principal modes.

Figure 24 represents the pressure amplitude versus frequency and phase plots for L/D = 0.5 and M_{∞} = 6.80. In the spectrum plots, there are pr4essure amplitude peaks of dominant frequency and multiples of the dominant frequency at different stations of the spike. The spectral analysis of the pressure reveals that the discrete frequencies of higher mode of oscillation are multiples of the principal modes. The vortex pattern inside the separated region is different for different spike lengths. Therefore the second and third mode frequencies are different for different location.

### 6.7. Self-excited oscillation for the spiked-blunt body

The fluid dynamics of the self-sustained oscillatory flow is analyzed using spring-mass analogy as well as the nonlinear oscillatory model. The self-excited oscillation is governed or autonomous and draws its energy from the external source by its own periodic motion. For small oscillations, energy is fed into the system and there is “negative damping” [50, 51]. For large flow oscillations, energy is taken from the system and therefore damped. The periodic pressure behavior is analogous with differential equation describing the self-sustained oscillation of Van der Pol equation [51]. Figure 25 shows [*d(Cp)/dt*] versus *Cp*. The phase plane portrait is analyzed to understand the characteristics of the oscillatory flow. The phase plane plots are computed using the time-dependent pressure data. The phase plots reveal the characteristics of the oscillatory pressure field. The motion tends to build up small oscillations and decrease for large oscillations, which indicates that the damping term is greater than zero, hence, after the initial transient, the motion becomes periodic, represented by a closed trajectory which is also called a limit cycle.

## 7. Conclusions

The numerical simulations are carried out over the hemisphere-cylinder model and the bulbous heat shield of a satellite launch vehicle at transonic Mach number. The time-dependent compressible axisymmetric Navier-Stokes equations are solved employing multi-stage Runge-Kutta time-stepping scheme with Baldwin-Lomax turbulence model. The pressure fluctuations are computed at different location on the model. The unsteady flowfield characteristics are analyzed using fast Fourier transform. The numerical analysis also includes the numerical flow visualization and comparison with the available experimental data.

The key concern of research into the area of protruding spikes ahead of blunt bodies is the unstable flow that has been observed to exist for particular families of blunt bodies in supersonic and hypersonic flow. The length of the spike had an impact on the frequency and mode of oscillations. In this study the spike length was the principle parameter of variation on both flat faced and hemispherical blunt bodies. To this point, the focus of numerical simulations had remained on the effects that variation to the spike length and blunt body profile had on the resulting flow. Research into spike tipped blunt bodies has typically focused on variations to two main design variables; the length of the spike relative to the diameter of the blunt body, and the geometric shape of the blunt body itself. This research has drawn conclusions about the different flow regimes and the relative spike lengths that this is observed to occur for specific flow conditions. It is the objective of this project to contribute to this understanding by analyzing the effect that variations to the spike diameter relative to the blunt body diameter have on the characteristics of the flow. The numerical analysis is extended to simulate the flow over spiked-blunt body under laminar condition. The numerical simulations captured pressure oscillations in the separation region. A limit cycle is obtained that describes the self-sustained oscillation of Van der Pol equation.

## Nomenclature

Cp = specific heat at constant pressure |

Cp = pressure coefficient |

D = diameter |

c = sound velocity |

d = booster diameter |

e = internal specific energy |

E = total specific energy, (e + 0.5(uiui)) |

F , G , H = flux vectors |

L = length of spike |

M = Mach number |

N = number of data points |

Pr = Prandtl number |

p = static pressure |

qj = heat flux components |

Re = Reynolds number |

SPL = sound pressure level |

T = temperature |

t = time |

ui = velocity components |

u, v = axial and radial velocity |

U = conservative variables in vector form |

xj = Cartesian coordinate |

x, r = axial and radial coordinate |

α = semi-cone angle |

β = angle of conical shock wave |

γ = ratio of specific heats |

δ = Kronecker delta |

μ = molecular viscosity |

ρ = density |

σij = viscous stress tensor |

τij = stress tensor |

Δ = increment |

Subscripts |

D = diameter |

rms = root-mean-square |

t = turbulent |

∞ = freestream |

Superscripts |

− = time mean |

' = turbulent fluctuation |

## Acknowledgments

The author is indebted to his parents and Vikram Sarabhai Space Centre, Trivandrum, India for their valuable encouragement, support and contributions to build the research career.