Week 3 - Adiabatic Flame Temperature calculation
1. Equivalence Ratio
One of the most important features of a combustion process is the highest temperature
of the combustion products that can be achieved. The temperature of the products
will be greatest when there are no heat losses to the surrounding environment and all of
the energy released from combustion is used to heat the products. In this project, a Python program was written to se the effect of equivalence ratio on AFT of Methane. An adiabatic constant-pressure analysis was used here to calculate the AFT. The task is finding the product temperature given the enthalpy of reactants.With an initial guess of flame temperature, the product temperature can be estimated by linear interpolation (Newton-Raphson Method). This method, although more accurate, still assumes complete combustion to the major products.
As such, the product (flame) temperature increases almost linearly with equivalence ratio,ϕ">ϕϕ, for lean combustion. As expected, the flame temperature peaks at the stoichiometric ratio. In rich combustion, the flame temperature decreases withϕ">ϕϕ.
The reason because AFT using Cantera is lower than using Python, it's that Cantera is based on The equilibrium state which determines the species concentrations and temperature under certain constraints such as constant enthalpy, pressure, or temperature, and it can be determined by minimizing the Gibbs free energy. Also, the dissociation of products at high temperature (T>1,500 K at ambient pressure) can take a significant portion of energy from combustion and hence the product temperature is lower than that calculated with only major components as products.
Code--->
"""
In this program we are going to calculate the AFT of methane at different equivalence ratios
Equivalence ratio phi=1/lambda
Stoichiometric CH4 + 2(O2+3.76N2)= CO2 +2H2O + 7.52N2
Lean mixture CH4 + (2/phi)(O2+3.76N2)= CO2 + 2H2O + (7.52/phi)N2 + (2/phi-2)O2
Rich misture CH4 + (2/phi)(O2+3.76N2)= (4/phi-3)CO2 + 2H2O + (7.52/phi)N2 + (4-4/phi)*CO
"""
import matplotlib.pyplot as plt
import math
"""
Function to evaluate enthalpy, constan pressure system
"""
def h(T,Coeff):
R= 8.314 #Jol/ 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
# Defining the coefficients
ch4_coeffs_l = [5.14987613E+00,-1.36709788E-02,4.91800599E-05,-4.84743026E-08,1.66693956E-11,-1.02466476E+04,-4.64130376E+00]
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]
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]
n2_coeffs_h = [0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04,0.05980528E+02]
o2_coeffs_h = [3.28253784E+00,1.48308754E-03,-7.57966669E-07,2.09470555E-10,-2.16717794E-14,-1.08845772E+03,5.45323129E+00]
co_coeffs_h = [2.71518561E+00,2.06252743E-03,-9.98825771E-07,2.30053008E-10,-2.03647716E-14,-1.41518724E+04,7.81868772E+00]
def f(T,phi):
Tstd =298.15 #k
R= 8.314 #Jol/ mol-k
h_ch4_r = h(Tstd,ch4_coeffs_l)
h_o2_r = h(Tstd,o2_coeffs_l)
h_n2_r = h(Tstd,n2_coeffs_l)
h_co2_p = h(T,co2_coeffs_h)
h_h2o_p = h(T,h2o_coeffs_h)
h_n2_p = h(T,n2_coeffs_h)
h_o2_p = h(T,o2_coeffs_h)
h_co_p = h(T,co_coeffs_h)
if phi<=1:
H_reactants = h_ch4_r + (2/phi)*h_o2_r + (7.52/phi)*h_n2_r
H_products = h_co2_p + 2*h_h2o_p + (7.52/phi)*h_n2_p + (2/phi-2)*h_o2_p
N_reactants = 1 + (2/phi) + (7.52/phi)
N_products = 1 + 2 + (7.52/phi) + (2/phi-2)
if phi>1:
H_reactants = h_ch4_r + (2/phi)*h_o2_r + (7.52/phi)*h_n2_r
H_products = (4/phi-3)*h_co2_p + 2*h_h2o_p + (7.52/phi)*h_n2_p + (4-4/phi)*h_co_p
N_reactants = 1 + (2/phi) + (7.52/phi)
N_products = (4/phi-3) + 2 + (7.52/phi) + (4-4/phi)
return H_products - H_reactants - R*(N_products*T - N_reactants*Tstd)
print(f)
def fprime(T, phi):
return ((f(T+ 1e-6, phi) - f(T,phi))/1e-6)
T_guess= [1500 for _ in range(14)] #K
tol = 1e-3
ct = 0
alpha = 0.5
i = 0
phi=[0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1,1.1,1.2,1.3,1.4]
for i in range(0,14):
while (abs(f(T_guess[i],phi[i]))>tol):
T_guess[i] = T_guess[i] - alpha*((f(T_guess[i],phi[i])/fprime(T_guess[i], phi[i])))
ct=ct+1
print(phi[i],T_guess[i])
T_aft=[681.6845,1009.5528,1301.8703,1567.8621,1812.8759,2040.4161,2252.9810,2452.4598,2640.3409,2817.8310,2704.1671,2592.9540,2484.0052,2377.1327]
plt.plot(phi,T_aft,'.',color='black')
"""
Using Cantera
"""
import cantera as ct
gas=ct.Solution('gri30.xml')
T_aft_cantera=[]
for i in range (0,len(phi)):
Total_mole = 1 + 2/(phi[i]) +7.52/(phi[i])
mole_fraction_ch4 = 1/Total_mole
mole_fraction_o2 = (2/(phi[i]))/Total_mole
mole_fraction_n2 = (7.52/(phi[i]))/Total_mole
gas.TPX= 298.15,101325,{'CH4':mole_fraction_ch4,'O2':mole_fraction_o2,'N2':mole_fraction_n2}
gas.equilibrate('UV','auto')
T_aft_cantera.append(gas.T)
plt.plot(phi,T_aft_cantera,'.',color='grey')
plt.title('Effect of Equivalence ratio on the AFT - Methane')
plt.xlabel('Equivalence Ratio')
plt.ylabel('Adiabatic Flame Temperature[K]')
plt.legend(['Python','Cantera'])
plt.grid()
plt.show()2. Heat loss at constant pressure
Python program was written to se the effect of Heat loss on FT of Methane. Heat loss was take a fraction of the "maximum possible heat loss". As higher the heat loss, as lower the FT. From previous graph AFT-Methane=2817.8310 K and in this one FT-Methane=1679.7053 K with het loss=35%
Code--->
"""
In this program we are going to calculate the heat loss on FT from a reactor with constant pressure
Methane
Stoichiometric CH4 + 2(O2+3.76N2)= CO2 +2H2O + 7.52N2
LHV=800900 J mol-1
"""
import matplotlib.pyplot as plt
import numpy as np
"""
Function to evaluate enthalpy, constan pressure system
"""
def h(T,Coeff):
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
# Defining the coefficients
ch4_coeffs_l = [5.14987613E+00,-1.36709788E-02,4.91800599E-05,-4.84743026E-08,1.66693956E-11,-1.02466476E+04,-4.64130376E+00]
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]
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]
n2_coeffs_h = [0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04,0.05980528E+02]
def f(T,Heat_loss):
"""
Function to evaluate enthalpy, constan pressure system
"""
Tstp =298.15 #k
LHV = 800950
h_ch4_r = h(Tstp,ch4_coeffs_l)
h_o2_r= h(Tstp,o2_coeffs_l)
h_n2_r = h(Tstp,n2_coeffs_l)
h_co2_p = h(T,co2_coeffs_h)
h_h2o_p = h(T,h2o_coeffs_h)
h_n2_p = h(T,n2_coeffs_h)
h_reactants = h_ch4_r+ 2*h_o2_r + 7.52*h_n2_r
h_products= h_co2_p+ 2*h_h2o_p + 7.52*h_n2_p
Loss = Heat_loss* LHV
return h_products - h_reactants + Loss
def fprime(T, Heat_loss):
return (f(T+ 1e-6, Heat_loss)- f(T,Heat_loss))/1e-6
T_guess=1500
tol = 1e-3
ct = 0
alpha = 1
Heat_loss=np.linspace(0,1,50)
T_ft=[]
for i in range (0,len(Heat_loss)):
while (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])))
ct=ct+1
T_ft.append(T_guess)
plt.grid()
plt.plot(Heat_loss,T_ft,'.',color='black')
plt.xlabel(' % Heat Loss')
plt.ylabel('Flame Temperature[K]')
plt.title('Effect of heat loss on FT - Methane')
plt.show()3. Alkanes, alkynes, and alkenes variants of 2 Carbon atom with heat loss=0.35
The effect of the number of hydrogen atoms in double hydrocarbons it can be seen in this program. AFT decrease while increase the number of hydrogen atoms. At lesser number of hydrogen atoms, at lesser water molecules of formed; which means less heat loss.
Code--->
"""
Estimation of FT for a Heat Loss Factor of 0.35 only for -ane, -ene and -yne variants of 2 Carbon atom
C2Hy + (2+ y/4)(O2 + 3.76N2)---->2CO2 + (y/2)H2O + 3.76*(2 + y/4)N2
"""
import matplotlib.pyplot as plt
import math
import numpy as np
"""
Function to evaluate enthalpy, constan pressure system
"""
def h(T,Coeff):
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
# Defining the coefficients
c2h2_coeffs_l = [8.08681094E-01,2.33615629E-02,-3.55171815E-05,2.80152437E-08,-8.50072974E-12,2.64289807E+04,1.39397051E+01]
c2h4_coeffs_l = [3.95920148E+00,-7.57052247E-03,5.70990292E-05,-6.91588753E-08,2.69884373E-11,5.08977593E+03,4.09733096E+00]
c2h6_coeffs_l = [4.29142492E+00,-5.50154270E-03,5.99438288E-05,-7.08466285E-08,2.68685771E-11,-1.15222055E+04,2.66682316E+00]
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]
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]
n2_coeffs_h = [0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04,0.05980528E+02]
def f(T,n):
Tstp = 298.15 #K
R= 8.314 #J/ mol-k
y=n
h_c2h2_r = h(Tstp,c2h2_coeffs_l)
h_c2h4_r = h(Tstp,c2h4_coeffs_l)
h_c2h6_r = h(Tstp,c2h6_coeffs_l)
h_o2_r = h(Tstp,o2_coeffs_l)
h_n2_r = h(Tstp,n2_coeffs_l)
h_loss = 0.35
if n==2: #Ethyne LHV=1258020 J/mol
h_reactants = h_c2h2_r + (2 + y/4)*h_o2_r + 3.76*(2 +y/4)*h_n2_r - (1258020*h_loss)
if n==4: #Ethene LHV=1323474 J/mol
h_reactants = h_c2h4_r + (2 + y/4)*h_o2_r + 3.76*(2 +y/4)*h_n2_r - (1323474*h_loss)
if n==6: #Ethane LHV=1428275 J/mol
h_reactants = h_c2h6_r + (2 + y/4)*h_o2_r + 3.76*(2 +y/4)*h_n2_r - (1428275*h_loss)
h_co2_p = h(T,co2_coeffs_h)
h_h2o_p = h(T,h2o_coeffs_h)
h_n2_p = h(T,n2_coeffs_h)
h_products= 2*h_co2_p + (y/2)*h_h2o_p + 3.76*(2+y/4)*h_n2_p
return h_products - h_reactants
def fprime(T, n):
return (f(T+ 1e-6, n)- f(T,n))/1e-6
T_guess= [1500 for _ in range(3)] #K
n=[2,4,6]
tol = 1e-3
ct = 0
alpha = 0.2
i = 0
for i in range(0,3):
while (abs(f(T_guess[i],n[i]))>tol):
T_guess[i] = T_guess[i] - alpha*(f(T_guess[i],n[i])/fprime(T_guess[i],n[i]))
ct=ct+1
print(n[i],T_guess[i])
plt.plot(n[i],T_guess[i],'o',color='grey')
plt.grid()
plt.title('Estimation of FT for a Heat Loss Factor of 0.35')
plt.xlabel('Number of hydrogen atoms in double hydrocarbons C2Hx')
plt.ylabel('Adiabatic Flame Temperature[K]')
plt.legend(['C2H2-Acetylene','C2H4-Ethylene','C2H6-Ethane'])
plt.show()
4. Effect of the number of carbon atoms (from 1 to 3) on the final temperature, with H_loss = 0
With this progrma we can see, that increasing the number of Carbon atoms result in the increase of AFT.
Code-->
"""
Effect of the number of carbon atoms (from 1 to 3) on the final temperature for a general alkane.
H_loss = 0
CxHy + (x+ y/4)(O2 + 3.76N2)---->xCO2 + (y/2)H2O + 3.76*(x + y/4)N2
"""
import matplotlib.pyplot as plt
import math
import numpy as np
"""
Function to evaluate enthalpy, constan pressure system
"""
def h(T,Coeff):
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
# Defining the coefficients
ch4_coeffs_l = [5.14987613E+00,-1.36709788E-02,4.91800599E-05,-4.84743026E-08,1.66693956E-11,-1.02466476E+04,-4.64130376E+00]
c2h6_coeffs_l = [4.29142492E+00,-5.50154270E-03,5.99438288E-05,-7.08466285E-08,2.68685771E-11,-1.15222055E+04,2.66682316E+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]
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]
n2_coeffs_h = [0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04,0.05980528E+02]
def f(T,n):
Tstp = 298.15 #K
R= 8.314 #J/ mol-k
x=n
y=2*n+2
h_ch4_r = h(Tstp,ch4_coeffs_l)
h_c2h6_r = h(Tstp,c2h6_coeffs_l)
h_c3h8_r = h(Tstp,c3h8_coeffs_l)
h_o2_r = h(Tstp,o2_coeffs_l)
h_n2_r = h(Tstp,n2_coeffs_l)
if n==1: #Ethyne LHV=1258020 J/mol
h_reactants = h_ch4_r + (x + y/4)*h_o2_r + 3.76*(x +y/4)*h_n2_r
if n==2: #Ethene LHV=1323474 J/mol
h_reactants = h_c2h6_r + (x + y/4)*h_o2_r + 3.76*(x +y/4)*h_n2_r
if n==3: #Ethane LHV=1428275 J/mol
h_reactants = h_c3h8_r + (x + y/4)*h_o2_r + 3.76*(x +y/4)*h_n2_r
h_co2_p = h(T,co2_coeffs_h)
h_h2o_p = h(T,h2o_coeffs_h)
h_n2_p = h(T,n2_coeffs_h)
h_products= x*h_co2_p + (y/2)*h_h2o_p + 3.76*(x+y/4)*h_n2_p
return h_products - h_reactants
def fprime(T, n):
return (f(T+ 1e-6, n)- f(T,n))/1e-6
T_guess= [1500 for _ in range(3)] #K
n=[1,2,3]
tol = 1e-3
ct = 0
alpha = 0.2
i = 0
for i in range(0,3):
while (abs(f(T_guess[i],n[i]))>tol):
T_guess[i] = T_guess[i] - alpha*(f(T_guess[i],n[i])/fprime(T_guess[i],n[i]))
ct=ct+1
print(n[i],T_guess[i])
plt.plot(n[i],T_guess[i],'o',color='grey')
plt.grid()
plt.title('Effect of alkane carbon atoms on the AFT')
plt.xlabel('Number of carbon atoms in alkane [CxHy]')
plt.ylabel('Flame Temperature[K]')
plt.legend(['n=1 CH4-Methane','n=2 C2H6-Ethene','n=3 C3H8-Propane'])
plt.show()