C.2: Equilibrium composition of mixture phases

This example will demonstrate what information can be obtained from ChemApp regarding the equilibrium composition of mixture phases. Along the way, the important concept of the stoichiometry matrix (see also Chapter E.3 for details) will be applied. The last calculation will demonstrate the relationship between partial and integral extensive property values for mixture phases, using a non-ideal liquid as an example.

Definition of goals

The aim is to show how ChemApp is used to retrieve information about the equilibrium state of a mixture phase, the gas phase of the system C-O-Si (data-file cosi.dat) being used as an example.

The information to be retrieved includes

The system Pb-Sn (data-file pbsn.dat) is then used to

Note that the items marked with an asterisk (*) indicate that the values retrieved for the properties involved are model-dependent. This is potentially always the case when getting information about phase constituents, as different mixture phase models might lead to different phase constituents. For example, when the Redlich-Kister formalism is applied to model a melt, a different set of phase constituents might be used to describe its thermochemical properties than when the quasichemical model is applied to the same melt. Thus also the information and variable values associated with the phase constituents will be different. Even values for variables like the amount of the phase might differ, if the amount is retrieved in mol and the number and/or type of phase constituents differ between two models.

Tools and concepts used

When using mixture phases in equilibrium calculations it is important to familiarise oneself with how the equilibrium state of these phases can be characterised. It is useful to understand which of the results that can be obtained from ChemApp are based on phase-internal information only, and which need phase-external information, relating to the equilibrium states of other phases and the whole system.

It will be shown that the stoichiometry matrix provides important information about the phase constituents, and how it is used to gain some of the results listed above. Familiarising oneself with the concept of the stoichiometry matrix is crucial to anybody who wants to understand the relationship between a mixture phase and its phase constituents. This is especially true for those who plan on assessing mixture phase data themselves, and thus editing thermochemical data-files.

ChemApp also allows for the retrieval of extensive properties of both phases and phase constituents. The partial and integral values for a phase and its constituents are of course related, and it is useful to see how ChemApp can be used to demonstrate this relationship.

Implementation


C Equilibrium composition of mixture phases
PROGRAM CAF101
IMPLICIT NONE
CHARACTER NAME*24
DOUBLE PRECISION STOI(30), VALS, VALUE, AGAS, ASYS, VALUES(50), * FRACSI, SITOT, WMASS, XSN, HINT1, HINT2, HPB, HSN
INTEGER NOERR, NUMCON, ICO, ISIO2, NPCON, I, NSCOM, ISI, * IPB, ISN, ILIQ
C Initialise ChemApp CALL TQINI(NOERR)
C Open the thermochemical data-file cosi.dat (system C-O-Si) C for reading CALL TQOPNA('cosi.dat', 10, NOERR)
C Read data-file CALL TQRFIL(NOERR)
C Close data-file CALL TQCLOS(10, NOERR)
C For the present example, the input composition consists of 2 mol of C CO/GAS/ and 3 mol of quartz CALL TQINPC('CO ', 1, ICO, NOERR) CALL TQSETC('IA ', 1, ICO ,2.D0, NUMCON, NOERR)
CALL TQINP('SiO2(quartz)', ISIO2, NOERR) CALL TQSETC('IA ', ISIO2, 1, 3.D0, NUMCON, NOERR)
C The temperature at which the equilibrium is to be calculated is set to C 1500 K CALL TQSETC('T ', 0, 0, 1500.D0, NUMCON, NOERR)
C Calculate the equilibrium and print a ChemSage result table CALL TQCEL(' ', 0, 0, VALS, NOERR)

Output:
 
 T = 1500.00 K
 P = 1.00000E+00 bar 
 V = 2.49275E+02 dm3

