f 

FINAL REPORT 

for 

SRB COMBUSTION DYNAMICS ANALYSIS 

COMPUTER PROGRAM (CDA-1) 

NAG8 - 6 2 7 

... ". 

T. J. Chung and 0. Y. Park 

Department of Mechanical Engineering 

The University of Alabama in Huntsville 

Huntsville, AL 35899 

prepared for 

George C. Marshall Space Flight Center 

National Aeronautics and Space Administration 

Marshall Space Flight Center, AL 35812 

(NASA-CR-179302) ) S8E COIBUSTXOY DYNAMICS N88-15032 

AUALYSIS COMPUTER PBOGRAB (CDA-1) Final 

Report [Alabaaa Unit.b 130 p CSCL 21B Unclas 

63/25 0118982 

February 1988 
c n 

ABSTRACT 

A two-dimensional numerical model is developed for the 

unsteady oscillatory combustion of the solid propellant 

flame zone. 

frequency responses across the long flame, 

double-base propellants, are accommodated. 

is based on a premixed, laminar flame with a one-step 

overall chemical reaction and the Arrhenius law of 

Variations of pressure with low and high 

such as in the 

The formulation 

.. . 

decomposition for the gaseous phase with no condensed phase 

reaction. Numerical calculations are carried out using the 

Galerkin finite elements, with perturbations expanded to 

the zeroth, first, and second orders. The numerical 

results indicate that amplification of oscillatory motions 

does indeed prevail in high frequency regions. 

second order system, the trend is similar to the first 

For the 

order system for low frequencies, but instabilities may 

appear at frequencies lower than those of the first order 

system. 

system is that the admittance is extremely oscillatory 

between moderately high frequency ranges. 

The most significant effect of the second order 

i 
J 

TABLE OF CONTENTS 

ABST.C ............................................ TABLE OF CONTENTS.o.ooo.ooo.o~o ..................... LIST OF FIGURES.....* ............................... LIST OF TAB.S ...................................... LIST OF S.BOLS.= ................................... I . INTRODUCTION ................................... I1 . PRELIMINARIES FOR COMBUSTION WAVES ............. 2.1 Premixed Laminar Flame .................... 2.2 

... 

Fluid Dynamically Excited Oscillations .... 2.3 Propellant Characteristics ................ 

Page 

i 

ii 

iv 

vii 

viii 

1 

3 

3 

6 

6 

2.4 Acoustic Admittance 7 ....................... I11 . COMBUSTION DYN.ICS ............................ 10 

3.1 General ................................... 10 

3.2 Governing Equations ....................... 10 

3.2.1 Gas Phase .......................... 11 

3.2.2 Solid Phase ........................ 13 

3.2.3 Perturbation Expansions ............ 15 

3.3 Finite Element Applications ............... 27 

IV . DISCUSSIONS .................................... 31 

V . CONCLUSIONS .................................... 43 

REFERENCES .......................................... 45 

ii 
I I ” 

APPENDICES 

A. DERIVATION OF NON-DIMENSIONAL GOVERNING 

8 

EQUATIONS......ee.....eoooooooea..e..e....e...e 79 

B. PERTURBATION OF MODEL EQUATIONS................ 88 

C. DERIVATIONS OF BOUNDARY EQUATIONS.............. 93 

De DERIVATIONS OF FINITE ELEMENT EQUATIONS........ 112 

... 

iii 
Figure 

3.1 

4.1 

4.2 

4.3 

4.4a 

4.4b 

4.5 

4.6 

4.7 

4.8 

4.9a 

LIST OF FIGURES 

Page 

Domain of study for one-dimensional case... 47 

Finite element geometry. ................... 48 

Steady-state distributions of field 

variables in flame zone.................... 49 

First. arde; pressure distributions 

vs. frequency .............................. 50 

First order distributions of field 

variables at w = 1 (real parts)............ 51 

First order distributions of field 

variables at w = 1 (imaginary parts).. ..... 52 

First order density distributions 

vs. frequency .............................. 53 

First order temperature distributions 

vs. frequency .............................. 54 

First order species (fuel) distributions 

vs. frequency. ............................. 55 

First order velocity distributions 

vs. frequency. ............................. 56 

First order acoustic admittance and 

burning rate vs. frequency ................. 57 

iv 
4.9b 

4.9c 

4.10 

4.11 

4.12 

4.13 

4.14 

4.15 

4.16 

4.17a 

4.17b 

4.18 

Steady oscillatory mode: typical dependence 

of amplitude on frequency of pressure 

oscillation... ............................. Real parts of acoustic admittance and 

burning rate vs. frequency.. ............... Second order time-independent pressure 

distributions vs. frequency ................ Second-order time-independent distributions 

of field variables at w = l................ 

Second order time-independent density 

distributions vs. frequency. ............... 

distributions vs. frequency. ............... 

(fuel) distributions vs. frequency ......... 

distributions vs. frequency ................ 

Second order time-independent temperature 

Second order time-independent species 

Second order time-independent velocity 

Second order time-dependent pressure 

distributions vs. frequency ................ Second order time-dependent distributions 

of field variables at w = 1 (real parts)... 

Second order time-dependent distributions 

of field variables at w = 1 (imaginary 

parts) ..................................... Second order time-dependent density 

distributions vs. frequency ................ 

sa 

sa 

59 

60 

61 

62 

63 

64 

65 

66 

67 

68 

V 
4.19 Second order time-dependent temperature 

distributions vs. frequency ................ 4.20 Second order time-dependent species 

(fuel) distributions vs. frequency.. ....... 4.21 Second order time-dependent velocity 

distributions vs. frequency..... ........... 4.22 Second order pressure distributions 

vs. frequency... ........................... 4.23 Second-order distributions of field 

variables at w = l............. 

4.24 Second order density distributions 

vs. frequency. ............................. 4.25 Second order temperature distributions 

vs. frequency.... .......................... 4.26 Second order species (fuel) distributions 

vs. frequency .............................. 4.27 Second order velocity distributions 

vs. frequency........ ...................... 4.28 Offset burning rate and acoustic admittances 

of the second order time-dependent and 

time-independent systems ................... 

69 

70 

71 

72 

73 

74 

75 

76 

77 

78 

vi 
Table 

1 

LIST OF TABLES 

Page 

Typical value and range of parameters 

used in numerical calculation............. 32 

.. . .,. 

vii 
LIST OF SYMBOLS 

a speed of sound 

B 

cP 

C8 

Da 

E 

G 

h 

H 

I 

k 

k 

R 

L 

frequency factor for gas phase reaction 

specific heat at constant pressure for gas 

specific heat of solid 

Damkohler number 

gas phase activation energy, acoustic energy in 

Eq. (2.14) 

surface activation energy 

eigenvector 

source term in numerical formulation 

heat of combustion per unit mass 

heat of reaction 

order of perturbation 

thermal conductivity, constant in Eq. 

wave number (2.29) 

reaction rate constant 

characteristic flame length 

latent heat of vaporization, chamber length in Eqs. 

(3.11), 

(2625) 

m mass flux 

Mb Mach number 

n order%of chemical reaction 

P pressure 

viii 
Pr 

Q 

r 

R 

RC 

t 

T 

U 

ui 

UC 

V 

V 

W 

W 

X 

X 

Y 

Y 

Prandtl number 

dimensionless heat of gas phase reaction 

burning rate 

gas constant 

dimensionless distance from surface in Eq. (3.22) 

time 

temperature 

axial gas velocity 

velocity ve&or .. 

core velocity in Eq. (3.22) 

normal gas velocity 

total volume 

reaction rate 

molecular weight 

axial distance parallel to the surface 

solution vector, mole fraction 

normal distance from the surface 

mass fraction of fuel species 

oxidizer-fuel ratio, z-axis in cylindrical coordinate 

thermal diffusivity 

ratio of solid to gas density 

ratio of heat capacity, surface area 

temperature exponent in Eq. (3.7) 

specific heat ratio 

perturbation parameter 

= k,"Cp*/k*C, * 

ix 
eigenvalue 

Lagrange multiplier 

viscosity coefficient 

stoichiometric coefficient 

= k*/k,* 

density 

scalar potential in Eq. (2.26) , 

function . 

frequency -. 

interior volume 

.,. 

Subscrhts and SuDerscriDts 

dimensional quantity 

steady-state mean value 

spatial variable 

it order perturbation coefficient 

gas side of interface 

solid side of interface 

propellant cold side 

time-dependent variable 

fuel 

homogeneous solution 

vector quantity 

time-independent variable 

steady-state, oxidizer 

particular solution, product 

reference value 

interpolation 

X 
S solid phase 

arsr 

Y,zI 

t complex conjugate 

finite element global node number 

... 

xi 
CHAPTER I 

INTRODUCTION 

Determination of oscillatory pressure and velocity 

variations from their mean values in solid propellant 

combustion chambers has been the subject of study for many 

years. Because of computational difficulties, however, 

rigorous and accurate methods of solution still remain in 

various stages of development. 

are obscured in oomplicated computational strategies 

view of the fact that reliable experimental measurements 

are difficult to obtain, it is crucial to possess an 

analytical tool to accurately predict combustion dynamics 

and combustion efficiency in designing successful solid 

propellant rocket motors. 

Often details of physics 

In 

Earlier works on combustion dynamics include Hart and 

McClure [l], Denison and Baum [Z], Culick [3] and T'ien [4] 

among others. These studies are centered around onedimensional, linear oscillatory burning. Culick [5,6], 

Baum and Levine [7], and Yang et al. [8] subsequently 

studied the nonlinear growth and limiting amplitude of 

acoustic waves in a combustion chamber. Nonlinear behavior 

as characterized by higher order perturbation expansions of 

governing equations has been investigated by Flandro [9] 

and Chung and Kim [lo] in recent years. They assumed that 

the pressure varies sinusoidally with time, the pressure 

1 
2 

amplitudes being a function of position. 

In the present approach, the pressure is assumed to 

vary in the unstable state. This allows the presence of a 

long flame such as in double-base solid propellants, in 

which high pressures and high frequencies may be 

accommodated. Galerkin finite elements [ll] are utilized 

to model numerically the geometries of burning surfaces and 

the flame thickness. The nonsteady governing equations are 

linearized by means of the zeroth, first, and second order 

perturbation expansions. The boundary conditions, .- . 

including the burning surfaces and flame edges, are imposed 

by means of the Lagrange multipliers. 

In the simple example problems demonstrated in this 

paper, the gaseous flame is assumed to follow the Arrhenius 

law with one-step forward chemical reactions. No condensed 

phase reaction is included in the formulation. The 

calculated results confirm the prediction reported by other 

investigators in the literature for the low frequency 

region. 

frequency region. 

oscillatory motions does prevail in the high frequency 

region. 

Extended studies are then carried out for the high 

It has been found that amplification of 

Discussions are presented in Section 4. 

The present study does not include turbulent boundary 

layers which may play an important role in erosive burning 

and combustion instability; this is the subject of the 

current stbdy. 
CHAPTER I1 

PRELIMINARIES FOR COMBUSTION WAVES 

2.1 Premixed Laminar Flame 

Combustion is the chemical process which gives rise to 

the conversion of reactants into products. 

combustion, chemical reactions are coupled with heat and 

mass transfer. 

rapid chemical reactions occurring exothermically in single 

or composite phase. 

chain sequence composed of four basic processes: 

initiation, branching, propagation, and termination of the 

radicals. 

in the form 

During 

Energy and species are transported with 

Most chemical reactions occur in a 

A chemical reaction to combustion can be written 

N N 

where Y, denotes the chemical species a and v,',~ and 

represent the stoichiometric coefficients of species Y, for 

the reactants and for the reaction products in the mth 

reaction, respectively. The stoichiometric coefficients 

and van;m must satisfy elemental material balances. 

The phenomenological principle of mass action states that 

the rate of reaction is proportional to the product of the 

concentration of the reactants. Thus, the reaction rate is 

3 
in which X0 is the mole fraction and k, is the reaction 

rate constant, with m indicating the reaction order. 

order of a reaction is determined by summing the exponents 

of the reactant mass fraction that express the reaction 

rate. 

empirically given by the Arrhenius law 

The 

The constant XA for a single step reaction is 

with R being the universal gas constant, E, the activation 

energy, and A,,, the frequency factor for the reaction step. 

According to the theory related to the collison process and 

potential energy barrier of the reaction, E, is the minimum 

barrier height which reactants must acquire before reaction 

process. The factor exp(-E,/RT) represents the fraction of 

collisions between reactant molecules in which products can 

be formed. Generally, E, is constant, but A, depends on T 

such that 

Am = B,T6m (2.4) 

where B, and 6, are constant. 

gives 

Combining those relations 
5 

(2 5) 

in which the detailed reaction kinetics are neglected. 

For a single step irreversible chemical reaction, the 

chemical process may be written in the form 

v'F + v'O + v"P 

f 0 P. 

with F, 0, and P-being the species representing the fuel, 

oxidizer, and product. In this case, the overall reaction 

rate is given by 

It is noted that the exponents r and s are not necessarily 

related to the stoichiometric coefficients vf' and yoB. In 

many cases, actual values of B, 6, r, s, and E are 

empirically determined. Based on this model, the local 

heat release rate Qf can be determined from Eq. (2.7) and 

the heat of combustion per unit mass of fuel h,, 

where Qf represents heat addition to the system. 

Eqs. (2.7) and (2.8) can be used in conjunction with the 

species and energy transport equations, respectively. 

Thus, 
6 

2.2 Fluid Pvna micallv Excited Oscillations 

The contribution of flow excitation is not usually 

distinguishable from that due to combustion. 

there are three rather distinct examples that indicate flow 

excitation by means other than combustion: (1) oscillations 

related to vortex shedding and transport from the segmenttransition regions of the motor cavity in large segmented 

motors at low frequencies [13]; (2) excitation by the 

interaction between transverse flow oscillation and axial 

mean flow [14]; and (3) excitation by turbulence [15]. The 

contribution of flow excitation to oscillations is usually 

small compared to the combustion oscillations, primarily 

because flow excitation results from the kinetic energy 

flux in the mean flow, while combustion oscillation results 

from the much larger thermal energy flux in the combustion 

zone. However, flow excitation may play a role in 

initiation of combustion-driven oscillations and 

correspondingly affect stability limits. 

However, 

... ”. 

2.3 ProDellant Characteristics 

Experience shows that combustion instability is highly 

dependent on propellant characteristics. It is clear that 

the effects can be due to either changes in combustion 

dynamics of the propellant or changes in damping. 

the stability is a property of the entire combustor and 

cannot generally be attributed to the propellant alone, the 

Although 
7 

variables involved in the response function are related to 

the characteristics of propellants. The growth constant in 

the response function shows that a high energy, fast 

burning propellant in a charge configuration with a large 

burning surface area will tend to be unstable unless the 

real part of the response function is small or the damping 

effect is large. It has been shown that the value of the 

response function for energetic propellants is high enough 

at some frequencies to cause motor instabilities. These 

propellants include ammonium perchlorate, double-base 

propellant, aluminized propellants, and combustion of 

propellants with other energetic oxidizers. 

2.4 Acoustic Admittance 

The acoustic wave equation for an adiabatic 

compressible inviscid flow is given as 

where a is the speed of sound and # is a scalar potential 

defined as ui = #,i. 

The corresponding acoustic pressure PI is 

a# 

p' 3 -p - 

at 

(2.10) 

Taking a partial derivative of Eq. (2.10) leads to 
8 

(2.11) 

The general solution of Eq. (2.9) in the one-dimensional 

case is 

(2.12) 

in which A and B-are complex quantities, k is the wave 

number, and w is-the frequency. If the positive going 

travelling wave is considered, it follows that 

= Cei (kx-ut 1 

Substituting Eq. (2.13) into Eq. (2.11) yields 

From Eqs. (2.13) and (2.14), we have 

P' 

U 

Pa -3 

(2.13) 

(2.14) 

(2.15) 

The quantity pa is called the specific impedance in 

acoustics and the reciprocal of the normalized specific 

impedance is denoted as the admittance in combustion 

instability. Thus, the acoustic admittance Y' is given as 

U 

y' 3.- - 

P' 

(2.16a) 
9 

Dividing by the steady-state value gives the nondimensional form of Y', 

Po 

u, P' 

y=-- (2.16b) 

where the subscript o indicates the mean flow. 

The so-called response function F is defined as 

F = Re(m'/m,)/(S'/Pb) (2.17) 

where m is the mass flow rate and Re is the real part of 

the complex quantity. 

response functions are representative of the combustion 

instability. In what follows, however, we shall use the 

acoustic admittance as a measure of combustion wave 

instabilities. 

Both acoustic admittances and 
CHAPTER I11 

COMBUSTION DYNAMICS 

3.1 General 

A numerical model is developed for the unsteady state 

combustion of a solid propellant subject to acoustic 

pressure oscillations in both low and high frequency 

ranges. 

second order and then solved by using the finite element 

method. After reproducing the analytical model given by 

Friedly and Petersen [16,17], the numerical results are 

verified for the one-dimensional case by comparison with 

the analytical solution. Thus, the work is extended for 

the multi-dimensional case where the critical assumptions 

required in the analytical method may be removed. 

result, the investigation of variable distributions in the 

flame zone is made possible and, consequently, the acoustic 

admittance is obtained to determine the combustion 

instability in a more realistic manner. 

The governing equations are perturbed to the 

As a 

Discussions of governing equations and finite element 

solutions will be presented in the following sections. 

3.2 Governina Emations 

The governing equations for a premixed laminar flame 

are derived under the following assumptions: (1) the Mach 

10 
11 

number is small, (2) body forces are negligible, 

(3) radiative heat transfer is negligible, (4) the Soret 

and Dufour effects are negligible, (5) diffusion caused by 

the pressure gradient is negligible, (6) the gas mixture is 

ideal and thermal coefficients are constant, and (7) the 

Lewis number is unity for the explicit relation between 

concentration and temperature. 

of the coordinate system for a one-dimensional model is 

shown in Fig. 3.1%. 

A schematic representation 

.,. 

3.2.1 Gas Phase 

For the simplicity of the combustion modeling of solid 

propellants, the gaseous flame is assumed to be multidimensional, premixed laminar, and calorically perfect, and 

a one-step forward chemical reaction occurs. The 

combustion of a solid propellant is approximated by the 

Arrhenius law. Thus, the conservation laws for the multicomponent reactive system in the gas-phase are represented 

in the non-dimensional forms as follows (see Appendix A for 

details) : 

Continuitv 
12 

Moment Urn 

I 3 I 1 

Pli - Pr ui,jj + - uj,ij = 0 

au 1 

+ pui,juj + - 

2 Pat rMb 

Enerav 

aT r-iap 

.. Species .. . 

ay 

State 

P = pT (3.5) 

where the commas denote partial derivatives, the repeated 

indices imply summing, Pr is the Prandtl number, and Y 

represents the fuel mass fraction. 

of three species equations (fuel, oxidizer, and product) is 

taken into account from the simple chemical reaction model 

[4]. 

render the above equations dimensionless: 

Note that only one out 

The following characteristic parameters are used to 

P = P*/PO * , P = P*/P,* , T = T*/T, * 

* ui = UI*/V, , t = t*V,*/j* , xi = Xi*/j* 
13 

, h = h*/Cp*To* I w = w ** a /Po * v, *2 * 

Mb 31 vo*/ao 

in which j* is the flame thickness given by u*/v,*, with U* 

being the thermal diffusivity, M, represents the Mach 

number; k*, the thermal conductivity; a,*, the speed of 

sound; h*, the combustion heat release; and w*, the 

reaction rate, whose dimensionless form is 

with z denoting the oxidizer-fuel ratio; n, the order of 

chemical reaction; E, the activation energy given by E = 

E*/RT,*; and B, the rate constant. 

represents dimensional quantities and subscript zero gives 

the mean value at the flame edge. 

The superscript * 

3.2.2 Solid Phase 

The solid propellant is assumed to be homogeneous, 

with no condensed phase chemical reaction. 

dimensionless form of the heat transfer equation in the 

solid phase is 

The 

(3.8a) 

Furthermore, assuming that the heat transfer is onedimensional gives 
. 14 

(3.8b) 

where 

13 - PS*/PO * , 1 = ks*Cp*/k*Cs* 

r = r*/r* with r* = po*vo*/Ps * 

Here, r denotes the burning rate at the solid surface and 

subscript s represents the solid phase. 

process of the solid propellant at the surface is assumed 

to follow an Arrhenius law; thus, 

The decomposition .-I 

r = ex{-Es[ - J] (3 9) 

in which E, is the dimensionless surface activation energy, 

E, = Es*/RTo*, and Fs is the mean temperature at the 

surface. 

determined by the dimensionless mass and energy balances 

across the interface, such that 

The solid-gas interface boundary conditions are 

(3.10) 

with 5: = k*/ks* and L = (H+* - H-*)/Cp*T0 * Here, L is the 

latent heat of vaporization of the propellant and H" 

denotes the enthalpy changes. 

represent the gas and solid side at the interface, 

The subscripts + and - 
respectively. 

15 

3.2.3 Perturbation ExDansions 

When a small pressure disturbance occurs in the 

combustion chamber, every field variable will be disturbed 

from its steady-state value and can be expressed in the 

form, 

F = F'O' + QF"' + Q F (*) + ... (3.11) 

... .. . where F = (p, ui, T, Y, P} and represents the 

perturbation parameter. 

parentheses indicate the perturbation order. 

assuming sinusoidal fluctuation of pressure with time 

The superscripts in the 

Furthermore, 

renders the variables in a different form: 

I I = 1,2,... (3.12) 

It is important to recognize that the sources in the 

second order consist of inhomogeneous terms that describe 

the nonlinearities in terms of the product of two first 

order variables. 

leads to the separation of each inhomogeneous term into a 

time-independent term and a term oscillating at the 

frequency of the second harmonic. 

only the real parts of F, 

quantities, the product indicates, in reality, that 

Considering the physical quantity of F' ' 

Because of the fact that 

and F,"' are the physical 
16 

Re(F,“’) (F,“’) = Re[$,(l’eiut ]Re [ P2 ( ’ e ‘ 1 

= Re[(F,, + iF,)(cos ut + i sin ut)]Re[(F2, + 

+ iF2) (cos ut + i sin ut)] 

1 1 

1 

2 

-- (F,,F2 + F,F,,)sin ut 

The right-hand s&de of the above equation equals the real 

part of (1/2)(F,(1’tF2‘1’ + F1(1)F2(1)), where the dagger 

represents the complex conjugate. Thus, the product of 

F,“’ and F2(’) may be replaced by 

(3.13) 

and each product of two first order terms can be replaced 

by a similar expression. 

constant term and a term oscillating at the frequency of 

the second harmonic, the dependent variable F2“’ may be 

rewritten as 

Since each nonlinearity yields a 

From the expression of Eq. (3.14), those equations for the 

second order system can be separated into two individual 
equations, one for the time-independent equations and the 

other for the time- dependent equations. 

coefficients of the time- independent equation set are 

exactly the same as those of the first order set, except 

the inhomogeneous terms which consist of each product of 

two first order variables. For a higher order, the same 

argument is applicable. 

All the 

Substituting Eqs. (3.11) ‘(3.14) into Eqs. (3.1) -( 3.5) 

and rearranging separately in the order of perturbation 

yield the final form of the governing equations 

corresponding to each order. 

Steadv-State Governina Emations 

The one-dimensional steady-state governing equations 

are given as follows: 

Continuitv 

Enerav 

Species 

(3.15) 

(3.16) 

(3.17) 
State 

18 

(3.18) 

Note that the uniform pressure is retained throughout the 

flame thickness (P' 

both the dependent variables, p"' and T"', are equal to 

unity at the flame edge, which is far from the origin on 

the scale of R*. . From Eqs. (3.15) and (3.18), we have 

= l), and from Eqs. (3.6) and (3.11), 

(3.19) 

and from Eqs. (3.16) and (3.17), Y'O' can be expressed as 

Therefore, 

(3.20) 

(3.21) 

where i = 0 and the second order chemical reaction 

is assumed. 

(n = 2) 

Now, the equations are expressed in terms of 

temperature in the steady-state; thus, only the solution of 

Eq. (3.16) is required. Flandro [9] suggests a simple 

analytical model of the mean flow field to facilitate 

multi-dimensional analysis in the higher order system. 

This model is given in the form 
19 

(3.22) 

Here, u, describes the flow speed along the local 

streamline and R, is a dimensionless distance from the 

solid surface. 

in the steady-state solution: 

At the flame edge, 

The following boundary conditions are used 

or 

dT‘ ’ 

=o 

dY 

At the solid phase, Eq. (3.8) gives 

(3.23) 

(3.24a) 

(3.24b) 

(3.25) 

where T, is the propellant cold side temperature. The 

continuous temperature condition at the interface requires 

that 

(3.26) 

The matching condition across the interface can be obtained 

by substituting Eq. (3.25) into Eq. (3.10) , which yields 
20 

(3.27) 

In solving Eq. (3.16), note that a correct rate constant B 

in Eq. (3.21) has to be determined such that the system 

satisfies the boundary conditions at the flame edge as well 

as at the interface. 

Hiqher Order Governins Eauations 

The higher order governing equations can be expressed 

in a single form because only the source terms are 

~.* 

-. 

different from each other. 

the right-hand side ,of the equations in terms of GS, the 

governing equations are represented as follows: 

Continuity 

Presenting the source terms on 

(3.28) 

Momentum 

Enerav 

(3.29) 

(3.30) 
21 

(3.31) 

(3.32) 

and the reaction-rate' is given by 

Here, G, is given as follows: for the first order system 

(I = l), G = 0; for the second order system (I = 2), 
22 

(3.34) 

Note that, including the pressure coupling, the velocity 

coupling is significant in the source terms (Appendix B). 

At the solid phase, Eq. (3.8) gives -I 

with 

Linearizing Eq. (3.9) results in 

where 

(3.36) 

Substituting Eq. (3.36) into Eq. (3.35) and solving 

analytically yield 
23 

in which 

c3c2 1 1 

iId[ 5 

G9 = [-XI + ;] + ; [X2c4 + cg - 2) - c,L 

and 

. .. 1 -. 

XI = - [l + 

25 

(1 + i41wS{)"2] 

-2 

=1 = E,& 

c2;( 1) 2 

cq = - 

2w213*5 

Equation (3.37) gives a boundary condition at the surface; 

other conditions are as follows: 

At the flame edge, 

;(I) 0 (3.38) 

Assumption of an isentropic flow near the flame edge gives 
24 

the temperature conditions. 

temperature gradient at the flame edge is almost zero, 

condition can be expressed in the form 

Noting that the steady-state 

this 

with 

Since the flux of each reactant species is always a 

fixed fraction of the total flux, the fuel mass flux 

fraction mf can be derived from Eq. (3.4) as 

For a one-dimensional expression, mf is given as 

(3.40) 

(3.41) 

where m is the mass flux equal to pv, and assumption of the 

constant burning rate at an instant has been used. 

surf ace, 

At the 

'A dY'O' 

+ GI1 (3.42) m' I) 

1 d$"' + A (I) = $(I) - 

dY m' 0 I2 dY m' 0 I m f 

in which 
25 

Using the relationship m "(I) = pS&(I) yields 

(3.43) 

with 

The normal velocity component is derived from the mass 

balance condition at the interface such that 

(3.44) 

where 

Moreover, the parallel velocity component may be obtained 

using the Taylor series expansion about the surface where 

the no-slip condition must be valid. Thus, 

(3.45) 
26 

with 

Note that G, - GI4 are valid only for I = 2. For higher 

order systems, Eqs. (3.37)-(3.39) and Eqs. (3.43)-(3.45) 

are used as boundary conditions to solve Eqs. (3.28)- 

(3.32). It is necessary to have more conditions for the 

density at both sides for better solutions, and the 

pressure fluctuation 'has to be forced at the flame edge. 

In the case of the second order time-independent system, 

care must be taken to use boundary equations (3.35) and 

(3.39). See Appendix C for details. 
3.3 Finite aement Amli cations 

The finite element method (FEM) is an approximate 

numerical procedure for solving partial differential 

equations of boundary and initial value problems. 

basic idea resides in the application of variational 

principles or equivalent concepts such as weighted residual 

The 

methods. 

to the finite element analysis is the Galerkin method, 

which is one of khe iethods of weighted residuals. 

It is well known that the most general approach 

In this 

scheme, test functions are the same as trial functions, 

known more popularly as interpolation functions. 

Typically, any scalar variable X is expressed as a linear 

combination of interpolation function BB and the nodal 

values of Xg such that 

(3.46) 

where B denotes the global nodes. The Galerkin finite 

element equations result from the inner product of the 

residual of a governing equation and the test function such 

as 

(3.47) 

Integrating this equation by parts gives rise to the 

simultaneous algebraic equations. 

is employed for the integrations of the spatial finite 

The Gaussian quadrature 
28 

element domain. Consequently, the global finite element 

equations are of the form 

[iIuA,, + B,,]X, = G, (3.48) 

in which i denotes the imaginary part of the complex 

number, X 

TB, Y,, P,], and G, is the inhomogeneous source term valid 

for the second order. 

represents the solution vector, X, = [p,, uRi, R 

The algebraic finite element equations are modified 

for the boundary conditions. 

boundary conditions is achieved by adopting the Lagrange 

multipliers method. 

boundary conditions have to be provided in the global form 

of the nodal values. 

in the form of the boundary matrix equations, 

Implementation of complicated 

The solution variables included in the 

The boundary conditions are expressed 

q,#, = b, (3.49) 

with Y = 1,2,...,m, with m being the number of equations. 

Thus, the modified algebraic finite element equations have 

the final form: 

(3.50) 

In combination with the Lagrange multipliers method, 

it is also possible to impose additional Dirichlet boundary 

conditions. Implementation of Dirichlet conditions is 
29 

. 

achieved by zeroing out the rows and columns of A af3 

Baa, placing 1 at the diagonal corresponding to the 

boundary nodes with all the boundary terms subtracted from 

G,. 

other equations. 

and 

The influence of this process is propagated to all 

It is seen that Eq. (3.50) lends itself to a standard 

eigenvalue problem, the solution of which determines the 

existence of excited natural frequencies and modes. 

Eq. (3.48) with zero ‘setting of G, gives the eigenvalue o 

which is the natural frequency of the given system. 

Using 

The solution of Eq. (3.50) at a given frequency is 

obtained by imposing the Dirichlet condition of the 

pressure at the flame edge. The finite element formulation 

contains two different kinds of test function used to 

represent the volume and surface integrals in the domain. 

The total number of field variables would be reduced by one 

if the density or pressure were eliminated using the 

perfect gas law. However, the stability of the matrix 

system is doubtful. 

Before calculations, the following is expected: since 

the Mach number is generally very small, the coefficient of 

the pressure term in Eq. (3.2) dominates the system unless 

the frequency considered is large enough such that other 

coefficients containing the frequency factor become 

comparable in magnitude. Therefore, in lower frequencies, 

the pressure gradient has to be negligible, thus resulting 
30 

in constant pressure. On the contrary, the gradient will 

become significant in higher frequencies; this results in 

pressure variance. In the latter case, severe pressure 

changes will occur if unbounded at the solid surface. Note 

that care must be exercised in expanding Eq. (3.9) due to 

the appearance of the exponential growth effect. 

For the steady-state case, as mentioned earlier, we 

calculate the eigenvalue B in which necessary initial 

conditions are satisfied. 

later, the eigenvalue is referred from the results of 

reference [4] for this study. The Galerkin finite element 

formulation of the governing equations in each system is 

described in Appendix D. 

However, as will be discussed 
CHAPTER IV 

DISCUSSION 

In the steady-state, the flame zone is governed by 

Eqs. (3.15) through (3.18). When an acoustic pressure wave 

at a certain frequency is imposed on the flame zone, the 

system will react according to Eqs. (3.28) through (3.32). 

The present computations on homogeneous solid propellants 

are performed for an adiabatic flame with a second order 

-. 

chemical reaction [Eq. (3.7)]. Instead of solving Eqs. 

(3.15) through (3.18), distributions of the field variables 

in the steady-state are calculated from the temperature 

given in Eq. (3.16). The coefficients used in this 

computation are chosen in such a way that the results may 

be compared with those in [4,9]. In view of this, the 

following parameters are utilized: z = 1, 6 = 0, T, = 0.15, 

T, = 0.35, I = = 1, E = 10, E, = 4, L = 0.15, R = 1000, 

- 

Y = 1.2, M, = 0.003, u, = 1.0, R, = 5.0, Pr = 1.0, and h = 

1.3. Representative dimensional parameters corresponding 

to the above dimensionless values are given in Table 1, and 

these are based on a single-base propellant [12]. 

To demonstrate the validity of the theory presented 

above, a two-dimensional rectangular geometry, shown in 

Fig. 4.1, is analyzed, in which a burning surface and flame 

edge can be established as boundaries. 

was chosen to resemble a one-dimensional behavior so that 

This configuration 

31 
32 

Table 1 Typical value and range of parameters 

Typical 

Parameter Value 

Physical 

Range Variable 

Typical 

Value 

1000 

1 

0.15 

0.35 

1.0 

10 

4 

- 

- 

- 

- 

1 

1 

2 

0.5 

- 

1.4 

0.15 

250-1000 

- 

- ... -. 

- 

4-15 

2-10 

0.5-2 

0.4-0.85 

10- 3-102 

0.05-0.3 

cal/gmole 

cal/gmole 

cal/gOk 

cal/gOk 

cal/cmOksec 

cal/cm ksec 

g/cm2sec 

atm 

cal/gmole 

1.5 

1.5 x 10'~ 

300 

700 

2000 

40 x io3 

16 x io3 

0.33 

0.33 

5 x 10'~ 

5 x 10'~ 

0.4 

9.5 

- 

- 

- 

- 

2.8 x io3 
33 

the results may be compared 

the literature [4,9]. 

The computations begin 

frequencies of the system. 

directly with those reported in 

with finding the natural 

The resonance frequencies of 

the system can be obtained by examining the acoustic 

admittance or response function at each natural frequency. 

As mentioned earlier, the natural frequencies of the system 

are given by the homogeneous solution of the total matrix 

equation (3.50), neglecting the boundary conditions. The 

frequency range of interest in this computation is shown to 

be between w = and w = lo2, corresponding to 

approximately 80 Hz and 8 MHz, respectively. Thus, sixteen 

different frequencies in this range are chosen to evaluate 

the present study, and these are w = 5 x 

-. 

5 x io-2, 5 x io-2, io-l, 5 x io-l, 1, 5, io, 20, 30, 

40, 50, and 100. 

Figure 4.2 shows a typical steady-state distribution 

of the field variables and the reaction rate along the 

flame thickness. In general, these data may be obtained by 

solving the energy equation (3.16) together with the 

reaction rate in Eq. (3.21). Ideally, a fourth order 

Runge-Kutta scheme may be used for this purpose. The 

result will be used as the basis of further calculations 

for combustion instability in the unsteady-state. 

Distributions of the field variables versus frequency 

for the first order system are shown in Figs. 4.3-4.8. 

These results are obtained by imposing an acoustic pressure 
amplitude of unity at the flame edge as the Dirichlet 

condition. Conventionally, the thickness of the burning 

zone has been assumed to be negligibly small compared with 

the wavelength of the acoustic oscillation in the 

intermediate frequency range. Thus, the uniform pressure 

is approximated throughout the domain of study and is 

assumed to vary with time only. On the contrary, since the 

oscillatory pressure is regarded as a spatially 

nonhomogeneous time-dependent source in the present study, 

it is possible to investigate the response of a specific 

... 

solid propellant at significantly high frequencies and to 

find the response even in the long flame of a double-base 

propellant. Figure 4.3 shows variations of the pressure 

distribution versus the acoustic frequencies. It is seen 

that the amplitude remains constant for w < 10. Random 

variations of these amplitudes occur for higher 

frequencies. Here, the imaginary parts representing the 

phase shift are zero. 

Figure 4.4 demonstrates distributions of various 

component fluctuations of the first order system at w = 1. 

A pressure wave with a certain amplitude striking the solid 

surface will give rise to a response in the burning rate. 

Generally, the response of the burning rate is nonlinear 

and has a complete Fourier series form that may not be 

expressed by a single frequency component. In Fig. 4.4, 

however, the pressure is uniform in the flame zone, 

although significant variations may occur at higher 

34 
35 

frequencies. 

The temperature amplification at the surface changes 

the burning rate in Eq. (3.36), while the velocity at the 

flame edge represents the acoustic admittance for pressure 

coupling, whose magnitude and sign indicate the instability 

of the system. 

implying that the acoustic wavelength is larger than the 

flame thickness in this case, the imaginary part of the 

velocity approaches a.constant slope at the edge by the 

assumption of an isentropic condition at the flame edge. 

It is interesting to note that the trend of variations of 

temperature and velocity appears to be similar, which is 

inversely proportional to the variations of species and 

density as expected. 

of T'ien [4] and Flandro [9]. 

Since the pressure remains constant, 

.. 

These results are comparable to those 

Density distributions for various frequencies are 

given in Fig. 4.5. At w < 0.5, changes of the 

distributions versus frequency are negligible; however, 

near w = 1 a significant reduction of the amplitude occurs, 

and then for w > 1 the amplitudes increase moderately. The 

distributions in the lower frequency region are very 

similar to those in the steady-state, from which the quasisteady assumption may be deduced. Note that the reverse 

peak moves in the flame zone, which implies a change of the 

reaction zone. Close to the surface, the effect of preheat 

in the solid phase increases, resulting in a higher burning 

rate. Positive amplitudes at the interface in each 
36 

frequency are caused by the rise in pressure and/or the 

increase of the burning rate. 

It is interesting to see that the amplitudes decrease 

while the frequency increases up to the unity and then 

increase along the frequency, from which the flame is 

expected to react sensitively at special frequencies. The 

negative amplitude is shown most notably at w = 1. 

may result from the net blowing effect, and the velocity at 

that point should be positive (see Fig. 4.8). 

This 

-. 

The normal velocity distributions are shown in Fig. 

4.8. 

negative slopes, the only exception being at o = 1. 

latter contains most of the distributions in the lower 

frequency region. 

interface, relating to an increase of the density. 

Although the pressure remains constant in lower 

frequencies, the velocity varies severely. At w > 20, the 

variation is more significant since from that frequency the 

pressure varies in the flame zone. At w = 0.01, note that 

the slope is positive even though it is a lower frequency. 

This result implies instability of the system for the 

quasi-steady case. 

interface results from adjusting between the density 

increase and the higher gasification rate. As indicated in 

Fig. 4.4, the slopes of the imaginary parts near the flame 

edge are constant. 

The profile may be classified into positive and 

The 

Negative amplitudes appear mostly at the 

At w = 1, a positive amplitude near the 

Rearranging the real part of the velocity at the flame 
37 

edge gives the distribution of the acoustic admittance, 

whose magnitude and sign indicate the amplifications or 

damping ability of the flame subject to the acoustic 

disturbance. Figure 4.9 shows the curve-smoothed acoustic 

admittance and burning rate from Eq. (3.36) at each 

frequency . At w = 1, there is a definite increase in the burning 

rate. The burning rate increases by increasing the heat 

transfer effect due to closer movement of the reaction zone 

to the surface. 

effects such as structural change in the solid phase and 

change of the chemical kinetics in the flame. These 

effects cause the feedback to the solid phase, resulting in 

a change of the burning rate. These effects, however, are 

not considered. 

-. 

There are actually several secondary 

Figure 9 also reveals a resonance in the condensed 

phase near w = 0.01, indicating that the system is 

unstable. 

and previous laboratory measurements which have shown that 

the response consists generally of a single peak in the 

range of frequency for which the quasi-steady approximation 

appears to be accurate. 

This verifies the result of Denison and Baum [Z] 

Some negative peaks exist at the 

other frequencies, implying that the resonance in the gas 

phase tends to damp the oscillatory motion. 

frequencies, the system seems to be unstable. 

The real part of the burning rate shows a trend 

At most higher 

similar to the acoustic admittance at the quasi-steady 
38 

region, although the magnitude is significantly different. 

But these trends differ from each other at the higher 

region. 

mode oscillation as the pressure varies in the flame zone. 

It must be emphasized that the admittance function alone is 

not sufficient to describe the stability of the system 

unless the velocity coupling is concerned. 

Over w = 100, the profile tends to have a second 

Temperature distributions are given in Fig. 4.6. 

Since the temperature is related to the density by Eq. 

(3.5), it can be analyzed easily by comparing the results 

with those in Fig. 4.5. Accordingly, the distribution of 

the temperature shows the profile opposed to that of the 

density. 

reaction region for w > 1. 

temperature changes occur at the interface as imposed by 

the boundary conditions. Negative amplitude near the flame 

edge may be caused by the stagnation phenomena in which the 

density increases. 

imaginary parts of the temperature are given as reciprocals 

of those of the density. 

.,. -. 

Temperature rises are predominant along the 

Note that no significant 

Because of the pressure unity, the 

The fuel consumption at a given frequency is shown in 

Fig. 4.7. 

those for temperature distributions in that negative 

maximum occurs along the reaction region. Note that 

linearly diminishing the fuel near the flame edge affects 

the temperature changes slowly (Fig. 4.6). 

Trends for species distributions are opposite to 

The results of the second order perturbation are shown 
39 

in Figs. 4.10-4.28. As previously indicated, each variable 

of the second order response to acoustic disturbance has 

two components: one time-dependent component that 

oscillates at twice the fundamental frequency and one that 

is time-independent. The latter represents a shift in the 

mean value, thus causing a shift of the mean burning rate. 

It should be noted that the variations for the second order 

system may be characterized by the sum of those two parts, 

with the factor ei 

indicating the case of time-independence. 

ranging between -1 and 1 and zero ..I 

Figures 4.10-4.15 show the distributions in the second 

order time-independent case. The calculations are also 

conducted by imposing the pressure of unity at the flame 

edge as the Dirichlet condition. General trends show that 

shifts in the mean values are evident. 

In Fig. 4.10, the pressure varies from o = 5, which is 

half of that in the first order system. 

frequency region, the oscillatory motion is significant, 

which may affect the velocity distribution. Figure 4.11 

shows all of the variables for w = 1, indicating that 

shifts in mean values may have the affect of damping; 

however, the opposite may be true for higher frequencies as 

demonstrated in Figs. 4.12-4.15. The trends are roughly 

similar to those of the first order illustrated in Fig. 

4.4a, but they have different amplitudes because of the 

nonlinearities in the higher order system. In Figs. 4.12- 

4.15, the distributions of field variables in lower 

For the higher 
frequencies appear to vary negligibly along the flame 

thickness, while the opposite is true for w > 1. In 

particular, since the pressure varies at w > 5, the 

velocity changes significantly at those frequencies, 

implying that the system is unstable at higher frequencies. 

As mentioned earlier, a shift of the mean burning rate is 

the most significant effect in the time-independent system. 

This 

case 

result will be discussed in Fig. 4.28. 

The distributions in the second order time-dependent 

are shown in Figs. 4.16-4.21. Here, the pressure also 

varies from w = 5 and oscillates for higher frequencies. 

Note that the amplitude variations are smaller than those 

in the time-independent case. The trends seem to reduce 

the total effect of the second order system. 

Figures 4.22-4.27 show the variations of field 

variables in the second order system by summing timeindependent and time-dependent effects, with the factor 

= 1 for time-dependent effects. Note that amplitude =iZot 

variations of each variable in the flame zone are 

significant at w = 1, which is the same result obtained for 

the case of the first order system. Finally, Fig. 4.28 

reveals that the shift in burning rate versus frequency is 

consistently negative, as asserted by other investigators. 

The offset is relatively small in the quasi-steady region, 

but increases with oscillatory motion along the 

frequencies. 

similar to those of the first order system, except that 

The variations of the acoustic admittance are 
41 

I. 

oscillatory deviations are more predominant at higher 

frequencies in both time-independent and time-dependent 

cases. 

Parameter studies are conducted for the first order 

and are summarized as follows. 

ratio I3 affects the variables shifted slightly to the 

negative direction, keeping the distribution profiles 

Decreasing the density 

constant. Changing the latent heat of solid L exerts a 

negligible effect on fhe variables, but a very small value 

of L shifts the system toward instability. Increasing the 

surface activation energy E, or the gas phase activation 

energy E reduces the magnitude of the variables, keeping 

the same profiles. 

-. 

Changes of the rate constant and 

viscosity effect coefficient strongly affect the system, 

such that every aspect discussed herein will change. 

The present study could be extended to the multidimensional case by introducing the appropriate axial mean 

flow field. It is well recognized that fluctuation of the 

gas velocity parallel to the propellant surface affects the 

burning rate in terms of velocity coupling; therefore, this 

quantity must be considered together with pressure coupling 

for a satisfactory measure of stability. 

A simple calculation has been accomplished using the 

artificial axial flow velocity in Eq. (3.22). However, 

difficulties of the boundary conditions could not be 

eliminated. A test run shows that the existence of a small 

amount of the axial flow reduces differences between the 
variable amplitudes in each frequency considered, thus 

leading the system toward stability. 

... ^. 
43 

CHAPTER V 

CONCLUSION 

A multi-dimensional numerical model for the premixed 

flame acoustic instability is proposed and solved using the 

finite element method. 

perturbed to the second order and formulated with the 

The governing equations are 

Galerkin finite elements. 

on the validity of published theories of solid propellant 

combustion instability at the lower frequency region where 

the uniform pressure assumption is valid. 

are made on the higher frequency region and some new 

results are obtained. 

The results have direct bearing -< 

Extended studies 

Under the restricted boundary conditions, the 

following conclusions, based on numerical computations, are 

reached : 

The stability characteristics for the low frequency 

range in the first order system have been verified to 

be the same as those reported in the literature. For 

example, the acoustic admittance is controlled by the 

burning rate, negative for low frequencies, whereas 

the opposite is true for high frequencies. 

Instabilities are likely to occur at high frequencies 

(w > 10). 

For the second order system, the trend is similar to 

the first order system for low frequencies, but 
44 

instabilities may appear at frequencies lower than 

those of the first order system. 

(4) The most significant effect of the second order system 

is that the admittance is extremely oscillatory 

between w - 1 and w = 10. 

... ”. 
45 

REFERENCES 

1. Hart, R.W. and McClure, F.T., ttCombustion Instability: 

Acoustic Interaction with a Burning Propellant 

Surfacett, J. Chem. Phvs., Vol. 30, Sept. 1959, pp. 

1501-1514. 

2. Denison, M.R. and Baum, E., "A Simplified Model of 

Unstable Burning in Solid Propellantstv, J., Vol. 

31, Aug. 1961, pp. 1112-1122. 

3. Culick, F.E;C., "A Review of Calculations for Unstable 

Burning of a Solid Propellanttt, AIM J., Vol. 6, Dec. 

1968, pp. 224102254. 

4. Ttien, J.S., tlOscillatory Burning of Solid Propellants 

Including Gas Phase Time Lagtt, Combustion Science and 

Technoloav, Vol. 5, 1972, pp. 47-54. 

5. Culick, F.E.C., ItNonlinear Behavior of Acoustic Waves 

in Combustion Chambers-IIt, Acta Astronautica, Vol. 3, 

1976, pp. 715-734. 

6. Culick, F.E.C., ItNonlinear Behavior of Acoustic Waves 

in Combustion Chambers-IItt, Acta Astronautica, Vol. 3, 

1976, pp. 735-757. 

7. Baum, J.D., Levine, J.N., and Lovine, R.L., IIPulseTriggered Instability in solid Rocket Motors", AIAA 

J., Vol. 22, No. 10, Oct. 1984, pp. 1413-1419. 

8. Yang, V., Kim, S.I. , and Culick, F.E.C. , "Third-order 

Nonlinear Acoustic Oscillations in Combustion 

Chambers, Part I", AIM paper 87-1873, June 1987. 

9. Flandro, G.A., I'Nonlinear Time-Dependent Combustion of 

a Solid Rocket Propellant", 19th JANNAF Combustion 

Meeting, Oct. 1982. 

10. Chung, T.J. and Kim, P.K., "Unsteady Response of 

Burning Surface in Solid Propellant Combustiontt, AIAA 

paper 85-0234, Jan. 1985. 

11. Chung, T.J., Finite Element Analysis in Fluid 

Dynamics, McGraw-Hill Book Co., 1978. 
46 

12. Cantrell, R.H., Hart, R.W.! and McClure, F.T., "Linear 

Acoustic Gains and Losses in Solid Propellant Rocket 

Motors", AIM (1., Vol. 2, No. 6, June 1964, pp. 1100- 

1105 . 13. Brown, R.S., Dunlap, R., Young, S.W., and Waugh, R.C., 

"Vortex Shedding as a Source of Acoustic Energy in 

Segmented Rockets", of SDacecraft and Rockets, Vol. 

18, July-August 1981, pp. 312-319. 

14. Flandro, G.A., "Rotating Flows in Acoustically 

Unstable Rocket Motors", AIM J., Vol. 2. July 1964, 

pp. 1305-1306. 

15. Hegde, U.G. and Strahle, W.C., "The Generation of 

Sound by Turbulence Inside Rocket Motor Cavities", 

AIAA paper 84-0287, .* 1984. 

Combustion Parameters on Instability in Solid 

Propellant Motors, Part I: Development of Model and 

Linear Analysis", AIM J., Vol. 4, No. 9, Sept. 1966, 

-. 

16. Friedly, J.C. and Petersen, E.E., "Influence of 

pp. 1605-1610. 

17. Friedly, J.C. and Petersen, E.E., IIInfluence of 

Combustion Parameters on Instability in Solid 

Propellant Motors, Part 11: Nonlinear Analysis", AIAA 

J., VOl. 4, No. 11, NOV. 1966, pp. 1932-1937. 
47 

* 

TC 

I I 

I I 

- -Y -OD 0 

CONDENSED $A&EE(aUS PHASE 

Fig. 3.1 Domain of study for one-dimensional case. 
48 

DIRICHLET B.C. SOLID SURFACE (Ts) 

SOLID PHASE 

* DEEP IN SOLID(Tc) 

-X 

Fig. 4.1 Finite element geometry. 
49 

? 

P) 

? 

(Y 

w 

3 

b 

I4 

n 

c 

4 

? 

4 

? 

0 

0.0 1.0 2.0 3.0 

FLAME LENGTH 

Fig. 4.2 Steady-state distributions of field 

variables in flame zone. Parameters 

used in the calculations are given 

in Table 1, including Pr=l and Mb= 

0.003. Reference values for non-dimensionalization are chosen from the 

flame edge. ( p :density,T:temperature, 

Y:species of fuel, v:velocity, P:pressure, and w:reaction rate). 
50 

Real 

Imaginary 

w=10 

~=O;O001-5 4 

-. 

0-0.0001-100 

I I 

0.0 I .o 2.0 3.0 

FLA;\lE LEXGTH 

Fig. 4.3 First order pressure distributions vs. 

frequency(Dirich1et condition is imposed 

at the flame edge). 
51 

\ T 

\ 

I I 

0.0 1 .o 2.0 

FLAME LENGTH 

3.0 

Fig. 4.4a First order distributions of field 

variables at w=l(real parts). 
52 

T 

P 

0.0 1.0 2.0 

FLAME LENGTH 

3.0 

Fig. 4.4b First order distributions of field 

variables atw= l(imaginary parts), 
Real 

----- Imaginary 

100 

0.0 1.0 2.0 

FLAME LENGTH 3.0 

Fig. 4.5 First order density distributions vs. 

frequency. 
54 

0 

r; 

0 

CIO . 

w 

Real 

Imaginary ------ 

0.0 1.0 2.0 

FLANE LENGTH 

3.0 

Fig. 4.6 First order temperature distributions 

vs . frequency. 
55 

Real ------ Imaginary 

- w=o.l\ 10, 

0.0 1 .o 2.0 

FLAME LENGTH 

3.0 

Fig. 4.7 First order species(fue1) distributions 

vs . frequency . 
56 

r; 

0 

NReal 

------ Imaginary 

Y 

w 

w 

n 

I 0=1.0 I 

I 

0 

r) 

I 

10 I I 

0.0 1.0 2.0 3.0 

FLAME LENGTH 

Fig. 4.8 First order velocity distributions 

vs. frequency. 
57 

Acoustic admittance 

----- Burning rate 

FREQUENCY 

Fig. 4.9a First order acoustic admittance and 

burning rate vs. frequency. Positive 

peak near u= 0.01 resembles the results 

of Denison and Baum(2) and T'ien(4). 
58 

0 

QD 

N 

X 

0 

u) 

In 

47n.r 

X 

h 

u 

c 

s 

I 

CUDs 

k 

rlE 

Y 

U 

0 
59 

0.0 1 .0 2.0 3.0 

FLAME LENGTH 

Fig. 4.10 Second order time-independent pressure 

distributions vs. frequency. 
u. 

2 

? 

Jn 

0 

0 

e 

0 

T 

0 

. 

z 

I 

P 4 

... 

\- 

0.0 1 .0 2.0 3.0 

FLAME LENGTH 

Fig. 4.11 Second order time-independent distributions of field variables at w = 1. 
61 

w = 1 (reduced by 1/2 

0.0 1.0 2.0 3.0 

FLAME LENGTH 

Fig. 4.12 Second order time-independent density 

distributions vs. frequency. 
62 

0.0 1.0 2.0 3.0 

FLAME LENGTH 

Fig. 4.13 Second order time-independent temperature 

distributions vs. frequency. 
63 

.- .-. 

I I 

0.0 1 .o 2.0 3.0 

FLAME LENGTH 

Fig. 4.14 Second order time-independent species 

(fuel) distributions vs. frequency. 
64 

u. 

s 

9 

In 

0 

0 

. 

0 

Y 

9 

z 

I 

by 1/21 

0.01 

0.0 1.0 2.0 

FLAME LENGTH 

3.0 

Fig. 4.15 Second order time-independent velocity 

distributions vs. frequency. 
I. 5 65 

real 

imaginary 

0.0 1 .o 2.0 3 .o 

FLAME LENGTH 

Fig. 4.16 Second order time-dependent pressure 

distributions vs. frequency. 
66 

I 1 

0.0 1.0 2.0 

FLAME LENGTH 

3 .0 

Fig. 4.17a Second order time-dependent distributions of field variables at w = 1 

(real parts). 
67 

? 

3 

9 

! 

9 

0 

d 

0 

e 

m 

s 

a 

9 

Y 

9 

5? 

I 

0 

e 

2 

I 

0 

0 

N 

I 

f " 

./ 

P 

0.0 1.0 2.0 

FLAME LENGTH 

3.0 

Fig. 4.17b Second order time-dependent distributions of field variables at w = 1 

(imaginary parts). 
68 

? 

5' 

9 

real 

---- imaginary 

a r' I 1 

0.0 1.0 2.0 

FLAME LENGTH 

3 .O 

Fig. 4.18 Second order time-dependent density 

distributions vs. frequency. 
69 

real 

---- imaginary 

w=l 

10 

0.0 1.0 2.0 3 .O 

FLAME LENGTH 

Fig. 4.19 Second order time-dependent temperature 

distributions vs. frequency. 
A 

GI 

w 

3 

Er 

VI 

w 

H 

u 

w 

PI 

VI 

Er 

0 

w 

3 

I3 

- 

n 

H 

Fig. 

L! 

0 

N 

9 

f3 

s 

s! 

0 

m' 

reduced by 1/2) 

0.0 1.0 2.0 3 -0 

FLAME LENGTH 

4.20 Second order time-dependent species 

(fuel) distributions vs. frequency. 
71 

’? 

0 

N 

---- imaginary 9 

2 

9 

0 

4 

0 

Y; 

0 

? ---eO.Ol’/ 

10 

(reduced by 1/1C 

q , , T 

0.0 1 .a 2.0 

FLAME LENGTH 

3 .O 

Fig. 4.21 Second order time-dependent velocity 

distributions vs. frequency. 
72 

I 0.1 0) .= 

0.0 1.0 2.0 3.0 

FLAME LENGTH 

Fig. 4.22 Second order pressure distributions 

vs. frequency. Both time-dependent 

and time-independent coefficients 

are added for denoting the maximum 

deviation from mean value. 
. 73 

0.0 1 .0 2.0 

FLAME LENGTH 

Fig. 4.23 Second order distributions of field 

variables at u = 1. 
74 

u . 

5 

? 

In 

9 

0 

? 

In 

I 

9 

z 

I 

0.0 

'1 (reduced by 1/21 

1.0 2.0 

FLAME LENGTH 

3.0 

Fig. 4.24 Second order density distributions 

vs . frequency. 
I .. 

w = 1.0 

1.0 2.0 3.0 

FLAME LENGTH 

Fig. 4.25 Second order temperature distributions vs. frequency. 
I. .. 76 

n 

I4 

w 

3 

L4 

v) 

w 

u 

w 

a 

v) 

L 

0 

w 

3 

b 

I4 

v 

l-4 

n 

I+ 

u. 

s! 

9 In 

9 

0 

9 

Y 

9 

0, 

I 

0.0 

.,. -. 

1.0 2.0 

FLAME LENGTH 

3.0 

Fig. 4.26 Second order species (fuel) distributions vs. frequency. 
c! 

5 

? 

In 

? 

0 

? 

Y 

9 

s! 

I 

0.0 

(reduced by 

1/10) 

1.0 2.0 

FLAME LENGTH 

3.0 

Fig. 4.27 Second order velocity distributions 

vs . frequency . 
I I. 

78 

r. 

-. 

PURNINF RATE 

--_ TIFE-DEPENDENT -. -.- TIME-I~DEPENDENT 

I 

F 

Fig. 4.'28 Offset burning rate and acoustic 

admittances of the second order time 

-dependent and time-independent systems. Negative value of the burning 

rate at the frequency region verifies 

the past investigations. 
APPENDIX A 

DERIVATION OF NON-DIMENSIONAL GOVERNING EQUATIONS 

1. Continuitv Emation 

The dimensional continuity equation is of the form 

ap* 

- + (P*Ui*),i = 0 

at* 

Substituting Eq. (3.6) yields 

PO*V,* ap 1 

- + - (po*vo*pui),i = 0 

j* at j* 

Dividing (A. 2) by po *vo */R * yields 

which is equivalent to Eq. (3.1). 

2. Momentum Eauation 

Neglecting the body forces, the dimensional governing 

equation is given by 

* 1, au 

at* I* 3 

* * * 

P* - + p*~i,j~j + ~,i - U* ui,jj + - uj,ij 

Substituting Eq. (3.6) gives 

PO * 

R* 

*2 

*- VO 

pUi,jUj + - P, i 

v,*~ au, 

PO* - P- + Po R* at R* 
80 

Dividing the equation by po*~o*2/~* results in 

I 3 I 1 

ui,jj + - uj,ij = 0 

au PO * U* 

at Po *vo * PO*VO*R* 

P - + PUi,jUj + P,i - 

Since 

PO*VO? 

*2 

Mb2 = [ >] = * 

YPO 

where M, is the Mach number at the reference point, then 

PO * 1 

For the last term, 

a* k* 

with a* being the thermal diffusivity. Thus, 

u* u* Po *CD*V0 * u*c,* = Pr 3 - - 

Po *vo *R * p *v, *k* k* 

Substituting Eqs. (A.7) and (A.8) yields 

au 1 f 1 \ 

P __+ PUi,jUj + - P,~ - Pr(ui,jj + - ujjijJ = o 

at YM, 3 
, . .. 81 

which is equivalent to Eq. (3.2). 

3. Enerw Eauation 

Neglecting the radiation effect, the dimensional governing 

equation is given by 

aT* ap* * n 

at* at a=l 

p*cp* - + p*Cp*T,i*ui* - - - k*T, i i - C w,*h,* = 0 (A.10) 

Substituting Eq. (3.6) yields 

~A 

1 3T vo*T0* po*vO* ap 

j* at R* j* at 

po*Cp*Vo*To* - P- + po*cp* PT,iUi - - 

* Po*vo*2 

Cp*T, C w,h, = 0 

TO * 

T,ii - - k* - 

l** a* a 

Dividing Eq. (A.ll) by po*Cp*~o*To*/~* gives 

" 

P- + pT,iUi - -- T,ii - C w,h, = O 

at pa*To*Cp* at U 

(A. 11) 

(A. 12) 

since p0* = P~*RT, * I 

R CV * 7-1 

--= l--= PO * 

Po*TlJ*Cp* CP * cP* Y 

where R = Cp * - Cv* and Y = Cp*/Cv* are used. Therefore, the 

non-dimensional governing equation is of the form 

aT r-iap 

P- + pT,iui - -- 

at Y at 

T,ii - Wh = O (A. 13) 

in which the last term is obtained from Eq. (A.30) Equation 
82 

(A.13) is equivalent to Eq. (3.3) 

4. Species Eau ation 

The dimensional governing equation is of the form 

a Ya k* 

P* - + P*Y,,iUi* - - ~,,ii - w,* = o , a = 1,2, ..., N 

(A. 14) cP at* 

where N denotes the total number of species considered. 

Substituting Eq. (3.6) yields 

VO * k* 1 

*- pYa,iUi - - - Ya, i i 

VO * a Ya 

R* CpR R** PO* - P- + Po R* at 

* *2 

Po vo - w, = 0 

a* 

Dividing Eq. (A.15) by po*v0*/R* gives the following 

aytz .. k* V,*P* u IP- + PYa,iUi - Ya,ii - w, = 0 

U* ** 

at cp R PO*%* 

(A. 15) 

(A. 16) 

Since 

U* k* 

VO 

1" a-= * 

Po vo*cp* * 

The coefficient of the third and the last terms are to be unity. 

Therefore, 

(A. 17) 
where Ya is defined by 

83 

(A. 18) 

The three mass fraction equations for the fuel, oxidizer, and 

product in Eq. (A.17) are similar to one another, allowing 

treatment of the fuel equation alone with the following relation 

Ya Yf 

- - - = constant- 

"a Yf 4 -. 

(A. 19) 

Therefore, Eq. (A.17) for the fuel only is enough to describe the 

whole system and is given as 

This is equivalent to Eq. (3.4). 

5. Reaction Rate of Fuel Species 

For a one-step forward chemical reaction, 

"f'[f] + YO'[O] + YP't[P] 

(A. 20) 

(A.21) 

where Y denotes the stoichiometric coefficient, and f,o,p 

represent fuel, oxidizer, and product, respectively. Using the 

Arrhenius law, the reaction rate of the fuel can be expressed by 

(A.22) 

in which B* is the frequency factor for the gas phase; R, the gas 
84 

constant: Wjt the molecular weight of the species; E*, the 

activiation energy: and 8 and n are constants. For the oxidizer 

and products, the reaction rate can be written as follows: 

vu" - VU 

W,* = w* I a = 0,p 

The the last term of Eq. (A.12) becomes 

** - 1 w,*hUo = w h 

(A.23) 

(A.24) 

where 

,* 

... 

c Wa(v,' - v,")ha0 

h* = 

Wf Vf 

is the heat of combustion per unit mass of fuel consumed. 

the definition in Eq. (A.18), 

Using 

and the following relation in Eq. (A.19) 

Ya y f 

- - - = constant 

"u Vf 

(A.25) 

(A. 26) 

the reaction rate term can be non-dimensionalized. 

Since both the fuel and oxidizer must vanish at the flame 

YO Yf ---=o (A. 27) 
Therefore, 

85 

YO 

Yo = - y f 

Vf 

From Eq. (A.25), vf = -1; thus, 

Yo = 'V0Yf 

or 

The non-dimensional reaction rate can then be written as 

where 

B*k*To*8po*" 

B= 

N -. 

po *2 V, *2 C, TI Wj"j 

j=1 

and this is eqyivalent to Eq. (3.7). 

(A. 28) 

(A.29) 

(A. 30) 
86 

6. Heat Tran sfer Emation in Solid Phase 

The dimensional governing equation is given as 

(A. 31) 

with c, * and k,* being the specific heat and the thermal 

conductivity of the propellant, respectively, and y = 0 is 

attached to the surface. Substituting Eq. (3.6) yields 

T0*2 

Vo*T,, iui - k,* - Ts,ii = 0 

vo* aT, TO * 

ps*c,*T0* - - .4 ps*'Ls* - 

j* at R* R*2 (A. 32) 

Dividing Eq. (A.32) by po*cs*To*vo*/~* results in 

Since R* = k*/po*vo*Cp*, the last term can be written as 

k,*To* k,*C,* 

Therefore, Eq. (A.33) becomes 

(A.33) 

(A.34) 

(A. 35) 

In the one-dimensional case, since po*vo* = ,os*:* if r is defined 

by r = r*/r* (actually r* is the same as ui* in Eq. (A.31)), then 

a2T, 

D- +r--(-- -0 

8% aTS 

at aY aY2 

(A.36) 
87 

which is equivalent to Eq. (3.8). 

7. Solid-Gas Tnterface Boundarv Condition 

The surface energy balance across the interface can be 

expressed as follows: 

aT* 

-k*( - 1 + p*v*H+* = -k, 

\ aY* I, 

(A. 37) 

with subscripts + and - wing used to indicate the gas and solid 

side of the interface, respectively, and H* denotes the enthalpy. 

Using Eq. (3.6) and defining [ = k/k,* and L = (H+* - H,*)/Cp*To* 

yield 

-. 

*-* TO 

R* 

+ ps r Cp*To*rL = -k,* - 

or 

*:* = * * and R* = k*/po*Vo*Cp * , then since ps Po vo 

or 

(A. 38) 

(A. 39) 

which is equivalent to Eq. (3.10). 
88 

APPENDIX B 

PERTURBATION OF MODEL EQUATIONS 

1. Perturbation of Reaction Rata 

The dimensionless reaction rate equation is given by Eq. 

(3.7). After taking 1 = 0 and n = 2, we have 

L 

1 E 1 E 

+€- T,' + - 6 '1 - E ]2T112][ 1+ c2 - (T2I - T112) T( 0 I T( 0 I 2 T( 0 I 
I. 89 

E Tl'](-2Tl' + 2P1') + 2c2 - E Y,'T,' + r2(2Y2' +- 

T'o) T'o I * 

-. 

I 1 E E 

+ Y,'2) + - 

2 

(T2' - T1l2) 

where F,' = $'lI/F'o' and F,' = $'2)/F'o) are used. Therefore, 

1 +'l) 2+'1I 291I 

+ + 

'O) [ [ 5-21 T" y'0, p' 0 1 

=w $1) 

(B.3) 

E E 

+ + 

p'0 I 

+ 3[1 - - 

T'oI 
90 

Since 

and 

we have another form of the reaction rate such that 

If we substitute T'O) for l/p(O) and (l/h)(l - T'O)) for Y'O), 

then 
91 

(B.9) 

Combining Eqs. (B.6) and (B.7) yields Eq. (3.33). 

2. Perturbation of Burninq Rate 

From Eq. (3.9), the burning rate r can be perturbed in the 

... 

f om 

(B.lO) 

Therefore, 

(B.ll) 
92 

where T, = ?, at y = 0 is used, 

(B.12) 

(B. 13) 

Combining Eqs. (B.12). and (B.13) yields Eq. (3.36). 
93 

APPENDIX C 

DERIVATIONS OF BOUNDARY CONDITIONS 

1. Temperature Condition & Interface 

(a) First Order Condition 

The first order solid phase energy equation (3.35) is 

rewritten as 

at aTy a2t 

iwat + r(O) - + ;f 1) -5-=o 

. aY aY 

-. 

aY 

where T denotes '?s ( for convenience. Knowing that r( = 1 and 

aT,(O)/ay = (i/5) (T, - T,)eY/l from Eq. (3.251, we have 

in which 

The general solution of Eq. (C.2) has the form 

with subscripts h and p denoting homogeneous and particular 

solution, respectively. The characteristic equation is then 

given by 

(C-i 2) 

SX' - x - iwR = o 
Thus, 

94 

which gives the homogeneous solution such as 

where th + 0 at deep-in golid (y + -a) is considered. 

The particular solution tP is of the form 

tP = ~,e~/f 

Since atplay = (i/l)tp and a2Tp/ay2 = (i/12)~p, we have 

Therefore, from Eqs. (C.3), (C.6), (C.7), and (C.8), 7 is given 

by 

From the faet that T = T, at y = 0, 

c2 

A, = T~ + 

id3 (C.10) 

which yields 
95 

ius 

and 

At y = 0, 

(C. 11) 

(C.12) 

(C. 13) 

Substituting Eq. (C.13) into Eq. (3.10) yields 

(C.14) 

where k'" = r(o)(~S/Tc2)+sr1) = cl+, is used. Therefore, 

(C.15) 
96 

(b) Second Ordel; Condition 

The second order solid phase energy equation is rewritten as 

or 

With the same argument of'the first order solution, the 

characteristic equation is given by 

1x2 - x - 1208 = 0 

of which solution is 

(1 + JX) 1 

1, = - 

21 

1 I 

(C.16) 

(C.17) 

(C.18) 

Therefore, the homogeneous solution is 

Th B,e XY 2 (C.19) 

Substituting Eq. (C.ll) into Eq. (C.16) yields the particular 

solution to be of the form 

(C.20) 

Since 
97 

XIB,eXiY + 

aY 

and 

I 

a27, 

-0 X 2BleX 1 y 

aY 

we have, from Eq. (C.16), the following equations: 

For the coefficient of the eXiY term, 

i2wAB1 5X12Bl = 0 

where 

+ X1B1 + F (C.21) 

F 

For the ey/l 

5 

the coefficient of term, 

c2 ;(1)2 

0 (C.22) 

Thus , 

F 

Bl 

B2 

(C.23) 

and 
98 e1 

(C.24) 

Since t = t, at y = 0, then 

which yields .A 

-. 

and 

Substituting Eq. (C.27) into Eq. (3.10) yields 
.. 99 

(C.28) 

Since 

(C.29) 
100 

where 

c2;r 1) 2 

c4 = - 

2wW5 

and 

1 

G, = c3c2 [-X2 + ;] + - 

i2~13k E 

are used. Combining Eqs. (C 15) and (C.30) yields Eq. (3.37). 
101 I. 

Note that for the time dependent condition, the term 9,' 

to be substituted by half of its value as indicated in Eq. 

has 

(3.13) 

If the second order time independent solution is required, 

then the boundary equation is different from the previous one and 

it is given as follows: From Eq. (3.35), 

The characteristic equation in this case is given by 

X(l - IX) = 0 

Therefore, the homogeneous solution is of the form 

and the particular solution is 

= D2eX1Y + D,yeY/' ?P 

((2.31) 

(C.32) 

(C.33) 

(C.34) 

Substituting Eq. (C.34) into Eq. (C.31) and investigating the 

coefficients gives the following relationships, i.e., 

F 

$(1)2, 

"(2) 2 

D, = r c2 - 

id{ 

(C.35) 

Now, the solution t is of the form 
102 

Since t = t, at y = 0 and T = 0 at y = -a, then 

Do = 0 

Finally, 

T = [$,(2 

F 

(+ e 

sX,2 - x, 

(C.36) 

(C.37) 

‘1Y 

(C.38) 

and 

At y = 0, 

(C.40) 

and then from Eq. (3.10) 
103 

I 

Therefore, 

1 

((2.41) 

... -. 

F 1 1 c2c3 

a+, 2) 

\+ I-- 1 -- 

-1 

ClC2 

- c3L] [ ft + - t + 4 (C.42) 

Note that in actual calculation the term GS ( 

substituted by (1/2)+$ ' ( as indicated in Eq. (3.13) . has to be 

s 
.* 

2. Temnerature Conditions & Flame Edae 

The temperature at the flame edge can be obtained by 

assuming that the flow very close to the flame is isentropic 

[21]). 

pressure yields 

Writing this expression in terms of temperature and 

P 

104 

(see 

Therefore, 

4 

-. 

(Y-l)/Y T P 

or 

Perturbing both sides yields 

(C.43) 

(C.44) 

(C.45) 

Since ln(l + X) = x - (1/2)~~ + (1/3)X3 - ... for I XI < 1, Eq. 

(C.45) can be rewritten as follows: 

1 

+ o(€~) 1 (C.46) 
105 

Thus, for the first order term, we have 

and for the second order term 

(C.47) 

(C.48) 

Since DS/Dt = 0, whe/re S denotes the entropy, which means 

-. 

that the entropy is conserved for each fluid particle convected 

away from the flame, Eq. (C.47) is of the form 

(C.49) 

Assuming that a+(l)/ay >> a+(l’/ax and $(I’ is constant at y = R 

(flame edge), then 

Substituting Eq. (C.47) into Eq. (C.48) shows 

(C.50) 

(C.51) 

and using the same argument about the entropy conservation, 
106 

Combining Eqs. (C.50) and (C.52) yields 

or 

which is equivalent to Eq. (3.39). 

3. Fuel Sizjecies Condition & Jnterface 

Since the mass flux fraction of the fuel and oxidizer is - 

fixed by the propellant composition, the species boundary 

condition at the gas-solid interface can be derived from Eq. 

(3.4). After neglecting the unsteady and reaction terms, the 

governing equation yields 

PUiY/i - Y/ii = 0 (C.54) 

For one-dimensional consideration, Eq. (C.54) can be integrated 

with respect to y and the results show 

(C.55) 

where m is the mass flux and mf is the fuel mass flux fraction. 

Perturbing Eq. (C. 55) yields 
107 

Therefore, 

((2.56) 

(C.57) 

Since Aft 1) = = 0, and using Eqs. (B.12) - (B.13) , we have 

(C.58) 
108 

Knowing that p,r(O’ = p0v(O) = 1, we have 

(C.60) 

(C.61) 

Combining the two eqttatigns yields Eq. (3.43). 

4. Normax Velocity Condition & Interface 

The boundary condition of the normal velocity component is 

decided using the mass balance condition at the interface. For 

the first order mass balance, we have 

Therefore, 

where 

(C.62) 

(C.63) 
109 

are used. For the second order condition, 

ps;12) (0+2) + "P1);Il) + ;;(2)v(o) P 

or 

(C.64) 

Combining Eqs. (C.63) and (C.64) yields Eq. (3.44). 

5. Parallel Velocitv Conditions & Interface 

Since the parallel component of the velocity has to be zero 

at the actual solid surface in order to satisfy the no-slip 

condition, we can derive the higher order condition using the 

Taylor series expansion about the origin such that 
au+ ( O 1 a2u+(0) 

aY 2 ay2 

ay + - ay2 + ... (0) + 

u+ = u+ 

110 

((2.65) 

Differentiating Eq. (C.65) with respect to t and multiplying p 

yields 

aut(o) 1 a 2~+ (0) aut 

at aY 2 aY 

P-= PV+ + - PV+ ay + ... 

For the first order condition, 

aY I+ 

or 

id3 aY It 

For the second order condition, 

1 aY I+ 

4w aY2 I+ 

or 

(C.66) 

(C.67) 
111 

aW2 ay2 

+ 

Combining Eqs. (C.67) and (C.68) yields Eq. (3.45). 

((2.68) 

... 4 
112 

APPENDIX D 

DERIVATIONS OF FINITE ELEMENT EQUATIONS 

1. Zeroth Order Analysis 

The formulation of the zeroth order governing 

equations (3.16) is 

(01 (ZAUB + ZBUB)Tb = ZG, 

4 

-. 

where 

P 

P 

and 

2. Piaher Order Analysis 

The finite element analogs for the higher order system 

are summarized below. 

Continuity 
(1) 

+ [FFuB + FHUB]vB = FG, 

where 

Momentum Eau ation. tX-comDonent) 
114 

where 

1 

XFUB = - 

3 

#a, i@fi, j dS1 , (j = 1) 

Pr s, 

XIuB = - Pr #u,i#g,j dS1 , (j = 2) 

XJUB = - J' klg, i dS1 

3 

1 

%* 51 
115 

Here, all i = 1. 

Momentum Euuation. (Y-comDonent) 

(I) (I) 

[YA,BI Pg + [YB,B + yc,&3 + [ (iIu)YDaB + YEaB 

(I) (I) 

+ YFaB + YHiB + YIaB]vB + YJ,& = YG, (D.4) 4 ”. 

where 

YA,B = xA,B I YBaB = XDaB 

YCaB = XFaB I YDaB = XBaB 

YEaB = XCaB I YFaB = XHaB 

YH,* = XEaB I YIaB = XIaB 

YJUB = XJaB f YG, = XG, 

Here, all i = 2. 
116 

... 

"51 
117 

SDecies Eauation 

where 
118 I