I
MURRAY
H. EDSON
Esso Research and Engineering
__ .
The Otto Cycle Engine
-
Co.,Linden, N. J. ...
,
...
A Mathematical Model for
COMBUSTION Rapid and detailed calculations of thermodynamic conditions prevailing during combustion in engines and other devices can be made, using these procedures
T H E EQUILIBRIUMthermodynamic
properties and composition of the working fluid are fundamental to a knowledge of combustion phenomena, and many methods for computing this information have been published. A useful procedure, applicable with any combustible system containing carbon, hydrogen, oxygen, and nitrogen, is due to Brinkley and others (7-6,9). Only a few methods have been published for computing this type of information in a time-dependent system-e.g., an Otto cycle engine. Hottel and others ( 8 ) and Lewis and others (70)have reported methods for computing temperature gradients for combustion occurring in a closed vessel.
With this method
c
I
' I I/ h, I
I
V
Hottel, advocating the use of thermodynamic property charts, also has treated combustion in the Otto cycle engine with piston motion. However, for the Otto cycle engine, the methods reported are lacking in generality and flexibility, do not yield all the information required, and do not lend themselves to repetitive computation, especially when the rate atwhich the unburnedmass is consumed depends on thermodynamic properties and composition of the working fluid. I n the Otto cycle engine, the thermodynamic properties and chemical composition of the working fluid change continuously with piston position. Recause
of this, the testing of any one relation put forth to describe the sequential burning of the unburned mass in a dynamic system is so complex as to be insurmountable except with the aid of a computer. I n this report, a new and improved technique is described, suitable for calculating with a digital computer the equilibrium conditions prevailing during the combustion process in the Otto cycle engine. The combustion process in the Otto cycle engine is envisioned as a sequence of i stages each of which contains a combustion step, a piston movement, and a mixing operation. The first stage starts at any given time (crank-angle position)
...
Any postulated relation to describe sequential burning of the unburned mass can be tested, even when burning depends on thermodynamic properties and composition of the working fluid. Further, combustion in other devices or containers can be mathematically simulated; for example, other engines or cycles, bombs, or underground combustion. Engine thermal efficiency and behavior such as detonation or knock can be studied. The model, programmed in a number of subroutines, is being used successfully with an IBM
704 digital computer, and cornbustion in an Otto cycle engine has been simulated, using assumptions regarding progressive burning of the unburned charge. About 45 seconds of computer time is needed for a single stage. The subroutines have been used successfully to determine unburned gas temperature and system pressure, together with other thermodynamic properties, in a spherical vessel with central ignition, using previously reported data ( 7 2) for initial conditions and burning velocity.
VOL. 52, NO. 12
0
DECEMBER 1960
1007
and repeats until all of the unburned mass contained in the system is consumed. The small segment of mass to be transformed from the unburned to the burned condition in a given increment of time is called the flame element. Before the first small segment of mass is burned, the system consists only of unburned material. After combustion of the first segment of mass, the system consists of burned and unburned fractions. As the process continues, the mass of unburned fraction decreases and the mass of burned fraction increases. In this study, a combustion step is defined as the constant pressue-burning of a small segment of the total mass contained in the combustion chamber followed by an isentropic compression of this segment and unburned and previously burned fractions to a final pressure such that the total energy and volume of the system are equal to the total energy and volume before the combustion step was initiated. In the Otto cycle engine the combustion step is followed by movement of the piston to the next crank angle position. This movement isentropically compresses (or expands) the burned and unburned fractions and flame element in the combustion chamber to a new volume determined by the physical characteristics of the engine. After compression (or expansion), the flame element and burned gas fraction are mixed into one
burned gas fraction. Each stage consisting of a combustion step, piston movement, and mixing is repeated until all unburned mass is consumed. Before the first combustion step is initiated, the composition and thermodynamic properties of the system are known. For the general case of a system consisting of a flame element and unburned and burned fractions, the mass of the total system at the ith stage of the process is given by : 11' =
WilU
+ wi,/+ i-1 Wn.b
i
1, 2, . . . . ( I )
n=O
w h e r e w ~= , ~0. At any ith stage piston position, 0, the mass, w ~ , is~ ,burned sure and then it and
corresponding to ith flame element at constant presthe burned mass,
i-1
=
M',,,~,
and unburned mass,
n=O
tropically compressed to a final pressure such that the total energy after the combustion step is equal to the total energy before the step began. The equilibrium conditions of each of the subsystems present are to be determined. Let T ( P , h, X)represent the temperature for given values of the pressure, spe-
The Mathematical Model In the model presented, each of the operations comprising a single stage (combustion step, piston movement, and mixing) is computed f o r the time interval corresponding to one crank-angle degree of crankshaft rotation. In addition, the following assumptions have been made:
b There i s no heat flow or mixing between adjacent elements of gas during combustion, compression, or expansion.
ing the flame element and burned fraction. Composition i s frozen at or below 1600" K.
N o pressure gradients exist in the combustion chamber.
b The mass of the system i s constant.
b No reactions Occur in the unburned fraction.
There are no external heat transfer effects.
b
Transport phenomena-i.e., viscosity, thermal conductivity, and diffusion-have been neglected.
At
temperatures
above
1680" K., the l a w of mass action is obeyed with respect to the chemical species compris-
1008
b
The equation of state, PV = nRT, holds.
b
INDUSTRIAL AND ENGINEERING CHEMISTRY
cific enthalpy, and the equilibrium concentration of the chemical species in the subsystem. Let X ( T , P, C) represent the equilibrium concentration of the chemical species in the subsystem for given values of the temperature, pressure, and the gross composition (gramatom fraction of each element available to a specified subsystem). Because the gram-atom fraction of each element available to any subsystem contained in the total system is invariant, temperature and pressure become the only independent variables governing the equilibrium concentration of the species. Thus, for any subsystem, the equilibrium concentration of the species may be written as X ( T , P). The equilibrium concentration is best obtained by the method developed by Brinkley and others. Let h ( T , P, X) represent the specific enthalpy for given values of the temperature, pressure, and concentration. Remembering that concentration also depends on T and P, specific enthalpy may be written as h ( T , P). Consider first the constant pressure burning of the subsystem consisting of the ith flame element. For combustion to chemical equilibrium at constant pressure, h(Tj,P)
- h(T,,P)
= 0
(2)
Specific enthalpy is best evaluated using Brinkley's method (4) which leads directly to the evaluation of a set of partial derivatives essential in computations to be described later. The solution of Equation 2 is quickly carried out by the Sewton Raphson (7) iterative method in which the specific enthalpy of the flame element at any temperature and constant pressure is approximated by the first few terms of a Taylor series expansion of the form: h(Tj,P)
+ (Tf h~(Ti,d'+ ) ... . . .
h(Tj,k,P)
- T/,ic)
(3)
I n view of Equation 2, and letting LIT^,^ = T,- T,,k,Equation 3 may be rewritten :
also T ~ , L +=I T I , , ,
+
ATi,k
(5)
The procedure consists of evaluating k P to determine the Equation 4 at T J , and correction, A T,,ii, which by Equation 5 is to be added to the kth estimate, T f , k ,to give the next best estimate, Ti,,+,, of T,. The process is continued until LIT^,^ approaches zero a t which time T,,L+l approaches T,. The specific heat, Cn=hr, is best evaluated using Brinkley's method, taking into account the concentration dependence on temperature and pressure. The specification of temperature, pressure, and composition of the flame element provides a means for obtaining the
M O D E L F O R COMBUSTION specific entropy, s ( T , P),and internal energy, e(T, P ) . T h e method for calculating specific entropy is given by Brinkley and corresponds to the procedure for calculating chemical equilibrium composition. For ideal gas mixtures, specific internal energy is given by e = h - - RT M
where the mean molecular weight is obtained in the process of calculating equilibrium composition. The density is computed from Equation 7. MP RT
p a -
(7)
Thus, after the constant pressure-burning of any ith flame element, the thermodynamic properties and composition of the flame element mass are known. From this, the remaining portion of the combustion step can be computed. The remaining portion consists of the simultaneous isentropic compression of all subsystems to a temperature and pressure such that the internal energy of the system is equal to the internal energy of the system before the constant pressure burning of the ith flame element. In a n isentropic process : sj(
i
T’, P’) -
sj(
given values of each subsystem temperature and pressure. Equations 8 and 9 must be satisfied for the isentropic process in which E‘S is specified. The solution of these equations completes a combustion step a t the ith stage. Equations 8 and 10 must be satisfied for the isentropic process in which V,+l is specified. The solution of these equations completes a combustion step with piston motion a t the ith stage. T h e simultaneous solution of either set of equations is quickly carried out by the Newton Raphson (7) iterative method in which for each subsystem the specific entropy, specific internal energy, and specific volume are approximated by the first few terms of a Taylor series expansion. The expansions are given by Equations 11, 12, and 13 for the flame element and with proper subscripts apply to the other subsystems.
+
4 T j , P) = ~ ( T f , k Pk) , ( T f - T f . k ) S T ( T ~Pk) ,~, ( P - P k ) S P ( T f , k , PA) ....
+
+
e ( T / , P) = e ( T f . k , Pk) ( T j - T j . k ) er(Tf.k,Pk) ( P - Pdep( T f , k , P k ) . .
+
4 T f , P) = V ( T/,ro P k )
T, P ) = 0
(Tf-
i j = b,f,
Here E’< represents the system internal energy before the ith flame element was burned and .ei,! represents the specific internal energy of each j subsystem at the ith stage for given values of each subsystem temperature and pressure. In the Otto cycle engine, movement of the piston to the next crank-angle position also is an isentropic process. Here, however, the process can be either expansion or compression. For an isentropic process at the ith stage in which the total volume V,+l of the system, after the combustion step is specified:
(11)
..
(12)
+
a =
(*) b log M b log T
p
=T
G
(5)
(21)
The procedure consists of solving Equation 14 a t T , and Pk to determine the corrections A T , and AP,, (APk = P Pk),which by Equations 19 and 20 are to be added to the kth estimate of T , and Pk to give the next best estimates, T , k + l and Pk+l of T , and P. The process is continued until A T , and APk approach zero a t which time T , , k + Iapproaches T , and P k f 1 approaches P. It is important to note that computation of the partial derivatives (bG/ b T ) p and (bG/bP)T are an integral part of the systematized procedure for computing the composition and thermodynamic properties, and that their evaluation permits the calculation of two additional important propertiesnamely, acoustical velocity and specific heat ratio.
-
-k
+ ... .
( 13)
In Equation 13, the Taylor series expansion for v ( T , P, M ) is written as v(T, P) because A4 in Equation 7 depends on concentration which in turn depends on T and P. Using matrix notation, the set of Equations { 8, 9 is written:
1
where
CP
T
(15)
-R(1 - a ) MP
(16)
R (1 - a ) --
(17)
ST
eT
sp
=
=
c p
=
M
+ @ ep = a P
T ~ . k + l=
pk+l
where v , , ~represents the specific volume of each j subsystem at the ith stage for
.
+
T f , k ) U ~ ( T f . k ,Pk)
(P - PA)u~( T/,k, Pk)
(8)
The specific entropy evaluated at T’ and P’ refers to the entropy value of each j subsystem before the isentropic process and the specific entropy evaluated at T and P refers to the entropy value of each j subsystem after the isentropic process. For an isentropic process at the ith stage in which the initial internal energy is specified :
+
ith stage of the process, the partial derivatives are evaluated at T j , k ,the kth estimate of temperature of subsystem j , and Pk the kth estimate of pressure of the system. Furthermore, in Equations 16, 17 and 18,
=
Ti.k pk
+ ATivk + APk
(18) (19 ) (20)
and where the summations extend over each subsystem ( j = f,b, and u ) ; a t the
Thus, after the combustion step a t the ith stage, the thermodynamic properties and composition of all subsystems are known. For the Otto cycle engine with piston motion, the combustion step is succeeded by an isentropic process in which the final volume, V , + is specified. For this case, the isentropic process calculations pertaining to the combustion step are combined with the isentropic process calculations pertaining to piston motion. Thus, the solution o f t h e set of equations IS, 10) completes a combustion step with piston motion a t the ith stage of the combustion process.
- -
VOL. 52, NO. 12
DECEMBER 1960
1009
Using matrix notation? the set of Equations { S , 10) is written:
iant. These values are easily computed from the known values of the specific
CR = compression ratio
D d
= displacement, cu. cm. = connecting rod length, cm.
E
=
e
G
= specific internal energy, cal./g. = the reciprocal of the number of moles; G = l / n
g
=
h
= =
A4 where (sT)?, (s,)~? A T 7 , k 3and AP, are given by Equations 15: 16. 19: and 20, respectively, and
The procedure and nomenclature referred to in the solution of Equation 14 are also applicable in the solution of Equation 25. In Equation 23 the total system volume a t any crank angle position 8 is computed ( 7 7 ) from the physical characteristics of the engine:
internal energy, density, and mass of the flame element and burned fraction before the mixing step was initiated. The pressure of the mixed system (now consisting of only burned material) is obtained using Equation 7. The specific entropy of the mixed burned fraction is computed from the calculated values of temperature, pressure, and equilibrium concentration. The remaining isentropic process calculations comprising the second step of the mixing operation are completed on application of Equation 25 (with omission of the flame element terms). Upon completion of the mixing operation each stage consisting of a combustion step, piston (1 -
\\.here 0 5 6 5 4n and 8 = 0,2n, and 471.radians a t BDC (bottom dead center). At the completion of one combustion step with piston motion, the thermodynamic properties and composition of each subsystem present have been determined at the ith stage. Before the i+l stage is begun, however, the flame element and burned fraction are mixed and the system is reduced to one containing two subsystems. This seems reasonable since mixing does occur in real systems. Furthermore, the computations are not aimed at determining a n instantaneous temperature gradient across the burned or unburned fractions but instead are aimed primarily a t describing the flame element and burned and unburned fractions a t stepped time intervals. The mixing operation is envisioned as a two-step process: first a mixing of the flame element and burned fraction a t constant volume (density) ; second a simultaneous isentropic compression-expansion of the mixed burned fraction and unburned fraction to a uniform system pressure and a t the specified total volume. For the first step, the temperature and equilibrium concentration of the mixed flame element and burned fraction are computed applying Equation 29 and using the Newton Raphson procedure previously outlined.
I n Equation 29, the mixture specific internal energy and density are invar-
1 010
f Sin2B )’”] 1 + -) 1
(28)
movement, and mixing is repeated and the repetition continued until all unburned mass is consumed. The computation procedure described assumed. a priori, a knowledge of the mass of flame element burned during the time interval corresponding to one crank angle degree of engine rotation. This mass may be computed from any postulated relation put forth to describe the rate a t which the unburned mass is consumed. The time interval corresponding ro one crank-angle degree is given by Equation 30 a t an engine speed of R revolutions per minute. t =
106 6R
- microsec.
Acknowledgment
The author is grateful to Esso Research and Engineering Co. for permission to publish this paper, and to J. F. Kunc, L. H. S. Roblee, Jr., H. L. Yowell, and others who contributed their help and gave valuable suggestions in preparing the manuscript.
Nomenclature
acoustical velocity, cm./sec. composition, gram-atom fraction of each element available to the system C, = specific heat at constant pressure, hT cal./g. K . ; C, C , = specific heat a t constant volume, cal. ’8. “ K . ; Cv E ( b e h T ) ,
a
=
C
= gross
INDUSTRIAL A N D ENGINEERING CHEMISTRY
n
P
= =
R
=
r S
=
=
V
= = = =
L,
=
W zw
= =
X
=
y
=
p
=
s
T t
total internal energy, cal.
dimensional constant, 980.665 (g.mtisAcm./sec.2 ( g + d specific enthalpy, cal./g. mean molecular weight moles pressure, 1 atm. = 1.03328 kg./ sq. cm. gas constant, 82.0567 cc. atm./ mole K . = 1.98719 cal./g. mole O K. cranking radius, cm. total entropy, cal./O K. specific entropy: cal./g. O K. temperature, K . time, microsec. total volume, cu. cm. specific volume, cu. cm.,’g.; L: = 1i P total mass, g. mass of subsystem, g. equilibrium concentration of chemical species, mole fraction C,,’C, density, g./’cu. cm.
Subscripts
b
= burned
f
= =
3 k
=
z
=
P
=
T
=
u
=
0
=
flame subsystem kth estimate stage partial derivative of s. e , h, or u with respect to pressure holding temperature constant partial derivative of s, e , h , or t with respect to tempei atur e holding pressure constant unburned crank angle position, radians
Literature Cited
(1) Brinkley: S. R., Jr., “Combustion Processes, p. 64, Princeton University Press, Princeton, N. J., 1956. (2) Brinkley, S. R., Jr., J . C h m . Phys. 14, 563, 686 (1946). (3) Zbid., 15, 197 (1947). (4) Brinkley, S. R., Jr., IND.ENG. CHEM. 43, 133, 2471 (1951). (5) Brinkley, S. R., Jr., Lewis, B., Bureau of Mines, Research Invest. 4806, 1952. (6) Brinkley, S. R., Jr., Smith, R. W.. Proc. Sem. of Sci. Comp., IBM, 5 8 , 1949. (7) Hildebrand, F. B., “Introduction to Numerical Analysis,” p. 447, McGraw Hill, New York, 1956. (8) Hottel. H. C., Eberhardt, J. E., C h m . l i e u . 21, 439 (1937). (9) Kandiner, H. J., Brinkley, S. R., Jr., IND.ENG.CHEM.42, 850, 1526 (1950). (10) Lewis, B., von Elbe, G., “Combustion, Flames, and Explosions,” p. 651, Academic Press, New York, 1951. (11) Lichty, L. C., “Internal Combustion Engines,” p. 480, McGraw Hill, New York, 1951. (12) Manton, J., von Elbe, G., Lewis, B., 4th Symposium (International) on Combustion, p. 358, T h e Williams and Wilkins Co., Baltimore, Md., 1953.
RECEIVED for review April 19, 1960 ACCEPTEDAugust 8, 1960
