Spatial linear analysis of the flow in a solid rocket motor with burning walls

Combustion and Flame 156 (2009) 865–888 Contents lists available at ScienceDirect Combustion and Flame www.elsevier.com/locate/combustflame Spatial...

Download PDF

806KB Sizes 1 Downloads 61 Views

Report

PDF Reader
Full Text

Combustion and Flame 156 (2009) 865–888

Contents lists available at ScienceDirect

Combustion and Flame www.elsevier.com/locate/combustflame

Spatial linear analysis of the ﬂow in a solid rocket motor with burning walls L. Massa ∗ Department of Mechanical & Aerospace Engineering, The University of Texas at Arlington, Arlington, TX 76019-0018, USA

a r t i c l e

i n f o

Article history: Received 26 June 2008 Received in revised form 30 September 2008 Accepted 11 November 2008 Available online 8 December 2008 Keywords: Solid propellant Erosive burning Rocket ﬂow

a b s t r a c t A theoretical investigation of the coupling between the core ﬂow process and combustion sublayer process in a rocket chamber ﬂow is presented. The focus of the investigation is on penetration of chamber disturbances into the sublayer and on the burn rate characteristic responses. Characteristic solutions to the local non-parallel instability problem supported by gasifying propellant in a double slab geometry are sought. In the fully decoupled limit the only characteristic solution response is a pressure response, thus comparisons between coupled and decoupled pressure responses are presented. The propellant is typiﬁed by a set of burning characteristics, and their effect on the coupling is explored. The goal of the analysis is to clarify the extent strand burning characteristics (i.e., in zero cross ﬂow conditions) can affect burn rate perturbations leading to the phenomenon of erosive burning. The ﬁndings of this research agree with experimental studies as it regards the effect of burning parameters on erosive behavior. Change in erosive sensitivity with propellant characteristics is thus correlated with a change in the core–sublayer interaction. © 2008 The Combustion Institute. Published by Elsevier Inc. All rights reserved.

1. Introduction This paper presents an investigation of the interaction between propellant burning and unstable core ﬂow in a solid propellant rocket chamber. Such interaction is most often considered responsible for the phenomenon of erosive burning (e.g., Razdan and Kuo [1]). Erosive burning denotes the burn rate augmentation of solid propellant in cross ﬂow conditions. The burn augmentation is an important phenomenon in rocket propulsion analysis because of the associated overpressure in the early stages after ﬁring. Due to the large difference in length scales between combustion sublayer and chamber core ﬂow, as discussed in more detail in Section 5 of this paper, the development of the instability and the ensuing turbulent ﬁeld are fairly independent of the burning characteristics. Therefore, the experimentally established link between erosive burning augmentation and propellant characteristics (e.g., King [2]) cannot be resolved without hypothesizing an effect of the instability on the micron-size combustion proﬁle; that is, hydro-dynamic considerations alone are not suﬃcient. Mukunda and Paul [3] argue that the effect of the burning characteristics can be taken into account by simply scaling the port mass ﬂow rate by the strand value. Still, the proposed effect of the Reynolds number on erosive burning cannot be explained based on the relationship between the same parameter and the instability growth rate. The most successful models that take into consideration the variation of the combustion proﬁle under cross ﬂow conditions are

*

Fax: +1 817 272 5010. E-mail address: [email protected].

those that postulate a turbulence increased transport in the combustion sublayer. Eminent models of this type are those of Most et al. [4] and King [2]. Although successful in matching experimental data, such models are based on a very simpliﬁed approach where the Prandtl mixing length is determined based upon ﬂat plate experiments, and a constant Prandtl number approximation is used to determine the turbulent thermal conductivity. Other computational analyses of erosive burning include those of Beddini [5] and Razdan and Kuo [6]. Although these analyses rely on different assumptions, they show a good agreement with experimental data. Beddini uses a simpliﬁed chemistry model where reactant mass fraction ﬂuctuations are ignored in the source term. Razdan and Kuo kinetic and diffusion models are also quite simpliﬁed, with assumptions that include the use of a constant Prandtl and Schmidt numbers, which do not allow for accurate resolution of the surface coupling effects. None of these studies considers in depth the interaction between the core and the sublayer, and the conditions that favor ﬂuctuation penetration in the ﬂame zone. The present paper attempts to provide a more rigorous analysis of the effect of burning parameters on the core ﬂow–sublayer coupling by seeking characteristic solutions to the linearized and localized solid–gas fully coupled problem. Although the linearized analysis does not provide a clear quantitative answer on the burning rate augmentation, burn rate eigenfunctions (ﬂuctuation intensity) quantify the importance of the coupling. The effects of burning parameters and motor geometry on the characteristic solutions are matched against experimental observations to reveal their importance in allowing penetration of ﬂuctuating eddys in the ﬂame structure. The approach provides insights on the role of parame-

0010-2180/$ – see front matter © 2008 The Combustion Institute. Published by Elsevier Inc. All rights reserved. doi:10.1016/j.combustﬂame.2008.11.007

866

L. Massa / Combustion and Flame 156 (2009) 865–888

2.2. Gas phase equations

Fig. 1. Double slab rocket motor geometry.

ters such as temperature sensitivity and pressure exponent on the penetration of core ﬂow disturbances in the sublayer, which were not included in previous studies. This paper also provides a basis for future non-linear analysis of erosive burning problem. Characteristic solutions are derived from a linear spatial normal mode analysis of the coupled solid–gas combustion problem. The analysis uses two principal assumptions. First, a local non-parallel analysis is used to represent the spatial development of the instability as detailed in Section 4.1. Second, ﬁnite rate chemistry is considered important only in the sublayer, while equilibrium conditions prevail in the (considerably hotter) core ﬂow. Consequently, the combustion model is composed of a decomposing reactant and an equilibrium product, the property of which are evaluated considering a large set of chemical species supported by the propellant atomic make-up, and minimizing the mixture Gibbs free energy, see Section 3. A set of burning characteristic is chosen to represent the combustion sublayer, as described in Section 3.1. The limiting case for which the difference in length scales is large enough that the core and sublayer can be considered separated is referred to as decoupled problem and described in Section 5. Homogeneous propellants of different chemical composition are analyzed as described in Section 6. The inﬂuence of the instability on the sublayer is studied by considering temperature and reactant mass fraction eigenfunctions, and comparing these to the decoupled problem analogs in Section 7. The effect of varying each burning characteristic on the coupling is presented in Section 8. 2. Governing equations 2.1. Geometry The double slab geometry motor is sketched in Fig. 1. This geometry has been widely used in computational and experimental rocket simulations, yields a set of governing equations simpler than the axisymmetric analog, and retains the fundamental attributes that concern the core–sublayer coupling. Mukunda and Paul [3] analyze various core ﬂow conﬁgurations and conclude that the erosive to non-erosive burn rate ratio is essentially a function of the Reynolds numbers and the ratio between port and propellant mass ﬂow rates. One of the ﬁrst erosive burning measurement on a double slab geometry was performed by Stokes et al. [7]. This geometry was most recently studied in an experimental investigation of erosive burning effects by Hasegawa et al. [8]. In the following the symbol · denotes a vector in the threedimensional Cartesian space, while capital letters denote general (directional) column vectors and matrices. The position vector is x . The internal energy is denoted with associated velocity vector u with U and is a non-linear function of temperature (T ) the pressure (p), and the reactant mass fraction (Y ) deﬁned by the minimization of the Helmholtz free energy as described in Section 3.2. The total energy is denoted as e 0 , the enthalpy as h (H is reserved for the enthalpy per unit mole), the density as ρ and the thermal conductivity as λ. The gas phase domain is (0, ∞) × (− D w , D w ) × (−∞, ∞), while the solid is (0, ∞) × (−∞, − D w ) × (−∞, ∞), and (0, ∞) × ( D w , ∞) × (−∞, ∞).

The governing equations for the gas phase are the reactive Navier–Stokes equations with a global gas phase reaction plus a surface zeroth order gasiﬁcation reaction. Consider the following vectors and matrices: the conservative variables vector Q = [ρ , ρ u 1 , ρ u 2 , ρ u 3 , ρ e 0 , ρ Y ] T ; the primitive variables vector P = [ p , u 1 , u 2 , u 3 , T , Y ] T ; the Jacobian matrix of the transformation Q ( P ), J = ∂∂ QP ; the advective ﬂux vectors

F 1 = Q 2 , Q 2 P 2 + P 1 , Q 3 P 2 , Q 4 P 2 , P 2 ( Q 5 + P 1 ), Q 6 P 2

T

,

F 2 (same as F 1 by interchanging u 1 with u 2 and rows 2 with 3), F 3 (same as F 1 by interchanging u 1 with u 3 and rows 2 with 4);

i the Jacobian of the transformations F i ( P , P ( Q )), A i = ∂∂FP ; the dif-

fusive V i and source term S vectors

∂ u1 ∂ u2 ∂ u1 ∂ u3 ∂ u1 , σ1 = μ 2 − 2/3∇ · u , + , + ∂ x1 ∂ x1 ∂ x2 ∂ x1 ∂ x3

(1a)

q = −(κ T , T ∇ T + κ T ,Y ∇ Y ),

(1b)

v = −(κY , T ∇ T + κY ,Y ∇ Y ),

(1c)

⎡

⎤

0

⎢ σ 1, 1 ⎢ ⎢ σ 1, 2 G =⎢ ⎢ σ 1, 3 ⎣ u · σ − q 1

1

1

⎥ ⎥ ⎥ ⎥, ⎥ ⎦

(1d)

˙ ]T , S = [0, 0, 0, 0, 0, −ω

(1e)

v1 V1 =

∂ G1 ∂ x1

and

the vector V 2 is obtained from V 1 by interchanging indexes 1 with 2 and rows 2 with 3, the vector V 3 is obtained from V 1 by interchanging indexes 1 with 3 and rows 2 with 4. In order to derive the linearized perturbation equations the vectors V i are here considered functions of the primitive variable vector and its derivative, V i = V i ( P , P xk , P xl ,xm ), for k, l, m = 1 . . . 3 and m l. The governing equation in the gas phase is

∂Fi ∂Q + − V i + S = 0. ∂t ∂ xi 3

3

i =1

i =1

(2)

2.3. Solid phase equations It is convenient to denote the solid temperature with a letter different from the gas phase one, T s . The three-dimensional heat conduction equation is the only equation considered in the solid phases x2 ∈ [−∞, − D w ] and x2 ∈ [ D w , ∞],

∂Ts ∂Ts + rb = K∇ 2 T s . ∂t ∂ x2

(3)

The burn rate rb is related to the surface temperature through a surface decomposition reaction, a zeroth-order Arrhenius equation,

rb = rb,0 exp −

Θ Ts

,

(4)

while K is the thermal diffusivity, which is taken to be a constant dependent on the chemical composition of the solid. 2.4. Connection conditions The connection conditions include non-slip velocity conditions plus ﬂux balances at the interfaces, x2 = − D w and x2 = D w ,

L. Massa / Combustion and Flame 156 (2009) 865–888

Table 1 Properties of solid ingredients. Material

λ (W/(m K))

ρ (kg/m3 )

AP HTPB

0.405 0.276

1950 920

u 1,3 = 0,

(5a)

T = T s,

(5b)

F 12

= ρs rb ,

(5c)

F 52 − G 25 = ρs rb h s,298

∂Ts , K + C s ( T s − 298 K) − λs ∂ x2

F 62 − G 26 = ρs rb Y s ,

(5d) (5e)

where C s is the solid phase speciﬁc heat capacity, taken constant and equal to 0.3 cal/g K; Y s denotes the percent of reactant in the solid propellant. This quantity is set to 1 and is clearly different from the mass fraction of reactant on the gas side of the interface, a result consistent with the ﬂux boundary conditions. The thermal conductivity and density of the solid phase λs and ρs are constant and related to the solid chemical composition through the homogenization formulae described by Chen et al. [9]. These formula are based on the oxidizer volumetric fraction, and on the properties of the pure ingredients discussed in Table 1. h298 K denotes the standard enthalpy of formation. The cold temperature propellant boundary condition completes the set of boundary conditions, T s = 0,

as x2 → −∞.

(6)

2.5. Non-dimensionalization Casalis et al. [10] report that for large values of the axial distance the local non-parallel eigenvalue problem tends to become independent of the Reynolds number. Moreover near the head end of a rocket motor the Mach number is very low due to the high gas temperature, and, as a result of the large difference in scale dimension between the core ﬂow and the combustion sublayer, the instability problem is essentially independent of propellant burning parameters when appropriately non-dimensionalized. Therefore the core ﬂow dimensions are used as references in the non-dimensionalization. The length scale that controls the inviscid, incompressible, non-burning ﬂow in a double-slab motor is a geometrical scale here taken to be half the port height, D w . The velocity scale, U w , is taken as the maximum normal velocity in the mean proﬁle. The use of the maximum velocity is justiﬁed by the normal acceleration in the combustion layer due to the strong temperature gradient, so that the velocity at the edge of the combustion layer is the maximum value. If the surface normal velocity were to be taken as scale, similarly to a non-burning analysis, nondimensional peak frequencies would be strongly dependent on the adiabatic ﬂame temperature. Results presented in Section 8.4 validate the assumption that the growth parameters scale with the core ﬂow dimensions and show that the combustion layer acceleration has little effect on the core ﬂuctuation dynamics. The temperature scale is taken as the adiabatic ﬂame temperature for strand (no-crossﬂow) conditions. The time scale is the convective value based on the previous indicated velocity and length scales. The pressure scale is twice the dynamic pressure based on the velocity scale and the density at the edge of the combustion layer. The edge of the combustion layer will be identiﬁed with the initialism ECL in the remainder of this paper. Besides the core ﬂow, two other phenomena control the propellant burning, the combustion sublayer and the solid-phase heat conduction. The scales associated with these two phenomena are

867

determined based upon the no cross ﬂow analysis. Strand burning propellant has no geometrical length scales, but a solid thermal , and a ﬂame standoff length, diffusion length, D t ≡ r ρλsc = K r b s p ,s

b

D f , representing a relationship between blowing velocity and chemical times. This research has focused on half port heights equal to 0.625, 1.25, 2.5, and 5 cm, and on several propellant speciﬁcations, each distinguished by a value for the solid thermal length. An average propellant speciﬁcation both in terms of thermal length and standoff distance is the 73% AP-HTPB mix. For this mix the thermal length at 30 atm is D t = 22 μm and the ﬂame stand-off measured as the maximum temperature derivative location for a strand is D f = 7.1 μm. The length scale analysis underscores a profound difference between core and sublayer dimensions, and suggests the two might considered independent. 2.6. Mean proﬁle Dunlap et al. [11] found the Culick–Taylor self-similar proﬁle to be an accurate approximation of a viscous laminar ﬂow in a rocket cold ﬂow model. This observation applies also to the core ﬂow of burning conﬁgurations because of the low Mach number, while diffusion contributions are important within the combustion wave. The application of a pseudo-time marching scheme to the solution of the steady, reactive Navier–Stokes equations is impractical because of the large difference in the problem length scales: the ﬁrst grid spacing should be approximately a tenth of a micron to resolve the combustion wave with a non-conservative 8th-order ﬁnite difference discretization. Therefore, the mean proﬁles for the spatial stability analysis are evaluated by solving the parabolic form of the Navier–Stokes equation. Stream-wise second derivatives and the normal pressure gradient are neglected to avoid elliptic contributions to the system of PDE. Laminar solutions are sought, thus no Reynolds stress is considered.

∇ · (ρ u ) = 0,

ρ u1

∂ u1 ∂ u1 ∂u ∂p ∂ + ρ u2 + = μ 1 , ∂ x1 ∂ x2 ∂ x1 ∂ x2 ∂ x2

∂U ∂U + ρ u2 + p ∇ · u ∂ x1 ∂ x2

2 ∂ u 1 ∂ u 22 ∂ u 1 ∂ u 2 ∂ u2 ∂ q2 4 +μ 1, =− + μ + − ∂ x2 3 ∂ x1 ∂ x2 ∂ x1 ∂ x2 ∂ x2

(7a) (7b)

ρ u1

ρ u1

∂Y ∂Y ∂ v2 + ρ u2 =− − ω˙ . ∂ x1 ∂ x2 ∂ x2

(7c) (7d)

The equations are solved using a steady one-dimensional approximation for the solid phase solution, which yields the following propellant boundary conditions,

ρ u 2 = ρs rb ( T ), 0

(8a)

ρ u 2 h − u 2 σ22 + q2 = ρs rb ( T )h s,298 K ,

(8b)

ρ u 2 Y + v 2 = ρs rb ( T ).

(8c)

The parabolic solution is initialized at the head-end by specifying a co-sinusoidal stream-wise velocity proﬁle u 1 with maximum value equal 1 cm/s, and a sinusoidal cross velocity proﬁle u 2 with wall velocity assigned based upon the strand combustion analysis. The initial temperature and species mass fractions are assigned based upon the strand propellant analysis. The initial pressure is set to 30 atm for all cases.

868

L. Massa / Combustion and Flame 156 (2009) 865–888

2.7. Linearized perturbation equations

3.1. Burning characteristics

Denoting the solution vector, X = [ P 1...6 , T s ], the perturbation expansion about the mean solution

The derivation of the chemical model here described was prompted by the exigence of a simple premixed-combustion model able to describe single-ﬂame composite propellant. Solid propellant combustion modeling contains a set of uncertainties so that the combustion model can seldom be obtained by reduction of a detailed model. Void fraction, solid phase reactions and solid–solid heterogeneous interface behavior are the most outstanding issues. A global reaction scheme is sought for which the parameters of the kinetic model are calibrated by matching a set of burning characteristics. At a given pressure, the set of burning characteristics used to describe propellant combustion in zero cross-ﬂow conditions are:

X = X¯ + X ,

(9)

leads to the linearized system of equations written in the general matrix form,

∂ X ∂ X ij ∂2 X + Z 00 X + Z i0 + Z ∂t ∂ xi ∂ xi x j 3

JX

3

i =1

3

i =1 j =i

+ Z v T s (x2 = − D w ) = 0.

(10)

The ﬁrst six rows and columns of the matrices Z i j are found based on the Navier–Stokes equations, while the (7, 7) element is found from the solid heat conduction equation. The vector Z v multiplies the interface (solid–gas) temperature and has compo∂r nents Z v = [0, . . . , 0, ∂ Tb ∂∂ xT s ] T . The solid occupies the domains s 2 x2 ∈ [−∞, − D w ] and x2 ∈ [ D w , ∞]. Due to the odd-even splitting of the perturbation modes as a consequence of symmetry, only the domain x2 ∈ [−∞, − D w ] is here considered. Furthermore,

JX = The Z

J 0

ij

...

(11)

sub-matrices based on the gas phase equations are

Z 100 ...6,1...6 =

−

0 . 1

3 ∂ S ∂ Ai ∂ P ∂ Ai ∂ P ∂ Ai ∂ P + , ,..., ∂P ∂ P 1 ∂ xi ∂ P 2 ∂ xi ∂ P 6 ∂ xi i =1

3 ∂V i i =1

∂P

3 ∂V i

j0

Z 1...6,1...6 = A j −

i =1

jk

Z 1...6,1...6 = −

(12a)

, P xk , P xl xm

∂ Pxj

3 ∂V i ∂ P x j ,xk P , P x i =1

(12b)

, P , P xk = j , P xl xm

1...3

.

(12c)

, P xl = j xm =k

While the heat conduction equation leads to Z 2, 0 = r b

and

Z i ,i = −K.

(13)

3. Chemistry and transport model The vast majority of rocket related propellants support heterogeneous combustion because a signiﬁcant portion of oxidizer particles have dimensions of the same order as the sublayer thermal scale. Nonetheless, oxidizer particles scales are much smaller than the core characteristic scales, so that there is no rigorous way to include ﬁnite grain effects in the linear analysis. Coarse grain propellants are analyzed by building the homogeneous combustion model on a set of burning characteristics that can be varied to match realistic formulations. By doing so the mutual interaction between chamber disturbances and diffusion limited combustion is a priori discarded, and the analysis addresses convective penetration by focusing only on the relations between core ﬂuctuation times and reaction times. The rationale for such approach is based on the modeling effort of Mukunda and Paul [3], and on data analysis of King [12]. They found virtually identical erosive augmentation values for mixes with different grain size but similar burning characteristics. Thus we can devise a homogeneous propellant that effectively reproduces the core–sublayer interaction of coarse grain propellant by varying its burning characteristics. Homogenization issues in the determination of propellant thermo-chemical properties are discussed in detail by Chen et al. [9].

(a) Burn rate, rb deﬁned as the normal regression rate. ∂r

(b) Cold temperature sensitivity,

σp ≡

( ∂ Tb ) p 0

rb

p ∂ rb ( ) . rb ∂ p T 0

.

(c) Pressure exponent, n p ≡ (d) Solid phase heat release (positive if exothermic), Qs ≡

λs ∂∂ xT s −λ g ∂∂xT 2

ρs rb

2

.

(e) Burn rate surface temperature sensitivity, Θ ≡ p

∂r

2 T surf ( ∂ ∂Trb ) p . rb surf

(f) Burn rate pressure sensitivity, nr ≡ r ( ∂ pb ) T surf ; this last paramb eter has been set to zero. A more detailed description on the value of this parameter in the context of solid propellant combustion is described by Massa et al. [13]. The above listed burning characteristics are assumed to fully describe the propellant burning, and will be varied in order to demonstrate their effect on the propellant chamber interaction. Besides this set, other parameters are important in describing the propellant burning under no-cross ﬂow conditions, most notably the adiabatic ﬂame temperature and its derivatives. The adiabatic ﬂame temperature can be accurately described using a premixed combustion model that includes equilibrium chemistry as main ingredient. The ﬂame temperature and its derivatives ( ∂ T∂ﬂame )T 0 and ( ∂ T∂ﬂame ) p (note in what follows the solid enthalpy p T0 is substituted to T 0 , the two variables are linearly related for a constant solid phase speciﬁc heat) are evaluated given the atomic composition of the propellant mix and its enthalpy of formation. 3.2. Kinetic model The smallest solid-propellant premixed chemistry model able to match the above described burning characteristics is composed of two reactions, (a) Surface decomposition (gasiﬁcation) reaction. (b) Global homogeneous irreversible reaction. The ﬁrst leads to the burn rate law, Eq. (4), and the previously described solid phase heat release. The product of the gasiﬁcation reaction is the gas phase reactant, here denoted with R. The second reaction, ω˙

R −→ P

(14)

leads to the formation of the ﬁnal products. The rate of mass decomposition of reactant is related to the local thermodynamic state through the source term

ω˙ = DY pnr exp(−θ/ T ),

(15)

where Y is the mass fraction of reactant. R is thus not a well deﬁned species but a variable representing the progress of the propellant decomposition.

L. Massa / Combustion and Flame 156 (2009) 865–888

The chemical composition of the products and consequently the heat release in Eq. (15) are determined by solving a chemical equilibrium problem similar to that described by Gordon and McBride in [14]. The product gas P is a mixture of μ gases with mass fractions Y k , k = 1, . . . , μ, molar fractions X k , and average (effective) molar mass M P . Two types of equilibrium computations are performed, those in which the pressure is given and Gibbs free energy is minimized, and those in which the density is given and the Helmholtz free energy is minimized. The major difference from what described in [14] is the presence of the reactant mass fraction (Y ) in both mass and enthalpy/internal energy balances as an (extra) independent variable. The reactant is characterized by assigning its effective molar mass, M R = 30 kg/kmol, the standard enthalpy of formation, hR,298 K , and the atomic molar fractions taken identical to the solid ones. The ﬁrst value is ﬁxed and similar to that considered by Jackson and Buckmaster [15]. The standard enthalpy of formation is determined as part of the calibration procedure to match the value of heat release. The chemical composition of the solid is speciﬁed through the mass percent of oxidizer χAP , in the AP-HTPB mix (blend). The atomic make-up of AP is taken to be NH4 ClO4 and the standard heat of formation of solid state AP, hAP,298 K = −601.67 cal/g. The atomic make-up of the polymeric HTPB binder is C7.075 H10.647 O0.223 N0.063 , and its standard heat of formation, hHTPB,298 K = −137.788 cal/g. The atomic make-up and enthalpy of formation for blends are evaluated by taking mass weighted averages, e.g., h s,298

K

= χAP hAP,298

K

+ (1 − χAP )hHTPB,298 K .

∂ T ﬂame ∂h ∂ T ﬂame ∂p

=

p

1 cp

=− h

(17a)

, 1

ρc p

1+

∂ ln ρ ∂ ln T

,

θ,

Θ and hR,298 K .

μ

(19)

From which the thermo-chemical transport coeﬃcients are obtained as

κT ,T = λ f −

qe

(20a)

,

∇T

∇ Y R =0

qe , ∇ Y R ∇ T =0

R MR W =− , ∇T ∇ Y R =0

R MR W =− , ∇ Y R ∇ T =0

κ T ,Y = −

(20b)

κY , T

(20c)

κY ,Y

(20d)

where λ f is the frozen thermal conductivity (calculated as detailed q

in [16, pp. 24–29]), and the symbol ( ∇eT )∇ Y R =0 indicates the scalar coeﬃcient of ∇ T in the vector expression qe . The remaining expression on the RHS of Eq. (20) are to be understood analogously. Also note that in this section the subscript R has been aﬃxed to the reactant mass fraction, to distinguish it from the other equilibrium species. Considering a stationary mixture at constant pressure, identifying a subset of ν < μ independent equilibrium species (cf. [18, p. 1638]), mass conservation considerations lead to the following relationship ν

i− ni , j W

ν

i =1

R, n i , j ψi + ψ j W

j = ν + 1, . . . , μ,

i =1

(21a) Yi

MR

ψ i = μ

j =1

Y j Mi

,

i = 1, . . . , μ.

(21b)

Denoting the binary diffusion coeﬃcients with Γ , and the associated scaled values,

(18)

3.3. Transport model

R HR. i Hi + W W

i =1

j =− W

∂ ln ρ

nr ,

qe =

(17b)

p ,Y

where the term ( ∂ ln T ) p ,Y is calculated as described in [14, p. 18], and c p is the equilibrium speciﬁc heat. The second derivative, Eq. (17b), would be identically zero if the frozen (perfect) gas model were considered. Finally, it is remarked that the model parameters calibrated to match the burning characteristic above listed are

D,

The thermo-chemical transport properties of the mixture are composed of two contributions, a frozen value and a reaction value. The latter is associated with the molecular ﬂux balances in equilibrium conditions. The theoretical background for the reactive contribution is discussed by Butler and Brokaw [18]. The presence of the non-equilibrium reactant changes the chemical equilibrium contribution in two basic ways, ﬁrst the contribution to the thermal conductivity, κ T , T , is altered; second equilibrium chemistry contributions affect the diffusive terms κ T ,Y , κY , T and κY ,Y introduced in Eq. (1). (W R identiﬁes the molar Denoting the molar ﬂuxes with W ﬂux of non-equilibrium reactant), and the enthalpy per unit mole with H , the heat ﬂux contribution due to chemical equilibrium is

(16)

The adiabatic ﬂame temperature is obtained by setting Y = 0 in the iterative Gibbs equation with pressure and enthalpy (set equal to h s,298 K ) assigned. It is therefore a non-linear function of the oxidizer weight fraction, an improvement to previously used propellant models [13]. For some of the numerical experiments later described the ﬂame temperature will be changed while maintaining the chemical composition ﬁxed. This will be achieved by arbitrarily increasing or decreasing the heat of formation of the mix. The derivatives of the ﬂame temperature are

869

Δ≡

RT

ΓP

,

the Stefan–Maxwell equations become, The two principal advantages of solving for the equilibrium composition are that the ﬂame temperature is exactly evaluated, and that the transport property of the mixture can be approximated in a manner consistent with the chemical composition. The transport properties for the μ equilibrium species forming the product gas, P, are obtained using the method detailed by Svehla and McBride [16]. The transport database is compiled based on the data reported by McBride et al. [17]. The transport cross sections of the (non-equilibrium) reactant, R, are estimated based on its effective molecular mass M R , using the empirical rules also suggested by Svehla and McBride (see [16, pp. 32–33]).

∇ Xk =

μ

j − X jW k) Δ j ,k ( X k W

j =1

R − XR W k ), + ΔR,k ( X k W ∇ XR =

μ

k = 1, . . . , μ,

j − X jW R ), Δ j ,R ( X R W

(22a) (22b)

j =1

∇ Y R = (1 − Y R )∇ X R

MR

ˆ M

−

μ YR

ˆ M

j =1

∇ X j M j.

(22c)

870

L. Massa / Combustion and Flame 156 (2009) 865–888

are shown in Fig. 2b. Both the Lewis and Prandtl number vary substantially across the combustion layer, an observation that suggests constant Lewis number approximations are inaccurate for propellant analysis. 4. Eigenvalue problem 4.1. Local non-parallel approach Following Casalis et al. [10], and, more recently, Abu-Irshaid et al. [19] we consider the local non-parallel instability problem as representative of perturbation growth in a rocket motor. Casalis et al. [10] and Griffond and Casalis [20] proved that the local non-parallel problem provides results consistent with more sophisticated approximations such as the parabolic stability equation, Herbert [21], where a weakly non-parallel perturbation expansion is considered. The wavelike perturbation expansion over a nonparallel local ﬂow approximation of the mean ﬂow,

X = Xˆ (x2 ) × exp i (α x1 + β x3 − ωt ) ,

(a)

(24)

leads to an eigenvalue problem for the complex number α , the negative imaginary part of which, −αi , is the spatial growth rate. 4.2. Discretization: Chebyshev tau method The perturbation equation are discretized using the Chebyshev tau method, Dongarra et al. [22]. The major difference from recent applications of the method, namely Abu-Irshaid et al. [19], is the treatment of the pressure and, more speciﬁcally, the pressure solid boundary condition. Because the highest order of pressure derivatives is one, and the continuity equation contains only ﬁrst-order derivative operators, the pressure expansion contains one fewer term than the other variables and the continuity equation yields one extra discrete equation, as shown in the following. Denoting with N the order of the maximum Chebyshev polynomial resolved in the truncation error, and C the vector of Chebyshev polynomials, (b) Fig. 2. Equilibrium contribution to the thermo-chemical transport properties for a 73% AP propellant strand. (a) Thermal conductivity. Lines: (—) single irreversible reaction with equilibrium product, (– –) chemical equilibrium. (b) Lewis and Prandtl numbers. Lines (single irreversible reaction model): (—) Lewis number, (– –) Prandtl number.

Finally the missing ν equations are provided by the equilibrium constant for the ν virtual reactions associated with each independent species. Using the van’t Hoff equation to determine temperature derivatives of the reaction constants,

Hk RT 2

∇T =

μ nkl l=ν +1

Xl

∇ Xl −

1 Xk

∇ Xk ,

k = 1, . . . , ν .

(23)

Algebraic substitution of Eq. (21) into Eq. (22), and inclusion of such results into Eq. (23), lead to the solution of a linear system of equation of size ν in the W i , i = 1, . . . , ν unknowns [compare to 8, Eq. (9)]. It is easily veriﬁed that the condition X R = 0 reduces the transport problem in Eq. (20) to that illustrated in [18] for an equilibrium mixture without ﬁnite rate reactions. The difference for X R > 0 is illustrated in Fig. 2a, which refers to the steady solution of a 73% AP propellant strand. The thermal conductivity obtained using the approach here presented is lower than the purely equilibrium analog (Butler and Brokaw [18] solution) at the propellant wall, where the reactant mass fraction is maximum. The Lewis and Prandtl numbers for the irreversible reaction case

⎡ N +1 ˇ k ⎤ k=0 X 1 C k ( z2 ) ⎢ N +2 Xˇ k C (z ) ⎥ ⎢ k=0 2 k 2 ⎥ ⎥. Xˆ = ⎢ ⎢ ⎥ .. ⎣ ⎦ . N +2 k ˇX C k (z2 ) 7 k=0

(25)

The Chebyshev expansion variable z2 is obtained by mapping the second Cartesian direction from [− D w , D w ] → [−1, 1] for the gas phase and from [−∞, − D w ] → [−1, 1] for the solid phase. Algebraic non-linear (tenth-order polynomial) mappings are used in order to obtain D1 ∂∂ xz2 = 1 × 10−3 at the surface. A small value w 2 of this derivative facilitates the task of capturing the variation of the solution in the combustion sublayer. Algebraic equations for the Chebyshev components Xˇ k are obtained by considering the associated coeﬃcients in the binary inner product, ·, ·, deﬁned as in Malik and Hooper [23]. The block coeﬃcient matrices of position k, l are

20 22 2 ¯ A = C T Z k00 ,l I + Z k,l D + Z k,l D , C i , i ,1...¯n

(26a)

k,l B¯

(26b)

k,l

T

Z k10,l I

Z k12,l D

= C + , Ci , T 30 23 = C Z k,l I + Z k,l D , C i , k , l ¯ D = C T Z k11,l I , C i , i ,1...¯n k,l E¯ i ,1...¯n = C T Z k33 ,l I , C i , k,l F¯ i ,1...¯n = C T Z k13,l I , C i , k,l G¯ i ,1...¯n = C T J kX,l I , C i i ,1...¯n

k,l C¯ i ,1...¯n

(26c) (26d) (26e) (26f) (26g)

L. Massa / Combustion and Flame 156 (2009) 865–888

¯ and k, l = 1 . . . 7. The row index limit m ¯ is equal for i = 1 . . . m to N + 2 for k = 1 and N + 1 otherwise, while the column index limit is equal to N + 2 for l = 1 and N + 3 otherwise. The last term in Eq. (10) is accounted for by adding to the ith row ¯ 7,7 the ith Chebyshev component of the term Z v . The of matrix A differentiation matrices, D and D 2 , are detailed in [22]. The inner product are easily evaluated analytically once the Z matrix entries are approximated with Chebyshev series. Boundary conditions provide the remaining twelve rows of the coeﬃcient matrices given in Eq. (26). Even and odd modes are sought independently owing to the symmetry of the problem. Separating the modes is a necessity in this analysis because the large length difference between combustion layer and chamber ﬂow length scales requires the use of a large set of Chebyshev polynomials to resolve the modes. The mode splitting leads to the solution of the (non-symmetric) solid heat conduction in only the lower part of the domain and implies that only seven boundary conditions are required. Equations (5) and (6) in order yield, ¯ A ¯ +1,1...¯n = 1, m

(27a)

¯ 4, 4 A ¯ +1,1...¯n m

(27b)

2, 2

= 1,

¯ A ¯ +1,1...¯n = −1 m 7, 7

¯ A ¯ +1,1...¯n = ι, m 7, 5

and

(27c)

∂ rb 2 ¯ 3,l ¯ 3, 7 A , ¯ +1,1...¯n = ιφl A 1,l and A m ¯ +1,1...¯n = −ρs m ∂Ts

∂ G 25 l 2 ¯ 5,l , A γ ¯ +1,1...¯n = ιφl A 5,l − m ∂ P l 1...¯n

∂ rb ¯ 5, 7 A = − ρ h + C ( T − 298 K ) + r C s s , 298 K s s s b ¯ m+1,1...¯n ∂Ts + λs γ17...¯n ,

∂ G 26 l 2 ¯ 6,l A = ι φ A − γ and l 6,l ¯ +1,1...¯n m ∂ P l 1...¯n ∂ rb ¯ 6, 7 A , ¯ +1,1...¯n = −ρs m ∂Ts ¯ 7, 7 ¯ − 1, A ¯ +2, j +1 = 1 − 2(i mod 2) for j = 0 . . . n m

(27d)

(27e)

2i 2 , 2(i + 1)2 ,

i = 0 . . . n¯ i = 0 . . . n¯

i B¯ − β F¯ I

(27g)

if φl × ι = 1 or l > 6, if φl × ι = −1 and l 6.

=α

¯ D 0

−i ω G¯ + A¯ + i β C¯ − β 2 E¯ 0 I

α Xˇ Xˇ

0

.

α Xˇ

The propellant burning problem has three set of scales, those associated with the chamber ﬂow, those associated with the combustion (sublayer) and those associated with the solid phase heat conduction. It has been shown in Section 2.5 that the length scales of the near wall combustion are much smaller than those controlling the chamber ﬂow. Given that the acceleration in the combustion layer is proportional to the ratio between ﬂame temperature and surface temperature, and that such ratio is in the range [2–3], one can anticipate that the same conclusion holds true for the time scales. On the other hand, due to the large density difference between solid and gas phases, the solid heat conduction time scale is of the same order as the characteristic time scales of the chamber ﬂow. For the 73% AP mix in strand conditions T t ≡ K2 = 2.94 × 10−3 s, while the chamber ﬂow time rb

scale is T w ≡ UD w = 2.2 × 10−2 s. Therefore, combustion and solid w heat conduction have comparable length scales, heat conduction and chamber ﬂow have comparable time scales, while combustion has both time and length scales considerably smaller than the chamber core ﬂow. In the limit for which the two latter phenomena can be considered independent, the combustion sublayer can be assumed at steady state and separated from the ﬂow chamber evolution by assuming zero ﬁrst derivative conditions at the sublayer–core boundary. This condition leads to a decoupled core ﬂow–sublayer problem. Decoupled combustion and chamber ﬁelds typify a propellant not supporting erosive burning augmentation, because the latter is identiﬁed with the additional heat feed-back caused by the chamber ﬂow convective penetration into the combustion layer. Convective penetration is determined as the ratio of temperature disturbance in the sublayer and pressure disturbance at the core–sublayer edge. The comparison between fully coupled and decoupled problems will provide information on the propensity to support erosive burning. 5.1. Decoupled problem

The eigenvalue problem is obtained after concatenating the coeﬃcient block matrices and using the wavelike representation of the perturbation, Eq. (24), to obtain derivatives with respect to x1 , x3 and t. The spatial modes are correlated by the quadratic form (see for example [24]),

5. Chamber propellant interaction

(27f)

where ι = 1 for even modes and ι = −1 for odd modes, l = 1 . . . 6. Moreover φ = [1, 1, −1, 1, 1, 1] is an index denoting whether each mean-proﬁle components of the primitive vector P is odd, only u 2 , or even, all other gas phase variables. Finally,

γl =

871

Xˇ

(28)

This generalized eigenvalue problem is solved using a QZ algorithm [25], by mean of which all the spatial eigenvalues and eigenvectors are determined. The determination of the complete spatial spectrum is signiﬁcantly more expensive than iterative methods for the solution of eigenvalue problems, especially for problems in which a dominant eigenvalue exists. Its major advantage is that it requires no iterations and eliminates the problem of converging to nondominant modes.

In the decoupled problem gas and solid phases are solved independently. The only way the chamber ﬂow affects the propellant burning is by virtue of the pressure perturbation at the wall,

p = p˜ × exp i (α x1 + β x3 − ωt ) . The linearized analysis of the propellant burning (including the combustion sublayer and the solid phase) subject to a wave-like pressure perturbation leads to the following problem,

∂ T s ∂T ∂ rb ∂ T¯ s + Ts + r¯b s = ∇ 2 T s , ∂t ∂ T s ∂ x2 ∂ x2 T s = 0 at x2 →= −∞, ∂ T s = G p p + G T T s at x2 →= − D w , ∂ x2

(29a) (29b) (29c)

where G p and G T are found by differentiating the steady, onedimensional, strand solution. Solution of Eq. (29) leads to the deﬁnition of the decoupled pressure response,

Θ , T2 surf

A = ( T surf − T 0 )

T

B=

1

, T σ p K 2 Ω =ω− α + β2 , i Λ = 1/2 1 + 1 − 4i ΩK/rb2 ,

(30a)

(30b) (30c) (30d)

872

L. Massa / Combustion and Flame 156 (2009) 865–888

R˜ p ≡

r˜b p p˜ rb

=

np A B

Λ + A /Λ + A B − (1 + A )

,

(30e)

where only the combustion parameters listed in Section 3.1 plus the surface temperature are involved. Note that in view of the relationship between gas and solid phase length scales and the fact (later discussed) that the perturbation wavelength, α and β , are only slightly changed by the presence of the combustion sublayer, neglecting the second term on the RHS of Eq. (30c) does not affect the pressure response. Therefore, R˜ p is identical to the pressure response of strand propellant to sinusoidal pressure waves. The gas phase for the decoupled problem is evaluated by solving a set of perturbation equations similar to those shown in Eq. (10), but without solid phase and with boundary conditions,

[u 1 , u 3 , Y ] = 0.

ρs rb ρ T u 2 = 2 R p − J 1,1 − J 1,5 p ,

ρ

T =

R p Cs

+

p

p

∂ T ﬂame ∂h

∂ T ﬂame ∂p

(31a)

T p

MIX Units

rb cm/s

np –

σ

K−1

Qs cal/g

Θ

73% AP 80% AP 100% AP

0.747 1.558 0.421

0.910 0.620 0.922

9 × 10−4 9 × 10−4 2 × 10−3

34.76 42.60 65.00

9462 9788 11,000

K

Table 3 Flame and ﬂow properties of pure propellant and homogeneous mixes at 30 atm. MIX Units

T ﬂame K

kg K/kJ

( ∂ T∂ﬂame )p h

( ∂ T∂ﬂame )h p K/atm

Re –

Mach –

73% AP 80% AP 100% AP

1631 2343 1397

0.512 0.499 0.659

0.0712 0.1532 0.2791

5.38 × 103 8.80 × 103 3.64 × 103

3.05 × 10−3 7.65 × 10−3 1.55 × 10−3

(31b)

p

Table 2 Burning characteristics of propellant mixes at 30 atm.

Λ − 1 + i ΩK/( Arb2 ) Λ

p ,

(31c)

h

where Eq. (31b) is obtained by linearizing ρ u 2 = ρs rb , while Eq. (31c) is found by integrating the perturbation equation, Eq. (29), in the constant property solid phase and imposing conservation of energy between the cold solid and the edge of the combustion layer (ECL),

hECL = −

1

ρs rb

ECL −∞

∂U Cs ρ dx2 = − ∂t rb

− D w −∞

∂Ts dx2 . ∂t

(32)

The integration over the gas phase in the second integral of Eq. (32) is neglected by virtue of the assumption,

ρgas ρsolid .

Fig. 3. Temperature mean proﬁles for the 73% AP mix at Dx = 20. Legend: (—) w D w = 2.5 cm, (– –) D w = 1.25 cm, (– ·) D w = 0.625 cm, (· · ·) D w = 5 cm, and (") corresponding strand proﬁles.

Note that − D w is the gas–solid interface location. 6. Propellant mixes Three ﬁne AP-HTPB propellant mixes have been considered. The experimental burn rate data for the model calibration were taken from Atwood et al. [26,27] (burn rate and cold temperature sensitivity, respectively) for the 100% propellant formulation, from Foster and Miller [28] for the 80% mix, and from King [29] for the 73% propellant mix. Experimental cold temperature sensitivity data for heterogeneous mono-propellants vary signiﬁcantly in the available literature. An increased ﬂame temperature associated with the oxidation of the fuel carbon yields a lower value for the burn rate sensitivity, thus the value suggested by Jackson and Buckmaster [15] is used. The surface heat release is determined by a mass weighted average, Q s = χAP Q AP + (1 − χAP ) Q HTPB ,

(33)

where the heat release of the pure solid components are chosen to be, Q AP = 65 cal/g, and Q HTPB = −47 cal/g. The fuel value matches that used in [15], while the oxidizer value was slightly increased from that suggested in [15], Q AP = 100 cal/g, to match the burn rate sensitivity. It is in fact impossible to obtain a high burn rate sensitivity with such a large part (100 cal/g) of the heat of combustion released in the combustion layer. The surface temperature sensitivity data, Θ , are taken from the analysis of Chen et al. [9], where issue concerning the properties of homogeneous monopropellant mixes are discussed. The chemical

composition and the burning characteristics at the nominal pressure of 30 atm are listed in Table 2, while the ﬂame temperature and its derivatives are listed in Table 3. The Reynolds and Mach numbers given therein are based on viscosity and speed of sound evaluated at the edge of the combustion layer. The above mentioned burning characteristics will be denoted as nominal conditions and correspond to the limit for the particle size to go to zero. Four of them, rb , n p , σ and T ﬂame will be arbitrarily varied in order to demonstrated their effect on the chamber-combustion interaction. Each time any of these variables is changed the combustion parameters listed in Eq. (18) are reevaluated. 7. Results 7.1. Mean proﬁle The mean temperature proﬁles are only slightly changed by a change in port size, and for all port sizes are very similar to the strand proﬁle. The non-dimensional mean proﬁles are shown in Fig. 3 at Dx1 = 20. The difference is mainly due to the use of w a non-dimensional coordinate in the abscissa. Plotting proﬁles at larger x1 values does not change the conclusion that the parabolic mean ﬂow has a weak inﬂuence on the burning. The burn rate along the axial distance is changed only because of the drop in pressure associated with the injecting ﬂow. The burn rate is plotted against the axial distance in Fig. 4 for the 73% AP mix. The

L. Massa / Combustion and Flame 156 (2009) 865–888

Fig. 4. Ratio between steady and strand (rb,s ) burn rate against the axial coordinate p for the 73% AP mix. Legend: (—) steady parabolic solution, and (") rb,s ( 30 atm )n p .

873

(a)

mean burn rate proﬁle agrees well with the steady burn rate law, p rb = rb,s ( 30 atm )n p , where rb,s is the strand burn rate at 30 atm. Therefore the pressure exponent is preserved along the streamwise direction. Lines drawn for different port diameters are indistinguishable on the scale of Fig. 4, which conﬁrms the weak correlation between steady parabolic ﬂow and burn rate. Temperature proﬁles at Dx1 = 20 (D w = 2.5 cm) are shown in w Fig. 5a for the three propellant cases, 73%, 80% and 100% AP. The differences between the parabolic proﬁles and the strand proﬁles are again marginal. The associated axial velocity proﬁles are shown in Fig. 5b. The difference in velocity proﬁles is very small, and a comparison with the Taylor–Culick proﬁle, shown with symbols, demonstrates that the Taylor–Culick proﬁle is reproduced by the combustion solver. 7.2. Eigenvalue veriﬁcation The solution to the eigenvalue problem was veriﬁed in three ways. First, the Orr–Sommerﬁeld modes supported by the twodimensional Couette ﬂow are compared to those reported by Dongarra et al. [22], second the local-non-parallel modes supported by the Taylor–Culick proﬁle for an incompressible ﬂow are compared to those reported by Casalis et al. [10], third the solution to the pressure response problem, Eq. (30), is obtained for zero cross ﬂow and imposed pressure perturbation in the combustion layer. This last veriﬁcation was necessary to test the propellant connection conditions and the combustion terms. All the solutions presented are carried out with 120 even or odd modes, and in all cases the method has been veriﬁed to be converged by comparing results for 90 and 120 polynomials. 7.3. Growth rate The number of unstable modes varies with the stream-wise coordinate, x1 . The number of odd and even modes is identical. The maximum growth rate for the even modes is plotted against the stream-wise coordinate in Fig. 6 for the 73%, 80% and 100% AP propellant cases. The ﬁgure shows that there is a maximum of six unstable even eigenmodes (there are as many odd modes) in the range 0 < x1 / D w < 40. For all modes, the growth rate eigenvalue −αi increases with the stream-wise direction. The modes become unstable more upstream for a faster burning propellant. This phenomenon is linked to the Reynolds number effect through the analysis that follows. There is a threshold value of the stream-

(b) Fig. 5. (a) Temperature and (b) stream-wise velocity mean proﬁles for three propellant mixes at Dx1 = 20. Legend: (—) 73% AP, (– –) 80% AP, (– ·) 100% AP and (") w strand proﬁle (panel a)/Taylor–Culick (panel b).

wise coordinate below which the ﬂow is stable. The threshold is fairly independent of the propellant composition. The pressure and velocity (u 1 and u 2 components) eigenfunctions for the six even and odd modes are plotted at x1 / D w = 36 (such value was chosen because the six modes appear at a single frequency) in Fig. 7, for the values of ω and β that maximize the most unstable mode growth rate. The temperature and reactant eigenfunctions will be investigated in detail in the next section. The modes are sorted according to their maximum growth rate −αi in descending order, starting from the least stable. The maximum velocity and pressure perturbation shift toward the burning surface as the mode number is decreased. The difference between even and odd eigenfunctions is accentuated for large mode numbers. Note that the even-odd mode splitting is based on all solution components but the normal velocity, for which such characterization is inverted: i.e., an even mode for pressure, tangential velocity, temperature, and mass fraction features an odd variation in the normal velocity. The remainder of this section focuses on results at Dx1 = 20 w and for the nominal propellants: i.e., the three mixes described

874

L. Massa / Combustion and Flame 156 (2009) 865–888

(a)

(b)

(c) Fig. 6. Maximally ampliﬁed even mode eigenvalues against the stream-wise coordinate. Lines of different types indicate different modes. (a) 73% AP. (b) 80% AP. (c) 100% AP.

(a)

(b)

Fig. 7. Even/odd modes eigenfunctions at x1 / D w = 36. Lines: (—) p, (– –) u 1 (multiplied by 5), and (– ·) u 2 (multiplied by 5). The panels are ordered left to right top to bottom in descending instability order. Only solutions in the bottom half of the domain are plotted. (a) Even modes. (b) Odd modes.

in Table 2. The most ampliﬁed even growth rates −αi in the β ∈ [0, ∞] range are plotted against the non-dimensional ﬂuctuation frequency ω in Fig. 8. Note that even though −αi in Fig. 8 is non-dimensional the length scale is identical for the three cases shown. All modes have a monomodal variation. The frequency at which the growth rate is maximum will be denoted peak fre-

quency in the remainder of this paper. The ﬁrst mode is strongly dominant over the remaining three, while differences among propellant mixes are more marked at high mode numbers. Therefore, the instability development is expected to be fairly independent of the propellant speciﬁcation and well represented by the ﬁrst mode. Even and odd modes eigenvalues are very similar. The difference

L. Massa / Combustion and Flame 156 (2009) 865–888

875

(a)

(b)

(c)

(d)

Fig. 8. Unstable even modes in nominal mixes; maximum growth rate eigenvalue −αi versus frequency (b) Second mode. (c) Third mode. (d) Fourth mode.

between even and odd growth rates is plotted in Fig. 9 against the non-dimensional frequency. The difference is more evident for low Reynolds numbers ﬂows, 100% AP, and low frequency ﬂuctuations; it is relatively more important for less unstable modes. At the maximum ampliﬁed frequency the difference between odd and even eigenmodes is negligible. This conclusion is true for all cases examined. The difference between associated eigenfunctions is also negligible. Thus only even modes will be considered in the reminder of this section. Three-dimensionality of the perturbation ﬁeld is investigated by considering the relationship between growth rate and β parameter at constant frequency. For the three cases discussed in this section, β = 0 identiﬁes the maximally ampliﬁed unstable modes. The variation of the ﬁrst even mode growth rate eigenvalue against β at its peak frequency is shown in Fig. 10 for the three propellant mixes. In summary the eigenvalue analysis reveals that the difference between propellant mixes is not marked. For all modes the fastest burning case 80% AP is also the most unstable. This observation is in apparent contrast with the experimental observation that erosive burning decreases with an increased burn rate, see for example Hasegawa et al. [8] and discussion therein.

ω . Lines: (—) 73% AP, (– –) 80% AP, and (– ·) 100% AP. (a) First mode.

7.4. Coupling effect on the growth rate The notion that the perturbation growth is weekly affected by the burning is reaﬃrmed by comparing coupled and decoupled (in the sense described in Section 5.1) solutions. The difference between the decoupled and coupled growth rates for the for the ﬁrst two eigenmodes is shown in Fig. 11. The decoupled solution is obtained using Eqs. (30) and (31) along with a mean ﬂow evaluated using the Taylor–Culick proﬁle for velocity and pressure, and a constant temperature equal to the adiabatic ﬂame temperature. The 80% AP decoupled solution differ the most from the coupled analog. The difference is due to the effect of the mean-proﬁle on the problem rather than the modiﬁcations to the perturbation equation listed in Eqs. (30) and (31). If the mean-proﬁle is kept ﬁxed and equal to the coupled case, the maximum difference in ﬁrst mode growth rate is,

|αi , I ,coupled − αi , I ,decoupled | < 2 × 10−4 , so the two are almost identical. A similar analysis carried out by selectively substituting components of the combustion meanproﬁle into the Taylor–Culick solution reveals that the reason why the deviation is most pronounced in the 80% case is the large change in temperature across the combustion layer. The surface

876

L. Massa / Combustion and Flame 156 (2009) 865–888

(a)

(b)

(c)

(d)

Fig. 9. Difference between even and odd modes growth rates in nominal mixes. Lines: (—) 73% AP, (– –) 80% AP, and (– ·) 100% AP. (a) First mode. (b) Second mode. (c) Third mode. (d) Fourth mode.

higher than for the other two cases, cf. Table 3. Note that a negative value of the variable plotted in Fig. 11 denotes a destabilizing effect of the combustion proﬁle. 7.5. Perturbation burn rate The eigenvalue analysis so far has indicated that the combustion–chamber ﬂow coupling has a weak inﬂuence on the growth rate, and the burning characteristics (see burn rate) inﬂuence on the perturbation growth rate is contrary to the observed relationship between erosive burning and combustion characteristics. A link between experimental measurements and linear analysis results is sought by considering burn rate eigenfunctions. The importance of coupling between core ﬂow and combustion sublayer on the perturbations is analyzed by comparing fully coupled and decoupled eigenfunctions. The analysis focuses on the perturbation response function Rˆ p ≡ Fig. 10. Growth rate against β for three nominal mixes, cf. Eq. (24). Lines: (—) 73% AP, (– –) 80% AP, and (– ·) 100% AP.

temperatures among the three formulations are very similar, while the adiabatic ﬂame temperature for the 80% case is substantially

rˆb p pˆ rb

,

and its relation with the decoupled analog deﬁned in Eq. (30e). The two are plotted against ω in Fig. 12. The thin lines refer to the coupled solution, while the thick lines represent the decoupled one, R˜ p . Note that for ω = 0 all response functions approach the pressure exponent n p , and that the lines in Fig. 12 cover only

L. Massa / Combustion and Flame 156 (2009) 865–888

877

(a) (a)

(b) Fig. 11. Difference between coupled and decoupled problem growth rates in nominal mixes. Lines: (—) 73% AP, (– –) 80% AP, and (– ·) 100% AP. (a) First mode. (b) Second mode.

the range of ω that supports unstable modes. The plots show that there is a signiﬁcant effect of the combustion parameters on the difference between coupled and decoupled eigenfunctions, that the most ampliﬁed eigenmodes have consistently stronger perturbation responses than weaker modes, and ﬁnally that for small value of ω the curves for all coupled modes approach rapidly the decoupled response. Large differences in response functions denote a strong interaction between chamber ﬂow and combustion processes in the combustion sublayer. The variation of the propellant response (at the peak frequency and for the most ampliﬁed mode) with the axial coordinate is shown in Fig. 13. The stream-wise decrease of Rˆ p is due to the increase of the peak frequency and the consequent loss of resonance between core ﬂow and solid phase modes. This phenomenon is more marked for the fastest burning mix, 80% AP. The ratio between coupled and decoupled responses increases monotonically with stream-wise distance, cf. Fig. 13b. Thus the interaction between core and sublayer processes becomes stronger as x1 increases. In the context of a local non-parallel approach, the streamwise coordinate can be interpreted as the port mass ﬂow rate; in fact a marginal variation of all solution components but u 1 was noted in the mean proﬁle. Therefore the results of Fig. 13b validate the link between coupling and erosive burning by comparison

(b)

(c) Fig. 12. Comparison of perturbation pressure responses of coupled ( Rˆ p ) and decoupled ( R˜ p ) problems, for nominal propellants. Thin lines: (—) ﬁrst mode, (– –) second mode, (– ·) third mode, and (· · ·) fourth mode. Thick solid line: R˜ p . (a) 73% AP. (b) 80% AP. (c) 100% AP.

878

L. Massa / Combustion and Flame 156 (2009) 865–888

(a)

(b)

Fig. 13. Perturbation pressure response against the stream-wise coordinate. Lines: (—) 73% AP, (– –) 80% AP, and (– ·) 100% AP. (a) | Rˆ p |. (b) | Rˆ p / R˜ p |.

with experimental measurements (see for example Marklund and Lake [30]). In order to clarify the nature and the extent of the coupling interaction, the temperature and reactant mass fraction eigenfunctions are analyzed for both the coupled and the decoupled models. The analysis considers only the maximally ampliﬁed even mode at the peak frequency, where, as mentioned above, even and odd solutions are identical. The decoupled perturbation problem, Eq. (29), can be solved analytically once the Jacobian matrix of the steady combustion sublayer solution is evaluated using a ﬁnite difference formula. The solution of the decoupled problem is composed of two parts, one representing perturbations in the combustion sublayer, the second in the chamber ﬂow. The decoupled sublayer solution is compared against the coupled analog in Figs. 14–16, and is extended to the whole domain by virtue of zero gradient boundary conditions. The solid solution occupies the region x2 −2.5 cm. The temperature eigenfunction in the core ﬂow obtained by solving Eq. (30) is compared to the coupled analog in Fig. 17. For this comparison the mean proﬁle was maintained ﬁxed to avoid the incongruence previously discussed. Fig. 17 demonstrates that the temperature eigenfunctions for the coupled problem can be split into two regions to which correspond two peaks in temperature magnitude; therefore the eigenfunctions have essentially a bimodal proﬁle. A ﬁrst peak is associated with the chamber ﬂow process, a second peak occurs in the region close to the burning interface and is associated with the combustion process. The ﬁrst peak is largest in magnitude for the fastest burning mix, 80% AP, while the second is largest for the 73% mix. The sublayer temperature is strongly dependent on the burning characteristics, and a comparison between coupled and decoupled problems links the phenomenon to convective penetration. The 73% and 100% formulations features a sublayer perturbation signiﬁcantly higher than the decoupled analogs. The interaction is much weaker for the fast burning mix, 80% AP, where the decoupled temperature and mass fraction eigenfunctions match well the coupled analogs, cf. Figs. 15b and 15d. The increased level of interaction leads to a larger perturbation of the surface temperature and consequently a substantially stronger burn rate eigenfunction; compare Figs. 12a and 12c to Fig. 12b. For all the cases, the decoupled chamber ﬂow solution is monomodal (the combustion peak is clearly absent) and matches well the coupled temperature eigenfunction away from the wall, cf. Fig. 17. The match is most accurate for the 80% AP. The fact that the decoupled analysis is most accurate both in the chamber ﬂow

and the combustion sublayer for the 80% AP is expected based on the fact that the combustion process thermal length scale D t is inversely proportional to the burning rate and that the decoupling is based on the combustion length scale being much smaller than the chamber ﬂow analog. 7.6. Normal velocity correlations Large temperature and mass fraction ﬂuctuations in the combustion layer noted for the slow burning propellants (see Figs. 14b and 16b) support the idea of disturbance penetration in the ﬂame structure. Convective penetration terms are related to non-linear velocity correlations, the most important of which were identiﬁed by Beddini [5] as u 1 u 2 and h u 2 . The form of these correlation terms is also important in the context of turbulence enhanced transport models [2,4], where the two are related by assuming a constant turbulent Prandtl number. The three works above mentioned ([2,4,5]) have focused on the fully development turbulent region. If we assume that the ﬂuctuation proﬁle remains similar to the linear regime, we can estimate the correlation terms using eigensolutions. The terms are calculated based upon linear ﬂuctuations by integrating over the spatial period, x1,0 + 2π /αr

u 1 u 2 dx1 x1,0

x1,0 + 2π /αr

h u 2 dx1 ,

and x1,0

where x1,0 is a reference location. The integral is proportional to the eigenvalue, which was shown to be almost identical for the three mixes, plus an x2 varying term of magnitude Re(uˆ 1 uˆ ∗2 ) (Re(hˆ uˆ ∗2 ) for the enthalpy term), where the superscript ∗ denotes the complex conjugate. These terms are plotted in 18 and 19. The abscissa extension of the insert plot is equal to 100 μm about four times D t , cf. Section 2.5. For the three propellants shown, the amount of both correlations in the sublayer is small compared to its value in the core ﬂow ∂ u

because of the small value of both u 2 and ∂ x 2 ; the latter is a result 2 of the almost incompressible conditions. The stream-wise velocity correlation u 1 u 2 does not change signiﬁcantly among the three cases, thus is substantially independent of burning parameters. The enthalpy correlation instead manifests the presence of a positive secondary peak in the near wall region only for the two slowest burning propellants (73% AP and 100% AP). The magnitude of the

L. Massa / Combustion and Flame 156 (2009) 865–888

879

(a)

(b)

(c)

(d)

Fig. 14. Comparison of coupled and decoupled temperature and mass fraction eigenfunctions for the 73% AP nominal mix. The solid–gas interface is at x2 = −2.5 cm. Lines: (—) coupled, (· · ·) decoupled. (a) Temperature in the lower half domain. (b) Temperature near the interface. (c) Reactant mass fraction in the lower half domain. (d) Reactant mass fraction near the interface.

peak increases with a decrease of nominal burning rate showing an increase of penetration in the ﬂame zone. The plots show a trend similar to those reported by Beddini [5] for a much simpler combustion model. The location of the secondary peak seems to be disconnected with the maximum in mean temperature gradient. 8. Non-nominal conditions The analysis presented in the previous section reveals a signiﬁcant change in the level of interaction between chamber ﬂow and combustion when the chemical composition of the mix is changed. In order to determine the role of each combustion characteristic in controlling the interaction, any of such characteristics are changed with the chemical composition ﬁxed. Doing so implies that the burning characteristics of the mix are not linked to experimental observations for ﬁne AP-HTPB propellant, such conditions will be denoted non-nominal. Non-nominal conditions can be experimentally obtained, for instance, by increasing the oxidizer particle size, while still ensuring that the particle diameter is much smaller

than the chamber ﬂow length. Povinelli et al. [31] use this approach to vary the propellant combustion characteristics in their experimental investigation of transverse mode instability in solid propellant rockets. Buckmaster et al. [32] compute the response function of heterogeneous propellant with particle sizes comparable to the thermal length but smaller than chamber scales. The same analysis leads to the deﬁnition of average burn rate and sensitivities of heterogeneous propellant. 8.1. Burn rate effect The effect of the burn rate on the chamber–combustion interaction is studied by modifying the burn rate of the 80% mix to the match the burn rate of the 73% mix (denoted as 800.749 ) and that of the 100% AP, 800.419 . The growth rates for the most unstable mode are compared in Fig. 20. This ﬁgure demonstrates that the larger growth rate previously noted in Fig. 8 for the 80% case when compared to the two other mixes is to be attributed to the higher blowing velocity, and thus to the increased Reynolds num-

880

L. Massa / Combustion and Flame 156 (2009) 865–888

(a)

(b)

(c)

(d)

Fig. 15. Comparison of coupled and decoupled temperature and mass fraction eigenfunctions for the 80% AP nominal mix. The solid–gas interface is at x2 = −2.5 cm. Lines: (—) coupled, (· · ·) decoupled. (a) Temperature in the lower half domain. (b) Temperature near the interface. (c) Reactant mass fraction in the lower half domain. (d) Reactant mass fraction near the interface.

ber. Such Reynolds number destabilization should be regarded as a minor effect and it will be shown not dependent on the coupling. Casalis et al. [10], point out that increasing the Reynolds number or moving the test section downstream result in a diminishing sensitivity of the growth rate to the Reynolds number. The instability mechanism is in effect mainly inviscid. On the other hand a decrease in burn rate has a very strong effect on the combustion chamber ﬂow interaction. The ratio between coupled and uncoupled pressure response,

| Rˆ p / R˜ p | is shown in Fig. 21. A comparison between the temperature eigenfunctions for coupled and decoupled combustion is shown in Fig. 22. These results reveal that for the lowest burn rate case, 800.419 , the peak temperature eigenfunction inside the combustion layer is larger than that associated with the chamber ﬂow process, while for the nominal case the opposite is true, cf. Fig. 15a. Note that temperature eigenfunctions are calculated according to the condition, | pˆ |x2 =− D w = 1, where the pressure perturbation is nondimensionalized by (ρ¯ u¯ 2 )2x =ECL . Because u 2,x2 =ECL ∝ rb , it is not 2

possible to compare peak perturbation values of Tˆ across cases, e.g., values in Figs. 22a, 22b, and 15a. Both Figs. 21 and 22 indicate that a decrease in burn rate has a very strong effect on the coupling, and thus on the burn rate perturbation at the wall. The analysis therefore agrees with the experimental observations (e.g., King [33] and Hasegawa et al. [8]) that link high erosive burning to low nominal burn rate. 8.2. Pressure exponent effect The effect of the pressure exponent is analyzed by varying the combustion wave parameters of the 73% mix. The ﬁrst nonnominal case, 730.62 , matches the 80% pressure exponent. Since the pressure exponents supported by the 73% and the 100% materials are very close together, the second case here introduced has n p = 1.2 and is denoted 731.2 . The pressure exponent has a substantially weaker inﬂuence on perturbation growth rate than nominal burn rate; lines of the type shown in Fig. 20 are virtually indistinguishable on the scale of the

L. Massa / Combustion and Flame 156 (2009) 865–888

881

(a)

(b)

(c)

(d)

Fig. 16. Comparison of coupled and decoupled temperature and mass fraction eigenfunctions for the 100% AP nominal solid. The solid–gas interface is at x2 = −2.5 cm. Lines: (—) coupled, (· · ·) decoupled. (a) Temperature in the lower half domain. (b) Temperature near the interface. (c) Reactant mass fraction in the lower half domain. (d) Reactant mass fraction near the interface.

same ﬁgure. To a smaller degree, the same is true for the effect of the pressure exponent on the burn rate eigenfunction; the ratio | Rˆ p / R˜ p | is shown in Fig. 23. A decreased pressure exponent increases the ratio between coupled and uncoupled pressure response. The relation between burn rate eigenfunction and pressure exponent explains the effect of particle size on erosive burn rate augmentation. It is well known that an increase in particle size decreases the pressure exponent. Experimental measurements of erosive burning conducted by King [12,33] and Razdan and Kuo [34] demonstrate that an increase in particle size has a positive effect on the burn rate augmentation. The phenomenon is strongly linked to a decrease in base burn rate. Comparison of formulations with similar burn rates but different particle size, e.g. 4685 and 4869 at p ≈ 80 atm in [12], shows that decreasing the pressure exponent (n p ,4869 ≈ 0.4 and n p ,4685 ≈ 0.8) slightly increases the burn rate sensitivity. This experimental observation therefore is in agreement with the results presented in Fig. 23.

8.3. Cold temperature sensitivity effect Two additional cases are considered to investigate the effect of σ p on the chamber–ﬂow combustion coupling. The ﬁrst case, 731.5×10−3 , is obtained from the 73% by changing the cold temperature sensitivity to 1.5 × 10−3 , while the second, 808×10−4 , is obtained by setting σ p = 8 × 10−4 K−1 for the 80% mix. A reduction in sensitivity affects the combustion wave by reducing the effective activation energy of the gas reaction, θ in Eq. (15). For instance, θ73 = 3.287 × 103 K and θ731.5×10−3 = 9.700 × 103 K, while

θ80 = 3.130 × 103 K and θ808×10−4 = 2.724 × 102 K. Low activation

energy ﬂames present thicker thermal gradients, so that one expects an augmented interaction with the chamber ﬂow. The effect of σ p on the growth rate eigenvalue is, as in the pressure exponent case, only marginal and indistinguishable if plotted. The ratio between coupled and uncoupled responses is shown in Fig. 24. An increase in propellant sensitivity yields a decreased pressure response for both propellant mixes. A comparison between coupled and decoupled temperature eigenfunctions is

882

L. Massa / Combustion and Flame 156 (2009) 865–888

(a)

(b)

(c) Fig. 17. Comparison of temperature eigenfunctions of coupled and decoupled problems for nominal condition propellants. Lines: (—) coupled problem, (– –) decoupled problem. The insert plot is a magniﬁed section of the interface region, x2 ≈ −2.5 cm. (a) 73% AP. (b) 80% AP. (c) 100% AP.

(a)

(b)

(c) Fig. 18. Comparison of u 1 u 2 non-linear correlation terms based on eigenfunctions for three nominal propellants. The correlations are non-dimensional. The insert plot is a magniﬁed section of the interface region, x2 ≈ −2.5 cm. (a) 73% AP. (b) 80% AP. (c) 100% AP.

L. Massa / Combustion and Flame 156 (2009) 865–888

(a)

883

Fig. 20. Growth rate against ω for three 80% AP mixes with decreasing burn rates. Lines: (—) 80% AP, (– –) 800.749 , and (– ·) 800.419 .

(b) Fig. 21. Pressure response against ω for three 80% AP mixes with decreasing burn rates. Lines: (—) 80% AP, (– –) 800.749 , and (– ·) 800.419 .

shown in Fig. 25. The reduced chamber ﬂow–combustion interaction associated with an increased sensitivity leads to a decreased difference between coupled and decoupled eigenfunctions in the combustion sublayer, cf. Fig. 14b. 8.4. Adiabatic ﬂame temperature effect

(c) Fig. 19. Comparison of h u 2 non-linear correlation terms based on eigenfunctions for three nominal propellants. The correlations are non-dimensional. The insert plot is a magniﬁed section of the interface region, x2 ≈ −2.5 cm. (a) 73% AP. (b) 80% AP. (c) 100% AP.

The adiabatic ﬂame temperature is not usually included among the burning characteristics because its value is ﬁxed given the chemical composition of the mix. Nonetheless, the ﬂame temperature can be greatly changed by adding a small quantity of metals (e.g., aluminum) to the propellant mix. In this section its value is varied at constant composition as discussed in Section 3.2. The solid mix heat of formation is increased of 300 cal/g to yield case 73300 and of 600 cal/g to yield case 73600 . In the ﬁrst case the adiabatic ﬂame temperature is T ﬂame = 2262 K, while in the second case is T ﬂame = 2815 K. The effect of an increase in ﬂame temperature on the growth rate eigenvalue is shown in Fig. 26. The increased temperature has a weak but noticeable stabilizing effect on the perturbation growth.

884

L. Massa / Combustion and Flame 156 (2009) 865–888

(a)

(b) Fig. 22. Comparison of coupled and decoupled temperature eigenfunctions for two 80% AP mixes with decreasing burn rates. Lines: (—) coupled, (· · ·) decoupled. The insert plot is a magniﬁed section of the interface region, x2 ≈ −2.5 cm. (a) 800.749 . (b) 800.419 .

The ratio between coupled and uncoupled pressure responses is shown in Fig. 27. The ﬂame temperature has a weak effect on the coupling. Both augmented ﬂame temperature mixes show an increased response, but the T ﬂame = 2262 K case has a pressure response larger than the T ﬂame = 2815 K case. The results in Fig. 27 showing that variations in adiabatic ﬂame temperature in the range 2000–3000 K have a much weaker effect on combustion layer perturbations than burn rate variations are in agreement with the experimental ﬁndings of King [12,33], Burick and Osborn [35], and Marklund and Lake [30]. These erosive burning measurements showed the ﬂame temperature effect to be weak for formulations with comparable burn rates, e.g. formulations 5565 and 4525 in [12]. 8.5. Half port height effect Three additional cases have been investigated to determine the inﬂuence of the half port height on the coupling. These are iden-

Fig. 23. Pressure response against ω for three 73% AP mixes with varying pressure exponents. Lines: (—) 73% AP, (– –) 730.62 , and (– ·) 731.2 .

Fig. 24. Pressure response against ω for two 73% AP mixes and two 80% AP mixes with varying cold temperature sensitivity. Lines: (—) 73% AP, (– –) 80% AP, (– ·) 731.5×10−3 , and (· · ·) 808×10−4 . The insert plot shows a close detail of the two 80% AP solutions close to the peak frequency marked on the abscissa with a “.”

tiﬁed as 730.625 , 731.25 and 735.0 , the subscript indicating the half port size in centimeters. The port height has a marked inﬂuence on the growth rate. Fig. 28 shows the growth rate for the most unstable even mode. Smaller heights lead to weaker instability; the gain in growth rate obtained doubling the height decreases signiﬁcantly as the height itself is increased. This implies that convective terms become quickly dominant in the instability equations. The effect of the height on the coupling is also very strong compared to the previously described parameters. The ratio of coupled and decoupled pressure response is shown in Fig. 29. The gain in wall eigenfunction magnitude is likely to off-set the decreased growth rate. These computations validate the experimental analysis of Hasegawa et al. [8], where the effect of port size on the erosive burning was shown to be very signiﬁcant. The eigenfunctions are shown in Fig. 30 for the two limiting half port heights D w = 0.625 cm and D w = 5 cm. As for previous cases the increased interaction with respect to the decoupled case is associated with a stronger ﬁrst eigenfunction peak, that associated with the combustion sublayer; the ﬁrst peak for the 730.625 is roughly 5 times the 735.0 analog.

L. Massa / Combustion and Flame 156 (2009) 865–888

885

Fig. 27. Pressure response against ω for three 73% AP mixes with increasing ﬂame temperature. Lines: (—) 73% AP, (– –) 73300 , and (– ·) 73600 . Fig. 25. Comparison of coupled and decoupled temperature eigenfunctions for the 731.5×10−3 case. Lines: (—) coupled, (· · ·) decoupled. The insert plot is a magniﬁed section of the interface region, x2 ≈ −2.5 cm.

Fig. 28. Growth rate against ω for four 73% AP mixes with increasing port height. Lines: (—) 73% AP, (– –) 731.25 , (– ·) 730.625 , and (· · ·) 735.0 . Fig. 26. Growth rate against ω for three 73% AP mixes with increasing ﬂame temperature. Lines: (—) 73% AP, (– –) 73300 , and (– ·) 73600 .

8.6. Comparison with experimental trends This section describes an order of magnitude comparison between linear analysis and King’s [12] experimental data. The experiments were performed in a rectangular channel of initial section 1.90 cm × 1.90 cm containing a composite HTPB-AP propellant. Hot combustion products were produced in a driver section, and erosive augmentation was correlated with the ﬂow speed in the test section. Eight propellant mixes of different stoichiometry and grain size were analyzed in the pressure range [10–100] atm. The comparison with the linear analysis focuses on propellant formulations 4685 (73/27 AP/HTPB, 5 μm AP), 4525 (73/27 AP/HTPB, 20 μm AP) and 5051 (73/27 AP/HTPB, 200 μm AP) at 30 atm and cross ﬂow velocity of 180 m/s (the minimum datum). The goal of the comparison is to support the argument that convective penetration in the sublayer is responsible for the burn rate augmentation, it is not to validate the analysis. The comparison is not rigorous for three reasons. First and foremost, measured augmentation is the results of non-linear interaction. In this section, experimental val-

Fig. 29. Pressure response against ω for four 73% AP mixes with increasing port height. Lines: (—) 73% AP, (– –) 731.25 , (– ·) 730.625 , and (· · ·) 735.0 .

886

L. Massa / Combustion and Flame 156 (2009) 865–888

(a)

(a)

(b)

(b)

Fig. 30. Comparison of coupled and decoupled perturbation temperature eigenfunctions for two 73% AP mixes with different half port heights. Lines: (—) coupled, (· · ·) decoupled. The insert plot is a magniﬁed section of the interface region. (a) 730.625 . (b) 735.0 .

Fig. 31. Comparison of linear analysis results with King’s [12] data. (a) Burn rates; symbols: (") experiments, (2) computations. (b) Temperature eigenfunctions; lines: (—) 5051, (– –) 4525, and (– ·) 4685.

ues of Δr ≡ (rb,180 − rb )/rb are scaled by the corresponding value for propellant 4685, and compared to the eigenfunction |ˆrb |/rb scaled in the same fashion. The scaling is necessary due to the linear nature of the computational results; rb,180 is the burn rate for 180 m/s cross-ﬂow. Second, only 4685 and to a lesser extent 4525 can be considered homogeneous mixes at 30 atm. This second issue is of lesser consequences than the ﬁrst one. Albeit the presence of diffusion effects has strong inﬂuence on linear regression speeds, it affects the erosive phenomenon mainly by changing the burn rate. The physical explanation lies in the large difference between characteristic ﬂuctuation lengths and the particle scales. This observation is supported by King’s comparison of different size/composition mixes, and by Mukunda and Paul [3] analysis of a vast pool of experimental results from different authors. Thus, changes in grain size affect the erosive behavior by changing burning characteristics. Third, the experimental set-up with driver section generated cross ﬂow is different from the double-slag geometry detailed in Section 2.1. In the linear analysis, mean ﬂows

are obtained for the three propellant formulations, and proﬁles are extracted at the stream-wise location where the centerline value of the stream-wise velocity is equal to 180 m/s. Such a location is obviously different for the three formulations. The difference in set-up is not likely to affect the order of magnitude comparison. In fact, the analysis presented in [3] shows that the augmentation can be essentially linked to the port ﬂow for a very various set of core-ﬂow conditions. The three propellant formulations are speciﬁed with different burn rates and pressure exponents, but are taken to have the same sensitivity, 9 × 10−4 K−1 . Increasing the AP size decreases the burn rate and the pressure exponent, both of which were shown to have a positive inﬂuence on the augmentation. Normalized burn rate augmentations Δr are plotted against the burn rate in Fig. 31a. The computed temperature eigenfunctions for the three formulations are compared in Fig. 31b. The comparison shows that the linear interaction analysis predicts combustion sublayer changes that scale with the no-crossﬂow burn rate in a manner consistent with King’s experiments. The use of burn rate eigenfunctions in

L. Massa / Combustion and Flame 156 (2009) 865–888

Fig. 31a is corroborated by the fact that the ratio between peak temperatures in Fig. 31b is of the same order of magnitude as the rˆb variations shown in Fig. 31a. Thus, the sublayer thermal ﬁeld features an increase in perturbation that is consistent in magnitude with measurements. 9. Conclusions The effect of propellant burning characteristics on the coupling between combustion sublayer and core ﬂow in rocket motors is investigated. At steady state, the interaction between the two processes is weak and a Taylor–Culick proﬁle describes accurately the core ﬂow. Unsteady coupling is analyzed by solving the threedimensional eigenproblem associated with a local-non-parallel expansion of the perturbation ﬁeld. The spatial linear stability analysis demonstrates a weak inﬂuence of combustion characteristics on the perturbation growth rate. The inﬂuence is essentially due to the Reynolds number variations among cases, rather than to the coupling. The burn rate and the adiabatic ﬂame temperature have the strongest effect on the growth rate among the analyzed burning characteristics; the ﬂow is stabilized as a result of an increase in ﬂame temperature, or a decrease in burn rate. Although the growth rate is not signiﬁcantly affected by, the thermo-chemical perturbations in the sublayer are markedly dependent on the burning characteristics. Penetration of core-ﬂow ﬂuctuations into the combustion layer is analyzed by comparing the coupled problem to a thermally decoupled problem where the interaction is limited to imposing the pressure at the edge of the core to the whole sublayer. The decoupled case typiﬁes conditions for which chamber length scales are much larger than sublayer ones. The comparison conﬁrms the weak inﬂuence of combustion scales on the core ﬂow instability: coupled and decoupled cases with ﬁxed mean proﬁle show very similar spatial growth rates. On the other hand, temperature and mass fraction eigenfunctions are strongly impacted by the coupling. Both solutions (coupled and decoupled) manifest a bimodal temperature magnitude distribution in the gas phase, with a peak in the sublayer and one in the core region. The thermo-chemical perturbations of the two problems closely agree for burning conditions non-conducive of erosive augmentation: high burn rate and large port height. Otherwise, the coupled problem features perturbations in the layer much stronger than the decoupled analog. For lower burn rates and port diameters the core temperature perturbation peak moves closer to the burning interface and the disturbance magnitude tend to be maximum in the combustion sublayer. Temperature and mass fraction ﬂuctuations are larger than velocity analogs in the sublayer; velocity correlations in the sublayer are marginally dependent on the burning characteristics, while enthalpy correlations are enhanced in slow burning propellants, showing penetration in the ﬂame zone. These considerations point to the importance of the reaction structure and scales mainly on thermo-chemical ﬂuctuations. A major effect of the coupling on the burn rate perturbation is the modiﬁed pressure response coeﬃcient, expressing the ratio between the burn rate perturbation and the pressure perturbation at the edge of the sublayer as a function of the non-dimensional frequency. An increased level of interaction leads to a signiﬁcantly stronger pressure response. The ratio between coupled and decoupled responses increases from unity with an increase in ﬂuctuation frequency. With ﬁxed combustion characteristics, less unstable eigenmodes have a consistently weaker pressure response than more unstable modes. Disturbance penetration is linked to burn rate augmentation through an increased value of second order correlations. Therefore, variations in eigenfunction responses with the burning char-

887

acteristics can be compared to experimental trends based on burn rate augmentation data. The ﬂame temperature has a substantially weaker effect on the sublayer perturbations than the burn rate. Linear analysis results are in line with experimental observations indicating the burn rate as the main factor in determining erosive behavior. Analysis of other combustion characteristics, besides burn rate and ﬂame temperature, reveals that an increase in either pressure exponent or cold temperature sensitivity has a negative effect on the pressure response. These observations suggest that a propellant with chemical composition closer to the stoichiometric conditions (lower cold temperature sensitivity) and with coarser oxidizer particles (lower pressure exponent) support a stronger interaction between core ﬂow and sublayer. Thus such propellants are more likely to be affected by erosive augmentation. Future work will focus on non-linear analysis of the interaction and on particle scale effects. Acknowledgments The author would like to thank Dr. Thomas L. Jackson for reading this manuscript. This work was supported by the US Department of Energy through the University of California under subcontract B523819. References [1] M.K. Razdan, K.K. Kuo, in: K.K. Kuo, M. Summerﬁeld (Eds.), Fundamentals of Solid Propellant Combustion, in: Progress in Astronautics and Aeronautics, vol. 90, AIAA, 1984, pp. 515–598. [2] M.K. King, J. Propul. Power 9 (6) (1993) 785–805. [3] H.S. Mukunda, P.J. Paul, Combust. Flame 109 (1997) 224–236. [4] J.M. Most, P. Joulain, G. Lengelle, J.C. Godon, Combust. Flame 105 (1996) 202– 210. [5] R.A. Beddini, AIAA J. 18 (1) (1980) 1346–1353. [6] M.K. Razdan, K.K. Kuo, AIAA J. 20 (1) (1982) 122–128. [7] B.B. Stokes, R.O. Hessler, L.H. Caveny, Erosive burning of non-metallized composite propellants, in: Proceedings of the 13th JANNAF Combustion Meeting, vol. II, CPIA publication 281, Monterey, CA, September 1976, pp. 437–451. [8] H. Hasegawa, M. Hanzawa, S. Tokudome, M. Kohno, J. Propul. Power 22 (5) (2006) 975–983. [9] M. Chen, J. Buckmaster, T.L. Jackson, L. Massa, Proc. Combust. Inst. 29 (2002) 2923–2929. [10] G. Casalis, G. Avalon, J.P. Pineau, Phys. Fluids 10 (10) (1998) 2558–2568. [11] R. Dunlap, P.G. Willoughby, R.W. Hermsen, AIAA J. 12 (10) (1974) 1440– 1442. [12] M.K. King, Proc. Combust. Inst. 18 (1981) 207–216. [13] L. Massa, T.L. Jackson, J.D. Buckmaster, Combust. Theory Model. 11 (4) (2007) 553–568. [14] S. Gordon, B.J. McBride, Computer program for calculation of complex chemical equilibrium compositions, rocket performance, incident and reﬂected shocks, and Chapman–Jouguet detonations, SP 273, NASA, 1976. [15] T.L. Jackson, J. Buckmaster, AIAA J. 40 (6) (2002) 1122–1130. [16] R.A. Svehla, B.J. McBride, FORTRAN IV computer program for calculation of thermodynamic and transport properties of complex chemical systems, TN D7056, NASA, 1973. [17] B.J. McBride, S. Gordon, M.A. Reno, Coeﬃcients for calculating thermodynamic and transport properties of individual species, TM 4513, NASA, 1993. [18] J.N. Butler, R.S. Brokaw, J. Chem. Phys. 26 (6) (1957). [19] E.M. Abu-Irshaid, J. Majdalani, G. Casalis, Phys. Fluids 19 (024101) (2007) 1– 11. [20] J. Griffond, G. Casalis, Phys. Fluids 13 (6) (2001) 1635–1644. [21] T. Herbert, Annu. Rev. Fluid Mech. 29 (1997) 245–283. [22] J.J. Dongarra, B. Straughan, D.W. Walker, Appl. Numer. Math. 22 (1996) 399– 434. [23] S.V. Malik, A.P. Hooper, Phys. Fluids 19 (052102) (2007) 1–18. [24] M. Borri, P. Mantegazza, Meth. Appl. Mech. Eng. 12 (1977) 19–31. [25] C.B. Moler, G.W. Stewart, SIAM J. Numer. Anal. 10 (1973) 241–256. [26] A.I. Atwood, T.L. Boggs, P.O. Curran, T.P. Parr, D.M. Hanson-Parr, J. Propul. Power 15 (6) (1999) 740–747. [27] A.I. Atwood, T.L. Boggs, P.O. Curran, T.P. Parr, D.M. Hanson-Parr, C.F. Price, J. Wiknich, J. Propul. Power 15 (6) (1999) 740–747.

888

L. Massa / Combustion and Flame 156 (2009) 865–888

[28] R.L. Foster, R.R. Miller, The burn rate temperature sensitivity of aluminized and non-aluminized HTPB propellants, in: 1980 JANNAF Propulsion Meeting, vol. IV, CPIA#315, 1981, pp. 667–693. [29] M.K. King, Model for steady state combustion of unimodal composite solid propellant, AIAA Paper 1978-216, AIAA, 1978. [30] T. Marklund, A. Lake, ARS J. 30 (1960) 173–178. [31] L.A. Povinelli, M.F. Heidmann, C.E. Feiler, Experimental investigation of transverse-mode solid-propellant combustion instability in a vortex burner, TN

D-3708, NASA, 1966. [32] J. Buckmaster, T.L. Jackson, L. Massa, M. Ulrich, Proc. Combust. Inst. 30 (2005) 2079–2086. [33] M.K. King, J. Spacecraft Rockets 16 (1979) 154–162. [34] M.K. Razdan, K.K. Kuo, AIAA J. 18 (1980) 669–677. [35] R.J. Burick, J.R. Osborn, Erosive combustion of double-base solid rocket propellant, in: 4th ICRPG Combustion Conference, vol. II, CPIA publication 162, December 1967, pp. 57–69.

Spatial linear analysis of the flow in a solid rocket motor with burning walls

Spatial linear analysis of the flow in a solid rocket motor with burning walls

Recommend Documents