Numerical Simulations of Combustion Instability and Pressure Oscillation in Solid Rocket Motor

Numerical Simulations of Combustion Instability and Pressure Oscillation in Solid Rocket Motor

L.G Praveen laws, S.M.Lakshmiganthan

Sree Sastha institute of engineering and technology, Chennai- 600 123.

Abstract

Solid rocket motor frequently experience unsteady gas motions and combustion instability. Pressure oscillations are a well- known problem of large solid rocket motors (SRMs), which lead to thrust oscillations and to significant dynamic loads that often cause the need for dampers on the payload, reducing the payload mass capacity and launcher usability. Instability in solid rocket motor is governed by the flow behaviour of combusted gas, pressure fluctuations and acoustic resonances in the Combustion chamber. In present investigation the computational analysis of high temperature air passing through a conventional solid rocket motor is studied, which describes vortex shedding, instabilities inside the core section of motor and thrust oscillation. The computation scale to scale model of

Shanbhogue[10] experimental setup has been developed (AIAA

paper 2003-4632, 39th AIAA/ASME/SAE/ASEE Joint propulsion Conference and Exhibit, 2003). The model is meshed to three different cells, i.e. 40000cells, 50000 cells and 60000 cells respectively. A commercial fluent software is used for the stimulating the hot flow inside the combustion chamber and the frequency of oscillation are predicted. The flow is assumed to be incompressible, and axi-symmetric. The analysis was carried out for both laminar and turbulent calculation .The observations are presented in the form of FFT plot which shows pressure-time data. The comparison of computational and experimental results has been investigated and found that there are significant changes in pressure wave at vortex region as flow changes from laminar to turbulent.

1. Introduction

   Combustion instability leads to unsteady thrust resulting in structural vibrations, uncomfortable ride for astronauts or payload, difficulty with guidance systems, and in extreme cases, erosion of the chamber wall, resulting in catastrophic motor failure. Most damaging frequencies occur in the low Hz to low kHz band. Periodic vortex shedding as a source of acoustic oscillations in solid rocket motors (SRM) has attracted the attention of several researchers during the last decade. The reason for this interest is an unpredicted oscillatory behaviour observed in several motors during test firing or during flight. The first evidence of vortex driven oscillations was encountered in motors with complex grain geometry, used in the upper stages of ballistic missiles. If the heat release during the combustion process fluctuates, it could couple or lock-on with other fluctuations, thereby initiating a feedback cycle.

   Oscillatory behaviour has been observed in the Space Shuttle. These oscillations have been attributed to a periodic vortex

   shedding due to the presence of an obstacle in the main flow constituted by the frontal thermal protections. When the vortex- shedding frequency is the same as the natural frequency of the chamber, the pressure oscillation reaches a maximum and this can reduce the overall performances and usability of the launcher. The rate of primary decomposition of solid propellant element and the rate of activation of intermediate products

   during the time lag are depended upon the gas oscillation [1].

   The concentration changes were governed by the rate processes of mixing, phase change, and chemical reaction in a rocket motor [2]. The unsteady flow in a combustion chamber could be

   governed by nonlinear partial differential equations which are difficult to study mathematically as well as computationally [3]. Computations were carried out for a cylindrical chamber in

   which the mean flow velocity, the stability of standing pressure oscillations were changed due to the distribution of combustion product during steady operation[4]. The vortices generated by

   acoustic oscillations were second order viscous effect, and tangential modes gave a single vortex swirling about their axis an axially perforated solid propellant rocket [5] and the vortex

   formation was affected oscillatory flow and sudden expansion of the chamber, vortex burning was a function of localized parameters like space and time, the vortex burning and corresponding vortex location were the function of chamber geometry, velocity field, characteristics of chemicals and

   hydrodynamic times[6]. The chamber pressure and

   characteristics velocity was lower for small characteristic lengths, than with larger L* in similar solid rocket motors[7].The investigated on cold flow effect on the instabilities

   in Ariane 5 solid propellant booster shows that the resonant pressure fluctuations which occurred at different velocities coincided to the existence of vortices[8]. Cold-Gas simulations

   would be easy and useful tools for analysis of vortex shedding and acoustics pressure oscillation [9]. Shanbhogue et al. [10] studied experimentally the aero acoustic instability in a rocket

   motor with the presence of a finocyl grain and cold flow. They observed that frequencies at which the pulsations occurred were depended on the upstream mean velocity, the upstream diameter and downstream length. The generated pulsations contained two discrete bands of Strouhal numbers and purely depended on the

   downstream length to step height ratio [11]. Utilizing CFD

   analysis, significant provided simplified understanding of internal flow of solid rocket motor flow and the computational results with experimental measured data at static tests of the reusable solid rocket motor. Detached eddy simulation (DES) provided accurate turbulence gradient properties compared to wall layer modeling (WLM). Modeling multiphase effect could be achieved more rigorously through the sub grid modeling in

   LES and would be more difficult in the DES. RANS approach for first order estimation of nozzle wall heat and momentum

   transfer and the statistical properties of turbulence could be better model by using RANS [12]. The main objective of the present work is to numerically simulate the flow instabilities in

   solid rocket motors. Combustion is not considered in this investigation. However, the basic understanding of the hot flow as well as the methodology developed in this work should hopefully provide a sound basis for extending the analysis to reacting flows as well. To this end, in the present work, the flow is assumed to be incompressible, axi-symmetric and turbulent. Numerical calculations have been carried out to determine the vortex shedding frequency and the effect of varying the flow rate on this frequency.

2. Solution methodology

   The experiment is to determine the frequency of the oscillations of the flow field and the impact of geometric parameters without the complicating effects of combustion. The flow is hot and non-reacting with air as the working fluid, incompressible, viscous flow and uniform injection of fluid normal to boundary and using classical Navier-Stokes equation in cylindrical co-ordinates.The maximum velocities seen in the experiment were of the order of about 75 m/s and the temperature range is of 1000 K. The flow is governed by the incompressible Navier-Stokes equations in cylindrical co- ordinates. In the present work, the realizable k – e model is used to model turbulence effects as it has been found to give good and reliable predictions over the entire range of flow rates under consideration and problem of vortex shedding is being modeled and compared with experiental results of

   Shanbhogue[10] experimental model and computational cold

   flow model, In addition for turbulence calculation, the inlet turbulence intensity is set to 3%. This value has been estimated on the relationship,

   I= 0.16(ReD)-0.125 for fully developed pipe flow.

   1. Modeling

      The objective of Shanbhogues[10] experimental work was to find out the frequency of oscillation on the flow field and geometrical parameters. In the experiments performed, the nozzle was not choked, thus allowing open-open boundary conditions. The computational analysis has been carried out using CFD, A commercial software GAMBIT & FLUENT .The model is generated in GAMBIT with complex geometry and meshed along boundary condition. Assumptionsconsidered for

      Shanbhogues[10] experimental model :

      • The model is generated in 2-D and 3-D, but the computational analysis is carried out in 2-D model only.

      • The nozzle inlet was not chocked, thus allowing open- open boundary condition.

      • Nozzle inlet was chocked, resulted an acoustic open- closed boundary condition.

   2. Configuration

      Shanbhogues experimental model and 2-D computational model has been shown in the fig.1 & fig.2. Scale to scale model of the experimental setup. Three different numbers of triangular and quadrilateral unstructured meshes are used for the

      discretization of the computational domain. Coarse and fine grids are used with 40000cells, 50000cells and 60000cells has been choosen.

   3. Solution and Computational method

   The computation through finite volume method is carried out for compressible, segregated solver with implicit and unsteady type scheme. The problem is evaluated with time stepping procedure and instability has been computed with FLUENT software, as the velocity and pressure distribution at the outlet are not known, outflow boundary condition is used here. No slip condition is imposed on the wall and symmetry conditions are imposed along the axis, The ReynoldsAverage Navier- Stokes has been used to find the turbulent intensity in complex flow in solid rocket motor. The computational model has been generated using GAMBIT software and imported to FLUENT using mesh file option, to solve the continuity, momentum and energy equation segregated and implicit solver formulation has been used .The grid is checked, material properties ,boundary conditions and operating condition have been assigned suitably. Initially steady state calculation has been carried out to remove the transient and unsteady calculation has been imposed when

   the residuals of all governing equations are of the order 10-5.

   For unsteady calculation the time step has been kept as 0.0001sec. Simple scheme has been used along with the segregated implicit solver. All the solution in this work has been obtain in the second order accuracy in space and time derivative. Turbulent calculations have been performed using the realizable k- model. The unsteady pressure data has been monitored at different location in computational domain and one of the locations of these pressure data has been taken in computational domain to compare the results with Shanbhogue experimental model. First Fourier Transfer (FFT) of these pressure data has been carried out to extract the dominant frequency of oscillation at different location in computational domain.

3. Results and Discussion

   In this section, numerical solutions obtained using the methodology outlined in the previous section are presented and discussed. The numerical solutions obtained using FLUENT

   have been compared with the experimental results reported by Shanbhogue et al[10] available, the vortex shedding frequency

   for various flow rates and downstream lengths have been compared. Results from laminar simulations for low flow rates are presented first, followed by turbulent calculations for higher flow rates.

   1. Laminar Calculations

      The laminar calculations have been carried out for lower flow rates (800-1000 lpm). Grid adaption has been done twice in the region of interest to capture the formation of small vortices, their growth and convection. Grid independence studies have been done using 40000, 50000 and 60000 cells.The vortex shedding frequency for two geometric configurations on the three grids. Estimated values of the wall Yplus value for this grid shows it to be less than 30 along the entire wall, indicating good resolution of the boundary layer. The laminar calculation is carried out for the different flow rates

      of 800lpm, 900lpm and 1000lpm because the velocities observed in the region of interest are quite low. First the laminar calculation is carried out for the flow rate of 800lpm with 40000 cells. Later grid adaptation has been done in the region of interest to capture the formation of small vortices, their growth and convection. By grid adaptation technic the laminar calculation is carried out for 50000 cells and 60000 cells. Finally datas such as contour of vorticity magnitude, and stream function and frequency curve data for the flow rate 800lpm with downstream length (Ld) 66mm and with 40000 cells, 50000 cells and 60000 cells has been collected .Datas collected for 800lpm with downstream length (Ld) 66mm show in fig.3,fig.4,fig.5,fig.6 & fig.7. Again the laminar calculation is carried out for 900lpm, with three different cells and with the same downstream length (Ld) 66mm has been carried out and frequency monitored due to the pressure variation, stream function and contour of vorticity has been collected. Finally laminar calculation for the flow rate 1000lpm at downstream length 66mm with 40000 cells, 50000 cells and 60000 cells has been absorbed and datas such as the contours of vorticity magnitude and stream function and frequency curve for 1000lpm has been collected, The highest fluctuation vortex shedding frequency resulted from computation has been presented in Table 1. It indicates that there is a significant pressure fluctuation and it can be related to the vorticity generation. The contours of vorticity and stream function have also represented this variation. The justification for this is based on the following observations. In this flow field, at the point of separation, the boundary layer is laminar .Also, the separated free shear layer is initially laminar. Based on the above reasons, it can be concluded that the wall bounded flow behind the backward facing step is an unsteady, quasi-laminar flow. It has been observed that fluctuations are higher with coarse grids than that of the fine grids. Grid independence is the term used to describe the improvement of result by using successive smaller cell size for calculation. The dominant frequency has been noted and it has been found that oscillation sustained for a longer period of time for the fine grids. so the fine grid values has been chosen for comparison with the experimental cold flow which are found much closer to the experimental analysis, the value of 50000cells fine grid have been taken for the comparison and shown in the Table 2.The comparison between computational and experimental values is found good.

   2. Turbulent calculations

      The turbulent calculation has been carried with flow rates from 1700 lpm to 2300lpm. Similar to laminar computations, the contour of velocity magnitude, stream function, vorticity magnitude, turbulence kinematic energy, vorticity magnitude and Fast Fourier transform of unsteady pressure data have been observed. The special interest has been taken to determine the frequency of vortex shedding at the region of interest for the different flow rates. The computational analysis has been carried out for 1700lpm (ReD = 97230) with coarse grid 40000 cells and fine grids (50000 cells , 60000 cells) and downstream length of 66mm.The stream function, vorticity magnitude, and turbulent kinematic energy, contours of velocity magnitude and Fast Furious Transformation of unsteady pressure data have been bserved and shown in fig.8 to fig.13 .The FFT plot for

      unsteady pressure data has been collected and the dominant frequency is 1475Hz with coarse grids and 502 with fine grids which is closer to the experimental data. Similarly, using 40000 cells, 5000cells and 60000 cells and volume flow rate of 1800 lpm (ReD =102950), the results has been obtained. Contours of vorticity magnitude for fully turbulent calculations (Q=1800 lpm, Re=102950) clearly show that the shed vortices are not growing and adverting. This is because the fully turbulent computations are excessively dissipative, causing rapid dissipation of the rolled up vortices within a short distance. Again adopting the same cells value and changing the volume flow rate into 1900lpm.The stream function, vorticity magnitude, turbulent kinematic energy, contours of velocity magnitude and also the FFT plot, the results has been collected. Similarly analysis has been carried out till 2300 lmp. Numerically obtained values for the vortex shedding frequency along with the experimental value have been shown in table 3. From the table it has been found that frequency is higher with 40000 cells when compared with 50000 cells and 60000 cells. From the comparisons it has been observed that the numerically calculated values for vortex shedding are closer to the experimental value. Out of the three grid values, the values of 50000 cells have been chosen for the comparison with experimental analysis. Since the values of 50000 cells are found more accurate with the experimental data. Turbulence arises from the flow at high Reynolds number and classified as shear layer instability. The shear layer instability is related to the sharp velocity gradient at separation point and downstream length. At large Reynolds number the shear layer may be idealized as surface of velocity discontinues. Shear layer are unstable and curl up into vortices when viscous forces are negligible. The velocity fluctuations are zero at the wall implies that the Reynolds stresses vanishes or Reynolds shear stress is zero and the only shear stress exerted directly at the wall is the viscous stress. But, away from the wall, fluctuating velocity components increase rapidly and Reynolds stresses are much larger than the viscous stresses. From the above observation, it is found that finer grids provide more accurate results in comparison with the coarse grids. The results from fine grids are very close to the experimental results. The vortex formation, growth and advection can be clearly seen from the contours. The vortices that are forming at the edge of the backward facing step are advected towards the nozzle and finally impinge on the wall. The effect of the grid refinement on the frequency of vortex shedding can be seen there is still a lot clutter in the spectrum. For low flow rates, the velocities in the region of interest are quite low; secondly, turbulence models tend to be excessively dissipative and suppress vortex formation, growth and shedding .It is clear that turbulent production is suppressing the vortex shedding, while the laminar calculations with a fine grid are able to predict the frequency quite accurately when compared with the experimental values.

      Wall

      Decoupler Backward

      Fig. 5 FFT plot for Q=1000lpm for 60000 cells.

      Step

      Step

      D1

      L L2 d

      Shanbhogues experimental model.

      Wall

      Fig. 1

      Fig. 6 Contour of stream function for Q=800lpm.

      Velocity inlet Axis Outflow

      Fig. 2 2-D Computational model.

      Fig. 3 FFT plot for Q=800lpm for 40000 cells

      Fig. 4 FFT plot for Q=900lpm for 50000 cells.

      Fig. 7 Contour of Vorticity Magnitude for Q=800lpm

      Fig. 8 FFT plot for Q = 1700lpm for 40000 cells.

      Fig. 9 FFT plot for Q = 1700lpm for 50000 cells.

      Fig. 10 FFT plot for Q = 1700lpm for 60000 cells.

      Fig. 11 Contour of stream function for Q =1700lpm.

      Fig.12 Contour of vorticity magnitude for

      Q =1700lpm.

      Fig. 13 Contour of turbulence kinetic energy for Q

      =1700lpm.

      Table 1 Computational results of vortex shedding frequency using different cells for laminar flow

      Downstream length Ld mm	Flow rate in lpm	Vortex shedding frequency (Hz)
      		Computational ( hot flow)
      		4000 0 cells	5000 0 cells	6000 0 cells
      66	800	393	373	393
      66	900	411	406	408
      66	1000	436	424	426

      Table 2 Comparison of vortex shedding frequency with computational hot flow value and experimental cold flow value for laminar calculation

      (50000 cells computational values are much closer to experimental data, so it is chosen for the comparison)

      Downstream length Ld mm	Flow rate in lpm	Vortex shedding frequency (Hz)
      		Comput ational ( hot flow)	Experi mental ( cold flow)
      		50000ce lls
      66	800	390	271
      66	900	406	268
      66	1000	424	274

      Downstr eam length Ld mm	Flow rate in lpm	Vortex shedding frequency (Hz)
      		Computational ( hot flow)
      		40000ce lls	50000ce lls	60000ce lls
      66	1700	1475	506	502
      66	1800	776	509	504
      66	1900	589	513	509
      66	2000	752	519	514
      66	2100	864	889	886
      66	2200	922	894	891
      66	2300	952	911	908

      Downstr eam length Ld mm	Flow rate in lpm	Vortex shedding frequency (Hz)
      		Computational ( hot flow)
      		40000ce lls	50000ce lls	60000ce lls
      66	1700	1475	506	502
      66	1800	776	509	504
      66	1900	589	513	509
      66	2000	752	519	514
      66	2100	864	889	886
      66	2200	922	894	891
      66	2300	952	911	908

      Table 3 Computational results of vortex shedding frequency using different grid for turbulent calculation

      Table 4 Computational results of vortex shedding frequency using fine grd (50000 cells) and their comparison with experimental cold flow value for turbulent calculation

      (50000 cells computational values are much closer to experimental data, so it is chosen for the comparison)

      Downstrea m length Ld mm	Flow rate in lpm	Vortex shedding frequency (Hz)
      		Computati onal ( hot flow)	Experimen tal ( cold flow)
      66	1700	506	511
      66	1800	509	515
      66	1900	513	519
      66	2000	519	525
      66	2100	889	891
      66	2200	894	898
      66	2300	911	913

