[Python/Examples] Update examples to use SolutionArray
This commit is contained in:
Ray Speth
parent
beb2056926
commit
59e0f5a6d9
6 changed files with 46 additions and 67 deletions
|
|
@ -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()
|
||||
|
|
|
|||
|
|
@ -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')
|
||||
|
|
|
|||
|
|
@ -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)
|
||||
|
|
|
|||
|
|
@ -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()
|
||||
|
|
|
|||
|
|
@ -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()
|
||||
|
|
|
|||
|
|
@ -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)')
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue