skyreflectedinmirrors
commented
6 years ago @chengdi123000
From the Paper (3.2. CEMA for the lifted flame):
. The least-negative eigenvalue aside from the zero eigenvalues associated with the conservation modes is displayed for non-explosive mixtures in figure 6
From this I infer they discarded zero-eigenvalues (resulting from conservation of mass / species / energy, as we discussed the other day).
Even after discarding eigenvalues identically equal to zero, I was not able to reproduce the plot you mention in the paper. However, by filtering eigenvalues by some small threshold (1e-5 in my case), I was able to achieve the expected behavior post-ignition, here's the updated example code:
%matplotlib inline
import matplotlib
matplotlib.rcParams['figure.figsize']=(12,10)
import matplotlib.pyplot as plt
#import matplotlib.animation as manimation
import pyjacob
from scipy.linalg import eig,eigvals,norm
from numpy import argmax,argmin,count_nonzero
import pyjac as pj
import cantera as ct
import numpy as np
import scipy as sp
from scipy.io import savemat
## problem setup
gas = ct.Solution('gri30.cti')
#reorder the gas to match pyJac
n2_ind = gas.species_index('N2')
specs = gas.species()[:]
specs2 = specs[:n2_ind] + specs[n2_ind + 1:] + [specs[n2_ind]]
gas = ct.Solution(thermo='IdealGas', kinetics='GasKinetics', species=specs2, reactions=gas.reactions())
Tin = 1100
P = ct.one_atm
gas.TPY = Tin, P, "H2:0.125, O2:1, N2:3.76"
# get the diabatic flame temperature Taf
oldstate = gas.TPY
gas.equilibrate('HP')
Taf = gas.T
gas.TPY=oldstate
r = ct.IdealGasConstPressureReactor(gas, name='R1')
sim = ct.ReactorNet([r])
## Henrici Measure
def HenriciMeasure(A):
return norm(jac.dot(jac.T)-jac.T.dot(jac),ord=2)/(norm(jac,ord=2)**2)
ts=[]
Ts=[]
num_abs=[]
HM=np.array([])
ermax=np.array([],dtype=complex)
ermin=np.array([],dtype=complex)
eimax=np.array([],dtype=complex)
ermax_nonzero=np.array([],dtype=complex)
Jacs=[]
g=[]
## simulation
t=0.0
print('Reactor state: time[s], temp[K], press[Pa], specific internal energy[J/kg]')
#sim.advance(1e-4)
for t in np.arange(0e-6,120e-6,1e-6):
sim.advance(t)
print('Reactor state: %10.3e %10.3f %10.3f %14.6e' %(sim.time, r.T, r.thermo.P, r.thermo.u))
#setup the state vector
y = np.zeros(gas.n_species)
y[0] = r.T
y[1:] = gas.Y[:-1]
#create a dydt vector
dydt = np.zeros_like(y)
pyjacob.py_dydt(0, P, y, dydt)
#create a jacobian vector
jac = np.zeros(gas.n_species * gas.n_species)
#evaluate the Jacobian
pyjacob.py_eval_jacobian(0, P, y, jac)
# reshape, and transpose
jac = jac.reshape((gas.n_species,gas.n_species)).T
# rescaling using (Taf - Tin)
jac[0,:] *= 1.0/(Taf -Tin)
jac[:,0] *= (Taf -Tin)
# delete argon
#jac=np.delete(jac,(ar_ind),axis=0)
#jac=np.delete(jac,(ar_ind),axis=1)
# get eigen
evals = eigvals(jac)
print('max real, min real, max imag, max nonzero, HenriciMeasure')
ts.append(t)
Ts.append(r.T)
HM = np.append(HM,HenriciMeasure(jac))
ermax=np.append(ermax,evals[argmax(evals.real)])
ermin=np.append(ermin,evals[argmin(evals.real)])
eimax=np.append(eimax,evals[argmax(evals.imag)])
# test nonzero
epsilon = 1e-5
nonzero = evals[np.where(np.abs(evals.real) > epsilon)]
ermax_nonzero = np.append(ermax_nonzero, nonzero[argmax(nonzero.real)])
num_abs.append(np.max(eigvals((jac+jac.T)/2)))
Jacs.append(jac)
g.append(gas.gibbs_mass)
print(ermax[-1],
ermin[-1],
eimax[-1],
ermax_nonzero[-1],
HM[-1])
matplotlib.rcParams.update({'font.size': 18})
from matplotlib.ticker import FormatStrFormatter
fig, ax1 = plt.subplots()
ax1.plot(ts, ermax_nonzero.real, 'b-')
ax1.set_xlabel('Residence time (s)')
# Make the y-axis label, ticks and tick labels match the line color.
ax1.set_ylabel('$\Re(\lambda)_{max} (Hz)$', color='b')
ax1.tick_params('y', colors='b')
ax1.set_xlim(60e-6,120e-6)
ax1.set_ylim(-1e3,1e3)
ax1.plot([0, ts[-1]], [0, 0], 'k-')
ax2 = ax1.twinx()
ax2.plot(ts, Ts, 'r-')
ax2.set_ylabel('Temperature (K)', color='r')
ax2.set_ylim(1000,2500)
ax2.tick_params('y', colors='r')
ax1.yaxis.set_major_formatter(FormatStrFormatter('%.2e'))
ax1.xaxis.set_major_formatter(FormatStrFormatter('%.2e'))
fig.tight_layout()
plt.show()This gives the 
Considering the scale of eigenvalues involved (maximum appears to be O(10^5)), this seems like a reasonable floating point precision threshold to me. You might try emailing the authors directly to see if they utilized something similar?
ghost
I doubt if CEMA is usable to other mechanism. Or the pyjac gave me wrong Jacobian.
I try to reproduce the figure 1-b of the following paper. I assumed that the GRI3.0 Mech in cantera is usable for H2-air autoignition as the mechanism from
Li, J., Zhao, Z. W., Kazakov, A. & Dryer, F. L. 2004 An updated comprehensive kinetic model of hydrogen combustion. Intl J. Chem. Kinet. 36, 566–575., which was used in Lu's paper.The problem is I did not detect the zero crossing. It is zero, but not crossing. And the transition is a very hard turn.
Paper: Lu, T. F., C. S. Yoo, J. H. Chen and C. K. Law (2010). "Three-dimensional direct numerical simulation of a turbulent lifted hydrogen jet flame in heated coflow: a chemical explosive mode analysis." Journal of Fluid Mechanics 652: 45-64.
Figure 1:Temperature and selected species profiles in auto-ignition of a stoichiometric hydrogen–air mixture under atmospheric pressure and an initial temperature of 1100 K. (a) Temperature and species concentrations and (b) time scale, λexp , of the explosive mode and the defectiveness index, dexp , of the chemical Jacobian. The dotted lines indicate the crossover point, where λexp =0.
code:
The
ermax.realis: