Multiphysics simulations of rocket engine combustion
Yen-Sen Chen
a,⇑, T.H. Chou
b, B.R. Gu
b, J.S. Wu
b, Bill Wu
a, Y.Y. Lian
a, Luke Yang
aa
National Space Organization, Hsinchu Science Park, Hsinchu, Taiwan
b
Department of Mechnical Engineering, National Chiao Tung University, Hsinchu, Taiwan
a r t i c l e
i n f o
Article history: Received 29 June 2010
Received in revised form 23 August 2010 Accepted 14 September 2010
Available online 19 September 2010 Keywords:
Hybrid rocket combustion N2O–HTPB system
Radiative heat transfer Real-fluid model Finite-rate chemistry
a b s t r a c t
Recently, the hybrid rocket propulsion has become attractive to the research community and has devel-oped the trend to become an alternative to the conventional liquid and solid rockets. The hybrid rocket is a combination of both the solid and liquid systems with half of the plumbing of the liquid rocket but retaining its operational flexibility and avoiding the explosive nature of the solid rocket. Among available hybrid systems, the N2O (Nitrous Oxide)–HTPB (Hydroxyl-Terminated PolyButadiene) hybrid propulsion
represents the simplest but sufficiently efficient design. Unfortunately, even until now, research in devel-oping hybrid N2O–HTPB propulsion system still strongly depends on trials-and-errors, which are
time-consuming and expensive. Thus, detailed understanding of the fundamental combustion processes that are involved in the N2O–HTPB propulsion system can greatly impact the research community in this field.
This may further facilitate the successful modeling of the combustion processes and help improving the design of N2O–HTPB propulsion system in the future. A comprehensive numerical model with real-fluid
properties and finite-rate chemistry was developed in this research to predict the combustion flowfield inside a N2O–HTPB hybrid rocket system. Good numerical predictions as compared to experimental data
are also presented.
Ó 2010 Elsevier Ltd. All rights reserved.
1. Introduction
In searching for payload mass fraction performance enhance-ment of human’s access to space technical capabilities, ramjet and scramjet (supersonic combustion ramjet) studies have been a long-term research effort in the aerospace community since the 1960s. Since the 1980s, computational modeling approaches have been gradually adopted in the aerospace community in the devel-opment of combustion devices and space launch systems. The advent of high performance computing infrastructures with im-proved computational techniques has turned complex numerical modeling methods, with substantial requirements in computa-tional resources, into practical daily design analysis tools. These advancements have allowed the researchers and designers of combustion devices to simulate the combustion processes and flow physics in detailed spatial and temporal resolutions that are very difficult and expensive to obtain using experimental means. Numerical models using computational fluid dynamics (CFD) methods have been applied to liquid and solid rocket combustion systems with successful supports to the technical programs
[1–6]. Recent successes in flight demonstrations of the HyShot of the University of Queensland and Hyper-X of NASA, plus the flight
test of Boeing X-51 in late 2009, have promoted greater research interests in computational modeling of advanced propulsion sys-tems. Numerical modeling efforts in simulating the supersonic combustion processes and the diagnoses of problem areas have contributed to the success of these programs.
In the recent development of space launch systems, hybrid rocket propulsion has drawn a lot of attention, especially in the civilian space tourism community, and has been demonstrated to become a viable alternative to the liquid and solid rockets, although it is still in its fledgling stage as compared to the other two types of rockets. The liquid rocket is an efficient system; how-ever, it requires very complex and expensive plumbing systems. The solid rocket premixes the fuel and oxidizer as a grainy solid, which is dense and compact. However, it is highly explosive and lack of thrust control or termination. The hybrid rocket is a combi-nation of both the solid and liquid systems with half of the plumb-ing of the liquid rockets, but retainplumb-ing their flexibility of operation and avoiding the explosive nature of the solid rockets. The major advantages of the hybrid rockets include[7]: (1) safety without explosive concern, (2) flexibility in throttling and termination, (3) grain robustness without the danger of explosion originating from cracks, (4) propellant versatility, (5) temperature insensitivity to the operating chamber pressure, and (6) low cost benefiting from the safety features. Nevertheless, there are several disadvantages which include[7]: (1) low regression rate (1 mm/s) due to the nature of the diffusion flame, (2) low bulk density (low volumetric 0045-7930/$ - see front matter Ó 2010 Elsevier Ltd. All rights reserved.
doi:10.1016/j.compfluid.2010.09.010
⇑ Corresponding author.
E-mail addresses: yenchen@nspo.org.tw (Y.-S. Chen), chongsin@faculty.nctu. edu.tw(J.S. Wu).
Contents lists available atScienceDirect
Computers & Fluids
fuel loading) due to many ports for increasing grain surface area, (3) low combustion efficiency, (4) O/F shift with burning time and location, and (5) slower transient response to throttling.
There are many types of hybrid combustion systems, in which fuel is classically a solid and the oxidizer is a liquid or gas. These include different solid fuels such as hydrocarbons (rubbers, plastics and even papers) and metals, and a wide variety of liquid oxidizers such as oxygen, hydrogen peroxide, nitric acid, nitrogen tetroxide, and nitrous oxide. Typical examples of combination of fuel and oxi-dizer with optimum O/F, Isp (184–326 s) and characteristic veloc-ity (1224–2118 m/s) are discussed in Ref.[7]. Among these, most experimental data were published for propellants using LOX and HTBP with or without aluminum powder. Although no detailed data were available in the literature for the hybrid propellants using nitrous oxide as the oxidizer and HTBP related materials as the fuel, recent development of most commercial companies tends to utilize this combination, e.g. the SpaceShip One (led by Burt Ru-tan) which won the X Prize in September 2004 by flying up to an altitude of 100 km and landing safely back to Earth. The maximum vacuum Isp of the N2O–HTPB propulsion system is only fair around
250 s; however, it represents the simplest one, mainly due to the self-pressurization feature of the nitrous oxide (60 bar at room temperature), which greatly simplifies the plumbing systems. Unfortunately, even until now, research in developing the hybrid N2O–HTPB propulsion system still strongly depends on
trials-and-errors, which are time-consuming and expensive. Very few
[8,9]have attempted to model complicated reactive flow phenom-ena of a realistic hybrid propulsion system, in which they employ an energy-balanced surface decomposition model with the assumption that decomposition of C4H6 (Iso-butadiene) starts at
820 K, and modeled the combustion process using some reduced reaction mechanisms with tuned rate constants.
The above modeling efforts have reached success to some ex-tent based on fitting of experiments and numerical simulations. However, the real-fluid effects were not considered in their mod-els. This can affect the overall flow structure in the combustion chamber, especially near the injectors, and affect the combustion processes and heat transfer characteristics, which is the key to the regression rates of the solid grain. Thus, a comprehensive numerical model is developed in this research to include the real-fluid property effects and finite-rate chemistry in the simula-tion of an N2O–HTPB hybrid rocket system.
2. Method of approach
The present numerical method solves a set of governing equa-tions describing the conservation of mass, momentum (Navier– Stokes equations), energy, species concentration and turbulence quantities. The governing equations are written as:
@
q
@tþ @ @xj ðq
ujÞ ¼ 0 ð1Þ @qai
@t þ @ @xjðq
ujajÞ ¼ @ @xjq
D þl
tr
a @ai
@xj þxi
ð2Þ @q
ui @t þ @ @xj ðq
ujuiÞ ¼ @p @xi þ@sij
@xj ð3Þ @q
H @t þ @ @xj ðq
ujHÞ ¼ @p @tþ Qrþ @ @xj K Cp þl
trH
r
H þ @ @xj ðl
þl
tÞ K Cp þl
trH
r
ðV2=2Þ þ @ @xj K Cp þl
trH
uk @uj @xk 2 3uj @uk @xk ð4Þ @q
k @t þ @ @xjðq
ujkÞ ¼ @ @xjl
þl
trk
@k @xj þq
ðP
e
Þ ð5Þ @qe
@t þ @ @xj ðq
ujeÞ ¼ @ @xjl
þl
tr
e @e
@xj þq
e
k C1P
C2eþ C3P
2 =e
ð6ÞThe convection terms of the governing equations are discretized with a second-order upwind scheme. Second-order central schemes are applied to the diffusion and source terms. For com-plete description of the thermal environment in the combustion chamber, a radiative heat transfer model with a finite-volume inte-gration method[1,4]is also employed in the present model. In the combustion chamber, the main participating species in the radia-tion model are carbon-dioxide and hot steam. For transient flow computations, an efficient second-order time-marching scheme, which has been validated for vortex shedding and transient start-up nozzle flows[6], is employed in the present study. An effi-cient method for comprehensive real-fluid equations of state and fluid properties is also tested for liquid propellant combustion flows. These numerical models are important for high fidelity sim-ulations of combustion physics. They are described in the following sections.
2.1. Radiative heat transfer model
Radiation effects are important for combustion simulations when radiative participating species appear in the combustion products and in the pyrolysis gas decomposed from the solid pro-pellant. Radiative transfer model using a finite-volume integration method is employed in this study[10,11]. Consider the radiative transfer equation (RTE) in a Cartesian coordinate system as shown inFig. 1a. The balance of energy passing in a specified directionX
through a small differential volume in an absorbing–emitting and scattering medium can be written as:
ð
X
r
ÞIkðr;X
Þ ¼ ðjk
þrk
ÞIkðr;X
Þ þjk
Ib;kðrÞ þr
k 4p
Z 0¼4p Ikðr;X
0ÞkðX
0 !X
Þd0 ð7Þwhere the subscript k represents the wave-number; Ikðr;XÞ is the
spectral radiative intensity, which is a function of position r and angular directionX; Ib;kðrÞ is the blackbody radiative intensity at
the temperature of the medium;
j
kandr
kare the spectral absorp-tion and scattering coefficients, respectively; andUkðX0!XÞ is thescattering phase function from the incomingX0direction to the out-going directionX; and d0 indicates integration over a 4
p
angulardirections. The term on the left-hand side represents the gradient of the intensity in the directionX. The three terms on the right hand side represent the changes in intensity due to absorption and out-scattering, emission and in-scattering, respectively. If the wall bounding the medium emits and reflects diffusely, then the radia-tive boundary condition for Eq.(7)is given by
Ikðrw;
X
þÞ ¼ek
Ib;kðrwÞ þ ð1ek
Þp
Z nX<0Ikðrw;X
ÞjnX
jd ð8ÞwhereXþandXdenote the leaving and arriving radiative intensity directions, respectively;
e
kis the spectral wall emissivity; nrepre-sents the unit normal vector on the wall; dindicates integration
over all angles of the arriving directions. Eq.(7)is a complex inte-gro-differential equation whose exact analytical solution is only possible for very simple and specific settings. This intrinsic diffi-culty has resulted in the development of several approximated models. In the present model, a finite volume method (FVM) is used
to solve the RTE and the numerical analysis procedures are briefly described here. In the following analysis, subscripts k are dropped for the sake of brevity and they are added whenever necessary.
In the FVM, the spatial and angular domains are divided into a finite number of control volumes and control angles, respectively. Then, Eqs.(7) and (8)are integrated over each control volume and control angle. Since the FVM shares the same computational grid as the employed CFD approach, the considered spatial domain is divided into MA discrete control volumes and surfaces by a grid generator. For numerical analysis of the FVM, a representative con-trol volume resulted from the spatial domain division is shown in
Fig. 1b. By referring to the division practice for the spatial domain, the angular domain (seeFig. 1c) at a node centered in a control vol-ume is divided into Nh N/= MB control angles with Nh and N/
representing numbers of control angle in polar angle h and azi-muthal angle / directions, respectively. These MB discrete solid angles are non-overlapping and their sum is 4
p
. Unlike the selec-tion of a quadrature scheme in the discrete ordinates method (DOM), there is no specific restriction in selecting control angles in the FVM. However, the control angles are usually chosen in a manner that best capture the physics of a given problem. This is analogous to the selection of control volumes.Multiplying Eq. (7) by a representative control volume DV (Fig. 1b) and a control angleDXm(Fig. 1c), carrying out the integra-tion, and transforming the left-hand side of the equation from the volume integral to a surface integral based on the divergence the-orem, Eq.(7)then becomes
X i¼e;n;t Im i Ai Z DXmð
X
m niÞdX
m X i¼w;s;b Im i Ai Z DXmðX
m niÞdX
m ¼ ðj
þr
ÞImþj
Ibþr
4p
XMB m0¼1DX
m0U
mm0 Im0 " # PD
VDX
m ð9Þwhere the subscripts e, w, etc. indicate the values on the eastern, western, etc. surfaces of the volume; The subscript P represents the value at the central node of the control volume and A represents the control volume surface area; Umm0
is the averaged scattering phase function from the control angleDXm0 to the control angle
DXm. Dividing Eq.(9)byDXm, we have
X i¼e;n;t Im i AiDmi X i¼w;s;b Im i AiDmi ¼ ½ð
j
þr
ÞI m þj
Ibþr
4p
X M m0¼1DX
m0U
mm0Im0 PD
V ð10Þ where Dm i ¼ niX
m¼ 1DX
m Z DXmðX
m niÞdX
m ð11ÞIn Eq.(11), Dimis the product of a surface unit normal vector
and the averaged intensity directionXm. To close the above
equa-tion, relations are needed between the intensities on the control volume surfaces and the nodal intensities. One appropriate closure relation for complicated geometries is based on the step scheme, Fig. 1. (a) Coordinate system for radiative transfer equation, (b) a representative control volume, and (c) a representative control angle.
which sets the downstream surface intensities equal to the up-stream nodal intensities. Use of the step scheme avoids the nega-tive intensities, overshoots, and undershoots which may occur in other radiation schemes such as the diamond scheme, and positive scheme. Furthermore it has much less connection with the neigh-boring nodes, and thus it is particularly suitable for parallel com-putation. The effect of the step scheme on communication costs is not significant. The final discretized equation for the FVM can be written as am PI m P ¼ a m EI m E þ a m WI m Wþ a m NI m Nþ a m SI m S þ a m TI m T þ a m BI m Bþ b m ð12Þ
where the intensities with the subscripts E, W, etc. denote the east-ern, westeast-ern, etc. nodal intensities, and
am P ¼ X i¼e;n;t maxðAiDmi ;0Þ X i¼w;s;b minðAiDmi ;0Þ þ
D
Vðj
þr
Þ ð13aÞ am I ¼ minðAiDmi ;0Þ i ¼ e; n; t and I ¼ E; N; T ð13bÞ am I ¼ maxðAiDmi ;0Þ i ¼ w; s; b and I ¼ W; S; B ð13cÞ bm¼D
Vj
Ibþr
4p
XMB m0¼1DX
m0U
mm0 Im0 " # P ð13dÞThe preceding discretization is carried out along only one con-trol angle at a node. The same procedure should be applied to all of the MB control angles at all of the MA nodes. This forms MAMB systems of non-symmetric algebraic equations. A solution of these equations only represents radiative contribution at a single wave number. The radiative divergence is the quantity used in the en-ergy equation and it should consist of radiative contributions from all wave numbers. The radiative divergence is expressed in terms of the radiative intensities as
r
qr¼ Z 1 0jk
Z 4p 0 Ikdx
þjk
Z 4p 0 Ib;kdx
dk ¼ Z 1 0 4pjk
Ib;kjk
XMB m¼1 Im kDX
m " # dk ð14ÞA typical radiatively participating gas consists of many lines whose absorption coefficients vary rapidly with wave number. Thus, it becomes a very difficult and time-consuming to evaluate the radiative properties over the actual band contour and include them into the RTE. To avoid this difficulty, the spectrum can be di-vided into MC bands and the radiative properties are assumed con-stant over each band. The integrated quantity in Eq.(14)is found as the summation over all bands of the individual contribution for each band, that is,
r
qr¼ X j¼MC j¼1 4p
kjjIbjjj
XMB m¼1 ImkDX
m " #D
kj ð15ÞThe above approach to the spectral problem essentially corre-sponds to the spectral discretization and it represents a good com-promise between accuracy and computational time. The number of MC can be changed from one (gray gas model) to several hundred (narrow band model). Obviously, the use of higher MC number provides more accurate results.
2.2. Transient time-marching scheme
In the present numerical model, the temporal terms of the transport equations are modeled with a second-order backward difference scheme. There are two ways of performing time-march-ing procedures. First, the momentum, energy and pressure-correc-tion equapressure-correc-tions are grouped together within an iterapressure-correc-tion loop and
drive to convergence for each time step. Although the transient sub-iteration method described in the previous section is accurate for transient flow applications, it requires many sub-iterations (more for highly dynamic cases) to advance one time step, espe-cially for low-speed flows. In order to save computational effort in transient flow applications, this sub-iteration algorithm can be replaced with a more efficient method, such as the operator split-ting technique[12–14]. This method consists of a predictor step plus two corrector steps to drive the discretization errors to sec-ond-order accuracy.
Testing of the transient model is described in the following sec-tions. Some test cases are identified for time accurate modeling and data comparisons. They are: (1) the laminar and turbulent vor-tex shedding flows behind a circular cylinder; and (2) unsteady multi-species flows past a circular cylinder with and without chemical reactions. The mesh and test conditions for these cases are prepared with initial uniform flow conditions. Proper time step sizes were selected for good solutions of these transient cases. Knowing the cycle time period from the experimental data, the computational time step size is then determined such that there are at least 200 time steps within one cycle. For cases without experimental data a priori, grid and time-step convergence tests are required numerically. This is a good guideline based on previ-ous experience although some adjustment may be needed for spe-cial cases with multiple dominant frequencies and/or frequency modulations.
2.3. Real-fluid combustion modeling
In order to perform analyses of phase change phenomena in cavitation, liquid sprays or cryogenic fluid flows, real-fluid thermal and caloric equations of state (EOS) were developed for use with the present CFD code. The HBMS equations of state[15–17]were chosen for this purpose. Thermal and caloric equations of state, va-por pressure, heat of vava-porization, surface tension, and transva-port properties are modeled with the equations of state proposed by Hirshfelder et al., as employed by Dooley [18] and Anon. [19]
(we term these the HBMS equations of state), and with conven-tional correlations[20]for the other properties. The property cor-relations used were not chosen for their absolute accuracy, but for their validity over a wide range of temperatures and pressures and for requiring a minimum of data to describe a particular species. These correlations are explicit in density and temperature. These equations are:
HBMS thermal equation of state:
P Pc¼ X4 j¼1 Tj2 r X6 i¼1 Bijqi2 r ; Tr¼ T Tc ;
q
r¼q
q
cHBMS caloric equation of state:
H H0 RT ¼ Zc Z qr 0 P Tr @P @Tr qr " #
q
2 r dqr
þ Zc Pq
rTr 1These equations are not only of sufficiently high order that properties are accurately predicted for a wide range of conditions, but component sub-models may be easily modified. In this in-stance, the vapor pressure curve and the liquid phase density cor-relations were improved over the original HBMS formulations. Other equations of state were considered, but were found not to be as satisfactory as the HBMS equations.
Multi-component mixtures were treated by adding partial spe-cific volumes or pressures. The partial volume methodology is essential to improve the accuracy of the prediction when a small amount of multi-component vapor and a large amount of liquid are present at the same point. The fluid property routines also
include correlations for the transport properties. This was conve-nient since much of the methodology was similar.
2.4. Real-fluid model computational speed
In assessing the computational overhead associated with using the real-fluid property sub-model, numerical experiments of simu-lating some converging–diverging pipe flows with ideal gas and li-quid oxygen working fluids were conducted to quantify the CPU time overhead when the real-fluid model is activated. The results show an overall CPU time increase of about 85% for the real-fluid cases. These results indicate that there are rooms for improvement in order to enhance the turnaround time of running multi-phase flow problems. To improve the computational efficiency, a real-fluid table lookup procedure was developed for the present CFD model. This method is much more efficient than the original real-fluid model which involves the inversion of curve-fit data describing the equation of state. In average, the real-fluid table lookup model is only 15% slower than its ideal-gas counterpart. This means about 70% saving in CPU time with the real-fluid table lookup method.
2.5. Multi-species mixture implementation
For general applications of the present table-lookup real-fluid model, one must consider the conditions of multi-species mixing processes. To streamline the modeling approach for better perfor-mance in terms of computational efficiency and program coding clarity, functions for calculating enthalpy, specific heat and density are created based on the current code structure of the present CFD model. These new functions now consist of extra sections for the real-fluid properties based on the local pressure, temperature and species concentrations. The table lookup section is activated when the temperature is lower than the value that divides the real-fluid and perfect gas properties. Within the table lookup sec-tion, the pressure level is first used to define a property curve based on the database provided. Then, the fluid property is ob-tained through interpolation based on this new curve.
3. Real-fluid model combustion test
A benchmark LOX/GH2 single-injector combustor test case of DLR, Lampoldshausen, Germany (RCM2), is simulated to validate the current real-fluid table lookup model. This is one of the bench-mark cases presented in the 2nd international workshop on rocket combustion modeling[21]. The test conditions and geometry of this case are summarized as:
chamber pressure = 1 MPa;
LOX mass flow rate = 0.05 kg/s (at 85 K); GH2mass flow rate = 0.0237 kg/s (at 287 K);
chamber geometry = 0.05 m 0.05 m square; chamber length (injector to throat) = 0.472 m; LOX inner tube diameter = 0.005 m;
GH2annulus duct width = 0.0032 m;
injector lip thickness = 0.0003 m; nozzle throat diameter = 0.015 m.
Computational model for this case includes an axisymmetric mesh with a grid size of 18,715, a time step size of 1.0E06 s, real-fluid table for LOX, a 6-species and 9-reaction chemistry sys-tem and a 1500 K spark region near the injector lip region. 35,000 time steps were required for a converged solution. Most of the CPU time was spent for chamber pressure buildup to balance the overall mass conservation of the system (i.e. chamber pressure
is reduced initially and then gradually increased after the flame zone is established). Methods for speeding up the pressure buildup process for this type of problem are important for better overall efficiency in computational efforts. Multi-grid matrix solver has been shown to be effective in this aspect through the reduction of long wave errors.
Another conclusion obtained in this case study is that the cur-rent real-fluid table model only cost about 15% more CPU time per step than the original ideal-gas model. The extra CPU time is spent for table lookup of the species properties and a searching procedure to obtain the mixture temperature. This also demon-strates a substantial improvement in computational efficiency over the original real-fluid model.Fig. 2shows the flowfield solutions of the RCM2 test case. The predicted shape of the flame zone is in close resemblance with the flame Abel-transformed emission im-age of the experiment, which is the same result as the original real-fluid model[21].
4. Hybrid rocket combustion experimental setup
A typical hybrid rocket hot-fire testing rig is schematically shown in Fig. 3. The facility includes a horizontal motor casing made of stainless steel, a pressurized vertical tank containing ni-trous oxide, and a solenoid or pneumatic throttling valve between the motor and the tank. This throttling valve is used to start, par-tially open, and stop of nitrous oxide flow. Preliminary design of the horizontal motor casing along with several locations of instru-mentation is incorporated. There is a pintle-type injector at the end of plumbing from the tank, which is used to inject the liquid nitrous oxide into the pre-combustion chamber. Note the liquid nitrous oxide will become mainly gaseous nitrous oxide as it enters the pre-combustion chamber because of the lower pressure as compared to the saturation pressure of the nitrous oxide.
Both the pressures in the pre- and post-combustion chambers are measured using pressure transducers. Measurement of the pressure in the post-combustion chamber requires cooling of the combusted gases through long metal tube to protect the pressure transducer. Temperatures in the pre- and post-combustion cham-bers are measured using S-type thermocouples. In addition, only temperatures inside the HTPB grain are measured since the temperature in the central port is too high to measure using ther-mocouples. All data are acquired and processed using a computer-ized data acquisition system (National Instruments, Inc.) through a well-known LabView software. This is an effective, reliable and low-cost experimental setup.
5. Hybrid rocket combustion testing and numerical modeling The present hybrid rocket combustion model include a 580 mm combustion chamber with the design of forward-end and aft-end mixing chambers. A simple convergent–divergent conical nozzle is attached to the end of the combustion chamber. A single-port simple solid grain of HTPB is cast in two sections and assembled into a single segment through bonding. To boost the mixing effi-ciency, a mixing enhancer (patent pending) is also installed near the forward corner of the solid grain. A pintle-type injector made of stainless steal is employed for steady injection of the N2O
oxi-dizer. The nitrous oxide tank upstream of the injector and control valve is setup in a vertical position. Due to the properties of the oxidizer and the absence of thermal control for the oxidizer tank and without using a pressurant tank upstream, the injection flow rates of this pressure-fed injection system depend directly on the temperature of the environment. A small pyro grain is attached to the forward face of the solid grain, which is ignited at engine
start-up to melt the HTPB solid and start the combustion process after the N2O control valve is commanded to open.
Typically, a hot-fire test of this experiment runs between 15 and 20 s. During the test, motor thrust, chamber pressures and HTPB port temperatures are measured. The oxidizer flow rates and the regression rates of the solid grain are obtained through post-pro-cessing after the hot-fire tests. Overall, the experimental data show that the averaged solid regression rate is around 1 mm/s for the current design. The measured specific impulse, sea-level Isp, of the motor with the mixing enhancer is around 213 s (or vacuum Isp of 222.18 s). From previous experimental study without the mixing enhancer, the measured sea-level Isp was around 178 s for a similar motor.
To model this test case, an axisymmetric mesh system with 91,960 elements is generated to model the geometry of the
com-bustion chamber, pintle injector, solid grain and the conical nozzle. The mixing enhancer is not modeled in the present simulation due to its 3-D complexity, which will be explored in future research. A real-fluid property database is created to represent the nitrous oxide fluid in the numerical model. A 13-species, 17-reaction fi-nite-rate chemistry model is employed to model the combustion process. This includes the thermal decomposition of the oxidizer, pyrolysis of the HTPB solid and combustion of the gaseous mixture. A radiation model with finite-volume integration method is used to calculate the radiative heat transfer effects. The turbulence fields are predicted with an extended two-equation model[12,13]. The total pressure (50 ATM) and total temperature (283 K) boundary conditions are imposed at the injector inlet with the consideration of estimated total pressure losses between the tank and the injec-tor inlet. However, this setup still can not represent the true test
O
)
d
(
)
K
(
e
r
u
t
a
r
e
p
m
e
T
)
a
(
2Mass Fraction
H
)
e
(
n
o
i
t
c
a
r
F
s
s
a
M
H
O
)
b
(
2Mass Fraction
(c)
H
2O Mass Fraction
(f) Flame Abel-Transformed Emission Image of Experiment
Fig. 2. Predicted flowfield of the RCM2 test case with the current real-fluid-table combustion model.
conditions because the real inlet conditions are transient in nature. During the test, the tank pressure and temperature are reduced rapidly due to high flow rate and the expansion effects. Therefore, in the present model, we compare the numerical predictions to the test data when the oxidizer injection stays in the liquid phase, i.e. between 2 and 7 s after ignition.
The numerical simulation takes 80,000 time steps with 1
l
s time step size to obtain a quasi steady-state solution of theflow-field.Fig. 4shows the predicted hybrid rocket motor combustion flowfield using the present numerical model. It is clear that well-organized shear-layer oscillation structure, similar to the Taylor– Goertler type instability effects of thin shear layers, is predicted for the injection system. This oscillation is then coupled with the oxidizer decomposition process and the flame along the solid grain surface. Clearly, the axisymmetric simplification of the present model can not represent the real physics of shear layer instabilities.
O
)
e
(
)
K
(
e
r
u
t
a
r
e
p
m
e
T
)
a
(
2Mass Fraction
H
)
f
(
r
e
b
m
u
N
h
c
a
M
)
b
(
2O Mass Fraction
O
C
)
g
(
)
M
T
A
(
e
r
u
s
s
e
r
P
)
c
(
2Mass Fraction
(d)
Density
(kg/m
3)
(
h
)
C
4H
6Mass Fraction
And, it is plausible that the oxidizer evaporation may also contrib-ute to the oscillation of the predicted flowfield. From the density plot of Fig. 4, a long liquid (or high density fluid) core is clearly shown in the predicted flowfield. This indicates that further improvement in the injector design or mixing enhancement method in the chamber is required in order to increase the overall performance of this combustion system. The predicted averaged sea-level Isp is 181.2 s (or vacuum Isp of 191.18 s), which is lower than the measured data as expected due to the fact that a mixing enhancer is not present in the numerical model.
To increase computational efficiency of the present complex model with multiphysics capabilities, parallel computation is em-ployed on an IBM cluster computer system. Classes of 8-CPU, 32-CPU and 64-32-CPU were tested in the study. Results show parallel efficiency of around 0.97 (7.76 times speed up), 0.86 (27.52 times speed up) and 0.71 (45.44 times speed up), respectively, for these three classes of runs. The sharp decrease in parallel efficiency for the 64-CPU cases is due to the reduced mesh size for each proces-sor and increased communication between procesproces-sors.
6. Conclusions
In this study, we have developed a comprehensive numerical model for predicting the combustion flowfield of a hybrid rocket motor. A real-fluid property model has been employed for accu-rate description of the thermal-fluid environment of the N2O–
HTPB hybrid rocket combustion system. In the real-fluid property modeling, a table lookup method has been implemented in the present study to improve the computational efficiency for realistic flow physics representation. The improved model has shown good solutions as it is applied to a liquid rocket engine combustion test case conducted at DLR, Germany. The overall improvement in computational efficiency for real-fluid combustion simulations is about 70%. This provides substantial saving in computational ef-forts for the simulations of combustion processes with real-fluid properties.
The experimental investigation of the N2O–HTPB hybrid rocket
combustion system has demonstrated the effective, reliable and low-cost setup of the present design. The measured motor performance data is provided to anchor the numerical model developed in the present study. The present numerical solutions have shown interesting flowfield in the hybrid motor combus-tion chamber. The predicted quasi steady-state specific impulse shows close comparison to the test data. Nevertheless, the simulation is still a 2-D axisymmetric representation of the experiment, which can not simulate the true mixing process of the shear layer of the oxidizer injection. For future study, full 3-dimensional simulations in transient process will be modeled to include the effects of realistic turbulence mixing downstream of the injector.
Acknowledgements
This work is supported by the Sounding Rocket Program of the National Space Organization of the National Applied Research Lab-oratory. The computational resources and services are provided by the National Center for High-performance Computing of the Na-tional Applied Research Laboratory.
References
[1] Chen YS, Liu J, Zhang S, Mallapragada P. An integrated tool for launch vehicle base-heating analysis. Final report, NAS8-00002. Huntsville: Engineering Sciences, Inc., AL; 2001.
[2] Liaw P, Shang HM, Chen YS. Particulate multi-phase flowfield calculation with combustion/breakup models for solid rocket motor. AIAA-94-2780, 30th AIAA/ ASME/SAE/ASEE joint propulsion conference, Indianapolis, IN, June 27–29, 1994.
[3] Chen YS, Zhang S, Liu J. Stage separation performance analysis project. Final report, H-34345D. Huntsville, AL: Engineering Sciences, Inc.; 2002. [4] Wang TS, Chen YS, Liu J, Myrabo LN, Mead Jr FB. Advanced performance
modeling of experimental laser lightcraft. J Propul Power 2002;18(6):1129–38. [5] Wang TS. Multidimensional unstructured-grid liquid rocket engine nozzle
performance and heat transfer analysis. J Propul Power 2005;22(1):78–84. [6] Wang T-S. Transient 3-D analysis of nozzle side load in regeneratively cooled
engines. AIAA paper 2005-3942, 41st AIAA/ASME/SAE/ASEE joint propulsion conference, Tucson, Arizona; 2005.
[7] Fundamentals of hybrid rocket combustion and propulsion. In: Chiaverini MI, Kuo KK, editors. Progress in astronautics and aeronautics, vol. 218; 2007. [8] Farmer RC, Cheng GC, Jones SH, Arves JP. A practical CFD model for simulating
hybrid motors. JANNAF 35th combustion, airbreathing propulsion, and propulsion systems hazards subcommittee joint meeting, Tucson, AZ, 1998. [9] Farmer RC, Cheng GC, Metshc TH. Analysis of steady and unsteady hybrid
combustion. JANNAF propulsion & joint subcommittee meeting, Albuquerque, NM, 1996.
[10] Liu JW, Chen YS. Parallel simulation of radiative heat transfer using an unstructured finite-volume method. Numer Heat Transf B 1999;36:115–37. [11] Liu JW, Chen YS. Development of an unstructured radiation model applicable
for two-dimensional planar, axisymmetric and three-dimensional geometries. J Quant Spectrosc Radiative Transf 2000;22:17–33.
[12] Kim YM, Chen CP, Ziebarth JP, Chen YS. Prediction of high frequency combustion instability in liquid propellant rocket engines. AIAA 92-3763, AIAA/SAE/ASME/ASEE 28th joint propulsion conference and exhibit, Nashville, TN, July 6–8, 1992.
[13] Chen YS, Cheng GC, Farmer RC. Reacting and non-reacting flow simulation for film cooling in 2-D supersonic flows. AIAA 92-3602, AIAA/SAE/ASME/ASEE 28th joint propulsion conference and exhibit, Nashville, TN, July 6–8, 1992. [14] Wang TS, Chen YS. Unified Navier–Stokes flowfield and performance analysis
of liquid rocket engines. J Propul Power 1993;9(5):678–85. September– October.
[15] Uyehara OA, Watson KM. A universal viscosity correlation. Natl Petrol News 1944:R-714–22.
[16] Maxwell JB. Data book on hydrocarbons. NY: Van Nostrand Co.; 1950. [17] Farmer RC, Cheng GC, Chen YS. Parallel processing of the 3D3P codes.
SECA-FR-02-02, SECA, Inc., Huntsville, AL, 2002.
[18] Dooley RB. The international association for the properties of water and steam. September 1997 release. See also 2001 and 2003 releases available on IAPWS web site.
[19] Anon. ASME steam properties for industrial use. NY: ASME Press; 1998. [20] Dean LE, Shurley LA. Characteristics of RP-1 rocket fuel. Aerojet general
corporation, technical report TCR-70, 1957.
[21] Cheng G, Farmer RC, Chen YS. CFD simulation of liquid rocket engine injectors, part 2: simulations of the RCM-2 experiment. 2nd international workshop, rocket combustion modeling, Lampoldshausen, Germany, March 25–27, 2001.
