From charlesreid1

This page contains gas mixing functions using Cantera, written in Python.

Many of these have been added to Pantera, a monkey-patch and convenience function library that wraps Cantera.

Types

The thermochemical of a gas state is determined by its temperature, pressure, and composition. The thermochemical state of a gas can be represented/stored using the Cantera library in a variety of ways:

  • Cantera Gas object - inherent to the Cantera Gas object is gas.saveState(), which can be used in conjunction with gas.restoreState() to make copies of gases, a la gas2.restoreState( gas.saveState() ). These methods transcieve all information required to specify the thermochemical state.
  • Composition string - Temperature and pressure specified, and composition specified using a string, e.g. "CH4:0.5, O2:1.0, N2:3.76"
  • Composition dict - Temperature and pressure specified, and composition specified using a Python dict type that would look something like this:
d = { 'CH4': 0.5,
 'N2': 3.76,
 'O2': 1.0 }

(As actually implemented, these amounts are normalized by the total moles. The representation ultimately doesn't matter, because gas objects are only representations of thermochemical states, and not of actual material amounts of gas (that kind of representation requires a Cantera Reactor object).

Functions that use Cantera Gas objects

It would be useful to create a set of functions that would turn a Cantera Gas object, which ultimately is simply a representation of a gas thermochemical state, into another representation that would be more useful (and easier to manipulate with Python).

I will start with composition, which is the primary (and most complicated) feature of the gas. The goal is to convert a Cantera gas object to a Python dict type, with the keys being species names and the values being species mole fractions.

Converting Gas to Composition Dict

Function

def convert_gas_to_composition_dict( gas ):
    d = {}
    for sp in gas.speciesNames():
        d[sp] = gas.moleFraction(sp)
    return d

Usage

Note that although total mole amounts are specified with the set() call, these amounts are normalized (because the convert_gas_to_composition_dict() method calls the moleFraction() method).

In [2]: g = GRI30()

In [3]: set(g,T=600.0,P=3*OneAtm,X="CH4:1.0, O2:0.5, N2:3.76")

In [4]: convert_gas_to_composition_dict
Out[4]: <function ocm.util.CanteraGasUtils.convert_gas_to_composition_dict>

In [5]: convert_gas_to_composition_dict(g)
Out[5]:
{'AR': 0.0,
 'C': 0.0,
...
 'CH3OH': 0.0,
 'CH4': 0.19011406844106463,
 'CN': 0.0,
 'CO2': 0.0,
...
 'HOCN': 0.0,
 'N': 0.0,
 'N2': 0.714828897338403,
 'N2O': 0.0,
 'NCO': 0.0
...
 'O': 0.0,
 'O2': 0.09505703422053231,
 'OH': 0.0}

Convert Gas to Composition String

This function can be used to convert a Cantera Gas object to a string representation like "CH4:1.0, O2:0.4" (etc.)

This prints 16 digits of precision because I was dealing with very small radical concentrations, and one of the limitations of this method is lack of flexibility. So, adjust as needed, if you don't need precision then discard the 16f in there.

def convert_gas_to_composition_string( gas ):
    """
    Converts a Cantera gas with a specified composition into
    a representative composition string of style "CH4:1.0, N2:4.0"
    """
    d = convert_gas_to_composition_dict(gas)
    s = convert_composition_dict_to_string(d)
    return s

def convert_composition_dict_to_string( d ):
    """
    Converts a composition dict of stype composition['CH4'] = 0.02,
    into a string like "CH4: 0.02"
    """
    X = ""
    for item in d:
        X += "%s:%.16f"%(item,d[item])
        X += ", "
    X=X[:-2]
    return X

Converting Composition Dicts and Composition Strings

Dict to String

def convert_composition_dict_to_string( d ):
    """
    Converts a composition dict of stype composition['CH4'] = 0.02,
    into a string like "CH4: 0.02"
    """
    X = ""
    for item in d:
        X += "%s:%.16f"%(item,d[item])
        X += ", "
    X=X[:-2]
    return X

String to Dict

def convert_composition_string_to_dict( X ):
    """
    Converts a composition string of style "CH4:0.02, N2:0.01, O2:0.45"
    into a dict, a la composition['CH4'] = 0.02
    """
    results = {}
    for sp in X.split(","):
        st = sp.strip()
        try:
            results[ st.split(":")[0].strip() ] = float( st.split(":")[1].strip() )
        except IndexError:
            # X is probably not a list
            # (why would we run split on a list?)
            # or is empty
            err = "ERROR: CanteraGasUtils: your X is probably specified incorrectly. Expected type list, got type "+type(X)
            raise Exception(err)
    return results


Mixing Gases

Finding Thermochemical State of Mixture

Mixing a gaggle of gases requires the ability to mix composition (made possible through the representation of gases as dicts) as well as temperature and pressure. The final pressure is based on partial pressures, which is based on the mole fractions. The temperature can be set directly, or it can be set indirectly through setting the mixture enthalpy, in which case Cantera iterates over mixture temperatures with a Newton solver until it finds a temperature corresponding to the user-specified enthalpy.

This method is illustrated in the function mixture_HPX. It often fails to converge.

A better and faster approximation is to treat heat capacities as constant over the range of gas temperatures specified (i.e., the change in heat capacity from the initial to the final gas temperatures is negligible). This reasonable as long as the temperature ranges aren't unrealistic (and this can be checked by evaluating the mixture enthalpies using the code.enthalpy_mole() method in conjunction with set(gas,T=Tinitial) and <gas>set(gas,T=Tfinal). One could even implement one's own numerical method in place of Cantera's Newton iterator.

This temperature method is illustrated in the function mixture_TPX.


def mixture_TPX( gases, Xs, verbosity=0):
    """
    Given a list of gases and their mole fractions,
    this returns the TPX info needed to create
    a Cantera Gas object that is 
    a mixture of these gases.
    """

    # --------------
    # X

    mixture_d = {}

    if verbosity > 0:
        print "Looping over all species of all gases to make a new mixture:"
    for gas,wx_i in zip(gases,Xs):
        for sp in gas.speciesNames():
            if verbosity > 0:
                if sp == 'CH3':
                    print "CH3 mole frac = %0.4g"%(gas.moleFraction(sp))
            if sp in mixture_d:
                mixture_d[sp] += wx_i * gas.moleFraction(sp)
            elif gas.moleFraction(sp) != 0.0: 
                mixture_d[sp] = wx_i * gas.moleFraction(sp)
            else:
                pass

    mixture_s = convert_composition_dict_to_string(mixture_d)
    if verbosity > 0:
        print "Mixing resulted in gas with X = "+mixture_s

    # --------------
    # T

    # Compute Tmix with molar heat capacities
    #
    # Define:
    # h_mix = C_pmix T_mix = ( sum_i n_i C_pi Ti )/( n_T )
    #
    # from which we get relationship:
    # T_mix = \sum_i [ (x_i C_pi )/( C_pmix ) ] T_i

    # first compute c_pmix
    cp_mix = 0
    for gas, wx_i in zip(gases,Xs):
        cp_mix += wx_i * gas.cp_mole()

    # next compute T_mix
    T_mix = 0
    for gas, wx_i in zip(gases,Xs):
        coeff = ( wx_i * gas.cp_mole() )/( cp_mix )
        T_mix += coeff * gas.temperature()

    # --------------
    # P

    press = 0.0
    for gas,wx_i in zip(gases,Xs):
        press += wx_i * gas.pressure() 

    # -------------------
    # Return TPX
    
    return T_mix, press, mixture_s



def mixture_HPX( gases, Xs ):
    """
    Given a mixture of gases and their mole fractions,
    this method returns the enthalpy, pressure, and 
    composition string needed to initialize 
    the mixture gas in Cantera.

    NOTE: The method of setting enthalpy 
    usually fails, b/c Cantera uses a Newton
    iterator to find the temperature that
    yields the specified enthalpy, and it 
    isn't very robust.
    Instead, approximate constant Cp's
    and find T_mix manually, as with the
    mixture_TPX() method.
    """

    # --------------
    # X

    mixture_d = {}

    for gas,wx_i in zip(gases,Xs):
        for sp in gas.speciesNames():
            if sp in mixture_d:
                mixture_d[sp] += wx_i * gas.moleFraction(sp)
            elif gas.moleFraction(sp) != 0.0: 
                mixture_d[sp] = wx_i * gas.moleFraction(sp)
            else:
                pass

    mixture_s = convert_composition_dict_to_string(mixture_d)

    # --------------
    # H

    # Compute Tmix with molar heat capacities
    #
    # Define:
    # h_mix = sum_i n_i h_i 
    #
    # where h is molar enthalpy

    # compute H_mix
    H_mix = 0
    for gas, wx_i in zip(gases,Xs):
        Hmix += wx_i * gas.enthalpy_mole()

    # --------------
    # P

    press = 0.0
    for gas,wx_i in zip(gases,Xs):
        press += wx_i * gas.pressure() 

    # -------------------
    # Return HPX
    
    return H_mix, press, mixture_s


Mixing the Gases: Providing a Master Function

In order to mix gases easily, the following wrapper function is defined:


def getGasMixture( gases,     # list of gases
                   Xs,        # list of moles/mole fractions of each gas
                   model_file=None, 
                   gas_phase_name=None ):
    """
    Convert a list of gases and mole fractions
    (or mass fractions, or volume fractions, 
    or partial pressures, eventually)
    into a single Gas mixture.
    Do this by manually computing 
    mass-averaged mixture enthalpy,
    mole-averaged pressure,
    and composition mixing via composition dict
    """

    if model_file == None or gas_phase_name == None:
        err = "ERROR: CanteraGasUtils: You must specify a model file (via params['model_file']) and gas phase name (via params['gas_phase_name']) for the gas mixture."
        raise Exception(err)

    # pass in a normalized Xs
    Xsarr = np.array(Xs)
    Xsnorm = Xsarr/sum(Xsarr) 

    # technically a little better, 
    # but setting HPX is expensive 
    # (Newton iterations over T)
    #try:
    #    H, P, X = mixture_HPX( gases, Xnorm )
    #except CanteraError:
    #    T, P, X = mixture_TPX( gases, Xsnorm )

    # this approximates constant Cp
    # (i.e, only works for limited T differences 
    #       among geses being mixed):
    T, P, X = mixture_TPX( gases, Xsnorm )
    gas = Cantera.importPhase(model_file,gas_phase_name)
    gas.set(T=T,P=P,X=X)

    return gas


How To Use These Functions

Copy and paste the block of all the functions, given in the section below. You can paste those into CanteraGasMixing.py (or whatever you want to name the file), and put that in your site-packages directory, or the directory in which you're working with Cantera.

These functions will then enable you to mix various Cantera Gases in various thermochemical states, all using a simple time-based ODE solver. I was able to use this functionality to build a simple Cantera-based mixing and transport model.

All The Functions In One File

def convert_gas_to_composition_dict( gas ):
    d = {}
    for sp in gas.speciesNames():
        d[sp] = gas.moleFraction(sp)
    return d


def convert_gas_to_composition_string( gas ):
    """
    Converts a Cantera gas with a specified composition into
    a representative composition string of style "CH4:1.0, N2:4.0"
    """
    d = convert_gas_to_composition_dict(gas)
    s = convert_composition_dict_to_string(d)
    return s


def convert_composition_dict_to_string( d ):
    """
    Converts a composition dict of stype composition['CH4'] = 0.02,
    into a string like "CH4: 0.02"
    """
    X = ""
    for item in d:
        X += "%s:%.16f"%(item,d[item])
        X += ", "
    X=X[:-2]
    return X



def convert_composition_string_to_dict( X ):
    """
    Converts a composition string of style "CH4:0.02, N2:0.01, O2:0.45"
    into a dict, a la composition['CH4'] = 0.02
    """
    results = {}
    for sp in X.split(","):
        st = sp.strip()
        try:
            results[ st.split(":")[0].strip() ] = float( st.split(":")[1].strip() )
        except IndexError:
            # X is probably not a list
            # (why would we run split on a list?)
            # or is empty
            err = "ERROR: CanteraGasUtils: your X is probably specified incorrectly. Expected type list, got type "+type(X)
            raise Exception(err)
    return results



# =======================================
# Gas mixing:

def getGasMixture( gases,     # list of gases
                   Xs,        # list of moles/mole fractions of each gas
                   model_file=None, 
                   gas_phase_name=None ):
    """
    Convert a list of gases and mole fractions
    (or mass fractions, or volume fractions, 
    or partial pressures, eventually)
    into a single Gas mixture.
    Do this by manually computing 
    mass-averaged mixture enthalpy,
    mole-averaged pressure,
    and composition mixing via composition dict
    """

    if model_file == None or gas_phase_name == None:
        err = "ERROR: CanteraGasUtils: You must specify a model file (via params['model_file']) and gas phase name (via params['gas_phase_name']) for the gas mixture."
        raise Exception(err)

    # pass in a normalized Xs
    Xsarr = np.array(Xs)
    Xsnorm = Xsarr/sum(Xsarr) 

    # technically a little better, 
    # but setting HPX is expensive 
    # (Newton iterations over T)
    #try:
    #    H, P, X = mixture_HPX( gases, Xnorm )
    #except CanteraError:
    #    T, P, X = mixture_TPX( gases, Xsnorm )

    # this approximates constant Cp
    # (i.e, only works for limited T differences 
    #       among geses being mixed):
    T, P, X = mixture_TPX( gases, Xsnorm )
    gas = Cantera.importPhase(model_file,gas_phase_name)
    gas.set(T=T,P=P,X=X)

    return gas


def mixture_TPX( gases, Xs, verbosity=0):
    """
    Given a list of gases and their mole fractions,
    this returns the TPX info needed to create
    a Cantera Gas object that is 
    a mixture of these gases.
    """

    # --------------
    # X

    mixture_d = {}

    if verbosity > 0:
        print "Looping over all species of all gases to make a new mixture:"
    for gas,wx_i in zip(gases,Xs):
        for sp in gas.speciesNames():
            if verbosity > 0:
                if sp == 'CH3':
                    print "CH3 mole frac = %0.4g"%(gas.moleFraction(sp))
            if sp in mixture_d:
                mixture_d[sp] += wx_i * gas.moleFraction(sp)
            elif gas.moleFraction(sp) != 0.0: 
                mixture_d[sp] = wx_i * gas.moleFraction(sp)
            else:
                pass

    mixture_s = convert_composition_dict_to_string(mixture_d)
    if verbosity > 0:
        print "Mixing resulted in gas with X = "+mixture_s

    # --------------
    # T

    # Compute Tmix with molar heat capacities
    #
    # Define:
    # h_mix = C_pmix T_mix = ( sum_i n_i C_pi Ti )/( n_T )
    #
    # from which we get relationship:
    # T_mix = \sum_i [ (x_i C_pi )/( C_pmix ) ] T_i

    # first compute c_pmix
    cp_mix = 0
    for gas, wx_i in zip(gases,Xs):
        cp_mix += wx_i * gas.cp_mole()

    # next compute T_mix
    T_mix = 0
    for gas, wx_i in zip(gases,Xs):
        coeff = ( wx_i * gas.cp_mole() )/( cp_mix )
        T_mix += coeff * gas.temperature()

    # --------------
    # P

    press = 0.0
    for gas,wx_i in zip(gases,Xs):
        press += wx_i * gas.pressure() 

    # -------------------
    # Return TPX
    
    return T_mix, press, mixture_s


def mixture_HPX( gases, Xs ):
    """
    Given a mixture of gases and their mole fractions,
    this method returns the enthalpy, pressure, and 
    composition string needed to initialize 
    the mixture gas in Cantera.

    NOTE: The method of setting enthalpy 
    usually fails, b/c Cantera uses a Newton
    iterator to find the temperature that
    yields the specified enthalpy, and it 
    isn't very robust.
    Instead, approximate constant Cp's
    and find T_mix manually, as with the
    mixture_TPX() method above.
    """

    # --------------
    # X

    mixture_d = {}

    for gas,wx_i in zip(gases,Xs):
        for sp in gas.speciesNames():
            if sp in mixture_d:
                mixture_d[sp] += wx_i * gas.moleFraction(sp)
            elif gas.moleFraction(sp) != 0.0: 
                mixture_d[sp] = wx_i * gas.moleFraction(sp)
            else:
                pass

    mixture_s = convert_composition_dict_to_string(mixture_d)

    # --------------
    # H

    # Compute Tmix with molar heat capacities
    #
    # Define:
    # h_mix = sum_i n_i h_i 
    #
    # where h is molar enthalpy

    # compute H_mix
    H_mix = 0
    for gas, wx_i in zip(gases,Xs):
        Hmix += wx_i * gas.enthalpy_mole()

    # --------------
    # P

    press = 0.0
    for gas,wx_i in zip(gases,Xs):
        press += wx_i * gas.pressure() 

    # -------------------
    # Return HPX
    
    return H_mix, press, mixture_s


Flags