- Email: [email protected]
Effect of swirl on the regression rate in hybrid rocket motors
Aerospace Science and Technology 29 (2013) 92–99
Contents lists available at SciVerse ScienceDirect
Aerospace Science and Technology www.elsevier.com/locate/aescte
Effect of swirl on the regression rate in hybrid rocket motors C. Palani Kumar 1 , Amit Kumar ∗,2 National Centre for Combustion Research and Development, Department of Aerospace Engineering, Indian Institute of Technology Madras, Chennai, Tamil Nadu-600036, India
a r t i c l e
i n f o
Article history: Received 7 July 2012 Received in revised form 29 January 2013 Accepted 30 January 2013 Available online 13 February 2013 Keywords: Hybrid rocket Regression rate Swirl flow Numerical simulation
a b s t r a c t In this work, the effect of swirl on regression rate (rb ) in a hybrid rocket motor is investigated numerically. The rb increased monotonically with inlet swirl number and was also found to depend on the inlet swirl profile. The swirl velocity profiles with the peak closer to the axis yielded higher rb . Parametric study on fuel grains of various lengths and diameters ( L / D 25) shows that swirl is more effective in improving the average rb for short grains ( L / D < 5) and large diameter grains. © 2013 Elsevier Masson SAS. All rights reserved.
1. Introduction Hybrid rockets have many distinct advantages over conventional chemical rockets, including simplicity, fuel costs, hardware costs, and safety. However, one concern with the operation of these rockets is low rb of the solid fuel [14,15]. While multi-port complex grain design can be used for obtaining reasonable thrust with low regression rates, this results in larger residuals and low fuel loading. Further the grain integrity also becomes a factor of concern [1]. Extensive research over the past few decades on hybrid rocket motors has contributed greatly towards the fundamental understanding of its working. One of the rb enhancement techniques established uses swirling flow inside the combustion chamber [9, 10,16,22,24]. Swirling flow in the combustion chamber can be achieved by either tangentially injecting oxidizer from the tail end [10], or using a swirl type injector [9,16,22,24] or using helical grain configuration of solid fuel [16,24]. Although rb improvements are reported in the literature [9,10, 16,22,24], the amount of swirl is most often not quantified. There is a lack of basic understanding of the effect of swirling flow on rb of the hybrid rocket motor. The present numerical study aims to systematically understand the underlying physics and relate rb to swirl strength, quantified by a non-dimensional swirl number.
* 1 2
Corresponding author. Tel.: +91 44 22574019; fax: +91 44 22574002. E-mail address: [email protected] (A. Kumar). Ph.D. Research Scholar. Associate Professor.
1270-9638/$ – see front matter © 2013 Elsevier Masson SAS. All rights reserved. http://dx.doi.org/10.1016/j.ast.2013.01.011
2. Numerical model and solution 2.1. Geometrical configuration and computational domain Fig. 1 presents the schematic of the hybrid rocket combustion chamber. Invariance in azimuthal direction is assumed to reduce the computational domain to a 2D axi-symmetric configuration. The shaded region in Fig. 1 is the computational domain. The domain consists of a combustion chamber and a CD nozzle. The dimensions of combustion chamber, nozzle and the details on grid are discussed later in Section 4. 2.2. Governing equations The processes occurring inside the hybrid rocket combustion chamber can be adequately described by basic flow equations of continuity, momentum, energy and species. The governing equations in 2D cylindrical coordinates are summarized below. Continuity equation
∂ ∂ ∂ ρ vr (ρ ) + (ρ v x ) + (ρ v r ) + =0 ∂t ∂x ∂r r
(1)
Axial momentum equation
∂ 1 ∂ 1 ∂ (ρ v x ) + (r ρ v x v x ) + (r ρ v r v x ) ∂t r ∂x r ∂r ∂p 1 ∂ ∂ vx 2 =− + rμ 2 − (∇ · v ) ∂x r ∂x ∂x 3 ∂ v x ∂ vr 1 ∂ rμ + + r ∂r ∂r ∂x
(2)
C.P. Kumar, A. Kumar / Aerospace Science and Technology 29 (2013) 92–99
93
Nomenclature A, As C C∗ D i ,m Di , D0 Dω Dt E a , E as G˜ k , G ω h, hi HR, HP J j k kr L p rb
Arrhenius pre-exponential factor for gas phase reaction and solid fuel pyrolysis Concentration of species Characteristic exhaust velocity . . . . . . . . . . . . . . . . . . . . m/s Diffusivity of species i in the mixture . . . . . . . . . . . m2 /s Initial and final port diameter . . . . . . . . . . . . . . . . . . . . . . . m Cross diffusion term Throat diameter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m Arrhenius activation energy for gas phase reaction, solid fuel pyrolysis Production of k and ω Sensible enthalpy of mixture and species i . . . . kJ/kg/s Heat of gas phase reaction and solid fuel pyrolysis Diffusion flux of species j Turbulent kinetic energy Reaction rate constant Fuel grain length . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m Pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . N/m2 Regression rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m/s
R0 Ru Sct T s, T∞ v x, vr , W , Wi x, r , z Yi Yk, Yω
α Γk , Γω λeff λs μ, μt
ρ ω
Radius of the combustion chamber . . . . . . . . . . . . . . . . . m Universal gas constant Turbulent Schmidt number Surface temperature and initial temperature of fuel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . K v z Axial, radial and swirl (tangential) velocity . . . . . m/s Global reaction rate and consumption rate of species i . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . mol/m3 /s Cylindrical coordinates in axial radial and tangential direction Mass fraction of species i Dissipation of k and ω Thermal diffusivity Effective diffusivities of k and ω Effective thermal conductivity Thermal conductivity of the solid fuel Molecular and turbulent viscosities . . . . . . . . . . . . N s/m2 Density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . kg/m3 Specific dissipation rate
The detailed turbulence models like RSM require rigorous closure strategies and grid quality requirements and involve complexities in computation and convergence [11]. Hence these are not suited for the parametric study intended here. Simpler models like lowRe k–ε turbulence models are widely popular in literature [3–5,15]. However, under intense injection the k–ε model fails in qualitative flow prediction [26]. Although the standard k–ω model overcomes this deficiency, it is sensitive to the free stream values of ω [20]. Hence, SST k–ω model was developed [21] which combines the advantages of both k–ε (free stream accuracy) and standard k–ω (near-wall accuracy). Thus SST k–ω was used to predict turbulence in the present study. Fig. 1. Top: Schematic of a single port cylindrical hybrid rocket motor showing solid fuel grain, combustion chamber and a convergent–divergent nozzle. The shaded region is the computational domain. Bottom: Schematic of the computational domain and the boundary types.
Radial momentum equation
∂ 1 ∂ 1 ∂ (ρ v r ) + (r ρ v x v r ) + (r ρ v r v r ) ∂t r ∂x r ∂r ∂p 1 ∂ ∂ v x ∂ vr =− + rμ + ∂r r ∂x ∂r ∂x 1 ∂ ∂ vr 2 + rμ 2 − (∇ · v ) r ∂r ∂r 3 − 2μ
vr r2
+
2μ 3 r
(∇ · v ) + ρ
v 2z r
(3)
∂ (ρ E ) + ∇ · v (ρ E + p ) ∂t h j J j + (τ¯¯ · v ) + H R = ∇ · λeff ∇ T −
(7)
(8)
(4)
Here E is defined as
E =h−
Tangential momentum equation
∂ 1 ∂ 1 ∂ (ρ v z ) + (r ρ v x v z ) + (r ρ v r v z ) ∂t r ∂x r ∂r 1 ∂ 1 ∂ 3 ∂ vz vr v z ∂ vz = + 2 −ρ rμ r μ r ∂x ∂x ∂r r r r ∂r
(6)
j
In the above equation the gradient of velocity vector is defined
∂ v x ∂ vr vr + + ∂x ∂r r
∂ ∂ ∂ ∂k Γk + G˜ k − Y k (ρ k) + (ρ ku i ) = ∂t ∂ xi ∂xj ∂xj ∂ ∂ ∂ ∂ω Γω + Gω − Yω + Dω (ρω) + (ρωu i ) = ∂t ∂ xi ∂xj ∂xj Energy equation
as
(∇ · v ) =
SST k–ω equations are given by
p
ρ
+
v2 2
(9)
Species transport equation
(5)
The flow in the combustion chamber and nozzle is likely to be turbulent, therefore, an appropriate turbulence model is required.
∂ (ρ Y i ) + ∇ · (ρ v Y i ) = −∇ · ( J i ) + W i ∂t
(10)
Here,
J i = −(ρ D j + μt /Sct )∇ Y j
(11)
94
C.P. Kumar, A. Kumar / Aerospace Science and Technology 29 (2013) 92–99
2.3. Gas phase chemistry model Combustion model with detailed chemistry is desired for accurate predictions. However, proper reaction mechanism and related kinetics are not readily available in the literature. Nevertheless, thermodynamic effect of multi-species can be readily incorporated [6]. Therefore, a global chemistry [8] is adopted here where heat of combustion is tuned to be consistent with thermodynamics following [12]. This is expected to give reasonable prediction of temperature field and thus heat transfer to the solid fuel. As butadiene is the main pyrolysis product of HTPB, the single step global reaction is taken as
C4 H6 + 5.5O2 → 4CO2 + 3H2 O
(12)
The reaction rate W is given by assuming a second order reaction as
W = kr C fuel C oxidizer = kr C C4 H6 C O2
(13)
kr can be expressed by Arrhenius law as kr = A (exp(− E a / R u / T )). Here, A = 8.8E+11 and E a = 1.2637E+08 J/kmol [4].
( S w ). Swirl number is defined as the ratio of axial fluxes of angular momentum to axial momentum, non-dimensionalized by inlet radius [7].
Sw =
ρ v x v z r 2 dr R 0 ρ v 2x r dr
(18)
Jiang et al. [7] also provided the analytical axial and swirl velocity profiles which can be considered as equivalent for swirling flow inlet boundary conditions in CFD as
vx =
−1 f x
r 2 − R 20
4 μ −1 f z 2 vz = r − R 0r 3 μ
(19a) (19b)
Substituting Eqs. (19a) and (19b) in Eq. (18) results in S w as a function of body force terms f x and f z . For a chosen v x ( f x = 1), v z is varied (different f z ) to achieve different inlet S w . The nozzle walls were considered adiabatic. Atmospheric pressure was specified at the nozzle exit. 2.6. Numerical method of solution
2.4. Pyrolysis model The pyrolysis of the solid fuel (HTPB grain) was modeled by the zeroth order Arrhenius equation following [23]. The mass of fuel released by pyrolysis is given by
˙ = ρs rb = A s exp(− E as / R u T s ) m
(14)
2.5. Boundary conditions The boundary types are mentioned in Fig. 1. The fuel inlet was defined as an interface where the mass and energy balance was applied. The reduced energy balance equation is given by
λeff (∂ T /∂ y ) g − λs (∂ T /∂ y )s + ρs rb ( H P ) = 0
(15)
The mass balance at the interface is given by
˙ = ρs rb = −ρ v m
(16)
A part of the thermal energy released in the reaction zone is convected to the solid–gas interface (1st term in Eq. (15); qconvection ). A part of this thermal energy is then transferred into solid fuel by means of conduction in the solid phase (2nd term; qconduction loss ). The remaining energy is utilized for solid fuel pyrolysis (3rd term; qpyrolysis ). The heat lost inside the solid fuel was determined by an analytical solution of 1D heat conduction inside the solid with moving hot boundary. The analytical temperature profile inside the solid fuel is given by
T ( y ) = ( T s − T ∞ ) exp(−rb y /α ) + T ∞
(17)
The radiation heat transfer is neglected in the present study. Simulations [12,15] have shown that the effect of radiation (compared to simulations without radiation) on change in net heat flux incident on the solid fuel surface is small (typically about 10%). This is because the radiation heat flux on the fuel surface increases the blowing velocity when radiation is included. This in turn reduces qconvection . Radiation from flame also reduces flame temperature which further reduces qconvection . The net result is only a marginal increase in net heat flux on the fuel surface (or a small increase in rb [19]). A fixed average mass flux was specified at the oxidizer inlet in all simulations. Swirl at the inlet was quantified by swirl number
The governing equations along with the boundary conditions were solved numerically using a pressure based, double precision, unsteady solver [2]. Second order upwind scheme was used for spatial discretization of the convective terms whereas second order central scheme was used for all other terms in the transport equations. The second order implicit formulation was used for temporal discretization. First order discretization schemes are not used because they are less accurate. Second order discretization schemes were opted as the velocity profiles and the local rb along the axis did not change appreciably between the second and third order discretization schemes. The discretized governing equations were iterated with the time step of 0.0002 s and 50 iterations per time step. The converged solution for short motors was attained typically in about 1000 time steps. Convergence was assessed by the condition that the scaled residuals reduce below a value of 1E−5 for all transport equations and become invariant with time. 3. Validation of numerical model Detailed experimental data on swirling reacting flow inside a hybrid rocket motor is not available in the literature. Hence the present numerical model was validated with cold turbulent swirling flow in a pipe [25] and 2D hybrid rocket experiments and simulations [15]. 3.1. Swirl flow in a pipe Steenbergen and Voskamp [25] reported experimental measurements on turbulent swirl velocity profiles at various axial locations (x/ D = 4.4, 7.6 and 14.5). Since the inlet velocity profiles were not known, the reported velocity profiles at x/ D = 4.4 were provided at the inlet in the present simulation. The computed profiles at the sections x/ D = 7.6 and 14.5 are compared with the experimental values in Fig. 2(A). The predicted velocity profiles match well with the measured profiles. 3.2. Two-dimensional motor Time averaged values of local rb along the axial length of a 2D slab motor and comparison with numerical simulations were reported in [15]. Here we discuss prediction of one representative experiment (test no. 11) with the present model. More detailed discussion on this is presented in [13].
C.P. Kumar, A. Kumar / Aerospace Science and Technology 29 (2013) 92–99
95
Fig. 2. (A) A comparison of axial and swirl velocity profiles predicted using present model with the experiments [25] in a long cylindrical pipe. (B) A comparison of predicted local regression rate with the experiments [15] for a two-dimensional slab motor.
Simulation was first carried out for a straight port at average port height based on half burn time [13]. Fig. 2(B) shows that the model under predicts rb and that the local rb increased more steeply than observed in the experiment. To make simulation more realistic an inclination was prescribed to the fuel surface based on the local rb at initial port height to account for non-uniform regression along the length of the grain [13]. The predicted local rb is seen to match well with experiment values over a large portion of fuel grain length. The mismatch near the head end of the motor is attributed to flat fuel surface assumption [13]. 4. Results and discussions Simulations were carried out primarily on a laboratory scale single port motor configuration which is easily amenable to experiments. A number of simulations were carried out to understand the influence of swirling flow on local and global regression rates, these include varying of inlet S w (0 to 2.5), variation in swirl profile for a fixed representative inlet S w (here 1) and varying the grain length (150 to 750 mm). Selected computations were also carried out on large diameter motors (of diameter up to 150 mm) at inlet S w = 1 and L / D = 25 to identify the effect of scaling. 4.1. Laboratory scale motor A laboratory scale hybrid rocket motor with port diameter of 30 mm and fuel grain length of 150 mm was taken as the base motor. The nozzle throat diameter was 14 mm and exit diameter was equal to the port diameter. Although the throat dimension of the nozzle affects the chamber pressure, the pressure is known to have negligible effect on rb of the hybrid rocket motor [18]. The average inlet oxidizer mass flux was maintained at 132 kg/m2 /s. The mesh comprised of quadrilateral grids which are clustered (with expansion ratio of about 1.1) towards fuel surface and oxidizer inlet to resolve the steep gradients and flame stabilization point accurately. Regression rate is dictated by the heat flux at the fuel surface. To ensure grid independence of the heat flux prediction, the first cell size from the fuel surface should be fine enough so as to capture the temperature gradient accurately. Grid independence study was conducted by reducing the size of smallest grid by an order till the difference between predicted average regression rates for two finest meshes was within 5%. A decrease in the first cell height from 5E−3 mm to 5E−4 mm, decreased regression rate by 10.43%. A further decrease in the cell height from 5E−4 mm to 5E−5 mm decreased regression rate by 4.86% (< 5%). Hence the
Fig. 3. Effect of swirl on the local heat transport at the interface along the length of the motor. The inset shows the corresponding local regression rate.
smallest cells of size 0.0005 mm (high) × 0.05 mm (wide) were used in flame stabilization zone. The largest cells of size 0.2 mm × 1 mm occupied most part of the combustion chamber. The grid size used in the combustion chamber was 176 × 111. The predicted average rb was found to be independent of grid. Grids for parametric study on other motors were decided such that the coarsest and the finest cell sizes were same as that of the base motor. 4.1.1. Effect of swirl on local heat transfer to the fuel and regression rate Since heat transport process primarily governs rb , we first look at the effect of swirl on the heat transfer and consequently the rb . Fig. 3 presents axial variation of components of heat flux at the fuel surface. The inset figure shows the local rb . The solid curves are for the case with swirl (inlet S w = 1) and the broken line curves present the base case without swirl. The gas phase convective transport, qconvection is the mechanism of heat transfer to the fuel and so determines the local rb . For the base case when there is no swirl in the flow (broken line curves), the local convective heat flux and therefore rb decreases initially and then increases along the axial length of the fuel grain, a trend observed previously in [15]. This occurs due to two competing phenomena. Near the head end, the flame is very close to the fuel surface and the heat flux is high but further downstream pyrolyzed fuel addition pushes the flame away from the fuel surface and the heat flux decreases. However, progressive addition of pyrolyzed fuel increases
96
C.P. Kumar, A. Kumar / Aerospace Science and Technology 29 (2013) 92–99
Fig. 4. Radial variation of velocity magnitude at various axial locations for the base case (left) and for the case with swirl (right).
Fig. 6. Effect of inlet swirl on the performance parameters. The inset shows local regression rate along the length of the motor.
4.1.2. Effect of inlet swirl number The inlet S w was varied by varying the inlet tangential velocity while maintaining the average inlet mass flux at 132 kg/m2 s. The variation of performance parameters with inlet S w is plotted in Fig. 6. The local rb along the length of the grain for selected inlet S w is illustrated as inset in Fig. 6. The local and average rb increase with increase in S w because of enhanced heat transfer to the fuel as discussed before. Consequently the motor that operated at oxidizer rich O / F ratio for low inlet S w approach near stoichiometric value at high S w . The thrust and the C ∗ also increase following rb at a decreasing rate with increase in S w . To carry out other parametric variations an inlet S w of reasonable magnitude (here S w = 1) was chosen. This choice will not affect the trend reported in this work as rb monotonically increases with increase in swirl strength. Fig. 5. Top: Temperature contours. Bottom: CO2 mass fraction, Y CO2 , contours for the base case (upper half) and for the case with swirl (lower half).
the local mass flux and thus the heat flux to the fuel surface. While the former effect is dominant near the head end, the latter effect is dominant downstream. For the case with swirl (solid line curves) the local rb follows the decreasing–increasing trend similar to the base case, however exhibit higher values (similar increase reported in [17]) and the minimum is shifted closer to the head end. As expected the heat flux components at interface are enhanced in presence of swirl. The reason for increased heat fluxes is explained next in Fig. 4 and Fig. 5. Velocity magnitude along the radius for the base and the swirl cases at four different axial locations are compared in Fig. 4. Swirl increases the magnitude of velocity in the combustion chamber. This results in larger velocity gradients at the fuel interface and thus increases heat and mass transport at the fuel surface. Another consequence of swirl in the flow is shown in Fig. 5. Fig. 5 presents contours of temperature (top) and CO2 mass fraction (bottom) for both the base case and swirl case (upper and lower half respectively in each figure). One can note that the flame (indicated by temperature contour 3000 K) shifts only slightly away from the fuel surface for the case with swirl. However, the core of the combustion chamber for the case with swirl gets heated much before the tail end. The products for the swirl case (indicated by CO2 mass fraction) also reach the core before the tail end. Thus introduction of swirl in flow also increases radial diffusive transport.
4.1.3. Effect of velocity profiles at the inlet Both swirl and axial velocity profiles at the inlet can vary depending upon the injector design. Therefore, computations were carried out to understand the influence of velocity profiles on rb . For this purpose two different axial velocity profiles viz., a constant velocity (plug flow) and a parabolic velocity profile as given by Eq. (19a) were chosen. The swirl velocity profile was prescribed by Eq. (20).
vz =
−1 3
r R0
n
fz 2 r − R 0r
μ
(20)
The value of the exponent n is chosen such that the peak swirl velocity is achieved at desired radial location (here three such locations are considered, 0.2R 0 , 0.5R 0 and 0.8R 0 ). In addition a radially linearly increasing swirl velocity profile (representing solid body rotation) was also investigated. Various combinations of above mentioned inlet axial and swirl velocity profiles at inlet S w = 1 were simulated. The average rb for these cases are shown in Fig. 7. One can note that the average rb is not affected axial velocity profile whereas it is sensitive to changes in swirl velocity profile. The swirl velocity profile with peak closest to the axis (here 0.2R 0 ) predicted highest average rb where as it was least for the solid body rotation profile. For completeness Fig. 7 also presents comparison of average rb over a range of S w , for two swirl profiles, one with peak at 0.5R 0 and other with solid body rotation. The inset in Fig. 7 presents the ratio of average rb for these two swirl cases. One can note that effect of swirl velocity profile is maximum at about S w = 1 which is seen to result in variation
C.P. Kumar, A. Kumar / Aerospace Science and Technology 29 (2013) 92–99
97
is higher when peak swirl is closer to axis we look at flow and temperature fields for selected swirl velocity profiles. Fig. 8 presents streamlines (upper half of the figures) and temperature contours (lower half) for simulations with parabolic axial velocity profile and swirl velocity profiles with peak at 0.2R 0 (top), 0.5R 0 (middle) and 0.8R 0 (bottom) for inlet S w = 1. The centrifugal force term in the radial momentum equation (ρ v 2z /r ) is inversely proportional to radius. With nearly same order for the magnitude of swirl velocity (not shown in the figure), the centrifugal force term is highest for peak location at 0.2R 0 . This results in larger radial velocity towards fuel surface at the head end. The temperature contours also reveal enhanced heat and mass transport from the reaction zone to the motor core. A closer look at the flame (here temperature contour of 3000 K) shows that flame is closest to the fuel surface for the case of peak swirl at 0.2R 0 where near-wall swirl velocity is small. As a result average rb is higher for the case peak swirl velocity at 0.2R 0 . Fig. 7. Effect of inlet velocity (swirl and axial) profiles at various swirl numbers ( S w ) on the average regression rate.
Fig. 8. Effect of inlet swirl velocity profile for S w = 1.0 on streamlines (drawn in upper half of subfigures) and gas temperature distribution (lower half). The peak inlet swirl velocity is located at 0.2R 0 (top), 0.5R 0 (middle), and 0.8R 0 (bottom).
in average rb by as much as 50%. At very small S w below 0.4, this variation is negligible. At high S w above 1.5, the swirl momentum becomes much higher than the axial momentum. Thus the effect of swirl velocity profile is less than 20%. To understand why rb
4.1.4. Effect of length of the combustion chamber Here we explore effectiveness of swirl in improving rb with fuel grain length. Fig. 9(A) shows the variation in local rb along the length of fuel grain of lengths (150 mm, 450 mm and 750 mm) at S w = 1. Figure also shows the base case rb profile represented by the solid line with symbols for motor with 750 mm long fuel grain ( L / D = 25). For the base case, the local rb peaks at around 500 mm and then decreases. This peak is because at this location all the oxygen is depleted and therefore, further downstream there is no more heat release and the bulk temperature drops due to dilution by fresh fuel addition. Consequently heat flux to the fuel surface and hence rb drops. For simulations with swirl, the local rb near the head end is higher compared to the base case. Hence the available oxygen is consumed at a shorter length (∼300 mm) and the peak rb also occurs at about this location. Interestingly the local rb for the short motors are identical to that of the long motor. One can note that the local regression rate curves for the base case and the case with swirl crossover. This tells that swirl is not effective in enhancing rb beyond the cross over grain length (here at about 450 mm). This is depicted in the inset of figure as the ratio of local rb of the swirl case to the base case. Fig. 9(B) shows comparison of performance parameters for the base case (dashed lines) with the swirl case (solid lines) for various fuel grain lengths. Longer motors have larger fuel surface area for a fixed oxidizer mass flow rate, hence have lower overall O / F
Fig. 9. Effectiveness of swirl with motor length. A comparison of (A) local regression rate and (B) performance parameters between the case with swirl ( S w = 1) and the base case for motors of various lengths.
98
C.P. Kumar, A. Kumar / Aerospace Science and Technology 29 (2013) 92–99
Fig. 10. Effect of diameter on the local regression rate along the length of the motor for motors with two different diameters at L / D = 25 for non-swirling inlet flow and swirling inlet flow with swirl number, S w = 1.
ratio. Owing to increased rb , thrust increases with grain length. However, the average rb exhibits a maxima as explained above. C ∗ strongly depends on the gas temperature at the nozzle inlet, therefore, a maximum is attained at a grain length required for complete oxygen depletion at the nozzle inlet (450 mm long fuel grain for base case compared to 300 mm long fuel grain for swirl case). In both the cases, maximum C ∗ was attained at an O / F ratio of 2 much lesser than stoichiometric value of 3.25. This is caused by presence of un-burnt fuel between the reaction zone and the fuel surface. This also results in low combustion efficiency of about 70% which can be improved if a post combustion chamber is used.
Fig. 11. A comparison of temperature contours in combustion chamber between non-swirling inlet flow (top) and swirling inlet flow (bottom) with S w = 1 for the diameters of 30 mm (upper half) and 150 mm (lower half). The L / D is fixed at 25.
4.2. Large scale motors To complete this study we finally explore effectiveness of swirl in large port diameter grain. Simulations were carried out for selected grain port diameters (D = 30, 45, 90, and 150 mm) at inlet swirl numbers, S w = 0 and 1 and average inlet oxidizer mass flux of 132 kg/m2 /s. Geometrical similarity was maintained with the long motor of L / D = 25 to include the effect of grain length. Fig. 10 presents local rb along the length of the grain for port diameters of 30 mm and 150 mm, for inlet S w = 0 and 1. One can note that rb is lower for simulation with large port diameter. This is because of less fuel vapor added per unit grain port area for the large diameter grain. Another consequence of large port diameter is that the core gets heated at a distance further downstream compared to the small port diameter grain as illustrated by temperature contours in Fig. 11. The enhancement of average regression for various X / D is plotted in Fig. 12. The enhancement over the no-swirl case is seen to be higher for fuel grains with large port diameter. This is due to lower values of average rb for the no-swirl case. It is also noted that higher enhancement in average rb is achievable for short motors ( L / D < 10). This is due to the decay of swirl in the combustion chamber. The variation of local swirl number along the length of the grain is plotted in the inset of Fig. 12. The swirl strength is seen to decay by about 90% at L / D of 10. Thus swirl can be effectively used to increase rb in short motors especially for large port diameters. 5. Conclusions The effect of swirling oxidizer flow at the inlet on the rb in hybrid rocket motors was numerically studied by varying inlet S w , inlet velocity profiles and fuel grain length. Effect of grain scal-
Fig. 12. Enhancement in the average regression rate (in percentage) along the length of motor for large port diameters and swirling flow. The inset figure shows axial variation of local swirl strength ( S w ).
ing was also investigated. The following conclusions may be drawn from this study. 1. The presence of swirl velocity component increases the velocity magnitude in the flow field and causes increased velocity gradient at the fuel surface. Heat and mass transport at the fuel surface increases which consequently increases the rb . 2. Increase in inlet S w increases average rb monotonically. 3. Variation in axial velocity profile did not affect average rb . However, a swirl velocity profile with peak near the axis of the motor results in maximum average rb because of enhanced transport processes and flame positioned closer to the fuel surface. The contribution of swirl velocity profile in increasing average rb is only a fraction of contribution of swirl itself, the maximum contribution being about 50% at about inlet S w = 1. 4. For long fuel grains, the local rb and C ∗ for base case attain a maximum value at a location where all the available oxidizer is consumed and the flame reaches the axis (occurs at L / D ∼ 15 for D = 30 mm). With swirling oxidizer flow, this maximum is achieved in shorter L / D (here 10) due to increased rb . 5. The local and average rb decrease with increase in grain port diameter (swirl or no-swirl) due to lower increase in local mass flux by addition of pyrolyzed fuel. Interestingly, the
C.P. Kumar, A. Kumar / Aerospace Science and Technology 29 (2013) 92–99
percentage enhancement in average rb from swirl is higher for large diameter grains. 6. The effect of swirl is high near the head end and reduces downstream. The swirl strength was seen to drop to about 5% at L / D = 10. References [1] D. Altman, A. Holzman, Overview and history of hybrid rocket propulsion, in: M.J. Chiaverini, K.K. Kuo (Eds.), Fundamentals of Hybrid Rocket Combustion and Propulsion, in: AIAA Progress in Astronautics & Aeronautics, vol. 218, AIAA, Reston, VA, 2007, pp. 1–36. [2] ANSYS-Fluent Software Release Version No. 6.3.26, 2006. [3] Y.-S. Chen, T.H. Chou, B.R. Gu, J.S. Wu, B. Wu, Y.Y. Lian, L. Yang, Multiphysics simulations of rocket engine combustion, Computers and Fluids 45 (1) (2011) 29–36. [4] G.C. Cheng, R.C. Farmer, H.S. Jones, J.S. McFarlane, Numerical simulation of internal ballistics of a hybrid rocket motor, AIAA Paper 1994-0554, 1994. [5] G. Gariani, F. Maggi, L. Galfetti, Numerical simulation of HTPB combustion in a 2D hybrid slab combustor, Acta Astronautica 69 (2011) 289–296. [6] S. Gordon, B.J. McBride, Computer program for calculation of complex chemical equilibrium compositions and applications, NASA Reference Publication 1311, 1996. [7] X. Jiang, G.A. Siamas, L.C. Wrobel, Analytical equilibrium swirling inflow conditions for computational fluid dynamics, Technical notes, AIAA Journal 46 (4) (2008) 1015–1019. [8] H.J. Kim, Y.M. Kim, Numerical modeling of combustion processes of hybrid rocket engine, in: 37th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, AIAA Paper 2001-4504, 2001. [9] K. Kitagawa, T. Mitsutani, T. Ro, S. Yuasa, Effects of swirling liquid oxygen flow on combustion of a hybrid rocket engine, in: 40th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, AIAA Paper 2004-3479, 2004. [10] W.H. Knuth, M.J. Chiaverini, Solid-fuel regression rate behavior of vortex hybrid rocket engines, Journal of Propulsion and Power 18 (3) (2002) 600–609. [11] M. Kumar, A.F. Ghoniem, Multiphysics simulations of entrained flow gasification. Part I: Validating the non-reacting flow solver and the particle turbulent dispersion model, Energy and Fuels 26 (2012) 451–463. [12] C. Kumar, A. Kumar, On the role of radiation and dimensionality in predicting flow opposed flame spread over thin fuels, Combustion Theory and Modeling 16 (3) (2012) 537–569.
99
[13] C.P. Kumar, A. Kumar, Effect of diaphragms on regression rate in hybrid rocket motors, Journal of Propulsion and Power (2013), in press, http://dx.doi.org/ 10.2514/1.B34671. [14] K.K. Kuo, M.J. Chiaverini, Challenges of hybrid rocket propulsion in the 21st century, in: M.J. Chiaverini, K.K. Kuo (Eds.), Fundamentals of Hybrid Rocket Combustion and Propulsion, in: AIAA Progress in Astronautics & Aeronautics, vol. 218, AIAA, Reston, VA, 2007, pp. 593–638. [15] K.K. Kuo, Y.-C. Lu, M.J. Chiaverini, D.K. Johnson, N. Serin, G.A. Risha, C.L. Merkle, S. Venkateshwaran, Fundamental phenomena on fuel decomposition and boundary layer combustion processes with applications to hybrid rocket motors, Semi-Annual Progress Report, NASA-CR-201843, 1996. [16] C. Lee, Y. Na, J.W. Lee, Y.H. Byun, Effect of induced swirl flow on regression rate of hybrid rocket fuel by helical grain configuration, Aerospace Science and Technology 11 (2007) 68–76. [17] T.-S. Lee, A. Potapkin, The performance of a hybrid rocket with swirling GOx injection, in: XI International Conference on the Methods of Aerophysical Research, 2002. [18] A. Lewin, J. Dennis, B. Conley, D. Suzuku, Experimental determination of performance parameters for a polybutadiene/oxygen hybrid rocket, in: 28th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, AIAA Paper 923590, 1992. [19] G.A. Marxman, Combustion in turbulent boundary layer on a vaporizing surface, in: Tenth International Symposium on Combustion, 1965, pp. 1337–1349. [20] F.R. Menter, Influence of free stream values on k–ω turbulence model predictions, AIAA Journal 30 (6) (1992) 1657–1659. [21] F.R. Menter, Two equation eddy-viscosity turbulence models for engineering applications, AIAA Journal 32 (8) (1994) 1598–1605. [22] A.V. Potapkin, T.S. Lee, Experimental study of thrust performance of a hybrid rocket motor with various methods of oxidiser injection, Combustion Explosion and Shock Waves 40 (4) (2004) 386–392. [23] P.A. Ramakrishna, P.J. Paul, H.S. Mukunda, Sandwich propellant combustion: Modeling and experimental comparison, Proceedings of the Combustion Institute 29 (2) (2002) 2963–2973. [24] K.H. Shin, C. Lee, Y.H. Yu, The experiments for the enhancement of regression rate of hybrid rocket fuel, KSME International 19 (10) (2005) 1926–1936. [25] W. Steenbergen, J. Voskamp, The rate of decay of swirl in turbulent pipe flow, Flow Measurement and Instrumentation 9 (2) (1998) 67–78. [26] E.P. Volchkov, V.I. Terekhov, V.V. Terekhov, Flow structure and heat and mass transfer in boundary layers with injection of chemically reacting substances (Review), Combustion Explosion and Shock Waves 40 (1) (2004) 1– 16.
Contents lists available at SciVerse ScienceDirect
Aerospace Science and Technology www.elsevier.com/locate/aescte
Effect of swirl on the regression rate in hybrid rocket motors C. Palani Kumar 1 , Amit Kumar ∗,2 National Centre for Combustion Research and Development, Department of Aerospace Engineering, Indian Institute of Technology Madras, Chennai, Tamil Nadu-600036, India
a r t i c l e
i n f o
Article history: Received 7 July 2012 Received in revised form 29 January 2013 Accepted 30 January 2013 Available online 13 February 2013 Keywords: Hybrid rocket Regression rate Swirl flow Numerical simulation
a b s t r a c t In this work, the effect of swirl on regression rate (rb ) in a hybrid rocket motor is investigated numerically. The rb increased monotonically with inlet swirl number and was also found to depend on the inlet swirl profile. The swirl velocity profiles with the peak closer to the axis yielded higher rb . Parametric study on fuel grains of various lengths and diameters ( L / D 25) shows that swirl is more effective in improving the average rb for short grains ( L / D < 5) and large diameter grains. © 2013 Elsevier Masson SAS. All rights reserved.
1. Introduction Hybrid rockets have many distinct advantages over conventional chemical rockets, including simplicity, fuel costs, hardware costs, and safety. However, one concern with the operation of these rockets is low rb of the solid fuel [14,15]. While multi-port complex grain design can be used for obtaining reasonable thrust with low regression rates, this results in larger residuals and low fuel loading. Further the grain integrity also becomes a factor of concern [1]. Extensive research over the past few decades on hybrid rocket motors has contributed greatly towards the fundamental understanding of its working. One of the rb enhancement techniques established uses swirling flow inside the combustion chamber [9, 10,16,22,24]. Swirling flow in the combustion chamber can be achieved by either tangentially injecting oxidizer from the tail end [10], or using a swirl type injector [9,16,22,24] or using helical grain configuration of solid fuel [16,24]. Although rb improvements are reported in the literature [9,10, 16,22,24], the amount of swirl is most often not quantified. There is a lack of basic understanding of the effect of swirling flow on rb of the hybrid rocket motor. The present numerical study aims to systematically understand the underlying physics and relate rb to swirl strength, quantified by a non-dimensional swirl number.
* 1 2
Corresponding author. Tel.: +91 44 22574019; fax: +91 44 22574002. E-mail address: [email protected] (A. Kumar). Ph.D. Research Scholar. Associate Professor.
1270-9638/$ – see front matter © 2013 Elsevier Masson SAS. All rights reserved. http://dx.doi.org/10.1016/j.ast.2013.01.011
2. Numerical model and solution 2.1. Geometrical configuration and computational domain Fig. 1 presents the schematic of the hybrid rocket combustion chamber. Invariance in azimuthal direction is assumed to reduce the computational domain to a 2D axi-symmetric configuration. The shaded region in Fig. 1 is the computational domain. The domain consists of a combustion chamber and a CD nozzle. The dimensions of combustion chamber, nozzle and the details on grid are discussed later in Section 4. 2.2. Governing equations The processes occurring inside the hybrid rocket combustion chamber can be adequately described by basic flow equations of continuity, momentum, energy and species. The governing equations in 2D cylindrical coordinates are summarized below. Continuity equation
∂ ∂ ∂ ρ vr (ρ ) + (ρ v x ) + (ρ v r ) + =0 ∂t ∂x ∂r r
(1)
Axial momentum equation
∂ 1 ∂ 1 ∂ (ρ v x ) + (r ρ v x v x ) + (r ρ v r v x ) ∂t r ∂x r ∂r ∂p 1 ∂ ∂ vx 2 =− + rμ 2 − (∇ · v ) ∂x r ∂x ∂x 3 ∂ v x ∂ vr 1 ∂ rμ + + r ∂r ∂r ∂x
(2)
C.P. Kumar, A. Kumar / Aerospace Science and Technology 29 (2013) 92–99
93
Nomenclature A, As C C∗ D i ,m Di , D0 Dω Dt E a , E as G˜ k , G ω h, hi HR, HP J j k kr L p rb
Arrhenius pre-exponential factor for gas phase reaction and solid fuel pyrolysis Concentration of species Characteristic exhaust velocity . . . . . . . . . . . . . . . . . . . . m/s Diffusivity of species i in the mixture . . . . . . . . . . . m2 /s Initial and final port diameter . . . . . . . . . . . . . . . . . . . . . . . m Cross diffusion term Throat diameter . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m Arrhenius activation energy for gas phase reaction, solid fuel pyrolysis Production of k and ω Sensible enthalpy of mixture and species i . . . . kJ/kg/s Heat of gas phase reaction and solid fuel pyrolysis Diffusion flux of species j Turbulent kinetic energy Reaction rate constant Fuel grain length . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m Pressure . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . N/m2 Regression rate . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . m/s
R0 Ru Sct T s, T∞ v x, vr , W , Wi x, r , z Yi Yk, Yω
α Γk , Γω λeff λs μ, μt
ρ ω
Radius of the combustion chamber . . . . . . . . . . . . . . . . . m Universal gas constant Turbulent Schmidt number Surface temperature and initial temperature of fuel . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . K v z Axial, radial and swirl (tangential) velocity . . . . . m/s Global reaction rate and consumption rate of species i . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . mol/m3 /s Cylindrical coordinates in axial radial and tangential direction Mass fraction of species i Dissipation of k and ω Thermal diffusivity Effective diffusivities of k and ω Effective thermal conductivity Thermal conductivity of the solid fuel Molecular and turbulent viscosities . . . . . . . . . . . . N s/m2 Density . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . kg/m3 Specific dissipation rate
The detailed turbulence models like RSM require rigorous closure strategies and grid quality requirements and involve complexities in computation and convergence [11]. Hence these are not suited for the parametric study intended here. Simpler models like lowRe k–ε turbulence models are widely popular in literature [3–5,15]. However, under intense injection the k–ε model fails in qualitative flow prediction [26]. Although the standard k–ω model overcomes this deficiency, it is sensitive to the free stream values of ω [20]. Hence, SST k–ω model was developed [21] which combines the advantages of both k–ε (free stream accuracy) and standard k–ω (near-wall accuracy). Thus SST k–ω was used to predict turbulence in the present study. Fig. 1. Top: Schematic of a single port cylindrical hybrid rocket motor showing solid fuel grain, combustion chamber and a convergent–divergent nozzle. The shaded region is the computational domain. Bottom: Schematic of the computational domain and the boundary types.
Radial momentum equation
∂ 1 ∂ 1 ∂ (ρ v r ) + (r ρ v x v r ) + (r ρ v r v r ) ∂t r ∂x r ∂r ∂p 1 ∂ ∂ v x ∂ vr =− + rμ + ∂r r ∂x ∂r ∂x 1 ∂ ∂ vr 2 + rμ 2 − (∇ · v ) r ∂r ∂r 3 − 2μ
vr r2
+
2μ 3 r
(∇ · v ) + ρ
v 2z r
(3)
∂ (ρ E ) + ∇ · v (ρ E + p ) ∂t h j J j + (τ¯¯ · v ) + H R = ∇ · λeff ∇ T −
(7)
(8)
(4)
Here E is defined as
E =h−
Tangential momentum equation
∂ 1 ∂ 1 ∂ (ρ v z ) + (r ρ v x v z ) + (r ρ v r v z ) ∂t r ∂x r ∂r 1 ∂ 1 ∂ 3 ∂ vz vr v z ∂ vz = + 2 −ρ rμ r μ r ∂x ∂x ∂r r r r ∂r
(6)
j
In the above equation the gradient of velocity vector is defined
∂ v x ∂ vr vr + + ∂x ∂r r
∂ ∂ ∂ ∂k Γk + G˜ k − Y k (ρ k) + (ρ ku i ) = ∂t ∂ xi ∂xj ∂xj ∂ ∂ ∂ ∂ω Γω + Gω − Yω + Dω (ρω) + (ρωu i ) = ∂t ∂ xi ∂xj ∂xj Energy equation
as
(∇ · v ) =
SST k–ω equations are given by
p
ρ
+
v2 2
(9)
Species transport equation
(5)
The flow in the combustion chamber and nozzle is likely to be turbulent, therefore, an appropriate turbulence model is required.
∂ (ρ Y i ) + ∇ · (ρ v Y i ) = −∇ · ( J i ) + W i ∂t
(10)
Here,
J i = −(ρ D j + μt /Sct )∇ Y j
(11)
94
C.P. Kumar, A. Kumar / Aerospace Science and Technology 29 (2013) 92–99
2.3. Gas phase chemistry model Combustion model with detailed chemistry is desired for accurate predictions. However, proper reaction mechanism and related kinetics are not readily available in the literature. Nevertheless, thermodynamic effect of multi-species can be readily incorporated [6]. Therefore, a global chemistry [8] is adopted here where heat of combustion is tuned to be consistent with thermodynamics following [12]. This is expected to give reasonable prediction of temperature field and thus heat transfer to the solid fuel. As butadiene is the main pyrolysis product of HTPB, the single step global reaction is taken as
C4 H6 + 5.5O2 → 4CO2 + 3H2 O
(12)
The reaction rate W is given by assuming a second order reaction as
W = kr C fuel C oxidizer = kr C C4 H6 C O2
(13)
kr can be expressed by Arrhenius law as kr = A (exp(− E a / R u / T )). Here, A = 8.8E+11 and E a = 1.2637E+08 J/kmol [4].
( S w ). Swirl number is defined as the ratio of axial fluxes of angular momentum to axial momentum, non-dimensionalized by inlet radius [7].
Sw =
ρ v x v z r 2 dr R 0 ρ v 2x r dr
(18)
Jiang et al. [7] also provided the analytical axial and swirl velocity profiles which can be considered as equivalent for swirling flow inlet boundary conditions in CFD as
vx =
−1 f x
r 2 − R 20
4 μ −1 f z 2 vz = r − R 0r 3 μ
(19a) (19b)
Substituting Eqs. (19a) and (19b) in Eq. (18) results in S w as a function of body force terms f x and f z . For a chosen v x ( f x = 1), v z is varied (different f z ) to achieve different inlet S w . The nozzle walls were considered adiabatic. Atmospheric pressure was specified at the nozzle exit. 2.6. Numerical method of solution
2.4. Pyrolysis model The pyrolysis of the solid fuel (HTPB grain) was modeled by the zeroth order Arrhenius equation following [23]. The mass of fuel released by pyrolysis is given by
˙ = ρs rb = A s exp(− E as / R u T s ) m
(14)
2.5. Boundary conditions The boundary types are mentioned in Fig. 1. The fuel inlet was defined as an interface where the mass and energy balance was applied. The reduced energy balance equation is given by
λeff (∂ T /∂ y ) g − λs (∂ T /∂ y )s + ρs rb ( H P ) = 0
(15)
The mass balance at the interface is given by
˙ = ρs rb = −ρ v m
(16)
A part of the thermal energy released in the reaction zone is convected to the solid–gas interface (1st term in Eq. (15); qconvection ). A part of this thermal energy is then transferred into solid fuel by means of conduction in the solid phase (2nd term; qconduction loss ). The remaining energy is utilized for solid fuel pyrolysis (3rd term; qpyrolysis ). The heat lost inside the solid fuel was determined by an analytical solution of 1D heat conduction inside the solid with moving hot boundary. The analytical temperature profile inside the solid fuel is given by
T ( y ) = ( T s − T ∞ ) exp(−rb y /α ) + T ∞
(17)
The radiation heat transfer is neglected in the present study. Simulations [12,15] have shown that the effect of radiation (compared to simulations without radiation) on change in net heat flux incident on the solid fuel surface is small (typically about 10%). This is because the radiation heat flux on the fuel surface increases the blowing velocity when radiation is included. This in turn reduces qconvection . Radiation from flame also reduces flame temperature which further reduces qconvection . The net result is only a marginal increase in net heat flux on the fuel surface (or a small increase in rb [19]). A fixed average mass flux was specified at the oxidizer inlet in all simulations. Swirl at the inlet was quantified by swirl number
The governing equations along with the boundary conditions were solved numerically using a pressure based, double precision, unsteady solver [2]. Second order upwind scheme was used for spatial discretization of the convective terms whereas second order central scheme was used for all other terms in the transport equations. The second order implicit formulation was used for temporal discretization. First order discretization schemes are not used because they are less accurate. Second order discretization schemes were opted as the velocity profiles and the local rb along the axis did not change appreciably between the second and third order discretization schemes. The discretized governing equations were iterated with the time step of 0.0002 s and 50 iterations per time step. The converged solution for short motors was attained typically in about 1000 time steps. Convergence was assessed by the condition that the scaled residuals reduce below a value of 1E−5 for all transport equations and become invariant with time. 3. Validation of numerical model Detailed experimental data on swirling reacting flow inside a hybrid rocket motor is not available in the literature. Hence the present numerical model was validated with cold turbulent swirling flow in a pipe [25] and 2D hybrid rocket experiments and simulations [15]. 3.1. Swirl flow in a pipe Steenbergen and Voskamp [25] reported experimental measurements on turbulent swirl velocity profiles at various axial locations (x/ D = 4.4, 7.6 and 14.5). Since the inlet velocity profiles were not known, the reported velocity profiles at x/ D = 4.4 were provided at the inlet in the present simulation. The computed profiles at the sections x/ D = 7.6 and 14.5 are compared with the experimental values in Fig. 2(A). The predicted velocity profiles match well with the measured profiles. 3.2. Two-dimensional motor Time averaged values of local rb along the axial length of a 2D slab motor and comparison with numerical simulations were reported in [15]. Here we discuss prediction of one representative experiment (test no. 11) with the present model. More detailed discussion on this is presented in [13].
C.P. Kumar, A. Kumar / Aerospace Science and Technology 29 (2013) 92–99
95
Fig. 2. (A) A comparison of axial and swirl velocity profiles predicted using present model with the experiments [25] in a long cylindrical pipe. (B) A comparison of predicted local regression rate with the experiments [15] for a two-dimensional slab motor.
Simulation was first carried out for a straight port at average port height based on half burn time [13]. Fig. 2(B) shows that the model under predicts rb and that the local rb increased more steeply than observed in the experiment. To make simulation more realistic an inclination was prescribed to the fuel surface based on the local rb at initial port height to account for non-uniform regression along the length of the grain [13]. The predicted local rb is seen to match well with experiment values over a large portion of fuel grain length. The mismatch near the head end of the motor is attributed to flat fuel surface assumption [13]. 4. Results and discussions Simulations were carried out primarily on a laboratory scale single port motor configuration which is easily amenable to experiments. A number of simulations were carried out to understand the influence of swirling flow on local and global regression rates, these include varying of inlet S w (0 to 2.5), variation in swirl profile for a fixed representative inlet S w (here 1) and varying the grain length (150 to 750 mm). Selected computations were also carried out on large diameter motors (of diameter up to 150 mm) at inlet S w = 1 and L / D = 25 to identify the effect of scaling. 4.1. Laboratory scale motor A laboratory scale hybrid rocket motor with port diameter of 30 mm and fuel grain length of 150 mm was taken as the base motor. The nozzle throat diameter was 14 mm and exit diameter was equal to the port diameter. Although the throat dimension of the nozzle affects the chamber pressure, the pressure is known to have negligible effect on rb of the hybrid rocket motor [18]. The average inlet oxidizer mass flux was maintained at 132 kg/m2 /s. The mesh comprised of quadrilateral grids which are clustered (with expansion ratio of about 1.1) towards fuel surface and oxidizer inlet to resolve the steep gradients and flame stabilization point accurately. Regression rate is dictated by the heat flux at the fuel surface. To ensure grid independence of the heat flux prediction, the first cell size from the fuel surface should be fine enough so as to capture the temperature gradient accurately. Grid independence study was conducted by reducing the size of smallest grid by an order till the difference between predicted average regression rates for two finest meshes was within 5%. A decrease in the first cell height from 5E−3 mm to 5E−4 mm, decreased regression rate by 10.43%. A further decrease in the cell height from 5E−4 mm to 5E−5 mm decreased regression rate by 4.86% (< 5%). Hence the
Fig. 3. Effect of swirl on the local heat transport at the interface along the length of the motor. The inset shows the corresponding local regression rate.
smallest cells of size 0.0005 mm (high) × 0.05 mm (wide) were used in flame stabilization zone. The largest cells of size 0.2 mm × 1 mm occupied most part of the combustion chamber. The grid size used in the combustion chamber was 176 × 111. The predicted average rb was found to be independent of grid. Grids for parametric study on other motors were decided such that the coarsest and the finest cell sizes were same as that of the base motor. 4.1.1. Effect of swirl on local heat transfer to the fuel and regression rate Since heat transport process primarily governs rb , we first look at the effect of swirl on the heat transfer and consequently the rb . Fig. 3 presents axial variation of components of heat flux at the fuel surface. The inset figure shows the local rb . The solid curves are for the case with swirl (inlet S w = 1) and the broken line curves present the base case without swirl. The gas phase convective transport, qconvection is the mechanism of heat transfer to the fuel and so determines the local rb . For the base case when there is no swirl in the flow (broken line curves), the local convective heat flux and therefore rb decreases initially and then increases along the axial length of the fuel grain, a trend observed previously in [15]. This occurs due to two competing phenomena. Near the head end, the flame is very close to the fuel surface and the heat flux is high but further downstream pyrolyzed fuel addition pushes the flame away from the fuel surface and the heat flux decreases. However, progressive addition of pyrolyzed fuel increases
96
C.P. Kumar, A. Kumar / Aerospace Science and Technology 29 (2013) 92–99
Fig. 4. Radial variation of velocity magnitude at various axial locations for the base case (left) and for the case with swirl (right).
Fig. 6. Effect of inlet swirl on the performance parameters. The inset shows local regression rate along the length of the motor.
4.1.2. Effect of inlet swirl number The inlet S w was varied by varying the inlet tangential velocity while maintaining the average inlet mass flux at 132 kg/m2 s. The variation of performance parameters with inlet S w is plotted in Fig. 6. The local rb along the length of the grain for selected inlet S w is illustrated as inset in Fig. 6. The local and average rb increase with increase in S w because of enhanced heat transfer to the fuel as discussed before. Consequently the motor that operated at oxidizer rich O / F ratio for low inlet S w approach near stoichiometric value at high S w . The thrust and the C ∗ also increase following rb at a decreasing rate with increase in S w . To carry out other parametric variations an inlet S w of reasonable magnitude (here S w = 1) was chosen. This choice will not affect the trend reported in this work as rb monotonically increases with increase in swirl strength. Fig. 5. Top: Temperature contours. Bottom: CO2 mass fraction, Y CO2 , contours for the base case (upper half) and for the case with swirl (lower half).
the local mass flux and thus the heat flux to the fuel surface. While the former effect is dominant near the head end, the latter effect is dominant downstream. For the case with swirl (solid line curves) the local rb follows the decreasing–increasing trend similar to the base case, however exhibit higher values (similar increase reported in [17]) and the minimum is shifted closer to the head end. As expected the heat flux components at interface are enhanced in presence of swirl. The reason for increased heat fluxes is explained next in Fig. 4 and Fig. 5. Velocity magnitude along the radius for the base and the swirl cases at four different axial locations are compared in Fig. 4. Swirl increases the magnitude of velocity in the combustion chamber. This results in larger velocity gradients at the fuel interface and thus increases heat and mass transport at the fuel surface. Another consequence of swirl in the flow is shown in Fig. 5. Fig. 5 presents contours of temperature (top) and CO2 mass fraction (bottom) for both the base case and swirl case (upper and lower half respectively in each figure). One can note that the flame (indicated by temperature contour 3000 K) shifts only slightly away from the fuel surface for the case with swirl. However, the core of the combustion chamber for the case with swirl gets heated much before the tail end. The products for the swirl case (indicated by CO2 mass fraction) also reach the core before the tail end. Thus introduction of swirl in flow also increases radial diffusive transport.
4.1.3. Effect of velocity profiles at the inlet Both swirl and axial velocity profiles at the inlet can vary depending upon the injector design. Therefore, computations were carried out to understand the influence of velocity profiles on rb . For this purpose two different axial velocity profiles viz., a constant velocity (plug flow) and a parabolic velocity profile as given by Eq. (19a) were chosen. The swirl velocity profile was prescribed by Eq. (20).
vz =
−1 3
r R0
n
fz 2 r − R 0r
μ
(20)
The value of the exponent n is chosen such that the peak swirl velocity is achieved at desired radial location (here three such locations are considered, 0.2R 0 , 0.5R 0 and 0.8R 0 ). In addition a radially linearly increasing swirl velocity profile (representing solid body rotation) was also investigated. Various combinations of above mentioned inlet axial and swirl velocity profiles at inlet S w = 1 were simulated. The average rb for these cases are shown in Fig. 7. One can note that the average rb is not affected axial velocity profile whereas it is sensitive to changes in swirl velocity profile. The swirl velocity profile with peak closest to the axis (here 0.2R 0 ) predicted highest average rb where as it was least for the solid body rotation profile. For completeness Fig. 7 also presents comparison of average rb over a range of S w , for two swirl profiles, one with peak at 0.5R 0 and other with solid body rotation. The inset in Fig. 7 presents the ratio of average rb for these two swirl cases. One can note that effect of swirl velocity profile is maximum at about S w = 1 which is seen to result in variation
C.P. Kumar, A. Kumar / Aerospace Science and Technology 29 (2013) 92–99
97
is higher when peak swirl is closer to axis we look at flow and temperature fields for selected swirl velocity profiles. Fig. 8 presents streamlines (upper half of the figures) and temperature contours (lower half) for simulations with parabolic axial velocity profile and swirl velocity profiles with peak at 0.2R 0 (top), 0.5R 0 (middle) and 0.8R 0 (bottom) for inlet S w = 1. The centrifugal force term in the radial momentum equation (ρ v 2z /r ) is inversely proportional to radius. With nearly same order for the magnitude of swirl velocity (not shown in the figure), the centrifugal force term is highest for peak location at 0.2R 0 . This results in larger radial velocity towards fuel surface at the head end. The temperature contours also reveal enhanced heat and mass transport from the reaction zone to the motor core. A closer look at the flame (here temperature contour of 3000 K) shows that flame is closest to the fuel surface for the case of peak swirl at 0.2R 0 where near-wall swirl velocity is small. As a result average rb is higher for the case peak swirl velocity at 0.2R 0 . Fig. 7. Effect of inlet velocity (swirl and axial) profiles at various swirl numbers ( S w ) on the average regression rate.
Fig. 8. Effect of inlet swirl velocity profile for S w = 1.0 on streamlines (drawn in upper half of subfigures) and gas temperature distribution (lower half). The peak inlet swirl velocity is located at 0.2R 0 (top), 0.5R 0 (middle), and 0.8R 0 (bottom).
in average rb by as much as 50%. At very small S w below 0.4, this variation is negligible. At high S w above 1.5, the swirl momentum becomes much higher than the axial momentum. Thus the effect of swirl velocity profile is less than 20%. To understand why rb
4.1.4. Effect of length of the combustion chamber Here we explore effectiveness of swirl in improving rb with fuel grain length. Fig. 9(A) shows the variation in local rb along the length of fuel grain of lengths (150 mm, 450 mm and 750 mm) at S w = 1. Figure also shows the base case rb profile represented by the solid line with symbols for motor with 750 mm long fuel grain ( L / D = 25). For the base case, the local rb peaks at around 500 mm and then decreases. This peak is because at this location all the oxygen is depleted and therefore, further downstream there is no more heat release and the bulk temperature drops due to dilution by fresh fuel addition. Consequently heat flux to the fuel surface and hence rb drops. For simulations with swirl, the local rb near the head end is higher compared to the base case. Hence the available oxygen is consumed at a shorter length (∼300 mm) and the peak rb also occurs at about this location. Interestingly the local rb for the short motors are identical to that of the long motor. One can note that the local regression rate curves for the base case and the case with swirl crossover. This tells that swirl is not effective in enhancing rb beyond the cross over grain length (here at about 450 mm). This is depicted in the inset of figure as the ratio of local rb of the swirl case to the base case. Fig. 9(B) shows comparison of performance parameters for the base case (dashed lines) with the swirl case (solid lines) for various fuel grain lengths. Longer motors have larger fuel surface area for a fixed oxidizer mass flow rate, hence have lower overall O / F
Fig. 9. Effectiveness of swirl with motor length. A comparison of (A) local regression rate and (B) performance parameters between the case with swirl ( S w = 1) and the base case for motors of various lengths.
98
C.P. Kumar, A. Kumar / Aerospace Science and Technology 29 (2013) 92–99
Fig. 10. Effect of diameter on the local regression rate along the length of the motor for motors with two different diameters at L / D = 25 for non-swirling inlet flow and swirling inlet flow with swirl number, S w = 1.
ratio. Owing to increased rb , thrust increases with grain length. However, the average rb exhibits a maxima as explained above. C ∗ strongly depends on the gas temperature at the nozzle inlet, therefore, a maximum is attained at a grain length required for complete oxygen depletion at the nozzle inlet (450 mm long fuel grain for base case compared to 300 mm long fuel grain for swirl case). In both the cases, maximum C ∗ was attained at an O / F ratio of 2 much lesser than stoichiometric value of 3.25. This is caused by presence of un-burnt fuel between the reaction zone and the fuel surface. This also results in low combustion efficiency of about 70% which can be improved if a post combustion chamber is used.
Fig. 11. A comparison of temperature contours in combustion chamber between non-swirling inlet flow (top) and swirling inlet flow (bottom) with S w = 1 for the diameters of 30 mm (upper half) and 150 mm (lower half). The L / D is fixed at 25.
4.2. Large scale motors To complete this study we finally explore effectiveness of swirl in large port diameter grain. Simulations were carried out for selected grain port diameters (D = 30, 45, 90, and 150 mm) at inlet swirl numbers, S w = 0 and 1 and average inlet oxidizer mass flux of 132 kg/m2 /s. Geometrical similarity was maintained with the long motor of L / D = 25 to include the effect of grain length. Fig. 10 presents local rb along the length of the grain for port diameters of 30 mm and 150 mm, for inlet S w = 0 and 1. One can note that rb is lower for simulation with large port diameter. This is because of less fuel vapor added per unit grain port area for the large diameter grain. Another consequence of large port diameter is that the core gets heated at a distance further downstream compared to the small port diameter grain as illustrated by temperature contours in Fig. 11. The enhancement of average regression for various X / D is plotted in Fig. 12. The enhancement over the no-swirl case is seen to be higher for fuel grains with large port diameter. This is due to lower values of average rb for the no-swirl case. It is also noted that higher enhancement in average rb is achievable for short motors ( L / D < 10). This is due to the decay of swirl in the combustion chamber. The variation of local swirl number along the length of the grain is plotted in the inset of Fig. 12. The swirl strength is seen to decay by about 90% at L / D of 10. Thus swirl can be effectively used to increase rb in short motors especially for large port diameters. 5. Conclusions The effect of swirling oxidizer flow at the inlet on the rb in hybrid rocket motors was numerically studied by varying inlet S w , inlet velocity profiles and fuel grain length. Effect of grain scal-
Fig. 12. Enhancement in the average regression rate (in percentage) along the length of motor for large port diameters and swirling flow. The inset figure shows axial variation of local swirl strength ( S w ).
ing was also investigated. The following conclusions may be drawn from this study. 1. The presence of swirl velocity component increases the velocity magnitude in the flow field and causes increased velocity gradient at the fuel surface. Heat and mass transport at the fuel surface increases which consequently increases the rb . 2. Increase in inlet S w increases average rb monotonically. 3. Variation in axial velocity profile did not affect average rb . However, a swirl velocity profile with peak near the axis of the motor results in maximum average rb because of enhanced transport processes and flame positioned closer to the fuel surface. The contribution of swirl velocity profile in increasing average rb is only a fraction of contribution of swirl itself, the maximum contribution being about 50% at about inlet S w = 1. 4. For long fuel grains, the local rb and C ∗ for base case attain a maximum value at a location where all the available oxidizer is consumed and the flame reaches the axis (occurs at L / D ∼ 15 for D = 30 mm). With swirling oxidizer flow, this maximum is achieved in shorter L / D (here 10) due to increased rb . 5. The local and average rb decrease with increase in grain port diameter (swirl or no-swirl) due to lower increase in local mass flux by addition of pyrolyzed fuel. Interestingly, the
C.P. Kumar, A. Kumar / Aerospace Science and Technology 29 (2013) 92–99
percentage enhancement in average rb from swirl is higher for large diameter grains. 6. The effect of swirl is high near the head end and reduces downstream. The swirl strength was seen to drop to about 5% at L / D = 10. References [1] D. Altman, A. Holzman, Overview and history of hybrid rocket propulsion, in: M.J. Chiaverini, K.K. Kuo (Eds.), Fundamentals of Hybrid Rocket Combustion and Propulsion, in: AIAA Progress in Astronautics & Aeronautics, vol. 218, AIAA, Reston, VA, 2007, pp. 1–36. [2] ANSYS-Fluent Software Release Version No. 6.3.26, 2006. [3] Y.-S. Chen, T.H. Chou, B.R. Gu, J.S. Wu, B. Wu, Y.Y. Lian, L. Yang, Multiphysics simulations of rocket engine combustion, Computers and Fluids 45 (1) (2011) 29–36. [4] G.C. Cheng, R.C. Farmer, H.S. Jones, J.S. McFarlane, Numerical simulation of internal ballistics of a hybrid rocket motor, AIAA Paper 1994-0554, 1994. [5] G. Gariani, F. Maggi, L. Galfetti, Numerical simulation of HTPB combustion in a 2D hybrid slab combustor, Acta Astronautica 69 (2011) 289–296. [6] S. Gordon, B.J. McBride, Computer program for calculation of complex chemical equilibrium compositions and applications, NASA Reference Publication 1311, 1996. [7] X. Jiang, G.A. Siamas, L.C. Wrobel, Analytical equilibrium swirling inflow conditions for computational fluid dynamics, Technical notes, AIAA Journal 46 (4) (2008) 1015–1019. [8] H.J. Kim, Y.M. Kim, Numerical modeling of combustion processes of hybrid rocket engine, in: 37th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, AIAA Paper 2001-4504, 2001. [9] K. Kitagawa, T. Mitsutani, T. Ro, S. Yuasa, Effects of swirling liquid oxygen flow on combustion of a hybrid rocket engine, in: 40th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, AIAA Paper 2004-3479, 2004. [10] W.H. Knuth, M.J. Chiaverini, Solid-fuel regression rate behavior of vortex hybrid rocket engines, Journal of Propulsion and Power 18 (3) (2002) 600–609. [11] M. Kumar, A.F. Ghoniem, Multiphysics simulations of entrained flow gasification. Part I: Validating the non-reacting flow solver and the particle turbulent dispersion model, Energy and Fuels 26 (2012) 451–463. [12] C. Kumar, A. Kumar, On the role of radiation and dimensionality in predicting flow opposed flame spread over thin fuels, Combustion Theory and Modeling 16 (3) (2012) 537–569.
99
[13] C.P. Kumar, A. Kumar, Effect of diaphragms on regression rate in hybrid rocket motors, Journal of Propulsion and Power (2013), in press, http://dx.doi.org/ 10.2514/1.B34671. [14] K.K. Kuo, M.J. Chiaverini, Challenges of hybrid rocket propulsion in the 21st century, in: M.J. Chiaverini, K.K. Kuo (Eds.), Fundamentals of Hybrid Rocket Combustion and Propulsion, in: AIAA Progress in Astronautics & Aeronautics, vol. 218, AIAA, Reston, VA, 2007, pp. 593–638. [15] K.K. Kuo, Y.-C. Lu, M.J. Chiaverini, D.K. Johnson, N. Serin, G.A. Risha, C.L. Merkle, S. Venkateshwaran, Fundamental phenomena on fuel decomposition and boundary layer combustion processes with applications to hybrid rocket motors, Semi-Annual Progress Report, NASA-CR-201843, 1996. [16] C. Lee, Y. Na, J.W. Lee, Y.H. Byun, Effect of induced swirl flow on regression rate of hybrid rocket fuel by helical grain configuration, Aerospace Science and Technology 11 (2007) 68–76. [17] T.-S. Lee, A. Potapkin, The performance of a hybrid rocket with swirling GOx injection, in: XI International Conference on the Methods of Aerophysical Research, 2002. [18] A. Lewin, J. Dennis, B. Conley, D. Suzuku, Experimental determination of performance parameters for a polybutadiene/oxygen hybrid rocket, in: 28th AIAA/ASME/SAE/ASEE Joint Propulsion Conference & Exhibit, AIAA Paper 923590, 1992. [19] G.A. Marxman, Combustion in turbulent boundary layer on a vaporizing surface, in: Tenth International Symposium on Combustion, 1965, pp. 1337–1349. [20] F.R. Menter, Influence of free stream values on k–ω turbulence model predictions, AIAA Journal 30 (6) (1992) 1657–1659. [21] F.R. Menter, Two equation eddy-viscosity turbulence models for engineering applications, AIAA Journal 32 (8) (1994) 1598–1605. [22] A.V. Potapkin, T.S. Lee, Experimental study of thrust performance of a hybrid rocket motor with various methods of oxidiser injection, Combustion Explosion and Shock Waves 40 (4) (2004) 386–392. [23] P.A. Ramakrishna, P.J. Paul, H.S. Mukunda, Sandwich propellant combustion: Modeling and experimental comparison, Proceedings of the Combustion Institute 29 (2) (2002) 2963–2973. [24] K.H. Shin, C. Lee, Y.H. Yu, The experiments for the enhancement of regression rate of hybrid rocket fuel, KSME International 19 (10) (2005) 1926–1936. [25] W. Steenbergen, J. Voskamp, The rate of decay of swirl in turbulent pipe flow, Flow Measurement and Instrumentation 9 (2) (1998) 67–78. [26] E.P. Volchkov, V.I. Terekhov, V.V. Terekhov, Flow structure and heat and mass transfer in boundary layers with injection of chemically reacting substances (Review), Combustion Explosion and Shock Waves 40 (1) (2004) 1– 16.