[Examples] Clean up NonIdealShockTube example
Eliminate pandas dependency and simplify some Matplotlib usage
This commit is contained in:
Ray Speth
parent
2a601c148f
commit
e78aac7b70
2 changed files with 98 additions and 135 deletions
|
|
@ -1,8 +1,16 @@
|
|||
"""
|
||||
Real gas n-dodecane-PAH mechanism.
|
||||
Mechanism reported in ‘Development of a reduced n-dodecane-PAH mechanism and its Application for n-dodecane Soot Predictions’
|
||||
Hu Wang, Youngchul Ra, Ming Jia and Rolf. D. Reitz. Fuel 136 (2014), p 25-36.
|
||||
Real gas n-dodecane-PAH mechanism.
|
||||
|
||||
Mechanism reported in ‘Development of a reduced n-dodecane-PAH mechanism and its
|
||||
Application for n-dodecane Soot Predictions’. Hu Wang, Youngchul Ra, Ming Jia
|
||||
and Rolf. D. Reitz. Fuel 136 (2014), p 25-36. doi:10.1016/j.fuel.2014.07.028
|
||||
|
||||
100 species and 432 reactions
|
||||
|
||||
Redlich-Kwong coefficients are based on tabulated critical properties or
|
||||
estimated according to the method of Joback and Reid, "Estimation of pure-
|
||||
component properties from group-contributions," Chem. Eng. Comm. 57 (1987)
|
||||
233-243
|
||||
"""
|
||||
|
||||
units(length='cm', time='s', quantity='mol', act_energy='cal/mol')
|
||||
|
|
@ -165,7 +173,7 @@ ideal_gas(name="nDodecane_IG",
|
|||
c6h5o A1- A1c2h- A1c2h
|
||||
A1c2h3 A2- A2r5 A3-
|
||||
A1 A2 A3 A4""",
|
||||
reactions='all',
|
||||
reactions='all',
|
||||
initial_state=state(temperature=300.0, pressure=OneAtm))
|
||||
|
||||
|
||||
|
|
@ -3113,5 +3121,3 @@ reaction('A3- + c2h2 <=> A4 + h', [6.600000e+24, -3.36, 17680.0])
|
|||
|
||||
# Reaction 553
|
||||
reaction('A4 + oh <=> A3- + ch2co', [2.000000e+13, 0.0, 41730.0])
|
||||
|
||||
|
||||
|
|
|
|||
|
|
@ -1,76 +1,79 @@
|
|||
|
||||
# coding: utf-8
|
||||
|
||||
# # Non-Ideal Shock Tube Example
|
||||
# Ignition delay time computations in a high-pressure reflected shock tube reactor
|
||||
#
|
||||
# In this example we illustrate how to setup and use a constant volume, adiabatic reactor to simulate
|
||||
# reflected shock tube experiments. This reactor will then be used to compute the ignition delay of
|
||||
# a gas at a specified initial temperature and pressure. The example is written in a general way,
|
||||
# i.e., no particular EoS is presumed and ideal and real gas EoS can be used equally easily.
|
||||
#
|
||||
# The reactor (system) is simply an 'insulated box,' and can technically be used for any number of
|
||||
# equations of state and constant-volume, adiabatic reactors.
|
||||
#
|
||||
# Other than the typical Cantera dependencies, plotting functions require that you have matplotlib
|
||||
# installed, and data storing and analysis requires pandas. See https://matplotlib.org/ and
|
||||
# http://pandas.pydata.org/index.html, respectively, for additional info.
|
||||
# Non-Ideal Shock Tube Example
|
||||
#
|
||||
# Ignition delay time computations in a high-pressure reflected shock tube
|
||||
# reactor
|
||||
#
|
||||
# In this example we illustrate how to setup and use a constant volume,
|
||||
# adiabatic reactor to simulate reflected shock tube experiments. This reactor
|
||||
# will then be used to compute the ignition delay of a gas at a specified
|
||||
# initial temperature and pressure. The example is written in a general way,
|
||||
# i.e., no particular EoS is presumed and ideal and real gas EoS can be used
|
||||
# equally easily.
|
||||
#
|
||||
# The reactor (system) is simply an 'insulated box,' and can technically be used
|
||||
# for any number of equations of state and constant-volume, adiabatic reactors.
|
||||
#
|
||||
# Other than the typical Cantera dependencies, plotting functions require that
|
||||
# you have matplotlib (https://matplotlib.org/) installed.
|
||||
|
||||
from __future__ import division
|
||||
from __future__ import print_function
|
||||
|
||||
# Dependencies: pandas, numpy, and matplotlib.pyplot
|
||||
import pandas as pd
|
||||
# Dependencies: numpy, and matplotlib
|
||||
import numpy as np
|
||||
import matplotlib.pyplot as plt
|
||||
from matplotlib import font_manager
|
||||
|
||||
import time
|
||||
|
||||
import cantera as ct
|
||||
print('Runnning Cantera version: ' + ct.__version__)
|
||||
print('Running Cantera version: ' + ct.__version__)
|
||||
|
||||
|
||||
|
||||
# Define the ignition delay time (IDT). This function computes the ignition delay from the occurence
|
||||
# of the peak concentration for the specified species.
|
||||
def ignitionDelay(df, species):
|
||||
return df[species].argmax()
|
||||
# Define the ignition delay time (IDT). This function computes the ignition
|
||||
# delay from the occurrence of the peak concentration for the specified
|
||||
# species.
|
||||
def ignitionDelay(states, species):
|
||||
i_ign = states(species).Y.argmax()
|
||||
return states.t[i_ign]
|
||||
|
||||
|
||||
# Define the reactor temperature and pressure:
|
||||
reactorTemperature = 1000 #Kelvin
|
||||
reactorPressure = 40.0*101325.0 #Pascals
|
||||
reactorTemperature = 1000 # Kelvin
|
||||
reactorPressure = 40.0*101325.0 # Pascals
|
||||
|
||||
# Define the gas
|
||||
# In this example we will choose a stoichiometric mixture of n-dodecane and air as the gas. For a
|
||||
# representative kinetic model, we use that developed by Wang, Ra, Jia, and Reitz
|
||||
# (https://www.erc.wisc.edu/chem_mech/nC12-PAH_mech.zip) by [H.Wang, Y.Ra, M.Jia, R.Reitz,
|
||||
# Development of a reduced n-dodecane-PAH mechanism. and its application for n-dodecane soot
|
||||
# predictions., Fuel 136 (2014) 25–36]
|
||||
# Define the gas: In this example we will choose a stoichiometric mixture of
|
||||
# n-dodecane and air as the gas. For a representative kinetic model, we use:
|
||||
#
|
||||
# H.Wang, Y.Ra, M.Jia, R.Reitz, Development of a reduced n-dodecane-PAH
|
||||
# mechanism. and its application for n-dodecane soot predictions., Fuel 136
|
||||
# (2014) 25–36. doi:10.1016/j.fuel.2014.07.028
|
||||
|
||||
# R-K constants are calculated according to their critical temperature (Tc) and pressure (Pc):
|
||||
# R-K constants are calculated according to their critical temperature (Tc) and
|
||||
# pressure (Pc):
|
||||
#
|
||||
# a = 0.4275*(R^2)*(Tc^2.5)/(Pc)
|
||||
# a = 0.4275*(R^2)*(Tc^2.5)/(Pc)
|
||||
#
|
||||
# and
|
||||
# and
|
||||
#
|
||||
# b = 0.08664*R*Tc/Pc
|
||||
# b = 0.08664*R*Tc/Pc
|
||||
#
|
||||
# where R is the gas constant.
|
||||
# where R is the gas constant.
|
||||
#
|
||||
# For stable species, the critical properties are readily available. For radicals and other
|
||||
# short-lived intermediates, the Joback method is used to estimate critical properties. See Joback
|
||||
# and Reid, "Estimation of pur-component properties from group-contributions," Chem. Eng. Comm. 57
|
||||
# (1987) 233-243, for details of the method.
|
||||
# For stable species, the critical properties are readily available. For
|
||||
# radicals and other short-lived intermediates, the Joback method is used to
|
||||
# estimate critical properties. For details of the method, see: Joback and Reid,
|
||||
# "Estimation of pure- component properties from group-contributions," Chem.
|
||||
# Eng. Comm. 57 (1987) 233-243, doi: 10.1080/00986448708960487
|
||||
|
||||
# There is a slight discontinuity in the thermo for three species at the mid-point temperatrue. We
|
||||
# are aware and okay, so we will suppress the warning statement (note: use this feature at your own
|
||||
# risk, in other codes!)
|
||||
# There is a slight discontinuity in the thermo for three species at the mid-
|
||||
# point temperature. We are aware and okay, so we will suppress the warning
|
||||
# statement (note: use this feature at your own risk in other codes!)
|
||||
ct.suppress_thermo_warnings()
|
||||
|
||||
|
||||
"""Real gas IDT calculation"""
|
||||
|
||||
# Load the real gas mechanism:
|
||||
real_gas = ct.Solution('nDodecane_Reitz.cti','nDodecane_RK')
|
||||
|
||||
|
|
@ -85,16 +88,12 @@ real_gas.set_equivalence_ratio(phi=1.0, fuel='c12h26',
|
|||
# In this example, this will be the only reactor in the network
|
||||
r = ct.Reactor(contents=real_gas)
|
||||
reactorNetwork = ct.ReactorNet([r])
|
||||
|
||||
# now compile a list of all variables for which we will store data
|
||||
stateVariableNames = [r.component_name(item) for item in range(r.n_vars)]
|
||||
|
||||
# Use the above list to create a DataFrame
|
||||
timeHistory_RG = pd.DataFrame(columns=stateVariableNames)
|
||||
#Tic
|
||||
timeHistory_RG = ct.SolutionArray(real_gas, extra=['t'])
|
||||
# Tic
|
||||
t0 = time.time()
|
||||
|
||||
# This is a starting estimate. If you do not get an ignition within this time, increase it
|
||||
# This is a starting estimate. If you do not get an ignition within this time,
|
||||
# increase it
|
||||
estimatedIgnitionDelayTime = 0.005
|
||||
t = 0
|
||||
|
||||
|
|
@ -104,13 +103,13 @@ while(t < estimatedIgnitionDelayTime):
|
|||
if (counter%20 == 0):
|
||||
# We will save only every 20th value. Otherwise, this takes too long
|
||||
# Note that the species concentrations are mass fractions
|
||||
timeHistory_RG.loc[t] = reactorNetwork.get_state()
|
||||
timeHistory_RG.append(r.thermo.state, t=t)
|
||||
counter+=1
|
||||
|
||||
# We will use the 'oh' species to compute the ignition delay
|
||||
tau_RG = ignitionDelay(timeHistory_RG, 'oh')
|
||||
|
||||
#Toc
|
||||
# Toc
|
||||
t1 = time.time()
|
||||
print('Computed Real Gas Ignition Delay: {:.3e} seconds. Took {:3.2f}s to compute'.format(tau_RG, t1-t0))
|
||||
|
||||
|
|
@ -128,14 +127,9 @@ ideal_gas.set_equivalence_ratio(phi=1.0, fuel='c12h26',
|
|||
|
||||
r = ct.Reactor(contents=ideal_gas)
|
||||
reactorNetwork = ct.ReactorNet([r])
|
||||
timeHistory_IG = ct.SolutionArray(ideal_gas, extra=['t'])
|
||||
|
||||
# now compile a list of all variables for which we will store data
|
||||
stateVariableNames = [r.component_name(item) for item in range(r.n_vars)]
|
||||
|
||||
# Use the above list to create a DataFrame
|
||||
timeHistory_IG = pd.DataFrame(columns=stateVariableNames)
|
||||
|
||||
#Tic
|
||||
# Tic
|
||||
t0 = time.time()
|
||||
|
||||
t = 0
|
||||
|
|
@ -146,56 +140,44 @@ while(t < estimatedIgnitionDelayTime):
|
|||
if (counter%20 == 0):
|
||||
# We will save only every 20th value. Otherwise, this takes too long
|
||||
# Note that the species concentrations are mass fractions
|
||||
timeHistory_IG.loc[t] = reactorNetwork.get_state()
|
||||
timeHistory_IG.append(r.thermo.state, t=t)
|
||||
counter+=1
|
||||
|
||||
# We will use the 'oh' species to compute the ignition delay
|
||||
tau_IG = ignitionDelay(timeHistory_IG, 'oh')
|
||||
|
||||
#Toc
|
||||
# Toc
|
||||
t1 = time.time()
|
||||
|
||||
print('Computed Ideal Gas Ignition Delay: {:.3e} seconds. Took {:3.2f}s to compute'.format(tau_IG, t1-t0))
|
||||
print('Ideal gas error: {:2.2f} %'.format(100*(tau_IG-tau_RG)/tau_RG))
|
||||
|
||||
# If you want to save all the data - molefractions, temperature, pressure, etc
|
||||
# uncomment the next line
|
||||
# timeHistory.to_csv("time_history.csv")
|
||||
|
||||
# Plot the result
|
||||
|
||||
plt.rcParams['axes.labelsize'] = 16
|
||||
plt.rcParams['xtick.labelsize'] = 12
|
||||
plt.rcParams['ytick.labelsize'] = 12
|
||||
plt.rcParams['figure.autolayout'] = True
|
||||
plt.rcParams['axes.labelsize'] = 14
|
||||
plt.rcParams['font.family'] = 'serif'
|
||||
|
||||
# Figure illustrating the definition of ignition delay time (IDT).
|
||||
|
||||
plt.figure()
|
||||
plt.plot(timeHistory_RG.index, timeHistory_RG['oh'],'-o',color='b',markersize=4)
|
||||
plt.plot(timeHistory_IG.index, timeHistory_IG['oh'],'-o',color='r',markersize=4)
|
||||
plt.xlabel('Time (s)',fontname='Times New Roman')
|
||||
plt.ylabel('$\mathdefault{OH\, mass\, fraction,}\, \mathdefault{y_{OH}}$',
|
||||
fontname='Times New Roman')
|
||||
plt.plot(timeHistory_RG.t, timeHistory_RG('oh').Y,'-o',color='b',markersize=4)
|
||||
plt.plot(timeHistory_IG.t, timeHistory_IG('oh').Y,'-o',color='r',markersize=4)
|
||||
plt.xlabel('Time (s)')
|
||||
plt.ylabel(r'OH mass fraction, $\mathdefault{Y_{OH}}$')
|
||||
|
||||
# Figure formatting:
|
||||
plt.xlim([0,0.00055])
|
||||
ax = plt.gca()
|
||||
font = plt.matplotlib.font_manager.FontProperties(family='Times New Roman',size=14)
|
||||
ax.annotate("",xy=(tau_RG,0.005), xytext=(0,0.005),
|
||||
arrowprops=dict(arrowstyle="<|-|>",color='r',linewidth=2.0),
|
||||
arrowprops=dict(arrowstyle="<|-|>",color='k',linewidth=2.0),
|
||||
fontsize=14,)
|
||||
plt.annotate('Ignition Delay Time (IDT)', xy=(0,0), xytext=(0.00008, 0.00525),
|
||||
family='Times New Roman',fontsize=16);
|
||||
fontsize=16);
|
||||
|
||||
for tick in ax.xaxis.get_major_ticks():
|
||||
tick.label1.set_fontsize(12)
|
||||
tick.label1.set_fontname('Times New Roman')
|
||||
for tick in ax.yaxis.get_major_ticks():
|
||||
tick.label1.set_fontsize(12)
|
||||
tick.label1.set_fontname('Times New Roman')
|
||||
|
||||
plt.legend(['Real Gas','Ideal Gas'],prop=font,frameon=0)
|
||||
plt.legend(['Real Gas','Ideal Gas'], frameon=False)
|
||||
|
||||
# If you want to save the plot, uncomment this line (and edit as you see fit):
|
||||
#plt.savefig('IDT_nDodecane_1000K_40atm.pdf',dpi=350,format='pdf')
|
||||
|
|
@ -207,21 +189,17 @@ plt.legend(['Real Gas','Ideal Gas'],prop=font,frameon=0)
|
|||
# with increasing temperature.
|
||||
|
||||
# Make a list of all the temperatures at which we would like to run simulations:
|
||||
T = [1250, 1225, 1200, 1150, 1100, 1075, 1050, 1025, 1012.5, 1000, 987.5, 975, 962.5, 950,
|
||||
937.5, 925, 912.5, 900, 875, 850, 825, 800]
|
||||
T = np.array([1250, 1225, 1200, 1150, 1100, 1075, 1050, 1025, 1012.5, 1000, 987.5,
|
||||
975, 962.5, 950, 937.5, 925, 912.5, 900, 875, 850, 825, 800])
|
||||
|
||||
# If we desire, we can define different IDT starting guesses for each temperature:
|
||||
estimatedIgnitionDelayTimes = np.ones(len(T))
|
||||
# But we won't, at least in this example :)
|
||||
estimatedIgnitionDelayTimes[:] = 0.005
|
||||
|
||||
# Now create a dataFrame for the real gas results:
|
||||
ignitionDelays_RG = pd.DataFrame(data={'T':T})
|
||||
ignitionDelays_RG['ignDelay'] = np.nan
|
||||
|
||||
|
||||
# Now, we simply run the code above for each temperature.
|
||||
"""Real Gas"""
|
||||
ignitionDelays_RG = np.zeros(len(T))
|
||||
for i, temperature in enumerate(T):
|
||||
# Setup the gas and reactor
|
||||
reactorTemperature = temperature
|
||||
|
|
@ -231,8 +209,8 @@ for i, temperature in enumerate(T):
|
|||
r = ct.Reactor(contents=real_gas)
|
||||
reactorNetwork = ct.ReactorNet([r])
|
||||
|
||||
# Create and empty data frame
|
||||
timeHistory = pd.DataFrame(columns=timeHistory_RG.columns)
|
||||
# create an array of solution states
|
||||
timeHistory = ct.SolutionArray(real_gas, extra=['t'])
|
||||
|
||||
t0 = time.time()
|
||||
|
||||
|
|
@ -241,7 +219,7 @@ for i, temperature in enumerate(T):
|
|||
while t < estimatedIgnitionDelayTimes[i]:
|
||||
t = reactorNetwork.step()
|
||||
if not counter % 20:
|
||||
timeHistory.loc[t] = r.get_state()
|
||||
timeHistory.append(r.thermo.state, t=t)
|
||||
counter += 1
|
||||
|
||||
tau = ignitionDelay(timeHistory, 'oh')
|
||||
|
|
@ -249,15 +227,11 @@ for i, temperature in enumerate(T):
|
|||
|
||||
print('Computed Real Gas Ignition Delay: {:.3e} seconds for T={}K. Took {:3.2f}s to compute'.format(tau, temperature, t1-t0))
|
||||
|
||||
ignitionDelays_RG.set_value(index=i, col='ignDelay', value=tau)
|
||||
|
||||
ignitionDelays_RG[i] = tau
|
||||
|
||||
|
||||
"""Repeat for Ideal Gas"""
|
||||
# Create a dataFrame for the ideal gas results:
|
||||
ignitionDelays_IG = pd.DataFrame(data={'T':T})
|
||||
ignitionDelays_IG['ignDelay'] = np.nan
|
||||
|
||||
ignitionDelays_IG = np.zeros(len(T))
|
||||
for i, temperature in enumerate(T):
|
||||
# Setup the gas and reactor
|
||||
reactorTemperature = temperature
|
||||
|
|
@ -267,8 +241,8 @@ for i, temperature in enumerate(T):
|
|||
r = ct.Reactor(contents=ideal_gas)
|
||||
reactorNetwork = ct.ReactorNet([r])
|
||||
|
||||
# Create and empty data frame
|
||||
timeHistory = pd.DataFrame(columns=timeHistory_IG.columns)
|
||||
# create an array of solution states
|
||||
timeHistory = ct.SolutionArray(ideal_gas, extra=['t'])
|
||||
|
||||
t0 = time.time()
|
||||
|
||||
|
|
@ -277,7 +251,7 @@ for i, temperature in enumerate(T):
|
|||
while t < estimatedIgnitionDelayTimes[i]:
|
||||
t = reactorNetwork.step()
|
||||
if not counter % 20:
|
||||
timeHistory.loc[t] = r.get_state()
|
||||
timeHistory.append(r.thermo.state, t=t)
|
||||
counter += 1
|
||||
|
||||
tau = ignitionDelay(timeHistory, 'oh')
|
||||
|
|
@ -285,19 +259,16 @@ for i, temperature in enumerate(T):
|
|||
|
||||
print('Computed Ideal Gas Ignition Delay: {:.3e} seconds for T={}K. Took {:3.2f}s to compute'.format(tau, temperature, t1-t0))
|
||||
|
||||
ignitionDelays_IG.set_value(index=i, col='ignDelay', value=tau)
|
||||
ignitionDelays_IG[i] = tau
|
||||
|
||||
|
||||
# Figure: ignition delay ($\tau$) vs. the inverse of temperature ($\frac{1000}{T}$).
|
||||
# Figure: ignition delay (tau) vs. the inverse of temperature (1000/T).
|
||||
fig = plt.figure()
|
||||
ax = fig.add_subplot(111)
|
||||
ax.plot(1000/ignitionDelays_RG['T'], 1e6*ignitionDelays_RG['ignDelay'],'-',
|
||||
linewidth=2.0,color='b')
|
||||
ax.plot(1000/ignitionDelays_IG['T'], 1e6*ignitionDelays_IG['ignDelay'],'-.',
|
||||
linewidth=2.0,color='r')
|
||||
ax.set_ylabel(r'$\mathdefault{Ignition\, Delay\, (\mu s)}$',fontname='Times New Roman',
|
||||
fontsize=16)
|
||||
ax.set_xlabel(r'$\mathdefault{1000/T\, (K^{-1})}$',fontname='Times New Roman', fontsize=16)
|
||||
ax.plot(1000/T, 1e6*ignitionDelays_RG, '-', linewidth=2.0, color='b')
|
||||
ax.plot(1000/T, 1e6*ignitionDelays_IG, '-.', linewidth=2.0, color='r')
|
||||
ax.set_ylabel(r'Ignition Delay ($\mathdefault{\mu s}$)', fontsize=14)
|
||||
ax.set_xlabel(r'1000/T (K$^\mathdefault{-1}$)', fontsize=14)
|
||||
|
||||
ax.set_xlim([0.8,1.2])
|
||||
|
||||
|
|
@ -307,26 +278,12 @@ ticks = ax.get_xticks()
|
|||
ax2.set_xticks(ticks)
|
||||
ax2.set_xticklabels((1000/ticks).round(1))
|
||||
ax2.set_xlim(ax.get_xlim())
|
||||
ax2.set_xlabel('Temperature (K)',fontname='Times New Roman',fontsize=16);
|
||||
ax2.set_xlabel('Temperature (K)', fontsize=14);
|
||||
|
||||
ticks_font = font_manager.FontProperties(family='Times New Roman', style='normal',
|
||||
size=12, weight='normal',
|
||||
stretch='normal')
|
||||
for label in ax.get_yticklabels():
|
||||
label.set_fontproperties(ticks_font)
|
||||
for label in ax.get_xticklabels():
|
||||
label.set_fontproperties(ticks_font)
|
||||
for label in ax2.get_xticklabels():
|
||||
label.set_fontproperties(ticks_font)
|
||||
|
||||
ax.legend(['Real Gas','Ideal Gas'],prop=font,frameon=0,loc=2)
|
||||
ax.legend(['Real Gas','Ideal Gas'], frameon=False, loc='upper left')
|
||||
|
||||
# If you want to save the plot, uncomment this line (and edit as you see fit):
|
||||
#plt.savefig('NTC_nDodecane_40atm.pdf',dpi=350,format='pdf')
|
||||
|
||||
# Show the plots.
|
||||
plt.show()
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue