Calculating auto-ignition time for methane combustion under various initial condition using Python and Cantera

INTRODUCTION:

Auto-ignition occurs when a mixture of gases or vapors ignites spontaneously in absence of any external ignition source such as a flame or a spark. The temperature at which auto-ignition occurs is called the auto-ignition temperature. It is the lowest temperature in a system where the rate of heat evolved from the gases or vapors increases beyond the rate of heat loss to the surroundings, resulting in the ignition. The auto-ignition temperature is not an intrinsic property of the gases or vapors and is affected by the pressure, vessel shape and volume, surface activity, contaminants, flow rate, reaction rate, gravity, reactant concentration, etc.

A Cantera reactor network represents the simplest form of a chemically reacting system. It corresponds to an extensive thermodynamic control volume `V`, in which all state variables are homogeneously distributed. Transient state changes due to chemical reactions are possible, however, thermodynamic (and not chemical) equilibrium is assumed to be present throughout the reactor at all instants of time. The ReactorNet() object is used to model the reactor class. 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 behavior of the reactor. Even for a single reaction, a reactor network must be created in Cantera to visualize the transient behavior. For this problem, we are defining a reactor network for an ideal gas reactor whose state variables are:

• `m [kg]`: the mass of the reactor's content.
• `V [m^3]`: the reactor volume.
• `T [K]`: the temperature (for describing the energy of the system).
• `Y_k`: the mass fraction of each species.

The governing equations for the reactor are modeled as:

1. Mass Conservation

`{dm}/{dt} = sum_{i n} dot{m}_{i n} - sum_{out} dot{m}_{out} + dot{m}_{wall}`

2. Species Conservation

`dot{m}_{k, g en} = V dot{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 + dot{m}_{k, g en}`

3. Energy Conservation

`U = m sum_{k} Y_k u_k(T)`

`{dU}/{dt} = u {dm}/{dt} + mc_v{dT}/{dt} + m sum_{k}u_k{dY_k}/{dt}`

IMPLEMENTATION IN PYTHON USING CANTERA MODULE:

The auto-ignition study of methane is performed as per the following steps:

• Defining the initial temperature and pressure condition and the mass fractions of each species `(CH_4, O_2` and `N_2)`.
• Defining the solution time (`10` `sec` for this case).
• Initializing the Solution object using the initial conditions.
• Defining the type of reactor (Ideal gas reactor for this case).
• Creating the reactor network from the single reactor.
• 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 to the output.
• Post-processing the output data.

The auto-ignition temperature for this case is defined as:

`T_text{auto-ignition} = T_text{init} + 400` `K`

The time at which this temperature is reached is called the ignition delay.

Part I

The auto-ignition of methane is studied under different initial conditions:

1. A constant temperature of `1250` `K` and varying pressure from `1` to `5` `atm`.

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

Part II

The rate of change of mole fractions of `H_2O, O_2` and `OH` is studied for `2` cases:

1. Temperature: `1000` `K` and Pressure: `5` `b ar`.

2. Temperature: `500` `K` and Pressure: `5` `b ar`.

A function is created to output the state of the reactor network for the specified initial conditions and time. This function is called through the main code which makes it easier to implement it for all the different cases. The entire code and the results are shown below:

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


Program by: Shrey Shah
Date: 08-Oct-19
"""

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

# 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 = np.linspace(950,1450,11)
P = 101325*np.linspace(1,5,9)
species_dict = {'CH4':0.095, 'O2':0.19, 'N2':0.715}
t = np.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 = plt.subplots()
for i in range(0,len(P)):
	plt.plot(t, T_ai_P[i], linewidth=2, label=str(P[i]/101325)+' atm')
plt.xlabel('Time [sec]', fontsize=14)
plt.ylabel('Temperature [K]', fontsize=14)
plt.title('Reaction Temperature with varying Pressure', fontsize=14, fontweight='bold')
plt.grid('both', linestyle='--', linewidth=0.8)
plt.legend(fontsize=12)

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

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

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

fig5, ax5 = plt.subplots(3)
ax5[0].plot(t, mole_frac1[0], linewidth=2)
plt.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)

fig6, ax6 = plt.subplots(3)
plt.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)

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

plt.show()

Part I:

From the above results, we can see that as we increase the initial temperature or pressure, the ignition delay reduces asymptotically to `0` `ms`. The final temperature of the reaction is also affected by the change in the initial conditions. The rate of change of ignition delay is higher with varying temperatures than varying pressure. A low temperature, say `950` `K`, will take almost `400` `ms` to ignite which is considerably larger compared to other ignition delays.

Part II:

Here, we can see that for an initial temperature of `1000` `K`, the auto-ignition occurs with some amount of ignition delay which increases the mole fractions of the products `H_2O` and `OH` and decreases the mole fraction of the reactant `O_2`. We can also see that the spike for the change in mole fraction is consistent for all the three species which indicates a spontaneous ignition of methane. However, for an initial temperature of `500` `K`, the mole fractions of the products `H_2O` and `OH` are nearly zero and the mole fraction of the reactant `O_2` does not budge from the initial mole fraction of `0.19`. This indicates that no reaction/ignition has taken place. This is because a temperature of `500` `K` is too low to cause any ignition of methane in the absence of any external spark or flame. This can also be seen by comparing the temperature plots for both conditions: