import logging
import pprint

from molmass import Formula

from atmodeller import (
    ChemicalSpecies,
    EquilibriumModel,
    SpeciesNetwork,
    ThermodynamicState,
    debug_logger,
)
from atmodeller.thermodata.core import CondensateActivity

logger = debug_logger()
logger.setLevel(logging.INFO)

# For more output use DEBUG
# logger.setLevel(logging.DEBUG)

Atmodeller initialized with double precision (float64)

Gas Mixing

This notebook is available at notebooks/gas_mixing.ipynb and is easiest to obtain by downloading the source code.

In gas mixing experiments, a predefined gaseous mixture is injected into a chamber and allowed to equilibrate thermodynamically under controlled temperature-pressure conditions. In this example, we consider a gas mixture that equilibrates in the presence of solid carbon (graphite).

# Enforce the stability of graphite
# Since in this example we do not provide carbon in the injected gas stream, we cannot solve for
# the stability of any carbon-bearing products because in order to do so requires specification of
# the mass of carbon in the system.
activity = CondensateActivity(1.0)
C_cr = ChemicalSpecies.create_condensed("C", activity=activity, solve_for_stability=False)

# Define allowable gas species at equilibrium
H2_g = ChemicalSpecies.create_gas("H2")
N2_g = ChemicalSpecies.create_gas("N2")
CH4_g = ChemicalSpecies.create_gas("CH4")
CHN_g = ChemicalSpecies.create_gas("CHN")
H_g = ChemicalSpecies.create_gas("H")

species = SpeciesNetwork((C_cr, H2_g, N2_g, CH4_g, CHN_g, H_g))

model = EquilibriumModel(species)

# Set the temperature and pressure
state = ThermodynamicState(temperature=1773.15, melt_fraction=0, pressure=1)

# Define the mole fractions of input gases
mole_fractions = {"H2": 0.5, "N2": 0.5}

# Define the composition of the input gas mixture by mass
mass_constraints = {key: value * Formula(key).mass for key, value in mole_fractions.items()}

# Solve
model.solve(state=state, mass_constraints=mass_constraints, solver="basic")

# The abundance of C cannot be solved for, so the output result for the number density of graphite
# should be ignored.
output = model.output
solution = output.quick_look()

logger.info("solution = %s", pprint.pformat(solution))

# Optionally dump the output to Excel
# output.to_excel("gas_mixing")

[19:34:59 - atmodeller.classes             - INFO     ] - species_network = ('C_cd: CondensateActivity, NoSolubility', 'H2_g: IdealGas, NoSolubility', 'N2_g: IdealGas, NoSolubility', 'CH4_g: IdealGas, NoSolubility', 'CHN_g: IdealGas, NoSolubility', 'H_g: IdealGas, NoSolubility')
[19:34:59 - atmodeller.classes             - INFO     ] - Thermodynamic data requires temperatures between 200 K and 6000 K
[19:34:59 - atmodeller.classes             - INFO     ] - reactions = {0: '1.0 C_cd + 2.0 H2_g = 1.0 CH4_g',
 1: '1.0 C_cd + 0.5 H2_g + 0.5 N2_g = 1.0 CHN_g',
 2: '0.5 H2_g = 1.0 H_g'}
[19:35:03 - atmodeller.classes             - INFO     ] - Solve (basic) complete: 1 (100.00%) successful model(s)
[19:35:03 - atmodeller.classes             - INFO     ] - Multistart summary: 1 (100.00%) models(s) required 1 attempt(s)
[19:35:03 - atmodeller.classes             - INFO     ] - Solver steps (max) = 53
[19:35:04 - atmodeller                     - INFO     ] - solution = {'CH4_g': array(0.000212698488103),
 'CH4_g_activity': array(0.000212698488103),
 'CHN_g': array(0.003057989685365),
 'CHN_g_activity': array(0.003057989685365),
 'C_cd': array(5.185286445067471e+18),
 'C_cd_activity': array(1.),
 'H2_g': array(0.498000976078994),
 'H2_g_activity': array(0.498000976078993),
 'H_g': array(0.000201308461558),
 'H_g_activity': array(0.000201308461558),
 'N2_g': array(0.49852702728598),
 'N2_g_activity': array(0.498527027285979)}