commit 1824abcb81f2b1b043c64e661b5d07f61a5335eb
Author: ignis <ignis@postech.ac.kr>
Date:   Wed Nov 17 21:13:31 2021 +0900

    initial coal caluator neglecting moisture and ash contents

diff --git a/coal.py b/coal.py
new file mode 100644
index 0000000..fc40600
--- /dev/null
+++ b/coal.py
@@ -0,0 +1,146 @@
+from functools import reduce
+
+import cantera as ct
+
+# Air stream Temperatures and Mass flow rates
+t = list(map(float, ''' 348.15  315.65  600.79  308.54   318.03  306.11  339.45 '''.split()))
+m = list(map(float, ''' 114.4   6.94    362.92  7.25     7.25    7.25    7.25 '''.split())) 
+
+coalT = 348.15
+coalMfr = 56.8813
+
+moist = 0.1551
+vm = 0.3047
+fc = 0.4347
+ash = 0.1055
+
+vm_fc = 1. - moist - ash
+fuelMfr = coalMfr * vm_fc
+
+gri_species = {S.name: S for S in ct.Species.listFromFile('gri30.cti')}
+nasa_species = {S.name: S for S in ct.Species.listFromFile('nasa_gas.cti')}
+nasa_species.update(gri_species)
+system_species = [nasa_species[S] for S in (str.split(' C S H2 O2 N2 CO CO2 H2O SO2 '))]
+
+
+"""#############################################################################
+Calculate averaged air property
+#############################################################################"""
+
+gases = [ct.Solution(thermo='IdealGas', species=system_species) for i in range(7)]
+for i, gas in enumerate(gases):
+    gas.TPY = t[i], ct.one_atm, "O2: 0.233, N2: 0.767"
+
+airs = [
+    ct.Quantity(gas, constant='HP', mass=m[i])
+    for i, gas in enumerate(gases)
+    ]
+
+for i, air in enumerate(airs):
+    pass # print("air", i, air.report())
+
+airmix = reduce(lambda a, b: a+b, airs)
+print("Total Air  flow rate    = ", airmix.mass)
+
+
+"""#############################################################################
+Dummy gaseous coal object containing elementary composition
+#############################################################################"""
+
+coal = ct.Solution(thermo='IdealGas', species=system_species)
+coal.TPY = 348.15, ct.one_atm, '''\
+    C:  81.41,
+    H2:  5.47,
+    O2: 10.83,
+    N2:  1.71,
+    S:   0.57
+    '''
+
+
+"""#############################################################################
+Coal + v O2 = (Y_C/W_C) CO2 + (Y_H/W_H/2) H2O + (Y_S/W_S) SO2 + (Y_N/W_N) N2
+Heating value = Sum(Product Enthalpy of Formation) - Sum(Reactant Enthalpy of Formation)
+#############################################################################"""
+
+# coal heating value = 5761 Kcal/kg * 4.184 kJ/kcal
+coalHV = 5761 * 4184.               # J/kg, with Moisture and Ash
+coalHVdaf = - 5761 * 4.184 / vm_fc  # kJ/kg, Dry and Ash Free, Negative since exothermic
+print("HV of Coal(daf), kJ/kg  = ", coalHVdaf)
+
+# sum product of 1kg coal enthalpy of formation - kJ/kg
+sum_product_hf = (
+        -  32762.45348 * coal.mass_fraction_dict()['C']
+        - 141887.60121 * coal.mass_fraction_dict()['H2']
+        -   8919.38250 * coal.mass_fraction_dict()['S'])
+print("Sum(Hf_product), kJ/kg  = ", sum_product_hf)
+
+sum_coal_hf = - coalHVdaf + sum_product_hf
+print("Sum(Hf_reactant), kJ/kg = ", sum_coal_hf)
+
+
+"""#############################################################################
+Coal Enthalpy at 348.15 K = \Delta H_f + (H(348.15) - H(298.15))
+ - Use graphite's Cp to calculate H difference (H(348.15) - H(298.15))
+#############################################################################"""
+
+gr = ct.Solution('graphite.cti')
+gr.TP = coalT, ct.one_atm
+coal_preheat_enthalpy = gr.enthalpy_mass / 1000.    # kJ/kg
+print("Coal preheat H  , kJ/kg = ", coal_preheat_enthalpy)
+
+coal_enthalpy = sum_coal_hf + coal_preheat_enthalpy
+print("Coal enthalpy   , kJ/kg = ", coal_enthalpy)
+print("Dummy Coal H    , kJ/kg = ", coal.enthalpy_mass/1000.)
+
+
+"""#############################################################################
+Can't Set coal object's H to coal_enthalpy since it is composed of gaseous
+    C S H2 O2 N2 and can't have H=coal_enthalpy with T >= 0 K.
+Therefore only difference between real coal enthalpy and dummy gas coal is 
+    caculated and enthalpy difference will be added later.
+#############################################################################"""
+
+enthalpy_added_after_mixing = (coal_enthalpy*1000 - coal.enthalpy_mass) * fuelMfr # J
+print("enthalpy to add later   = ", enthalpy_added_after_mixing)
+
+
+################################################################################
+## Adiabatic Flame Temperature
+##  - also adjust enthalpy difference btw dummy gas fuel and real solid fuel
+################################################################################
+
+reactant = airmix + ct.Quantity(coal, constant='HP', mass=fuelMfr)
+print("Total Mass flow rate    = ", reactant.mass)
+
+
+################################################################################
+## Adiabatic Flame Temperature
+##  - also adjust enthalpy difference btw dummy gas fuel and real solid fuel
+################################################################################
+
+reactant.equilibrate("HP")
+reactant.HP = reactant.enthalpy/reactant.mass + enthalpy_added_after_mixing/reactant.mass, ct.one_atm
+reactant.equilibrate("HP")
+
+print("")
+print("Chemical Equilibrium")
+reactant.phase()
+
+
+################################################################################
+## Heat Transfer to pipes and walls
+################################################################################
+
+heatXferPipe = 677656911.111 # J/s
+heatXferWall = 592070866.667 # J/s
+
+reactant.HP = reactant.enthalpy/reactant.mass - (heatXferWall+heatXferPipe)/reactant.mass, ct.one_atm
+print("After Heat Transfer")
+reactant.phase()
+
+note_ref='''\
+Reference for Estimation of Coal Enthalpy of Formation
+
+Eqn (14) and (15) of
+Sciazko, M. (2013). Rank-dependent formation enthalpy of coal. Fuel, 114, 2-9. doi:10.1016/j.fuel.2012.06.099
+'''