4. Conclusion

It has been found that the vortex shedding frequency is more oscillating at initial stage and then the oscillation gradually dies out. The values from fine grids are more accurate and closer to experimental data than that of course grids. Vortex shedding frequency has been found to be dependent on the flow rate. Vortex shedding frequency which is obtained by the turbulent calculation with coarse and fine grids has been found to be increased with the flow rate. The vortices that are forming at the edge of the backward facing step are adverted towards the nozzle and finally impinge on the wall. The increment quantity of the frequency using fine grid has been found almost linear and this has not been found with coarse grid. Again, it has been found that only one dominating peak and few low secondary peaks noted from the FFT plot.The fluctuation is irregular, and it is also observed that fluctuations are higher with fine grids than that of the coarse grid. The frequency of oscillation has been found to be more prominent in the initial stage with coarse grid, after that this dies out. Vortex formation and its growth has been seen in the contours and these vortices finally impinge on the wall of the nozzle. Numerically predicted values for the frequency of the vortex shedding agree quite well with the experimental data.

Acknowledgements

The author expresses his heartfelt thanks to Prof. Dr. P.K. Dash, for guiding to work on this particular topic and his continued assistance provided for this work. The author expresses his sincere thanks to Prof. Dr. D.G. Roychowdhury, Dean (Research), Hindustan University, for guiding and encouraging during the thesis work. And his keen interest, constant help and discussion during the work.

