Cleanup of upper lvl dir.
This commit is contained in:
Harry Moffat
parent
50e325563c
commit
5b30416449
9 changed files with 13 additions and 761 deletions
|
|
@ -1,7 +1,7 @@
|
|||
#!/bin/sh
|
||||
|
||||
PY_DEMOS = combustor_sim functors_sim mix1_sim mix2_sim piston_sim reactor1_sim \
|
||||
reactor2_sim sensitivity_sim surf_prf_sim
|
||||
reactor2_sim sensitivity_sim surf_pfr_sim
|
||||
|
||||
all:
|
||||
@(for py in $(PY_DEMOS) ; do \
|
||||
|
|
|
|||
|
|
@ -1,84 +0,0 @@
|
|||
# Mixing two streams.
|
||||
|
||||
# Since reactors can have multiple inlets and outlets, they can be
|
||||
# used to implement mixers, splitters, etc. In this example, air and
|
||||
# methane are mixed in stoichiometric proportions. Due to the low
|
||||
# temperature, no reactions occur. Note that the air stream and the
|
||||
# methane stream use *different* reaction mechanisms, with different
|
||||
# numbers of species and reactions. When gas flows from one reactor or
|
||||
# reservoir to another one with a different reaction mechanism,
|
||||
# species are matched by name. If the upstream reactor contains a
|
||||
# species that is not present in the downstream reaction mechanism, it
|
||||
# will be ignored. In general, reaction mechanisms for downstream
|
||||
# reactors should contain all species that might be present in any
|
||||
# upstream reactor.
|
||||
#
|
||||
#-----------------------------------------------------------------------
|
||||
|
||||
from Cantera import *
|
||||
from Cantera.Reactor import *
|
||||
|
||||
|
||||
# Use air for stream a. Note that the Air() function does not set the
|
||||
# composition correctly; thus, we need to explicitly set the
|
||||
# composition to that of air.
|
||||
gas_a = Air()
|
||||
gas_a.set(T = 300.0, P = OneAtm, X = 'O2:0.21, N2:0.78, AR:0.01')
|
||||
rho_a = gas_a.density()
|
||||
|
||||
|
||||
# Use GRI-Mech 3.0 for stream b (methane) and for the mixer. If it is
|
||||
# desired to have a pure mixer, with no chemistry, use instead a
|
||||
# reaction mechanism for gas_b that has no reactions.
|
||||
gas_b = GRI30()
|
||||
gas_b.set(T = 300.0, P = OneAtm, X = 'CH4:1')
|
||||
rho_b = gas_b.density()
|
||||
|
||||
|
||||
# Create reservoirs for the two inlet streams and for the outlet
|
||||
# stream. The upsteam reservoirs could be replaced by reactors, which
|
||||
# might themselves be connected to reactors further upstream. The
|
||||
# outlet reservoir could be replaced with a reactor with no outlet, if
|
||||
# it is desired to integrate the composition leaving the mixer in
|
||||
# time, or by an arbitrary network of downstream reactors.
|
||||
res_a = Reservoir(gas_a)
|
||||
res_b = Reservoir(gas_b)
|
||||
downstream = Reservoir(gas_b)
|
||||
|
||||
|
||||
# Create a reactor for the mixer. A reactor is required instead of a
|
||||
# reservoir, since the state will change with time if the inlet mass
|
||||
# flow rates change or if there is chemistry occurring.
|
||||
mixer = Reactor(gas_b)
|
||||
|
||||
|
||||
# create two mass flow controllers connecting the upstream reservoirs
|
||||
# to the mixer, and set their mass flow rates to values corresponding
|
||||
# to stoichiometric combustion.
|
||||
mfc1 = MassFlowController(upstream = res_a, downstream = mixer,
|
||||
mdot = rho_a*2.5/0.21)
|
||||
|
||||
mfc2 = MassFlowController(upstream = res_b, downstream = mixer,
|
||||
mdot = rho_b*1.0)
|
||||
|
||||
|
||||
# connect the mixer to the downstream reservoir with a valve.
|
||||
outlet = Valve(upstream = mixer, downstream = downstream, Kv = 1.0)
|
||||
|
||||
sim = ReactorNet([mixer])
|
||||
|
||||
# Since the mixer is a reactor, we need to integrate in time to reach
|
||||
# steady state. A few residence times should be enough.
|
||||
t = 0.0
|
||||
for n in range(30):
|
||||
tres = mixer.mass()/(mfc1.massFlowRate() + mfc2.massFlowRate())
|
||||
t += 0.5*tres
|
||||
sim.advance(t)
|
||||
print '%14.5g %14.5g %14.5g %14.5g %14.5g' % (t, mixer.temperature(),
|
||||
mixer.enthalpy_mass(),
|
||||
mixer.pressure(),
|
||||
mixer.massFraction('CH4'))
|
||||
|
||||
# view the state of the gas in the mixer
|
||||
print mixer.contents()
|
||||
|
||||
|
|
@ -1,103 +0,0 @@
|
|||
# Mixing two streams.
|
||||
|
||||
# Since reactors can have multiple inlets and outlets, they can be
|
||||
# used to implement mixers, splitters, etc. In this example, air and
|
||||
# methane are mixed in stoichiometric proportions. Due to the low
|
||||
# temperature, no reactions occur. Note that the air stream and the
|
||||
# methane stream use *different* reaction mechanisms, with different
|
||||
# numbers of species and reactions. When gas flows from one reactor or
|
||||
# reservoir to another one with a different reaction mechanism,
|
||||
# species are matched by name. If the upstream reactor contains a
|
||||
# species that is not present in the downstream reaction mechanism, it
|
||||
# will be ignored. In general, reaction mechanisms for downstream
|
||||
# reactors should contain all species that might be present in any
|
||||
# upstream reactor.
|
||||
#
|
||||
#-----------------------------------------------------------------------
|
||||
|
||||
from Cantera import *
|
||||
from Cantera.Reactor import *
|
||||
|
||||
|
||||
# Use air for stream a. Note that the Air() function does not set the
|
||||
# composition correctly; thus, we need to explicitly set the
|
||||
# composition to that of air.
|
||||
gas_a = Air()
|
||||
gas_a.set(T = 300.0, P = OneAtm, X = 'O2:0.21, N2:0.78, AR:0.01')
|
||||
rho_a = gas_a.density()
|
||||
|
||||
|
||||
# Use GRI-Mech 3.0 for stream b (methane) and for the mixer. If it is
|
||||
# desired to have a pure mixer, with no chemistry, use instead a
|
||||
# reaction mechanism for gas_b that has no reactions.
|
||||
gas_b = GRI30()
|
||||
gas_b.set(T = 300.0, P = OneAtm, X = 'CH4:1')
|
||||
rho_b = gas_b.density()
|
||||
|
||||
|
||||
# Create reservoirs for the two inlet streams and for the outlet
|
||||
# stream. The upsteam reservoirs could be replaced by reactors, which
|
||||
# might themselves be connected to reactors further upstream. The
|
||||
# outlet reservoir could be replaced with a reactor with no outlet, if
|
||||
# it is desired to integrate the composition leaving the mixer in
|
||||
# time, or by an arbitrary network of downstream reactors.
|
||||
res_a = Reservoir(gas_a)
|
||||
res_b = Reservoir(gas_b)
|
||||
downstream = Reservoir(gas_b)
|
||||
|
||||
|
||||
# Create a reactor for the mixer. A reactor is required instead of a
|
||||
# reservoir, since the state will change with time if the inlet mass
|
||||
# flow rates change or if there is chemistry occurring.
|
||||
mixer = Reactor(gas_b)
|
||||
|
||||
|
||||
# create two mass flow controllers connecting the upstream reservoirs
|
||||
# to the mixer, and set their mass flow rates to values corresponding
|
||||
# to stoichiometric combustion.
|
||||
mfc1 = MassFlowController(upstream = res_a,
|
||||
downstream = mixer,
|
||||
mdot = rho_a*2.5/0.21)
|
||||
|
||||
mfc2 = MassFlowController(upstream = res_b,
|
||||
downstream = mixer,
|
||||
mdot = rho_b*1.0)
|
||||
|
||||
|
||||
# add an igniter to ignite the mixture. The 'igniter' consists of a
|
||||
# stream of pure H.
|
||||
gas_c = IdealGasMix('h2o2.cti')
|
||||
gas_c.set(T = 300.0, P = OneAtm, X = 'H:1')
|
||||
igniter = Reactor(gas_c)
|
||||
|
||||
mfc3 = MassFlowController(upstream = igniter, downstream = mixer,
|
||||
mdot = 0.05)
|
||||
|
||||
|
||||
# connect the mixer to the downstream reservoir with a valve.
|
||||
outlet = Valve(upstream = mixer, downstream = downstream, Kv = 1.0)
|
||||
|
||||
sim = ReactorNet([mixer])
|
||||
|
||||
# Since the mixer is a reactor, we need to integrate in time to reach
|
||||
# steady state. A few residence times should be enough.
|
||||
t = 0.0
|
||||
for n in range(30):
|
||||
tres = mixer.mass()/(mfc1.massFlowRate() + mfc2.massFlowRate())
|
||||
t += 0.5*tres
|
||||
sim.advance(t)
|
||||
|
||||
# if ignited, turn the igniter off.
|
||||
# We also need to restart the integration in this case.
|
||||
if mixer.temperature() > 1200.0:
|
||||
mfc3.set(mdot = 0.0)
|
||||
sim.setInitialTime(t)
|
||||
|
||||
print '%14.5g %14.5g %14.5g %14.5g %14.5g' % (t, mixer.temperature(),
|
||||
mixer.enthalpy_mass(),
|
||||
mixer.pressure(),
|
||||
mixer.massFraction('CH4'))
|
||||
|
||||
# view the state of the gas in the mixer
|
||||
print mixer.contents()
|
||||
|
||||
|
|
@ -1,99 +0,0 @@
|
|||
"""
|
||||
Gas 1: a stoichiometric H2/O2/Ar mixture
|
||||
Gas 2: a wet CO/O2 mixture
|
||||
|
||||
-------------------------------------
|
||||
| || |
|
||||
| || |
|
||||
| gas 1 || gas 2 |
|
||||
| || |
|
||||
| || |
|
||||
-------------------------------------
|
||||
|
||||
The two volumes are connected by an adiabatic free piston.
|
||||
The piston speed is proportional to the pressure difference
|
||||
between the two chambers.
|
||||
|
||||
Note that each side uses a *different* reaction mechanism
|
||||
|
||||
"""
|
||||
from Cantera import *
|
||||
from Cantera.Reactor import *
|
||||
import sys
|
||||
|
||||
fmt = '%10.3f %10.1f %10.4f %10.4g %10.4g %10.4g %10.4g'
|
||||
print '%10s %10s %10s %10s %10s %10s %10s' % ('time [s]','T1 [K]','T2 [K]',
|
||||
'V1 [m^3]', 'V2 [m^3]',
|
||||
'V1+V2 [m^3]','X(CO)')
|
||||
|
||||
gas1 = importPhase('h2o2.cti')
|
||||
gas1.set(T = 900.0, P = OneAtm, X = 'H2:2, O2:1, AR:20')
|
||||
|
||||
gas2 = GRI30()
|
||||
gas2.set(T = 900.0, P = OneAtm, X = 'CO:2, H2O:0.01, O2:5')
|
||||
|
||||
r1 = Reactor(gas1, volume = 0.5)
|
||||
r2 = Reactor(gas2, volume = 0.1)
|
||||
w = Wall(left = r1, right = r2, K = 1.0e3)
|
||||
|
||||
reactors = ReactorNet([r1, r2])
|
||||
|
||||
tim = []
|
||||
t1 = []
|
||||
t2 = []
|
||||
v1 = []
|
||||
v2 = []
|
||||
v = []
|
||||
xco = []
|
||||
xh2 = []
|
||||
|
||||
for n in range(30):
|
||||
time = (n+1)*0.002
|
||||
reactors.advance(time)
|
||||
print fmt % (time, r1.temperature(), r2.temperature(),
|
||||
r1.volume(), r2.volume(), r1.volume() + r2.volume(),
|
||||
r2.moleFraction('CO'))
|
||||
|
||||
tim.append(time)
|
||||
t1.append(r1.temperature())
|
||||
t2.append(r2.temperature())
|
||||
v1.append(r1.volume())
|
||||
v2.append(r2.volume())
|
||||
v.append(r1.volume() + r2.volume())
|
||||
xco.append(r2.moleFraction('CO'))
|
||||
xh2.append(r1.moleFraction('H2'))
|
||||
|
||||
|
||||
# plot the results if matplotlib is installed.
|
||||
# see http://matplotlib.sourceforge.net to get it
|
||||
args = sys.argv
|
||||
if len(args) > 1 and (args[1] == '-plot' or
|
||||
args[1] == '-p' or
|
||||
args[1] == '--plot'):
|
||||
try:
|
||||
from matplotlib.pylab import *
|
||||
clf
|
||||
subplot(2,2,1)
|
||||
plot(tim,t1,'-',tim,t2,'r-')
|
||||
xlabel('Time (s)');
|
||||
ylabel('Temperature (K)');
|
||||
subplot(2,2,2)
|
||||
plot(tim,v1,'-',tim,v2,'r-',tim,v,'g-')
|
||||
xlabel('Time (s)');
|
||||
ylabel('Volume (m3)');
|
||||
subplot(2,2,3)
|
||||
plot(tim,xco);
|
||||
xlabel('Time (s)');
|
||||
ylabel('CO Mole Fraction (right)');
|
||||
subplot(2,2,4)
|
||||
plot(tim,xh2);
|
||||
xlabel('Time (s)');
|
||||
ylabel('H2 Mole Fraction (left)');
|
||||
show()
|
||||
except:
|
||||
print """matplotlib required.
|
||||
http://matplotlib.sourceforge.net"""
|
||||
|
||||
else:
|
||||
print """To view a plot of these results, run this script with the option -plot"""
|
||||
|
||||
|
|
@ -1,71 +0,0 @@
|
|||
"""
|
||||
|
||||
Constant-pressure, adiabatic kinetics simulation.
|
||||
|
||||
"""
|
||||
import sys
|
||||
|
||||
from Cantera import *
|
||||
from Cantera.Reactor import *
|
||||
from Cantera.Func import *
|
||||
from Cantera import rxnpath
|
||||
|
||||
gri3 = GRI30()
|
||||
|
||||
gri3.set(T = 1001.0, P = OneAtm, X = 'H2:2,O2:1,N2:4')
|
||||
r = Reactor(gri3)
|
||||
|
||||
env = Reservoir(Air())
|
||||
|
||||
# Define a wall between the reactor and the environment, and
|
||||
# make it flexible, so that the pressure in the reactor is held
|
||||
# at the environment pressure.
|
||||
w = Wall(r,env)
|
||||
w.set(K = 1.0e6) # set expansion parameter. dV/dt = KA(P_1 - P_2)
|
||||
w.set(A = 1.0)
|
||||
|
||||
sim = ReactorNet([r])
|
||||
time = 0.0
|
||||
tim = zeros(100,'d')
|
||||
data = zeros([100,5],'d')
|
||||
|
||||
for n in range(100):
|
||||
time += 1.e-5
|
||||
sim.advance(time)
|
||||
tim[n] = time
|
||||
data[n,0] = r.temperature()
|
||||
data[n,1] = r.moleFraction('OH')
|
||||
data[n,2] = r.moleFraction('H')
|
||||
data[n,3] = r.moleFraction('H2')
|
||||
print '%10.3e %10.3f %10.3f %14.6e' % (sim.time(), r.temperature(),
|
||||
r.pressure(), r.intEnergy_mass())
|
||||
|
||||
|
||||
# plot the results if matplotlib is installed.
|
||||
# see http://matplotlib.sourceforge.net to get it
|
||||
args = sys.argv
|
||||
if len(args) > 1 and args[1] == '-plot':
|
||||
try:
|
||||
from matplotlib.pylab import *
|
||||
clf
|
||||
subplot(2,2,1)
|
||||
plot(tim,data[:,0])
|
||||
xlabel('Time (s)');
|
||||
ylabel('Temperature (K)');
|
||||
subplot(2,2,2)
|
||||
plot(tim,data[:,1])
|
||||
xlabel('Time (s)');
|
||||
ylabel('OH Mole Fraction');
|
||||
subplot(2,2,3)
|
||||
plot(tim,data[:,2]);
|
||||
xlabel('Time (s)');
|
||||
ylabel('H Mole Fraction');
|
||||
subplot(2,2,4)
|
||||
plot(tim,data[:,3]);
|
||||
xlabel('Time (s)');
|
||||
ylabel('H2 Mole Fraction');
|
||||
show()
|
||||
except:
|
||||
pass
|
||||
else:
|
||||
print """To view a plot of these results, run this script with the option -plot"""
|
||||
|
|
@ -1,127 +0,0 @@
|
|||
"""
|
||||
|
||||
This script simulates the following situation. A closed cylinder with
|
||||
volume 2 m^3 is divided into two equal parts by a massless piston that
|
||||
moves with speed proportional to the pressure difference between the
|
||||
two sides. It is initially held in place in the middle. One side is
|
||||
filled with 1000 K argon at 20 atm, and the other with a combustible
|
||||
500 K methane/air mixture at 0.1 atm (phi = 1.1). At t = 0 the piston
|
||||
is released and begins to move due to the large pressure difference,
|
||||
compressing and heating the methane/air mixture, which eventually
|
||||
explodes. At the same time, the argon cools as it expands. The piston
|
||||
is adiabatic, but some heat is lost through the outer cylinder walls
|
||||
to the environment.
|
||||
|
||||
Note that this simulation, being zero-dimensional, takes no account of
|
||||
shock wave propagation. It is somewhat artifical, but nevertheless
|
||||
instructive.
|
||||
|
||||
"""
|
||||
|
||||
import sys
|
||||
from Cantera import *
|
||||
from Cantera.Reactor import *
|
||||
from Cantera.Func import *
|
||||
|
||||
#-----------------------------------------------------------------------
|
||||
# First create each gas needed, and a reactor or reservoir for each one.
|
||||
#-----------------------------------------------------------------------
|
||||
|
||||
# create an argon gas object and set its state. This function is
|
||||
# defined in module Cantera.gases, as are functions 'Air()', and
|
||||
# 'GRI30()'
|
||||
|
||||
ar = Argon()
|
||||
ar.set(T = 1000.0, P = 20.0*OneAtm, X = 'AR:1')
|
||||
|
||||
# create a reactor to represent the side of the cylinder filled with argon
|
||||
r1 = Reactor(ar)
|
||||
|
||||
|
||||
# create a reservoir for the environment, and fill it with air.
|
||||
env = Reservoir(Air())
|
||||
|
||||
|
||||
# use GRI-Mech 3.0 for the methane/air mixture, and set its initial state
|
||||
gri3 = GRI30()
|
||||
|
||||
gri3.set(T = 500.0, P = 0.2*OneAtm, X = 'CH4:1.1, O2:2, N2:7.52')
|
||||
|
||||
# create a reactor for the methane/air side
|
||||
r2 = Reactor(gri3)
|
||||
|
||||
|
||||
#---------------------------------------------------------------------
|
||||
# Now couple the reactors by defining common walls that may move (a piston)
|
||||
# or conduct heat
|
||||
#----------------------------------------------------------------------
|
||||
|
||||
# add a flexible wall (a piston) between r2 and r1
|
||||
w = Wall(r2, r1)
|
||||
w.set(area = 1.0, K=0.5e-4, U = 100.0)
|
||||
|
||||
|
||||
# heat loss to the environment. Heat loss always occur through walls,
|
||||
# so we create a wall separating r1 from the environment, give it a
|
||||
# non-zero area, and specify the overall heat transfer coefficient
|
||||
# through the wall.
|
||||
w2 = Wall(r2, env)
|
||||
w2.set(area = 1.0, U=500.0)
|
||||
|
||||
sim = ReactorNet([r1, r2])
|
||||
|
||||
# Now the problem is set up, and we're ready to solve it.
|
||||
print 'finished setup, begin solution...'
|
||||
|
||||
time = 0.0
|
||||
f = open('piston.csv','w')
|
||||
writeCSV(f,['time (s)','T1 (K)','P1 (Bar)','V1 (m3)',
|
||||
'T2 (K)','P2 (Bar)','V2 (m3)'])
|
||||
temp = zeros([300, 2], 'd')
|
||||
pres = zeros([300, 2], 'd')
|
||||
vol = zeros([300, 2], 'd')
|
||||
tm = zeros(300,'d')
|
||||
for n in range(300):
|
||||
time += 4.e-4
|
||||
print time, r2.temperature(),n
|
||||
sim.advance(time)
|
||||
tm[n] = time
|
||||
temp[n,:] = [r1.temperature(), r2.temperature()]
|
||||
pres[n,:] = [1.0e-5*r1.pressure(), 1.0e-5*r2.pressure()]
|
||||
vol[n,:] = [r1.volume(), r2.volume()]
|
||||
writeCSV(f, [tm[n], temp[n,0], pres[n,0], vol[n,0],
|
||||
temp[n,1], pres[n,1], vol[n,1]])
|
||||
f.close()
|
||||
import os
|
||||
print 'Output written to file piston.csv'
|
||||
print 'Directory: '+os.getcwd()
|
||||
|
||||
args = sys.argv
|
||||
if len(args) > 1 and args[1] == '-plot':
|
||||
try:
|
||||
from matplotlib.pylab import *
|
||||
clf
|
||||
subplot(2,2,1)
|
||||
plot(tm, temp[:,0],'g-',tm, temp[:,1],'b-')
|
||||
legend(['Reactor 1','Reactor 2'],2)
|
||||
xlabel('Time (s)');
|
||||
ylabel('Temperature (K)');
|
||||
|
||||
subplot(2,2,2)
|
||||
plot(tm, pres[:,0],'g-',tm, pres[:,1],'b-')
|
||||
legend(['Reactor 1','Reactor 2'],2)
|
||||
xlabel('Time (s)');
|
||||
ylabel('Pressure (Bar)');
|
||||
|
||||
subplot(2,2,3)
|
||||
plot(tm, vol[:,0],'g-',tm, vol[:,1],'b-')
|
||||
legend(['Reactor 1','Reactor 2'],2)
|
||||
xlabel('Time (s)');
|
||||
ylabel('Volume (m^3)');
|
||||
|
||||
show()
|
||||
except:
|
||||
pass
|
||||
else:
|
||||
print """To view a plot of these results, run this script with the option -plot"""
|
||||
|
||||
12
Cantera/python/examples/reactors/reactor2_sim/.cvsignore
Normal file
12
Cantera/python/examples/reactors/reactor2_sim/.cvsignore
Normal file
|
|
@ -0,0 +1,12 @@
|
|||
Makefile
|
||||
air.xml
|
||||
argon.xml
|
||||
ct2ctml.log
|
||||
diff_csv.txt
|
||||
diff_out_0.txt
|
||||
gri30.xml
|
||||
output_0.txt
|
||||
piston.csv
|
||||
runtest
|
||||
transport_log.xml
|
||||
|
||||
|
|
@ -1,97 +0,0 @@
|
|||
"""
|
||||
|
||||
Constant-pressure, adiabatic kinetics simulation with sensitivity analysis
|
||||
|
||||
"""
|
||||
import sys
|
||||
|
||||
from Cantera import *
|
||||
from Cantera.Reactor import *
|
||||
from Cantera.Func import *
|
||||
|
||||
gri3 = GRI30()
|
||||
temp = 1500.0
|
||||
pres = OneAtm
|
||||
|
||||
gri3.set(T = temp, P = pres, X = 'CH4:0.1, O2:2, N2:7.52')
|
||||
r = Reactor(gri3)
|
||||
|
||||
air = Air()
|
||||
air.set(T = temp, P = pres)
|
||||
env = Reservoir(air)
|
||||
|
||||
# Define a wall between the reactor and the environment, and
|
||||
# make it flexible, so that the pressure in the reactor is held
|
||||
# at the environment pressure.
|
||||
w = Wall(r,env)
|
||||
w.set(K = 1.0e6) # set expansion parameter. dV/dt = KA(P_1 - P_2)
|
||||
w.set(A = 1.0)
|
||||
|
||||
# enable sensitivity with respect to the rates of the first 10
|
||||
# reactions (reactions 0 through 9)
|
||||
r.addSensitivityReaction(reactions = range(10))
|
||||
|
||||
sim = ReactorNet([r])
|
||||
|
||||
# set the tolerances for the solution and for the sensitivity
|
||||
# coefficients
|
||||
sim.setTolerances(rtol = 1.0e-6, atol = 1.0e-15,
|
||||
rtolsens = 1.0e-5, atolsens = 1.0e-5)
|
||||
time = 0.0
|
||||
np = 400
|
||||
tim = zeros(np,'d')
|
||||
data = zeros([np,6],'d')
|
||||
|
||||
for n in range(np):
|
||||
time += 5.0e-6
|
||||
sim.advance(time)
|
||||
tim[n] = time
|
||||
data[n,0] = r.temperature()
|
||||
data[n,1] = r.moleFraction('OH')
|
||||
data[n,2] = r.moleFraction('H')
|
||||
data[n,3] = r.moleFraction('CH4')
|
||||
|
||||
# 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)
|
||||
|
||||
print '%10.3e %10.3f %10.3f %14.6e' % (sim.time(), r.temperature(),
|
||||
r.pressure(), r.intEnergy_mass())
|
||||
#sim.sensitivity("OH",0))
|
||||
|
||||
|
||||
# plot the results if matplotlib is installed.
|
||||
# see http://matplotlib.sourceforge.net to get it
|
||||
args = sys.argv
|
||||
if len(args) > 1 and args[1] == '-plot':
|
||||
try:
|
||||
from matplotlib.pylab import *
|
||||
clf
|
||||
subplot(2,2,1)
|
||||
plot(tim,data[:,0])
|
||||
xlabel('Time (s)');
|
||||
ylabel('Temperature (K)');
|
||||
subplot(2,2,2)
|
||||
plot(tim,data[:,1])
|
||||
xlabel('Time (s)');
|
||||
ylabel('OH Mole Fraction');
|
||||
subplot(2,2,3)
|
||||
plot(tim,data[:,2]);
|
||||
xlabel('Time (s)');
|
||||
ylabel('H Mole Fraction');
|
||||
subplot(2,2,4)
|
||||
plot(tim,data[:,3]);
|
||||
xlabel('Time (s)');
|
||||
ylabel('H2 Mole Fraction');
|
||||
figure(2)
|
||||
plot(tim,data[:,4],'-',tim,data[:,5],'-g')
|
||||
legend([r.sensParamName(2),r.sensParamName(3)],'best')
|
||||
xlabel('Time (s)');
|
||||
ylabel('OH Sensitivity');
|
||||
show()
|
||||
except:
|
||||
print 'could not make plots'
|
||||
else:
|
||||
print """To view a plot of these results, run this script with the option -plot"""
|
||||
|
|
@ -1,179 +0,0 @@
|
|||
# This example solves a plug flow reactor problem, where the chemistry
|
||||
# is surface chemistry. The specific problem simulated is the partial
|
||||
# oxidation of methane over a platinum catalyst in a packed bed
|
||||
# reactor.
|
||||
|
||||
from Cantera import *
|
||||
from Cantera.Reactor import *
|
||||
from Cantera import rxnpath
|
||||
import math
|
||||
import sys
|
||||
|
||||
|
||||
#######################################################################
|
||||
|
||||
# unit conversion factors to SI
|
||||
|
||||
cm = 0.01
|
||||
minute = 60.0
|
||||
|
||||
|
||||
|
||||
#######################################################################
|
||||
#
|
||||
# Input Parameters
|
||||
#
|
||||
#######################################################################
|
||||
|
||||
tc = 800.0 # Temperature in Celsius
|
||||
|
||||
length = 0.3 * cm # Catalyst bed length
|
||||
area = 1.0 * cm * cm # Catalyst bed area
|
||||
cat_area_per_vol = 1000.0 / cm # Catalyst particle surface area
|
||||
# per unit volume
|
||||
velocity = 40.0 * cm / minute # gas velocity
|
||||
porosity = 0.3 # Catalyst bed porosity
|
||||
|
||||
# input file containing the surface reaction mechanism
|
||||
cti_file = 'methane_pox_on_pt.cti'
|
||||
|
||||
# The PFR will be simulated by a chain of 'NReactors' stirred
|
||||
# reactors.
|
||||
NReactors = 200
|
||||
dt = 1.0
|
||||
|
||||
|
||||
#####################################################################
|
||||
|
||||
|
||||
t = tc + 273.15 # convert to Kelvin
|
||||
|
||||
# import the gas model
|
||||
gas = importPhase(cti_file,'gas')
|
||||
|
||||
# set the initial conditions
|
||||
gas.set(T = t, P = OneAtm, X = 'CH4:1, O2:1.5, AR:0.1')
|
||||
rho0 = gas.density()
|
||||
nsp = gas.nSpecies()
|
||||
g_names = gas.speciesNames()
|
||||
|
||||
# import the surface model
|
||||
surf = importInterface(cti_file,'Pt_surf', [gas])
|
||||
surf.setTemperature(t)
|
||||
s_names = surf.speciesNames()
|
||||
nsurf = surf.nSpecies()
|
||||
|
||||
rlen = length/NReactors
|
||||
rvol = area * rlen * porosity
|
||||
|
||||
names = gas.speciesNames()
|
||||
|
||||
f = open('surf_pfr_output.csv','w')
|
||||
writeCSV(f, ['Distance (mm)', 'T (C)', 'P (atm)'] + g_names + s_names)
|
||||
|
||||
# catalyst area in one reactor
|
||||
cat_area = cat_area_per_vol*rvol
|
||||
|
||||
mass_flow_rate = velocity * rho0 * area
|
||||
|
||||
# The plug flow reactor is represented by a linear chain of
|
||||
# zero-dimensional reactors. The gas at the inlet to the first one has
|
||||
# the specified inlet composition, and for all others the inlet
|
||||
# composition is fixed at the composition of the reactor immediately
|
||||
# upstream. Since in a PFR model there is no diffusion, the upstream
|
||||
# reactors are not affected by any downstream reactors, and therefore
|
||||
# the problem may be solved by simply marching from the first to last
|
||||
# reactor, integrating each one to steady state.
|
||||
|
||||
for n in range(NReactors):
|
||||
|
||||
# create a new reactor
|
||||
r = Reactor(contents = gas, energy = 'off', volume = rvol)
|
||||
|
||||
# create a reservoir to represent the reactor immediately
|
||||
# upstream. Note that the gas object is set already to the
|
||||
# state of the upstream reactor
|
||||
upstream = Reservoir(gas, name = 'upstream')
|
||||
|
||||
# create a reservoir for the reactor to exhaust into. The
|
||||
# composition of this reservoir is irrelevant.
|
||||
downstream = Reservoir(gas, name = 'downstream')
|
||||
|
||||
# use a 'Wall' object to implement the reacting surface in the
|
||||
# reactor. Since walls have to be installed between two
|
||||
# reactors/reserviors, we'll install it between the upstream
|
||||
# reservoir and the reactor. The area is set to the desired
|
||||
# catalyst area in the reactor, and surface reactions are
|
||||
# included only on the side facing the reactor.
|
||||
w = Wall(left = upstream, right = r, A = cat_area, kinetics = [None, surf])
|
||||
# We need a valve between the reactor and the downstream reservoir.
|
||||
# This will determine the pressure in the reactor. Set Kv large
|
||||
# enough that the pressure difference is very small.
|
||||
v = Valve(upstream = r, downstream = downstream, Kv = 3.0e-6)
|
||||
|
||||
# The mass flow rate into the reactor will be fixed by using a
|
||||
# MassFlowController object.
|
||||
m = MassFlowController(upstream = upstream,
|
||||
downstream = r, mdot = mass_flow_rate)
|
||||
|
||||
sim = ReactorNet([upstream, r, downstream])
|
||||
|
||||
# set relative and absolute tolerances on the simulation
|
||||
sim.setTolerances(rtol = 1.0e-6, atol = 1.0e-15)
|
||||
|
||||
time = 0
|
||||
while 1 > 0:
|
||||
time = time + dt
|
||||
sim.advance(time)
|
||||
|
||||
# check whether surface coverages are in steady
|
||||
# state. This will be the case if the creation and
|
||||
# destruction rates for a surface (but not gas) species
|
||||
# are equal.
|
||||
alldone = 1
|
||||
|
||||
# Note: netProduction = creation - destruction. By
|
||||
# supplying the surface object as an argument, only the
|
||||
# values for the surface species are returned by these
|
||||
# methods
|
||||
sdot = surf.netProductionRates(surf)
|
||||
cdot = surf.creationRates(surf)
|
||||
ddot = surf.destructionRates(surf)
|
||||
for ks in range(nsurf):
|
||||
ratio = sdot[ks]/(cdot[ks] + ddot[ks])
|
||||
if ratio < 0.0: ratio = -ratio
|
||||
if ratio > 1.0e-11 or time < 10*dt:
|
||||
alldone = 0
|
||||
|
||||
if alldone: break
|
||||
|
||||
# set the gas object state to that of this reactor, in
|
||||
# preparation for the simulation of the next reactor
|
||||
# downstream, where this object will set the inlet conditions
|
||||
gas = r.contents()
|
||||
|
||||
dist = n*rlen * 1.0e3 # distance in mm
|
||||
|
||||
# write the gas mole fractions and surface coverages
|
||||
# vs. distance
|
||||
writeCSV(f, [dist, r.temperature() - 273.15,
|
||||
r.pressure()/OneAtm] + list(gas.moleFractions())
|
||||
+ list(surf.coverages()))
|
||||
|
||||
f.close()
|
||||
|
||||
# make a reaction path diagram tracing carbon. This diagram will show
|
||||
# the pathways by the carbon entering the bed in methane is convered
|
||||
# into CO and CO2. The diagram will be specifically for the exit of
|
||||
# the bed; if the pathways are desired at some interior point, then
|
||||
# put this statement inside the above loop.
|
||||
#
|
||||
# To process this diagram, give the command on the command line
|
||||
# after running this script:
|
||||
# dot -Tps < carbon_pathways.dot > carbon_pathways.ps
|
||||
# This will generate the diagram in Postscript.
|
||||
|
||||
element = 'C'
|
||||
rxnpath.write(surf, element, 'carbon_pathways.dot')
|
||||
|
||||
|
||||
Loading…
Add table
Reference in a new issue