mirror of
https://github.com/Cantera/cantera.git
synced 2025-02-25 18:55:29 -06:00
c49c95d666
Introduce a monkey patch to avoid executing specific examples that won't run under sphinx-gallery.
77 lines
2.3 KiB
Python
77 lines
2.3 KiB
Python
"""
|
|
Constant-pressure, adiabatic kinetics simulation with sensitivity analysis
|
|
==========================================================================
|
|
|
|
Requires: cantera >= 2.5.0, matplotlib >= 2.0
|
|
|
|
.. tags:: Python, combustion, reactor network, sensitivity analysis, plotting
|
|
"""
|
|
|
|
import sys
|
|
import numpy as np
|
|
|
|
import cantera as ct
|
|
|
|
gas = ct.Solution('gri30.yaml')
|
|
temp = 1500.0
|
|
pres = ct.one_atm
|
|
|
|
gas.TPX = temp, pres, 'CH4:0.1, O2:2, N2:7.52'
|
|
r = ct.IdealGasConstPressureReactor(gas, name='R1')
|
|
sim = ct.ReactorNet([r])
|
|
|
|
# enable sensitivity with respect to the rates of the first 10
|
|
# reactions (reactions 0 through 9)
|
|
for i in range(10):
|
|
r.add_sensitivity_reaction(i)
|
|
|
|
# set the tolerances for the solution and for the sensitivity coefficients
|
|
sim.rtol = 1.0e-6
|
|
sim.atol = 1.0e-15
|
|
sim.rtol_sensitivity = 1.0e-6
|
|
sim.atol_sensitivity = 1.0e-6
|
|
|
|
states = ct.SolutionArray(gas, extra=['t', 's2', 's3'])
|
|
|
|
for t in np.arange(0, 2e-3, 5e-6):
|
|
sim.advance(t)
|
|
s2 = sim.sensitivity('OH', 2) # sensitivity of OH to reaction 2
|
|
s3 = sim.sensitivity('OH', 3) # sensitivity of OH to reaction 3
|
|
states.append(r.thermo.state, t=1000*t, s2=s2, s3=s3)
|
|
|
|
print('{:10.3e} {:10.3f} {:10.3f} {:14.6e} {:10.3f} {:10.3f}'.format(
|
|
sim.time, r.T, r.thermo.P, r.thermo.u, s2, s3))
|
|
|
|
# plot the results if matplotlib is installed.
|
|
# see http://matplotlib.org/ to get it
|
|
if '--plot' in sys.argv:
|
|
import matplotlib.pyplot as plt
|
|
plt.subplot(2, 2, 1)
|
|
plt.plot(states.t, states.T)
|
|
plt.xlabel('Time (ms)')
|
|
plt.ylabel('Temperature (K)')
|
|
plt.subplot(2, 2, 2)
|
|
plt.plot(states.t, states('OH').X)
|
|
plt.xlabel('Time (ms)')
|
|
plt.ylabel('OH Mole Fraction')
|
|
plt.subplot(2, 2, 3)
|
|
plt.plot(states.t, states('H').X)
|
|
plt.xlabel('Time (ms)')
|
|
plt.ylabel('H Mole Fraction')
|
|
plt.subplot(2, 2, 4)
|
|
plt.plot(states.t, states('CH4').X)
|
|
plt.xlabel('Time (ms)')
|
|
plt.ylabel('CH4 Mole Fraction')
|
|
plt.tight_layout()
|
|
|
|
plt.figure(2)
|
|
plt.plot(states.t, states.s2, '-', label=sim.sensitivity_parameter_name(2))
|
|
plt.plot(states.t, states.s3, '-g', label=sim.sensitivity_parameter_name(3))
|
|
plt.legend(loc='best')
|
|
plt.xlabel('Time (ms)')
|
|
plt.ylabel('OH Sensitivity')
|
|
plt.tight_layout()
|
|
plt.show()
|
|
else:
|
|
print("""To view a plot of these results, run this script with the option '--plot""")
|