BlogFor BusinessTalk to a counsellor
CoursesIIT ProgramsPlacementsLog in
BlogFor Business
Talk to a counsellor
Challenge

Adiabatic Flame Temperature calculation

MMahesh U13 Mar 2021

Aim:

  • To calculate Adiabatic flame temperature for a given fuel using python programming.

Adiabatic Flame Temperature is a benchmark temperature for any combustion system. AFT is the highest temperature acheived when a fuel reacts with an oxidizer and approaches equilibrium by means of an adiabatic and isobaric process. In a combustion chamber, the highest temperature the material can be exposed to is limited by metallurgical considerations. Hence, AFT is an important consideration in the design of combustion chambers, nozzles and gas turbines.

Objectives:

  • To calculate AFT for methane at constant volume with a range of equivalence ratios and compare the results obtained from python and cantera.
  • To calculate AFT for methane with heat loss of 35% from the chamber.
  • To calculate AFT for an alkane,alkene and alkyne and to compare the results.
  • To calculate AFT for different hydrocarbons with varying carbon atoms(1,2 &3) and to compare the results.

Case-1:To calculate AFT for methane at constant volume with a range of equivalence ratios using python progamming and cantera.

The equivalence ratio, 'phi' is the ratio of actual air/fuel ratio to stoichiometric air/fuel ratio. At an equvalence ratio of one, the combustion is said to be stoichiometric. If it is less than one, the combustion is lean with excess air, and if it is greater than one, the combustion is rich with incomplete combustion.
Here the fuel used is methane. Stoichiometric combustion of methane is given by;

CH4 + 2(O2 + 3.76N2) => CO2 + 2H20 +7.52N2

If the equivalence ratio, `phi` is less than one or in other words, if the combustion is lean. In lean combustion, air is in excess and hence there will be extra oxygen in the products.

CH4 + (2/phi)(O2 + 3.77N2) => CO2 + 2H20 + (7.52/phi)N2 + (2/phi-2)O2

If `phi` is greater than one, or ifthe combustion is rich(excess fuel). The unburnt fuel molecules result in the formation of carbon monoxide.

CH4 + (2/phi)(O2+3.77N2) => (4/phi-3)CO2 + 2H20 + (7.52/phi)N2 + (4-4/phi)CO

The adiabatic flame temperature had to be calculated using two methods,

1) Python programming: The enthalpies of all participating reactants and products are calculated using NASA polynomials.

`H =( a_1+(a_2/2)*T+(a_3/3)*T^2+(a_4/4)*T^3+(a_5/5)*T^4+(a_6/T))*R*T`

Where,

H = enthalpy

`a_1,a_2,...,a_6` = NASA polynomials

T = temperature of the reactants/products.

The reactants are considered to be at standard temperature of T = 298.15 K. In case of reactants, temperature is the adiabatic flame temperature and is unknown.

To calculate AFT, Newton Raphson's method of approximation of the root of a real values function is used.

`T_(AFT) = T_(guess)- alpha*(f(x)/(f'(x)))`

Since this is an iterative technique, an initial value for AFT had to be guessed. To calculate AFT, `f(x)` had to be set as change in enthalpy(`DeltaH`). But here, the adiabatic flame temperature was be calculated in a constant volume chamber. In a constant volume chamber, workdone and heat transfer are zero.

Hence, change in internal energy, `U_p-U_r = 0`

Or, `(H_p - N_p*R*T_p) - (H_r-N_r*R*T_r)` = 0

The Newton Raphson's method is used to approximate the root of above equation.

2) Cantera: Cantera Python module is a tool for problems involving chemical kinetics, thermodynamics, and transport processes. In cantera, the temperature,pressure and molefractions of the reactants had to be set to perform the combustion.

Using both these methods, AFT for a range of equivalence ratios were calculated.

import matplotlib.pyplot as plt
import cantera as ct
import numpy as np
'''
stoichiometric combustion
CH4 + 2(O2 + 3.76N2) => CO2 + 2H20 +7.52N2

Lean combustion(excess air)
CH4 + (2/phi)(O2 + 3.77N2) => CO2 + 2H20 + (7.52/phi)N2 + (2/phi-2)O2

Rich combustion(Lack of air)
CH4 + (2/phi)(O2+3.77N2) => (4/phi-3)CO2 + 2H20 + (7.52/phi)N2 + (4-4/phi)CO

Equation to find Enthalpy of reactants/products using NASA polynomials
H/(R*T) = a1 + a2*T/2 + a3*(T^2)/3 + a4*(T^3)/4 + a5*(T^4)/5 + a6/T

'''
#NASA polynomials for required reactants/products from Berkeley gri-mech file

CH4_cff = [5.14987613E+00, -1.36709788E-02, 4.91800599E-05, -4.84743026E-08, 1.66693956E-11, -1.02466476E+04]
O2_cff_r = [3.78245636E+00, -2.99673416E-03, 9.84730201E-06, -9.68129509E-09, 3.24372837E-12, -1.06394356E+03]
N2_cff_r = [0.03298677E+02, 0.14082404E-02, -0.03963222E-04, 0.05641515E-07, -0.02444854E-10, -0.10208999E+04]
CO2_cff = [3.85746029E+00, 4.41437026E-03, -2.21481404E-06, 5.23490188E-10, -4.72084164E-14, -4.87591660E+04]
H2O_cff = [3.03399249E+00, 2.17691804E-03, -1.64072518E-07, -9.70419870E-11, 1.68200992E-14, -3.00042971E+04]
N2_cff_p= [0.02926640E+02, 0.14879768E-02,-0.05684760E-05, 0.10097038E-09, -0.06753351E-13, -0.09227977E+04]
O2_cff_p =[3.28253784E+00, 1.48308754E-03, -7.57966669E-07, 2.09470555E-10, -2.16717794E-14, -1.08845772E+03]
CO_cff = [2.71518561E+00, 2.06252743E-03, -9.98825771E-07, 2.30053008E-10, -2.03647716E-14, -1.41518724E+04]

