Adiabatic Flame Temperature calculation
Introduction
Adiabatic flame temperature is the highest temperature that can be obtained in a combustion process during a constant pressure/ constant volume cycle.
Since there is no heat loss, finding out the adiabatic flame temperature can be useful to find the maximum possible temperature that a fuel can reach during combustion.
Using the Newton-Raphson techqnique, when the enthalpy of products and reactants are equal, the root can be found which gives the AFT.
Aim
1. Adiabatic flame temperature variation with Equivalence ratio (In constant volume-chamber)
2. Adiabatic flame temperature variation with Heat loss (In constant pressure-chamber)
3. Adiabatic flame temperature variation with respect to Alkane, Alkene, Alkyne
4. Adiabatic flame temperature variation with respect to a number of Carbon atom chains
Code, explaination and results
Case 1
import cantera as ct
import math
import numpy
import matplotlib.pyplot as plt
#creating the equivalence ratio
er = numpy.linspace(0.1,2.0,200)
#creating a function that returns enthalpy
def h(T,co_eff):
R = 8.314
a1 = co_eff[0]
a2 = co_eff[1]
a3 = co_eff[2]
a4 = co_eff[3]
a5 = co_eff[4]
a6 = co_eff[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_coeff_l = [5.14987613E+00,-1.36709788E-02,4.91800599E-05 ,-4.84743026E-08,1.66693956E-11,-1.02466476E+04,-4.64130376E+00]
o2_coeff_l = [3.78245636E+00, -2.99673416E-03,9.84730201E-06,-9.68129509E-09,3.24372837E-12,-1.06394356E+03,3.65767573E+00]
n2_coeff_h = [0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04,0.05980528E+02]
n2_coeff_l = [0.03298677E+02,0.14082404E-02,-0.03963222E-04,0.05641515E-07,-0.02444854E-10,-0.10208999E+04,0.03950372E+02]
co2_coeff_h = [3.85746029E+00,4.41437026E-03,-2.21481404E-06,5.23490188E-10,-4.72084164E-14,-4.87591660E+04,2.27163806E+00]
h20_coeff_h = [3.03399249E+00,2.17691804E-03,-1.64072518E-07,-9.70419870E-11,1.68200992E-14,-3.00042971E+04,4.96677010E+00]
o2_coeff_h = [3.28253784E+00, 1.48308754E-03,-7.57966669E-07, 2.09470555E-10,-2.16717794E-14 ,-1.08845772E+03 ,5.45323129E+00]
co_coeff_h = [2.71518561E+00, 2.06252743E-03,-9.98825771E-07, 2.30053008E-10,-2.03647716E-14 ,-1.41518724E+04, 7.81868772E+00]
#creating a function that returns LHV
def f(T,er):
R = 8.314
h_co2_p = h(T,co2_coeff_h)
h_n2_p = h(T,n2_coeff_h)
h_h2o_p = h(T,h20_coeff_h)
h_co_p = h(T,co_coeff_h)
h_o2_p = h(T,o2_coeff_h)
if er > 1:
H_products = ((4/er)-3)*h_co2_p + (4-(4/er))*h_co_p + 2*h_h2o_p + (7.52/er)*h_n2_p
N_products = (((4/er)-3) + 2 + (4-4/er) + (7.52/er))*R*T
else:
H_products = h_co2_p + ((2/er)-2)*h_o2_p + 2*h_h2o_p + (7.52/er)*h_n2_p
N_products = (((2/er)-2) + 3 + (7.52/er))*R*T
Tstd = 298.15
h_ch4_r = h(Tstd,ch4_coeff_l)
h_o2_r = h(Tstd,o2_coeff_l)
h_n2_r = h(Tstd,n2_coeff_l)
H_reactants = h_ch4_r + (2/er)*h_o2_r + ((2/er)*3.76)*h_n2_r
N_reactants = (1+(2/er)+((2/er)*3.76))*R*Tstd
return (H_reactants - H_products) - (N_reactants-N_products)
def fprime(T,er):
return(f(T+1e-6,er)-f(T,er))/1e-6
T_guess = 1500
tol = 1e-3
ct = 0
alpha = 0.2
temp = []
for j in er:
while(abs(f(T_guess,j))>tol):
T_guess = T_guess - alpha*((f(T_guess,j)/fprime(T_guess,j)))
temp.append(T_guess)
plt.plot(er,temp,color = 'red')
plt.xlabel('phi')
plt.ylabel('AFT (K)')
plt.show()
It can be seen that after an equivalence ratio of 1, there is a reduction in flame temperature indicating that a rich mixture allows lower flame temperatures to be produced.
case 2
import cantera as ct
import math
import numpy
import matplotlib.pyplot as plt
#creating the equivalence ratio
er = 1
#creating a function that returns enthalpy
def h(T,co_eff):
R = 8.314
a1 = co_eff[0]
a2 = co_eff[1]
a3 = co_eff[2]
a4 = co_eff[3]
a5 = co_eff[4]
a6 = co_eff[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_coeff_l = [5.14987613E+00,-1.36709788E-02,4.91800599E-05 ,-4.84743026E-08,1.66693956E-11,-1.02466476E+04,-4.64130376E+00]
o2_coeff_l = [3.78245636E+00, -2.99673416E-03,9.84730201E-06,-9.68129509E-09,3.24372837E-12,-1.06394356E+03,3.65767573E+00]
n2_coeff_h = [0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04,0.05980528E+02]
n2_coeff_l = [0.03298677E+02,0.14082404E-02,-0.03963222E-04,0.05641515E-07,-0.02444854E-10,-0.10208999E+04,0.03950372E+02]
co2_coeff_h = [3.85746029E+00,4.41437026E-03,-2.21481404E-06,5.23490188E-10,-4.72084164E-14,-4.87591660E+04,2.27163806E+00]
h20_coeff_h = [3.03399249E+00,2.17691804E-03,-1.64072518E-07,-9.70419870E-11,1.68200992E-14,-3.00042971E+04,4.96677010E+00]
o2_coeff_h = [3.28253784E+00, 1.48308754E-03,-7.57966669E-07, 2.09470555E-10,-2.16717794E-14 ,-1.08845772E+03 ,5.45323129E+00]
co_coeff_h = [2.71518561E+00, 2.06252743E-03,-9.98825771E-07, 2.30053008E-10,-2.03647716E-14 ,-1.41518724E+04, 7.81868772E+00]
H_loss = numpy.linspace(0.1,1,200)
max_heatloss = 80200
def f(T,H_loss):
R = 8.314
h_co2_p = h(T,co2_coeff_h)
h_n2_p = h(T,n2_coeff_h)
h_h2o_p = h(T,h20_coeff_h)
h_co_p = h(T,co_coeff_h)
h_o2_p = h(T,o2_coeff_h)
H_products = h_co2_p + ((2/er)-2)*h_o2_p + 2*h_h2o_p + (7.52/er)*h_n2_p
Tstd = 298.5
h_ch4_r = h(Tstd,ch4_coeff_l)
h_o2_r = h(Tstd,o2_coeff_l)
h_n2_r = h(Tstd,n2_coeff_l)
H_reactants = h_ch4_r + (2/er)*h_o2_r + ((2/er)*3.76)*h_n2_r
return (H_products - H_reactants) + (H_loss*max_heatloss)
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
alpha = 0.9
temp = []
for j in H_loss:
while(abs(f(T_guess,j))>tol):
T_guess = T_guess - alpha*((f(T_guess,j)/fprime(T_guess,j)))
temp.append(T_guess)
p = numpy.searchsorted(H_loss,0.35)
print(p)
temp_aft = temp[p]
print(temp_aft)
plt.plot(H_loss,temp,color = 'red')
plt.xlabel('Heat loss coefficient')
plt.ylabel('AFT (K)')
plt.grid('on')
plt.show()
It can be seen that there is a linear variation as the heat loss increases, which in turn reduces the AFT.
AFT at 0.35 = 
case 3
import cantera as ct
import math
import numpy
import matplotlib.pyplot as plt
#creating a function that returns enthalpy
def h(T,co_eff):
R = 8.314
a1 = co_eff[0]
a2 = co_eff[1]
a3 = co_eff[2]
a4 = co_eff[3]
a5 = co_eff[4]
a6 = co_eff[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
c2h2_coeff_l = [8.08681094E-01, 2.33615629E-02,-3.55171815E-05, 2.80152437E-08,-8.50072974E-12 ,2.64289807E+04, 1.39397051E+01]
c2h4_coeff_l = [3.95920148E+00,-7.57052247E-03, 5.70990292E-05 ,-6.91588753E-08 ,2.69884373E-11, 5.08977593E+03, 4.09733096E+00 ]
c2h6_coeff_l = [4.29142492E+00,-5.50154270E-03, 5.99438288E-05, -7.08466285E-08, 2.68685771E-11,-1.15222055E+04, 2.66682316E+00 ]
ch4_coeff_l = [5.14987613E+00,-1.36709788E-02,4.91800599E-05 ,-4.84743026E-08,1.66693956E-11,-1.02466476E+04,-4.64130376E+00]
o2_coeff_l = [3.78245636E+00, -2.99673416E-03,9.84730201E-06,-9.68129509E-09,3.24372837E-12,-1.06394356E+03,3.65767573E+00]
n2_coeff_h = [0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04,0.05980528E+02]
n2_coeff_l = [0.03298677E+02,0.14082404E-02,-0.03963222E-04,0.05641515E-07,-0.02444854E-10,-0.10208999E+04,0.03950372E+02]
co2_coeff_h = [3.85746029E+00,4.41437026E-03,-2.21481404E-06,5.23490188E-10,-4.72084164E-14,-4.87591660E+04,2.27163806E+00]
h20_coeff_h = [3.03399249E+00,2.17691804E-03,-1.64072518E-07,-9.70419870E-11,1.68200992E-14,-3.00042971E+04,4.96677010E+00]
o2_coeff_h = [3.28253784E+00, 1.48308754E-03,-7.57966669E-07, 2.09470555E-10,-2.16717794E-14 ,-1.08845772E+03 ,5.45323129E+00]
co_coeff_h = [2.71518561E+00, 2.06252743E-03,-9.98825771E-07, 2.30053008E-10,-2.03647716E-14 ,-1.41518724E+04, 7.81868772E+00]
Tstd = 298.15
h_c2h2_r = h(Tstd,c2h2_coeff_l)
h_c2h4_r = h(Tstd,c2h4_coeff_l)
h_c2h6_r = h(Tstd,c2h6_coeff_l)
H_alkane = h_c2h6_r
H_alkene = h_c2h4_r
H_alkyne = h_c2h2_r
H = [H_alkane,H_alkene,H_alkyne]
def f(T,H):
h_co2_p = h(T,co2_coeff_h)
h_n2_p = h(T,n2_coeff_h)
h_h2o_p = h(T,h20_coeff_h)
h_o2_r = h(Tstd,o2_coeff_l)
h_n2_r = h(Tstd,n2_coeff_l)
if H == H_alkane:
x=2
y=6
a = x + y/4
b = y/2
H_reactants = H_alkane + a*h_o2_r + 3.76*a*h_n2_r
H_products = x*h_co2_p + b*h_h2o_p + 3.76*a*h_n2_p
elif H == H_alkene:
x=2
y=4
a = x + y/4
b = y/2
H_reactants = H_alkene + a*h_o2_r + 3.76*a*h_n2_r
H_products = x*h_co2_p + b*h_h2o_p + 3.76*a*h_n2_p
elif H == H_alkyne:
x=2
y=6
a = x + y/4
b = y/2
H_reactants = H_alkyne + a*h_o2_r + 3.76*a*h_n2_r
H_products = x*h_co2_p + b*h_h2o_p + 3.76*a*h_n2_p
return (H_products - H_reactants)
def fprime(T,H):
return(f(T+1e-6,H)-f(T,H))/1e-6
T_guess = 1500
tol = 1e-3
ct = 0
alpha = 0.1
Temp = []
for i in H:
while(abs(f(T_guess,i))>tol):
T_guess = T_guess - alpha*((f(T_guess,i)/fprime(T_guess,i)))
Temp.append(T_guess)
plt.ylabel('AFT (K)')
plt.bar(['Alkane','Alkene','Alkyne'],Temp)
plt.show()
There is a lack of h20 produced during alkyene combustion, thus leading to high flame temperatures.
case 4
import cantera as ct
import math
import numpy
import matplotlib.pyplot as plt
#creating a function that returns enthalpy
def h(T,co_eff):
R = 8.314
a1 = co_eff[0]
a2 = co_eff[1]
a3 = co_eff[2]
a4 = co_eff[3]
a5 = co_eff[4]
a6 = co_eff[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_coeff_l = [5.14987613E+00,-1.36709788E-02,4.91800599E-05 ,-4.84743026E-08,1.66693956E-11,-1.02466476E+04,-4.64130376E+00]
c2h6_coeff_l = [4.29142492E+00,-5.50154270E-03,5.99438288E-05, -7.08466285E-08, 2.68685771E-11,-1.15222055E+04,2.66682316E+00 ]
c3h8_coeff_l = [0.93355381E+00, 0.26424579E-01, 0.61059727E-05,-0.21977499E-07, 0.95149253E-11,-0.13958520E+05, 0.19201691E+02]
o2_coeff_l = [3.78245636E+00, -2.99673416E-03,9.84730201E-06,-9.68129509E-09,3.24372837E-12,-1.06394356E+03,3.65767573E+00]
n2_coeff_h = [0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04,0.05980528E+02]
n2_coeff_l = [0.03298677E+02,0.14082404E-02,-0.03963222E-04,0.05641515E-07,-0.02444854E-10,-0.10208999E+04,0.03950372E+02]
co2_coeff_h = [3.85746029E+00,4.41437026E-03,-2.21481404E-06,5.23490188E-10,-4.72084164E-14,-4.87591660E+04,2.27163806E+00]
h20_coeff_h = [3.03399249E+00,2.17691804E-03,-1.64072518E-07,-9.70419870E-11,1.68200992E-14,-3.00042971E+04,4.96677010E+00]
Tstd = 298.15
h_ch4_r = h(Tstd,ch4_coeff_l)
h_c2h6_r = h(Tstd,c2h6_coeff_l)
h_c3h8_r = h(Tstd,c3h8_coeff_l)
C1 = h_ch4_r
C2 = h_c2h6_r
C3 = h_c3h8_r
C = [C1,C2,C3]
def f(T,C):
h_co2_p = h(T,co2_coeff_h)
h_n2_p = h(T,n2_coeff_h)
h_h2o_p = h(T,h20_coeff_h)
h_o2_r = h(Tstd,o2_coeff_l)
h_n2_r = h(Tstd,n2_coeff_l)
if C == C1:
x = 1
y = 2*(x+1)
a = x + y/4
b = y/2
H_reactants = C1 + a*h_o2_r + 3.76*a*h_n2_r
H_products = x*h_co2_p + b*h_h2o_p + 3.76*a*h_n2_p
elif C == C2:
x = 2
y = 2*(x+1)
a = x + y/4
b = y/2
H_reactants = C2 + a*h_o2_r + 3.76*a*h_n2_r
H_products = x*h_co2_p + b*h_h2o_p + 3.76*a*h_n2_p
elif C == C3:
x = 3
y = 2*(x+1)
a = x + y/4
b = y/2
H_reactants = C3 + a*h_o2_r + 3.76*a*h_n2_r
H_products = x*h_co2_p + b*h_h2o_p + 3.76*a*h_n2_p
return (H_products-H_reactants)
def fprime(T,C):
return(f(T+1e-6,C)-f(T,C))/1e-6
T_guess = 1900
tol = 1e-3
ct = 0
alpha = 0.1
Temp = []
for i in C:
while(abs(f(T_guess,i))>tol):
T_guess = T_guess - alpha*((f(T_guess,i)/fprime(T_guess,i)))
Temp.append(T_guess)
plt.ylabel('AFT (K)')
plt.bar(['1 Carbon','2 carbon','3 carbon'],Temp)
plt.show()
A slight increase in AFT can be observed due to an increase of carbon atoms.