Week 3 - Adiabatic Flame Temperature calculation

                                         ADIABATIC FLAME TEMPERATURE CALCULATION(using Cantera and Python)

 

Objective;

1. Adiabatic flame temperature variation with Equivalence ratio (In constant volume-chamber).

2. Adiabatic flame temperature variation with Heat loss (In constant pressure-chamber).

3. Adiabatic flame temperature variation with respect to Alkane, Alkene, Alkyne.

4. Adiabatic flame temperature variation with respect to a number of Carbon atom chains.

 

Theory;

The adiabatic flame temperature (AFT) is defined as the temperature attained when all of the chemical reaction heat released heats combustion products or the maximum temperature at constant volume or pressure that is obtained after the fuel is completely burned.This temperature is obtained from the system where there is no heat transfer.

                                                         

Adiabatic Flame Temperature can be splited to two cases ;Constant volume and Constant Pressure

Constant volume:Combustion leads to the rise in entropy.

Constant Pressure: Combustion leads to the rise in enthalpy.

both of which describe temperature that combustion products theoretically can reach if no energy is lost to the outside environment.Inother words The constant volume adiabatic flame temperature is the temperature that results from a complete combustion process that occurs without any work,heat transferr or changes in kinetic or potential energy. Its temperature is higher than the constant pressure process because no energy is utilized to change the volume of the system. 

 

 In constant pressure:

Hreactants(Tinitial,P)=Hproducts(Tadiabatic,P)

but in case of constant volume the combustion has extra terms to be considered:

here Ureactants(Tinitial,Pinitial)=Uproducts(Tadiabatic,Pfinal)

rearranging the equation gives us:

Hreactants-Hproducts-V*(Pinitial-Pfinal)=0------1

using real gas equation PV=NRT ,since V is constant equation 1 can be arranged as:

Hreactants-Hproducts-(Nreactants*Tinitial-Nproducts*Tadiabatic)=0-----2

 

LHV:

Lower Heat Value can be said as how much work or energy can be extracted from a fuel.

LHV can be determined by subtracting Enthalpy of reactants and products

 

Newton-Raphson method

Iterative method,Objective is tho find the root of the function F(x)=0

Newton-Raphson method is a way to quickly find a good approximation for the root of a real-valued function $f(x)=0$. It uses the idea that a continuous and differentiable function can be approximated by a straight line tangent to it.

Suppose you need to find the root of a continuous, differentiable function $f(x)$, and you know the root you are looking for is near the point x = x_0x=x0​. Then Newton's method tells us that a better approximation for the root is

                               

This process may be repeated as many times as necessary to get the desired accuracy. In general, for any $x$-value x_nxn​, the next value is given by

                                           

Note: the term "near" is used loosely because it does not need a precise definition in this context. However, x_0x0​ should be closer to the root you need than to any other root (if the function has multiple roots).

 

 

 

Limitation of newton raphson is that the method may not work if there are points of inflection, local maxima or minima around X0​ or the root.

 

Equivalance Ratio:

The equivalence ratio is defined as the ratio of the actual fuel/air ratio to the stoichiometric fuel/air ratio. Stoichiometric combustion occurs when all the oxygen is consumed in the reaction, and there is no molecular oxygen(O2) in the products.

If the equivalence ratio is equal to one, the combustion is stoichiometric. If it is < 1, the combustion is lean with excess air, and if it is >1, the combustion is rich with incomplete combustion.

                                                       

 

Skeleton of the Process/code:

• In the following process tools like numpy(for numerical simulation),Matplotlib.pyplot and cantera are used.
• The use of linespace is to create a sequence of equivalence ratio between 0.1-2 with specific steps
• The enthalpy of products and reactants are found through NASA Polynomial equation in Python
• 
• Low and High coefficients of the reactants and products are defined from GRI mech 3.0
• Enthalpy is calculated using a function
• Using forward difference method derivative of the function is defined
• Apply Newton raphson method by guessing a temperature value
• Plot the respective result with matplotlib
• If so compare the results with Cantera 

Solutions:

 

 1. Adiabatic flame temperature variation with Equivalence ratio (In constant volume-chamber).

we have general Hydrocarbon stochiometric equation as

General eqn:CxHy+ ((x+y)/4)(O2+3.76N2)----xCO2+(y/2)H2O+((X+Y)/4)N2

for lean mixture:CH4+(2/e)(O2+3.76 N2)----CO2+2H20+((2/e)-2)O2+(7.52/e)N2

for rich mixture:CH4+(2/e)(O2+3.76 N2)----((4/e)-3)CO2 +2H2O+(4-(4/e))CO+(7.52/e)N2

where e=equivalance ratio.

 

________________________________________________________________________________________________________________________________________________________________________

#Assume that methane is contained in a constant volume chamber.
import cantera as ct
import matplotlib.pyplot as plt
import numpy as np

#General eqn:CxHy+ ((x+y)/4)(O2+3.76N2)----xCO2+(y/2)H2O+((X+Y)/4)N2

#for lean mixture:CH4+(2/e)(O2+3.76 N2)----CO2+2H20+((2/e)-2)O2+(7.52/e)N2

#for rich mixture:CH4+(2/e)(O2+3.76 N2)----((4/e)-3)CO2 +2H2O+(4-(4/e))CO+(7.52/e)N2

#function to calculate enthalpy of reactants and Products
def h(T,co_effs):
R=8.314 #J/Mol
a1=co_effs[0]
a2=co_effs[1]
a3=co_effs[2]
a4=co_effs[3]
a5=co_effs[4]
a6=co_effs[5]
#nasa Polynomial eqn
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 emperature and High Temperature coefficients of species

ch4_coeffs_l = [5.14987613E+00,-1.36709788E-02, 4.91800599E-05,-4.84743026E-08, 1.66693956E-11,-1.02466476E+04]

o2_coeffs_l = [3.78245636E+00,-2.99673416E-03, 9.84730201E-06,-9.68129509E-09, 3.24372837E-12,-1.06394356E+03]

o2_coeffs_h = [3.28253784E+00,1.48308754E-03, -7.57966669E-07,2.09470555E-10, -2.16717794E-14,-1.08845772E+03]

n2_coeffs_l = [0.03298677E+02, 0.14082404E-02,-0.03963222E-04, 0.05641515E-07,-0.02444854E-10,-0.10208999E+04]

n2_coeffs_h = [0.02926640E+02, 0.14879768E-02,-0.05684760E-05, 0.10097038E-09,-0.06753351E-13,-0.09227977E+04]

co2_coeffs_l = [2.35677352E+00, 8.98459677E-03,-7.12356269E-06, 2.45919022E-09,-1.43699548E-13,-4.83719697E+04]

co2_coeffs_h = [3.85746029E+00, 4.41437026E-03,-2.21481404E-06, 5.23490188E-10,-4.72084164E-14,-4.87591660E+04]

h2o_coeffs_l = [4.19864056E+00,-2.03643410E-03, 6.52040211E-06,-5.48797062E-09, 1.77197817E-12,-3.02937267E+04]

h2o_coeffs_h = [3.03399249E+00, 2.17691804E-03,-1.64072518E-07,-9.70419870E-11, 1.68200992E-14,-3.00042971E+04]

co_coeffs_l = [3.57953347E+00,-6.10353680E-04, 1.01681433E-06, 9.07005884E-10,-9.04424499E-13,-1.43440860E+04]

co_coeffs_h = [2.71518561E+00, 2.06252743E-03,-9.98825771E-07, 2.30053008E-10,-2.03647716E-14,-1.41518724E+04]

#fn to find enthalpy of product and reactants

def f(T,eq):

Tstd =298.15 #k
R= 8.314 #Jol/ mol-k

h_ch4_r = h(Tstd,ch4_coeffs_l)
h_o2_r = h(Tstd,o2_coeffs_l)
h_n2_r = h(Tstd,n2_coeffs_l)
#choosing the coefficients using T>1000 condition

