Adiabatic Flame Temperature calculation using Python and Cantera
Objective:
PART_1- To find out the effect of equivalence ratio on AFT(Adiabatic Flame temperature) assuming METHANE(CH4) in a constant volume chamber.
PART-2- To find out the effect of equivalence ratio on AFT by considering Constant pressure reactor with heat loss.
- Test it for -ane, -ene, and -yne variants of 2 Carbon atom chain and plot a trend.
- To plot the effect of the number of carbon atoms (from 1 to 3) on the final temperature.
Theory:
Adiabatic Flame Temperature
This temperature is defined as the temperature that is achieved by allowing fuel to react with an oxidizer and approach equilibrium by means of an isobaric and adiabatic process. The adiabatic flame temperature is the highest temperature a system can achieve at constant Pressure and also at constant volume. This temperature is achieved by a combustion system at constant pressure when there are no heat losses from the system. In other words, all of the heat from the combustion reaction is used to heat up the products of combustion.
A benchmark temperature for any combustion system is the “Adiabatic Flame Temperature”. It gives the essence of the maximum possible temperature rises by burning the fuel.
Here, we can define two adiabatic flame temperatures: one for constant pressure combustion and one for a constant volume.
In Constant-Pressure combustion, if a fuel-air mixture burns adiabatically at constant pressure, the absolute enthalpy of the reactant at the initial state will be equal to the enthalpy of the products at the final state.
`H_("reactant")(T_i,P)=H_("product")(T_(ad),P)`_____equation 1
Above equation can be written as
`h_("reactant")(T_i,P)=h_("product")(T_(ad),P)`______per mass of mixture basis
This equation will be called it a constant-pressure adiabatic flame temperature. This equation helps us to deal with a gas-turbine combustor or a furnace.
In Constant-Volume combustion, this we will basically be required for analysis of otto-cycle analysis
`U_("reactant")(T_("initial"),P_("initial"))=U_("product")(T_(ad),P_f)`
we can rearrange the equation n following form,
`H("reactant")-H("product")-V(P_("init")-P_f)=0`_____equation 2
According to ideal gas law, PV= nRT
`P_("initial")*V= N_("reactants")*R_U*T_("initial")`
`P_("f")*V= N_("products")*R_U*T_("ad")`
Thus, equation 2 will be written as
`H_("react")-H_("prod")-R_U*(N_("rect")*T_("init")-N_("prod")*T_("ad"))=0`
Newton-Raphson Method
The Newton-Raphson method is one of the most widely used methods for root finding. It can be easily generalized to the problem of finding solutions of a system of non-linear equations, which is referred to as Newton's technique. Moreover, it can be shown that the technique is quadratically convergent as we approach the root.
Unlike the bisection and false position methods, the Newton-Raphson methods require only one initial value `x_0`, which we will refer to as the initial guess for the root. To see how the N-R method works, we can rewrite the function f(x) using a Taylor series expansion in (x-x0):
`f(x)= f(x_0)+f(x_0)(x-x_0)+(1/2)f'(x_0)(x-(x_0))^2+....=0`
where, `f(x)`denotes the first derivative of `f(x)` with respect to x, `f'(x)`is the second derivative, and so forth. Now, suppose the initial guess is pretty close to the real root. Then (x-x0) is small, and only the first few terms in the series are important to get an accurate estimate of the true root, given x0. By truncating the series at the second term (linear in x), we obtain the N-R iteration formula for getting a better estimate of the true root:
`x_1=x_0-(f(x_0))/(f'(x_0))`
Thus the N-R method finds the tangent to the function f(x) at x=x0 and extrapolates it to intersect the x-axis to get x1. This point of intersection is taken as the new approximation to the root and the procedure is repeated until convergence is obtained whenever possible. Mathematically, given the value of x = xi at the end of the ith iteration, we obtain xi+1 as
`x_(i+1)=x_i-(f(x_i))/(f'(x_i))`
Graphically we can represent this method is as:
Equivalence ratio:
The premixing of the fuel and air prior to the flame zone creates a mixture strength that is characterized by an equivalence ratio designated by the Greek letter phi (φ). The equivalence ratio is defined as the ratio of the stoichiometric air/fuel ratio to the actual air/fuel ratio. That is
`phi=("Fuel"/("Air"))_"Actual condition"/("Fuel"/("Air"))_"Stoichiometric condition"`
as we can say,
φ = 1 Stoichiometric mixture
φ < 1 Fuel lean mixture (excess air)
φ > 1 Fuel rich mixture (excess fuel)
SOLUTIONS:
PART-1: Effect of equivalence ratio on AFT of fuel METHANE(CH4) assuming constant volume chamber:
The general form of a stoichiometric equation
`CxHy + (x+y/4)*(O2+3.76*N2) -> x*Co2 + (y/2)*H2O + (x+y/4)*N2`
For Lean mixture
`CH4 + (2/e_r)*(O2+3.76*N2) -> CO2 + 2H2O + ((2/e_r)-2)O2 + (7.52/e_r)*N2`
For Rich Mixture
`CH4 + (2/e_r)*(O2+3.76*N2) -> ((4/e_r)-3)CO2 + 2H2O + (4-(4/e_r))CO + (7.52/e_r)*N2`
where, `e_r` represents equivalnece ratio.
Python code:
NumPy- To perform numerical analysis
Matplotlib.pyplot- To create charts for better visualization of the results.
Cantera- It is a standard chemical kinetics library that contains standard chemical mechanisms GRI mechanisms, thermodynamic properties of several species that are included in a reaction.
STEP-1: I have created the range of equivalence ratio that will be participating in the reaction using linspace in-build function ranging from 0.1 till 2.0 where 200 points will be taken into calculations.
STEP-2: Next we have defined the functions from which we will help to get enthalpy of reactants and products by using the equation from NASA polynomial form
`H/(R*T)=a1+a2*(T/2)+a3*(T^2/3)+a4*(T^3/4)+a5*(T^4/5)+(a6)/T`
STEP-3: Defined the low and high-temperature coefficients for participating species.
STEP-4: Function defined which will actually go for finding root iteratively, in this function based on the selection of a range of equivalence ratio enthalpy of the product will be calculated. As reactants are at initial state enthalpy of the reactant will be calculated by considering low-temperature polynomial coefficients.
STEP-5: In this step, for the second derivative numerical derivative of function defined using FORWARD-DIFFERENCING.
STEP-6: Next, we have set up the loop for performing the N-R method for root finding by considering all values of equivalence ratio from which we will get AFT of methane at all values of e_r.
"""
Programm to find out the effect of equivalence ratio on AFT of CH4 assuming constant volume chamber
General form of stoichiometric equation
CxHy + (x+y/4)*(O2+3.76*N2) -> x*Co2 + (y/2)*H2O + (x+y/4)*N2
For Lean mixture
CH4 + (2/e_r)*(O2+3.76*N2) -> CO2 + 2H2O + ((2/e_r)-2)O2 + (7.52/e_r)*N2
For Rich Mixture
CH4 + (2/e_r)*(O2+3.76*N2) -> ((4/e_r)-3)CO2 + 2H2O + (4-(4/e_r))CO + (7.52/e_r)*N2
By- GAURAV V. KHARWADE | Skill-lync-2020
"""
import cantera as ct
import math
import numpy
import matplotlib.pyplot as plt
# To define the range of equivalence ratio (0.1,0.2,0.3.....2.0)
i = numpy.linspace(0.1,2.0,200)
e_r = i # e_r represents range of equivalence ratio
def h(T,co_effs): # Function to calculate enthalpy of all reactants and products
R = 8.314 # J/mol-K
a1 = co_effs[0]
a2 = co_effs[1]
a3 = co_effs[2]
a4 = co_effs[3]
a5 = co_effs[4]
a6 = co_effs[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
# Low temperature coefficients (Reactants) at STP condition
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 (Products)
N2_coeffs_h = [0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04,0.05980528E+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 ]
CO_coeffs_h = [2.71518561E+00,2.06252743E-03,-9.98825771E-07,2.30053008E-10,-2.03647716E-14,-1.41518724E+04,7.81868772E+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]
def f(T,e_r):
# To calculate enathlpy of reactants and products and to perform root finding problem
# For calculation of enthalpy of Products in Jouls
R = 8.314 #J/mol-K
h_N2_p = h(T,N2_coeffs_h) # J/Mol-K
h_CO2_p = h(T,CO2_coeffs_h) # J/Mol-K
h_H2O_p = h(T,H2O_coeffs_h) # J/Mol-K
h_CO_p = h(T,CO_coeffs_h) # J/Mol-K
h_O2_p = h(T,O2_coeffs_h) # J/Mol-K
if e_r > 1: # Conditional loop based on equivalence ratio enthalpy of products will be find out
H_products = ((4/e_r)-3)*h_CO2_p + (4-(4/e_r))*h_CO_p + 2*h_H2O_p + (7.52/e_r)*h_N2_p # Enthalpy in J basis
N_products = (((4/e_r)-3) +2 + (4-4/e_r) + (7.52/e_r)) * R * T
else:
H_products = h_CO2_p + ((2/e_r)-2)*h_O2_p + 2*h_H2O_p + (7.52/e_r)*h_N2_p # Enthalpy in J basis
N_products = (((2/e_r)-2) + 3 + (7.52/e_r)) * R * T
print(H_products)
# For calculation of enthalpy of Reactants in Jouls
Tstd = 298.15 # in K at STP condition
h_CH4_r = h(Tstd,CH4_coeffs_l) # J/Mol-K
h_O2_r = h(Tstd,O2_coeffs_l) # J/Mol-K
h_N2_r = h(Tstd,N2_coeffs_l) # J/Mol-K
H_reactants = h_CH4_r + (2/e_r)*h_O2_r + ((2/e_r)*3.76)*h_N2_r # Enthalpy in J basis
N_reactants = (1 + (2/e_r) + ((2/e_r)*3.76)) * R * Tstd
print(H_reactants)
return (H_reactants - H_products) - (N_reactants - N_products)
def fprime(T,e_r):
"""
To define numerical derivative of function (forward-differencing)
"""
return (f(T+1e-6,e_r)-f(T,e_r))/1e-6
# Using Newton-Raphson to solve iteratively
tol = 1e-3
alpha = 0.9
Temp= [] # Defined arrays to store the values of temperature
for j in e_r: # Loop to perform a root-finding procedure for all values of equivalence ratio between 0.1 to 2
T_guess = 1500
while(abs(f(T_guess,j)) > tol):
T_guess = T_guess - alpha*((f(T_guess,j)/fprime(T_guess,j))) # N-R iteration formula for getting better estimate of true root
Temp.append(T_guess)
print(Temp)
plt.plot(e_r,Temp,color='red')
plt.xlabel('Equivalence ratio (e_r)')
plt.ylabel('Adiabatic Flame temperature (K)')
plt.title('Equivalence Ratio VS Adiabatic Flame temperature for CH4')
plt.grid('on')
#####################################################################################################
# AFT calculation using CANTERA
gas = ct.Solution('gri30.xml')
n= 200
phi= numpy.linspace(0.1,2.0,n)
tad = numpy.zeros(n)
for i in range(n):
gas.set_equivalence_ratio(phi[i], 'CH4', 'O2:1.0, N2:3.76')
gas.TP = 298.15,101325
gas.equilibrate('UV','auto') # UV where we have fixed the volume
tad[i] = gas.T
print('At phi = {0:12.4g}, tad = {1:12.4g}'.format(phi[i], tad[i]))
print(phi,tad)
plt.plot(phi,tad)
plt.legend(('Python','Cantera'))
plt.show()Results:
- We can observe here, we will get the high adiabatic flame temperature at e_r=1 using python but if we take a look at the curve of AFT using Cantera high AFT is in between the equivalence ratio of 1 to 1.25.
- Nos. of participating species at the product side using python is less compare to Cantera hence we get a somewhat high adiabatic flame temperature than that of CANTERA. In CANTERA, since it is a standard chemical kinetics library the nos. of species will be more at the product side.
- Practically we operate at an equivalence ratio of 0.5 to 1.5.
- Cantera gives slightly more value of AFT after the equivalence ratio increases more than 1.3 approximately.
PART-2: Constant pressure reactor with heat loss.
2.1- To find out the effect of the heat loss coefficient on AFT of methane(CH4) for constant pressure reactor.
In this part, we are going to look at the effect of heat loss from the constant pressure reactor on AFT where stoichiometric combustion of methane (CH4) fuel will occur in presence of O2 and N2.
`CH_4 + 2*(O_2+3.76*N_2) -> CO_2 + 2*H_2O + 7.52*N_2`
Here, we have assumed the equivalence ratio as 1.
Fuel utilization efficiency will be written as,
`eta=dotQ/(dot(m_f)*LHV)`
for simplicity of calculation of maximum possible heat loss, assume efficiency is 100% or 1.
The equation will be written as
`dotQ=(dot(m_f)*LHV)`
for methane (CH4),
mf = 16.04 g/mol
LHV = 50 MJ/Kg
Python code:
I have followed the same steps as mentioned in the earlier program. Instead of a range of equivalence ratios, I have defined the range of the Heat loss coefficient from 0.1 to 1.
"""
Programm to find out the effect of heat loss on AFT of methane(CH4) for constant pressure reactor
The general form of the stoichiometric equation
CH4 + 2*(O2+3.76*N2) -> Co2 + 2*H2O + 7.52*N2
By- GAURAV V. KHARWADE | Skill-Lync-2020
"""
import math
import numpy
import matplotlib.pyplot as plt
"""
Maximum possible heat loss will be calculated as
Q = mass of fuel * LHV
Q = mf * LHV where mf - (g/mol)
mf = 16.04 g/mol
LHV = 50 MJ/kg
HHV = 55 MJ/kg
Q = 16.04e-3 * 50e+6
Q = 80200
"""
# To define the range of heat loss coefficient (0.1,0.2,0.3.....1.0)
H_loss = numpy.linspace(0.1,1.0,200)
Qmax = 80200 # Maximum possible heat loss
e_r = 1 # Equivalence ratio is assumed to be 1
def h(T,co_effs): # Function to calculate the enthalpy of all reactants and products
R = 8.314 # J/mol-K
a1 = co_effs[0]
a2 = co_effs[1]
a3 = co_effs[2]
a4 = co_effs[3]
a5 = co_effs[4]
a6 = co_effs[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
# Low-temperature coefficients (Reactants) at STP condition
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 (Products)
N2_coeffs_h = [0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04,0.05980528E+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 ]
CO_coeffs_h = [2.71518561E+00,2.06252743E-03,-9.98825771E-07,2.30053008E-10,-2.03647716E-14,-1.41518724E+04,7.81868772E+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]
def f(T,H_loss):
#To calculate the enthalpy of reactants and products and to perform root-finding problem
# For calculation of enthalpy of Products in Jouls
R = 8.314 #J/mol-K
h_N2_p = h(T,N2_coeffs_h) # J/Mol-K
h_CO2_p = h(T,CO2_coeffs_h) # J/Mol-K
h_H2O_p = h(T,H2O_coeffs_h) # J/Mol-K
h_CO_p = h(T,CO_coeffs_h) # J/Mol-K
h_O2_p = h(T,O2_coeffs_h) # J/Mol-K
if e_r > 1:
H_products = ((4/e_r)-3)*h_CO2_p + (4-(4/e_r))*h_CO_p + 2*h_H2O_p + (7.52/e_r)*h_N2_p # Enthalpy in J basis
N_products = (((4/e_r)-3) +2 + (4-4/e_r) + (7.52/e_r)) * R * T
else:
H_products = h_CO2_p + ((2/e_r)-2)*h_O2_p + 2*h_H2O_p + (7.52/e_r)*h_N2_p # Enthalpy in J basis
N_products = (((2/e_r)-2) + 3 + (7.52/e_r)) * R * T
print(H_products)
# For calculation of enthalpy of Reactants in Jouls
Tstd = 298.15 # in K at STP condition
h_CH4_r = h(Tstd,CH4_coeffs_l) # J/Mol-K
h_O2_r = h(Tstd,O2_coeffs_l) # J/Mol-K
h_N2_r = h(Tstd,N2_coeffs_l) # J/Mol-K
H_reactants = h_CH4_r + (2/e_r)*h_O2_r + ((2/e_r)*3.76)*h_N2_r # Enthalpy in J basis
N_reactants = (1 + (2/e_r) + ((2/e_r)*3.76)) * R * Tstd
print(H_reactants)
return (H_products - H_reactants) + (H_loss*Qmax)
def fprime(T,H_loss):
"""
To define numerical derivative of function (forward-differencing)
"""
return (f(T+1e-6,H_loss)-f(T,H_loss))/1e-6
# Using Newton-Raphson to solve iteratively
tol = 1e-3
alpha = 0.9
Temp= [] # To define array of temperature
for j in H_loss:
T_guess = 1500
while(abs(f(T_guess,j)) > tol):
T_guess = T_guess - alpha*((f(T_guess,j)/fprime(T_guess,j)))
Temp.append(T_guess)
print(Temp)
# To Find the AFT of CH4 with H_loss = 0.35
h_l = numpy.searchsorted(H_loss, 0.35) # To find the index of element H_loss = 0.35
print('Index value is:',h_l)
t = Temp[h_l] # To find out the value of temperature from array Temp at index no 56.
print('The AFT of CH4 with H_loss @ 0.35 is:',t)
plt.plot(H_loss,Temp,color='red')
plt.xlabel('Heat Loss coefficients (H_loss)')
plt.ylabel('Adiabatic Flame temperature (K)')
plt.title('Heat Loss coefficients VS Adiabatic Flame temperature for CH4')
plt.grid('on')
plt.show()Result:
- From the graph, we can say increasing in the loss coefficient will increase in heat rejection rate from the constant pressure reactor which leads to a decrease in adiabatic flame temperature.
- AFT variation is linear with the Heat loss coefficient as we can see in the above plot.
- AFT @ heat loss coefficient= 0.35 is found out to be 2261.5395K
2.2: Test it only for -ane, -ene and -yne variants of 2 Carbon atom chain and plot a trend.
Alkanes, Alkenes, & Alkynes are all open-chain structures, while cyclanes and aromatics exhibit ring structures.
Since we have fixed carbon atoms to 2, different fuel species will participate in the combustion process.
Python code:
"""
Programme to calculate AFT of alkane(CnH2n+2), alkene(CnH2n) & alkyne(CnH2n-2) using python programming.
Assuming carbon atoms are 2
The stoichiometric reaction of alkane(C2H6-Ethane), alkene(C2H4-Ethylene) & alkyne(C2H2-Acetylene)
General representation of stoichiometric reaction written as
CxHy + (x+y/4)*(O2+3.76*N2) -> x*Co2 + (y/2)*H2O + (x+y/4)*3.76*N2
(Reactants) -> (Products)
By- Gaurav V. Kharwade || Skill-Lync-2020
"""
import cantera as ct
import numpy as np
import matplotlib.pyplot as plt
import math
def h(T,co_effs): # Function to calculate the enthalpy of all reactants and products
R = 8.314 # J/mol-K
a1 = co_effs[0]
a2 = co_effs[1]
a3 = co_effs[2]
a4 = co_effs[3]
a5 = co_effs[4]
a6 = co_effs[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
# Low-temperature coefficients of Reactants at STP condition
C2H6_coeffs_l = [4.29142492E+00,-5.50154270E-03,5.99438288E-05,-7.08466285E-08,2.68685771E-11,-1.15222055E+04,2.66682316E+00]
C2H4_coeffs_l = [3.95920148E+00,-7.57052247E-03,5.70990292E-05,-6.91588753E-08,2.69884373E-11,5.08977593E+03,4.09733096E+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]
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 of Products
N2_coeffs_h = [0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04,0.05980528E+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 ]
Tstd = 298.15 # in K at STP condition
h_C2H6_r = h(Tstd,C2H6_coeffs_l)
h_C2H4_r = h(Tstd,C2H4_coeffs_l)
h_C2H2_r = h(Tstd,C2H2_coeffs_l)
H_alkane = h_C2H6_r
H_alkene = h_C2H4_r
H_alkyne = h_C2H2_r
H = [H_alkane,H_alkene,H_alkyne]
def f(T,H):
"""
To calculate the enthalpy of reactants and products and to perform root-finding problem
"""
# For calculation of enthalpy of Products in Jouls
h_N2_p = h(T,N2_coeffs_h) # J/Mol-K
h_CO2_p = h(T,CO2_coeffs_h) # J/Mol-K
h_H2O_p = h(T,H2O_coeffs_h) # J/Mol-K
# For calculation of enthalpy of Reactants in Jouls
h_O2_r = h(Tstd,O2_coeffs_l) # J/Mol-K
h_N2_r = h(Tstd,N2_coeffs_l) # J/Mol-K
if H == H_alkane:
x = 2
y = 6
a = x + y/4
b = y/2
H_reactants = H_alkane + a*h_O2_r + 3.76*a*h_N2_r # Enthalpy in J basis
print(H_reactants)
H_products = x*h_CO2_p + b*h_H2O_p + 3.76*a*h_N2_p # Enthalpy in J basis
print(H_products)
elif H == H_alkene:
x = 2
y = 4
a = x + y/4
b = y/2
H_reactants = H_alkene + a*h_O2_r + 3.76*a*h_N2_r # Enthalpy in J basis
print(H_reactants)
H_products = x*h_CO2_p + b*h_H2O_p + 3.76*a*h_N2_p # Enthalpy in J basis
print(H_products)
elif H == H_alkyne:
x = 2
y = 2
a = x + y/4
b = y/2
H_reactants = H_alkyne + a*h_O2_r + 3.76*a*h_N2_r # Enthalpy in J basis
print(H_reactants)
H_products = x*h_CO2_p + b*h_H2O_p + 3.76*a*h_N2_p # Enthalpy in J basis
print(H_products)
return (H_products - H_reactants)
def fprime(T,H):
"""
To define numerical derivative of function (forward-differencing)
"""
return (f(T+1e-6,H)-f(T,H))/1e-6
T_guess = 1900
tol = 1e-3
alpha = 0.1
iter = 0
Temp= []
for i in H:
while(abs(f(T_guess,i)) > tol):
T_guess = T_guess - alpha*((f(T_guess,i)/fprime(T_guess,i)))
Temp.append(T_guess)
print(Temp)
print('Nos. of iterations=',iter)
print('AFT of CH4 is =',T_guess)
plt.ylabel('AFT in K')
plt.title('Effect of -ane,-ene & -yne on AFT using PYTHON')
plt.bar(['Alkane','Alkene','Alkyne'],Temp,color='red')
plt.show()
#####################################################
# AFT calculation using CANTERA
gas = ct.Solution('gri30.xml')
gas.TPX = 298.15,101325,{'C2H6':1,'O2':3.5,'N2':3.5*3.76}
gas.equilibrate('HP','auto')
T_alkane = gas.T
gas.TPX = 298.15,101325,{'C2H4':1,'O2':3,'N2':3*3.76}
gas.equilibrate('HP','auto')
T_alkene = gas.T
gas.TPX = 298.15,101325,{'C2H2':1,'O2':2.5,'N2':2.5*3.76}
gas.equilibrate('HP','auto')
T_alkyne = gas.T
T = (T_alkane,T_alkene,T_alkyne)
T= np.array(T)
print(T)
plt.bar(['Alkane','Alkene','Alkyne'],T,color='red')
plt.title('Effect of -ane,-ene & -yne on AFT using CANTERA')
plt.show()
Results:
- We can observe here is that lowering the hydrogen atoms leads to increasing the AFT of CH4 because it will produce less H2O on the product side.
- In Cantera, for all fuel types, AFT is lower because of an increase in species at the product's side.
2.3: For a general alkane, plot the effect of the number of carbon atoms (from 1 to 3) on the final temperature.
Assume Heat loss is 0.
Python code:
"""
Program to calculate AFT of alkane(CnH2n+2) by varying nos. of carbon atoms from 1 to 3 using python programming.
General representation of stoichiometric reaction written as
CxHy + (x+y/4)*(O2+3.76*N2) -> x*Co2 + (y/2)*H2O + (x+y/4)*3.76*N2
(Reactants) -> (Products)
By- Gaurav V. Kharwade || Skill-Lync-2020
"""
import cantera as ct
import numpy as np
import matplotlib.pyplot as plt
import math
def h(T,co_effs): # Function to calculate the enthalpy of all reactants and products
R = 8.314 # J/mol-K
a1 = co_effs[0]
a2 = co_effs[1]
a3 = co_effs[2]
a4 = co_effs[3]
a5 = co_effs[4]
a6 = co_effs[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
# Low-temperature coefficients of Reactants at STP condition
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 ]
# High-temperature coefficients of Products
N2_coeffs_h = [0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04,0.05980528E+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 ]
Tstd = 298.15 # in K at STP condition
h_CH4_r = h(Tstd,CH4_coeffs_l)
h_C2H6_r = h(Tstd,C2H6_coeffs_l)
h_C3H8_r = h(Tstd,C3H8_coeffs_l)
C_1 = h_CH4_r
C_2 = h_C2H6_r
C_3 = h_C3H8_r
C = [C_1,C_2,C_3]
def f(T,C):
"""
To calculate the enthalpy of reactants and products and to perform root-finding problem
"""
# For calculation of enthalpy of Products in Jouls
h_N2_p = h(T,N2_coeffs_h) # J/Mol-K
h_CO2_p = h(T,CO2_coeffs_h) # J/Mol-K
h_H2O_p = h(T,H2O_coeffs_h) # J/Mol-K
# For calculation of enthalpy of Reactants in Jouls
h_O2_r = h(Tstd,O2_coeffs_l) # J/Mol-K
h_N2_r = h(Tstd,N2_coeffs_l) # J/Mol-K
if C == C_1:
x = 1
y = 2*(x+1)
a = x + y/4
b = y/2
H_reactants = C_1 + a*h_O2_r + 3.76*a*h_N2_r # Enthalpy in J basis
print(H_reactants)
H_products = x*h_CO2_p + b*h_H2O_p + 3.76*a*h_N2_p # Enthalpy in J basis
print(H_products)
elif C == C_2:
x = 2
y = 2*(x+1)
a = x + y/4
b = y/2
H_reactants = C_2 + a*h_O2_r + 3.76*a*h_N2_r # Enthalpy in J basis
print(H_reactants)
H_products = x*h_CO2_p + b*h_H2O_p + 3.76*a*h_N2_p # Enthalpy in J basis
print(H_products)
elif C == C_3:
x = 3
y = 2*(x+1)
a = x + y/4
b = y/2
H_reactants = C_3 + a*h_O2_r + 3.76*a*h_N2_r # Enthalpy in J basis
print(H_reactants)
H_products = x*h_CO2_p + b*h_H2O_p + 3.76*a*h_N2_p # Enthalpy in J basis
print(H_products)
return (H_products - H_reactants)
def fprime(T,C):
"""
To define numerical derivative of function (forward-differencing)
"""
return (f(T+1e-6,C)-f(T,C))/1e-6
T_guess = 1900
tol = 1e-3
alpha = 0.1
Temp= []
for i in C:
while(abs(f(T_guess,i)) > tol):
T_guess = T_guess - alpha*((f(T_guess,i)/fprime(T_guess,i)))
Temp.append(T_guess)
print(Temp)
plt.ylabel('AFT in K')
plt.title('Effect of carbon atomes of ALKANE on AFT using PYTHON')
plt.bar(['1-carbon(CH4)','2-carbon(C2H6)','3-carbon(C3H8)'],Temp,color='red')
plt.show()
Results:
- More nos. of carbon atoms means it will produce more heat, hence fuel with 3 carbon atoms gives us slightly higher AFT than that of fuel with lower carbon atoms of the alkane family
In this way, I have shown through this study that what is AFT, how this AFT will be impacted with changing equivalence ratio, type of fuel, nos. of carbon atoms of alkane, and also I have demonstrated the effect of heat loss coefficient on AFT of constant pressure chamber where combustion is being done by using METHANE.