Feature request: Add function to estimate Hc based on molecular formula

[yalinli2]

Is your feature request related to a problem? Please describe. I was using the boilerturbogenerator unit, and I realized it was pulling Hc of chemicals to calculate the amount of steam/electricity produced. But many chemicals in the database do not have a default Hc.

Describe the solution you'd like I'm wondering if we can add an option to implement functions to estimate Hc. In the past I've used Dulong's formula (empirical correlation initially developed to estimate energy value of coal, I used that to calculate Hc of biomass, which is a common approach in literature) to estimate Hc based on the molecular formula.

Additional context I'll see if I can get this work. I think it'll be helpful to add molecular formula an attribute for of the chemicals, this can make it possible to automatically populate ParallelReaction for combustion that can be used in boilerturbogenerator, molecular formula should be easily extracted from InChI.

[yalinli2]

Here are the functions I made to do what I mentioned above:

To retrieve molecular formula from database: It's so darn hard parsing something like CHN12O12SCa10Fe, I wanted to make it capable of reading the floats but gave up, spent like 5 hours coming up with very verbose codes

First you need import re

def parse_formula(chemical):
    if not chemical.InChI:
        return chemical.ID + ' does not have molecular formula'
    split_InChI = chemical.InChI.split('/')
    formula = split_InChI[0]
    list_formula = list(formula)
    if '.' in list_formula:
        return 'Molecular formula of ' + chemical.ID + ' has float, cannot be parsed'

    formula_dict = {}
    complicated_elements = []
    elements = re.findall(r'\D+', formula)
    counts = re.findall(r'\d+', formula)
    # In case that formula ends with elements
    if len(counts) < len(elements): counts.append(1)

    for element, count in zip(elements, counts):
        formula_dict.update({element: int(count)})

    # In case that the formula has neighbouring elements with one of more of their number(s) being one
    for i in elements:
        if len(list(i))==2 and list(i)[-1].islower(): pass
        elif len(list(i))>1:
            complicated_elements.append(i)
            complicated_element_number = formula_dict[i]
            atom = ''
            atom_count = 0
            atom_dict = {}

            for x in i:
                if x.isupper():
                    if atom != '':
                        atom_dict[atom] = 1
                        atom = ''
                    atom = x

                elif x.islower():
                    atom += x

                else:
                    atom_count = str(atom_count) + str(x)
                    atom_dict[atom] = atom_count
                    atom_count = 0
                    atom = ''

            if atom != '': atom_dict[atom] = 1

            atom_dict[list(atom_dict.keys())[-1]] = complicated_element_number

            formula_dict.update(atom_dict)

    for i in complicated_elements:
        formula_dict.pop(i)

    return formula_dict

[yalinli2]

Set molecular weight based on molecular formula First you need a dictionary of molecular weight of all elements:

MW_of_elements = {'H': 1.00794, 'He': 4.002602, 'Li': 6.941, 'Be': 9.012182, 'B': 10.811, 'C': 12.0107, 'N': 14.0067,
              'O': 15.9994, 'F': 18.9984032, 'Ne': 20.1797, 'Na': 22.98976928, 'Mg': 24.305, 'Al': 26.9815386,
              'Si': 28.0855, 'P': 30.973762, 'S': 32.065, 'Cl': 35.453, 'Ar': 39.948, 'K': 39.0983, 'Ca': 40.078,
              'Sc': 44.955912, 'Ti': 47.867, 'V': 50.9415, 'Cr': 51.9961, 'Mn': 54.938045,
              'Fe': 55.845, 'Co': 58.933195, 'Ni': 58.6934, 'Cu': 63.546, 'Zn': 65.409, 'Ga': 69.723, 'Ge': 72.64,
              'As': 74.9216, 'Se': 78.96, 'Br': 79.904, 'Kr': 83.798, 'Rb': 85.4678, 'Sr': 87.62, 'Y': 88.90585,
              'Zr': 91.224, 'Nb': 92.90638, 'Mo': 95.94, 'Tc': 98.9063, 'Ru': 101.07, 'Rh': 102.9055, 'Pd': 106.42,
              'Ag': 107.8682, 'Cd': 112.411, 'In': 114.818, 'Sn': 118.71, 'Sb': 121.760, 'Te': 127.6,
              'I': 126.90447, 'Xe': 131.293, 'Cs': 132.9054519, 'Ba': 137.327, 'La': 138.90547, 'Ce': 140.116,
              'Pr': 140.90465, 'Nd': 144.242, 'Pm': 146.9151, 'Sm': 150.36, 'Eu': 151.964, 'Gd': 157.25,
              'Tb': 158.92535, 'Dy': 162.5, 'Ho': 164.93032, 'Er': 167.259, 'Tm': 168.93421, 'Yb': 173.04,
              'Lu': 174.967, 'Hf': 178.49, 'Ta': 180.9479, 'W': 183.84, 'Re': 186.207, 'Os': 190.23, 'Ir': 192.217,
              'Pt': 195.084, 'Au': 196.966569, 'Hg': 200.59, 'Tl': 204.3833, 'Pb': 207.2, 'Bi': 208.9804,
              'Po': 208.9824, 'At': 209.9871, 'Rn': 222.0176, 'Fr': 223.0197, 'Ra': 226.0254, 'Ac': 227.0278,
              'Th': 232.03806, 'Pa': 231.03588, 'U': 238.02891, 'Np': 237.0482, 'Pu': 244.0642, 'Am': 243.0614,
              'Cm': 247.0703, 'Bk': 247.0703, 'Cf': 251.0796, 'Es': 252.0829, 'Fm': 257.0951, 'Md': 258.0951,
              'No': 259.1009, 'Lr': 262, 'Rf': 267, 'Db': 268, 'Sg': 271, 'Bh': 270, 'Hs': 269, 'Mt': 278,
              'Ds': 281, 'Rg': 281, 'Cn': 285, 'Nh': 284, 'Fl': 289, 'Mc': 289, 'Lv': 292, 'Ts': 294, 'Og': 294}

