Challenge

Week 3 - Adiabatic Flame Temperature calculation

Adiabatic flame temperature variation with Equivalence ratio (In constant volume-chamber) :

AIM: To calculate the Adiabatic Flame Temperature.

OBJECTIVES:

1) Plot the effect of equivalence ratio on the final adiabatic flame temperature and further compare the results with Cantera.

2) Find out the Adiabatic Flame Temperature of Methane Air Combustion in constant pressure with heat loss.

3) To examine the Adiabatic Flame Temperature of Alk-ane, Alk-ene, Alk-yne with no heat loss.

4) Find out the Adiabatic flame temperature variation with respect to a number of carbon atom chains.

import cantera as cnt
import matplotlib.pyplot as plt
import math
import numpy as np

# CH4 + 2(O2 + 3.76N2) = CO2 + 2H2O + 7.52N2
def h(T,co_effs):
	
	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

R = 8.314 #j/mol-k
CH4_coeff_lw = [5.14987613E+00,-1.36709788E-02,4.91800599E-05,-4.84743026E-08,1.66693956E-11,-1.02466476E+04,-4.64130376E+00]
O2_coeff_lw = [3.78245636E+00,-2.99673416E-03,9.84730201E-06,-9.68129509E-09,3.24372837E-12,-1.06394356E+03,3.65767573E+00]
N2_coeff_lw = [0.03298677E+02,0.14082404E-02,-0.03963222E-04,0.05641515E-07,-0.02444854E-10,-0.10208999E+04,0.03950372E+02]
CO2_coeff_hg = [3.85746029E+00,4.41437026E-03,-2.21481404E-06,5.23490188E-10,-4.72084164E-14,-4.87591660E+04,2.27163806E+00]
N2_coeff_hg = [0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04,0.05980528E+02]
H2O_coeff_hg = [3.03399249E+00,2.17691804E-03,-1.64072518E-07,-9.70419870E-11,1.68200992E-14,-3.00042971E+04,4.96677010E+00]
O2_coeff_hg = [3.28253784E+00,1.48308754E-03,-7.57966669E-07,2.09470555E-10,-2.16717794E-14,-1.08845772E+03,5.45323129E+00]
CO_coeff_hg = [2.71518561E+00,2.06252743E-03,-9.98825771E-07,2.30053008E-10,-2.03647716E-14,-1.41518724E+04,7.81868772E+00]


def f(T,psi):

	H_CO2_p = h(T,CO2_coeff_hg)
	H_H2O_p = h(T,H2O_coeff_hg)
	H_N2_p = h(T, N2_coeff_hg)
	H_O2_p = h(T, O2_coeff_hg)
	H_CO_p = h(T,CO_coeff_hg)

	T_std = 298.15 # for reactants 
	H_CH4_r = h(T_std,CH4_coeff_lw)
	H_O2_r = h(T_std,O2_coeff_lw)
	H_N2_r = h(T_std,N2_coeff_lw)


	# for lean mixture psi < 1
	if psi <= 1:

		H_reactants = H_CH4_r+(2/psi)*(H_O2_r+3.76*H_N2_r)
		H_products = H_CO2_p+2*H_H2O_p+(7.52/psi)*H_N2_p+((2/psi)-2)*H_O2_p

		N_reactants = 1+(2/psi)*(1+3.76)
		N_productants = 1+2+(7.52/psi)+((2/psi)-2)

	else: # for rich mixture 

		H_reactants = H_CH4_r+(2/psi)*(H_O2_r+(3.76*H_N2_r))
		H_products = ((4/psi)-3)*H_CO2_p+2*H_H2O_p+(7.52/psi)*H_N2_p+(4-(4/psi))*H_CO_p

		N_reactants = 1+(2/psi)*(1+3.76)
		N_productants = ((4/psi)-3)+2+(7.52/psi)+(4-(4/psi))
	
	return(-H_products+H_reactants-R*(-N_productants*T+N_reactants*T_std))

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

# NRM


alfa = 1
ct = 0
T_guess = 1500

x = np.linspace(0.2,2,50)
T_ad = []
for psi in x:
	tol = 1e-6
	T_guess = 1500
	while (abs(f(T_guess,psi))>tol):
		T_guess = T_guess-alfa*((f(T_guess,psi))/(fprime(T_guess,psi)))
		#print(T_guess)
		ct = ct+1
	T_ad.append(T_guess)
plt.plot(x,T_ad,color = 'red',label='Python')
print(ct)
#plt.show()

# cantera 

gas = cnt.Solution('gri30.yaml')
gas.TPX = 298.15,101325,{'CH4':1,'O2':2,'N2':7.75}
#print(gas.X)
gas.equilibrate('UV','auto')
print(gas.T)
T_cantera=[]
for phi in x:
	gas = cnt.Solution('gri30.yaml')
	gas.TPX = 298.15,101325,{'CH4':1,'O2':2/phi,'N2':7.75/phi}
	gas.equilibrate('UV','auto')
	T_cantera.append(gas.T)
	
plt.plot(x,T_cantera,color='blue', label='cant')
plt.legend()
plt.show()

AFT with Heat loss

import cantera as cnt
import matplotlib.pyplot as plt
import math
import numpy as np

# CH4 + 2(O2 + 3.76N2) = CO2 + 2H2O + 7.52N2
def h(T,co_effs):
	
	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

R = 8.314 #j/mol-k
lower_heating_value = 50E6
Mfuel = 16.04E-3
heat_supplied = lower_heating_value*Mfuel
# low tempareures coefficients
CH4_coeff_lw = [5.14987613E+00,-1.36709788E-02,4.91800599E-05,-4.84743026E-08,1.66693956E-11,-1.02466476E+04,-4.64130376E+00]
O2_coeff_lw = [3.78245636E+00,-2.99673416E-03,9.84730201E-06,-9.68129509E-09,3.24372837E-12,-1.06394356E+03,3.65767573E+00]
N2_coeff_lw = [0.03298677E+02,0.14082404E-02,-0.03963222E-04,0.05641515E-07,-0.02444854E-10,-0.10208999E+04,0.03950372E+02]

# high temparetures coefficients
CO2_coeff_hg = [3.85746029E+00,4.41437026E-03,-2.21481404E-06,5.23490188E-10,-4.72084164E-14,-4.87591660E+04,2.27163806E+00]
N2_coeff_hg = [0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04,0.05980528E+02]
H2O_coeff_hg = [3.03399249E+00,2.17691804E-03,-1.64072518E-07,-9.70419870E-11,1.68200992E-14,-3.00042971E+04,4.96677010E+00]
O2_coeff_hg = [3.28253784E+00,1.48308754E-03,-7.57966669E-07,2.09470555E-10,-2.16717794E-14,-1.08845772E+03,5.45323129E+00]
CO_coeff_hg = [2.71518561E+00,2.06252743E-03,-9.98825771E-07,2.30053008E-10,-2.03647716E-14,-1.41518724E+04,7.81868772E+00]


def f(T,heat_loss):

	H_CO2_p = h(T,CO2_coeff_hg)
	H_H2O_p = h(T,H2O_coeff_hg)
	H_N2_p = h(T, N2_coeff_hg)
	H_O2_p = h(T, O2_coeff_hg)
	H_CO_p = h(T,CO_coeff_hg)

	T_std = 298.15 # for reactants 
	H_CH4_r = h(T_std,CH4_coeff_lw)
	H_O2_r = h(T_std,O2_coeff_lw)
	H_N2_r = h(T_std,N2_coeff_lw)


	
	H_reactants = H_CH4_r+2*(H_O2_r+3.76*H_N2_r)
	H_products = H_CO2_p+2*H_H2O_p+7.52*H_N2_p
					
	
	return(H_products-H_reactants)+(heat_loss*heat_supplied)

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

# NRM


alfa = 1
ct = 0
T_guess = 1500
x = np.linspace(0,1,21)


T_ad = []
for heat_loss in x:
	tol = 1e-6
	T_guess = 1500
	while (abs(f(T_guess,heat_loss))>tol):
		T_guess = T_guess-alfa*((f(T_guess,heat_loss))/(fprime(T_guess,heat_loss)))
		ct = ct+1
		print(T_guess)
	T_ad.append(T_guess)
plt.plot(x,T_ad,color = 'red',label='Python')
plt.xlabel('Heat loss')
plt.ylabel('Adiabatic flame temp')
plt.title('variation AFT with Heat loss')
plt.show()
  • AFT varies linearly with heat loss coefficient

Adiabatic flame temperature variation with respect to Alkane, Alkene, Alkyne :

import matplotlib.pyplot as plt
import numpy as np 
import cantera as cnt 

