Reference

NBS

Publications

'AirWJii—2-s-ferffETr

NAT'L INST. OF STANO & TECH R.I.C.

A11104 ESTfibE

NBSIR 78-1509

An Algorithm and Basic

Computer Program for

Calculating Simple Coal

Gasification Equilibria

... /

William S. Horton

Chemical Stability and Corrosion Division

Center for Materials Science

National Measurement Laboratory

National Bureau of Standards

Washington, D C. 20234

NATIONAL BUREAU OF STANDARDS

a*mJ Wof Standards APR 17 197§

NBSIR 78-1509

» *• AN ALGORITHM AND BASIC

COMPUTER PROGRAM FOR

CALCULATING SIMPLE COAL

GASIFICATION EQUILIBRIA

William S. Horton

Chemical Stability and Corrosion Division

Center for Materials Science

National Measurement Laboratory

National Bureau of Standards

Washisngton, D C. 20234

August 1978

MCl&UaU, OSc U.S. DEPARTMENT OF COMMERCE, Juanita M. Kreps, Secretary

Dr. Sidney Harman, Under Secretary

Jordan J. Baruch, Assistant Secretary for Science and Technology

NATIONAL BUREAU OF STANDARDS, Ernest Ambler. Director
1II '
TABLE OF CONTENTS

Page

Introduction 1

Thermodynamic Basis and Assumptions 4

Thermodynamic and Mathematical Treatment 7

General Statement of the Algorithm 17

Application to a Simplified Coal Gasification System 25

Thermodynamic Data 30

Detailed Algorithm 32

Discussion of the Listing 45

Information for Using the Program 47

The Phase Rule 50

Examples 55

Appendices 53

1. Symbols - Text 53

2. Symbols - Listing 55

3. Listing 57

References 75

An Algorithm and BASIC Computer Program for Calculating Simple

Coal Gasification Equilibria

William S. Horton

Center for Materials Science

National Bureau of Standards

Washington, DC 20234

INTRODUCTION

Designing plants for the gasification of coal is being aided by

testing materials expected to be capable of withstanding the required

conditions. A variety of physical and chemical tests are performed

under these possible conditions and in the same gaseous environments.

Generally, however, to fully understand the observed effects, tests are

made in a variety of environments related to those of the gasification

plants and at a variety of conditions of temperature and pressure. The

gas mixtures used as the environments are generally not at chemical

equilibrium, and depending upon the conditions and the possible presence

of catalyzing substances, chemical reactions will take place [2] at some

unpredictable speed. The tester is concerned about the true nature of

the environment, but this can only be found by carefully arranged chemical

analyses. On the other hand it is of value to know the maximum change

that can occur during the test. This can be determined if the equilibrium composition for the mixture can be calculated.
2

This report is concerned with a simple procedure for calculating

the equilibrium composition of such testing mixtures. The aims in

devising this procedure were to provide an interactive computer program

that was easy to use, not inherently complex, that exploited the relative

simplicity of the BASIC programming language, and was confined to the

coal gasification problem alone. Such a program would be usable by nonprogrammers, who merely need to sign on an interactive timesharing

terminal and type answers to a few questions about the testing conditions.

There are batch programs already available for solving simultaneous

chemical equilibria, e.g., CEC71 [8], however these have many disadvantages

for the problem at hand. They are often complex, handle many more

problems than are needed, require extensive storage, have excessive

print-out, and require some learning prior to use. There is an interesting algorithm [14] with a computer program written in ALGOL, which

treats multiphase chemical equilibria involved in coal combustion. This

includes formation of slag.

The development of a program fulfilling these aims is presented

here. Included are the mathematical and chemical thermodynamic treatment.

Instructions for use and a listing are provided. Also included are some

results obtained using the program which illustrate the use of the phase

rule to predict the potential condensation of solid carbon during a

materials test under coal gasification related conditions.

The program, called COLGAS, is capable of calculating equilibrium

gas compositions containing CH^, CO, CO

2

, and h^O as well as the

quantity of solid carbon and/or liquid water when these are present.
3

Constant temperature is assumed, but the operator may choose to hold

either pressure or volume constant.

The exposition may appear more detailed than is warranted for the

simplicity of the problem to which the development is applied. The

author found this to be helpful in seeing the relationship of the parts

and the exact nature of the assumptions used. The detail has been

retained in the report with the hope that some readers may benefit

thereby.

Wherever possible the recommendations of The International Union of

Pure and Applied Chemistry (IUPAC) for physicochemical symbols and units [3]

have been followed. This includes the SI units. However, the atmosphere, liter, and degree Celsius have been made usable. The SI units

used are the megapascal (MPa) and the kelvin (K). Because of the desire

to keep symbolism the same in the text and in the computer program

wherever possible, occasionally the IUPAC symbol has not been used. For

example, k is used for equilibrium constant rather than K.
4

THERMODYNAMIC BASIS AND ASSUMPTIONS

The literature dealing with calculation of chemical equilibrium

compositions is abundant, and most of it up to 1970 is given by van

Zeggeren and Storey [9]. Other useful reviews are by Zeleznik and

Gordon [13] and Klein [18]. The two known approaches depend upon either

solving simultaneous equilibrium constant equations or minimizing the

Gibbs energy of the system (for a constant pressure problem). Both

include the constraints for conserving the chemical elements. The first

of these methods is rather old, but Brinkley [10] is credited with

developing a generalized scheme for a systematic treatment using equilibrium

constants. This procedure is amenable to the use of computers. The

first reference to minimizing the Gibbs energy as a route to solving

simultaneous chemical equilibria was probably by White, Johnson, and

Dantzig [11]. In principal, however, these two approaches are equivalent.

This becomes obvious if minimization is effected by setting the partial

derivatives with respect to the constituents equal to zero. The result

is a set of equations involving the logarithms of the species concentrations. If logarithms are taken of the Brinkley equilibrium constant

equations a set of similar equations is obtained. With a certain

choice of "constituents" in the Brinkley sense, these sets of equations

are identical. A real difference in the methods can occur if minimization is effected by other methods such as steepest descent or linear

programming [11].

Minimization of the Gibbs energy by use of the partial derivatives

has at least one advantage over the Brinkley method. There is no need
5

to choose substances likely to be present in larger abundance than

others and to treat these differently. All substances in the gas phase

are treated alike and therefore the algorithm can be a bit less complex.

Furthermore, the addition of a new substance is always done in the same

way and does not lead to rewriting the equations already in the algorithm.

Only new equations are added. For these reasons minimizing the Gibbs

energy was the chosen method. Minimization is accomplished by the use

of Lagrangian multipliers and setting the partial derivatives equal to

zero followed by the Newton-Raphson method of solution. Minimization of

the Helmholtz energy for constant volume problems leads to equations only

slightly different from those for constant pressure and to minor additions

to the program.

For the testing of coal gasification plant materials a reference

environment has been chosen [12] and is shown in Table 1.

Table 1. Standard Gas Environment

Volume Percent

CH

4

CO co

2

H

2

H

2

0 h

2

s NH

3

5 18 12 24 39.5 0.5 1

The presence of H^S can lead to S0

2

and the NH^ can give N 2< However, all four of these will be in relatively low concentration and it

was decided to consider the first five species only. If desired, the

concentrations of the others can be simply calculated once these five

are known, as well as the amount of sulfur and/or nitrogen.
6

All gases are assumed to be ideal. Condensed water and carbon are

assumed to be pure. For example, the solubility of carbon dioxide in

water is ignored. For constant volume problems, the volume of condensed

species is ignored in comparison to the gas volumes. That this is

generally acceptable is evident from the fact that when 90% of the

original gas has been equilibrated as condensed species, these latter

contribute only about 0.5% to the total volume in the worst case.
7

THERMODYNAMIC AND MATHEMATICAL TREATMENT

Assume a system of s substances, g of which are gases, at pressure

P, with volume V, and thermodynamic temperature T. If the system is not

at chemical equilibrium, chemical reactions will adjust the number of

moles of each substance, n-, until equilibrium is attained. At this

point the Gibbs energy, G, will have become a minimum if the pressure

and temperature are kept constant (ref 1, p 209). If, instead, the

volume and temperature are kept constant, the Helmholtz energy, A, will

have become a minimum, (ref 1, p 208). The s-g > 0 substances are

condensed phases and are assumed to be pure and immiscible with one

another. The gases are assumed to be ideal and insoluble in the condensed

phases.

The Gibbs energy may be expressed as the sum of the contributions

from each substance.

The p. are Gibbs energies per mole and are known as the partial molar

Gibbs energies or the chemical potentials (ref 1, pp 214,215). For the

g gases in an ideal mixture the chemical potentials are a function of

the partial pressures, p • (ref 1, p 261).

s

( 1 )

p. = p? + RT In p. ( 2 )
8

The p° are constants at fixed temperature and represent the Gibbs

energy of one mole of pure ideal gas at 0.101325 megapascals (one atmosphere) pressure. For pure condensed phases the p. are merely the p°

for these substances. Consequently, equation (1) can now be rewritten.

s g G = H n.p° + RT I n.ln p. (3)

1 1 1 1

1 1

Note that the first g subscripts have been assigned to the gases. The partial

pressures are equal to the product of the mole fractions, n^./n^, anc* the

total pressure.

p.p = (n./n ) P

i i g'

g

n = \ ' n.

g L*t i

1

Thus, equation (3) may be rewritten.

s g

G = X"iM