Then it's easy

def set_MW(chemical):
    if chemical.MW: 
        return chemical.ID + ' already has MW'
        pass
    if not chemical.InChI:
        return chemical.ID + ' does not have molecular formula to set MW'
        pass

    formula_dict = parse_formula(chemical)
    MW = 0
    for element in formula_dict:
        MW += formula_dict[element] * MW_of_elements[element]

    chemical.MW = MW

[yalinli2]

Set heat of combustion based for organics

Note that I edited this as the previous one used Dulong's formula incorrectly, I checked both estimation methods for ethanol and octane against data in the database, they give like 10±5% errors, I consider them acceptable.

def set_organics_Hc(chemical):
    # Set heat of combustion (J/mol) based on existing properties,
    # this is only applicable for organics containing only C/H/O/N/S,
    # assume molecular formula of the chemical is C(n1)H(n2)O(n3)N(n4)S(n5)
    # does not consider latent heat of water
    # (a) If chemical has molecular formula and heat of formation Hf,
    #     then calculate Hc based on stoichiometric combustion as:
    #     C(n1)H(n2)O(n3)N(n4)S(n5) + (n1+n2/4+n4+n5-n3/2) O2 
    #                               -> n1 CO2 +  n2/2 H2O + n4 NO2 + n5 SO2
    # (b) If chemical has molecular formula and heat of formation but not Hf, 
    #     then calculate Hc based on Dulong's equation as:
    #     Hc (J/g) = 338*C% + 1428(H%-O%/8)+ 95*S%
    #     Ref:  Brown et al., Energy Fuels 2010, 24 (6), 3639–3646.
    #     Note: this technically only good for <10% O_content

    if chemical.Hc: 
        return chemical.ID + ' already has Hc'
        pass
    if not chemical.InChI: 
        return chemical.ID + ' does not have molecular formula to set Hc'
        pass
    formula_dict = parse_formula(chemical)
    combustable_elements = ['C', 'H', 'O', 'N', 'S']
    non_combustable_elements = []

    for i in formula_dict.keys():
        if not i in combustable_elements:
            non_combustable_elements.append(i)
        else: pass

    if len(non_combustable_elements) > 0:
        return chemical.ID + ' has inorganic elements, Hc cannot be set'
    else:
        if not chemical.MW: set_MW(chemical)
        try: 
            n1 = formula_dict['C']
            C_content = n1*12.0107/chemical.MW*100
        except: n1 = C_content = 0

        try: 
            n2 = formula_dict['H']
            H_content = n2*1.00794/chemical.MW*100
        except: n2 = H_content = 0

        try: 
            n3 = formula_dict['O']
            O_content = n3*15.9994/chemical.MW*100
        except: n3 = O_content = 0

        try: 
            n4 = formula_dict['N']
            N_content = n4*14.0067/chemical.MW*100
        except: n4 = N_content = 0

        try: 
            n5 = formula_dict['S']
            S_content = n5*32.065/chemical.MW*100
        except: n5 = S_content = 0

        if chemical.Hf:
            chemical.Hc = (n1*-393530.0  # Hf of CO2
                           +n2/2*-241820.0  # Hf of H2O
                           +n4*33100.0  # Hf of NO2
                           +n5*-296850.0) - chemical.Hf # Hf of SO2 
        else:
            if not chemical.MW: set_MW(chemical)
            chemical.Hc = - (338*C_content + 1428*(H_content-O_content/8)+ 95*S_content)*chemical.MW

[yoelcortes]

