Week 3 - Adiabatic Flame Temperature calculation

AIM -  

To calculate the adiabatic flame temperature using python and cantera

 

OBJECTIVES - 

To calculate the adiabatic flame temperature in a constant volume chamber with variation in rquivalence ratio.

To calculate the flame temperature with variation in the heat loss

To calcualte the adiabatic flame temperature variation with respect to alkane, alkene and alkyne

To calcualte the adaibatic flame temperature variation with respect to number of carbon atoms

 

INTRODUCTION -

Adiabatic flame temperature is the temperature of the flame that is generated when a given amount of fuel in burnt completely inside an adiabatic chamber either with constant volume or constant pressure.

To calculate the adiabatic flame temperature we will be using a method called as Newton Raphson Method. In this method good approximation for the roots of a real valued function can be found out quickly. It uses the idea that a continuous and differentiable function can be approximated by a straigh 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 underset(o)()`​. Then Newton's method tells us that a better approximation for the root is 

`x underset(1)() = x underset(o)() - f (x underset(o)()) / (f'(x underset(o)()))`

This process is repeated until the desired accuracy is reached , in general for any x having value `x underset(n)()` the next value can be given by 

`x underset(n+1)() = x underset(n)() - f (x underset(n)()) / (f'(x underset(n)()))`

the graphical representation of the Newton-Raphson method is shown below 

 

NASA Polynomials to calulate the thermodynamic properties :

To calculate the thermodynamic properties like the specific heat, enthalpy and entropy we will be using the nasa polynomials that are shown below

in this polynomials the coefficients a1, a2, a3 untill a7 can be found in the Thermodynamic CHEMKIN file from the website of the university of Berkeley. this file contains the data for the coefficients for alot of chemical species , a small number of species can be seen in the following image 

there are two sets of coefficients for each species one set is for temperature above 1000K and another set is for temperature below 1000K this coefficients are listed in the row 2,3,4 as numbered on the right side. 

 

 

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

In this part we will be :

• Assuming that methane is contained in a constant volume chamber, we will write a program in Python to plot the effect of equivalence ratio on the final adiabatic flame temperature. 
• We will also calculate the results with Cantera for varying equivalence ratios using the equilibrate method assuming that N2 is present. 
• Finally we will compare the results of AFT  (at different equivalence ratio) for python and cantera by plotting in a single graph. and explain why AFT of both differs.

For Methane as fuel the chemical reation will be 

CH4 + 2 (O2 + 3.76 N2) => CO2 + 2 H2O + 7.52 N2

 

here we have a constant volume chamber so `Delta V = 0`

enthalpy : `H = U + PV`

differenetiating enthalpy : `Delta H = Delta U + Delta P. V`...............(1)

also from the first law : `d U = dQ - P Delta V`

since we have a constant volume and adiabatic chamber `d U = 0`............(2)

from (1) and (2) we get 

`Delta H - V Delta P =0`........(3)

now we know 

`PV = n Runderset(u)() T`

so for the reactants we have `P underset(i)() V = n underset(r)() R underset(u)() T underset(i)()`

and for the products `P underset(f)()V = n underset(p)() R underset(u)() T underset(f or ad)`

So the equation (3) becomes

`H underset(r)() - H underset(p)() - R underset(u)() .(n underset(r)() Tunderset(i)() - n underset(p)() T underset(ad)() ) = 0`

So we have our equation that we will be using in the Newton Raphson method to find the roots that is the adiabatic flame temperature  

also in this case we will be having a range of equivalence ratio so that we can observe that variation in the adiabatic flame temperature with change in the equivalence ratio.

`sf "equivalence ratio "(phi) = (sf"fuel" / sf"air") / ((sf"fuel" / sf"air")| underset(sf"stoichiometric")())`

 

When thare is excess of air present in the mixture , the equivalence ratio is lower than 1 and the mixture of fuel and air is called lean mixture and oxygen is present in the products part as well. 

For Lean Mixture : CH4 + x (O2 + 3.76 N2) => a CO2 + b H2O + 3.76*x N2 + c O2

 `phi = (1/x)/(1/2)`

therefore , `x = 2/ phi`

carbon : a = 1,

hydrogen: b = 2, 

oxygen: 2x = 2a + b + 2c     =>    c = x - 2    =>   c = `2/phi`- 2

 

the equation for lean mixture (`phi < 1`) becomes :

CH4 + `2/phi` ( O2 + 3.76 N2)    =>   CO2 + 2 H2O + `7.52/phi` N2 +  `2/phi -2` O2

 

When there is excess of fuel present in the mixture , the equivalence ratio becomes greater than 1 and the mixture is called rich mixture and we get incompletely burnt fuel as cabon monooxide on the products side 

For Rich Mixture : CH4 + x (O2 + 3.76 N2)   =>   a CO2 +  b H2O + 3.76*x N2 + c CO

from the equivalence ratio: x = `2/phi`

carbon: 1 = a+c

hydrogen: 4=2b  =>  b=2 

oxygen: 2x = 2a + b + c  =>  2*`2/phi` = a + b + (a+c)  =>  a = `4/phi`- 3

c  = 4 - `4/phi`

 

the equation for rich mixture (`phi >1`) becomes :

CH4 + `2/phi` ( O2 + 3.76 N2)    => `4/phi`- 3  CO2 + 2 H2O + `7.52/phi` N2 +  4 - `4/phi` O2

 

To calculate the AFT using cantera first we will import the cantera library than use the 'gri30.xml' Solution to assign gas as methane here we will give the STP conditions and the mole factions as the inputs and use the equilibrate method to find the AFT and for the parameter we will use UV as the internal energy and the volume is constant in ur combustion of constant volume chamber and use the solver parameter as 'auto'.

 

The python code for the calculation of Adiabatic flame temperature is given below:

""" 
	In  this code we are going to write the program to calculate the 
	Adiabatic Flame Temperature(AFT) at constant volume at different 
	equivalence ratios

	Chemical equation
	For Stiochometric ratio :: CH4 + 2 (O2 + 3.76 N2) => CO2 + 2 H2o + 7.52 N2
	For lean mixture:: CH4 + 2/eq (O2 + 3.76 N2) => CO2 + 2 H2o + (7.52/eq) N2 + (2/eq - 2) O2
	For Rich mixture:: CH4 + 2/eq (O2 + 3.76 N2) => (4/eq - 3) CO2 + 2 H2O + (7.52/eq) N2 + (4-4/eq) CO

	Program by - Amol Patel

"""

import matplotlib.pyplot as plt
import math
import numpy as np 

# universal gas constant
R = 8.314 #[J/mol-K]

def h (T, co_effs):
	""" 
		Function to calculate the Enthalpy 
	"""
	a1 = co_effs[0]
	a2 = co_effs[1]
	a3 = co_effs[2]
	a4 = co_effs[3]
	a5 = co_effs[4]
	a6 = co_effs[5]
	a7 = co_effs[6]

	return (a1 + a2*T/2 + a3*pow(T,2)/3 + a4*pow(T,3)/4 + a5*pow(T,4)/5 + a6/T)*R*T

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]

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

