- Email: [email protected]
Analysis of hybrid rocket combustion
Acta Astronautica Vol.7, pp. 619--631 PergamonPressLtd., 1980. Printedin GreatBritain
Analysis of hybrid rocket combustion D. L. C H E R N G
AND C. C. T A O
Chung Shan Institute of Science and Technology, Taiwan, Republic of China (Received 10 December 1979)
Abstract--The problem involving flow of a turbulent boundary layer with mass injection and combustion is studied to investigate the hydrid rocket combustion. Because of the presence of combustion, the momentum equation is coupled to the energy and species equations. Particular attention is given in this paper to determine the eddy viscosity, and a combined model for it is proposed as follows: (1) the Cebeci's model for eddy viscosity is modified to allow for mass injection and combustion in the inner region; (2) for the eddy viscosity in the outer region of boundary layer, the Wilcox-Traci's two-equation model is effectively applied. The turbulent boundary layer equations are solved numerically in a finite-difference method developed by Patankar and Spalding. The input data are taken from Jones. According to the proposed model, the distribution of the viscosity parameter is obtained and similar in shape to Jones' data. Therefore, the proposed combined model of eddy viscosity has been developed to the stage where it can be used to predict quantitatively the behavior of compressible turbulent boundary layer with mass injection and combustion. However, the discrepancy between prediction and experimental data still exists, which may be due to the flame sheet assumption. Further investigation is needed to improve the prediction.
Introduction RECENT interest has been generated in turbulent b o u n d a r y layer with combustion p h e n o m e n a for the h y b r i d - r o c k e t - m o t o r investigation. A n u m b e r of studies were m a d e . t Considering the hybrid rocket c o m b u s t i o n in this paper, the problem proposing flow of a turbulent b o u n d a r y layer with mass injection and c o m bustion is studied. Because of the p r e s e n c e of the combustion, the density is related to the t e m p e r a t u r e and composition of the mixture in the b o u n d a r y layer. Thus, the m o m e n t u m equation is coupled to the energy and species equations. In this study, the following a s s u m p t i o n s are made: (1) the Lewis n u m b e r is unity, (Brokaw, 1962, Wooldridge and Muzzy, 1965). (2) The c o m b u s t i o n near the wall occurs in a zone that occupies a small portion of the total b o u n d a r y layer thickness, so that the "flame s h e e t " model is applied (Chen, 1963; Wooldridge and Muzzy, 1965). With the a b o v e assumptions, the conservation equations of species and energy are simplified. ,As to the e d d y viscosity which is a dominant portion of effective viscosity a p p e a r e d in all governing equations, a c o m b i n e d model is p r o p o s e d in this paper. tSee Marxman and Gilbert, 1963; Marxman a aL, 1964; Wooldridge and Muzzy, 1965; Jones and Isaacson, 1970; Jones et al., 1971; Helman and Manheimer-Timnat, 1973; Helman et al., 1974; and Gany et al., 1976. 619
620
D.L. Cherng and C. C. Tao
The viscosity parameter which appeared in Jones' analysis (Jones and Isaacson, 1970; Jones et al., 1971) is markedly reduced by combustion; thus, the Cebeci's model (Cebeci, 1971, Cebeci and Smith, 1974) for eddy viscosity is modified to allow for mass injection and combustion in the inner region. For the eddy viscosity in the outer region of boundary layer, the Wilcox-Traci's (1976) two-equation model is effectively applied. All governing differential equations and the Wilcox-Traci's K - Wt equations are parabolic equations, which can be solved in a finite-difference method developed by Patankar and Spalding (1970). The computer program is also considered to include Cebeci's model (1973) for variable turbulent Prandtl number, Touloukian's viscosity expressions, and variable specific heats (McBride et al., 1963). The input data are taken from Jones and Isaacson (1970). According to the proposed combined model, the distribution of viscosity parameter is obtained and similar in shape to Jones' et ai. data (1971). The prediction results of velocity profile, boundary layer thickness, flame sheet position and flame temperature are in good agreement with Jones and Isaacson' experimental data (1970).
Problem formulation To simplify the analysis of hybrid rocket combustion, the problem proposing flow of a compressible turbulent boundary layer with mass injection and combustion is shown in Fig. 1. A noncombustion turbulent flow extends to x = x2 on a flat plate; for x > x2, a mixture of hydrogen (H2) and nitrogen (N2) is injected through the plate into the boundary layer, and combustion flame front is formed between the wall and the edge of the boundary layer after ignition.
Governing equations This problem can be described as two-dimensional compressible, turbulent flow with mass injection and combustion. Assuming a Lewis number of unity and "flame sheet" model, the governing equations for the case of no pressure gradient in boundary layer can be expressed as: Conservation of mass
a(pu)lax +
(1)
a ( p v ) / a y = o.
s Ue
Air T
r
j._.~- - / j /
Edge of turbulent boundary layer
j..tJ-
8. ~ /
~ _ / -
////z.ZZlZ/////z//14-~
X I
-
4
~ X2 -
~
I
Flame front
~
,
'
--
4
T
~ Vw
X3
i
'
~4 .....................
Fig. 1. Illustration of the problem.
[ _
....
T W
_
~
Analysis o[ hybrid rocket combustion
621
Conservation of momentum
pu( Ou/Ox) + ov( Ou/Oy) = O{iz~( Ou/Oy)}/ Oy
(2)
where ~.n =/.t + ~t. Conservation of energy
ou( OH/ Ox) + ov( OHI Oy) = O{(t~e~lPr~n)(OH/ Oy)}/ Oy + O{tZe~(1-- l/Pr, n)O(u2/2)/Oy}/Oy
(3)
where H = h + u2/2. Conservation of species Ysz
Ou( OYs2/ Ox) + Ov( OYN2/Oy) = O{(I.~efdPren)(OYs2/ Oy)}/ Oy
(4)
with Lewis number: Le=l;
Scen=Pr~.
Conservation of species Y
pu( OY/ ax) + pv( d Y[ ay) = a{(lzen/PrefO(a Y/ ay)}/ Oy
(5)
where Y = Yo2-acY.2. The above boundary conditions are: Momentum u(x, O) = 0
v(x,0) = {0, /)w,
x<-x2 X > X2
(6)
lira u(x, y) = u,. y--~
Energy
H(x, O) = Hw(x) lim H(x, y) = H,. y~
(7)
Species YN~ for x > x:
YN~(x, 0) = YN~.,, lira YN~(X,y ) = YN2.e" y~
(8)
622
D.L. Cherng and C. C. Tao
Species Y for x > x2 Y ( x , O) = - acYH2,~,,
(9)
lira Y ( x , y ) = Yo2,e and y = 6 I at Y = 0 . The initial condition (at x = x 0 are: = ;~(x;). Initial velocity profile u is calculated from Van Driest's equation (1956): , ( - y,/A,)]2} 112 u , = foy. 1 + {1 + 4k,.Zy,2[12 -d yexp
(10)
where U , = Id/Ur
y, = yudvw u, = ( r d p , , ) u2 rw = c~p,u, 2
cI = 0.0592/Relx/5 A , = 26 k,, = 0.435 and initial temperature profile is approximated by Schlichting (1968): (T-
(12)
T w ) / ( T , - Tw) = UlUe.
A combined model
Effective viscosity/xo~ in all governing equations is defined as tzen =/.L +/~t.
(13)
The laminar viscosity/x, which is a function of temperature, can be expressed as (14)
IX = ~'. Y~Ui
where/z~ is the laminar viscosity for species i and is given as lxi • IIY = bl + b z T + b 3 T 2 + b 4 T 3 b s T 4 + b 6 T 5
[cp]
(15)
where T is in degrees Kelvin. The exponent j and the coefficients b l - b6 are taken from Touloukian's viscosity expressions.
Analysis of hybrid rocket combustion
623
The turbulent viscosity/z, is /x, = pe.
(16)
where e, is the eddy viscosity. In Jones' similar solution (Jones and Isaacson, 1970; Jones et al., 1971), a "viscosity parameter" 4' which appeared in transformed momentum equation was defined as ,b = (pt~lp, t~,)(1 +
(17)
e.lv).
Jones indicated that the viscosity parameter was markedly reduced by combustion. Thus for the eddy viscosity which is a dominant portion of effective viscosity appeared in all governing equations, a combined model is proposed as shown in Fig. 2, and described as a combination of modified Cebeci's model and Wilcox-Traci's two-equation model. Modified Cebeci's model. The Cebeci's model (Cebeci, 1971; Cebeci and Smith, 1974) for eddy viscosity is modified to allow for mass injection and combustion in the inner region of boundary layer. In this model, the eddy viscosity e.,., and emi.. are given by e.,., =/3c(k.y)2{1 - e x p (-ylA)}21OulSyJ, 8
y <- 8/
2
(18) ,,9,
where Bc = (I -
cay)"
(20)
c a = A~8i[82
(21)
A = 26(v/u,)(p/pw)m(l/N)
(22)
N = {exp [11.8(tzdlx)(vdu~)]} 'n.
(23)
!
[ A PPlied b Y Wilcox " Traci's" Outer region
IfY
Z
Inner region Applied by
i/I//11/i/////i/////////////////////////fl////
Fig. 2. Illustration of combined model.
~m, ~,% X
624
D.L. Cherng and C. C. Tao
The exponent n and constants A s in eqns (20) and (21) are supposed to be universal constants. They are obtained here by finding the values that give the best fit to experimental data. It is noted that the proposed model is reduced to Cebeci's model for noncombustion turbulent flow, because, in the absence of combustion, 8: = 0 and c~ = 0. Wilcox-Traci's two-equation model. Because the above modified Cebeci's model does not take into account the processes of convection or diffusion of turbulence which play an important part away from the wall for turbulent flow in the outer region of boundary layer, Wilcox and Traci (1976) suggested that the local eddy viscosity can be calculated from
(24)
emo = K I Wt
where the turbulent kinetic energy K and dissipation rate Wt are governed by following partial differential equations:
r = Tyt, 0 : ~' , ~'°~-~-] a r ~ + ~*p [IWyl 0 . I t - ~*pW, K p u ~0 -r + pv 00y (convection)
(diffusion)
(25)
(generation) (dissipation)
aWt2+ OWt2 _ a / aw, zx_ au 2 [K a' m~] 2]~Wt3 (26) OuT pv ay - ~-~t cru ¢f-~-y ] -,- a P l~y l Wt - p { [3 + 2or ][~y] [t--~-t (convection)
(diffusion)
(generation)
(dissipation)
where a*, a, /3*, B, or*, cr are supposed to be universal constants. The initial conditions are:
K = (u,2/a*)cos2(~*rylS) { Wt = K/(kmu,y), WI = K/(0.09u,8),
(27)
for y18 <-O.09/km for y/8 > O.09/km,
(28)
The boundary conditions are: K=0,
at y = 0 (29)
K = O.O15u~z/a*, at y = 8 Wt = 2500vda*k 2, Wt = 0.015ud(0.09a*8),
at y = 0
(3O) at y = &
Chamber and Wilcox (1977) discovered that the Wilcox-Traci's two-equation model admitted straight forward integration through the viscous sublayer which
Analysis of hybrid rocket combustion
625
is the source of difficulty in other two-equation models as Jones-Lauder's (1972) and NG-Spalding's (1972) models. Thus, for the eddy viscosity in the outer region of boundary layer, the Wilcox-Traci's model can be effectively applied to connect easily with the modified Cebeci's model in the inner region. The line separating the inner and outer regions is determined from the relation era, e m o at y = y,.. =
Turbulent P r a n d t l n u m b e r
An expression for the turbulent Prandtl number Prt was given by Cebeci (1973) as: Prt = er,/eh = kin{1 - exp ( - y/A)}/kh{1 - exp ( - y/B)}
(3 I)
where, for Z --- 0.3 (Z = Re x 10-3, Ro = u,O/p, and 0 = momentum thickness), values of kin, kh, A and B are chosen as km= 0.4 + 0.19/(1 + 0.49Z 2) kh = 0.44 + 0.22/(1 + 0.42Z 2)
A + = 26+ 14/(1 + Z 2) (32) A = A+t,(%lp) - m
B + = 35 + 25/(1 + 0.55Z 2) B = B+v(%/p) - m
and (%/p)112 = (~.dp,,)U2(pdp)mN.
It is noted that high Reynolds number (Re >5000), values of A ÷ = 2 6 and ks = 0.4 are used. While the turbulent Prandtl number is evaluated from eqn (31), the effective Prandtl number, by definition, is Prel = (v + em)/(a + eh) = {1 + (e,dv)}/Pr -I + Pril(em/~)}
(33)
where e,, and v are kinematic eddy viscosity and kinematic viscosity, respectively.
Fluid properties The fluid properties that appear in the momentum, energy and species
626
D.L. Cherng and C. C. Tao
equations are density (p) and enthalpy (h). The enthalpy h is calculated as:
h = ~, Yihi
(34)
where hi's are calculated from the relations (McBride et al., 1963) as:
hi = RT(al + a2T/2 + 83T2/3 + 84T3/4 + 85T4/5 + a6]T)
(35)
where R is universal gas constant and al - a6 are the coefficients which depend on the particular species. The density p is obtained from the equation of state and from the assumption that the static pressure Pe remains constant within the boundary layer, thus
p = (Pd T)/(R/M)
(36)
where M is the average molecular weight, i.e. M = {Y.(Yi/Mi)}-1. All properties, density, viscosity and enthalpy are related to temperature T. which follows from h. The enthalpy h is related to the total enthalpy H and velocity u and these quantities are obtained from the energy and momentum equation.
Numerical
solution
Transformation of equations All governing differential equations (2)-(5) and the Wilcox-Traci's K - Wt equations (25) and (26) are parabolic equations, and they have the general form: pu( cgflOx) + pv( Ofl Oy) = O{(l~enl~enj )( Ofl Oy)}/ Oy + $I
(37)
where f is equal to u, H, YN2, Y, K and Wt2. Partankar and Spalding (1970) developed a finite-difference method, using a marching, tri-diagonal matrix technique, and cast eqn (37) into the form: aflax +
f th" + ~o(rh"- rh~)] . . . .
1:
~,-~-~
jo110co=
{#utz•
}
a (~e - tk.,)2o'eft,i Of/am / a ~ + S'f
(38)
where the cross-stream coordinate, oJ, is introduced by o) =
(~b - ~#w)/(~b, - ~bw),
0 -< o~ _ 1
(39)
with stream function ~b defined as
dO = pu dy d~,ldx = - m ., . . .
(40) (pv ),
(41)
Analysis of hybrid rocket combustion
627
m""w-- - (pv)w.
d~bw/dx = -
(42)
The particular expressions for source form S~ and tre~j in eqn (38) are given in Table 1. Table 1. Expressions for S~ and tren~,
f
S~
u
~e~j
0
H
O{~(l-
Ys2 Y
1
1/Pre~)O(u2/2)/Oto}/Oto 0 0
K Wt~
a*JIpu/(~b,- ~bw)}Ou/&olK/u-
B*
WtK/u
- 4~,)}aulaw[ W H u - (l/u) x{B+ 2o'[~a(KIW,2)"21&o ]2} W~3 a[{pul(~b,
Pre~ Pr,~ Prea lhr*
l/tr
Full information of the numerical method can be found in the program developed by Patankar and Spalding. The flow diagram for this study is presented in Fig. 3. Recommended
constants
in t h e c o m b i n e d
model
Equations (18)-(26) for the combined model contains nine constants n, A,, kin, a, fl, tr, a*, fl* and tr*, which by hypothesis are universal constants. Based on the Saffman's arguments (1970), the six constants which appeared in the WilcoxTraci's model could all be estimated quite closely by some requirements on the quantitative structure of the solutions and a universal observation on the XO = Xu + ,~X Initial Xu
y, u, H, YN2 Y,K, WL T. p /~, r,,, u . 8 at xv I I I L _ _ m
Xu = XD y, U, H, YN2, Y K, Wt2, r, o, /z &. 7w, u. at Xu are stored in computer memory for next step.
Fig. 3. Flow diagram for computer program.
] Print out
628
D.L. Cherng and C. C. Tao
structure of turbulence. The other three constants which a p p e a r in the modified Cebeci's model can be also estimated by the c o m p a r i s o n of the numerical results with J o n e s ' experimental data. By the a b o v e techniques, the nine constants appearing in the combined model equations are given in Table 2, where the values of a * , / 3 " , or, ~r* are the same as the W i l c o x - T r a c i ' s model. Table 2. The values of the constants in the combined model n
A
k,.
ct
13
6-
o'*
a*
/3*
1.56
ll.3
0.435
0.171
0.165
0.5
0.5
0.3
0.09
Results and comparisons
The input data are t a k e n f r o m Jones and I s a a c s o n (1970) and listed in Table 3. According to the p r o p o s e d c o m b i n e d model and r e c o m m e n d e d constants, the distribution of the viscosity p a r a m e t e r is obtained as s h o w n in Fig. 4 and it is similar in shape to Jones et al. data (1971). Table 3. Input data Input data xl x2
x3 8~1 u, P~ T~ Tw Yo2.~ YN:.,
Y.2.w YN2.w F = (pv)wl(pu)e
at
Value 63.5 cm 76.2 cm 118.1 cm 3.0734 cm 1024.13 cm/sec 1.0334 kg/cm 302.8°K 350.0°K 0.232 0.768 0.038 0.962 3.12 × 10-3 8.0
The predicted results of velocity profile at location x3 along the fiat plate is c o m p a r e d with the m e a s u r e d data b y Jones and I s a a c s o n (1970) as shown in Fig. 5. Also, the results of b o u n d a r y layer thickness & flame sheet position 3t and flame t e m p e r a t u r e TI are listed in Table 4. The c o m p a r i s o n s show that numerical results of this study are in good a g r e e m e n t with J o n e s ' experimental data.
Analysis of hybrid rocket combustion
629
I00 /
60
Prediction Experiment
\
~
80
\\
40
20
o
i
0
0.2
i
i
0.4
0.6
1.0
0.8
9 /r/e
Fig. 4. Comparison of computed and measured viscosity parameter. I 0
0.8
0.6
/
0.4
0.2
//~
Prediction Experiment -
_
_
0 0
1.0
2.0
3.0 y~
4.0
5.0
6.0
70
cm
Fig. 5. Comparison of computed and measured velocity profile.
Table 4. Comparison of computed and measured 6, 8t and Tf
Jones' experimental data Predicted results
8
8~
r~
6.731 cm 6.662 cm
1.524 cm 1.386 cm
1185°K 1184.9°K
Concluding remarks Based on this investigation for the hybrid rocket combustion, the proposed c o m b i n e d model of eddy viscosity has been developed to the stage where it can be used to predict quantitatively the behavior of compressible turbulent boundary layer with mass injection and combustion. H o w e v e r , the discrepancies b e t w e e n prediction and experimental data still exist in Figs. 4 and 5, which m a y be due to the flame sheet assumption. Further investigation is needed to improve the prediction.
630
D . L . Cherng and C. C. Tao
Acknowledgement--The co-author of this paper would like to express his sincere appreciation to Dr. H. A. Hassan, Professor, North Carolina State University, U.S.A., for guidance and encouragement during the course in part of this study. References Brokaw R. S. (1962) The Lewis number. Progress in International Research on Thermodynamic and Transport Properties, pp. 271-278. Academic Press, New York. Cebeci T. (1971) Calculation of compressible turbulent boundary layers with heat and mass transfer. AIAA J. 9, 1091-1097. Cebeci T. (1973) A model for eddy conductivity and turbulent Prandtl number. Z Heat Trans. (ASME) 95, 227-234. Cebeci T. and Smith A. M. O. (1974) Analysis of Turbulent Boundary Layers. Academic Press, New York. Chen T. N. and Toong T. Y. (1963) Laminar boundary layer wedge flows with evaporation and combustion. Prog. in Astron. Aeron. 15, 643--664. Chamber T. L. and Wilcox D. C. (1977) Critical examination of two-equation turbulence closure models for boundary layer. AIAA 3". 15, 821-828. Gany A., Manheimer-Timnat Y. and Wolfshtein M. (1976) Two-phase flow effect on hybrid combustion. Acta Astronautica 3, 241-263. Helman D. and Manheimer-Timnat Y. (1973) Investigation of hybrid rocket combustion. Israel Z Technol II, 17-25. Helman D., Wolfshtein M. and Manheimer-Timnat Y. (1974) Theoretical investigation of hybrid rocket combustion by numerical methods. Combustion and Flame 22, 171-190. Jones J. W. and Isaacson L. K. (1970) A Turbulent Boundary Layer with Mass Addition, Combustion, and Pressure Gradient. AFOSR Scientific Pep. 70-1428TR. Jones J. W., lsaacson L. K and Vreeke S. (1977) A turbulent boundary layer with mass addition, combustion and pressure gradient. AIAA J. 9, 1762-1768. Jones W. P. and Lauder B. E. (1972) The prediction of laminarization with a two-equation model of turbulence. Int. J. Heat and Mass Trans. 15, 301-314. Marxman G. and Gilbert M. (1963) Turbulent boundary layer combustion in the hybrid rocket. 9th Syrup. (Int.) on Combustion, pp. 371-383. Marxman G., Wooldridge C. E. and Muzzy R. J. (1964) Fundamental of hybrid boundary-layer combustion. Prog. in Astron. Aeron 15, 485-522. McBride B. J., Heimel S., Ehlers J. G. and Gordon S. (1963) Thermodynamic Properties to 6000 K for 210 Substances Involving the First 18 Elements. NASA-SP-3001. NG K. H. and Spalding D. B. (1972) Turbulence model for boundary layer near walls. Physics Fluids 15, 20-30. Patankar S. V. and Spalding D. B. (1970) Heat and Mass Transfer in Boundary Layer. lntertext Books, London. Schlichting H. (1968) Boundary Layer Theory. McGraw-Hill, New York. Saffman P. G. (1970) A model for inhomogeneous turbulent flow. Pro. R. Sac. A. 317, 417--433. Touloukian Y. S. Viscosity. TRPC Data Book, Vol. 2. Purdue University, Lafayette, Indiana. Van Driest E. R. (1956) On turbulent flow near a wall. J. Aeronaut. Sci. 23, 1007-1011. Wooldridge C. E. and Muzzy R. J. (1965) Measurements in a turbulent boundary layer with porous wall injection and combustion, lOth Syrup. (Int.) on Combustion, pp. 1351-1362. Wilcox D. C. and Traci R. M. (1976) A Complete Model of Turbulence. AIAA paper 76--351, San Diego, California.
Appendix
Nomenclature a thermal diffusivity, cm21sec A B constant in eqn (21)
c~ coefficient in eqn (20) cf local surface skin friction coefficient
Analysis of hybrid rocket combustion
[ F h H k k,. K Le rh" M
dependent variables mass injection rate, (pv)w/(pu), enthalpy, cal/mole total enthalpy, h + u2/2, cal/mole roughness height of the wall, cm Von Karman constant turbulent kinetic energy, cm2/secz Lewis number mass transfer rate across a boundary, kg/cm2 • sec average of molecular weight,
{X(YJM~)}-' n N P Pr Prf Pren R Re Scefr T
an exponent in eqn (20) defined in eqn (23) pressure, kg/cm z laminar Prandtl number turbulent Prandti number effective Prandtl number universal gas constant, call°K • kg-mole Reynolds number effective Schmidt number temperature, °K u,v x and y component of velocity, cm/sec u , dimensionless velocity, u/u, u, friction velocity, (rw/pw)tt2, cm/sec x,y rectangular coordinates, cm y,. value of y separating two model, cm y , dimensionless distance, yudv¢ Y species, Yo2-acYH2 Ya2 mass fraction of hydrogen Ys2 mass fraction of nitrogen Yo~ mass fraction of oxygen a,a* constants in Wilcox-Traci's two-equation model ac stoichiometric oxygen/hydrogen mass ratio
631
/3,/3* constants in Wilcox-Traci's two equation model /3c defined in eqn (20) 0 momentum thickness, cm stream function, kg/cm • sec tz laminar viscosity, kg/cm, sec /Lt turbulent viscosity, kg/cm, sec tze~ effective viscosity, tz +/~,, kg/cm • sec p density, kg/cm3 v kinematic viscosity, cm2/sec e,. eddy kinematic viscosity, e,, = Izt/p, cmZ/sec 8 boundary layer thickness cm 8r flame front location, cm or,or* constants in Wilcox-Traci's two-equation model cre~,r effective Prandtl number for variable / z shear stress, kg/cm 2 ~o dimensionless stream function viscosity parameter, defined in eqn (17) rt
"o = 0 , ~ , / ( 2 ~ )
,~ l~ =
'2
fo'
(pip,) dy
O~u~, dx, kg 2 seC/cm 4
Subscripts D e .f i o U w
downstream outer edge of boundary layer flame front inner region; species outer region upstream wall
Analysis of hybrid rocket combustion D. L. C H E R N G
AND C. C. T A O
Chung Shan Institute of Science and Technology, Taiwan, Republic of China (Received 10 December 1979)
Abstract--The problem involving flow of a turbulent boundary layer with mass injection and combustion is studied to investigate the hydrid rocket combustion. Because of the presence of combustion, the momentum equation is coupled to the energy and species equations. Particular attention is given in this paper to determine the eddy viscosity, and a combined model for it is proposed as follows: (1) the Cebeci's model for eddy viscosity is modified to allow for mass injection and combustion in the inner region; (2) for the eddy viscosity in the outer region of boundary layer, the Wilcox-Traci's two-equation model is effectively applied. The turbulent boundary layer equations are solved numerically in a finite-difference method developed by Patankar and Spalding. The input data are taken from Jones. According to the proposed model, the distribution of the viscosity parameter is obtained and similar in shape to Jones' data. Therefore, the proposed combined model of eddy viscosity has been developed to the stage where it can be used to predict quantitatively the behavior of compressible turbulent boundary layer with mass injection and combustion. However, the discrepancy between prediction and experimental data still exists, which may be due to the flame sheet assumption. Further investigation is needed to improve the prediction.
Introduction RECENT interest has been generated in turbulent b o u n d a r y layer with combustion p h e n o m e n a for the h y b r i d - r o c k e t - m o t o r investigation. A n u m b e r of studies were m a d e . t Considering the hybrid rocket c o m b u s t i o n in this paper, the problem proposing flow of a turbulent b o u n d a r y layer with mass injection and c o m bustion is studied. Because of the p r e s e n c e of the combustion, the density is related to the t e m p e r a t u r e and composition of the mixture in the b o u n d a r y layer. Thus, the m o m e n t u m equation is coupled to the energy and species equations. In this study, the following a s s u m p t i o n s are made: (1) the Lewis n u m b e r is unity, (Brokaw, 1962, Wooldridge and Muzzy, 1965). (2) The c o m b u s t i o n near the wall occurs in a zone that occupies a small portion of the total b o u n d a r y layer thickness, so that the "flame s h e e t " model is applied (Chen, 1963; Wooldridge and Muzzy, 1965). With the a b o v e assumptions, the conservation equations of species and energy are simplified. ,As to the e d d y viscosity which is a dominant portion of effective viscosity a p p e a r e d in all governing equations, a c o m b i n e d model is p r o p o s e d in this paper. tSee Marxman and Gilbert, 1963; Marxman a aL, 1964; Wooldridge and Muzzy, 1965; Jones and Isaacson, 1970; Jones et al., 1971; Helman and Manheimer-Timnat, 1973; Helman et al., 1974; and Gany et al., 1976. 619
620
D.L. Cherng and C. C. Tao
The viscosity parameter which appeared in Jones' analysis (Jones and Isaacson, 1970; Jones et al., 1971) is markedly reduced by combustion; thus, the Cebeci's model (Cebeci, 1971, Cebeci and Smith, 1974) for eddy viscosity is modified to allow for mass injection and combustion in the inner region. For the eddy viscosity in the outer region of boundary layer, the Wilcox-Traci's (1976) two-equation model is effectively applied. All governing differential equations and the Wilcox-Traci's K - Wt equations are parabolic equations, which can be solved in a finite-difference method developed by Patankar and Spalding (1970). The computer program is also considered to include Cebeci's model (1973) for variable turbulent Prandtl number, Touloukian's viscosity expressions, and variable specific heats (McBride et al., 1963). The input data are taken from Jones and Isaacson (1970). According to the proposed combined model, the distribution of viscosity parameter is obtained and similar in shape to Jones' et ai. data (1971). The prediction results of velocity profile, boundary layer thickness, flame sheet position and flame temperature are in good agreement with Jones and Isaacson' experimental data (1970).
Problem formulation To simplify the analysis of hybrid rocket combustion, the problem proposing flow of a compressible turbulent boundary layer with mass injection and combustion is shown in Fig. 1. A noncombustion turbulent flow extends to x = x2 on a flat plate; for x > x2, a mixture of hydrogen (H2) and nitrogen (N2) is injected through the plate into the boundary layer, and combustion flame front is formed between the wall and the edge of the boundary layer after ignition.
Governing equations This problem can be described as two-dimensional compressible, turbulent flow with mass injection and combustion. Assuming a Lewis number of unity and "flame sheet" model, the governing equations for the case of no pressure gradient in boundary layer can be expressed as: Conservation of mass
a(pu)lax +
(1)
a ( p v ) / a y = o.
s Ue
Air T
r
j._.~- - / j /
Edge of turbulent boundary layer
j..tJ-
8. ~ /
~ _ / -
////z.ZZlZ/////z//14-~
X I
-
4
~ X2 -
~
I
Flame front
~
,
'
--
4
T
~ Vw
X3
i
'
~4 .....................
Fig. 1. Illustration of the problem.
[ _
....
T W
_
~
Analysis o[ hybrid rocket combustion
621
Conservation of momentum
pu( Ou/Ox) + ov( Ou/Oy) = O{iz~( Ou/Oy)}/ Oy
(2)
where ~.n =/.t + ~t. Conservation of energy
ou( OH/ Ox) + ov( OHI Oy) = O{(t~e~lPr~n)(OH/ Oy)}/ Oy + O{tZe~(1-- l/Pr, n)O(u2/2)/Oy}/Oy
(3)
where H = h + u2/2. Conservation of species Ysz
Ou( OYs2/ Ox) + Ov( OYN2/Oy) = O{(I.~efdPren)(OYs2/ Oy)}/ Oy
(4)
with Lewis number: Le=l;
Scen=Pr~.
Conservation of species Y
pu( OY/ ax) + pv( d Y[ ay) = a{(lzen/PrefO(a Y/ ay)}/ Oy
(5)
where Y = Yo2-acY.2. The above boundary conditions are: Momentum u(x, O) = 0
v(x,0) = {0, /)w,
x<-x2 X > X2
(6)
lira u(x, y) = u,. y--~
Energy
H(x, O) = Hw(x) lim H(x, y) = H,. y~
(7)
Species YN~ for x > x:
YN~(x, 0) = YN~.,, lira YN~(X,y ) = YN2.e" y~
(8)
622
D.L. Cherng and C. C. Tao
Species Y for x > x2 Y ( x , O) = - acYH2,~,,
(9)
lira Y ( x , y ) = Yo2,e and y = 6 I at Y = 0 . The initial condition (at x = x 0 are: = ;~(x;). Initial velocity profile u is calculated from Van Driest's equation (1956): , ( - y,/A,)]2} 112 u , = foy. 1 + {1 + 4k,.Zy,2[12 -d yexp
(10)
where U , = Id/Ur
y, = yudvw u, = ( r d p , , ) u2 rw = c~p,u, 2
cI = 0.0592/Relx/5 A , = 26 k,, = 0.435 and initial temperature profile is approximated by Schlichting (1968): (T-
(12)
T w ) / ( T , - Tw) = UlUe.
A combined model
Effective viscosity/xo~ in all governing equations is defined as tzen =/.L +/~t.
(13)
The laminar viscosity/x, which is a function of temperature, can be expressed as (14)
IX = ~'. Y~Ui
where/z~ is the laminar viscosity for species i and is given as lxi • IIY = bl + b z T + b 3 T 2 + b 4 T 3 b s T 4 + b 6 T 5
[cp]
(15)
where T is in degrees Kelvin. The exponent j and the coefficients b l - b6 are taken from Touloukian's viscosity expressions.
Analysis of hybrid rocket combustion
623
The turbulent viscosity/z, is /x, = pe.
(16)
where e, is the eddy viscosity. In Jones' similar solution (Jones and Isaacson, 1970; Jones et al., 1971), a "viscosity parameter" 4' which appeared in transformed momentum equation was defined as ,b = (pt~lp, t~,)(1 +
(17)
e.lv).
Jones indicated that the viscosity parameter was markedly reduced by combustion. Thus for the eddy viscosity which is a dominant portion of effective viscosity appeared in all governing equations, a combined model is proposed as shown in Fig. 2, and described as a combination of modified Cebeci's model and Wilcox-Traci's two-equation model. Modified Cebeci's model. The Cebeci's model (Cebeci, 1971; Cebeci and Smith, 1974) for eddy viscosity is modified to allow for mass injection and combustion in the inner region of boundary layer. In this model, the eddy viscosity e.,., and emi.. are given by e.,., =/3c(k.y)2{1 - e x p (-ylA)}21OulSyJ, 8
y <- 8/
2
(18) ,,9,
where Bc = (I -
cay)"
(20)
c a = A~8i[82
(21)
A = 26(v/u,)(p/pw)m(l/N)
(22)
N = {exp [11.8(tzdlx)(vdu~)]} 'n.
(23)
!
[ A PPlied b Y Wilcox " Traci's" Outer region
IfY
Z
Inner region Applied by
i/I//11/i/////i/////////////////////////fl////
Fig. 2. Illustration of combined model.
~m, ~,% X
624
D.L. Cherng and C. C. Tao
The exponent n and constants A s in eqns (20) and (21) are supposed to be universal constants. They are obtained here by finding the values that give the best fit to experimental data. It is noted that the proposed model is reduced to Cebeci's model for noncombustion turbulent flow, because, in the absence of combustion, 8: = 0 and c~ = 0. Wilcox-Traci's two-equation model. Because the above modified Cebeci's model does not take into account the processes of convection or diffusion of turbulence which play an important part away from the wall for turbulent flow in the outer region of boundary layer, Wilcox and Traci (1976) suggested that the local eddy viscosity can be calculated from
(24)
emo = K I Wt
where the turbulent kinetic energy K and dissipation rate Wt are governed by following partial differential equations:
r = Tyt, 0 : ~' , ~'°~-~-] a r ~ + ~*p [IWyl 0 . I t - ~*pW, K p u ~0 -r + pv 00y (convection)
(diffusion)
(25)
(generation) (dissipation)
aWt2+ OWt2 _ a / aw, zx_ au 2 [K a' m~] 2]~Wt3 (26) OuT pv ay - ~-~t cru ¢f-~-y ] -,- a P l~y l Wt - p { [3 + 2or ][~y] [t--~-t (convection)
(diffusion)
(generation)
(dissipation)
where a*, a, /3*, B, or*, cr are supposed to be universal constants. The initial conditions are:
K = (u,2/a*)cos2(~*rylS) { Wt = K/(kmu,y), WI = K/(0.09u,8),
(27)
for y18 <-O.09/km for y/8 > O.09/km,
(28)
The boundary conditions are: K=0,
at y = 0 (29)
K = O.O15u~z/a*, at y = 8 Wt = 2500vda*k 2, Wt = 0.015ud(0.09a*8),
at y = 0
(3O) at y = &
Chamber and Wilcox (1977) discovered that the Wilcox-Traci's two-equation model admitted straight forward integration through the viscous sublayer which
Analysis of hybrid rocket combustion
625
is the source of difficulty in other two-equation models as Jones-Lauder's (1972) and NG-Spalding's (1972) models. Thus, for the eddy viscosity in the outer region of boundary layer, the Wilcox-Traci's model can be effectively applied to connect easily with the modified Cebeci's model in the inner region. The line separating the inner and outer regions is determined from the relation era, e m o at y = y,.. =
Turbulent P r a n d t l n u m b e r
An expression for the turbulent Prandtl number Prt was given by Cebeci (1973) as: Prt = er,/eh = kin{1 - exp ( - y/A)}/kh{1 - exp ( - y/B)}
(3 I)
where, for Z --- 0.3 (Z = Re x 10-3, Ro = u,O/p, and 0 = momentum thickness), values of kin, kh, A and B are chosen as km= 0.4 + 0.19/(1 + 0.49Z 2) kh = 0.44 + 0.22/(1 + 0.42Z 2)
A + = 26+ 14/(1 + Z 2) (32) A = A+t,(%lp) - m
B + = 35 + 25/(1 + 0.55Z 2) B = B+v(%/p) - m
and (%/p)112 = (~.dp,,)U2(pdp)mN.
It is noted that high Reynolds number (Re >5000), values of A ÷ = 2 6 and ks = 0.4 are used. While the turbulent Prandtl number is evaluated from eqn (31), the effective Prandtl number, by definition, is Prel = (v + em)/(a + eh) = {1 + (e,dv)}/Pr -I + Pril(em/~)}
(33)
where e,, and v are kinematic eddy viscosity and kinematic viscosity, respectively.
Fluid properties The fluid properties that appear in the momentum, energy and species
626
D.L. Cherng and C. C. Tao
equations are density (p) and enthalpy (h). The enthalpy h is calculated as:
h = ~, Yihi
(34)
where hi's are calculated from the relations (McBride et al., 1963) as:
hi = RT(al + a2T/2 + 83T2/3 + 84T3/4 + 85T4/5 + a6]T)
(35)
where R is universal gas constant and al - a6 are the coefficients which depend on the particular species. The density p is obtained from the equation of state and from the assumption that the static pressure Pe remains constant within the boundary layer, thus
p = (Pd T)/(R/M)
(36)
where M is the average molecular weight, i.e. M = {Y.(Yi/Mi)}-1. All properties, density, viscosity and enthalpy are related to temperature T. which follows from h. The enthalpy h is related to the total enthalpy H and velocity u and these quantities are obtained from the energy and momentum equation.
Numerical
solution
Transformation of equations All governing differential equations (2)-(5) and the Wilcox-Traci's K - Wt equations (25) and (26) are parabolic equations, and they have the general form: pu( cgflOx) + pv( Ofl Oy) = O{(l~enl~enj )( Ofl Oy)}/ Oy + $I
(37)
where f is equal to u, H, YN2, Y, K and Wt2. Partankar and Spalding (1970) developed a finite-difference method, using a marching, tri-diagonal matrix technique, and cast eqn (37) into the form: aflax +
f th" + ~o(rh"- rh~)] . . . .
1:
~,-~-~
jo110co=
{#utz•
}
a (~e - tk.,)2o'eft,i Of/am / a ~ + S'f
(38)
where the cross-stream coordinate, oJ, is introduced by o) =
(~b - ~#w)/(~b, - ~bw),
0 -< o~ _ 1
(39)
with stream function ~b defined as
dO = pu dy d~,ldx = - m ., . . .
(40) (pv ),
(41)
Analysis of hybrid rocket combustion
627
m""w-- - (pv)w.
d~bw/dx = -
(42)
The particular expressions for source form S~ and tre~j in eqn (38) are given in Table 1. Table 1. Expressions for S~ and tren~,
f
S~
u
~e~j
0
H
O{~(l-
Ys2 Y
1
1/Pre~)O(u2/2)/Oto}/Oto 0 0
K Wt~
a*JIpu/(~b,- ~bw)}Ou/&olK/u-
B*
WtK/u
- 4~,)}aulaw[ W H u - (l/u) x{B+ 2o'[~a(KIW,2)"21&o ]2} W~3 a[{pul(~b,
Pre~ Pr,~ Prea lhr*
l/tr
Full information of the numerical method can be found in the program developed by Patankar and Spalding. The flow diagram for this study is presented in Fig. 3. Recommended
constants
in t h e c o m b i n e d
model
Equations (18)-(26) for the combined model contains nine constants n, A,, kin, a, fl, tr, a*, fl* and tr*, which by hypothesis are universal constants. Based on the Saffman's arguments (1970), the six constants which appeared in the WilcoxTraci's model could all be estimated quite closely by some requirements on the quantitative structure of the solutions and a universal observation on the XO = Xu + ,~X Initial Xu
y, u, H, YN2 Y,K, WL T. p /~, r,,, u . 8 at xv I I I L _ _ m
Xu = XD y, U, H, YN2, Y K, Wt2, r, o, /z &. 7w, u. at Xu are stored in computer memory for next step.
Fig. 3. Flow diagram for computer program.
] Print out
628
D.L. Cherng and C. C. Tao
structure of turbulence. The other three constants which a p p e a r in the modified Cebeci's model can be also estimated by the c o m p a r i s o n of the numerical results with J o n e s ' experimental data. By the a b o v e techniques, the nine constants appearing in the combined model equations are given in Table 2, where the values of a * , / 3 " , or, ~r* are the same as the W i l c o x - T r a c i ' s model. Table 2. The values of the constants in the combined model n
A
k,.
ct
13
6-
o'*
a*
/3*
1.56
ll.3
0.435
0.171
0.165
0.5
0.5
0.3
0.09
Results and comparisons
The input data are t a k e n f r o m Jones and I s a a c s o n (1970) and listed in Table 3. According to the p r o p o s e d c o m b i n e d model and r e c o m m e n d e d constants, the distribution of the viscosity p a r a m e t e r is obtained as s h o w n in Fig. 4 and it is similar in shape to Jones et al. data (1971). Table 3. Input data Input data xl x2
x3 8~1 u, P~ T~ Tw Yo2.~ YN:.,
Y.2.w YN2.w F = (pv)wl(pu)e
at
Value 63.5 cm 76.2 cm 118.1 cm 3.0734 cm 1024.13 cm/sec 1.0334 kg/cm 302.8°K 350.0°K 0.232 0.768 0.038 0.962 3.12 × 10-3 8.0
The predicted results of velocity profile at location x3 along the fiat plate is c o m p a r e d with the m e a s u r e d data b y Jones and I s a a c s o n (1970) as shown in Fig. 5. Also, the results of b o u n d a r y layer thickness & flame sheet position 3t and flame t e m p e r a t u r e TI are listed in Table 4. The c o m p a r i s o n s show that numerical results of this study are in good a g r e e m e n t with J o n e s ' experimental data.
Analysis of hybrid rocket combustion
629
I00 /
60
Prediction Experiment
\
~
80
\\
40
20
o
i
0
0.2
i
i
0.4
0.6
1.0
0.8
9 /r/e
Fig. 4. Comparison of computed and measured viscosity parameter. I 0
0.8
0.6
/
0.4
0.2
//~
Prediction Experiment -
_
_
0 0
1.0
2.0
3.0 y~
4.0
5.0
6.0
70
cm
Fig. 5. Comparison of computed and measured velocity profile.
Table 4. Comparison of computed and measured 6, 8t and Tf
Jones' experimental data Predicted results
8
8~
r~
6.731 cm 6.662 cm
1.524 cm 1.386 cm
1185°K 1184.9°K
Concluding remarks Based on this investigation for the hybrid rocket combustion, the proposed c o m b i n e d model of eddy viscosity has been developed to the stage where it can be used to predict quantitatively the behavior of compressible turbulent boundary layer with mass injection and combustion. H o w e v e r , the discrepancies b e t w e e n prediction and experimental data still exist in Figs. 4 and 5, which m a y be due to the flame sheet assumption. Further investigation is needed to improve the prediction.
630
D . L . Cherng and C. C. Tao
Acknowledgement--The co-author of this paper would like to express his sincere appreciation to Dr. H. A. Hassan, Professor, North Carolina State University, U.S.A., for guidance and encouragement during the course in part of this study. References Brokaw R. S. (1962) The Lewis number. Progress in International Research on Thermodynamic and Transport Properties, pp. 271-278. Academic Press, New York. Cebeci T. (1971) Calculation of compressible turbulent boundary layers with heat and mass transfer. AIAA J. 9, 1091-1097. Cebeci T. (1973) A model for eddy conductivity and turbulent Prandtl number. Z Heat Trans. (ASME) 95, 227-234. Cebeci T. and Smith A. M. O. (1974) Analysis of Turbulent Boundary Layers. Academic Press, New York. Chen T. N. and Toong T. Y. (1963) Laminar boundary layer wedge flows with evaporation and combustion. Prog. in Astron. Aeron. 15, 643--664. Chamber T. L. and Wilcox D. C. (1977) Critical examination of two-equation turbulence closure models for boundary layer. AIAA 3". 15, 821-828. Gany A., Manheimer-Timnat Y. and Wolfshtein M. (1976) Two-phase flow effect on hybrid combustion. Acta Astronautica 3, 241-263. Helman D. and Manheimer-Timnat Y. (1973) Investigation of hybrid rocket combustion. Israel Z Technol II, 17-25. Helman D., Wolfshtein M. and Manheimer-Timnat Y. (1974) Theoretical investigation of hybrid rocket combustion by numerical methods. Combustion and Flame 22, 171-190. Jones J. W. and Isaacson L. K. (1970) A Turbulent Boundary Layer with Mass Addition, Combustion, and Pressure Gradient. AFOSR Scientific Pep. 70-1428TR. Jones J. W., lsaacson L. K and Vreeke S. (1977) A turbulent boundary layer with mass addition, combustion and pressure gradient. AIAA J. 9, 1762-1768. Jones W. P. and Lauder B. E. (1972) The prediction of laminarization with a two-equation model of turbulence. Int. J. Heat and Mass Trans. 15, 301-314. Marxman G. and Gilbert M. (1963) Turbulent boundary layer combustion in the hybrid rocket. 9th Syrup. (Int.) on Combustion, pp. 371-383. Marxman G., Wooldridge C. E. and Muzzy R. J. (1964) Fundamental of hybrid boundary-layer combustion. Prog. in Astron. Aeron 15, 485-522. McBride B. J., Heimel S., Ehlers J. G. and Gordon S. (1963) Thermodynamic Properties to 6000 K for 210 Substances Involving the First 18 Elements. NASA-SP-3001. NG K. H. and Spalding D. B. (1972) Turbulence model for boundary layer near walls. Physics Fluids 15, 20-30. Patankar S. V. and Spalding D. B. (1970) Heat and Mass Transfer in Boundary Layer. lntertext Books, London. Schlichting H. (1968) Boundary Layer Theory. McGraw-Hill, New York. Saffman P. G. (1970) A model for inhomogeneous turbulent flow. Pro. R. Sac. A. 317, 417--433. Touloukian Y. S. Viscosity. TRPC Data Book, Vol. 2. Purdue University, Lafayette, Indiana. Van Driest E. R. (1956) On turbulent flow near a wall. J. Aeronaut. Sci. 23, 1007-1011. Wooldridge C. E. and Muzzy R. J. (1965) Measurements in a turbulent boundary layer with porous wall injection and combustion, lOth Syrup. (Int.) on Combustion, pp. 1351-1362. Wilcox D. C. and Traci R. M. (1976) A Complete Model of Turbulence. AIAA paper 76--351, San Diego, California.
Appendix
Nomenclature a thermal diffusivity, cm21sec A B constant in eqn (21)
c~ coefficient in eqn (20) cf local surface skin friction coefficient
Analysis of hybrid rocket combustion
[ F h H k k,. K Le rh" M
dependent variables mass injection rate, (pv)w/(pu), enthalpy, cal/mole total enthalpy, h + u2/2, cal/mole roughness height of the wall, cm Von Karman constant turbulent kinetic energy, cm2/secz Lewis number mass transfer rate across a boundary, kg/cm2 • sec average of molecular weight,
{X(YJM~)}-' n N P Pr Prf Pren R Re Scefr T
an exponent in eqn (20) defined in eqn (23) pressure, kg/cm z laminar Prandtl number turbulent Prandti number effective Prandtl number universal gas constant, call°K • kg-mole Reynolds number effective Schmidt number temperature, °K u,v x and y component of velocity, cm/sec u , dimensionless velocity, u/u, u, friction velocity, (rw/pw)tt2, cm/sec x,y rectangular coordinates, cm y,. value of y separating two model, cm y , dimensionless distance, yudv¢ Y species, Yo2-acYH2 Ya2 mass fraction of hydrogen Ys2 mass fraction of nitrogen Yo~ mass fraction of oxygen a,a* constants in Wilcox-Traci's two-equation model ac stoichiometric oxygen/hydrogen mass ratio
631
/3,/3* constants in Wilcox-Traci's two equation model /3c defined in eqn (20) 0 momentum thickness, cm stream function, kg/cm • sec tz laminar viscosity, kg/cm, sec /Lt turbulent viscosity, kg/cm, sec tze~ effective viscosity, tz +/~,, kg/cm • sec p density, kg/cm3 v kinematic viscosity, cm2/sec e,. eddy kinematic viscosity, e,, = Izt/p, cmZ/sec 8 boundary layer thickness cm 8r flame front location, cm or,or* constants in Wilcox-Traci's two-equation model cre~,r effective Prandtl number for variable / z shear stress, kg/cm 2 ~o dimensionless stream function viscosity parameter, defined in eqn (17) rt
"o = 0 , ~ , / ( 2 ~ )
,~ l~ =
'2
fo'
(pip,) dy
O~u~, dx, kg 2 seC/cm 4
Subscripts D e .f i o U w
downstream outer edge of boundary layer flame front inner region; species outer region upstream wall