Thanks for request and the code! Thermosteam already calculates heat of combustion based on the formula/atoms (see function below). ~Many chemicals do not have the formula but do have the InChI key, so we could parse it and get the formula and atoms (just as you suggested). We can also use Dulong's equation too to set Hc for those chemicals that do not have it. As for estimating the MW, we can include that as part of a user method to set formula/atoms when creating a new chemical.~

EDIT: Every chemical available in Thermosteam has the formula, so no need to parse to get the formla. The atom counts are found using Thermosteam's built-in parser. The atom counts and the formula are used internally when initializing a chemical, but are not saved as attributes. The reason Hc is not computed for some chemicals is because Hf is missing. Another issue is that Hc is vague as it could mean the either lower or higher heating value. So hear are the main changes that I could make:

• Save the formula as a chemical property formula. Add a get_atom_counts method that returns a dictionary of atom counts.
• Deprecate the Hc attribute in favor of LHV and HHV
• Add a built-in function interface to estimate Hf, LHV and HHV given the atoms, and optionally the heat of formation, and a desired method (e.g. Dulong, Tillman, Corbitt, Demirbas, etc.). Non-combustibles will have a heat of combustion of 0.
• Make Hf, LHV, and HHV required chemical properties that must be present before compiling a Chemicals object.
• Define combustion reactions when initializing chemical as a string. To maintain the mass balance, non-combustible atoms will form a chemical called "ash".
• When a Chemicals object is compiled, it will create a ParallelReaction object for combustion using the combustion reaction strings. A Stream object may then use the ParallelReaction object for things like emissions.

I will go ahead and include these proposed features :-)

Just FYI, this is currently what has been implemented:

def estimate_Hc(atoms, Hf=None, HfH2O=-285825, HfCO2=-393474,
                HfSO2=-296800, HfBr2=30880, HfI2=62417, HfHCl=-92173,
                HfHF=-272711, HfP4O10=-3009940, HfO2=0, HfN2=0):
    '''
    Calculates the heat of combustion, in J/mol.
    Value non-hydrocarbons is not correct, but still calculable.

    Parameters
    ----------
    atoms : dict
        Dictionary of atoms and their counts, []
    Hf : float
        Heat of formation of given chemical, [J/mol]
    HfH2O : float, optional
        Heat of formation of water, [J/mol]
    HfCO2 : float, optional
        Heat of formation of carbon dioxide, [J/mol]
    HfSO2 : float, optional
        Heat of formation of sulfur dioxide, [J/mol]
    HfBr2 : float, optional
        Heat of formation of bromine, [J/mol]
    HfI2 : float, optional
        Heat of formation of iodine, [J/mol]
    HfHCl : float, optional
        Heat of formation of chlorine, [J/mol]
    HfHF : float, optional
        Heat of formation of hydrogen fluoride, [J/mol]
    HfP4O10 : float, optional
        Heat of formation of phosphorus pentoxide, [J/mol]
    HfO2 : float, optional
        Heat of formation of oxygen, [J/mol]
    HfN2 : float, optional
        Heat of formation of nitrogen, [J/mol]

    Returns
    -------
    Hc : float
        Heat of combustion of chemical, [J/mol]

    Notes
    -----
    Default heats of formation for chemicals are at 298 K, 1 atm.

    Examples
    --------
    Liquid methanol burning

    >>> Hcombustion({'H': 4, 'C': 1, 'O': 1}, Hf=-239100)
    -726024.0
    '''
    if not Hf or not atoms: return None
    if 'C' not in atoms: return None
    nC = atoms['C']
    if not nC: return None
    nH = atoms.get('H', 0)
    nN = atoms.get('N', 0)
    nO = atoms.get('O', 0)
    nS = atoms.get('S', 0)
    nBr = atoms.get('Br', 0)
    nI = atoms.get('I', 0)
    nCl = atoms.get('Cl', 0)
    nF = atoms.get('F', 0)
    nP = atoms.get('P', 0)
    nO2_req = nC + nS + nH/4. + 5*nP/4. - (nCl + nF)/4. - nO/2.
    nCO2 = nC
    nBr2 = nBr/2.
    nI2 = nI/2.
    nHCl = nCl
    nHF = nF
    nSO2 = nS
    nN2 = nN/2.
    nP4O10 = nP/4.
    nH2O = (nH - nCl - nF)/2.
    Hc = (nBr2*HfBr2 + nI2*HfI2) + (nHCl*HfHCl + nHF*HfHF) + nSO2*HfSO2 + \
         nN2*HfN2 + nP4O10*HfP4O10 + nH2O*HfH2O - nO2_req*HfO2 + nCO2*HfCO2 - Hf
    return Hc

NOTE: The heats of combustion should be negative if releasing heat. There is a bug in the BoilerTurbogenerater unit operation and the biorefineries that assumes in should be positive. I will also fix this up in the next update.

