Simulating auto ignition using Cantera

Objective:

To find the autoignition time of methane under varying conditions of pressure & temperature and analyze the evolution of molar concentrations of reactants with combustion progress.

For this we will consider:

1. Auto ignition time with
   1. Constant temperature of 1250 K and pressure variation from 1 to 5 atmospheres.
   2. Constant pressure of 5 atmospheres and temperature variation from 950 K to 1450 K.
2. Molar concentrations at
   1. 1000 K and 5 atm.
   2. 500 K and 5 atm.

Introduction:

When a mixture of gases combust spontaneously without the presence of an ignition source like a flame or spark the phenomenon is called autoignition.

Autoignition depends mainly on the mixture, temperature, and pressure. The lowest temperature at a given pressure at which the mixture combusts spontaneously is called autoignition temperature or kindling point. Autoignition temperature is an extrinsic property and can depend on flow rate, vessel geometry, pressure, reaction rate, concentration gradients, etc.

Autoignition temperature is of prime importance in combustion engines, gas turbines, rocket propulsion, etc

Approach:

Combustion will be simulated in python using reactor network available in cantera package. A cantera reactor represents the simplest form of a chemically reacting system and has the following characteristics:

• It is a thermodynamic control volume V, in which all state variables are homogeneously distributed.
• The system is generally unsteady -- that is, all states are functions of time so that transient state changes due to chemical reactions are possible.
• Thermodynamic (but not chemical) equilibrium is assumed to be present throughout the reactor at all instants of time.

The state of this reactor will be changed through a series of many equilibrium states each separated by a small time step so that the process (combustion) can be closely approximated as a quasi-static process.

CODE

• For simulating combustion, a function combust() is defined which takes initial temperature, pressure, and species concentration as user inputs. Optionally, the time step size can be changed & also whether to simulate only till ignition.
• Inside combust(), a reactor object is made which mimics an ideal gas combustion mixture at a given state.
• Next, the reactor network object is made from the reactor object which will simulate the system.
• A solution array is made using cantera so that the state at each time step can be saved.
• The time loop is then run. Each iteration changes time by a small amount then, the reactor state is computed at that time using state variables at the previous time step.
• As soon as reactor temperature increases above 400 K of initial value auto-ignition has occurred and the temperature is saved.
• The gas object, states (Solution array), and ignition_delay are returned by the function.
• The results are plotted after providing appropriate arguments to combust().

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

def combust(T: float, P: float, spcies_dict: dict, reactor: str, dt: float= 0.0005, till_ignition: bool=True):
    gas = ct.Solution('gri30.yaml')
    gas.TPX = T, P, spcies_dict
    # Selecting reactor
    if reactor == 'constant temperature':
        r = ct.IdealGasReactor(gas)
    
    elif reactor == 'constant pressure':
        r = ct.IdealGasConstPressureReactor(gas)
    
    else:
        raise TypeError(f"Improper reactor type '{reactor}'")

    # Making reactor network
    sim = ct.ReactorNet([r])
    # Making a solution array  object
    states = ct.SolutionArray(gas, extra=['time_in_ms', 'time_in_sec'])

    # Simulates for  10 sec = dt*n_steps
    simulation_time = 10
    # time step is dt
    n_steps = int(simulation_time/dt)
    ignition_delay = 0.0
    
    # A flag so that auto-ignition temperature is computed only 
    # once, the first time temperatue spikes by 400 [K]
    ignited = False
    
    # Time loop
    time = 0.0
    for n in range(n_steps):
        time += dt
        sim.advance(time)
        states.append(r.thermo.state, time_in_sec= time, time_in_ms= time*1e3)
        
        if ignited == False:
            if states.T[n] >= (T + 400):
                ignition_delay = time
                ignited = True
                if till_ignition == True:
                    break

    return gas, states, ignition_delay

# %%
# PART I-A
fig1, axs = plt.subplots(2,1, figsize=(10, 10), layout='tight')

igd_lst = []
P = np.linspace(1,5,9)
for t in P:
    _, states, igd = combust(1250, t*ct.one_atm, {'CH4':1, 'O2':2, 'N2': 7.52}, 'constant temperature', dt= 2e-4)
    axs[0].plot(states.time_in_ms, states.T, '.-', label= f'{t: .3} atm, I.D. = {igd*1000: .5} ms')
    igd_lst.append(igd)

axs[0].set_xlabel('Time [ms]$\longrightarrow$')
axs[0].set_ylabel("Flame Temperature [K]$\longrightarrow$")
axs[0].grid(linestyle='-.')
axs[0].legend()

axs[1].plot(np.array(igd_lst)*1000, P, '.-')
axs[1].set_xlabel('Auto-ignition time [ms]$\longrightarrow$')
axs[1].set_ylabel("Pressure [atm]$\longrightarrow$")
axs[1].grid(linestyle='-.')

fig1.suptitle(f"Auto-ignition of methane at 1250 Kn and pressure variation")
# %%
# PART I-B
fig2, axs = plt.subplots(2,1, figsize=(10, 10), layout='tight')

igd_lst = []
T = np.linspace(950, 1450, 9)
for t in T:
    _, states, igd = combust(t, 5*ct.one_atm, {'CH4':1, 'O2':2, 'N2': 7.52}, 'constant pressure')
    axs[0].plot(states.time_in_ms, states.T, '.-', label= f'{t} K, I.D. = {igd*1000: .5} ms')
    igd_lst.append(igd)

axs[0].set_xlabel('Time [ms]$\longrightarrow$')
axs[0].set_ylabel("Flame Temperature [K]$\longrightarrow$")
axs[0].grid(linestyle='-.')
axs[0].legend()

axs[1].plot(np.array(igd_lst)*1000, T, '.-')
axs[1].set_xlabel('Auto-ignition time [ms]$\longrightarrow$')
axs[1].set_ylabel("Initial temperature [K]$\longrightarrow$")
axs[1].grid(linestyle='-.')

fig2.suptitle(f"Auto-ignition of methane at 5 atmn and initial temperature variation")
#%%
# PART II-A & B
gasA, statesA, igdA = combust(1000, 5*ct.one_atm, {'CH4':1, 'O2':2, 'N2': 7.52}, 'constant temperature', till_ignition= False)
gasB, statesB, igdB = combust(500, 5*ct.one_atm, {'CH4':1, 'O2':2, 'N2': 7.52}, 'constant temperature', dt= 0.01, till_ignition= False)
#%%
fig3, axs = plt.subplots(3,1, figsize=(10, 10), layout='tight')

axs[0].plot(statesA.time_in_sec, statesA.X[:,gasA.species_index('H2O')], '.-', color='orange', label= f'1000 K, I.D. = {igdA*1000: .5} ms')
axs[0].plot(statesB.time_in_sec, statesB.X[:,gasB.species_index('H2O')], '.-', label= f'500 K, I.D. = {igdB*1000: .5} ms')
axs[0].set_ylabel("$[H_2O]$")
axs[0].grid(linestyle='-.')
axs[0].legend()

axs[1].plot(statesA.time_in_sec, statesA.X[:,gasA.species_index('O2')], '.-', color='orange', label= f'1000 K, I.D. = {igdA*1000: .5} ms')
axs[1].plot(statesB.time_in_sec, statesB.X[:,gasB.species_index('O2')], '.-', label= f'500 K, I.D. = {igdB*1000: .5} ms')
axs[1].set_ylabel("$[O_2]$")
axs[1].grid(linestyle='-.')
axs[1].legend()

axs[2].plot(statesA.time_in_sec, statesA.X[:,gasA.species_index('Oh')], '.-', color='orange', label= f'1000 K, I.D. = {igdA*1000: .5} ms')
axs[2].plot(statesB.time_in_sec, statesB.X[:,gasB.species_index('Oh')], '.-', label= f'500 K, I.D. = {igdB*1000: .5} ms')
axs[2].set_ylabel("$[OH]$")
axs[2].grid(linestyle='-.')
axs[2].legend()

fig3.supxlabel('Time [sec]$\longrightarrow$')
fig3.supylabel("Molar concentration$\longrightarrow$")

plt.show()

Results & Inferences:

PART I-A

• The first plot shows the evolution of temperature only till auto-ignition is deemed to occur i.e. until the initial temperature increases by 400K.
• In the first plot, despite a small time step, ignition causes the temperature to spike so rapidly that the temperature rise curve is not captured adequately. For this, a very small time step will be needed which puts a load on computational resources.
• Clearly, auto-ignition time reduces with an increase in pressure.
• An increase in pressure increases the collision frequency leading to an increase in reaction rate.

PART I- B

• Clearly, auto-ignition time reduces with increasing initial temperature.
• This is because as per the Arrhenius equation, reaction rate varies exponentially with temperature.
• The trend line shows that auto-ignition time reduces asymptotically with temperature.

PART II-A & II-B

• At 1000K & 5 atm, `OH` & `H_2O` concentration rapidly increases while molar concentration of `O_2` reduces.
• This indicates the spontaneous combustion of methane.
• However, at 500K & 5 atm molar concentrations do not change. This is because the pressure & temperature are too low to auto-ignite methane.