Week 3 - Adiabatic Flame Temperature calculation

AIM:To calculate the abiabatic flame temperature for the combustion of fuel using pyhton.

Objective: 1. Use the concept of newton raphsons method to solve the enthapy equation of products and reactants at constant pressure.

2.To monitor the effect of equivalence ratio on adiabatic flame temperature

Adiabatic flame temperature:

When a combustion reaction takes place energy is released to the combustion products. If no heat is lost in this process, the temperature of the combustion products is known as the "Adiabatic Flame Temperature."

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.

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

In constant volume system the internal energy of the reaction is conserved unlike to cinstant pressure system where enthalphy is conserved. Hence the internal energy of reactant is equal to internal energy of product at initial and final state.

H_products-H_reactant-(n_product)*R*T+n_reactant*R*Tstd

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

p=np.arange(0.1,2,0.1)
R=8.314 #JOL/M-K
def h(T,coeff):
	
	a1=coeff[0]
	a2=coeff[1]
	a3=coeff[2]
	a4=coeff[3]
	a5=coeff[4]
	a6=coeff[5]
	return(a1+a2*T/2+a3*pow(T,2)/3+a4*pow(T,3)/4+a5*pow(T,4)/5+a6/T)*R*T
ch4_coeff_l=[5.14987613E+00,-1.36709788E-02,4.91800599E-05,-4.84743026E-08,1.66693956E-11,-1.02466476E+04,-4.64130376E+00]
o2_coeff_l=[3.78245636E+00,-2.99673416E-03,9.84730201E-06,-9.68129509E-09,3.24372837E-12,-1.06394356E+03,3.65767573E+00]
n2_coeff_l=[0.03298677E+02,0.14082404E-02,-0.03963222E-04,0.05641515E-07,-0.02444854E-10,-0.10208999E+04,0.03950372E+02]
n2_coeff_h=[ 0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04,0.05980528E+02]
co2_coeff_h=[3.85746029E+00,4.41437026E-03,-2.21481404E-06,5.23490188E-10,-4.72084164E-14,-4.87591660E+04,2.27163806E+00]
h2o_coeff_h=[3.03399249E+00,2.17691804E-03,-1.64072518E-07,-9.70419870E-11,1.68200992E-14,-3.00042971E+04,4.96677010E+00]
o2_coeff_h=[3.28253784E+00,1.48308754E-03,-7.57966669E-07,2.09470555E-10,-2.16717794E-14,-1.08845772E+03,5.45323129E+00]
co_coeff_h=[2.71518561E+00,2.06252743E-03,-9.98825771E-07,2.30053008E-10,-2.03647716E-14,-1.41518724E+04,7.81868772E+00]
def f(T,phi):
	h_co2_p=h(T,co2_coeff_h)
	h_n2_p=h(T,n2_coeff_h)
	h_h2o_p=h(T,h2o_coeff_h)
	h_o2_p=h(T,o2_coeff_h)
	h_co_p=h(T,co_coeff_h)

	if phi<1:
		H_products=h_co2_p+2*h_h2o_p+(7.52/phi)*h_n2_p+(2/phi-2)*h_o2_p
		n_product=1+(7.52/phi)+2+((2/phi)-2)
	if phi==1:
		H_products=h_co2_p+2*h_h2o_p+7.52*h_n2_p
		n_product=1+2+7.52
	if phi>1:
		H_products=((4/phi)-3)*h_co2_p+2*h_h2o_p+(7.52/phi)*h_n2_p+(4-4/phi)*h_co_p
		n_product=((4/phi)-3)+2+(7.52/phi)+(4-(4/phi))
	
	#print(H_products)

	Tstd=298.15
	h_ch4_r=h(Tstd,ch4_coeff_l)
	h_o2_r=h(Tstd,o2_coeff_l)
	h_n2_r=h(Tstd,n2_coeff_l)

	if phi<1:
		H_reactant=h_ch4_r+(2/phi)*(h_o2_r+3.76*h_n2_r)
		n_reactant=1+(2/phi)+3.76*(2/phi)

	elif phi==1:
		H_reactant=h_ch4_r+2*h_o2_r+7.52*h_n2_r
		n_reactant=1+2+7.52
	elif phi>1:
		H_reactant=h_ch4_r+(2/phi)*(h_o2_r+3.76*h_n2_r)
		n_reactant=1+(2/phi)+3.76*(2/phi)

	#print(H_reactant)

	return H_products-H_reactant-R*(n_product)*T+R*n_reactant*Tstd

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

T_guess=1500
alpha=0.2
tol=1e-3
T_g=[]
for phi in p:
	while(abs(f(T_guess,phi))>tol):
		T_guess=T_guess-alpha*(f(T_guess,phi)/fprime(T_guess,phi))
		print(T_guess)
	T_g.append(T_guess)

plt.plot(p,T_g,color='blue')
plt.xlabel('Equivalance Ratio')
plt.ylabel('Temperature')
plt.show()

Plot:

1.The plot shows clear that maximum adiabatic flame temperature is obtained is 2817.29 at equivalance ratio of 1 

2. The adiabatic flame temperature increases with equivalance ratio as more energy is produces while burning the fuel in excess air. As we move ahead in x axis the aft remperature reduces as excess amount of fuel results in incomplete combustion.

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

when Heat loss at constant pressure is considered the equation to find aft is obtained by equating the enthalpy of product and reactants

H_products-H_reactant+(H_LOSS*Q)

where H_LOSS=heat loss factor and Q=Maximum possible heat loss out of 100

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

HL=np.arange(0.1,1.1,0.1)

def h(T,coeff):
	R=8.314 #JOL/M-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
ch4_coeff_l=[5.14987613E+00,-1.36709788E-02,4.91800599E-05,-4.84743026E-08,1.66693956E-11,-1.02466476E+04,-4.64130376E+00]
o2_coeff_l=[3.78245636E+00,-2.99673416E-03,9.84730201E-06,-9.68129509E-09,3.24372837E-12,-1.06394356E+03,3.65767573E+00]
n2_coeff_l=[0.03298677E+02,0.14082404E-02,-0.03963222E-04,0.05641515E-07,-0.02444854E-10,-0.10208999E+04,0.03950372E+02]
n2_coeff_h=[ 0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04,0.05980528E+02]
co2_coeff_h=[3.85746029E+00,4.41437026E-03,-2.21481404E-06,5.23490188E-10,-4.72084164E-14,-4.87591660E+04,2.27163806E+00]
h2o_coeff_h=[3.03399249E+00,2.17691804E-03,-1.64072518E-07,-9.70419870E-11,1.68200992E-14,-3.00042971E+04,4.96677010E+00]
o2_coeff_h=[3.28253784E+00,1.48308754E-03,-7.57966669E-07,2.09470555E-10,-2.16717794E-14,-1.08845772E+03,5.45323129E+00]
co_coeff_h=[2.71518561E+00,2.06252743E-03,-9.98825771E-07,2.30053008E-10,-2.03647716E-14,-1.41518724E+04,7.81868772E+00]
def f(T,H_LOSS):
	h_co2_p=h(T,co2_coeff_h)
	h_n2_p=h(T,n2_coeff_h)
	h_h2o_p=h(T,h2o_coeff_h)
	h_o2_p=h(T,o2_coeff_h)
	h_co_p=h(T,co_coeff_h)
	H_products=h_co2_p+2*h_h2o_p+7.52*h_n2_p
	
#print(H_products)

	Tstd=298.15
	h_ch4_r=h(Tstd,ch4_coeff_l)
	h_o2_r=h(Tstd,o2_coeff_l)
	h_n2_r=h(Tstd,n2_coeff_l)
	H_reactant=h_ch4_r+2*h_o2_r+7.52*h_n2_r
#print(H_reactant)
#Q=maximum possible heat loss
#H_LOSS=heat loss factor
	Q=h(Tstd,ch4_coeff_l)+2*h(Tstd,o2_coeff_l)-h(Tstd,co2_coeff_h)-2*h(Tstd,h2o_coeff_h)
	return H_products-H_reactant+(H_LOSS*Q)

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

T_guess=1500
alpha=0.2
tol=1e-3
T_g=[]
for H_LOSS in HL:
	while(abs(f(T_guess,H_LOSS))>tol):
		T_guess=T_guess-alpha*(f(T_guess,H_LOSS)/fprime(T_guess,H_LOSS))
		print(T_guess)
	T_g.append(T_guess)

plt.plot(HL,T_g,color='blue')
plt.xlabel('Heat Loss factor')
plt.ylabel('Temperature')
plt.grid('on')
plt.show()

Plot:

The AFT is high as lower heat loss factor as heat loss, it eventually reduces the temperature reached by combustion.

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

The AFT is obtained for Ethane, ethene and ethyne to understand the effect of saturated and unsaturated hydrocarbons on AFT. 

Ethane:

C2H6+3.5(O2+3.75N2) ---> 2CO2+3H2O+13.16N2

Ethene:

C2H4+3(O2+3.75N2) ---> 2CO2+2H2O+11.28N2

Ethyne:

C2H2+2.5(O2+3.75N2) ---> 2CO2+H2O+9.4N2

import math
import matplotlib.pyplot as plt
ane=1
ene=2
yne=3
a=[1,2,3]
def h(T,coeff):
	R=8.314 #JOL/M-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
c2h6_coeff_l=[4.29142492E+00,-5.50154270E-03,5.99438288E-05,-7.08466285E-08,2.68685771E-11,-1.15222055E+04,2.66682316E+00]
c2h4_coeff_l=[3.95920148E+00,-7.57052247E-03,5.70990292E-05,-6.91588753E-08,2.69884373E-11,5.08977593E+03,4.09733096E+00]
c2h2_coeff_l=[8.08681094E-01,2.33615629E-02,-3.55171815E-05,2.80152437E-08,-8.50072974E-12,2.64289807E+04,1.39397051E+01]
o2_coeff_l=[3.78245636E+00,-2.99673416E-03,9.84730201E-06,-9.68129509E-09,3.24372837E-12,-1.06394356E+03,3.65767573E+00]
n2_coeff_l=[0.03298677E+02,0.14082404E-02,-0.03963222E-04,0.05641515E-07,-0.02444854E-10,-0.10208999E+04,0.03950372E+02]
n2_coeff_h=[ 0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04,0.05980528E+02]
co2_coeff_h=[3.85746029E+00,4.41437026E-03,-2.21481404E-06,5.23490188E-10,-4.72084164E-14,-4.87591660E+04,2.27163806E+00]
h2o_coeff_h=[3.03399249E+00,2.17691804E-03,-1.64072518E-07,-9.70419870E-11,1.68200992E-14,-3.00042971E+04,4.96677010E+00]

def f(T,alk):
	h_co2_p=h(T,co2_coeff_h)
	h_n2_p=h(T,n2_coeff_h)
	h_h2o_p=h(T,h2o_coeff_h)

	if alk==ane:
		x=2
		y=6
		H_products=x*h_co2_p+(y/2)*h_h2o_p+13.16*h_n2_p
	if alk==ene:
		x=2
		y=4
		H_products=x*h_co2_p+(y/2)*h_h2o_p+11.28*h_n2_p
	if alk==yne:
		x=2
		y=2
		H_products=x*h_co2_p+(y/2)*h_h2o_p+9.4*h_n2_p

	Tstd=298.15
	h_c2h6_r=h(Tstd,c2h6_coeff_l)
	h_c2h4_r=h(Tstd,c2h4_coeff_l)
	h_c2h2_r=h(Tstd,c2h2_coeff_l)
	h_o2_r=h(Tstd,o2_coeff_l)
	h_n2_r=h(Tstd,n2_coeff_l)

	if alk==ane:
		H_reactant=h_c2h6_r+3.5*h_o2_r+13.16*h_n2_r
	if alk==ene:
		H_reactant=h_c2h4_r+3*h_o2_r+11.28*h_n2_r
	if alk==yne:
		H_reactant=h_c2h2_r+2.5*h_o2_r+9.4*h_n2_r

	return H_products-H_reactant

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

T_guess=1500
alpha=0.2
tol=1e-3
ct=1
for alk in a:
	while(abs(f(T_guess,alk))>tol):
		T_guess=T_guess-alpha*(f(T_guess,alk)/fprime(T_guess,alk))
		ct=ct+1
		plt.plot(alk,T_guess,'*',color='blue')
	print(f(T_guess,alk))
	print(T_guess)
plt.xlabel('Hydrocarbon')
plt.ylabel('Temperature')
plt.show()

Plot:

The Adiabatic flame temperature of Ethene and ethyne is more as compared to ethane, which can be explained by the double and triple bond between carbon atoms which releases requires more energy to break hence increses the temperature.

