import os
import cantera as ct
import numpy as np
%matplotlib inline
from matplotlib import pyplot as plt
import re
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 = 'heptane.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':7, 'H':16}:
print (species.name)
def get_ignition_delay(model, temperature, pressure=10, phi=1.0, plot=False):
"""
Get the ignition delay at temperature (K) and pressure (bar) and equivalence ratio (phi),
for n-heptane. Oxidizer/bath is 20% O2 and 80% Nitrogen.
"""
fuel = 'NC7H16'
model.set_equivalence_ratio(phi, fuel, 'O2:1.0, N2:4.0')
model.TP = temperature, pressure*1e5
reactor=ct.IdealGasReactor(model)
reactor_network=ct.ReactorNet([reactor])
time=0.0
end_time=5.0 # second
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}".format(temperature,pressure,ignition_time_ms, fuel))
return ignition_time_ms
else:
print("At {0} K {1} bar, no ignition is detected for {2}" .format(temperature, pressure, fuel))
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, 600, 10, plot=True)
temperatures=1./np.linspace(1./600,1./2000,15)
pressures = [1,50]
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_pressure = dict()
for pressure_atm in pressures:
ignition_delay_times = np.zeros_like(temperatures)
for i,T in enumerate(temperatures):
P_bar = pressure_atm * 1.01325
ignition_delay_times[i] = get_ignition_delay(original, T, P_bar, phi=phi)
original_idt_by_pressure[pressure_atm] = ignition_delay_times
return original_idt_by_pressure
original_idt_by_pressure = get_original_delays(phi=1.0)
plt.figure(figsize=(6,5))
for pressure, color in zip(pressures, 'rb'):
plt.semilogy(1000./temperatures, original_idt_by_pressure[pressure],'o-'+color, label='{0} atm'.format(pressure))
plt.legend(loc='best')
plt.xlabel("1000K / temperature")
plt.ylabel("Ignition delay time (ms)")
plt.ylim()
plt.show()
This is figure 8 from Curran, H. J.; Gaffuri, P.; Pitz, W. J.; Westbrook, C. K. A Comprehensive Modeling Study of n-Heptane Oxidation. Combust. Flame 1998, 114 (1-2), 149–177 DOI: 10.1016/S0010-2180(97)00282-4.
def get_roaming_delays(alpha=0.1, 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,'heptane.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,'heptane.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_pressure = dict()
for pressure in pressures:
ignition_delay_times = np.zeros_like(temperatures)
for i,T in enumerate(temperatures):
P_bar = pressure * 1.01325
ignition_delay_times[i] = get_ignition_delay(roaming, T, P_bar, phi=phi)
roaming_idt_by_pressure[pressure] = ignition_delay_times
return roaming_idt_by_pressure
#temp_file_path = os.path.join(cantera_files_directory,'heptane.temporary.cti')
#roaming = ct.Solution(temp_file_path)
#get_ignition_delay(roaming, 1000/.85,10, plot=True)
# Check there is absolutely no change when ALPHA=0
roaming00_idt_by_pressure = get_roaming_delays(alpha=0)
def plot_comparison(temperatures, original_idt_by_pressure, roaming_idt_by_pressure):
print "For temperatures of {!s} K".format(temperatures)
for pressure in pressures:
change = (roaming_idt_by_pressure[pressure]-original_idt_by_pressure[pressure])/original_idt_by_pressure[pressure]
print "At {!s} atm there's a {!s} percent change in IDT".format(pressure, change*100)
from matplotlib import gridspec
plt.figure(figsize=(5,5))
gs = gridspec.GridSpec(2, 1,height_ratios=[2, 1])
colors, markers = 'rb', '^o'
ax0 = plt.subplot(gs[0])
for pressure, color, marker in zip(pressures, colors, markers):
plt.semilogy(1000./temperatures,
original_idt_by_pressure[pressure],
color=color,
marker=marker,
linestyle='-',
label='{0} atm'.format(pressure)
)
plt.semilogy(1000./temperatures,
roaming_idt_by_pressure[pressure],
color=color,
marker=None,
linestyle=':',
label=None)
plt.legend(loc='best', frameon=False)
plt.xlabel("1000K / temperature")
plt.ylabel("Ignition delay time(ms)")
plt.ylim(1e-4,2000)
ax1 = plt.subplot(gs[1])
max_change = 0
for pressure, color, marker in zip(pressures, colors, markers):
changes = (roaming_idt_by_pressure[pressure]-original_idt_by_pressure[pressure]) / original_idt_by_pressure[pressure]
ax1.plot(1000./temperatures,
100.*changes,
color=color,
linestyle=':',
marker=None,
)
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 <= 12:
plt.ylim(-0.2,12)
plt.yticks([0,5,10])
plt.tight_layout()
mprint("## Original model ($-$) vs 0% roaming ($\cdots$)\n(should be identical)")
plot_comparison(temperatures, original_idt_by_pressure, roaming00_idt_by_pressure)
plt.savefig('ignition-heptane-00pc.pdf')
# See how much change when ALPHA = 0.1 (reasonable upper limit)
roaming10_idt_by_pressure = get_roaming_delays(alpha=0.1)
# Redo with alpha=0.1 but phi = 0.5
original_idt_by_pressure_phi05 = get_original_delays(phi=0.5)
roaming10_idt_by_pressure_phi05 = get_roaming_delays(alpha=0.1, phi=0.5)
# Redo with alpha=0.1 but phi = 2.0
original_idt_by_pressure_phi20 = get_original_delays(phi=2.0)
roaming10_idt_by_pressure_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_pressure, original_idt_by_pressure_phi05)
mprint("## 10% roaming with $\phi = 0.5$")
plot_comparison(temperatures, original_idt_by_pressure_phi05, roaming10_idt_by_pressure_phi05)
plt.savefig('ignition-heptane-10pc-phi05.pdf')
mprint("## 10% roaming with $\phi = 1$")
plot_comparison(temperatures, original_idt_by_pressure, roaming10_idt_by_pressure)
plt.savefig('ignition-heptane-10pc.pdf')
mprint("## 10% roaming with $\phi = 2$")
plot_comparison(temperatures, original_idt_by_pressure_phi20, roaming10_idt_by_pressure_phi20)
plt.savefig('ignition-heptane-10pc-phi20.pdf')
# See how much change when ALPHA = 0.05 (moderate?)
roaming05_idt_by_pressure = get_roaming_delays(alpha=0.05)
# See how much change when ALPHA = 0.5 (excessive)
roaming50_idt_by_pressure = get_roaming_delays(alpha=0.5)
mprint("## 5% roaming ($\phi=1$)")
plot_comparison(temperatures, original_idt_by_pressure, roaming05_idt_by_pressure)
plt.savefig('ignition-heptane-05pc.pdf')
mprint("## 50% roaming ($\phi=1$)")
plot_comparison(temperatures, original_idt_by_pressure, roaming50_idt_by_pressure)
plt.savefig('ignition-heptane-50pc.pdf')