Numerical solution technique for transient, two-dimensional combustion with multi-step kinetics

COMPUTER METHODS IN APPLIED MECHANICS AND ENGINEERING 83 (1990) 9-31 NORTH-HOLLAND

NUMERICAL SOLUTION TECHNIQUE FOR TRANSIENT, TWODIMENSIONAL COMBUSTION WITH MULTI-STEP KINETICS T.M. KIEHNE Defense Advanced Research Projects Agency, U.S.A

D.E. WILSON, R.D. MATI'HEWS University of Texas, Austin, TX 78712, U.S.A. Received 21 April 1989 Revised manuscript received 2 November 1989 For numerical combustion problems, computational time and storage requirements scale approximately with the square of the number of chemical species. Therefore, the most difficult cases that can be solved using elementary chemical kinetics have been one-dimensional transient and two-dimensional steady formulations. This limitation and the established shortcomings of global reactions have encouraged recent development of simplified but accurate multi-step kinetics schemes. The present investigation uses one of these simplified reaction mechanisms to determine if solution of a multidimensional, transient, reacting flow problem is possible without resort to overly simplified, global kinetics. A laminar, two-dimensional, transient, viscous flow model has been developed which incorporates a chemical kinetics mechanism for propane consisting of four reversible reactions and only seven species. The model is used to investigate the initiation and propagation of a flame subjected to heat transfer, wall motion and shear flow. The mathematical model simulates the physical and chemical processes occurring near the centerline of a thin..gap formed by two cold, solid surfaces as one surface moves first toward and then away from the other. The governing conservation equations for mass, moraentum, energy and chemical species are simplified using the thin-gap approximation. The problem is transformed to eliminate the continuity equation by using a compressible stream function formulation. The resulting mixed order system of nonlinear, boundary value problems is solved by a semi-discrete numerical procedure which employs a finite element collocation method and adaptive grid technique. Both non-reacting and reacting cases are solved to demonstrate the effect of the shear flow and cold surfaces on the flow field and on the developing laminar flame. The results reveal a complex flow with numerous substructures caused by the interaction of the developing flame with cold fluid trapped in the 0oundary layer adjacent to the stationary wall. The reacting flow results represent the first calculation of a laminar, reacting flow in a moving wall geometry using nonequilibrium, multi-step chemistry while simultaneously resolving the momentum and thermal boundary layers.

1. Introduction The use of detailed, elementary chemical reaction mechanisms is the established m e t h o d for accurately representing the chemistry occuring in flames. However, for hydrocarbon fuels such mechanisms usually include 100 or more reversible, interdependent individual reaction steps and 25 or more chemical species. Because computer time and storage requirements scale roughly with the square of the n u m b e r of species considered [1], the most sophisticated ~a~_-7~9~,n,~,$03.50 © 1990- Elsevier Science Publishers B.V. (North-Holland)

10

T.M. Kiehne et al., Numerical solution technique for combustion with multi-step kinetics

problems that generally can be modeled using detailed chemistry are one-dimensional transient or two-dimensional steady formulations. Because most practical combustion systems are not encompassed by either one of these formulations, models of such systems have necessarily relied upon use of global, simplified kinetics schemes to model the chemistry. Simplified reaction mechanisms and their limitations are discussed in the review by Westbrook and Dryer [1]. The objective of the various types of simplified chemical mechanisms is to reproduce the primary features of the detailed kinetics while using a manageably small number of reaction steps involving a limited number of chemical species. The major advantage of overall kinetics modeling is the tremendous reduction in the number of ordinary differential equations necessary to predict the evolution of the primary reactants and products as well as the energy release as a function of time. This results directly in a reduction of the computational problems associated with the coupling of the chemistry and the other physical processes. The major disadvantage is the loss of detailed information due to reduction in the number of species included in the model. Typically, such schemes consist of a single irreversible reaction, but in some instances two or more steps are used. One reason for the general inadequacy of all of these simplifed reaction mechanisms is their failure to fully account for the sequential nature of hydrocarbon oxidation in that they do not include reactions representing the formation and oxidation of intermediate hydrocarbons. Inclusion of the intermediate hydrocarbons is necessary to accurately represent the delay in energy release and wall interactions of a flame. For these reasons, there has been a recent effort to develop simplified but accurate multi-step kinetics mechanisms [2-4]. The goal of the present investigation is to determine whether the computational time and storage savings associated with use of one of these recently developed simplified but accurate multi-step kinetics mechanisms are sufficient to allow solution of a combustion problem which is not presently possible using detailed, elementary kinetics. The kinetics mechanism chosen for this purpose involves 4 reversible steps and only 7 species, as shown in Table 1. This mechanism, to be refe~ed to as the 8-step mechanism, was developed by Kiehne et al. [4] following experimental research by Dryer and coworkers at Princeton [5, 6]. It has been shown to adequately represent the major features (flame speed, quench distance, flame temperature, total energy release, postflame composition and ignition delay) of both freely propagating and quenching flames over a wide range of temperatures, pressures and equivalence ratios. Although propane is used for the present work, research at Princeton [6] indicates that development of similar mechanisms for higher molecular weight fuels is possible. The specific problem chosen for the present investigation is a simplified simulation of the physics of combustion in a homogeneous charge spark ignition engine. This problem was chosen for study because the long-term goal of the authors is to implement the 8-step kinetics in an engine simulation code. However, a complete engine model requires a fully elliptic representation and additional equations to account for the turbulent flow field. These requirements obviously add to the computational complexity of the model. If the time savings associated with the use of the 8-step kinetics is not first shown to be sufficient to allow solution for a less demanding problem, then any attempt to solve the more complex problem is obviously premature. Further, if success is shown for the present problem, it is anticipated that other investigators will be encouraged to apply one of the recently available simplified but accurate kinetics mechanisms for solution of combustion problems which previously could be solved only by resorting to inaccurate global chemistry.

T.M. Kiehne et al., Numerical solution technique for combustion with multi-step kinetics

11

Table 1 Simplified but accurate kinetics m e c h a n i s m and rate expressions used in the numerical simulation [4] ( D i m e n s i o n a l units used are cal, gmole, cm, sec and K) 2C3H s ~ 3C2H 4 + 2 H 2 C2H 4 + 0 2 ~-- 2 H 2 + 2 C O 2CO + 0 2 ~ 2CO2 2H 2 + 0 2 ~ 2 H 2 0 with forward rate expressions R, R2R3S 4-

f , ( P ) 2 . 0 8 9 x 10 '7 exp(-49,600/RT)[C3Hs]°'5°[O2]"°7[C2H4] °'4° f 2 ( P ) 2 x 10 '3 exp(-50,000/RT)[C2H4]°'9°[O2]'"S[C3Hs] -°'37 S ( 0 ) 1.5 x 10 '3 e x p ( - 4 0 , 0 0 0 / R T ) [ C O ] " ° [ O 2 ] ° ' 2 5 [ H 2 0 ] °'5° 3.311 x 10 '3 e x p ( - 38,100/RT)[H2]°'85[O2] 1"42[C2H4] -0"56

and reverse rate expressions S 5 - 4 . 9 2 0 x l 0 s exp(-49,600/SY)[C3Hs]°"27[O2]t'°7[C2H4] °'4° R 6 - 2.25 x 109 exp(-50,000/RT)[C2H4]°'52S[O2]L's[C3Hs]-0"37 R 7 - 4.16 x 10 '6 T-1/2 e x p ( - 106,950/RY)[CO2]l'°[O2]-°'25[H20] 0"50 S s - 6 . 1 2 x 10 '5 T -1/2 exp(-lOO,586/RT)[H2]-°'z53[O2]°'9'6[C2H4]-°'s63[H20] ''°

where f t ( P ) ~ 6.434P -°.s'16 f 2 ( e ) - 1.115 - 1.125 exp -°'25P S ( 0 ) - min[1.0, 16.0 e x p ( - 2 . 4 8 O ) ]

Regardless of the level of complexity modeled, three features of an engine-like combustion process are identifiable as significant physical characteristics. First, the process occurs in a relatively thin-gap near Top Dead Center (TDC). Second, the combustible medium is subjected to momentum interactions with a moving wall while confined within a closed container. Third, thermal interactions with the walls have significant influence on the chemistry of the combustion process. The physics of these processes impact upon one another in an interplay which is not well understood, primarily because accurate chemistry has not been used in such models and because of uncertainty regarding the turbulence. Thus, the problem chosen for study is the transient propagation of a laminar flame away from a spark kernel on the centerline near a stationary w~l as the opposing wall moves first toward and then away from the stationary wall. The thin-gap feature is employed to develop the simulation used for the present research. During the combustion process the ratio of the gap near TDC (8) to the cylinder radius (R) is small (~e~, 8/R <0.1). This allows elimination of terms in the governing equations resulting in transformation of the character of the governing equations from fully elliptic to parabolic in the radial coordinate. For a typical engine, the model is then valid for approximately 20° before to 20° after TDC. The 40° of total crank angle is characteristic of the combustion event while the 20° before TDC is within the range of spark advance for many spark ignition engines. The motion of the moving wall serves the same function as that of a piston in an engine. The moving wall makes this problem more complex and interesting than would be the case if a simple laminar flame propagation process in a fixed geometry were studied. Both planar and axisymmetric geometries are studied. The transient,

12

T.M. Kiehne et al., Numerical solution technique for combustion with multi-step kinetics

two-dimensional formulation includes moving wall, shear flow and boundary layer effects with the temperature specified independently at each wall. The model allows for a fundamental investigation of flow, wall temperature and heat transfer effects on evolution of the chemistry in the near-wall regions. It should be emphasized that this investigation is fundamental in nature and is not intended to be a multi-dimensional engine combustion model. The intent is to advance the level of sophistication in modeling of combustion and fluid dynamic interactions with stationary and moving walls. The primary goal of this research is to demonstrate that the use of the 8-step kinetics mechanism allows solution of combustion problems that were previously intractable without resort to gross oversimplification of the chemistry.

2. Mathematical model

The combustible medium is assumed to be a uniform premixture of ideal gases near stoichiometric conditions. The fuel is taken to be propane since it exhibits characteristics (flame speed, minimum ignition energy and flame temperature) similar to higher saturated hydrocarbons. The walls are assumed to be noncatalytic. Thus, only homogeneous gas phase reactions are considered. Further assumptions involved in formulating the model are addressed below. 2.1. Basic equations The most general conservation and transport equations which apply can be found in e.g. [7, 8]. If the azimuthal swirl aspects of the flow are ignored, a simplified, two.dimensional, transient set of governing equations results. Letting a = 0 denote the planar case and a - 1 the axisymmetric case, and with superscript • denoting a dimensional quantity, the set of governing conservation equations is Conservation of mass: Op * + 1 0 ( p , r , ~ v , ) + 0 , Ot* r *~ Or" ~ (p*vz)=0 ;

(1)

Conservation of momentum: r-direction:

o,(o~,* \'~

oo:

a~*~

+ v* Or* + v* - ~ / =

ap* +[ 1 Or*

r*" Or*

¢ , , ~ , , ) _ ~ "oo 7

+ Oz* J

(2)

z-direction:

- ~ + v* - ~ + v* oz*/=

~;

+

r ;~ Or*

"~

'

(3)

T.M. Kiehne et al., Numerical solution technique for combustion with multi-step kinetics

13

Conservation of energy: , _,[ aT*

_

1

r;"

aT*

aT*

0 A,r,= Or* Or*

+

s

[ ,

"~

Oz* / J

, 0T*] O (r*=T * ) + v . 'aFz* J r*= Or*

, , 1 Y1 Cpi [}ir

- P* ~

i=1

$ - ~ w •i h •i + i=1

(op*

OO~ , D *r , Ov* rrr* Or* + @T0O -;z + Tzz az* + ¢'z~T~r* + oz* /

ap* +

+ ~+v*

Or*

op*].

(4)

V*az,/,

Conservation of species"

Op*i

Ot* =

+

1

a

O

(p*,*=.*) +

, , (p,u,)

r *= Or* ' - -- • 1 a **=, O **] r ;= Or* (p ir Dir) "~" " ~ ( p i o i z )

, + W i"

(5)

The components of the Newgonian stress tensor which appear above are given by *

Trr :

/z*

[2 Or* Or* -

T:=* = /Z* 2 aZ*

]

coo =/z

3 (V" V*) ,

T,, =/Z* ~

where

v. v* = 1

a (r,=v,) +

r *= Or*

r

0v: OZ*'

,

,

[2a v,r*'

2 "3 (V. V"~) ,

2 (V.V*)

3

+ aZ* J ,

]

(6)

(7)

(8)

In the above equations, the symbols are those in conventional usage with A* being the mixture thermal conductivity,/z* the mixture viscosity, h* the specific total enthalpy of species i, o]' the species diffusion velocity, W* the species molecular weight, w~' the mass rate of production of species i calculated from the rate expressions in Table 1, and the upper limits S and R on the summations indicate sums over the species and reactants, respectively. The effects of body forces, radiative transport, Soret effects, Dufour effects, pressure diffusion, forced diffusion and photon emission or absorption have been neglected. The assumption of ideal gas behavior has been made to write the energy equation in the form shown. The required equation of state (the multicomponent Ideal Gas Law in terms of the universal gas constant and Y*, the species mass fractions) is given by p* = p*lT*

i=IE W? ] "

(9)

14

T.M. Kiehne et al., Numerical solution technique for combustion with multi-step kinetics

In addition, sub-models for the diffusive and thermal transport properties of the fluid are necessary. The above equations then form a tractable set for the 5 + S dependent variables (p*, p*, T*, V* (v*, v*) and Y*) in terms of the independent variables (r*, z*, t*). The solution is found by solving the 5 + S equations ((1)-(4) and (9) along with S-equations with the form of (5)) with appropriate boundary and initial conditions. 2.2. Development o f the thin-gap equations

The system of governing equations can be further simplified by first nondimensionalizing with appropriate time, velocity and length scales, and then invoking the thin-gap approximation. Denoting the initially uniform mixture conditions with a subscript (o), letting ~o be the maximum crankshaft angular velocity and Vp the maximum piston linear velocity, allows definition of the nondimensional variables: r*

z*

r = -R- '

T*

T=--

p*

z = - 8- '

To ' C*

Cp ~- C?p~,

Yi--

t = o~t*,

Y*

h'=Cporo

*

O r

Po '

v,= (VpR/8)'

p*

Yio ' *

Cpi "~ C pi

/z*

P = poO;,R Wi =

C.o' h~'

p=

Y '

W*

'

,~__ ao '

(10)

D i] - D ij

n,o' 01,*

w~'

Cporo'

A*

*

Wo'

Cpl T*

~o, m.~

V*7 vz= Vp'

w, WjR, poo,'

o,, (D,o/R)'

• v~z

o,, (O,o/8)"

Using the above, the equations are then simplified by using the thin-gap approximations resulting from constraining the model to those cases for which the ratio of the gap (8) to the cylinder radius (R) is small, i.e., 8/R <0.1. Substituting (10) into (1)-(9) and neglecting terms of O(8/R) or higher, the thin-gap equations are Conservation o f mass: Op 1 0 (pr=o,) + a [3 - ~ + r-~ 0"~ "~ ( P Vz ) = 0;

(11)

Conservation o f momentum: r-direction:

-ff+VrT;+v'

oz

0 Or

ore

#z

" "~'z/ '

(12)

z-direction:

ap

0z =0;

(13)

T.M. Kiehne et al., Numerical solution technique for combustion with multi-step kinetics

15

Conservation of energy: 0T 0T 0T 13 T i + v" -g; + v" oz

s

1

1 [0-~ (A-~zT)] - Cp PeLe oCpPe +

(Y-1)M2 (R)2 [/~ t~P pCp

i=1

s y~oYiCp,V,~ OToz

op

~Or, ~' W~R,hi

P,-,p i=1 ;

+ Vr " ~

(14)

Conservation of species: OYi OYi OY~_ [3 -~i" + v" "~" + vz Oz -

1 pPeLe

[ O

"0-'z'(pY,

+ 7 ,o w, ;

(15)

Equation of state: =

s (Y, oY,

P 'yWoM2(R)2[p/(CpTi~=I ,'W-~i ,] •

(!6)

Where the nondimensional groups are defined as Gap Reynolds number = Re = poVpSIizo ,

(17)

Prandtl number = Pr =/~oCpo/Ao,

(18)

Peclet number = Pe = Re Pr,

(19)

Lewis number = Le = Ao/(poDaoCpo),

(20)

Mach number= M= Vp/V?-RTo = VplV(7- 1)CpoTo,

(21)

Nondimensional frequency parameter =/3 = 8oo/Vp.

(22)

These governing equations are essentially unsteady, thin-shear layer equations for a reacting flow with each transient term multiplied by the nondimensional frequency parameter. Application of the thin-gap approximation has changed the character of the governing equations from fully elliptic to parabolic in the radial coordinate. With the imposition of suitable boundary conditions, the above equations allow for the determination of the importance of shear flow and wall thermal effects on laminar flame development in a thin-gap moving wall geometry when the solution is restricted to regions near the centerline of the computational domain where end-wall effects are unimportant.

2.3. Initial and boundary conditions The governing equations are first order in time, first order in the radial coordinate and second order in the axial coordinate. Therefore, the equations require an initial condition, a

16

T.M. Kiehne et al., Numerical solution technique for combustion with multi-step kinetics

starting condition from the centerline and two conditions at each surface for all dependent variables except axial velocity, v~, which appears only to the first order in (11). However, physically the model contains two conditions on v~ one at the moving piston and one at the stationary cylinder head. This apparent paradox arises from the reduction of the axial momentum equation and is resolved in the next section when a compressible stream function formulation is introduced. The boundary conditions can be stated in dimensionless variables where at the stationary wall v,(r, O, t) = O,

vz(r, O, t) = O,

T(r, O, t) = T h ,

oY,(r,O, t) Oz

=0;

(23)

at the moving wall v,(r, zp, t) = O,

vz(r, Zp, t) = Vp(t),

T(r, Zp, t) = Tp ,

0 Yi(r, zp, t) Oz =0; (24)

and at the centerline

v,(0, z, t)= 0,

7"(0, z, t)= To,,

Y,(0, z, t) = Y,o,.

(as)

The radial domain is open and unspecified at infinity since the problem is parabolic in the radial direction. The problem allows for either Dirichlet or Neumann boundary conditions on T which may vary with time and radial position. Thus, in lieu of the specification of temperature at each wall, the wall heat flux could be specified. Only fixed wall temperature conditions were imposed in the present study. The Neumann condition on mass fraction at each wall implies no mass flux and the absence of wall reactions or transpiration. The initial conditions are based on the state of the combustible mixture prior to ignition. The temperature and composition of the mixture before ignition are reasonably expected to be uniform and homogeneous. Therefore, stipulation of a constant temperature and mass fraction distribution in the mixture is justified. In an engine the development of the flow field prior to ignition is extremely complex and is influenced by the intake and compression processes which have occurred prior to ignition. However, for the present purposes, the initial velocity field in the region of the centerline is assumed to be stagnant. With this as a base condition, complexities can be added to the initial velocity field at will. The model is ostensibly capable of accepting any reasonable initial conditions on velocity and could conceivably be supplied from a separate code. Employing the above assumption, the initial conditions in dimensionless form are v,(r, z, O)= O,

T(r, z, O)= 1,

Yi(r, z, O)= 1.

(26)

2.4. Stream function transformation

The governing equations are essentially one-dimensional flame equations in the axial direction coupled to the radial momentum equation which allows for a shear flow to be

T.M. Kiehne et al., Numerical solution technique for combustion with multi-step kinetics

17

imposed. It is possible to retain the momentum, energy and species equations in a common form while eliminating the conservation of mass equation through use of a general, compressible stream function transformation. The transformation is found under various names in the literature. The original reference appears to be attributable to Dorodnitzyn [9], whose intent was to avoid the explicit calculation of the density in the three-dimensional, transient boundary layer equations. The Dorodnitzyn transformation, written in nondimensional variables and modified to accomodate either the planar or axisymmetric geometries and the unsteady frequency parameter, is

~: = r, vr =

~/= 1

pr ~

a¢ az

'

t

p dz,

• = - •

(27)

/3'

oz = - -

1(1 p

r~

a¢ 07/) +/3 Or -~ .

(28)

The transformation from physical coordinates (r, z, t) to stream function coordinates (~, 7/, ~-) is depicted graphically in Fig. 1. In terms of the transformed variables, v, and v z are given by

1

111

]

(29)

where the 'z' derivatives are

1

°I;1 Z

V////~ ~

I

/ / / / / / / / / / / / / / / / / / / / / ~

Dependent Verloblee: v r-- v r (r,z,t) v z = v z (r,z.t) Yi = Yi (roZ.t) T = T (r,z,t)

hill

I

r

ki

PHYSICAL PLANE

Dependent Variables: ,,,

Yi T

•

,,

= =

c~,,~,~)

Yl (~'q"~) T (F,oq.~)

T

"qp(~.'t) $

*

TRANSFORHEDPLANE

Fig. 1. Graphical depiction of the physical and transformed coordinate systems.

18

T.M. Kiehne et al., Numerical solution techniquefor combustion with multi-step kinetics

Here, and in subsequent equations, the transformed variables when used as subscripts indicate differentiation with respect to the variable specified. When the transformation is applied to the governing thin-gap equations, the result is Conservation of momentum:

o

¢.

~" on

=

p 8~: + R e On / ~ p ~

;

(31)

Conservation of energy: 8,

~:" 8~

~:" on = CpP-'-"e

AO

('y-- 1)M 2

s + CpPe Le 8n ~=1 p2

aY,

8 T E YioDijCpi an

'~2~"ap

,__, w,n,h, +

+#¢': 0# + ~

N

;

(32)

Conservation of species: "0,r + ~'" 0~

~'~ 0n = Ve L-"--e ~

p2Dij ~

+~io

Wi;

(33)

Equation of state:

w, //J'

(34)

The conservation of mass equation is satisfied identically by the transformation and therefore need no longer be directly considered. The axial momentum equation simply implies that the axial pressure gradient is essentially zero and, while used repeatedly during the transformation, is no longer stated explicitly. In the above, an expression based on Fick's Law of Diffusion is used to relate the axial diffusion velocity to the mass fraction gradients and binary diffusion coefficients. In transformed variables, the expression is

1 aY~

v~ = -PDii ~

On '

(35)

where D~ remains to be modeled. The use of Fick's law is the simplest approach to modeling of the diffusion velocity. The work of Heimerl and Coffee [I0, 11] for hydrogen and methane flames indicates that use of a diffusivity model based on Fick's law is computationally economical and totally adequate in comparison to more involved methods. It is particularly appropriate when thermally induced mass diffusion effects are being ignored as is the case here. The boundary conditions for the transformed problem are

19

T.M. Kiehne et al., Numerical solution technique for combustion with multi-step kinetics

at the moving wall

,/,,,(~:, rip, ~-) =o,

l~lr~(~:, rip, T) ~--"~ a [ p ( Z . r -- ~p)](~. Wp.,r),

(36)

o Y,( ~, ,~, ,)

T(*,,,,, ~') = T,,,

=0;

at the stationary wall

,/,,,(~:, o, ~-) = o

¢d¢,o, ~-) =o

'

'

T(~', o, ~-) = r,,

'

or,(~', o, ~-) o~/

=0;

(37) and at the centerline

,/,,,(o,,1, ,,-)= o,

T(0,,~, ,,-)= L,,

Y,(0,,1, ~-)= Y,o,.

(38)

In (36), ~/p is the axial location of the moving wall in the transformed plane as given by ~p =

o(r,

z, 0 dz.

(39)

This integral represents a new unknown introduced through the transformation. The transformed piston location is not known until the density field is calculated. Therefore, the point at which the moving wall boundary condition needs to be applied is tied directly to an unknown in the problem--the density field. In essence, the above equation requires that the moving wall be positioned in the transformed coordinates such that mas,~; conservation is satisfied. When transformed, the initial conditions become

,/,,,(~, ,1, o) = o,

7"(~, ,~, o) = t,

Y,(~, ,~, o) = 1.

(40)

The final step in the problem formulation involves eliminating the radia~Lpressure gradient, ap/O~, from both the momentum and energy equations. This is accomplished for the momentum equation by differentiating (31) with respect to ~/and then using the fact that ap/O~ =0. This step is similar to that used in the traditional strear:a function-vorticity formulation. The final result after grouping of terms and rearranging is 1

0e (¢,,) - p

+

a,la6

Re OT/2

/f" 017

Re 0,1 0,1 ~p -~'JJ"

(41)

20

T.M. Kiehne et al., Numerical solution technique for combustion with multi-step kinetics

The density gradients which appear in (41) are determined by differentiating (34) with respect to 7. The final form of the momentum equation is fourth order and now allows for the four velocity (now stream function) boundary conditions stated previously. The pressure no longer appears explicitly in the momentum equation. However, it enters implicitly through the transformation. The effect of the pressure gradient on the flow field has not been neglected but has been incorporated into the inertial, viscous and pressure forces which act in establishing the flow. The effect of the longitudinal pressure gradient in the energy equation, however, is henceforth neglected. The temporal variation of pressure in an actual engine typically varies by a factor of fifty while the spatial gradients are nearly zero. This is true since spatial pressure variations are equilibrated on a time scale which is much shorter than the time scale associated with a propagating flame (i.e., the acoustic waves which equilibrate spatial pressure variations propagate much more rapidly than the flame). Therefore, any local pressure fluctuations which may arise are equilibrated almost immediately. With this justification, spatial pressure variations for the energy equation are assumed to be zero. The derivative of pressure with respect to time which remains in the energy equation accounts for the work interaction of the moving piston with the compressible gas mixture. The temporal behavior of the fluid pressure must be known to evaluate this derivative. Since the governing equations are parabolic in the radial direction, the presence of the confining end-waUs which establish this pressure are not included. A model which predicts the correct pressure behavior would require formulation of the fully elliptic problem. This difficulty is circumvented by using a curve fit to an actual pressure-,time profile for a typical spark ignition engine to evaluate the needed derivative and absolute pressure. 3. Numerical solution

The governing equations are solved by a variation of the standard semi-discrete approach or what is called the Method of Lines. There is an extensive literature--mostly Russian--for the Method of Lines. A comprehensive but dated review was done by Liskovets [12]. More recent work in the combustion area [13-15] has focused on the one-dimensional flame equations. In these works, the discretization is performed in the space dimension and the equations are maintained continuous in time. The formulation used here differs from this established Method of Lines approach. In this work, discretization is performed on two independent variables. Specifically, the governing equations are discretized in the parabolic or marching variables--~ and ~-. The resulting equations are then continuous in the axial direction, 7, and consist of a coupled set of two-point, nonlinear boundary value problems. This approach permits expeditious handling of the axial flame propagation on a manageable adaptive grid. Without this numerical scheme, it would be extremely difficult to resolve the strong axial temperature and velocity gradients. While this procedure has been alluded to in the literature [16, 17], no actual use of the technique has been found. Differencing of the governing equations in the time and radial coordinates is relatively straightforward. Since there is a principle flow direction, upwind or donor cell differencing is used in the radial direction. This technique avoids coupling of grid values for dependent variables ahead of the current continuous line of interest. The temporal differencing is

21

T.M. Kiehne et al., Numerical solution technique for combustion with multi-step kinetics

backward in time and naturally yields to a time marching formulation. Since the boundary value problem code which will be used requires isolation of values for the highest derivative in each of the governing conservation equations, the appropriate difference equations are given in the form shown below.

Conservation of momentum:

Of,,.=

Re [ (dh-)-/~p h ( 0 " - O~) + d ~*' +

(,.

*" (~, - ~'b) -- 2a dh -,-'~, " I ] - ~ ) - d -~-

1 ell Re (--~-) [h(O' - ~,~) + d "~-~ (0' 0~,) d "~---a (0 - 0b ) -- a dh --~-j (dh)/~p

/£P [02(,U,p)+ (._~)O(/~,p)] .

072

07

(42)

'

Conservation of energy: CpPe

T W# . . ~ (dh)Ap

~'

T'

[h(T- T~) + d -~--a( T - Tb) -- d ~-'a ( , - *b)

]

s , /3Pe s P r' 2 YioDoCpiY, + Ap2 i=1 E WiRih, A Le i=z

(Y-1)M2pe(R) Ap2

Conservation of species (i = 1-7):

Pe Le [

Y'i

¢'

Y'; = (d~-~-Do h(Y,- Y,,) + d ~ ( Y , - Y,b)- d ~ ( , /3RiPe Le p Do o W,

y[ p2Dij

O(p2Do) 07

.

*b)

] (44)

The notation used in these equations is shown in Fig. 2. In each expression, primes indicate differentiation with respect to 7, subscript c indicates the value of a dependent variable at the current radial location and previous time step, subscript b indicates the value of a dependent variable backward in radial location but at the current time, h is the radial space step (A~), and d is the time step (A,.). The differencing used is accurate to first order and written so that only the solution for the previous time at each line need be retained in computer memory. The property derivatives with respect to the transformed axial coordinate arise from chain rule differentiation involved in isolating the highest derivative of each dependent variable. The boundary conditions which involve radial gradients of the stream function must also be finite differenced to allow for their application on a discrete radial grid. The condition on the

22

T.M. Kiehne et al., Numerical solution technique for combustion with multi-step kinetics PREVIOUS SOLUTION IN SPACE (SUBSCRIPT"B')

~ I ]

ITI ~-I

/

,7,

•

,

f

LINE"ALONG WHICH I IiI [~ I i k I ' 1 CURRENT

/

I Iii I~ l i 1 " 4 I SOLUTION I I~1 [ I~1 I~-:SSOUGtn

/ /V

.

EVIOUS SOLUTION (SUBSCRIPT""C')

Fig. 2. Finite difference notation used in discretising the governing equations.

radial gradient of the stream function at the stationary wall in (37) becomes O(/j, O, ~) = O,

(45)

where an arbitrary choice of zero stream function at this wall has been made. This choice establishes a reference for the stream function to which all subsequent values are now tied. When differenced, the companion condition on radial gradient of the stream function at the moving wall from (37) becomes

Pc

~p, 1')

,

(46)

where the notation is the same as that used for the governing conservation equations. Again, the differencing is first order accurate. The above expression indicates that the stream function at the moving wall is established by an integrated compressibility effect as well as by the wall motion. It also ties the value of the stream function at the moving wall to the density field both directly through the integral and indirectly through another integral, (39), which establishes the transformed piston location, 'l~p. The transformed and differenced governing equations are a coupled set of nonlinear, two-point, boundary value problems in the axial coordinate. The equations are coupled to each other through the dependent variables and to previous solutions in both time and space by differencing. The set consists of a fourth order equation for the stream function along with eight second order equations for the temperature and mass fractions of the seven reacting chemical species. The density of the nonreacting species, N2, is employed to assure conservation of mass. Thus, the system is of order 20. In addition to being nonlinear, the equations are expected to be 'stiff' due to the presence of chemical production terms which have an exponential dependence on the temperature. The appearance of the higher derivatives of the dependent variables in the governing equations imposes the added difficulty of accurately determining these derivatives. All of these complications place demanding requirements on the solution algorithm to be employed. COLSYS [18], a collocation solver for mixed order systems of boundary value problems,

T.M. Kiehne et al., Numerical solution technique for combustion with multi-step kinetics

23

was selected to accomplish the solution of the governing equation set. COLSYS implements a method spline collocation at Gaussian points. Nonlinear problems are solved using a damped Newton's method of quasilinearization. Jacobian evaluations to support the Newton iteration are capable of being managed with a finite difference formulation. One of the most attractive features of COLSYS for the present problem is the fact that the equations are solved and all derivatives evaluated on a sequence of refined meshes until a user specified set of tolerances is satisfied. Thus, adaptive grid selection is accomplished automatically. Mesh redistribution is performed to roughly equidistribute the error in the collocation solution. Further details on this package are available in the documentation. An array of submodels is required to support calculation of a solution to the problem posed above. The differenced equations and boundary conditions require knowledge of/~, A, Dij and Vp. In the interest of facilitating the initial development of the present concepts, extremely simple submodels primarily based on absolute temperature are used for the property variations. Thermochemical calculations and chemical rate equation evaluations are performed using a version of CHEMKIN [19] modified to accomodate the nonstandard rate equations of the simplfied propane kinetics scheme. The thermochemical data base used is that of Gordon and McBride [20]. The instantaneous position and velocity of the moving wall are based upon the kinematics of a 4.1 liter, Buick V6 engine.

4. Results and discussion

4.1. Nonreacting flow results To assess the viability of the numerical approach, a series of nonreacting flow cases were computed. These cases used a fixed temperature, at 500 K, and constant pressure. Therefore, the flow is incompressible since no density or composition changes occur. Calculations were performed for fixed pressures of I, 10 and 40 atmospheres with the wall motion simulating an engine speed of 200 rpm for both planar and axisymmetric geometries. The rationale for choosing this speed is explained in the next subsection. Under these conditions, the piston velocity at 20° before TDC (BTDC) was approximately 40cm/sec. These conditions correspond to gap Reynolds numbers of 72, 720 and 2880, respectively. For the I and 10 atmosphere calculations, the time step was fixed and corresponded to one-half degree of crankshaft rotation. The time step for the 40 atmosphere case was likewise fixed but at one-quarter degree of crankshaft rotation to ensure computational stability. The initial nonreacting flow calculation produced an anomaly in the velocity profiles during the first few time steps. Specifically, sharp peaks in the velocity profile near the moving wall occurred at the beginning of each calculation. This problem was traced to the impulsive motion that the initially stagnant fluid experienced when the moving wall was suddenly accelerated. The behavior was very similar to that of an unsteady shear wave in the classical Stoke's problem of impulsive planar wall motion. The effect was insignificant at low piston velocities (low rpm) but quite pronounced at higher speed. Although the velocity peaks disappear with time due to viscous diffusion, a modification in the initial velocity field was incorporated into all calculations in order to more accurately simulate conditions which might exist physically prior to ignition in an engine. The method selected was to linearly ramp the

24

T.M. Kiehne et al., Numerical solution technique for combustion with multi-step kinetics

piston velocity from zero at 24 ° BTDC to its actual value at 20 ° BTDC. This scheme allowed the velocity field to develop gradually and removed all vestiges of the anomalous velocity peaks near the moving wall. The nonreacting flow calculations, interesting in their own right, provided valuable lessons in refining both the mathematical model and the computer code. The radial and axial velocity profiles for the 10 atmosphere, axisymmetric case are shown in Figs. 3 and 4 in increments of 4 ° of crankshaft angle beginning at 18 ° before TDC (piston velocity of - 3 5 . 6 c m / s e c ) to 6 ° after TDC (piston velocity 12.1 cm/sec). The planar results (not shown) are qualitatively similar in behavior. A key point to remember throughout the following discussion is that prior

18°

00J 14° BTD¢

I0 ° BTDC

6° BTDO

_t--

~~

~,

2 ° BTI)C

i

o

w

i

i

i

q

I

•

il

18° BTDC

2 ° ATDC

"~ I0 • ,., -

~

6° ATDC

' L 3 -_ _ _ J 0-30

;~0 -I0

0

I0

20

30

RADIAL VELOCITY(CM/S)

Fig. 3. Radial velocity distributions for nonreacting axisymmetric flow. The lines represent radial locations of 1, 2, 3, 4 and 5 % of the cylinder radius away from the centerline.

•

I

-40,0 -300 -200 -I0.0 O0 IOO AXIAL VELOCITY (CMIS)

20,0

Fig, 4. Axial velocity distributions for nonreacting axisymmetric flow at a radial location of 5 % of the cylinder radius away from the centerline.

T.M. Kiehne et al., Numerical solution technique for combustion with multi-step kinetics

25

to TDC the moving wall is continually decelerating. Conversely, after TDC the moving wall is continuously accelerating away from the fixed surface (z = 0). In Fig. 3, the ordinate is the nondimensional gap which changes with time as the upper surface moves. Each line in this figure is an increment in radial distance corresponding to 1% of the radial dimension. Therefore, the last line is approximately 2.5 mm removed from the centerline for the reference conditions employed. Except for computer time and storage limitations of the CDC Dual Cyber 170/750 used for these calculations, there is no reason why numerical predictions could not have been continued to larger radial distances. At 18° BTDC, the uniformity of the profiles over the lower two thirds of the gap indicates that viscous effects have not yet diffused from the walls to influence the fluid momentum in this region. At 14° BTDC, viscous effects have further infringed upon the uniform central core and the peak velocities have increased accordingly. These basic observations continue to apply at 10° BTDC. However, the peak velocities at each line are now essentially unchanged from those at the previous time. Between 14° and 10° BTDC the profiles have peaked and the effects of the decelerating wall are now being felt. This is even more apparent at 6° BTDC where the shear stress at the moving wall and the peak velocities have obviously begun to decrease when compared to the profiles from the previous time. Near the stationary wall the decelerations have not yet been felt due to a viscous time lag. At 6° B'I~DC, the uniform velocity region near the stationary wall has essentially vanished. At 2° BTDC, a major change in the nature of the flow field is apparent. The piston is now rapidly decelerating as the approach to TDC continues. The fluid near the moving wall has responded to this deceleration by reversing its direction of flow. This flow reversal begins at approximately 3° BTDC. Development of the profiles following flow reversal is of particular interest. The effect of the surface deceleration up to and then acceleration away from the fixed surface is rapidly transmitted through the fluid. At 2° after TDC (ATDC), complete reversal has occurred except at a saddle-point or cusp which is a remnant of the peak upstream velocity at 2° BTDC. At 6° ATDC, the flow field has been entirely reversed in direction and the maximum gradients are now at the moving wall as the surface velocity away from TDC increases rapidly. The corresponding planar case exhibits essentially the same behavior except that the development of the flow and the onset of flow reversal are advanced by approximately one-half degree of crankshaft rotation. A valid issue with the flow reversal is the advisability of continuing the calculation in light of the upwind differencing utilized. The issue is one of properly reflecting the convective forces which the flow experiences. In a strict sense, calculation of the derivatives in the finite difference formulation should involve upstream information upon flow reversal. However, the boundary conditions and geometry in the incompressible thin-gap case stipulate that these derivatives are unchanged at a given time as a function of radial location. This is most easily seen by considering the axial velocity profiles shown in Fig. 4. Although the profiles are shown for a radial position corresponding to the last radial line of Fig. 3, the axial velocity varied little from one radial position to another at a fixed time. For example, the maximum change in the profiles from the first line to the last at 10 atmospheres is approximately 2% and occurs near the stationary wall. Therefore, the axial gradient of the axial velocity is nearly constant with radial location. Consideration of the mass conservation equation for the incompressible case then dictates that the radial gradient of the radial velocity is also a constant. If both of these gradients are constant, then it is immaterial whether information is passed from upstream or downstream upon flow reversal since the information is unchanged with radial

26

T.M. Kiehne et al., Numerical solution technique for combustion with multi-step kinetics

location. This result follows directly from the fixed boundary condition at both surfaces. Therefore, for the incompressible nonreacting flow situation, continuation of the calculation past flow reveral is justified. Of course, these conditions do not apply in the reacting flow case since variations in density and composition must be considered. Therefore, reacting flow calculations were not performed past flow reversal. In addition to the above case, constant pressure nonreacting flow calculations at 1 and 40 atmospheres were made to assess the Reynolds number effect on the flow field and upon numerical stability. Also, variable pressure calculations using engine pressure-time histories were performed to test the effect of a varying density field as a prelude to the reacting flow calculations. Although these results are not shown graphically, a brief discussion is of interest. In the 1 atmosphere case, viscous effects influence the development of the central core more rapidly. As a consequence, the flow reversal is delayed as compared to the 10 atmosphere case. In contrast, the 40 atmosphere case has a gap Reynolds number which is four times that of the 10 atmosphere case. Therefore, inertial forces dominate the viscous forces. This feature is apparent in the profiles as a uniform central core persists all the way to the onset of flow reversal. As a direct result, the shear layer at the stationary wall is very thin and the flow is forced to turn upstream more rapidly. The lack of influence of viscous effects also leads to an earlier flow reversal as the fluid reacts to the convective decelerations imposed by the decelerating surface. By 2° ATDC, the flow has not yet entirely reversed but remains qualitatively similar to the 10 atmosphere case. In general, the wall shear layers at the higher pressures are extremely thin which explains the need for a smaller time step to resolve the evolution. The 40 atmosphere, axisymmetric results exhibit a characteristic which may have been present at the lower pressures but is now more prominent. Specifically, there is a slight concavity to the profiles near the moving wall which persists, though shifte,d downward, at 2° BTDC. More importantly, from a computational perspective, the high Reynolds number calculations exhibited none of the numerical instabilities which are common in fully discrete calculations. The 1 atmosphere, 40 atmosphere, and variable pressure nonreacting flow results are discussed in further detail in [21].

#.2. Reacting flow results Before discussion of the reacting flow results, issues related to the ignition of the mixture and selection of an appropriate surface velocity are addressed. Successful promotion of ignition in a combustible mixture is dependent upon the physical and thermodynamic properties of both the mixture and ignition kernel, the oxidation kinetics and the mode of ignition energy release. It was not the purpose of this study to exhaustively investigate the ignition phenomenon; although this issue could conceivably be addressed to a useful degree with the model. The conceptual objective is to deposit energy in the mixture in a direct manner while retaining some semblance of the physics of the ignition process. Thus, in what follows, the ignition kernel is modeled as a hot pocket of burned gases along 10% of the centerline adjacent to the stationary wall. This gas pocket is assigned properties corresponding to the equilibrium postflame region at the pressure of interest. Therefore, the kernel temperature is the adiabatic flame temperature and the mixture is dominated by CO2, H20 and N 2. Deposition of energy in the mixture then occurs by the process of convection as energy and chemical species which promote reaction are transported into the fluid domain.

7'. M. Kiehne et al., Numerical solution technique for combustion with multi-step kinetics

27

The rapidity of this transport is governed by the velocity of the moving surface. Therefore, a surface velocity was used which was low enough to allow for the flame-wall interactions to be adequately observed in the presence of a slowly developing laminar flame. That is, since the turbulent flame speed in an engine is of the order of a factor of ten higher than the laminar flame speed, the model used a piston velocity corresponding to 200rpm to simulate the duration of combustion in an engine operating at approximately 2000 rpm. The reacting flow case which is discussed below is for a stoichiometric mixture initially at 10 atmospheres and 500 K. The wall motion is identical to that of the nonreacting flow case just presented and the surface temperatures are fixed at 500 K. Of course, chemical reactions resulting in density and composition changes are now allowed. Figures 5 and 6 display the temperature, mixture composition and velocity profiles at a sequence of times (in crankangle degrees) following ignition of the mixture as described above at 18° BTDC. Forty subintervals

~..[ . . . . . . ~..413.304-13.289\ I . . . i ~ L 13.64-,

',

"

1

.

.

.

.

1

r ' - - ' -'--'J 13304 -'13.289-1~ / / / - - - - -~13.64 --7 // ii I

"~

~" ""'" "" "-

i~1

,

i

" "

"5 " " i

""

i

0

N

2

TEMPERATURE (K) 0~,(%)

Z '~' ~

~" rq

....

/

,,L.- l - - -I 15304-13.289-~

,--I---] 13.64-- ~ / /[15.18--I . . . . ;, ~\ // / / j I 6 4 5 ~ . - ~ / \

. . .a. . .

i o

4

i

CsH e (%)

4

6

0

2

~

I I I

,,-:

/

- - - ~3:3o~,-b3Z89-

/---

//"

,

,

6

8

:-, . . . .

13.64--I'-~

.-15.18,

~

4

Hz0 (%)

C02(%)

~

°o

/

-

m

zo

- 17.5\

,

I

\~

\

. . . . . .

,

,

,

,

,

,

0.00 0.02 004 0.06 0,00 0,08 0,15 0.23 0.30

C2H 4 (%)

H e (%)

CO (%)

Fig. 5, Temperature and species mass percent distributions at a radial location of 1 % of the cylinder radius away from the centerline for a sequence of times (in crankshaft angle degrees before T D C ) following ignition at 18° before TDC.

~ cJ ~

N N

O

'

,,,-"

,,.-17.5° BTDC---.,.. ,,,.16.450-....... ''" "---~,...---

I

II

I

! /

/ ~-- -13.304~ ~)~-' :15.295 °

°0.0 1.5 3.0 4.5 s.o ?.5 910 10.s 12.0-~-40 -32-24-16 -s 0

LONGITUDINAL TRANSVERSE VELOCITY (CM/SEC) Fig. 6. Velocity distributions for reacting flow at a radial location of ) % of the cylinder radius away from the centerline for a sequence of times following ignition at 18° BTDC.

28

T.M. Kiehne et al., Numerical solution technique for combustion with multi-step kinetics

with 200 collocation points were used to resolve the gradients in the axial direction. Computer storage limitatio,s of the CDC Dual Cyber 170/750 permitted calculations at only one radial position and prevented calculation of the problem through the rapid energy release associated with oxidation of the CO and H e in the reacting mixture. Therefore, the results shown are restricted to conditions near the centerline and for initial pyrolysis of the fuel along with oxidation of the products ,~f pyrolysis. Nonetheless, the model permits calculation of the effects of the flame-wall interaction on the early development of the flow field and chemistry in a wall dominated viscous environment. These results represent the first calculation of a transient, laminar, reacting flow in a moving wall geometry using nonequilibrium, multi-step chemistry while resolving the momentum and thermal boundary layers. The temperature profiles shown in Fig. 5 reveal the convection of energy from the original pocket of hot gases into the combustible mixture where diffusion in the radial direction occurs. The fixed wall temperature at the stationary wall leads to increasingly severe temperature gradients and rapidly increasing heat transfer near this wall. With the passage of time, automatic time step controls built into the computer code at first increase and then rapidly decrease the time step. At nearly 6.6 msec, the chemical activity has increased sufficiently to reduce the time step to levels commensurate with those previously experienced with a one-dimensional laminar flame code at 10 atmospheres [21]. As a result, the further development of the temperature profiles becomes so small as to be indistinguishable at the last time steps depicted. As rapid, exothermic CO and H 2 oxidation commences, the temperature would be expected to rise rapidly in a narrow zone near the peak of the existing temperature profile. The species profiles, expressed in mass percent in Fig. 5, show the development of the chemistry. The CO, and HeO profiles are dominated by convection from the centerline followed by axial diffusion. The fixed wall temperature produces a thermal boundary layer which causes the species concentrations at the wall to lag behind the peak upstream value due to the effect of temperature on the diffusion coefficients. The fuel initially undergoes pyrolysis very slowly. However, as the peak temperature rises to nearly 1200 K, very rapid fuel consumption occurs. Ignition delay studies for this kinetics mechanism [21] indicate that approximately 0.5 msec of delay prior to rapid energy release is involved at this temperature. This value is further confirmed by these fuel disappearance profiles. After the initial delay period, production of ethylene due to propane pyrolysis is also seen to be extremely rapid. As the ethylene is produced by pyrolysis of the fuel, its conversion to CO and H2 through oxidation is accelerated. Hydrogen is produced and oxygen consumed rapidly at the latest times shown. The hydrogen profile has begun to develop a shape similar to that observed in one-dimensional flame calculations. Specifically, an undershoot in the developing postflame region is becoming apparent. A series of minor peaks at the saddle-point in the hydrogen profiles is due to the onset of rapid conversion of fuel and ethylene to hydrogen. This feature is even more obvious in the carbon monoxide plot since all the production of CO is due to the conversion of ethylene. The early development in both the H 2 and CO profiles is due to convection but the onset of rapid chemical production is obvious. Even at the final very small time steps, the conversion to hydrogen and carbon monoxide is proceeding rapidly. The radial and axial velocity profiles shown in Fig. 6 are plotted against the entire axial domain to highlight the surface motion. The early development of the radial velocity, except

T.M. Kiehne et al., Numerical solution technique for combustion with multi-step kinetics

29

for an increased magnitude and a more severe gradient near the moving wall, is similar to that for the nonreacting flow. Near the stationary wall where the pocket of hot gases is established at the centerline, a wave in the profile develops early in the calculation and moves away from the wall until engulfed in the developing reaction region. This structure is believed to be caused by an interaction of the hot fluid with the cooler fluid trapped next to the wall. At the final times depicted, numerous substructures develop near the wall. The fluid on the opposite side of the developing hot region does not develop these structures which reinforces the belief that the origin of these substructures is an interaction with the wall. The axial velocity profiles indicate an enhanced flow of cold gases into the hot wall region. This explains the increased magnitude of the radial velocity as the additional flow is forced to turn downstream. The carbon monoxide profile reveals the origin of the storage limitation problem encountered on the computer available when these calculations were done. Insufficient computer memory was available to allocate additional subintervals for resolution of the rapidly developing chemistry. The present calculation was terminated when the need arose to further subdivide the grid to resolve the developing gradients. The fundamental problem is shown in Fig. 7, which shows grid point spacing versus axial location from the stationary wall. The grid shown is for the energy and species equations at the last time calculated and therefore has been established and refined by the adaptive grid algorithm over numerous time steps. The lack of truly severe gradients prior to the last short time steps in the calculation has led to a grid spacing which is small but relatively uniform in the near-wall region. The grid then expands to very large spacing in the region where no changes are occuring. The smallest spacing is in the zone of the peak ethylene gradient which is relatively slow in developing. However, the conversion of ethylene to hydrogen and carbon monoxide is very rapid following fuel disappearance. The grid selection algorithms have inadequate warning of this impending development and simply do not have sufficient time to prepare for it. The fundamental problem is that the grid must be halved which implies at least double the storage presently being used.

i 1,4

0.00

0.02

0.04

0.06

GRID SPACING

0.08

O.0

.12

(NONDIM)

Fig. 7. Grid spacing allocation at the last computed timestep (13.289 ° BTDC).

30

T.M. Kiehne et al., Numerical solution technique for combustion with multi-step kinetics

5. Summary A laminar, two-dimensional, transient, viscous flow model which includes finite rate chemistry has been developed to simulate the physical and chemical processes occurring during flame initiation and propagation near the centerline of a thin-gap formed by two cold, solid surfaces as one surface moves first toward and then away from the other. The primary motivation of this research was to determine whether the computational time and storage savings associated with use of a recently developed simplified but accurate kinetics mechanism were sufficient to allow solution of combustion problems which were previously too complex to model without resort to oversimplification of the chemistry. The model includes the effects of shear flows and boundary layers on each surface, accounts for nonisobaric conditions and uses a recently developed and verified chemical kinetics scheme which includes four reversible steps with 7 chemical species and accounts for the role of intermediate hydrocarbons. The practical considerations which have led to solution of the governing, thin-gap equations have been addressed. Two general purpose computer subprograms--COLSYS and CHEMKINm have been employed in an extensive driver program to solve the 20th order set of nonlinear, boundary value problems which result from semi-discrete numerical formulation of the problem. The development of the model and major features of the solution procedure have been discussed. The results of a series of nonreacting flow calculations at various fixed pressures for an axisymmetric geometry are presented. The major effect of pressure change is to alter the fluid density and therefore the Reynolds number. The flow then develops subsequent to the onset of surface motion in a manner which is understandable using standard Reynolds number considerations. Variable pressure, nonreacting flow calculations were also performed. The major effect of the non-isobaric condition is to alter the rate at which the shear flow develops. A reacting flow case computed up to the point of rapid production of carbon monoxide and hydrogen due to ethylene oxidation is also presented. Although the computer memory requirement to further resolve the spatial nature of the chemistry then exceeded the very limited capabilities of the computer used for this investigation, the authors are convinced that further calculations could be achieved using any of the more recently designed computers since these, in general, have a much greater storage capacity than the CDC Dual Cyber 170/750. Nevertheless, this reacting flow case is the first calculation of a transient, laminar, reacting flow in a time-dependent geometry using accurate, multi-step, finite rate kinetics while resolving the momentum and thermal boundary layers. It is believed that the success of this investigation should encourage other investigators to adopt one of the recently available simplified but accurate kinetics mechanisms to develop models for complex combustion systems that previously had to rely on oversimplified chemistry. Another facet of this investigation involved the application of the thin-gap approximations to a transient reacting flow. The, various future situations which may ultimately lend themselves to investigation with the thin-gap model are extensive. The impact of prescribed heat transfer, variable wall temperature, differing initial flow conditions, various initial mixture conditions and different chemical formulations are all possible. Conduct of ignition delay studies which include the physics of convection and diffusion are also of future interest.

T.M. Kiehne et al., Numerical solution technique for combustion with multi-step kinetics

31

References

[1] C.K. Westbrook and F.L. Dryer, Chemical kinetic modeling of hydrocarbon combustion, Progr. Energy Combust. Sci. 10 (1984) 1-57.

N N. Peters, Numerical and asymptotic analysis of systematically reduced reaction schemes for hydrocarbon

flames, in: R. Olowinski, B. Larrouturou and R. Temam, eds., Numerical Simulation of Combustion Phenomena, Lecture Notes in Physics 241 (Springer, Berlin, 1985) 90-109. [3] O. Paczko, P.M. Lefdal and N. Peters, Reduced reaction schemes for methane, methanol, and propane flames, Twenty-First Symposium (International) on Combustion, The Combustion Institute, Pittsburgh (1986) 739-748. [41 T.M. Kiehne, R.D. Matthews and D.E. Wilson, An eight-step kinetics mechanism for high temperature propane flames, Combust. Sci. Techn. 54 (1987) 1-23. [51 D.J. Hautman, F.L. Dryer, K.P. Schug and L. Glassman, A multiple-step overall kinetic mechanism for the oxidation of hydrocarbons, Combust. Sci. Techn. 25 (1981) 219-235. [6] W.M. Proscia, High temperature flow experiments and multiple-step overall kinetic mechanisms for the oxidation of high carbon number paraffin hydrocarbons, Masters Thesis 1625-T, Dept. of Mechanical and Aerospace Engineering, Princeton University, 1983. [7] R.A. Strehlow, Combustion Fundamentals (McGraw-Hill, New York, 1984). [8] F.A. Williams, Combustion Theory (Benjamin/Cummings, Menlo Park, CA, 1985). [9] A.A. Dorodnitzyn, Prikl. Mat. Mekh. 6 (1942). [10] T.P. Coffee and J.M. Heimerl, Transport algorithms for premixed, laminar, steady state flames, Combustion and Flame 43 (1981) 273-289. [11] J.M. Heimerl and T.P. Coffee, Transport algorithms for methane flames, Combust. Sci. Techn. 34 (1983) 31-43. [121 O.A. Liskovets, The Method of Lines (review), Differentiarnye Uravneniya (Differential Equations) 1 (12) (1965) 1662-1678. [131 G.R. Grey and H.A. Dwyer, Numerical study of the interaction of fast chemistry and diffusion, AIAA J. 17 (1979) 606-613. [14] S. Galant, Numerical simulation of unsteady laminar flame propagation via the Method of Lines: Further mathematical refinements and results obtained on halogen inhibition, Combust. Sci. Techn. 34 (1983) 111-148. [15] R.J. Kee, L.R. Petzold, M.D. Smooke and J.F. Grear, Implicit methods in combustion and chemical kinetics modeling, in: Multiple Time Scales (Academic Press, New York, 1985). [16] S.G. Mikhlin and K.L. Smolitskiy, Approximate Methods for Solution of Differential and Integral Equations (Elsevier, New York, 1967). [17] D.J. Jones, J.C. South and E.B. Klunker, On the numerical solution of elliptic partial differential equations by the Method of Lines, J. Comput. Phys. 9 (1972) 496-527. [18] U. Ascher, J. Christiansen and R.D. Russell, COLSYS--A collocation code for boundary value problems, Lecture Notes in Computer Science, Vol. 76 (Springer, Berlin, 1978) 164. R.J. Kee, J.A. Miller and T.H. Jefferson, CHEMKIN: A general-purpose, problem-independent, transport[19] able Fortran chemical kinetics code package, Report SAND-80-8003, Sandia Laboratories, New Mexico, March 1980. [2o] S. Gordon and B.J. McBride, Computer program for calculation of complex chemical equilibrium compositions, rocket performance, incident and reflected shocks, and Chap~nan-Jouquet detonations, NASA-SP-273, 1971. [21] T.M. Kiehne, Development of an eight-step kinetics mechanis,~',~for propane and its application in a transient, thin-gap, engine model, Ph.D. Dissertation, Dept. of Mechanical Engineering, University of Texas at Austin, 1985. [22] T.M. Kiehne, R.D. Matthews and D.E. Wilson, The significance of intermediate hydrocarbons during wall quench of propane flames, Twenty-first Symposium (International) on Combustion, The Combustion Institute, Pittsburgh (1986) 481-489.