Using the friendly Module

This chapter provides instructions on how to use ChemApp for Python’s friendly module, which provides a friendlier Python interface to the ChemApp library than the tq* interface in the basic module.

The purpose of the friendly module is to make is quicker and easier to set up and do thermochemical calculations, and process calculation results. It also aims to make the code you write less error prone, and easier to read and understand.

ChemApp for Python’s friendly vs basic Module

The friendly module differs from the basic module in the following ways, which are discussed further below:

  • Descriptive naming

  • Groups of functions

  • More specific functions

  • More powerful functions

  • Sophisticated data types

Descriptive Naming

This aspect is explained in the Introduction chapter of this manual, and it is especially relevant to the friendly module. We use clear names for function groups, and clearly defined abbreviations to create function, parameter and variable names. This makes the code that you write with ChemApp for Python much easier to read and understand than code written with ChemApp, or with the ChemApp for Python’s basic module.

Functions Organised into Groups

The ChemApp API consists of 75 functions that can be used for a variety of tasks related to doing thermochemical calculations. It may be a daunting prospect to get started with this interface.

The friendly module organises functions into groups. Each group contains functions related to a specific aspect of thermochemical calculations. This makes them easier to find, learn, and remember. The groups are listed and explained in the Function Groups section.

More Specific Functions

Many of ChemApp’s tq* functions can be used to accomplish a variety of different tasks. They are very compact, but this can make them more difficult to decipher and use. As an example, the tqsetc function can be used for 21 different things.

The friendly module takes a different approach. Rather than having a small number of functions that can each do many things, we create more functions that each does one specific thing. This makes it easier to learn and understand what a function does and how to use it. It also makes the code you write easier to read.

More Powerful Functions

The friendly module also provides functions that do more work with less coding. For example, to open, read, and close a thermochemical data file with ChemApp requires at least three lines of code. With friendly you can do this with one.

There are also functions in the friendly module that help you to print lists of system components, phases, phases constituents, etc. with a single line of code. With ChemApp this would requires more effort.

Using Sophisticated Data Types

The friendly module uses the sophisticated data types provided by The core Module. This means that you can access objects like system components, phases, and phase constituents, knowing that these objects contain all their relevant information. Here are two examples illustrating the benefits. They both print out the details of a specific system component. The first uses only primitive data types, and the second uses a core module class.

from chemapp.friendly import ThermochemicalSystem as cats

# load a thermochemical data file
cats.load("fec.dat")

# get the configuration details of a system component
index = 0
name = cats.get_name_sc(index)
mm = cats.get_mm_sc(index)

# print the system component configuration
print("================")
print("SYSTEM COMPONENT")
print("================")
print(f"Index:      {index}")
print(f"Name:       {name}")
print(f"Molar mass: {mm}")
print("================")

When executing this code, the output is as follows:

================
SYSTEM COMPONENT
================
Index:      0
Name:       Fe
Molar mass: 55.847
================

Here is the alternative, employing the SystemComponentConfig class. When calling the ThermochemicalSystem.get_config_sc function, a SystemComponentConfig object is returned.

from chemapp.friendly import ThermochemicalSystem as cats

# load a thermochemical data file
cats.load("fec.dat")

# get the configuration of a system component
sc = cats.get_config_sc("Fe")

# print the system component configuration
print(sc)

The output of the second example is:

=========================================================================
System Component Configuration
-------------------------------------------------------------------------
Class:      SystemComponentConfig
Index:      0
Name:       Fe
Molar mass: 55.847
=========================================================================

With the core module classes you get more done with less effort. The same principle demonstrated here for system components applies to phases and phase constituents, and to the thermochemical system as a whole as well.

Function Groups

ChemApp for Python uses Python static classes to divide all its functions into easy-to-understand groups, so that they are easier to find and use. The following groups are available:

  • Info

  • Units

  • ThermochemicalSystem

  • EquilibriumCalculation

  • PhaseMapCalculation

  • StreamCalculation

Each group is explained briefly below. For more information about the functions and how to use them, please consult the friendly module’s reference documentation, which forms part of the Package Reference chapter.

Info Group

The purpose of this group is to provide you with functions to get information about ChemApp itself. This includes things like a copyright message, version number, the user that ChemApp is licensed to, and the sizes of arrays that ChemApp uses.

Units Group

The Units group is used to manage the units of measure used in thermochemical calculations. It contains functions to get the current units, change units, and print out the set of units currently in use as nicely formatted strings.

Upon loading ChemApp for Python, the default units will be the same as the default units of ChemApp:

  • Temperature - K

  • Pressure - bar

  • Amount - mol

  • Energy - J

  • Volume - dm3

ThermochemicalSystem Group

This group is used to load and interact with thermochemical systems. A thermochemical system is described by a thermochemical data file, which can be in text format with a .dat extension, or in encrypted format with a .cst extension. The encrypted files are also called transparent data files.

The functions in this group are used to get information about system components, phases, phase constituents, etc. It is also used to extract configuration objects from the thermochemical system if you want to work with the more sophisticated class data types provided by the core module.

EquilibriumCalculation Group

This group is used to do absolute equilibrium calculations that provide extensive property results as \text{C}_p, \text{H}, \text{S}, \text{G}, rather than as \Delta\text{C}_p, \Delta\text{H}, \Delta\text{S}, and \Delta\text{G}.

StreamCalculation Group

This group is used to do relative equilibrium calculations that can provide extensive property results as both absolute (e.g. \text{H}) and relative values (\Delta\text{H}). This is achieved by defining streams that enter an equilibrium calculation. The final extensive properties of the system are then calculated relative to those of the incoming streams.

PhaseMapCalculation Group

This group allows one-dimensional phase map calculations, which is another type of absolute equilibrium calculation. You can use these calculations to search for phase boundaries as a function of temperature, pressure, composition, chemical potential, etc. You can even use it to construct binary and pseudo-binary phase diagrams.

Base Classes

When you look a bit deeper into the friendly module, you will find more static classes that represent function groups. They are EquilibriumCalculationBase and AbsoluteEquilibriumCalculationBase. The word Base at the end of the names indicates that they are base classes that we use internally in ChemApp for Python. You should not use them in your code, not even import them. They are only used internally in ChemApp for Python.

Worked Examples

This section provides the same set of worked examples that is contained in the ChemApp documentation provided by GTT-Technologies, but using the friendly interface instead of the basic one. When you compare the basic and friendly worked examples, the benefits of using the friendly interface should be clear.

Equilibrium Composition of Mixture Phases

"""
Equilibrium composition of mixture phases.
"""

# Import the ChemApp for Python modules that we will use in our calculations.
from chemapp.core import (TemperatureUnit,
                            PressureUnit, VolumeUnit, AmountUnit,
                            EnergyUnit)
from chemapp.friendly import Units as cau
from chemapp.friendly import ThermochemicalSystem as cats
from chemapp.friendly import EquilibriumCalculation as caec
from chemapp.friendly import StreamCalculation as casc

# Determine the data file path.
from pathlib import Path
import chemapp
data_dir = Path(chemapp.__file__).parents[0] / 'data'

# Initialise ChemApp for Python.
# Open, read and close the thermochemical data-file cosi.dat (system
# C-O-Si).
cats.load(str(data_dir / 'cosi.dat'))

# Change default units to desired units.
cau.set(T = TemperatureUnit.K,
        P = PressureUnit.bar,
        V = VolumeUnit.dm3,
        A = AmountUnit.mol,
        E = EnergyUnit.J)

# For the present example, the input composition consists of 2 mol of
# CO/GAS/ and 3 mol of quartz.
ico = cats.get_index_pc(0, 'CO')  # index of gas phase is always 0.
caec.set_IA_pc(0, ico, 2.0)

