ADIABATIC FLAME TEMPERATURE CALCULATION FOR COMBUSTION REACTION OF HYDROCARBONS AS FUEL IN CONSTANT PRESSURE AND CONSTANT VOLUME SYSTEMS

INTRODUCTION

In the study of combustion, there are two types of adiabatic flame temperature (AFT)

depending on how the process is completed: the constant volume and constant

pressure; both of which describe temperature that combustion products theoretically can

reach if no energy is lost to the outside environment.

The AFT is calculated from the enthalpies of reactants and products. We use the Newton

Raphson iterative solver to calculate the AFT by equating the enthalpies of the reactants

and the products in a chemical reaction.

The enthalpies of the reactants and products are calculated by the NASA polynomial

functions.

SCOPE OF CURRENT STUDY

We investigate four different cases for the variation of AFT:

Case 1: Adiabatic flame temperature of methane-air combustion in a constant volume

chamber by varying the equivalence ratio.

Case 2: Flame temperature of methane-air combustion under constant pressure in the

presence of a heat loss.

Case 3: Adiabatic flame temperature for general alkane-air combustion with 1, 2 and

Carbon atom chains.

Case 4: Adiabatic flame temperature of -ane, -ene, and -yne variants of Carbon atom chains

combustion in air.

The enthalpy of the reactants and products at a particular temperature are calculated using

the NASA polynomial function as:

`H(T)=(a1+a2*T/2+a3*T^2/3+a4*T^3/4+a5*T^4/5+ a6*T/6)RT`

, where `a_i` are the coefficients of each component at either low or high temperatures

obtained from the GRI-Mech 3.0 Thermodynamics file. Generally, the reactant coefficients are

taken at low temperatures (assuming STP) and product coefficients are taken at high

temperatures.

The AFT is calculated by finding the root of the equilibrium condition. The Newton-Raphson

method is an iterative method to find the root of a function by starting at a guess value for

and iterating till the value of the function `f(x)` reaches below a tolerance as:

NOTE: `x_n` stands for the new value of \'x\'

`x_n =x_(guess) - α*[f(x_(guess) )/(f^\' (x_(guess))]]`

Here, `alpha` is the relaxation factor which determines the intensity of the iterative solver.

For alpha < 1` the solution is under-relaxed and may take slightly a greater number of

iterations than higher values, but it also tries to avoid divergence in a solution. For `alpha >

1` the solution is rigorous and may lead to divergence.

CASE 1 - VARIATION W.R.T. EQUIVALENCE RATIO OF REACTION at CONSTANT VOLUME:

The equivalence ratio is defined as the ratio of fuel-to-air ratio in the mixture to the

stoichiometric fuel-to-air ratio. Depending on the equivalence ratio, the combustion may take

place in a controlled manner or may result in an explosion. Consider the stoichiometric

combustion of methane and air:

`CH_4+2(O_2+3.76N_2 )→CO_2+2H_2 O+7.52N_2`

Here the stoichiometric fuel-to-air ratio is 1/2. For an equivalence ratio of , the reaction can

be rewritten as:

`CH_4+(2/phi)(O_2+3.76N_2 )→xCO_2+yH_2 O +` other products

The products can also differ with different equivalence ratios. Typically, for `phi > 1`, the

mixture is termed as a rich mixture containing an excess fuel and for `phi < 1` , the mixture

is termed as a lean mixture containing an excess of air.

Rich Mixture `(phi > 1)`

`CH_4+(2/phi)(O_2+3.76N_2 )→(4/phi-3)CO_2+2H_2 O+(7.52/phi) N_2 +(4-4/phi)CO`

Lean Mixture `(phi < 1)`

`CH_4+(2/phi)(O_2+3.76N_2 )→CO_2+2H_2 O+(7.52/phi) N_2 +(2/phi-2) O_2`

Since, the reaction is occurring in a constant volume combustion chamber, in addition to

enthalpy balance, the number of moles also need to be balanced as:

`〖(H〗_products-H_reactants)-R(N_products T-N_reactants T_STP)=0`

At STP, the temperature is 298.15 K. This temperature is used to calculate the enthalpies of

the reactants and the equilibrium temperature for the products is calculated iteratively

through the Newton-Raphson method to find the AFT. The relaxation parameter for the

Newton-Raphson method is set as `alpha = 1` and the tolerance value is set as `1 *

