Modifying custom ODE example to model constant volume but pressure off

[George Eappachan]

Hi Cantera users,

I am working on modifying the custom ODE example available at https://cantera.org/examples/python/reactors/custom.py.html in Python. The example provided models a constant pressure reactor. I am modifying the custom ODE example to model a constant volume reactor.

As I will have to work with modifying the governing equations (in the future) I am not using the IdealGasReactor. Just to make sure I am doing right, I compared the results with the ODE integration approach and the IdealGasReactor (the MWE to do this is included). However, the pressure values are underestimated by the ODE approach (comparison given below), due to some error on my side probably, but I cannot figure out what. Shouldn't the pressure be calculated using ideal gas law if the energy and species equations are correct, am I missing something?

The ignition point and OH mole fraction profile appears to match. What could I be possibly doing wrong that the pressure is underestimated and temperature is slightly off?

Any help in this regard would be highly appreciated.

Thank you.

Regards,
George

MWE

"""
Requires: cantera >= 2.5.0, matplotlib >= 2.0
"""


import cantera as ct
import numpy as np
import scipy.integrate


class ReactorOde:
    def __init__(self, gas):
        # Parameters of the ODE system and auxiliary data are stored in the
        # ReactorOde object.
        self.gas = gas

    def __call__(self, t, y):
        """the ODE function, y' = f(t,y) """

        # State vector is [T, Y_1, Y_2, ... Y_K]
        self.gas.set_unnormalized_mass_fractions(y[1:])
        self.gas.TP = y[0], self.gas.P
        rho = self.gas.density

        wdot = self.gas.net_production_rates
        dTdt = - (np.dot(self.gas.partial_molar_enthalpies, wdot) /
                (rho * self.gas.cp))
        dYdt = wdot * self.gas.molecular_weights / rho

        return np.hstack((dTdt, dYdt))



P = 1e5
T = 990

gas_reactor = ct.Solution('gri30.yaml')
gas_ode     = ct.Solution('gri30.yaml')

gas_reactor.TPX = T, P, 'H2:2,O2:1,N2:4'
gas_ode.TPX     = T, P, 'H2:2,O2:1,N2:4'

r = ct.IdealGasReactor(gas_reactor)

sim = ct.ReactorNet([r])
# sim.verbose = True

# limit advance when temperature difference is exceeded
delta_T_max = 20.
r.set_advance_limit('temperature', delta_T_max)

dt = 1.e-5
t_end = 100 * dt
states_reactor = ct.SolutionArray(gas_reactor, extra=['t'])

# print('{:10s} {:10s} {:10s} {:14s}'.format(
#   't [s]', 'T [K]', 'P [Pa]', 'u [J/kg]'))
while sim.time < t_end:
    sim.advance(sim.time + dt)
    states_reactor.append(r.thermo.state, t=sim.time)
#   print('{:10.3e} {:10.3f} {:10.3f} {:14.6f}'.format(
#           sim.time, r.T, r.thermo.P, r.thermo.u))



y0 = np.hstack((gas_ode.T, gas_ode.Y))

# Set up objects representing the ODE and the solver
ode = ReactorOde(gas_ode)
solver = scipy.integrate.ode(ode)
solver.set_integrator('vode', method='bdf', with_jacobian=True)
solver.set_initial_value(y0, 0.0)

# Integrate the equations, keeping T(t) and Y(k,t)
# t_end = 1e-3
states_ode = ct.SolutionArray(gas_ode, 1, extra={'t': [0.0]})

while solver.successful() and solver.t < t_end:
    solver.integrate(solver.t + dt)
    gas_ode.TPY = solver.y[0], gas_ode.P, solver.y[1:]
    states_ode.append(gas_ode.state, t=solver.t)

states_reactor.t = states_reactor.t * 1e3
states_ode.t = states_ode.t * 1e3

import matplotlib.pyplot as plt
plt.clf()
plt.subplot(2, 2, 1)
plt.plot(states_reactor.t, states_reactor.T, label='Reactor')
plt.plot(states_ode.t, states_ode.T, label='ODE')
plt.xlabel('Time (ms)')
plt.ylabel('Temperature [K]')
plt.legend()
plt.subplot(2, 2, 2)
plt.plot(states_reactor.t, states_reactor.P/1e5, label='Reactor')
plt.plot(states_ode.t, states_ode.P/1e5, label='ODE')
plt.xlabel('Time (ms)')
plt.ylabel('Pressure [bar]')
plt.legend()
plt.subplot(2, 2, 3)
plt.plot(states_reactor.t, states_reactor.X[:, gas_reactor.species_index('OH')], label='Reactor')
plt.plot(states_ode.t, states_ode.X[:, gas_ode.species_index('OH')], label='ODE')
plt.xlabel('Time (ms)')
plt.ylabel('H Mole Fraction')
plt.legend()
plt.subplot(2, 2, 4)
plt.plot(states_reactor.t, states_reactor.X[:, gas_reactor.species_index('H2')], label='Reactor')
plt.plot(states_ode.t, states_ode.X[:, gas_ode.species_index('H2')], label='ODE')
plt.xlabel('Time (ms)')
plt.ylabel('H2 Mole Fraction')
plt.legend()
plt.tight_layout()
plt.show()

[Bryan Weber]

Hi,

Based on the line:

self.TP = y[0], self.gas.P

I think you're setting the T and P of the gas, which does not allow the P to float. If you want a constant volume and constant mass, then T and rho should be used to set the state, I think. You'll also need to change how the energy equation is computed, to use c_v instead of c_p. Check out the similar Matlab example: https://cantera.org/examples/matlab/conuv.m.html

Best,

Bryan

[George Eappachan]

Thanks, Bryan. This fixed the difference in results. The profiles agree well now.


--G