Week 3 - Adiabatic Flame Temperature calculation
Adiabatic Flame Temperature: For a combustion process that takes place adiabatically with no shaft work, the temperature of the products is referred to as the adiabatic flame temperature. This is the maximum temperature that can be achieved for given reactants. Heat transfer, incomplete combustion, and dissociation all result in lower temperature. The maximum adiabatic flame temperature for a given fuel and oxidizer combination occurs with a stoichiometric mixture (correct proportions such that all fuel and all oxidizer are consumed). The amount of excess air can be tailored as part of the design to control the adiabatic flame temperature.
Adiabatic Flame Temperature
The Air fuel ratio is defined as the ratio of air and fuel of a mixture prepared for combustion. The ideal (theoretical) air-fuel ratio, for complete combustion, is called the stoichiometric air-fuel ratio. For a gasoline (petrol) engine, the stoichiometric air fuel ratio is around 14.7:1. This means that, in order to burn completely 1 kg of fuel, we need 14.7 kg of air. The combustion is possible even is the AFR is different than stoichiometric. For the combustion process to take place in a gasoline engine, the minimum AFR is around 6:1 and the maximum can go up to 20:1.
Air fuel ratio (AF or AFR) is defined as the ratio between the mass of air ma and mass fuel mf,
`AFR = m_a/m_f`
Equivalence Ratio is defined as the ratio of the actual fuel/air ratio to the stoichiometric fuel/air ratio. Stoichiometric combustion occurs when all the oxygen is consumed in the reaction, and there is no molecular oxygen(O2) in the products. If the equivalence ratio is equal to one, the combustion is stoichiometric. If it is < 1, the combustion is lean with excess air, and if it is >1, the combustion is rich with incomplete combustion.
The adiabatic flame temperature is dependent on two types of combustion processes:
A.) Constant Volume Process: Systems such as enclosed explosions and otto cycle deal with constant volume combustion. The Adiabatic flame temperature can be calculated by using
`H_r-H_p-R_u(N_rT_i-N_pT_ad)=0`
B.) Constant Pressure Process: Systems such as gas turbine combustion and rocket engines deal with constant pressure combustion. The adiabatic flame temperature can be calculated as:
`DeltaH=0`
1. Adiabatic flame temperature variation with Equivalence ratio (In constant volume-chamber) :
Lean mixture equation:
`CH_4 + 2/phi * (O_2 +3.76*N_2) -> CO_2 + 2*H_2O + 7.5/phi*N_2 + (2/phi - 2)*O_2`
Rich Mixture Equation:
`CH_4 + 2/phi*(O_2+3.76*N_2) -> (4/phi-3)*CO +2*H_2O +7.5/phi * N_2 + (4-4/phi)O_2`
import matplotlib.pyplot as plt
import math
import numpy as np
import cantera as ct
def h(T,coeff):
# To evaluate enthalpy
R = 8.314 #J/mol-K
a1 = coeff[0]
a2 = coeff[1]
a3 = coeff[2]
a4 = coeff[3]
a5 = coeff[4]
a6 = coeff[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]
O2_coeff_h = [3.28253784E+00, 1.48308754E-03, -7.57966669E-07, 2.09470555E-10, -2.16717794E-14, -1.08845772E+03, 5.45323129E+00]
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]
CO_coeff_h = [2.71518561E+00, 2.06252743E-03, -9.98825771E-07, 2.30053008E-10, -2.03647716E-14, -1.41518724E+04, 7.81868772E+00]
H2O_coef_h = [3.03399249E+00, 2.17691804E-03, -1.64072518E-07,-9.70419870E-11, 1.68200992E-14, -3.00042971E+04, 4.96677010E+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]
"""
Lean mixture equatiom
CH4 + 2/phi * (O2 +3.76N2) --> CO2 + 2*H2O + 7.5/phi*N2 + (2/phi - 2)*O2
Rich Mixture Equation
CH4 + 2/phi*(O2+3.76*N2) --> (4/phi-3)*CO +2*H2O +7.5/phi * N2 + (4-4/phi)O2
"""
def f(T,phi):
R = 8.314 #J/mol-K
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/phi)*h_O2_r + (7.52/phi)*h_N2_r
#print(H_reactants)
N_reactants = 1 + 9.52/phi
h_CO2_p = h(T,CO2_coeff_h)
h_H2O_p = h(T, H2O_coef_h)
h_N2_p = h(T, N2_coeff_h)
h_O2_p = h(T,O2_coeff_h)
h_CO_p = h(T,CO_coeff_h)
#Lean Mixture products
if phi<1:
H_products = h_CO2_p + 2*h_H2O_p + (7.5/phi)*h_N2_p + (2/phi - 2)*h_O2_p
N_products = (1+ (2/phi -2) + 2 + (7.52/phi))
# Rich mixture
elif phi>1:
H_products =(4/phi - 3)*h_CO2_p + 2*h_H2O_p + (7.5/phi)*h_N2_p + (4-4/phi)*h_CO_p
N_products = (4/phi -3) + 2 +(4- 4/phi) +(7.52/phi)
return (H_products-H_reactants) - (N_products*T - N_reactants*Tstd)*R
def fprime(T,phi):
return (f(T+1e-6,phi) - f(T,phi))/1e-6 #Backward differencing
T_guess = 1500
phi = np.linspace(0.1,2,100)
tol = 1e-3
alpha = 0.2
iter = 0
Temp=[]
for i in range(0,len(phi)):
while np.abs(f(T_guess,phi[i])) > tol :
T_guess = T_guess - alpha*((f(T_guess,phi[i]))/(fprime(T_guess,phi[i])))
iter += 1
Temp.append(T_guess)
print(T_guess)
#Cantera
gas = ct.Solution('gri30.xml')
T_cant =[]
for i in range(0,len(phi)):
m_net = 1 + (2/phi[i]) + (7.52/phi[i])
m_ch4 = 1/ m_net
m_o2 = (2/phi[i])/ m_net
m_n2 = (7.52/phi[i])/ m_net
gas.TPX = 298.15, 101325, {'CH4':m_ch4, 'O2':m_o2, 'N2':m_n2}
gas.equilibrate('UV','auto')
T_cant.append(gas.T)
print(T_cant)
print(T_guess)
plt.plot(phi,T_cant,'-o',color='blue',label='Cantera')
plt.plot(phi,Temp,'-o',color='red',label='Pyhton')
plt.xlabel('Equivalence ratio')
plt.ylabel('Adiabatic flame temperature [K]')
plt.title('AFT v/s Equivalence ratio for Methane Combustion')
plt.legend()
plt.show()
Output:
The reason that there is difference in AFT calculation by Python logic and Cantera is that python only takes into account a few of by products and cantera on the whole accounts for more by products.
2. Flame temperature variation with Heat loss (In constant pressure-chamber) : Next Flame temperature is calculated by constant pressure method
import matplotlib.pyplot as plt
import math
import numpy as np
def h(T,coeff):
# To evaluate enthalpy
R = 8.314 #J/mol-K
a1 = coeff[0]
a2 = coeff[1]
a3 = coeff[2]
a4 = coeff[3]
a5 = coeff[4]
a6 = coeff[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]
O2_coeff_h = [3.28253784E+00, 1.48308754E-03, -7.57966669E-07, 2.09470555E-10, -2.16717794E-14, -1.08845772E+03, 5.45323129E+00]
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]
CO_coeff_h = [2.71518561E+00, 2.06252743E-03, -9.98825771E-07, 2.30053008E-10, -2.03647716E-14, -1.41518724E+04, 7.81868772E+00]
H2O_coef_h = [3.03399249E+00, 2.17691804E-03, -1.64072518E-07,-9.70419870E-11, 1.68200992E-14, -3.00042971E+04, 4.96677010E+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]
"""
stoichiometric Eqn
#CH4 + 2(O2 + 3.76N2) = CO2 + 2H2O + 7.5N2
"""
def f(T,heat_loss):
Tstd = 298.15
"""
Qmax= h(Tstd,ch4_coeff_l)+h(Tstd, O2_coeff_l)-h(Tstd,CO2_coeff_h)-h(Tstd, H2O_coef_h)
"""
mass = 1 #Kg
LHV = 802300
Qmax = mass*LHV
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*h_O2_r + 7.52*h_N2_r
h_CO2_p = h(T,CO2_coeff_h)
h_H2O_p = h(T, H2O_coef_h)
h_N2_p = h(T, N2_coeff_h)
H_products = h_CO2_p + 2*h_H2O_p + 7.5*h_N2_p
return ((H_products-H_reactants) + Qmax*heat_loss)
def fprime(T,heat_loss):
return (f(T+1e-6,heat_loss) - f(T,heat_loss))/1e-6 #Backward differencing
T_guess = 1500
heat_loss = np.linspace(0,1,10)
tol = 1e-3
alpha = 0.2
iter = 0
Temp=[]
for i in range(0,len(heat_loss)):
while np.abs(f(T_guess,heat_loss[i])) > tol :
T_guess = T_guess - alpha*((f(T_guess,heat_loss[i]))/(fprime(T_guess,heat_loss[i])))
iter += 1
Temp.append(T_guess)
print(Temp)
plt.plot(heat_loss,Temp,'-o',color='blue')
plt.xlabel('Heat loss Percentage')
plt.grid('on')
plt.ylabel('Adiabatic flame temperature [K]')
plt.title('Heat loss and AFT plot')
plt.show()
It can be deduced from the above plot that as the heat loss is increased from the container the Temperature inside will be decreasing
3. Adiabatic flame temperature variation with respect to Alkane, Alkene, Alkyne : The Adiabatic flame temperature distribution with respect to different bonds in carbon atom is thus calculated here.
import matplotlib.pyplot as plt
import math
import numpy as np
def h(T,coeff):
# To evaluate enthalpy
R = 8.314 #J/mol-K
a1 = coeff[0]
a2 = coeff[1]
a3 = coeff[2]
a4 = coeff[3]
a5 = coeff[4]
a6 = coeff[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
C2H6_coeff_l = [4.29142492E+00, -5.50154270E-03, 5.99438288E-05,-7.08466285E-08, 2.68685771E-11, -1.15222055E+04, 2.66682316E+00]
C2H4_coeff_l = [3.95920148E+00, -7.57052247E-03, 5.70990292E-05, -6.91588753E-08, 2.69884373E-11, 5.08977593E+03, 4.09733096E+00]
C2H2_coeff_l = [8.08681094E-01, 2.33615629E-02, -3.55171815E-05, 2.80152437E-08, -8.50072974E-12, 2.64289807E+04, 1.39397051E+01]
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]
O2_coeff_h = [3.28253784E+00, 1.48308754E-03, -7.57966669E-07, 2.09470555E-10, -2.16717794E-14, -1.08845772E+03, 5.45323129E+00]
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]
CO_coeff_h = [2.71518561E+00, 2.06252743E-03, -9.98825771E-07, 2.30053008E-10, -2.03647716E-14, -1.41518724E+04, 7.81868772E+00]
H2O_coef_h = [3.03399249E+00, 2.17691804E-03, -1.64072518E-07,-9.70419870E-11, 1.68200992E-14, -3.00042971E+04, 4.96677010E+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]
"""
stoichiometric Eqn
#ethane
#C2H6 + 3.5(O2 + 3.76N2) = 2CO2 + 3H2O + 13.16N2
#ethene
#C2H4 + 3(O2 + 3.76N2) = 2CO2 + 2H2O + 11.28N2
#ethyne
#C2H2 + 2.5(O2 + 3.76N2) = 2CO2 + H2O + 9.4N2
"""
#Ethane
def f(T):
Tstd = 298.15
h_C2H6_r = h(Tstd,C2H6_coeff_l)
h_O2_r = h(Tstd, O2_coeff_l)
h_N2_r = h(Tstd, N2_coeff_l)
H_reactants = h_C2H6_r + 3.5*h_O2_r + 13.16*h_N2_r
h_CO2_p = h(T,CO2_coeff_h)
h_H2O_p = h(T, H2O_coef_h)
h_N2_p = h(T, N2_coeff_h)
H_products = 2*h_CO2_p + 3*h_H2O_p + 13.16*h_N2_p
return (H_products-H_reactants)
def fprime(T):
return (f(T+1e-6) - f(T))/1e-6 #Backward differencing
T_guess = 1500
tol = 1e-3
alpha = 0.2
ct = 0
while(abs(f(T_guess)) > tol):
T_guess = T_guess - alpha*((f(T_guess))/(fprime(T_guess)))
ct += 1
print(T_guess)
#Ethene
def g(T):
Tstd = 298.15
h_C2H4_r = h(Tstd,C2H4_coeff_l)
h_O2_r = h(Tstd, O2_coeff_l)
h_N2_r = h(Tstd, N2_coeff_l)
H_reactants = h_C2H4_r + 3*h_O2_r + 11.28*h_N2_r
h_CO2_p = h(T,CO2_coeff_h)
h_H2O_p = h(T, H2O_coef_h)
h_N2_p = h(T, N2_coeff_h)
H_products = 2*h_CO2_p + 2*h_H2O_p + 11.28*h_N2_p
return (H_products-H_reactants)
def gprime(T):
return (g(T+1e-6) - g(T))/1e-6 #Backward differencing
T_guess_1 = 1500
tol = 1e-3
alpha = 0.2
a = 0
while(abs(g(T_guess_1)) > tol):
T_guess_1 = T_guess_1 - alpha*((g(T_guess_1))/(gprime(T_guess_1)))
a += 1
print(T_guess_1)
#Ethyne
#C2H2 + 2.5(O2 + 3.76N2) = 2CO2 + H2O + 9.4N2
def k(T):
Tstd = 298.15
h_C2H2_r = h(Tstd,C2H2_coeff_l)
h_O2_r = h(Tstd, O2_coeff_l)
h_N2_r = h(Tstd, N2_coeff_l)
H_reactants = h_C2H2_r + 2.5*h_O2_r + 9.4*h_N2_r
h_CO2_p = h(T,CO2_coeff_h)
h_H2O_p = h(T, H2O_coef_h)
h_N2_p = h(T, N2_coeff_h)
H_products = 2*h_CO2_p + h_H2O_p + 9.4*h_N2_p
return (H_products-H_reactants)
def kprime(T):
return (k(T+1e-6) - k(T))/1e-6 #Backward differencing
T_guess_2 = 1500
tol = 1e-3
alpha = 0.2
a = 0
while(abs(k(T_guess_2)) > tol):
T_guess_2 = T_guess_2 - alpha*((k(T_guess_2))/(kprime(T_guess_2)))
a += 1
print(T_guess_2)
x = ["Ethane", "Ethene", "Ethyne"]
y = [T_guess,T_guess_1,T_guess_2]
plt.bar(x,y)
plt.title('AFT Variation')
plt.show()
It is clear that as we increase the number of bonds between carbon atoms the Flame temperature will also rise.
4. Adiabatic flame temperature variation with respect to number of Carbon atom chains : A variation on AFT with respect to number of carbon atoms is calculated as:
import matplotlib.pyplot as plt
import math
import numpy as np
def h(T,coeff):
# To evaluate enthalpy
R = 8.314 #J/mol-K
a1 = coeff[0]
a2 = coeff[1]
a3 = coeff[2]
a4 = coeff[3]
a5 = coeff[4]
a6 = coeff[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
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]
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]
O2_coeff_h = [3.28253784E+00, 1.48308754E-03, -7.57966669E-07, 2.09470555E-10, -2.16717794E-14, -1.08845772E+03, 5.45323129E+00]
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]
CO_coeff_h = [2.71518561E+00, 2.06252743E-03, -9.98825771E-07, 2.30053008E-10, -2.03647716E-14, -1.41518724E+04, 7.81868772E+00]
H2O_coef_h = [3.03399249E+00, 2.17691804E-03, -1.64072518E-07,-9.70419870E-11, 1.68200992E-14, -3.00042971E+04, 4.96677010E+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]
"""
stoichiometric Eqn
#Methane
#CH4 + 2(O2 + 3.76N2) = CO2 + 2H2O + 7.5N2
#Ethane
#C2H6 + 3.5(O2 + 3.76N2) = 2CO2 + 3H2O + 13.16N2
#Propane
#C3H8 + 5(O2 + 3.76N2) = 3CO2 + 4H2O + 18.8N2
"""
#Methane
def f(T):
h_CO2_p = h(T,CO2_coeff_h)
h_H2O_p = h(T, H2O_coef_h)
h_N2_p = h(T, N2_coeff_h)
H_products = h_CO2_p + 2*h_H2O_p + 7.5*h_N2_p
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*h_O2_r + 7.52*h_N2_r
return H_products-H_reactants
def fprime(T):
return (f(T+1e-6) - f(T))/1e-6 #Backward differencing
T_guess = 1500
tol = 1e-3
alpha = 0.2
ct = 0
while(abs(f(T_guess)) > tol):
T_guess = T_guess - alpha*((f(T_guess))/(fprime(T_guess)))
ct += 1
print(T_guess)
#Ethane
def q(T):
Tstd = 298.15
h_C2H6_r = h(Tstd,C2H6_coeff_l)
h_O2_r = h(Tstd, O2_coeff_l)
h_N2_r = h(Tstd, N2_coeff_l)
H_reactants = h_C2H6_r + 3.5*h_O2_r + 13.16*h_N2_r
h_CO2_p = h(T,CO2_coeff_h)
h_H2O_p = h(T, H2O_coef_h)
h_N2_p = h(T, N2_coeff_h)
H_products = 2*h_CO2_p + 3*h_H2O_p + 13.16*h_N2_p
return (H_products-H_reactants)
def qprime(T):
return (q(T+1e-6) - q(T))/1e-6 #Backward differencing
T_guess_1 = 1500
tol = 1e-3
alpha = 0.2
a = 0
while(abs(q(T_guess_1)) > tol):
T_guess_1 = T_guess_1 - alpha*((q(T_guess_1))/(qprime(T_guess_1)))
a += 1
print(T_guess_1)
#Propane
#C3H8 + 5(O2 + 3.76N2) = 3CO2 + 4H2O + 18.8N2
def k(T):
Tstd = 298.15
h_C3H8_r = h(Tstd,C3H8_coeff_l)
h_O2_r = h(Tstd, O2_coeff_l)
h_N2_r = h(Tstd, N2_coeff_l)
H_reactants = h_C3H8_r + 5*h_O2_r + 18.8*h_N2_r
h_CO2_p = h(T,CO2_coeff_h)
h_H2O_p = h(T, H2O_coef_h)
h_N2_p = h(T, N2_coeff_h)
H_products = 3*h_CO2_p + 4*h_H2O_p + 18.8*h_N2_p
return (H_products-H_reactants)
def kprime(T):
return (k(T+1e-6) - k(T))/1e-6 #Backward differencing
T_guess_2 = 1500
tol = 1e-3
alpha = 0.2
a = 0
while(abs(k(T_guess_2)) > tol):
T_guess_2 = T_guess_2 - alpha*((k(T_guess_2))/(kprime(T_guess_2)))
a += 1
print(T_guess_2)
x = ["Methane", "Ethane", "Propane"]
y = [T_guess,T_guess_1,T_guess_2]
plt.bar(x,y)
plt.title('AFT Variation')
plt.show()