° + RT XX ln (n

i P/V 1 1

The p. may be evaluated by considering compounds. For example, the

compound A

x B^ may be made, in principle, from its elements, A and B,

according to the following chemical reaction.

xA + yBJ

-> A B x y

(4)

(5)
I 9

The Gibbs energy change for the reversible formation of one mole of

A^By is given by AG° with all the substances in the standard states (ref

1, p 283).

AG

f = m A B

’ x m A

" y m

B x y

But for an equilibrium among these three substances, the AG^ is related

to the equilibrium constant (ref 1, p 283).

AG° = - RT In k.

He re k. is the equilibrium constant for formation of one mole of substance and R is the molar gas constant. Consequently, one can

express the chemical potential of a compound in terms of the absolute

chemical potentials of its elements and the equilibrium constant of

formati on.

m A B = X m A

+ y P

B

" RT 1n k

f x y

For the general case let k. represent the equilibrium constant of

formation for the i - th substance.

'i S V

i J j=l

RT In k. (7)

The summation is over all e elements present in the system, with

v

. . = 0 for any element, j, not present in substance i. Furthermore,

^ J this expression may be used for substances which are themselves chemical elements by recognizing that for these In k. = 0. Using these

ideas, the first term in equation (6) may be rewritten.

RT In k.

i i=l j=l
10

It can be shown that the right hand side here can be arranged as

fol lows.

h Vi = i; n

i v ij- RT E ni inki («)

1 j=l i=l i=l

Because the elements are neither destroyed nor created by shifts among

the n.j, S n_j v — - E j , the number of gram atoms of element j originally

introduced into the system. Consequently, the first term on the righthand side of equation (8) is given by

H E. = G°

j=l J J

where G° is the absolute Gibbs energy due to the elements, a constant.

Therefore, the sum of Gibbs energies for all the substances in standard

states is given by the following.

s s £ n.q? = G° - RT En. In k. (9)

1 1 1

1

1 1

This may be used in (6).

s g

G = G° - RT S n. In ^

+ RT S n. In (n. P/n

g

) (10)

The right-hand member may be rearranged to give a new function,
11

9 s

r G = TT = £ "iln {n

i P/ V ' ?

n

i ln ki <")

I 3

1

Now, because G°, R, and T are constants, minimizing will minimize G.

Turning at this point to the problem of equilibrium at constant

volume and temperature, use is made of the relationship between G and A

(ref 1, p 202).

A = G - PV (12)

Use of equation (10) gives an expression which will lead to one that is

si mi 1 ar to that in (11).

s g

A = G° - RT y n. In k. + RT y n. In (n. P/n ) - PV (13)

i i Ar* i i q

1 1 y

Assuming that the volume of the condensed phases is negligible compared

to that of the gases permits using the (ideal) gas volume for V.

PV = n gRT (14)

At the same time it is possible to replace the variable P.

s g A = G° - RT E n. ln k. + RT Vn. In (n.RT/V) - n RT (15)

l i ' ii g

1 1

In view of equation (5) the last term may be put inside the gas summation and the expression rearranged to define the function =

9 S

r A

= ^1 = E n

{

[In (n. RT/VJ-1]- L n. ln k

i (16)
12

Minimization of Valso minimizes A.

The problem of equilibration has now been shown to require minimization of r

G or PA respectively for constant pressure or constant

volume. However, these minimizations with respect to the n. must be

carried out subject to constraints: 1. The chemical elements are

conserved as previously noted, and, 2. No n- may be less than zero.

The constraints may be expressed as follows.

s

L = E n, ,, - E, = 0 j = 1(1 )e (17)

J

1

1 1J J

n. > 0 i = 1 (1 )s (18)

There are e functions, L . , each of which expresses the conservation of

element j such that the E- equals the number of gram atoms of that

1 element in the system. The E. are determined from the initial composition. As before, . . is the subscript for the element j in subJ stance i. Constraint (18) is handled in this development by replacing negative estimates of gaseous constituents by small positive numbers. A negative n. for a condensed phase assumed to be present indicates the absence of that phase and is retained as an indicator of that

fact.

The method of minimization adopted here is the classic one using

Lagrangian multipliers, A., [6]. To minimize V form the function F by

J multiplying each of the L. by A. and adding to V.

J J
13

F r = Z n. In (n. P/n ) - Z n. In k. + Z Ab 1 1 9 -| 1 1 i J

1 1 In,

F. = Z n. In (n. RT/V) - Z n. In k. + Z A.

" 1

1 l -jiT-iJ1

, 'I v if E

j

£ n. v. .-E.

1

T 1J J

(19)

( 20 )

J

Minimization of F^ will be treated first in detail. It will be seen

later that a few simple changes in the equations will suffice to minimize

F A- Setting the appropriate partial derivatives of F^ equal to zero

results in three sets of equations, one set for the gases, one for the

condensed phases, and one for the elements.

e

f . = In n. - In n. +

g

ln (PA

i ) + £ «ij

= 0

j=l

i=l(l)g (2 1)

f . = In k.. ++ V A •

. vi„

• - — 0

1 2- J 1

J

j=l

i=(g+i)(i)s (22)

f_- = V n. v -

-• " E,

i=l

t Z i v ij j

j=l (1 )e (23)

It can be seen that there are s + e equations, sufficient to completely

determine the s n. and the e A.. The variable n is a defined function

i J 9

which is calculated when the n. have been estimated.

Equations (21) are nonlinear and an iterative procedure, NewtonRaphson [7], is used. In summary, this method iterates from initial

estimates of the n- and the A.- to a final set which satisfies some

J stopping criterion. The iterations may be expressed by use of a matrix

equation [13].
14

( 24 )

Here x and f are column matrices (vectors) representing the unknowns and

the functions (21) through (23). Subscript t indicates that the values

used are from the t-th trial. J is the Jacobian of the f's, a square

matrix with elements equal to the partial derivatives of the functions,

f, with respect to the and the A .

.

There are a variety of types of elements in the Jacobian. Equations

(21) yield the following entries.

3f

.

l n

i i=l(l)g (25)

3f

3n

1 = - 1/n

|i=KDg

jk=l(l)s (26)

' k^i

V

i j

Equations (22) yield the following.

ji=KDg

/j=l(l)e (27)

3f

.

l 3n,

k = 0

(i=g+l(l)s

| k=l (1 )s (28)

3f

.

3A • J

v

i j

j i=g+l(l)s

j j=l(l)e ( 29 )

Equations (23) yield the following.
15

8n

.

i

jj=l(l)e

Vij f i=l(l )s (30)

8f -

k

j— 1 ( 1 )e

k=l (1 )e (31)

Expression (25) gives nonzero values to the first g Jacobian elements

along the diagonal. Expressions (27) and (29) may be combined.

Expressions (26), (30), and (32) show that these elements are symmetrically disposed with respect to the diagonal. Since all other elements

are zero, J is symmetric. This is an advantage for some methods of

calculating the inverse.

These results are easily modified for equilibria at constant volume.

The functions (21) are replaced by (33), but (22) and (23) are retained.

i=l (1 )s

j=l(l)e (32)

e

( 33 )

Of the Jacobian elements, those in (25) and (26) are changed, replaced

by (34) and (35).
16

3f • 57T = 1/n9n^. i i=KDg (34)

af

.

-- = 0 3n

k

i=Ki)g

k=l (1 )s (35)

k*i
17

GENERAL STATEMENT OF THE ALGORITHM

In principle there is no difficulty obtaining equilibrium compositions from the development of the previous section. A problem is defined

by giving values for the constant temperature and the constant pressure

(or volume) and for the composition of the test gas in moles of each

chemical element may be calculated. From the formulae of the substances

be obtained from the literature. One needs only initial estimates of

the n . and the A. . in order to proceed.

1. Choose initial estimates.

2. Evaluate the functions, equations (21) to (23), using the

initial estimates.

3. Calculate the elements of the Jacobian, equations (25) to (32).

4. Invert the Jacobian.

5. Calculate the new estimates from equation (24).

6. Check the estimates against the stopping criterion.

7. Depending upon the result of step 6, either stop or repeat

steps 2 through 6 using the new estimates for evaluation

purposes.

* In the later section about using the program the instructions ask for

mole fractions or percentages although these data are treated as moles.

Experience with COLGAS has shown that when n is much less or much

9

greater than one, difficulty with over- or underflow may occur.

substance.

the v.. are also determined. The equilibrium constants of formation may
18

As straightforward as this procedure is, there are a number of

practical points associated with these steps which are worth commenting

on. However, the handling of the practical problems depends upon the

general aims of the program. Therefore, although some general remarks

will be made, it must be understood that the particular procedure

adopted in this report is based upon the restricted nature of the

problems being solved.

Choosing the Estimates . The program to be developed was to give a "one

shot" solution. That is, a single run would produce only a single

answer. A different temperature and/or a different pressure (or volume)

or test gas composition will require a new run. Consequently, the

program has no prior experience to build upon as is used in some programs

[8, 14]. It was decided to use the test gas composition as the initial

estimates for the substances, and this has turned out to be satisfactory.

For each A., zero has been a satisfactory initial estimate.

There are two versions of the program in the same listing. The

commonly used version has a short printout. Either the short or the

long printout is obtained by an appropriate response to the initial

query at the terminal. The long version permits the operator to suggest

which of the four combinations of condensed species situations is believed

to be correct. It then allows input of a set of initial estimates

different from that used automatically by the program. These alternatives

have been provided so that on occasions of non-convergence an operator

may try to redo the iterations with a better starting situation.
19

Evaluating the Functions . Equations (21) to (23) present no problems

beyond obtaining values for the equilibrium constants of formation for

compounds in the system. (The total number of moles of gases, n , is

evaluated anew at each iteration using equation (5)). The equilibrium

constants used in the program, except for liquid water, come from the

JANAF Tables [15] and are discussed in a separate section.

Calculating the Elements of the Jacobian . Examination of equations (25)

to (32) shows that the only elements varying in value from iteration to

iteration are the elements of the upper left 5x5 submatrix. These will

require recalculation during the iteration. All other elements are

constant integers and need not be re-evaluated. However, if the presence

or absence of a condensed phase turns out to be different than assumed,

the presence or absence of the appropriate constant elements must be

adjusted.

Inverting the Jacobian . There is an extensive literature on the inversion

of matrices, particularly symmetric ones [16]. It is tempting to seek

a direct analytical solution with perhaps the possibility of avoiding

some round-off error in the computer. In fact, at constant volume and

for the case where solid carbon is present, but liquid water is not, it

is relatively simple to do this using either an LP or an LDR decomposition

of the Jacobian [17]. There is a minor advantage to this beyond the

possibility of reducing the round-off error. Some of the expressions

for the n . contain functions of them which may be replaced by the Ej

constants for hydrogen and oxygen. This should improve the accuracy of
20

the estimates and perhaps reduce the number of required iterations.

However, when the condensed phase situation is reversed, some of the

expressions for the estimates become quite complicated and hardly worth

evaluating directly. Although not essential, it is desirable to handle

all varieties of condensed phase situations in the same way for simplicity of programming. This is particularly true when using the BASIC

programming language because of the one-statement matrix inversion

procedure provided. A direct analytical solution was, therefore, not

used.

Calculating the New Estimates . As is well known, the Newton-Raphson

procedure often overcorrects in the early stages, particularly if the

initial estimates are somewhat far from the mark. This results in

divergence. To avoid this a common procedure is to modify the correction term in the estimating equation, (24), by a weighting factor, k,

less than unity.

x

t+i = xt

k

t V ft (36)

In many nonlinear equation problems it is adequate to use k = 0.5 in the

early stages and to change it to unity later on. In running this

program it was found that this was not always sufficient. This difficulty was overcome by making k a variable function of t, the number of

the trial. A sinusoidal function of t sets k equal to about 0.01 for t

= 1, to 0.5 for t = 8, and to 1.0 for t = 16. The value of k remains at

unity for any further iterations. This is done only for the gaseous
21

constituents. Experience with the algorithm has shown that the estimates

of the Langrangian multipliers converge easily with k = 1.0 at all

stages.

Sometimes new estimates of the gaseous moles become negative. This

not only violates the constraint, equation (18), but it cannot be tolerated when evaluating the functions, f . , which contain In n.. Trans2 formations may avoid this problem, y

= n [4] or y

= In n [5], but introduce numerous back and forth transformation steps in the algorithm which

were thought to be undesirable. The algorithm reported here replaces

such numbers with a very small number, originally 1 x 10 ^

.

It was

found that this occasionally led to difficulty when the size of the

correction was larger than this figure. A loop of repeated values could

result with continually negative estimates for the same n- leading to a

lack of convergence. This difficulty was partially overcome by replacing

the 7 in the exponent by a random number distributed uniformly between 7

and 10 inclusive. This choice occasionally leads to a new difficulty

with some processors. During the matrix inversion step, usually a

Gauss-Jordan routine, the product of the five (1/n.j-l/n ) is evaluated.

If for a given iteration a number of n^ estimates become negative the n^

product may cause an overflow. This has happened with a compiler for

which the exponent of numerical constants is restrained to be within ±38.

For a processor with a larger range, say greater than ±55, this difficulty

may not arise.
22

Checking the Estimates for Stopping . The criterion used here for

stopping depends upon the magnitude of the unweighted correction to the

n. and the A.. The vector c. is the correction and it comes from (24).

i J t

c t

= a

1 f t

The iterations are considered complete when

c —t -5 — < 5 x 10

3

(38)

n t

for all n.j and Aj. This particular stopping criterion has been criticized, (ref 18, p 520) as not leading necessarily to the same solution

as a criterion based upon the derivatives: f^ being sufficiently close

to zero. However, because the interest here is in the gas composition

the criterion given above is preferable. To assure adequate minimization, the pertinent independent chemical equilibria constants are

checked for agreement when calculated from Gibbs energies of formation

and from the calculated mole fractions. According to Brinkley [23]

there are at most R = N - C reactions for N chemical species and C components. In this discussion C = 3 and N equals 5, 6, or 7 depending

upon the case. For case IV, all gas, (see below), N = 5 and R = 2.

The reactions chosen were the following.

A: C0

2 + H

? = CO + H2

0

B: CO + 3H

2 = CH

4

+ H

2

0

When C(s) is present, reaction C is included, and for liquid water

reaction D.
23

C: CG

2 + C(s) = 2C0

D: H

2 0(1) = H

2 0(g)

Terminating the Iterations . When the stopping criterion, equation (38)

has been satisfied for all n. and A. one set of iterations is complete.

J The problem may not be solved, however. In order to perform the calculations a choice had to be made among the possible combinations of

condensed species situations as shown in Table 2.

Table 2. Condensed Phase Cases

H 2 0(1)

YES NO

YES 1-8 II - 7

C(s)

NO III - 9 IV - 8

Here the yes-no tabulations refer to the presence or absence of the

phase at the top or at the left. The Roman numerals are used to identify

the cases, and the Arabic numerals show the number of functions, f.,

and correspond to the order of the relevant Jacobian.

Solving the seven-equation system correspond!' ng to case II, say,

has assumed the presence of solid carbon and the absence of liquid

water. However, this assumption may not be correct. Therefore, after

going through steps 1 to 7 and deciding to stop, it is necessary to

check the assumption of case II. If liquid water is really present
24

and/or solid carbon is really absent as suggested by the final n. , then

the appropriate case must be chosen and the iterations redone. There

are, therefore, two hierarchies of iterations, one on the n. and one on

J the cases. When the final n. are found to be consistent with the case

assumed, the calculation is finally finished. Details of checking case

consistency are covered in the section on DETAILS OF THE ALGORITHM.
25

APPLICATION TO A SIMPLIFIED COAL GASIFICATION SYSTEM

The work shown in this report was undertaken to provide results

with atmospheres used for testing materials to be used in coal gasification

plants. For simplicity the problem was limited to one in which an

initial gas mixture was made from two or more of the following: CH^, CO,

CC^, and F^O at some temperature, T. The test may be carried out

dynamically in which case constant pressure is assumed. A static test

implies constant volume. Constant pressure will be assumed here. From

what has gone before, extension to constant volume will be obvious.

The first step is to set down the pertinent functions, according to

equations (21) through (23) of the second section, but it is necessary

to decide upon allowable condensed states. Because conditions of test

will not permit the existence of ice, only liquid water need be considered.

The conditions will be above the critical temperatures of the other four

substances so no other solid or liquid phase for these need be considered.

Flowever, decomposition of CH^ and disproportionation of CO can produce

solid carbon. Consequently, there are two condensed phases to consider,

liquid water and solid carbon. Examination of the simultaneous equations

in the light of this shows that n , the number of moles of solid carbon,

appears only once, and that in the equations (23) corresponding to that

element. Allowing indices one through five to refer to the gases in the

order mentioned above

n n + n 0 + n 0 + n_ - 12 3c C = 0 (39)
26

Here C is the E. for the element carbon. It is unnecessary, therefore,

J to solve for n

c iteratively. Instead the equation for the constraint

due to the element carbon will be deleted. When the other n. have all

been determined, n

c may be calculated from (39). For case II (solid

carbon only) the functions to be used are the following.

f , = In n, - In n - In k n + 4A,

1 1 g 1 h

f 0 - In n 0

- In n - In k„ + A

2 2 g 2 o

f, - In n. - In n - In k 0 + 2A

3 3 g 3 o

f 4 = 1n n4

" ln n g

~ 1n k

4 + 2*

h (40)

fr = In n c

~ In n * In L + 2A, + A

5 5 g 5 h o

f 6 = 4n

] + 2n

4 + 2n5

“ h

f

7 ~ n 2

+ 2n

3

+ n 5

' 0

Here h and o are the E. for hydrogen and oxygen.

Equations (25) to (31) of the second section lead simply to the

Jacobian.
27

J =

1 _ J_

n, n

1 9

1 1 n 2

n g

1 1 n 0 n

3 g

1 1 n- n 4 g

1 1 n c n 5 g

-1

0 |(41)

The functions of equation (40) and the inverse of the Jacobian (41)

are the needed entries for the right hand side of equation (24) of the

second section.
28

The development of case II gives the smallest set of equations to

be solved by iteration and give rise to the "basic" Jacobian. This is a

7x7 upper left submatrix common to the cases. Consequently, the elements

in it may be defined without regard to case.

The development here was for a particular case. The one with

condensation, Case I, with both condensed phases present, follows by

adding terms for liquid water to fg and f^: 2n^ and n, respectively. A

function for liquid water is needed.

The pertinent Jacobian is equation (41) with another row and column

added. The row is given by jg.

The column is the transpose of this.

Case III, with liquid water but no solid carbon, requires fg

and n but has no n . Because there is no n equation (39) becomes

w c c

This function must be included as well as a A . A new row and column

c are added to the Jacobian consisting of ones in the first three elements

and zeros elsewhere. This is the largest matrix and has nine rows and

(42)

j 8 - {0,0,0,0,0,2,1,01 (43)

(44)

columns

.
29

Case IV, with neither liquid water nor solid carbon, is derivable

from Case III by elimination of the eighth row and column. These correspond to liquid water and A^. Also the function, fg, is deleted. The

problem of determining equilibrium compositions in the system is then

resolved to choosing the appropriate form of the Jacobian, (41), depending upon the case. Deciding the appropriate case for any given problem

without any prior knowledge is difficult. As was seen in the previous

y section, the algorithm given here assumes a starting case and checks the

results for consistency with this assumption.
30

THERMODYNAMIC DATA

For the five gases, CH^, CO, CO^, H

2

, H

2

O, the tabular data of the

JANAF Tables [15] between 300 and 1400 K were fitted by least squares to

equations of the form

m

In k. = a^ + a

2

In T + X] a i+2 (45)

Here, as stated earlier, k. is the equilibrium constant of formation of

the i-th compound in its standard state and from its elements in their

standard states. The constant for hydrogen is, of course, unity so that

the corresponding In k. equals zero. The value of m was determined by

trials with m = 2(1)5. The smallest m > 2 so that the standard error of

the residuals < .001 was chosen. The form of this equation is similar

to that of Barron, Porter, and Hammond [19] but the latter authors used

a fixed highest power equal to five. The least squares coefficients

were determined with the computer language OMNITAB [20] which has an

accurate least squares routine [21].

Liquid water is not included in the JANAF Tables, so resort was had

to combining the JANAF data for gaseous H

2

O, k^, with the vapor pressure,

?s1n k

£ = 1n k5

" ln p

s (46)

For the vapor pressure data, Schmidt's tabular values [24] were used and

fitted by least squares to an equation of the following form.
31

m

In p

s = b1 + b

2

/T + E b i+2T

1 (47)

i=l

In this case the data were sufficiently detailed that a tighter criterion

could be used for m. The smallest m was chosen for which the maximum

residual < 0.0005. The results, In for the gases and In p

& for

H

2 0(£) are shown in the table.

Table 3. Empirical Coefficients for the In and Ln p

$

In k . a

i a

2

a

3 xl0

3

a^xlO^

9

a^xlO3 in 13 a^xlO

:H4

8212.35 -0.816063 -10.0801 7.73804 -3.06981 5.08552

:o 13632.3 1.8055 - 2.3936 0.372836 — —

,°

2

47326.7 0.0777872 - 0.280706 0.0348747 — —

H 2 °(g) 28786.6 -0.697311 - 1.45957 0.83769 -0. 180286 —

b

l b2

b

3

b^xlO

4

brXlO

7

5

, n

l 0 b^xlO

! " p

s 55.5049 -9430.74 -0.159684 3.07499 -3.0549 1.24325

1 The maximum residual for any of the nonzero ln k. to be used is

0.00186 for ln p . Program trials with In k differing by 0.01 showed

absolute values of changes in the computed n-j < 0.0001. The accuracy

of the coefficients in the table is quite adequate for the purpose.
32

DETAILED ALGORITHM

Master Flowchart . The details of the algorithm will be described in

terms of flowcharts for some sections and with an overall chart showing

the relationship of the pieces. It will be noticed that although some

sections of the program could be subroutines, they are not. There are

no "gosub" and "return" statements. There is no advantage to having

formal subroutines here. Replacement of a section is no more difficult

than replacing a subroutine because BASIC statements must all have line

numbers which must be adjusted in either case. Furthermore, because one

cannot fruitfully use dummy variables, one of the advantages of a

FORTRAN type subroutine is lost.

The master flowchart is shown as Figure 1. Most of the items are

self-explanatory, however, some remarks about others will be useful.

The first major branch point reflects the impossibility of H

2 OO)

existing above its critical temperature, 647.3 K, and below that at

total pressures less than its vapor pressure. This is a useful distinction

because it rules out two of the four cases, I and III. The choice of II

and IV as starting cases for the respective branches resulted from the

experience that these led more often to convergence of the first set of

iterations. A longer experience may suggest alternate choices. Note

that the variable c-j is assigned the value of the case number assumed at

the start of a set of iterations. This is required for the case-consistency

checking section.

The Newton-Raphson iterations follow in either case. These will be

treated later by a separate flowchart. Upon successful convergence of
33

Colgas

0

Fiaure la
34

Figure lb
35

the iterations, the next step is to check the case suggested by the

results. Upon unsuccessful convergence, the " lack-of-convergence

output" results. The latter not only states lack of convergence but

also gives the case number that had been assumed and the last set of

iterated results. These outputs are useful if it is desired to try

again with a different set of initial estimates.

The case-checking section will be presented in a detailed flowchart.

Essentially it computes the case number, c^, suggested by the successfully converged results and compares this with c-j . If they are equal

then a "successful output" results which gives the equilibrium gas

composition, quantities of liquid water and/or solid carbon if either or

both are present, and the changed volume (or pressure). If c

^ does not

equal c^ then a counter is compared with the number four. If four cases

have been tried then an "unsuccessful output" is produced stating this

fact. This kind of result has not yet appeared in practice (April 11,

1978), however, the possibility has been allowed for. If four cases

have not yet been tried c-j is set equal to c? and the program returns to

that case for another set of iterations.

Newton-Raphson Iterations . The flowchart for this section appears in

Figure 2. The iteration counter, t-j , allows up to 30 trials. It has

been the usual experience with this program that if convergence has not

been reached in 30 iterations, it will not be reached by going beyond

this point. Commonly, convergence is attained with between 10 and 20

trials. It would take less, but, for the cautiously slow start enforced

by the weighting factor, k.
36
37

Figure 2b
38

I

Figure 2c
39

Figure 2d
40

Initially, the seven functions which, with slight variations, are

common to all four cases, are evaluated using the previously determined

n.. The variations required for the different cases are achieved by

using two constant multipliers b and d, which are set independently

equal to zero or unity depending upon the case. A third constant, e, is

the subscript for A , the Lagrangian multiplier for the element carbon

in cases III and IV when solid carbon is not present. Since A is

always last in these cases, e is set equal to n, the size of the appropriate

Jacobian. Next, via an "on c^" statement, the appropriate remaining

functions for cases I, III, and IV are assigned. This is accomplished

with the user-defined BASIC functions fnw and fnc.

Following the assignment of the functions, the elements of the

upper left 5x5 submatrix of the Jacobian are evaluated. As stated

earlier, these depend upon whether the problem is for constant pressure

or constant volume.

Once evaluated, the Jacobian is inverted and unweighted corrections

calculated. The next part of the diagram introduces the weight factor,

k, which is applied only to the five gases and then up to only trial 15.

Following the correction of all n^ , the number of moles of gas is

calculated.

The next operation is a loop testing all n^ against the stopping

criterion. If successful for all, s, the counter of iteration sets is

incremented by one and the calculation proceeds to testing the condensed

phase case found against the one assumed.
41

Case-Consistency Check . In Figure 3, this section begins by calculating

mole fractions for the five gases. These are used to calculate equilibrium constants for the three independent reactions and the partial

pressure of water. As stated earlier these are used to check the final

solution found by minimization of the reduced relative Planck function.

(See the subsection on Checking the Estimates for Stopping.) Although

only k^, for the di sproporti onation reaction, is used in this section

they are all calculated together.

The determination of the case found by the converged iterations

depends upon the case that was assumed. If that was case III or IV then

the conservation of the element carbon was enforced. In these circumstances, by comparing the equilibrium constant for disproportionation

(chemical equation C) calculated with the mole fractions with that

calculated from the Gibbs energies of formation the need for solid

carbon can be determined. On the other hand, if carbon is assumed

present then the defined function fnc indicates whether its presence has

been verified.

For liquid water cases I and III go together by assuming its presence.

For those cases the sign of n^ is the indicator. For cases II and IV,

n-j is not calculated, therefore, the partial pressure of water is checked

against the vapor pressure.

The value of c

2

is built up by an initial assignment and addition

as the program works its way through this section. Upon reaching its

final value, c

^ is compared with c-j as stated earlier. Before returning

to the appropriate case, if that is necessary, c-j is set equal to c
42

Case- consistency

check

k4 = *2 *4

k§ = x s x s / x 2 x| p

2

\ d! = p-exp|i7 )

