diff --git a/Cantera/python/examples/reactors/surf_pfr.py b/Cantera/python/examples/reactors/surf_pfr.py new file mode 100644 index 000000000..6a0529c29 --- /dev/null +++ b/Cantera/python/examples/reactors/surf_pfr.py @@ -0,0 +1,168 @@ +# 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 Centigrade + +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 +# reactors, 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 obbject. + 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. + alldone = 1 + sdot = surf.netProductionRates(surf) + cdot = surf.creationRates(surf) + ddot = surf.destructionRates(surf) + for ks in range(nsurf): + ratio = sdot[ks]/(cdot[ks] + ddot[ks]) + #print ks, ratio + if ratio < 0.0: ratio = -ratio + if ratio > 1.0e-11 or time < 10*dt: + alldone = 0 + if alldone: break + + 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') + +