1# -*- coding: utf-8 -*- 2""" 3Simulation of a (gaseous) Diesel-type internal combustion engine. 4 5The simulation uses n-Dodecane as fuel, which is injected close to top dead 6center. Note that this example uses numerous simplifying assumptions and 7thus serves for illustration purposes only. 8 9Requires: cantera >= 2.5.0, scipy >= 0.19, matplotlib >= 2.0 10""" 11 12import cantera as ct 13import numpy as np 14 15from scipy.integrate import trapz 16import matplotlib.pyplot as plt 17 18######################################################################### 19# Input Parameters 20######################################################################### 21 22# reaction mechanism, kinetics type and compositions 23reaction_mechanism = 'nDodecane_Reitz.yaml' 24phase_name = 'nDodecane_IG' 25comp_air = 'o2:1, n2:3.76' 26comp_fuel = 'c12h26:1' 27 28f = 3000. / 60. # engine speed [1/s] (3000 rpm) 29V_H = .5e-3 # displaced volume [m**3] 30epsilon = 20. # compression ratio [-] 31d_piston = 0.083 # piston diameter [m] 32 33# turbocharger temperature, pressure, and composition 34T_inlet = 300. # K 35p_inlet = 1.3e5 # Pa 36comp_inlet = comp_air 37 38# outlet pressure 39p_outlet = 1.2e5 # Pa 40 41# fuel properties (gaseous!) 42T_injector = 300. # K 43p_injector = 1600e5 # Pa 44comp_injector = comp_fuel 45 46# ambient properties 47T_ambient = 300. # K 48p_ambient = 1e5 # Pa 49comp_ambient = comp_air 50 51# Inlet valve friction coefficient, open and close timings 52inlet_valve_coeff = 1.e-6 53inlet_open = -18. / 180. * np.pi 54inlet_close = 198. / 180. * np.pi 55 56# Outlet valve friction coefficient, open and close timings 57outlet_valve_coeff = 1.e-6 58outlet_open = 522. / 180 * np.pi 59outlet_close = 18. / 180. * np.pi 60 61# Fuel mass, injector open and close timings 62injector_open = 350. / 180. * np.pi 63injector_close = 365. / 180. * np.pi 64injector_mass = 3.2e-5 # kg 65 66# Simulation time and parameters 67sim_n_revolutions = 8 68delta_T_max = 20. 69rtol = 1.e-12 70atol = 1.e-16 71 72##################################################################### 73# Set up IC engine Parameters and Functions 74##################################################################### 75 76V_oT = V_H / (epsilon - 1.) 77A_piston = .25 * np.pi * d_piston ** 2 78stroke = V_H / A_piston 79 80 81def crank_angle(t): 82 """Convert time to crank angle""" 83 return np.remainder(2 * np.pi * f * t, 4 * np.pi) 84 85 86def piston_speed(t): 87 """Approximate piston speed with sinusoidal velocity profile""" 88 return - stroke / 2 * 2 * np.pi * f * np.sin(crank_angle(t)) 89 90 91##################################################################### 92# Set up Reactor Network 93##################################################################### 94 95# load reaction mechanism 96gas = ct.Solution(reaction_mechanism, phase_name) 97 98# define initial state and set up reactor 99gas.TPX = T_inlet, p_inlet, comp_inlet 100cyl = ct.IdealGasReactor(gas) 101cyl.volume = V_oT 102 103# define inlet state 104gas.TPX = T_inlet, p_inlet, comp_inlet 105inlet = ct.Reservoir(gas) 106 107# inlet valve 108inlet_valve = ct.Valve(inlet, cyl) 109inlet_delta = np.mod(inlet_close - inlet_open, 4 * np.pi) 110inlet_valve.valve_coeff = inlet_valve_coeff 111inlet_valve.set_time_function( 112 lambda t: np.mod(crank_angle(t) - inlet_open, 4 * np.pi) < inlet_delta) 113 114# define injector state (gaseous!) 115gas.TPX = T_injector, p_injector, comp_injector 116injector = ct.Reservoir(gas) 117 118# injector is modeled as a mass flow controller 119injector_mfc = ct.MassFlowController(injector, cyl) 120injector_delta = np.mod(injector_close - injector_open, 4 * np.pi) 121injector_t_open = (injector_close - injector_open) / 2. / np.pi / f 122injector_mfc.mass_flow_coeff = injector_mass / injector_t_open 123injector_mfc.set_time_function( 124 lambda t: np.mod(crank_angle(t) - injector_open, 4 * np.pi) < injector_delta) 125 126# define outlet pressure (temperature and composition don't matter) 127gas.TPX = T_ambient, p_outlet, comp_ambient 128outlet = ct.Reservoir(gas) 129 130# outlet valve 131outlet_valve = ct.Valve(cyl, outlet) 132outlet_delta = np.mod(outlet_close - outlet_open, 4 * np.pi) 133outlet_valve.valve_coeff = outlet_valve_coeff 134outlet_valve.set_time_function( 135 lambda t: np.mod(crank_angle(t) - outlet_open, 4 * np.pi) < outlet_delta) 136 137# define ambient pressure (temperature and composition don't matter) 138gas.TPX = T_ambient, p_ambient, comp_ambient 139ambient_air = ct.Reservoir(gas) 140 141# piston is modeled as a moving wall 142piston = ct.Wall(ambient_air, cyl) 143piston.area = A_piston 144piston.set_velocity(piston_speed) 145 146# create a reactor network containing the cylinder and limit advance step 147sim = ct.ReactorNet([cyl]) 148sim.rtol, sim.atol = rtol, atol 149cyl.set_advance_limit('temperature', delta_T_max) 150 151##################################################################### 152# Run Simulation 153##################################################################### 154 155# set up output data arrays 156states = ct.SolutionArray( 157 cyl.thermo, 158 extra=('t', 'ca', 'V', 'm', 'mdot_in', 'mdot_out', 'dWv_dt'), 159) 160 161# simulate with a maximum resolution of 1 deg crank angle 162dt = 1. / (360 * f) 163t_stop = sim_n_revolutions / f 164while sim.time < t_stop: 165 166 # perform time integration 167 sim.advance(sim.time + dt) 168 169 # calculate results to be stored 170 dWv_dt = - (cyl.thermo.P - ambient_air.thermo.P) * A_piston * \ 171 piston_speed(sim.time) 172 173 # append output data 174 states.append(cyl.thermo.state, 175 t=sim.time, ca=crank_angle(sim.time), 176 V=cyl.volume, m=cyl.mass, 177 mdot_in=inlet_valve.mass_flow_rate, 178 mdot_out=outlet_valve.mass_flow_rate, 179 dWv_dt=dWv_dt) 180 181 182####################################################################### 183# Plot Results in matplotlib 184####################################################################### 185 186def ca_ticks(t): 187 """Helper function converts time to rounded crank angle.""" 188 return np.round(crank_angle(t) * 180 / np.pi, decimals=1) 189 190 191t = states.t 192 193# pressure and temperature 194xticks = np.arange(0, 0.18, 0.02) 195fig, ax = plt.subplots(nrows=2) 196ax[0].plot(t, states.P / 1.e5) 197ax[0].set_ylabel('$p$ [bar]') 198ax[0].set_xlabel(r'$\phi$ [deg]') 199ax[0].set_xticklabels([]) 200ax[1].plot(t, states.T) 201ax[1].set_ylabel('$T$ [K]') 202ax[1].set_xlabel(r'$\phi$ [deg]') 203ax[1].set_xticks(xticks) 204ax[1].set_xticklabels(ca_ticks(xticks)) 205plt.show() 206 207# p-V diagram 208fig, ax = plt.subplots() 209ax.plot(states.V[t > 0.04] * 1000, states.P[t > 0.04] / 1.e5) 210ax.set_xlabel('$V$ [l]') 211ax.set_ylabel('$p$ [bar]') 212plt.show() 213 214# T-S diagram 215fig, ax = plt.subplots() 216ax.plot(states.m[t > 0.04] * states.s[t > 0.04], states.T[t > 0.04]) 217ax.set_xlabel('$S$ [J/K]') 218ax.set_ylabel('$T$ [K]') 219plt.show() 220 221# heat of reaction and expansion work 222fig, ax = plt.subplots() 223ax.plot(t, 1.e-3 * states.heat_release_rate * states.V, label=r'$\dot{Q}$') 224ax.plot(t, 1.e-3 * states.dWv_dt, label=r'$\dot{W}_v$') 225ax.set_ylim(-1e2, 1e3) 226ax.legend(loc=0) 227ax.set_ylabel('[kW]') 228ax.set_xlabel(r'$\phi$ [deg]') 229ax.set_xticks(xticks) 230ax.set_xticklabels(ca_ticks(xticks)) 231plt.show() 232 233# gas composition 234fig, ax = plt.subplots() 235ax.plot(t, states('o2').X, label='O2') 236ax.plot(t, states('co2').X, label='CO2') 237ax.plot(t, states('co').X, label='CO') 238ax.plot(t, states('c12h26').X * 10, label='n-Dodecane x10') 239ax.legend(loc=0) 240ax.set_ylabel('$X_i$ [-]') 241ax.set_xlabel(r'$\phi$ [deg]') 242ax.set_xticks(xticks) 243ax.set_xticklabels(ca_ticks(xticks)) 244plt.show() 245 246###################################################################### 247# Integral Results 248###################################################################### 249 250# heat release 251Q = trapz(states.heat_release_rate * states.V, t) 252output_str = '{:45s}{:>4.1f} {}' 253print(output_str.format('Heat release rate per cylinder (estimate):', 254 Q / t[-1] / 1000., 'kW')) 255 256# expansion power 257W = trapz(states.dWv_dt, t) 258print(output_str.format('Expansion power per cylinder (estimate):', 259 W / t[-1] / 1000., 'kW')) 260 261# efficiency 262eta = W / Q 263print(output_str.format('Efficiency (estimate):', eta * 100., '%')) 264 265# CO emissions 266MW = states.mean_molecular_weight 267CO_emission = trapz(MW * states.mdot_out * states('CO').X[:, 0], t) 268CO_emission /= trapz(MW * states.mdot_out, t) 269print(output_str.format('CO emission (estimate):', CO_emission * 1.e6, 'ppm')) 270