Figure 3a
43

I 1 hf

Figure 3b
44

Also it is necessary to clear the eigth and ninth rows and columns of

the Jacobian preparatory to inserting the appropriate values for the

next case.
45

DISCUSSION OF THE LISTING

Appendix 1 gives the list ‘of symbols for the mathematical development in the text and Appendix 2 is for the source program. Wherever

possible, the same variable symbol has been used in the text and in the

computer program. To avoid cluttering the text, variables used there

often have fewer subscripts than the corresponding symbols of the

program.

The listing, which is Appendix 3, is given in the full ASCII

character set. As stated earlier, it is dependent upon the BASIC

* compiler used by the Computer Sciences Corporation. At this writing

(April 11, 1978) there is no standard for the BASIC language as there is

for FORTRAN [22]. On the other hand almost all BASIC compilers provide

reasonably useful diagnostics and internal editing commands so that

translation to another compiler should not be a great chore.

Use of the full ASCII set permits use of upper and lower case

letters. This has been used to differentiate vectors and matrices.

The former, even when the unity dimension is explicitly indicated, are

given by lower case letters, as are unsubscripted variables. Matrices

are in upper case as are elements of matrices.

Sections are separated by remark statements with leading single

quote characters. Some statements are amplified by remarks after a

* A particular commerical source is identified in this paper in order

