Week 3 - Adiabatic Flame Temperature calculation
Aim
To determine the adiabatic flame temperature for combustion of hydrocarbon fuels and study the trend for variation in equivalence ratio, number of carbon atoms, type of carbon bond using Python and Cantera
- Introduction
In a chemical reaction, the conversion of reactants to products involves heat transfer in the system. Considering the system to be adiabatic, that is, there is no heat transfer into the environment; the temperature attained inside the system is defined as the Adiabatic Flame Temperature (AFT). AFT is an ideal case temperature and the actual flame temperature observed for a particular fuel is always less than AFT.
In this study, we use programming to calculate AFT for methane combustion and study the influence of equivalence ratio on AFT. The results of programming using Python and Cantera are compared. The study is further extended to calculate flame temperature for a particular heat loss and the influence of number of carbon atoms, type of carbon bonds are studied.
- Theory
2.1 Newton Rhapson
It is one of the approximation method to find roots of a real-valued function f(x) =0. The method is based on finding the minima and maximum values by approximating using a tangent to a continuous and differentiable function. Using x_guess as the initial solution approximation, the new value is calculated as follows:
X_new = x_guess – f(x) / f’(x_guess)
Where f’(x_guess) is the differential value of f(x) at point x_guess. The derivative can be calculated using finite difference method (Taylor series approximation).
2.2 Enthalpy equation for constant volume combustion
Using the enthalpy definition and first law of thermodynamics, we can write the equations as:
h = U + pV
Applying differentials on both sides, we get
dh = dU + Vdp (dV=0 for constant volume process)
For a constant volume process, there is no work done by the system due to no change in the volume. Hence, the heat generated q is completely used in increasing the internal energy U of the system. That is,
dU = q
Therefore, dh = q +Vdp
For an adiabatic process, q=0. Also, by substituting pV=nRT, we get the equation as,
dh = Rd(nT)
That is, hproduct – hreactant = R (nproduct x Tproduct – nreactant x Treactant)
2.3 Enthalpy equation for constant pressure combustion
For the constant pressure combustion, the enthalpy change is zero. This can be explained by the work done by the system to expand the volume of the system. Consider the system to have a heat loss factor, we can write the equation as
dh = 0
hproduct +hloss – hreactant = 0
Where, hloss = loss factor x Lower Heating Value of fuel
- Effect of equivalence ratio on AFT in constant volume system
In this section, the AFT calculation for combustion of methane fuel with varying equivalence ration is studied. For a rich mixture, we consider only CO as intermediary product and for a lean mixture, formation of O2. The equations are balanced in terms of equivalence ratio and the algorithm is generated to calculate AFT for a range of equivalence ratio from 0.1 to 2.
The results are plotted using matplotlib script commands. The same case is solved using Cantera and a comparison plot of the results is generated.
'''
Program to calculate ADT for different equivalence ratio
Stoichiometric Reaction
CH4 + 2(O2+ 3.76 N2) <=> CO2 + 2H2O + 7.52 N2+ 0O2+ 0CO
Lean mixture
CH4 + (2/phi)(O2+ 3.76 N2) <=> CO2 + 2H2O + (7.52/phi) N2 + (2*(1-phi)/phi)O2+ 0CO
Rich mixture
CH4 + (2/phi)(O2+ 3.76 N2) <=> (4/phi-3)CO2 + 2H2O + (7.52/phi) N2 + 0O2 + (4-4/phi)CO
Program written by Vignesh V
'''
import matplotlib.pyplot as plt
import math
#Function to evaluate enthalpy
R=8.314 #J/kmol-K
def h(T, coeff):
a1 = coeff[0]
a2 = coeff[1]
a3 = coeff[2]
a4 = coeff[3]
a5 = coeff[4]
a6 = coeff[5]
a7 = coeff[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
#Species coefficient from NASA polynomials
CH4_coeff = [5.14987613E+00,-1.36709788E-02,4.91800599E-05,-4.84743026E-08,1.66693956E-11,-1.02466476E+04,-4.64130376E+00]
O2react_coeff= [3.78245636E+00,-2.99673416E-03,9.84730201E-06,-9.68129509E-09, 3.24372837E-12,-1.06394356E+03,3.65767573E+00]
N2react_coeff=[0.03298677E+02,0.14082404E-02,-0.03963222E-04,0.05641515E-07,-0.02444854E-10,-0.10208999E+04,0.03950372E+02]
CO2_coeff=[3.85746029E+00,4.41437026E-03,-2.21481404E-06,5.23490188E-10,-4.72084164E-14,-4.87591660E+04,2.27163806E+00]
H2O_coeff=[3.03399249E+00,2.17691804E-03,-1.64072518E-07,-9.70419870E-11,1.68200992E-14,-3.00042971E+04,4.96677010E+00]
N2product_coeff=[0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04,0.05980528E+02]
O2product_coeff=[3.28253784E+00,1.48308754E-03,-7.57966669E-07,2.09470555E-10,-2.16717794E-14,-1.08845772E+03,5.45323129E+00]
CO_coeff=[2.71518561E+00,2.06252743E-03,-9.98825771E-07,2.30053008E-10,-2.03647716E-14,-1.41518724E+04,7.81868772E+00]
#Function to define the root finding problem
def f(T,mol_react,mol_product):
T_std= 298.15
h_CH4= h(T_std,CH4_coeff)
h_O2_react= h(T_std,O2react_coeff)
h_N2react= h(T_std,N2react_coeff)
h_CO2= h(T,CO2_coeff)
h_H2O= h(T,H2O_coeff)
h_N2product= h(T,N2product_coeff)
h_O2_product= h(T,O2product_coeff)
h_CO= h(T,CO_coeff)
N_react= mol_react[0]+mol_react[1]+mol_react[2]
N_product= mol_product[0]+mol_product[1]+mol_product[2]+mol_product[3]+mol_product[4]
h_reactant= mol_react[0]*h_CH4 + mol_react[1]*h_O2_react+ mol_react[2]*h_N2react
h_product= mol_product[0]*h_CO2+ mol_product[1]*h_H2O+ mol_product[2]*h_N2product+mol_product[3]*h_O2_product+mol_product[4]*h_CO
return h_product-h_reactant - R*(N_product*T-N_react*T_std)
def fprime(T,mol_react,mol_product):
return (f(T+1e-6,mol_react,mol_product)-f(T,mol_react,mol_product))/1e-6
T_guess=1500
tol=1e-3
x=0
#loop for phi from 0 to 2
import numpy
phi = numpy.zeros(50)
for x in range (0,50):
phi[x]=(2/50)*(x+1)
x=x+1
#loop for ADT for different phi
AFT = numpy.zeros(50)
i=0
for i in range (0,50):
#define moles for species based on equivalence ratio
if phi[i]==1:
mol_reactant=[1,2,7.52]
mol_products=[1,2,7.52,0,0]
elif phi[i]<1:
mol_reactant=[1,2/phi[i],7.52/phi[i]]
mol_products=[1,2,7.52/phi[i],((2/phi[i])-2),0]
elif phi[i]>1:
mol_reactant=[1,2/phi[i],7.52/phi[i]]
mol_products=[(4/phi[i])-3,2,7.52/phi[i],0,4-(4/phi[i])]
#Newton Rhapson method
while (abs(f(T_guess,mol_reactant,mol_products))>tol):
T_guess=T_guess-(f(T_guess,mol_reactant,mol_products)/fprime(T_guess,mol_reactant,mol_products))
AFT[i]= T_guess
i=i+1
'''
Code for AFT calculation using Cantera
'''
AFT_Cantera = numpy.zeros(50)
import cantera as ct
gas = ct.Solution('gri30.xml')
i=0
for i in range (0,50):
#define moles for species based on equivalence ratio
if phi[i]==1:
mol_reactant=[1,2,7.52]
mol_products=[1,2,7.52,0,0]
elif phi[i]<1:
mol_reactant=[1,2/phi[i],7.52/phi[i]]
mol_products=[1,2,7.52/phi[i],((2/phi[i])-2),0]
elif phi[i]>1:
mol_reactant=[1,2/phi[i],7.52/phi[i]]
mol_products=[(4/phi[i])-3,2,7.52/phi[i],0,4-(4/phi[i])]
gas.TPX = 298.15, 101325,{'CH4':mol_reactant[0],'O2':mol_reactant[1],'N2':mol_reactant[2]}
gas.equilibrate('UV','auto')
AFT_Cantera[i]=gas.T
i=i+1
plt.plot(phi,AFT,'*',color='blue',label='Python')
plt.plot(phi,AFT_Cantera,'.',color='red',label='Cantera')
plt.xlabel('Equivalence ratio')
plt.ylabel('Adiabatic flame temperature (in K)')
plt.title('Effect of equivalence ratio on AFT')
plt.legend()
plt.show()
It can be observed from the plot at lower equivalence ratio, the results obtained from Python and Cantera are very similar. However, at higher equivalence ratio, the results diverge. The values of AFT in the region close to stoichiometric region are overestimated. This can be attributed to the heat lost in heating the inert N2 resulting in formation of intermediary NOx compounds. Since the python is unable to capture this, the AFT values observed are more than that calculated using Cantera.
At rich mixtures, the Python program underestimates the Adiabatic Flame Temperature as compared to Cantera. The complex nature of combustion which involves several hundred intermediary species is accurately estimated in Cantera whereas the python program assumes presence of one additional species. This assumption is the reason for the deviation in the plots generated.
- Flame temperature for heat loss factor of 35% in constant pressure system
The flame temperature for a heat loss of 35% of Lower Heating Value (LHV) is calculated using programming in Python. For this case, we consider a stoichiometric reaction of methane combustion in constant pressure chamber.
'''
Program to calculate ADT for different hydrocarbons in a system with heat loss
Stoichiometric Reaction - Alkane
CxHy + (y/2)(O2+ 3.76 N2) <=> xCO2 + (y/2)H2O + (y/2)*3.76 N2
Lower Heating Value
LHV= 55 MJ/kg = 800000 J/mol
Program written by Vignesh V
'''
import matplotlib.pyplot as plt
import math
#Function to evaluate enthalpy
def h(T, coeff):
R=8.314 #J/kmol-K
a1 = coeff[0]
a2 = coeff[1]
a3 = coeff[2]
a4 = coeff[3]
a5 = coeff[4]
a6 = coeff[5]
a7 = coeff[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
#Species coefficient from NASA polynomials
CH4_coeff = [5.14987613E+00,-1.36709788E-02,4.91800599E-05,-4.84743026E-08,1.66693956E-11,-1.02466476E+04,-4.64130376E+00]
O2react_coeff= [3.78245636E+00,-2.99673416E-03,9.84730201E-06,-9.68129509E-09, 3.24372837E-12,-1.06394356E+03,3.65767573E+00]
N2react_coeff=[0.03298677E+02,0.14082404E-02,-0.03963222E-04,0.05641515E-07,-0.02444854E-10,-0.10208999E+04,0.03950372E+02]
CO2_coeff=[3.85746029E+00,4.41437026E-03,-2.21481404E-06,5.23490188E-10,-4.72084164E-14,-4.87591660E+04,2.27163806E+00]
H2O_coeff=[3.03399249E+00,2.17691804E-03,-1.64072518E-07,-9.70419870E-11,1.68200992E-14,-3.00042971E+04,4.96677010E+00]
N2product_coeff=[0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04,0.05980528E+02]
#Function to define the root finding problem
def f(T,mol_react,mol_product,h_loss):
T_std= 298.15
LHV = 800000 #J/mol-K
h_CH4= h(T_std,CH4_coeff)
h_O2_react= h(T_std,O2react_coeff)
h_N2react= h(T_std,N2react_coeff)
h_CO2= h(T,CO2_coeff)
h_H2O= h(T,H2O_coeff)
h_N2product= h(T,N2product_coeff)
h_reactant= mol_react[0]*h_CH4 + mol_react[1]*h_O2_react+ mol_react[2]*h_N2react
h_product= mol_product[0]*h_CO2+ mol_product[1]*h_H2O+ mol_product[2]*h_N2product
return h_product-h_reactant + h_loss*LHV
def fprime(T,mol_react,mol_product,h_loss):
return (f(T+1e-6,mol_react,mol_product,h_loss)-f(T,mol_react,mol_product,h_loss))/1e-6
#AFT for methane with heat loss of 35%
T_guess=1500
tol=1e-3
h_loss = 0.35
x=1
y=2*x+2
mol_reactant= [1, y/2,y/2*3.76]
mol_products= [x,y/2,y/2*3.76]
#Newton Rhapson method
while (abs(f(T_guess,mol_reactant,mol_products,h_loss))>tol):
T_guess=T_guess-(f(T_guess,mol_reactant,mol_products,h_loss)/fprime(T_guess,mol_reactant,mol_products,h_loss))
AFT_CH4_loss35 = T_guess
print(AFT_CH4_loss35)
It is observed that Flame temperature obtained is 1680 K. This is much lower than the AFT obtained in previous case, 2325 K. The heat lost from the system results in less enthalpy of products which explains the drop in temperature.
- Flame temperature for alkane, alkene and alkyne hydrocarbons
This section involves calculating the flame temperature for ethane, ethane and ethyne for a heat loss factor of 35%. The program algorithm involves defining moles of reactants and moles of products for each case. Using the NASA polynomials and the Newton Raphson method, Flame temperature is calculated at thermal equilibrium.
'''
Comparison of AFT for C2H8, C2H6, C2H4
C2H6 + 3.5*(O2 + 3.76 N2) = 2CO2 + 3H2O + (3.5*3.76) N2
C2H4 + 3*(O2 + 3.76 N2) = 2CO2 + 2H2O + (3*3.76) N2
C2H2 + 2.5*(O2 + 3.76 N2) = 2CO2 + H2O + (2.5*3.76) N2
'''
import matplotlib.pyplot as plt
import math
#Function to evaluate enthalpy
def h(T, coeff):
R=8.314 #J/kmol-K
a1 = coeff[0]
a2 = coeff[1]
a3 = coeff[2]
a4 = coeff[3]
a5 = coeff[4]
a6 = coeff[5]
a7 = coeff[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
#Species coefficient from NASA polynomials
C2H6_coeff = [4.29142492E+00,-5.50154270E-03,5.99438288E-05,-7.08466285E-08,2.68685771E-11,-1.15222055E+04,2.66682316E+00]
C2H4_coeff = [3.95920148E+00,-7.57052247E-03,5.70990292E-05,-6.91588753E-08,2.69884373E-11,5.08977593E+03,4.09733096E+00]
C2H2_coeff = [8.08681094E-01,2.33615629E-02,-3.55171815E-05,2.80152437E-08,-8.50072974E-12,2.64289807E+04,1.39397051E+01]
O2react_coeff= [3.78245636E+00,-2.99673416E-03,9.84730201E-06,-9.68129509E-09, 3.24372837E-12,-1.06394356E+03,3.65767573E+00]
N2react_coeff=[0.03298677E+02,0.14082404E-02,-0.03963222E-04,0.05641515E-07,-0.02444854E-10,-0.10208999E+04,0.03950372E+02]
CO2_coeff=[3.85746029E+00,4.41437026E-03,-2.21481404E-06,5.23490188E-10,-4.72084164E-14,-4.87591660E+04,2.27163806E+00]
H2O_coeff=[3.03399249E+00,2.17691804E-03,-1.64072518E-07,-9.70419870E-11,1.68200992E-14,-3.00042971E+04,4.96677010E+00]
N2product_coeff=[0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04,0.05980528E+02]
#Function to define the root finding problem
def f(T,h_loss,x):
T_std= 298.15
LHV = 800000 #J/mol-K
h_C2H6= h(T_std,C2H6_coeff)
h_C2H4= h(T_std,C2H4_coeff)
h_C2H2= h(T_std,C2H2_coeff)
h_O2_react= h(T_std,O2react_coeff)
h_N2react= h(T_std,N2react_coeff)
h_CO2= h(T,CO2_coeff)
h_H2O= h(T,H2O_coeff)
h_N2product= h(T,N2product_coeff)
if x==0:
h_fuel = h_C2H6
mol_reactant= [1,3.5,3.5*3.76]
mol_products= [2,3,3.5*3.76]
elif x==1:
h_fuel = h_C2H4
mol_reactant= [1,3,3*3.76]
mol_products= [2,2,3*3.76]
elif x==2:
h_fuel = h_C2H2
mol_reactant= [1,2.5,2.5*3.76]
mol_products= [2,1,2.5*3.76]
h_reactant= mol_reactant[0]*h_fuel + mol_reactant[1]*h_O2_react+ mol_reactant[2]*h_N2react
h_product= mol_products[0]*h_CO2+ mol_products[1]*h_H2O+ mol_products[2]*h_N2product
return h_product-h_reactant + h_loss*LHV
def fprime(T,h_loss,x):
return (f(T+1e-6,h_loss,x)-f(T,h_loss,x))/1e-6
#AFT for ethane, ethene and ethyne
T_guess=1500
tol=1e-3
h_loss = 0.35
import numpy
AFT= numpy.zeros(3)
i=0
for i in range (0,3):
#Newton Rhapson method
while (abs(f(T_guess,h_loss,i))>tol):
T_guess=T_guess-(f(T_guess,h_loss,i)/fprime(T_guess,h_loss,i))
AFT[i] = T_guess
i=i+1
print(AFT)
plt.plot(['C2H6'],AFT[0],'*')
plt.plot(['C2H4'],AFT[1],'*')
plt.plot(['C2H2'],AFT[2],'*')
plt.title('AFT for ethane,ethene,ethyne')
plt.ylabel('Adiabatic flame temperature (in K)')
plt.show()From the plot, it can be observed that with the flame temperature increases with an increase in number of carbon bonds. An explanation for this can be the reduction of Hydrogen atoms. Due to this, less number of H2O molecules is formed for alkyne as compared to alkene and alkanes. The result of this is the high flame temperature at thermal equilibrium for combustion of alkyne hydrocarbon compared to alkene and alkane.
- Flame temperature trend for increase in number of carbon atoms
The section determines the flame temperature for three hydrocarbons- methane, ethane and propane with stoichiometric combustion in a constant pressure system with zero heat loss. The resultant adiabatic flame temperature are plotted to perform a comparative study.
'''
Program to calculate ADT for different hydrocarbons in a system with heat loss
Stoichiometric Reaction - Alkane
CxHy + (x+y/4)(O2+ 3.76 N2) <=> xCO2 + (y/2)H2O + (x+y/4)*3.76 N2
Lower Heating Value
LHV= 55 MJ/kg = 800000 J/mol
Program written by Vignesh V
'''
import matplotlib.pyplot as plt
import math
#Function to evaluate enthalpy
def h(T, coeff):
R=8.314 #J/kmol-K
a1 = coeff[0]
a2 = coeff[1]
a3 = coeff[2]
a4 = coeff[3]
a5 = coeff[4]
a6 = coeff[5]
a7 = coeff[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
#Species coefficient from NASA polynomials
CH4_coeff = [5.14987613E+00,-1.36709788E-02,4.91800599E-05,-4.84743026E-08,1.66693956E-11,-1.02466476E+04,-4.64130376E+00]
C2H6_coeff = [4.29142492E+00,-5.50154270E-03,5.99438288E-05,-7.08466285E-08,2.68685771E-11,-1.15222055E+04,2.66682316E+00]
C3H8_coeff = [0.93355381E+00,0.26424579E-01,0.61059727E-05,-0.21977499E-07,0.95149253E-11,-0.13958520E+05,0.19201691E+02]
O2react_coeff= [3.78245636E+00,-2.99673416E-03,9.84730201E-06,-9.68129509E-09, 3.24372837E-12,-1.06394356E+03,3.65767573E+00]
N2react_coeff=[0.03298677E+02,0.14082404E-02,-0.03963222E-04,0.05641515E-07,-0.02444854E-10,-0.10208999E+04,0.03950372E+02]
CO2_coeff=[3.85746029E+00,4.41437026E-03,-2.21481404E-06,5.23490188E-10,-4.72084164E-14,-4.87591660E+04,2.27163806E+00]
H2O_coeff=[3.03399249E+00,2.17691804E-03,-1.64072518E-07,-9.70419870E-11,1.68200992E-14,-3.00042971E+04,4.96677010E+00]
N2product_coeff=[0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04,0.05980528E+02]
#Function to define the root finding problem
def f(T,mol_react,mol_product,h_loss,x):
T_std= 298.15
LHV = 800000 #J/mol-K
h_CH4= h(T_std,CH4_coeff)
h_C2H6= h(T_std,C2H6_coeff)
h_C3H8= h(T_std,C3H8_coeff)
h_O2_react= h(T_std,O2react_coeff)
h_N2react= h(T_std,N2react_coeff)
h_CO2= h(T,CO2_coeff)
h_H2O= h(T,H2O_coeff)
h_N2product= h(T,N2product_coeff)
if x==1:
h_fuel = h_CH4
elif x==2:
h_fuel = h_C2H6
elif x==3:
h_fuel = h_C3H8
h_reactant= mol_react[0]*h_fuel + mol_react[1]*h_O2_react+ mol_react[2]*h_N2react
h_product= mol_product[0]*h_CO2+ mol_product[1]*h_H2O+ mol_product[2]*h_N2product
return h_product-h_reactant + h_loss*LHV
def fprime(T,mol_react,mol_product,h_loss,x):
return (f(T+1e-6,mol_react,mol_product,h_loss,x)-f(T,mol_react,mol_product,h_loss,x))/1e-6
#AFT for carbon atoms from 1 to 3 with no loss
T_guess=1500
tol=1e-3
h_loss = 0
import numpy
AFT= numpy.zeros(3)
i=0
x=[1,2,3]
for i in range (0,3):
y=2*x[i]+2
mol_reactant= [1, x[i]+y/4,(x[i]+y/4)*3.76]
mol_products= [x[i],y/2,(x[i]+y/4)*3.76]
#Newton Rhapson method
while (abs(f(T_guess,mol_reactant,mol_products,h_loss,x[i]))>tol):
T_guess=T_guess-(f(T_guess,mol_reactant,mol_products,h_loss,x[i])/fprime(T_guess,mol_reactant,mol_products,h_loss,x[i]))
AFT[i] = T_guess
i=i+1
print(AFT)
plt.plot([1],AFT[0],'*',label='CH4')
plt.plot([2],AFT[1],'*',label='C2H6')
plt.plot([3],AFT[2],'*',label='C3H8')
plt.title('AFT trend for hydrocarbon with increasing C atoms')
plt.ylabel('Adiabatic flame temperature (in K)')
plt.xlabel('Number of carbon atoms')
plt.legend()
plt.show()From the plot, it can be concluded that AFT increases with an increase in carbon atoms. With an increase in number of carbon atoms, mass of fuel is more which combines with more number of air molecules. That is, the combustion mixture is enriched for the same volume of air when more number of carbon atoms are present. The energy released due to this increases the product enthalpies and the net enthalpy difference increases. In turn, AFT observed is more than of lesser number of carbon atoms in case of hydrocarbon fuel combustion.
- Conclusion
The following are the conclusions of the current study:
- Adiabatic Flame Temperature plot for a combustion process with variation in equivalence ratio is more accurate when programmed using Cantera
- Heat loss from the combustion system results in a drop in the flame temperature at thermal equilibrium
- The flame temperature is highest for alkyne amongst the alkane, alkene and alkyne hydrocarbons
- With an increase in number of carbon atoms in a hydrocarbons, the flame temperature increases.
The link to program scripts are below: