Numerical Simulation of a Dual Pulse Solid Rocket Motor Flow Field

Numerical simulations are carried out for the internal flow field of a dual pulse solid rocket motor port to understand the flow behaviour. Three dimensional Reynolds Averaged Navier Stokes equations are solved alongwith shear stress transport turbulence model using commercial code. The combustion gas is assumed as a mixture of alumina and gases and single phase flow calculations are done with the thermo chemical properties provided for the mixture. The simulation captures all the essential features of the flow field. The flow accelerates through the pulse separation device (PSD) port and high temperature and high velocity gas is seen to impinge the motor wall near the PSD port. The overall total pressure drop through motor port and through PSD is found to be moderate.

For the past few decades, solid rocket motors have been extensively developed for tactical purposes and for civilian applications. Currently, different types of solid propellants rocket motors are used for different missions. Dual pulse rocket motor is one such variant of solid rocket where more than one pulse is used in the motor burn time. The first pulse of the rocket motor launches and accelerates the aerospace vehicle and after the burn out of first pulse, when the vehicle needs high manoeuvrability and power during its end game, the second pulse is fired to accelerate the missile for a second time. The better performance and higher manoeuvrability is achieved by using the same nozzle at the vehicle’s rear end for the two pulses. The aerodynamic performance of the missile is not affected by the operation of the second pulse motor.

The dual pulse solid rocket motor design consists of two burning chambers, separated by a bulkhead, designated as pulse separation device (PSD) and a nozzle. The pulse separation device protects the propellant grain in the second pulse chamber against high temperature and pressure impact during the first pulse operation. At initiation of the second pulse, the PSD reliably opens for the gas flow and the combustion gas from the second pulse chamber passes through the empty first pulse chamber and the nozzle. This design allows the initiation of the second pulse at any time after burn out of the first pulse. The use of one central nozzle for both pulses and the avoidance of lateral nozzles help the missile to show outstanding aerodynamic stability in manoeuvres during the second pulse phase. Number of dual pulse rocket motor was designed, manufactured and successfully tested in missile flights1-6 and the utility of this technology is demonstrated.

Schematic configuration of a dual pulse rocket motor is shown in Fig. 1. The first pulse chamber is filled with a finocyl-shaped aluminized composite propellant. Its burn rate is moderate. In the second pulse chamber a star-shaped low aluminized composite propellant is cast. Its burn rate is high. Both chambers are screwed together. Between both chambers the PSD is jammed. A nozzle is attached to the rear of the first pulse chamber. Typical thrust time curve of a dual pulse rocket motor is shown in Fig. 2.

Figure 1. Schematic of a dual pulse rocket motor.

Figure 2. Typical thrust-time history of a dual pulse rocket motor.

Although, it is reported in the literature that the motion of reactive gases ejected from the regressing wall of a solid propellant rocket can be explained through Taylor–Culick equation7 and Majdalani8, et al. have adapted the theories for predicting the nozzle performances in both nozzle-adapted and nozzleless configurations, the application of CFD methods to simulate internal aerodynamics inside solid rocket motors (SRM) is advantageous. The numerical simulation of the internal flow in a motor usually involves a flow that ranges from the low subsonic in the head end of the motor to the supersonic regime of nozzle exit. Starting from internal three-dimensional inviscid flow in a Titan solid rocket motor by Johnston9, much progress has been made in the CFD simulation of SRM10-15 which includes consideration of complicated geometry, viscous simulation, two phase calculation with solid particles, chemical reaction both gas and condensed phases, combustion instability, etc

In the present work, three dimensional viscous simulations are carried out in a dual pulse solid rocket motor port to understand the flow behaviour and to generate necessary thermo dynamical parameters for the design of rocket motor liner. The flow field is also analysed to estimate the pressure drop across the PSD.

Three dimensional Reynolds averaged Navier stokes (RANS) equations are solved in a fully implicit manner using a commercial CFD software CFX-1116 which is an integrated software system capable of solving diverse and complex multidimensional fluid flow problems. It is a finite volume method and is based on a finite element approach to represent the geometry. The method retains much of the geometric flexibility of finite element methods as well as the important conservation properties of the finite volume method.

