research/cantera
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity Actions

[Python/Examples] Update examples to use new SolutionArray features

Browse source
This commit is contained in:
Ray Speth 2016-07-01 23:02:42 -04:00
parent 8ffbbea93f
commit faa4c36e1d
5 changed files with 34 additions and 59 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

10
interfaces/cython/cantera/examples/reactors/custom.py
View file

@ -52,23 +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]
states = ct.SolutionArray(gas, 1)
states = ct.SolutionArray(gas, 1, extra={'t': [0.0]})
dt = 1e-5
while solver.successful() and solver.t < t_end:
solver.integrate(solver.t + dt)
t_out.append(solver.t)
gas.TPY = solver.y[0], P, solver.y[1:]
states.append(gas.state)
states.append(gas.state, t=solver.t)
# Plot the results
try:
import matplotlib.pyplot as plt
L1 = plt.plot(t_out, states.T, color='r', label='T', lw=2)
L1 = plt.plot(states.t, states.T, color='r', label='T', lw=2)
plt.xlabel('time (s)')
plt.ylabel('Temperature (K)')
plt.twinx()
L2 = plt.plot(t_out, states.Y[:,gas.species_index('OH')], label='OH', lw=2)
L2 = plt.plot(states.t, states('OH').Y, 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()

11
interfaces/cython/cantera/examples/reactors/ic_engine.py
View file

@ -223,10 +223,10 @@ plt.savefig('ic_engine_Q_W.png')
# gas composition
plt.figure()
plt.clf()
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.plot(t, states('O2').X, label='O2')
plt.plot(t, states('CO2').X, label='CO2')
plt.plot(t, states('CO').X, label='CO')
plt.plot(t, states('C3H8').X * 10, label='C3H8 x10')
plt.legend(loc=0)
plt.ylabel('$X_i$ [-]')
plt.xlabel('$\phi$ [deg]')
@ -245,8 +245,7 @@ Q = trapz(heat_release_rate, t)
W = trapz(d_W_v_d_t, t)
eta = W / Q
MW = states.mean_molecular_weight
CO_emission = trapz(MW * mdot_out * states.X[:, gas.species_index('CO')], t) \
/ trapz(MW * mdot_out, t)
CO_emission = trapz(MW * mdot_out * states('CO').X, t) / trapz(MW * mdot_out, t)
print('Heat release rate per cylinder (estimate):\t' +
format(Q / t_sim / 1000., ' 2.1f') + ' kW')
print('Expansion power per cylinder (estimate):\t' +

9
interfaces/cython/cantera/examples/reactors/periodic_cstr.py
View file

@ -75,21 +75,18 @@ network = ct.ReactorNet([cstr])
t = 0.0
dt = 0.1
tm = []
y = []
states = ct.SolutionArray(gas, extra=['t'])
while t < 300.0:
t += dt
network.advance(t)
tm.append(t)
y.append(cstr.thermo['H2','O2','H2O'].Y)
states.append(cstr.thermo.state, t=t)
if __name__ == '__main__':
print(__doc__)
try:
import matplotlib.pyplot as plt
plt.figure(1)
plt.plot(tm, y)
plt.legend(['H2','O2','H2O'])
plt.plot(states.t, states('H2','O2','H2O').Y)
plt.title('Mass Fractions')
plt.show()
except ImportError:

26
interfaces/cython/cantera/examples/reactors/piston.py
View file

@ -48,12 +48,8 @@ w = ct.Wall(r1, r2, velocity=v)
net = ct.ReactorNet([r1, r2])
tim = []
v1 = []
v2 = []
v = []
states1 = ct.SolutionArray(r1.thermo)
states2 = ct.SolutionArray(r2.thermo)
states1 = ct.SolutionArray(r1.thermo, extra=['t','v'])
states2 = ct.SolutionArray(r2.thermo, extra=['t','v'])
for n in range(200):
time = (n+1)*0.001
@ -62,31 +58,27 @@ for n in range(200):
print(fmt % (time, r1.T, r2.T, r1.volume, r2.volume,
r1.volume + r2.volume, r2.thermo['CO'].X[0]))
tim.append(time * 1000)
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)
states1.append(r1.thermo.state, t=1000*time, v=r1.volume)
states2.append(r2.thermo.state, t=1000*time, v=r2.volume)
# 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, states1.T, '-', tim, states2.T, 'r-')
plt.plot(states1.t, states1.T, '-', states2.t, states2.T, 'r-')
plt.xlabel('Time (ms)')
plt.ylabel('Temperature (K)')
plt.subplot(2,2,2)
plt.plot(tim,v1,'-',tim,v2,'r-',tim,v,'g-')
plt.plot(states1.t, states1.v,'-', states2.t, states2.v, 'r-',
states1.t, states1.v + states2.v, 'g-')
plt.xlabel('Time (ms)')
plt.ylabel('Volume (m3)')
plt.subplot(2,2,3)
plt.plot(tim, states2.X[:,states2.species_index('CO')])
plt.plot(states2.t, states2('CO').X)
plt.xlabel('Time (ms)')
plt.ylabel('CO Mole Fraction (right)')
plt.subplot(2,2,4)
plt.plot(tim, states1.X[:,states1.species_index('H2')])
plt.plot(states1.t, states1('H2').X)
plt.xlabel('Time (ms)')
plt.ylabel('H2 Mole Fraction (left)')
plt.tight_layout()

37
interfaces/cython/cantera/examples/reactors/sensitivity1.py
View file

@ -26,52 +26,41 @@ sim.atol = 1.0e-15
sim.rtol_sensitivity = 1.0e-6
sim.atol_sensitivity = 1.0e-6
states = ct.SolutionArray(gas, extra=['t','s2','s3'])
n_times = 400
tim = np.zeros(n_times)
data = np.zeros((n_times,6))
time = 0.0
for n in range(n_times):
time += 5.0e-6
sim.advance(time)
tim[n] = 1000 * time
data[n,0] = r.T
data[n,1:4] = r.thermo['OH','H','CH4'].X
# sensitivity of OH to reaction 2
data[n,4] = sim.sensitivity('OH',2)
# sensitivity of OH to reaction 3
data[n,5] = sim.sensitivity('OH',3)
for t in np.arange(0, 2e-3, 5e-6):
sim.advance(t)
s2 = sim.sensitivity('OH', 2) # sensitivity of OH to reaction 2
s3 = sim.sensitivity('OH', 3) # sensitivity of OH to reaction 3
states.append(r.thermo.state, t=1000*t, s2=s2, s3=s3)
print('%10.3e %10.3f %10.3f %14.6e %10.3f %10.3f' %
(sim.time, r.T, r.thermo.P, r.thermo.u, data[n,4], data[n,5]))
(sim.time, r.T, r.thermo.P, r.thermo.u, s2, s3))
# plot the results if matplotlib is installed.
# see http://matplotlib.org/ to get it
if '--plot' in sys.argv:
import matplotlib.pyplot as plt
plt.subplot(2,2,1)
plt.plot(tim,data[:,0])
plt.plot(states.t, states.T)
plt.xlabel('Time (ms)')
plt.ylabel('Temperature (K)')
plt.subplot(2,2,2)
plt.plot(tim,data[:,1])
plt.plot(states.t, states('OH').X)
plt.xlabel('Time (ms)')
plt.ylabel('OH Mole Fraction')
plt.subplot(2,2,3)
plt.plot(tim,data[:,2])
plt.plot(states.t, states('H').X)
plt.xlabel('Time (ms)')
plt.ylabel('H Mole Fraction')
plt.subplot(2,2,4)
plt.plot(tim,data[:,3])
plt.plot(states.t, states('CH4').X)
plt.xlabel('Time (ms)')
plt.ylabel('H2 Mole Fraction')
plt.ylabel('CH4 Mole Fraction')
plt.tight_layout()
plt.figure(2)
plt.plot(tim,data[:,4],'-',tim,data[:,5],'-g')
plt.plot(states.t, states.s2, '-', states.t, states.s3, '-g')
plt.legend([sim.sensitivity_parameter_name(2),sim.sensitivity_parameter_name(3)],'best')
plt.xlabel('Time (ms)')
plt.ylabel('OH Sensitivity')