Update ic_engine example using new cantera capabilities

interfaces/cython/cantera/examples/reactors/ic_engine.py

@@ -2,27 +2,36 @@
 """
 Simulation of a (gaseous) Diesel-type internal combustion engine.
 
-The use of pure propane as fuel requires an unrealistically high compression
+The simulation uses n-Dodecane as fuel, which is injected close to top dead
-ratio.
+center. Note that this example uses numerous simplifying assumptions and
-
+thus serves for illustration purposes only.
 """
 
 import cantera as ct
 import numpy as np
 
-#######################################################################
+from scipy.integrate import trapz
+import matplotlib.pyplot as plt
+
+#########################################################################
 # Input Parameters
-#######################################################################
+#########################################################################
+
+# reaction mechanism, kinetics type and compositions
+reaction_mechanism = 'nDodecane_Reitz.xml'
+phase_name = 'nDodecane_IG'
+comp_air = 'o2:1, n2:3.76'
+comp_fuel = 'c12h26:1'
 
 f = 3000. / 60.  # engine speed [1/s] (3000 rpm)
 V_H = .5e-3  # displaced volume [m**3]
-epsilon = 50.  # compression ratio [-]
+epsilon = 20.  # compression ratio [-]
 d_piston = 0.083  # piston diameter [m]
 
 # turbocharger temperature, pressure, and composition
 T_inlet = 300.  # K
 p_inlet = 1.3e5  # Pa
-comp_inlet = 'O2:1, N2:3.76'
+comp_inlet = comp_air
 
 # outlet pressure
 p_outlet = 1.2e5  # Pa
@@ -30,15 +39,12 @@ p_outlet = 1.2e5  # Pa
 # fuel properties (gaseous!)
 T_injector = 300.  # K
 p_injector = 1600e5  # Pa
-comp_injector = 'C3H8:1'
+comp_injector = comp_fuel
 
 # ambient properties
 T_ambient = 300.  # K
 p_ambient = 1e5  # Pa
-comp_ambient = 'O2:1, N2:3.76'
+comp_ambient = comp_air
-
-# Reaction mechanism name
-reaction_mechanism = 'gri30.xml'
 
 # Inlet valve friction coefficient, open and close timings
 inlet_valve_coeff = 1.e-6
@@ -54,200 +60,203 @@ outlet_close = 18. / 180. * np.pi
 injector_open = 350. / 180. * np.pi
 injector_close = 365. / 180. * np.pi
 injector_mass = 3.2e-5  # kg
-injector_t_open = (injector_close - injector_open) / 2. / np.pi / f
 
-# Simulation time and resolution
+# Simulation time and parameters
-sim_n_revolutions = 8.
+sim_n_revolutions = 8
-sim_n_timesteps = 100000.
+delta_T_max = 20.
+rtol = 1.e-12
+atol = 1.e-16
 
-###################################################################
+#####################################################################
+# Set up IC engine Parameters and Functions
+#####################################################################
+
+V_oT = V_H / (epsilon - 1.)
+A_piston = .25 * np.pi * d_piston ** 2
+stroke = V_H / A_piston
+
+def crank_angle(t):
+    """Convert time to crank angle"""
+    return np.remainder(2 * np.pi * f * t, 4 * np.pi)
+
+def piston_speed(t):
+    """Approximate piston speed with sinusoidal velocity profile"""
+    return - stroke / 2 * 2 * np.pi * f * np.sin(crank_angle(t))
+
+#####################################################################
+# Set up Reactor Network
+#####################################################################
 
 # load reaction mechanism
-gas = ct.Solution(reaction_mechanism)
+gas = ct.Solution(reaction_mechanism, phase_name)
 
-# define initial state
+# define initial state and set up reactor
 gas.TPX = T_inlet, p_inlet, comp_inlet
-r = ct.IdealGasReactor(gas)
+cyl = ct.IdealGasReactor(gas)
+cyl.volume = V_oT
+
 # define inlet state
 gas.TPX = T_inlet, p_inlet, comp_inlet
 inlet = ct.Reservoir(gas)
+
+# inlet valve
+inlet_valve = ct.Valve(inlet, cyl)
+inlet_delta = np.mod(inlet_close - inlet_open, 4 * np.pi)
+inlet_valve.valve_coeff = inlet_valve_coeff
+inlet_valve.set_time_function(
+    lambda t: np.mod(crank_angle(t) - inlet_open, 4 * np.pi) < inlet_delta)
+
 # define injector state (gaseous!)
 gas.TPX = T_injector, p_injector, comp_injector
 injector = ct.Reservoir(gas)
