In [1]:
import pandas as pd
In [2]:
# Pyrolysis reactions in the Debiagi 2018 reaction scheme
# THE SPECIES "G{COH2}STIFF" AND "G{COH2}LOOSE" ARE HERE NAMED "COH2S" AND "CH2OS"
rxn1 = {"reactants":[ 1, "CELL"], "products": [1, "CELLA"]}
rxn2 = {"reactants":[1, "CELLA"], "products": [0.4, "CH2OHCHO",
0.03, "CHOCHO",
0.17, "CH3CHO",
0.25, "C6H6O3",
0.35, "C2H5CHO",
0.2, "CH3OH",
0.15, "CH2O",
0.49, "CO",
0.05, "COS",
0.43, "CO2",
0.13, "H2",
0.93, "H2O",
0.05, "CH2OS",
0.02, "HCOOH",
0.05, "CH2OHCH2CHO",
0.05, "CH4",
0.1, "H2S",
0.66, "CHAR"]}
rxn3 = {"reactants":[ 1, "CELLA"], "products": [1, "C6H10O5"]}
rxn4 = {"reactants":[ 1, "CELL"], "products": [4.45, "H2O",
5.45, "CHAR",
0.12, "COH2S",
0.18, "CH2OS",
0.25, "COS",
0.125, "H2S",
0.125, "H2"]}
rxn5 = {"reactants":[ 1, "GMSW"], "products": [0.7, "HCE1", 0.3, "HCE2"]}
rxn6 = {"reactants":[ 1, "XYHW"], "products": [0.35, "HCE1", 0.65, "HCE2"]}
rxn7 = {"reactants":[ 1, "XYGR"], "products": [0.12, "HCE1", 0.88, "HCE2"]}
rxn8 = {"reactants":[ 1, "HCE1"], "products": [0.25, "C5H8O4",
0.25, "C6H10O5",
0.16, "FURFURAL",
0.13, "C6H6O3",
0.09, "CO2",
0.1, "CH4",
0.54, "H2O",
0.06, "CH2OHCH2CHO",
0.1, "CHOCHO",
0.02, "H2",
0.1, "CHAR"]}
rxn9 = {"reactants":[ 1, "HCE1"], "products": [0.4, "H2O",
0.39, "CO2",
0.05, "HCOOH",
0.49, "CO",
0.01, "COS",
0.51, "CO2S",
0.05, "H2S",
0.4, "CH2O",
0.43, "CH2OS",
0.3, "CH4",
0.325, "CH4S",
0.1, "C2H4",
0.075, "C2H4S",
0.975, "CHAR",
0.37, "COH2S",
0.1, "H2",
0.2, "C2H6S"]}
rxn10 = {"reactants":[ 1, "HCE2"], "products": [0.3, "CO",
0.5125, "CO2",
0.1895, "CH4",
0.5505, "H2",
0.056, "H2O",
0.049, "C2H5OH",
0.035, "CH2OHCHO",
0.105, "CH3CO2H",
0.0175, "HCOOH",
0.145, "FURFURAL",
0.05, "CH4S",
0.105, "CH3OH",
0.1, "C2H4S",
0.45, "CO2S",
0.18, "CH2OS",
0.7125, "CHAR",
0.21, "H2S",
0.78, "COH2S",
0.2, "C2H6S"]}
rxn11 = {"reactants":[ 1, "LIGH"], "products": [1, "LIGOH",
0.5, "C2H5CHO",
0.4, "C2H4",
0.2, "CH2OHCHO",
0.1, "CO",
0.1, "C2H6"]}
rxn12 = {"reactants":[ 1, "LIGO"], "products": [1, "LIGOH", 1, "CO2"]}
rxn13 = {"reactants":[ 1, "LIGC"], "products": [0.35, "LIGCC",
0.1, "VANILLIN",
0.1, "C6H5OCH3",
0.27, "C2H4",
1, "H2O",
0.17, "CH2OS",
0.4, "COH2S",
0.22, "CH2O",
0.21, "CO",
0.1, "CO2",
0.36, "CH4S",
5.85, "CHAR",
0.2, "C2H6S",
0.1, "H2S"]}
rxn14 = {"reactants":[ 1, "LIGCC"], "products": [0.25, "VANILLIN",
0.15, "CRESOL",
0.15, "C6H5OCH3",
0.35, "CH2OHCHO",
0.7, "H2O",
0.45, "CH4",
0.3, "C2H4",
0.7, "H2",
1.15, "CO",
0.4, "COS",
6.80, "CHAR",
0.4, "C2H6"]}
rxn15 = {"reactants":[ 1, "LIGOH"], "products": [0.9, "LIG",
1, "H2O",
0.1, "CH4",
0.6, "CH3OH",
0.3, "CH3OHS",
0.05, "CO2",
0.65, "CO",
0.6, "COS",
0.05, "HCOOH",
0.45, "CH2OS",
0.4, "COH2S",
0.25, "CH4S",
0.1, "C2H4S",
0.15, "C2H6S",
4.25, "CHAR",
0.025, "C24H28O4",
0.1, "C2H3CHO"]}
rxn16 = {"reactants":[ 1, "LIG"], "products": [1, "VANILLIN",
0.1, "C6H5OCH3",
0.5, "C2H4",
0.6, "CO",
0.3, "CH3CHO",
0.1, "CHAR"]}
rxn17 = {"reactants":[ 1, "LIG"], "products": [ 0.6, "H2O",
0.3, "CO",
0.1, "CO2",
0.2, "CH4",
0.4, "CH2O",
0.2, "COS",
0.4, "CH4S",
0.5, "C2H4S",
0.4, "CH3OHS",
1.25, "CH2OS",
0.65, "COH2S",
6.1, "CHAR",
0.1, "H2S"]}
rxn18 = {"reactants":[ 1, "LIG"], "products": [0.6, "H2O",
2.6, "CO",
0.6, "CH4",
0.4, "CH2O",
0.75, "C2H4",
0.4, "CH3OH",
4.5, "CHAR",
0.5, "C2H6"]}
rxn19 = {"reactants":[ 1, "TGL"], "products": [1, "C2H3CHO",
2.5, "MLINO",
0.5, "U2ME12"]}
rxn20 = {"reactants":[ 1, "TANN"], "products": [0.85, "C6H5OH",
0.15, "C6H5OHS",
1, "COS",
1, "H2O",
1, "ITANN"]}
rxn21 = {"reactants":[ 1, "ITANN"], "products": [5, "CHAR",
2, "CO",
1, "H2O",
0.55, "CH2OS",
0.45, "COH2S"]}
rxn22 = {"reactants":[ 1, "CO2S"], "products": [1, "CO2"]}
rxn23 = {"reactants":[ 1, "COS"], "products": [1, "CO"]}
rxn24 = {"reactants":[ 1, "CH3OHS"], "products": [1, "CH3OH"]}
rxn25 = {"reactants":[ 1, "CH2OS"], "products": [0.2, "CO",
0.2, "H2",
0.8, "H2O",
0.8, "CHAR"]}
rxn26 = {"reactants":[ 1, "C2H6S"], "products": [1, "C2H6"]}
rxn27 = {"reactants":[ 1, "CH4S"], "products": [1, "CH4"]}
rxn28 = {"reactants":[ 1, "C2H4S"], "products": [1, "C2H4"]}
rxn29 = {"reactants":[ 1, "C6H5OHS"], "products": [1, "C6H5OH"]}
rxn30 = {"reactants":[ 1, "COH2S"], "products": [0.8, "CO",
0.8, "H2",
0.2, "H2O",
0.8, "CHAR"]}
rxn31 = {"reactants":[ 1, "H2S"], "products": [1, "H2"]}
rxn32 = {"reactants":[ 1, "ACQUA"], "products": [1, "H2O"]}
rxn_list = [rxn1, rxn2, rxn3, rxn4, rxn5, rxn6, rxn7, rxn8, rxn9, rxn10,
rxn11, rxn12, rxn13, rxn14, rxn15, rxn16, rxn17, rxn18, rxn19,
rxn20, rxn21, rxn22, rxn23, rxn24, rxn25, rxn26, rxn27, rxn28,
rxn29, rxn30, rxn31, rxn32]
In [3]:
def nasa_polynomial(species, thermo_df, R):
"""
Calculates heat of formation of a species from nasa polynomial
coefficients in the `thermo_df` dataframe.
Input arguments:
`species` should be a string containing the species name
`species` must be present in `thermo_df`
`R` is the ideal gas constant. Its units determine
the units of the heat of formation.
Returns the heat of formation of `species`.
"""
row = thermo_df.loc[thermo_df['Name'] == species]
if not all('' == s or s.isspace() for s in row["Hº(298.15)/R"]):
species_has_heat_data = True
else:
species_has_heat_data = False
if species_has_heat_data and T == 298.15:
heat = float(row.iloc[0]["Hº(298.15)/R"]) * R
else:
if T >= 1000:
a1 = row.iloc[0]["a1 T>= 1000 K"]
a2 = row.iloc[0]["a2 T>= 1000 K"]
a3 = row.iloc[0]["a3 T>= 1000 K"]
a4 = row.iloc[0]["a4 T>= 1000 K"]
a5 = row.iloc[0]["a5 T>= 1000 K"]
b1 = row.iloc[0]["b1 T>= 1000 K"]
b2 = row.iloc[0]["b2 T>= 1000 K"]
else:
a1 = row.iloc[0]["a1 T<= 1000 K"]
a2 = row.iloc[0]["a2 T<= 1000 K"]
a3 = row.iloc[0]["a3 T<= 1000 K"]
a4 = row.iloc[0]["a4 T<= 1000 K"]
a5 = row.iloc[0]["a5 T<= 1000 K"]
b1 = row.iloc[0]["b1 T<= 1000 K"]
b2 = row.iloc[0]["b2 T<= 1000 K"]
# Heat of reaction polynomial from NASA
p1 = float(a1)
p2 = float(a2)*(T/2)
p3 = float(a3)*((T**2)/3)
p4 = float(a4)*((T**3)/4)
p5 = float(a5)*((T**4)/5)
p6 = float(b1)/T
heat = (R*T) * (p1 + p2 + p3 + p4 + p5 + p6)
return heat
def h_of_rxn(rxn, thermo_df, T, R):
"""
Use with `nasa_polynomials()` function to calculate heat
of reaction.
Input arguments:
`rxn` is the reaction of interest, from the Reactions cell.
for example, rxn22
`thermo_df` is the dataframe containing nasa polynomial coefficients
for each species.
`T` is the temperature for use in the calculation, its units should be
consistent with the units for `R`
`R` is the ideal gas constant. Its units determine the units of the
heat of reaction that is returned from this function.
Returns a list containing the reaction, sum of the heats of formations
for products, sum of heats of formations for reactants, and the heat
of reaction calculated.
"""
reactant_stoich_list = []
reactant_species_list = []
product_stoich_list = []
product_species_list = []
reactant_heats = []
reactant_mw = []
product_heats = []
# Gather list of all stoichiometric coefficients and
# species names of reactants
for i in range(0, len(rxn["reactants"]), 2):
reactant_stoich = rxn["reactants"][i]
reactant_stoich_list.append(reactant_stoich)
reactant_species = rxn["reactants"][i+1]
reactant_species_list.append(reactant_species)
# Gather list of all stoichiometric coefficients and
# species names of products
for i in range(0, len(rxn["products"]), 2):
product_stoich = rxn["products"][i]
product_stoich_list.append(product_stoich)
product_species = rxn["products"][i+1]
product_species_list.append(product_species)
for i in range(len(reactant_species_list)):
np_reactant = nasa_polynomial(reactant_species_list[i], thermo_df, R)
reactant_heats.append(np_reactant)
for i in range(len(product_species_list)):
np_product = nasa_polynomial(product_species_list[i], thermo_df, R)
product_heats.append(np_product)
# Multiply heats by stoichiometric cofficients
for i in range(len(reactant_heats)):
reactant_heats[i] = reactant_heats[i] * reactant_stoich_list[i]
for i in range(len(product_heats)):
product_heats[i] = product_heats[i] * product_stoich_list[i]
total_heat_products = sum(product_heats)
total_heat_reactants = sum(reactant_heats)
heat_of_reaction = total_heat_products - total_heat_reactants
return [rxn, total_heat_products, total_heat_reactants, heat_of_reaction]
def convert_heat_g(rxn, thermo_df, heat):
"""
Converts heat from molar basis to mass basis
using molecular weight of reactant in g/mol.
For now, it assumes only one reactant.
Input arguments:
`rxn` is the reaction number of interest
`thermo_df` is the dataframe containing molecular weights
of each species
`heat` is the heat of reaction in energy/mol
Returns the heat of reaction on an energy/g basis
"""
# multiply heat in kJ/mol by molecular weight of reactant
# to convert to kJ/g
# assumes only one reactant for the time being
reactant_name = rxn["reactants"][1]
reactant_stoich = rxn["reactants"][0]
row = thermo_df.loc[thermo_df['Name'] == reactant_name]
mw = row.iloc[0]["MW (g/mol)"]
heat_g = (heat * reactant_stoich) / mw
return heat_g
In [4]:
T = 298.15 # Kelvin
thermo_file = "thermo.tdc"
# Column names from nasa-polynomial-doc.pdf p. 9
column_names =["Name",
"a1 T>= 1000 K",
"a2 T>= 1000 K",
"a3 T>= 1000 K",
"a4 T>= 1000 K",
"a5 T>= 1000 K",
"b1 T>= 1000 K",
"b2 T>= 1000 K",
"a1 T<= 1000 K",
"a2 T<= 1000 K",
"a3 T<= 1000 K",
"a4 T<= 1000 K",
"a5 T<= 1000 K",
"b1 T<= 1000 K",
"b2 T<= 1000 K",
"Hº(298.15)/R"]
# Initialize dataframe with column names
thermo_df = pd.DataFrame(columns=column_names)
# Read in thermo.tdc and populate dataframe with species names &
# nasa polynomial coefficients
with open(thermo_file, "r") as f:
# skip first 2 lines (header)
for i in range(2):
line = next(f)
for line in f:
if not line.isspace():
species_name = line[:11].rstrip()
species_name = species_name.lstrip()
try:
line = next(f)
# Column formatting from nasa-polynomial-doc.pdf p. 9
a1_gt = line[:15]
a2_gt = line[15:30]
a3_gt = line[30:45]
a4_gt = line[45:60]
a5_gt = line[60:75]
except Exception as e:
logger.error(
'Failed to get coefficients a1-a5 for' \
f'T > 1000 K: for {species_name}' + str(e))
try:
line = next(f)
b1_gt = line[:15]
b2_gt = line[15:30]
a1_lt = line[30:45]
a2_lt = line[45:60]
a3_lt = line[60:75]
except Exception as e:
logger.error(
'Failed to get coefficients b1-b2 for' \
'T > 1000 K or a1-a3 for T < 1000K: for' \
f'{species_name}' + str(e))
try:
line = next(f)
a4_lt = line[:15]
a5_lt = line[15:30]
b1_lt = line[30:45]
b2_lt = line[45:60]
# Not all entries have heat of reaction data
h_rxn_298 = line[60:79]
row_data = [species_name, a1_gt, a2_gt, a3_gt,
a4_gt, a5_gt, b1_gt, b2_gt, a1_lt,
a2_lt, a3_lt, a4_lt, a5_lt, b1_lt,
b2_lt, h_rxn_298]
species_df = pd.DataFrame([row_data], columns=column_names)
thermo_df = thermo_df.append(species_df, ignore_index=True)
except Exception as e:
logger.error(
'Failed to get coefficients a4-b2 for' \
f'T < 1000K: for {species_name}' + str(e))
f.close()
# Example row for species containing nasa polynomial coefficients
species_df
Out[4]:
In [5]:
# Load in molecular weight data
mw_df = pd.read_csv("./molecular-weights.csv", header=0)
mw_list = []
index_list = []
name_list = []
for index1, row1 in thermo_df.iterrows():
name1 = row1[0].rstrip()
name1 = name1.lstrip()
for index2, row2 in mw_df.iterrows():
name2 = row2[0].rstrip()
name2 = name2.lstrip()
if name1 == name2 and name1 not in name_list:
name_list.append(row1[0])
mw = pd.DataFrame([row2[1]], columns=["MW g/mol"])
index_list.append(index1)
mw_list.append(row2[1])
# Check that all species have mws (ash has NaN for now)
for index, row in thermo_df.iterrows():
if row[0] not in name_list:
print(f"Species {row[0]} is missing molecular weight data.")
# Add molecular weights for each species to dataframe
thermo_df["MW (g/mol)"] = mw_list
# Double check molecular weights for each species is correct
# by uncommenting the line below
#thermo_df[['Name','MW (g/mol)']]
In [6]:
# Write C_p equation with coefficients for each species to .txt file for use with COMSOL
# T < 1000 K
Cp_string_list = []
for index, row in thermo_df.iterrows():
name = row[0]
a1 = row[8]
a2 = row[9]
a3 = row[10]
a4 = row[11]
a5 = row[12]
Cp_string = f"Cp_{name} = R_const*({a1} + {a2}*T + {a3}*T^2 + {a4}*T^3 + {a5}*T^4)\n"
Cp_string_list.append(Cp_string)
with open("heat-capacities.txt", "w") as f:
for string in Cp_string_list:
f.write(string)
f.close()
In [7]:
# Write heats of formation for reaction products, reactant, and total heat of reaction
# in kJ/mol and kJ/g to heats-25.csv
R = 8.314/1000 #ideal gas constant in kJ/(mol*K)
# Delete or rename this file if it already exists
# to avoid appending to end
outfile = "heats-25C.csv"
with open(outfile, "a") as f:
i = 1
f.write("Reaction Number, H_products (kJ/mol), H_reactants (kJ/mol), deltaH_reaction (kJ/mol), deltaH_reaction (kJ/g)\n")
for reaction in rxn_list:
data_to_write = h_of_rxn(reaction, thermo_df, T, R)
rxn = data_to_write[0]
h_prod = data_to_write[1]
h_react = data_to_write[2]
h_rxn = data_to_write[3]
h_rxn_kJg = convert_heat_g(reaction, thermo_df, h_rxn)
f.write(f"Reaction {i},{h_prod},{h_react},{h_rxn}, {h_rxn_kJg}\n")
i +=1
f.close()
In [8]:
# Testing reaction 22 CO2S -> CO2
# 3 Test Scenarios:
# Compare values returned from function, calculations by manually querying dataframe, and
# calculations by hard coding coefficient values
# They should all be the same value
T = 298.15
R = 8.314/1000 #ideal gas constant in kJ/(mol*K)
# Scenario 1: Manually queru dataframe to calculate heat of reaction
reactant = "CO2S"
product = "CO2"
mw_r = 44.01
r_row = thermo_df.loc[thermo_df['Name'] == reactant]
p_row = thermo_df.loc[thermo_df['Name'] == product]
r_a1 = r_row.iloc[0]["a1 T<= 1000 K"]
r_a2 = r_row.iloc[0]["a2 T<= 1000 K"]
r_a3 = r_row.iloc[0]["a3 T<= 1000 K"]
r_a4 = r_row.iloc[0]["a4 T<= 1000 K"]
r_a5 = r_row.iloc[0]["a5 T<= 1000 K"]
r_b1 = r_row.iloc[0]["b1 T<= 1000 K"]
p_a1 = p_row.iloc[0]["a1 T<= 1000 K"]
p_a2 = p_row.iloc[0]["a2 T<= 1000 K"]
p_a3 = p_row.iloc[0]["a3 T<= 1000 K"]
p_a4 = p_row.iloc[0]["a4 T<= 1000 K"]
p_a5 = p_row.iloc[0]["a5 T<= 1000 K"]
p_b1 = p_row.iloc[0]["b1 T<= 1000 K"]
# calculating reactant heat of formation
r_p1 = float(r_a1)
r_p2 = float(r_a2)*(T/2)
r_p3 = float(r_a3)*((T**2)/3)
r_p4 = float(r_a4)*((T**3)/4)
r_p5 = float(r_a5)*((T**4)/5)
r_p6 = float(r_b1)/T
r_heat = (R*T) * (r_p1 + r_p2 + r_p3 + r_p4 + r_p5 + r_p6)
# calculating product heat of formation
p_p1 = float(p_a1)
p_p2 = float(p_a2)*(T/2)
p_p3 = float(p_a3)*((T**2)/3)
p_p4 = float(p_a4)*((T**3)/4)
p_p5 = float(p_a5)*((T**4)/5)
p_p6 =float(p_b1)/T
p_heat = (R*T) * (p_p1 + p_p2 + p_p3 + p_p4 + p_p5 + p_p6)
h_rxn = p_heat - r_heat # kJ/mol
h_rxn_kJg = h_rxn/mw_r
print("Scenario 1 Results")
print(h_rxn)
print(h_rxn_kJg)
print("-------------------")
# Scenario 2: Heats of reaction returned from functions
from_func = h_of_rxn(rxn22, thermo_df, T, R)
kJg_from_func = convert_heat_g(rxn22, thermo_df, from_func[3])
print("Scenario 2 Results")
print(from_func[3])
print(kJg_from_func)
print("-------------------")
# Scenario 3: Calculate heats of reaction from hard coded nasa polynomial values
r_a1 = 7.93959940e+00
r_a2 = 5.29306627e-03
r_a3 = 0.00000000e+00
r_a4 = 0.00000000e+00
r_a5 = 0.00000000e+00
r_b1 = -4.00155073e+04
p_a1 = 2.20664321E+00
p_a2 = 1.00970086E-02
p_a3 = -9.96338809E-06
p_a4 = 5.47155623E-09
p_a5 = -1.27733965E-12
p_b1 = -4.83529864E+04
# calculating reactant heat of formation
r_p1 = float(r_a1)
r_p2 = float(r_a2)*(T/2)
r_p3 = float(r_a3)*((T**2)/3)
r_p4 = float(r_a4)*((T**3)/4)
r_p5 = float(r_a5)*((T**4)/5)
r_p6 = float(r_b1)/T
r_heat = (R*T) * (r_p1 + r_p2 + r_p3 + r_p4 + r_p5 + r_p6)
# calculating product heat of formation
p_p1 = float(p_a1)
p_p2 = float(p_a2)*(T/2)
p_p3 = float(p_a3)*((T**2)/3)
p_p4 = float(p_a4)*((T**3)/4)
p_p5 = float(p_a5)*((T**4)/5)
p_p6 =float(p_b1)/T
p_heat = (R*T) * (p_p1 + p_p2 + p_p3 + p_p4 + p_p5 + p_p6)
h_rxn = p_heat - r_heat # kJ/mol
h_rxn_kJg = h_rxn/mw_r
print("Scenario 3 Results")
print(h_rxn)
print(h_rxn_kJg)
print("-------------------")
In [9]:
# Test to make sure calculation from nasa polynomial coefficients is accurate
# Methane
# delHº(298.15)/R should equal -8.97226656E+03
R = 8.314/1000 #ideal gas constant in kcal/(mol*K)
T = 298.15
species = "CH4"
heat_from_fn = nasa_polynomial(species, thermo_df, R) / R
print(f"Heat from `nasa_polynomial` function: {heat_from_fn}")
a1 = 5.14911468E+00
a2 = -1.36622009E-02
a3 = 4.91453921E-05
a4 = -4.84246767E-08
a5 = 1.66603441E-11
b1 = -1.02465983E+04
b2 = -4.63848842E+00
p1 = a1
p2 = a2*(T/2)
p3 = a3*((T**2)/3)
p4 = a4*((T**3)/4)
p5 = a5*((T**4)/5)
p6 = b1/T
heat = (T) * (p1 + p2 + p3 + p4 + p5 + p6)
print(f"Heat calculated from hard coded nasa polynomial coefficients: {heat}")
In [ ]: