THERMODYNAMIC PROPERTIES FOR ENGINE COMBUSTION SIMULATIONS Yun X'ao Department of Mechanical Engineering and Applied Mechanics University of Michigan, Ann Arbor July 9, 1989 UM-.IEAM-90-07 The properties of gases in combustion environments are complex issues dealing with chemical kinetics and equilibrium. A large number of species are present and the detailed chemistry is not always known particularly for fuel rich mixtures. The following discussion presents a suitable approach to model compositions and properties of hydrocarbon fuel air mixtures and their combustion products for engine combustion simulations. The most accurate and detailed analysis up to now is to use a full thermodynamic equilibrium model including a huge number of species and reactions. The NASA equilibrium programs [4] and [11] are readily available for this purpose and are well documented. But the complexities and the substantial computer time involved in such models make it prohibitive to use them in engine simulation models, where computational time limits are critical. On the other hand, even those detailed equilibrium programs cannot predict some kinetics controlled engine combustion phenomena such as knock, particulate formation and flame speed with a reasonable accuracy. For such processes a simpler equilibrium model coupled with the submodels describing those particular kinetics is preferred which presents a compromise between accuracy and simplicity. The properties of premixed fuel air mixtures needed for gasoline engine simulations have been computed within a limited equivalence ratio range using simpler models as shown by Hires [7] and Olikara [10]. The simplest model is to assume a complete combustion for lean mixtures and include the consideration of equilibrium of the water gas reaction for rich mixtures. At higher temperatures effects of other dissociation reactions become important and they must therefore be included. Non-homogeneous mixtures in spray combustion processes could have a wide range of fuel air ratio variations. It may change from pure vaporized fuel to pure air from one location to another. There may be over rich or lean mixtures beyond the combustible limits and partial burn is sometimes encountered. Soot formation often occurs in fuel rich regions. In order to cover all those possibilities, the present model extends the fuel 1

Xc i\R Q-0 -3 9-

air ratio to the full range of 0 < q < oo in contrast with the previous models [7] and [10]. The model also includes the consideration of combustion progress limit and allows partially burned mixtures to exist. These are done by a simplified model without going into detailed complex chemistry for fuel rich mixtures. It also includes the possibility to take into account of soot formation effects on the properties. 1 Hydrocarbon Air Combustions Mixtures of hydrocarbon fuels and air are considered as they are applicable for many combustion processes. The general reactant mixture can be described on the basis of one mole of oxygen e ( CxHyOzNw + 02 + 1N2 where the coefficients are Vb The molar N:O ratio of the air. X Equivalence ratio, (mf /ma)/(mf /ma)o. E Stoichiometric mole fuel oxygen ratio. X The number of C atoms in fuel molecules. Y The number of H atoms in fuel molecules. Z The number of O atoms in fuel molecules. W The number of N atoms in fuel molecules. A limited number of species are considered for combustion products. According to the computation results [3] using the NASA program [11], the important species presented in products are listed as 1 2 3 4 5 6 7 8 9 10 11 C02, H20, CO, H2, 02, N2, OH, NO, O, H, Fuel. provided the equivalence ratio q < 3. At low temperatures the list can be simplified further to neglect OH, NO, 0, H radicals since they would have significant effects only when T > 2000[1K]. The list has been kept as small as possible in order to simplify the computation. A more comprehensive analysis must be done to include other radicals, soot and cracked fuels. Assume a complete combustion for a stoichiometric mixture (q = 1), the reaction is then can be described as E CxHyOZNW + 02 + bN2 -- nlC02 + n2H20 + n6N2 (1) 2

The mole fuel oxygen ratio E can be determined from the balances of C, H, O atoms 4(2) 4X+Y-2Z (2) The sections below discuss the models for low temperature range and high temperature range separately. Comparisons between them are also made to show their differences and the applicable range for each method. 2 Low Temperature Range 2.1 Chemistry at low temperature At low temperature ( T < 2000[K] ), chemical dissociation effects on the properties are small, only the. water gas reaction will be considered. Combustion processes could be described as e ~ CxHyOzNw + 02 + bN2 nlCO2 + n2H20 + n3CO + n4H2 + n502 + n6N2 + (1 - B)ecqCxHyOzNw (3) where ni Moles of species i per mole of 02 reactant. B The burned mass fraction of the fuel, or the progress variable. Consider the burned and unburned gases separately, ni can be determined by using following assumptions * For lean mixtures (0 < b < 1), CO,H2 can be neglected. * For rich mixtures (1 < _< Cmax), 02 can be neglected. and n3 = C will be determined from the water gas reaction. * For over-rich mixtures (q > qmax), CO2 and H20 burn to form CO, 12, extra fuel is left without burning. The results for unburned gas B = 0 and burned gas B = 1 are listed in Table 1 For a mixture of the unburned and burned gases (0 < B < 1), the mole fraction yi can be expressed using the mass fraction xi as =i = i/M (4) ZEx/Mi where the Mi is the molecular weight of ith species. 3

Species B = 1 B - 0O Ii 0_ < 1 1 < O_< OmOma < X O < 0 < oo O02 OeX erX - C 0 0 H20 2 COY - Y- 2(O- -1)+C 0 0 CO 0 C 2x 0 X-Z H2 0 2( - 1) - C 0 02 1 - 0 0 1 N2 ~+ 2W +-+ ~,W + x-z,/ Fuel 0 0 ex-z2 c Total (X+ lY+_ W)ec (X+ Y- 1W)EqX 2X+Y+W-2 + 2 T2 2 2 ~r JCYI X-J nT +(1- )+1) + b +| +c + l Table 1: Number of Moles ni for Mixtures with B = 1 and B = 0 Let us consider a mixture composed of two components: unburned gas'U' and burned gas'B', the mole fraction of each component can be expressed as a function of burned fuel mass fraction B B 1-B YB= 1-B B YU 1B B (5) MU MB U MB The mole fraction of each species can be computed as Yi = Yi Yu + Yi Y= (1-B)(niu) + B(nig)=ni (6) (1- B)(niu) + B(niB) nT It can also be derived, from another point of view, by dividing the total reactant per mole of 02 in equation 3 into two parts. Let the amount B(eqCxHyOzNw + 02 + bN2) of reactant be completely burned, the remaining (1 - B)(EqCxHyOzNw + 02 + 4N2) stay completely unburned, then mix them as two non-reacting components. The water gas reaction mentioned above is the reaction C02 + H2 CO + H20 (7) 4

