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

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):
    print("air {}, T = {}, mass flow rate = {}".format(i+1, air.T, air.mass))

airmix = reduce(lambda a, b: a+b, airs)
print("Total Air  flow rate    = ", airmix.mass)


"""#############################################################################
Dummy gaseous coal object containing elementary composition
    mass of moisture and ash contents is added to N2
#############################################################################"""

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: 36.9647930755,
    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 * 4.184          # 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     , kJ/kg  = ", coalHV)
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 = - coalHV + 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
'''