Adiabatic Flame Temperature calculation   (Cantera, Skill Lync, Easy)
Aim:
1) To calculate adiabatic flame temperature for given species in constant volume combustion chamber.
2) To calculate adiabatic flame temperature for given species in constant pressure chamber (considering effect of product at initial reactant temperature).
3) Observing the effect of bond strength and number of carbon atoms on AFT.
Introduction:
For designing any combustion chamber, whether a constant volume type or a constant pressure (steady flow) type as in gas turbines, maximum temperature attained has to kept below keeping in mid the metallurgical considerations of parts in chamber. The maximum temperature can be attained only when all of heat released is used to increase temperature of products. This maximum temperature is called Adiabatic Flame temperature. This temperature is dependent on the air fuel ratio, type of fuel used and type of combustion chamber.
Following report has 4 cases:
1) Calculating AFT in Methane-Air combustion for different equivalence ratios, in a constant volume chamber using Python and Cantera
2) Calculating AFT in Methane-Air combustion in a constant pressure chamber. (considering effect of product at initial reactant temperature).
3) Effect of -ane, -ene, -yne variants of a hydrocarbon chain (Ethane, Ethene, Ethyne) on AFT.
4) Effect of number of carbon atoms in hydrocarbon chain (Methane, Ethane, Propane) on AFT.
Theory:
1) Moles of reactants and products for different combustion equivalence ratios:
Assuming no dissociation, the combustion of an arbitrary hydrocarbon with our simplified air can be represented as:
Equivalence Ratio for any combustion is:
For lean or stoichiometric conditions (φ<=1)
For rich conditions (φ>1)
According to calculation shown in Ref [1]:
For Lean and Stoichiometric condition:
Since there is no product species that would dissociate leading to chemical equilibrium conditions, the calculation of number of moles is pretty straightforward method of balancing the atoms.
For Rich conditions:
In these conditions, the CO2 molecules have tendency to react with H2 to form CO and H2O. This reaction is reversible and thus CO and H2O can also form CO2 and H2. Due to this, a chemical equilibrium condition is used to evaluate the number of moles. The equilibrium constant Kp is calculated using the relation given below.
The values of Kp can be found using equation:`Kp=0.3754*x^(-0.719)` (Ref [1])
2) Constant Volume chamber
The volume of chamber before and after combustion remains same. So, we apply first law of thermodynamics:
Using Ideal Gas Law:
Since there is no path for heat to be released, the H20 vapours will not condense. Thus, we will only get Lower Heating Value (LHV)
3) Constant Pressure Chamber
This case is useful for combustion gas turbines. The flow of air and fuel is continuous and there change in pressure due to combustion and expansion. The governing equation is:
In a constant pressure chamber, some part of heat released is used to perform work of expansion of volume to keep pressure constant. This is the reason AFT is less in constant pressure chamber than in constant volume chamber.
Higher and lower heating values:
If during combustion the reactants at STP are burned to form products, there is release of heat. This heat value is called Lower Heating Value. When a gas is cooled, it releases heat. If the gaseous products are cooled, we can release more heat along with LHV. This total heat released is called Higher Heating Value (HHV). This heat is captured in preheating or can be captured by achieving more expansion.
The maximum heat loss in a combustion reaction is HHV. In case 2 of report, the plot of AFT is obtained for different proportions of cooling of products.
4) Newton Raphson Solver
The equations written above are solved using Newton Raphson solver. This solver starts from a guess value and using curve’s slope at guess value it guesses new value which has less curve slope. In this way we reach at a point which is a Minima. This iterative process continues till the difference in slope becomes less than tolerance provided.
How does the code work?
1) Defining all the constants needed for finding enthalpy of a species at a given temperature. (These constants are taken from NASA’s Thermodynamic data files)
2) Then the values of equivalence ratio for which AFT is to be found are placed in a variable ‘phi’
3) The Newton Raphson solver starts using the parameters given. It calculates the LHS of constant volume chamber’s equation and tries to equate it to 0.
4) The LHS is calculated by function ‘E_N’. It takes in values of phi and temperature for which it is asked to solve by the solver. The function is able to calculate for different species of hydrocarbon, even if it has multiple carbon atoms or double or triple bond.
5) The enthalpy calculations are solved using enthalpy values of different species. This enthalpy is calculated by function ‘h’. It uses the NASA polynomial and coefficients to calculate the enthalpy.
6) The slope of the curve at given temperature is found using function ‘E_N_prime’.
7) Using all these values the solver finds the AFT for a single equivalence ratio. This whole procedure is repeated for all the values required.
8) Similar code is used for all cases with minor changes.
For case 3 & 4: The x and y values in function 'E_N' are changed according to the fuel used and 'H_loss' value is kept as 0. The value of 'T_hloss' at index 0 is used as AFT.
Code 1: (for case 1 of constant volume chamber)
#Program to calculate AFT for methane combustion in constant volume chamber for different equivalence ratios
# CxHy + a(O2 + 3.76N2) --> bCO2 + cCO + dH2O + eH2 + fO2 + 3.76aN2 (General Reaction)
# Date: 22/05/2020
import matplotlib.pyplot as plt
import math
Tstp = 298.15 #K
#function to return the value of enthalpy in J/mol for a given species at a given temperature
def h(T, co_effs):
a_co = co_effs
enthal_py = (a_co[0] + (a_co[1]*T)/2 + (a_co[2]*(T**2))/3 + (a_co[3]*(T**3))/4 + (a_co[4]*(T**4))/5 + (a_co[5]/T)) * 8.314 * T
return enthal_py
# Low temperature 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]
# High Temperature coefficients
co2_coeffs_h = [3.85746029E+00, 4.41437026E-03, -2.21481404E-06, 5.23490188E-10, -4.72084164E-14, -4.87591660E+04, 2.27163806E+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]
h2o_coeffs_h = [3.03399249E+00, 2.17691804E-03, -1.64072518E-07, -9.70419870E-11, 1.68200992E-14, -3.00042971E+04, 4.96677010E+00]
h2_coeffs_h = [3.33727920E+00, -4.94024731E-05, 4.99456778E-07, -1.79566394E-10, 2.00255376E-14, -9.50158922E+02, -3.20502331E+00]
o2_coeffs_h = [3.28253784E+00, 1.48308754E-03, -7.57966669E-07, 2.09470555E-10, -2.16717794E-14, -1.08845772E+03, 5.45323129E+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]
# Function to find enthalpy in J/mol for a given temperature and equivalence ratio
def E_N(T, phi):
x = 1
y = 4
if phi<=1:
a = (x + (y/4))/phi
b = x
c = 0
d = y/2
e = 0
f = (1-phi) * a
g = 3.76*a
elif phi>1:
a = (x + (y/4))/phi
kp = 0.3754 * (T**(-0.719))
b = (2*a*(kp-1) +x+ (y/2))/(2*(kp-1))
t1 = (2*a*(kp-1) +x+ (y/2))**2
t2 = 4*kp*(kp-1)*(2*a*x - x**2)
b = b - ((t1-t2)**0.5)/(2*(kp-1))
c = x - b
d = 2*a -b-x
e = -2*a +b+x+ (y/2)
f = 0
g = 3.76*a
ch4_h_r = h(Tstp, ch4_coeffs_l)
o2_h_r = h(Tstp, o2_coeffs_l)
n2_h_r = h(Tstp, n2_coeffs_l)
H_reac = ch4_h_r + (2/phi)*o2_h_r + (7.52/phi)*n2_h_r # Enthalpy of reactants at STP
NRT_reac = (1 + (2/phi) + (7.52/phi))*8.314*Tstp # NRT term of reactants at STP
co2_h_p = h(T, co2_coeffs_h)
co_h_p = h(T, co_coeffs_h)
h2o_h_p = h(T, h2o_coeffs_h)
h2_h_p = h(T, h2_coeffs_h)
o2_h_p = h(T, o2_coeffs_h)
n2_h_p = h(T, n2_coeffs_h)
H_prod = b*co2_h_p + c*co_h_p + d*h2o_h_p + e*h2_h_p + f*o2_h_p + g*n2_h_p # Enthalpy of products at Temperature T
NRT_prod = (b+c+d+e+f+g)*8.314*T # NRT term of products at temperature T
return (H_prod - H_reac) - (NRT_prod - NRT_reac)
# Defining the derivative of the function
def E_N_prime(T, phi):
return (E_N(T+1e-4, phi) - E_N(T, phi))/(1e-4)
# Setting values of equivalence ratio to find
alpha = 1
tol = 1e-4
phi_min = 0.6
phi_max = 1.4
phi_diff = 0.05
phi_num = int(round((phi_max - phi_min)/phi_diff))
phi = []
for count_er in range(0,phi_num+1):
phi.append(round((((count_er)*phi_diff)+phi_min), 4))
# Newton-Raphson solver
AFT_i = []
for i in range(0,len(phi)):
T_guess = 2000 # K
while abs(E_N(T_guess, phi[i])) >= tol:
T_guess = T_guess - alpha*((E_N(T_guess, phi[i]))/(E_N_prime(T_guess, phi[i])))
AFT_i.append(T_guess)
print(AFT_i)
#Solution using Cantera
import cantera as ct
gas = ct.Solution('gri30.xml')
AFT_cant = []
for i in range(0,len(phi)):
gas.TPX = 298.15, 101325, {'CH4':1, 'O2':2/phi[i], 'N2':7.52/phi[i]}
gas.equilibrate('UV','auto')
AFT_cant.append(gas.T)
print(AFT_cant)
plt.plot(phi,AFT_i,'-o',color='red',label='Python')
plt.plot(phi,AFT_cant,'-s',color='blue',label='Cantera')
plt.xlabel('Equivalence Ratio (phi)', fontsize=10)
plt.ylabel('Adiabatic Flame Temperature [K]', fontsize=10)
plt.grid('on')
plt.minorticks_on()
plt.grid(which ='minor',alpha = 0.2)
plt.legend(fontsize=12)
plt.title('AFT vs Equivalence Ratio for Methane-Air Combustion' 'n' 'in constant volume chamber', fontweight='bold')
plt.show()Code 2: (For case 2, 3, 4)
import matplotlib.pyplot as plt
import math
Tstp = 298.15 #K
phi = 1
#function to return the value of enthalpy in J/mol for a given species at a given temperature
def h(T, co_effs):
a_co = co_effs
enthal_py = (a_co[0] + (a_co[1]*T)/2 + (a_co[2]*(T**2))/3 + (a_co[3]*(T**3))/4 + (a_co[4]*(T**4))/5 + (a_co[5]/T)) * 8.314 * T
return enthal_py
# Low temperature 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]
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]
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
co2_coeffs_h = [3.85746029E+00, 4.41437026E-03, -2.21481404E-06, 5.23490188E-10, -4.72084164E-14, -4.87591660E+04, 2.27163806E+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]
h2o_coeffs_h = [3.03399249E+00, 2.17691804E-03, -1.64072518E-07, -9.70419870E-11, 1.68200992E-14, -3.00042971E+04, 4.96677010E+00]
h2_coeffs_h = [3.33727920E+00, -4.94024731E-05, 4.99456778E-07, -1.79566394E-10, 2.00255376E-14, -9.50158922E+02, -3.20502331E+00]
o2_coeffs_h = [3.28253784E+00, 1.48308754E-03, -7.57966669E-07, 2.09470555E-10, -2.16717794E-14, -1.08845772E+03, 5.45323129E+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]
# Function to find enthalpy in J/mol for a given temperature and equivalence ratio
def E_N(T, phi,hloss):
x = 1
y = 4
if phi<=1: # for lean and stoichiometric mixtures
a = (x + (y/4))/phi
b = x
c = 0
d = y/2
e = 0
f = (1-phi) * a
g = 3.76*a
elif phi>1: # for rich mixtures
a = (x + (y/4))/phi
kp = 0.3754 * (T**(-0.719))
b = (2*a*(kp-1) +x+ (y/2))/(2*(kp-1))
t1 = (2*a*(kp-1) +x+ (y/2))**2
t2 = 4*kp*(kp-1)*(2*a*x - x**2)
b = b - ((t1-t2)**0.5)/(2*(kp-1))
c = x - b
d = 2*a -b-x
e = -2*a +b+x+ (y/2)
f = 0
g = 3.76*a
if (x == 1) and (y == 4): # for selecting Methane as fuel
fuel_h_r = h(Tstp, ch4_coeffs_l)
if (x==2):
if (y == 2): # for selecting Ethyne as fuel
fuel_h_r = h(Tstp, c2h2_coeffs_l)
if (y == 4): # for selecting Ethene as fuel
fuel_h_r = h(Tstp, c2h4_coeffs_l)
if (y == 6): # for selecting Ethane as fuel
fuel_h_r = h(Tstp, c2h6_coeffs_l)
if (x == 3) and (y==8): # for selecting Propane as fuel
fuel_h_r = h(Tstp, c3h8_coeffs_l)
o2_h_r = h(Tstp, o2_coeffs_l)
n2_h_r = h(Tstp, n2_coeffs_l)
H_reac_STP = fuel_h_r + (2/phi)*o2_h_r + (7.52/phi)*n2_h_r # Enthalpy of reactants at STP
co2_h_p_stp = h(Tstp, co2_coeffs_h)
co_h_p_stp = h(Tstp, co_coeffs_h)
h2o_h_p_stp = h(Tstp, h2o_coeffs_h)
h2_h_p_stp = h(Tstp, h2_coeffs_h)
o2_h_p_stp = h(Tstp, o2_coeffs_h)
n2_h_p_stp = h(Tstp, n2_coeffs_h)
H_prod_STP = b*co2_h_p_stp + c*co_h_p_stp + d*h2o_h_p_stp + e*h2_h_p_stp + f*o2_h_p_stp + g*n2_h_p_stp # Enthalpy of products at STP
co2_h_p = h(T, co2_coeffs_h)
co_h_p = h(T, co_coeffs_h)
h2o_h_p = h(T, h2o_coeffs_h)
h2_h_p = h(T, h2_coeffs_h)
o2_h_p = h(T, o2_coeffs_h)
n2_h_p = h(T, n2_coeffs_h)
H_prod_AFT = b*co2_h_p + c*co_h_p + d*h2o_h_p + e*h2_h_p + f*o2_h_p + g*n2_h_p # Enthalpy of products at Temp T
return (H_prod_AFT - H_reac_STP) - hloss*(H_prod_STP - H_reac_STP)
# Defining the derivative of the function
def E_N_prime(T, phi, hloss):
return (E_N(T+1e-4, phi, hloss) - E_N(T, phi, hloss))/(1e-4)
# Setting values of heat loss to find
alpha = 1
tol = 1e-4
H_Loss_min = 0
H_Loss_max = 1
H_Loss_diff = 0.05
H_Loss_num = int(round((H_Loss_max - H_Loss_min)/H_Loss_diff))
H_Loss = []
for count_er in range(0,H_Loss_num+1):
H_Loss.append(round((((count_er)*H_Loss_diff)+H_Loss_min), 4))
# Newton-Raphson solver
T_hloss = []
for i in range(0,len(H_Loss)):
T_guess = 2000 # K
while abs(E_N(T_guess, phi, H_Loss[i])) >= tol:
T_guess = T_guess - alpha*((E_N(T_guess, phi, H_Loss[i]))/(E_N_prime(T_guess, phi, H_Loss[i])))
T_hloss.append(T_guess)
plt.plot(H_Loss,T_hloss,'-o',color='red')
plt.xlabel('Heat loss', fontsize=10)
plt.ylabel('Adiabatic Flame Temperature [K]', fontsize=10)
plt.grid('on')
plt.minorticks_on()
plt.grid(which ='minor',alpha = 0.2)
plt.title('AFT vs heat loss for Methane-Air Combustion' 'n' 'in constant pressure chamber', fontweight='bold')
plt.show()
print(T_hloss[0])Results:
Case 1: Combustion of Methane in air in a constant volume chamber for different equivalence ratios:
As it can be seen from the graph, the maximum AFT is achived at stiochiometric condition if python program is used. The cantera program gives maximum AFT at some point in between 1.05 and 1.1. This is because cantera uses an approach of minimization of gibbs free energy, which gives us a practical value of AFT. It can also be stated that since cantera uses this approach, the maximum AFT is lower in cantera program than in python program.
2) Combustion of methane in air in constant pressure chamber (considering effect of product at initial reactant temperature):
The AFT value on the left of the graph is obtained when the products are at AFT, i.e. it is value when LHV is heat output of reaction.
The AFT value on the right of the graph is obtained when all of the products are at reactant temperature, i.e. it is value when HHV is heat output of reaction.
3) Effect of type of bond on AFT:
The AFT increases as the bond strength between carbon increases. This is due to more energy released when higher order bonds break. The temperature difference is quite distinct (~500K) between ethane and ethyne.
4) Effect of number of carbon atoms in carbon chain on AFT:
The AFT increases as the carbon chain length increases. The increase is small (~70K) between methane and propane. The values approximately match with Ref [2] and Ref [3]. These values of Ref [2] and Ref [3] are less because they used cantera code to find the temperatures.
Refrences:
1) "An Introduction to Combustion" by Stephen Turns.
2) http://elearning.cerfacs.fr/combustion/tools/adiabaticflametemperature/index.php
3) https://www.engineeringtoolbox.com/adiabatic-flame-temperature-d_996.html