Numerical Exploration of Solid Rocket Motor Blast Tube Flow Field

The blast tube flowfield of a solid rocket motor is explored numerically by solving 3-D RANS equations with SST Turbulence model using a commercial computational fluid dynamics (CFD) software CFX-10. Parametric studies are carried out to find out the effect of the blast tube diameter on the total pressure loss in the rocket motor. It is observed that the total pressure loss in the rocket motor is less than 4 per cent and the blast tube is contributing less than 1 per cent. It is also found out that higher the blast tube diameter, the lesser the drop in the total pressure. Blast tube geometry is not found to contribute significantly in the overall thrust and specific impulse in the rocket motor.

Keywords:  /b> Blast tube  total pressure loss  computational fluid dynamics

The combustion chamber of a solid rocket motor has solid propellant grain consisting of a composition that generates energy and increases chamber pressure while undergoing chemical reactions. This thermal energy is converted to mechanical energy by accelerating the products of combustion through a nozzle that generates a thrust force. The nozzles generally employed in solid rocket motors, depending upon their applications and mission requirements, can be generally classified in five categories 1 : (1) Fixed, (2) Movable, (3) Submerged, (4) Extendible, and (5) Blast tube mounted nozzles. The Blast tube mounted nozzle is generally used in tactical missiles with diameter constraints to allow space for different subsystems. It also allows the centre of gravity (CG) of rocket motor to be close to or ahead of the complete vehicle centre of gravity. This minimizes the CG travel as the motor burns, making vehicle stabilization and control much easier. However, there is no benefit of the straight duct in blast tube mounted nozzle, as per as the motor performance is concerned. There is a decrease in motor performance because of unnecessary friction losses. The designers accept this penalty in order to achieve other design objectives in the form of availability of packing space for other non-propulsion subsystems and management of CG travel.

Internal flow in the blast tube is preferably subsonic, otherwise there will be shock reflections and losses, flow changes, failure to achieve correct exhaust conditions etc. leading to more loss of thrust. In some cases where aluminized propellant is used with small throat diameter and longer operation time, alumina deposition on the throat leads to problem of increased chamber pressure in the case of subsonic blast tube. In order to avoid the alumina particles to cool and deposit on the throat, a ‘supersonic blast tube’ is utilised2. Except for such specialised applications, the blast tube internal flow field is kept subsonic. The flow through blast tube is associated with losses in total pressure leading to drop in Isp. These losses are dependent on the tube geometry viz. the diameter and length of the tube. The quantification of the thrust losses in the blast tube are of great importance to the designer. A small diameter of the blast tube allows more space for other subsystems while the flow losses are expected to be higher due to small cross sectional area. Detailed flow simulation in the blast tube is not reported adequately in the literature. A numerical study was made on the performance for different configurations of blast tubes by Sinha and Javed3. In this study, the performance of different configurations of supersonic and subsonic blast tubes was studied in terms of the thrust delivered.

Tahsini and Ebrahimi4 have carried out a parametric study for the effect of blast tube on the internal ballistics of a solid rocket motor by solving quasi one-dimensional unsteady Euler equations with a model for wall friction. In absence of quantitative data for the losses in blast tube, designers often consider some adhoc losses in the blast tube leading to the conservative motor design. In the present work high fidelity flow simulations are carried out for different blast tube geometries to quantify the losses in blast tube. Three dimensional Navier Stokes equations along with SST turbulence model have been solved in the flow field using a commercial CFD software CFX 10.05. Parametric studies are carried out with different blast tube diameters to find its effect on total pressure loss. The computed total pressure loss and specific impulse (Isp) loss for these tubes at sea level are compared and relative merits are discussed.

Four different geometries of blast tube with different diameters (70 mm, 75 mm, 80 mm, and 85 mm) are considered for this study. The inlet diameters of all the geometries are same and the divergent portions of the convergent divergent nozzles attached at the end of the tail pipe are also same. The length of the tail pipe including the nozzle is also kept constant (400 mm) only the diameter of the tubular portion is varied, and accordingly the interface region of blast tube and the nozzle also gets changed especially up to the throat of the nozzle. A typical geometry of the rocket motor with blast tube is shown in Fig.1.

Figure 1. The penetration test system.