Species 0 _< < 1 1 < < q5 mas | qmax < I CO2 ccXB (eOX - C)B 0 H20 cOYB [-Iq5Y-2( - 1)+C]B 0 CO 0 CB 2XB H2 O [2(o - 1)- C]B YB x-z 02 1-B 1 - B 1-B N2 0 + ~cOWB + WB + Fuel e(1 — B) cq(1-B) 2BX-Z Total (X+ 1Y+ 1W)e4B (X+.Y+ ~W)e4B 2X+Y+W-2 B nT +1+[e-(1+E)B]q+4 +(1-B)(1+eq)++b +1- B + 1 E + | Table 2: Number of Moles ni for a Mixture with 0 < B < 1 which is assumed to be in equilibrium at the given T. The equilibrium constant K is a function of temperature and it is approximated to fit the JANAF Table [9] as 1000 K(T) = exp(2.743 - 1.761Z - 1.611Z2 + 0.2803Z3); Z 1000 (8) The relation between the equilibrium constant and partial pressures of species (T) PCOPH2 _ C[ EY - 2(q - 1) + C] PC02 PH2 (eqX - C)[2(q5 - 1) - C] gives a quadratic equation for C (1 - K)C2 + {-COY - 2(0 - 1) + [2(0 - 1) + cEX]K}C - 2,EqX(o - 1)K = 0 (10) Introducing following transformations for convenience 6 =.EqY - 2($ - 1) a= 2(q5- 1) + CX a = 2(1- K) (11) = + aKi 7 = 2cEqX(q- 1) 5

equation 10 can be simplified as C2 + pC- -yK = 0 (12) 2 To be physically meaningful, C is solved by taking only the positive root as C = +2 +2a7K (13) In order to determine the maximum equivalence ratio qmax in the products, the physical constraint that the number of moles for each species cannot be negative is applied, which results following conditions for C C _< EX [CO2] C > 2(q$-1)-~ 2cY [H20] C>[GO] (14) C > O [CO] C < 2(k-1) [H2] The results for octane C8H18, which has c = 0.08, are plotted in the Figure 1. The possible C - q solutions have to be located within the shaded area. When q approaches qmax, C has to be under the [CO2] line and above the [H20] line. In the limit, we have EpmaxX = Cmax = (2X - Z)E4max - 2 (15) which gives 2 2X ma= (X ); Ca= X (16) E(X - Z) X - z When q >_ qmax, the reaction can be described as EcmaxCXHYOZNW + 02 +,N2 + E(4 - qmax)CxHy~OzNw CmaxCO + [2(rmax - 1)- Cmax]H2 + (ib + cEOWBmax)N2 (17) +(1- Bmax,)ECxHyOzNw where the underlines denote the extra fuel left over. Noting the analysis also imposes a constraint for B at the large fuel-air ratio range 1 ~ < Omax Bmax =- maa (18)

OCTANE 3[ H2 ]/ [ HO ]i Cm.1 [co2] 2- ------ - Temperature T = 600 K T = 750 K 1-:!;.~ A,,_./ _ _ I T = 1000 K ~ r ~/ /"~ I T = 600- - - [o ]' max 0 0 1 2 3 4 Figure 1: Physical Constraint for C and q. 0.8O 1 2 3 4 5 a 7 8 9 10 Figure 2: Valid B Range 7