isio2 = cats.get_index_ph('SiO2(quartz)')
caec.set_IA_pc(isio2, 0, 3.0)

# The temperature at which the equilibrium is to be calculated is set to
# 1500 K.
caec.set_eq_T(1500.0)

# Calculate the equilibrium and print a ChemSage result table.
casc.calculate_eq(print_results=True)

# First, get the equilibrium amount of the gas phase, which is the only
# mixture phase in the system. Note that since we changed the units, the
# amount unit (mol) is used.
agas = caec.get_eq_A_ph(0)

print(f'\nEquilibrium amount of gas phase in mol: {agas:.4f}')

# To determine the mole fraction of the gas phase, relative to the whole
# system, we first need to get the amount contained in the system.
asys = caec.get_eq_A()

print(f'\nMole fraction of the gas phase: {agas / asys:.5f}')

# The equilibrium amounts of all phase constituents of the gas phase can
# be obtained by using the get_eq_A_pcs_in_ph function and the phase
# index.
values = caec.get_eq_A_pcs_in_ph(0)

# Print the equilibrium amount of each phase constituent.
npcon = cats.get_count_pcs_in_ph(0)
print('\nEquilibrium amounts of gas phase constituents:')

names = cats.get_name_pcs_in_ph(0)
for i in range(npcon):
    print(f'{i:>5} : {names[i]:<20} {values[i]:.5E} mol')

# The mole fractions of the phase constituents can be calculated using
# the values just retrieved.

# Print the mole fractions of each phase constituent.
print(f'\nMole fractions of gas phase constituents:')
for i in range(npcon):
    print(f'{i:>5} : {names[i]:<20}  {values[i] / agas:.5E}')

# Information relating to the system components (in the present case
# C, O, and Si) to the phases can also be obtained.

# Get the mole fraction of the system components in the gas phase.
nscom = cats.get_count_scs()
print(f'\nMole fractions of system components in the gas phase:')
names = cats.get_name_scs()
values = caec.get_eq_X_scs_in_ph(0)
for i in range(nscom):
    print(f'{i:>5} : {names[i]:<20} {values[i]:.5}')

# To get the fraction of a system component that reports to the gas
# phase, we have to make use of the information contained in the
# stoichiometry matrix.
isi = cats.get_index_sc('Si')
fracsi = 0.0

# Get the equilibrium amount of each gas phase constituent, then use the
# stoichiometry matrix to determine how much Si is contained in the
# constituent, and sum it all up.
for i in range(npcon):
    wmass = cats.get_mm_pc(0, i)
    stoi = cats.get_stoich_pc(0, i)
    value = caec.get_eq_A_pc(0, i)
    fracsi += stoi[isi] * value

# Get the total amount of Si in the system.
sitot = caec.get_eq_A_sc(isi)

print(f'\nMole fraction of system component Si that reports to the gas '
      f'phase: {fracsi / sitot:.5E}')

# To demonstrate how ChemApp for Python can be used to check the relationship
# between partial and integral values of extensive properties, a
# different data-file for a system that contains non-ideal mixture
# phases is used.

# Open a new thermochemical data-file (pbsn.dat, system Pb-Sn) for
# reading.
cats.load(str(data_dir / 'pbsn.dat'))

# Find out the system component index numbers for Pb and Sn, as well as
# the phase index number for the liquid phase:
ipb = cats.get_index_sc('Pb')
isn = cats.get_index_sc('Sn')
iliq = cats.get_index_ph('LIQUID')

# Using the index numbers just determined, we can enter our incoming
# amounts, in this case an alloy consisting of 0.7 mol Sn and 0.3 mol
# Pb:
xsn = 0.7
caec.set_IA_sc(isn, xsn)
caec.set_IA_sc(ipb, 1.0 - xsn)

# Set the temperature to 600 K. Under these conditions, the liquid phase
# is the only stable phase in this system (see the Pb-Sn phase diagram
# provided in worked example 3).
caec.set_eq_T(600)

# Calculate the equilibrium.
caec.calculate_eq()

# Get the integral enthalpy of the liquid phase.
hint1 = caec.get_eq_H_ph(iliq)

# Get the partial molar enthalpies of Sn/LIQUID/ and Pb/LIQUID/.
hsn = caec.get_eq_HM_pc(iliq, isn)
hpb = caec.get_eq_HM_pc(iliq, ipb)

# Calculate the integral enthalpy of the liquid phase from the two
# partial values.
hint2 = xsn * hsn + (1.0 - xsn) * hpb

# Print both values.
print(f'\nhint1 = {hint1:.3f} J')
print(f'hint2 = {hint2:.3f} J')

# As expected, the two values are the same.

Equilibrium Composition of Mixture Phases

"""
Equilibrium composition of mixture phases.
"""

# Import the ChemApp for Python modules that we will use in our calculations.
from chemapp.core import (TemperatureUnit,
                            PressureUnit, VolumeUnit, AmountUnit,
                            EnergyUnit)
from chemapp.friendly import Units as cau
from chemapp.friendly import ThermochemicalSystem as cats
from chemapp.friendly import EquilibriumCalculation as caec
from chemapp.friendly import StreamCalculation as casc

# Determine the data file path.
from pathlib import Path
import chemapp
data_dir = Path(chemapp.__file__).parents[0] / 'data'

# Initialise ChemApp for Python.
# Open, read and close the thermochemical data-file cosi.dat (system
# C-O-Si).
cats.load(str(data_dir / 'cosi.dat'))

# Change default units to desired units.
cau.set(T = TemperatureUnit.K,
        P = PressureUnit.bar,
        V = VolumeUnit.dm3,
        A = AmountUnit.mol,
        E = EnergyUnit.J)

# For the present example, the input composition consists of 2 mol of
# CO/GAS/ and 3 mol of quartz.
ico = cats.get_index_pc(0, 'CO')  # index of gas phase is always 0.
caec.set_IA_pc(0, ico, 2.0)

isio2 = cats.get_index_ph('SiO2(quartz)')
caec.set_IA_pc(isio2, 0, 3.0)

# The temperature at which the equilibrium is to be calculated is set to
# 1500 K.
caec.set_eq_T(1500.0)

# Calculate the equilibrium and print a ChemSage result table.
casc.calculate_eq(print_results=True)

# First, get the equilibrium amount of the gas phase, which is the only
# mixture phase in the system. Note that since we changed the units, the
# amount unit (mol) is used.
agas = caec.get_eq_A_ph(0)

print(f'\nEquilibrium amount of gas phase in mol: {agas:.4f}')

# To determine the mole fraction of the gas phase, relative to the whole
# system, we first need to get the amount contained in the system.
asys = caec.get_eq_A()

print(f'\nMole fraction of the gas phase: {agas / asys:.5f}')

# The equilibrium amounts of all phase constituents of the gas phase can
# be obtained by using the get_eq_A_pcs_in_ph function and the phase
# index.
values = caec.get_eq_A_pcs_in_ph(0)

# Print the equilibrium amount of each phase constituent.
npcon = cats.get_count_pcs_in_ph(0)
print('\nEquilibrium amounts of gas phase constituents:')

names = cats.get_name_pcs_in_ph(0)
for i in range(npcon):
    print(f'{i:>5} : {names[i]:<20} {values[i]:.5E} mol')

# The mole fractions of the phase constituents can be calculated using
# the values just retrieved.

# Print the mole fractions of each phase constituent.
print(f'\nMole fractions of gas phase constituents:')
for i in range(npcon):
    print(f'{i:>5} : {names[i]:<20}  {values[i] / agas:.5E}')

# Information relating to the system components (in the present case
# C, O, and Si) to the phases can also be obtained.