n2_coeffs_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]

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

# eq is the notation for equivalence ratio

def f(T,eq):
	""" 
		Function that represents root finding problem 
	"""
	# temp at STP
	Tstd = 298.15 #[K]

	# enthaly for each reactant
	h_ch4_r = h(Tstd,ch4_coeffs_l)
	h_o2_r = h(Tstd,o2_coeffs_l)
	h_n2_r = h(Tstd,n2_coeffs_l)

	# total enthalpy of the reactants
	H_reactants = h_ch4_r + 2/eq*(h_o2_r + 3.76*h_n2_r)
	
	# number of moles of reactant
	N_reactants = 1+2/eq*(1+3.76) 


	# enthaly for each product 
	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)


	#for lean mixture
	if eq <= 1:
		# total enthalpy of products for lean mixture
		H_products = h_co2_p + 2*h_h2o_p + 7.52/eq*h_n2_p +(2/eq -2)*h_o2_p
		
		# total number of moles of product 
		N_products = 1 + 2 + 7.52/eq + 2/eq- 2

	#for rich mixture 	
	else:
		# total enthalpy of products for rich mixture
		H_products = (4/eq - 3)*h_co2_p + 2*h_h2o_p + 7.52/eq*h_n2_p + (4 - 4/eq)*h_co_p

		# total number of moles of product 
		N_products = 4/eq -3 + 2 + 7.52/eq + 4 - 4/eq

	return (H_reactants - H_products - R*(N_reactants*Tstd - N_products*T))


