COMPUTING COMBUSTION EQUILIBRIA ON A HIGH SPEEDDIGITAL COMPUTER
April, 1959
chemisorbed CO toward subsequent physisorption. The initial rapid chemisorption of CO increases the electron work function of the metal considerably so that the subsequent chemisorption of CO occurs with less heat of adsorption. It is probable that the heat of chemisorption may be of the same order as that of physical adsorption. Therefore at - 1 9 4 O , CO adsorbed subsequently will be bound to the surface in a quasi-chemical way without significantly contributing to the magnitude of the positive charge on the metal. A further assumption that the extent of physisorption is the same over the weakly and strongly chemisorbed CO gives an explanation for the result that the presence of weakly chemisorbed CO does not add any more to the suppression already caused by the strong CO chemisorption. Since the transition from physisorption to chemisorption is continuous and the heat of chemisorption a t high coverages may be as low as that of physisorption, it is not possible to visualize any
521
temperature a t which only the physisorbed molecules are removed by evacuation. Nor is there any theoretical justification for assuming tacitly that the amount of such weak chemisorption is inconsiderable in comparison with the amount held more strongly. Strictly speaking, therefore, it would not be possible to estimate accurately the true extent of the CO-chemisorption by the technique of evacuation a t any temperature higher than that of liquid nitrogen and readsorption a t the same temperature. Evacuation a t liquid nitrogen temperature will not be any better, because it is bound to leave behind a considerable quantity of the physisorbed gas. I n view of the uncertainties caused on the one hand by the suppression effect and on the other hand by the desorbability of the chemisorbed gas over a range of low temperatures, neither of the procedures suggested in the paper of Emmett and Skau can be expected to give an accurate estimate of the surface metal sites on cobalt catalysts.
A METHOD OF SUCCESSIVE APPROXIMATIONS FOR COMPUTING COMBUSTION EQUILIBRIA ON A HIGH SPEED DIGITAL COMPUTER1 BY D. S. VILLARS Contributionfrom the Research Department,
U.S. Naval Ordnance Test Station, China Lake, California
Receiued JuZy 18, 1968
A rocedure has been developed for ra idly solving complicated thermodynamic equilibria by a flexible iteration method whici can be readily extended to inclug additional chemical elements. Changes in composition are computed for only one reaction at a time, neglectin the interaction of such changes in composition on the other equilibria. Current values of concentrations are used to calcaate back all equilibrium constants. For the next computation the program selects that reaction showin the greatest discrepancy between calculated and given equilibrium constant. Discrepancies are then recalculated andranother determination made of the reaction showing the greatest discrepancy. The process is repeated until the maximum discrepancy is reduced to a value less than an error, E, specified as a parameter of the problem. Speed of convergence may be maximized by expressing all species in terms of components existing in the largest concentrations at equilibrium. A subiteration procedure is utilized for solving the individual e uations. This converges upon the solution by halving successive tentative intervals. For errors of 0.001% the JPN equilhum at 2 atmos heres is computable with 28 iterations for 2500OK. and 43 for 4000’K. A composition involving excess solid carbon has feen computed with 14 iterations for 2500OK. and 20 for 4000’K. Each iteration involves about 0.2 second machine time.
Introduction proportion to the number of moles at equilibrium of It is well known2 that the specific impulse Ispof a each of the different species present, it is highly depropellant is calculable from flame temperature T,, sirable to have an easy method of calculating the the average ratio of specific heats y , and the average molecular weight M of the gaseous mixture =i
-
( l / g ) [ ( R T e / W I ~ Y / ( Y1 ) ) {1
- ( P ~ / P & v - ~ ) / v ]]’’*
(A)
if equilibrium is frozen during expansion through the nozzle, or from the change in enthalpy a t constant entropy AHs, on adiabatic expansion between chamber and exit 1.p
(-2AHs/M)’h/g
(B)
if the reaction mixture continually re-equilibrates during expansion. In the above equations g is the acceleration of gravity, R the molar gas constant, P e the external pressure, and Pc the chamber pressure. Since the above averages are weighted in (1) Preaented in part before the Division of Physioal and Inorganic Chemistry 131st Meetine of the American Chemical Society, Miami, April 12, 1957. (2) 6. 8. Penner, “Chemistry Problema in Jet Propubion,” Pergsmon Press. New York, N. Y., 1957, see Chapter 14.
equilibrium composition. I n order to discuss a concrete problem let us consider the typical propellant, JPN,* which is a double base ballistite consisting of carbon, hydrogen, oxygen and nitrogen. At flame temperatures in a rocket chamber six equilibria produce 10 reaction species in appreciable concentrations. To solve for 10 mole fractions one requires 10 equations. Four of the 10 necessary equations are atom balances in the four elements and are linear. However, the six equilibrium constant equations are quite non-linear, several of them being cubic and others quadratic. The only practical way of solving such a system of equations is by successive approximations. Such procedures are likely to be very laborious for desk machine calculations but with the recent availability of high speed digital computers which thrive on (3) R. N. Wimpreas. “Internal Ballistics of Solid-Fuel Rockets.” McGraw-Hill Book Co., New York, N. Y., 1950, Table 2-1, p. 4.
D. S. VILLARS
522
multiple repetitions of simple arithmetical operations, the possibility is opened up of writing programs which can complete one of these problems in a few seconds machine time. The calculations of the equilibrium concentrations of gases which result from the combustion of materials of varying composition thus becomes quite feasible. As to method of approximation to be used there is probably an infinite variety available. I n the long run those methods will survive which combine simplicity and speed of operation with accuracy of results. The purpose of the present paper is to disclose an iterative procedure which rapidly gives results accurate to within the accuracy of the given equilibrium constants, yet is so simple that it is readily capable of extension to include additional chemical elements. Space is not available to cite all approximation methods which have been developed for this problem. One important class of schemes involves linearization; ref. 4 one of the most elegant.4 Taylor’s series expansions are made of the equations around differences between the solutions and “guestimates.” On the assumption that the initial guestimates are good these differences will be small and higher powers than the first may be dropped. This enables one to solve the system as a set of linear equations. The writer independently applied the linearization method to the present problem in early 1950. However, although theoretically beautiful, the method is subject to breakdown under two sets of circumstances. For certain types of functions, the guestimates have to be quite good for the method to converge (in other words, higher powers which have been discarded often are not negligible). Also, if the matrix of coefficients approaches singularity, the solution is unstable. Because of these complications the linearization method was not considered foolproof enough for the present work. A different scheme used in this La%oratory was developed by McEwan6 and consisted of setting up a rigorous algebraic solution of the ten equations in terms of a parameter whose value had to be guessed and improved by trial and error. The procedure was coded for the IBM 7016 but the technique was such that a completely new derivation and reprogramming was required for extension to cases involving additional chemical elements. Equilibrium Constant Equations I n order to distinguish between smoky and nonsmoky flame we include
c+;o2=co
(0)
For the reacting elements, C , H, 0 and N, equili(4) H. J. Kandiner and S. R. Brinkley, Jr., Ind. Eng. Chem., 42, 850 (1950). (5) U.S. Naval Ordnance Test Station: “Equilibrium Composition and Thermodynamic Properties of Combustion Gases,” b y William 6. MoEwan, China Lake, Calif., N O T s 26 May 1950. (NAVORD Report 1229,P a r t 1, N O T s 289). (6) U. S. Naval Ordnance Test Station: “A Thermodynamic Investigation of the Product Gases of Carbon-Hydrogen-Oxygen-Nitrogen-Propellant Systems,” b y Mary M. Williams and William S. McEwan, China Lake, Calif., N O T s , 12 January 1955 (NAVORD Report 3421,Part 1, NOT8 1020).
Vol. 63
bria (1) to (6) are generally accepted as being of importance coz = co (‘/do2 (1) H2O = Hz HzO = OH (‘/z)Hz
+ + (’/dOz + (‘/2)H2
(‘/2)02
NO
=
(‘/z)Nz
=
H
(2)
(3) (4)
= 0 (‘j2)Oz
+
(5) (6)
Consider the general chemical equation
I
(C)
bjXj = 0 j
I
where Xj is t h e j t h chemical species, products having positive coefficients bj, reactants negative. According to the principles of conventional physical chemistry, changes in composition necessary to achieve equilibrium may be calculated by a formula of the general type Kp[(S
+ @)/PId
(Xio
=
i
+ bjyYj
(D)
I
I
-
-
1
1
where
I
x
=
cx,o
I
(E)
I
j
(the superscript” referring to initial concentrations) and d = z b j
(F)
I
I
j
The mathematical problem now is: given certain initial values of C, H, 0 and N, and the equilibrium constants a t constant pressure, the pressure and temperature, find the final concentrations of the 10 reaction species that will satisfy all six equilibrium constant equations corresponding to reactions 1 through 6 when there * no excess carbon, or find the final amounts of the 11 reaction species when solid carbon is in excess. A further condition on the solution is that the total sum of atoms of each element is to be maintained unchanged. As stated a t the outset, it is proposed to do this by solving the equations one at a time. In view of the fact that the number of individual types of equations which (D) represents is small, their solution may be conveniently programmed as subroutines. Or, bether, the general type (D) may be programmed as a single subroutine. This allows us to capitalize upon the high speed potentiality of the digital computer which works best for repeated routine operations. Iteration Procedure.-The proposed iterative procedure is to start immediately with the initial composition specified by the propellant under consideration, and calculate the complete change in composition for one reaction a t a time, neglecting the interaction of such changes in composition on the other equilibria. Each time an equation is solved, the amounts listed for the different species are revised to accord with the new results and used for solving the next equation. Different variants of this procedure have been tried. I n the first, a test sum, T, of the absolute values of the changes of all six reactions computed in succession was recorded and the computations repeated again and again until this test sum of changes per round became less than an error, E, specified as a parameter of, the problem. This scheme gave results satisfactory
I
- 1 I I.
April, 1959
COMPUTING COMBUSTION EQUILIBRIA ON A HIGHSPEEDDIGITALCOMPUTER
for specific impulse calculation, but was shown by the referee of this paper to be defective in that it did not converge to correct answers as higher accuracy was ~pecified.~This procedure was modified by determining the fractional discrepancy between calculated and given equilibrium constants Di = (calcd. Kijgiven K,) - 1 and computing that reaction next which shows the greatest discrepancy. After concentrations are adjusted to the new values, the discrepancies are recomputed and that reaction again chosen for solution which corresponds to the greatest current discrepancy. This process is repeated until the maximum discrepancy is less than the error, E. (Each repetition will be referred to herein as an iteration). Such a procedure converges satisfactorily as smaller and smaller errors are specified. Improvement of Rate of Convergence.-In the first iteration scheme used, over two thousand equilibria had to be computed to achieve an (absolute) error of E = for JPN at 2 atmospheres and 2500°K. To speed up convergence, use of an additional reaction, the water gas equilibrium, was suggested by Prof. K. S. Pitzer and ihdependently by L. S. Kassel. The chemical equation corresponding to the latter equilibrium is the linear combination, (2)-(1). Inclusion of this reaction reduced the number of computations of equilibria to 66. The reason that convergence was so markedly affected, as pointed out by L. 5. Kassel, is that (1) and (2) involve an intermediate (molecular oxygen) which at 2500°K. is present in such small amounts that even 1 0 0 ~ consumption o in each reaction computation accounts for such a very small change in CO, COzt Hz and H20that a large number of calculations is required for any progress to be made. The water gas reaction does not involve this intermediate, so the complete change can be effected at once. This same effect shows its importance also in the computation of NO by reaction (6). It became of preponderant importance in calculations involving excess carbon where the oxygen mole fraction is of the order 10-13. It was ultimately realized (during discussions with Prof. Verner Schomaker) that the principle introduced by Kandiner and Brinkley4 to aid the convergence of their linearized equations is just as beneficial to the present iteration method and represents the key to rapid convergence. This is to guess which reaction species are likely to be present in highest concentrations in the final equilibrium mixture (call them "components") and express all other species in terms of reactions of formation from these components. By such a technique there is no question of trying to compute the formation of relatively high concentration species through intermediates existing only in low concentrations. It may perhaps be objected that use of the Kandiner and Brinkley suggestion requires one to "know the answer before he gets it." While success with the linearization procedure often depends upon one providing the problem with fairly close quantita(7) The author is indebted t o the referue, L. S. Kassel, for the suggestion that this may have been due t o internal round-off errors in the fixed point arithmetic used.
523
tive estimates of the final equilibrium concentrations, the present method only utilizes the attributive (qualitative) estimate of which constituents are to be considered as the comporients. Furthermore, the method does not require this information. Its provision is recommended only to reduce the number of iterations required for arriving at the final solution. Indeed, all of the components used below were obtained from preliminary runs which were carried out without benefit of use of components. For JPN at 2500°K. it turns out that the species existing in greatest concentrations are CO, H20, COz and NP. Calling these the components, the reactions of formation of the other species are the following linear combinations J P N 2500°K. Components: CO, HzO, COn, Nz
coz = co + ( ' / z 1 0 2
+ CO = Hz + CO HzO + COz = CO + 20H HZO + CO = COz + 2 H co, = co + 0 COz + ('/z)Nz = NO + CO HzO
(1)
(12) (13) (14)
(15) (16)
On the other hand the components indicated for 4000°K. are CO, H, 0 and N,. The reactions of formation are accordingly J P N 4000'K. Components: CO, H, 0, Na ('/?)HI = H (1/2)0? = 0 Hz0 = 0 2H CH=O+H coz = CO 0 NO = 0 ('/z)Nz
+ +
+
(4) (5)
(22) (23) (15)
(26)
To test a composition giving excess carbon, JPN was modified by adding carbon. Preliminary results indicated the components are CO, Hz, Nz, and C for 2500" and CO, H, Nz and C for 4000". Since Hz was still of the same order of magnitude (0.0626) as H (0.2826) at 4000°,the same set of components was used for both temperatures. The reactions of formation are J P N modified, 2500'K. Components: CO, I&, Nz, C
c+ = co c + coz = 2co
(0)
(1/2)02
C C
+ H20 CO + H, + OH = CO + ('/z)H?
('/z)Hz
=
H
C+O=CO C NO = ('/2)N2
+
(31)
(32) (33) (4)
+ CO
(35) (36)
Results For input compositions and equilibrium constants listed, respectively, in Tables I and I1 results were obtained as given in Table 111. It may be seen that the discrepancies listed at the bottom of Table 111 with which the final mole fractions reproduce the respective equilibrium constants are all less than Specification of the specified error, E = higher accuracy is meaningless in view of the number of significant figures of the given equilibrium constants. It may be seen that use of the Kandiner and Brinkley principle has brought the number of
D.S.VILLARS
524
iterations to achieve this error down to 28 at 2500"K.,43 at 4000" for JPN; and to 14 at 2500" and 20 at 4000" for the soot case. In terms of machine time thede problems are computed in the order of 3 to 9 seconds. TABLE I
VOl. 63
tual calculation of the equilibrium is carried out when carbon is in excess. Fig. - 1.
Input, print: com ositions, equif: consts.
In ut:
+ E, F, $P,
+
R
H2 0 2 N2
JPN
Soot case
0,34523 .23075 .32847 .09556
0.70000 .23075 .32847 .09556
Test excess for t-C Soot
TABLE I1 EQUILIBRIUM CONSTANTS 2500°K.
KPO KPI KP~ Kpa Kp4 KPS KP3
c -
4000OK.
25OOOK.
800-
No
Calc. RO equil.
0.69375 X IO8 4.998 0.60159 .92132 1.5933 1.5528 3.3585
Burn all to co
.1
1
J
0 2
Compute discrepancies route 3
TABLE rII RF~ULTS FOR 2 ATMOSPHERES PRESSURE c-JPN2500'K.
I
soot
4000OK.
0.71531 X IO7 .037818 .0062008 .0047528 .02509 ,015575 16.983
s
i
COMPOSITIONS, MOLEFRACTIONS C
Calcd.: K,/P~,
Compute discrepancies route, R
Find largest discr.
~
1
4000'K.
Equilibrium mole fraction of ges 0.000258 0.022272 0.173 X 10-1: 0.137 X 10-11 .000177 .I63861 ,145 X 10-8 .000001 .072434 .053543 .229434 .061793 ,260696 .008498 H .004775 ,280059 OH .003316 .052885 .483 X IO-' ,446 X 10'8 Ha0 .265538 ,009843 .000007 .170 X 10'6 co .319592 .316061 ,665282 .574568 cox .I92100 .013346 .000003 ,190 x 10-1 NI .I41455 .084707 .096773 .OS3577 . l o 1 X 10-6 NO .000356 .012940 ,241 X 10-8 Residual C (.043081) ( ,043074)
Or 0 Ht
Discrepancies 0 1 2 3 4 5
0.031 .016
.884 .020 .010 ,025
x
106 for error specification, E = 10-6 0.153 .I59 .007 .536 .060 .328 43
0.004 .025 .096 .I22 .869 .185 .064 14
ro RO
Increase I by 1 Compute next equil.
I
0.004 .639 .006 ,057 .498 ,153 ,745 20
Convert to mole fractions Print final results Read in new problem
Three methods of solution of the individual equationa were coded before a practical one was devised. The third method which turned out to be quite Noitms. 28 satisfactory makes use of the knowledge that, in The program was coded so that by means of a order to avoid negative concentrations and still sense switch it is possible to print out the order in have finite positive equilibrium constants, physiwhich the different reactions are selected. In the cally acceptable solutions for x = -y in equation first run they went in the order: 1, 12, 1, 13, 14, (D) have to lie between 15, 16, 12, 1, 16, 14, 15, 13, 1, 16, 12, 1, 16, 14, 15, (Xj0/b,)*R < x < (Xj'/bj)*p 13,1, 16,12,1, 16,14 and 15. Program.-Principal steps of one variant of the where the lower limit of x is the quotient having program are shown in the flow diagram, Fig. 1. smallest absolute value among the reactants (negaOnly general functions are outlined. As indicated tive bj's) and the upper limit is the quotient having above, each solution of an equilibrium equation is smallest absolute value among the products (posiachieved by invoking a subroutine. All calcula- tive b.'s). Furthermore, the value of the right side of (D) changes monotonically as x is changed betions are in floating point arithmetic. Before each computation of discrepancies a tween these limiting values. Half their difference check is made of reaction (0) to determine whether is therefore taken to define a starting value of a comor not carbon is in excess. Since the equilibrium putation variable z and an initial value of x is set (Xjo/bj)*~. Values of Xjo are reconstant for this reaction is of the order 10+6 up to UP equal to x 4000"K., this reaction has been coded for excess placed with Xj' bjy = XjO - b,x and the discrepoxygen on the approximation of complete combus- ancy, Di,is computed (see Iteration Procedure). tion of oxygen and carbon. Such an assumption If Di comes out the next value of x is set to be is valid only for oxidizer-rich propellants. An ac- +z (which in themeanwhile has been divided by 6
++ +,
I ,
April, 1959
CELLSHAVINGA FUSEDZINC CHLORIDE ELECTROLYTE
2); if Di comes out -, x is set to equal -2. The procedure is repeated as many times as is necessary to get \Oil- F < O . In the present work, solutions of the individual equilibria were carried out using the simplifying assumption that y in equation (D) is negligible compared with 8. It is believed that avoidance of this assumption would have shortened the calculations only to a minor extent. It is recommended that F,the subiteration error, be set an order of magnitude smaller than E, else there may be danger of an infinite loop occurring. Furthermore, E should not be specified with too small a value such that it represents digits beyond those carried by floating point numbers. Extension to Additional Elements.-The exten-
525
sion of the present procedure to include additional (chemical) elements is quite simple. The first step, of course, is to decide which are the important reactions. Usually an estimate will be available from previous work as to the most reasonable species to choose as components. Then, after rewriting the chemical equations to correspond to reactions of formation of the other species from the components, one need merely provide for inclusion of the additional data and code in appropriate calling sequences for the new reactions. It is proposed to carry out such extensions in subsequent work. Availability of Program.-The author will be pleased to reproduce a copy of the program for any one wishing to see it.
EXPERIMENTS WITH THE CELL Bi, Bi203,ZnC12/ZnC12/ZnC12,Zn BETWEEN 450 AND 510' BY REUBENE. WOOD National Bureau of Standards, Washington
D. C.
Rsosived J u l y SI, 1968
As a p v t ,of study of cells having a fused zinc chloride electrolyte, experiments were carried out on such a cell in which the negative electrode was liquid zinc and the positive electrode was liquid bismuth on top of which had been put some BiZOa. The cell potentials are quite rapidly and accurately reproducible after temperature changes and after the passage of current. The cell Dotentials and temperature coefficient are reported. Phase equilibria in this system may be rather complicated and have not yet been studied.
Introduction The calomel electrode has long been of interest and utility sts a reference electrode in aqueous electrochemistry, but it is difficult to devise a close, fused-salt analog of this electrode because of the relatively high mutual solubility of the metal halides. As R possible basis for a satisfactory fused-salt reference electrode, similar if not exactly analogous to the aqueous calomel electrode, the use of bismuth and bismuth oxychloride in a suitable fused halide electrolyte suggested itself. The cell that was assembled and studied to test this idea may be represented as Bi, BizOsIZnCln/ZnClz/ZnCln, Zn
(1)
What solid phases are present at the bismuth electrode when this cell is operating reversibly is not yet known. It seems unlikely that they include Bi2O3, more likely that they include bismuth oxychloride and ainc oxide or zinc oxychloride. In fact, although he did not report studies of systems containing ainc, SillBn's work' on two-metal oxyhalide compounds suggests the possibility that reaction between zinc chloride and bismuth and zinc oxides could produce a series of double oxychlorides. However, the above representation is used for the cell to indicate that BizOaand the other substances shown were the only ones put into the cell. This cell was studied over the temperature range 450 to 510" and over a period of several weeks. (1)
L. G . SillBn, Nalurwissenscha/lsn, 80, 318 (1942).
Apparatus and Procedures The cell system used was similar t o the conventional H cell except that there were three vertical tubes and two cross connect,ions between each outside vertical tube and the middle vertical tube. The two lower cross connections provided the electrolyte path between the electrodes, the u per cross connections permitted the passage of helium wiich was maintained as an inert atmosphere above the cells. The pu ose of having three vertical tubes was to allow the use two Zn-ZnCls reference electrodes against which to measure the potential of the bismuth electrode. One of these zinc reference electrodes was in the center tube, the other wah in one of the outside tubes. At the bottom of each vertical tube were ap roximately 4 ml. of molten metal, zinc in two tubes, bismutlin the other. 2.4 g. of B&Oswas laced above the bismuth and then about 100 g. of liquid Zn& was introduced into the cell. A glassenclosed thermocouple junction dipped into the molten metal of each electrode. Electrical connection was made to the molten metals through graphite rods introduced into the cells through side arms so that these rods were not in contact with the electrolvte. Reagent-grade chemicals were used exclusively except that the zinc was an especially pure product asserted by the manufacturer to be 99.99% The zinc chloride was dried by bubbling anhydrous !%?through it at about 500' for 10 hours, then dried oxygen for about two hours. The cell was kept hot in a cylindrical electric furnace. This furnace was supplied with current from a constantvoltage transformer through a variable resistor system. No automatic temperature regulator was used but because of the constancy of power input and of the room temperature, the temperature conditions within the furnace remained quite static &R long as it was undisturbed and the resistance was not changed. Under such conditions, gradual temperature variations at the rate of 0.5" per hour would be typical. All potentials were measured with a T e K potentiometer using a standard cell calibrated in g e s e laboratories. The thermocouples were calibrated only at the melting
2