The software has four major modules:
(a) Pre-processor - sets up the boundary condition and initial field condition,
(b) Solver manager - solves the flow field based on the grid and the boundary condition, and
(c) Post-processor - visualizes and extracts the results.

The software has a provision to use different numerical upwind schemes for and it ensures global convergence of mass, momentum, energy and species. In the present study, the discretisation of the convective terms is done by second order upwind difference scheme. SST turbulence model17 is used for the turbulence closure. To find out the accuracy and the range of applications, the software has been validated extensively for various complex flows including supersonic base flows18, transverse injection in supersonic flow for missile control19 etc and very good quantitative agreement with the experimental results are obtained. Recently, the same software was applied to predict the flow field of reacting gases coming out of a regressing wall of a solid rocket motor with composite propellant20 and the computed nozzle damping coefficients match with the experimental results for different throat to port area ratio. The details of governing equations, turbulence models and the discretisation schemes are given in the following subsections.

2.1 Governing Equation

The appropriate system of equations governed the turbulent compressible gas may be written as
Continuity equation:

$\frac{\partial \rho }{\mathrm{\partial t}}+\frac{\partial }{\partial {\chi }_{\kappa }}\left(\rho {\mu }_{\kappa }\right)=0$     k = 1,2,3

Momentum equation :

$\frac{\partial }{\partial t}\left(\rho {u}_{i}\right)+\frac{\partial }{\partial {x}_{\mathit{\kappa }}}\left(\rho {u}_{i}{u}_{k}\right)+\frac{\partial \rho }{\partial {x}_{\mathrm{i}}}=\frac{\partial \left({\mathit{\tau }}_{\mathrm{ik}}\right)}{\partial {x}_{\mathrm{i}}}$     ,i,k =1,2,3

Energy equation :

$\frac{\mathit{\partial }}{\mathit{\partial }t}\mathit{\left(}\mathit{\rho }E\mathit{\right)}\mathit{+}\frac{\mathit{\partial }}{\mathit{\partial }{x}_{k}}\mathit{\left(}\mathit{\rho }{u}_{k}H\mathit{\right)}\mathit{=}\mathit{-}\frac{\mathit{\partial }}{\mathit{\partial }{x}_{k}}\mathit{\left(}{u}_{j}{\mathit{\tau }}_{\mathrm{jk}}\mathit{\right)}\mathit{+}\frac{\mathit{\partial }{q}_{k}}{\mathit{\partial }{x}_{k}}$ ,j,k = 1,2,3

where, ρ,ui,p, E and H are the density, velocity components, pressure and total energy and enthalpy respectively Turbulent shear stress is defined as

${\tau }_{ik}=\mu \left(\frac{\partial {u}_{i}}{\partial {x}_{k}}+\frac{\partial {u}_{k}}{\partial {x}_{i}}\right)$

μ= μl + μt is the total viscosity; μl , μt being the laminar and turbulent viscosity

Laminar viscosity (μl ) is calculated from Sutherland law as

${\mu }_{l}={\mu }_{ref}\text{\hspace{0.17em}}{\left(\frac{T}{{T}_{ref}}\right)}^{3/2}\left(\frac{{T}_{ref}+S}{T+S}\right)$ where, T is the temperature and μref , Tref and S are known coefficient.

In eddy viscosity models, the stress tensor is expressed as a function of turbulent viscosity (μt ). Based on dimensional analysis, few variables (k, ε, ω) are defined as given below,

Turbulent kinetic energy k,

$\mathrm{k}=\overline{{{\mathrm{u}}_{\mathrm{i}}}^{\prime }{{\mathrm{u}}_{\mathrm{i}}}^{\prime }}/2$

Turbulent dissipation rate ε,

$\mathit{ϵ}\mathit{\equiv }\mathit{\nu }\overline{\frac{\mathit{\partial }{{u}_{i}}^{\mathit{\prime }}}{\mathit{\partial }{x}_{j}}\left(\frac{\mathit{\partial }{u}_{i}}{\mathit{\partial }{u}_{j}},+,\frac{{{\mathit{\partial }u}^{\mathit{\prime }}}_{j}}{{\mathit{\partial }x}_{i}}\right)}$