STREAM CONSTITUENTS AMOUNT/mol CO/GAS/ 2.0000E+00 SiO2(quartz) 3.0000E+00
EQUIL AMOUNT MOLE FRACTION FUGACITY PHASE: GAS mol bar CO 1.9974E+00 9.9935E-01 9.9935E-01 CO2 1.2879E-03 6.4438E-04 6.4438E-04 SiO 1.9134E-06 9.5730E-07 9.5730E-07 SiO2 3.9280E-12 1.9652E-12 1.9652E-12 Si 7.1155E-14 3.5600E-14 3.5600E-14 O 2.6014E-14 1.3015E-14 1.3015E-14 C 3.4466E-17 1.7244E-17 1.7244E-17 O2 2.0479E-17 1.0246E-17 1.0246E-17 C3 8.9044E-18 4.4551E-18 4.4551E-18 C2 2.6909E-19 1.3463E-19 1.3463E-19 Si2C 5.6221E-20 2.8129E-20 2.8129E-20 SiC 9.2988E-21 4.6524E-21 4.6524E-21 Si2 3.6182E-22 1.8103E-22 1.8103E-22 Si3 2.2100E-28 1.1057E-28 1.1057E-28 O3 1.8757E-34 9.3844E-35 9.3844E-35 TOTAL: 1.9987E+00 1.0000E+00 1.0000E+00 mol ACTIVITY SiO2(tridymite) 3.0000E+00 1.0000E+00 C 1.2860E-03 1.0000E+00 SiO2(cristobalite) 0.0000E+00 9.9878E-01 SiO2(quartz) 0.0000E+00 9.4548E-01 SiO2(liquid) 0.0000E+00 8.2521E-01 SiC 0.0000E+00 5.2057E-04 Si 0.0000E+00 3.5247E-06 ******************************************************************** Cp_EQUIL H_EQUIL S_EQUIL G_EQUIL V_EQUIL J.K-1 J J.K-1 J dm3 ******************************************************************** 2.88120E+02 -2.62627E+06 9.37547E+02 -4.03259E+06 2.49275E+02
Mole fraction of system components: GAS C 4.9984E-01 O 5.0016E-01 Si 4.7850E-07

C First, get the equilibrium amount of the gas phase, which is the only C mixture phase in the system. Note that since we did not change any C units, the default amount unit (mol) is used. CALL TQGETR('A ', 1, 0, AGAS, NOERR)
WRITE(*,FMT='(A,G12.5)') 'Equilibrium amount of gas phase ' // * 'in mol: ', AGAS

Output:
Equilibrium amount of gas phase in mol:   1.9987

C To determine the mole fraction of the gas phase, relative to the whole C system, we first need to get the amount contained in the system CALL TQGETR('A ', 0, 0, ASYS, NOERR)
WRITE(*,FMT='(A,G12.5)') 'Mole fraction of the gas phase: ', * AGAS/ASYS

Output:
Mole fraction of the gas phase:  0.39974

C The equilibrium amounts of all phase constituents of the gas phase can C be obtained by a call to TQGETR which passes an array together with C the proper INDEX/INDEXP combination CALL TQGETR('A ', 1, -1, VALUES, NOERR)
C Print the equilibrium amount of each phase constituent CALL TQNOPC(1, NPCON, NOERR) WRITE(*,*) 'Equilibrium amounts of gas phase constituents' DO I=1, NPCON CALL TQGNPC(1, I, NAME, NOERR) WRITE(*,FMT='(I4,A,G12.5,A)') I, ' : ' // NAME, * VALUES(I), ' mol' ENDDO

Output:
 Equilibrium amounts of gas phase constituents
   1 : C                        0.34466E-16 mol
   2 : C2                       0.26909E-18 mol
   3 : C3                       0.89044E-17 mol
   4 : CO                        1.9974     mol
   5 : CO2                      0.12879E-02 mol
   6 : O                        0.26014E-13 mol
   7 : O2                       0.20479E-16 mol
   8 : O3                       0.18757E-33 mol
   9 : Si                       0.71155E-13 mol
  10 : Si2                      0.36182E-21 mol
  11 : Si2C                     0.56221E-19 mol
  12 : Si3                      0.22100E-27 mol
  13 : SiC                      0.92988E-20 mol
  14 : SiO                      0.19134E-05 mol
  15 : SiO2                     0.39280E-11 mol

C The mole fractions of the phase constituents can be calculated using C the values just retrieved
C Print the mole fractions of each phase constituent WRITE(*,*) 'Mole fractions of gas phase constituents' DO I=1, NPCON CALL TQGNPC(1, I, NAME, NOERR) WRITE(*,FMT='(I4,A,G12.5)') I, ' : ' // NAME, * VALUES(I)/AGAS ENDDO

Output:
 Mole fractions of gas phase constituents
   1 : C                        0.17244E-16
   2 : C2                       0.13463E-18
   3 : C3                       0.44551E-17
   4 : CO                       0.99935    
   5 : CO2                      0.64438E-03
   6 : O                        0.13015E-13
   7 : O2                       0.10246E-16
   8 : O3                       0.93844E-34
   9 : Si                       0.35600E-13
  10 : Si2                      0.18103E-21
  11 : Si2C                     0.28129E-19
  12 : Si3                      0.11057E-27
  13 : SiC                      0.46524E-20
  14 : SiO                      0.95730E-06
  15 : SiO2                     0.19652E-11

C Information relating the system components (in the present case C, O, C and Si) to the phases can also be obtained.
C Get the mole fraction of the system components in the gas phase CALL TQNOSC(NSCOM, NOERR) WRITE(*,*) 'Mole fractions of system components in the gas phase' DO I=1, NSCOM CALL TQGNSC(I, NAME, NOERR) CALL TQGETR('XP ', 1, I, VALUE, NOERR) WRITE(*,FMT='(I4,A,G12.5)') I, ' : ' // NAME, * VALUE ENDDO