def h(T, cff): #Function to calculate enthalpy
	a1 = cff[0]
	a2 = cff[1]
	a3 = cff[2]
	a4 = cff[3]
	a5 = cff[4]
	a6 = cff[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

R = 8.314 #J/molK Gas constant
T_std = 298.15 #K Temperature at standard condition

def f(T, mol_p,mol_r): #fuction to calculate difference in enthapies of products and reactants
	h_CH4 = h(T_std, CH4_cff)
	h_O2_r = h(T_std, O2_cff_r)
	h_N2_r = h(T_std,N2_cff_r)

	h_reactants = mol_r[0]*h_CH4 + mol_p[1]*h_O2_r + mol_r[2]*h_N2_r
	N_r = mol_r[0] + mol_r[1] + mol_r[2]

	h_CO2 = h(T, CO2_cff)
	h_H20 = h(T, H2O_cff)
	h_N2_p = h(T,N2_cff_p)
	h_O2_p = h(T, O2_cff_p)
	h_CO = h(T, CO_cff)

	h_products = mol_p[0]*h_CO2 + mol_p[1]*h_H20 + mol_p[2]*h_N2_p + mol_p[3]*h_CO +mol_p[4]*h_O2_p
	N_p = mol_p[0] + mol_p[1] + mol_p[2] + mol_p[3] + mol_p[4]

	return h_products- (R*N_p*T_guess) - h_reactants - (R*N_r*T_std)


def fprime(T,mol_p,mol_r): #forward differencing
	return (f(T+1e-6,mol_p,mol_r) - f(T,mol_p,mol_r))/1e-6


T_guess = 1500 #K initial guess value of AFT
alpha = 0.2 #relaxation factor
tol = 1e-3 #tolerance
AFT= []
phi= np.zeros(20) #creating an array for equivalence ratios
for i in range(0,20):
	phi[i] = (i*0.1)+0.1

	if phi[i] <= 1:
		mol_r = [1,2/phi[i],7.54/phi[i]]
		mol_p = [1,2,7.54/phi[i],0,2/phi[i]-2]
	elif phi[i] > 1:
		mol_r = [1,2/phi[i],7.54/phi[i]]
		mol_p = [4/phi[i]-3, 2, 7.54/phi[i], 4-4/phi[i], 0]

  #Newton Raphson method of root finding
	while abs(f(T_guess,mol_p,mol_r)) > tol:
		T_guess = T_guess - alpha*(f(T_guess,mol_p,mol_r)/fprime(T_guess,mol_p,mol_r))
	AFT.append(T_guess)
	i+=1
plt.plot(phi,AFT, label = 'Python')
plt.plot(phi[9],AFT[9], 'b*')
plt.text(phi[9],AFT[9],str(AFT [9]) +' K')
plt.grid('on')


#AFT using Cantera
gas = ct.Solution('gri30.cti')
AFT_ct= []
for i in range(0,20):
	phi[i] = (i*0.1)+0.1
	mol_r = [1,2/phi[i],7.54/phi[i]]
	gas.TPX = 298.15,101325, { 
		'CH4':mol_r[0], 'O2':mol_r[1], 'N2':mol_r[2]
	}#setting the temperature,pressure and mole fractions of reactants.
	gas.equilibrate('UV', 'auto') #constant volume combustion.
	AFT_ct.append(gas.T)
	i+=1

plt.plot(phi,AFT_ct,'r', label = 'Cantera')
plt.plot(phi[10],AFT_ct[10], 'rs')
plt.text(phi[10],AFT_ct[10],str(AFT_ct[10])+' K')
plt.title('Equivalence ratio Vs AFT', fontweight = 'bold')
plt.xlabel('Equivalence ratio')
plt.ylabel('Adiabatic flame temperature(K)')
plt.legend()
plt.grid('on')
print('Peak adiabatic flame temperature using python= ' + str(AFT [9])+ ' K'+ 'nPeak adiabatic flame temperature using Cantera= ' + str(AFT_ct[10])+ ' K')
plt.show()

Results and observations:

The adiabatic flame temperatures for equivalence ratios varying from 0.1 to 2 were obtained as follows:

Plot.1 Equivalence ratio vs AFT

  • The peak temperature for both cases were obtained at an equivalence ratio of roughly one. Equivalence ratio of one indicates stochiometric combustion.
  • A significant difference in the peak temperature was observed while comparing the results. The AFT obtained from python programming was higher than that of cantera.
  • The difference in AFT can be attributed to the fact that, the python code did not take into consideration the numerous products that would be formed apart from O2,CO2 and N2 in actual combustion of Methane.
  • Although the other products maybe in smaller concentrations, the presence of these products will definitely have an effect on the adiabatic flame temperature. In the python method of calculating AFT, this effect is not considered and hence the results are inaccurate.
  • The adiabatic flame temperature is lower for cantera method because it takes into consideration all the product concentrations and their effects on flame temperature before calculating the AFT and hence is the accurate result.


Case-2:To calculate AFT for methane with heat loss of 35% from the chamber.

For this case, adiabatic temperature had to be calculated with an assumption that heatloss of 35% takes place from the combustion chamber. The heat loss was calculated interms of the fuel's lower heating value. LHV is the amount of heat released by the fuel. It can be calculated from the difference in enthalpies of reactants and products at 298.15 K. The heat loss factor was then multiplied by the LHV to determine heat loss in J/mol-K. 

import matplotlib.pyplot as plt
import numpy as np
'''
stoichiometric combustion
CH4 + 2(O2 + 3.76N2) => CO2 + 2H20 +7.52N2

Equation to find Enthalpy of reactants/products using NASA polynomials
H/(R*T) = a1 + a2*T/2 + a3*(T^2)/3 + a4*(T^3)/4 + a5*(T^4)/5 + a6/T

heat loss = 35% = (Lower heating value*heat loss factor)
'''
#NASA polynomials for required reactants/products from Berkeley gri-mech file

CH4_cff_l = [5.14987613E+00, -1.36709788E-02, 4.91800599E-05, -4.84743026E-08, 1.66693956E-11, -1.02466476E+04]
O2_cff_l = [3.78245636E+00, -2.99673416E-03, 9.84730201E-06, -9.68129509E-09, 3.24372837E-12, -1.06394356E+03]
N2_cff_l = [0.03298677E+02, 0.14082404E-02, -0.03963222E-04, 0.05641515E-07, -0.02444854E-10, -0.10208999E+04]
CO2_cff_h = [3.85746029E+00, 4.41437026E-03, -2.21481404E-06, 5.23490188E-10, -4.72084164E-14, -4.87591660E+04]
H2O_cff_h = [3.03399249E+00, 2.17691804E-03, -1.64072518E-07, -9.70419870E-11, 1.68200992E-14, -3.00042971E+04]
N2_cff_h= [0.02926640E+02, 0.14879768E-02,-0.05684760E-05, 0.10097038E-09, -0.06753351E-13, -0.09227977E+04]
CO2_cff_l = [2.35677352E+00, 8.98459677E-03, -7.12356269E-06, 2.45919022E-09, -1.43699548E-13, - 4.83719697E+04]
H2O_cff_l = [4.19864056E+00, -2.03643410E-03, 6.52040211E-06, -5.48797062E-09, 1.77197817E-12, -3.02937267E+04]


def h(T, cff):  #Function to calculate enthalpy
	a1 = cff[0]
	a2 = cff[1]
	a3 = cff[2]
	a4 = cff[3]
	a5 = cff[4]
	a6 = cff[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

R = 8.314 #J/molK Gas constant
T_std = 298.15 #K Temperature at standard condition
#To calculate LHV in J/mol-K, LHV = enthalpy of reactants - Enthalpy of products (both at standard temperatures)
LHV = (h(T_std,CH4_cff_l) + 2*(h(T_std,O2_cff_l) + 3.77*h(T_std,N2_cff_l))) - (h(T_std, CO2_cff_l) + 2*h(T_std, H2O_cff_l) + 7.54*h(T_std, N2_cff_l))

def f(T, heatloss): #fuction to calculate difference in enthapies of products and reactants

	h_CH4 = h(T_std, CH4_cff_l)
	h_O2_r = h(T_std, O2_cff_l)
	h_N2_r = h(T_std,N2_cff_l)
	h_reactants = h_CH4 + 2*(h_O2_r + 3.77*h_N2_r)

	h_CO2 = h(T, CO2_cff_h)
	h_H20 = h(T, H2O_cff_h)
	h_N2_p = h(T,N2_cff_h)
	h_products = h_CO2 + 2*h_H20 + 7.54*h_N2_p

	return h_products - h_reactants + heatloss

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


T_guess = 1500 #K initial guess value of AFT
alpha = 0.2 #relaxation factor
tol = 1e-3 #tolerance
h_loss = np.zeros(20)
AFT= []
for i in range(0,20):
	h_loss[i] = (i*0.05) +.05
	heatloss = LHV*h_loss[i]
  #Newton Raphson method of root finding
	while abs(f(T_guess,heatloss)) > tol:
		T_guess = T_guess - alpha*(f(T_guess,heatloss)/fprime(T_guess,heatloss))
	AFT.append(T_guess)
	i+=1
plt.plot(h_loss, AFT, '-b*')
plt.plot(h_loss[6], AFT[6],'ro', label = '0.35 Heat loss factor')
plt.text(0.38, 1676.3,'1676.3 K', fontsize='8')
plt.xlim(0,1)
plt.title('Heat Loss factor Vs AFT', fontweight = 'bold')
plt.xlabel('Heat loss factor')
plt.ylabel('Adiabatic flame temperature(K)')
plt.grid('on')
plt.legend()
print('Adiabatic Flame temperature with 35% Heat loss = ' +str(AFT[6]) +' K')
plt.show()

Results and observations:

Plot.2 Heat loss factor Vs AFT

The adiabatic flame temperature with heat loss of 35% was obtained as 1676.3 K. As expected, the flame temperature can be seen as decreasing with increase in heat loss.


Case-3: To calculate AFT for an alkane,alkene and alkyne.

The adiabatic flame temperatures of hydrocarbons with two carbon atom chains were found using cantera.

import matplotlib.pyplot as plt
import cantera as ct
'''
stoichiometric combustion of: 
Ethane
C2H6 + 3.5(O2+3.76N2) => 2CO2 + 3H2O + 13.16N2
Ethene
C2H4 + 3(O2+3.76N2) => 2CO2 + 2H20 + 11.28N2
Ethyne
C2H2 + 2.5(O2+3.76N2) => 2CO2 + 1H20 + 9.4N2 
'''
#Calculating AFT using Cantera

gas = ct.Solution('gri30.cti') 
H_num = ['C2H6','C2H4','C2H2'] #Hydrocarbon variants
it=0
AFT=[]
for i in H_num:
	if i== H_num[0]:
		gas.TPX = 298.15,101325,{'C2H6':1,'O2':3.5,'N2':13.16}
		gas.equilibrate('HP','auto')
	elif i== H_num[1]:
		gas.TPX = 298.15,101325,{'C2H4':1,'O2':3,'N2':11.28}
		gas.equilibrate('HP','auto')
	elif i== H_num[2]:	
		gas.TPX = 298.15,101325,{'C2H2':1,'O2':2.5,'N2':9.4}
		gas.equilibrate('HP','auto')
	AFT.append(gas.T)
	it+=1
	
plt.bar([1,2,3],[AFT[0], AFT[1], AFT[2]], width=0.4)
plt.title('AFT for Ethane/Ethene/Ethyne', fontweight = 'bold')
plt.ylabel('Adiabatic flame temperature (K)')
plt.xticks([1,2,3], ['Ethane','Ethene','Ethyne'])
print('Adiabatic Flame temperature obtained: nEthane = '+str(AFT[0])+' K'+'nEthene = '+str(AFT[1])+' K'+'nEthyne = '+str(AFT[2])+' K' )
plt.show()

Results and observations:

Plot.3 AFT variation for ethane/ethene/ethyne.

Hydrocarbon Adiabatic flame temperature (K)
Ethane 2258.73
Ethene 2368.64
Ethyne 2539.67

The adiabatic flame temperature increased for hydrocarbons with fewer hydrogen atoms.

`AFT_(ALKANE)< AFT _(ALKENE) < AFT_( ALKYNE)`

Case-4: To calculate AFT with respect to number of carbon atom chains:

The adiabatic flame temperature was calculated for methane, ethane and propane using cantera.

import matplotlib.pyplot as plt
import cantera as ct
'''
stoichiometric combustion of: 
Methane
CH4 + 2(O2 + 3.76N2) => CO2 + 2H20 +7.52N2 
Ethane
C2H6 + 3.5(O2+3.76N2) => 2CO2 + 3H2O + 13.16N2 
Propane
C3H8 + 5(O2 + 3.76N2) => 3CO2 + 4H20 + 18.8N2 
'''

#Calculating AFT using Cantera
gas = ct.Solution('gri30.cti')
H_num = ['CH4','C2H6','C3H8'] #Hydrocarbon variants
it=0
AFT=[]
for i in H_num:
	if i== H_num[0]:
		gas.TPX = 298.15,101325,{'CH4':1,'O2':2,'N2':7.52}
		gas.equilibrate('HP','auto')
	elif i== H_num[1]:
		gas.TPX = 298.15,101325,{'C2H6':1,'O2':3.5,'N2':13.16}
		gas.equilibrate('HP','auto')
	elif i== H_num[2]:	
		gas.TPX = 298.15,101325,{'C3H8':1,'O2':5,'N2':18.8}
		gas.equilibrate('HP','auto')
	AFT.append(gas.T)
	it+=1
	
plt.bar([1,2,3],[AFT[0], AFT[1], AFT[2]], width=0.4)
plt.title('AFT for Methane/Ethane/Propane', fontweight = 'bold')
plt.ylabel('Adiabatic flame temperature (K)')
plt.xticks([1,2,3], ['Methane','Ethane','Propane'])
print('Adiabatic Flame temperature obtained: nMethane = '+str(AFT[0])+' K'+'nEthane = '+str(AFT[1])+' K'+'nPropane = '+str(AFT[2])+' K' )
plt.show()

Results and observations:

Hydrocarbon Adiabatic flame temperature(K)
Methane 2224.61
Ethane 2258.73
Propane 2265.70

The adiabatic flame temperature increased for increase in carbon atom chains.

`AFT_(METHANE)< AFT _(ETHANE) < AFT_(PROPANE)`

Conclusion:

Successfully calculated adiabatic flame temperature for multiple cases using python coding and cantera. Python coding approach of solving combustion was found to be very inefficient in obtaining accurate results. Cantera gave accurate results with fewer lines of codes.

Your portfolio, next

Ship work like this
in your own course.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.

Talk to a career counsellor

Free 20-minute call. We'll map the right course to where you want to work.

By submitting you agree to our Terms and Privacy Policy.