From e78aac7b70265705aef4d62cde8c675003dc50e5 Mon Sep 17 00:00:00 2001 From: Ray Speth Date: Mon, 18 Sep 2017 20:17:05 -0400 Subject: [PATCH] [Examples] Clean up NonIdealShockTube example Eliminate pandas dependency and simplify some Matplotlib usage --- data/inputs/nDodecane_Reitz.cti | 18 +- .../examples/reactors/NonIdealShockTube.py | 215 +++++++----------- 2 files changed, 98 insertions(+), 135 deletions(-) diff --git a/data/inputs/nDodecane_Reitz.cti b/data/inputs/nDodecane_Reitz.cti index 9990fb8e3..a9b583742 100644 --- a/data/inputs/nDodecane_Reitz.cti +++ b/data/inputs/nDodecane_Reitz.cti @@ -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]) - - diff --git a/interfaces/cython/cantera/examples/reactors/NonIdealShockTube.py b/interfaces/cython/cantera/examples/reactors/NonIdealShockTube.py index d2796126b..e986f2a04 100644 --- a/interfaces/cython/cantera/examples/reactors/NonIdealShockTube.py +++ b/interfaces/cython/cantera/examples/reactors/NonIdealShockTube.py @@ -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() - - - -