Output:
 Mole fractions of system components in the gas phase
   1 : C                        0.49984    
   2 : O                        0.50016    
   3 : Si                       0.47850E-06

C To get the fraction of a system component in the gas phase, we have to C make use of the information contained in the stoichiometry matrix.
CALL TQINSC('Si ', ISI, NOERR) FRACSI = 0.D0
C Get the equilibrium amount of each gas phase constituent, then use the C stoichiometry matrix returned by TQSTPC to determine how much Si is C contained in the constituent, and sum it all up DO I = 1, NPCON CALL TQSTPC(1, I, STOI, WMASS, NOERR) CALL TQGETR('A ', 1, I, VALUE, NOERR) FRACSI = FRACSI + STOI(ISI)*VALUE ENDDO
C Get the total amount of Si in the system CALL TQGETR('a ',0,ISI,SITOT,NOERR)
WRITE(*,FMT='(A,G12.5)') 'Mole fraction of system component ' // * 'Si in the gas phase: ', FRACSI/SITOT

Output:
Mole fraction of system component Si in the gas phase:  0.63779E-06


C To demonstrate how ChemApp can be used to check the relationship C between partial and integral values of extensive properties, a C different data-file for a system that contains non-ideal mixture C phases is used.
C Open a new thermochemical data-file (pbsn.dat, system Pb-Sn) C for reading CALL TQOPNA('pbsn.dat', 10, NOERR)
C Read data-file CALL TQRFIL(NOERR)
C Close data-file CALL TQCLOS(10, NOERR)
C Find out the system component index numbers for Pb and Sn, as well as C the phase index number for the liquid phase: CALL TQINSC('Pb ', IPB, NOERR) CALL TQINSC('Sn ', ISN, NOERR) CALL TQINP('LIQUID ', ILIQ, NOERR)
C Using the index numbers just determined, we can enter our incoming C amounts, in this case an alloy consisting of 0.7 mol Sn and 0.3 mol Pb: XSN = 0.7D0 CALL TQSETC('IA ', 0, ISN, XSN, NUMCON, NOERR) CALL TQSETC('IA ', 0, IPB,1.D0 - XSN, NUMCON, NOERR)
C Set the temperature to 600 K. Under these conditions, the liquid phase C is the only stable phase in this system (see Figure 7). CALL TQSETC('T ', 0, 0, 600.D0, NUMCON, NOERR)
C Calculate the equilibrium CALL TQCE(' ', 0, 0, VALS, NOERR)
C Get the integral enthalpy of the liquid phase CALL TQGETR('H ', ILIQ, 0, HINT1, NOERR)
C Get the partial molar enthalpies of Sn/LIQUID/ and Pb/LIQUID/ CALL TQGETR('HM ', ILIQ, ISN, HSN, NOERR) CALL TQGETR('HM ', ILIQ, IPB, HPB, NOERR)
C Calculate the integral enthalpy of the liquid phase from the two C partial values HINT2 = XSN * HSN + (1.D0 - XSN) * HPB
C Print both values WRITE(*,FMT='(A,F12.3,A)') 'HINT1 = ', HINT1,' J' WRITE(*,FMT='(A,F12.3,A)') 'HINT2 = ', HINT2,' J'

Output:
HINT1 =    16066.466 J
HINT2 =    16066.466 J
C As expected, the two values are the same

END

Suggestions for enhancements

  1. Before adding to the code, answer the following questions for yourself:

    1. In the above example, the fraction of system component Si in the gas phase has been calculated twice. Since the values are different, what Si fraction exactly has been calculated in each case?

    2. Which of the results obtained using the calls to TQGETR can also be read from the ChemSage result table printed by TQCEL?

    3. Looking at the ChemSage result table, why does the column labelled 'fugacity' contain the same values as the column labelled 'mole fraction'?

    4. When retrieving the integral enthalpy for the Pb-Sn melt, the 'H' option of TQGETR is used. Will it make a difference for the example used if the molar enthalpy (option 'HM') is used instead?

    5. What happens when a different total amount in mol is entered into the equilibrium calculation instead, will the options 'H' and 'HM', as used in the example, still result in both integral enthalpies to be equal?

  2. Extend the part of the program which gets the fraction of the system component Si in the gas phase (relative to the total amount of Si in the whole system) so that a list of fractions for all system components is printed.

  3. Calculate the integral Gibbs energy of the Pb-Sn melt without using the options 'G' or 'GM' for TQGETR. Verify it in a similar fashion as it is done for the enthalpy in the above example.


ChemApp Programmer's Manual, Edition 3.6© GTT-Technologies, 2003