Analysis of Auto-ignition temperature at different operating conditions using Cantera.

Title:- Autoignition using Cantera

Problem Statement:-

Using the inbuilt reactor network code in Cantera, the following objectives to be attained.

1. Variation of Auto Ignition time of Methane with a constant temperature of 1250 K and pressure to be varied from 1 to 5 atm.

2. Variation of Auto Ignition time of Methane with a constant pressure of 5 atm and temperature to be varied from 950K to 1450K.

3. Rate of change of molar concentrations of ' H2O, O2, OH' w.r.t time for 10 sec simulation time.

4. The behavior of the following graphs:-

• When Temperature is 1000K and pressure 5 bar

• When Temperature is 500K and pressure 5 bar

Analysis:-

Introduction:-

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. The system is generally unsteady, i.e., all states are functions of time. In particular, transient state changes due to chemical reactions are possible. However, thermodynamic (but not chemical) equilibrium is assumed to be present throughout the reactor at all instants of time.

Reactors can interact with the surrounding environment in multiple ways:

• Expansion/compression work.
• Heat transfer.
• Mass transfer.
• Surface interaction.

All of these interactions do not have to be constant but can vary as a function of time or state. Interactions of the reactor with the environment are defined in one or more walls, inlets, and outlets.

In addition to single reactors, Cantera is also able to interconnect reactors into a Reactor Network. Each reactor in a network may be connected so that the contents of one reactor flow into another. Reactors may also be in contact with one another or the environment via walls that conduct heat or move to do work.

As the above objectives to be fulfilled by a single reactor, the use of a single reactor to simulate the system will be enough. Here ideal gas reactor is used, so the equations to be solved are:

In the case of the Ideal Gas Reactor Model, the reactor temperature 'T' used instead of the total internal energy 'U' as a state variable. For an ideal gas, we can rewrite the total internal energy in terms of the mass fractions and temperature. Following energy equations will be solved in addition to mass and species concentration governing equations.

`U = m sum_k Y_k u_k(T)`

`(dU)/dt = u*((dm)/(dt)) + m*C_v*((dT)/(dt)) + m*sum_k u_k * (dY_k)/(dt`

`mc_v*(dT)/( dt) =−p (dV)/(dt) − Q+ m_(i) sum_ m (h_i − sum_k u_k Y_ki )- (pV)/m *sum_(out) m_(out) - sum_k m_k,g*u_k`

While this form of the energy equation is somewhat more complicated, it significantly reduces the cost of evaluating the system Jacobian, since the derivatives of the species equations are taken at constant temperature instead of constant internal energy.

Steps in Cantera to solve reactor network:

In Cantera, the following steps are typically necessary to investigate a reactor network:

1. Define Solutionobjects for the fluids to be flowing through your reactor network.
2. Define the reactor type(s) and reservoir(s) that describe your system. Chose Ideal Gas (Constant Pressure) Reactor(s) if you only consider ideal gas phases.
3. Optional: Set up the boundary conditions and flow devices between reactors or reservoirs.
4. Define a reactor network which contains all the reactors previously created.
5. Advance the simulation in time, typically in a for- or while-loop. Note that only the current state is stored in Cantera by default. If you want to observe the transient states, you manually have to keep track of them.
6. Analyze the data.

Case I: Autoignition time when pressure varied from 1bar to 5bar.

The code for solving the above objective is:

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

gas = Solution('gri30.xml')
pressure  = linspace(1,5,9)
T_int = 1250	
time_start =0.0
time_end = 10.0
time_step =1e-3
n = int((time_end - time_start)/time_step)
t_autoI_P = []
states1 = []


for i in pressure:
	gas.TPX = T_int, i*101325, {'CH4':1,'O2':2,'N2':7.52}

	T_autoI = 0
	t_autoI = 0
	time = 0
	r = IdealGasReactor(gas)
	sim = ReactorNet([r])
	states = SolutionArray(gas, extra=['time_in_ms','t'])
	print('{:10s} {:10s} {:10s}'.format('t [s]','T [K]','P [Pa]'))

	for n in range(n):
		time += 1e-3
		sim.advance(time)
		states.append(r.thermo.state,time_in_ms= time*1e3, t=time)
		print('{:10.3e} {:10.3f} {:10.3f} '.format(sim.time, r.T, r.thermo.P))
	states1.append(states.T)

	for i in range(n):
		T_error = states.T[i] - T_int -400
		if T_error > 0:
			T_autoI = states.T[i]
			t_autoI = states.time_in_ms[i]
			break
	t_autoI_P.append(t_autoI)


states1 = np.array(states1)
for k in range(len(pressure)):
	figure(1)
	plot(states1[k])
	xlabel('Number of timesteps')
	ylabel('Temperature')
	title('Variation of Auto-ignition Temperature with inlet Pressure')
	legend(['1 Bar','1.5 Bar', '2 Bar', '2.5 Bar','3 Bar','3.5 Bar', '4 Bar', '4.5 Bar', '5 Bar'])

figure(2)
plot(pressure,t_autoI_P,'r-*')
xlabel('Pressure [bar]')
ylabel('Auto-ignition Time [ms]')
title('Ignition time delay Vs Pressure')
show()

The results of the above codes are:

Case II: Autoignition time when temperature varied from 950 K to 1450 K.

The code for solving the above objective is:

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

gas = Solution('gri30.xml')
Temperature  = linspace(950,1450,11)
time_start =0.0
time_end = 10.0
time_step =1e-3
n = int((time_end - time_start)/time_step)
t_autoI_T = []
states1 = []


for i in Temperature:

	gas.TPX = i, 5*one_atm, {'CH4':1,'O2':2,'N2':7.52}

	T_autoI = 0
	t_autoI = 0
	time = 0
	r = IdealGasReactor(gas)
	sim = ReactorNet([r])
	states = SolutionArray(gas, extra=['time_in_ms','t'])
	print('{:10s} {:10s} {:10s}'.format('t [s]','T [K]','P [Pa]'))

	for j in range(n):
		time += 1e-3
		sim.advance(time)
		states.append(r.thermo.state,time_in_ms= time*1e3, t=time)
		print('{:10.3e} {:10.3f} {:10.3f} '.format(sim.time, r.T, r.thermo.P))
	states1.append(states.T)

	for j in range(n):
		T_error = states.T[j] - i -400
		if T_error > 0:
			T_autoI = states.T[j]
			t_autoI = states.time_in_ms[j]
			break
	t_autoI_T.append(t_autoI)

states1 = array(states1)
for k in range(len(Temperature)):
	figure(1)
	plot(states1[k])
	xlabel('Number of timesteps')
	ylabel('Temperature')
	title('Variation of Auto-ignition Temperature with inlet Temperature')
	legend(['950 K','1000 K', '1050 K', '1100 K','1150 K','1200 K', '1250 K', '1300 K', '1350 K', '1400 K', '1450 K'])

figure(2)
plot(Temperature,t_autoI_T,'b-*')
xlabel('Temperature [K]')
ylabel('Auto-ignition Time [ms]')
title('Ignition time delay Vs Initial Temperature')
show()

The results of the above codes are:

Case III: Rate of change of molar concentrations of ' H2O, O2, OH' w.r.t time for 10 sec simulation time.

The code for the finding out molar concentration is:

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

gas = Solution('gri30.xml')
Temperature  = linspace(950,1450,11)
time_start =0.0
time_end = 10.0
time_step =1e-3
n = int((time_end - time_start)/time_step)
t_autoI_T = []
states1 = []


for i in Temperature:

	gas.TPX = i, 5*one_atm, {'CH4':1,'O2':2,'N2':7.52}

	T_autoI = 0
	t_autoI = 0
	time = 0
	r = IdealGasReactor(gas)
	sim = ReactorNet([r])
	states = SolutionArray(gas, extra=['time_in_ms','t'])
	print('{:10s} {:10s} {:10s}'.format('t [s]','T [K]','P [Pa]'))

	for j in range(n):
		time += 1e-3
		sim.advance(time)
		states.append(r.thermo.state,time_in_ms= time*1e3, t=time)
		print('{:10.3e} {:10.3f} {:10.3f} '.format(sim.time, r.T, r.thermo.P))
	states1.append(states.T)

	for j in range(n):
		T_error = states.T[j] - i -400
		if T_error > 0:
			T_autoI = states.T[j]
			t_autoI = states.time_in_ms[j]
			break
	t_autoI_T.append(t_autoI)

figure(1)

subplot(3, 1, 1)
plot(states.t, states.X[:,gas.species_index('H2O')])
xlabel('Time (ms)')
ylabel('H2O Mole Fraction')
subplot(3, 1, 2)
plot(states.t, states.X[:,gas.species_index('O2')])
xlabel('Time (ms)')
ylabel('O2 Mole Fraction')
subplot(3, 1, 3)
plot(states.t, states.X[:,gas.species_index('OH')])
xlabel('Time (ms)')
ylabel('OH Mole Fraction')
title('Molar concentrations of H2O, O2, OH for 10 sec simulation time')
tight_layout()

show()

The concentration after 10 sec is:

a) When Temperature varied from 950K to 1450K.

b) When Pressure varied from 1 bar to 5 bar.

Case IV: Comparison of temperature and auto-ignition time at 500K & 1000K.

from sys import*
from numpy import*
from cantera import*
from matplotlib.pyplot import*

gas = Solution('gri30.xml')

Temp = [500,1000]
states1 = []

for i in Temp:
	gas.TPX = i, 506625, {'CH4':1, 'O2':2,'N2':7.52}
	r = IdealGasReactor(gas)
	sim = ReactorNet([r])
	time =0.0
	states = SolutionArray(gas)

	for n in range(1000):
		time += 1e-3
		sim.advance(time)
		states.append(r.thermo.state)

	states1 = np.array(states.T)
	for k in range(1000):
		figure(1)
		plot(states1)
		xlabel('Number of timesteps')
		ylabel('Temperature')
		title('Temperature variation at 500 K & 1000 K at 5 bar')
show()

Results:

1. The initial temperature of 500K is not enough to auto-ignite the mixture even if the pressure is 5 bar.

2. The temperature 1000k will start igniting after 0.2 sec and it will reach a maximum temperature of around 2900K.

3. The temperature 500K will not ignite mixture even if pressure is increased to any magnitude.