10^-3`.

The solution from Python programming is compared to the solution obtained by Cantera. The

AFT obtained by Cantera is usually lower than the AFT calculated by Python as Cantera

minimizes the Gibbs\' Free Energy to obtain equilibrium. The solution obtained by Cantera is

more accurate than the solution obtained by Python even though there is not muchdifference

between the two. The mixture in Cantera is equilibrated at constant internal energy and

volume to obtain the AFT.

CODE to calculate AFT for case 1 in Python and Cantera

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

\"\"\"
Program to study the dependence of Adiabatic Flame Temperature on Equivalence Ratio
Assumption: Constant Volume System
Stoichiometric (phi = 1): CH4 + 2(O2 + 3.76N2) <=> CO2 + 2H2O + 7.52N2
Rich Mixture (phi > 1): CH4 + (2/phi)(O2 + 3.76N2) <=> (4/phi - 3)CO2 + 2H2O + (7.52/phi)N2 + (4 - 4/phi)CO
Lean Mixture (phi < 1): CH4 + (2/phi)(O2 + 3.76N2) <=> CO2 + 2H2O + (7.52/phi)N2 + (2/phi - 2)O2
Program by: Shouvik Bandopadhyay
\"\"\"

R = 8.314 # J/mol-K
# Function to calculate enthalpy of a compound through NASA polynomial function
def h(T, coeff):
	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

# Coefficients for the NASA polynomial for low temperatures
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]

# Coefficients for the NASA polynomial for high temperatures
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]

# Standard temperature condition
Tstp = 298.15 # K
# Defining the function for root-solving by Newton-Raphson method
def f(T, phi):
	h_ch4_r = h(Tstp, ch4_coeff_l)
	h_o2_r = h(Tstp, o2_coeff_l)
	h_n2_r = h(Tstp, n2_coeff_l)
	
	H_reactants = h_ch4_r + (2/phi)*h_o2_r + (7.52/phi)*h_n2_r
	N_reactants = (1 + (2/phi) + (7.52/phi))*R*Tstp
	
	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)
	h_o2_p = h(T, o2_coeff_h)
	
	if phi >= 1:
		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_products = ((4/phi - 3) + 2 + (7.52/phi) + (4-4/phi))*R*T

	elif phi < 1:
		H_products = h_co2_p + 2*h_h2o_p + (7.52/phi)*h_n2_p + (2/phi-2)*h_o2_p
		N_products = (1 + 2 + (7.52/phi) + (2/phi-2))*R*T

	return (H_products - H_reactants) - (N_products - N_reactants)

# Defining the derivative of the function
def fprime(T, phi):
	return (f(T+1e-6, phi) - f(T, phi))/(1e-6)

# Newton-Raphson solver
alpha = 1
tol = 1e-3
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, 1.5, 1.6, 1.7, 1.8, 1.9, 2]
T_aft_py = []

for i in range(0,len(phi)):
	T_guess = 1500 # K
	while abs(f(T_guess, phi[i])) >= tol:
		T_guess = T_guess - alpha*((f(T_guess, phi[i]))/(fprime(T_guess, phi[i])))
	T_aft_py.append(T_guess)

print(T_aft_py)
print()

\"\"\"
Solution obtained by Cantera
\"\"\"

gas = ct.Solution(\'gri30.xml\')
T_aft_ct = []

for i in range(0,len(phi)):
# Computing the mole fractions of each reactant
	n_tot = 1 + (2/phi[i]) + (7.52/phi[i])
	n_ch4 = 1/n_tot
	n_o2 = (2/phi[i])/n_tot
	n_n2 = (7.52/phi[i])/n_tot
# Defining the gas mixture

	gas.TPX = 298.15, 101325, {\'CH4\':n_ch4, \'O2\':n_o2, \'N2\':n_n2}
# Calculating the equilibrium conditions

	gas.equilibrate(\'UV\',\'auto\')
	T_aft_ct.append(gas.T)

print(T_aft_ct)

# Plotting the results
plt.plot(phi,T_aft_py,\'-o\',color=\'green\',label=\'Python\')
plt.plot(phi, T_aft_ct, \'-o\',color=\'yellow\',label=\'Cantera\')
plt.xlabel(\'Equivalence Ratio ($phi$)\', fontsize=12)
plt.ylabel(\'Adiabatic Flame Temperature [K]\', fontsize=12)
plt.yticks(np.arange(600, 3000, step=200))
plt.xticks(np.arange(0, 2.1, step=0.2))
plt.grid(\'on\')
plt.title(\'AFT vs. Equivalence Ratio for Methane-Air Combustion\', fontsize=13, fontweight=\'bold\')
plt.legend(fontsize=12)
plt.show()

Output from the above code (Equivalence ratios and Corrosponding AFTs)

`phi` PYTHON CANTERA

0.1 681.68 679.62

0.2 1009.55 1009.43

0.3 1301.87 1300.75

0.4 1567.86 1563.56

0.5 1812.86 1802.09

0.6 2040.42 2018.37

0.7 2252.98 2211.72

0.8 2452.46 2377.51

0.9 2640.34 2506.67

1.0 2817.83 2585.88

1.1 2704.17 2600.54

1.2 2592.95 2556.49

1.3 2484.01 2484.06

1.4 2377.13 2403.26

1.5 2272.15 2321.03

1.6 2168.86 2239.55

1.7 2067.07 2159.52

1.8 1966.59 2081.14

1.9 1867.22 2004.48

2.0 1768.75 1929.53

• We see a linear curve for Python data points while the curve for the Cantera data points give us the actual solution.
• AFT increases consistently from equivalence ratio 0.1 to 1 and declines after the value of 1.
• Hence it is safe to conclude that stoichiometric combustion is the most aggressive combustion reaction when compared to either lean or rich mixture combustion reaction.

CASE 2 - VARIATION W.R.T. HEAT LOSS AT CONSTANT PRESSURE:

In this case, we will be analyzing the variation of AFT with heat loss from a constant pressure

chamber. Due to the heat loss, the reaction is no longer adiabatic and hence, we will refer to

the AFT as just FT (flame temperature). The reaction is stoichiometric combustion of

methane and air:

`CH_4+2(O_2+3.76N_2 )→CO_2+2H_2 O+7.52N_2`

The heat loss from the reaction is programmed as a factor (H_loss) of the maximum possible

heat loss from the reaction. The maximum possible heat loss from the reaction is calculated

as the lower heating value (LHV) of the reaction. LHV is calculated from the difference in the

enthalpies of reactants and products when both are considered at STP. Since the `N_2` from

both the reactants and the products side is considered at `T_(STP) = 298.15 K`, it will

cancel each other out and can effectively be neglected from the calculations altogether as:

`LHV = H_(CH_4)^298.15+2H_(O_2)^298.15- H_(CO_2)^298.15-2H_(H_2 O)^298.15`

The heat loss factor (H_loss) is varied from (no heat loss) to (all heat lost) as a multiplier of

the LHV. Since the heat is lost from the reaction, lesser heat is available for the products to

form and the equilibrium condition (root-finding function) for this reaction is:

`H_(pro ducts) - H_(reactants) + H_(loss) = 0`

CODE to calculate AFT for case 2 in Python and Cantera

import matplotlib.pyplot as plt
import numpy as np

\"\"\"
Program to study the dependence of Adiabatic Flame Temperature on Heat Loss Ratio
Assumpton: Constant pressure system
Methane: CH4 + 2(O2 + 3.76N2) <=> CO2 + 2H2O + 7.52N2
Program by: Shouvik Bandopadhyay
\"\"\"

R = 8.314 # J/mol-K

# Function to calculate enthalpy of a compound through NASA polynomial function
def h(T, coeff):
	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

# Coefficients for the NASA polynomial for low temperatures
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]

# Coefficients for the NASA polynomial for high temperatures
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]

# Standard temperature condition
Tstp = 298.15 # K

# Computing the maximum possible heat loss
LHV = h(Tstp,ch4_coeff_l) + 2*h(Tstp,o2_coeff_l) - h(Tstp,co2_coeff_h) - 2*h(Tstp,h2o_coeff_h)

# Defining the function for root-solving by Newton-Raphson method
def f(T, H_loss):
	h_ch4_r = h(Tstp, ch4_coeff_l)
	h_o2_r = h(Tstp, o2_coeff_l)
	h_n2_r = h(Tstp, 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_coeff_h)
	h_n2_p = h(T, n2_coeff_h)
	H_products = h_co2_p + 2*h_h2o_p + 7.52*h_n2_p
	Loss = H_loss*LHV
	return H_products - H_reactants + Loss

# Defining the derivative of the function
def fprime(T, H_loss):
	return (f(T+1e-6, H_loss) - f(T, H_loss))/(1e-6)

# Newton-Raphson solver
alpha = 0.5
tol = 1e-3
H_loss = np.arange(0, 1.05, step=0.05)
T_aft_py = []
for i in range(0,len(H_loss)):
	T_guess = 1500 # K
	while abs(f(T_guess, H_loss[i])) >= tol:
		T_guess = T_guess - alpha*((f(T_guess, H_loss[i]))/(fprime(T_guess, H_loss[i])))
	T_aft_py.append(T_guess)
print(T_aft_py[7])

# Plotting the results
plt.plot(H_loss,T_aft_py,\'-o\',color=\'green\')
plt.plot(H_loss[7],T_aft_py[7],\'o\',color=\'black\',label=\'H_loss = 0.35\')
plt.xlabel(\'Heat Loss Ratio\', fontsize=12)
plt.ylabel(\'Flame Temperature [K]\', fontsize=12)
plt.yticks(np.arange(250, 2500, step=250))
plt.xticks(np.arange(0, 1.05, step=0.1))
plt.grid(\'on\')
plt.legend(fontsize=12)
plt.title(\'FT vs. Heat Loss Ratio for Methane-Air Combustion\', fontsize=13, fontweight=\'bold\')
plt.show()

Output from the above code

• It can be clearly seen that as the heat loss from the reaction increases, the flame temperature of the products consistently reduces as lesser heat is available
  for the reaction.

• At `H_(loss) = 0`, no heat is lost from the reaction and hence FT = AFT = 2325.6 K. At `H_(loss) = 1`, all the heat is lost from the reaction and no heat is left for the products to form, hence the flame temperature is `T_(STP) = 298.15 K`.

• From the graph, the flame temperature of methane combustion with `H_(loss) = 0.35` is calculated as 1678.81 K.

CASE 3 - VARIATION W.R.T. C-ATOM BONDS FOR N = 2 (Constant Pressure):

In this case, we will study the AFT of ethane, ethene, and ethyne to see how the AFT varies

with the type of hydrocarbon. The stoichiometric combustion reactions for the three

hydrocarbons are:

Ethane

`C_2 H_(6 )+3.5(O_2+3.76N_2 )→2CO_2+3H_2 O+13.16N_2`

Ethene

`C_2 H_(4 )+3(O_2+3.76N_2 )→2CO_2+2H_2 O+11.28N_2`

Ethyne

`C_2 H_(2 )+2.5(O_2+3.76N_2 )→2CO_2+H_2 O+9.4N_2`

The combustion takes place in a constant pressure chamber without any heat loss. Hence the

equilibrium condition (root-finding function) for all the above reactions is:

`H_(pro ducts) - H_(reactants) = 0`

CODE to calculate AFT for case 3 in Python and Cantera

import matplotlib.pyplot as plt
import cantera as ct

\"\"\"
Program to study the dependence of Adiabatic Flame Temperature on the Type of Hydrocarbon (n = 2)
General Hydrocarbon: CxHy + [x+(y/4)](O2 + 3.76N2) <=> xCO2 + (y/2)H2O + 3.76[x+(y/4)]N2
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
Program by: Shouvik Bandopadhyay
\"\"\"

R = 8.314 # J/mol-K

# Function to calculate enthalpy of a compound through NASA polynomial function
def h(T, coeff):
	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

# Coefficients for the NASA polynomial for low temperatures
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]
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]

# Coefficients for the NASA polynomial for high temperatures
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]
# Standard temperature condition
Tstp = 298.15 # K

# Defining the function for root-solving by Newton-Raphson method
def f(T, ntype):
	h_c2h6_r = h(Tstp, c2h6_coeff_l)
	h_c2h4_r = h(Tstp, c2h4_coeff_l)
	h_c2h2_r = h(Tstp, c2h2_coeff_l)
	h_o2_r = h(Tstp, o2_coeff_l)
	h_n2_r = h(Tstp, n2_coeff_l)
	h_co2_p = h(T, co2_coeff_h)	
	h_h2o_p = h(T, h2o_coeff_h)
	h_n2_p = h(T, n2_coeff_h)
	if ntype==1:
		x = 2
		y = 6
		h_fuel = h_c2h6_r
	elif ntype==2:
		x = 2
		y = 4
		h_fuel = h_c2h4_r
	elif ntype==3:
		x = 2
		y = 2
		h_fuel = h_c2h2_r
	H_reactants = h_fuel + (x+y/4)*h_o2_r + 3.76*(x+y/4)*h_n2_r
	H_products = x*h_co2_p + (y/2)*h_h2o_p + 3.76*(x+y/4)*h_n2_p
	return H_products - H_reactants

# Defining the derivative of the function
def fprime(T, ntype):
	return (f(T+1e-6, ntype) - f(T, ntype))/(1e-6)

# Newton-Raphson solver
alpha = 1
tol = 1e-3
ntype = [1, 2, 3]
T_aft_py = []
for i in range(0,len(ntype)):
	T_guess = 1500 # K
	while abs(f(T_guess, ntype[i])) >= tol:
		T_guess = T_guess - alpha*((f(T_guess, ntype[i]))/(fprime(T_guess, ntype[i])))
	T_aft_py.append(T_guess)
	print(T_aft_py)
print()

# Plotting the Python results
plt.bar([\'Ethane ($C_2H_6$)\',\'Ethene ($C_2H_4$)\',\'Ethyne ($C_2H_2$)\'],T_aft_py)
plt.ylabel(\'Adiabatic Flame Temperature [K]\',fontsize=12)
plt.title(\'AFT vs. Type of Hydrocarbon (Python)\',fontsize=13,fontweight=\'bold\')
plt.show()

\"\"\"
Solution obtained by Cantera
\"\"\"

gas = ct.Solution(\'gri30.xml\')
T_aft_ct = []
for i in range(0,len(ntype)):

# Computing the mole fractions of each reactant
	if ntype[i]==1:
		x = 2
		y = 6
		fuel = \'C2H6\'
	elif ntype[i]==2:
		x = 2
		y = 4
		fuel = \'C2H4\'
	elif ntype[i]==3:
		x = 2
		y = 2
		fuel = \'C2H2\'
	n_tot = 1 + (x+y/4) + 3.76*(x+y/4)
	n_fuel = 1/n_tot
	n_o2 = (x+y/4)/n_tot
	n_n2 = (3.76*(x+y/4))/n_tot

# Defining the gas mixture
gas.TPX = 298.15, 101325, {fuel:n_fuel, \'O2\':n_o2, \'N2\': n_n2}

# Calculating the equilibrium conditions
gas.equilibrate(\'HP\',\'auto\')
T_aft_ct.append(gas.T)
print(T_aft_ct)

# Plotting the Cantera results
plt.bar([\'Ethane ($C_2H_6$)\',\'Ethene ($C_2H_4$)\',\'Ethyne ($C_2H_2$)\'],T_aft_py)
plt.ylabel(\'Adiabatic Flame Temperature [K]\',fontsize=12)
plt.title(\'AFT vs. Type of Hydrocarbon (Cantera)\',fontsize=13,fontweight=\'bold\')
plt.show()

Output from the above code

HYDROCARBON FUEL	PYTHON AFT IN KELVIN	CANTERA AFT IN KELVIN
C2H6	2738.85	2259.74
C2H4	2565.53	2368.65
C2H2	2909.35	2538.68

• We can see that the AFT increases considerably for ethene and ethyne compared to ethanebecause of their multiple C-atom bonds.

• For multiple C-atom bonds, more energy is required to break the bond and cause the reaction, and hence higher temperatures are generated.

• Also, we can see that the difference in the solution obtained from Python and Cantera differ by more amounts as the hydrocarbon becomes more complex. This is because the combustion reaction of ethene and ethyne can produce some additional products than anticipated.

CASE 4 - VARIATION W.R.T. NUMBER OF C-ATOM CHAINS:

In this case, we will study the AFT of methane, ethane, and propane to see how the AFT

varies with the number of C-atoms in a hydrocarbon. The stoichiometric combustion

reactions for the three hydrocarbons are:

Methane

`CH_4+2(O_2+3.76N_2 )→CO_2+2H_2 O+7.52N_2`

Ethane

`C_2 H_(6 )+3.5(O_2+3.76N_2 )→2CO_2+3H_2 O+13.16N_2`

Propane

`C_3 H_(8 )+5(O_2+3.76N_2 )→3CO_2+4H_2 O+18.8N_2`

The combustion takes place in a constant pressure chamber without any heat loss. Hence the

equilibrium condition (root-finding function) for all the above reactions is:

`H_(pro ducts) - H_(reactants) = 0`

CODE to calculate AFT for case 4 in Python and Cantera

import matplotlib.pyplot as plt
import cantera as ct

\"\"\"
Program to study the dependence of Adiabatic Flame Temperature on No. of Carbon Atoms for Alkanes
Assumption: Constant Pressure System
General Hydrocarbon: CxHy + [x+(y/4)](O2 + 3.76N2) <=> xCO2 + (y/2)H2O + 3.76[x+(y/4)]N2
Methane: CH4 + 2(O2 + 3.76N2) <=> CO2 + 2H2O + 7.52N2
Ethane: C2H6 + 3.5(O2 + 3.76N2) <=> 2CO2 + 3H2O + 13.16N2
Propane: C3H8 + 5(O2 + 3.76N2) <=> 3CO2 + 4H2O + 18.8N2
Program by: Shouvik Bandopadhyay
\"\"\"

R = 8.314 # J/mol-K

# Function to calculate enthalpy of a compound through NASA polynomial function
def h(T, coeff):
	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

# Coefficients for the NASA polynomial for low temperatures
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]
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]

# Coefficients for the NASA polynomial for high temperatures
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]

# Standard temperature condition
Tstp = 298.15 # K

# Defining the function for root-solving by Newton-Raphson method
def f(T, n):
	h_ch4_r = h(Tstp, ch4_coeff_l)
	h_c2h6_r = h(Tstp, c2h6_coeff_l)
	h_c3h8_r = h(Tstp, c3h8_coeff_l)
	h_o2_r = h(Tstp, o2_coeff_l)
	h_n2_r = h(Tstp, n2_coeff_l)
	h_co2_p = h(T, co2_coeff_h)
	h_h2o_p = h(T, h2o_coeff_h)
	h_n2_p = h(T, n2_coeff_h)
	if n==1:
		x = 1
		y = 4
		h_fuel = h_ch4_r
	elif n==2:
		x = 2
		y = 6
		h_fuel = h_c2h6_r
	elif n==3:
		x = 3
		y = 8
		h_fuel = h_c3h8_r
	H_reactants = h_fuel + (x+y/4)*h_o2_r + 3.76*(x+y/4)*h_n2_r
	H_products = x*h_co2_p + (y/2)*h_h2o_p + 3.76*(x+y/4)*h_n2_p
	return H_products - H_reactants

# Defining the derivative of the function
def fprime(T, n):
	return (f(T+1e-6, n) - f(T, n))/(1e-6)

# Newton-Raphson solver
alpha = 1
tol = 1e-3
n = [1, 2, 3]
T_aft_py = []
for i in range(0,len(n)):
	T_guess = 1500 # K
	while abs(f(T_guess, n[i])) >= tol:
		T_guess = T_guess - alpha*((f(T_guess, n[i]))/(fprime(T_guess, n[i])))
	T_aft_py.append(T_guess)
	print(T_aft_py)
print()

# Plotting the Python results
plt.bar([\'Methane ($CH_4$)\',\'Ethane ($C_2H_6$)\',\'Propane ($C_3H_8$)\'],T_aft_py)
plt.ylabel(\'Adiabatic Flame Temperature [K]\',fontsize=12)
plt.title(\'AFT vs. No. of C atoms - Alkanes (Python)\',fontsize=13,fontweight=\'bold\')
plt.show()

\"\"\"
Solution obtained by Cantera
\"\"\"
gas = ct.Solution(\'gri30.xml\')
T_aft_ct = []
for i in range(0,len(n)):

# Computing the mole fractions of each reactant
	if n[i]==1:
		x = 1
		y = 4
		fuel = \'CH4\'
	elif n[i]==2:
		x = 2
		y = 6
		fuel = \'C2H6\'
	elif n[i]==3:
		x = 3
		y = 8
		fuel = \'C3H8\'
	n_tot = 1 + (x+y/4) + 3.76*(x+y/4)
	n_fuel = 1/n_tot
	n_o2 = (x+y/4)/n_tot
	n_n2 = (3.76*(x+y/4))/n_tot

# Defining the gas mixture
gas.TPX = 298.15, 101325, {fuel:n_fuel, \'O2\':n_o2, \'N2\': n_n2}

# Calculating the equilibrium conditions
gas.equilibrate(\'HP\',\'auto\')
T_aft_ct.append(gas.T)
print(T_aft_ct)

# Plotting the Cantera results
plt.bar([\'Methane ($CH_4$)\',\'Ethane ($C_2H_6$)\',\'Propane ($C_3H_8$)\'],T_aft_py)
plt.ylabel(\'Adiabatic Flame Temperature [K]\',fontsize=12)
plt.title(\'AFT vs. No. of C atoms - Alkanes (Cantera)\',fontsize=13,fontweight=\'bold\')
plt.show()

Output from the above code

HYDROCARBON	PYTHON AFT IN KELVIN	CANTERA AFT IN KELVIN
CH4	2325.8	2224.65
C2H6	2379.84	2258.76
C3H8	2392.11	2268.7

• We can see that the AFT does not vary much `(∆ <100 K)` with a consequent increase in the number of C-atom chains in a hydrocarbon.

• Also, in this case, the difference between the solution obtained from Python and Cantera do not differ by a huge amount.

• The reaction mechanism for alkanes is fairly simple and it does not produce any unanticipated products as opposed to the previous case with the combustion of alkenes and alkynes.