Week 3 - Adiabatic Flame Temperature calculation
AIM
* To write a program in python to plot the effect of equivalence ratio on Adiabatic flame temperature (AFT) of methane contained in a constant volume chamber assuming N2 is present.
* To find out AFT of methane contained in a constant pressure chamber with `phi` = 1, `lamda` = 1 and H_loss = 0.35.
* To find out AFT of 2 carbon variants of Alkane, Alkene and Alkyne in a constant pressure chamber.
* To plot the effect of number of carbon atoms of a general alkane (1 to 3) on AFT with H_loss = 0.
INTRODUCTION
Before starting with the codes for the above problem statements it will be good to provide brief explaination about what is a combustion reaction, adiabatic flame temperature equivalence ratio.
Combustion reaction - Any reaction in which a substance reacts with oxygen in the air and hence produces heat energy. This is a simple definition of a combustion reaction, when a fuel burns it combines with oxygen, because of this reaction the heat energy is produced. In this problem solving codes methane (CH4) is used as fuel. Along with oxygen nitrogen is also considered as reactant. A combustion results in formation of CO2 as well as H2O as products along with many other minor products. All the chemical kinetics solver consider all the minor products as oxides of sulphur (SO_x), oxides of nitrogen (NO_x), carbon monoxide (CO) while computing combustion calculations but for this calculations under this challenge only CO2, H2O and N2 are considered on the product side.
Chemical equation of a combustion reaction with CH4 as fuel is given as `CH_(4) + ar(O_2 + 3.76N_2) = aCO_2 + bH_2O + ar*3.76N_2`, coefficients ar, a, b is determined by balancing the equation. Balancing the above chemical equation gives ar = 2, a = 1, b = 2 and hence equation can be written as `CH_(4) + 2(O_2 + 3.76N_2) = CO_2 + 2H_2O + 7.52N_2` this equation is also called as stoichometric equation.
Equivalence ratio (`phi`) - It is defined as fuel to air ratio of given mixture to that of a stoichometric mixture. Equivalence ratio can be better explained by considering an example. Consider a mixture `CH_(4) + 1.8(O_2 + 3.76N_2)`in which the reactant side is shown, from the stoichometric equation of the mixture on reactant side we have `CH_(4) + 2(O_2 + 3.76N_2) =`
so fuel to air ratio, `phi = (F/A)_(mix)/(F/A)_(st)`
hence we have `phi = (1/1.8)/(1/2) = 2/1.8 = 1.11`,
similarly we may have mixture where `phi` is less than 1, so mixtures with `phi < 1` is called lean mixture and `phi > 1` is called rich mixture.
Balanced chemical equation for different value of `phi` is mentioned below,
for rich mixture (`phi > 1`) - `CH_4 + (2/phi)*(O_2 + 3.76N_2) = (4/phi-3)CO_2 + 2H_2O + (7.52/phi)N_2 + (4-4/phi)CO`
for lean mixture (`phi < 1`) - `CH_4 + (2/phi)*(O_2 + 3.76N_2) = CO_2 + 2H_2O + (7.52/phi)N_2 + (2/phi - 2)O_2`
Adiabatic flame temperature (AFT) - when a combustion reaction takes place heat is released to the combustion product if there is no heat loss the max temperature attained by the combustion product is called Adiabatic flame temperature (AFT). In this challenge the AFT is computed using newton raphson method.
1. To write codes in python to plot the effect of equivalence ratio on AFT, fuel used CH4 in a constant volume chamber:
Logic:
* AFT can be determined using newton raphson method.
* For a constant volume combustion reaction difference of entalphy of product to that of reactant is zero. i.e, `H_(prouct) - H_(reactant) = 0`. Entalphy (H) can be calculated using NASA polynomials. A snip of the NASA polynomials is shown below`
In which a1, a2,....a7 are the coefficients, T is the temperature and R is universal gas constant. If there is no H_loss then the temperature of the product will be AFT. so this AFT is determined using newton raphson method.
* Coefficients changes from species to species and that too for a given species the coefficients varies for a range of temperature. This can be better understood by looking at thermo data file of GRI 3.0 mechanism file. A snip of the coefficients of `O_2` is shown below
First 7 coefficients can be used for entalphy calculation when the temperature T > 1000 K and less than 3500 K next 7 coefficients can be used when the temperature T > 200 K and less than 1000 K. In the highlighted temperature range 1000 K is the mid range. For each species in the themo30.dat file a specific temperature range is mentioned. Coefficients of reactants at 298.15 K and that of the product (at higher temp. range) is copied from thermo30 dat folder.
* A function is h(T,coeff) is defined to calculate the entalphy.
* Criteria used for a constant volume vessel is `H_p - H_r + 8.314*(n_r*T_(std) - n_p*T)`, where `H_p -`entalphy of product, `H_r -`entalphy of reactant
* T_aft = T_aft - (alpha*f(T_aft)/fprime(T_aft)), this formulation obtained using newton raphson method is applied to calculate AFT. T_aft value keeps updating from first iteration till the last iteration, so a function f(T) is defined such that the entalphy value of products is calculated for every iteration as T_aft value changes also it returns the criteria `H_(prouct) - H_(reactant) = 0`
* Since aim is to plot the variation of AFT with `phi`, to store value of AFT for each of `phi` value 2 matrices are used viz AFT_value and phi_value.
* A if loop is defined inside the function f(T) such that for lean mixture (`phi < 1`) `CH_4 + (2/phi)*(O_2 + 3.76N_2) = CO_2 + 2H_2O + (7.52/phi)N_2 + (2/phi - 2)O_2` is used also for rich mixture (`phi > 1`) `CH_4 + (2/phi)*(O_2 + 3.76N_2) = (4/phi-3)CO_2 + 2H_2O + (7.52/phi)N_2 + (4-4/phi)CO`is used.
Below is the python code for the same,
"""
Calculating variation of AFT with equivalence ratio
"""
import matplotlib.pyplot as plt
"""
Defining function to calculate entalphy
"""
def h(T,coeffs):
R = 8.314 # J/mol*K
a1 = coeffs[0]
a2 = coeffs[1]
a3 = coeffs[2]
a4 = coeffs[3]
a5 = coeffs[4]
a6 = coeffs[5]
return (a1 + a2/2*T + a3/3*pow(T,2) + a4/4*pow(T,3) + a5/5*pow(T,4) + a6/T)*R*T
"""
Defining coefficients of reactants
"""
coeff_CH4_r = [5.14987613E+00,-1.36709788E-02,4.91800599E-05,-4.84743026E-08,1.66693956E-11,-1.02466476E+04]
coeff_O2_r = [3.78245636E+00,-2.99673416E-03,9.84730201E-06,-9.68129509E-09,3.24372837E-12,-1.06394356E+03]
coeff_N2_r = [0.03298677E+02,0.14082404E-02,-0.03963222E-04,0.05641515E-07,-0.02444854E-10,-0.10208999E+04]
"""
Defining coefficients of product
"""
coeff_CO2_p = [3.85746029E+00,4.41437026E-03,-2.21481404E-06,5.23490188E-10,-4.72084164E-14,-4.87591660E+04]
coeff_H2O_p = [3.03399249E+00,2.17691804E-03,-1.64072518E-07,-9.70419870E-11,1.68200992E-14,-3.00042971E+04]
coeff_N2_p = [0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04]
coeff_O2_P = [3.28253784E+00,1.48308754E-03,-7.57966669E-07,2.09470555E-10,-2.16717794E-14,-1.08845772E+03]
coeff_CO_p = [2.71518561E+00,2.06252743E-03,-9.98825771E-07,2.30053008E-10,-2.03647716E-14,-1.41518724E+04]
"""
Defining matrix for phi
"""
phi_value = []
"""
Defining matrix for AFT
"""
AFT_value = []
""""
Defining a function to calculate entalphy of products
"""
def f(T):
h_CO2_p = h(T,coeff_CO2_p)
h_H2O_p = h(T,coeff_H2O_p)
h_N2_p = h(T,coeff_N2_p)
h_O2_p = h(T,coeff_O2_P)
h_CO_p = h(T,coeff_CO_p)
if phi <= 1:
h_prod = h_CO2_p + (2*h_H2O_p) + (7.52/phi*h_N2_p) + (((2/phi)-2)*h_O2_p)
n_p = 1 + 2 + (7.52/phi) + ((2/phi)-2)
else:
h_prod = ((4/phi)-3)*h_CO2_p + (2*h_H2O_p) + (7.52/phi*h_N2_p) + ((4-(4/phi))*h_CO_p)
n_p = ((4/phi)-3) + 2 + (7.52/phi) + (4-(4/phi))
return (h_prod - h_react + 8.314*((n_r*T_std) - (n_p*T)))
"""
Defining fprime
"""
def fprime(T):
return (f(T+1e-6) - f(T))/(1e-6)
for i in range (1,18):
phi = i/10
phi_value.append(phi)
"""
Calculating entalphy of reactants
"""
T_std = 298.15 # Known value of temperature
h_CH4_r = h(T_std,coeff_CH4_r)
h_O2_r = h(T_std,coeff_O2_r)
h_N2_r = h(T_std,coeff_N2_r)
# CH4 + (2/phi)*(O2 + 3.76N2) = CO2 + 2H2O + (7.52/phi)N2 + (2/phi - 2)O2 - Stoichometric equation for lean mixture
# CH4 + (2/phi)*(O2 + 3.76N2) = (4/phi-3)CO2 + 2H2O + (7.52/phi)N2 + (4-4/phi)CO - Stoichometric equation for rich mixture
h_react = (h_CH4_r*1) + (2/phi*h_O2_r) + (7.52/phi*h_N2_r)
n_r = 1 + (2/phi) + (7.52/phi)
tol = 1e-4 # Tolerance value
T_aft = 3000 # Initial guess value of temperature
alpha = 0.2 # Relaxation factor
while(abs(f(T_aft)) > tol):
T_aft = T_aft - (alpha*f(T_aft)/fprime(T_aft))
AFT_value.append(T_aft)
print(AFT_value)
print(phi_value)
plt.plot(phi_value,AFT_value)
plt.xlabel('Phi values')
plt.ylabel('AFT values')
plt.title('AFT v/s phi plot')
plt.show()
executing the above code, prints the AFT value and phi value a snip of the same is provided below,
Variation of `phi` with AFT is plotted, which is as shown below,
So from the above plot the variation of AFT with `phi`is obtained, but this plot needs to be verified for which cantera module is used.
Before proceeding to the codes using cantera modules it will be better to understand a bit about cantera.
Cantera is a chemical kinetics library used in many chemical kinetics solver. It has a rich collection of information pretaining to different species. AFT can be calculated using cantera module in python, for this first cantera library needs to be imported. 'import cantera as ct' is the command line used for the same
Logic:
* Import cantera library as ct, using the command 'import cantera as ct'.
* define a class named as gas, which can be done using 'gas = ct.Solution('gri30.xml').
* Executing the command help(gas), prints all the information and different function that can be used to calculate different recation aspects related to different species. One such parameter that is available is TPX.
* Stoichometric combustion of CH4 with N2 can be given as `CH_(4) + ar(O_2 + 3.76N_2) = aCO_2 + bH_2O + ar*3.76N_2`AFT of this stoichometric combustion can be determined from cantera module using 2 functions called TPX and equilibriate.
gas.TPX = 298.15,101325{'CH4':1,'O2':2,'N2':7.52}
gas.equilibriate('UV','auto')
executing the above 2 command's gives us the AFT of CH4
* Since the reaction is considered as a constant volume reaction the 'UV' is mentioned inside .equilibriate, it represents that `U_2 =U_1`is used to derive the criteria and the reaction is constant volume.
* As `phi`value changes the molar coefficients of species on reactant side also changes so, from the stoichometric equation for lean and rich mixture the molar coefficients of CH4 is always 1, molar coefficient of o2 is (2/`phi`) and molar coefficient of N2 is `(2/phi)*3.76`
* So accorindly the gas.TPX is modified as gas.TPX = 298.15,101325{'CH4':1,'O2':(2/`phi`),'N2':(2/phi)*3.76}
* Values of AFT for different phi values is stored in 2 different matrices, AFT_value and phi_value
* For plotting using a python code it is necessary that we import appropriate library. In here matplotlib.pyplot is imported as plt
Codes for the same is shown below:
"""
Determining variation of AFT with equivalence ratio using cantera module
"""
import cantera as ct
import matplotlib.pyplot as plt
# CH4 + (2/phi)*(O2 + 3.76N2) = CO2 + 2H2O + (7.52/phi)N2 + (2/phi - 2)O2 - Stoichometric equation for lean mixture
# CH4 + (2/phi)*(O2 + 3.76N2) = (4/phi-3)CO2 + 2H2O + (7.52/phi)N2 + (4-4/phi)CO - Stoichometric equation for rich mixture
gas = ct.Solution('gri30.xml')
"""
Defining a matrix for phi
"""
phi_value = []
"""
Defining a matrix for AFT
"""
AFT_value = []
for i in range(1,18):
phi = i/10
phi_value.append(phi)
gas.TPX = 298.15,101325,{'CH4':1,'O2':2/phi,'N2':7.52/phi}
gas.equilibrate('UV','auto')
AFT_value.append(gas.T)
print(AFT_value)
print(phi_value)
plt.plot(phi_value,AFT_value)
plt.xlabel('Phi values')
plt.ylabel('AFT values')
plt.title('AFT v/s phi using cantera module')
plt.show()
Executing the above command prints the value of AFT obtained at different values of `phi`. A snip of the same is shown below,
plot obtained as follows,
Observation: From the plot obtained using code and cantera module matches and hence the variation of AFT with `phi`is verified.
2. To find AFT of methane conatined in a constant pressure chamber with H_loss = 0.35*H_max
Logic:
* Logic is same as that of the logic mentioned above, all the functions defined while determining the variation of AFT with equivalence ratio need's to be defined in here also.
* Only difference is in the criteria, for a constant pressure combustion the criteria used for newton raphson method is `H_(p) - H_(r) = 0`, since a `H_(loss) = 0.35*H_(max)`the criteria is `H_(p) -H_(r) + 0.35*H_(max)`
* `H_(max)` mentioned above is the max entalphy generated. Max entalphy, `H_(max) = H_(r) - H_(p)`, for this computation the coefficients of both reactant species as well as product species is taken at 298.15 K i.e, coefficients in the low temperature range is considered.
Codes for the same is shown below,
"""
Determining the AFT of CH4 considering entalphy loss
"""
"""
CH4 + 2(O2 + 3.76N2) = CO2 + 2H2O + 7.52N2 - Stoichometric equation for combustion of CH4
"""
import matplotlib.pyplot as plt
"""
Defining a function to calculate entalphy
"""
def h(T,coeff):
a1 = coeff[0]
a2 = coeff[1]
a3 = coeff[2]
a4 = coeff[3]
a5 = coeff[4]
a6 = coeff[5]
R = 8.314 # J/mol*K, Universal gas constant
return (a1 + a2/2*T + a3/3*pow(T,2) + a4/4*pow(T,3) + a5/5*pow(T,4) + a6/T)*R*T
"""
Defining coefficients for reactants
"""
coeff_CH4_r = [5.14987613E+00,-1.36709788E-02,4.91800599E-05,-4.84743026E-08,1.66693956E-11,-1.02466476E+04]
coeff_O2_r = [3.78245636E+00,-2.99673416E-03,9.84730201E-06,-9.68129509E-09,3.24372837E-12,-1.06394356E+03]
coeff_N2_r = [0.03298677E+02,0.14082404E-02,-0.03963222E-04,0.05641515E-07,-0.02444854E-10,-0.10208999E+04]
"""
Defining coefficients for products
"""
coeff_CO2_p = [3.85746029E+00,4.41437026E-03,-2.21481404E-06,5.23490188E-10,-4.72084164E-14,-4.87591660E+04]
coeff_H2O_p = [3.03399249E+00,2.17691804E-03,-1.64072518E-07,-9.70419870E-11,1.68200992E-14,-3.00042971E+04]
coeff_N2_p = [0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04]
"""
Defining coefficients for products at low temp. range
"""
coeff_CO2_l = [2.35677352E+00,8.98459677E-03,-7.12356269E-06,2.45919022E-09,-1.43699548E-13,-4.83719697E+04]
coeff_H2O_l = [4.19864056E+00,-2.03643410E-03,6.52040211E-06,-5.48797062E-09,1.77197817E-12,-3.02937267E+04]
coeff_N2_l = [0.03298677E+02,0.14082404E-02,-0.03963222E-04,0.05641515E-07,-0.02444854E-10,-0.10208999E+04]
"""
Calculating entalphy of reactants
"""
T_std = 298.15 # Initial temperature before combustion
h_CH4_r = h(T_std,coeff_CH4_r)
h_O2_r = h(T_std,coeff_O2_r)
h_N2_r = h(T_std,coeff_N2_r)
h_react = h_CH4_r*1 + h_O2_r*2 + 7.52*h_N2_r
"""
Defining a function to calculate entalphy of products
"""
def f(T):
h_CO2_p = h(T,coeff_CO2_p)
h_H2O_p = h(T,coeff_H2O_p)
h_N2_p = h(T,coeff_N2_p)
h_prod = h_CO2_p + 2*h_H2O_p + 7.52*h_N2_p
h_max = h_react - h(T_std,coeff_CO2_l) - (2*h(T_std,coeff_H2O_l)) - (7.52*h(T_std,coeff_N2_l))
return h_prod - h_react + 0.35*h_max
"""
Defining a function for fprime
"""
def fprime(T):
return (f(T + 1e-6) - f(T))/1e-6
tol = 1e-4 # Tolerance
T_aft = 3000 # Initial guess temperature
alpha = 0.2 # Relaxation factor
iter = 0
while(abs(f(T_aft)) > tol):
T_aft = T_aft - (alpha*f(T_aft)/fprime(T_aft))
iter = iter + 1
plt.plot(iter,T_aft,'*',color = 'r')
plt.xlabel('No. of iterations')
plt.ylabel('T_aft')
plt.title('AFT variation with h_loss')
print(T_aft)
print(iter)
plt.show()
Executing the above command gives AFT of CH4 with h_loss = 0.35*h_max. A snip of the AFT obtained is shown below,
so the AFT = 1678.411 K and it took 101 iterations using newton raphson method to get a converged solution. A plot of variation of AFT value with the Number of iteration is plotted which is as shown below,
This result is verified, like without H_loss the AFT of CH4 = 2325.49 K so with a H_loss = 0.35 * H_max giving AFT = 1678.41 K is promising.
3. To find out AFT of 2 carbon atom variants of Alkane, Alkene and Alkyne considering a constant pressure reaction.
Logic:
* As far as defining the functions are concerned it is same as mentioned above.
* A for loop is defined in the range of 1 to 4, a variable c is used to store the loop counter. When c = 1 AFT calculations are carried for ethane. Similarly when c = 2 and c = 3, AFT of ethene and ethyne is calculated.
Codes for the same is shown below,
"""
Determining AFT for 2 carbon variants of general hyrdrocarbons
"""
"""
C2H6 + 3.5(O2 + 3.76N2) = 2CO2 + 3H2O + 13.16N2 - Stoichometric combustion of ethane
C2H4 + 3(O2 + 3.76N2) = 2CO2 + 2H2O + 11.28N2 - Stoichomatric combustion of ethene
C2H2 + 2.5(O2 + 3.76N2) = 2CO2 + H2O + 9.4N2 - Stoichometric combustion of ethyne
"""
"""
Defining a function to calculate entalphy
"""
def h(T,coeff):
a1 = coeff[0]
a2 = coeff[1]
a3 = coeff[2]
a4 = coeff[3]
a5 = coeff[4]
a6 = coeff[5]
R = 8.314 #J/mol*K - Universal gas constant
return (a1 + a2*T/2 + a3/3*pow(T,2) + a4/4*pow(T,3) + a5/5*pow(T,4) + a6/T)*R*T
"""
Defining coefficients of reactants
"""
coeff_C2H6_r = [4.29142492E+00,-5.50154270E-03,5.99438288E-05,-7.08466285E-08,2.68685771E-11,-1.15222055E+04]
coeff_O2_r = [3.78245636E+00,-2.99673416E-03,9.84730201E-06,-9.68129509E-09,3.24372837E-12,-1.06394356E+03]
coeff_N2_r = [0.03298677E+02,0.14082404E-02,-0.03963222E-04,0.05641515E-07,-0.02444854E-10,-0.10208999E+04]
coeff_C2H4_r = [3.95920148E+00,-7.57052247E-03,5.70990292E-05,-6.91588753E-08,2.69884373E-11,5.08977593E+03]
coeff_C2H2_r = [8.08681094E-01,2.33615629E-02,-3.55171815E-05,2.80152437E-08,-8.50072974E-12,2.64289807E+04]
"""
Defining coefficients of products
"""
coeff_CO2_p = [3.85746029E+00,4.41437026E-03,-2.21481404E-06,5.23490188E-10,-4.72084164E-14,-4.87591660E+04]
coeff_H2O_p = [3.03399249E+00,2.17691804E-03,-1.64072518E-07,-9.70419870E-11,1.68200992E-14,-3.00042971E+04]
coeff_N2_p = [0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04]
"""
Defining a matrix for AFT values
"""
AFT_value = []
"""
Defining a function to calculate AFT
"""
def f(T):
# Calculating entalphy of reactants
T_std = 298.15 # Ambient temperature
h_C2H6_r = h(T_std,coeff_C2H6_r)
h_O2_r = h(T_std,coeff_O2_r)
h_N2_r = h(T_std,coeff_N2_r)
h_C2H4_r = h(T_std,coeff_C2H4_r)
h_C2H2_r = h(T_std,coeff_C2H2_r)
# Calculating entalphy of products
h_CO2_p = h(T,coeff_CO2_p)
h_H2O_p = h(T,coeff_H2O_p)
h_N2_p = h(T,coeff_N2_p)
if c == 1:
h_react = h_C2H6_r + 3.5*h_O2_r + 13.16*h_N2_r
h_prod = 2*h_CO2_p + 3*h_H2O_p + 13.16*h_N2_p
elif c == 2:
h_react = h_C2H4_r + 3*h_O2_r + 11.28*h_N2_r
h_prod = 2*h_CO2_p + 2*h_H2O_p + 11.28*h_N2_p
else:
h_react = h_C2H2_r + 2.5*h_O2_r + 9.4*h_N2_r
h_prod = 2*h_CO2_p + 1*h_N2_p + 9.4*h_N2_p
return (h_prod - h_react)
"""
Defining fprime
"""
def fprime(T):
return (f(T + 1e-6) - f(T))/1e-6
tol = 1e-4 # Tolerance
T_aft = 3000 # Initial guess temperature
alpha = 0.2 # relaxation factor
for i in range(1,4):
c = i # c=1 --> Ethane, c=2 --> Ethene, c=3 --> ethyne
while(abs(f(T_aft)) > tol):
T_aft = T_aft - (alpha*f(T_aft)/fprime(T_aft))
AFT_value.append(T_aft)
print(AFT_value)Printed AFT value is :
where AFT of ethane = 2379.849 K
AFT of ethene = 2564.533 K
AFT of ethyne = 2497.267 K
Above result can be verified using the cantera module. The validation can be done in a similar way as we verified the variation of AFT of CH4 with `phi`.
Logic:
* First the cantera module needs to be imported as ct using the command 'import cantera as ct'
* Same .TPX and .equilibriate command is used to calculate the AFT
* A for loop in range of 1 to 4 is used and a variable c is used to store the loop counter. At c = 1 AFT of ethane is calculated similarly at c = 2 and c = 3 AFT of ethene and ethyne is calculated.
Codes for the same is shown below,
"""
Determining the AFT of 2 carbon variants of general alkane using cantera module
"""
"""
C2H6 + 3.5(O2 + 3.76N2) = 2CO2 + 3H2O + 13.16N2 - Stoichometric combustion of ethane
C2H4 + 3(O2 + 3.76N2) = 2CO2 + 2H2O + 11.28N2 - Stoichomatric combustion of ethene
C2H2 + 2.5(O2 + 3.76N2) = 2CO2 + H2O + 9.4N2 - Stoichometric combustion of ethyne
"""
import cantera as ct
gas = ct.Solution('gri30.xml')
"""
Defining a matrix for AFT values
"""
AFT_value = []
for i in range(1,4):
c = i/1 # Defining a variable to store the loop counter
# C = 1 --> Ethane, C = 2 --> Ethene, C = 3 --> Ethyne
if c == 1:
gas.TPX = 298.15,101325,{'C2H6':1,'O2':3.5,'N2':13.16}
gas.equilibrate('HP','auto') # Reaction is constant pressure with criteria H_p = H_r
AFT_value.append(gas.T)
elif c == 2:
gas.TPX = 298.15,101325,{'C2H4':1,'O2':3,'N2':11.28}
gas.equilibrate('HP','auto')
AFT_value.append(gas.T)
else:
gas.TPX = 298.15,101325,{'C2H2':1,'O2':2.5,'N2':9.4}
gas.equilibrate('HP','auto')
AFT_value.append(gas.T)
print(AFT_value)
Executing the above command prints the AFT. A snip of the same given below,
* Comparing the AFT values obtained from code as well as cantera there is a difference and this difference is because the cantera uses many minor product species.
4. To plot the effect of No. of carbon atom of general alkane (1 to 3) for H_loss = 0.
Logic:
* Defined functions for AFT calculation and entalphy calculation is same.
* A for loop is defined in the range 1 to 4. A variable c is defined to store the loop counter. At c =1 AFT calculation is done for Methane, at c = 2 and c = 3 AFT calculations are done for ethane and propane.
Codes for the same as follows,
""""
Determining effect of carbon atoms of a general alkane on AFT
"""
"""
C_n H_2n+2 + (3*n+1)/2(O2 + 3.76N2) = nCO2 + (n+1)H2O + (3*n+1)/2*3.76N2 # General stoichometric equation for alkane
"""
import matplotlib.pyplot as plt
"""
Defining a function to calculate entalphy
"""
def h(T,coeffs):
a1 = coeffs[0]
a2 = coeffs[1]
a3 = coeffs[2]
a4 = coeffs[3]
a5 = coeffs[4]
a6 = coeffs[5]
R = 8.314 # J/mol*K
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 coefficients of reactants
"""
coeff_CH4_r = [5.14987613E+00,-1.36709788E-02,4.91800599E-05,-4.84743026E-08,1.66693956E-11,-1.02466476E+04]
coeff_O2_r = [3.78245636E+00,-2.99673416E-03,9.84730201E-06,-9.68129509E-09,3.24372837E-12,-1.06394356E+03]
coeff_N2_r = [0.03298677E+02,0.14082404E-02,-0.03963222E-04,0.05641515E-07,-0.02444854E-10,-0.10208999E+04]
coeff_C2H6_r = [4.29142492E+00,-5.50154270E-03,5.99438288E-05,-7.08466285E-08,2.68685771E-11,-1.15222055E+04]
coeff_C3H8_r = [0.93355381E+00,0.26424579E-01,0.61059727E-05,-0.21977499E-07,0.95149253E-11,-0.13958520E+05]
"""
Defining coefficients of products
"""
coeff_CO2_p = [3.85746029E+00,4.41437026E-03,-2.21481404E-06,5.23490188E-10,-4.72084164E-14,-4.87591660E+04]
coeff_H2O_p = [3.03399249E+00,2.17691804E-03,-1.64072518E-07,-9.70419870E-11,1.68200992E-14,-3.00042971E+04]
coeff_N2_p = [0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04]
"""
Calculating entaphy of reactants
"""
T_std = 298.15 # Known temperature value
h_CH4_r = h(T_std,coeff_CH4_r)
h_O2_r = h(T_std,coeff_O2_r)
h_N2_r = h(T_std,coeff_N2_r)
h_C2H6_r = h(T_std,coeff_C2H6_r)
h_C3H8_r = h(T_std,coeff_C3H8_r)
"""
Defining a function to calculate AFT
"""
def f(T):
h_CO2_p = h(T,coeff_CO2_p)
h_H2O_p = h(T,coeff_H2O_p)
h_N2_p = h(T,coeff_N2_p)
if c == 1:
h_react = h_CH4_r + ((3*c)+1)/2*h_O2_r + ((3*c)+1)/2*3.76*h_N2_r
h_prod = c*h_CO2_p + (c+1)*h_H2O_p + ((3*c)+1)/2*3.76*h_N2_p
elif c == 2:
h_react = h_C2H6_r + ((3*c)+1)/2*h_O2_r + ((3*c)+1)/2*3.76*h_N2_r
h_prod = c*h_CO2_p + (c+1)*h_H2O_p + ((3*c)+1)/2*3.76*h_N2_p
else:
h_react = h_C3H8_r + ((3*c)+1)/2*h_O2_r + ((3*c)+1)/2*3.76*h_N2_r
h_prod = c*h_CO2_p + (c+1)*h_H2O_p + ((3*c)+1)/2*3.76*h_N2_p
return (h_prod - h_react)
"""
Defining fprime
"""
def fprime(T):
return (f(T + 1e-6) - f(T))/1e-6
tol = 1e-4 # Tolerance
T_aft = 3000 # Initial guess value of temperature in kelvin
alpha = 0.2 # Relaxation factor
""""
Defining matrix for No. of carbon atoms
"""
Number_of_C_atom = []
"""
Defining matrix for AFT values
"""
AFT_value = []
for i in range (1,4):
c = i/1
Number_of_C_atom.append(c)
while(abs(f(T_aft)) > tol):
T_aft = T_aft - (alpha*f(T_aft)/fprime(T_aft))
AFT_value.append(T_aft)
print(AFT_value)
print(Number_of_C_atom)
plt.plot(Number_of_C_atom,AFT_value)
plt.xlabel('Number of carbon atoms')
plt.ylabel('AFT_value')
plt.title('No. of carbon atoms v/s AFT')
plt.show()
when above codes are executed the AFT values are printed
the AFT value is plotted against the No. of carbon atom. The plot is shown below,
The above result can be verified using cantera module. .TPX and .equilibrate functions are used along with a for loop codes for the same is shown below,
"""
Determining variation of AFT with No. of carbon atoms for general alkane using cantera module
"""
"""
C_n H_2n+2 + (3n+1)/2*(O2 + 3.76N2) = nCO2 + (n+1)H2O + (3n+1)/2*3.76*N2
"""
import cantera as ct
import matplotlib.pyplot as plt
gas = ct.Solution('gri30.xml')
"""
Defining a matrix for No.of carbon atoms
"""
Number_of_carbon_atoms = []
"""
Defining a matrix for AFT values
"""
AFT_value = []
for i in range(1,4):
c = i/1 # Defining a variable for No. of carbon atoms
Number_of_carbon_atoms.append(c)
if c == 1:
gas.TPX = 298.15,101325,{'CH4':1,'O2':((3*c)+1)/2,'N2':((3*c)+1)/2*3.76}
gas.equilibrate('HP','auto')
AFT_value.append(gas.T)
elif c == 2:
gas.TPX = 298.15,101325,{'C2H6':1,'O2':((3*c)+1)/2,'N2':((3*c)+1)/2*3.76}
gas.equilibrate('HP','auto')
AFT_value.append(gas.T)
else:
gas.TPX = 298.15,101325,{'C3H8':1,'O2':((3*c)+1)/2,'N2':((3*c)+1)/2*3.76}
gas.equilibrate('HP','auto')
AFT_value.append(gas.T)
print(AFT_value)
print(Number_of_carbon_atoms)
plt.plot(Number_of_carbon_atoms,AFT_value)
plt.xlabel('No. of carbon atoms')
plt.ylabel('AFT')
plt.title('No. of carbon atoms of alkane v/s AFT using cantera module')
plt.show()
from the above codes we get AFT values as,
plot obtained as follows,
Both the plot looks same and hence the result obtained from codes can be verified.