Following are the AFT for ethane, ethene and ethyne

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

Let us consider Methane, Ethane and Propane for the case having 1,2 and 3 number of carbon atoms

Methane:

CH4+2(O2+3.76N2) ---> CO2+2H2O+7.56N2

Ethane:

C2H6+3.5(O2+3.76N2) ---> 2CO2+3H2O+13.16N2

Propane:

C2H8+5(O2+3.76N2) ---> 3CO2+4H2O+18.8N2

import math
import matplotlib.pyplot as plt
methane=1
ethane=2
propane=3
a=[1,2,3]
def h(T,coeff):
	R=8.314 #JOL/M-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
ch4_coeff_l=[5.14987613E+00,-1.36709788E-02,4.91800599E-05,-4.84743026E-08,1.66693956E-11,-1.02466476E+04,-4.64130376E+00]
c2h6_coeff_l=[4.29142492E+00,-5.50154270E-03,5.99438288E-05,-7.08466285E-08,2.68685771E-11,-1.15222055E+04,2.66682316E+00]
c3h8_coeff_l=[0.93355381E+00,0.26424579E-01,0.61059727E-05,-0.21977499E-07,0.95149253E-11,-0.13958520E+05,0.19201691E+02]
o2_coeff_l=[3.78245636E+00,-2.99673416E-03,9.84730201E-06,-9.68129509E-09,3.24372837E-12,-1.06394356E+03,3.65767573E+00]
n2_coeff_l=[0.03298677E+02,0.14082404E-02,-0.03963222E-04,0.05641515E-07,-0.02444854E-10,-0.10208999E+04,0.03950372E+02]
n2_coeff_h=[ 0.02926640E+02,0.14879768E-02,-0.05684760E-05,0.10097038E-09,-0.06753351E-13,-0.09227977E+04,0.05980528E+02]
co2_coeff_h=[3.85746029E+00,4.41437026E-03,-2.21481404E-06,5.23490188E-10,-4.72084164E-14,-4.87591660E+04,2.27163806E+00]
h2o_coeff_h=[3.03399249E+00,2.17691804E-03,-1.64072518E-07,-9.70419870E-11,1.68200992E-14,-3.00042971E+04,4.96677010E+00]

def f(T,alk):
	h_co2_p=h(T,co2_coeff_h)
	h_n2_p=h(T,n2_coeff_h)
	h_h2o_p=h(T,h2o_coeff_h)

	if alk==methane:
		x=1
		y=4
		H_products=x*h_co2_p+(y/2)*h_h2o_p+7.52*h_n2_p
	if alk==ethane:
		x=2
		y=6
		H_products=x*h_co2_p+(y/2)*h_h2o_p+13.16*h_n2_p
	if alk==propane:
		x=3
		y=8
		H_products=x*h_co2_p+(y/2)*h_h2o_p+18.8*h_n2_p

	Tstd=298.15
	h_ch4_r=h(Tstd,ch4_coeff_l)
	h_c2h6_r=h(Tstd,c2h6_coeff_l)
	h_c3h8_r=h(Tstd,c3h8_coeff_l)
	h_o2_r=h(Tstd,o2_coeff_l)
	h_n2_r=h(Tstd,n2_coeff_l)

	if alk==methane:
		H_reactant=h_ch4_r+2*h_o2_r+7.52*h_n2_r
	if alk==ethane:
		H_reactant=h_c2h6_r+3.5*h_o2_r+13.16*h_n2_r
	if alk==propane:
		H_reactant=h_c3h8_r+5*h_o2_r+18.8*h_n2_r

	return H_products-H_reactant

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

T_guess=1500
alpha=0.2
tol=1e-3
ct=1
for alk in a:
	while(abs(f(T_guess,alk))>tol):
		T_guess=T_guess-alpha*(f(T_guess,alk)/fprime(T_guess,alk))
		ct=ct+1
		plt.plot(alk,T_guess,'*',color='blue')
	print(f(T_guess,alk))
	print(T_guess)
plt.xlabel('No of C atom')
plt.ylabel('Temperature')
plt.show()

Plot:

Following are the AFT for Methane,Ethane and Propane

 The AFt for methane is lowest as it has lower carbon atoms as compared to propane.