+
+# injector is modeled as a mass flow controller
+injector_mfc = ct.MassFlowController(injector, cyl)
+injector_delta = np.mod(injector_close - injector_open, 4 * np.pi)
+injector_t_open = (injector_close - injector_open) / 2. / np.pi / f
+injector_mfc.mass_flow_coeff = injector_mass / injector_t_open
+injector_mfc.set_time_function(
+    lambda t: np.mod(crank_angle(t) - injector_open, 4 * np.pi) < injector_delta)
+
 # define outlet pressure (temperature and composition don't matter)
 gas.TPX = T_ambient, p_outlet, comp_ambient
 outlet = ct.Reservoir(gas)
+
+# outlet valve
+outlet_valve = ct.Valve(cyl, outlet)
+outlet_delta = np.mod(outlet_close - outlet_open, 4 * np.pi)
+outlet_valve.valve_coeff = outlet_valve_coeff
+outlet_valve.set_time_function(
+    lambda t: np.mod(crank_angle(t) - outlet_open, 4 * np.pi) < outlet_delta)
+
 # define ambient pressure (temperature and composition don't matter)
 gas.TPX = T_ambient, p_ambient, comp_ambient
 ambient_air = ct.Reservoir(gas)
 
-# set up connecting devices
+# piston is modeled as a moving wall
-inlet_valve = ct.Valve(inlet, r)
+piston = ct.Wall(ambient_air, cyl)
-injector_mfc = ct.MassFlowController(injector, r)
-outlet_valve = ct.Valve(r, outlet)
-piston = ct.Wall(ambient_air, r)
-
-# convert time to crank angle
-def crank_angle(t):
-    return np.remainder(2 * np.pi * f * t, 4 * np.pi)
-
-# set up IC engine parameters
-V_oT = V_H / (epsilon - 1.)
-A_piston = .25 * np.pi * d_piston ** 2
-stroke = V_H / A_piston
-r.volume = V_oT
 piston.area = A_piston
-def piston_speed(t):
-    return - stroke / 2 * 2 * np.pi * f * np.sin(crank_angle(t))
 piston.set_velocity(piston_speed)
 
-# create a reactor network containing the cylinder
+# create a reactor network containing the cylinder and limit advance step
-sim = ct.ReactorNet([r])
+sim = ct.ReactorNet([cyl])
+sim.rtol, sim.atol = rtol, atol
+cyl.set_advance_limit('temperature', delta_T_max)
+
+#####################################################################
+# Run Simulation
+#####################################################################
 
 # set up output data arrays
-states = ct.SolutionArray(r.thermo)
+states = ct.SolutionArray(cyl.thermo,
-t_sim = sim_n_revolutions / f
+                          extra=('t', 'ca', 'V', 'm',
-t = (np.arange(sim_n_timesteps) + 1) / sim_n_timesteps * t_sim
+                                 'mdot_in', 'mdot_out',
-V = np.zeros_like(t)
+                                 'dWv_dt', 'heat_release_rate'))
-m = np.zeros_like(t)
-test = np.zeros_like(t)
-mdot_in = np.zeros_like(t)
-mdot_out = np.zeros_like(t)
-d_W_v_d_t = np.zeros_like(t)
-heat_release_rate = np.zeros_like(t)
 
-# set parameters for the automatic time step refinement
+# simulate with a maximum resolution of 1 deg crank angle
-n_last_refinement = -np.inf  # for initialization only
+dt = 1. / (360 * f)
-n_wait_coarsening = 10
+t_stop = sim_n_revolutions / f
+while sim.time < t_stop:
 
-# do simulation
+    # perform time integration
-for n1, t_i in enumerate(t):
+    sim.advance(sim.time + dt)
-    # define opening and closing of valves and injector
-    if (np.mod(crank_angle(t_i) - inlet_open, 4 * np.pi) <
-            np.mod(inlet_close - inlet_open, 4 * np.pi)):
-        inlet_valve.valve_coeff = inlet_valve_coeff
-        test[n1] = 1
-    else:
-        inlet_valve.valve_coeff = 0.
-    if (np.mod(crank_angle(t_i) - outlet_open, 4 * np.pi) <
-            np.mod(outlet_close - outlet_open, 4 * np.pi)):
-        outlet_valve.valve_coeff = outlet_valve_coeff
-    else:
-        outlet_valve.valve_coeff = 0.
-    if (np.mod(crank_angle(t_i) - injector_open, 4 * np.pi) <
-            np.mod(injector_close - injector_open, 4 * np.pi)):
-        injector_mfc.set_mass_flow_rate(injector_mass / injector_t_open)
-    else:
-        injector_mfc.set_mass_flow_rate(0)
 
