Week 7 - Auto ignition using Cantera

Aim: Calculation of auto-ignition time of methane using Python and Cantera

Objective:

• Plotting auto-ignition for constant temperature and varing pressure and for constant pressure and varying temperature.
• Observing and explaining the behaviour of change of molar concentration of given species.

Theory:

The auto-ignition temperature or kindling point of a substance is the lowest temperature in which it spontaneously ignites in normal atmosphere without an external source of ignition, such as flame or spark. This temperature is required to supply the activation energy needed for combustion. This temperature reduces as the pressure or oxygen concentration increases.

A cantera network represents the simplest form of a chemically reacting system. It corresponds to an extensive thermodynamic control volume `V`, where all state variables are homogeneously distributed. Transient state changes due to chemical reactions are possible, however, thermal equilibrium is assumed to be present throughout the reactor at all instants of time. The ReactorNet() object is used to model the reactor to the specified time. The advance(time) function is used to advance the reactor to the specified time. This thermodynamic state of the reactor can then be extracted at each time-step to study the transient behaviour of the reactor. Even for a single reaction, a reactor network must be created in Cantera to visualize the transient behaviour. For this problem, we are defining a reactor network for an ideal gas reactor whose state variables are:

• `m` : the mass of reactor's content (in `kg`)
• `V` : the reactor volume (in `m^3`)
• `T` : temperature (in `K`)
• `Y_k` : mass fraction of each species

We apply the following conservation principles to the reactor model:

1. Mass Conservation

`(dm)/(dt) = sum_(i n)dot(m_{i n}) - sum_(out)dot(m_{out}) + dot(m_{wall})`

2. Energy Conservation

`U = msum_(k)Y_ku_k(T)`

`(dU)/(dt) = u(dm)/(dt) + mc_v(dT)/(dt) + msum_(k)u_k(dY_k)/(dt)`

3. Species Conservation

`dot(m_{k,g e n}) = Vdot(omega_k)W_k + dot(m_{wall})`

`(d(mY_k))/(dt) = sum_(i n)dot(m_{i n})Y_{k,i n} - sum_(out)dot(m_{out})Y_{k,out} + dot(m_{k,g e n})`

Code:

The study is done by implemeting the following steps:

• Defining the initial temperature and pressure of the species and the solution time.
• Initializing Solution object for the initial conditions.
• Defining and creating the reactor network and initializing the thermodynamic state array of the reactor for output.
• Advancing the reactor network through the desired time and adding the thermodynamic state at each time step and post processing the output data.

The auto-ignition, here, is defined as:

`T_{a u t o - i g n i t i o n} = T_{i n i t} + 400K`

The time required to reach this temperature is called ignition delay.

Part-I

The auto-ignition of methane is done as:

1. A constant temperature of 1250K and varying pressure from 1 to 5 atm.
2. A constant pressure 5 atm and varying temperature from 950to 1450K.

Part-II

The rate of change of mole fractions of `H_2O, O_2` and `OH`:

1. Temperature: 1000K and Pressure: 5bar
2. Temperature: 500K and Pressure: 5 bar

A function is used to determine the output state of the reactor network for specified initial condition. This function is called in the code which makes it easy to gather ouptuts for different cases:

"""
Program to study the auto-ignition of methane under various conditions.

"""

import cantera as ct
from matplotlib.pyplot import *
from numpy import *

# Defining the function to output the thermodynamic states at each time-step
def auto_ignition(T, P, species_dict, t):
	gas = ct.Solution('gri30.xml')
	gas.TPX = T, P, species_dict

	r = ct.IdealGasReactor(gas)

	sim = ct.ReactorNet([r])

	states = ct.SolutionArray(gas, extra=['time_in_ms'])

	for tp in t:
		sim.advance(tp)
		states.append(r.thermo.state, time_in_ms = tp*1e3)

	return states

# Defining the various initial conditions
T = linspace(950,1450,11)
P = 101325*linspace(1,5,9)
species_dict = {'CH4':0.095, 'O2':0.19, 'N2':0.715}
t = linspace(0,10,10001)

# Initializing the output parameters
T_ai_T = [[0 for time in t] for temp in T]
T_ai_P = [[0 for time in t] for pres in P]
T_mole = [1000, 500]
T_comp = [[0 for time in t] for i in range(0,2)]
time_ai_T = []
time_ai_P = []
mole_frac1 = [[0 for time in t] for i in range(0,3)]
mole_frac2 = [[0 for time in t] for i in range(0,3)]

# Auto-ignition for constant temperature (1250 K)
for i in range(0,len(P)):
	states = auto_ignition(T[6], P[i], species_dict, t)
	T_ai_P[i] = states.T
	for j in range(0,len(states.T)):
		if states.T[j] >= states.T[0] + 400:
			time_ai_P.append(states.time_in_ms[j])
			break

# Auto-ignition for constant pressure (5 atm)
for i in range(0,len(T)):
	states = auto_ignition(T[i], P[-1], species_dict, t)
	T_ai_T[i] = states.T
	for j in range(0,len(states.T)):
		if states.T[j] >= states.T[0] + 400:
			time_ai_T.append(states.time_in_ms[j])
			break

# Mole Fractions for temperature of 1000 K
states = auto_ignition(T_mole[0], 5*100000, species_dict, t)
mole_frac1[0] = states.X[:, states.species_index('H2O')]
mole_frac1[1] = states.X[:, states.species_index('O2')]
mole_frac1[2] = states.X[:, states.species_index('OH')]
T_comp[0] = states.T