to explain the procedure adquately. In no case does such identification

imply recommendation or endorsement by the National Bureau of Standards,

nor does it imply that the service is necessarily the best for the

purpose.
46

single quote character. Loops have been doubly indented in order to

catch the eye, and inner, nested loops are each indented again. The

listing is packed, i.e., blanks have been omitted in the program statements except for explanatory matter and output text. The experienced

BASIC programmer reads this with ease, and a some storage space is

saved.
47

INFORMATION FOR USING THE PROGRAM

Input . After signing on at the terminal the program COLGAS is LOADed

and RUN. The following queries will result.

The program first asks whether the problem is to be for constant

pressure or volume and whether a short or long printout is desired. The

long printout is generally only used when attempting to rerun a nonconverging problem. The answer to these joined questions is therefore

generally "p,s" or "v,s". For the long printout "s" is replaced by "1".

The next question concerns the temperature and pressure (or volume).

The answers, again separated by a comma, are in kelvins and atmospheres

(or liters). Alternatively, degrees Celsius and/or megapascals may be

used if the values are preceded by an asterisk without separation by

commas or blanks.

The concentrations of the assumed gases are requested next, and the

response is to be in mole or volume (decimal) fractions. The answers

are typed in the order CH^, CO, C0

2

, H

2 0 with separating commas. If

any of these are not present, zeros must be entered in the appropriate

places. If the input concentrations do not satisfy certain requirements,

the operator will be requested to try again.

The computer program does not consider problems in which all three

elements, C, H, and 0 are not present, i.e., if CH^ + H

2

, H

2 + H2 0, or

CO + C0

2 are used as starting mixtures. In the first two cases only

decomposition of the compounds can occur; in the third only disproportionation. (It is assumed that the decomposition of H

2

and CO are

negligible.) Such equilibria are calculated by hand directly and
48

without iteration. It is also true that a mixture of CO^ and F^O does

not undergo appreciable reaction.

At this stage, input for the short output is complete. If the long

output was chosen the following extra input is requested.

The opportunity to choose the starting case will be next presented.

The response is "yes" or "no". If "yes" is typed, the (Arabic) number

of the preferred case is requested.

Finally, the operator is to decide if a preferred set of initial

estimates is to be used; "yes" or "no". If "yes", the preferred values

are requested and these are typed in the same fashion as the assumed gas

concentrations. Estimates for only the five gases are to be entered.

This is the end of the operator input.

Output . All input data are reprinted as an echo check. After the

input gas concentrations have been echoed, the unchosen quantity, volume

or pressure, is given. Because the number of moles of gas almost

always changes in the final equilibrium, all extensive quantities are

quoted on the basis of the number of moles of the assumed gas. Although

pressure is not an extensive quantity it too can only be quoted on this

basis when a constant volume problem is calculated because the number of

moles of gas in the input volume has not been requested.

A long printout gives at this point the assumed case and the results

of each iteration until the program terminates successfully or unsuccessfully. In the latter situation a statement about unsuccessful convergence

is followed by the last set of estimates. In any case the mole fractions
49

are next and these are followed by the number of moles of gas and the

new volume (or pressure).

If liquid water or solid carbon are produced a pertinent statement

is made at this point. If the temperature is above the critical value

for water, such a statement is made next.

Normally, that is the end of the printout. However, if the run was

unsuccessful the independent equilibria will not have been satisfied and

statements to this effect will be the last output.
50

THE PHASE RULE

For the restricted sets of problems considered here, the phase rule

is simple to apply. It is given by equation (48).

F = C - P + 2 (48)

The number of degrees of freedom, F, is the number of independent

choices of variables possible with a given system which has P phases and

C components. In order to use the traditional symbols for the phase

rule, P has been used here, although it has been used before. Because

this applies only to this section, confusion should not result. There

are three possible phases, gas, liquid water, and solid carbon. The

number of components in these problems is equal to the number of chemical

elements, three. This simplifies equation (48).

F = 5 - P (49)

In case IV (see the section on General Statement of the Algorithm)

there is only the gas phase, and there must be four variables available

for choice. Since two are assumed to be temperature and pressure (or

volume), there are two left for concentrations . Thus, if a triangular

diagram is used for the concentrations of the elements carbon, hydrogen,

and oxygen, the possible elemental concentrations would be represented

by an area. See Figure 4. If two atomic compositions are chosen, the

third is fixed by the requirement that the percentages add to 100.

On the other hand for an equilibrium that has either solid carbon

or liquid water (but not both) present, there is only one degree of
Figure
52

freedom, and this must be represented by a line. Choosing the percentage

for just one element fixes both of the others if one condensed phase is

present. Also, because the condensed phase always has a fixed composition, the variable composition refers to the gas phase. Referring

to Figure 3, the curved line shows the gas compositions in equilibrium

with solid carbon at 700 kelvins and 1 atmosphere. The curve has been

drawn by connecting points calculated with COLGAS which turned out to be

case II. If a gas mixture contains 15 atomic percent carbon and is in

equilibrium with solid carbon the the gas could contain also about 48

atomic percent hydrogen and 37 percent oxygen. Alternatively, it could

contain 81 percent hydrogen and 4 percent oxygen. The two possibilities

are at the intersections of the 15 percent carbon line with the calculated

curve.

A valuable use of the phase diagram in Figure 4 is to show potential

carbon precipitation for an assumed test gas mixture. One needs merely

to calculate the elemental composition from the compound composition and

locate that point on the diagram. If it falls in the area above the

curve, carbon will precipitate. Furthermore, extending the line through

the apex (for c = 100%) and the point to the curve will give the elemental,

but, not the compound composition of the gas in equilibrium with the

carbon. It should be noted that gas mixtures made from the five postulated

substances will have elemental compositions within the region bounded by

connecting the compositions of all of the substances with straight

lines. The lines should connect CH^, CO, CO

2

, F^O, F^, and return to

CH^ in that order.
53

Figure 4 represented a situation above the critical temperature of

water. Figure 5 illustrates the situation for 500K and 50 atmospheres.

The diagram is rather more complicated. Although not always indicated,

a gas phase is present for all regions of the diagram. There are two