Specific dissipation rate ω,
ω = ε /k
The turbulent viscosity μt is calculated as

${\mathit{\mu }}_{t}\mathit{=}{c}_{\mathit{\mu }}\frac{\mathit{\rho }{k}^{\mathit{2}}}{\mathit{ϵ}}$
The heat flux qk is calculated as ${q}_{k}\mathit{}\mathit{=}\mathit{-}\mathit{\lambda }\frac{\mathit{\partial }T}{\mathit{\partial }{x}_{k}}$, λ is the thermal conductivity

2.1.1 k-ω Turbulence Model

The turbulent viscosity21 is calculated as function of k and ω ${\mu }_{\mathrm{t}}=\mathrm{f}\left(\frac{{\rho }^{\mathrm{k}}}{\omega }\right)$

Turbulent kinetic energy (k) equation:

$\frac{\partial }{\partial \mathrm{t}}\left(\rho \mathrm{k}\right)+\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{i}}}\left(\rho \mathrm{k}{\mathrm{u}}_{\mathrm{i}}\right)=\frac{\partial }{\partial {\mathrm{x}}_{\mathrm{j}}}\left({\Gamma }_{\mathrm{k}}\frac{\partial \mathrm{k}}{\partial {\mathrm{x}}_{\mathrm{j}}}\right)+{\mathrm{G}}_{\mathrm{k}}-{\mathrm{Y}}_{\mathrm{k}}$

Specific dissipation rate (ω) equation:

$\frac{\mathit{\partial }}{\mathit{\partial }t}\mathit{\left(}\mathit{\rho }\mathit{\omega }\mathit{\right)}\mathit{+}\frac{\mathit{\partial }}{\mathit{\partial }{x}_{i}}\mathit{\left(}\mathit{\rho }\mathit{\omega }{u}_{i}\mathit{\right)}\mathit{=}\frac{\mathit{\partial }}{\mathit{\partial }{x}_{j}}\left({\mathit{\Gamma }}_{\mathit{\omega }}\frac{\mathit{\partial }\mathit{\omega }}{\mathit{\partial }{x}_{j}}\right)\mathit{+}{G}_{\mathit{\omega }}\mathit{-}{Y}_{\mathit{\omega }}$
where Gk, Yk, Гk and Gω, Yω Гω are the production, dissipation and diffusion terms of the k and ω equations respectively

2.1.2 Shear Stress Transport Turbulence Model

To retain the robust and accurate formulation of Wilcox’s k-ω model in the near wall region, and to take advantage of the freestream independence of the k-ε model in the outerpart of the boundary layer, Menter 17. blended both the models through a switching function. k-ε model was transformed into Wilcox’s k-ω formulation and was multiplied by (1–F1) and added to original k-ω model multiplied by F1. The blending function F1 will be one in the near wall region and zero away from the surface.
In the second step, the definition of eddy viscosity was modified in the following way to account for the
${\nu }_{\mathrm{t}}=\frac{{a}_{\mathit{1}}k}{\mathrm{max}\left({\mathrm{a}}_{1}\omega ;\Omega {\mathrm{F}}_{2}\right)}$

transport of the principal turbulent shear stress ($\mathit{\tau }\mathit{}\mathit{=}\mathit{\text{}}\overline{\mathit{-}\mathit{\rho }{u}^{\mathit{\prime }}{v}^{\mathit{\prime }}}$ )

2.2 Discretisation of Governing Equations

