Computer calculations of multicomponent phase diagrams

Scripta METALLURGICA

Vol. 4, pp. 685-692, 1970 Printed in the United States

Pergamon Press,

Inc.

C o m p u t e r Calculations of Multicomponent Phase D i a g r a m s by Henri Gaye and C. H. P. Lupis

Department of Metallurgy and Materials Science Carnegie-Mellon University Pittsburgh, Pennsylvania, U.S.A. 15213 (Received

July

1,

1970)

The r e t r i e v a l of m u c h t h e r m o d y n a m i c information is l o s t without an adequate c o m p u t e r i z e d mapping of p h a s e d i a g r a m s .

M o r e o v e r , as the s t a t e of a t h e r m o d y n a m i c s y s t e m may be r e p r e s e n t e d

by a point in a m u l t i d i m e n s i o n a l phase d i a g r a m , this mapping allows the r e p r e s e n t a t i o n and s i m u l a tion of v a r i o u s t r a n s f o r m a t i o n s .

The c o m p u t e r can thus s e r v e as an invaluable guide in the design

of p r o c e s s e s , provided, h o w e v e r , that the c a l c u l a t i o n s can be made with sufficient efficiency to b e a r a reasonable cost.

The method we developed is outlined below.

It is valid for any n u m b e r of

c o m p o n e n t s , but h a s been t e s t e d , so f a r , only for b i n a r y and t e r n a r y s y s t e m s . The c r i t e r i o n of e q u i l i b r i u m for the c o e x i s t e n c e of two p h a s e s ~ and ~ in an m - c o m p o n e n t s y s t e m at c o n s t a n t t e m p e r a t u r e and p r e s s u r e may be e x p r e s s e d in two a l t e r n a t e ways: 1./

by the equality of the c h e m i c a l potential of e a c h component in the two p h a s e s : .(;)

= .(f)

(i= 1,2...m)

(1)

o r , in t e r m s of activity coefficients: Xl c') X(fl)

T1

Ti

=

-

( i = 1,2 . . . . m)

(2)

RT

i

o w h e r e Di r e p r e s e n t s the m o l a r f r e e e n e r g y of the p u r e c o m p o n e n t i and X i the mole f r a c t i o n of i. 2./

by m i n i m i z i n g , at c o n s t a n t t e m p e r a t u r e and p r e s s u r e , the Gibbs f r e e energy of the

s y s t e m , i . e . the sum of the f r e e e n e r g i e s of the two phases: G s = nzv G (~) + nflG (fl)

(3)

The m o l a r f r e e e n e r g y of e a c h phase can be e x p r e s s e d as the sum of the Raoultian s t a n d a r d f r e e energy, the ideal f r e e e n e r g y of mixing and the e x c e s s f r e e energy: G = Z X i bL°

+ RT r Xi ~nXi + G E .

685

(4)

686

Multicomponent

Phase Diagrams

Vol.

4,

No. 9

At a given temperature and p r e s s u r e , the activity coefficients and excess free energies are functions of the concentrations, and usually, can be satisfactorily represented by polynomials [1 ] . The coefficients of these polynomials are related by the Gibbs-Duhem equation [2 ] and the number of n e c e s s a r y p a r a m e t e r s (i.e. the degree of the polynomial) depends on the nature of the system and the accuracy desired [1 ].

These p a r a m e t e r s may be determined from isothermal activity data, from

points on the experimental phase diagram, or by combination of both; they can also be estimated through statistical models of solutions. The difference of the standard chemical potentials is a function of the temperature and it may be expressed as: A /~i o(~-a)

= AHI~-'a)/1 - Tt)

+

Aa i IT - T: ) - T4m ~T]

Abi 2

(T_Ti)t 2

l

-

2

T

-

(5)