if T>1000:
h_co2_p = h(T,co2_coeffs_h)
h_h2o_p = h(T,h2o_coeffs_h)
h_n2_p = h(T,n2_coeffs_h)
h_o2_p = h(T,o2_coeffs_h)
h_co_p = h(T,co_coeffs_h)
else:
h_co2_p = h(T,co2_coeffs_l)
h_h2o_p = h(T,h2o_coeffs_l)
h_n2_p = h(T,n2_coeffs_l)
h_o2_p = h(T,o2_coeffs_l)
h_co_p = h(T,co_coeffs_l)

if eq<1: # Lean Mixture with reactant product enthalpy and no.of moles
h_reactants = h_ch4_r+ (2/eq)*h_o2_r + (7.52/eq)*h_n2_r
h_products = h_co2_p+ 2*h_h2o_p +((2/eq)-2)*h_o2_p + (7.52/eq)*h_n2_p
N_reactants = 1+ (2/eq)+ (7.52/eq)
N_products = 1+2 +((2/eq)-2)+ (7.52/eq)

else: # Rich Mixture (eq > 1)
h_reactants = h_ch4_r+ (2/eq)*h_o2_r + (7.52/eq)*h_n2_r
h_products = ((4/eq) -3)*h_co2_p+ 2*h_h2o_p +(4-(4/eq))*h_co_p + (7.52/eq)*h_n2_p
N_reactants = 1+ (2/eq)+ (7.52/eq)
N_products = ((4/eq) -3)+ 2+(4-(4/eq))+ (7.52/eq)
return h_reactants-h_products-R*(N_reactants*Tstd-N_products*T)
#numercal derivative of the function
def fprime(T,eq):
return (f(T+1e-6, eq)-f(T, eq))/1e-6
"""" _________________________newton raphson method"""
T_guess=1500
tol=1e-3
alpha=0.8
iter=0
i=np.linspace(0.1,2,200)#creating euivalance ratio series
for x in i:#loop for performing root finding for all equivalance series
while (abs(f(T_guess, x))>tol):

T_guess=T_guess-(alpha*(f(T_guess, x)/fprime(T_guess,x)))
iter=iter+1
plt.plot(x,T_guess,'.',color='black')
print('T_guess=',T_guess)
############################################################CANTERA
gas=ct.Solution('gri30.xml')

j=np.linspace(0.1,2.0,200)

temp = []#creating an empty array to store temperature

for k in j:
gas.set_equivalence_ratio(k, 'CH4', 'O2:1, N2:3.76')
gas.TP=298.15,101325
gas.equilibrate('UV','auto')#const volume v and entropy u
temp.append(gas.T)
print(gas.T)
plt.plot(j,temp,'.',color='red')
plt.xlabel('Equivalence Ratio')
plt.ylabel('Adiabatic Flame Temperature(K)')
plt.title('Effect of Equivalence ratio on the Adiabatic flame temperature at Constant Volume')
plt.legend(['Python=black','Cantera=red'])
plt.grid()
plt.show()

 

Code access:  https://drive.google.com/drive/u/1/folders/1MphgMS7MkOfPHiMZhcTE77WsYTZhTwGE

________________________________________________________________________________________________________________________________________________________                                                        

Observation:

• High AFT at equivalance ratio ==1 in python 
• High AFT at the range of  1 to 1.25 in Cantera
• In one scenario cantera exceeds AFT reading after 1.3
• Since species participation is less in python we can expect higher temperature results than cantera with build in library offers

 

 

________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________________

2. Adiabatic flame temperature variation with Heat loss (In constant pressure-chamber) :

 Consider a constant pressure reactor with heat loss. The reactor is initially set at STP conditions.

The heat loss should be input in your program and must be a fraction of the "maximum possible heat loss". You will control heat loss using a variable called H_loss. For example, H_loss = 0.5 would mean a heat loss that is equal to 50% of the maximum possible heat loss.

The fuel type can be a general hydrocarbon of the form CxHy. Phi and lambda = 1.0. Find the AFT of CH4 with H_loss = 0.35

To find out heat loss we assume that the hydrocarbon is CH4 with equivalance ratio 1.

 

import matplotlib.pyplot as plt
import numpy as np
#function to calculate enthalpy of reactants and Products
def h(T,Coeff):

