Adiabatic Flame Temperature calculation by using Python

Combustion of Methane Gases

A common example of the burning of gaseous fuel is the combustion of natural gas which mainly contains methane. One molecule of methane reacts with two molecules of oxygen to produce one molecule of carbon dioxide and two molecules of water vapor.

Equation:-

CH4 + 2O2 == CO2 + 2H2O

In a combustion process, the complete burning cannot be accomplished with the exact theoretical amount of air. Therefore it is essential to supply an excessive amount of air for the combustion of fuel. The volume of twice the theoretical air may be used; otherwise, it is decided based on the nature of the fuel used and the method of the burning process. The valid reason for incomplete combustion of fuel in presence of the theoretical amount of air may be due to the fact that each particle of oxygen in the air does not have the intimate contact with the fuel particles in the combustion process. Therefore an excess air is to be supplied.

CH4 +2O2+7.524 N2 = CO2+2H2O+7.524 N2

The stoichiometric equation is CH4 +2O2+7.524 N2 = CO2+2H2O+7.524 N2 ( 2 moles of oxygen required 7.524 moles of nitrogen, that means 9.524 moles of air.

1 mole of methane requires 9.524 moles of stoichiometric air or 9.524 ml of air at STP per ml of methane.

1. Program for Effect of equivalence ratio on AFT

\"\"\"

@author: Chocolate

Calculating the AFT of Methane at different equivalence ratios

CH4+2(O2+3.76N2)=CO2+2h2o+7.52N2
\"\"\"
import matplotlib.pyplot as plt
import math

def h(T, co_effs):
R =8.314 #j/mol-k
a1 = co_effs[0]
a2 = co_effs[1]
a3 = co_effs[2]
a4 = co_effs[3]
a5 = co_effs[4]
a6 = co_effs[5]

return (a1+a2*T/2+a3*pow(T,2)/3+a4*pow(T,3)/4+a5*pow(T,4)/5+a6/T)*R*T

ch4_coeffs_l =[5.14987613E+00,-1.36709788E-02,4.91800599E-05,-4.84743026E-08,1.66693956E-11,-1.02466476E+04,-4.64130376E+00]
o2_coeffs_l = [3.78245636E+00,-2.99673416E-03,9.84730201E-06,-9.68129509E-09,3.24372837E-12,-1.06394356E+03,3.65767573E+00]
n2_coeffs_h = [0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04,0.05980528E+02]
n2_coeffs_l= [0.03298677E+02,0.14082404E-02,-0.03963222E-04,0.05641515E-07,-0.02444854E-10,-0.10208999E+04,0.03950372E+02]
co2_coeffs_h = [3.85746029E+00,4.41437026E-03,-2.21481404E-06,5.23490188E-10,-4.72084164E-14,-4.87591660E+04,2.27163806E+00]
h2o_coeffs_h = [3.03399249E+00,2.17691804E-03,-1.64072518E-07,-9.70419870E-11,1.68200992E-14,-3.00042971E+04,4.96677010E+00]

def f(T):
\"\"\" Products of reactions\"\"\"

Tstd = 298.15

h_co2_p = h(T,co2_coeffs_h)
h_n2_p = h(T,n2_coeffs_h)
h_h2o_p = h(T,h2o_coeffs_h)

H_products = h_co2_p + 7.52*h_n2_p + 2*h_h2o_p
print(H_products)

\"\"\"Reactants of reactions\"\"\"

h_ch4_r = h(Tstd,ch4_coeffs_l)
h_o2_r = h(Tstd,o2_coeffs_l)
h_n2_r = h(Tstd,n2_coeffs_l)

H_reactants = h_ch4_r + 2*h_o2_r + 7.52*h_n2_r
print(H_reactants)

return H_products - H_reactants

def fprime(T):
return (f(T+1e-6) - f(T))/1e-6

T_guess = 1500

tol = 1e-3
ct = 0
alpha = 0.2

while(abs(f(T_guess))> tol):
T_guess = T_guess - alpha*((f(T_guess))/(fprime(T_guess)))
plt.plot(ct,T_guess,\'*\',color=\'red\')
ct = ct+1
print(T_guess)

plt.show()
print(ct)
print(T_guess)

Use of Cantera

import matplotlib.pyplot as plt
import cantera as ct
gas=ct.Solution(\'gri30.xml\')
phi=0.05
for i in range(0,15):
gas.TPX= 298,101325,{\'CH4\':1,\'O2\':2/phi,\'N2\':7.52/phi}
gas.set_equivalence_ratio(phi,\'CH4\',\'O2:1, N2:3.76\')
gas.equilibrate(\'UV\',\'auto\')
print(phi, gas.T)
plt.plot(phi,gas.T,\'*\',color=\'red\')
phi=phi+0.1
plt.xlabel(\'Equivalence ratio\')
plt.ylabel(\'Adiabatic Flame Temperature\')
plt.show()

2. Program for Constant pressure reactor with heat loss

A. Effect of Heat loss on the Flame Temperature

import matplotlib.pyplot as plt
import numpy as np

def h(T,coeffs):

R = 8.314 #j/mol-k

a1=coeffs[0]
a2=coeffs[1]
a3=coeffs[2]
a4=coeffs[3]
a5=coeffs[4]
a6=coeffs[5]

return (a1 + (a2*T/2) + (a3*pow(T,2)/3) + (a4*pow(T,3)/4) + (a5*pow(T,4)/5) + (a6/T))*R*T

#low temperature coefficients
CH4_coeffs_l=[5.14987613E+00, -1.36709788E-02, 4.91800599E-05, -4.84743026E-08, 1.66693956E-11, -1.02466476E+04]
O2_coeffs_l=[3.78245636E+00, -2.99673416E-03, 9.84730201E-06, -9.68129509E-09, 3.24372837E-12, -1.06394356E+03]
N2_coeffs_l=[0.03298677E+02, 0.14082404E-02, -0.03963222E-04, 0.05641515E-07, -0.02444854E-10, -0.10208999E+04]

#high temperature coefficients
N2_coeffs_h=[0.02926640E+02, 0.14879768E-02, -0.05684760E-05, 0.10097038E-09, -0.06753351E-13, -0.09227977E+04 ]
CO2_coeffs_h=[3.85746029E+00, 4.41437026E-03, -2.21481404E-06, 5.23490188E-10, -4.72084164E-14, -4.87591660E+04]
H20_coeffs_h=[3.03399249E+00,2.17691804E-03,-1.64072518E-07,-9.70419870E-11,1.68200992E-14,-3.00042971E+04]

def f(T,H_loss):

R = 8.314 # J/mol.k
Tstd = 298.15
LHV = 800950

h_CH4_r=h(Tstd,CH4_coeffs_l)
h_O2_r=h(Tstd,O2_coeffs_l)
h_N2_r=h(Tstd,N2_coeffs_l)

h_CO2_p=h(T,CO2_coeffs_h)
h_N2_p=h(T,N2_coeffs_h)
h_H2O_p=h(T,H20_coeffs_h)

h_reactants = h_CH4_r + 2*h_O2_r + 7.52*h_N2_r
h_products = h_CO2_p + 2*h_H2O_p +7.52*h_N2_p

loss = H_loss*LHV

return h_products - h_reactants + loss

def fprime(T,H_loss):
return (f(T+1e-6,H_loss)-f(T,H_loss))/1e-6

T_guess = 1500
tol = 1e-3
ct = 0
H_loss = np.linspace(0,1,10)
alpha=1
T_ft = []

for i in range (0,len(H_loss)):
while(abs(f(T_guess,H_loss[i])) > tol ):
T_guess = T_guess - alpha*((f(T_guess,H_loss[i]) / fprime(T_guess,H_loss[i])))
ct = ct+1
T_ft.append(T_guess)

print(T_ft)

plt.plot(H_loss,T_ft, \'-.*\')
plt.xlabel(\'Heat Loss factor\')
plt.ylabel(\'Flame Temperature(K)\')
plt.title(\'Effect of Heat loss on the Flame Temperature\')
plt.show()

B. Estimation of FT for a Heat Loss Factor of 0.35

\"\"\"
THIS CODE IS TO CALCULATE THE AFT OF ALKANE, ALKENE, ALKYNE FOR 2 CARBON ATOM AND TO COMPARE WITH EACH OTHER

Alkane (ETHANE)
C2H6 + 3.5O2 +3.5 *3.76 N2 = 2CO2 + 3H2O + 3.5 *3.76 N2 + 0.35* LHV

Alkene (ETHENE)
C2H4 + 3O2 + 3 *3.76 N2 = 2CO2 + 2H2O + 3 *3.76 N2 + 0.35* LHV

Alkyne (ETHYNE)
C2H2 + 2.5O2 + 2.5 *3.76 N2 = 2CO2 + H2O + 3 *3.76 N2 + 0.35* LHV

\"\"\"
import matplotlib.pyplot as plt
import numpy as np

def h(T,coeffs):

R = 8.314 #j/mol-k

a1=coeffs[0]
a2=coeffs[1]
a3=coeffs[2]
a4=coeffs[3]
a5=coeffs[4]
a6=coeffs[5]
a7=coeffs[6]

return (a1 + (a2*T/2) + (a3*pow(T,2)/3) + (a4*pow(T,3)/4) + (a5*pow(T,4)/5) + (a6/T))*R*T

c2h6_coeffs_l =[4.29142492E+00,-5.50154270E-03 ,5.99438288E-05 ,-7.08466285E-08 ,2.68685771E-11,-1.15222055E+04 ,2.66682316E+00]
c2h4_coeffs_l =[3.95920148E+00,-7.57052247E-03 ,5.70990292E-05 ,-6.91588753E-08 ,2.69884373E-11 ,5.08977593E+03 ,4.09733096E+00]
c2h2_coeffs_l =[8.08681094E-01 ,2.33615629E-02,-3.55171815E-05 ,2.80152437E-08,-8.50072974E-12 ,2.64289807E+04 ,1.39397051E+01]
n2_coeffs_l = [0.03298677E+02, 0.14082404E-02,-0.03963222E-04, 0.05641515E-07,-0.02444854E-10,-0.10208999E+04, 0.03950372E+02]

n2_coeffs_h = [0.02926640E+02, 0.14879768E-02,-0.05684760E-05, 0.10097038E-09,-0.06753351E-13,-0.09227977E+04, 0.05980528E+02]
o2_coeffs_l = [3.78245636E+00,-2.99673416E-03, 9.84730201E-06,-9.68129509E-09, 3.24372837E-12,-1.06394356E+03, 3.65767573E+00]
co2_coeffs_h = [3.85746029E+00, 4.41437026E-03,-2.21481404E-06, 5.23490188E-10,-4.72084164E-14,-4.87591660E+04, 2.27163806E+00]
h2o_coeffs_h = [3.03399249E+00, 2.17691804E-03,-1.64072518E-07,-9.70419870E-11, 1.68200992E-14,-3.00042971E+04, 4.96677010E+00]

\"\"\"
Alkane (ETHANE)
C2H6 + 3.5O2 +3.5 *3.76 N2 = 2CO2 + 3H2O + 3.5 *3.76 N2 + 0.35* LHV
\"\"\"
def f_alkane(T):

R=8.314
Tstd = 298.15
LHV = 1428275 #J/mole

h_c2h6_r = h(Tstd,c2h6_coeffs_l)
h_o2_r = h(Tstd,o2_coeffs_l)
h_n2_r = h(Tstd,n2_coeffs_l)

h_co2_p = h(T,co2_coeffs_h)
h_h2o_p = h(T,h2o_coeffs_h)
h_n2_p = h(T,n2_coeffs_h)

H_reactants = h_c2h6_r + 3.5*h_o2_r + (3.5*3.76*h_n2_r)
H_products = 2*h_co2_p+ 3*h_h2o_p + (3.5*3.76*h_n2_p)

return H_products - H_reactants + (0.35*LHV)

def fprime_alkane(T):
return (f_alkane(T+1e-6)-f_alkane(T))/1e-6

\"\"\"
Alkene (ETHENE)
C2H4 + 3O2 + 3 *3.76 N2 = 2CO2 + 2H2O + 3 *3.76 N2 + 0.35* LHV

\"\"\"
def f_alkene(T):

R=8.314
Tstd = 298.15
LHV = 1323474 #J/mole

h_c2h4_r = h(Tstd,c2h4_coeffs_l)
h_o2_r = h(Tstd,o2_coeffs_l)
h_n2_r = h(Tstd,n2_coeffs_l)

h_co2_p = h(T,co2_coeffs_h)
h_h2o_p = h(T,h2o_coeffs_h)
h_n2_p = h(T,n2_coeffs_h)

H_reactants = h_c2h4_r + 3*h_o2_r + (3*3.76*h_n2_r)
H_products = 2*h_co2_p+ 2*h_h2o_p + (3*3.76*h_n2_p)

return H_products - H_reactants + (0.35*LHV)

def fprime_alkene(T):
return (f_alkene(T+1e-6)-f_alkene(T))/1e-6

\"\"\"
Alkyne (ETHYNE)
C2H2 + 2.5O2 + 2.5 *3.76 N2 = 2CO2 + H2O + 3 *3.76 N2 + 0.35* LHV

\"\"\"

def f_alkyne(T):

R=8.314
Tstd = 298.15
LHV = 1258020 #J/mole

h_c2h2_r = h(Tstd,c2h2_coeffs_l)
h_o2_r = h(Tstd,o2_coeffs_l)
h_n2_r = h(Tstd,n2_coeffs_l)

h_co2_p = h(T,co2_coeffs_h)
h_h2o_p = h(T,h2o_coeffs_h)
h_n2_p = h(T,n2_coeffs_h)

H_reactants = h_c2h2_r + 2.5*h_o2_r + (2.5*3.76*h_n2_r)
H_products = 2*h_co2_p+ h_h2o_p + (3*3.76*h_n2_p)

return H_products - H_reactants + (0.35*LHV)

def fprime_alkyne(T):
return (f_alkyne(T+ 1e-6)-f_alkyne(T))/1e-6

T_alkane_guess = 1500
T_alkene_guess = 1500
T_alkyne_guess = 1500

tol = 1e-3
alpha = 1

T_aft = []

#alkane

while(abs(f_alkane(T_alkane_guess)) > tol):
T_alkane_guess = T_alkane_guess - alpha*(f_alkane(T_alkane_guess) / fprime_alkane(T_alkane_guess))

T_aft.append(T_alkane_guess)

#alkene

while(abs(f_alkene(T_alkene_guess)) > tol):
T_alkene_guess = T_alkene_guess - alpha*(f_alkene(T_alkene_guess) / fprime_alkene(T_alkene_guess))

T_aft.append(T_alkene_guess)

#alkyne

while(abs(f_alkyne(T_alkyne_guess)) > tol):
T_alkyne_guess = T_alkyne_guess - alpha*(f_alkyne(T_alkyne_guess) / fprime_alkyne(T_alkyne_guess))

T_aft.append(T_alkyne_guess)

fuel=[\'Ethane\',\'Ethene\',\'Ethyne\']
print(T_aft)
plt.plot(fuel,T_aft, \'-.*\')
plt.xlabel(\'Type of Hydrocarbon\')
plt.ylabel(\' Flame Temperature(K)\')
plt.title(\'Estimation of FT for a Heat Loss Factor of 0.35\')
plt.show()

C. VARIATION IN AFT WITH CHANGE IN NUMBER OF CARBON ATOM

THIS CODE IS TO SEE THE VARIATION IN AFT WITH CHANGE IN NUMBER OF CARBON ATOM

Methane
CH4 + 2O2 +2 *3.76 N2 = CO2 + 2H2O + 2 *3.76 N2

Ethane
C2H6 + 3.5O2 +3.5*3.76 N2 = 2CO2 + 3H2O + 3.5 *3.76 N2

Propane
C3H8 + 5O2 + 5*3.76 N2 = 3CO2 + 4H2O + 5*3.76 N2

import matplotlib.pyplot as plt
import numpy as np

def h(T,co_effs):

R = 8.314 #j/mol-k

a1=co_effs[0]
a2=co_effs[1]
a3=co_effs[2]
a4=co_effs[3]
a5=co_effs[4]
a6=co_effs[5]
a7=co_effs[6]

return (a1 + (a2*T/2) + (a3*pow(T,2)/3) + (a4*pow(T,3)/4) + (a5*pow(T,4)/5) + (a6/T))*R*T

ch4_coeffs_l = [5.14987613E+00,-1.36709788E-02, 4.91800599E-05,-4.84743026E-08, 1.66693956E-11,-1.02466476E+04,-4.64130376E+00]
c2h6_coeffs_l =[4.29142492E+00,-5.50154270E-03, 5.99438288E-05,-7.08466285E-08 ,2.68685771E-11,-1.15222055E+04 ,2.66682316E+00]
c3h8_coeffs_l =[0.93355381E+00, 0.26424579E-01,0.61059727E-05,-0.21977499E-07 ,0.95149253E-11 ,-0.13958520E+05 ,0.19201691E+02]
n2_coeffs_l = [0.03298677E+02, 0.14082404E-02,-0.03963222E-04, 0.05641515E-07,-0.02444854E-10,-0.10208999E+04, 0.03950372E+02]

n2_coeffs_h = [0.02926640E+02, 0.14879768E-02,-0.05684760E-05, 0.10097038E-09,-0.06753351E-13,-0.09227977E+04, 0.05980528E+02]
o2_coeffs_l = [3.78245636E+00,-2.99673416E-03, 9.84730201E-06,-9.68129509E-09, 3.24372837E-12,-1.06394356E+03, 3.65767573E+00]
co2_coeffs_h = [3.85746029E+00, 4.41437026E-03,-2.21481404E-06, 5.23490188E-10,-4.72084164E-14,-4.87591660E+04, 2.27163806E+00]
h2o_coeffs_h = [3.03399249E+00, 2.17691804E-03,-1.64072518E-07,-9.70419870E-11, 1.68200992E-14,-3.00042971E+04, 4.96677010E+00]

# Methane
# CH4 + 2O2 +2 *3.76 N2 = CO2 + 2H2O + 2 *3.76 N2

def f_methane(T):

R=8.314
Tstd = 298.15
LHV = 1428275 #J/mole

h_ch4_r = h(Tstd,ch4_coeffs_l)
h_o2_r = h(Tstd,o2_coeffs_l)
h_n2_r = h(Tstd,n2_coeffs_l)

h_co2_p = h(T,co2_coeffs_h)
h_h2o_p = h(T,h2o_coeffs_h)
h_n2_p = h(T,n2_coeffs_h)

H_reactants = h_ch4_r + 2*h_o2_r + 7.52*h_n2_r
H_products = h_co2_p+ 2*h_h2o_p + 7.52*h_n2_p

return H_products - H_reactants

def fprime_methane(T):
return (f_methane(T+1e-6)-f_methane(T))/1e-6

#Ethane
#C2H6 + 3.5O2 +3.5*3.76 N2 = 2CO2 + 3H2O + 3.5 *3.76 N2

def f_ethane(T):

R=8.314
Tstd = 298.15
LHV = 1428275 #J/mole

h_c2h6_r = h(Tstd,c2h6_coeffs_l)
h_o2_r = h(Tstd,o2_coeffs_l)
h_n2_r = h(Tstd,n2_coeffs_l)

h_co2_p = h(T,co2_coeffs_h)
h_h2o_p = h(T,h2o_coeffs_h)
h_n2_p = h(T,n2_coeffs_h)

H_reactants = h_c2h6_r + 3.5*h_o2_r + (3.5*3.76*h_n2_r)
H_products = 2*h_co2_p+ 3*h_h2o_p + (3.5*3.76*h_n2_p)

return H_products - H_reactants

def fprime_ethane(T):
return (f_ethane(T+1e-6)-f_ethane(T))/1e-6

# Propane
# C3H8 + 5O2 + 5*3.76 N2 = 3CO2 + 4H2O + 5*3.76 N2

def f_propane(T):

R=8.314
Tstd = 298.15
LHV = 1428275 #J/mole

h_c3h8_r = h(Tstd,c3h8_coeffs_l)
h_o2_r = h(Tstd,o2_coeffs_l)
h_n2_r = h(Tstd,n2_coeffs_l)

h_co2_p = h(T,co2_coeffs_h)
h_h2o_p = h(T,h2o_coeffs_h)
h_n2_p = h(T,n2_coeffs_h)

H_reactants = h_c3h8_r + 5*h_o2_r + (5*3.76*h_n2_r)
H_products = 3*h_co2_p+ 4*h_h2o_p + (5*3.76*h_n2_p)

return H_products - H_reactants

def fprime_propane(T):
return (f_propane(T+1e-6)-f_propane(T))/1e-6

T_methane_guess = 1500
T_ethane_guess = 1500
T_propane_guess = 1500

tol = 1e-3
alpha = 1

T_aft = []

#methane

while(abs(f_methane(T_methane_guess)) > tol):
T_methane_guess = T_methane_guess - alpha*(f_methane(T_methane_guess) / fprime_methane(T_methane_guess))

T_aft.append(T_methane_guess)

#ethane

while(abs(f_ethane(T_ethane_guess)) > tol):
T_ethane_guess = T_ethane_guess - alpha*(f_ethane(T_ethane_guess) / fprime_ethane(T_ethane_guess))

T_aft.append(T_ethane_guess)

#propane

while(abs(f_propane(T_propane_guess)) > tol):
T_propane_guess = T_propane_guess - alpha*(f_propane(T_propane_guess) / fprime_propane(T_propane_guess))

T_aft.append(T_propane_guess)

fuel=[\'Methane\',\'Ethane\',\'Propane\']
print(T_aft)
plt.plot(fuel,T_aft, \'-.*\')
plt.xlabel(\'Type of Hydrocarbon\')
plt.ylabel(\' Flame Temperature(K)\')
plt.title(\'Effect of alkane carbon atoms on the AFT\')
plt.show()