Reference

1. Cheng, Sin-I(1954),High frequency combustion instability in solid propellant rockets, (Part -1 and Part- 2), Jr. of Jet Propulsion.

2. Wise, H., and Lorell, J. (1954), Excitation of oscillation by chemical reaction, Jr. of Jet Propulsion, Technical Reports.

3. Malhotra, S. (2004), On combustion instability in solid rocket, Ph.D. thesis, California Institute of Technology, California, USA.

4. Culic, F.E.C. (1963), Stability of high frequency pressure oscillation in rocket combustion chamber, Jr. of AIAA, Vol.1,

   No.5

5. Sotter, G., and Swithenbank, J. (1964), Vortex generation in solid propellant rockets.

6. Matveev, K.I., and Culic, F.E.C. (2003), A model for combustion instability involving vortex shedding, Jr. of Combustion Science and Technology, Vol. 175, No. 6, pp. 1059-1083.

7. Schoyer, H.F.R (1983), Incomplete combustion: a possible cause of combustion instability, Jr. of AIAA, Vol. 21, No. 8.

8. Anthoine , J., and Olivari, D. , and Plaquart, P., (1998), Cold flow investigation of the flow acoustic coupling in solid propellant booster, 36 th AIAA Aerospace Science Meeting and Exhibit, Reno, Nevada, Paper 98-0475.

9. Vetel, J., Plourde, F., and DoanKim, S. (2005), Cold gas simulations of the influence of inhibitor shape in combustor combustion.

10. Shanbhogue, S.J., Sujith, R.I., Chakravarthy, S.R. (2003), Aero acoustics of rocket motors with finocyl grain.39thAIAA Joint Propulsion Conference and Exhibit, Alabama, USA.

11. Culic, F.E. (2004) Combustion instabilities in solid propellant rocket motors, AIAA Conference Paper, Report No.

    A291524

12. Wasistho, B., Najjar, F.M., and Haselbacher, A. (2005), Effect of turbulence and particle combustion on solid rocket motor instabilities, 41st AIAA/ASME/ SAE/ ASEE/ Joint Propulsion Conference and Exhibit, Jucson, Arizona.