R= 8.314 #Jol/ mol-k
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
#low emperature and High Temperature coefficients of species
ch4_coeffs_l = [5.14987613E+00,-1.36709788E-02, 4.91800599E-05,-4.84743026E-08, 1.66693956E-11,-1.02466476E+04]

o2_coeffs_l = [3.78245636E+00,-2.99673416E-03, 9.84730201E-06,-9.68129509E-09, 3.24372837E-12,-1.06394356E+03]

o2_coeffs_h = [3.28253784E+00,1.48308754E-03, -7.57966669E-07,2.09470555E-10, -2.16717794E-14,-1.08845772E+03]

n2_coeffs_l = [0.03298677E+02, 0.14082404E-02,-0.03963222E-04, 0.05641515E-07,-0.02444854E-10,-0.10208999E+04]

n2_coeffs_h = [0.02926640E+02, 0.14879768E-02,-0.05684760E-05, 0.10097038E-09,-0.06753351E-13,-0.09227977E+04]

co2_coeffs_l = [2.35677352E+00, 8.98459677E-03,-7.12356269E-06, 2.45919022E-09,-1.43699548E-13,-4.83719697E+04]

co2_coeffs_h = [3.85746029E+00, 4.41437026E-03,-2.21481404E-06, 5.23490188E-10,-4.72084164E-14,-4.87591660E+04]

h2o_coeffs_l = [4.19864056E+00,-2.03643410E-03, 6.52040211E-06,-5.48797062E-09, 1.77197817E-12,-3.02937267E+04]

h2o_coeffs_h = [3.03399249E+00, 2.17691804E-03,-1.64072518E-07,-9.70419870E-11, 1.68200992E-14,-3.00042971E+04]

co_coeffs_l = [3.57953347E+00,-6.10353680E-04, 1.01681433E-06, 9.07005884E-10,-9.04424499E-13,-1.43440860E+04]

co_coeffs_h = [2.71518561E+00, 2.06252743E-03,-9.98825771E-07, 2.30053008E-10,-2.03647716E-14,-1.41518724E+04]

#finding LHV
Tstd = 298.15
LHV=(h(Tstd,ch4_coeffs_l)+2*h(Tstd,o2_coeffs_l)+7.52*h(Tstd,n2_coeffs_l))-(1*h(Tstd,co2_coeffs_l)+2*h(Tstd,h2o_coeffs_l)+7.52*h(Tstd,n2_coeffs_l))

#fn to find enthalpy of product and reactants
def f(T,H_loss):
h_ch4_r = h(Tstd,ch4_coeffs_l)
h_o2_r = h(Tstd,o2_coeffs_l)
h_n2_r = h(Tstd,n2_coeffs_l)

if T>1000:
h_co2_p = h(T,co2_coeffs_h)
h_h2o_p = h(T,h2o_coeffs_h)
h_n2_p = h(T,n2_coeffs_h)
else:
h_co2_p = h(T,co2_coeffs_l)
h_h2o_p = h(T,h2o_coeffs_l)
h_n2_p = h(T,n2_coeffs_l)

h_reactants = h_ch4_r+ 2*h_o2_r + 7.52*h_n2_r
h_products= h_co2_p+ 2*h_h2o_p + 7.52*h_n2_p
#we have
Loss = H_loss* LHV
return h_products - h_reactants + Loss
#numercal derivative of the function
def fprime(T, H_loss):
return (f(T+ 1e-6, H_loss)- f(T,H_loss))/1e-6
#####newton raphson root finding method
T_guess=1500
tol = 1e-3
alpha = 0.5

H_loss=np.linspace(0,1,20)

Tl=[]#creating an empty array to store temperature

for i in H_loss:
while (abs(f(T_guess,i))>tol):
T_guess = T_guess - alpha*((f(T_guess,i)/fprime(T_guess, i)))
Tl.append(T_guess)

plt.plot(H_loss,Tl,color='green')
plt.xlabel('Heat Loss')
plt.ylabel('Flame Temperature(K)')
plt.title('Effect of Heat loss on the Flame Temperature')
plt.grid()
plt.show()
#To find the heat loss at 0.35 we need to access the (N-1)th element in list i.e (7-1)th
print('Hloss at 0.35=',Tl[6],'k')

Code access:  https://drive.google.com/drive/u/1/folders/1MphgMS7MkOfPHiMZhcTE77WsYTZhTwGE

_______________________________________________________________________________________________________________________________________________________________________

                                                                             

Obseravation:

• Increase in Heat loss decreases the Adiabatic flame temperature
• AFT v/s Heat loss is linear and the heat loss at 0.35 is found out to be 1743.1330271543134 k

______________________________________________________________________________________________________________________________________________________________________________________________________

 

 

 

3. Adiabatic flame temperature variation with respect to Alkane, Alkene, Alkyne :

""Test it only for Alk-ane, -ene and -yne variants of 2 Carbon atom chain and plot a trend.""  Assume no heat loss.

#Assume no heat loss.
#the stochiometric equation for combustion of a hydrocarbon can be expressed as
#''''CxHy+((x+y/4))(O2+3.76 N2)------> xCO2 + (y/2)H2o +((x+y/4))N2'''
import cantera as ct
import matplotlib.pyplot as plt
import numpy as np
#function to calculate enthalpy of reactants and Products
def h(T,co_effs):
R=8.314 #J/Mol
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)

#lower and higher heat coefficients;

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]

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 #k'

H_alkane=h(Tstd,c2h6_coeffs_l)
H_alkene=h(Tstd,c2h4_coeffs_l)
H_alkyne=h(Tstd,c2h2_coeffs_l)

#and let H be the list containing these hydrocarbons
H=[H_alkane,H_alkene,H_alkyne]
#defining the function with H
def f(T,H):
h_n2_p=h(T,n2_coeffs_h)
h_co2_p=h(T,co2_coeffs_h)
h_h20_p=h(T,h2o_coeffs_h)

#defining reactants
h_n2_r=h(Tstd,n2_coeffs_l)
h_o2_r=h(Tstd,o2_coeffs_l)
#now for the hydrocarbons we can use the loop condition to choose

if H==H_alkane:
#giving stochiometric coefficients;2 carbon chain so x=2|let x+y/4=a and y/2=b
x=2
y=6
a=(x+y/4)
b=y/2
#hence the reactants will be
H_reactants=H_alkane+(a*(h_o2_r))+(a*3.76*h_n2_r)
H_products=(x*h_co2_p)+(b*h_h20_p)+(a*3.76*h_n2_p)
print('Alkane react|prod=',H_reactants,'|',H_products)
elif H==H_alkene:
x=2
y=4
a=(x+y/4)
b=y/2
#hence the reactants will be
H_reactants=H_alkene+(a*(h_o2_r))+(a*3.76*h_n2_r)
H_products=(x*h_co2_p)+(b*h_h20_p)+(a*3.76*h_n2_p)
print('Alkene react|prod=',H_reactants,'|',H_products)
elif H==H_alkyne:
x=2
y=2
a=(x+y/4)
b=y/2
#hence the reactants will be
H_reactants=H_alkyne+(a*(h_o2_r))+(a*3.76*h_n2_r)
H_products=(x*h_co2_p)+(b*h_h20_p)+(a*3.76*h_n2_p)
print('Alkyne react|prod=',H_reactants,'|',H_products)
return (H_reactants-H_products)
#numercal derivative of the function
def fprime(T, H_loss):
return (f(T+ 1e-6, H_loss)- f(T,H_loss))/1e-6
#####newton raphson root finding method
T_guess=1500
tol = 1e-3
alpha = 0.2
iter=0
temp=[]#creating an empty array to store temperature
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)
fuel=['Alkane','Alkene','Alkyne',]
y=np.array(temp)#to plot on bar graph using numpy array
x=np.array(fuel)
plt.ylabel('AFT')

plt.bar(x,y)
plt.show()

##########################################Cantera

gas=ct.Solution('gri30.xml')
gas.TPX=298.15,101325,{'C2h6':1,'O2':3.5,'N2':3.5*3.72}
gas.equilibrate('HP','auto')
T_alkane=gas.T
#___________
gas.TPX=298.15,101325,{'C2h4':1,'O2':3,'N2':3*3.72}
gas.equilibrate('HP','auto')
T_alkene=gas.T
#------------------
gas.TPX=298.15,101325,{'C2h2':1,'O2':2.5,'N2':2.5*3.72}
gas.equilibrate('HP','auto')
T_alkyne=gas.T

T=[T_alkane,T_alkene,T_alkyne]
print(T)
y1=np.array(T)
x=fuel
plt.bar(x,y1,color='orange')
plt.show()

Code access:  https://drive.google.com/drive/u/1/folders/1MphgMS7MkOfPHiMZhcTE77WsYTZhTwGE

_______________________________________________________________________________________________________________________________________________________________________

in python:(AFT)

['Alkane','Alkene','Alkyne',]=[2379.849457826236k, 2564.5336879576425k, 2909.3545659758656k]

                                                      

 

In cantera:

['Alkane','Alkene','Alkyne']=[2266.9015630814106k, 2376.377705459795k, 2546.813043128299k]

                                               

 

Observation:

• AFT is lower in cantera since it had the advantage of  having more species.
• less carbon atoms more AFT i.e more hydrogen atoms  

 

______________________________________________________________________________________________________________________________________________________________________________________________________

4. Adiabatic flame temperature variation with respect to a number of Carbon atom chains :

 

""For a general alkane, plot the effect of the number of carbon atoms (from 1 to 3) on the final temperature. Assume no heat loss.""

 

Similar to above code just changed the alkane,alkene and alkyne to carbon 1,2 and 3:

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

def h(T,co_effs):
R=8.314 #J/Mol
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)

"""CxHy+((x+Y/4)(O2+3.76 N2)------> xCO2 + (y/2)H2o +((x+y/4)N2 |
The Hydrocarbon in alkanes are written as CnH2n+2
Since assignment gives us 3 carbon atom:
The alkanes include CH4, C2H6 and C3H8"""

#coefficients of rectants and products:
ch4_coeffs_l=[5.14987613E+00,-1.36709788E-02, 4.91800599E-05,-4.84743026E-08,1.66693956E-11,-1.02466476E+04]
c2h6_coeffs_l=[ 4.29142492E+00,-5.50154270E-03,5.99438288E-05,-7.08466285E-08,2.68685771E-11,-1.15222055E+04]
c3h8_coeffs_l=[0.93355381E+00,0.26424579E-01,0.61059727E-05,-0.21977499E-07,0.95149253E-11,-0.13958520E+05 ]

n2_coeffs_h=[0.02926640E+02 ,0.14879768E-02,-0.05684760E-05, 0.10097038E-09,-0.06753351E-13,-0.09227977E+04]
h2o_coeffs_h=[3.03399249E+00,2.17691804E-03,-1.64072518E-07,-9.70419870E-11, 1.68200992E-14,-3.00042971E+04]
co2_coeffs_h=[3.85746029E+00 ,4.41437026E-03,-2.21481404E-06,5.23490188E-10,-4.72084164E-14,-4.87591660E+04]

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]

Tstd=298.15#k
c1=h(Tstd,ch4_coeffs_l)
c2=h(Tstd,c2h6_coeffs_l)
c3=h(Tstd,c3h8_coeffs_l)

carbon=[c1,c2,c3]
def f(T,carbon):
h_n2_p=h(T,n2_coeffs_h)
h_co2_p=h(T,co2_coeffs_h)
h_h20_p=h(T,h2o_coeffs_h)

#defining reactants
h_n2_r=h(Tstd,n2_coeffs_l)
h_o2_r=h(Tstd,o2_coeffs_l)
#now for the hydrocarbons we can use the loop condition to choose