invariant points and two regions within which gas compositions will

change to one or the other resulting in condensed water and carbon.

There are two regions within which the change will result in carbon

only, and two more from which only water will condense. Note the two

small regions within which no condensation will occur, although five

gaseous substances may change their relative amounts.
*3" ln < <

00 «-> <

GAS

ONLY
55

EXAMPLES

Examples A through D are representative of cases I through IV and

are successful short printout results. Example B, for 700K and 1

atmosphere, corresponds to the conditions for the phase diagram of

Figure 3. The original elemental proportions are .6:1. 6:. 8 and the

atomic fractions are .2, .533, and .267. The line extended from the

apex through this point on the diagram intersects the curve for solid

carbon at approximately .13, .58, and .29. The equilibrium mole fractions

for the five substances give for the atomic fractions in the order C,

H, 0: .1288, .5809, and .2904.

Example E is for a problem which did not converge during the shortprintout run. The last set of estimates suggests that the difficulty

involved the amounts of CO and CO^ which appeared to be continually

going negative. Often in a situation of this kind it was found that the

initial estimates for such substances were not small enough. Consequently,

for the long printout run, example F, initial estimates were chosen

equal to the last set from the short run except for the CO and CO

2

which

were made smal ler.
56

Example A

col gas 13:14 04/07/78

Is this a constant pressure (p), or a constant volume (v) problem, and

do you want a short (s) or a long (1) printout?p,s

What are the temperature (kelvins) and pressure (atmospheres)?500 , 1 00

Temperature = 500 K = 226.85 C Pressure = 100 atm = 10.1325 MPa

What are the assumed concentrations in mole fractions of the gases in

the order: CH^, CO, CO^, H^, and H^O? . 1, . 1, . 1

, . 4 , .

3

.1 .1 .1 .4 .3

The initial volume = .410284 liters per mole of original gas

Gas contains C, H, and 0 atoms in the proportions .3 , 1.8 , and .6

Equilibrium mole fractions in the order CH^, CO, CO

2

, H^O are: .7301 .0000 .0077 .0016 .2606

The number of moles of gas = .2074 per mole of original gas

The volume = 8.50928e~2 liters per mole of original gas

Solid carbon is present in the amount of .1470 moles per mole of original

gas

Liquid water is present in the amount of .5428 mole per moles of original

gas now at end
57

Example B

colgas 13:16 04/07/78

Is this a constant pressure (p), or a constant volume (v) problem, and

do you want a short (s) or a long (1) printout?p,s

What are the temperature (kelvins) and pressure (atmospheres)?700 ,

1

Temperature = 700 K = 426.85 C Pressure = 1 atm = .101325 MPa

What are the assumed concentrations in mole fractions of the gases in

the order: CH

4

, CO, C0

2

, H

2

, and H

2

0?. 2, . 2, . 2, . 2, .

2

.2 .2 .2 .2 .2

The initial volume = 57.4398 liters per mole of original gas

Gas contains C, H, and 0 atoms in the proportions .6 , 1.6 , and .8

Egui librium mole fractions in the order CH^, CO, C0

2

, H^, H

2

0 are: .1693

.0076 .2344 .1376 .4511

The number of moles of gas = .8626 per mole of original gas

The volume = 49.5486 liters per mole of original gas

Solid carbon is present in the amount of .2452 moles per mole of

original gas

Above the critical temperature for water

now at end
58

Example C

col gas 13:18 04/07/78

Is this a constant pressure (p), or a constant volume (v) problem, and

do you want a short (s) or a long (1) printout?p,s

What are the temperature (kelvins) and pressure (atmospheres)?500 ,50

Temperature = 500 K = 226.85 C Pressure = 50 atm = 5.06625 MPa

What are the assumed concentrations in mole fractions of the gases in

the order: CH^, CO, C02> H^, and h^O?. 05, . 1, . 05, . 5, .

3

.05 .1 .05 .5 .3

The initial volume = .820569 liters per mole of original gas

Gas contains C, H, and 0 atoms in the proportions .2 , 1.8 , and .5

Equilibrium mole fractions in the order CH^, CO, CO^, are: .4730

.0000 .0012 .0047 .5211

The number of moles of gas = .4218 per mole of original gas

The volume = .34611 liters per mole of original gas

Liquid water is present in the amount of .2792 moles per mole of original

gas now at end
59

Example D

col gas 14:22 05/03/78

Is this a constant pressure (p), or a constant volume (v) problem, and

do you want a short (s) or a long (1) printout?v,s

What are the temperature (kelvins) and volume ( 1 i ters )? 1255,1.514

Temperature = 1255 K = 981.85 C Volume = 1.514 1 iters

What are the assumed concentrations in mole fractions of the gases in

the order: CH^, CO, CO

2

, , and H^O?. 05, . 18, . 12, . 25, .

4

.05 .18 .12 .25 .4

The initial pressure - 68.0194 atm per mole of original gas

Gas contains C, H, and 0 atoms in the proportions .35 , 1.5 , and .82

Equilibrium mole fractions in the order CH^, CO, C02 # E^O are: .0169 .1902 .1218 .3342 .3368

The number of moles of gas = 1.0640 per mole of original gas

The pressure = 72.37 atm = 7.33289 MPa per mole of original gas

Above the critical temperature for water

now at end
60

Example E

colgas 14:14 05/03/78

Is this a constant pressure (p), or a constant volume (v) problem, and

do you want a short (s) or a long (1) printout?p,s

What are the temperature (kelvins) and pressure (atmospheres)?500,50

Temperature = 500 K = 226.85 C Pressure = 50 atm = 5.06625 MPa

What are the assumed concentrations in mole fractions of the gases in

the order: CH

4

, CO, C0

2

, H

2

, and H

2

0?. 12,0, .01 ,.75,. 12

.12 0 .01 .75 .12

The initial volume = .820569 liters per mole of original gas

Gas contains C, H, and 0 atoms in the proportions .13 , 2.22 , and .14

30 iterations did not converge!

The results at this stage are as follows

Case IV

The last set of estimates is, moles and multipliers:

.13 l.e-9 l.e-9 .71 .141.79487 54.3204 13.1832

Equilibrium mole fractions in the order CH^, CO, C0 2> H

2> H

2

0 are: .1327

.0000 .0000 .7245 .1429

The number of moles of gas = .9800 per mole of original gas

The volume = .804157 liters per mole of original gas

Equilibrium constants do not check out

First equilibrium constant = 7.24886e~3 vs. .197184

Second equilibrium constant = 1.20341e+10 vs. 19534.9

Third equilibrium constant = 1.61 777e- 9 vs. 5.10204e~8

now at end
61

Example F

col gas 14:16 05/03/78

Is this a constant pressure (p), or a constant volume (v) problem, and

do you want a short (s) or a long (1) printout?p,l

What are the temperature (kelvins) and pressure (atmospheres)?500,50

Temperature = 500 K = 226.85 C Pressure = 50 atm = 5.06625 MPa

What are the assumed concentrations in mole fractions of the gases in

the order: CH^, CO, CO

2

, H^, and F^O?. 12,0, .01 ,.75,. 12

.12 0 .01 .75 .12

The initial volume = .820569 liters per mole of original gas

Gas contains C, H, and 0 atoms in the proportions .13 , 2.22 , and .14

Do you want to choose the starting case?no

no

Do you want to provide your own initial estimates?yes

yes

What are they in the same order as above? .13 l.e-15 l.e-14 .71

.14

0 .13 l.e-15 l.e-14 .71 .14

0 0 0
62

Case IV

1 .13 1 . 00465e- 1

5

1.01427e-14 .71 . 14 -1. 79487 54.3204 13. 1832

2 . 13 1 . 023e- 1

5

1.071053-14 .71 . 14 -1. 79487 54.3204 13. 1832

3 .13 1 . 0628e- 1

5

1. 1989e-14 .71 .14 -1. 79487 54.3204 13. 1832

4 .13 1 . 1 2872e~ 1

5

1 .4278e-14 .71 . 14 -1. 79487 54.3204 13. 1832

5 .13 1 . 21 986e~ 1

5

1.78602e-14 .71 . 14 -1. 79487 54.3204 13. 1832

6 .13 1.32744e-15 2. 28501 e~ 14 .71 . 14 -1. 79487 54.3204 13. 1832

7 .13 1 .43493e-15 2.89084e-14 .71 . 14 -1. 79487 54.3204 13.1832

8 .13 1.52342e-15 3. 5031 3e- 1

4

.71 .14 -1. 79487 54.3204 13. 1832

9 .13 1 . 581 22e- 1

5

3. 98772e- 14 .71 . 14 -1. 79487 54.3204 13. 1832

10 . 13 1 . 60992e- 1

5

4.26874e-14 .71 . 14 -1. 79487 54.3204 13. 1832

11 .13 1 . 62027e- 1

5

4.38108e-14 .71 .14 -1. 79487 54.3204 13. 1832

12 .13 1 . 62284e- 1

5

4.41047e-14 .71 . 14 -1. 79487 54.3204 13. 1832

13 .13 1 . 62325e- 1

5

4.41 521 e- 1

4

.71 . 14 -1. 79487 54.3204 13. 1832

14 .13 1 . 62328e- 1

5

4.41 564e~14 .71 . 14 -1. 79487 54.3204 13. 1832

15 .13 1 . 62329e- 1

5

4. 41 565e- 14 .71 .14 -1. 79487 54.3204 13. 1832

Equi librium mole fractions in the order CH^, CO, co

2

, h

2

. H 2

0 are: .132653 1 . 65641 e- 1 5 4.50576e-14 .72449 .142857

The number of moles of gas = 1. per mole of original gas

The volume = .820569 liters per mole of original gas now at end
APPENDIX 1

SYMBOLS - TEXT

Helholtz energy of the system

Number of components

Number of gram atoms of element j

Number of degrees of freedom

Function actually minimized for constant volume problems

Function actually minimized for constant pressure problems

Gibbs energy of the system

Gibbs energy of the chemical elements

Jacobian of the functions set equal to zero for minimization

Inverse of J

Function expressing conservation of element j

Number of chemical species

Number of phases (discussion of phase rule only)

Pressure of the system

Number of independent chemical reactions (subsection Checking the

Estimates for Stopping Only)

Universal gas constant

Thermodynamic temperature

Volume of the system

Coefficients in the empirical equation for In kCoefficients in the empirical equation for the natural logarithm

of the vapor pressure of water

Number of gram atoms of carbon in the system

Number of chemical elements in the system

Value at the t-th iteration of the functions set equal to zero

for minimization

Number of gaseous substances

Number of gram atoms of hydrogen in the system

Equilibrium constant of formation of A B Equilibrium constant of formation of tne^i-th substance

Value of weighting factor k for t-th trial

Natural logarithm

Maximum exponent of temperature in equations for equilibrium

constants of formation and vapor pressure of water

Number of moles of solid carbon

Total number of moles of gas

Number of moles of the i-th substance

Number of gram atoms of oxygen in the system

Partial pressure of substance i Saturated vapor pressure of water

