A quadratic finite difference approach to phase equilibria calculations

CALJ’HALI Vol. 14, No. 3, pp. 283-288, 7990 Printed in the USA.

A QUADRATIC

FINITE

DIFFERENCE

0364-5916t90 $3.00 + .OO (c) 1990 Pergamon Press plc

APPROACH

TO PEASE

EQUILIBRIA

CALCULATIONS

R. 0. Williams 906 W. Outer Drive Oak Ridge, Tenn. 31030

ABSTRACT Over the yeara we have routinely used an iterative method based on a quadratic representation derived from finite differences for solving nonlinear problems in many variables. As the method addresses one variable at a time no differentiation is required which simplifies the coding and reduces the chance of errors and makes the use of different defining equations particularly simple. Besides using the method for phase equilibria problems we.hrrve also used it for the treatment of problems in coherency, the cluster variation method and the Bragg Williams modeling with many sublattices. The success of the method depends to some degree on how the variables are chosen but the use of more variables than required to define a problem causes no problems whatsoever and can expedite the solution in some cases.

INTRODUCTION Following Gibbs it is well known that the equilibrium in a multiphase system of N components is determined by the minimization of the Gibbs free energy and that this corresponds to a plane that is tangent to the free energy surfaces of the different phases. The tangent plane and the free energy surfaces are in fact N-l dimensional figures and thus are two dimensional only for a ternary system. This construction also corresponds to the equality of slopes of the free energy surfaces at the tangent points and the slope of the tangent plane and also corresponds to the activity of each component being identical in all the phases. Each of these three equivalent formulations have been used to calculate equilibria. It is possible to calculate the equilibrium in closed form only for very simple cases so iterative solutions are almost universally required. A graphical construction can be used for a binary system but even in this case the precision may not be adequate. The solution becomes complex as the number of components and phases increase. A representative collection of the work published in the metallurgical field is contained in references 1-12. In most of the references the mathematics and computer program are rather complex because the solution is based on estimating new values for all the variables during each iteration. The approach is generally based on the Newton Raphson method although in some cases the simplex method of Nelder and Mead(l3) is used. The present approach is based on adjusting the variables individually in succession. The papers by Counsel1 et al(2), Pelton et al(10) and Talley and Pelton(l1) relate most ____________________-----Received 22 September 1989.

283

R.O. WILLIAMS