if carbon==c1:
#giving stochiometric coefficients;2 carbon chain so x=2|let x+y/4=a and y/2=b
x=1
y=(x+1)*2
a=(x+y/4)
b=y/2
#hence the reactants will be
H_reactants=c1+(a*(h_o2_r))+(a*3.76*h_n2_r)
H_products=(x*h_co2_p)+(b*h_h20_p)+(a*3.76*h_n2_p)
print('Alkane react|prod=',H_reactants,'|',H_products)
elif carbon==c2:
x=2
y=(x+1)*2
a=(x+y/4)
b=y/2
#hence the reactants will be
H_reactants=c2+(a*(h_o2_r))+(a*3.76*h_n2_r)
H_products=(x*h_co2_p)+(b*h_h20_p)+(a*3.76*h_n2_p)
print('Alkene react|prod=',H_reactants,'|',H_products)
elif carbon==c3:
x=2
y=(x+1)*2
a=(x+y/4)
b=y/2
#hence the reactants will be
H_reactants=c3+(a*(h_o2_r))+(a*3.76*h_n2_r)
H_products=(x*h_co2_p)+(b*h_h20_p)+(a*3.76*h_n2_p)
print('Alkyne react|prod=',H_reactants,'|',H_products)
return (H_reactants-H_products)

def fprime(T, carbon):
return (f(T+ 1e-6,carbon)- f(T,carbon))/1e-6

#Newton Raphson

T_guess=1500
tol = 1e-3
alpha = 0.2
iter=0
temp=[]
for i in carbon:
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)
fuel=['carbon1','carbon2','carbon3',]
y=np.array(temp)
x=np.array(fuel)
plt.ylabel('AFT')

plt.title('AFT v/s No, of Carbon atoms in Hydro carbon in python')

plt.bar(x,y,color='red')
plt.show()

##########################################Cantera

gas=ct.Solution('gri30.xml')
gas.TPX=298.15,101325,{'Ch4':1,'O2':2,'N2':2*3.72}
gas.equilibrate('HP','auto')
T_CH4=gas.T
#___________
gas.TPX=298.15,101325,{'C2h6':1,'O2':3.5,'N2':3.5*3.72}
gas.equilibrate('HP','auto')
T_C2H6=gas.T
#------------------
gas.TPX=298.15,101325,{'C3h8':1,'O2':5,'N2':5*3.72}
gas.equilibrate('HP','auto')
T_C3H8=gas.T

T=[T_CH4,T_C2H6,T_C3H8]
print(T)
y1=np.array(T)
x=fuel
plt.bar(x,y1,color='red')

plt.title('AFT v/s No, of Carbon atoms in Hydro carbon in cantera')
plt.show()

Code access: https://drive.google.com/drive/u/1/folders/1MphgMS7MkOfPHiMZhcTE77WsYTZhTwGE

________________________________________________________________________________________________________________________________________________________________________

plots:

python:

fuel=['carbon1','carbon2','carbon3']=[2325.59812786638k, 2379.8494579179182k, 2353.7066176808203k]

note: here carbon 1 means hydrocarbon with 1 carbon in the alkane family i.e CH4.Respectively 

carbon 1 =CH4

carbon 2=C2H6

carbon 3= C3H8

                                                                                  

 

Cantera:

fuel=['carbon1','carbon2','carbon3']=[2232.8777130628387k, 2266.9015630814106k, 2273.8525645061973k]

                                                                                       

 

Observation

• Similar trend Python gives higher AFT values than Cantera
• As carbon number increases the AFT increases gives a hint about the hydrogens influence

 

 

 

 Concusion:

From this Project we can conclude that Combustion simulation through cantera is more to physical results.Although python using Newton raphson model can approximate the result, for more close and valid results we should truly use cantera,Since it have more data about the species .From this experiment we can also say that representation of species in the combustion matters a lot.This can significantly influence the results.

From this experiment we observed the trends of AFT with various parameters like Heat loss,Constant volume situations,Hydrocarbons and no.of carbon atoms from all these we saw a lesser value for cantera results making us think the influence of forgotten species formed during the combustion recation.

In general for any chemical system the reactions are gonna take place till the point of minimum value i.e we get a minimize gibbs free energy function.This is how Cantera works 

 code access:  https://drive.google.com/drive/u/1/folders/1MphgMS7MkOfPHiMZhcTE77WsYTZhTwGE