Total number of substances

Number of trial or interation

Coefficients and subscripts in chemical equation of formation
64

AG°

r A

a

g

A B

x y

with all

Vector representing the variables being estimated at

iteration number t Gibbs energy of formation of one mole of

substances in the standard states

Function of A which is minimized

Function of G which is minimized

Lagrangian multiplier

Chemical potential, partial molar Gibbs energy of the

i-th substance
65

I ( i , j ) J(i , j)

N(i , tl

)

a$

a,b,d,e

b$

cc$

c(i,D

c(i )$

cl

c2

d$

dl

f(i,D

fnc

fnw

9

g(ti

)

h

k

k k k 4*5*6

Hi)

n

Pp$

p( tl

)

pO

tt$

tO , t9

APPENDIX 2

SYMBOLS - PROGRAM

Inverse of the Jacobian,

J

Jacobian

Number of moles of the i-th substance for the

t-th trial

"Yes" or "no" to query about choosing case Constants used to differentiate functions for the cases

"Yes" or "no" to query about choosing initial estimates

Number of gram atoms of carbon in the system

"p" or "v" for constant pressure or constant volume

problem

Unweighted correction to the i-th variable

Case number i as a string variable

Case number during an iteration set

Case number found as a result of iterations

"1" or "s" for length of printout desired

Relative difference between saturated vapor pressure

of water and partial pressure of water calculated from

mole fraction found by iteration

The i-th function being set equal to zero Carbon function

Water function

Original number of moles of gas

Number of moles of gas

Number of gram atoms of hydrogen in the system

Weighting factor for applying corrections

Equilibrium constants for the three independent

reactions calculated from the logarithms of the

equilibrium constants of formation

Equilibrium constants for the three independent

reactions calculated from mole fractions found

by iteration

Natural logarithm of equilibrium constant of

formation of the i-th substance

Size of the pertinent Jacobian

Number of gram atoms of oxygen in the system

Pressure in atmospheres

Pressure as a string variable

Pressure of a constant volume system at the tl-th

trial

Initial pressure in a constant volume system

Universal gas constant

Number of sets of iterations performed

Temperature in kelvins

Temperature as a string variable

Starting and ending indexes for the iteration

loop
cm

ro

<3-

66

t-. Number of an iteration

* Square of t

Cube of t

Fourth power of t

Volume

vO Initial volume in a constant pressure system

x(i) Mole fraction of the i-th gaseous substance
67

APPENDIX 3

LISTING

TOO' colgas WSHorton 9/7/77

110' Calculates equilibrium composition for the CH4, CO, C02, H2, H20

system

120' by minimizing the Gibbs energy for a constant pressure problem or the

130' Helmolz energy for a constant volume problem.

140'

150' Initial Statements

160'

170width80

180randomize

190dimg( 150) ,1(15,15) ,J(15,15) ,N(15,150) ,p(150)

200deffnw=2*n(6,tl-l )+n(7,tl-l)-l(6)

21 Odef fnc=c-n( 1 ,tl-l)-n(2,tl-l)-n(3,tl-l)

220r=. 08205686 ' V-zero = 24. 41383e-3mt3molet-l

230c ( 1 )$=" I" 240c(2)$="II"

250c(3)$=" III"

260c(4)$=" IV"

270'

280' Input

290'

300print"Is this a constant pressure (p), or a constant volume (v) problem,"

310print"and do you want a short (s) or a long (1) printout";

320 i nputc$ ,d$

330print"What are the temperature (kelvins) and";

340 i fc$="p"orc$="P"then380

350print" volume (liters)";

360 i nputt$ ,

v

370goto400

380print" pressure (atmospheres)";

390 i nputt$ ,p$

400 i f t$( 1, 1 )="*" then430

41 0t=val (t$)

420goto440

430t=val (t$ (2)) +273. 15

440pri nt"Temperature ="t;"K =" ; t- 273 . 15;"C"

450 i fc$="p"orc$="P"then 480

460pri nt"Vol ume ="v;"l iters"

470goto530

480i fp$( 1 ,l)="*"then510

490p=val (p$)

500goto520

510p=val(p$(2))/. 10325

520pri nt"Pressure ="p;"atm =" ;p*0. 101325;"MPa"

530pri nt

540pri nt"What are the assumed concentrations in mole fractions of the gases"
68

550print"in the order: CH4, CO, C02, H2, and H20"

;

560 i nputn( l,0),n(2,0),n(3,0),n(4,0),n(5,0)

570g=g(0)=n( 1

, 0)+n(2 ,0)+n(3 ,0)+n(4,0)+n(5 ,0)

580c=n(l ,0)+n(2,0)+n(3,0)

590h=4*n(l ,0)+2*n(4,0)+2*n(5,0)

600o=n(2 ,0)+2*n(3 ,0)+n(5 ,0)

61 Oi fc<t<>0andh<t<>0ando<t<>0then650

620print"It is necessary to have compounds containing carbon, oxygen,"

630print"and hydrogen. Please try again."

640goto540

650 i f n ( 1 , 0)+n(2 ,0)+n(4 ,0)<t<>0then680

660pri nt"Carbon dioxide and water alone do not react significantly"

670goto3850' print, end

680fori=l to5

690 printn(i,0);

700 i fn( i ,0)<>0then720

710 n(i,0)=le-7

720nexti

730pri nt

740pri nt

750 i fc$="p"orc$="P"then810

760a=0

770p=p0=r*t/v

780p(0)=g(0)*p0 ' Initial estimate of pressure

790print"The initial pressure ="p(0);"atm per mole of original gas"

800goto840

810a=l

820v0=r*t/p' Initial volume

830print"The initial volume ="v0;"l iters per mole of original gas"

840pri nt

850print"Gas contains C, H, and 0 atoms in the proportions'^;" ,"h;" , and"o

860ifd$="s"ord$="S"then970

870print"Do you want to choose the starting case";

880 i nputa$

890pri nta$

900 i fa$="no"then960

91 Opri nt"What case number";

920 i nputcl

930pri nt"Case number"cl

940pri nt

950goto990

960pri nt

970c 1 =0 1 Causes bypass of line 1530

980 i fd$="s"ord$="S"thenl 180

990print"Do you want to provide your own initial estimates";

1 000 i nputb$

1 01 Opri ntb$

1 020pri nt

1 030i fb$="yes" then 1 050

1 040gotol 090
69

1050pri nt"What are they in the same order as above";

1 060 i nputn(l ,0),n(2,0),n(3,0),n(4,0),n(5,0)

1070printn(l ,0);n(2,0);n(3,0);n(4,Q);n(5,0)

1080print

1090pri nt" 0";

1 1 00g(0)=0

1 1 10fori=l to5

1120 printn(i,0);

1130 ifn(i ,0)<>0thenl 150

1140 n(i ,tl)=10i(-7-int(4*rnd))

1150 g(0)=g(0)+n(i ,0)

1 160nexti

1 170print

1 180n(6,0)=n(7,0)=n(8,0)=n(9 ,0)=0

1190'

1200' Thermodynamic Data

1210 ' 1220t2=t*t

1 230t3=t*t2

1 240t4=t2*t2

1 250£( 1 )=821 2. 35/ t- . 81 6063*1 og(t)~ 1 . 00801 e~2*t+7. 73804e~6*t2

1260£(l)=£(l)-3. 06981 e-9*t3+5. 08552e-13*t4

1 270£(2)=1 3632. 3/t+l . 8055*1 og(t)- . 0023936* t+3. 72836e~7*t2

1280£(3)=47326.7/t+. 0777872*1 og(t)-2. 80706e-4*t+3. 48747e-8*t2

1 290£(4)=0

1 300£(5)=28786. 6/t-. 69731 I*log(t)-.00145957*t+8.3769e-7*t2-1.80285e-10*t3

1 31 0£(7)=55. 5049-9430. 74/t- . 1 59684*t+3. 07499e-4*t2-3. 0549e-7*t3+l . 24325e~ 1 0*t4

1320£(6)=1(5)-1(7)

1330kl=exp(£(2)+£(5)-£(3)-£(4))

1340k2=exp(£(l)+£(5)-£(2)-3*£(4))

1 350k3=exp(2*£(2)-£(3)

)

1 360fori=l to5

1370 m(i )=log(p)-£(i

)

1 380nexti

1390'

1400' Assign max. size and constant elements of basic Jacobian

1410'

1420n=9

1430matJ=zer(n , n)

1 440J (2,7 )=J (7,2 )=J (5,7 )=J ( 7 , 5 )=1

1450J(3 , 7)=J(7 ,3)=J(4 ,6)=J(6 ,4)=J(5 ,6)=J(6 ,5)=2

1 460J ( 1 ,6)=J(6,1)=4

1470'

1480' Assign constants for first set of iterations

1490'

1 500s=0

1 51 0t0=l

1 520t9=30

1 530onc 1 goto 1 630 , 1 790 , 1 920 , 2080

1540'
70

1550' Test for possibility of liquid water"

1560'

1 5 7 0 i ft<647. 3andp>exp(£(7)) then 1600

1 580cl=2

1 590gotol 790

1 600cl=4

1610goto2080

1620‘

1630' Case I 1640'

1 650 i f s>0thenl 680

1 660 i fd$=" s"ord$="S" then 1 680

1670pri ntn(6 ,0) ;n(7,0);n(8,0)

1 680b=e=0

1 690d=l

1 700n=8

1710J(6,8)=J(8,6)=2

1 720J(7 ,8)=J(8,7)=1

1 730i fd$="s"ord$="S"thenl 770

1 740pri nt

1 750pri nt"Case I"

1 760pri nt

1 770goto2220 1 iterate

1780'

1790' Case II 1800'

1810i fs>0thenl840

1 820 i fd$="s"ord$="S" then 1840

1830printn(6,0) ;n(7,0)

1840b=d=e=0

1850n=7

1 860 i fd$="s"ord$="S" then 1 900

1870pri nt

1 880pri nt"Case II"

1890pri nt

1900goto2220' iterate

1910'

1920' Case III

1930'

1 940 i f s>0thenl 970

1 950 i fd$="s"ord$="S" then 1 970

1 9 6 0 p r i ntn(6 ,0) ;n(7,0);n(8,0);n(9,0)

1 970b=d=l

1 980e=n=9

1 990J(6 ,8)=J(8,6)=2

2000J(7 ,8)=J(8,7)=1

201 0J( 1 ,9)=J(9,1)=J(2,9)=J(9,2)=J(3,9)=J(9,3)=1

2020 i fd$=" s"ord$="S" then2060

2030pri nt

2040pri nt"Case III"
71

2 0 5 0 p r i nt

2060goto2220' iterate

2070'

2080' Case IV

2090'

2100ifs>0then2130

21 1 Oi fd$="s"ord$="S"then21 30

2120pri ntn(6 ,0) ;n(7,0);n(8,0)

21 30b=l

21 40d=0

2150e=n=8

2160J(1 ,8)=J(8,1)=J(2,8)=J(8,2)=J(3,8)=J(8,3)=1

21 70i fd$="s"ord$="S"then2220

2180pri nt

2190pri nt"Case IV"

2200pri nt

2210“

2220' Newton-Raphson Iterations

2230'

2240fortl=t0tot9

2250 f(l ,l)=log(n(l ,tl-l))-a*log(g(tl-l))+4*n(6,tl-l)+b*n(e,tl-l)+m(l

)

2260 f(2,l)=log(n(2,tl-l))-a*log(g(tl-l))+n(7,tl-l)+b*n(e,tl-l)+m(2)

2270 f(3,l)=log(n(3,tl-l))-a*log(g(tl-l))+2*n(7,tl-l)+b*n(e,tl-l)+m(3)

2280 f(4,l)=log(n(4,tl-l))-a*log(g(tl-l))+2*n(6,tl-l)+m(4)

2290 f(5,l)=log(n(5,tl-l))-a*log(g(tl-l))+2*n(6,tl-l)+n(7,tl-l)+m(5)

2300 f (6 , 1 )=4*n( 1 ,tl-l)+2*n(4,tl-l)+2*n(5,tl-l)+d*2*n(8,tl-l)-h

2310 f(7, 1 )=n( 2 , 1

1 - 1 )+2*n(3,tl-l )+n(5,tl-l )+d*n(8,tl-l )-o

2320 onclgoto2330, 2390, 2350, 2380

2330 f(8,l)=fnw

2340 goto2390

2350 f (8 , 1 )=f nw 2360 f(9,l)=-fnc

2370 goto2390

2380 f (8 , 1 )=-fnc

2390 i fc$="p“orc$="P"then2440

2400 fori=lto5

2410 J(i,i)=l/n(i,tl-l)

2420 nexti

2430 goto2520

2440 fori=lto5

2450 forj=lto5

2460 i f i=jthen2490

2470 J(i,J)=-l/g(tl-l)

2480 goto2500

2490 J(i , i )=l/n(i , tl-1 )-l/g(tl-l

)

2500 nextj

2510 nexti

2520 matI(n,n)=inv(J(n,n))

2530 matc(n,l )=I(n,n)*f(n,l)

2540 ifd$="s"ord$="S"then2560
72

2550 printtl;

2560 fori=lton

2570 ifi>5then261Q

2580 iftl>15then2610

2590 k=sin((t.l-8)*&pi/16)/2+l/2

2600 got.o2620

2610 k=l

2620 n(i ,tl)=n(i ,tl-l )-k*c(i , 1

)

2630 ifn(i ,tl )>0then2660

2640 i f i>5then2660

2650 n(i ,tl )=10t(-7-int(4*rnd))

2660 i fd$="s"ord$="S"then2680

2670 printn(i ,tl )

;

2680 nexti

2690 ifd$="s"ord$="S"then2710

2700 print

2710 g(tl)=n(l ,tl)+n(2,tl)+n(3,tl)+n(4,tl)+n(5,tl)

2720 i fc$="p"orc$="P ll then2740

2730 p=p(tl)=g(tl)*pO

2740 fori=lton

2750 ifabs(c(i ,1 )/n(i , t 1 ))>5e~5then2780

2760 nexti

2770 goto2880

2780nexttl

2790tl=tl-l

2800pri nttl ; " iterations did not converge!"

2810fori=l to5

2820 x(i )=n(i , tl )/g(tl

)

2830nexti

2840k4=x(2)*x(5)/(x(3)*x(4))' C02 + H2 = CO + H20

2850k5=x( 1 )*x(5)/(x(2)*x(4)i3*p*p) 1 CO + 3H2 - CH4 + H20

2860k6=x(2)*x(2)*p/x(3) 1 C02 + C(s) = 2C0

2870goto3270

2880s=s+l

2890'

2900' Case-consistency checking

2910“

2920fori=l to5

2930 x(i )=n(i ,tl )/g(tl

)

2940nexti

2950k4=x(2)*x(5)/(x(3)*x(4))' C02 + H2 = CO + H20

2960k5-x(l)*x(5)/(x(2)*x(4)i3*p*p)' CO + 3H2 = CH4 + H20

2970k6=x(2)*x(2)*p/x(3) 1 C02 + C(s) = 2C0

2980d]=(x(5)*p~exp( 1 (7)))/exp(l (7)) ' H20: (p.p. - v.p.)/v.p.