The CFX-TASCflow solver utilizes a finite volume approach, in which the conservation equations in differential form are integrated over a control volume described around a node, to obtain an integral equation. The pressure integral terms in momentum integral equation and the spatial derivative terms in the integral equations are evaluated using finite element approach. An element is described with eight neighbouring nodes. The advective term is evaluated using upwind differencing with physical advection correction. The set of discretised equations form a set of algebraic equations: $\stackrel{}{A}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\stackrel{\to }{x}\text{\hspace{0.17em}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}b$ where $\stackrel{\to }{x}$ is the solution vector. The solver uses an iterative procedure to update an approximated ${x}_{n}$ (solution of x at nth time level) by solving for an approximate correction ${x}^{\prime }$ from the equation $A\text{\hspace{0.17em}}\text{\hspace{0.17em}}\stackrel{\to }{{x}^{\prime }}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\stackrel{\to }{R}$ , where $\stackrel{\to }{R}\text{\hspace{0.17em}}=\stackrel{\to }{b}\text{\hspace{0.17em}}-\text{\hspace{0.17em}}\text{\hspace{0.17em}}A\text{\hspace{0.17em}}{\stackrel{\to }{x}}_{n}$ is the residual at nth time level. The equation $A\text{\hspace{0.17em}}\text{\hspace{0.17em}}\stackrel{\to }{{x}^{\prime }}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\stackrel{\to }{R}$ is solved approximately using an approach called Incomplete Lower Upper factorization method.

2.3 Computational Grid

Taking the advantage of the symmetry, only one-eighth sector the symmetry in geometry. The origin has been chosen at the centre of the head-end of pulse-2 motor. Two unstructured grids of 1.1 million and 1.4 million points have been made for the 1/8th geometry using ICEM-CFD22 grid generator to demonstrate the grid independence of the results (shown later). The grids are clustered near the critical sections of the geometry where large gradient of flow properties are expected like PSD region, near the steps provided at the entry empty pulse-1 motor and nozzle throat regions. Hexahedral grids are used near solid boundary to capture boundary layer effects. Typical computational grid along with with the computational boundary is shown in Fig. 3. It has two symmetry planes namely symmetry-1 (the plane passing middle of two adjacent PSD) symmetry-2 (the plane passing the middle of PSD).

Figure 3. Computational geometry along with grid in the plane of symmetry 2

CFX-11.0 software16 has been used for the numerical simulation of this problem. The inflow parameters used for the simulation is given in Table 1. Mass flow rate and total temperature inlet conditions are provided at the inlet (grain surface) of pulse-2 motor. The exit of the nozzle is been kept as a supersonic outlet. There are two symmetry boundary conditions. Adiabatic and no slip conditions have been used on the walls. The exhaust gas constitutes solid/liquid alumina also due to combustion of aluminum particles. The thermo chemical properties are given for the mixture of the alumina and gases. The flow simulation is carried out assuming a single phase flow with the thermo chemical properties provided for the mixture. As mentioned earlier, Shear Stress Transport (SST) turbulence model is used to model the turbulence. The simulations have been carried out till a convergence level of 5.0x10-5 for maximum residuals.

Table 1. Inflow parameter for the simulation

The axial distribution of area average Mach number along the motor-PSD-nozzle assembly for two grids is presented in Fig. 4 to demonstrate the grid independence of the solution. The area average Mach number remains low throughout the motor but increases to supersonic value at the nozzle end. The increase in the Mach number at X/L=0.36 is due to the presence of PSD. A close agreement between the area averaged Mach numbers with two grids clearly depicts the grid independence of the solution. Static temperature at symmetry-2 plane along the length of the combustor is shown in Fig. 5. Mach number distributions at various axial planes in pulse-2 motor, PSD, pulse-1 motor and nozzle have been shown in Fig. 6.

Figure 4. Axial distribution of average Mach number for 2 different grids.

Figure 5. Static temperature distribution in the symmetry plane 2.

Figure 6. Mach number distribution at various axial locations in the rocket motor.

Moderately high Mach number and velocity have been found in the PSD port and adjacent core region of pulse-1 motor. Flow is found to expand at the convergent and divergent nozzle attached with pulse-1 motor due to which both Mach number and velocity is very high at the exit of the nozzle. To study the flow field around PSD, Mach number and velocity distribution at various axial planes in PSD and adjacent regions of pulse-1 motor have been presented in Fig.7. The flow velocity is seen to expand through the PSD and reduces again due to increase in the area in 1st pulse motor. Velocity distribution on symmetry-1 and symmetry-2 planes depicted in Fig. 8 reveals that high velocity flow (~ 500 m/s) generated from the PSD port is hitting the walls of the various steps adjacent to the PSD before turning towards the centre of the pulse-1 motor. However, velocity on the inner wall of the pulse-1 motor adjacent to PSD regions is about 35-45 m/s. The area-average axial velocity, static pressure and static temperature are evaluated at different axial locations and their axial distributions are presented in Fig.9. Flow has been found to slightly expand in the PSD regions due to which slightly higher values of Mach number and velocity are observed in that region. However, static temperature has been found to be almost constant in pulse-2 motor, PSD and pulse-1 motor. Flow expanded towards the exit of the nozzle.