-    # perform time integration, refine time step if necessary
+    # calculate results to be stored
-    solved = False
+    dWv_dt = - (cyl.thermo.P - ambient_air.thermo.P) * A_piston * \
-    for n2 in range(4):
+        piston_speed(sim.time)
-        try:
+    heat_release_rate = - cyl.volume * ct.gas_constant * cyl.T * \
-            sim.advance(t_i)
+        np.sum(gas.standard_enthalpies_RT * cyl.thermo.net_production_rates, 0)
-            solved = True
-            break
-        except ct.CanteraError:
-            sim.set_max_time_step(1e-6 * 10. ** -n2)
-            n_last_refinement = n1
-    if not solved:
-        raise ct.CanteraError('Refinement limit reached')
-    # coarsen time step if too long ago
-    if n1 - n_last_refinement == n_wait_coarsening:
-        sim.set_max_time_step(1e-5)
 
-    # write output data
+    # append output data
-    states.append(r.thermo.state)
+    states.append(cyl.thermo.state,
-    V[n1] = r.volume
+                  t=sim.time, ca=crank_angle(sim.time),
-    m[n1] = r.mass
+                  V=cyl.volume, m=cyl.mass,
-    mdot_in[n1] = inlet_valve.mdot(0)
+                  mdot_in=inlet_valve.mdot(sim.time),
-    mdot_out[n1] = outlet_valve.mdot(0)
+                  mdot_out=outlet_valve.mdot(sim.time),
-    d_W_v_d_t[n1] = - (r.thermo.P - ambient_air.thermo.P) * A_piston * \
+                  dWv_dt=dWv_dt,
-        piston_speed(t_i)
+                  heat_release_rate=heat_release_rate)
-    heat_release_rate[n1] = - r.volume * ct.gas_constant * r.T * \
-        np.sum(gas.standard_enthalpies_RT * r.thermo.net_production_rates, 0)
 
-
+#######################################################################
-#####################################################################
 # Plot Results in matplotlib
-#####################################################################
+#######################################################################
 
-import matplotlib.pyplot as plt
+def ca_ticks(t):
+    """Helper function converts time to rounded crank angle."""
+    return np.round(crank_angle(t) * 180 / np.pi, decimals=1)
+
+t = states.t
 
 # pressure and temperature
-plt.clf()
+fig, ax = plt.subplots(nrows=2)
-plt.subplot(211)
+ax[0].plot(t, states.P / 1.e5)
-plt.plot(t, states.P / 1.e5)
+ax[0].set_ylabel('$p$ [bar]')
-plt.ylabel('$p$ [bar]')
+ax[0].set_xlabel('$\phi$ [deg]')
-plt.xlabel('$\phi$ [deg]')
+ax[0].set_xticklabels([])
-plt.xticks(plt.xticks()[0], [])
+ax[1].plot(t, states.T)
-plt.subplot(212)
+ax[1].set_ylabel('$T$ [K]')
-plt.plot(t, states.T)
+ax[1].set_xlabel('$\phi$ [deg]')
-plt.ylabel('$T$ [K]')
+ax[1].set_xticklabels(ca_ticks(ax[1].get_xticks()))
-plt.xlabel('$\phi$ [deg]')
-plt.xticks(plt.xticks()[0], crank_angle(plt.xticks()[0]) * 180 / np.pi,
-           rotation=17)
 plt.show()
-plt.savefig('ic_engine_t_p_T.png')
 
 # p-V diagram
-plt.clf()
+fig, ax = plt.subplots()
-plt.plot(V[t > 0.04] * 1000, states.P[t > 0.04] / 1.e5)
+ax.plot(states.V[t > 0.04] * 1000, states.P[t > 0.04] / 1.e5)
-plt.xlabel('$V$ [l]')
+ax.set_xlabel('$V$ [l]')
-plt.ylabel('$p$ [bar]')
+ax.set_ylabel('$p$ [bar]')
 plt.show()
-plt.savefig('ic_engine_p_V.png')
 
 # T-S diagram