2990 i fcl>2then3030

3000 i f f nc>0then3050

301 0c2=2

3020goto3060

3030ifk6>k3then3050

3040goto301

0
73

3050c2=0

3060onclgoto3070, 3090, 3070, 3090

30701 f n(8 , tl )>0then31 20

3080goto3100

30901 fx(5)*p>exp(£(7))then3120

3100c2=c2+2

31 1 0goto31 30

31 20c2=c2+l

31 301 fc2=cl then3340 ' successful exit

3 1 4 0 i fs>3then3250 ' unsuccessful exit; four inconsistent sets of iterations

3150cl=c2

3160fori=8ton

3170 forj=lton

3180 J(i,J)=J(J,i)=0

3190 nextj

3200nexti

3210onclgotol630, 1790, 1920,2080

3220'

3230' Output

3240'

3250pri nt"Coul d not determine correct case with four sets of iterations"

3260pri nt

3270pri nt"The results at this stage are as follows"

3280pri nt"Case "c(cl)$

3290pri nt"The last set of estimates is, moles and multipliers:"

3300fori=l ton

3310 pri ntn( i

, tl )

;

3320nexti

3330pri nt

3340pri nt

3350pri nt"Equi 1 ibri urn mole fractions in the order CH4, CO, C02, H2, H20 are:"

3360fori=l to5

3370 i fd$="£"ord$=" L"then3400

3380 printusing".#### ",x(i);

3390 goto341

0

3400 printx(i);

341 Onexti

3420pri nt

3430pri nt

3440 i fd$=" 1 "ord$="L"then3470

3450printusing"The number of moles of gas = #.#### per mole of original gas",g(t

1 )/g

3460goto3480

3470pri nt"The number of moles of gas ="g(tl )/g;"per mole of original gas"

3480 i fc$="p"orc$="P" then3510

3490pri nt"The pressure ="g(tl )*p0/g;"atm = ". 101325*g(tl )*p0/g;

"MPa per mole of original gas"

3500goto3520

3510print"The volume ="g(tl)*v0/g;" liters per mole of original gas"

3520pri nt
74

3530oncgoto3540 , 3540,3590 , 3630

3540 i fd$=ll £"ord$=" L"then3570

3550pri ntusi ng"Sol id carbon is present in the amount of #.#### moles per mole of

original gas",fnc/g

3560goto3580

3570print"Sol id carbon is present in the amount of "fnc/g;"moles per mole of original

gas"

3580 i fCl=2then3630

3 5 9 0 i fd$="£" or d$="L"then3620

3600printusing"Liquid water is present in the amount of #.#### moles per mole of

original gas" ,n(8,tl )/g(0)

361 0goto3630

3620print"Liquid water is present in the amount of "n(8,tl )/g;"moles

per mole of original gas"

3630pri nt

3640ifabs(k4-kl)/kl>le-4then3710

3650ifabs(k5-k2)/k2>le-4then3710

3660onc 1 go to3670 , 3670 , 3690 , 3690

3670ifabs(k6-k3)/k3>le-4then3710

3680goto3750

3690if(k3-k3)/k3>le-4then3710

3700goto3750

37 1 Opr i nt"Equi 1 i bri um constants do not check out"

3720print"First equilibrium constant ="kl;"vs."k4

3730pri nt"Second equilibrium constant ="k2;"vs."k5

3740pri nt"Thi rd equilibrium constant =" k3

; "vs. 11 k6

3 7 5 0 i ft>647. 3then3840

3760onclgoto3770, 3790, 3770, 3790

3 7 7 0 i fabs(dl )>5e-5then3800

3780goto3850

3790 i fdl <5e-5then3850

3800print"Relationship of H20 partial pressure and vapor pressure is wrong"

3810pri nt"Vapor pressure of water ="exp(£(7))

3820print"Partial pressure of water ="p*x(5)

3830goto3850

3840pri nt"Above the critical temperature for water"

3850pri nt

3860' end
75

REFERENCES

[1] Glasstone, S. , "Thermodynamics for Chemists," D. van Nostrand

Company Inc. , New York, 1947.

[2] Stanczyk, N. H. , "Research Support on Materials for Coal

Conversion Process Task II," Fourth Quarterly Report, JulySeptember 1977, ERDA Contract E(49-18)-2219, Tuscaloosa

Metallurgical Research Center, Bureau of Mines.

[3] McGlashan, M. L. , Ed., "Manual of Symbols and Terminology for

Physiochemical Quantities and Units," International Union of Pure

and Applied Chemistry, Butterworths , London, 1970.

[4] Hamming, R. W. , "Introduction to Applied Numerical Analysis," McGraw

Hill Book Company, New York, 1971, p 238.

[5] Krieger, F. J. and White, W. B. , J. Chem. Phys. , ]_6, 358-60 (1948).

[6] Sokolnikoff, I. S. , Advanced Calculus, McGraw Hill Book Company,

New York, 1939, p 327-34.

[7] Scarborough, J. B. , Numerical Mathematics and Analysis, 5th ed.

,

John Hopkins Press, Baltimore, 1962, p 213-17.

[8] Gordon, S. and McBride, B. J. , Computer Program for Calculation of

Complex Chemical Equilibrium Compositions, Rocket Performance,

Incident and Reflected Shocks, and Chapman-Jouget Detonations,

NASA SP-273, NTIS No. N71-37775, 1971, p 245.

[9] Van Zeggeren, F. and Storey, S. H. , The Computation of Chemical

Equilibria, Cambridge University Press, Cambridge, 1970.

[10]

Brinkley, S. R. , J. Chem. Phys., JJ5, 107-10 (1947).
76

[11] White, W. B. , Johnson, S. M. and Dantzig, G. B. , J. Chem. Phys.

28, 751-5 (1958).

[12] Schaefer, A. 0., Samans, C. H. , Howes, M. A., Bhattacharyga , S.

Bangs, E. R. , Hill, V. L. and Chang, F. C. , Annual Report for

1975, No. FE 1784-12, Metals Properties Council, Inc., January 20, 1976.

[13] Zeleznik, F. J. and Gordon, S. , Ind. Eng. Chem. 60 27-57 (1968).

[14] Spencer, F. E. and Orning, A. A., An Algorithm for Calculating

Multiphase Chemical Equilibrium, Bureau of Mines Report RI 7786,

1973, p 42.

[15] Stull, 0. R. and Prophet, H. , JANAF Thermochemical Tables, 2nd ed.

Nat. Bur. Stand. Publ. NSRDS-NBS37, U.S. Gov't. P.0. Cat. No.

C 13.48:37.

[16] Jennings, A., Matrix Computation for Engineers and Scientists, John Wiley

and Sons, New York, 1977.

[17] Householder, A. S. , The Theory of Matrices in Numerical Analysis, Dover,

p 10, 1975.

[18] Klein, M. , Practical Treatment of Coupled Gas Equilibrium, Chapter 7

in Physical Chemistry: An Advanced Treatise in Eleven Volumes, Ed.

H. Eyring, W. Jost and D. Henderson, Academic Press, New York, Volume

I, Thermodynamics, 1971.

[19] Baron, R. E. , Porter, J. H. and Hammond, 0. H. , Jr., Chemical

Equilibria in Carbon-Hydrogen-Oxygen Systems, MIT Press, (Cambridge

USA), 1976.

'

[20] Hogben, D. , Peavy, S. T. and Varner, R. N. , OMNITAB II Users

Reference Manual, Nat. Bur. Stand. Tech. Note 552, October 1971.
77

[21] Wampler, R. H. , J. Am. Stat. Assoc., 65, 549-65 (1970).

[22] American Standard FORTRAN, X3. 9-1966, American Standards Association,

New York, March 7, 1966, p 38.

[23] Brinckley, S. R. , Jr., J. Chem. Phys. 1_4, 563-4 (1946).

[24] Schmidt, E. , Properties of Water and Steam in Si-Units, SpringerVerlag New York, Inc. , New York, 1969.

FEDERAL INFORMATION PROCESSING STANDARD SOFTWARE SUMMARY

01. Summary date

Yr. Mo. Day 7 8 0 8 0 4

04. Software date

Yr. Mo. Day

7 8 0 4 1 4

02. Summary prepared by (Name and Phone)

W. S. Horton (301) 921-2859

05. Software title BASIC Computer Program for Calculating Simple

Coal Gasification Equilibria

06. Short title COLGAS

03. Summary action New Replacement Deletion

Previous EDInternal Software ID

07. Internal Software ID N.A.

08. Software type

Automated Data Q System

fc] Computer Program

| [ Subroutine/Module

09. Processing mode

[Xj Interactive Batch

I I Combination

10. GeneralApplication area

Computer Systems

j~| Support/Utility

[jg Scientific/ Engineering

11. Submitting organization and address 12. Technical contact(s) and phone

Center for Materials Science

National Bureau of Standards W. S. Horton

Washington, D. C. 20234 (301) 921-2859

13. Narrative

Management/

Business

Process Control Other

Specific

Chemical Equilibria

COLGAS calculates the concentrations of gases and amounts of solid carbon and/or liquid

water at equilibrium for constituents found in coal gasification plants. The gases

are restricted to methane, carbon monoxide, carbon dioxide, hydrogen, and water and are considered to be ideal. Minimization of Gibbs energy is achieved by setting partial

derivatives equal to zero and solving the resulting nonlinear equations by the NewtonRaphson technique.

14. Keywords

BASIC; carbon; carbon dioxide; carbon monoxide; coal gasification; COLGAS; hydrogen;

interactive computer program; methane; minization of Gibbs energy; phase rule;

15. Computer manufr and model 16. Computer operating system 17. Programing language(s) 18. Number of source program state- , UNIVAC 1108 CSTS-Release 18 BASIC ments 3??

19. Computer memory requirements 20. Tape drives 21. Disk/ Drum units 22. Terminals Program only: 10K None None

one, time-sharing,

110 or 300 baud

23. Other operational requirements

None

24. Software availability

Available Limited In-house only

25. Documentation availability

Available Inadequate

m

In-house only

26. FOR SUBMITTING ORGANIZATION USE

Documentation available from NTIS or the author: W. S. Horton, "An Algorithm and

BASIC Computer Program for Calculating Simple Coal Gasification Equilibria", NBS-IR

78-1509

185-101 Standard Form 185

1974 July

U.S. Dept, of Commerce—NBS

(FIPS. Pub. 30)
INSTRUCTIONS

01. Summary Date. Enter date summary prepared. Use Year, Month, Day format: YYMMDD.

02. Summary Prepared By. Enter name and phone number (including area code) of individual who prepared this summary.

03. Summary Action. Mark the appropriate box for new summary, replacement summary or deletion of summary. If this software summary is a replace ment, enter under “Previous Internal Software ID" the internal software identification as reported in item 07 of the original summary, and enter the new internal software identification in item 07 of this form; complete all other items as for a new summary. If a software summary is to be deleted, enter under “Previous Internal Software ID" the internal software identification as reported in item 07 of the original summary; complete only items

01, 02, 03 and 11 on this form.

04. Software Date. Enter date software was completed or last updated. Use Year, Month, Day format: YYMMDD.

05. Software Title. Make title as descriptive as possible.

06. Short Title. (Optional) Enter commonly used abbreviation or acronym which identifies the software.

07. internal Software ID. Enter a unique identification number or code.

08. Software Type. Mark the appropriate box for an Automated Data System (set of computer programs), Computer Program, or Subroutine/Module,

whichever best describes the software.

09. Processing Mode. Mark the appropriate box for an Interactive, Batch, or Combination mode, whichever best describes the software.

10. Application Area. General: Mark the appropriate box which best describes the general area of application from among: Computer Systems Support/Utility Process Control Management/Business Bibliographic/Textual

Scientific/Engineering Other

Specific: Specify the sub-area of application; e.g.: “COBOL optimizer” if the general area is "Computer Systems Support/Utility”; “Payroll” if the general area is "Management/Business"; etc. Elaborate here if the general area is “Other.”

11. Submitting Organization and Address. Identify the organization responsible for the software as completely as possible, to the Branch or Division level, but including Agency, Department (Bureau/Administration), Service, Corporation, Commission, or Council. Fill in complete mailing address, including mail code, street address, city, state, and ZIP code.

12. Technical Contact(s) and Phone: Enter person(s) or office(s) to be contacted for technical information on subject matter and/or operational aspects

of software, include telephone area code. Provide organization name and mailing address, if different from that in item 11.

13. Narrative. Describe concisely the problem addressed and methods of solution. Include significant factors such as special operating system modifications, security concerns, relationships to other software, input and output media, virtual memory requirements, and unique hardware features.

Cite references, if appropriate.

14. Keywords. List significant words or phrases which reflect the functions, applications and features of the software. Separate entries with semicolons.

15. Computer Manufacturer and Model. Identify mainframe computer(s) on which software is operational.

16. Computer Operating System. Enter name, number, and release under which software is operating. Identify enhancements in the Narrative (item 13).

17. Programing Language(s). Identify the language(s) in which the software is written, including verata6££&, ANSI COBOL, FORTRAN V, SUBSCRIPT 11.5, SLEUTH II.

18. Number of Source Program Statements. 'Include statements in this software, separate macros, called su^-nutinea, e*c.

19. Computer Memory Requirements. Enter minimum internal memory necessary to execute software, exclusive of memory required for the operating system. Specify words, bytes, characters, etc., and number of bits per unit. Identify virtual memory requirements In the Narrative (item 13).

20. Tape Drives. Identify number needed to operate software. Specify, if critical, manufacturer, model, tracks, recording density, etc.

21. Disk/Drum Units. Identify number and size (in same units as “Memory”—item 19) needed to operate software. Specify, if critical, manufacturer, model, etc.

22. Terminals. Identify number of terminals required. Specify, if critical, type, speed, character set, screen/line size, etc.

23 Other Operational Requirements. Identify peripheral devices, support software, or related equipment not indicated above, e.g., optical character devices,

facsimile, computer-output microfilm, graphic plotters.

24. Software Availability. Mark the appropriate box which best describes the software availability from among: Available to the Public, Limited Availability

(e.g.: for government use only), and For-ln-house Use Only. If the software is "Available”, include a mail or phone contact p>oint, as well as the pnee

and form in which the software is available, if possible.

25. Documentation Availability. Mark the appropriate box which best describes the documentation availability from among: Available to the Public, In- adequate for Distribution, and For In-house Use Only. If documentation is “Available”, include a mail or phone contact point, as well as the price and

form in which the documentation is available, if possible. If documentation is presently “Inadequate”, show the expected availability date.

26. For Submitting Organization Use. Thi$ area is provided for the use of the organization submitting this summary. It may contain any information deemed useful for internal operation.

GPO : 1974 O - 554-655
NBS-114A (REV. 7-73)

U.S. DEPT. OF COMM.

BIBLIOGRAPHIC DATA

SHEET

1. PUBLICATION OR REPORT NO.

NBSIR 78-1509

2. Gov’t Accession

No.

3. Recipient’s Accession No.

4. TITLE AND SUBTITLE

An Algorithm and BASIC Computer Program for Calculating

Simple Coal Gasification Equilibria

5. Publication Date

6. Performing Organization Code

7. AUTHOR(S)

William S. Horton

8. Performing Organ. Report No.

9. PERFORMING ORGANIZATION NAME AND ADDRESS

NATIONAL BUREAU OF STANDARDS

DEPARTMENT OF COMMERCE

WASHINGTON, D.C. 20234

10. Project/Task/Work Unit No.

11. Contract/Grant No.

12. Sponsoring Organization Name and Complete Address (Street, City, State, ZIP)

Same as Item 9

13. Type of Report & Period

Covered

14. Sponsoring Agency Code

15. SUPPLEMENTARY NOTES

16. ABSTRACT (A 200-word or less factual summary of most significant information. If document includes a significant

bibliography or literature survey, mention it here.)

Calculation of the equilibrium composition for the gases CHi* , CO, CO2 , H2 , and H2 O

is treated by minimizing the Gibbs energy, G. Minimization is constrained by the

conservation of chemical elements. With the use of Lagrangian multipliers, the

minimum is found by setting the partial derivatives of G with respect to the amount

of each substance equal to zero. The resulting non-linear equations axe solved

iteratively by the Newton-Raphson method. This algorithm is implemented with an interactive computer program written in the BASIC language and named COLGAS. The

aim of this work was to provide people who test materials in coal-gasification-like

atmospheres an easy way to obtain the equilibrium composition of their gas mixtures.

A knowledge of computer programming is not required in order to use the program.

A listing of the program is given and also six sample computer calculations. The

phase rule is applied to the C-H-0 system and two ternary diagrams are shown

illustrating the condensation of solid carbon and liquid water.

17. KEY WORDS (six to twelve entries; alphabetical order; capitalize only the first letter of die first key word unless a proper name; separated by semicolons)

BASIC; carbon; carbon dioxide; carbon monoxide; coal gasification; COLGAS; hydrogen;

interactive computer program; methane; minimization of Gibbs energy; phase rule;

simultaneous chemical equilibria; water

18. AVAILABILITY [X^ Unlimited 19. SECURITY CLASS 21. NO. OF PAGES

(THIS REPORT)

Q For Official Distribution. Do Not Release to NTIS 81

UNCLASSIFIED

f 1 Order From Sup. of Doc., U.S. Government Printing Office 20. SECURITY CLASS 22. Price

Washington, D.C. 20402, SD Cat. No. 03 (THIS PAGE)

|

‘ Order From National Technical Information Service (NTIS) $6. 00

Springfield, Virginia 22151 UNCLASSIFIED

USCOMM.DC 29042- P 7 4
.