h0
is our enthalpy, this calculation is done internally, and we would simply write:gas.HPY = h0, P, y
T = gas.T
Now divide through by $dt$:
$$\frac{dh}{dt} = c_p\frac{dT}{dt} + \sum_ih_i\frac{dy_i}{dt} = 0,$$Use $dy_i/dt = \dot{m}_i^{\prime\prime\prime}/\rho$, and solve for $dT/dt$:
$$\frac{dT}{dt} = -\frac{1}{\rho c_p}\sum_ih_i\dot{m}_i^{\prime\prime\prime}.$$Also,
$$c_p = \sum_iy_ic_{p,i},$$ $$h_i = h_{f,i}(T_{ref}) + \int_{T_{ref}}^{T}c_{p,i}dT.$$Also,
$$c_v = \sum_iy_ic_{v,i},$$ $$u_i = u_{f,i}(T_{ref}) + \int_{T_{ref}}^{T}c_{v,i}dT.$$Hence,
$$\frac{d}{dt} \rightarrow \frac{\dot{m}^{\prime\prime}}{\rho}\frac{d}{dz}$$If we nondimensionalize this equation, then we need a timescale $\tau$, a lengthscale $L$, and a scale $y_{i,ref}$ (which cancels).
The nondimensional form is
as the characteristic diffusion timescale.
But this is just our PSR mixing term. The point is, that the PSR mixing term can be used as an analogy to more complex diffusive mixing in terms of a mixing timescale, where mixing and reaction processes are in a kind of balance.
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
from scipy.integrate import odeint
import cantera as ct
gas = ct.Solution("gri30.yaml")
#------------- Set batch reactor equation
def rhsf(y,t):
gas.HPY = h0, gas.P, y
return gas.net_production_rates * gas.molecular_weights / gas.density
#------------- set initial condition and enthalpy
gas.TPX = 1400, 101325, "CH4:1, O2:2, N2:7.52"
h0 = gas.enthalpy_mass
y0 = gas.Y
#------------- solve the ODE system
nt = 1000
times = np.linspace(0,0.01,nt)
y = odeint(rhsf, y0, times)
#------------- recover the temperature
T = np.zeros(nt)
for i in range(nt):
gas.HPY = h0, gas.P, y[i,:]
T[i] = gas.T
#------------- equilibrium state
gas.equilibrate("HP")
Teq = gas.T
yeq = gas.Y
#------------- plot results
plt.rc("font", size=14)
plt.plot(times*1000,T, label="T")
plt.plot(times[-1]*1000,Teq, 'kd', label="Teq")
plt.xlabel("time (ms)")
plt.ylabel("T (K)")
plt.legend(frameon=False);
plt.figure()
species = ["CH4", "O2", "CO2", "CO", "OH"]
for sp in species:
plt.plot(times*1000,y[:,gas.species_index(sp)], label=sp)
plt.plot(times[-1]*1000,yeq[gas.species_index(sp)], 'kd', label="eq" if sp==species[-1] else "")
plt.xlabel("time (ms)")
plt.ylabel("y")
plt.legend(frameon=False);
plt.plot(times*1000, y[:,gas.species_index("CO2")]/y[-1,gas.species_index("CO2")])
plt.plot(times*1000, y[:,gas.species_index("NO")]/y[-1,gas.species_index("NO")])
[<matplotlib.lines.Line2D at 0x106d81be0>]