1 where T t is the temperature of the transformation ( ~ -* a) for the pure component i, AH (~ " O0 is the corresponding changeo in enthalpy per mole of i and Aai, Ab I and/$c i are coefficients describing the difference AC(W.~ a ) in the specific heats of i in the two phases: p,l AC ( ~ a ) . p, 1

= A a i + Ab.T !

+ AClT-2

(6)

Although the theoretical concepts a t e relatively simple, the actual solution of either method is a difficult problem because of the nature of the transcendental equations involved. Except in the very s i m p l e s t cases, analytical solutions are impossible, and the careful selection of efficient numerical methods becomes a prerequisite.

The establishment of a stable or metastable phase

diagram also requires that the relative stabilities of various possible equilibria be compared. For convenience, the calculations are presented for the equilibrium between two phases.

.Firstmethod ibinary systems, A solution to equations (2) was investigatedin the case of binary systems only, the method becoming cumbersome when the number of components m is greater than 2. The equilibrium

:11

~(a)n

R

A2

(LI n

n

RT

R

conditions, at a given temperature, constitute a system of two equations with two unknowmx~a'and 1.n l - X 2 (a)

+

N~ (L(nl)a

1 - X~ fl)

n=0 \

x(a) 2

~

/L(2) a n

n=0

RT

6n X~ ~)

+

RT

K .

~ - nffi0

X~P!

X RT (Ta)

(7b)

Vol.

4, No.

9

Multicomponent

Phase Diagrams

687

where the L's and K's r e p r e s e n t enthalpy and entropy interaction p a r a m e t e r s [1]. Hiskes and T i l l e r [3 ] have used two s i m i l a r equations to deduce from the phase diagram the dependence of the chemical potentials on composition and temperature.

In equations (7), this consists

in fixing the values of --X(:) and X (fl) at several temperatures to the values read from the phase diagram, 2 and in calculating then the L and K unknowns. Since the equations are linear in the L's and K's, the solution is easily found. In the p r e s e n t study however, the L's and K's are assumed known and, at a given temperature, the variables are .x 2

and X~B);- the equations are then transcendental and, in

general, difficult to solve. Kaufman and Bernstein [4 ] solved them by the Newton-Raphson iteration teclmique in the simple cases of ideal solutions (all the K and L coefficients are zero) and of regular solutions (all the K coefficients are zero and only one L coefficient is different from zero).

However, in many systems

of metallurgical i n t e r e s t these assumptions are unrealistic and more powerful procedures than the ones they used are also n e c e s s a r y . of values

Rudman [5] used a trial and e r r o r method to determine the couple

of X(2a)- and X~~)- that best fits the equations (7). However, it appears that this approach

is rather time consuming; its most obvious weakness is thatD when computing the equilibrium compositions at some temperatureD no use is made of previous calculations performed at different neighboring t e m p e r a t u r e s .

Stepwise methods such as that described below are much more powerful.

The equilibrium concentrations at the temperature T being known, the equilibrium c o n c e n t r a tions at the temperature (T + AT) are computed by the relation:

X2 (T( + VAT) )

(v) d X ~ u ) \ A T + (d2X2V)~ = X 2 (T) + ~--~TT ~T

where

d~)

t

dT

dX~) dT

~

.2 v(¢0

O

A2

~AT)2

(v=a,~)

(8)

d2 X ~ ) are calculated by differentiation of the system of

dT 2

dT 2

equations (7), and the temperature i n c r e m e n t is fixed at each step so that: d2X~)

I

(AT) 2 <

1

(v) dX 2

l

dT 2