def fprime(T,eq):
	"""
	fuction of find the derivative of the Function f 
	"""
	return (f(T+1e-6,eq) - f(T,eq))/1e-6

# temperature array for the python part
T_python = []

# assigning a range for equivalence ratio
eq = np.linspace(0.2,2.0,50)

for i in eq:
	#guess value for temp
	T_guess = 1000

	alpha = 0.5 # relaxation factor
	iter = 0 # iteration number
	tol = 1e-6 # tollerence

	# Newton_Raphson Method
	while (abs(f(T_guess,i) > tol)):
		T_guess = T_guess - alpha*(f(T_guess,i)/(fprime(T_guess,i)))

		iter = iter + 1
	# adding the temp for each value of equivalence ratio in the temperature array
	T_python.append(T_guess)
	
	# plotting the variation of temperature with change in equivalence ratio
	plt.plot(i,T_guess,"*", color = "red")
	
print("Temperature by Python: \n",T_python)



################################## CANTERA-PART #######################################################

import cantera as ct

#temperature array for the cantera part
T_ct = []

for j in eq:

	# assigning gas a methane
	gas = ct.Solution("gri30.xml")
	# assigning the STP value and the mole fraction for the gas
	gas.TPX = 298.15, 101325, {'CH4':1 , 'O2':2/j, 'N2':7.52/j}

	# using the gas equilibrate method with constant internal energy and volume along with "auto" solver
	gas.equilibrate('UV', 'auto')
	# adding the value of temperature in the array 
	T_ct.append(gas.T)
	
	# plotting the variation of temperature with change in the equivalence ratio
	plt.plot(j,gas.T,"o", color = 'blue')

print("Temperature by Cantera: \n",T_ct)

# adding labels to plot 	
plt.xlabel("equivalence-ratio")
plt.ylabel("temperature")
plt.grid("on")
plt.show()	

 

after runnig the code we get the output of the plot as shown below 

also the temperature array will be printed on the output window

 

 

Conclusion for part 1:

1. we see that there is maximum devaiton in the temperature for the python and cantera is when the equivalence ratio is equal to 1.
2. the value of temperature is lower in case of cantera due to the application of gibbs minimization that isn't done in the python code 
3. the value for the temperatures in cantera exceeds the pthyon temperature after equivalence ratio of almost 1.3 

 

 

 

PART 2 - Flame temperature variation with Heat loss (In constant pressure-chamber) :

In this part we will :

• Consider a constant pressure reactor with heat loss. The reactor is initially set at STP conditions.
• We will vary the heat loss term by multiplying the “maximum heat loss” with a factor known as H_loss_factor. Example, H_loss_factor = 0.5 would mean a heat loss for that reaction that is equal to 50% of the maximum possible heat loss
• The maximum heat loss for the reaction is given below    

• By varying the H_loss_factor from 0 to 1. Plot variation in Flame temperature w.r.t H_loss_factor

 

For Methane as fuel the chemical reation will be 

CH4 + 2 (O2 + 3.76 N2) => CO2 + 2 H2O + 7.52 N2

 

here we have a constant pressure chamber so `Delta P = 0`

enthalpy : `H = U + PV`

differenetiating enthalpy : `Delta H = Delta U + Delta P. V + P. Delta V`

`therefore Delta H = Delta U + P .Delta V`................(1)

also from the first law : `d U = dQ - P Delta V`...............(2)

from (1) and (2) we get 

`Delta H = Delta Q - P. Delta V + P. Delta V`

`Delta H = Delta Q`

`therefore Delta H - Delta Q= 0`............(3)

in the above equation it is assumed all the heat is lost but for our calcualtion we will have to vary the heat loss factor that is the amount of heat lost from the chamber.

`Delta Q = sf"Loss"`

for maximum heat loss for the reaction is given in the above image has the enthalpy at the STP temperature and it is known as Lower Heating Value (LHV) 

 

The total loss will be the product of heat loss factor and the maximum heat loss 

`sf"Loss" = sf"H_loss_factor" * sf"Maximum heat loss(LHV)"`  

 

now our equation (3) becomes 

`H_r - H_p - sf"Loss"` 

Now we have our equation which is a fucntion of temperature that will be used to calculate the flame temperture using the Newton-Raphson method

 

here we will be varying the heat loss factor from 0 to 1 to see the variation in the flame temperature 

 

The Python code for the same is shown below: 

"""
Program to calculate the Flame Temperature variation with Heat loss in constant pressure chamber

Chemical equation
	CH4 + 2 (O2 + 3.76 N2) => CO2 + 2 H2o + 7.52 N2

Program by - Amol Patel
"""

import matplotlib.pyplot as plt
import math 
import numpy as np

def h(T,co_effs):
	"""
		function to evaluate enthalpy
	"""
	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]
	a7 = co_effs[6]

	return (a1 + a2*T/2 + a3*pow(T,2)/3 + a4*pow(T,3)/4 + a5*pow(T,4)/5 + a6/T )*R*T

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 ]

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]

co2_coeffs_l = [ 2.35677352E+00 ,8.98459677E-03,-7.12356269E-06,2.45919022E-09,-1.43699548E-13,-4.83719697E+04 ,9.90105222E+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]

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

def f(T,H_loss_factor):
	"""
		Function that represents the root finding problem
	"""
	# enthalpy of each products
	h_co2_p = h(T,co2_coeffs_h)
	h_h2o_p = h(T,h2o_coeffs_h)
	h_n2_p = h(T,n2_coeffs_h)
	# total enthalpt of products
	H_products = h_co2_p + 2*h_h2o_p + 7.52*h_n2_p
	
	# temperature at STP
	Tstd = 298.15

	# enthalpy for each reactants
	h_ch4_r = h(Tstd,ch4_coeffs_l)
	h_o2_r = h(Tstd,o2_coeffs_l)
	h_n2_r = h(Tstd,n2_coeffs_l)
	# total enthalpy of the reactants
	H_reactants = h_ch4_r + 2*h_o2_r + 7.52*h_n2_r
	
	# Lower Heating Values
	LHV = (h(Tstd,co2_coeffs_l) + 2*h(Tstd,h2o_coeffs_l) + 7.52*h(Tstd,n2_coeffs_l)) - (h(Tstd,ch4_coeffs_l) + 2*h(Tstd,o2_coeffs_l) + 7.52*h(Tstd,n2_coeffs_l)) 

	# loss in heat genrated 
	Loss = H_loss_factor * LHV 

	return (H_products - H_reactants) - Loss



def fprime(T,H_loss_factor):
	"""
	function to get the derivative of the function f
	"""
	return (f(T+1e-6,H_loss_factor) -f(T,H_loss_factor))/1e-6

# guess values for temperature
T_guess = 1500

