- Email: [email protected]
Numerical modeling for the combustion of simulated solid rocket motor propellant
Computers & Fluids 89 (2014) 29–37
Contents lists available at ScienceDirect
Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d
Numerical modeling for the combustion of simulated solid rocket motor propellant A.M. Hegab a,⇑, H.H. Sait a, A. Hussain b, A.S. Said c a
Mechanical Engineering Department, Faculty of Engineering at Rabigh King Abdulaziz University, Saudi Arabia Department of Nuclear Engineering, King Abdulaziz University, Jeddah 21589, Saudi Arabia c Deanship of Scientific Research, King Abdulaziz University, Jeddah 21589, Saudi Arabia b
a r t i c l e
i n f o
Article history: Received 7 April 2013 Received in revised form 2 October 2013 Accepted 22 October 2013 Available online 30 October 2013 Keywords: Solid rocket propellants Moving interfaces Level set method Sandwich propellant AP/HTPB Alternating-Direction-Implicit (ADI) solver
a b s t r a c t The numerical procedure for the burning of Ammonium Perchlorate (AP) with a Fuel-Binder (Hydroxyl Terminated Polybutadience HTPB) heterogeneous propellant is presented. This model accounts for the two-step reaction mechanism for the primary diffusion flame at the interface between Binder (B) and the oxidizer AP and the AP premixed flame. The complete coupling between the gas-phase, the condensed-phase, and the unsteady non-uniform regression of the propellant surface is considered. The parameters used in this model are chosen to fit experimental data for the burning of the composite AP/Binder. The propagation of the unsteady non-planer regression surface is described, using the Essentially-Non-Oscillatory (ENO) scheme with the aid of the level set strategy. The Alternating-DirectionImplicit (ADI) solver is employed to solve the full Navier–Stokes equations in the gas phase for the variable density model. The results show the effect of various parameters on the surface propagation speed, flame structure, and the burning surface geometry. A comparison between the computational and experimental results is presented. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction The complex flame structure that is generated by burning of a heterogeneous solid rocket propellant is proposed by Beckstead, Derr, and Price (the BDP model [1]). Three separate flames can be identified in the gas phase: 1. a primary flame between the decomposition products of the binder and the oxidizer; 2. a premixed oxidizer flame; and 3. a final diffusion flame between the products of the other two flames. In spite of the BDP model is one-dimensional and necessarily omits or fails to properly account for important physics, but attempts to account for many of the significant feature of the combustion field. The influence of this work (published in 1970) still endures [2], and 1D models are still used [3]. Several improvements to BDP model of steady-state burning have been conducted. Lee et al. [4] presented a modified picture for the flame structure for AP-Binder-AP sandwich as in Fig. 1. This sketch show the principles of the combustion zone, in which the oxidizer-fuel flames consists of a leading-Edge Flame (LEF) that stands in the mixing region of the oxidizer and fuel vapors, and a diffusion flame that trails from the LEF up to a point where the fuel vapor is all consumed. The LEF is a region of very high heat release as com⇑ Corresponding author. Address: P.O. Box 344, Rabigh 21911, Saudi Arabia. Tel.: +966 548834995. E-mail addresses: [email protected], [email protected] (A.M. Hegab). 0045-7930/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compfluid.2013.10.029
pared to the rest of the diffusion flames and contributes most of the heat transfer back of the propellant surface. This edge occurs because the diffusion flame cannot extend all the way to the surface, the temperature there being too low. The theoretical studies for the combustion of heterogeneous solid rocket propellant have faced a lot of difficulties because of the chemical and physical complexity of the propellant and the microscopic scale of the combustion zone. Therefore, few experimental studies have been performed for the simplest model of the combustion of ammonium perchlorate sandwiches [4,5]. The propellant was made from sheets of AP-binder-AP. The AP formed by dry pressing ultra pure AP powder. Four different binders were used. Observations for the combustion were made by high-speed photography and microscopic examination of quenched samples. In addition, Lee et al. [4] illustrated the effect of inclusion of particulate AP in the binder on the combustion surface and the flame structure. The effect of three types of fuel binder and oxidizer particle diameter on the decomposition and combustion behavior of ammonium perchlorate is studies by Al-Harthi and Williams [6]. Few decades ago, several theoretical studies on the combustion field of the burning of the heterogeneous propellant have been conducted. These researches are divided into two main categories. The first one is concentrated on the gas phase modeling without consideration for the condensed phase process, for example [7– 11] The second one is studied the condensed phase reaction as
30
A.M. Hegab et al. / Computers & Fluids 89 (2014) 29–37
Nomenclature reaction rate constant overall regression speed specific heat diffusivity Damkohler numbers activation energies characteristic lengths Lewis number Mach number mass flux Peclet number Prandtl number pressure heats of reactions burning rate reaction rate Schmitt number
t T, X, Y
time temperature, mass fractions of AP and binder
Greek symbols b AP/binder stoichiometric ratio g surface function k heat conductivity m fractional binder thickness q density / equivalence ratio w function negative in binder, positive in AP Subscripts AP ammonium perclorate B binder g gas phase s solid phases
Periodic B. C.
Gas Phase Periodic B. C.
Vn
y Moving Interface at
AP
−ν
-L
t=0
x
rb
FUEL BINDER
the most important factor, for example [12,13]. Recently, several computational models [14–17] have been conducted to account the complex coupling between the solid-phase and gas-phase process. In particularly, the complexity that arises from the consideration of the unsteady non-planar regression surfaces. In the present paper a complete numerical strategy to examine what is perhaps, the simplest model is developed and account the following ingredients: the diffusion flame located at the interface between the decomposition products of the binder and the oxidizer (AP), the primary premixed AP flame, different properties (density, conductivity) of the AP and binder, temperature-dependent gas-phase transport properties, an unsteady non-planer regression surface; and a proper accounting of the fluid-mechanics in the gas-phase. These ingredients are applied to the problem of Periodic Sandwich Propellant (PSP).
-2L
AAP,B C Cp D Da1,2 E1,2 L Le M m Pe Pr P Qg1,2 rb Rr Sc
AP
ν
L
Fig. 2. Periodic sandwich propellant configuration.
2. The physical and mathematical models The physical model of the PSP model is shown in Fig. 2. This model consists of a sheet of fuel-binder of thickness ‘‘m’’, layered between two sheets of ammonium perchlorate (AP). The gas phase, consisting of a mixture of the decomposition products of the solid oxidizer and fuel is located above the combustion surface. Periodic boundary conditions are applied at x = ±L. The AP-binder-AP sandwich geometry has been chosen as a useful framework in which to gain fundamental insights into propellant combustion (e.g. [16–18]) and a notable experimental program has been pursued for some years by Price and his colleagues [5]. In addition, combustion behavior of the simpler sandwiches is much easier to observe and is enabled us to describe the combustion fields of the complex random packing propellant.
Diffusion flame
Leading edge flame (LEF) Stoichiometric surface
2.1. Constant and variable density models It is useful to summarize the formulation of the constant density model before addressing the complete problem, as this enable us to introduce most of the model ingredients together with various convenient scaling in the context of a model set of equations. The specific details of the constant density model for our problem are as follows: the density is set equal to constant (so that the equation of state, Charles law, is jettisoned); and a uniform velocity field u = 0 and v(y) = constant is adopted, which satisfies both the continuity and momentum equations. The one-step kinetics that include the primary flame is examined separately in order to achieve a good understanding of the unsteady burning of periodic sandwich propellant with complete coupling between the solid and gas phases. Thus R1
APðXÞ ! decomposition products ðZÞ R2
ð1Þ
Z þ binderðYÞ ! final products
Oxidizer-fuel mixing fan
R1 and R2 are assumed to have the forms;
AP flame Oxidizer (AP)
B
Oxidizer (AP)
Fig. 1. Flame structure for an AP-Binder-AP sandwich [4].
R1 ¼ B1 PX expðE1 =Ru TÞ and R2 ¼ B2 P ng YZ expðE2 =Ru TÞ; where B1 and B2 are the exponential prefatory, E1 and E2 are the activation energy in the gas phase, P is the pressure with exponent
31
A.M. Hegab et al. / Computers & Fluids 89 (2014) 29–37
ng, Ru is the universal gas constant, and (T, X, Y, and Z) are the temperature, oxidizer, fuel, and decomposition products respectively. The corresponding gas phase equations are;
equal to that in the gas phase for simplicity. The possibility of differing densities and thermal properties in the solid phase is allowed and setting by;
D/ ~ kg ~ qg ¼r r/ þ W Dt cp
qs ¼
ð2Þ
R1
R2
R1 bR2
T
When full fluid-mechanics coupling (variable density model) is accounted for, the system of Eq. (2) is replaced by;
@Q @F @G þ þ ¼H @t @x @y where
2
qg
3
ð3Þ
3 qg u 2 7 6 qg u þ P sxx 7 6 7 6 q u v s xy g 7 6 6 F ¼ 6 ðqg e þ PÞu ðusxx þ v sxy qx Þ 7 7; 7 6 qg uX qg Dg X x 7 6 5 4 qg uY qg Dg Y x qg uZ qg Dg Z x 3 2
qg v
7 6 qg uv sxy 7 6 2 7 6 q v þ P s yy g 7 6 6 ¼ 6 ðqg e þ PÞv ðusxy þ v syy qy Þ 7 7; 7 6 qg v X qg Dg X y 7 6 5 4 qg v Y qg Dg Y y qg v Z qg Dg Z y and
e ¼ qP þ 12 ðu2 þ v 2 Þ; g
ks ¼
kAP
wP0
kB
w<0
ð7Þ
;
gt þ gx
dx dy þ gy ¼ 0: dt dt
ð8Þ
and the final equation that controls the moving of the gas/solid interface as in Fig. 3 is derived by Hegab et al. [16,17] and may be written as follows;
~ gj ¼ 0; gt rb jr 2
6 qg u 7 7 6 6 qg v 7 7 6 6 Q ¼ 6 qg e 7 7; 6 qg X 7 7 6 4 qg Y 5 qg Z 2
qAP qB
The function w(x, y) is a level set function which demarks the regions of AP from binder (B) within the solid, so that a point (x, y) lies in the AP if w(x, y) P 0, and in the binder if w(x, y) < 0. Suppose the solid/gas interface defined by g(x(t), y(t), t) = 0. Then;
where
/ ¼ ½ T X Y Z T ; W ¼ ðQ g1 R1 þ Q g2 R2 Þ=cp
sxx ¼ l
where rb is defined as the speed of the front which moves in the directions of the solid. In general rb is a function of x and t and is given by the following simple pyrolysis law;
(
G rb ¼
3 0 7 6 0 7 6 7 6 0 7 6 6 H ¼ 6 þQ g1 R1 þ Q g2 R2 7 7 7 6 R1 7 6 5 4 R2 þR1 bR2
4
u 23 v y ; 3 x
ð9Þ
4
syy ¼ l 3 v y 23 ux ;
sxy = l(vx + uy); qx = kgTx; qy = kgTy Lewis number is taken to be unity, then;
qg Dg ¼ kg =cp
ð10Þ
Note that pressure dependence has been added to the pyrolysis law for generality. In the study, the propellant surface is not flat and its shape changes with time. Therefore, the following mapping function is used;
g ¼ y f ðx; tÞ
ð11Þ
and the front of Eq. (9) reduces to the simple Hamilton–Jacobi equation;
qffiffiffiffiffiffiffiffiffiffiffiffiffi ft þ r b ðx; tÞ 1 þ fx2 ¼ 0;
ð12Þ
ð5Þ
2.3. Boundary/jump conditions
2.2. Solid-phase and solid/gas interface equations
The appropriate jump conditions across the gas/solid interface are;
½qð~ v ~n þ rb Þ ¼ 0;
ð13Þ
½T ¼ 0; ½k~ n rT ¼ Q s m;
In the solid-phase, the following heat equation is used;
m½vi ¼ ½qD~ n rvi ; ð6Þ
Here, qs is the density of the solid, T the temperature, and ks is the solid thermal conductivity. The specific heat cp is assumed to be
ð14Þ i ¼ 1; 2; 3;
ð15Þ ð16Þ
where [/] = /g /s denotes the jump in the quantity / across the interface, vI refers to X, Y, and Z, m is the mass flux. ~ n is the unit normal pointing in the direction of the gas; ~ n ¼ rg=jrgj, Fig. 3. Qs is the solid phase heat release term defined by
Qs ¼
Fig. 3. Sketch showing the coordinate system for the moving surface.
w<0
Further information about the non-planar moving of the gas/solid interfaces using the level set strategy is mentioned in details in [16,17].
where ax and ay are the stoichiometric coefficients. Here there are six unknowns in the gas-phase, (u, v, T, X, Y) and one unknown in the solid-phase (Ts)
k qs T t ¼ s r2 T cp
r B ¼ AB ðP=P0 ÞnB expfEB =Ru T B;s g
ð4Þ
With the aid of the equation of state;
P ¼ qRT
r AP ¼ AAP ðP=P0 ÞnAP expfEAP =Ru T AP;s g w P 0
Q AP
wP0
Qb
w60
ð17Þ
For an exothermic surface reaction, Q s > 0, and for an endothermic reaction, Q s < 0. Typically the AP is considered an exothermic reaction, while the binder an endothermic one. Studies by Hegab et al. [16,17] proved that the length and time scales for the front and the solid are the same order of magnitudes. On the other side, the ratio of the gas to solid or the ratio of the gas to front is of the order of 103. Thus for the present purpose, the quasi-steady approximation for the gas phase is employed. Note
32
A.M. Hegab et al. / Computers & Fluids 89 (2014) 29–37
that disturbances with time scales of order 103 s would effect the solid phase, but not the gas phase; changes on time scales of order 105 s are needed to generate an unsteady gas phase and changes of this nature have been discussed in [11]. 2.4. Nondimensionalization
2
t s ¼ L=r b;ref ;
t g ¼ L=V g ;
e ¼ tg =ts 1 V g ¼ qAP rb;ref =qg ; T ref ¼ 2700 K; Q ref ¼ cp 2700 j=g
Q s;AP =ðcp T ref Þ w P 0 Q s;B =ðcp T ref Þ
kratio ¼
qratio ¼
w < 0;
1
wP0
kB =kAP
w < 0;
1
wP0
qB =qAP w < 0:
3 0 7 6 aux 7 6 7 6 4=3 a v x 7 6 2 7 V 2 ðQ ; Q y Þ ð c 1ÞM a fuu þ 4=3 vv g þ b T W 2 ðQ ; Q y Þ ¼ 6 y y y7 6 7 6 l =P X e y 7 6 5 4 l=Pe Y y l=Pe Z y 3 2 0 7 6 2=3av y 7 6 7 6 u y 7 6 2 6 ¼ 6 ðc 1ÞM af2=3uv y þ v uy g 7 7W 1 ðQ ; Q x Þ 7 6 0 7 6 5 4 0
ð18Þ
ð19Þ
ð20Þ
ð21Þ
7 6 av x 7 6 7 6 2=3 a u x 7 6 2 7 ð c 1ÞM a fu v 2=3 v u g ¼6 x x 7 6 7 6 0 7 6 5 4 0 0
ft þ r b
2
R2 ¼ Da2 P YZ expðhg2 =TÞ; but temperature-dependent transport is accounted for, viz; kg ¼ kg;ref ~ kðT; T ref Þ where kg,ref is a reference heat conduction coefficient. The value of k at the reference temperature Tref, specially kg = 1.08104 T + 0.0133 is choosing with dimensions W/m K when T is assigned in degrees Kelvin, so that
u ¼ 0;
1:08 104 T ref þ 0:0133
v ¼ qratio ft ;
R ¼ Da PXY expðhg =TÞ
ð22Þ ð23Þ
The non-dimensional equations of motion can be re-written in the following form:
@Q @FðQ Þ @GðQ Þ @V 1 ðQ; Q x Þ @V 2 ðQ ; Q y Þ þ þ ¼ þ @t @x @y @x @x @W 1 ðQ ; Q x Þ @W 2 ðQ ; Q y Þ þ þ þH @y @y
where
ð25Þ
qffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ fx2 ¼ 0;
ð26Þ
Tðx; 0þ ; tÞ ¼ Tðx; 0 ; tÞ
ð27Þ
1
2 f T þ 1 þ f T x x g C
x ^kB qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi @ A
1 þ fx2
0
0þ
g > 0 for Variable Density Model
e
where a ¼ PSce l; b ¼ PSc k, and h ¼ e þ c1 c p=q e Pr g < 0 (Solid Phase)
g = 0 (Moving Interface)
R1 ¼ Da1 PX expðhg1 =TÞ
~kðT; T ref Þ ¼
3
0
qratio ðT t ft T g Þ ¼ ðkratio =Pec ÞDT
R1 and R2 have the following forms;
1:08 104 T ref T þ 0:0133
0
2
In non-dimensional form the equations and boundary/connections form conditions are: g > 0 for Constant Density Model
~ ð Kr ~ TÞ þ Q g1 R1 þ Q g2 R2 eT t þ uT x þ v T g ¼ ð1=Peg Þr ~ ~ XÞ R1 eX t þ uX x þ v X g ¼ ð1=Peg Þr ðKr ~ ðKr ~ YÞ R2 eY t þ uY x þ v Y g ¼ ð1=Peg Þr ~ ðKr ~ ZÞ þ R1 bR2 eZ t þ uZ x þ v Z g ¼ ð1=Peg Þr
3
2
Pressure po (bar), surface speed rb,ref (cm/s), and mass flux mref = qAPrb,ref. Length L (half of the computational domain, which is the sum of the binder and the AP thickness). Time t = L/rb,ref. Then the following non-dimensional parameters are defined: Peclet numbers Peg = qgVg L cp/kg,ref, and Pec = qAP rb,ref L cp/kAP, Activation energy h = E/(RuTref)
bs ¼ Q
2
3 0 7 6 4=3aux 7 6 7 6 a v x 7 6 2 6 ¼ 6 ðc 1ÞM af4=3uux þ vv x g þ b T x 7 7 7 6 l=Pe X x 7 6 5 4 l=Pe Y x l=Pe Z x
T ¼ T=T ref ; X ¼ X=X s ; Y ¼ Y=Y s ; Z ¼ Z=Z s ; P ¼ P=Po ; q ¼ q=qo ; ðu ; v Þ ¼ ðu; v Þ=V g ; f ¼ f =L; rb ¼ rb =rb;ref ; t ¼ t=t s ;
3
2
2
The following reference values are taken to nondimensionalize the equations;
ðx ; g Þ ¼ ðx; gÞ=L;
3
q qu qv 6 qu 7 6 qu2 þ p 6 quv 7 ~7 7 7 7 6 6 6 2 6 qv 7 6 quv 7 6 qv þ p ~7 7 7 7 6 6 6 7 7 7 6 6 Q ¼6 6 qe 7 FðQ Þ ¼ 6 quh 7 GðQÞ ¼ 6 qv h 7 V 1 ðQ ; Q x Þ 6 qX 7 6 quX 7 6 qv X 7 7 7 7 6 6 6 4 qY 5 4 quY 5 4 qv Y 5 qZ quZ qv Z
ð24Þ
1
2 f T þ 1 þ f T x x g B C
x ^ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kkratio @ A
1 þ fx2
0
b s mP e ¼ Q g
ð28Þ
8 91
>
> k^ <fx X x þ 1 þ fx2 X g =C
B qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi @mX A
> P eg > : ;
1 þ fx2 0
8 91
> > < 2 ^ k fx Y x þ 1 þ fx Y g =C
B q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi mY @ A
> P eg > : ;
1 þ fx2
¼
m wP0 0
w<0
0
wP0
ð29Þ
0þ
0
8 91
> > < 2 ^ k fx Z x þ 1 þ fx Z g =C
B qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi
@mZ A > > P eg : ;
1 þ fx2
0
¼
m w<0
0þ
0
¼ 0þ
0
wP0
m w<0
ð30Þ
A.M. Hegab et al. / Computers & Fluids 89 (2014) 29–37
g ! 1 T ¼ T o g ! þ1
ð31Þ
@ ðÞ ¼ 0 @g
ð33Þ
3.1. Modified conservation equations To avoid pressure singularity at low Mach numbers, the numerical strategy outlines in [20,21] is used. The pressure is rescaled in the momentum equation since it is the pressure gradient, not the actual pressure that is involved in the momentum balance. The rescaled pressure is applied to retain accuracy in calculating the momentum conservation. As a result, the pressure is divided into constant and fluctuating parts as follow;
b y; tÞ Pðx; y; tÞ ¼ 1 þ cM 2 Pðx;
ð34Þ
and substitute into the equations of motion given above. Where M is the Mach number (M = Vo/Co). The equation of state now becomes;
b y; tÞ 1 þ cM 2 Pðx; T
ð35Þ
and is used to update the density. Since the quasi-steady state solution in the gas phase (e 1) is required, the physical time term in (24) is eliminated. In addition an accelerate convergence technique is employed [20]. The technique begins by adding a pseudo-time derivative to the conservation Eq. (24). So the equations to be solved for the gas phase become;
e @FðQÞ @GðQ Þ @V 1 ðQ ; Q Þ @V 2 ðQ ; Q y Þ @Q x þ þ ¼ þ @x @x @y @x @s @W 1 ðQ ; Q x Þ @W 2 ðQ ; Q y Þ þ þ þH @y @y
ð36Þ
where s represents the pseudo-time domain. Because the pseudotime derivative vanishes at convergence to the steady state solution in the gas phase, a certain amount of liberty is given in choosing the e . A new scaled pressure term P=b e is added as a pseuvariables in Q do-time derivative term to the continuity equation. The remaining e are then fixed by rewriting the momentum, energy, variables in Q and species equations from their neoconservative form by means of the modified equation of mass conservation. As a result, the pseue and its associated preconditioning matrix do-time variable vector Q C take the form;
2 3 2 b 1=b0 0 0 0 0 0 P 6 7 6 0 q 0 0 0 0 6u7 6 u=b 6 7 6 6 7 6 v =b0 0 q 0 0 0 6v 7 6 6 7 6 e ¼ 6 h 7C ¼ 6 h=b0 1 0 0 q 0 0 Q 6 7 6 6 7 6 6X 7 6 X=b0 0 0 0 q 0 6 7 6 6 7 6 0 0 0 0 0 q 4Y 5 4 Y=b Z
Z=b0
0
0
0
0
0
0
3
7 07 7 07 7 7 07 7 7 07 7 7 05
ð37Þ
q
b Its selecFactor b’ is used to properly scale the time derivative of P. tion is proportional to the dynamics pressure of the flow field. 3.2. ADI algorithm development The first step for advancing the solution of (36) is using the Delta form scheme [19] as follows;
ð38Þ
and
DQ n ¼ Q nþ1 Q n
3. Numerical methods
C
DQ n ¼
ð32Þ
jxj ¼ 1 Periodic boundary conditions
q¼
hDs @ Ds @ n ð DQ n Þ þ ðQ n Þ þ DQ n1 1þn 1 þ n @s 1 þ n @s þ O½ðh 1=2 nÞDs2 þ Ds3
33
ð39Þ
The time-difference formulas (38) and (39), with the appropriate choice of the parameters n and h reproduce many familiar two and three level, implicit schemes. The three level, second order accuracy implicit schemes (n = 1/2 and h = 1) with special treatment for the cross derivative terms at the level (n 1) is applied. By inserting the temporal derivative of Eq. (36) in (38) and by straightforward derivation, the resulting approximated form can be splitted into two-tridiagonal systems;
" I þ hDsC1 " I þ hDsC1
!# @ @2 ðA P þ Rx Þn 2 ðRn Þ DQ ¼ C1 RHS @x @x !# 2 @ e þ Sy Þn @ ðSn Þ DQ ¼ DQ ðB Q @y @y2
ð40Þ ð41Þ
e n ; Sn , and Rn are the jacobian matrices for where An ; Pn ; Rn ; Rnx Bn ; Q y the implicit x- and y-sweep (40) and (41) respectively. I am a unit matrix (6 6) and C1 am the inverse of the preconditioning matrix. The three-points second order accuracy central difference approximation is used for the spatial differences in the LHS of x-sweep form (40). This approximation along with the applied periodic boundary conditions at jxj = 1 produce a system of Periodic Block-Tridiagonal Equations (PBTE). After the computation of DQ⁄ at the interior points by solving this system of PBTE with each block having dimensions 6 6 components, the code is ready for the implicit y-sweep form (41). Here again, the three-point second order accuracy central difference approximation are used for the spatial differences in the LHS of (41). This approximation along with the rigid wall boundary conditions at y = 0 produce a system of NonPeriodic Block-Tridiagonal Equations (NPBTE). The final delta form DQ can be computed by solving the NPBTE system. Then the solution at new time step (n + 1) can be determined from Eq. (40). Note that, the cross-derivatives terms at (n 1) are treated explicitly to avoid the implicit coupling of adjacent boundary points. The surface Eq. (12) is solved in order to follow the non-planar regression surface by the first order temporal scheme [22]. Besides the non-flat regression surface mapping as in (10), another transformation is applied for clustering grid points adjacent to the wall, where most of the flow parameters changes rapidly. The solution of the final mapped equations is advanced in the solid phase using physical time (t). Simultaneously the solution in the gas phase using pseudo-time (s) to the local steady state at the first physical time step (t) is advanced. The boundary/jump conditions are continually updated as in Eqs. (25)–(33). Then the Hamilton–Jacobi Eq. (12) is advanced at the physical time by a third order ENO and a fifth-order WENO (weighted essentially non-oscillatory) solver [23,24]. All numerical calculations were performed on a 140 70 grid, uniform in the x-direction and stretched in the y-direction. At each physical time, the solution in the gas phase is advanced until the relative difference between each two different pseudo-time values is less than some prescribed tolerance, taken here to be 106. Convergence tests where carried out and it was determined that any further refinement resulted in less than 1% relative error.
34
A.M. Hegab et al. / Computers & Fluids 89 (2014) 29–37
Fig. 4. Reaction rate contours at different times in the gas-phase for hg1 = 3.5, hg2 = 6.0, b = 7.51, Peg = 6, P1s = 2.9, Da1 = 5104, Da1 = 1.7105, m = 0.22, Qg1 = 0.833, Qg2 = 3.86, Qs,AP = 0.4, and Qs,B = 0.06.
4. Results and discussion The understanding of the complex combustion structure of APBinder-AP sandwich, as a simple model to the heterogeneous solid rocket propellant, is studied in details by two different models. The first model is the constant density model and the second one is the variable density model or the Navier–Stokes model. Initially, the solution starts from a flat surface f(x, t = 0) = 0. Then the solution is advanced in the solid phase, gas phase with simultaneously moving interface. The first set of the results are for the constant density model. In this model, it has been taken l = K(T, Tref) which is one of the more realistic choices rather l = constant in earlier studies [7].
In the gas phase, the contours plots for the reaction rate contours (R) at different times are shown in Fig. 4 for hg1 = 3.5, hg2 = 6.0, b = 7.51, Peg = 6, P1s = 2.9, Da1 = 5104, Da1 = 1.7105, m = 0.22, Qg1 = 0.833, Qg2 = 3.86, Qs,AP = 0.4, and Qs,B = 0.06, as benchmark computation. The upper portion represents the gas phase and the lower one refers to the solid phase. The dark gray region in the latter represents the binder layer between the two AP sheets (light gray). The contour values are written in the upper part of each plot. The plot shows the location and shape of the generating flame. These reaction rate contours show a combined of two different flames. The first one is the diffusion flame that formed at the interface between the oxidizer and the fuel, while the second one is the premixed flame that formed above the oxidizer regions.
Burning surface
Stoichiometric Contour Fig. 5. Comparison of experimental image for the burning of AP-Binder-AP sandwich[26] {left} with the numerical model {right}.
A.M. Hegab et al. / Computers & Fluids 89 (2014) 29–37
Fig. 6. Variations in the global regression (c) rate with pressure (p).
The contours of the diffusion flame show that reaction rate is characterized by two strong mixing structure each centered at jxj 0.223. The successive curved shapes through the solid phase show that, the surface is initially flat and then as the solution is advanced, the combustion surface retreats in an unsteady fashion and finally retreats at a fixed speed with unchanged shape by t = 2.49. These contours reveal not only the significant effect of the surface profiles and the burning rate on the shape of the flame but also on the maximum reaction rate values which decrease as time increase. The differences in profiles with time advanced reflect the behavior of the burning rate at the propellant surface and reconstruct the corresponding flame structures and identify the parts of the flame structure that dominate the sandwich burning rate and the surface heat flux. In addition there is a sharp discontinuity in slope at the interface between the binder and AP regions. This interesting phenomenon is predicted in [25]. Fig. 5 shows a comparison between the experimental image (the left) [26] and the current computation model (the middle for the reaction rate and right for the fuel vapor contours) for the structure of the flame shape and the burning surface. The stoichiometric level surface shows the location at which the fuel and oxidizer meet together in shoichiometric proportions. The inner
35
region for the stoichiometry contour is the fuel-rich region (see the fuel vapor image on the left side), while the area out of the stoichiometry envelop represent the fuel-lean composition. The two stoichiometric surface in both the computational and experimental images represent the two strong mixing structure (at the AP-Binder interface) each centered at jxj 0.223 and forming the flame envelop in addition to the AP-premixed flame that control the shape of the combustion surface above the AP region. This comparison show how the computational steadily surface regressing profile (at t > 5) is qualitatively consistent with the experimental emission-transmission composite image by Brewster et al. [26], [26] and also with the experimental photography for the spontaneous quench samples of AP/BPAN/AP sandwich that reported by Lee et al. [4]. Both theoretical and experimental results verified that conditions that gave the highest burning rates resulted in narrow smooth bands of binder in the middle with little ‘‘protrusion’’ of AP at the interface. The variation of the burning rate c with the pressure is presented in Fig. 6 for the benchmark data as in Fig. 4. The response of the burning rate to the exposed pressure is one of the important and essential characteristics of the propellant. This figure shows that the burning rate is directly proportional to the pressure with exponential component about (n = 0.6). This exponential value has great interest to the rocket designer, in particular for the modeling of the composite propellant rather the sandwich one. The effect of the exposed pressure on the burning surface profiles is presented in Fig. 7 at three pressure values P = 1, 10, 15 atm in the gas phase for the same parameters in Fig. 4. A comparison between the experimental image and the computational model for the response of the combustion surface profiles to the pressure change (experimental at 7 atm in the left and the computational at 1, 10, and 15 atm in the right) is presented in Fig. 7. Both images illustrate that increasing the pressure causes the binder to stick out and the surface consist of a ‘‘trough’’ centered on the binder lamina. In addition a comparison of these clear pictures with the experimental results by Price and et al. [5] is possible but caution is appropriate as the sandwiches that Price examined are isolated, not periodic. Furthermore, the effect of the pressure on the surface heat flux reveal that the heat flow vector, near the gas–solid interface and AP-binder solid interface, is from the binder toward the AP in the solid phase. Consequently there is a weak (hot) portion of the
Fig. 7. Effect of pressure on the burning surface structure. Left figure is experimental [5] and the right figures are the computational model.
36
A.M. Hegab et al. / Computers & Fluids 89 (2014) 29–37
Fig. 8. Total reaction rate contours (a) constant density model and (b) variable density model.
Fig. 9. The axial velocity component for the variable density model (2) at different times for h = 7, b = 7, Pe = 6, Dag = 5106, P = 10, M = 0.02.
binder slightly up from the corner interface and into the binder. This might account for the appearance of the V-shape (notches) in the binder. These phenomena have been seen in the experimental quenched samples by Handley and Strahle [25]. The last set of the results are for the variable density model (2). In this section the results obtained from the constant density model (1) with those obtained using the Navier–Stokes equations, model (2) are compared, and thus to validate the simpler strategy. Fig. 8 shows a comparison between the two models for the reaction rate contours for the parameters in Fig. 4. The reaction rate values
varied from 12 to 2 with incremental 2 for the AP-premixed flame, while varied from 1 to 0.2 with incremental 0.2 for the AP-Binder diffusion flame that centered near the interface between the fuel and oxidizer biased little bite toward the oxidizer surface. It is noted that the flame is setting at the same transverse location with slight differences in the maximum reaction rate contours. RRmax is 11.88678 for model (2), compared to 11.66972 for model (1) and in the adiabatic flame temperature, where T 0max is 2448.09 K for model (2) compared to 2446.57 for model (1). These interesting values reveal that there is only a small difference in the far-field for the
A.M. Hegab et al. / Computers & Fluids 89 (2014) 29–37
two models. These very slight differences are related to the slow flow accelerated since the axial velocity component u is not zero when the Navier–Stokes equations are used and the transverse velocity component becomes more or less in the neighborhood of the flame for model (2). Fig. 9 shows the axial velocity contours at different times for the Navier–Stokes model. It is noted that an axial velocity is generated near the combustion surface due to the surface morphology, where the velocity cells are positive on the left hand side, since the flow goes down hill and the right hand side cells are negative, since the flow goes up hill. Moreover, the absolute value of the axial velocity increases as time increases, since the curvature in the combustion surface profiles become more deeper than at earlier time. The small differences between the two models reveal that, a useful exploration calculations can be carried out using the constant density model, since the generated axial velocity component is very small, 0.0 6 juj P 0.04 and also the variation in v is consistent with the mass conservation. In general, the Navier–Stokes model may have a significant effect when the current solution to the modeling of randomly packed heterogeneous propellant, particularly in 3-D solution is advanced. 5. Conclusion Few decades ago the details of the heterogeneous propellant combustion has been discussed. The theoretical framework has consists of a 2D physical pictures, but most modeling efforts have been 1D mathematical model as in [1–4]. This has obviously been of value and important insights have been achieved, but a multidimensional numerical framework is required. Here, for the first time, the 2D calculations to the combustion of heterogeneous solid propellant, accounting for the gas phase physics, the solid phase physics and an unsteady non-planar description of the regressing propellant surface is developed. There are a number of issues that have discussed. The speed within which the combustion surface recedes depends on the exposed pressure in the gas phase, the effect of several parameters on the combustion and shape of the flame. In addition a variety of steady-state surface shapes are achieved. At higher pressure values, the binder is tending to stick out and the surface consists of a ‘‘trough’’ centered on the binder lamina. These trends were also recognized in an experimental observation [5]. A comparison between the computational steadily surface regressing profile with the experimental emission-transmission composite image by Brewster et al. [26], [24] and also with the experimental photography for the spontaneous quench samples of AP/BPAN/AP sandwich that reported by Lee et al. [4] is qualitatively consistent. In this study, the Navier–Stokes equations were used rather the constant density model in earlier work [7,10,16,18]. A comparison between the constant density model and the Navier–Stokes solutions reveals very small differences. More intensive work for the solution of the full Navier–Stokes as in [27–32] is required to perform more complex computation for the reactive whole scale of the solid rocket motor chamber. Acknowledgments This work was funded by the Deanship of Scientific Research (DSR), King Abdulaziz University, Jeddah, under Grant No. (829003-D1433). The authors, therefore, acknowledge with thanks DSR technical and financial support. The first author would like to acknowledge Prof. John Buckmaster and Prof. Tomas Jackson at University of Illinois at Urbana–Champaign (UIUC), USA for their support for developing the variable density model.
37
References [1] Beckstead MW, Derr RL, Price CF. A model of composite solid-propellant combustion based on multiple flames. AIAA J 1970;8(12):2200–7. [2] Rasmussen Jr B, Frederick RA. A nonlinear heterogeneous model of composite solid propellant combustion. In: AIAA paper no. 99-2228. 35th AIAA/ASME/ SAE/ASEE joint propulsion conference; June 1999. [3] Murphy JJ, Krier H. Heterogeneous effects on dynamic burning in composite solid propellant. In: Proceeding of the combustion institute, vol. 28; 2000. [4] Lee S, Price E, Sigman R. Effect of multidimensional flamelets in composite propellant combustion. J Propul Power 1994;10(6):761–8. [5] Price E, Handley J, Panyam R, Sigman R, Ghosh A. Combustion of ammonium perchlorate-polymer sandwiches. AIAA J 1981;19(3):380–6. [6] Al-Harthi A, Williams A. Effect of fuel binder and oxidizer particle diameter on the combustion of ammonium perchlorate based propellants. Fuel 1998;77(13):1451–68. [7] Jackson TL, Buckmaster J, Hegab A. Periodic propellant flames and fluidmechanical effects. J Propul Power 2001;17(2):371–9. [8] Buckmaster J, Jackson TL, Yao J. An elementary diffusion of propellant flame geometry. Combust Flame 1999;117:541–52. [9] Jackson TL, Buckmaster J, Hoeflinger J. Three-dimensional flames supported by heterogeneous propellants. In: JANNAF paper, JANNAF CS/PSHS/APS joint meeting, Cocoa Beach, FL; October 1999. [10] Jackson TL, Buckmaster J. Nonpremixed periodic flames supported by heterogeneous propellant. J Propul Power 2000;16(3):498–504. [11] Cai W, Yang V. A model of AP/HTPB composite propellant combustion. In: AIAA paper 2000-0311, 38th aerospace science meeting; 2000. [12] Waesche RH, Wenograd J. Calculation of solid propellant burning rate from condensed-phase decomposition kinetics. AIAA paper 69-145, NY; January 1969. [13] Hermance CE. A model of composite propellant combustion including surface heterogeneity and heat generation. AIAA J 1960;4:1629–37. [14] Hegab AM. An overview of solid propellant combustion. In: Proceeding of the 12th aerospace science & aviation technology conference, (ASAT-12), 29–31 May, 2007, BAL-03, Military Technical College, Cairo, Egypt; 2007. p. 1–24. [15] Hegab AM, Balabel A. Effect of ammonium perchlorate grain sizes on the combustion of solid rocket propellant. In: Proceeding of the 12th aerospace science & aviation technology conference, (ASAT-12), 29–31 May, 2007, CHA1, Military Technical College, Cairo, Egypt; 2007. p. 1–20. [16] Hegab A, Jackson T, Buckmaster J, Stewart S. The burning of periodic sandwich propellants. In: 36th AIAA/ASME/SAE/ASEE joint propulsion conference, AIAA paper 2000-3459; 2000. [17] Hegab A, Jackson T, Buckmaster J, Stewart S. Nonsteady burning of periodic sandwich propellant with complete coupling between the solid and gas phases. Combust Flame 2001;125(1/2):1055–70. [18] Knott GM, Brewster MQ. Two-dimensional combustion modeling of heterogeneous propellants with finite Peclet number. Combust Flame 2000;121(1/2):91–106. [19] Tseng IS, Yang V. Combustion of a double-base homogeneous propellant in a rocket motor. Combust Flame 1994;96:325–42. [20] Roh TS, Yang V. Transient combustion responses of solid propellants to acoustics disturbances in rocket motors. AIAA paper 95-0602; 1995. [21] Beam RM, Warming RF. An implicit factored scheme for the compressible Navier–Stokes equations. AIAA J 1978;16(4). [22] King MK. Examination of chemical approaches to stabilizing composite propellant combustion. J Propuls Power 1996;12:554–63. [23] Harten A. Uniformly high order accurate essentially non-oscillatory schemes, III. J Comput Phys 1987;71:231–303. [24] Jiang G, Shu C. Efficient implementation of weighted ENO scheme. J Comput Phys 1996;126:202–28. [25] Handley JC, Strahle WC. The behaviour of several catalysts in the combustion of solid propellant sandwiches. AIAA paper 74-122, Washington (DC); 1974. [26] Brewster MQ, Knott GM, Chorpening BT. Combustion of AP/HTPB laminate propellants. In: AIAA paper 2001-4501, 37th AIAA/ASME/SAE/ASEE joint propulsion conference, Salt Lake City, UT, USA; 2001. [27] Hegab AM, Sait HH, Hussien A, Balabel A, Almutawa M, Djouider F. Modeling of microscale solid rocket composite propellant. Life Sci J 2013;10(4). [28] Hegab AM, Hasanien SA, Mostafa HE, Maraden AM. Experimental and numerical investigations for the combustion of selected composite solid propellant. Eng Res J 2012;35(3). ISSN:1110-1180. [29] Maraden A, Hossam E, Hasanien S, Hegab A. Effect of packing parameters on the combustion behaviour of ap based composite solid propellant. In: 6th international conference on chemical & environmental engineering: energetic material I, Military Technical College; 29–31 May, 2012. [30] El-Askary WA, Balabel A, El-behery SM, Hegab AM. Computational of a compressible turbulent flow in a rocket-motor chamber configuration with symmetric and asymmetric injection. J Comput Model Eng Sci 2011;82(1):29–54. [31] El-Askary WA, Wilson SA, Hegab AM. Simulation of reactive fluid flow in a solid rocket motor combustion-chamber with/without nozzle. J Comput Model Eng Sci 2011;76(4):235–66. [32] Balabel A, Hegab A, Nasr M, El-Behery S. Assessment of turbulence models for gas flow in convergence–divergence rocket nozzle. Appl Math Model 2011;35:3408–22.
Contents lists available at ScienceDirect
Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d
Numerical modeling for the combustion of simulated solid rocket motor propellant A.M. Hegab a,⇑, H.H. Sait a, A. Hussain b, A.S. Said c a
Mechanical Engineering Department, Faculty of Engineering at Rabigh King Abdulaziz University, Saudi Arabia Department of Nuclear Engineering, King Abdulaziz University, Jeddah 21589, Saudi Arabia c Deanship of Scientific Research, King Abdulaziz University, Jeddah 21589, Saudi Arabia b
a r t i c l e
i n f o
Article history: Received 7 April 2013 Received in revised form 2 October 2013 Accepted 22 October 2013 Available online 30 October 2013 Keywords: Solid rocket propellants Moving interfaces Level set method Sandwich propellant AP/HTPB Alternating-Direction-Implicit (ADI) solver
a b s t r a c t The numerical procedure for the burning of Ammonium Perchlorate (AP) with a Fuel-Binder (Hydroxyl Terminated Polybutadience HTPB) heterogeneous propellant is presented. This model accounts for the two-step reaction mechanism for the primary diffusion flame at the interface between Binder (B) and the oxidizer AP and the AP premixed flame. The complete coupling between the gas-phase, the condensed-phase, and the unsteady non-uniform regression of the propellant surface is considered. The parameters used in this model are chosen to fit experimental data for the burning of the composite AP/Binder. The propagation of the unsteady non-planer regression surface is described, using the Essentially-Non-Oscillatory (ENO) scheme with the aid of the level set strategy. The Alternating-DirectionImplicit (ADI) solver is employed to solve the full Navier–Stokes equations in the gas phase for the variable density model. The results show the effect of various parameters on the surface propagation speed, flame structure, and the burning surface geometry. A comparison between the computational and experimental results is presented. Ó 2013 Elsevier Ltd. All rights reserved.
1. Introduction The complex flame structure that is generated by burning of a heterogeneous solid rocket propellant is proposed by Beckstead, Derr, and Price (the BDP model [1]). Three separate flames can be identified in the gas phase: 1. a primary flame between the decomposition products of the binder and the oxidizer; 2. a premixed oxidizer flame; and 3. a final diffusion flame between the products of the other two flames. In spite of the BDP model is one-dimensional and necessarily omits or fails to properly account for important physics, but attempts to account for many of the significant feature of the combustion field. The influence of this work (published in 1970) still endures [2], and 1D models are still used [3]. Several improvements to BDP model of steady-state burning have been conducted. Lee et al. [4] presented a modified picture for the flame structure for AP-Binder-AP sandwich as in Fig. 1. This sketch show the principles of the combustion zone, in which the oxidizer-fuel flames consists of a leading-Edge Flame (LEF) that stands in the mixing region of the oxidizer and fuel vapors, and a diffusion flame that trails from the LEF up to a point where the fuel vapor is all consumed. The LEF is a region of very high heat release as com⇑ Corresponding author. Address: P.O. Box 344, Rabigh 21911, Saudi Arabia. Tel.: +966 548834995. E-mail addresses: [email protected], [email protected] (A.M. Hegab). 0045-7930/$ - see front matter Ó 2013 Elsevier Ltd. All rights reserved. http://dx.doi.org/10.1016/j.compfluid.2013.10.029
pared to the rest of the diffusion flames and contributes most of the heat transfer back of the propellant surface. This edge occurs because the diffusion flame cannot extend all the way to the surface, the temperature there being too low. The theoretical studies for the combustion of heterogeneous solid rocket propellant have faced a lot of difficulties because of the chemical and physical complexity of the propellant and the microscopic scale of the combustion zone. Therefore, few experimental studies have been performed for the simplest model of the combustion of ammonium perchlorate sandwiches [4,5]. The propellant was made from sheets of AP-binder-AP. The AP formed by dry pressing ultra pure AP powder. Four different binders were used. Observations for the combustion were made by high-speed photography and microscopic examination of quenched samples. In addition, Lee et al. [4] illustrated the effect of inclusion of particulate AP in the binder on the combustion surface and the flame structure. The effect of three types of fuel binder and oxidizer particle diameter on the decomposition and combustion behavior of ammonium perchlorate is studies by Al-Harthi and Williams [6]. Few decades ago, several theoretical studies on the combustion field of the burning of the heterogeneous propellant have been conducted. These researches are divided into two main categories. The first one is concentrated on the gas phase modeling without consideration for the condensed phase process, for example [7– 11] The second one is studied the condensed phase reaction as
30
A.M. Hegab et al. / Computers & Fluids 89 (2014) 29–37
Nomenclature reaction rate constant overall regression speed specific heat diffusivity Damkohler numbers activation energies characteristic lengths Lewis number Mach number mass flux Peclet number Prandtl number pressure heats of reactions burning rate reaction rate Schmitt number
t T, X, Y
time temperature, mass fractions of AP and binder
Greek symbols b AP/binder stoichiometric ratio g surface function k heat conductivity m fractional binder thickness q density / equivalence ratio w function negative in binder, positive in AP Subscripts AP ammonium perclorate B binder g gas phase s solid phases
Periodic B. C.
Gas Phase Periodic B. C.
Vn
y Moving Interface at
AP
−ν
-L
t=0
x
rb
FUEL BINDER
the most important factor, for example [12,13]. Recently, several computational models [14–17] have been conducted to account the complex coupling between the solid-phase and gas-phase process. In particularly, the complexity that arises from the consideration of the unsteady non-planar regression surfaces. In the present paper a complete numerical strategy to examine what is perhaps, the simplest model is developed and account the following ingredients: the diffusion flame located at the interface between the decomposition products of the binder and the oxidizer (AP), the primary premixed AP flame, different properties (density, conductivity) of the AP and binder, temperature-dependent gas-phase transport properties, an unsteady non-planer regression surface; and a proper accounting of the fluid-mechanics in the gas-phase. These ingredients are applied to the problem of Periodic Sandwich Propellant (PSP).
-2L
AAP,B C Cp D Da1,2 E1,2 L Le M m Pe Pr P Qg1,2 rb Rr Sc
AP
ν
L
Fig. 2. Periodic sandwich propellant configuration.
2. The physical and mathematical models The physical model of the PSP model is shown in Fig. 2. This model consists of a sheet of fuel-binder of thickness ‘‘m’’, layered between two sheets of ammonium perchlorate (AP). The gas phase, consisting of a mixture of the decomposition products of the solid oxidizer and fuel is located above the combustion surface. Periodic boundary conditions are applied at x = ±L. The AP-binder-AP sandwich geometry has been chosen as a useful framework in which to gain fundamental insights into propellant combustion (e.g. [16–18]) and a notable experimental program has been pursued for some years by Price and his colleagues [5]. In addition, combustion behavior of the simpler sandwiches is much easier to observe and is enabled us to describe the combustion fields of the complex random packing propellant.
Diffusion flame
Leading edge flame (LEF) Stoichiometric surface
2.1. Constant and variable density models It is useful to summarize the formulation of the constant density model before addressing the complete problem, as this enable us to introduce most of the model ingredients together with various convenient scaling in the context of a model set of equations. The specific details of the constant density model for our problem are as follows: the density is set equal to constant (so that the equation of state, Charles law, is jettisoned); and a uniform velocity field u = 0 and v(y) = constant is adopted, which satisfies both the continuity and momentum equations. The one-step kinetics that include the primary flame is examined separately in order to achieve a good understanding of the unsteady burning of periodic sandwich propellant with complete coupling between the solid and gas phases. Thus R1
APðXÞ ! decomposition products ðZÞ R2
ð1Þ
Z þ binderðYÞ ! final products
Oxidizer-fuel mixing fan
R1 and R2 are assumed to have the forms;
AP flame Oxidizer (AP)
B
Oxidizer (AP)
Fig. 1. Flame structure for an AP-Binder-AP sandwich [4].
R1 ¼ B1 PX expðE1 =Ru TÞ and R2 ¼ B2 P ng YZ expðE2 =Ru TÞ; where B1 and B2 are the exponential prefatory, E1 and E2 are the activation energy in the gas phase, P is the pressure with exponent
31
A.M. Hegab et al. / Computers & Fluids 89 (2014) 29–37
ng, Ru is the universal gas constant, and (T, X, Y, and Z) are the temperature, oxidizer, fuel, and decomposition products respectively. The corresponding gas phase equations are;
equal to that in the gas phase for simplicity. The possibility of differing densities and thermal properties in the solid phase is allowed and setting by;
D/ ~ kg ~ qg ¼r r/ þ W Dt cp
qs ¼
ð2Þ
R1
R2
R1 bR2
T
When full fluid-mechanics coupling (variable density model) is accounted for, the system of Eq. (2) is replaced by;
@Q @F @G þ þ ¼H @t @x @y where
2
qg
3
ð3Þ
3 qg u 2 7 6 qg u þ P sxx 7 6 7 6 q u v s xy g 7 6 6 F ¼ 6 ðqg e þ PÞu ðusxx þ v sxy qx Þ 7 7; 7 6 qg uX qg Dg X x 7 6 5 4 qg uY qg Dg Y x qg uZ qg Dg Z x 3 2
qg v
7 6 qg uv sxy 7 6 2 7 6 q v þ P s yy g 7 6 6 ¼ 6 ðqg e þ PÞv ðusxy þ v syy qy Þ 7 7; 7 6 qg v X qg Dg X y 7 6 5 4 qg v Y qg Dg Y y qg v Z qg Dg Z y and
e ¼ qP þ 12 ðu2 þ v 2 Þ; g
ks ¼
kAP
wP0
kB
w<0
ð7Þ
;
gt þ gx
dx dy þ gy ¼ 0: dt dt
ð8Þ
and the final equation that controls the moving of the gas/solid interface as in Fig. 3 is derived by Hegab et al. [16,17] and may be written as follows;
~ gj ¼ 0; gt rb jr 2
6 qg u 7 7 6 6 qg v 7 7 6 6 Q ¼ 6 qg e 7 7; 6 qg X 7 7 6 4 qg Y 5 qg Z 2
qAP qB
The function w(x, y) is a level set function which demarks the regions of AP from binder (B) within the solid, so that a point (x, y) lies in the AP if w(x, y) P 0, and in the binder if w(x, y) < 0. Suppose the solid/gas interface defined by g(x(t), y(t), t) = 0. Then;
where
/ ¼ ½ T X Y Z T ; W ¼ ðQ g1 R1 þ Q g2 R2 Þ=cp
sxx ¼ l
where rb is defined as the speed of the front which moves in the directions of the solid. In general rb is a function of x and t and is given by the following simple pyrolysis law;
(
G rb ¼
3 0 7 6 0 7 6 7 6 0 7 6 6 H ¼ 6 þQ g1 R1 þ Q g2 R2 7 7 7 6 R1 7 6 5 4 R2 þR1 bR2
4
u 23 v y ; 3 x
ð9Þ
4
syy ¼ l 3 v y 23 ux ;
sxy = l(vx + uy); qx = kgTx; qy = kgTy Lewis number is taken to be unity, then;
qg Dg ¼ kg =cp
ð10Þ
Note that pressure dependence has been added to the pyrolysis law for generality. In the study, the propellant surface is not flat and its shape changes with time. Therefore, the following mapping function is used;
g ¼ y f ðx; tÞ
ð11Þ
and the front of Eq. (9) reduces to the simple Hamilton–Jacobi equation;
qffiffiffiffiffiffiffiffiffiffiffiffiffi ft þ r b ðx; tÞ 1 þ fx2 ¼ 0;
ð12Þ
ð5Þ
2.3. Boundary/jump conditions
2.2. Solid-phase and solid/gas interface equations
The appropriate jump conditions across the gas/solid interface are;
½qð~ v ~n þ rb Þ ¼ 0;
ð13Þ
½T ¼ 0; ½k~ n rT ¼ Q s m;
In the solid-phase, the following heat equation is used;
m½vi ¼ ½qD~ n rvi ; ð6Þ
Here, qs is the density of the solid, T the temperature, and ks is the solid thermal conductivity. The specific heat cp is assumed to be
ð14Þ i ¼ 1; 2; 3;
ð15Þ ð16Þ
where [/] = /g /s denotes the jump in the quantity / across the interface, vI refers to X, Y, and Z, m is the mass flux. ~ n is the unit normal pointing in the direction of the gas; ~ n ¼ rg=jrgj, Fig. 3. Qs is the solid phase heat release term defined by
Qs ¼
Fig. 3. Sketch showing the coordinate system for the moving surface.
w<0
Further information about the non-planar moving of the gas/solid interfaces using the level set strategy is mentioned in details in [16,17].
where ax and ay are the stoichiometric coefficients. Here there are six unknowns in the gas-phase, (u, v, T, X, Y) and one unknown in the solid-phase (Ts)
k qs T t ¼ s r2 T cp
r B ¼ AB ðP=P0 ÞnB expfEB =Ru T B;s g
ð4Þ
With the aid of the equation of state;
P ¼ qRT
r AP ¼ AAP ðP=P0 ÞnAP expfEAP =Ru T AP;s g w P 0
Q AP
wP0
Qb
w60
ð17Þ
For an exothermic surface reaction, Q s > 0, and for an endothermic reaction, Q s < 0. Typically the AP is considered an exothermic reaction, while the binder an endothermic one. Studies by Hegab et al. [16,17] proved that the length and time scales for the front and the solid are the same order of magnitudes. On the other side, the ratio of the gas to solid or the ratio of the gas to front is of the order of 103. Thus for the present purpose, the quasi-steady approximation for the gas phase is employed. Note
32
A.M. Hegab et al. / Computers & Fluids 89 (2014) 29–37
that disturbances with time scales of order 103 s would effect the solid phase, but not the gas phase; changes on time scales of order 105 s are needed to generate an unsteady gas phase and changes of this nature have been discussed in [11]. 2.4. Nondimensionalization
2
t s ¼ L=r b;ref ;
t g ¼ L=V g ;
e ¼ tg =ts 1 V g ¼ qAP rb;ref =qg ; T ref ¼ 2700 K; Q ref ¼ cp 2700 j=g
Q s;AP =ðcp T ref Þ w P 0 Q s;B =ðcp T ref Þ
kratio ¼
qratio ¼
w < 0;
1
wP0
kB =kAP
w < 0;
1
wP0
qB =qAP w < 0:
3 0 7 6 aux 7 6 7 6 4=3 a v x 7 6 2 7 V 2 ðQ ; Q y Þ ð c 1ÞM a fuu þ 4=3 vv g þ b T W 2 ðQ ; Q y Þ ¼ 6 y y y7 6 7 6 l =P X e y 7 6 5 4 l=Pe Y y l=Pe Z y 3 2 0 7 6 2=3av y 7 6 7 6 u y 7 6 2 6 ¼ 6 ðc 1ÞM af2=3uv y þ v uy g 7 7W 1 ðQ ; Q x Þ 7 6 0 7 6 5 4 0
ð18Þ
ð19Þ
ð20Þ
ð21Þ
7 6 av x 7 6 7 6 2=3 a u x 7 6 2 7 ð c 1ÞM a fu v 2=3 v u g ¼6 x x 7 6 7 6 0 7 6 5 4 0 0
ft þ r b
2
R2 ¼ Da2 P YZ expðhg2 =TÞ; but temperature-dependent transport is accounted for, viz; kg ¼ kg;ref ~ kðT; T ref Þ where kg,ref is a reference heat conduction coefficient. The value of k at the reference temperature Tref, specially kg = 1.08104 T + 0.0133 is choosing with dimensions W/m K when T is assigned in degrees Kelvin, so that
u ¼ 0;
1:08 104 T ref þ 0:0133
v ¼ qratio ft ;
R ¼ Da PXY expðhg =TÞ
ð22Þ ð23Þ
The non-dimensional equations of motion can be re-written in the following form:
@Q @FðQ Þ @GðQ Þ @V 1 ðQ; Q x Þ @V 2 ðQ ; Q y Þ þ þ ¼ þ @t @x @y @x @x @W 1 ðQ ; Q x Þ @W 2 ðQ ; Q y Þ þ þ þH @y @y
where
ð25Þ
qffiffiffiffiffiffiffiffiffiffiffiffiffi 1 þ fx2 ¼ 0;
ð26Þ
Tðx; 0þ ; tÞ ¼ Tðx; 0 ; tÞ
ð27Þ
1
2 f T þ 1 þ f T x x g C
x ^kB qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi @ A
1 þ fx2
0
0þ
g > 0 for Variable Density Model
e
where a ¼ PSce l; b ¼ PSc k, and h ¼ e þ c1 c p=q e Pr g < 0 (Solid Phase)
g = 0 (Moving Interface)
R1 ¼ Da1 PX expðhg1 =TÞ
~kðT; T ref Þ ¼
3
0
qratio ðT t ft T g Þ ¼ ðkratio =Pec ÞDT
R1 and R2 have the following forms;
1:08 104 T ref T þ 0:0133
0
2
In non-dimensional form the equations and boundary/connections form conditions are: g > 0 for Constant Density Model
~ ð Kr ~ TÞ þ Q g1 R1 þ Q g2 R2 eT t þ uT x þ v T g ¼ ð1=Peg Þr ~ ~ XÞ R1 eX t þ uX x þ v X g ¼ ð1=Peg Þr ðKr ~ ðKr ~ YÞ R2 eY t þ uY x þ v Y g ¼ ð1=Peg Þr ~ ðKr ~ ZÞ þ R1 bR2 eZ t þ uZ x þ v Z g ¼ ð1=Peg Þr
3
2
Pressure po (bar), surface speed rb,ref (cm/s), and mass flux mref = qAPrb,ref. Length L (half of the computational domain, which is the sum of the binder and the AP thickness). Time t = L/rb,ref. Then the following non-dimensional parameters are defined: Peclet numbers Peg = qgVg L cp/kg,ref, and Pec = qAP rb,ref L cp/kAP, Activation energy h = E/(RuTref)
bs ¼ Q
2
3 0 7 6 4=3aux 7 6 7 6 a v x 7 6 2 6 ¼ 6 ðc 1ÞM af4=3uux þ vv x g þ b T x 7 7 7 6 l=Pe X x 7 6 5 4 l=Pe Y x l=Pe Z x
T ¼ T=T ref ; X ¼ X=X s ; Y ¼ Y=Y s ; Z ¼ Z=Z s ; P ¼ P=Po ; q ¼ q=qo ; ðu ; v Þ ¼ ðu; v Þ=V g ; f ¼ f =L; rb ¼ rb =rb;ref ; t ¼ t=t s ;
3
2
2
The following reference values are taken to nondimensionalize the equations;
ðx ; g Þ ¼ ðx; gÞ=L;
3
q qu qv 6 qu 7 6 qu2 þ p 6 quv 7 ~7 7 7 7 6 6 6 2 6 qv 7 6 quv 7 6 qv þ p ~7 7 7 7 6 6 6 7 7 7 6 6 Q ¼6 6 qe 7 FðQ Þ ¼ 6 quh 7 GðQÞ ¼ 6 qv h 7 V 1 ðQ ; Q x Þ 6 qX 7 6 quX 7 6 qv X 7 7 7 7 6 6 6 4 qY 5 4 quY 5 4 qv Y 5 qZ quZ qv Z
ð24Þ
1
2 f T þ 1 þ f T x x g B C
x ^ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi kkratio @ A
1 þ fx2
0
b s mP e ¼ Q g
ð28Þ
8 91
>
> k^ <fx X x þ 1 þ fx2 X g =C
B qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi @mX A
> P eg > : ;
1 þ fx2 0
8 91
> > < 2 ^ k fx Y x þ 1 þ fx Y g =C
B q ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi mY @ A
> P eg > : ;
1 þ fx2
¼
m wP0 0
w<0
0
wP0
ð29Þ
0þ
0
8 91
> > < 2 ^ k fx Z x þ 1 þ fx Z g =C
B qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ffi
@mZ A > > P eg : ;
1 þ fx2
0
¼
m w<0
0þ
0
¼ 0þ
0
wP0
m w<0
ð30Þ
A.M. Hegab et al. / Computers & Fluids 89 (2014) 29–37
g ! 1 T ¼ T o g ! þ1
ð31Þ
@ ðÞ ¼ 0 @g
ð33Þ
3.1. Modified conservation equations To avoid pressure singularity at low Mach numbers, the numerical strategy outlines in [20,21] is used. The pressure is rescaled in the momentum equation since it is the pressure gradient, not the actual pressure that is involved in the momentum balance. The rescaled pressure is applied to retain accuracy in calculating the momentum conservation. As a result, the pressure is divided into constant and fluctuating parts as follow;
b y; tÞ Pðx; y; tÞ ¼ 1 þ cM 2 Pðx;
ð34Þ
and substitute into the equations of motion given above. Where M is the Mach number (M = Vo/Co). The equation of state now becomes;
b y; tÞ 1 þ cM 2 Pðx; T
ð35Þ
and is used to update the density. Since the quasi-steady state solution in the gas phase (e 1) is required, the physical time term in (24) is eliminated. In addition an accelerate convergence technique is employed [20]. The technique begins by adding a pseudo-time derivative to the conservation Eq. (24). So the equations to be solved for the gas phase become;
e @FðQÞ @GðQ Þ @V 1 ðQ ; Q Þ @V 2 ðQ ; Q y Þ @Q x þ þ ¼ þ @x @x @y @x @s @W 1 ðQ ; Q x Þ @W 2 ðQ ; Q y Þ þ þ þH @y @y
ð36Þ
where s represents the pseudo-time domain. Because the pseudotime derivative vanishes at convergence to the steady state solution in the gas phase, a certain amount of liberty is given in choosing the e . A new scaled pressure term P=b e is added as a pseuvariables in Q do-time derivative term to the continuity equation. The remaining e are then fixed by rewriting the momentum, energy, variables in Q and species equations from their neoconservative form by means of the modified equation of mass conservation. As a result, the pseue and its associated preconditioning matrix do-time variable vector Q C take the form;
2 3 2 b 1=b0 0 0 0 0 0 P 6 7 6 0 q 0 0 0 0 6u7 6 u=b 6 7 6 6 7 6 v =b0 0 q 0 0 0 6v 7 6 6 7 6 e ¼ 6 h 7C ¼ 6 h=b0 1 0 0 q 0 0 Q 6 7 6 6 7 6 6X 7 6 X=b0 0 0 0 q 0 6 7 6 6 7 6 0 0 0 0 0 q 4Y 5 4 Y=b Z
Z=b0
0
0
0
0
0
0
3
7 07 7 07 7 7 07 7 7 07 7 7 05
ð37Þ
q
b Its selecFactor b’ is used to properly scale the time derivative of P. tion is proportional to the dynamics pressure of the flow field. 3.2. ADI algorithm development The first step for advancing the solution of (36) is using the Delta form scheme [19] as follows;
ð38Þ
and
DQ n ¼ Q nþ1 Q n
3. Numerical methods
C
DQ n ¼
ð32Þ
jxj ¼ 1 Periodic boundary conditions
q¼
hDs @ Ds @ n ð DQ n Þ þ ðQ n Þ þ DQ n1 1þn 1 þ n @s 1 þ n @s þ O½ðh 1=2 nÞDs2 þ Ds3
33
ð39Þ
The time-difference formulas (38) and (39), with the appropriate choice of the parameters n and h reproduce many familiar two and three level, implicit schemes. The three level, second order accuracy implicit schemes (n = 1/2 and h = 1) with special treatment for the cross derivative terms at the level (n 1) is applied. By inserting the temporal derivative of Eq. (36) in (38) and by straightforward derivation, the resulting approximated form can be splitted into two-tridiagonal systems;
" I þ hDsC1 " I þ hDsC1
!# @ @2 ðA P þ Rx Þn 2 ðRn Þ DQ ¼ C1 RHS @x @x !# 2 @ e þ Sy Þn @ ðSn Þ DQ ¼ DQ ðB Q @y @y2
ð40Þ ð41Þ
e n ; Sn , and Rn are the jacobian matrices for where An ; Pn ; Rn ; Rnx Bn ; Q y the implicit x- and y-sweep (40) and (41) respectively. I am a unit matrix (6 6) and C1 am the inverse of the preconditioning matrix. The three-points second order accuracy central difference approximation is used for the spatial differences in the LHS of x-sweep form (40). This approximation along with the applied periodic boundary conditions at jxj = 1 produce a system of Periodic Block-Tridiagonal Equations (PBTE). After the computation of DQ⁄ at the interior points by solving this system of PBTE with each block having dimensions 6 6 components, the code is ready for the implicit y-sweep form (41). Here again, the three-point second order accuracy central difference approximation are used for the spatial differences in the LHS of (41). This approximation along with the rigid wall boundary conditions at y = 0 produce a system of NonPeriodic Block-Tridiagonal Equations (NPBTE). The final delta form DQ can be computed by solving the NPBTE system. Then the solution at new time step (n + 1) can be determined from Eq. (40). Note that, the cross-derivatives terms at (n 1) are treated explicitly to avoid the implicit coupling of adjacent boundary points. The surface Eq. (12) is solved in order to follow the non-planar regression surface by the first order temporal scheme [22]. Besides the non-flat regression surface mapping as in (10), another transformation is applied for clustering grid points adjacent to the wall, where most of the flow parameters changes rapidly. The solution of the final mapped equations is advanced in the solid phase using physical time (t). Simultaneously the solution in the gas phase using pseudo-time (s) to the local steady state at the first physical time step (t) is advanced. The boundary/jump conditions are continually updated as in Eqs. (25)–(33). Then the Hamilton–Jacobi Eq. (12) is advanced at the physical time by a third order ENO and a fifth-order WENO (weighted essentially non-oscillatory) solver [23,24]. All numerical calculations were performed on a 140 70 grid, uniform in the x-direction and stretched in the y-direction. At each physical time, the solution in the gas phase is advanced until the relative difference between each two different pseudo-time values is less than some prescribed tolerance, taken here to be 106. Convergence tests where carried out and it was determined that any further refinement resulted in less than 1% relative error.
34
A.M. Hegab et al. / Computers & Fluids 89 (2014) 29–37
Fig. 4. Reaction rate contours at different times in the gas-phase for hg1 = 3.5, hg2 = 6.0, b = 7.51, Peg = 6, P1s = 2.9, Da1 = 5104, Da1 = 1.7105, m = 0.22, Qg1 = 0.833, Qg2 = 3.86, Qs,AP = 0.4, and Qs,B = 0.06.
4. Results and discussion The understanding of the complex combustion structure of APBinder-AP sandwich, as a simple model to the heterogeneous solid rocket propellant, is studied in details by two different models. The first model is the constant density model and the second one is the variable density model or the Navier–Stokes model. Initially, the solution starts from a flat surface f(x, t = 0) = 0. Then the solution is advanced in the solid phase, gas phase with simultaneously moving interface. The first set of the results are for the constant density model. In this model, it has been taken l = K(T, Tref) which is one of the more realistic choices rather l = constant in earlier studies [7].
In the gas phase, the contours plots for the reaction rate contours (R) at different times are shown in Fig. 4 for hg1 = 3.5, hg2 = 6.0, b = 7.51, Peg = 6, P1s = 2.9, Da1 = 5104, Da1 = 1.7105, m = 0.22, Qg1 = 0.833, Qg2 = 3.86, Qs,AP = 0.4, and Qs,B = 0.06, as benchmark computation. The upper portion represents the gas phase and the lower one refers to the solid phase. The dark gray region in the latter represents the binder layer between the two AP sheets (light gray). The contour values are written in the upper part of each plot. The plot shows the location and shape of the generating flame. These reaction rate contours show a combined of two different flames. The first one is the diffusion flame that formed at the interface between the oxidizer and the fuel, while the second one is the premixed flame that formed above the oxidizer regions.
Burning surface
Stoichiometric Contour Fig. 5. Comparison of experimental image for the burning of AP-Binder-AP sandwich[26] {left} with the numerical model {right}.
A.M. Hegab et al. / Computers & Fluids 89 (2014) 29–37
Fig. 6. Variations in the global regression (c) rate with pressure (p).
The contours of the diffusion flame show that reaction rate is characterized by two strong mixing structure each centered at jxj 0.223. The successive curved shapes through the solid phase show that, the surface is initially flat and then as the solution is advanced, the combustion surface retreats in an unsteady fashion and finally retreats at a fixed speed with unchanged shape by t = 2.49. These contours reveal not only the significant effect of the surface profiles and the burning rate on the shape of the flame but also on the maximum reaction rate values which decrease as time increase. The differences in profiles with time advanced reflect the behavior of the burning rate at the propellant surface and reconstruct the corresponding flame structures and identify the parts of the flame structure that dominate the sandwich burning rate and the surface heat flux. In addition there is a sharp discontinuity in slope at the interface between the binder and AP regions. This interesting phenomenon is predicted in [25]. Fig. 5 shows a comparison between the experimental image (the left) [26] and the current computation model (the middle for the reaction rate and right for the fuel vapor contours) for the structure of the flame shape and the burning surface. The stoichiometric level surface shows the location at which the fuel and oxidizer meet together in shoichiometric proportions. The inner
35
region for the stoichiometry contour is the fuel-rich region (see the fuel vapor image on the left side), while the area out of the stoichiometry envelop represent the fuel-lean composition. The two stoichiometric surface in both the computational and experimental images represent the two strong mixing structure (at the AP-Binder interface) each centered at jxj 0.223 and forming the flame envelop in addition to the AP-premixed flame that control the shape of the combustion surface above the AP region. This comparison show how the computational steadily surface regressing profile (at t > 5) is qualitatively consistent with the experimental emission-transmission composite image by Brewster et al. [26], [26] and also with the experimental photography for the spontaneous quench samples of AP/BPAN/AP sandwich that reported by Lee et al. [4]. Both theoretical and experimental results verified that conditions that gave the highest burning rates resulted in narrow smooth bands of binder in the middle with little ‘‘protrusion’’ of AP at the interface. The variation of the burning rate c with the pressure is presented in Fig. 6 for the benchmark data as in Fig. 4. The response of the burning rate to the exposed pressure is one of the important and essential characteristics of the propellant. This figure shows that the burning rate is directly proportional to the pressure with exponential component about (n = 0.6). This exponential value has great interest to the rocket designer, in particular for the modeling of the composite propellant rather the sandwich one. The effect of the exposed pressure on the burning surface profiles is presented in Fig. 7 at three pressure values P = 1, 10, 15 atm in the gas phase for the same parameters in Fig. 4. A comparison between the experimental image and the computational model for the response of the combustion surface profiles to the pressure change (experimental at 7 atm in the left and the computational at 1, 10, and 15 atm in the right) is presented in Fig. 7. Both images illustrate that increasing the pressure causes the binder to stick out and the surface consist of a ‘‘trough’’ centered on the binder lamina. In addition a comparison of these clear pictures with the experimental results by Price and et al. [5] is possible but caution is appropriate as the sandwiches that Price examined are isolated, not periodic. Furthermore, the effect of the pressure on the surface heat flux reveal that the heat flow vector, near the gas–solid interface and AP-binder solid interface, is from the binder toward the AP in the solid phase. Consequently there is a weak (hot) portion of the
Fig. 7. Effect of pressure on the burning surface structure. Left figure is experimental [5] and the right figures are the computational model.
36
A.M. Hegab et al. / Computers & Fluids 89 (2014) 29–37
Fig. 8. Total reaction rate contours (a) constant density model and (b) variable density model.
Fig. 9. The axial velocity component for the variable density model (2) at different times for h = 7, b = 7, Pe = 6, Dag = 5106, P = 10, M = 0.02.
binder slightly up from the corner interface and into the binder. This might account for the appearance of the V-shape (notches) in the binder. These phenomena have been seen in the experimental quenched samples by Handley and Strahle [25]. The last set of the results are for the variable density model (2). In this section the results obtained from the constant density model (1) with those obtained using the Navier–Stokes equations, model (2) are compared, and thus to validate the simpler strategy. Fig. 8 shows a comparison between the two models for the reaction rate contours for the parameters in Fig. 4. The reaction rate values
varied from 12 to 2 with incremental 2 for the AP-premixed flame, while varied from 1 to 0.2 with incremental 0.2 for the AP-Binder diffusion flame that centered near the interface between the fuel and oxidizer biased little bite toward the oxidizer surface. It is noted that the flame is setting at the same transverse location with slight differences in the maximum reaction rate contours. RRmax is 11.88678 for model (2), compared to 11.66972 for model (1) and in the adiabatic flame temperature, where T 0max is 2448.09 K for model (2) compared to 2446.57 for model (1). These interesting values reveal that there is only a small difference in the far-field for the
A.M. Hegab et al. / Computers & Fluids 89 (2014) 29–37
two models. These very slight differences are related to the slow flow accelerated since the axial velocity component u is not zero when the Navier–Stokes equations are used and the transverse velocity component becomes more or less in the neighborhood of the flame for model (2). Fig. 9 shows the axial velocity contours at different times for the Navier–Stokes model. It is noted that an axial velocity is generated near the combustion surface due to the surface morphology, where the velocity cells are positive on the left hand side, since the flow goes down hill and the right hand side cells are negative, since the flow goes up hill. Moreover, the absolute value of the axial velocity increases as time increases, since the curvature in the combustion surface profiles become more deeper than at earlier time. The small differences between the two models reveal that, a useful exploration calculations can be carried out using the constant density model, since the generated axial velocity component is very small, 0.0 6 juj P 0.04 and also the variation in v is consistent with the mass conservation. In general, the Navier–Stokes model may have a significant effect when the current solution to the modeling of randomly packed heterogeneous propellant, particularly in 3-D solution is advanced. 5. Conclusion Few decades ago the details of the heterogeneous propellant combustion has been discussed. The theoretical framework has consists of a 2D physical pictures, but most modeling efforts have been 1D mathematical model as in [1–4]. This has obviously been of value and important insights have been achieved, but a multidimensional numerical framework is required. Here, for the first time, the 2D calculations to the combustion of heterogeneous solid propellant, accounting for the gas phase physics, the solid phase physics and an unsteady non-planar description of the regressing propellant surface is developed. There are a number of issues that have discussed. The speed within which the combustion surface recedes depends on the exposed pressure in the gas phase, the effect of several parameters on the combustion and shape of the flame. In addition a variety of steady-state surface shapes are achieved. At higher pressure values, the binder is tending to stick out and the surface consists of a ‘‘trough’’ centered on the binder lamina. These trends were also recognized in an experimental observation [5]. A comparison between the computational steadily surface regressing profile with the experimental emission-transmission composite image by Brewster et al. [26], [24] and also with the experimental photography for the spontaneous quench samples of AP/BPAN/AP sandwich that reported by Lee et al. [4] is qualitatively consistent. In this study, the Navier–Stokes equations were used rather the constant density model in earlier work [7,10,16,18]. A comparison between the constant density model and the Navier–Stokes solutions reveals very small differences. More intensive work for the solution of the full Navier–Stokes as in [27–32] is required to perform more complex computation for the reactive whole scale of the solid rocket motor chamber. Acknowledgments This work was funded by the Deanship of Scientific Research (DSR), King Abdulaziz University, Jeddah, under Grant No. (829003-D1433). The authors, therefore, acknowledge with thanks DSR technical and financial support. The first author would like to acknowledge Prof. John Buckmaster and Prof. Tomas Jackson at University of Illinois at Urbana–Champaign (UIUC), USA for their support for developing the variable density model.
37
References [1] Beckstead MW, Derr RL, Price CF. A model of composite solid-propellant combustion based on multiple flames. AIAA J 1970;8(12):2200–7. [2] Rasmussen Jr B, Frederick RA. A nonlinear heterogeneous model of composite solid propellant combustion. In: AIAA paper no. 99-2228. 35th AIAA/ASME/ SAE/ASEE joint propulsion conference; June 1999. [3] Murphy JJ, Krier H. Heterogeneous effects on dynamic burning in composite solid propellant. In: Proceeding of the combustion institute, vol. 28; 2000. [4] Lee S, Price E, Sigman R. Effect of multidimensional flamelets in composite propellant combustion. J Propul Power 1994;10(6):761–8. [5] Price E, Handley J, Panyam R, Sigman R, Ghosh A. Combustion of ammonium perchlorate-polymer sandwiches. AIAA J 1981;19(3):380–6. [6] Al-Harthi A, Williams A. Effect of fuel binder and oxidizer particle diameter on the combustion of ammonium perchlorate based propellants. Fuel 1998;77(13):1451–68. [7] Jackson TL, Buckmaster J, Hegab A. Periodic propellant flames and fluidmechanical effects. J Propul Power 2001;17(2):371–9. [8] Buckmaster J, Jackson TL, Yao J. An elementary diffusion of propellant flame geometry. Combust Flame 1999;117:541–52. [9] Jackson TL, Buckmaster J, Hoeflinger J. Three-dimensional flames supported by heterogeneous propellants. In: JANNAF paper, JANNAF CS/PSHS/APS joint meeting, Cocoa Beach, FL; October 1999. [10] Jackson TL, Buckmaster J. Nonpremixed periodic flames supported by heterogeneous propellant. J Propul Power 2000;16(3):498–504. [11] Cai W, Yang V. A model of AP/HTPB composite propellant combustion. In: AIAA paper 2000-0311, 38th aerospace science meeting; 2000. [12] Waesche RH, Wenograd J. Calculation of solid propellant burning rate from condensed-phase decomposition kinetics. AIAA paper 69-145, NY; January 1969. [13] Hermance CE. A model of composite propellant combustion including surface heterogeneity and heat generation. AIAA J 1960;4:1629–37. [14] Hegab AM. An overview of solid propellant combustion. In: Proceeding of the 12th aerospace science & aviation technology conference, (ASAT-12), 29–31 May, 2007, BAL-03, Military Technical College, Cairo, Egypt; 2007. p. 1–24. [15] Hegab AM, Balabel A. Effect of ammonium perchlorate grain sizes on the combustion of solid rocket propellant. In: Proceeding of the 12th aerospace science & aviation technology conference, (ASAT-12), 29–31 May, 2007, CHA1, Military Technical College, Cairo, Egypt; 2007. p. 1–20. [16] Hegab A, Jackson T, Buckmaster J, Stewart S. The burning of periodic sandwich propellants. In: 36th AIAA/ASME/SAE/ASEE joint propulsion conference, AIAA paper 2000-3459; 2000. [17] Hegab A, Jackson T, Buckmaster J, Stewart S. Nonsteady burning of periodic sandwich propellant with complete coupling between the solid and gas phases. Combust Flame 2001;125(1/2):1055–70. [18] Knott GM, Brewster MQ. Two-dimensional combustion modeling of heterogeneous propellants with finite Peclet number. Combust Flame 2000;121(1/2):91–106. [19] Tseng IS, Yang V. Combustion of a double-base homogeneous propellant in a rocket motor. Combust Flame 1994;96:325–42. [20] Roh TS, Yang V. Transient combustion responses of solid propellants to acoustics disturbances in rocket motors. AIAA paper 95-0602; 1995. [21] Beam RM, Warming RF. An implicit factored scheme for the compressible Navier–Stokes equations. AIAA J 1978;16(4). [22] King MK. Examination of chemical approaches to stabilizing composite propellant combustion. J Propuls Power 1996;12:554–63. [23] Harten A. Uniformly high order accurate essentially non-oscillatory schemes, III. J Comput Phys 1987;71:231–303. [24] Jiang G, Shu C. Efficient implementation of weighted ENO scheme. J Comput Phys 1996;126:202–28. [25] Handley JC, Strahle WC. The behaviour of several catalysts in the combustion of solid propellant sandwiches. AIAA paper 74-122, Washington (DC); 1974. [26] Brewster MQ, Knott GM, Chorpening BT. Combustion of AP/HTPB laminate propellants. In: AIAA paper 2001-4501, 37th AIAA/ASME/SAE/ASEE joint propulsion conference, Salt Lake City, UT, USA; 2001. [27] Hegab AM, Sait HH, Hussien A, Balabel A, Almutawa M, Djouider F. Modeling of microscale solid rocket composite propellant. Life Sci J 2013;10(4). [28] Hegab AM, Hasanien SA, Mostafa HE, Maraden AM. Experimental and numerical investigations for the combustion of selected composite solid propellant. Eng Res J 2012;35(3). ISSN:1110-1180. [29] Maraden A, Hossam E, Hasanien S, Hegab A. Effect of packing parameters on the combustion behaviour of ap based composite solid propellant. In: 6th international conference on chemical & environmental engineering: energetic material I, Military Technical College; 29–31 May, 2012. [30] El-Askary WA, Balabel A, El-behery SM, Hegab AM. Computational of a compressible turbulent flow in a rocket-motor chamber configuration with symmetric and asymmetric injection. J Comput Model Eng Sci 2011;82(1):29–54. [31] El-Askary WA, Wilson SA, Hegab AM. Simulation of reactive fluid flow in a solid rocket motor combustion-chamber with/without nozzle. J Comput Model Eng Sci 2011;76(4):235–66. [32] Balabel A, Hegab A, Nasr M, El-Behery S. Assessment of turbulence models for gas flow in convergence–divergence rocket nozzle. Appl Math Model 2011;35:3408–22.