Figure 7. Blown up view of (a) Mach number and (b) velocity distribution in PSD and adjacent pulse-1 motor.

Figure 8. Velocity distribution in (a) symmetry-1 (b) symmetry-2 at PSD and pulse 1 motor

Figure 9.Axial distribution of area average (a) Axial velocity, (b) Static pressure, and (c) Static temperature.

The details of the total pressure drop of the individual components are shown in the Fig. 10. It can be observed that the pressure drop across the pulse-2 motor along with the interface with PSD is about 4.11 per cent of the total pressure of pulse-2 motor (P02m). The drops of total pressure in Pulse Separation Device (PSD) and its interface with step-1 are 2.18 per cent and 3.27 per cent, respectively. The drop in the 3-steps (1st, 2nd, and 3rd) placed in between the PSD and pulse-1 motor is about 0.6 per cent. Whereas, the total pressure drop in the empty casing of the pulse-1 motor is 0.12 per cent. This low total pressure drop can be explained in the view of low velocities inside the pulse-1 motor casing. The total pressure drop inside the nozzle is about 6.66 per cent. The overall total pressure drop from the head end of pulse-2 motor to entry of the nozzle is found to be around 10.3 per cent.

Figure 10. Area averaged total pressure distribution in the axial direction

Three dimensional viscous simulations are presented to analyse the internal flow field of a dual pulse solid rocket motor. Reynolds averaged Navier Stokes (RANS) equations are solved along with SST turbulence model using commercial code. Although, the combustion products contain solid/liquid alumina due to combustion of aluminium, single phase flow calculations are done with the thermo chemical properties provided for the mixture. The simulation captures all the essential features of the flow field. Grid independence of the simulation is demonstrated by comparing the flow variables with two different grids. The flow velocity is seen to expand through the PSD and reduces again due to increase in the area in 1st pulse motor. High temperature combustion gas with velocity 500 m/s generated from the PSD port is seen to hit the walls of the various steps adjacent to the PSD. Flow velocity reduces significantly near the inner wall of the pulse-1 motor adjacent to PSD regions due to area increase. The computed flow field is analysed to obtain the pressure drop through various section of the motor geometry. It was found that out of total pressure drop of 10.3 per cent across the whole motor length, PSD and adjacent steps contribute about 6.05 per cent. The computed thermochemical variables in the motor chamber provide valuable input for the liner design of motor casing.

1. Stadler, L.J.; Hoffmann, S.; Niedermaier, H.; Hacker, A.; Stingl, R.; Bénayon, G. Trouillot, P. & Naumann, K.W. Testing and verification of the LFK NG dual pulse motor. AIAA Paper-2006-4765, 2006.

2. Naumann, K.W.; Stadler, L.J.; Trouillot, P.; Weigand, A.; Schilling, S. & Niedermaier, H. Double pulse solid rocket motor technology at Bayern-Chemie /Protac. AIAA Paper-2006-4761, 2006.

3. Trouillot, P.; Audri, D.; Rozière, J.M..; Cohen, Y.; Mialon, L.; Niedermaier, H.; Stadler, L.J. & Naumann, K.W. Design of internal thermal insulation and structures for the LFK NG double pulse motor. AIAA Paper-2006-4763, 2006.

4. Hacker, A.; Niedermaier, H. & Naumann, K.W. The safety and delay device for the LFK NG double pulse motor. AIAA Paper-2006-4764, 2006.

5. Schilling, S.; Trouillot, P. & Weigand, A. On the development and testing of a 120 mm calibre double pulse motor (DPM). AIAA Paper-2004-3387, 2004.

6. Stadler, L.J.; Trouillot, P.; Rienäcker, C.; Niedermaier, H.; Audri, D.; Ruiz, S.; Hacker, A. & Naumann, K.W. The dual pulse motor For LFK NG. AIAA Paper-2006-4762, 2006.