# range of heat loss factor is between 0 and 1 
H_loss_factor = np.linspace(0,1,101)

# assigning a empty array for the temperature values to be stored after calculation for each value of heat loss factor 
temp = []
for i in H_loss_factor: 

	tol = 1e-6 # tollerence
	ct = 0 # number of iteration
	alpha = 0.5 # relaxation factor 

	# Newton Raphson Method
	while (abs(f(T_guess,i)) > tol):
		T_guess = T_guess - alpha*((f(T_guess,i)/fprime(T_guess,i)))
		ct = ct + 1 

	# plotting for variation of temperature with the heat loss factor
	plt.plot(i,T_guess,"o", color = 'blue')
	temp.append(T_guess) # adding the the calculated value of temperature to the array

#output as temperature array 
print("temperature: \n",temp)

plt.xlabel("Heat-loss-factor")
plt.ylabel("temperature")
plt.grid("on")
plt.show()

 

the output for the temperature in the output window is shown below

also the output plot for the variation in flame temperature with variation in heat loss factor is shown below

  

 

Conclusioin for part 2 :

1. there is a reduction in the flame temperature with increase in the heat loss factor 
2. when the heat loss factor is 0 the condition is adiabatic
3. when the heat loss factor is 1 all the is removed from the system and we see that the temperature is less than STP condition.

 

 

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

In this part we will

• Test it only for Alk -ane, -ene and -yne variants of 2 Carbon atom chains and plot a trend assuming no heat loss.

For this part, the Chemical reaction for alkane, alkene and alkyne for 2 carbon atoms can be given as 

Alkane Reaction: C2H6 + 7/2 ( O2 + 3.76 N2) => 2 CO2 + 3 H2O + 7/2*3.76 N2

Alkene Reaction: C2H4 + 3 (O2 + 3.76 N2) => 2 CO2 + 2 H2O + 3*3.76 N2

Alkyne Reaction: C2H2 + 5/2 (O2 + 3.76 N2) => 2 CO2 + H2O + 5/2*3.76 N2

Using this above reation assuming constant pressure and the equivalence ratio is 1, we will calculate the Adiabatic Flame Temerature

the equation this is the function of temperature for this case will be as shown below

`Delta H = 0`

`therefore H underset(sf"reactant")() - H underset(sf"product")() = 0`

 

The Python code for the same calculation is shown below:

"""
In this code we are going to write a program to calculate the AFT of alkane, alkene and alkyne 
that is for a 2 carbon atom chain with differnt number of bonds

Chemical equations: 
ALkanes
	C2H6 + 7/2 ( O2 + 3.76 N2 ) => 2 CO2 + 3 H2O + 7/2*3.76 N2

Alkenes 
	C2H4 + 3 ( O2 + 3.76 N2 ) => 2 CO2 + 2 H2O + 3*3.76 N2

Alkynes
	C2H2 + 5/2 ( O2 +3.76 N2 ) => 2 CO2 + H2O + 5/2*3.76 N2


Program by - Amol Patel
"""

import matplotlib.pyplot as plt
import math
import numpy as np

def h(T,co_effs):
	"""
		function to evaluate enthalpy
	"""
	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]
	a7 = co_effs[6]

	return (a1 + a2*T/2 + a3*pow(T,2)/3 + a4*pow(T,3)/4 + a5*pow(T,4)/5 + a6/T )*R*T


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]


def f(T,n_bonds):
	"""
		Function that represents the root finding problem
	"""
	Tstd = 298.15
	# calculating the enthalpy of each species on the reactants side
	h_o2_r = h(Tstd,o2_coeffs_l)
	h_n2_r = h(Tstd,n2_coeffs_l)
	h_c2h6_r = h(Tstd,c2h6_coeffs_l)
	h_c2h4_r = h(Tstd,c2h4_coeffs_l)
	h_c2h2_r = h(Tstd,c2h2_coeffs_l)
	# calculating the enthalpy of each species on the products side 
	h_co2_p = h(T,co2_coeffs_h)
	h_h2o_p = h(T,h2o_coeffs_h)
	h_n2_p = h(T,n2_coeffs_h)
	

	# n_bonds is the number of bonds between carbon atoms
	if n_bonds == 1:
		# for alkane
		# enthalpy of reactant and products
		H_reactants = h_c2h6_r + 7/2*(h_o2_r + 3.76*h_n2_r)
		H_products = 2*h_co2_p + 3*h_h2o_p + 7/2*3.76*h_n2_p


	elif n_bonds == 2:
		# for alkene		
		# enthalpy of reactant and products
		H_reactants = h_c2h4_r + 3*(h_o2_r + 3.76*h_n2_r)
		H_products = 2*h_co2_p + 2*h_h2o_p + 3*3.76*h_n2_p

	else:
		#for alkyne		
		# enthalpy of reactant and products
		H_reactants = h_c2h2_r +5/2*(h_o2_r + 3.76*h_n2_r)
		H_products = 2*h_co2_p + h_h2o_p + 5/2*3.76*h_n2_p


	return H_products - H_reactants


def fprime(T,n_bonds):
	"""
	function to get the derivative of the function f
	"""
	return (f(T+1e-6,n_bonds) -f(T,n_bonds))/1e-6

# guess value of temperature
T_guess = 1500

# giving value for tollerance, iteration(ct) and ralaxation-factor(alpha)
tol = 1e-6
ct = 0 
alpha = 0.5

# giving the values for the number of bonds as 1 , 2 and 3
n_bonds = np.linspace(1,3,3)

# defining temperature array to store the values of temperature corresponding to the number of bonds
Temp = []

for i in n_bonds:
	# newton-raphson method 
	while (abs(f(T_guess,i)) > tol):
		T_guess = T_guess - alpha*((f(T_guess,i)/fprime(T_guess,i)))
		ct = ct +1

	plt.plot(i,T_guess,"*",color = "red")	
	# adding value to the temperature array
	Temp.append(T_guess)

# getting the output as derired 
print("Temp_alkane =" , Temp[0] ,"\nTemp_alkene =", Temp[1], "\nTemp_alkyne =" ,Temp[2])
plt.xlabel("number of bonds")
plt.ylabel("temperature")
plt.show()

  

 after running the code we get the following output for the plot

and in the output window the temperature values are shown as below

 

Conclusion for part 3 :

1. The Adiabatic Flame Temperature increase with increase in the number of bonds between the carbon atoms.
2. As we know that the strength of the triple bond > double bond > single bond the enregy released in alkyne is lot more than that released in alkanes

 

PART 4 - Adiabatic flame temperature variation with respect to number of Carbon atom chains :

In this part we will be

• plotting  the effect of the number of carbon atoms (from 1 to 3) for a general alkane, on the final temperature assuming no heat loss

The chemical reaction for the general alkanes from 1 to 3 carbon atom assuming that we have a constant pressure chamber and the equivalence ratio is 1 , we will calculate the Adiabatic Flame Temperature

The Pyhton code for this calculation is shown below:

"""
In this code we are going to write a program to calculate the AFT of general alkane with carbon atoms from 1 to 3

Chemical equations
C1 
	CH4 + 2 (O2 + 3.76 N2) => CO2 + 2 H2o + 7.52 N2

C2 
	C2H6 + 7/2 ( O2 + 3.76 N2 ) => 2 CO2 + 3 H2O + 7/2*3.76 N2

C3
	C3H8 + 5 (O2 + 3.76 N2) => 3 CO2 + 4 H2O + 5*3.76 N2


Program by - Amol Patel
"""

import matplotlib.pyplot as plt
import math
import numpy as np