[

2

100

AT

(v = ~, ~)

(9)

dT

The validity of this expansion is then tested at temperature intervals of 1 degree by an iteration at constant temperature: (10)

where 8 ( X ') and i

8tx(~h 2 .j

are calculated by differentiation of equations (7). This provides an

accuracy of + 10 -6 on the equilibrium mole fractions.

688

Multicomponent

Phase Diagrams

Vol.

4, No. 9

Three subroutines were written to describe a eutectic, a peritectic, and a system with complete miscibility. They can be called together to describe a complex phase diagram, as illustrated by fig. 1. for the aluminum-zinc system. The number of parameters to describe any phase is arbitrary and the efficiency of the program is practically insensitive to this number. In the present example it was found that four parameters for the liquid phase, three for a and two for ~ were sufficient to represent adequately both the phase diagram and the activity data.

l~.'Ot~ [ -- . . . .

-1

~--

T

........

I

!~ ~ . ~ . _

T °K 6

f

0

oo AI

Liquid

0

~

...................................

400

200

I

u +

i

I

I

I

02

04

06

08

XZ.

Zn

FIG. I. Phase diagram of the aluminum-zinc system. The dashed line is reported by Hultgren et al [9 ] and the continuous line is calculated and plotted by the computer.

Second method: multicomponent systems. Kaufmsn and Bernsteln [4 ] generalized to isothermal sections of regular ternary solutions the method they used for binaries; it is essentially identical to that used by Hurle and Pike [6], also for regular ternary solutions. The technique outlined above for binaries, although very efficient, was found to become somewhat cumbersome for s y s t e m with more than two components. Iusteadw a more satisfactory method based on the minimization of the free energy was adopted. The free energy of a system of one mole of composition X° ( i = l , . . . m ) decomposing into L

a phase ~ and a phase ~, can be expressed, in terms of the variables na~ n ~ x_ i(~) and X(i~l(for i=213,,

o. m)

as: Gs

=

...

o,

Vol.

4,

No. 9

Multicomponent

Phase Diagrams

689

The 2 m variables are bound by m independent relations expressing the conservation of the different species: na + n~ ffi 1

(12)

and:

na x(a) + n~ x(i~) i

= Xoi

(i=2,3,...m)

(13)

G contains logarithmic and polynomial t e r m s and the constraints (13) a r e nonlinear; time-consuming S

search techniques would therefore be required to obtain a direct solution. Instead, a stepwise procedure was again developed. The free energy change when going from a known equilibrium state (Xla)" , Xi~)) of-

overall

o to a neighboring state of overall composition (Xi + dX l ) can be approximated by a composition Xi, second order T a y l o r ' s s e r i e s .

Noting that dn~ = -dnct ,

+ dna

l iffi2 L

+

t

m~

Ina

-

~

~ x(a) i

dX ~) i

b 2 G ( a ) . (a)

dxi(a) dX!a)j + n~

52 G (~) i

dXlt (14)

(-)

At equilibrium, the variables d n ~ dX i

and dXl~) (for 1=2,3, " . . m ) minimize AG s and obey

the constraints obtained by differentiation of equation (13):

na

i

This may be recognized as a quadratic programming problem.

/ However 0 it can be solved by

a simpler technique than the usual one (based on the simplex algorithm) because of the form of the .__(,v) constraints which do not include inequalities. Replacing in equation (14) the ( m - l ) variables ctXi in t e r m s of dn and dX!"~) (for i=2, 3 , . . . m), the solution is obtained by writing that the derivative of AG with r e s p e c t to the l a t t e r m variables are zero. S

The calculation is therefore reduced to the

*

O

solution of a system of m linear equations . The i n c r e m e n t in overall composition, dXl , is fixed at

When the coefficients in these equations are very d i s s i m i l a r , as in the ease where the system is dilute with r e s p e c t to one of its components, conventional techniques are very inadequate and, consequently, Cholesky's method [7 ] has been adopted.

690

Multicomponent

Phase

Diagrams

4, No. 9

Vol.

each step so that the new equilibrium compositions are in a p r e d e t e r m i n e d neighborhood of the p r e c e d ing ones: dX ~) <_ R "~t

(l=l,2,...m;l~=Ct,~)

(16)

where R is a p a r a m e t e r adjusted at each step to ensure the validity of equations (14) and (15). Near a o c r i t i c a l point~ the value of dX i(Is)"ts generally much l a r g e r than that of d X ; , so that a test of completion of the calculations is provided by imposing a minimum value to the length of the vector dX i

e

In o r d e r to obtain the complete isothermal section of the phase boundaries, at each step the O

vector dXi is chosen normal to the tie-line.

The starting points for the procedure a r e selected as a

couple of equilibrium compositions in one of the binaries (calculated by the f i r s t method), or, if that is not possible (e.g. in the c a s e of an "island" miscibility gap), they a r e determined in the neighborhood of a c r i t i c a l point, the p r o p e r t i e s of which are exploited analytinally [8 ] . An example of the r e s u l t s is given in fig. 2 for the c a s e of the lead - zinc - bismuth system. Twelve p a r a m e t e r s were used to d e s c r i b e the composition dependence of the excess f r e e energy: nine

U E

U q

U B

U,e

~lZa

FIG. 2. Isoactivity c u r v e s and isothermal section of the l e a d - z i n c - b i s m u t h system at 520UC. The experimental points on the boundaries of the miscibility gap were r e p o r t e d by Hansen e t al [11] (triangles) s Seith et al [12 ] (black c i r c l e s ) and Valenti et al [10 ] (white c i r c l e s ) . The dotted line is ~dyS~tr e p o r t e ~ by Valenti et al [10] and the continuous line is calculated me computer.

Vol.

4, No. 9

Multicomponent

Phase Diagrams

691

were obtained from the data r e p o r t e d by Hultgren et al ~9 ] on the binaries, and three f r o m the m e a s u r e m e n t s of the activity of zinc in the ternary by Valenti, Oleari and Fiorani E10~. The isoactivity curves calculated on that basis agree very well with those calculated by Valenti, et al El0 1. Fig. 2 shows that the computed phase boundaries a r e consistent with the data reported by Hansen and Anderko ~11~, Valenti et al ~10~ and Seith, Johnen and Wagner ~12~. The input of the computer p r o g r a m s , for binary systems as well as multicomponent ones, consists of the small set of p a r a m e t e r s describing each phase.

The output is given either as a table

of the equilibrium compositions, or, m o r e directly, as actual plots, of the complete phase diagram in the case of a binary and of isothermal sections in the case of a ternary. The software is written in FORTRAN V and used on Carnegie-Mellon University's UNIVAC 1108 computer.

The plot is provided by the CALCOMP PLOTTER.

The typical execution time for the

calculation of a complete binary phase diagram by the f i r s t method v a r i e s f r o m ten to twenty seconds, depending on the complexity of the system and the extent of the t e m p e r a t u r e range studied.

The

p r o g r a m for isothermal sections of t e r n a r i e s has an execution time of the s a m e o r d e r of magnitude. It was equal to five seconds for the section in fig. 2. A detailed description of the p r o g r a m s will appear in subsequent a r t i c l e s . The authors wish to acknowledge the financial support of the Mellon Foundation and the National Science Foundation (Grant GK - 11103).

References 1. C . H . P . Lupis, Acta Met., 16, 1365-1375 (1968). 2. C . H . P . Lupis and J . F . Elltott, Acta M e t . , 14~ 529-538 (1966). 3. R. Hlskes and W.A. T i l l e r , Materials SCi. and Eng.,2~ 320-330 (1968) and 4, 163-184 (1969). 4. L. Kaufman and H. Bernstein, "Computer Calculations of Phase Dlagrams"~ Academic P r e s s , New York s 1970. 5. P.S. Rudman, J. Mater. Scl. E n g . , to be published. 6. D . T . J . Hurle and E . R . Pike, J. Mater. Sci. E n g . , 1, 399-402, (1966). 7. L. Fox, U. S. Nat'1. Bur. Stand., Appl. Math. Series, 39 t 1-54, (1954). 8. H. Gaye and C. H. P. L u l l s , unpublished r e s u l t s . 9. R. Hultgren, L . R . O r r , P.D. Anderson and K.K. Kelley, "Selected Values of Thermodynamic P r o p e r t i e s of Metals and Alloys", John Wiley & Sons, Inc., New York, 1963. 10. V. Valenti, L. Oleari and M. F i o r a n i , Gazz. chim. Ital., 86, 930-941, (1956). 11. M. Hansen and K. Anderko, "Constitution of Binary Alloys", McGraw-Hill, New York, 1958. 12. W. Seith, H. Johnen and J. Wagner, Z. Metallkunde, 46, 773-779, (1955).