import os
import cantera as ct
import numpy as np
%matplotlib inline
from matplotlib import pyplot as plt
import re
import IPython
from IPython.display import display, Markdown, HTML
def mprint(s):
"A convenience to format things in Jupyter notebooks via MarkDown syntax"
display(Markdown(s))
cantera_file_name = 'butanol.original.cti'
cantera_files_directory = 'CanteraFiles'
cantera_file_path = os.path.join(cantera_files_directory,cantera_file_name)
print os.path.abspath(cantera_file_path)
assert os.path.exists(cantera_file_path)
original = ct.Solution(cantera_file_path)
# Find the fuel
for species in original.species():
if species.composition == {'C':4, 'H':10, 'O':1}:
print (species.name)
def get_ignition_delay(model, temperature, pressure=10, phi=1.0, plot=False, isomer='n'):
"""
Get the ignition delay at temperature (K) and pressure (bar) and stochiometry (phi),
for the butanol isomer (n,s,t,i)
"""
assert isomer in ['n','s','t','i'], "Expecting isomer n,s,t, or i not {}".format(isomer)
oxygen_mole = 1.0
argon_mole = 96./4.*oxygen_mole # 4% O2 in Ar
butanol_mole = phi * oxygen_mole/6.
X_string = isomer + 'c4h9oh:{0}, o2:{1}, ar:{2}'.format(butanol_mole, oxygen_mole, argon_mole)
model.TPX = temperature, pressure*1e5, X_string
reactor=ct.IdealGasReactor(model)
reactor_network=ct.ReactorNet([reactor])
time=0.0
end_time=25e-3
times=[]
concentrations=[]
pressures=[]
temperatures=[]
print_data = True
while time < end_time:
time=reactor_network.time
times.append(time)
temperatures.append(reactor.T)
pressures.append(reactor.thermo.P)
concentrations.append(reactor.thermo.concentrations)
reactor_network.step()
print("reached end time {0:.4f} ms in {1} steps ". format(times[-1]*1e3, len(times)))
concentrations = np.array(concentrations)
times = np.array(times)
pressures = np.array(pressures)
temperatures = np.array(temperatures)
dTdt = (temperatures[1:] - temperatures[:-1]) / (times[1:] - times[:-1])
if plot:
plt.subplot(2,1,1)
plt.plot(times,temperatures)
plt.subplot(2,1,2)
plt.plot(times[1:], dTdt)
plt.show()
step_with_fastest_T_rise = dTdt.argmax()
if step_with_fastest_T_rise > 1 and step_with_fastest_T_rise < len(times)-2:
ignition_time_ms = 1e3 * times[step_with_fastest_T_rise]
print("At {0} K {1} bar, ignition delay time is {2} ms for {3}-butanol".format(temperature,pressure,ignition_time_ms,isomer))
return ignition_time_ms
else:
print("At {0} K {1} bar, no ignition is detected for {2}-butanol" .format(temperature, pressure, isomer))
return np.infty
"""
# For excited OH emission
i_ch=gas.species_index('ch')
i_o2=gas.species_index('o2')
excited_oh_generation=concentrations[:,i_o2] * concentrations[:,i_ch]
step_with_highest_oh_gen = excited_oh_generation.argmax()
if step_with_highest_oh_gen > 1 and excited_oh_generation.max() > 1e-20:
ignition_time_ms = 1e3 * times[step_with_highest_oh_gen]
print("At {0} K {1} bar, ignition delay time is {2} ms".format(temperature,pressure,ignition_time_ms))
return ignition_time_ms
else:
print("At {0} K {1} bar, no ignition is detected" .format(temperature, pressure))
return np.infty
"""
get_ignition_delay(original, 1000/.85,43.56, plot=True, isomer='t')
temperatures=1./np.linspace(1./1000,1./1400,8)
def get_original_delays(phi=1.0):
"""
Get the ignition delay times, sorted by pressure then temperature,
for the requested phi, using the original model.
"""
original_idt_by_isomer = dict()
for isomer in 'tnis':
ignition_delay_times = np.zeros_like(temperatures)
for i,T in enumerate(temperatures):
P = 43 * 1.01325
ignition_delay_times[i] = get_ignition_delay(original, T, P, isomer=isomer, phi=phi)
original_idt_by_isomer[isomer] = ignition_delay_times
return original_idt_by_isomer
original_idt_by_isomer = get_original_delays(phi=1.0)
plt.figure(figsize=(6,5))
for isomer, color in zip('tnis', 'mgbr'):
plt.semilogy(1000./temperatures, original_idt_by_isomer[isomer],'o-'+color, label='{0}-butanol'.format(isomer))
plt.legend(loc='best')
plt.xlabel("1000K / temperature")
plt.ylabel("Ignition delay time (ms)")
plt.ylim(1e-2,10)
plt.show()
Compare that with Sarathy et al's Figure 27b and it looks pretty good (especially now that I'm using Ar not N2 bath). They were a bit vague about what pressure they ran at, saying "near 43 atm" and "The initial pressure was varied according to the average pressure of the measurements, thereby better representing the actual conditions of the experiments". Makes it hard to reproduce exactly. Probably the pressures reported in the figure legend, rather than 43 atm as used here. (Also reports "CHEMKIN PRO using constant volume homogeneous batch reactor simulations at $\phi$ = 1"). Anyway, good enough to be confident we're using the model correctly.
def get_roaming_delays(alpha=0.1, P_atm=43, phi=1.0):
"""
For the given alpha, and phi, find the new ignition delay times
and return them for comparison to the previously calculated
ones.
"""
cantera_file_path = os.path.join(cantera_files_directory,'butanol.roaming.cti')
assert os.path.exists(cantera_file_path)
cti_def = open(cantera_file_path).read()
cti_def, substitutions = re.subn('ALPHA = [0-9.]+', 'ALPHA = {0:f}'.format(alpha), cti_def)
assert substitutions == 1
for line in cti_def.splitlines():
if line.startswith('ALPHA'): print line
# Write to temporary file and read from disk, see https://github.com/Cantera/cantera/issues/416
temp_file_path = os.path.join(cantera_files_directory,'butanol.temporary.cti')
with open(temp_file_path,'w') as output_file:
output_file.write(cti_def)
roaming = ct.Solution(temp_file_path)
roaming_idt_by_isomer = dict()
for isomer in 'tnis':
ignition_delay_times = np.zeros_like(temperatures)
for i,T in enumerate(temperatures):
P_bar = P_atm * 1.01325
ignition_delay_times[i] = get_ignition_delay(roaming, T, P_bar, isomer=isomer, phi=phi)
roaming_idt_by_isomer[isomer] = ignition_delay_times
return roaming_idt_by_isomer
# Check there is absolutely no change when ALPHA=0
roaming_idt_by_isomer = roaming00_idt_by_isomer = get_roaming_delays(alpha=0.0)
def plot_both(temperatures, original_idt_by_isomer, roaming_idt_by_isomer):
plt.figure(figsize=(6,6))
for isomer, color in zip('tnis', 'mgbr'):
plt.semilogy(1000./temperatures, original_idt_by_isomer[isomer],'o-'+color, label='{0}-butanol'.format(isomer))
plt.semilogy(1000./temperatures, roaming_idt_by_isomer[isomer],'o:'+color, label=None)
plt.legend(loc='best')
plt.xlabel("1000K / temperature")
plt.ylabel("Ignition delay time(ms)")
plt.ylim(1e-2,)
plt.show()
plot_both(temperatures, original_idt_by_isomer, roaming_idt_by_isomer)
def plot_comparison(temperatures, original_idt_by_isomer, roaming_idt_by_isomer):
print "For temperatures of {!s} K".format(temperatures)
for isomer in 'tnis':
change = (roaming_idt_by_isomer[isomer]-original_idt_by_isomer[isomer])/original_idt_by_isomer[isomer]
print "For {!s}-butanol there's a {!s} percent change in IDT".format(isomer, change*100)
from matplotlib import gridspec
plt.figure(figsize=(5,5))
gs = gridspec.GridSpec(2, 1,height_ratios=[2, 1])
isomers, colors, markers = 'tsin', 'krbg', 'o^sv'
ax0 = plt.subplot(gs[0])
for isomer, color, marker in zip(isomers, colors, markers):
plt.semilogy(1000./temperatures,
original_idt_by_isomer[isomer],
marker=marker,
color=color,
linestyle='-',
label='{0}-butanol'.format(isomer),
)
plt.semilogy(1000./temperatures,
roaming_idt_by_isomer[isomer],
marker=None,
color=color,
linestyle=':',
label=None,
)
plt.legend(loc='best', frameon=False)
plt.xlabel("1000K / Temperature")
plt.ylabel("Ignition delay time (ms)")
plt.ylim(0.01, 20)
ax1 = plt.subplot(gs[1])
max_change = 0
for isomer, color, marker in zip(isomers, colors, markers):
changes = (roaming_idt_by_isomer[isomer]-original_idt_by_isomer[isomer]) / original_idt_by_isomer[isomer]
ax1.plot(1000./temperatures,
100.*changes,
linestyle=':',
marker=None,
color=color,
)
max_change = max(max_change, max(changes*100))
plt.ylabel("Change")
plt.xlabel("Temperature")
plt.xlim(ax0.get_xlim())
from matplotlib.ticker import FuncFormatter
def reciprocal(x, pos):
'The two args are the value and tick position'
return '{0:.0f}K'.format(1000./x)
def percentage(x, pos):
'The two args are the value and tick position'
return '{0:g}%'.format(x)
formatter = FuncFormatter(reciprocal)
ax1.xaxis.set_major_formatter(formatter)
ax1.yaxis.set_major_formatter(FuncFormatter(percentage))
ax1.spines["top"].set_visible(False)
ax1.spines["right"].set_visible(False)
ax1.spines["bottom"].set_visible(False)
plt.axhline(y=0, color='k', linewidth=1)
if max_change <= 8:
plt.ylim(-1,8)
plt.yticks([0,4,8])
plt.tight_layout()
mprint("## Original model ($-$) vs 0% roaming ($\cdots$)\n(should be identical)")
plot_comparison(temperatures, original_idt_by_isomer, roaming00_idt_by_isomer)
plt.savefig('ignition-butanol-00pc.pdf')
# See how much change when ALPHA = 0.1 (reasonable upper limit)
roaming10_idt_by_isomer = get_roaming_delays(alpha=0.10)
# Redo with alpha=0.1 but phi = 0.5
original_idt_by_isomer_phi05 = get_original_delays(phi=0.5)
roaming10_idt_by_isomer_phi05 = get_roaming_delays(alpha=0.1, phi=0.5)
# Redo with alpha=0.1 but phi = 2.0
original_idt_by_isomer_phi20 = get_original_delays(phi=2.0)
roaming10_idt_by_isomer_phi20 = get_roaming_delays(alpha=0.1, phi=2.0)
mprint("## Difference between reference cases for $\phi=1 (-)$ and $\phi=0.5 (\cdots)$ ")
plot_comparison(temperatures, original_idt_by_isomer, original_idt_by_isomer_phi05)
mprint("## 10% roaming with $\phi = 0.5$")
plot_comparison(temperatures, original_idt_by_isomer_phi05, roaming10_idt_by_isomer_phi05)
plt.savefig('ignition-butanol-10pc-phi05.pdf')
mprint("## 10% roaming with $\phi = 1$")
plot_comparison(temperatures, original_idt_by_isomer, roaming10_idt_by_isomer)
plt.savefig('ignition-butanol-10pc.pdf')
mprint("## 10% roaming with $\phi = 2$")
plot_comparison(temperatures, original_idt_by_isomer_phi20, roaming10_idt_by_isomer_phi20)
plt.savefig('ignition-butanol-10pc-phi20.pdf')
# See how much change when ALPHA = 0.05 (moderate?)
roaming05_idt_by_isomer = get_roaming_delays(alpha=0.05)
# See how much change when ALPHA = 0.5 (excessive)
roaming50_idt_by_isomer = get_roaming_delays(alpha=0.5)
mprint("## 5% roaming ($\phi=1$)")
plot_comparison(temperatures, original_idt_by_isomer, roaming05_idt_by_isomer)
plt.savefig('ignition-butanol-05pc.pdf')
mprint("## 50% roaming ($\phi=1$)")
plot_comparison(temperatures, original_idt_by_isomer, roaming50_idt_by_isomer)
plt.savefig('ignition-butanol-50pc.pdf')