B B < Bmax B Bmax B > Bma (19) The valid B range is shown as the shaded area on Figure 2. 2.2 Partial derivatives of composition at low temperature The partial derivative of ~ is nonzero when 1 ~< k < qmax, where it is OT = TAB[i] ( )B (20) with TAB = [ -1, 1, 1,-1, 0, 0, 0] The partial derivative of C with respect to T can be obtained from equation 8 and equation 13 using the chain rule dC OC dK (21) AT OK dT The first partial derivative on the right hand side is dC C2 - aC + (22) OK ac + / noting d =-2 d=, d = 0 and equation 12 have been used to simplify the result. Getting the second part is straight forward dK = K(-1.761 - 3.222Z + 0.8409Z2)d (23) dT dT And also using equation 13 and d = 0, 2 = (2+EX)K-(2-E1Y), -- 2EX(2q0- 1) OC 2KeX(24 - 1) - [(2 + eX)K - (2 - ~cY)]C (24) The other partial derivatives needed in the computation are listed in Tables 3 and 4. Obviously, these tables are obtained by operating on Table 2. The derivative of B with respect to q in Table 3 can be derived from equations 18 and 19 dB _ [ +2 B > Bmax and 0 > Xmax -l 0 otherwise (25) 8

Species 0 < < 1 1 _< < ~ max | max < _ C02 |XB (eX- C)B 0 CO 0 a2XB (X - B H2 0 (2 - ac)B a 02 -B 0 aB N2 WB WB x-z ) Fuel c(1 - B) c(1- B) -(x ) anT (X + -Y + ~W)eB (X + ~Y + -W)cB (2x+Y+W-2).B 2 2 ~ "' 2 2 X-Z |0 | + - (1 + e)B +E(1 - B) + e - Table 3: Partial Derivative of n with respect to 4 Species 0 < 0< 1 1 -< < _maz | max < q CO2 cEX c4X - C 0 H20 ~ 1cY icY- 2(4- 1)+C 0 CO 0 C 2X H2 0 2(~-1)-C 02 -1 -1 N2 l"EOirr',zcW w N2 2EqSW ~EcqW x-Z Fuel - kE -c4 -xz OnT (X+ 2Y+ LW)c4 (1 ( + (X+Y+W)EX 2XY+W-2 _ B 4: Partial Deritive of n with respect to B Table 4: Partial Derivative of n with respect to B 9

2.3 Basic thermodynamic properties at low temperature Combustible mixtures and their products inside engine combustion chambers could be treated as an ideal gas. On a mole basis (denoted by an over-bar) R- 8.31434 [kJ/kmol K] (26) h = (yihi) [kJ/kmol] (27) s = (yisi) [kJ/kmol K] (28) On a mass basis R = R ( )R [kJ/kg K] (29) ~MT~ 1 h = = -:(ni ) [kJ/kg] (30) S = M Z(ni:i) = - [ni( -rln )] [kJ/kgK] (31) where mT is the total mass per mole of 02. mT = EqCMf + M02 + OMN2 (32) It is a conserved quantity for both reactant and product and is a function of d only for given fuel and air compositions. It has a derivative 9MT -= Mf (33) The molecular weight of a fuel in above equations is Mf = 12.01X + 1.008Y + 16.00Z + 14.007W (34) 10

2.4 Partial derivatives of properties at low temperature From expressions in section 2.3 it is easy to show 1 9R - o R aT 1 MR - R P =1 AR 1 IfT ( _ 1 amT) (35) Rt8 +fT\ 7X / mT a3 1AR - 1 1 (()T R B - nT B B ah _ [( 1 ( (31 1 aT - mT aT mT aC[n oh [) and ah 1 E[( dan; h r (amTT ) n~y1(36) T = mT S( Po) m S _ )[ -l- m as1a as {( )[, o _ Rln(yi o)]} T -- mT aB The last terms in equations 36.1 and 37.1 are related to the frozen specific heat which will be given in equation 84.1. The first terms in those equations are are nonzero only when 1 < ~ < oma,, and'9ni are given in equation 20 The other required partial derivatives can be found in section 2.2. Proceeding up to now, we have had all necessary information to calculate the required thermodynamic properties and their derivatives at low temperature for computations. The calculated results from the low temperature model are given later together with the results from the high temperature model and comparisons are made there. 3 High Temperature Range 3.1 Chemistry at high temperature At high temperature (T > 2000[K] ), chemical dissociation effects have to be taken into account. The basic approach is based on the method developed by Olikara and Borman [10] where a rapid means for thermodynamic property calculation is presented compared to the extensive NASA program [4] in which a large number of species are considered. In order to simplify, the matter further, only the species of importance because of dissociation are included. These species are OH, NO, O, H (see [3]), provided 11

the equivalence ratio q < 3, beyond which solid carbon will be formed in the products. Thus the species list in the program is terminated at i = 11. The hydrocarbon-air combustion process described in equation 3 could be expanded as ecfCxHyOzNw + 02 + bN2 - n1C02 + n2H20 + n3CO + n4H2 + n502 + n6N2 + n70H + n8NO + n90 + nloH +(1 - B)EeCxHyOzNw (38) where the nomenclature is the same as before and the stoichiometric fuel oxygen mole ratio e is given in equation 2. Atom balancing yields C: BcX = (Y1 + Y3)nT H: BY = (2y2 + 2y4 + y7 + Ylo)nT 0: 2 + BecZ = (2y + Y2 + Y3+ 2y + y7 + ys+ y)nT N: 2'2 + BeqW = (2y6 + ys8)nT Also the constraint that the mole fractions of all species add up to unity must be satisfied 10 (B - 1)c= (y - 1)nT (40) i=l In order to solve for the 11 unknowns ( nT, yi,yi = 1,...,10 ), 6 more equations are needed which may be derived from the consideration of equilibriums among the products. Following 6 major dissociation reactions are considered in the program 2 H2 = H r=, 2~2=~ K /2 1 1/2 Y4 102 - 0 K2 - 1/2. Y5 2N2 + 202 NO K3 = 1 (41) Y6 Y5 (41) 1H2 + 1202 = OH KI4 =,YT, Y4 Y5 H2 + 102 = H20 K5 = 1/2 2Y4Y5pl/2 CO + 02 C02 y6 1/2 Equilibrium constants in equation 41 are curve fitted to the JANAF Thermochemical Table [9] in the form T B2 log i = Ain00 + T + Ci + DiT + EiT2 (42) 12

where T is in [K], P in equation 41 is in [atm]. The coefficients Ai, Bi, Ci, Di, and Ei are listed in Table 7 for the temperature range 300 < T < 5000K. Rearranging the expressions. for equilibrium constants, the mole fractions of the species i = 1,2,7,8,9,10 can be expressed as functions of the mole fractions of the species i = 3, 4, 5, 6 Y1 = ClY3Y/Y C1 = IC6P2 Y2 = C2Y4Y5 C2 = IK 5P,, 1/2 1/2 C3 = K4 Y7 -= 03Y4 Y5 C3 - 4() (43) Ys = C4y6 Y5 C4 = 13 1/2 Y9 = C5Y5 C5 = 0K2P- Ylo = C6y41 C6 = K1P-1 fining D = Bcq d, = 2b+ DW d2 =2+DZ (44) d3 = DX d4 = 1 + B(X - 1) d5 = BX equation 39 and equation 40 can be manipulated into a homogeneous form so that the nT term is cancelled out Y(yl + y3)- X(2y2 + 2y4 + 7 + Y10o) = 0 dl(2yl + Y2 + y3 + 2y5 + y7 + Ys + y) - d2(2y6 + Ys) = 0 dl(yl + y3) - d3(2y6 + Y8) = 0 d4(y1 + Y3) + d5(y2 + y4 + Y5 + Y6 + Y7 + Y8 + Y9 + Y10 - 1) = 0 Utilizing equation 43, equation 45 can be expressed symbolically as functions of y3, Y4, Y5, Y6 only fj(Y3, Y4, Y5s, Y6) 0; (j = 1,2,3,4) (46) If we have an approximate solution [y3n, y4, y5, y6], where n denotes the approximation after the nth correction, reasonably close to the real solution of equation 46, a further correction can be made by yl+l = y + Ayn, (i = 3,4,5,6) (47) until the correction terms are smaller than a prescribed error allowance, which is the so called Newton Raphson iteration scheme. 13

The correction terms could be generated by a Taylor expansion of equation 46 about a known vector. Neglecting the second and higher order partial derivatives we get a set of linear equations for Ay Ofj a fj a fj Of y fj + ( >yY + ( )/Y4 + (Y )Ay5 + (AY6 oY, 0, (j = 1, 2, 3,4) (48) Or writing in matrix notation [A] [Y] = [F] (49) where [F] = [-f], [Y] = [Ay] and [A] = []. For convenience, defining following partial derivatives Oyi Dij - Oy (i= 1,2, 7, 8,9, 10; j 3, 4,5,6) ay3 so that D13 = Cly/2 D15 = =ClY3Y/ 1_2 -1/2 D24 = C2Y D25 = 1C2y4 /2 D74 = C31/ 1/2 75 = I3y 1/2y -1/2 (50) D6- 4Y6 Y D = C Ys Elements in the [A] matrix can then be expressed as All = Y(D13 + 1) A12 = -X(2D24 + 2 + D74 + D104) A13 = YD15 - X(2D25 + D75) A14 = 0 A21 = dl(2D13 + 1) A22 = d1(D24 + D74) A23 = di(2D15 + D25 + 2 + D75 + D85 + D95) - d2D85 A24 = dlD6 - d2(2 + D86) 1) A31 = dl(D13 + 1) A32 = 0 A33 = d1D15 - d3D85 A34 = -d3(2 + D86) A41 = d4(D13 + 1) A42 = d5(D24 + 1 + D74 + D104) A43 = d4D15 + d5(D25 + 1 + D75 + D85 + D95) A44 = d5(1 + Ds6) 14

The solution of the matrix equation 49 by means of the Gaussian elimination gives the Ay required in equation 47. After the independent variables Y3,Y4, Y5,Y6 in the iteration are corrected, the dependent variables are updated using equation 43 at every step until the convergence criteria are satisfied. 3.2 Initial estimate It is necessary to make an initial estimation for the iteration described above since it is not self starting. A low temperature model, where only the 6 species C02, H20, CO, H2, 02, N2 are considered, is found to give a remarkably good initial estimation. The simplified equation reads BeqbCxHyOzNw + 02 + 4N2 -- nlC02 + n2H20 + n3CO + n4H2 + ns502 + n6N2 (52) The C balance BecX = (Y1 + Y3)nT with equations 43 and 44 give DX - (+ C1Y5/2 )nT (53) The H balance BeVY = (2y2 + 2y4)nT with equations 43 and 44 give DY Y4 = DY (54) 2(1 + C2y/2 )nT Insert above relations into the O balance BEqZ + 2 = (2yl + Y2 + Y3 + 2y5)nT, a nonlinear equation for Y5 is derived f5 = D[X(1 + g9) + Yg2 - Z] + 2(ysnT - 1) = 0 (55) with 1 X2 Y 1 = D[-1 + 2] 3/2 + (56) 2 C 2C2 y5 so that Y5 can be iterated to a given accuracy n+l yn _ (f5)n (57) nfY5 5 (57) The function g is defined as and the total mole number nT can be taken from Table 2 as a first approximation. 51/2 15 -3 A L -. -. -I -- 1 5 - - l - -'.-I - -I 0 -, - - A - -

The mole fraction for nitrogen Y6 can be estimated as V/ + ~DW Y6= 2 (59) nT The results of computed composition are shown in Figures 3 - 6. Having the composition available, we pursue further to calculate the partial derivatives of mole fractions. 3.3 Partial derivatives of composition with respect to T and P at high temperature The independent derivatives are obtained by solution of the matrix equation that results from differentiating equation 48 Ofj Ofj Oy4 Ofj y O 5 af j Oya 6 3 ) 6 aai+f (h 0 +f(if fiY + (0ifj +\('9 Y8if = o, Y6- (j = 1,2,3,4) (60) + (y- + (Ty4) - + (jy5) +( )y0- -0 In matrix notation [A] [Y,j] = [F, ] (61) where ~ could be T,P,k, or B, [F,] = [-f ], [Y,] = [X] and [A] is identical to that used in equation 49. The elements ~ on the right hand side can be evaluated from differentiating equation 45 2,- = yL( a) - x[2 U(a~) + a-U( )] -d2( aT )U -lCl aT ) ( C2, z aT ( ) + IC( )+ -() + -C4 (aT) + C ( aT + C6 ( )](62) af3 = Y.L( a) X[322 C2 ( )+ (T a) )] PL = dl[2y L(OC) \+ UaC2 )+ Ju(ac-5)] a = dl C1 ( ap a~f = d4L() - x[r2( C) + &s(aP)+C 6(aP)] ap C1 aP C2 aP C5 a C6 aP The required partial derivatives of C with respect to T and P are 16

T =1250[ K ] T = 2250 [ K ] O 0- -o - Ha OHH Co 00-3-3H 0 00 NO 0 -4- 4- aO -5 - -5 NO 0, OH 0 -6 o o 0 2 3 0 1 23 N, ~~~~~~~~~~~~~N. Figure 3: ~~ComutoPrdc ComoiinasFctnsf I(P30[m] HgO HO NO cog OP~~~~~O -2~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ NO -2 OH H, 0 r H 00 -3 I 1 2 O o 0~N OH -4 -3 0 0 U I I CO OH 0 II NO -4 - -4 H, OH -51- - - 0 -5 -6 - 0 1 2 3 0 2 Figure 3: Combustion Product Compositions as Functions of c (P = 30 [atm]) 17

P =10 [atm] P= 50 [atm] 0-' | R, 1, - - - D H I(D 1. 0 ]:0. 1 1.0l l h -2- 3-2 O O -3- -3 -4 411 -41000 2000 3000 4000 5000 1000 2000 3000 4000 5000ooo T [K] T (D= 0.8 (1)= 1.2 SaO H HO 0 HO Figure 4: Combustion Product Equilibrium Compositions as Functions of T con C 0 NO OH NoHOH -4 -4- -4O 1000 2000 3000 4000 5000 0 1000 2000 3000 4000 5000 T [K] T [K] 18

T =15O0 [K] T = 3000 [K] 0 0 ~~~~~~~~~NN OO 02 o 0 co I I H -4$~~~~~ F....s His =0.8 = 1.2IP 30[atm]l 0 ~~~~~~~~~~~~~~N, N Hc0 Hc II II \I II I- l CO, -3 -1 N~~~~~~~~~COO Oa~~~~~~~~~ OH -4 - - I 1 -4 0 0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 F F 0 0 Iii ~ ~ ~ ~ ~ ~ ~ ~ ~~~N — 3 H~~~~~~~~~~~T = 3000 [K] = 00 -4'' 0 20 ~~4'0 6~0 8'0 10 0 0 2 0 4 0 6 0 8 0 10 0 P [atm] P [atmn] 19 H20 H20 -1- -1 cog NO H hD-2- OH hO-2 H2 ~~~~~~~~~~NOHo~~~~~~~~~~~~ H 02 IT 30 00 [K] III IT 3000 [K]l -4 20 410 610 810 160 4 210 410 6I0 810 160 P [atm] P [atm]

= 0.8, T = 1000 [K] D= 1.2, T = 1000 [K] N~ Ng2 0. 2H0 o o -4I I I I I I 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 B B = 0.8, T = 2000 [K] = 1.2, T= 2000 [K] 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 1 c og~~~~B Figure 6: Combustion Product Equilibrium Compositions as Functions of B 20

ac. _ dKs pi ac. _ c1 aT - dT P -2 P ac2: _ dK p ac _ 1 C2 aT - dT aP - 2 P aC3 _ dK4 aC3 _ o aT dT aP - aC4 _ dKE oC4 _ 0 (63) aT - dTP -a ac_ _ dK' p-~ ac5 _ _. l. aT dT PP 2 P ac& _ dK p-1 ac _ _1c aT dT aP - 2 P where dKi Ai Bi dT = 2.302585Ki( - + Di + 2EiT) (64) The remaining derivatives can be obtained in terms of the independent ones.yl _ _ ay_ yl _1 yl( )_D13(_ ) _D15(a_ ) aT = C, ( a9') +D13((93)+Dj5(aY'5) ap =c(& ap) D3( aT C1 T J aT \ aT aP - C1 P J P 5 P T —.C2,(%\)....D24(o#)_..D25(_.) __...(p )+D24(09')_.. D25(2 p') T _.C.3(aCT) +D74(aa_)+D75() D74(a )+D5( ap) aT C T T T - a aT C3 aT T) +D95('\ aT J AP \ aP J )+D95( P 2.a — 05 ( aP D86 Csj)') + D85 Ol) TaY1Y C6 ) — +D z 10 _4( {T c6 \)+D O4( aT C6 derivatives of composition with respect to P C P a nd B 3.4 Partial derivatives of composition with respect to 0 and B at high temperature The partial derivatives of the mole fractions with respect to q and B are similar, except with 9O/O, &/0B replacing 0/dT, a/aP. Noticing only the parameter D is a function of ~ but none of Ci is, the expressions are therefore much simpler 21

3fi ( + a - -, -(D )[W(2yl + Y2 + y3 + 2y5 + y7 + Y8 + Y9) - Z(2y6 + Y8)] -a - ()[W(yl + y3) - X(2y6 + Ys)] 0~ = ()(-Y - Y3- XY11) (66) aB - aB = (~)[W(2yl + Y2 + y3 + 2Ys + Y7 + Ys + y)- Z(2Y6 + Ys)] aB = (-) [W(yl + Y3) - X(2y6 + Y(s)] AB =-Y1 - Y3 - XY11 where OD OB aD Oc = E(B + ); (67) and the partial derivative B is given in equation 25, where the equivalence ratio limit qmax is set to be approximately equal to 3 to prevent the presence of negative mole fractions. The other partial derivatives of mole fractions are all D13(8~,)+ D15(8y5) Oyl = D13( 8y3 ) + D15(ay) 4 a — 15 a aB aB - =2 _ D24(9 )+ D25(2-) 2 = D24(=) + D2s(_a) a a — aB aB - B =Y7 D74(aY4) + D75('5 ) y = D74 (_ + D75 ( ) ao ao a O aB &B aB (68) ay9 _ D95s(() Yay5) aYo = D104(Oa_4 ) a, _ D104( a=..) 4 aB a3.5 Basic thermodynamic properties at high temperature The required thermodynamic properties in engine combustion simulations are R = R [kJ/kg.K] (69) M h = (yihi) [kJ/kg] (70) 22

1 s = M- (yii) [kJ/kg K] (71) The equivalent molecular weight M of the mixture in above equations is M = (yiMi) (72) and it has derivatives aM Oyi =( )Mi] (73) The entropy 7i in equation 71 should include both temperature contribution and pressure contribution under the partial pressure Pi = yiP of the species "i". -o Pi Si = S - Rln(P (74) So that equation 71 can be written as s = M [yi(si - Rln yi)] - Rln( In ) (75) 3.6 Partial derivatives of properties at high temperature The partial derivatives of properties can be derived from section 3.5 as 1 OR 1,aM ~1 iR _ 1O(8M = -(76) R O6 M dS _h h OM ) + I ay- w, + y Ohi (77) h - -( )+ -[(')i+ Yi(' A)] Os s om 1 y yi 1 Ogg R aP X - = -M( t )+M RI[( ~)(~-Rln Y i)]+ M Yi )- ( y) (78) Enthalpy hi and entropy so are only functions of T and they can be found from equation 84 as Cp and -. Therefore, the last term in equation 77 and the third term in equation 78 will vanish except when ( = T. The last term in equation 78 is finite only when = P. 23

3.7 Relations between F and q The representations of the fuel mass fraction F and the equivalence ratio q are exchangeable. Their relations can be derived from the definitions mf (m —t) F -= ma (79) ma +mr (ml )o maSmf ~ c~fm A simple manipulation gives F 1-F F. (80) 1-Fo and F Fo 0 (81) (1- Fo) + Foq where Fo is the stoichiometric fuel mass fraction EMf o eMf + Mo2 + bMN2 (82) The derivatives have a correlation d 1-Fo 1 a (83) OF Fo (1 - F)2 (83) 4 Results and Discussion The calculated results for properties using the both methods are shown in Figures 7 12. 1. Unlike other available programs which will predict a unrealistic negative species concentration when the equivalence ratio 4 > 3, the present program is capable of covering the entire fuel mass fraction range 0 < F < 1 ( 0 _< < oo ) at an expense of neglecting the solid carbon formation under fuel rich conditions (Figure 10 and 11). It preserves the original fuel structure in order to simplify the calculation. If one is interested in modeling fuel cracking or soot formation, an additional subroutine is needed to carry out the calculations of corresponding reaction equations. This capability is crucial in stochastic Diesel combustion modeling since we are often confronted with computing the properties of fuel rich elements due to the heterogeneous nature of the combustion process. Suppose a pure air element could always be found to mix with, the sequence of the resulting equivalence ratio from a pure octane fuel element would be 4 = oo, 15.03, 5.01,2.15, 1.00,'". 24

2. The dissociation effects on properties become important starting about T = 2000 [K], (Figure 7, 8) The rapid increases in h and Cp are mainly due to the shifting of composition from C02, H20 to high enthalpy containing radicals, such as OH, NO, O, H (Figure 4, 13.b). The maximum deviation of the calculated specific heat from the simple routine happens at around T = 4000[K]. Beyond which, the effect of the decrease of CO2, H20 (Figure 4), which have large specific heats (Figure 13.a), dominates, that is, the decrease of the frozen specific heat dominates, although the concentrations of radicals O, H still keep increasing. This provides us an opportunity to save some computer time by using the simple routine if one is only interested in energy balance with the exception of near the TDC range where the peak temperature is too high to neglect the dissociation effects. Also, because of the lacking of the initial conditions to start calculation, simulation programs usually run for several cycles using guessed initial conditions until some convergence criteria are reached. The convergence process can be sped up by using the simple routine in pre-convergence cycles. 3. The pressure effects are small compared with effects of the other three independent variables T, F, B. Increasing the pressure actually reduces the effects of dissociation reactions (Figure 7, 8). It helps in some degree when the simple routine is used in pre-convergence cycle calculations because high gas temperature near the TDC is always associated with the peak pressure. The simple routine is totally independent of pressure since the pressure dose not appear in the water-gas reaction equilibrium equation, and all derivatives with respect to pressure therefore vanish in the simple routine. 4. In some range of fuel rich region, increasing the progress variable B gives an endothermic reaction instead of a exothermic one (Figure 9.b, 11.a). The reason for those phenomena is that the maximum utilization of the oxidizer is allowed in the program. The definition of the B is the mass fraction of the burned fuel over the total available fuel, and there is no restriction on the available oxidizer. Therefore in a fuel rich case, small B allows fuel to find more oxidizer and the actual air fuel ratio of the products will be closer to the stoichiometry until the oxidizer is used up. This in turn gives more energy release. As the B increases, more and more fuel is put into reaction, this will cause a shift of the actual air fuel ratio towards to the richer input value, the reaction will be endothermic. If a totally unmixed model is desirable, such as in SI engine simulations, multizone modeling, the program can be called twice setting B = 1 and B = 0, then combining the results using the mixing rule according to the fractions of the both burned and unburned gases. The corresponding results are shown in Figure 12. 25

9z L JO suojound se sepnpoid uoijsnquwo3 Jo speeH DipladS pus saidIevllpu:L aInl'i [X] JL 0009 000 O 000: 0001 0001 0 0, C) -u popp. 018 aIdtur..S........... —.. _. E] L 0009 000' 000~C 000, 0001 0 8 0 = C 6 I - 5"

40l 30- Simple Method o P E < 20CD 0 = 0.8................ = 1.0 CI =_1.2 100 II 0 1000 2000 3000 4000 5000 T [K] 1.4 1.3~ ~- ~..,..,.Simple Method O, ~ \ l. 1.2eD = 1.o (D = 1.2 " "1.1 *I I 0 1000 2000 3000 4000 5000 T[K] Figure 8: M and'y as Functions of T for Combustion Products 27

X 4 - 0-4- 0 0.2 0.4 0.6 0.8 F 5 3 1IT =2000 /K], P 50 [a _ "~B = 1.0 -1 0 0.2 0.4 0.6 0.8 F Figure 9: h as Functions of F for Combustion products 28

120 100T = 2000 [K], P = 50 [atm] 80, 0 0)0 0.2 0.4 0.6 0.8 1 400.5 04 B= 1.0 T = 2000 [K], P = 50 [atm] 0.1 0.0- I 0 0.2 0.4 0.6 0.8 1 F Figure 10: M and R as Functions of F for Combustion products 29

-2 = 1.2 -4 0 0.2 0.4 0.6 0.8 B 0.31. IP=50[atm] 0.30 = 0.8................, = 1.0 0.2 0.4 0.6 0.8 B 0.31 Figure 11: h and R as Functions of B 30 0 0.20.4 06 0.

0.5 B= 1.0 0.4 IT = 2000 [K], P = 50 [atm] AB = 0.2 0.2 0 0.2 0.4 0.6 0.8 F X - 2 AB = 0.2 IT =2000 [K], P = 50 [atm] cB=1.0 0 0.2 0.4 0.6 0.8 F Figure 12: Results of the Unmixed Model 31

SPECIFIC HEAT ENTHALPY ENTROPY Cp [ kJ/kmol K] h [ kJ/klX0 k/kmol K ] CO 3~~5 ThermHo Hf 35 CO H ~15 ~-400 150 2000 3000 4000 600 0 1000 2000 3000 4000 6000 0 1000 2000 3000 4000 5000 T[K] TK] T [K] Figure 13: Curve Fitting of the JANAF Table Data 5 Thermodynamic Properties for Elements A simultaneous polynomial fitting of the specific heat Cp, enthalpy h and entropy s~ to the JANAF Table [9] data is used to evaluate those thermodynamic properties as functions of temperature at the reference pressure P, = 1 [atm]. The functions employed are -= al + a2T + a3T2 + a4T3 + a5T4 R RT = 2 a3l T2+ T + -3- + 4 (84) -a,+ T+. T as T8 RT 2 3 4 5 T 80. a3T2 + a4T3 a5 4 -R = al InT + a2T + 3T + 4T +a7 R ~2 3 4 the enthalpy h and entropy s~ are related to Cp as h = JCpdT; s-0 =J TPdT 32

SPECIFIC HEAT ENTHALPY ENTROPY c[ kJ/kmol K h kJ/kmol ]X10-3 S [ kJ/kmolK] 600 1200 1400 500- CH, 12004,~~~~~~~~~~~~oo'~"l~~~~~~~~~~~~~~~~~~~~~~~~~~~~, ool ~~l /~//I { ooo11 / / / 200 6000 CIII. 200 C Ht 100 400 0 on 0a ~~~~~~~~-200 200 0 100 1400 0100 2000 2500 0 700 1400 2100 2000 35000 0 700 1400 2100 2000 3500 T[K. TIK] T K] Figure 14: Fuel Thermodynamic Properties SPECIFIC HEAT ENTHALPY Cp [kJ/kmol K] h [kJ/kmol] X10-3 600 1200 1000 500 800 400.ool 3001 CSH 400`v oo1- 1 /' ~5~ ~. ~~Hla 2~~00 20' 200_sCH30H~~3 100- C H 4: - 0~~~~~~~~~~~~~~~~~0 1 -200 0 700 1400 2100 2800 3600 o 700 1400 2100 2800 3500 T[K] T[K] Dashed lines are from Ferguson [ ], dotted lines are from Heywood [ ]. Figure 15' Comparison between Different Polynomials 33

The coefficients are listed in Tables 5 and 6. There are two sets of coefficients to cover both low temperature range (300- 1000[K]) and high temperature range (1000-5000[K]) because it was found long ago that the thermal properties tend to have characteristic knees around 1000 [K]. Thus a single polynomial fit may not give a good reproduction of the original data, the pinned polynomials used here give relative low deviations. The results of the fitting are shown in the Figure 13. Other alternate polynomials for fuel properties are found in Ferguson [3] or Heywood [6]. Since most thermodynamic data tables cover temperature ranges too low for combustion calculations, those polynomials, when extrapolated to higher temperatures, will give improper results, especially for Cp curves (see Figure 15). To overcome this inconvenience, the coefficients in the Table 6 were calculated using the Woihoit's method [13] which forces Cp curves to approach asymptotically the correct classical upper values at Too. 6 Transport Properties A simple approach for computing transport properties for working fluids is assume that they can be represented by the values for air. Approximate correlations fitted to air data could be used T1.5 p = 1.468 x 10-6 T+ 110.4 [Pa s] (85) A = 3.65 x 10-7 T0.748 [ ] (86) m K The principal advantage of these correlations is computation speed. For such approximations the properties vary as function of temperature only, but not vary with gas composition and pressure. The measured property dependence can only be explained by more sophisticated models. A more rigorous treatment of gas transport properties based on more realistic intermolecular potential energy models can be found in Hirschfelder el at. [8], which also presents methods for computing the transport properties of mixtures of gases. The procedures used in our program to compute the transport properties are basically based on those techniques. According to the Chapman-Enskog theory [1], the viscosity and thermal conductivity of a pure monatomic gas may be written as Hi = 2.6693 x 10-6` [Pa s] (87) 34

100 150 1001 1~~o X,0 YW 50 l< 50 0 I 0 0 1000 2000 3000 0 1000 2000 3000 T [K] T [K] Figure 16: Transport Properties of Air 35

43\ C2 0 0.1 1 10 100 T* Figure 17: Integral Function for Transport Properties Ai = 8.3224 x 10-5 aiTi mKW (88) a2 [m K] where M is the molecular weight, a is the molecular diameter in [A], and 2 is the dimensionless integral of order unity dependent on the details of the force interaction during a collision. For the Lennard-Jones interaction, which assumes the potential energy' is a function of the distance r separating two molecules, 4i = 4 Ei [(ni)12 _ (1i)6], the integral function is given by Hattikudur [5] and may be curve fitted as Qi = 1.155T*-0.1462 + 0.3945e-0'6672T* + 2.05e-2 168T* (89) where T* is a dimensionless temperature scaled by the Boltzmann's constant k and the maximum energy of interaction during the collision E T kT (90) A partial list of the interaction parameters is given in Table 8. When values of a and E are not known, they may be estimated from the properties of the fluid at the critical point, 0.77T-; aa 8.41A3 3 113.75( p) 36

From equation 87 and 88, the monatomic gas relation between A and lu is 15 5 Ai = Ripi = 5CviPi (91) Although the viscosity formula was derived for monatomic gases, measurements show it is remarkably good for polyatomic gases as well. The thermal conductivity model dose not take into account the vibrational and rotational energy exchange in collisions between polyatomic molecules which contribute to energy transport in gases of interest in engines. Comparisons with experimental data show some discrepancy for polyatomic gases. Eucken [2] proposed an empirical correction expression Ai = (Cpi + 4Ri)Ii (92) The monatomic formula is a special case of equation 92 because Cp = 5R for monatomic gases. Equation 92 also gives a simple method of estimating the Prandtl number Pr=P- C7C Pr, _ (93) A Cp + 1.25R which is in good agreement for monatomic and diatomic gases but less satisfactory for more complex gases. A semiempirical formula for the viscosity of a multicomponent gas mixture was proposed by Wilke [14] n yili n Yi Ai IL el S A (94):= E..~j l Yij' Yj~ij (9) in which [1+ () 2 ( )4 ]2, ~=Mi In the case of reacting mixtures, that are mixtures in which, compositions shift with temperature and pressure, as opposed to frozen mixtures, there is also a contribution due to the dissociation and recombination as molecules fluctuate back and forth in the temperature gradient driving the heat flux. In order to take into account this effect following approximation could be used Se = P (96) CPf where Cp is the equilibrium specific heat and Cpf is the frozen specific heat. The results are given in Figures 18 and 19. 37

B=0 B=1 100 100 P=50[atm]] F0 IP = 50 [atm]l F=0 75- 75 A F = 0.2 A F = 0.2 Io..o II o X - X IX ro 50- W 5so0 0 I I 0 0 0.2 00 2000 000 0 1000 2000 3000 T [K] T [K] P = 50 [atm] P = 50,[atm] 0 0 T 1000 [K] 25-1 1 I1 225 I = 1.0 Calculated Viscosity 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 F B 38

P = 50 [atm] =1 300 1 ~Q. =0 0.9 =0 B=O0 200- D=1 Burned x X -OO [ ] A 10.80 B = -~.* ~~~0.7 1000 0.6 0 - I 0.5 0 1000 2000 3000 0 1000 2000 3000 T [K] T [K] P =50 [atm] = 400 400 B=O.......... IP= 50 [atm]l 300 B = 1 300 T3000 [K] T = 2000 [K] _ _ 2000 100- 100T ==10001~~~~~~~~~~~ = K]T1000 [K] 0 0.2 0.4 0.6 0.8 1 0 0.2 0.4 0.6 0.8 x F B Figure 19: Calculated Thermal Conductivity 39

References [1] Chapman, S. and Cowling, T. G.: "Mathematical Theory of Non-Uniform Gases", Cambridge University Press, 1951. [2] Eucken, A., Physik. Z., 14, pp 324-332, 1913. [3] Ferguson, C. R.: "Internal Combustion Engines", John Wiley & Sons, 1986. [4] Gordon, S. and Mcbrige, B. J.: "Computer Program for the Calculation of Complex Chemical Equilibrium Composition, Rocket Performance, Incident and reflected Shocks, and Chapman-Jouguet Detonations", NASA publication SP-273, 1971. [5] Hattikudur, V. R. and Thodos, G.: "Equations for the Collision Integral", Journal of Chemical Physics, 52 (8), pp 4313. [6] Heywood, J. B.: "Internal Combustion Engine Fundamentals", McGraw-Hill, 1988. [7] Hires, S. D., Ekchain, A., Heywood, J. B., Tabaczynski, R. J. and Wall, J. C., "Performance and NOx Emissions Modeling of a Jet Ignition Prechamber Stratified Charge Engine", SAE Paper 760161. [8] Bird, R. B., Stewart, W. E. and Lightfoot, E. N.: "Transport Phenomena", John Wiley & Sons, 1960. [9] "JANAF Thermochemical Tables", Second Edition, The Dow Chemical Co., Midland, Michigan, 1971. [10] Olikara, C. and Borman, G. L.: "A Computer program for Calculating Properties of Equilibrium Combustion Products with Some Applications to I.C. Engines", SAE paper 750468. [11] Svehla, R. A. and Mcbrige, B. J.: "Fortran IV Computer Program for Calculation of Thermodynamic and Transport Properties of Complex Chemical System", NASA technical notes TND-7056, 1973. [12] Van Wylen, G. J. and Sonntag, R. E.: "Fundamentals of Classical Thermodynamics", JOHN WILEY & SONS, 1986. [13] Wilhoit, R. C.: "Thermodynamics Research Center Current Data News", Vol. 3, No. 2, Texas A & M University, 1975. 14] Wilke, C. R., Journal of Chemical Physics, 18, pp 517-519, 1950. 40

A APPENDIX: COEFFICIENTS 41

_-emperature Range 300 < T < 1000 [ K ] a3 a4 a5: 6 a7 -' —-2t.83465E-05 -.48251E-08.83841E-13 -.10027E+05.44047E+01 ~-t e.57805E-05 -.12578E-07.45863E-11 -.11239E+05.14432E+02 _ L.54314E-05 -.21260E-07.92433E-11 -.13955E+05.19355E+02 -1.81007E-05 -.29215E-07.12673E-10 -.17126E+05.18349E+02 I —1t.12532E-04 -.37015E-07.15256E-10 -.20038E+05.18773E+02 5-_-t'.12596E-04 -.44284E-07.18722E- 10 -.22927E+05.20881E+02....I.21813E-04 -.54234E-07.20809E- 10 -.26003E+05.17509E+02 — 41t.13469E-04 -.58321E-07.26347E-10 -.28873E+05.21120E+02 —:9.71701E-05 -.87932E-08.23906E- 11 -.25353E+05.11233E+02 - 3.32672E-05 -.16375E-07.73522E- 11 -.29674E+05.18271E+02 (-::-I -.21823E-05 -.34845E-.07.18421E-10.89179E+04.37036E+02: demperature Range 1000 < T < 5000 [ K | a3 _a4 a5 a6 a7 — 2!-.27888E-05.40232E-09 -.20206E- 13 -.10162E+05.75069E+01 - 4.-_1 -.45573E-05.67250E-09 -.35982E-13 -.12718E+05 -.52395E+01 — il4 -.62839E-05.91794E-09 -.48124E-13 -.16465E+05 -.17844E+02 — AI -.78760E-05.11502E-08 -.60251E-13 -.20519E+05 -.32140E+02 -— 3 | -.35333E-05 -.57422E-09.15159E-12 -.25537E+05 -.63729E+02 --. 35307E-05 -.54662E- 09.14789E- 12 -.30738E+05 -.95832E+02 __-4 — -.11120E-04.17132E-08 -.96212E-13 -.33679E+05 -.94335E+02 _=i... -.12878E-04.17911E-08 -.88109E-13 -.36113E+05 -.89323E+02 -— 2 -.30503E-05.43588E-09 -.22247E-13 -.26158E+05.23782E+01: —'x_ -.35144E-05.42572E-09 -.15468E-13 -.31944E+05 -.16427E+02: —— ~-.65489E-05.98388E-09 -.53417E-13.40879E+04 -.43938E+02 _Eicients for Fuel Thermodynamic Properties 43

i A B C D E K1.432168E+00 -.112464E+05.267269E+01 -.745744E-04.242484E-08 K2.310805E+00 -.129540E+05.321779E+01 -.738336E-04.344645E-08 K3.150879E-01 -.470959E+04.646096E+00.272805E-05 -.154444E - 08 K4 -.141784E+00 -.213308E+04.853461E+00.355015E- 04 -.310227E-08 K5 -.752364E+00.124210E+05 -.260286E+01.259556E-03 -.162687E-07 K6 -.415302E-02.148627E+05 -.475746E+01.124699E-03 -.900227E-08 Table 7: Coefficients for the Equilibrium Constants Species cr [A] c/k [K] Species a [A] c/lk [K] C02 3.996 190.0 CH4 3.822 137.0 H20 2.641 809.1 C2H6 4.418 230.0 CO 3.590 110.0 C3H8 5.061 254.0 H2 2.915 38.0 C4H10 5.341 313.0 02 3.433 113.0 C5H12 5.769 345.0 N2 3.681 91.5 C6H14 5.909 413.0 OH 3.147 79.8 C8Hs1 7.451 320.0 NO 3.470 119.0 C6H6 5.270 440.0 0 3.050 106.7 H 2.708 37.0 Table 8: Collision Parameters 44

UNIVERSITY OF MICHIGAN 3 9015 02527 80141111111 3 9015 02527 8014