diff --git a/interfaces/cython/cantera/examples/reactors/custom.py b/interfaces/cython/cantera/examples/reactors/custom.py index 51f3e6c48..d0412e63d 100644 --- a/interfaces/cython/cantera/examples/reactors/custom.py +++ b/interfaces/cython/cantera/examples/reactors/custom.py @@ -52,25 +52,21 @@ solver.set_initial_value(y0, 0.0) # Integrate the equations, keeping T(t) and Y(k,t) t_end = 1e-3 t_out = [0.0] -T_out = [gas.T] -Y_out = [gas.Y] +states = ct.SolutionArray(gas, 1) dt = 1e-5 while solver.successful() and solver.t < t_end: solver.integrate(solver.t + dt) t_out.append(solver.t) - T_out.append(gas.T) - Y_out.append(gas.Y) - -Y_out = np.array(Y_out).T + states.append(gas.state) # Plot the results try: import matplotlib.pyplot as plt - L1 = plt.plot(t_out, T_out, color='r', label='T', lw=2) + L1 = plt.plot(t_out, states.T, color='r', label='T', lw=2) plt.xlabel('time (s)') plt.ylabel('Temperature (K)') plt.twinx() - L2 = plt.plot(t_out, Y_out[gas.species_index('OH')], label='OH', lw=2) + L2 = plt.plot(t_out, states.Y[:,gas.species_index('OH')], label='OH', lw=2) plt.ylabel('Mass Fraction') plt.legend(L1+L2, [line.get_label() for line in L1+L2], loc='lower right') plt.show() diff --git a/interfaces/cython/cantera/examples/reactors/ic_engine.py b/interfaces/cython/cantera/examples/reactors/ic_engine.py index 6ad369549..315ae5f75 100644 --- a/interfaces/cython/cantera/examples/reactors/ic_engine.py +++ b/interfaces/cython/cantera/examples/reactors/ic_engine.py @@ -105,20 +105,16 @@ piston.set_velocity(piston_speed) sim = ct.ReactorNet([r]) # set up output data arrays +states = ct.SolutionArray(r.thermo) t_sim = sim_n_revolutions / f t = (np.arange(sim_n_timesteps) + 1) / sim_n_timesteps * t_sim -p = np.zeros_like(t) V = np.zeros_like(t) -T = np.zeros_like(t) -s = np.zeros_like(t) m = np.zeros_like(t) test = np.zeros_like(t) mdot_in = np.zeros_like(t) mdot_out = np.zeros_like(t) -MW = np.zeros_like(t) d_W_v_d_t = np.zeros_like(t) heat_release_rate = np.zeros_like(t) -species_X = np.zeros((t.size, gas.n_species)) # set parameters for the automatic time step refinement n_last_refinement = -np.inf # for initialization only @@ -158,15 +154,11 @@ for n1, t_i in enumerate(t): sim.set_max_time_step(1e-5) # write output data - p[n1] = r.thermo.P + states.append(r.thermo.state) V[n1] = r.volume - T[n1] = r.T - s[n1] = r.thermo.s m[n1] = r.mass mdot_in[n1] = inlet_valve.mdot(0) mdot_out[n1] = outlet_valve.mdot(0) - MW[n1] = r.thermo.mean_molecular_weight - species_X[n1] = r.thermo.X d_W_v_d_t[n1] = - (r.thermo.P - ambient_air.thermo.P) * A_piston * \ piston_speed(t_i) heat_release_rate[n1] = - r.volume * ct.gas_constant * r.T * \ @@ -183,12 +175,12 @@ import matplotlib.pyplot as plt plt.figure() plt.clf() plt.subplot(211) -plt.plot(t, p / 1.e5) +plt.plot(t, states.P / 1.e5) plt.ylabel('$p$ [bar]') plt.xlabel('$\phi$ [deg]') plt.xticks(plt.xticks()[0], []) plt.subplot(212) -plt.plot(t, T) +plt.plot(t, states.T) plt.ylabel('$T$ [K]') plt.xlabel('$\phi$ [deg]') plt.xticks(plt.xticks()[0], crank_angle(plt.xticks()[0]) * 180 / np.pi, @@ -199,7 +191,7 @@ plt.savefig('ic_engine_t_p_T.png') # p-V diagram plt.figure() plt.clf() -plt.plot(V[t > 0.04] * 1000, p[t > 0.04] / 1.e5) +plt.plot(V[t > 0.04] * 1000, states.P[t > 0.04] / 1.e5) plt.xlabel('$V$ [l]') plt.ylabel('$p$ [bar]') plt.show() @@ -208,7 +200,7 @@ plt.savefig('ic_engine_p_V.png') # T-S diagram plt.figure() plt.clf() -plt.plot(m[t > 0.04] * s[t > 0.04], T[t > 0.04]) +plt.plot(m[t > 0.04] * states.s[t > 0.04], states.T[t > 0.04]) plt.xlabel('$S$ [J/K]') plt.ylabel('$T$ [K]') plt.show() @@ -231,10 +223,10 @@ plt.savefig('ic_engine_Q_W.png') # gas composition plt.figure() plt.clf() -plt.plot(t, species_X[:, gas.species_index('O2')], label='O2') -plt.plot(t, species_X[:, gas.species_index('CO2')], label='CO2') -plt.plot(t, species_X[:, gas.species_index('CO')], label='CO') -plt.plot(t, species_X[:, gas.species_index('C3H8')] * 10, label='C3H8 x10') +plt.plot(t, states.X[:, gas.species_index('O2')], label='O2') +plt.plot(t, states.X[:, gas.species_index('CO2')], label='CO2') +plt.plot(t, states.X[:, gas.species_index('CO')], label='CO') +plt.plot(t, states.X[:, gas.species_index('C3H8')] * 10, label='C3H8 x10') plt.legend(loc=0) plt.ylabel('$X_i$ [-]') plt.xlabel('$\phi$ [deg]') @@ -252,7 +244,8 @@ from scipy.integrate import trapz Q = trapz(heat_release_rate, t) W = trapz(d_W_v_d_t, t) eta = W / Q -CO_emission = trapz(MW * mdot_out * species_X[:, gas.species_index('CO')], t) \ +MW = states.mean_molecular_weight +CO_emission = trapz(MW * mdot_out * states.X[:, gas.species_index('CO')], t) \ / trapz(MW * mdot_out, t) print('Heat release rate per cylinder (estimate):\t' + format(Q / t_sim / 1000., ' 2.1f') + ' kW') diff --git a/interfaces/cython/cantera/examples/reactors/pfr.py b/interfaces/cython/cantera/examples/reactors/pfr.py index 1ec110009..982da873a 100644 --- a/interfaces/cython/cantera/examples/reactors/pfr.py +++ b/interfaces/cython/cantera/examples/reactors/pfr.py @@ -54,17 +54,14 @@ dt = t_total / n_steps t1 = (np.arange(n_steps) + 1) * dt z1 = np.zeros_like(t1) u1 = np.zeros_like(t1) -T1 = np.zeros_like(t1) -X_H2_1 = np.zeros_like(t1) +states1 = ct.SolutionArray(r1.thermo) for n1, t_i in enumerate(t1): # perform time integration sim1.advance(t_i) # compute velocity and transform into space u1[n1] = mass_flow_rate1 / area / r1.thermo.density z1[n1] = z1[n1 - 1] + u1[n1] * dt - # write output data - T1[n1] = r1.T - X_H2_1[n1] = r1.thermo['H2'].X + states1.append(r1.thermo.state) ##################################################################### @@ -115,10 +112,9 @@ sim2 = ct.ReactorNet([r2]) # define time, space, and other information vectors z2 = (np.arange(n_steps) + 1) * dz t_r2 = np.zeros_like(z2) # residence time in each reactor -t2 = np.zeros_like(z2) u2 = np.zeros_like(z2) -T2 = np.zeros_like(z2) -X_H2_2 = np.zeros_like(z2) +t2 = np.zeros_like(z2) +states2 = ct.SolutionArray(r2.thermo) # iterate through the PFR cells for n in range(n_steps): # Set the state of the reservoir to match that of the previous reactor @@ -132,8 +128,7 @@ for n in range(n_steps): t_r2[n] = r2.mass / mass_flow_rate2 # residence time in this reactor t2[n] = np.sum(t_r2) # write output data - T2[n] = r2.T - X_H2_2[n] = r2.thermo['H2'].X + states2.append(r2.thermo.state) ##################################################################### @@ -145,8 +140,8 @@ for n in range(n_steps): import matplotlib.pyplot as plt plt.figure() -plt.plot(z1, T1, label='Lagrangian Particle') -plt.plot(z2, T2, label='Reactor Chain') +plt.plot(z1, states1.T, label='Lagrangian Particle') +plt.plot(z2, states2.T, label='Reactor Chain') plt.xlabel('$z$ [m]') plt.ylabel('$T$ [K]') plt.legend(loc=0) @@ -154,8 +149,8 @@ plt.show() plt.savefig('pfr_T_z.png') plt.figure() -plt.plot(t1, X_H2_1, label='Lagrangian Particle') -plt.plot(t2, X_H2_2, label='Reactor Chain') +plt.plot(t1, states1.X[:, gas1.species_index('H2')], label='Lagrangian Particle') +plt.plot(t2, states2.X[:, gas2.species_index('H2')], label='Reactor Chain') plt.xlabel('$t$ [s]') plt.ylabel('$X_{H_2}$ [-]') plt.legend(loc=0) diff --git a/interfaces/cython/cantera/examples/reactors/piston.py b/interfaces/cython/cantera/examples/reactors/piston.py index f40824e31..f48b90ea6 100644 --- a/interfaces/cython/cantera/examples/reactors/piston.py +++ b/interfaces/cython/cantera/examples/reactors/piston.py @@ -49,13 +49,11 @@ w = ct.Wall(r1, r2, velocity=v) net = ct.ReactorNet([r1, r2]) tim = [] -t1 = [] -t2 = [] v1 = [] v2 = [] v = [] -xco = [] -xh2 = [] +states1 = ct.SolutionArray(r1.thermo) +states2 = ct.SolutionArray(r2.thermo) for n in range(200): time = (n+1)*0.001 @@ -65,20 +63,18 @@ for n in range(200): r1.volume + r2.volume, r2.thermo['CO'].X[0])) tim.append(time * 1000) - t1.append(r1.T) - t2.append(r2.T) + states1.append(r1.thermo.state, v=r1.volume) + states2.append(r2.thermo.state) v1.append(r1.volume) v2.append(r2.volume) v.append(r1.volume + r2.volume) - xco.append(r2.thermo['CO'].X[0]) - xh2.append(r1.thermo['H2'].X[0]) # plot the results if matplotlib is installed. if '--plot' in sys.argv: import matplotlib.pyplot as plt plt.subplot(2,2,1) - plt.plot(tim,t1,'-',tim,t2,'r-') + plt.plot(tim, states1.T, '-', tim, states2.T, 'r-') plt.xlabel('Time (ms)') plt.ylabel('Temperature (K)') plt.subplot(2,2,2) @@ -86,11 +82,11 @@ if '--plot' in sys.argv: plt.xlabel('Time (ms)') plt.ylabel('Volume (m3)') plt.subplot(2,2,3) - plt.plot(tim,xco) + plt.plot(tim, states2.X[:,states2.species_index('CO')]) plt.xlabel('Time (ms)') plt.ylabel('CO Mole Fraction (right)') plt.subplot(2,2,4) - plt.plot(tim,xh2) + plt.plot(tim, states1.X[:,states1.species_index('H2')]) plt.xlabel('Time (ms)') plt.ylabel('H2 Mole Fraction (left)') plt.tight_layout() diff --git a/interfaces/cython/cantera/examples/reactors/reactor1.py b/interfaces/cython/cantera/examples/reactors/reactor1.py index 91cc6a274..da15962b2 100644 --- a/interfaces/cython/cantera/examples/reactors/reactor1.py +++ b/interfaces/cython/cantera/examples/reactors/reactor1.py @@ -14,15 +14,14 @@ r = ct.IdealGasConstPressureReactor(gri3) sim = ct.ReactorNet([r]) time = 0.0 times = np.zeros(100) -data = np.zeros((100,4)) +states = ct.SolutionArray(gri3, 100) print('%10s %10s %10s %14s' % ('t [s]','T [K]','P [Pa]','u [J/kg]')) for n in range(100): time += 1.e-5 sim.advance(time) times[n] = time * 1e3 # time in ms - data[n,0] = r.T - data[n,1:] = r.thermo['OH','H','H2'].X + states[n].TDY = r.thermo.TDY print('%10.3e %10.3f %10.3f %14.6e' % (sim.time, r.T, r.thermo.P, r.thermo.u)) @@ -32,19 +31,19 @@ if '--plot' in sys.argv[1:]: import matplotlib.pyplot as plt plt.clf() plt.subplot(2, 2, 1) - plt.plot(times, data[:,0]) + plt.plot(times, states.T) plt.xlabel('Time (ms)') plt.ylabel('Temperature (K)') plt.subplot(2, 2, 2) - plt.plot(times, data[:,1]) + plt.plot(times, states.X[:,gri3.species_index('OH')]) plt.xlabel('Time (ms)') plt.ylabel('OH Mole Fraction') plt.subplot(2, 2, 3) - plt.plot(times, data[:,2]) + plt.plot(times, states.X[:,gri3.species_index('H')]) plt.xlabel('Time (ms)') plt.ylabel('H Mole Fraction') plt.subplot(2, 2, 4) - plt.plot(times,data[:,3]) + plt.plot(times, states.X[:,gri3.species_index('H2')]) plt.xlabel('Time (ms)') plt.ylabel('H2 Mole Fraction') plt.tight_layout() diff --git a/interfaces/cython/cantera/examples/reactors/reactor2.py b/interfaces/cython/cantera/examples/reactors/reactor2.py index b96be88c0..c713cdec4 100644 --- a/interfaces/cython/cantera/examples/reactors/reactor2.py +++ b/interfaces/cython/cantera/examples/reactors/reactor2.py @@ -67,8 +67,8 @@ outfile = open('piston.csv', 'w') csvfile = csv.writer(outfile) csvfile.writerow(['time (s)','T1 (K)','P1 (Bar)','V1 (m3)', 'T2 (K)','P2 (Bar)','V2 (m3)']) -temp = np.zeros((n_steps, 2)) -pres = np.zeros((n_steps, 2)) +states1 = ct.SolutionArray(ar) +states2 = ct.SolutionArray(gri3) vol = np.zeros((n_steps, 2)) tm = np.zeros(n_steps) @@ -77,11 +77,11 @@ for n in range(n_steps): print(n, time, r2.T) sim.advance(time) tm[n] = time - temp[n,:] = r1.T, r2.T - pres[n,:] = 1.0e-5*r1.thermo.P, 1.0e-5*r2.thermo.P + states1.append(r1.thermo.state) + states2.append(r2.thermo.state) vol[n,:] = r1.volume, r2.volume - csvfile.writerow([tm[n], temp[n,0], pres[n,0], vol[n,0], - temp[n,1], pres[n,1], vol[n,1]]) + csvfile.writerow([tm[n], r1.thermo.T, r1.thermo.P, vol[n,0], + r2.thermo.T, r2.thermo.P, vol[n,1]]) outfile.close() print('Output written to file piston.csv') print('Directory: '+os.getcwd()) @@ -90,13 +90,13 @@ if '--plot' in sys.argv: import matplotlib.pyplot as plt plt.clf() plt.subplot(2,2,1) - h = plt.plot(tm, temp[:,0],'g-',tm, temp[:,1],'b-') + h = plt.plot(tm, states1.T, 'g-', tm, states2.T, 'b-') #plt.legend(['Reactor 1','Reactor 2'],2) plt.xlabel('Time (s)') plt.ylabel('Temperature (K)') plt.subplot(2,2,2) - plt.plot(tm, pres[:,0],'g-',tm, pres[:,1],'b-') + plt.plot(tm, states1.P / 1e5, 'g-', tm, states2.P / 1e5, 'b-') #plt.legend(['Reactor 1','Reactor 2'],2) plt.xlabel('Time (s)') plt.ylabel('Pressure (Bar)')