7. Saad, T. & Majdalani, J. On the lagrangian optimization of wall-injected flows: from the Hart-McClure Potential to the Taylor-Culick Rotational Motion. In the Proceedings of the Royal Society, London, Series A, 2010, 466(2114), pp. 331-362.

8. Majdalani, J., Vyas, A. B. & Flandro, G.A. Higher-order mean-flow solution for a rocket motor with radially regressing walls. AIAA Journal, 2002, 40(9), 1780-788.

9. Johnston, W.A. A computational fluid dynamics analysis of the internal flow in a Titan SRMU. AIAA Paper-90-2079, 1990.

10. Bruno, C. Flow analysis of a solid propellant rocket motor with aft fins. J. Propuls. Power, 1997, 13(2), page no.

11. Shimada, T.; Masumi, S. & Sekino, N. Flow inside a solid rocket motor with relation to nozzle inlet ablation. AIAA Journal, 2007, 45(6), page no.

12. Shimada, T.; Masumi, S. & Mihoko, F. Numerical investigation of roll torque induced by solid rocket motor internal flow. J. Propuls. Power, 2009, 25(6), 1300-310.

13. Dupays, J. Two-phase unsteady flow in solid rocket motors. Aerosp. Sci. Technol., 2002, 6, 413-22.

14. Fabignon, Y.; Dupays, J.; Avalon, G.; Vuillot, F.; Lupoglazoff, N.; Casalis, G. & Michel, P. Instabilities and pressure oscillations in solid rocket motors. Aerosp. Sci. Technol., 2003, 7, 191–200

15. Lee, S.N.; Seung, W.B. & Kim, K.M. Numerical analysis of quasi-steady-combustion characteristics in a solid rocket motor. J. Propuls. Power, 2010, 26(5), 980-86.

16. ANSYS CFX, Release 11, Installation and Overview, January - 2007.

17. Menter, F.R. Performance of popular turbulence models for attached and separated adverse pressure gradient flows. AIAA Journal, 1992, 30(8), 2066-072.

18. Dharavath, M.; Sinha, P.K. & Chakraborty, D. Simulation of supersonic base flow – Effect of computational grid and turbulence model. ASME J. Aerosp. Eng., 2010, 224(3), 311-19.

19. Aswin, G. & Chakraborty, D. Numerical simulation of transverse side jet interaction with supersonic free stream. Aerosp. Sci. Technol. J., 2010, 14, 295-301.

20. Javed, A & Chakraborty, D. Prediction of solid rocket motor nozzle damping coefficient using CFD technique. In the 26th Convention of Aerospace Engineers. Hyderabad, 24-25 November 2012.

21. Wilcox, D.C. Multiscale model for turbulent flows. AIAA Journal, 1988, 26(11), 1311-320.

22. ANSYS, ICEM-CFD -11: Installation and Overview. January - 2007.

 Mr Afroz Javed did his MTech from Indian Institute of Technology (IIT) Madras, in 1996. Presently he is working in the Directorate of Computational Dynamics, Defence Research & Development Laboratory (DRDL), Hyderabad. He is pursuing his PhD from Indian Institute of Science (IISc), Bengaluru. He has published 2 journal papers and 3 conference papers. His research interests include: Computational fluid dynamics (CFD), propulsion, combustion instability, and turbulent compressible mixing. Dr P Manna obtained his PhD (Thermal Science and Engineering) from IIT, Kharagpur. Presently he is working as a Scientist in the Directorate of Computational Dynamics, DRDL, Hyderabad. He has published 10 journal and 20 conference publications in the national and international level. His research interests include: CFD, propulsion, heat transfer, and high speed reacting flow. Dr Debasis Chakraborty obtained his PhD (Aerospace Engineering) from IISc, Bengaluru. Presently, he is working as Technology Director, Computational Dynamics, DRDL, Hyderabad. He was the guest editor of special Issue of Defence Science Journal on computational fluid dynamics. He has about 46 journal and 58 conference publications to his credit. His research interests include: CFD, aero-dynamics, high-speed combustion, and propulsion.