NOx Model
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
SOCIETY OF AUTOMOTIVE ENGINEERS, INC.
400 Commonwealth Drive, Warrendale, Pa. 15096
A Computer Program for
Calculating Properties of
Equilibrium Combustion
Products with Some
Applications to I.C. Engines
Cherian Olikara and Gary L. Borman
University of Wisconsin
Automotive Engineering Congress and Exposition
Detroit, Michigan
February 24-28, 1975
750468
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
Copyright © Society of Automotive Engineers, Inc. 1975
All rights reserved.
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
750468
A Computer Program for
Calculating Properties of
Equilibrium Combustion
Products with Some
Applications to I.C. Engines
Cherian Olikara and Gary L. Borman
University of Wisconsin
MANY MODELS FOR ENGINE COMBUSTION PROBLEMS use
the First Law of Thermodynamics as applied to
either the entire cylinder contents or to sub
systems. Examples include models for spark ig nition engines (1)*, diesel engines (2) and stratified charge engines (3). Typically the major species of the products of combustion may be assumed to follow a shifting equilibrium process for thermodynamic purposes. For pur
are given by References 6-9. The NASA-Lewis program (9) is very extensive and includes ther modynamic data for hundreds of species.
The purpose of the present program is very narrow by comparison to the NASA program. Its purpose is to specifically deal only with the gas phase products of combustion of hydrocarbon
fuels (containing, C, H, O, N atoms) and air.
The program is however extended to calculate the
partial derivatives of internal energy and molec ular weight which are helpful when numerically at the equilibrium concentration (4). Because solving the first law as a differential equation of this wide spread use it is important to have in time using the method of Reference 10. Most a rapid means of calculating the equilibrium in importantly, the program is to provide a rapid ternal energy, molecular weight and species mole means of calculation comparable to the use of fractions of the products of hydrocarbon-air com the regression analysis equations given in Ref poses of chemical kinetics calculations, many of the major species may also be assumed to be
bustion.
erence 11.
Calculation of the equilibrium composition and internal energy of combustion products for
ENERGY EQUATION FORMULATION
engine cycles goes back to the early work of
Tizard and Pye (5) and has been the subject of numerous computer studies, a sampling of which
*Numbers in parentheses designate References at end of paper
In order to treat the cylinder gas system
thermodynamically the system may be divided into
cells such that the composition and temperature
of each cell is uniform. This implies that the molar average temperature of the cell may be
ABSTRACT
bon fuel and air is described.
A subroutine is
A computer program which rapidly calculates the equilibrium mole fractions and the partial
also given which calculates the gas constant,
for the products of combustion of any hydrocar
pressure and equivalence ratio. Some examples of the uses of the programs are also given.
derivatives of the mole fractions with respect to temperature, pressure and equivalence ratio
enthalpy, internal energy and the partial deriv atives of these with respect to temperature,
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
2 used to obtain the internal energy of the cell.
As an example, some models used for the spark ignition engine products divide the products in
to subsystems (1) but assume the same pressure
for each subsystem while others (11) assign one temperature for the entire product mass. The first law for each system of products is
necessarily be true for the reactants (12).
The specific gas constant, where M is the mean molecular weight, is a function of p, T,
F. The derivatives dp/dt in Equation 2 may be obtained by differentiation of the ideal gas law and substitution of:
for the dR/dt term.
This gives:
The ideal gas equation cannot be explicitly solved for p in terms of density, temperature and composition. p, T, F and ? are known
If
at state "n" then the pressure at state n+1 may
be obtained by a first order approximation to
In order to solve Equation 1 it is necessary R. This gives:
to supply an equation of state and the internal
energy as a function of the composition and two state variables such as pressure and temperature.
The volume must be determined from the fact that
the sum of all system volumes is equal to the known total volume. A mass balance plus appro
priate flow equations determine the mass in the
system. The heat transfer may be modeled by the use of empirical convection and radiation equa
tions which generally will be functions of the temperature, density and composition as well as some characteristic velocity.
For a numerical solution of Equation 1 it
is convenient to introduce the internal energy as the sum of the chemical energy and sensible energy. For equilibrium the energy then can be expressed as a function of the variables p, T and F. The variable, F, is the fuel/air equiva lence ratio. For products formed directly from the reactants (fuel plus air) the value of F is the equivalence ratio of the reactants. For products of such a reaction mixed with air or additional products of the same fuel the value
Equation 4 may be solved numerically using
Equation 5 to evaluate p and appropriate differ
ential conservation equations for m and F. An example of this technique is given later, but first we shall consider the program to evaluate u and R and their derivatives.
EQUILIBRIUM THERMODYNAMICS
Let the fuel CnHmOlNk and air (13) at equiv
alence ratio F react and the products subject to
temperature T and pressure p attain equilibrium.
The numbers n and m should be nonzero; l and k
may or may not be zero.
of F may be computed for a hypothetical reactant
F which would give the same atomic properties as the resulting mixture. The term du/dt may be written The partial derivative terms can be calculated
directly from the equilibrium thermodynamics as functions of T, p and F.
Little error is made for normal engine con ditions if the products are assumed to be an ideal gas mixture although the same may not
where x1 through x12 are mole fractions of the product species. The number, x13, represents
the moles of fuel that will give one mole of
products.
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
3
The left side of the equation may be written as Atom balances for the various elements give
REDUCTION IN NUMBER OF EQUATIONS
Equation (8-1) gives
Equation (8-5) gives
These are now used to eliminate X12 and
x13 in Equations 8 and 9.
The constraint that the mole fraction of all the
products add up to unity requires that
To solve for the 13 unknowns, we need 7 more
Denote the ratios
equations which are provided by the criteria of equilibrium among the products, expressed by the following 7 non-redundant hypothetical reactions. Now Equations 11 can be used in Equations 12 to
eliminate all variables except x4, x6, x8 and
X11. The resulting four nonlinear equations with four unknowns may be expressed symbolically by Assume we know a vector
Where p is the pressure in atmospheres.
The equilibrium constants are curve fitted from
data in JANAF Thermochemical Tables, Second Edi
tion (1971). Details of fitting are given in
Appendix A.
The expressions for equilibrium
which is reasonably close to the true solution vector constants can be rearranged to express mole
fractions of all the products of combustion in terms of x4, x6, x8 and X11 the mole-fractions
of H2, CO, O2 and N2 respectively.
The functions on the left side of the equations
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
4 can then be expanded around the known vector as a Taylor 's series.
Substituting for x6, x9 and X10, rearranging terms and dividing through by x13,
Neglecting the partial derivatives of second or
der and higher we get a set of linear equations
where ?xi are the approximate corrections.
Thus
The quantity X13 can be estimated with good ac
curacy from the condition that the mole fractions should add up to unity.
where functions fj and their partial derivatives are evaluated at the known vector.
This set of linear equations can be solved for ?x4, ?x6, ?x8 and ?x11 using a Gaussian
For F = 1 a good estimate of X13 obtained be can from
elimination method. The improved values are then
For F > 1 an estimate of x13 can be obtained
The improved vector can now be used to evaluate the partial derivatives and functions.
from
Further corrections can then be calculated and
applied. The iterative procedure is to be con tinued until the relative changes in all the com ponents are less than a specified value.
Details
of the computation of the elements of Eq. 14 are given in Appendix B.
method, the other unknowns may be obtained di
rectly by substitution into Equations 16-1, 16-2,
INITIAL ESTIMATION OF MOLE FRACTIONS
The Newton Raphson iteration is not self
starting. Unless an approximate solution is known from a previous execution of the subroutine at conditions close to those presently encoun
tered, it is necessary to make an initial esti mation to get the iteration started. The follow
ing method of estimation developed through judi
cious simplifying assumptions is found to be re markably good.
Assume that the products are only H2, CO,
O2, H2O, CO2, N2 and Ar.
Substitution of the estimated value of x13
into Eq. 16-6 gives an equation in the single un known x8. After solving Eq. 16-6 by Newton 's
Equation 6 becomes
and 16-3.
These initial estimates can then be
used to start the solution of Equations 13 by
use of Equations 14 and 15.
PARTIAL DERIVATIVES OF THE MOLE FRACTIONS
Starting with Equations 13, the functions
and hence the mole fractions are dependent upon
temperature, pressure and equivalence ratio.
For example, taking the total derivatives of each equation with respect to temperature we get four linear simultaneous equations in four un knowns. From the C balance and Eq. 11-7
It is seen that the coefficient matrix is
the Jacobian that appeared in Newton Raphson
From the H balance and Eq. 11-6
iteration. Similarly, it is seen that the matrix equations for solving for the partial derivatives with respect to pressure and equivalence ratio would also be identical except that T is replaced by p or F. Equations for the computations of the elements of the coefficient matrix and the con
and from the N, Ar, and O balances
stant vector are given in Appendix C. Solution of the equations sets (Equations 18 and similar ones for p and F) can be accomplished in a straight forward manner using the Gaussian elim ination method. As shown in Appendix C the re maining partial derivatives of x1,x2, etc. are easily obtained by use of Equations 11 and 8-5 and direct substitution.
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
5
The computation of the equilibrium mole fractions and partial derivatives ofthe mole
fractions with respect to T, p, F iscarried
out by a FORTRAN coded program EQMD.A guide for users is given in Appendix D anda program listing is given in Appendix E.
ENTHALPY INTERNAL
AND
ENERGY
The absolute enthalpy used here is defined as The values of the specific heats and heats of
formation were taken from Ref. 14 and were used
with a table look up and linear interpolation program. Using the results of Program EQMD the en
thalpy of the mixture, h, is easily calculated as is the average molecular weight, M, and gas
constant, R.
we consider the expansion of a closed system of products. For the adiabatic constant pressure reaction we first calculate the reactant mixture enthalpy.
If the fuel heat of formation is not known it can
be estimated using the heating value and H/C ra
tio (15). For higher pressures the non-ideal ef fects (12) should be included. Having computed
the reactant enthalpy, hr, we hold p and F con
stant and solve for T
The partial derivatives are also easily computed. Using Newton 's method
For example,
The other expressions for the partiais of h, u
and R with respect to T, p and F are given in
Appendix C.
The computation of u, h, R and their par
where Tn is the first estimate and Tn+1 is the improved estimate. The first estimate can be any reasonable number but preferably larger than the expected solution. For our example, we take methane at 536.4°R, one atmosphere and F = 0.8,
0.9, 1.0, 1.1, 1.2, Figure 1 shows a plot of the tials is carried out by FORTRAN coded program PER. result and Table 1 shows the progressive esti
A users guide is given in Appendix D and a mates of T for the F = 1.1 case. A plot of program listing is given in Appendix E.
?(h-hr)/?F at the calculated F and T values of Figure 1 gives the dashed curve of Figure 1. The zero point determines the F at the maximum flame
EXAMPLE CALCULATIONS temperature. The following example of expansion of a closed system of products of combustion is given to show the use of the programs. First we con sider the classical problem of computing the con only to illustrate the application of the pro stant pressure adiabatic flame temperature. Second, grams to such calculations. The kinematic and
Two simple cases are given here as examples
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
6
Total mass of gas inside the system is
0.001244 lbm and there is no flow in or out of
the system. The fuel is C8H16 and equivalence
ratio F is 1.2.
Initial conditions are T =
5600°R and p = 533 psi at ? = 20°. The expan sion is assumed to follow shifting equilibrium.
Improved Euler method is used for solving the differential equations. The steps taken at each interval are listed below.
Let the temperature Tn and pressure Pn be
known at time tn where the derivative, dT/dt, is calculated using Eq. 4.
heat transfer models used below are extremely simplistic. are used as starting condi
tions for the subsequent time interval.
The results of the program for expansion from ? = 20° to 180° is plotted in Figure 2.
In addition, the integrated values of pdV and Q were also accumulated during the calculation.
The change in internal energy was calculated from initial and final conditions.
For a step
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
7 size of 1° crank angle:
"Engine Simulation Studies Conducted at the Uni versity of Wisconsin." Tech. Report No. 11398,
U.S. Army Tank Automotive Command, March 1972.
3.L.W. Evers, P.S. Myers, and O.A. Uyehara,
"A Search For a Low Nitric Oxide Engine." Cen tral States Combustion Institute Meeting,
Madison, Wisconsin, Spring 1974.
4.H.K. Newhall, "Kinetics of Engine-Generated
ACKNOWLEDGMENT
Nitric Oxides and Carbon Monoxide."
Twelfth
Symposium (International) on Combustion, 1969.
5.Tizard and Pye, "Empire Motor Fuels Commit tee Report." Proc. Inst. of Automotive Engrs.,
This report was submitted in partial ful fillment of Grant Number R-802589-01-0 by the
University of Wisconsin-Madison under the partial Vol. 18, No. 1, 1923. sponsorship of the Environmental Protection
6.H.N. Powell, A. Schaffer and S.N. Suciu,
Agency. Work was completed as of July 1, 1974.
"Thermodynamic Properties--Properties of Combus tion Gases, System: CnH2n-Air." General Elec
NOMENCLATURE
Aij afi/axj dQ/dt - rate of heat
-
fi
-
transfer nonlinear function of mole fractions
F - fuel/air equivalence ratio h - enthalpy per unit mass
tric Co., Cincinnati, Ohio, 1955.
7.H.K. Newhall and E.S. Starkman, "Thermody
namic Properties of Octane and Air for Engine
Performance Calculations." SAE Paper 633G, 1963.
8.R.J. Steffensen, "A FORTRAN IV Program of
Thermochemical Calculations Involving the Ele
hi hr -
enthalpy per mole of species i except
ments Al, B, Be, C, F, H, Li, Mg, N and 0 and
Their Compounds," PhD Thesis, Purdue University,
-
reactant enthalpy per unit mass
9.S. Gordon, "Complex Chemical Equilibrium
Calculations." NASA Sp-239, Kinetics and Thermo
m
-
Ki
where otherwise defined
-
partial pressure equilibrium constant
mass or number of H atoms in fuel molecule M
-
n
-
? -
molecular weight of mixture
number of carbon atoms in fuel mole cule pressure
r,r ',r" - proportions of N2, O2, Ar in air
R - specific gas constant of mixture
Ro - universal gas constant t - time
T - temperature
?xi/?xj
Tij - gas side wall surface temperature
Tw
-
u
-
internal energy per unit mass
V - volume
xi
? overdot -
mole fraction
-
gas constant correction factor indicates d/dt
-
REFERENCES
1966.
dynamics in High Temperature Gases, March 1970.
10.G.L. Borman, "Mathematical Simulation of
Internal Combustion Engine Processes and Perform
ance Including Comparisons with Experiments."
PhD Thesis, University of Wisconsin, Madison,
1964.
11.R.B. Krieger, and G.L. Borman, "The Compu
tation of Apparent Heat Release for Internal Com
bustion Engines." ASME Proc. Diesel Gas Power,
ASME Paper 66-WA/DGP-4, 1966.
12.J. Maniotes, "Ideal and Non-ideal Theoreti
cal Equilibrium Calculations for Mixtures of
Methane and Air." PhD Thesis, Purdue University
1962.
13."CRC Handbook of Chemistry and Physics."
71st Edition, Sea level atmospheric composition
of dry atmosphere, pp. F-147, 1970-71.
14."JANAF Thermochemical Tables." Second Edi
tion, The Dow Chemical Co., Midland, Michigan,
1971.
15.H.N. Powell, "Applications of an Enthalpy
1.P. Blumberg and J.T. Kummer, "Prediction of Fuel/Air Ratio Diagram to 'First Law ' Combustion
NO Formation in Spark-Ignited Engines--An Analy Problems." Trans. ASME, pp. 1129-1138, July
1957.
sis of Methods of Control." C.S.T., Vol. 4, pp. 73-95, 1971.
16.I. Klotz, "Introduction to Chemical Ther
2.G.L. Borman, P.S. Myers, and O.A. Uyehara, modynamics." W.A. Benjamin Co., New York, 1964.
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
8
APPENDIX A
CURVE FITTING OF EQUILIBRIUM CONSTANTS
The equilibrium constants data was taken
from JANAF Thermochemical Tables (14), where lNew Text og10Kp (formation) for all species are tabulated as functions of the absolute temperature (°K).
Equilibrium constants for the reactions consid ered in the equilibrium thermodynamics were cal culated using the relation
where T is the absolute temperature and A, B, C,
D, E are constants.
This model was used to fit the tabulated
data by means of a least squares fitting program.
A trade-off between the conflicting interests of
obtaining either high accuracy in a narrow range of temperature or lower accuracy in a wider range
led to choosing 600 to 4000°K (1080 to 7200°R) as
the range most relevant to the study of combus tion phenomena in engines.
The log Kp predicted by the equations were
compared with the original data and the devia tions were less than 0.0009. (The original data is tabulated only to the third decimal place and so there is an inherent uncertainty of 0.0005.
Hence the deviations were not considered signi
Theoretical studies (16) suggest a function relationship of the type
ficant.)
A transformed temperature TA defined as
0.005 T/9 where T is in °R was used for fitting.
The constants A, B, C, D and E are listed in
Table A-1.
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
9
APPENDIX B
MOLE FRACTION COMPUTATION EQUATIONS
In the specific problem under consideration, we define for convenience, the following partial derivatives From Eqs. 11,
we get
With the above notations, the elements of
Eq. 14 written as a matrix equation [A][?x] = [B] can be expressed as
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
10
APPENDIX C
PARTIAL DERIVATIVE COMPUTATIONEQUATIONS
Define
From Equations 11, Yi are seen to be functions of x4, X6, x8 and x11 only.
Substituting for xi from Equation Cl in
Equation 12 and differentiating with respect to
T, we derive
Define
The model used for fitting equilibrium constants
Partial derivatives with respect to p are similar
KP is
to above except that ?/?T is replaced by ?/?p.
Note that C5 and C7 are not functions of p and hence terms involving ?C5/?p and ?C7/?p drop out.
The parameters d2, d3 and d4 are functions of F, but none of the Ci are. Hence
From which we derive
CALCULATION OF PARTIAL DERIVATIVES OF REMAINING
MOLE FRACTIONS
From Equations 11
From the expression for
Ci and di
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
11
CALCULATION OF PARTIAL DERIVATIVES of R, h and u
From Equations 19 to 22
For partial derivatives with respect to p, substitute ?/?p where ?/?T appears above.
Example:
For partial derivatives with respect to F, substitute zero where ?Ci/?T appears and ?/?F elsewhere. Example:
Note that
the specific heats of the
species at constant pressure.
For partial derivatives with respect to p and F, replace ?/T with ?/?p and ?/?F respec tively. The exceptions are
Exception:
APPENDIX D
GUIPE TO USING SUBROUTINES EQMD AND PER
(where AN, AM, AL are the number of C, H and O atoms in fuel molecules) free carbon will defi
PURPOSE OF THE SUBROUTINES - Subroutine EQMD nitely be formed. For C8H16, this maximum F is can calculate the mole fractions xi of the prod 3. ucts of combustion of any hydrocarbon fuel and
The products of combustion are assumed to air, at equilibrium under specified conditions of be ideal gases. This assumption is not valid at extremely high pressures (12). temperature T, pressure p and equivalence ratio
HOW TO ACCESS THE SUBROUTINES - Transfer of
F. It can also optionally calculate the partial
derivatives of the xi with respect to T, p and F. information between the calling program and the
subroutines are through labelled COMMON areas
Subroutine PER can calculate the average and hence every program that calls the subrou molecular weight M, gas constant R, enthalphy h and internal energy u of the equilibrium products tines must contain of combustion. It can also optionally calculate the partial derivatives of R, h and u with re spect to t, p and F.
In addition to carbon and hydrogen, the fuel may or may not contain oxygen and nitrogen atoms.
The product species considered are H, O, N, H2,
OH, CO, NO, O2, H2O, CO2, N2 and Ar in gas phase.
LIMITATIONS OF THE SUBROUTINE - The equilib
rium constants used in the subroutines were
fitted as a function of temperature in the range
1080°R to 7200°R. The subroutines cannot be used outside this range.
The subroutine cannot handle the formation
of free carbon. lence ratio,
It can be shown that for equiva
The above variables are input to the subroutine
and should be defined at the time of calling.
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
12
The following are output variables.
IERR - Error code with meaning as specified
tered during calculation of
partial derivatives of mole
later.
XEQ(12) - Mole fractions of the product species.
XPA - The numbering code is given later.
A parameter used in the subroutine.
It represents the modes of fuel that
will give one mole of products.
fractions.
Immediately after every call to the subrou tine, the value of TERR must be checked. It should be zero. If it is not, the numbers
stored in the other output variable locations are meaningless and should not be used. Prefer rably print out the value of IERR and also the
condition (AN,AM,AL,AK,F,T,P,KLO)at which the failure occurred.
AVM - Average molecular weight of products of combustion
R - Gas constant BTU
H - Enthalpy BTU
°R-1
U - Internal energy BTU
The partial derivatives defined earlier
have units of the variable per °R, per psi , per unit fuel air equivalence ratio (dimensionless) respectively. For example,
DXT(1) has units of °R-1
DRP has units of BTU
DHF has units of BTU
psi-1
°R-1
REFERENCE STATES - For computation of en thalpies and internal energies, the following states were assigned zero value of enthalpies.
H2 (ideal gas) at 0°R
O2 (ideal gas) at 0°R
N2 (ideal gas) at 0°R
C (Graphite) at 0°R
Ar (ideal gas) at 0°R
HOW TO CALL SUBROUTINE PER - Both the
$INCLUDE cards and all three of the COMMON state
ments are required. The subroutine PER will be executed every time the following statement is encountered. CONTROL VARIABLE K-]OL KLO = 0
The subroutine will make a fresh
estimate before iterating
KLO = 1
Closeness criteria built into the subroutine will decide whether or not to make a fresh estimate. THIS
IS RECOMMENDED.
KLO = 2
The subroutine will use input values of H2, CO, O2 and N2 mole fractions as initial estimates
ERROR CODE I-]ERR
IERR 0
=
= 1
= 2
No error.
Specified temperature out side 1080 to 7200°R range.
Specified equivalence ratio too high as to cause forma tion of free carbon.
= 3
Attempted to calculate an estimate; but none found within
= 4
outputs of EQMD are also available after a call
to PER.
It is not necessary to call EQMD again.
However if only the mole fractions and/or their
partial derivatives are needed and not the prop
erties, some computation time could be saved by directly calling EQMD as explained below.
HOW TO CALL SUBROUTINE EQMD ALONE - In this case # INCLUDE CO*PER and COMMON/PROP/ cards are not needed. The subroutine EQMD will be executed each time the following statement is encountered.
A singular matrix was encoun
= 5
The subroutine PER internally calls EQMD
with the same value of JDR and all the relevant
mole fractions
Few mole fractions were cal
reasonable limits.
tered during calculation of
culated as negative during iteration. = 6
=7
Newton Raphson iteration did not converge in 25 attempts.
A singular matrix was encoun
SAMPLE PROGRAM - Listing of a sample program which calls subroutine PER is given on the fol lowing page.
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
13
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
14
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
15
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
16
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
17
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
18
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
19
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
20
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
21
This paper is subject to revision. Statements and opinions advanced in papers or discussion are the author 's and are
in Transactions. For permission to publish this paper in full or in part,
SAE
been edited by SAE for uniform styling and format. Dis cussion will be printed with the paper if it is published
tion through SAE should send the manuscript or a 300 word abstract of a pro posed manuscript to: Secretary, Engineering Activities Board, SAE.
his responsibility, not the Society 's; however, the paper has
contact the SAE Publications Division.
Persons wishing to submit papers to be considered for presentation or publica
page booklet.
Printed in U.S.A.
References: Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM 2
SOCIETY OF AUTOMOTIVE ENGINEERS, INC.
400 Commonwealth Drive, Warrendale, Pa. 15096
A Computer Program for
Calculating Properties of
Equilibrium Combustion
Products with Some
Applications to I.C. Engines
Cherian Olikara and Gary L. Borman
University of Wisconsin
Automotive Engineering Congress and Exposition
Detroit, Michigan
February 24-28, 1975
750468
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
Copyright © Society of Automotive Engineers, Inc. 1975
All rights reserved.
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
750468
A Computer Program for
Calculating Properties of
Equilibrium Combustion
Products with Some
Applications to I.C. Engines
Cherian Olikara and Gary L. Borman
University of Wisconsin
MANY MODELS FOR ENGINE COMBUSTION PROBLEMS use
the First Law of Thermodynamics as applied to
either the entire cylinder contents or to sub
systems. Examples include models for spark ig nition engines (1)*, diesel engines (2) and stratified charge engines (3). Typically the major species of the products of combustion may be assumed to follow a shifting equilibrium process for thermodynamic purposes. For pur
are given by References 6-9. The NASA-Lewis program (9) is very extensive and includes ther modynamic data for hundreds of species.
The purpose of the present program is very narrow by comparison to the NASA program. Its purpose is to specifically deal only with the gas phase products of combustion of hydrocarbon
fuels (containing, C, H, O, N atoms) and air.
The program is however extended to calculate the
partial derivatives of internal energy and molec ular weight which are helpful when numerically at the equilibrium concentration (4). Because solving the first law as a differential equation of this wide spread use it is important to have in time using the method of Reference 10. Most a rapid means of calculating the equilibrium in importantly, the program is to provide a rapid ternal energy, molecular weight and species mole means of calculation comparable to the use of fractions of the products of hydrocarbon-air com the regression analysis equations given in Ref poses of chemical kinetics calculations, many of the major species may also be assumed to be
bustion.
erence 11.
Calculation of the equilibrium composition and internal energy of combustion products for
ENERGY EQUATION FORMULATION
engine cycles goes back to the early work of
Tizard and Pye (5) and has been the subject of numerous computer studies, a sampling of which
*Numbers in parentheses designate References at end of paper
In order to treat the cylinder gas system
thermodynamically the system may be divided into
cells such that the composition and temperature
of each cell is uniform. This implies that the molar average temperature of the cell may be
ABSTRACT
bon fuel and air is described.
A subroutine is
A computer program which rapidly calculates the equilibrium mole fractions and the partial
also given which calculates the gas constant,
for the products of combustion of any hydrocar
pressure and equivalence ratio. Some examples of the uses of the programs are also given.
derivatives of the mole fractions with respect to temperature, pressure and equivalence ratio
enthalpy, internal energy and the partial deriv atives of these with respect to temperature,
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
2 used to obtain the internal energy of the cell.
As an example, some models used for the spark ignition engine products divide the products in
to subsystems (1) but assume the same pressure
for each subsystem while others (11) assign one temperature for the entire product mass. The first law for each system of products is
necessarily be true for the reactants (12).
The specific gas constant, where M is the mean molecular weight, is a function of p, T,
F. The derivatives dp/dt in Equation 2 may be obtained by differentiation of the ideal gas law and substitution of:
for the dR/dt term.
This gives:
The ideal gas equation cannot be explicitly solved for p in terms of density, temperature and composition. p, T, F and ? are known
If
at state "n" then the pressure at state n+1 may
be obtained by a first order approximation to
In order to solve Equation 1 it is necessary R. This gives:
to supply an equation of state and the internal
energy as a function of the composition and two state variables such as pressure and temperature.
The volume must be determined from the fact that
the sum of all system volumes is equal to the known total volume. A mass balance plus appro
priate flow equations determine the mass in the
system. The heat transfer may be modeled by the use of empirical convection and radiation equa
tions which generally will be functions of the temperature, density and composition as well as some characteristic velocity.
For a numerical solution of Equation 1 it
is convenient to introduce the internal energy as the sum of the chemical energy and sensible energy. For equilibrium the energy then can be expressed as a function of the variables p, T and F. The variable, F, is the fuel/air equiva lence ratio. For products formed directly from the reactants (fuel plus air) the value of F is the equivalence ratio of the reactants. For products of such a reaction mixed with air or additional products of the same fuel the value
Equation 4 may be solved numerically using
Equation 5 to evaluate p and appropriate differ
ential conservation equations for m and F. An example of this technique is given later, but first we shall consider the program to evaluate u and R and their derivatives.
EQUILIBRIUM THERMODYNAMICS
Let the fuel CnHmOlNk and air (13) at equiv
alence ratio F react and the products subject to
temperature T and pressure p attain equilibrium.
The numbers n and m should be nonzero; l and k
may or may not be zero.
of F may be computed for a hypothetical reactant
F which would give the same atomic properties as the resulting mixture. The term du/dt may be written The partial derivative terms can be calculated
directly from the equilibrium thermodynamics as functions of T, p and F.
Little error is made for normal engine con ditions if the products are assumed to be an ideal gas mixture although the same may not
where x1 through x12 are mole fractions of the product species. The number, x13, represents
the moles of fuel that will give one mole of
products.
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
3
The left side of the equation may be written as Atom balances for the various elements give
REDUCTION IN NUMBER OF EQUATIONS
Equation (8-1) gives
Equation (8-5) gives
These are now used to eliminate X12 and
x13 in Equations 8 and 9.
The constraint that the mole fraction of all the
products add up to unity requires that
To solve for the 13 unknowns, we need 7 more
Denote the ratios
equations which are provided by the criteria of equilibrium among the products, expressed by the following 7 non-redundant hypothetical reactions. Now Equations 11 can be used in Equations 12 to
eliminate all variables except x4, x6, x8 and
X11. The resulting four nonlinear equations with four unknowns may be expressed symbolically by Assume we know a vector
Where p is the pressure in atmospheres.
The equilibrium constants are curve fitted from
data in JANAF Thermochemical Tables, Second Edi
tion (1971). Details of fitting are given in
Appendix A.
The expressions for equilibrium
which is reasonably close to the true solution vector constants can be rearranged to express mole
fractions of all the products of combustion in terms of x4, x6, x8 and X11 the mole-fractions
of H2, CO, O2 and N2 respectively.
The functions on the left side of the equations
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
4 can then be expanded around the known vector as a Taylor 's series.
Substituting for x6, x9 and X10, rearranging terms and dividing through by x13,
Neglecting the partial derivatives of second or
der and higher we get a set of linear equations
where ?xi are the approximate corrections.
Thus
The quantity X13 can be estimated with good ac
curacy from the condition that the mole fractions should add up to unity.
where functions fj and their partial derivatives are evaluated at the known vector.
This set of linear equations can be solved for ?x4, ?x6, ?x8 and ?x11 using a Gaussian
For F = 1 a good estimate of X13 obtained be can from
elimination method. The improved values are then
For F > 1 an estimate of x13 can be obtained
The improved vector can now be used to evaluate the partial derivatives and functions.
from
Further corrections can then be calculated and
applied. The iterative procedure is to be con tinued until the relative changes in all the com ponents are less than a specified value.
Details
of the computation of the elements of Eq. 14 are given in Appendix B.
method, the other unknowns may be obtained di
rectly by substitution into Equations 16-1, 16-2,
INITIAL ESTIMATION OF MOLE FRACTIONS
The Newton Raphson iteration is not self
starting. Unless an approximate solution is known from a previous execution of the subroutine at conditions close to those presently encoun
tered, it is necessary to make an initial esti mation to get the iteration started. The follow
ing method of estimation developed through judi
cious simplifying assumptions is found to be re markably good.
Assume that the products are only H2, CO,
O2, H2O, CO2, N2 and Ar.
Substitution of the estimated value of x13
into Eq. 16-6 gives an equation in the single un known x8. After solving Eq. 16-6 by Newton 's
Equation 6 becomes
and 16-3.
These initial estimates can then be
used to start the solution of Equations 13 by
use of Equations 14 and 15.
PARTIAL DERIVATIVES OF THE MOLE FRACTIONS
Starting with Equations 13, the functions
and hence the mole fractions are dependent upon
temperature, pressure and equivalence ratio.
For example, taking the total derivatives of each equation with respect to temperature we get four linear simultaneous equations in four un knowns. From the C balance and Eq. 11-7
It is seen that the coefficient matrix is
the Jacobian that appeared in Newton Raphson
From the H balance and Eq. 11-6
iteration. Similarly, it is seen that the matrix equations for solving for the partial derivatives with respect to pressure and equivalence ratio would also be identical except that T is replaced by p or F. Equations for the computations of the elements of the coefficient matrix and the con
and from the N, Ar, and O balances
stant vector are given in Appendix C. Solution of the equations sets (Equations 18 and similar ones for p and F) can be accomplished in a straight forward manner using the Gaussian elim ination method. As shown in Appendix C the re maining partial derivatives of x1,x2, etc. are easily obtained by use of Equations 11 and 8-5 and direct substitution.
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
5
The computation of the equilibrium mole fractions and partial derivatives ofthe mole
fractions with respect to T, p, F iscarried
out by a FORTRAN coded program EQMD.A guide for users is given in Appendix D anda program listing is given in Appendix E.
ENTHALPY INTERNAL
AND
ENERGY
The absolute enthalpy used here is defined as The values of the specific heats and heats of
formation were taken from Ref. 14 and were used
with a table look up and linear interpolation program. Using the results of Program EQMD the en
thalpy of the mixture, h, is easily calculated as is the average molecular weight, M, and gas
constant, R.
we consider the expansion of a closed system of products. For the adiabatic constant pressure reaction we first calculate the reactant mixture enthalpy.
If the fuel heat of formation is not known it can
be estimated using the heating value and H/C ra
tio (15). For higher pressures the non-ideal ef fects (12) should be included. Having computed
the reactant enthalpy, hr, we hold p and F con
stant and solve for T
The partial derivatives are also easily computed. Using Newton 's method
For example,
The other expressions for the partiais of h, u
and R with respect to T, p and F are given in
Appendix C.
The computation of u, h, R and their par
where Tn is the first estimate and Tn+1 is the improved estimate. The first estimate can be any reasonable number but preferably larger than the expected solution. For our example, we take methane at 536.4°R, one atmosphere and F = 0.8,
0.9, 1.0, 1.1, 1.2, Figure 1 shows a plot of the tials is carried out by FORTRAN coded program PER. result and Table 1 shows the progressive esti
A users guide is given in Appendix D and a mates of T for the F = 1.1 case. A plot of program listing is given in Appendix E.
?(h-hr)/?F at the calculated F and T values of Figure 1 gives the dashed curve of Figure 1. The zero point determines the F at the maximum flame
EXAMPLE CALCULATIONS temperature. The following example of expansion of a closed system of products of combustion is given to show the use of the programs. First we con sider the classical problem of computing the con only to illustrate the application of the pro stant pressure adiabatic flame temperature. Second, grams to such calculations. The kinematic and
Two simple cases are given here as examples
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
6
Total mass of gas inside the system is
0.001244 lbm and there is no flow in or out of
the system. The fuel is C8H16 and equivalence
ratio F is 1.2.
Initial conditions are T =
5600°R and p = 533 psi at ? = 20°. The expan sion is assumed to follow shifting equilibrium.
Improved Euler method is used for solving the differential equations. The steps taken at each interval are listed below.
Let the temperature Tn and pressure Pn be
known at time tn where the derivative, dT/dt, is calculated using Eq. 4.
heat transfer models used below are extremely simplistic. are used as starting condi
tions for the subsequent time interval.
The results of the program for expansion from ? = 20° to 180° is plotted in Figure 2.
In addition, the integrated values of pdV and Q were also accumulated during the calculation.
The change in internal energy was calculated from initial and final conditions.
For a step
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
7 size of 1° crank angle:
"Engine Simulation Studies Conducted at the Uni versity of Wisconsin." Tech. Report No. 11398,
U.S. Army Tank Automotive Command, March 1972.
3.L.W. Evers, P.S. Myers, and O.A. Uyehara,
"A Search For a Low Nitric Oxide Engine." Cen tral States Combustion Institute Meeting,
Madison, Wisconsin, Spring 1974.
4.H.K. Newhall, "Kinetics of Engine-Generated
ACKNOWLEDGMENT
Nitric Oxides and Carbon Monoxide."
Twelfth
Symposium (International) on Combustion, 1969.
5.Tizard and Pye, "Empire Motor Fuels Commit tee Report." Proc. Inst. of Automotive Engrs.,
This report was submitted in partial ful fillment of Grant Number R-802589-01-0 by the
University of Wisconsin-Madison under the partial Vol. 18, No. 1, 1923. sponsorship of the Environmental Protection
6.H.N. Powell, A. Schaffer and S.N. Suciu,
Agency. Work was completed as of July 1, 1974.
"Thermodynamic Properties--Properties of Combus tion Gases, System: CnH2n-Air." General Elec
NOMENCLATURE
Aij afi/axj dQ/dt - rate of heat
-
fi
-
transfer nonlinear function of mole fractions
F - fuel/air equivalence ratio h - enthalpy per unit mass
tric Co., Cincinnati, Ohio, 1955.
7.H.K. Newhall and E.S. Starkman, "Thermody
namic Properties of Octane and Air for Engine
Performance Calculations." SAE Paper 633G, 1963.
8.R.J. Steffensen, "A FORTRAN IV Program of
Thermochemical Calculations Involving the Ele
hi hr -
enthalpy per mole of species i except
ments Al, B, Be, C, F, H, Li, Mg, N and 0 and
Their Compounds," PhD Thesis, Purdue University,
-
reactant enthalpy per unit mass
9.S. Gordon, "Complex Chemical Equilibrium
Calculations." NASA Sp-239, Kinetics and Thermo
m
-
Ki
where otherwise defined
-
partial pressure equilibrium constant
mass or number of H atoms in fuel molecule M
-
n
-
? -
molecular weight of mixture
number of carbon atoms in fuel mole cule pressure
r,r ',r" - proportions of N2, O2, Ar in air
R - specific gas constant of mixture
Ro - universal gas constant t - time
T - temperature
?xi/?xj
Tij - gas side wall surface temperature
Tw
-
u
-
internal energy per unit mass
V - volume
xi
? overdot -
mole fraction
-
gas constant correction factor indicates d/dt
-
REFERENCES
1966.
dynamics in High Temperature Gases, March 1970.
10.G.L. Borman, "Mathematical Simulation of
Internal Combustion Engine Processes and Perform
ance Including Comparisons with Experiments."
PhD Thesis, University of Wisconsin, Madison,
1964.
11.R.B. Krieger, and G.L. Borman, "The Compu
tation of Apparent Heat Release for Internal Com
bustion Engines." ASME Proc. Diesel Gas Power,
ASME Paper 66-WA/DGP-4, 1966.
12.J. Maniotes, "Ideal and Non-ideal Theoreti
cal Equilibrium Calculations for Mixtures of
Methane and Air." PhD Thesis, Purdue University
1962.
13."CRC Handbook of Chemistry and Physics."
71st Edition, Sea level atmospheric composition
of dry atmosphere, pp. F-147, 1970-71.
14."JANAF Thermochemical Tables." Second Edi
tion, The Dow Chemical Co., Midland, Michigan,
1971.
15.H.N. Powell, "Applications of an Enthalpy
1.P. Blumberg and J.T. Kummer, "Prediction of Fuel/Air Ratio Diagram to 'First Law ' Combustion
NO Formation in Spark-Ignited Engines--An Analy Problems." Trans. ASME, pp. 1129-1138, July
1957.
sis of Methods of Control." C.S.T., Vol. 4, pp. 73-95, 1971.
16.I. Klotz, "Introduction to Chemical Ther
2.G.L. Borman, P.S. Myers, and O.A. Uyehara, modynamics." W.A. Benjamin Co., New York, 1964.
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
8
APPENDIX A
CURVE FITTING OF EQUILIBRIUM CONSTANTS
The equilibrium constants data was taken
from JANAF Thermochemical Tables (14), where lNew Text og10Kp (formation) for all species are tabulated as functions of the absolute temperature (°K).
Equilibrium constants for the reactions consid ered in the equilibrium thermodynamics were cal culated using the relation
where T is the absolute temperature and A, B, C,
D, E are constants.
This model was used to fit the tabulated
data by means of a least squares fitting program.
A trade-off between the conflicting interests of
obtaining either high accuracy in a narrow range of temperature or lower accuracy in a wider range
led to choosing 600 to 4000°K (1080 to 7200°R) as
the range most relevant to the study of combus tion phenomena in engines.
The log Kp predicted by the equations were
compared with the original data and the devia tions were less than 0.0009. (The original data is tabulated only to the third decimal place and so there is an inherent uncertainty of 0.0005.
Hence the deviations were not considered signi
Theoretical studies (16) suggest a function relationship of the type
ficant.)
A transformed temperature TA defined as
0.005 T/9 where T is in °R was used for fitting.
The constants A, B, C, D and E are listed in
Table A-1.
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
9
APPENDIX B
MOLE FRACTION COMPUTATION EQUATIONS
In the specific problem under consideration, we define for convenience, the following partial derivatives From Eqs. 11,
we get
With the above notations, the elements of
Eq. 14 written as a matrix equation [A][?x] = [B] can be expressed as
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
10
APPENDIX C
PARTIAL DERIVATIVE COMPUTATIONEQUATIONS
Define
From Equations 11, Yi are seen to be functions of x4, X6, x8 and x11 only.
Substituting for xi from Equation Cl in
Equation 12 and differentiating with respect to
T, we derive
Define
The model used for fitting equilibrium constants
Partial derivatives with respect to p are similar
KP is
to above except that ?/?T is replaced by ?/?p.
Note that C5 and C7 are not functions of p and hence terms involving ?C5/?p and ?C7/?p drop out.
The parameters d2, d3 and d4 are functions of F, but none of the Ci are. Hence
From which we derive
CALCULATION OF PARTIAL DERIVATIVES OF REMAINING
MOLE FRACTIONS
From Equations 11
From the expression for
Ci and di
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
11
CALCULATION OF PARTIAL DERIVATIVES of R, h and u
From Equations 19 to 22
For partial derivatives with respect to p, substitute ?/?p where ?/?T appears above.
Example:
For partial derivatives with respect to F, substitute zero where ?Ci/?T appears and ?/?F elsewhere. Example:
Note that
the specific heats of the
species at constant pressure.
For partial derivatives with respect to p and F, replace ?/T with ?/?p and ?/?F respec tively. The exceptions are
Exception:
APPENDIX D
GUIPE TO USING SUBROUTINES EQMD AND PER
(where AN, AM, AL are the number of C, H and O atoms in fuel molecules) free carbon will defi
PURPOSE OF THE SUBROUTINES - Subroutine EQMD nitely be formed. For C8H16, this maximum F is can calculate the mole fractions xi of the prod 3. ucts of combustion of any hydrocarbon fuel and
The products of combustion are assumed to air, at equilibrium under specified conditions of be ideal gases. This assumption is not valid at extremely high pressures (12). temperature T, pressure p and equivalence ratio
HOW TO ACCESS THE SUBROUTINES - Transfer of
F. It can also optionally calculate the partial
derivatives of the xi with respect to T, p and F. information between the calling program and the
subroutines are through labelled COMMON areas
Subroutine PER can calculate the average and hence every program that calls the subrou molecular weight M, gas constant R, enthalphy h and internal energy u of the equilibrium products tines must contain of combustion. It can also optionally calculate the partial derivatives of R, h and u with re spect to t, p and F.
In addition to carbon and hydrogen, the fuel may or may not contain oxygen and nitrogen atoms.
The product species considered are H, O, N, H2,
OH, CO, NO, O2, H2O, CO2, N2 and Ar in gas phase.
LIMITATIONS OF THE SUBROUTINE - The equilib
rium constants used in the subroutines were
fitted as a function of temperature in the range
1080°R to 7200°R. The subroutines cannot be used outside this range.
The subroutine cannot handle the formation
of free carbon. lence ratio,
It can be shown that for equiva
The above variables are input to the subroutine
and should be defined at the time of calling.
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
12
The following are output variables.
IERR - Error code with meaning as specified
tered during calculation of
partial derivatives of mole
later.
XEQ(12) - Mole fractions of the product species.
XPA - The numbering code is given later.
A parameter used in the subroutine.
It represents the modes of fuel that
will give one mole of products.
fractions.
Immediately after every call to the subrou tine, the value of TERR must be checked. It should be zero. If it is not, the numbers
stored in the other output variable locations are meaningless and should not be used. Prefer rably print out the value of IERR and also the
condition (AN,AM,AL,AK,F,T,P,KLO)at which the failure occurred.
AVM - Average molecular weight of products of combustion
R - Gas constant BTU
H - Enthalpy BTU
°R-1
U - Internal energy BTU
The partial derivatives defined earlier
have units of the variable per °R, per psi , per unit fuel air equivalence ratio (dimensionless) respectively. For example,
DXT(1) has units of °R-1
DRP has units of BTU
DHF has units of BTU
psi-1
°R-1
REFERENCE STATES - For computation of en thalpies and internal energies, the following states were assigned zero value of enthalpies.
H2 (ideal gas) at 0°R
O2 (ideal gas) at 0°R
N2 (ideal gas) at 0°R
C (Graphite) at 0°R
Ar (ideal gas) at 0°R
HOW TO CALL SUBROUTINE PER - Both the
$INCLUDE cards and all three of the COMMON state
ments are required. The subroutine PER will be executed every time the following statement is encountered. CONTROL VARIABLE K-]OL KLO = 0
The subroutine will make a fresh
estimate before iterating
KLO = 1
Closeness criteria built into the subroutine will decide whether or not to make a fresh estimate. THIS
IS RECOMMENDED.
KLO = 2
The subroutine will use input values of H2, CO, O2 and N2 mole fractions as initial estimates
ERROR CODE I-]ERR
IERR 0
=
= 1
= 2
No error.
Specified temperature out side 1080 to 7200°R range.
Specified equivalence ratio too high as to cause forma tion of free carbon.
= 3
Attempted to calculate an estimate; but none found within
= 4
outputs of EQMD are also available after a call
to PER.
It is not necessary to call EQMD again.
However if only the mole fractions and/or their
partial derivatives are needed and not the prop
erties, some computation time could be saved by directly calling EQMD as explained below.
HOW TO CALL SUBROUTINE EQMD ALONE - In this case # INCLUDE CO*PER and COMMON/PROP/ cards are not needed. The subroutine EQMD will be executed each time the following statement is encountered.
A singular matrix was encoun
= 5
The subroutine PER internally calls EQMD
with the same value of JDR and all the relevant
mole fractions
Few mole fractions were cal
reasonable limits.
tered during calculation of
culated as negative during iteration. = 6
=7
Newton Raphson iteration did not converge in 25 attempts.
A singular matrix was encoun
SAMPLE PROGRAM - Listing of a sample program which calls subroutine PER is given on the fol lowing page.
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
13
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
14
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
15
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
16
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
17
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
18
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
19
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
20
Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM
21
This paper is subject to revision. Statements and opinions advanced in papers or discussion are the author 's and are
in Transactions. For permission to publish this paper in full or in part,
SAE
been edited by SAE for uniform styling and format. Dis cussion will be printed with the paper if it is published
tion through SAE should send the manuscript or a 300 word abstract of a pro posed manuscript to: Secretary, Engineering Activities Board, SAE.
his responsibility, not the Society 's; however, the paper has
contact the SAE Publications Division.
Persons wishing to submit papers to be considered for presentation or publica
page booklet.
Printed in U.S.A.
References: Downloaded from SAE International by University Of Wales Trinity Saint, Friday, January 17, 2014 09:29:36 AM 2