[yalinli2]

Sounds great! I'll look forward to the updates.

I started doing this because the current BoilerTurbogenerater gives negativemakeup_water flow when Hc in the ins is smaller than heat needed for to evaporate all water, then I realized a part of that was because I forgot to define Hc for many chemicals. So I was trying to make it easier to estimate Hc, maybe let the BoilerTurbogenerater estimate Hc in case the user forgot to define it. I also want to improve the _design function for the unit because now it gives invalid (negative) cost in this situation.

Another thing I want to improve is how the BoilerTurbogenerater unit makes outs. Right now the solid outs just copies the solid ins, so I want the unit to automatically create combustion reactions. The idea of including the combustion reaction into the reaction string is super cool and will make this very easy! For now I'll try to include this in the unit and you can look at the codes and evaluate whether it's better to include this in the Stream object or in the unit.

I'm making a new BoilerTurbogenerater unit for the orgacids package, so I'll let you know when I get this done (my time scale of code development is at least 10X larger than yours 😂).

Also ah correct on the Hc - thanks for the notes!!!

PS, I'll leave this open until the new features are added and I finish making the new BoilerTurbogenerater unit.

[yalinli2]

Along this line, another thing I think will be helpful to the users is for units to have:

1. A design_requirement attribute (or a section in the unit help information) which includes what properties of the streams are needed in unit simulation (like Hc in the case of the BoilerTurbogenerator.
2. Automatically calls functions to calculate the needed properties if users forget to define and it's possible to do so.
3. In raising errors, explicitly tell the user that properties prevent simulation. This can be dumb but helpful for non-chemical engineers (like me!). For example, it took me a while to figure out why the Distillation unit kept telling me I had intermediate chemicals between light and heavy key, it'll be much more intuitive if the error can have something in parentheses like Tb of (IntermediateChemical) is between Tb of (LightKey) and Tb of (HeavyKey). I know this one can be a lot to ask, but I think will be very helpful 😝

[yoelcortes]

I believe Sarang can help you out with the negative cost for the boiler-turbogenerator. It has to do with wrong values of heats of combustion. I helped him debug it.

As for the suggestions:

1. The notes section of a docstring should include details on how the unit operation works (e.g., the need for Hc in the BoilerTurbogenerator). However, many times an implementation is only temporary just to get the job done, so it is not well documented. In these cases, best is look at the source code and feel free to ask any questions :)
2. That's the goal!
3. Thanks for pointing these out. Specific descriptions for each possible error is the right mindset.

Note: Sometimes the code might error for an unforeseen reason and spit out an error that doesn't guide you well. For these kinds of bugs, you could also try going to the source code and start printing stuff (in addition to posting the issue). Also, as you may have noticed already, I made thermosteam prevent users from compiling chemicals when key properties are missing. I'll include Tb, LHV, HHV, as a key properties in the next version. That way missing property issues are cached early on.

[yalinli2]

Cool! All makes sense and thanks for incorporating those features!

[yalinli2]

NOTE: The heats of combustion should be negative if releasing heat. There is a bug in the BoilerTurbogenerater unit operation and the biorefineries that assumes in should be positive. I will also fix this up in the next update.

If you are going to update biorefineries, then I just noticed another bug: in the biorefineries.cornstover.chemicals module, you added Hc to a bunch of chemicals using Hc of lignin and cellulosic, the numbers you used were 21000 and 17000 J/g, respectively, but the unit of Chemical.Hc is J/mol, so you need to add the conversion.

[yoelcortes]

Thanks for noting this, Hc_cellulosic was only used once for the "Cellulose" chemical, while Hc_lignin was not used. However, all the heat of combustions was replaced in "system.py" using the combustion formulas in Humbird's report.

That said, I'll try to simplify/refactor the code. I started working on the new features, it should take about a week to be integrated and tested. :)

[yalinli2]

Thanks for doing this!

BTW, I just finished rewriting the chemicals.py module for orgacids. I'll certainly update the script after you finish integrating the new features, many functions won't be needed then. But I put it here for your references (compressed it, GitHub won't let me upload ".py"), let me know any comments! No hurries though, I probably won't be touching this in at least a week, thanks!

chemicals.py.zip

The priority hierarchy I used in defining the chemicals:

1. If information available in the database, then follow the database
2. If no information in database, then follow literature (mostly Humbird et al.)
3. Use functions to estimate missing properties if possible

[yoelcortes]

The new features are out! For now you can look at the thermosteam.functors.combution module if you would like details on how LHV and HHV are calculated. I am planning to move some modules around, so I won't include the module in readthedocs just yet.

[yalinli2]

COOL!!! I'm fine without the docs now, I'm getting used to looking at the codes now 😄