# Mole Fractions for temperature of 500 K
states = auto_ignition(T_mole[1], 5*100000, species_dict, t)
mole_frac2[0] = states.X[:, states.species_index('H2O')]
mole_frac2[1] = states.X[:, states.species_index('O2')]
mole_frac2[2] = states.X[:, states.species_index('OH')]
T_comp[1] = states.T

# Plotting all the results
fig1, ax1 = subplots()
for i in range(0,len(P)):
	plot(t, T_ai_P[i], linewidth=2, label=str(P[i]/101325)+' atm')
xlabel('Time [sec]', fontsize=14)
ylabel('Temperature [K]', fontsize=14)
title('Reaction Temperature with varying Pressure', fontsize=14, fontweight='bold')
grid('both', linestyle='--', linewidth=0.8)
legend(fontsize=12)
savefig('Reaction_Temperature_with_varying_Pressure.png')
show()

fig2, ax2 = subplots()
plot(P/101325, time_ai_P, 'o-', linewidth=2, markersize=10, color='r')
xlabel('Pressure [atm]', fontsize=14)
ylabel('Ignition Delay [ms]', fontsize=14)
title('Ignition Delay with varying Pressure', fontsize=14, fontweight='bold')
grid('both', linestyle='--', linewidth=0.8)
savefig('Ignition_Delay_with_varying_Pressure.png')
show()

fig3, ax3 = subplots()
for i in range(0,len(T)):
	plot(t, T_ai_T[i], linewidth=2, label=str(int(T[i]))+' K')
xlabel('Time [sec]', fontsize=14)
ylabel('Initial Temperature [K]', fontsize=14)
title('Reaction Temperature with varying Initial Temperature', fontsize=14, fontweight='bold')
grid('both', linestyle='--', linewidth=0.8)
legend(fontsize=12)
savefig('Reaction_Temperature_with_varying_Initial_temperature.png')
show()

fig4, ax4 = subplots()
plot(T, time_ai_T, 'o-', linewidth=2, markersize=10, color='r')
xlabel('Initial Temperature [K]', fontsize=14)
ylabel('Ignition Delay [ms]', fontsize=14)
title('Ignition Delay with varying Initial Temperature', fontsize=14, fontweight='bold')
grid('both', linestyle='--', linewidth=0.8)
savefig('Ignition_Delay_with_varying_Initial_Temperature.png')
show()

fig5, ax5 = subplots(3)
ax5[0].plot(t, mole_frac1[0], linewidth=2)
suptitle('Mole Fractions for Reaction at 1000 K and 5 bar', fontsize=14, fontweight='bold')
ax5[0].grid('both', linestyle='--', linewidth=0.8)
ax5[0].set_xticklabels([])
ax5[1].plot(t, mole_frac1[1], linewidth=2)
ax5[1].grid('both', linestyle='--', linewidth=0.8)
ax5[1].set_xticklabels([])
ax5[2].plot(t, mole_frac1[2], linewidth=2)
ax5[0].set_ylabel('$X_{H_2O}$', fontsize=14)
ax5[1].set_ylabel('$X_{O_2}$', fontsize=14)
ax5[2].set_ylabel('$X_{OH}$', fontsize=14)
ax5[2].set_xlabel('Time [sec]', fontsize=14)
ax5[2].grid('both', linestyle='--', linewidth=0.8)
savefig('Mole_Fractions_for_Reaction_at_1000K_and_5bar.png')
show()

fig6, ax6 = subplots(3)
suptitle('Mole Fractions for Reaction at 500 K and 5 bar', fontsize=14, fontweight='bold')
ax6[0].plot(t, mole_frac2[0], linewidth=2)
ax6[0].grid('both', linestyle='--', linewidth=0.8)
ax6[0].set_xticklabels([])
ax6[1].plot(t, mole_frac2[1], linewidth=2)
ax6[1].grid('both', linestyle='--', linewidth=0.8)
ax6[1].set_xticklabels([])
ax6[2].plot(t, mole_frac2[2], linewidth=2)
ax6[0].set_ylabel('$X_{H_2O}$', fontsize=14)
ax6[1].set_ylabel('$X_{O_2}$', fontsize=14)
ax6[2].set_ylabel('$X_{OH}$', fontsize=14)
ax6[2].set_xlabel('Time [sec]', fontsize=14)
ax6[2].grid('both', linestyle='--', linewidth=0.8)
savefig('Mole_Fractions_for_Reaction_at_500K_and_5bar.png')
show()

fig7, ax7 = subplots()
for i in range(0,2):
	plot(t, T_comp[i], linewidth=2, label=str(T_mole[i])+' K')
xlabel('Time [sec]', fontsize=14)
ylabel('Temperature [K]', fontsize=14)
title('Reaction Temperature with varying Initial Temperature', fontsize=14, fontweight='bold')
grid('both', linestyle='--', linewidth=0.8)
legend(fontsize=12)
savefig('Reaction_Temperature_with_varying_Initial_Temperature.png')
show()

Output:

Part-I

Part-II

Result&Conclusion:

• We observe that, with increase in initial temperature or pressure, the ignition delay approaches zero (and then beccomes an asymptote).
• The final temperature also increases with increase in initial temperature. Also, rate of decrease of ignition delay is greater in varying temperature as compared to varying pressure.
• At initial condition of 1000K and 5bar, we can observe after certain duration of time, there is a spike in formation of `H_2O` and `OH` and reduction of `O_2`. This signifies spontaneous ignition of methane. But, for 500K and 5bar condition the molar fraction of `O_2` remains constant, i.e. combustion did not take place as the temperature, 500K is too low for ignition to take place without a flame.