x Pkys. Chem. Sotidf. Vol. 42, pp. 297-301.19si Printed in Grtti Britain. Per^mori Pre** Ltd. 1 A REGULAR SOLUTION MODEL FOR PHASES WITH SEVERAL COMPONENTS AND SUBLATTICES,, SUITABLE FOR COMPUTER APPLICATIONS Bo Sundman and John Agren Division of Physical Metallurgy, Roya! Institute of Technology, S 100 44 Stockholm 70, Sweden (Received 24 March 1980; accepted 18 August J980) Abstract—A thermodynamic^ model, based on the regular solution approximation is presented and a formalism, suitable for phases with an arbitrary number of sublattices, is developed. A new concept, the component array, is introduced in order to simplify the analytical expressions for the integral Gibbs energy. The definition of the component array allows a straightforward procedure for the derivation of expressions describing the Gibbs energy for any kind of phase. Expressions for the partial Gibbs energy are derived. The implementation of the model on a computer is discussed. I. INTKtmUCTION The term regular solution was introduced by HildebrandfJ] to describe mixtures, whose behaviour show some experimental regularities. Later Guggenheim {2] proposed the term "strictly-regular-solutions" to cover mixtures that show an idea! entropy of mixing but have a non-zero interchange energy. The enthalpy of mixing was obtained by counting the number of nearest neighbours of different kinds. As long as the various species, atoms or molecules, are sufficiently similar in size, shape, electronegativity etc., it may be a good approximation to assume random mixing. If the species are more different, they will not mix randomly and as a consequence the solution will not have an ideal entropy of mixing. If they are sufficiently different, they may even occupy different sublattices and only those which are similar go into the same sublattice. It may then be possible to use a regular-solution model which assumes random mixing within each sublattice. Recently, Hillerf and Staffansson [3] proposed a model that describes mixing in two different sublattices. In the first version the model covered the case with 4 components that mixed pair-wise in two sublattices. An extended version, capable of treating an arbitrary number of components on both sublattices, was later presented by Harvig[4]. More recently, Hiliert and Waldenstrom [5] introduced concentration dependent interaction parameters. In the future there may be a need for further extensions, e.g. by the introduction of more sublattices and higher-order interaction parameters. The mathematical description of the model grows more complex as more components and more sublattices are added. However, very complex systems can be handled by the application of computer techniques. As such techniques are now being developed for the purpose of calculating phase equilibria, it is interesting to generalize the regular solution model to an arbitrary number of sublattices and components. This will be done in the present report. The result will be presented in such a form that it can be applied to any special case without further mathematical derivations. Furthermore, a fully automatic, computer-operated procedure wilt be described which can handle all such cases. t. the model 2.1 Definition of site fractions The mole fractions in a phase are generally defined by the following relation xt~ni/%nl = nh (1) where n; and nt are the number of moles of component / and ; respectively and n is the total number of moles in the phase. Xhe_number of components is c. It is convenient to introduce a"sTmilar composition parameter for each sublattice if the phase is composed of two or more sublattices and to include the number of moles of vacant sites, nv„ in the summation = v.+2«,'^ (2) The superscript -s denotes the sublattice under consideration, y* thus represents the site fraction of component i on sublattice s. One can then define the site fraction of vacant sites and obtain the relation v^=i-2yA (3) k It is often convenient to regard the vacant sites in a sublattice as a component and include them in the summation for the sublattice in eqn (2) but not in the summation for the whole phase in eqn (1). We also introduce a' as the number of sites on the sublattice s per mole of formula units of the phase. The site fractions can be arranged in a / x c matrix if 291 B. Sundman and J. Acren / there are / sublattices and c components y.' ^2 y2J y> y>2 superscript and subscript can be represented by a component array /, which defines one component for each sublattice. Using this notation one can write eqn (7) in the following generalized form (8) yi — (4) Each row represents a sublattice and each column a component. Since most components do not enter into all sublattices, many components y' are zero. It is always possible to calculate the mole fraction from a given Y matrix as "- —7-v- (5) 2.2 The ideal entropy of mixing \ The assumption of random mixing on each sublattice yields an expression for the ideal molar entropy of mixing in the phase. It can be written in the following -s(nw",//? = 2"*2.v*iny*. (6) Sometimes it is possible to express all the y values in terms of the x values and thus to express the ideal entropy as a function of the x parameters. However, such an expression is not always possible to formulate and when possible it is always more complicated. " 2.3 The state of reference for the Gibbs energy We need to define a frame of reference for the Gibbs energy of a solution which contains two or more sublattices. This problem was considered by Hillert and Staffansson[3]. For the simple case where there are two sublattices and two components on each one, A, B and C, D respectively, they suggested the following expression . = y^cGA.^y^GBt^y^G^t^ - +yB>-D°GB.1oOJ ^y.- • - . :.; (7) where "Ga^q, is the Gibbs energy of one mole of formula units of the pure compound Aa\Cai and the other quantities have an equivalent definition. An alternative notation of these quantities would be °GJ' where / is the component on the s sublattice and j the component on the t sublattice. More sublattices can easily be added. When eqn (7) is generalized, the principle suggested by Hillert and Staffansson makes it necessary to consider every possible compound and multiply its Gibbs energy by the fraction of that compound, i.e. the product of the corresponding site fractions, e.g. y'yf in the case of two sublattices. As another alternative the information contained in the "Gi represents the Gibbs energy of the compound defined by J which may of course depend on temperature and pressure, and Pi(Y) represents the corresponding product of site fractions from the Y matrix. The frame of reference defined by eqn (8) is composed of up to c'-.terms if there arg t sqblatticesj.and c com-. ; ponents and this wouTo*"yieIa''a ma^munTnumber of 16. terms for the case of 4 components and 2 sublattices. However, only 4 terms were included in eqn (7) because A and B only entered the first sublattice and C and D the second. Four y values were thus zero and a large number of hypothetical compounds could be omitted. 2.4 The excess Gibbs energy ; , . It is now possible to define an excess Gibbs energy, EGm, of a solution from the following expression" ■ ' . •" . ,-r\- ,•>,• .:;•. Gm = GJ«.-TS„'a<«+EGm (?) fG„ is mainly composed of the interaction energies between different components in the same sublattice. The interactions between neighbouring atoms in different sublattices are essentially described by the C^"' terms. According to the original regular solution model the contribution from the pair-wise interaction between atoms of different components can be expressed by the terms. /■--'-■:;"-;.'.: •;>■;;;.;•' - -lit (10) where Ka represents the interaction energies. By generalizing this model to phases with two or more sublattices we obtain the following expression''' v ■ ;r L .1 I I - (H) By definition K« =0. The values of the K quantities should normally depend upon what components' are present in the other sublattices and Hillert and Staffans^ son suggested the following expression for the case of two sublattices, ' : '• '-: ' ' '■k EGm=|S22Sy'y^'^ " cm where j denotes one sublattice and t the other. The superscript sst gives the; sublattices where the components i, j and Jt are located. Each such L quantity can be identified by giving a component array I, which tells that the components < and k are in the sublattices s and /, respectively, and in addition giving the information that there is another component; m the j sublattice. The A Tegular solution mode? for phases for computer applications 30) The partial Gibbs energy with respectto a component array is We first examine how D,t will operate on the sum of infraction defined as . ■ ■ , ■ (A3) We mm- derive some useful quantities that will be used to formulate eqn (A3) in terms "of derivatives of G„. First we express the number of moles of component; on sublatticc s This summation is taken over all 70 with i in sublatticc s. Of course, the following relation holds between the quantities for all the components and for all the component arrays we obtain From (A6) we find and from have used (A14) Notice that if Lk2 is composition dependent, it can only depend on the components given by KZ. If we assume that 2*, is (A5) composition independent, then eqn (Al 3) can be simplified by the observation thai the last term will run through all components / in all sublattices. As the component array contains l+Z com-ponents, this term will be non zero 1+2 times and its value wilt be Pkz(Y)Lkz each time. As before, I is the number of sublattices. Thus ..; x^-d+P^PumLKz (A1S) This simplification can be applied to any case where Lxz is composition independent. An important case is Z = 0 where lxz is replaced by 6Clt> which is composition independent by definition. We now let D,t operate on the ideal entropy of mixing, eqn (6) : Sfr1-HA""1- ~RDK(2 c'2 In>/) =-x{2 a* ^ y/ln>/+2 o'O+ln yU -X-'^W + y/In,/)} (A16) By using 2 y/ = 1 we obtain ■ (A17) By combining (A13) and (A17) the general expression for the We now apply (A10) on the expression of Gm given by eqn partial Gibbs energy can be obtained. For the simpler case (18). First we define "the partial operator" Thus ; G/o'AqC. * ;- :' (A12) composition independent interaction parameters we combine (A15) and (A17) to 4*rX<.'lny5.4^gPtf(y)L«(]+2^-(f+Z)). (A18) A regular solution model for phases for computer applications 299 following notation could be used Is}, However, this notation is slightly arbitrary since one could just as well include the / component in the component array instead of i and add the further information that there is another component / on the s sublattice. This can be avoided by simply generalizing the component array by allowing more than one component to be given for each sublattice.; \,-; ;. ''■ ■ lyl:':^■y.-..;"'', .;:•-':'"..;/.. From now on, we shall refer to a component .array, that in one sublattice contains two components but only one in each of the remaining sublattices, as a component array of the first order and it will be denoted by II. The type of component array that was introduced when defining the reference state could then be denoted /Oand be referred to as a component array of the zeroth order. Using this notation we could write «qn (12) as ; (13) where the summation is taken over all different Jl. It is easy to extend to higher order interaction by simply introducing higher order component arrays IZ. However, in order to make the definition of the interaction parameters unambiguous, we must impose the restriction that a componenl array must not contain any component more than once In each sublattice. In this way we automatically exclude parameters such as L'£ In eqn (12). ... :;;„-,v The excess Gibbs energy can then be written ~... ".; the binary A-B solution with only one sublattice is G- = JaGa + yB"GB + RT{yA In yA + yB In yB] + i'Ayi>LA:B. (16a) For the case with two sub-lattices with a1 and a7 sites respectively the model gives Ga = yA'yAi0GA.A + yA'yt?DGAB + yBWGB.A ; +ys1YB^GB.B + l^^{flWlnyA,.+VinyB,) + a2(yA*lnyAi + yB1tlnyB1)} i.+yA'yB'yAiLA.,B,A+yA,yB,ya2LA:B,B +XiVx»V^;i..X:e.: (16b) In order to be able to handle the currently available data on binary systems, interaction parameters may be given different kinds of composition dependency. The most common expression for binary systems is z>olz (14) The interaction parameters L« iriayiof course depend on temperature and pressure. As an example of this formalism one may consider a term used by Blander 16] for the interaction between all four components in a reciprocal system, i.e. two components A and B on sublattice land C and J) on sublattice 2 yAyyBycyol-A:B.ciT> 05) which is included in eqn (14) for Z=2. Notice that a component array is written with a comma between components in different sublattices and a colon between components in the same sublattice. .:■ ' As IZ must not contain a component more than once for the same sublattice, the value of Z is limited, A binary system.with one sublattice and no vacancies will thus have just one interaction parameter. Such systems are by far the most commonly investigated [7] and it is a general experience that a single parameter usually is not sufficient to represent the experimental data. A composition dependency is often introduced. However, if more than one sublattice is used to describe the system, the number of parameters may be large enough even if the Liz parameters are composition independent. As an example, consider an A-B system first having only one sublattice and then two with A and B mixing randomly on each one but not between them. The Gibbs energy for (17) where i', and »" are the two components which occupy sublattice s according to the component 11. The number b is called the degree of the parameter *Gj i. Any such expression can be incorporated in the model. We can now write the Gibbs energy per mole PUti'Gjo+RTXa-J, y'lny/ + 2 SPizOO-I* z>0 iz . 3. the partial gibbs energy The partial Gibbs energy with respect to a component i is defined in any textbook as (19) For a phase with several sublattices one can define the following-partial energy with respect to a component array of zeroth order, 70 ' XdniJr.p.tfi, (20) where »/0 is the number of moles of formula units of 70. This quantity is, in contrast to the partial Gibbs energy for the individual components, always possible to formulate explicitly. In the appendix it is shown that the well-known expression for the calculation of partial quantities can be generalized as follows • ■ t Gio=Gm + £*E it. V / 9G„ rKdyl (21) lliiiSiiliISi / 300 B. Sundman and J. Agren where i, is the component in sublattice s of the component array 10. By applying eqn (18) for Gm, the partial derivatives are obtained as follows, dGm To y; , + 2 2/6^(1* + *%). (22) z>o jz . . % \ 3y\ I The quantity S'jn is a special form of the Kronecker delta. Shi is unity if / is one of the components in sublattice s, of the component array JZ and zero otherwise. It is shown in the appendix that in the case of composition independent interaction parameters we obtain by inserting eqn (22) into eqn (21), where / is the number of sublattices. '■: ~ / • 4. implementation of the model on computer . There is a considerable need for thermodynamic evaluations in many types of calculations relating to problems of chemical equilibrium or kinetics. In order to simplify such calculations on computer it would be helpful to have a library of subroutines for the evaluation of the various thermodynamic quantities involved. In order to keep the number of such subroutines at a minimum, it is necessary to make such subroutines; as general as possible. It was thus felt that the present thermodynamic model should be programmed in a general way, such that the person working on an application program could always use the same subroutine and would only have to. specify the phase by giving the number of components and sites on each sublattice. Such a program has now been constructed in the FORTRAN language for a Nord-, 10 minicomputer. It can be modified for usage on other computers. • V-. i-^c.t'^-: — .. :-.. Since there may be a need for calculating various quantities from the integral molar Gibbs energy Gm, it was necessary to structure the data in some special way. The method chosen was inspired by. the monograph, "■The Art of Computer Programming" by: Knuth[8]. Some programming languages like PL/I, PASCAL and SIMULA are designed to handle structural data.' However, it turned out to be quite straight-forward to implement dynamic storage allocation and list processing for the present purpose in FORTRAN. The data was structured in such a way that both Gm and partial derivatives with respect to fractions can be calculated by the same subroutine. This is achieved by using the same techniques as in programs for symbolic manipulations of formulas, but in the present case the procedure is taken one step further and a numerical value rather than a formula is the final result. Another benefit of this tech- nique is that only non-zero parameter values need to be stored which reduces the space required for storage. As an extra feature the possibility of evaluating the molar volume and its derivatives was included in the program, as well as a composition dependent ferromagnetic contribution to the Gibbs energy. ■: The program has an interactive command monitor written as a separate subroutine. With this monitor the user can, from his terminal, specify his system and give all the necessary data and also inspect it and obtain listings.: There is a special Users. Guide for this monitor[9]. '.•:. ■ -•;>■•"-.'■ v The program has been tested in several applications. At present there are two different application programs available for phase diagram calculations, one for optimizing parameters in the thermodynamic mode! to fit experimental tie-lines or activities and one for simulating diffusion controlled transformations in solid or solid-liquid systems. Some results from the application programs have been reported [10,11]. v The experience has shown that the program has some advantages in addition to saving time for the person writing the' application program. In particular, it yields safer, application programs since -the highly structural data storage eliminates many errors. In view of the fact that the program derives all quantities analytically from the integral quantities, it is not possible for the user to make logical errors in his handling of thermodynamics. Among the drawbacks of the program it should be mentioned that k yields execution times which are unnecessary long for very simple cases. Furthermore, it uses a rather large part of the primary memory of the Nord-10 minicomputer. It may thus be necessary to segment an application program or to use overlay techniques. On the other hand, such drawbacks will grow less and less important due to the rapid progress of computer handware. Acknowledgement—The authors wish to thank Prof. Mats Hillert for many critical and helpful discussions.: .': references "'".'"■'' 1. Hildebrand J. H., J/Am. Chem.Soc.Sl, 66 (1929). ! 2. Guggenheim E. A. Pnc. Roy. Soc. A148,304 (1935). 3. Hitlert M. and Staffansson L.-I., Acta Chem. Scand. 24,3618 . . (1970). :;■„-.^ ':;o . - V' ,:"'- . ,■ 4. HarvigX Acta Chem, Scand. 2S, 3199 (1971). ... 5. Waldenstrom M., An Experimental and Thermodynamic •;.;>• Study of Phase Equilibria in Complex Alloyed Steels. Dok- torsavhandling. KTH(1976). :' h- - 6. Blander M., / Chem. Phys. 39,2610 (1963). : ', ■ 7. Ansara L, /. Metals Rev. 24,20 (1979). - .8. Knulh D. E.,The Art of Computer Programming, 2nd Edn, • p. 288 ff. Addision-Wtsley (1976). 9. Sundman B.( ft/. Rep. D5,1978, Rev. (1979). - . : 10. Agren S., Met. Trans A, 10A, 1847 (1979). : It, Agren J., TR/TA-AMC-0154 (1979). ./: . - - . " " APPENDIX Calculation of partial quantities Gibbs energy for n) moler of formula unils is : ■ . G = ruGm ,- (Al) where G„ is a function of temperature, pressure and the fractions of ail components given by the Y matrix Gm = Gm(T,P,Y). (A2) The par defined as We now -formulate express lh> This sumn course, the the compo and for all From (A6) and from (. Equaiior where £„• i: eqn (A3Vf We now (18). First* Thus ( 5 - T 1 lit