# Get the mole fraction of the system components in the gas phase.
nscom = cats.get_count_scs()
print(f'\nMole fractions of system components in the gas phase:')
names = cats.get_name_scs()
values = caec.get_eq_X_scs_in_ph(0)
for i in range(nscom):
    print(f'{i:>5} : {names[i]:<20} {values[i]:.5}')

# To get the fraction of a system component that reports to the gas
# phase, we have to make use of the information contained in the
# stoichiometry matrix.
isi = cats.get_index_sc('Si')
fracsi = 0.0

# Get the equilibrium amount of each gas phase constituent, then use the
# stoichiometry matrix to determine how much Si is contained in the
# constituent, and sum it all up.
for i in range(npcon):
    wmass = cats.get_mm_pc(0, i)
    stoi = cats.get_stoich_pc(0, i)
    value = caec.get_eq_A_pc(0, i)
    fracsi += stoi[isi] * value

# Get the total amount of Si in the system.
sitot = caec.get_eq_A_sc(isi)

print(f'\nMole fraction of system component Si that reports to the gas '
      f'phase: {fracsi / sitot:.5E}')

# To demonstrate how ChemApp for Python can be used to check the relationship
# between partial and integral values of extensive properties, a
# different data-file for a system that contains non-ideal mixture
# phases is used.

# Open a new thermochemical data-file (pbsn.dat, system Pb-Sn) for
# reading.
cats.load(str(data_dir / 'pbsn.dat'))

# Find out the system component index numbers for Pb and Sn, as well as
# the phase index number for the liquid phase:
ipb = cats.get_index_sc('Pb')
isn = cats.get_index_sc('Sn')
iliq = cats.get_index_ph('LIQUID')

# Using the index numbers just determined, we can enter our incoming
# amounts, in this case an alloy consisting of 0.7 mol Sn and 0.3 mol
# Pb:
xsn = 0.7
caec.set_IA_sc(isn, xsn)
caec.set_IA_sc(ipb, 1.0 - xsn)

# Set the temperature to 600 K. Under these conditions, the liquid phase
# is the only stable phase in this system (see the Pb-Sn phase diagram
# provided in worked example 3).
caec.set_eq_T(600)

# Calculate the equilibrium.
caec.calculate_eq()

# Get the integral enthalpy of the liquid phase.
hint1 = caec.get_eq_H_ph(iliq)

# Get the partial molar enthalpies of Sn/LIQUID/ and Pb/LIQUID/.
hsn = caec.get_eq_HM_pc(iliq, isn)
hpb = caec.get_eq_HM_pc(iliq, ipb)

# Calculate the integral enthalpy of the liquid phase from the two
# partial values.
hint2 = xsn * hsn + (1.0 - xsn) * hpb

# Print both values.
print(f'\nhint1 = {hint1:.3f} J')
print(f'hint2 = {hint2:.3f} J')

# As expected, the two values are the same.

Simple Process Using Streams

"""
A simple process using streams.
"""
# Note that this program contains phase target calculations, which
# cannot be performed with the 'light' version of ChemApp for Python.

from chemapp.core import(ResultVariable, ConditionVariable,
                           StreamVariable, Quantity,
                           AmountUnit, TemperatureUnit, PressureUnit,
                           VolumeUnit, EnergyUnit, TargetVariable)
from chemapp.friendly import Info as cai
from chemapp.friendly import Units as cau
from chemapp.friendly import ThermochemicalSystem as cats
from chemapp.friendly import EquilibriumCalculation as caec
from chemapp.friendly import StreamCalculation as casc

# import packages for plotting graphs.
import matplotlib.pyplot as plt
import pandas as pd

# import package for 3D-plot.
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.mlab import griddata
import scipy.interpolate

# Determine the data file path.
from pathlib import Path
import chemapp
data_dir = Path(chemapp.__file__).parents[0] / 'data'

# The following function calculates the sum of all equilibrium amounts
# of the condensed phases. It starts with phase #1, since the gas phase
# is #0. The index number of the last phase to be considered is passed
# as an argument to the function and is retrieved in the main program.

def getacp(lastph):
    getacp = 0.0
    # Add up the equilibrium amounts of all condensed phases.
    for i in range(lastph):
        amount = caec.get_eq_A_ph(i)
        getacp += amount
    return amount

# Main program.

# Initialise ChemApp for Python.
# Open data-file for reading. In this case, a data-file containing the
# thermochemical data for the system carbon-nitrogen-oxygen is selected.
# Read data-file.
# Close data-file.

cats.load(str(Path(data_dir) / 'cno.dat'))

# Since we will perform an extensive property target calculation later,
# we will use streams. This way we can assign an initial state regarding
# temperature and pressure to the incoming substances, which will be
# taken into account when ChemApp for Python calculates the extensive property
# (e.g. heat) balance.

# Since we are going to use degrees Celsius as temperature unit
# throughout the calculations, we will change the unit now.
cau.set(T = TemperatureUnit.C,
        V = VolumeUnit.dm3,
        P = PressureUnit.bar,
        E = EnergyUnit.J,
        A = AmountUnit.mol)

# According to the process scheme shown above, define the first stream.
# It is the incoming stream of mixed CO2 and CO at 200 C and 1 bar,
# which we will call 'CO2CO'. The amount of CO2 and CO in this stream
# is set to 0.5 and 2.5 mol, respectively.
casc.create_st('CO2CO', 200.0, 1.0)

# Set the constituent amounts of CO2 and CO. To do this, get the index
# number of these constituents of the gas phase, which is phase number
# 0.
ico2 = cats.get_index_pc(0, 'CO2')
ico = cats.get_index_pc(0, 'CO')

# Use the index numbers just obtained to set the constituent amount of
# CO2 in the stream 'CO2CO' to 0.5 mol.
casc.set_IA_pc('CO2CO', 0, ico2, 0.5)

# Similarly, set the amount of CO in the same stream to 2.5 mol.
casc.set_IA_pc('CO2CO', 0, ico, 2.5)

# Define the second stream ('O2N2'), which is air supplied at room
# temperature and 1 bar for the combustion. The air is assumed to
# consist of O2 and N2 at a ratio of 1:4. The amount of O2 is
# initially set to 0.5 mol (thus N2 = 2.0 mol), and will later
# be varied.
casc.create_st('O2N2', 25.0, 1.0)

# Get the index numbers of O2 and N2.
in2 = cats.get_index_pc(0, 'N2')
io2 = cats.get_index_pc(0, 'O2')

# Set the amount of O2.
casc.set_IA_pc('O2N2', 0, io2, 0.5)

# Set the amount of N2.
casc.set_IA_pc('O2N2', 0, in2, 2.0)

# Set the temperature of the reactor at which the equilibrium reaction
# (i.e. the combustion) will take place. It is initially set to 600 C
# at first, but will be varied later.
casc.set_eq_T(600.0)

# In the result table to be printed out subsequently, the amount unit
# for the condensed phases formed should be 'gram', while the input
# was given in 'mol'.
cau.set_A_unit(AmountUnit.gram)

# Calculate the equilibrium and print out the result table for a
# quick overview on the equilibrium state.
casc.calculate_eq(print_results = True)

# As can be seen from the above table, both the gas phase and carbon
# are stable under the conditions chosen. It can also be seen that
# N2, CO2, and CO are the only majority gas phase species. All
# other gas phase species are produced in negligible amounts only.

# The values for the amounts of condensed phases and the gas phase
# can be retrieved using get_count_phs(). Since the gas phase, if
# present in the data-file, has always the phase index number 0,
# this means that we have to add up all the amounts of the phases
# from phase #2 to the last one in order to retrieve the total
# amount of condensed phases.

# Find out how many phases we have in the data-file, and use this
# number in all subsequent calls to getacp.
nphase = cats.get_count_phs()

# Get the total equilibrium amount of the condensed phases, using
# the function getacp declared above.
acond = getacp(nphase)

