Adiabatic Flame temperature for different Hydrocarbons using Python and Cantera
Aim –
- Determination of Adiabatic Flame temperature of Methane in a Constant Volume process at different equivalence ratio using Python and comparing it with the results obtained from Cantera.
- Determination of Flame temperature of different Hydrocarbons (Ethane, Ethene and Acetylene) with heat loss in a constant pressure process.
- Determination of Adiabatic Flame temperatures of different alkanes (Methane, Ethane and Propane)
- Studying the variation in Flame temperature of Methane with increase in heat loss in a constant pressure combustion
Theory –
The maximum temperature obtained when a fuel is burnt adiabatically in air is termed as Adiabatic Flame temperature.
The combustion can be achieved either in a constant pressure process or in a constant volume process.
Assuming no dissociation, the combustion of a general hydrocarbon, in terms of equivalence ratio (`phi`) can be written as –
`ϕ≤1`
`C_x H_y+a(O_2+3.76N_2 )→b CO_2+cH_2 O+dO_2+eN_2`
By performing a simple mass balance of each element, we get
|
a |
(x+y/4)/`phi` |
|
b |
x |
|
c |
y/2 |
|
d |
(1-`phi`)* a |
|
e |
3.76a |
`phi>1`
`C_x H_y+a(O_2+3.76N_2 )→b CO_2+cCO+dH_2 O+eN_2`
Again by performing the mass balance,
|
a |
(x+y/4)/`phi` |
|
b |
2a-y/2-x |
|
c |
2x+y/2-2a |
|
d |
y/2 |
|
e |
3.76a |
When the mixture is rich, H2 formation starts when a<0. This scenario has been ignored in this study to make it simpler.
Constant Pressure Combustion
For combustion in a constant pressure process, the adiabatic flame temperature is the temperature when
`ΔH=0`
`H_p=H_r`
When heat loss is involved, the temperature is called as Flame temperature and is achieved when
`H_p+Heat loss=H_r`
Constant Volume Combustion
In a constant volume process, by first law the work done is zero. Thus, for an adiabatic process,
`ΔU=0`
`U=H-PV`
Thus
`ΔU=ΔH-ΔPV=0`
For an ideal gas
PV = nRT
Thus, Adiabatic Flame temperature is obtained using an iterative process when
`H_p-H_r-R(N_p T_(ad)-N_r T_r)=0`
Here R = 8.314 J/molK
`N_p`is the total number of moles of all the products
`N_r`is the total number of moles of all the reactants
Solution Algorithm
- Input the variables, namely, x and y of the hydrocarbons, initial temperature and pressure, equivalence Ratio
- Calculate the enthalpy of reactants and products using the NASA polynomials and the coefficients from a GRI Mechanism File available online, depending on the type of process (constant volume or pressure).
- Determine the first derivative of the appropriate equation using forward difference formula
- To solve the nonlinear equation, utilize the Newton Raphson method.
Results
From the below graph, it can be concluded that the AFT values obtained from first law analysis are significantly higher than those obtained from Cantera. There can be two possible reasons for this. First, the number of species considered by Cantera is significantly higher and all these extra species absorb a certain part of the energy. Second, Cantera also incorporates the second law and minimizes the gibbs free energy. This is an important consideration in chemical systems which govern when the combustion process stops.
It can be observed from the below graph that Flame temperature increases with the increase in carbon atoms in an alkane. This is because with the increase in carbon atoms the number of water moles on the product side also increases, thus increasing the enthalpy of the products.
With the increase in the number of bonds (alkane, alkene and alkyne), the amount of energy released after combustion due to breaking of bonds is higher. Thus the Flame temperature also increases when transition from single bonded carbon atoms to triple bonded carbon atoms take place.
CODE FOR AFT OF CH4 VS EQUIVALENCE RATIO (PART 1 OF PROBLEM STATEMENT)
https://drive.google.com/file/d/1EzudYpE2DWQeoZqnEREkAWUjMptpTZMl/view?usp=sharing
"""
THIS IS A GENERAL CODE WHICH CAN BE USED FOR AFT CALCULATION OF ANY HYDROCARBON PROVIDED THE NASA POLYNOMIAL COEFFICIENT
OF THAT HC IS AVAILABLE
CH4 combustion in a constant volume chamber INITIALLY at STP conditions
Effect of equivalence ratio(phi) on AFT has been studied
phi = (F/A)/(F/A)_stoich
From second edition SR Turns Pg - 48
Assuming no dissociation, combustion of an arbitrary hydrocarbon is represented as
when phi<=1 i.e. Lean combustion
CxHy + a(O2+3.76N2) =>bCO2 + cH2O + dO2 + eN2
a = (x+y/4)/phi
b = x
c = y/2
d = ((1-phi)/phi)*(x+(y/4))
e = 3.76*a
when phi>1 i.e. Rich combustion;H2 has not been included to reduce the complexity. This is only valid till b=0. For CH4
the limit of phi is 4/3. Beyond this b becomes negative which is erroneous.
CxHy + a(O2+3.76N2) =>bCO2 + cCO + dH2O + eN2
a = (x+y/4)/phi
d = y/2
b = 2*((x+y/4)/phi)-(y/2)-x
c = 2*x + (y/2) - 2*((x+y/4)/phi)
e = 3.76*a
|Sharique Iqbal|24 Oct 2020|
"""
import math
import cantera as ct
import numpy as np
import matplotlib.pyplot as plt
#CxHy -> Arbitrary Hydrocarbon
x = 1
y = 4
T_reac = 298.15 #Initial temperature in Kelvin
P = 101325 #Inital pressure in pascals
phi = np.linspace(0.1,4/3,20) #phi must not be zero else the stoichio coeff of air would become not defined
T_aft_py = []
T_aft_ct = []
#AFT calculated using cantera
gas = ct.Solution('gri30.xml')
for i in range(len(phi)):
a = (x+y/4)/phi[i] #Stoichiometric coeff of reactant side air
n_ch4 = 1
n_o2 = a
n_n2 = a*3.76
ntot = n_ch4+n_o2+n_n2
x_ch4 = n_ch4/ntot
x_o2 = n_o2/ntot
x_n2 = n_n2/ntot
gas.TPX = T_reac,P,{'CH4':x_ch4,'O2':x_o2,'N2':x_n2}
gas.equilibrate('UV','auto') #Constant Volume Combustion
T_aft_ct.append(gas.T)
#AFT Calculated using first law principles only
#function to calucalte enthalpy at a given temp using NASA Polynomial coefficients
def h(T,coeffs):
a1 = coeffs[0]
a2 = coeffs[1]
a3 = coeffs[2]
a4 = coeffs[3]
a5 = coeffs[4]
a6 = coeffs[5]
a7 = coeffs[6]
R = 8.314 #J/molK
return (a1 + a2*T/2 + a3*pow(T,2)/3 + a4*pow(T,3)/4 + a5*pow(T,4)/5 + a6/T)*R*T #J/mol
#NASA Polynomial coefficient from GRI Mech file. Special Attention must be paid while copying this data from the GRI mech file
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]
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]
h2o_coeff_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]
co_coeff_h = [2.71518561E+00,2.06252743E-03,-9.98825771E-07,2.30053008E-10,-2.03647716E-14,-1.41518724E+04,7.81868772E+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]
#Newtin Raphson Method Functions f(T) and fprime(T) for root finding problem for Constant volume
#f(T) = H_p-H_r - R*(N_reac*T_reac - N_prod*T_ad)
def f(T,x,y,phi):
h_ch4_r = h(T_reac,ch4_coeff_l)
h_o2_r = h(T_reac,o2_coeff_l)
h_n2_r = h(T_reac,n2_coeff_l)
a = (x+y/4)/phi
H_r = h_ch4_r + a*h_o2_r + 3.76*a*h_n2_r
N_reac = 1+a+(3.76*a)
if (phi<=1):
h_co2_p = h(T,co2_coeff_h)
h_h2o_p = h(T,h2o_coeff_h)
h_n2_p = h(T,n2_coeff_h)
h_o2_p = h(T,o2_coeff_h)
b = x
c = y/2
d = ((1-phi)/phi)*(x+(y/4))
e = 3.76*a
H_p = b*h_co2_p + c*h_h2o_p + d*h_o2_p + e*h_n2_p
N_prod = b+c+d+e
else:
h_co2_p = h(T,co2_coeff_h)
h_h2o_p = h(T,h2o_coeff_h)
h_n2_p = h(T,n2_coeff_h)
h_co_p = h(T,co_coeff_h)
d = y/2
b = 2*((x+y/4)/phi)-(y/2)-x
c = 2*x + (y/2) - 2*((x+y/4)/phi)
e = 3.76*a
H_p = b*h_co2_p + c*h_co_p + d*h_h2o_p + e*h_n2_p
N_prod = b+c+d+e
R = 8.314 #J/molK
return H_p-H_r-R*(N_prod*T-N_reac*T_reac)
def fprime(T,x,y,phi):
return (f(T+1e-6,x,y,phi)-f(T,x,y,phi))/1e-06
T_guess = 1500
ct = 0
tol = 1e-06
alpha = 0.5
for j in range(len(phi)):
while(abs(f(T_guess,x,y,phi[j]))>tol):
T_guess = T_guess - alpha*(f(T_guess,x,y,phi[j])/fprime(T_guess,x,y,phi[j]))
T_aft_py.append(T_guess)
"""
print('Adiabatic Flame Temperature = ',T_guess)
print('Number of iterations = ',ct)
"""
plt.plot(phi,T_aft_ct,'r-o')
plt.plot(phi,T_aft_py,'b-o')
plt.xlabel('Equivalence ratio')
plt.ylabel('Adiabatic Flame temperature in Kelvin')
plt.xticks(phi,rotation=45)
plt.legend(["Using Cantera","Solved the AFT equation in Python"])
plt.title('AFT vs Equivalence ratio for Methane-Air combustion in a Constant volume process')
plt.grid('on')
plt.show()CODE FOR AFT AS A FUNCTION OF HEAT LOSS FOR DIFFERENT HCs(PART 2 OF PROBLEM STATMENT)
https://drive.google.com/file/d/1P3PbbUKqPI7YsEiQiYlcBZcDJrG5Ljlw/view?usp=sharing
"""
AFT Calculator for a Constant pressure system with heat loss as a function of LHV for a general HC CxHy
The code has been tested for 2 arbon atom alkane alkene and alkynes.
Stoichiometric reactions have only been considered
LHV Values have been taken from Book by SR Turns
THIS IS A VERY GENERIC CODE AND CAN BE MODIFIED EASILIY FOR ANY HYDROCARBON WITH A FEW CHANGES
AUTHOR : Sharique Iqbal
Date : 25 Oct 2020
"""
import math
import matplotlib.pyplot as plt
import cantera as ct
import numpy as np
#CxHy -> Arbitrary Hydrocarbon
x = [1,2,3] #Carbon atoms
y = [4,6,8] #Corresponding Hydrogen atoms
T_reac = 298.15 #Initial temperature in Kelvin
P = 101325 #Inital pressure in pascals
H_loss = 0.0 #Heat loss factor. For 0 use 0.0.This must have length of 1 if this code is to be modifed for multiple hydrocarbons
check = isinstance(H_loss,float) #Len function doesn't work for float type objects. This is a way around.
#function to calucalte enthalpy at a given temp using NASA Polynomial coefficients
def h(T,coeffs):
a1 = coeffs[0]
a2 = coeffs[1]
a3 = coeffs[2]
a4 = coeffs[3]
a5 = coeffs[4]
a6 = coeffs[5]
a7 = coeffs[6]
R = 8.314 #J/molK
return (a1 + a2*T/2 + a3*pow(T,2)/3 + a4*pow(T,3)/4 + a5*pow(T,4)/5 + a6/T)*R*T #J/mol
#NASA Polynomial coefficient from GRI Mech file
ch4_coeff_l = [5.14987613E+00,-1.36709788E-02,4.91800599E-05,-4.84743026E-08,1.66693956E-11,-1.02466476E+04,-4.64130376E+00]
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]
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]
c2h4_coeff_l = [3.95920148E+00,-7.57052247E-03,5.70990292E-05,-6.91588753E-08,2.69884373E-11,5.08977593E+03,4.09733096E+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]
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]
h2o_coeff_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]
#Functions f(T) and fprime(T) for root finding problem
#f(T) = H_p-H_r+H_loss*LHV
#With heat loss H_p+Hloss = Hr
#LHV values have been taken from SR Turns. Their units were kJ/kg. Converted to J/mol by multiplying with MW of the fuel
if(check==1):
def f(T,x,y,H_loss):
a = x+y/4
MW = 12*x + 1*y #MW of the fuel in g/mol
if(y==((2*x)+2)): #Check for alkane
if (x==1):
h_fuel = h(T_reac,ch4_coeff_l)
LHV = 50016*MW
elif(x==2):
h_fuel = h(T_reac,c2h6_coeff_l)
LHV = 47489*MW
else:
h_fuel = h(T_reac,c3h8_coeff_l)
LHV = 46357*MW
if(y==2*x): #Check for alkene. Only one alkyne has been considered i.e. C2H4
h_fuel = h(T_reac,c2h4_coeff_l)
LHV = 47161*MW
if(y==((2*x)-2)): #Check for alkyne.Only acetylene has been considered
h_fuel = h(T_reac,c2h2_coeff_l)
LHV = 48225*MW
h_o2_r = h(T_reac,o2_coeff_l)
h_n2_r = h(T_reac,n2_coeff_l)
H_r = h_fuel + a*h_o2_r + 3.76*a*h_n2_r
#Enthalpy of products
h_co2_p = h(T,co2_coeff_h)
h_h2o_p = h(T,h2o_coeff_h)
h_n2_p = h(T,n2_coeff_h)
H_p = x*h_co2_p + (y/2)*h_h2o_p + 3.76*a*h_n2_p
return H_p-H_r+H_loss*LHV
def fprime(T,x,y,H_loss):
return (f(T+1e-6,x,y,H_loss)-f(T,x,y,H_loss))/1e-06
T_guess = 1500
ctr = 0
tol = 1e-06
alpha = 0.8
T_aft = []
if (len(y)>1 and len(x)==1):
for i in range(len(y)):
while (abs(f(T_guess,x,y[i],H_loss))>tol):
T_guess = T_guess - alpha*(f(T_guess,x,y[i],H_loss)/fprime(T_guess,x,y[i],H_loss))
T_aft.append(T_guess)
elif(len(x)>1 and len(y)>1):
for i in range(len(x)):
while (abs(f(T_guess,x[i],y[i],H_loss))>tol):
T_guess = T_guess - alpha*(f(T_guess,x[i],y[i],H_loss)/fprime(T_guess,x[i],y[i],H_loss))
T_aft.append(T_guess)
elif(len(x)>1 and len(y)==1):
for i in range(len(x)):
while (abs(f(T_guess,x[i],y,H_loss))>tol):
T_guess = T_guess - alpha*(f(T_guess,x[i],y,H_loss)/fprime(T_guess,x[i],y,H_loss))
T_aft.append(T_guess)
else:
while (abs(f(T_guess,x,y,H_loss))>tol):
T_guess = T_guess - alpha*(f(T_guess,x,y,H_loss)/fprime(T_guess,x,y,H_loss))
T_aft = T_guess
x = ['Methane(CH4)','Ethane(C2H6)','Ethylene(C3H8)'] #Need to be changed dependin on x and y values
plt.plot(x,T_aft,'ro-')
print('Adiabatic Flame Temperature = ',T_aft)
plt.xlabel('Hydrocarbon Name')
plt.ylabel('Adiabatic Flame Temperature in Kelvin')
plt.title('Adiabatic Flame temperature for different HCs in a constant pressure combustion without heat loss')
plt.grid('on')
plt.show()
elif(len(H_loss)>1):
def f(T,x,y,H_loss):
a = x+y/4
MW = 12*x + 1*y #MW of the fuel in g/mol
if(y==((2*x)+2)): #Check for alkane
if (x==1):
h_fuel = h(T_reac,ch4_coeff_l)
LHV = 50016*MW
elif(x==2):
h_fuel = h(T_reac,c2h6_coeff_l)
LHV = 47489*MW
else:
h_fuel = h(T_reac,c3h8_coeff_l)
LHV = 46357*MW
if(y==2*x): #Check for alkene. Only one alkyne has been considered i.e. C2H4
h_fuel = h(T_reac,c2h4_coeff_l)
LHV = 47161*MW
if(y==((2*x)-2)): #Check for alkyne.Only acetylene has been considered
h_fuel = h(T_reac,c2h2_coeff_l)
LHV = 48225*MW
h_o2_r = h(T_reac,o2_coeff_l)
h_n2_r = h(T_reac,n2_coeff_l)
H_r = h_fuel + a*h_o2_r + 3.76*a*h_n2_r
#Enthalpy of products
h_co2_p = h(T,co2_coeff_h)
h_h2o_p = h(T,h2o_coeff_h)
h_n2_p = h(T,n2_coeff_h)
H_p = x*h_co2_p + (y/2)*h_h2o_p + 3.76*a*h_n2_p
return H_p-H_r+H_loss*LHV
def fprime(T,x,y,H_loss):
return (f(T+1e-6,x,y,H_loss)-f(T,x,y,H_loss))/1e-06
T_guess = 1500
ctr = 0
tol = 1e-06
alpha = 0.8
T_aft = []
for k in range(len(H_loss)):
while (abs(f(T_guess,x,y,H_loss[k]))>tol):
T_guess = T_guess - alpha*(f(T_guess,x,y,H_loss[k])/fprime(T_guess,x,y,H_loss[k]))
T_aft.append(T_guess)
plt.plot(H_loss,T_aft,'ro-')
print('Adiabatic Flame Temperature = ',T_aft)
plt.xlabel('Heat Loss Factor')
plt.ylabel('Flame Temperature in Kelvin')
plt.title('Flame temperature of CH4 in constant pressure combustion with heat loss')
plt.xticks(H_loss)
plt.grid('on')
plt.show()