directly to the present method. In particular, Counsel1 et alta‘) expressed the problem in terms of the total free energy of the system, a procedure that we also follow. They choose to use the simplex method(l3) to obtain the solution but it is clear from their Fig. 2(b) that for the particular problem they are illustrating the two variables can be treated alternately and the solution could be reached very quickly. The way in which Pelton et al{101 set up the problem is different from by Counsel1 et al(a) but it appears to represent the same basic formulation. They adjust the variables in succession just as we do. In the Talley and Pelton(llf paper they transfer individual components in succession from one phase to another based on the activity. When our formulation is set up in terms of the individual components it is approximately equivalent to this last approach since any adjustment always results in a decrease in the free energy of the system for both methods. THE FOBMULATZON OF THE SOLUTION We express the total free energy of a fixed amount of material distributed between two or more phases and assume that it can be approximated by an expression of the form G = A + BX + CXP

11)

where X is one of possibly many variables of the system. In particular, X may be the amount of one component in one phase of a two phase mixture or, alternately, the fraction of one of the phases. The value of X corresponding to a minimum is calculated by setting the derivative with respect to X to zero to give X=

- B/x

(31

provided that C is positive in which case the value of X may now be set to this new value. If C is negative we are in a poorly behaved region where no minimum exists. The solution is continued by similarly examining the additional variables in order and then the entire process is repeated until no further decrease in the total free energy can be obtained. The way in which this approach is implemented is as follows. We calculate the total free energy for the values of X - dX, X and X + dX and implement the above concept in the form GA=A-

BdX + CdXp

(3a)

GB =A

(3bf

GC = A + BdX + CdXz

(3c)

The value of the increment, dX, corresponding to the minimum is given as dX' = .5(GA - GC)dX/(GA + GC - 2GB)

(4)

provided that 2GB < GA + GC. Otherwise this is a naximhm so the sign of dX' must be changed but no estimate of the required magnitude to move into a well behaved region is possible. There are always physical limits on the permissible magnitude of dX’. In particular the compositions must never become negative. Zero compositions are permissible but the ideal entropy contribution for such components must be set to zero since the standard expression will now give an error. The first point about this approach is that it is so simple than it can be implemented readily from memory or else Eq. 4 can be rederived when needed. A

PHASE EQUILIBRIA CALCULATIONS

285

Taylor expansion is certainly a more precise representation but this requires more mathematical operations, computer coding and more involved calculations. The particular variables that are chosen will depend on the problem but for phase equilibria in systems of three or more components the fraction of a given component in a pair of phases examined two at a time seems to work well enough. But when only two components are present this does not work since we now are using only a single variable as the two components are exactly the choice is to treat the volume fraction related. Thus, for a binary system of each phase as the two variables. As Counsel1 et al(a) show these two variables are nearly independent so the solution converges vary rapidly. An important advantage of this approach is that one can use extra variables without causing problems. In particular, a three component system of three phases requires six variables whereas we might represent this as a nine variable problem in our approach. The inclusion of extra variables can lead to simpler coding and possibly a more efficient solution. If the increment is too small this can lead to round off errors. This may be solved by using a larger increment. One could also calculate the term GA + GC -2GB in Eq. 4 in a way to improve the accuracy or one can go to double precision calculations. Our experience has been that it is always been possible to choose the increment size so as to eliminate this problem. When the volume fraction of one phase becomes small or when the amount of one component becomes small in a phase a variable size increment may be required. There is the obvious point that one does not make an adjustment so large as to lead to negative compositions but further, one uses some care not to make adjustments appreciable larger than the increments used in the calculations as one can get into regions where the expression is no longer of adequate precision. We routinely use adjustments no larger in magnitude than twice the increment. In the appendix we give a BASIC program that implements this approach for the calculation of an asymmetric miscibility gap. The free energy of the solution is defined by the DEF FNG statement that may be changed for gaps for different free energy functions. Statements DEF FNA and DEF FNB give the activity of the A and B components within the present formulation where X is the fraction of A in each case. Other minor changes, particularly the temperature, may be required. Additional phases can be included by additional statements with the required changes made to call them. The bulk composition, X0, should correspond approximately to a 50/5Q mixture at the start. Once the solution is obtained for each temperature the activities of both components are calculated for both phases to verify the correctness of the solution. For dilute solutions the activities can change very rapidly with composition such that the activities may differ by a small amount. This program does not blow up above the miscibility gap but the compositions shift until the fraction of one of the phases becomes negative. The program detects this condition and such results are excluded from the final results. It is possible to change the size of the adjustment from that given by Eq. 4 by changing the magnitude of FAC . In the present program a value of 0.95 was used but we believe that in some problems appreciably higher values will speed the search without causing problems. If the value is too high the compositions will oscillate and a solution may not be reached. For the present parameters the temperature of 1070 k is very near the critical temperature and the behavior of the program is a bit erratic but the correct solution is reached. Execution times on a PC should be about a minute. DISCUSSION As already noted our formulation normally uses the total free energy of the

286

R.O. WILLIAMS

representation and they further adjust one variable at a time. We believe that our approach is more efficient in arriving at the solution. There is some similarity in our approach and that due to Talley and Pelton(l0) but ours is simpler in that only free energy calculations are required and the quadratic representation should reduce the number of iterations required. We certainly expect that others have used the method described here but it is not common as we have found no reference. Certainly it has not been commonly used to calculate phase equilibria. The brevity of the included program supports our claim of the utility of this approach. In nonlinear problems it is always possible that more than one minimum exists in which case it is possible to arrive at a solution that is not the absolute minimum. However, we are not aware that this has occurred in our work with the kinds of problems discussed here. Another well recognized problem is that one or more phases may not have been included so one may be treating a constrained or metastable situation. This may be intentionally. We have used this approach to solve a number of related problems that include the effects of coherency on phase equilibria(l4), the cluster variation method following Kukuchi and Murray(l5) and the Bragg Williams modeling using 8 sublattices for the FCC structure following Buth and Inden(l6). A section through the miscibility gap in the U-Nb-Zr system calculated by this method is given in reference 17. The Bragg-Williams problem is interesting. It is a seven variable problem in which one could interact one sublattice with the other seven in turn or, using our approach, we could interact each sublattice with all the others giving a total of 56 steps per iteration. Suppose that l/3 of the sublattices are lean in B, l/3 have the correct B content and l/3 are B rich. With the later formulation about l/3 of the time little transfer of B would be possible, 4/9 of the time a modest transfer could take place and 2/9 of the time conditions exist for a very favorable transfer. We are not sure that for this particular problem the use of the extended variable set expedites the solution but in phase equilibria problems it certainly happens that certain combinations of phases will be unfavorable and the use of an extended variable set will expedite the solution. This occurs when one is trying to transfer appreciable amounts of a component through a phase that contains very little of this component. In dealing with complex problems there may well be a desire to understand the structure in some detail and this approach may be found deficient in this respect but the simplicity and expediency make up for this in some measure. One reason why the method works so well is that the behavior of the total free energy around a minimum is well represented by this quadratic representation. Had we used only a linear representation things would work less well since we then we would have no direct method of estimating the size of the corrections. CONCLUSIONS We have presented a general approach to solving complex problems based on the repeated use of a quadratic expression based on finite differences. Some care in how the variables are chosen will improve the efficiency of the solution but extra variables can be included that reduce problems arising from poor choices. The efficiency of the method is so great that there is little incentive to reduce the variables to a more or less optimum set. ACKNOWLEDGEMENT We appreciate the constructive comments of Drs. T. B. Lindemer and E. C. Beahm of the Oak Ridge National Laboratory where we were also employed before retirement.

287

PHASE EQUILIBRIA CALCULATIONS

REFERENCES 1.

I. Ansara. J. N. Barbier, P. Y. Chevalier and M. Ii. Rand, CHEMICAL METALLURGY, A Tribute to Carl Wagner, ed. N. A. Gokcen, p. 310, AIME, Warrendale (1981).

2.

J. F. Counsell, E. B. Lees and P. J. Spencer, Metal Science J., 5, 210 (1971).

3.

P. Dorner, L. J. Gauckler, K. Krieg, Ii. L. Lukas, G. Petzow and J. Weiss, CALPHAD, 3, 241 (1979).

4.

G. Eriksson, Acta Chem. Scan. 25, 2651 (1971).

5.

B. Gale and J. M. Davis, Metal Sci. J. 5, 25 (1971).

6.

Ii. Gaye and C. H. P. Lupis, Metall. Trans., 6A. 1049 (1975).

7.

H. Gaye and C. H. P. Lupis, Metall. Trans., 6A, 1057 (1975).

8.

H. L. Lukas, J. Weiss and E. T. Henig, CALPHAD, 6, 229 (1982).

9.

L. Kaufman and H. Bernstein, COMPUTER CALCULATION Academic Press, New York, (1970).

OF

PHASE DIAGRAMS,

10.

A. D. Pelton, C. W. Bale and M. Rigaud, 2. Metallk., 68, 135 (1977).

11.

P. K. Talley and A. D. Pelton, CALPHAD, 12, 171 (1988).

12.

J. S. van Duijneveldt, F. 13, 133 (1989).

13.

J. A. Nelder and R. Mead, Computer J., 7, 308 (1965).

14.

R. 0. Williams, CALPHAD, 8, 1 (1984).

15.

R. Kikuchi and J. L. Murray, CALPHAD, 9, 311 (1985).

16.

J. Buth and G. Inden, Acta Metall., 30, 213 (1982).

17.

R. 0. Williams and J. Brynestad, CALPHAD, 6, 25 (1982).

S.

A. Baas and H. A. J. Oonk, CALPHAD,

R.O. WILLIAMS

APPENDIX A BASIC

Program

to Calculate

Asymmetric

Miscibility

Gaps

10 'a-24-89, Robin Williams 20 DEF ~G(X,Y)=X*Y*(16OOO!+4OOO!*~Y-X~~+8.314*T*~X*LOG~X~+Y*LOG~Y~~ 30 DEF FNA(X)=X*EXP((1-X~-2*(16000+4000*~1-4*X~~/~8.314*T~~ 40 DEF FNB(X~=~1-X~*EXP~X-2*~16000+4000*~3-4*X~~/~8.314*T~~ 50 DIM T(20),X1(20),X2(20),1T(20),G(20) 60 XO=.5:X1=.01:X2=.99:DX=.OO4:Yl=O:Y2=l!:FAC=.95 70 FOR I=1 TO 12:T=470!+1*50 80 PRINT:PRINT" Xl x2 G T,k = ",T:PRINT 90 G2=FNG(X2,1-X2):FOR J=l TO 50 100 Zl=X1-DX:V1=~X2-XO)/~X2-Zl~:GA=Vl*FNG~Z1,l-Zl~+~l-Vl~*G2 110 Zl=Xl:Vl=(X2-XO)/(X2-Zl):GB=Vl*FNG(Zl,l-Zl~+~l-Vl~*G2 120 Zl=X1+DX:V1=~X2-XO)/~X2-Zl~:GC=Vl*FNG~Zl,l-Zl~~~l-Vl~*G2 130 GOSUB 290:X1=X1+DX*BB*FAC:Gl=FNG(Xl,l-Xl) 140 Z2=X2-DX:Vl=(Z2-XO)/~Z2-Xl):GA=Vl*Gl+~1-Vl~*FNG~Z2,l-Z2~ 150 Z2=X2:V1=(Z2-XO)/(Z2-Xl):GB~Vl*Gl+(l-Vl~*FNG~Z2,l-Z2~ 160 Z2=X2+DX:Vl=~Z2-XO)/~Z2-Xl~:GC=Vl*Gl+~1-Vl~*FNG~Z2,l-Z2~ 170 GOSUB 290:X2=X2+DX*BB*FAC:G2=FNG(X2,l-X2) 180 V1=(X2-XO)/(X2-X1):GO=Vl*Gl+~l-Vl~*G2 190 PRINT USING "######.####";Xl,X2,GO 200 IF ABS(Yl-Xl)+ABS(YZ-X2) < .00004 GOT0 220 210 Yl=Xl:Y2=X2:NEXT J:II=I-l:GOTO 260 220 II=I:IF Xl > X0 OR X2 < X0 THEN GOT0 260 230 X1~I~=X1:X2~I~=X2:T~I~=T:IT(I~=J:XO=~X1+X2~/2:G~I~=GO:PRINT 240 PRINT"Activity of X ";:PRINT USING "######.###";FNA(Xl),FNA(X2) 250 PRINT"Activity of X-l ";:PRINT USING "######.###";FNB(Xl),FNB(X2):NEXT I Iterations, Fat = ",FAC:PRINT 260 PRINT:PRINT" T,k Xl x2 270 FOR I=1 TO 1I:PRINT USING "###t##.";T(II; 280 PRINT USING "#####.###";Xl(I),X2(I);:PRINT w ";IT(I):NEXT 1:END 290 CC=(GA+GC_2*GB)*2:IF CC < 0 THEN CC=-CC 300 BB=(GA-GC)/CC:IF ABS(BB) > 2 THEN BB=SGN(BB)*2! 310 RETURN