#For alkane
#C2H6 + 3.5(O2 + 3.76N2) <-> 2CO2 + 3H20 + 13.16N2
gas  = cnt.Solution('gri30.xml')
gas.TPX = 298.15, 101325, {'C2H6':1, 'O2':3.5, 'N2':13.16}
gas.equilibrate('HP','auto')	
T_alkane = gas.T
print(T_alkane)

#For alkene
#C2H4 + 3(O2 + 3.76N2) <-> 2CO2 + 2H20 + 11.28N2

gas  = cnt.Solution('gri30.xml')
gas.TPX = 298.15, 101325, {'C2H4':1, 'O2':3, 'N2':11.28}
gas.equilibrate('HP','auto')	
T_alkene = gas.T
print(T_alkene)


#For alkyne
#C2H6 + 2.5(O2 + 3.76N2) <-> 2CO2 + H20 + 9.4N2

gas  = cnt.Solution('gri30.xml')
gas.TPX = 298.15, 101325, {'C2H2':1, 'O2':2.5, 'N2':9.4}
gas.equilibrate('HP','auto')	
T_alkyne = gas.T
print(T_alkyne)


Temp = [T_alkane, T_alkene, T_alkyne]
plt.figure()
plt.bar(['Ethane', 'Ethene', 'Ethyne'],Temp,color='red')
plt.ylabel('Adiabetic Flame Temperatur (AFT) in K')
plt.title('AFT for various hydrocarbon fuels')
plt.show()

Adiabatic flame temperature variation with respect to a number of Carbon atom chains :

"""
Adiabatic flame temperature variation with respect to Alkane, Alkene and Alkyne
"""

# general hydrocarbon: CxHy + [x + (y/4)](O2 + 3.76N2) => xCO2 + (y/2)H2O + 3.76[x+(y/4)]N2
import matplotlib.pyplot as plt
import cantera as ct

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

#low temperature coefficients
c2h6_coeffs_l = [4.29142492E+00,-5.50154270E-03,5.99438288E-05,-7.08466285E-08,2.68685771E-11,-1.15222055E+04,2.66682316E+00]
ch4_coeffs_l = [5.14987613E+00,-1.36709788E-02,4.91800599E-05,-4.84743026E-08,1.66693956E-11,-1.02466476E+04,-4.64130376E+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]
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_l = [0.03298677E+02,0.14082404E-02,-0.03963222E-04,0.05641515E-07,-0.02444854E-10,-0.10208999E+04,0.03950372E+02]

#high temperature coefficients
n2_coeffs_h = [0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04,0.05980528E+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]

#standard temperature
Tstd = 298.15
ntype = [1,2,3]

# Part A - Newton Raphson Method
def f(T,ntype):
    # Products
    h_co2_p = h(T,co2_coeffs_h)
    h_n2_p = h(T,n2_coeffs_h)
    h_h2o_p = h(T,h2o_coeffs_h)
    
    # Reactants
    h_ch4_r = h(Tstd,ch4_coeffs_l)
    h_c2h6_r = h(Tstd,c2h6_coeffs_l)
    h_c3h8_r = h(Tstd,c3h8_coeffs_l)
    h_o2_r = h(Tstd,o2_coeffs_l)
    h_n2_r = h(Tstd,n2_coeffs_l)
    
    if ntype == 1:
        x = 1
        y = 4
        h_fuel = h_ch4_r
    elif ntype == 2:
        x = 2
        y = 6
        h_fuel = h_c2h6_r
    elif ntype == 3:
        x = 3
        y = 8
        h_fuel = h_c3h8_r
        
    H_products = x*h_co2_p + (y/2)*h_h2o_p + 3.76*(x + y/4)*h_n2_p
    H_reactants = h_fuel + (x + y/4)*h_o2_r + 3.76*(x + y/4)*h_n2_r
    return (H_products - H_reactants)

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

T_guess = 1500
tol = 1e-3
alpha = 0.2
ntype = [1,2,3]
temp = []
 
    
for i in range(0,len(ntype)):
    while(abs(f(T_guess,ntype[i])) > tol):
        T_guess = T_guess - alpha*((f(T_guess,ntype[i])/fprime(T_guess,ntype[i])))
    temp.append(T_guess)

print(temp)
# Plots
plt.figure
plt.bar(['Methane','Ethane','Propane'],temp)
#plt.ylabel('Adiabatic Flame Temperature')
#plt.title('Effect of AFT with respect to no carbon atom chain')
plt.show()
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.