These geometries are imported as CAD models for mesh generation. Taking the advantage of symmetry of the geometry, a 60o sector is taken for numerical simulation. Unstructured grids are generated using CFX5.66 software. The grid contains hexahedral elements near the wall in order to capture boundary layer. The grid distribution for 70 mm diameter blast tube is shown in Fig 2. It can be seen from these Figs. that grids are clustered sufficiently near the wall so that viscous losses can be properly estimated.

Figure 2.Unstructured grids in symmetry plane

Three dimensional Reynolds Averaged Navier Stokes (RANS) equations are solved using CFX-10 code5, which is an integrated software capable of solving diverse and complex multidimensional fluid flow problems. The code is fully implicit, finite volume method with finite element based discretization of geometry. The convective terms are discretized through 2nd order scheme and SST turbulence model is used. Log-normalized maximum residue of -04 is considered as the convergence criteria. The details of the governing equations, thermodynamics and the discretization schemes are given in the following subsections. To find out the accuracy and the range of applications, the software has been validated for various internal flow fields in the rectangular duct behind backward facing step7,8, base flow9, free jets10, free stream and jet interaction11-12, dual pulse rocket motor13, air intakes14, etc. and good quantitative agreement has been obtained between experimental and computational results. In order to model turbulence SST turbulence model is used. SST turbulence model is derived by blending k-ω and k-ε turbulence models through a blending function, to retain the robust and accurate formulation of Wilcox’ k-ω model in the near wall region, and to take advantage of the free stream independence of the k-ε model in the outer part of the boundary layer. This turbulence model has shown very good pressure drop predictions for both non-reacting13 and reacting15 flow simulations dealing with internal flow situations carried out by the authors.

3.1 Governing Equations

The appropriate system of equations governs the turbulent flow of a compressible gas may be written as

Continuity

$\frac{\partial \overline{\rho }}{\partial t}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\frac{\partial }{\partial {x}_{j}}\left(\overline{\rho }\text{\hspace{0.17em}}{\stackrel{˜}{u}}_{j}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}0$

Momentum (Navier – Stokes)

$\frac{\partial }{\partial t}\left(\overline{\rho }\text{\hspace{0.17em}}{\stackrel{˜}{u}}_{i}\right)\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\frac{\partial }{\partial {x}_{j}}\left(\overline{\rho }\text{\hspace{0.17em}}{\stackrel{˜}{u}}_{i}{\stackrel{˜}{u}}_{j}\text{\hspace{0.17em}}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\text{\hspace{0.17em}}-\frac{\partial \overline{P}}{\partial {x}_{j}}\text{\hspace{0.17em}}+\frac{\partial }{\partial {x}_{j}}\left(\begin{array}{l}\text{\hspace{0.17em}}{\mu }_{eff}\left(\frac{\partial {\stackrel{˜}{u}}_{i}}{\partial {x}_{j}}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\frac{\partial {\stackrel{˜}{u}}_{j}}{\partial {x}_{i}}\right)\\ -\text{\hspace{0.17em}}\frac{2}{3}{\mu }_{eff}\left(\frac{\partial {\stackrel{˜}{u}}_{l}}{\partial {x}_{l}}\right){\delta }_{ij}\end{array}\right)\text{\hspace{0.17em}}$
where
${\mu }_{eff}\text{\hspace{0.17em}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\text{\hspace{0.17em}}\mu \text{\hspace{0.17em}}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{\mu }_{t}$
${\mu }_{t}$ is a modelling constant and is known as eddy viscosity. It is further discussed in turbulence modelling.

Energy

$\begin{array}{l}\frac{\partial }{\partial t}\left(\overline{\rho }\stackrel{˜}{H}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}-\frac{\partial \overline{P}}{\partial t}+\text{\hspace{0.17em}}\frac{\partial }{\partial {x}_{j}}\left(\overline{\rho }\stackrel{˜}{H}{\stackrel{˜}{u}}_{j}\text{\hspace{0.17em}}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{\partial }{\partial {x}_{j}}\left(\text{\hspace{0.17em}}\lambda \frac{\partial \stackrel{˜}{T}}{\partial {x}_{j}}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\frac{{\mu }_{t}}{{\mathrm{Pr}}_{t}}\frac{\partial \stackrel{˜}{h}}{\partial {x}_{j}}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\\ \frac{\partial }{\partial {x}_{j}}\left({\stackrel{˜}{u}}_{i}\text{\hspace{0.17em}}\left(\text{\hspace{0.17em}}{\mu }_{eff}\left(\frac{\partial {\stackrel{˜}{u}}_{i}}{\partial {x}_{j}}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\frac{\partial {\stackrel{˜}{u}}_{j}}{\partial {x}_{i}}\right)-\text{\hspace{0.17em}}\frac{2}{3}{\mu }_{eff}\left(\frac{\partial {\stackrel{˜}{u}}_{l}}{\partial {x}_{l}}\right){\delta }_{ij}\right)\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\mu \frac{\partial k}{\partial {x}_{j}}\right)\end{array}$
Apart from these equations, equation of state is used to close the system of equations.
$\overline{P}\text{\hspace{0.17em}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\overline{\rho }\text{\hspace{0.17em}}R\left({\stackrel{˜}{Y}}_{m}\right)\stackrel{˜}{T}$

Turbulence Modelling

Turbulence is modelled using Menter’s Shear Stress Transport (SST) model16. The SST turbulence model is derived by blending k - ω model applied to the inner portion of the turbulent boundary layer with a high Reynolds number form of the k - ε turbulence model transformed into the k and ω variables being applied to the outer portion of the turbulent boundary layer. A parameter F1 is defined so as be unity for the near wall region and to vary smoothly to zero as the outer region of the turbulent boundary layer is reached. By assigning a weight of F1 to the inner k - ω model and a weight of (1-F1) to the outer transformed k - ε model, advantages of both models are incorporated into the new SST model. Additionally, the eddy viscosity relation is modified to provide a lag in development of the eddy viscosity for strong interaction flows. Both the model equations are given as follows.

Original k-ω model

$\frac{\partial \left(\overline{\rho }k\right)}{\partial t}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\text{\hspace{0.17em}}\frac{\partial \left(\overline{\rho }{\stackrel{˜}{u}}_{j}k\right)}{\partial {x}_{j}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\text{\hspace{0.17em}}{P}_{k}\text{\hspace{0.17em}}-\text{\hspace{0.17em}}{\beta }^{*}\overline{\rho }k\omega \text{\hspace{0.17em}}+\text{\hspace{0.17em}}\frac{\partial }{\partial {x}_{j}}\text{\hspace{0.17em}}\left[\left(\mu \text{\hspace{0.17em}}+\text{\hspace{0.17em}}\frac{{\mu }_{t}}{{\sigma }_{k1}}\right)\text{\hspace{0.17em}}\frac{\partial k}{\partial {x}_{j}}\right]$
$\frac{\partial \left(\overline{\rho }\omega \right)}{\partial t}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\text{\hspace{0.17em}}\frac{\partial \left(\overline{\rho }{\stackrel{˜}{u}}_{j}\omega \right)}{\partial {x}_{j}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\alpha }_{1}\frac{\omega }{k}{P}_{k}\text{\hspace{0.17em}}-\text{\hspace{0.17em}}{\beta }_{1}\overline{\rho }{\omega }^{2}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\frac{\partial }{\partial {x}_{j}}\text{\hspace{0.17em}}\left[\left(\mu \text{\hspace{0.17em}}+\text{\hspace{0.17em}}\frac{{\mu }_{t}}{{\sigma }_{\omega 1}}\right)\text{\hspace{0.17em}}\frac{\partial \omega }{\partial {x}_{j}}\right]$

Transformed k - ε model

$\frac{\partial \left(\overline{\rho }k\right)}{\partial t}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\text{\hspace{0.17em}}\frac{\partial \left(\overline{\rho }{\stackrel{˜}{u}}_{j}k\right)}{\partial {x}_{j}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\text{\hspace{0.17em}}{P}_{k}\text{\hspace{0.17em}}-\text{\hspace{0.17em}}{\beta }^{*}\overline{\rho }k\omega \text{\hspace{0.17em}}+\text{\hspace{0.17em}}\frac{\partial }{\partial {x}_{j}}\text{\hspace{0.17em}}\left[\left(\mu \text{\hspace{0.17em}}+\text{\hspace{0.17em}}\frac{{\mu }_{t}}{{\sigma }_{k2}}\right)\text{\hspace{0.17em}}\frac{\partial k}{\partial {x}_{j}}\right]$
$\begin{array}{l}\frac{\partial \left(\overline{\rho }\omega \right)}{\partial t}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\text{\hspace{0.17em}}\frac{\partial \left(\overline{\rho }{\stackrel{˜}{u}}_{j}\omega \right)}{\partial {x}_{j}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\alpha }_{2}\frac{\omega }{k}{P}_{k}\text{\hspace{0.17em}}-\text{\hspace{0.17em}}{\beta }_{2}\overline{\rho }{\omega }^{2}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\frac{\partial }{\partial {x}_{j}}\text{\hspace{0.17em}}\left[\left(\mu \text{\hspace{0.17em}}+\text{\hspace{0.17em}}\frac{{\mu }_{t}}{{\sigma }_{\omega 2}}\right)\text{\hspace{0.17em}}\frac{\partial \omega }{\partial {x}_{j}}\right]\text{\hspace{0.17em}}+\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}2\overline{\rho }{\sigma }_{\omega 2}\frac{1}{\omega }\text{\hspace{0.17em}}\frac{\partial k}{\partial {x}_{j}}\text{\hspace{0.17em}}\frac{\partial \omega }{\partial {x}_{j}}\end{array}$

Now the first two equations for k-ω model are multiplied by function F1 and the equations for k - ε model are multiplied by a function 1-F1 and the corresponding k and ω-equations are added to give the following equations.

$\frac{\partial \left(\overline{\rho }k\right)}{\partial t}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\text{\hspace{0.17em}}\frac{\partial \left(\overline{\rho }{\stackrel{˜}{u}}_{j}k\right)}{\partial {x}_{j}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\text{\hspace{0.17em}}{P}_{k}\text{\hspace{0.17em}}-\text{\hspace{0.17em}}{\beta }^{*}\overline{\rho }k\omega \text{\hspace{0.17em}}+\text{\hspace{0.17em}}\frac{\partial }{\partial {x}_{j}}\text{\hspace{0.17em}}\left[\left(\mu \text{\hspace{0.17em}}+\text{\hspace{0.17em}}\frac{{\mu }_{t}}{{\sigma }_{k3}}\right)\text{\hspace{0.17em}}\frac{\partial k}{\partial {x}_{j}}\right]$
$\begin{array}{l}\frac{\partial \left(\overline{\rho }\omega \right)}{\partial t}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\text{\hspace{0.17em}}\frac{\partial \left(\overline{\rho }{\stackrel{˜}{u}}_{j}\omega \right)}{\partial {x}_{j}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\text{\hspace{0.17em}}{\alpha }_{3}\frac{\omega }{k}{P}_{k}\text{\hspace{0.17em}}-\text{\hspace{0.17em}}{\beta }_{3}\overline{\rho }{\omega }^{2}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\frac{\partial }{\partial {x}_{j}}\text{\hspace{0.17em}}\left[\left(\mu \text{\hspace{0.17em}}+\text{\hspace{0.17em}}\frac{{\mu }_{t}}{{\sigma }_{\omega 3}}\right)\text{\hspace{0.17em}}\frac{\partial \omega }{\partial {x}_{j}}\right]\text{\hspace{0.17em}}+\\ \text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\left(1-{F}_{1}\right)2\overline{\rho }{\sigma }_{\omega 2}\frac{1}{\omega }\text{\hspace{0.17em}}\frac{\partial k}{\partial {x}_{j}}\text{\hspace{0.17em}}\frac{\partial \omega }{\partial {x}_{j}}\end{array}$

The production term is given as

${P}_{k}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}{\mu }_{t}\left(\frac{\partial {\stackrel{˜}{u}}_{i}}{\partial {x}_{j}}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\frac{\partial {\stackrel{˜}{u}}_{j}}{\partial {x}_{i}}\right)\text{\hspace{0.17em}}\frac{\partial {\stackrel{˜}{u}}_{i}}{\partial {x}_{j}}\text{\hspace{0.17em}}-\text{\hspace{0.17em}}\frac{2}{3}{\delta }_{ij}\left(\overline{\rho }k\text{\hspace{0.17em}}+\text{\hspace{0.17em}}{\mu }_{t}\frac{\partial {\stackrel{˜}{u}}_{k}}{\partial {x}_{k}}\right)\text{\hspace{0.17em}}\frac{\partial {\stackrel{˜}{u}}_{i}}{\partial {x}_{j}}$

The coefficients of the new model are a linear combination of the corresponding coefficients of the underlying models.

${\Phi }_{3}\text{\hspace{0.17em}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\text{\hspace{0.17em}}{F}_{1}{\Phi }_{1}\text{\hspace{0.17em}}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\left(1-{F}_{1}\right){\Phi }_{2}$

All coefficients are listed below.

$\beta *=0.09$, ${\alpha }_{1}=5/9$, ${\beta }_{1}=3/40$, ${\sigma }_{k1}=2$, ${\sigma }_{\omega 1}=2$, ${\alpha }_{2}=0.44$, ${\beta }_{2}=0.0828$, ${\sigma }_{k2}=1$, ${\sigma }_{\omega 2}=0.856$

The proper transport behaviour to predict separation is obtained by a limiter to the formulation of eddy viscosity. This is given as follows.

${\nu }_{t}\text{\hspace{0.17em}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\frac{{a}_{1}k}{\mathrm{max}\left({a}_{1}\omega ,\text{\hspace{0.17em}}S\text{\hspace{0.17em}}{F}_{2}\right)}$  and   ${\mu }_{t}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\overline{\rho }{\nu }_{t}$

Again F2 is a blending function similar to F1, which restricts the limiter to the wall boundary layer, as the underlying assumptions are not correct for free shear flows. S is an invariant measure of the strain rate as given below. The constant a1 is equal to 0.31.

$S\text{\hspace{0.17em}}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\text{\hspace{0.17em}}\sqrt{\left(\frac{\partial {\stackrel{˜}{u}}_{i}}{\partial {x}_{j}}\text{\hspace{0.17em}}+\text{\hspace{0.17em}}\frac{\partial {\stackrel{˜}{u}}_{j}}{\partial {x}_{i}}\right)\text{\hspace{0.17em}}\text{\hspace{0.17em}}\frac{\partial {\stackrel{˜}{u}}_{i}}{\partial {x}_{j}}}$

Blending Functions

The blending functions are critical to the success of the method. Their formulation is based on the distance to the nearest surface and on the flow variables.

${F}_{1}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\mathrm{tanh}\left({\mathrm{arg}}_{1}^{4}\right)$with ${\mathrm{arg}}_{1}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\mathrm{min}\left\{\mathrm{max}\left(\frac{\sqrt{k}}{{\beta }^{*}\omega y};\text{\hspace{0.17em}}\frac{500\nu }{{y}^{2}\omega }\right);\text{\hspace{0.17em}}\frac{4\rho {\sigma }_{\omega 2}k}{C{D}_{k\omega }{y}^{2}}\right\}$
and
$C{D}_{k\omega }\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\mathrm{max}\left(2\rho {\sigma }_{\omega 2}\frac{1}{\omega }\text{\hspace{0.17em}}\frac{\partial k}{\partial {x}_{j}}\frac{\partial \omega }{\partial {x}_{j}},\text{\hspace{0.17em}}1.0×{10}^{-10}\right)$
${F}_{2}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\mathrm{tanh}\left({\mathrm{arg}}_{2}^{2}\right)$ with ${\mathrm{arg}}_{2}\text{\hspace{0.17em}}=\text{\hspace{0.17em}}\mathrm{max}\left(2\frac{\sqrt{k}}{{\beta }^{*}\omega y};\text{\hspace{0.17em}}\frac{500\nu }{{y}^{2}\omega }\right)$

3.2 Thermodynamic Model

The products of combustion are assumed to be a perfect gas for the present study. The properties of this exhaust gas is evaluated from NASA CEA 40017,18program, which gives values of molecular weight, specific heats, coefficient of dynamic viscosity and thermal conductivity for the given initial temperature of the propellant composition and chamber pressure. For the propellant combination considered at ambient temperature, for a nominal chamber pressure of 100 × 105 Pa the temperature of the exhaust gases comes out to be 3500 K. These values for rocket exhaust gases are summarised in Table 1.. Within the constant diameter tube, across which losses need to be evaluated, the velocities are subsonic and only little changes in temperature and pressure is expected. Due to these small changes, assumption of constant properties is justified.

Table 1. Thermochemical properties of rocket exhaust

3.3 Discretization of Governing Equations

The CFX-10 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 the momentum integral equation and the spatial derivative terms in the integral equations are evaluated using finite element approach. The advective term is evaluated using a high resolution scheme explained in detail in reference 5. Following the standard finite element approach, shape functions are used to evaluate spatial derivatives for all the diffusion terms and pressure gradient terms. The set of discretized 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. An algebraic multigrid method is implemented to reduce low frequency errors in the solution of the algebraic equations.

3.4 Boundary Conditions and Initial Conditions

Typical computational domain is shown in Fig. 3 with various boundaries marked on it. At the computational domain inlet, a subsonic boundary condition has been used with total pressure and total temperature specified. A no slip adiabatic boundary condition is employed on the wall and outlet is modeled with supersonic boundary condition. Automatic fluid time step option has been chosen for time stepping. Thermochemical properties of the rocket exhaust used in the simulation are summarised in 1. The simulations are carried out using 32-cluster distributed computing platform. It takes around 12-14 hours for the convergence.

Figure 3.Locations of the applied boundary conditions for the simulations.

Grid independence study has been carried out for 70 mm diameter tailpipe with 248715 nodes for coarse grid and 408939 nodes for fine grid. The total pressure is normalised by the inlet total pressure (100 × 105 Pa), and the axial location is normalised by the length of the tube and nozzle assembly. The normalised total pressure variation with the normalised axial location across the tailpipe is compared and shown in Fig. 4. The total pressure drop from inlet to the exit of the nozzle is found to be 3.9 per cent for Coarse Grid and 3.7 per cent for fine grid. It can be observed that most of the difference in the total pressure drop occurs in the nozzle portion while in the tube portion the drops for both the grids are almost identical. Further the total pressure drop seems to reduce marginally with the refinement of grids. With these observations no further refinement has been carried out and the fine grid has been used for the simulations.

Figure 4.Total pressure with coarse and fine grids in 70 mm tailpipe.

The total pressure distribution is shown in Fig. 5, it can be noticed that there is higher total pressure in the tube compared to the nozzle. This indicates negligible loss in total pressure in tube whereas in the nozzle total pressure loss is significant. Also the nozzle contains some very weak wave structures arising due to the viscous effects near the wall in the throat region. The losses observed in the nozzle are the combined effect of viscous losses and losses through waves.

The total pressure drops across the pipe, nozzle and complete assembly have been calculated for all the blast tube geometries by considering area average total pressures at appropriate cross sections. These losses are presented in Fig. 6.

It may be observed from the Fig. 6 that only a small fraction of total pressure drop takes place in the tube. The loss is minimum in higher diameter tube because of less restriction of the flow leading to less viscous loss. A major portion of the total pressure drop takes place in the nozzle

Figure 5. Total pressure distribution in different tail pipes.

Figure 6.Total pressure losses in the nozzle, tube, and complete assembly for different tube diameters.

Table 2. Computed sea level thrust, Isp, and Total pressure drop

part. Also the variation in the total pressure drop in all the tail pipe geometry does not seem to differ considerably. For the increase of diameter from 70 mm to 85 mm, the cross section area of the blast tube increases by 47 per cent, while the total pressure loss in complete assembly is merely 0.16 per cent less for largest diameter tube than the smallest diameter tube. The computed sea level thrust and specific impulse are presented in Table 2 along with the numerical values of the total pressure drops. The variation in thrust is observed to be small for all the four geometries. Also the variations in the sea level specific impulses are found to be minimal.

The flow of hot gases through four different blast tube geometries are computed by solving 3-D RANS equations with SST Turbulence model using a commercial CFD software CFX-10. It is observed that the losses in total pressure in the rocket motor are less than 4 per cent and the blast tube is contributing less than 1 per cent. It is also found out that higher the blast tube diameter, the lesser the drop in the total pressure. Also a large increase (47 per cent) in the cross section area of the blast tube offers only a small improvement in the total pressure loss (0.16 per cent). The thrust and sea level specific impulse almost remain unchanged. It is concluded that blast tube geometry variations is not contributing significantly in the overall thrust and specific impulse in the rocket motor.

1.     Sutton, G.P. & Biblarz, O. Rocket propulsion elements.John Wiley & Sons Inc., 7th Edn, 2001, pp. 550-552.

2.     Mohan Reddy, T.; Subanandha Rao, A. & Sambasiva Rao, M. Design and development of propulsion system for antitank guided missile. Def. Sci. J., 1995, 45(3),  213-219.

3.     Sinha, P.K. & Javed, Afroz. Performance analysis of propulsive blast tube configurations using CFD. In the Proccedings of the seminar on Aerospace Technology: Challenges in the Millenium. December 15-16, 2003, AeSI, Hyderabad, pp.178-188.

4.     Tahsini, A. M. & Ebrahimi, M. Blast tube effects on internal ballistics of SRM, 2006, AIAA paper No. 2006-4958.

5.     ANSYS CFX, Release 10.0 :Installation and Overview, July 2005.

6.     CFX 5, Version 5.6, Installation guide, Introduction and New Feature Tutorials, 2003.

7.     Manna, P. & Chakraborty, D. Numerical investigation of transverse sonic injection in a nonreacting supersonic combustor. Proc. IMechE Part G: J. Aerospace Eng., 2005, 219(3), 205-215.

8.     Manna, P. & Chakraborty, Debasis. Numerical investigation of confinement effect on supersonic turbulent flow past backward facing step with and without transverse injection. J. Aerospace Sci. Techno., 2009, 61(2), 283-294.

9.     Dharavath, Malsur; Sinha, P.K. & Chakraborty, Debasis. Simulation of supersonic base flow – effect of computational grid and turbulence model. ASME J. Aerospace Eng., 2010, 224(3), 311-319.

10.   Dharavath, Malsur & Chakraborty, Debasis. Numerical simulation of underexpanded sonic jet. J. Aerospace Sci. Techno.2012, 14(4), 259-267.

11.   Aswin, G. & Chakraborty, Debasis. Numerical simulation of transverse side jet interaction with supersonic free stream.
J. Aerospace Sci. Techno., 2010, 14(5), 295-301.

12.   Saha, S.; Rathod, S.; Chandra Murty, M. S. R.; Sinha, P.K. & Chakraborty, Debasis. Numerical simulation of base flow of a long range flight vehicle. Acta Astronautica, 2012, 74(3), 112-119.

13.   Javed, A.; Manna, P. & Chakraborty, Debasis. Numerical simulation of dual pulse rocket motor flow field. Def. Sci. J., 2012, 62(6), 369-374.

14.   Saha, Soumyajit; Sinha, P.K. & Chakraborty, Debasis. CFD prediction of ramjet intake characteristics at angle of attack. In Proceedings of 8th National Conference on Air breathing engines and Aerospace Propulsion, held on December 12-14, 2006, at DIAT, Pune, pp. 97-107

15.   Javed, A. & Chakraborty, D. Numerical simulation of supersonic combustion of pylon injected hydrogen fuel in scramjet combustor. J. Institution Eng. (India), 2006, 87, 1-6.

16.   Menter, F.R. Two – Equation Eddy – Viscosity turbulence models for engineering applications. AIAA J., 1994, 32(8), 1598-1605.

17.   Gordon, S. & McBride, B. J. Computer program for calculation of complex chemical equilibrium compositions and applications - I. Analysis, 1994, NASA RP-1311

18.   Gordon, S. & McBride, B.J. Computer program for calculation of complex chemical equilibrium compositions and applications - II. Users manual and program description, 1996, NASA RP-1311.

 Dr Afroz Javed Dr Afroz Javed did his MTech from Indian Institute of Technology (IIT) Madras, in 1996 and PhD from Indian Institute of Science (IISc), Bengaluru in 2013. Presently he is working in the Directorate of Computational Dynamics, Defence Research & Development Laboratory (DRDL), Hyderabad. He has published 5 journal papers and 3 conference papers. His research interests include: Computational fluid dynamics (CFD), propulsion, combustion instability, and turbulent compressible mixing. Mr Pradip Kumar Sinha has completed his BE from NIT, Durgapur and M.E. in Mechanical Engineering from Jadavpur University, Kolkata. He is working in the area of Computational Fluid Dynamics at Defence Research and Development Laboratory, Hyderabad from 1990. He has contributed significantly in the development and application of indigenous Euler codes for missile aerodynamic characterization. He has received DRDO Young scientist award 1998. He has nearly 30 publications in various conferences and journals. Dr Debasis Chakraborty obtained his PhD (Aerospace Engineering) from Indian Institute of Science (IISc), Bengaluru. Presently, he is working as Technology Director, Computational Fluid Dynamics (CFD), DRDL, Hyderabad. He was the guest editor of special issue of Defence Science Journal on CFD. His research areas includes: CFD, aerodynamics, high speed combustion and propulsion. He has more than 120 technical papers in referred journals and conferences.