def h(T,co_effs):
	"""
		function to evaluate enthalpy
	"""
	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]
	a7 = co_effs[6]

	return (a1 + a2*T/2 + a3*pow(T,2)/3 + a4*pow(T,3)/4 + a5*pow(T,4)/5 + a6/T )*R*T

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]

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]


def f(T,n_carbons):
	"""
		Function that represents the root finding problem
	"""
	Tstd = 298.15
	# calculating the enthalpy of each species on the reactant side
	h_o2_r = h(Tstd,o2_coeffs_l)
	h_n2_r = h(Tstd,n2_coeffs_l)
	h_ch4_r = h(Tstd,ch4_coeffs_l)
	h_c2h6_r = h(Tstd,c2h6_coeffs_l)
	h_c3h8_r = h(Tstd,c3h8_coeffs_l)

	# calculating the enthalpy of each species on the product side
	h_co2_p = h(T,co2_coeffs_h)
	h_h2o_p = h(T,h2o_coeffs_h)
	h_n2_p = h(T,n2_coeffs_h)


	if n_carbons == 1:
		# for C1
		# Total enthalpy of reactants and products
		H_reactants = h_ch4_r + 2*(h_o2_r + 3.76*h_n2_r)
		H_products = h_co2_p + 2*h_h2o_p + 2*3.76*h_n2_p

	elif n_carbons == 2:
		# for C2
		# Total enthalpy of reactants and products
		H_reactants = h_c2h6_r + 7/2*(h_o2_r + 3.76*h_n2_r)
		H_products = 2*h_co2_p + 3*h_h2o_p + 7/2*3.76*h_n2_p

	else:
		# for C3
		# Total enthalpy of reactants and products
		H_reactants = h_c3h8_r +5*(h_o2_r + 3.76*h_n2_r)
		H_products = 3*h_co2_p + 4*h_h2o_p + 5*3.76*h_n2_p


	return H_products - H_reactants


def fprime(T,n_carbons):
	"""
	function to get the derivative of the function f
	"""
	return (f(T+1e-6,n_carbons) -f(T,n_carbons))/1e-6

# guess value for temerature
T_guess = 1500

# giving tollerence, iteration number (ct) and relaxation factor (alpha)
tol = 1e-6
ct = 0 
alpha = 0.5

# number of carbon atoms
n_carbons = np.linspace(1,3,3)

# assigning empty array for temperature values
Temp = []

for i in n_carbons:

	# Newton-Raphson Method
	while (abs(f(T_guess,i)) > tol):
		T_guess = T_guess - alpha*((f(T_guess,i)/fprime(T_guess,i)))
		#plt.plot(ct,T_guess,"*", color = "red")
		#print(T_guess)
		ct = ct +1
	plt.plot(i,T_guess,"*",color = "red")	
	# addign calcualted temperature values in the temp array
	Temp.append(T_guess)

# getting the desired output
print("Temp_c1 = ",Temp[0],"\nTemp_c2 =",Temp[1], "\nTemp_c3 =", Temp[2])
plt.xlabel("number of carbons")
plt.ylabel("temperature")
plt.show()

 

the output in the from of plot can be seen as given below

also the temperature values in the output window is shown below

 

Conclusion for part 4 :

1. we can see as the number of carbon atom increases the AFT values also increases but the growth in temperature values is not linear.

 

CONCLUSIONS -

• We can see that using the python code we get higher value of Adiabatic Flame temperature which is because the gibbs minimization has not been applied in our code while in cantera it is already applied , this cause our result to higher than actual
• We can see the behavior of Flame temperature due to the variation in the heat loss factor, when the heat loss factor is 1 we get a temperature value lower than the STP condition, it is beacuse we havent used the intermidiate reaction that take place during combustion and also the species formed in this reactions are not taken into account so we have this error.
• We can observe that the adiabatic flame temperature is  highest for the alkyne and lowest for the alkane.
• We can observe that the higher temperature is attained with increase in the number of carbon atom for the general alkane series.