-plt.clf()
+fig, ax = plt.subplots()
-plt.plot(m[t > 0.04] * states.s[t > 0.04], states.T[t > 0.04])
+ax.plot(states.m[t > 0.04] * states.s[t > 0.04], states.T[t > 0.04])
-plt.xlabel('$S$ [J/K]')
+ax.set_xlabel('$S$ [J/K]')
-plt.ylabel('$T$ [K]')
+ax.set_ylabel('$T$ [K]')
 plt.show()
-plt.savefig('ic_engine_T_S.png')
 
 # heat of reaction and expansion work
-plt.clf()
+fig, ax = plt.subplots()
-plt.plot(t, heat_release_rate, label='$\dot{Q}$')
+ax.plot(t, 1.e-3 * states.heat_release_rate, label='$\dot{Q}$')
-plt.plot(t, d_W_v_d_t, label='$\dot{W}_v$')
+ax.plot(t, 1.e-3 * states.dWv_dt, label='$\dot{W}_v$')
-plt.ylim(-1e5, 1e6)
+ax.set_ylim(-1e2, 1e3)
-plt.legend(loc=0)
+ax.legend(loc=0)
-plt.ylabel('[W]')
+ax.set_ylabel('[kW]')
-plt.xlabel('$\phi$ [deg]')
+ax.set_xlabel('$\phi$ [deg]')
-plt.xticks(plt.xticks()[0], crank_angle(plt.xticks()[0]) * 180 / np.pi,
+ax.set_xticklabels(ca_ticks(ax.get_xticks()))
-           rotation=17)
 plt.show()
-plt.savefig('ic_engine_Q_W.png')
 
 # gas composition
-plt.clf()
+fig, ax = plt.subplots()
-plt.plot(t, states('O2').X, label='O2')
+ax.plot(t, states('o2').X, label='O2')
-plt.plot(t, states('CO2').X, label='CO2')
+ax.plot(t, states('co2').X, label='CO2')
-plt.plot(t, states('CO').X, label='CO')
+ax.plot(t, states('co').X, label='CO')
-plt.plot(t, states('C3H8').X * 10, label='C3H8 x10')
+ax.plot(t, states('c12h26').X * 10, label='n-Dodecane x10')
-plt.legend(loc=0)
+ax.legend(loc=0)
-plt.ylabel('$X_i$ [-]')
+ax.set_ylabel('$X_i$ [-]')
-plt.xlabel('$\phi$ [deg]')
+ax.set_xlabel('$\phi$ [deg]')
-plt.xticks(plt.xticks()[0], crank_angle(plt.xticks()[0]) * 180 / np.pi,
+ax.set_xticklabels(ca_ticks(ax.get_xticks()))
-           rotation=17)
 plt.show()
-plt.savefig('ic_engine_t_X.png')
 
-
+######################################################################
-#####################################################################
 # Integral Results
-#####################################################################
+######################################################################
 
-from scipy.integrate import trapz
+# heat release
-Q = trapz(heat_release_rate, t)
+Q = trapz(states.heat_release_rate, t)
-W = trapz(d_W_v_d_t, t)
+print('{:45s}{:>4.1f} kW'.format('Heat release rate per cylinder (estimate):',
+                                 Q / t[-1] / 1000.))
+
+# expansion power
+W = trapz(states.dWv_dt, t)
+print('{:45s}{:>4.1f} kW'.format('Expansion power per cylinder (estimate):',
+                                 W / t[-1] / 1000.))
+
+# efficiency
 eta = W / Q
+print('{:45s}{:>4.1f} %'.format('Efficiency (estimate):',
+                                eta * 100.))
+
+# CO emissions
 MW = states.mean_molecular_weight
-CO_emission = trapz(MW * mdot_out * states('CO').X[:,0], t) / trapz(MW * mdot_out, t)
+CO_emission = trapz(MW * states.mdot_out * states('CO').X[:, 0], t)
-print('Heat release rate per cylinder (estimate):\t' +
+CO_emission /= trapz(MW * states.mdot_out, t)
-      format(Q / t_sim / 1000., ' 2.1f') + ' kW')
+print('{:45s}{:>4.1f} ppm'.format('CO emission (estimate):',
-print('Expansion power per cylinder (estimate):\t' +
+                                  CO_emission * 1.e6))
-      format(W / t_sim / 1000., ' 2.1f') + ' kW')
-print('Efficiency (estimate):\t\t\t' + format(eta * 100., ' 2.1f') + ' %')
-print('CO emission (estimate):\t\t' + format(CO_emission * 1.e6, ' 2.1f') +
-      ' ppm')