print(f'\nAmount of condensed material: {acond:.4f} gram.')


# Retrieve the amount of gas phase. Note: The current amount unit
# is still 'gram'. If the amount of the gas phase is desired in
# 'mol' instead, tqcsu would need to be called again.
agas = casc.get_eq_A()

print(f'\nAmount of gas phase: {agas:.2f} gram.')

# Retrieve values for some extensive properties.  Since we used
# streams for this example, these extensive property values are
# calculated relative to the incoming streams, when the whole
# system is selected.

# Get the enthalpy change associated with this reaction.
result = casc.get_eq_H()

print(f'\nDelta H: {result:.4E} J')

# The negative sign of this value indicates that energy is released
# by the system.

# Get the change in volume of the process, which is the volume of
# the equilibrium gas phase minus the sum of the volumes of the
# incoming streams.
result =  casc.get_eq_V()
print(f'Delta V: {result:.2f} dm^3')

# Note: using the option 'V' with tqgetr normally would retrieve
# the _total_ volume of the system. There is one exception, in
# which case it retrieves the _change_ in volume, and that is the
# case when streams are used _and_ the whole system is selected.
# To retrieve the _total_ volume instead in the present case, use
# the option 'VT':
result =  casc.get_eq_VT()
print(f'Total V: {result:.2f} dm^3')

# Compare the two values for the volume with the corresponding
# ones in the output table produced earlier. The volume printed
# at the top of such an output table is always the _total_ volume,
# whereas the line specifying the extensive property balance toward
# the end of the table shows the _change_ in volume.

# tqstxp can be used to retrieve extensive property values for the
# incoming streams. Get the volume of the incoming CO2/CO stream:
result =  casc.get_st_V('CO2CO')
print(f'\nV(CO2CO): {result:.2f} dm^3')

# Get the value of the incoming air stream:
result =  casc.get_st_V('O2N2')
print(f'V(O2N2): {result:.3f} dm^3')

# For users of ChemApp versions prior to 3.2.0: Please note that
# the values returned by tqstxp always er to the _last_ calculated
# equilibrium state. Even if one uses tqstxp to retrieve properties
# of _incoming_ streams, one thus needs to perform an equilibrium
# calculation before the correct values are obtained.

# As can be seen from the above results, this process, under the
# conditions we have chosen, results in solids to be formed in the
# reactor. We will now further investigate this systems by looking
# at how the amount of condensed material changes when
# a) the reaction temperature is varied, and
# b) the amount of air supplied is varied.

# *********************************************************************
# Varying the reaction temperature.
# *********************************************************************

# Since several individual equilibria will be calculated, the results
# (temperature and amount of condensed phases) will be stored to a
# file for further processing.
file_11 = open('caf103_1.rst', 'w')

# Vary the temperature between 200 C and 1000 C in steps of 5 C.
for j in range(200, 1001, 5):

    # Use the variable temp to hold the current temperature.
    temp = j
    casc.set_eq_T(temp)

    # Calculate the equilibrium.
    casc.calculate_eq()


    # Get the amount of condensed phases.
    acond = getacp(nphase)

    # Write the temperature and the amount of condensed phases to the
    # result file.
    file_11.write(f'{temp},{acond}\n')

# Close the result file.
file_11.close()

# When plotted, the resulting graph looks like the following Figure.

# Read the created caf103_1.rst file, and name the first column
# 'temp' and name the second 'acond'.
data_1 = pd.read_csv('caf103_1.rst', names = ['temp', 'acond'])
plt.axis([200, 1000, 0, 10])
plt.ylabel('Amount/gram')
plt.xlabel('Temperature/C')
plt.plot(data_1['temp'], data_1['acond'], 'r-')
plt.show()

# *********************************************************************
# Varying the amount of air supplied.
# *********************************************************************

# Instead of the temperature, the amount of incoming air will be varied
# now. The temperature will be kept constant at the initial value of
# 600 C. To plot the results later on, we will again store them in a
# separate file.
file_11 = open('caf103_2.rst', 'w')

# Set the temperature for the equilibrium calculation (600 C).
casc.set_eq_T(600.0)

# Calculate 201 equilibria...
for j in range(201):

    # ... with the O2 amount varying between 0 and 1 mol...
    iao2 = j * 0.005
    # ... and the N2 amount between 0 and 4 mol.
    ian2 = iao2 * 4.0

    # To make sure ChemApp for Python interprets these values as 'mol', set the
    # amount unit accordingly.
    cau.set_A_unit(AmountUnit.mol)

    # Set the amount of O2 in the air stream.
    casc.set_IA_pc('O2N2', 0, io2, iao2)

    # Set the amount of N2 in the air stream.
    casc.set_IA_pc('O2N2', 0, in2, ian2)

    # Calculate the equilibrium.
    casc.calculate_eq()

    # Set the amount unit to 'gram'.
    cau.set_A_unit(AmountUnit.gram)

    # Get the amount of condensed phases.
    acond = getacp(nphase)

    # Write the amounts of O2 and condensed phases to the result file.
    file_11.write(f'{iao2},{acond}\n')

# Close the result file.
file_11.close()

# When plotted, the resulting graph looks like the following Figure:
data_2 = pd.read_csv('caf103_2.rst', names = ['iao2', 'acond'])
plt.axis([0, 1, 0, 12])
plt.ylabel('Amount/gram')
plt.xlabel('n(O2)/mol')
plt.plot(data_2['iao2'], data_2['acond'], 'r-')
plt.show()

# *********************************************************************
# Varying the amount of air supplied _and_ the temperature.
# *********************************************************************

# To investigate the combined influence of temperature and amount of air
# supplied on the resulting amount of condensed phases, both parameters
# are varied in nested loops. The results are again stored in a file for
# further processing.
file_11 = open('caf103_3.rst', 'w')

# The outer loop will vary the temperature between 200 C and 1000 C.
for j in range(200, 1001, 20):

    temp = j
    casc.set_eq_T(temp)

    # The inner loop varies the amount of incoming air as above.
    for k in range(0, 201, 5):
        iao2 = k * 0.005
        ian2 = iao2 * 4.0

        cau.set_A_unit(AmountUnit.mol)
        # Set the incoming amount of O2.
        casc.set_IA_pc('O2N2', 0, io2, iao2)

        # Set the incoming amount of N2.
        casc.set_IA_pc('O2N2', 0, in2, ian2)

        casc.calculate_eq()
        cau.set_A_unit(AmountUnit.gram)

        # Get the amount of condensed phases
        acond = getacp(nphase)

        # Write the temperature and the amounts of O2 and condensed
        # phases to the result file.
        file_11.write(f'{temp},{iao2},{acond}\n')

# Close the result file.
file_11.close()

# When plotted, the resulting graph looks like the following Figure:

data_3 = pd.read_csv('caf103_3.rst', names = ['temp', 'iao2', 'acond'])
fig = plt.figure()

ax = fig.add_subplot(111, projection='3d')
ax.set_xlim3d(200, 1000)
ax.set_ylim3d(0, 1)
ax.set_zlim3d(-10, 20)
ax.set_xlabel('Temperature/C')
ax.set_ylabel('n(O2)/mol')
ax.set_zlabel('Amount/gram')

x = data_3['temp']
y = data_3['iao2']
z = data_3['acond']
ax.scatter(x, y, z, color = 'red')
plt.show()

# Check if we are working with the 'light' version. If we do, omit
# the following phase target calculation(s).
islite = cai.is_light()
if islite == True:
    print(f'*** Phase target calculations have been omitted here,\n '
          f'*** since they are not possible with the\n '
          f'*** ''light''version of ChemApp for Python.')

else:

    # *****************************************************************
    # Varying the amount of air supplied _and_ the temperature,
    # searching for the explicit condition where the amount of condensed
    # phase is zero.
    # *****************************************************************

    # Assuming that one is interested in the above process _not_
    # producing any condensed phases (i.e. complete combustion), it would
    # be helpful to have a diagram which corresponds to the bottom plane
    # of the above figure: A graph which, as a function of the incoming
    # amount of air, shows the temperature at which the amount of
    # condensed phases formed is zero.

    # For the calculation procedure with ChemApp for Python, this results in the
    # following steps:
    # 1) We are looking for the temperature at which condensed phases
    # _just_ start to form. In terms of ChemApp for Python, we thus have to define
    # a precipitation target with respect to the gas phase.
    # 2) For each value of the amount of air supplied, define the 'O2N2'
    # stream with the proper amounts of O2 and N2.
    # 3) ChemApp for Python is then supposed to vary the temperature until the
    # target (a phase precipitates from the gas phase) is met.
    # Temperature is thus the target variable.
    # 4) The equilibrium is calculated for each of the amounts of air
    # supplied, then the equilibrium temperature is retrieved and stored
    # in the result file.
    file_11 = open('caf103_4.rst', 'w')

    # Reset the amount unit to 'mol'.
    cau.set_A_unit(AmountUnit.mol)

    # Define the target as a precipitation target for the gas
    # phase (phase index 0).
    casc.set_target_precipitation_from_ph(0)

    # The initial estimate for the temperature is passed via the first
    # element of vals.
    vals = 600.0

    # Vary the amount of supplied air.
    for j in range(201):
        iao2 = j * 0.005
        ian2 = iao2 * 4.0

        # Set the incoming amount of O2.
        casc.set_IA_pc('O2N2', 0, io2, iao2)

        # Set the incoming amount of N2.
        casc.set_IA_pc('O2N2', 0, in2, ian2)

        # Define the target variable (temperature), with the initial
        # estimate in the first element of vals, and calculate the
        # equilibrium.
        casc.calculate_eq_T(vals)

        # Retrieve the calculated equilibrium temperature.
        temp = casc.get_eq_T()

        # Store the amount of O2 supplied and the calculated
        # equilibrium temperature to the result file.
        file_11.write(f'{iao2},{temp}\n')
    # Close data-file.
    file_11.close()

    # When plotted, the resulting graph looks like the following Figure.
    data_4 = pd.read_csv('caf103_4.rst', names = ['iao2', 'temp'])
    plt.axis([0, 1, 500, 800])
    plt.ylabel('Temperature/C')
    plt.xlabel('n(O2)/mol')
    plt.plot(data_4['iao2'], data_4['temp'], 'r-')
    plt.show()

One-dimensional Phase Mapping

"""
One-dimensional phase mapping calculations.

Note that this program contains phase target calculations, which cannot be
performed with the 'light' version of ChemApp for Python.
"""

from chemapp.core import (Status,
    ConditionVariable, TargetVariable, ResultVariable, PhaseMapVariable,
    AmountUnit)
from chemapp.friendly import Units as cau
from chemapp.friendly import Info as cai
from chemapp.friendly import ThermochemicalSystem as cats
from chemapp.friendly import EquilibriumCalculation as caec
from chemapp.friendly import PhaseMapCalculation as capmc

# Determine the data file path.
from pathlib import Path
import chemapp
data_dir = Path(chemapp.__file__).parents[0] / 'data'

# import packages for plotting graphs.
import matplotlib.pyplot as plt

# The function psp prints the value of the 'search variable' (temperature,
# pressure, or composition), the names of the stable phases, and their amounts.
# It is used after every call to calculate_ph_txs.
# The variable_str is set to 'Temperature/K' for temperature calls,
# 'Pressure/bar' for pressure etc. This function also prints a header for the result
# table that is populated by psp.
def psp(variable_str, variable, results):
    print(f'\n{variable_str:<26}{"Stable phases":<18}{"Amount/mol"}')

    eps = 1E-8

    for r in results:
        print(f'{getattr(r, variable):.2f}')

        for ph in r.phs.values():
            if ph.A > eps:
                print(f'{"":<26}{ph.name:<24}{ph.A:.2f}')


# Since the calculate_ph_txs_IA is different, a psp_ia function was written and
# it works on the same principles as psp.
def psp_ia(results):
    print(f'\n{"X(Fe2SiO4)":<26}{"X(Mg2SiO4)":<26}{"Stable phases":<18}'
          f'{"Amount/mol"}')

    eps = 1E-8

    for r in results:
        print(f'{r.IAs[0].A:<26.3f}{r.IAs[1].A:<26.3f}')

        for ph in r.phs.values():
            if ph.A > eps:
                print(f'{"":<52}{ph.name:<24}{ph.A:.2f}')


asize = 100
eps = 1E-3
nused = 0

# Initialise ChemApp for Python, open, read and close the thermochemical data-file
# pbsn.dat (system Pb-Sn).
cats.load(str(data_dir / 'pbsn.dat'))


# Check if we are working with the 'light' version.
# If we are, omit the following phase target calculation(s).
islite =  cai.is_light()

if islite == True:
    print(f'*** Phase target calculations have been omitted here,\n'
          f'*** since they are not possible with the\n'
          f'*** ''light'' version of ChemApp for Python.')

    exit(0)


# The data-file for system Pb-Sn contains a second copy of the BCT_A5 phase,
# because of the inherent metastable miscibility gap. Since this miscibility
# gap will stay metastable under the conditions used in this example, we
# will eliminate the second copy of that phase, which means it will be
# disregarded completely in the following calculations.
bct2 = cats.get_index_ph('BCT_A5#2')
cats.set_status_ph(bct2, Status.ELIMINATED)

# Determine the index numbers for the system components Pb and Sn.
ipb = cats.get_index_sc('Pb')
isn = cats.get_index_sc('Sn')

# As indicated in the Pb-Sn phase diagram, label "1", the
# first phase map should be done with the temperature as the search
# variable, for X(Sn)=0.2.

# Set the incoming amounts for the Sn and Pb.
capmc.set_IA_sc(isn, 0.2)
capmc.set_IA_sc(ipb, 0.8)

tarray = list(range(410, 600, 10))

# The purpose of the array tarray is to hold the temperature values at which
# the equilibria will be calculated later.
# This set of temperature values consists of,
# 1) a number of evenly spaced temperature intervals,
# 2) and the temperatures at which calculate_ph_txs_T finds phase boundaries.

# One call is made to ChemApp for Python friendly, and all of the transitions are
# calculated and sorted.
results = capmc.calculate_ph_txs_T(400.0, 600.0)
psp('Temperature/K', 'T', results)

# We are now ready to calculate the equilibria for all the temperatures in
# tarray. These equilibrium calculations now being performed are _not_
# target calculations of any kind any more, that's why we reset _all_
# previous conditions including all targets and target variables.
caec.remove_eq_conditions_all()

# The following equilibria are to be calculated for the same alloy
#  concentration of course, so we use the same incoming amounts again.
caec.set_IA_sc(isn, 0.2)
caec.set_IA_sc(ipb, 0.8)

# The graphical representation of the one-dimensional phase map will be
# done using two x/y graphs: one displaying the phase amounts vs. the
# temperature, the other one the phase activities vs. the
# temperature.
# For every temperature in tarray, calculate the equilibrium:
for t in tarray:
    caec.set_eq_T(t)

    # Note again that these are regular equilibrium calculations, no
    # targets are involved:
    caec.calculate_eq()
    results.append(caec.get_result_object())

# Sort the result list by temperature. This is required because results are appended
# to the list.
results.sort(key=lambda r: r.T)

# When plotted, the two graphs, displaying the equilibrium phase
# amounts and phase activities vs. temperature look like Figures
# 3 and 4.

T_array = [r.T for r in results]
A_liq_array = [r.phs['LIQUID'].A for r in results]
A_bct_array = [r.phs['BCT_A5#1'].A for r in results]
A_fcc_array = [r.phs['FCC_A1'].A for r in results]

# Print the sorted values in T_array in ascending order.
for i in range(0, len(T_array)):
    print(f'TARRAY({i:>2}) = {T_array[i]:>8.5f}')

# Phase amounts vs. temperature.
plt.axis([400, 600, 0, 1])
plt.ylabel('Amount/mol')
plt.xlabel('Temperature/K')
plt.plot(T_array, A_liq_array, 'r-+',
         linewidth=0.5, markersize=5, label='LIQUID')
plt.plot(T_array, A_bct_array, 'g-x',
         linewidth=0.5, markersize=5, label='BCT_A5#1')
plt.plot(T_array, A_fcc_array, 'm-s',
         linewidth=0.5, markersize=5, markerfacecolor='none',
         label='FCC_A1')
plt.legend(loc='best')
# plt.show()

AC_liq_array = [r.phs['LIQUID'].AC for r in results]
AC_bct_array = [r.phs['BCT_A5#1'].AC for r in results]
AC_fcc_array = [r.phs['FCC_A1'].AC for r in results]

# Activities vs. temperature.
plt.axis([400, 600, 0.55, 1])
plt.ylabel('Activity')
plt.xlabel('Temperature/K')
plt.plot(T_array, AC_liq_array, 'r-+',
         linewidth=0.5, markersize=5, label='LIQUID')
plt.plot(T_array, AC_bct_array, 'g-x',
         linewidth=0.5, markersize=5, label='BCT_A5#1')
plt.plot(T_array, AC_fcc_array, 'm-s',
         linewidth=0.5, markersize=5, markerfacecolor='none',
         label='FCC_A1')
plt.legend(loc='best')
# plt.show()

# For the next one-dimensional phase mapping calculations we will only
# use tabular output. First, remove all previously set conditions.
capmc.remove_eq_conditions_all()

# We will now calculate another one-dimensional phase map in the same
# system, again with temperature as the search variable. The new
# alloy composition will cause calculate_ph_txs_T to find the eutectic.
capmc.set_IA_sc(isn, 0.5)
capmc.set_IA_sc(ipb, 0.5)

# For this, we call calculate_ph_txs_T to calculate the phase transitions
# according to temperature.
results = capmc.calculate_ph_txs_T(400.0, 600.0)
psp('Temperature/K', 'T', results)

# Note that the above table lists a 3-phase equilibria now, which for the
# current system means we found the eutectic temperature (compare this
# with the results for the same equilibrium determined in worked example
# 3: Phase equilibria and phase target calculations).

# So far we have calculated one-dimensional phase maps for a
# temperature/composition phase diagram, in our case keeping the
# composition (and implicitly the pressure) constant and using the
# temperature as the search variable. The next phase map we will
# calculate at a constant temperature, this time varying the
# pressure. For this purpose, we will use a different thermochemical
# data-file for the system Fe2SiO4-Mg2SiO4 which includes
# pressure-dependent Gibbs-energy data.

# Open, Read and close the thermochemical data-file femgsio4.dat
# (system Fe2SiO4-Mg2SiO4).
cats.load(str(data_dir / 'femgsio4.dat'))

# Determine the index numbers for the stoichiometric compounds Fe2SiO4
# and Mg2SiO4. These two phases will be used to define the incoming
# amounts for the next calculations. Note that they are included in the
# thermochemical data-file only for this purpose (to allow for a
# somewhat easier definition of incoming amounts), since the Gibbs
# energy functions for these two stoichiometric compounds are only
# dummies.

# Select a composition of 2 mol-% Fe2SiO4, the rest is Mg2SiO4:
capmc.set_IA_pc('Fe2SiO4', 0, 0.02)
capmc.set_IA_pc('Mg2SiO4', 0, 0.98)

# Select the range for the search variable, which is the pressure in the
# current case. We will search for any phase boundaries between 120 kbar
# and 150 kbar (see also the line marked '1' in Figure 2.)

# The temperature is set to 1000 K.
capmc.set_eq_T(1000.0)

# Call calculate_ph_txs_P to calculate the transitions pressure.
results = capmc.calculate_ph_txs_P(120000.0, 150000.0)

# Call psp (print stable phases), passing the value for the pressure
# we just retrieved. psp will print a formatted entry for the result
# table.
psp('Pressure/bar', 'P', results)

# The result table is complete now. Compare the phase boundaries found
# with the phase diagram in Figure 2. Note that the total
# amount in the system is always calculated to be 2 mol. This is due to
# the fact that the phase constituents of all mixture phases in this
# data-file (Spinell, Olivine, and \DF-phase) have been scaled with a
# factor of 0.5 in order to obtain the correct ideal entropy between Fe and
# Mg.

# The last one-dimensional phase map will use composition as the search
# variable, keeping both temperature and pressure constant. As with
# other phase target calculations which use composition as the target
# variable, the input in this case is a bit more complicated (see also
# worked example 3: Phase equilibria and phase target calculations).
# This is due to the fact that when one is interested in comparing the
# result for this type of calculation with graphical phase diagram
# information, it is necessary to keep the total substance amount in the
# system constant (i.e. use concentrations rather than absolute
# amounts).

# First, reset all conditions and targets again:
capmc.remove_eq_condition(-1)

# Set the pressure to 135 kbar (see Figure 2):
capmc.set_eq_P(135000.0)

# The composition range for Fe2SiO4 has a minimum of 0 mol and a maximum of 0.2
# mol, compare again to Figure 2. The inputs required by calculate_ph_txs_IA
# are defined in the list "ranges". Here the ph, pc (None for this example),
# range_start and range_end are defined.

# The results are printed by psp_ia.
ranges = [capmc.IA_range('Fe2SiO4', '', 0.0, 0.2),
          capmc.IA_range('Mg2SiO4', '', 1.0, 0.8)]
results = capmc.calculate_ph_txs_IA(ranges)
psp_ia(results)

Tutorials

This section provides examples to show you how to use the ChemApp for Python friendly module function groups. The examples are small bits of code that you can use to build a full size investigation.

Working with Info Functions

The easiest way to get all the information about ChemApp, is by requesting a ChemAppInfo object, and printing it. Here is an example, code: 30-friendly/Info/01-basic-use.py

"""
This script gives examples of how to get a info object
"""
from chemapp.friendly import Info as cai

# extract all information related to chemapp into a ChemAppInfo object
info = cai.get_object()

# print the object
print(info)

Here is the result:

========================================================================
 ChemApp Information
------------------------------------------------------------------------
Class:           ChemAppInfo

Program ID:      CAFU
Version:         816
Is Light:        False
Expiration date: 2022/7

User ID:         1234
User Name:       Acme Ltd.

Dongle Type:
Dongle ID:       0

Copyright:
This program contains ChemApp.
Copyright GTT-Technologies, Kaiserstrasse 103, D-52134 Herzogenrath, Germany.
http://www.gtt-technologies.de
-------------------------------------------------------------------------
Array Sizes:
Item   Maximum   Used  Description
na       3000    3000  Number of constituents
nb         48      48  Number of system components
nc        150     150  Number of mixture phases
nd       2000    2000  Number of excess Gibbs energy coefficients for a
                       mixture phase
ne        500     500  Number of excess magnetic coefficients for a
                       mixture phase
nf          5       5  Number of sublattices for a mixture phase
ng       3000    3000  Number of constituents of a sublattice
nh         48      48  Number of oxide constituents of a phase described
                       by the Gaye-Kapoor-Frohberg or modified
                       quasichemical formalisms
ni         20      20  Number of Gibbs energy/heat capacity equations
                       for a constituent
nj      20000   20000  Number of Gibbs energy/heat capacity equations
nk       3000    3000  Number of constituents with P,T-dependent molar
                       volumes
===========================================================================

You can also use individual functions in the Info group, and you can use individual properties of ChemAppInfo objects that you obtain from the Info group. This provides ample flexibility for doing what you need. You can use the help function in python to find other functions of the Info functions.

Working with Units Functions

Checking Units

It is good practise to check the current units in ChemApp for Python before doing calculations. This can be done by printing them. Here is an example based on this example code: 30-friendly/Units/01-basic-use.py

"""
This script gives examples of how to use the core unit
"""
from chemapp.friendly import Units

# check the current units in ChemAppPy
units = Units.get_str()

# print the units
print(units)

Here is the result:

Units Currently Used by ChemApp
===============================
Pressure:    bar
Volume:      dm3
Temperature: K
Energy:      J
Amount:      mol

Changing Units

It is simple to change any of the units in ChemApp for Python. To change the pressure unit, here is an example, the code is found here: 30-friendly/Units/02-change-unit.py

"""
This script gives examples of how to use the core unit
"""
from chemapp.friendly import Units, PressureUnit

# set the pressure units to atmospheres
Units.set_P_unit(PressureUnit.atm)

# print to check that the units have been changed
print(Units.get_P_unit())

To check that the units are now correct you can print the unit again. The result is:

PressureUnit.atm

Working with Thermochemical Functions

Load thermochemical datafile

Before being able to do equilibrium calculations with ChemApp for Python you need to create a thermochemical system. You do this by loading a thermochemical data file, see the example below. This refers to 30-friendly/ThermochemicalSystem/01-basic-use.py

"""
ThermochemicalSystem Example 1
==============================

This example demonstrates the most basic steps involved in working with the
``friendly`` module's ``ThermochemicalSystem`` group of functions.
"""

from chemapp.friendly import (
    ThermochemicalSystem)


# load a thermochemical system from a data file
ThermochemicalSystem.load('../../data/cosi.dat')

# convert the thermochemical system to a core.ThermochemicalSystem object
ts = ThermochemicalSystem.get_object()

# print the core.ThermochemicalSystem object
print(ts)

To check if the thermochemical data file was loaded you can print it. It is useful to save the printed information in a text file so that you can refer to it later.

Getting the molar mass

You can get the molar mass of a particular phase, phase constituent or system component. This example is for a phase constituent, and the example code can be found here: 30-friendly/ThermochemicalSystem/055-get-molar-mass.py.

"""
This script gets the molar mass of a phase constituent
"""

from chemapp.friendly import ThermochemicalSystem, Units, AmountUnit


# load a thermochemical system from a data file
ThermochemicalSystem.load('../../data/fec.dat')

Units.set_A_unit(AmountUnit.gram)

# get the molar mass of O2 in the gas phase
mm_c = ThermochemicalSystem.get_mm_sc('C')

# print the molar mass
print(mm_c*1000.0)

The molar mass is:

12.010700000000002

Get Y (mass fraction) of a system component or phase constituent

You can get the mass fraction of a phase constituent in a phase or the mass fraction of a system component in a phase constituent. The example below shows how to calculate the mass fraction of carbon in carbon monoxide in the gas phase. This example refers to 30-friendly/ThermochemicalSystem/03-get-mass-fraction.py:

"""
This script gives an example of how to get the mass fraction of a system
component in a phase constituent
"""

from chemapp.friendly import ThermochemicalSystem

# load a thermochemical system from a data file
ThermochemicalSystem.load('../../data/cosi.dat')

# get mass fraction of C in CO
y_c = ThermochemicalSystem.get_y_sc_pc('GAS', 'CO', 'C')

# print result
print(y_c)

Result:

0.42880501528003884

Set the status of a phase

Before running a calculation, you can set the status of a particular phase or phase constituent. The options are ENTERED, ELIMINTATED or DORMANT. Here is an example. The code is: 30-friendly/ThermochemicalSystem/04-set-phase-status.py

"""
This scripts sets the status of a phase in ChemAppPy
"""

from chemapp.friendly import ThermochemicalSystem, Status

# load a thermochemical system from a data file
ThermochemicalSystem.load('../../data/fec.dat')

# set the status of the FCC phase as ENTERED
ThermochemicalSystem.set_status_ph('FCC_A1', Status.ENTERED)

# check if the status is now ENTERED
status = ThermochemicalSystem.get_status_ph('FCC_A1')

# print status
print(status)

The result:

Status.ENTERED

Working with EquilibriumCalculation Functions

Equilibrium calculation

One can do simple equilibrium calculation in ChemApp for Python. The example below shows how to calculate the mass of liquid alloy formed when heating up cementite with pure carbon. Here is an example, the code is here: 30-friendly/EquilibriumCalculation/01-basic-use.py.

"""
This scripts shows how to do a simple equilibrium calculation with
EquilibriumCalculation
"""

from chemapp.friendly import (
    ThermochemicalSystem as cats,
    EquilibriumCalculation as caec,
    Status, Units)

# load a thermochemical system from a data file
cats.load('../../data/fec.dat')
Units.set()

# convert the thermochemical system to a core.ThermochemicalSystem object
ts = cats.get_object()

# set input parameters
T_feed = 25.00  # C, feed temperature
T_equilib = 1600.00 # C, equilibrium temperature
P_feed = 1.0  # atm, feed pressure
P_equilib = 1.0  # atm, equilibrium pressure
m_fec = 500.0  # kg, mass of cementite added
m_c = 1.0  # kg, mass of graphite added

# set input parameters in ChemApp for Python with EquilibriumCalculation
caec.set_IA_pc('Fe3C_CEMENTITE', 'Fe3C_CEMENTITE', m_fec)
caec.set_IA_pc('C_GRAPHITE', 'C_GRAPHITE', m_c)

# set equilibrium conditions
caec.set_eq_T(T_equilib)
caec.set_eq_P(P_equilib)

# calculate equilibrium
caec.calculate_eq(print_results=False)

# get amount of liquid alloy formed
m_liq = caec.get_eq_A_ph('LIQUID')

# print result
print('Liquid-alloy: %f kg' % m_liq )

The mass of liquid alloy formed is:

Liquid-alloy: 493.640494 kg

Target equilibrium calculation

It is also possible to set the target formation phase for the calculation. Here is an example to do this, code: 30-friendly/EquilibriumCalculation/02-set-target-eq-condition.py In this example we calculate the temperature at which liquid alloy starts forming:

"""
This scripts shows how to do a simple equilibrium calculation with
EquilibriumCalculation
"""

from chemapp.friendly import (
    ThermochemicalSystem as cats,
    EquilibriumCalculation as caec,
    Info as cai,
    Units)

# load a thermochemical system from a data file
cats.load('../../data/fec.dat')
Units.set()

# Check if we are working with the 'light' version.
# If we do, omit the following phase target calculation(s).
if cai.is_light():
    raise ValueError(f'Target calculations are not supported when using ChemApp Light.'
                     ' This example is included in the documentation to illustrate the'
                     ' abilities available with the full version of ChemApp.')


# convert the thermochemical system to a core.ThermochemicalSystem object
ts = cats.get_object()

# set input parameters
T_feed = 25.00  # C, feed temperature
T_equilib = 1600.00 # C, equilibrium temperature
P_feed = 1.0  # atm, feed pressure
P_equilib = 1.0  # atm, equilibrium pressure
m_fec = 500.0  # kg, mass of cementite added
m_c = 1.0  # kg, mass of graphite added

# set input parameters in ChemApp for Python with EquilibriumCalculation
caec.set_IA_pc('Fe3C_CEMENTITE', 'Fe3C_CEMENTITE', m_fec)
caec.set_IA_pc('C_GRAPHITE', 'C_GRAPHITE', m_c)

# set equilibrium conditions
caec.set_target_formation_of_ph('LIQUID')
caec.set_eq_P(P_equilib)

# calculate equilibrium
caec.calculate_eq_T(T_equilib)

# get amount of liquid alloy formed
T_formation = caec.get_eq_T()

# print result
print(T_formation)

Result:

1153.4357094806123 C

Working with StreamCalculation Functions

Equilibrium calculation with two streams

An example of a simple stream calculation is to calculate the amount of energy required or change in energy when smelting iron oxide with graphite at a high temperature. Here is an example, see code: 30-friendly/StreamCalculation/01-basic-use.py

"""
This script is a basic example of how an stream equilibrium calculation is done.
"""

from chemapp.friendly import \
    ThermochemicalSystem as cats, \
    StreamCalculation

# load a thermochemical system from a data file
cats.load('../../data/fec.dat')

# get thermochemical system as object
ts = cats.get_object()

# set parameters
T_amb = 25.00 # C
P_amb = 1 # atm
T_equilib = 1500.00 # C
P_equilib = 1 # atm
m_red = 1.0 # kg
m_ore = 500.0 # kg

y_feo = 0.4 # kg/kg, mass fraction feo in ore
y_fe2o3 = 0.6 # kg/kg, mass fraction fe2o3 in ore

m_feo = m_ore*y_feo
m_fe2o3 = m_ore*y_fe2o3

# create ore stream
StreamCalculation.create_st('ore', T_amb, P_amb)
StreamCalculation.set_IA_pc('ore', 'Fe3C_CEMENTITE', 'Fe3C_CEMENTITE',
                            m_feo)


# create reductant stream
StreamCalculation.create_st('reductant', T_amb, P_amb)
StreamCalculation.set_IA_pc('reductant', 'C_GRAPHITE', 'C_GRAPHITE',
                            m_red)

# set equilibrium conditions
StreamCalculation.set_eq_T(T_equilib)
StreamCalculation.set_eq_P(P_equilib)

# calculate equilibrium
StreamCalculation.calculate_eq(print_results=False)

# calculate change in enthalpy
dH = StreamCalculation.get_eq_H()

# print result
print(f'dH reaction: {dH} kWh')

Results:

dH reaction: 73.290423 kWh

Target calculation with stream calculation

You can also calculate the equilibrium by seting a formation or precipitation phase target. This is an example of a target precipitation phase, 30-friendly/StreamCalculation/02-set-target-eq-condition.py

"""
This script is a basic example of how an stream equilibrium calculation is
done with a target precipitation phase.
"""

from chemapp.friendly import \
    ThermochemicalSystem as cats, \
    StreamCalculation

# load a thermochemical system from a data file
cats.load('../../data/fec.dat')

# get thermochemical system as object
ts = cats.get_object()

# set parameters
T_amb = 1600.00  # C
P_amb = 1.0  # atm
T_equilib = 500.00  # C
P_equilib = 1  # atm
m_alloy = 500.0  # kg

y_fe = 0.8  # kg/kg, mass fraction fe in alloy
y_c = 0.05  # kg/kg, mass fraction carbon in alloy

m_fe = m_alloy*y_fe
m_c = m_alloy*y_c

# create alloy stream
StreamCalculation.create_st('alloy', T_amb, P_amb)
StreamCalculation.set_IA_pc('alloy', 'LIQUID', 'Fe',
                            m_fe)
StreamCalculation.set_IA_pc('alloy', 'LIQUID', 'C',
                            m_c)

# set equilibrium conditions
# set precipitation target phase as the liquid alloy
StreamCalculation.set_target_formation_of_ph('LIQUID')
StreamCalculation.set_eq_P(P_equilib)

# calculate equilibrium
StreamCalculation.calculate_eq_T(T_equilib)

# get the temperature at which precipitation occurs
T_precept = StreamCalculation.get_eq_T()

# print result
print(T_precept)

The result:

1153.4358038309547

Working with PhaseMap Functions

These functions allow one-dimensional phase mapping. To use these functions, one of three variables can be adjusted; pressure, temperature or composition.

Equilibrium calculation with two streams

Phase transitions can be determined by varying the temperature. See the example, 30-friendly/PhaseMapCalculation/01-basic-use.py.

from chemapp.friendly import (
    ThermochemicalSystem as cats,
    PhaseMapCalculation as pmc,
    Units)
from chemapp.core import (
    AmountUnit,
    TemperatureUnit)
import chemapp

from pathlib import Path

# Local ChemApp for Python data directory
cap_path = Path(chemapp.__file__).parent
data_path = cap_path / "data"

# Load the PbSn thermochemical system
cats.load(data_path / 'pbsn.dat')

# Set amount units to kg and temperature units to K
Units.set_A_unit(AmountUnit.kg)
Units.set_T_unit(TemperatureUnit.K)

# Get and print the thermochemical system
ts = cats.get_object()
print(ts)

# Set incoming amounts to 0.8 mol Pb and 0.2 mol Sn
pmc.set_IA_sc("Pb", 0.8)
pmc.set_IA_sc("Sn", 0.2)

# Calculate the one dimensional phase boundaries when varying temperature
results = pmc.calculate_ph_txs_T(300, 700, False)

# Print the results in table format
template = "{TEMP:<30}{PHASES:<30}{AMOUNT:<30}"
print(template.format(
    TEMP="Temperature(K)",
    PHASES="Stable phases",
    AMOUNT="Amount(kg)"))

for result in results:
    print(template.format(TEMP=result.T,
                          PHASES="",
                          AMOUNT=""))
    stable_phases = result.get_eq_phs()
    for stable_phase in stable_phases:
        print(template.format(TEMP="",
                              PHASES=stable_phase.name,
                              AMOUNT=stable_phase.A))
Temperature (K)    Stable Phases      Amount (mol)
300.000000
                    BCT_A5#1             0.178784
                    FCC_A1             0.821216
454.561924
                    BCT_A5#1             0.038190
                    FCC_A1             0.961810
541.876513
                    LIQUID             1.000000
700.000000
                    LIQUID             1.000000

Vary incoming amounts to get phase transition

You can use composition to find the phase transitions as well. In this case you will keep temperature and pressure constant. See this example, 30-friendly/PhaseMapCalculation/02-incoming-amounts-variation.py

"""
This script gives an example of how to do a phase map calculation with ChemApp for Python
"""

from chemapp.friendly import \
    ThermochemicalSystem as cats, \
    PhaseMapCalculation as pmc

# load a thermochemical system from a data file
cats.load('../../data/femgsio4.dat')

# set parameters
P_equilib = 1  # atm

# set the equilibrium pressure
pmc.set_eq_P(P_equilib)

# create ranges of incoming amounts
fe2sio4_range = pmc.IA_range('Fe2SiO4', '', 0.0, 0.2)
mg2sio4_range = pmc.IA_range('Mg2SiO4', '', 1.0, 0.8)
ranges = [fe2sio4_range, mg2sio4_range]

# calculate phase map results
results = pmc.calculate_ph_txs_IA(ranges, print_results=False)

# Print the results in table format
print('X(Fe2SiO4)     X(Mg2SiO4)       Stable Phases        Amount (mol)')
for result in results:
    print('%f             %f' % (result.IAs[0].A, result.IAs[1].A))
    stable_phases = result.get_eq_phs()
    for stable_phase in stable_phases:
        print('                                %s               %f' % (
            stable_phase.name, stable_phase.A))

The result:

X(Fe2SiO4)     X(Mg2SiO4)       Stable Phases        Amount (mol)
0.000000             1.000000
                                Olivine               1.000000
0.200000             0.800000
                                Olivine               1.000000