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

[Cython] Translated more samples to use the new API

Browse source
This commit is contained in:
Ray Speth 2013-02-01 23:40:41 +00:00
parent eb01aee6f5
commit 2be392e6e1
9 changed files with 1072 additions and 1 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

36
interfaces/cython/cantera/examples/dusty_gas.py Normal file
View file

@ -0,0 +1,36 @@
"""
Dusty Gas transport model.
The Dusty Gas model is a mulicomponent transport model for gas transport
through the pores of a stationary porous medium. This example shows how to
create a transport manager that implements the Dusty Gas model and use it to
compute the multicomponent diffusion coefficients.
"""
import cantera as ct
# create a gas-phase object to represent the gas in the pores, with a
# dusty gas transport manager
g = ct.DustyGas('h2o2.cti')
# set the gas state
g.TPX = 500.0, ct.one_atm, "OH:1, H:2, O2:3, O:1.0E-8, H2:1.0E-8, H2O:1.0E-8, H2O2:1.0E-8, HO2:1.0E-8, AR:1.0E-8"
# set its parameters
g.porosity = 0.2
g.tortuosity = 4.0
g.mean_pore_radius = 1.5e-7
g.mean_particle_diameter = 1.5e-6 # lengths in meters
# print the multicomponent diffusion coefficients
print(g.multi_diff_coeffs)
# compute molar species fluxes
T1, rho1, Y1 = g.TDY
g.TP = g.T, 1.2 * ct.one_atm
T2, rho2, Y2 = g.TDY
delta = 0.001
print(g.molar_fluxes(T1, T1, rho1, rho1, Y1, Y1, delta))
print(g.molar_fluxes(T1, T2, rho1, rho2, Y1, Y2, delta))

41
interfaces/cython/cantera/examples/multiphase/plasma_equilibrium.py Normal file
View file

@ -0,0 +1,41 @@
"""
An equilibrium example with charged species in the gas phase
and multiple condensed phases.
"""
import cantera as ct
import csv
# create objects representing the gas phase and the condensed phases. The gas
# is a mixture of multiple species, and the condensed phases are all modeled
# as incompressible stoichiometric substances. See file KOH.cti for more
# information.
phases = ct.import_phases('KOH.cti', ['K_solid', 'K_liquid', 'KOH_a', 'KOH_b',
'KOH_liquid', 'K2O2_solid', 'K2O_solid',
'KO2_solid', 'ice', 'liquid_water',
'KOH_plasma'])
# create the Mixture object from the list of phases
mix = ct.Mixture(phases)
csvfile = open('equil_koh.csv', 'w')
writer = csv.writer(csvfile)
writer.writerow(['T'] + mix.species_names)
# loop over temperature
for n in range(100):
t = 350.0 + 50.0*n
print('T = {}'.format(t))
mix.T = t
mix.P = ct.one_atm
mix.species_moles = "K:1.03, H2:2.12, O2:0.9"
# set the mixture to a state of chemical equilibrium holding
# temperature and pressure fixed
# mix.equilibrate("TP",maxsteps=10000,loglevel=1)
mix.equilibrate("TP", max_steps=10000, log_level=0)
# write out the moles of each species
writer.writerow([t] + list(mix.species_moles))
csvfile.close()

116
interfaces/cython/cantera/examples/onedim/flame_fixed_T.py Normal file
View file

@ -0,0 +1,116 @@
"""
FIXED_T_FLAME - A burner-stabilized, premixed methane/air flat flame with
multicomponent transport properties and a specified temperature profile.
"""
import cantera as ct
# read temperature vs. position data from a file.
# The file is assumed to have one z, T pair per line, separated by a comma.
def getTempData(filename):
# open the file containing the temperature data for reading
lines = open(filename).readlines()
# check for unix/Windows/Mac line ending problems
if len(lines) == 1:
print('Warning: only one line found.')
print('Possible text file line-ending problem?')
print('The one line found is: ', lines[0])
z = []
T = []
for line in lines:
if line[0] == '#': # use '#' as the comment character
continue
try:
zval, tval = line.split(',')
z.append(float(zval))
T.append(float(tval))
except Exception:
pass
print('read {} temperature values.'.format(len(z)))
# convert z values into non-dimensional relative positions.
n = len(z)
zmax = z[n-1]
for i in range(n):
z[i] = z[i]/zmax
return z,T
################################################################
# parameter values
p = ct.one_atm # pressure
tburner = 373.7 # burner temperature
mdot = 0.04 # kg/m^2/s
comp = 'CH4:0.65, O2:1, N2:3.76' # premixed gas composition
# The solution domain is chosen to be 1 cm, and a point very near the
# downstream boundary is added to help with the zero-gradient boundary
# condition at this boundary.
initial_grid = [0.0, 0.0025, 0.005, 0.0075, 0.0099, 0.01] # m
tol_ss = [1.0e-5, 1.0e-9] # [rtol atol] for steady-state problem
tol_ts = [1.0e-5, 1.0e-4] # [rtol atol] for time stepping
loglevel = 1 # amount of diagnostic output (0 to 5)
refine_grid = True # 'True' to enable refinement
################ create the gas object ########################
#
# This object will be used to evaluate all thermodynamic, kinetic, and
# transport properties. It is created with two transport managers, to enable
# switching from mixture-averaged to multicomponent transport on the last
# solution.
gas = ct.Solution('gri30.xml', 'gri30_mix')
# set its state to that of the unburned gas at the burner
gas.TPX = tburner, p, comp
# create the BurnerFlame object.
f = ct.BurnerFlame(gas=gas, grid=initial_grid)
# set the properties at the burner
f.burner.mdot = mdot
f.burner.X = comp
f.burner.T = tburner
# read in the fixed temperature profile
[zloc, tvalues] = getTempData('tdata.dat')
# set the temperature profile to the values read in
f.flame.set_fixed_temp_profile(zloc, tvalues)
f.flame.set_steady_tolerances(default=tol_ss)
f.flame.set_transient_tolerances(default=tol_ts)
# show the initial estimate for the solution
f.show_solution()
# don't solve the energy equation
f.energy_enabled = False
# first solve the flame with mixture-averaged transport properties
f.set_refine_criteria(ratio=3.0, slope=0.3, curve=1)
f.set_max_jac_age(50, 50)
f.set_time_step(1.0e-5, [1, 2, 5, 10, 20])
f.solve(loglevel, refine_grid)
f.save('ch4_flame_fixed_T.xml','mixav',
'solution with mixture-averaged transport')
print('\n\n switching to multicomponent transport...\n\n')
f.transport_model = 'Multi'
f.set_refine_criteria(ratio=3.0, slope=0.1, curve=0.2)
f.solve(loglevel, refine_grid)
f.save('ch4_flame_fixed_T.xml','multi',
'solution with multicomponent transport')
# write the velocity, temperature, density, and mole fractions to a CSV file
f.write_csv('flame_fixed_T.csv', quiet=False)
f.show_stats()

74
interfaces/cython/cantera/examples/onedim/tdata.dat Normal file
View file

@ -0,0 +1,74 @@
#
# This data file lists temperature vs. height values for a burner-stabilized flame.
# This file is used by example 'fixed_T_flame.py'.
#
0, 373.7
0.00015625, 465.4070428
0.000234375, 510.4311676
0.000390625, 599.5552837
0.00046875, 643.8342938
0.000507813, 665.9335545
0.000546875, 688.0122338
0.000625, 732.1284327
0.000664062, 754.1744755
0.000703125, 776.2170662
0.000742188, 798.2588757
0.00078125, 820.3020011
0.000820313, 842.348001
0.000859375, 864.3979228
0.000898437, 886.4523159
0.0009375, 908.5112198
0.001015625, 952.6396629
0.001054688, 974.7018199
0.00109375, 996.7515831
0.001132813, 1018.777651
0.001171875, 1040.765863
0.001210938, 1062.69948
0.00125, 1084.558639
0.001289062, 1106.320078
0.001328125, 1127.956918
0.001367187, 1149.438472
0.00140625, 1170.730129
0.001445313, 1191.793309
0.001484375, 1212.585506
0.001523438, 1233.060477
0.0015625, 1253.168589
0.001601563, 1272.857384
0.001640625, 1292.072391
0.00171875, 1328.859767
0.001757812, 1346.323998
0.001796875, 1363.101361
0.001835937, 1379.147594
0.001875, 1394.425274
0.001914063, 1408.905834
0.001953125, 1422.569115
0.001992188, 1435.40408
0.00203125, 1447.410648
0.002070313, 1458.597668
0.002109375, 1468.982722
0.002148438, 1478.590978
0.0021875, 1487.453914
0.002226563, 1495.607879
0.002265625, 1503.092709
0.002304688, 1509.950449
0.00234375, 1516.224147
0.002382813, 1521.956853
0.002421875, 1527.19079
0.002460938, 1531.966722
0.0025, 1536.32348
0.002578125, 1543.891739
0.00265625, 1550.203579
0.002734375, 1555.480771
0.0028125, 1559.908135
0.002890625, 1563.637879
0.00296875, 1566.794144
0.003046875, 1569.477867
0.003125, 1571.77099
0.00328125, 1575.385829
0.0034375, 1578.108169
0.00359375, 1580.194856
0.00375, 1581.820666
0.00390625, 1583.106578
0.0087, 1589.51315
0.01, 1589.578955

132
interfaces/cython/cantera/examples/surface_chemistry/catalytic_combustion.py Normal file
View file

@ -0,0 +1,132 @@
"""
CATCOMB -- Catalytic combustion of methane on platinum.
This script solves a catalytic combustion problem. A stagnation flow is set
up, with a gas inlet 10 cm from a platinum surface at 900 K. The lean,
premixed methane/air mixture enters at ~ 6 cm/s (0.06 kg/m2/s), and burns
catalytically on the platinum surface. Gas-phase chemistry is included too,
and has some effect very near the surface.
The catalytic combustion mechanism is from Deutschman et al., 26th
Symp. (Intl.) on Combustion,1996 pp. 1747-1754
"""
import numpy as np
import cantera as ct
# Parameter values are collected here to make it easier to modify them
p = ct.one_atm # pressure
tinlet = 300.0 # inlet temperature
tsurf = 900.0 # surface temperature
mdot = 0.06 # kg/m^2/s
transport = 'Mix' # transport model
# We will solve first for a hydrogen/air case to use as the initial estimate
# for the methane/air case
# composition of the inlet premixed gas for the hydrogen/air case
comp1 = 'H2:0.05, O2:0.21, N2:0.78, AR:0.01'
# composition of the inlet premixed gas for the methane/air case
comp2 = 'CH4:0.095, O2:0.21, N2:0.78, AR:0.01'
# the initial grid, in meters. The inlet/surface separation is 10 cm.
initial_grid = [0.0, 0.02, 0.04, 0.06, 0.08, 0.1] # m
# numerical parameters
tol_ss = [1.0e-5, 1.0e-9] # [rtol, atol] for steady-state problem
tol_ts = [1.0e-4, 1.0e-9] # [rtol, atol] for time stepping
loglevel = 1 # amount of diagnostic output (0 to 5)
refine_grid = True # enable or disable refinement
################ create the gas object ########################
#
# This object will be used to evaluate all thermodynamic, kinetic, and
# transport properties. The gas phase will be taken from the definition of
# phase 'gas' in input file 'ptcombust.cti,' which is a stripped-down version
# of GRI-Mech 3.0.
gas = ct.Solution('ptcombust.cti', 'gas')
gas.TPX = tinlet, p, comp1
################ create the interface object ##################
#
# This object will be used to evaluate all surface chemical production rates.
# It will be created from the interface definition 'Pt_surf' in input file
# 'ptcombust.cti,' which implements the reaction mechanism of Deutschmann et
# al., 1995 for catalytic combustion on platinum.
#
surf_phase = ct.Interface('ptcombust.cti', 'Pt_surf', [gas])
surf_phase.TP = tsurf, p
# integrate the coverage equations in time for 1 s, holding the gas
# composition fixed to generate a good starting estimate for the coverages.
surf_phase.advance_coverages(1.0)
# create the object that simulates the stagnation flow, and specify an initial
# grid
sim = ct.ImpingingJet(gas=gas, grid=initial_grid, surface=surf_phase)
# Objects of class StagnationFlow have members that represent the gas inlet
# ('inlet') and the surface ('surface'). Set some parameters of these objects.
sim.inlet.mdot = mdot
sim.inlet.T = tinlet
sim.inlet.X = comp1
sim.surface.T = tsurf
# Set error tolerances
sim.flame.set_steady_tolerances(default=tol_ss)
sim.flame.set_transient_tolerances(default=tol_ts)
# Show the initial solution estimate
sim.show_solution()
# Solving problems with stiff chemistry coulpled to flow can require a
# sequential approach where solutions are first obtained for simpler problems
# and used as the initial guess for more difficult problems.
# start with the energy equation on (default is 'off')
sim.energy_enabled = True
# disable the surface coverage equations, and turn off all gas and surface
# chemistry.
sim.surface.coverage_enabled = False
surf_phase.set_multiplier(0.0)
gas.set_multiplier(0.0)
# solve the problem, refining the grid if needed, to determine the non-
# reacting velocity and temperature distributions
sim.solve(loglevel, refine_grid)
# now turn on the surface coverage equations, and turn the chemistry on slowly
sim.surface.coverage_enabled = True
for mult in np.logspace(-5, 0, 6):
surf_phase.set_multiplier(mult)
gas.set_multiplier(mult)
print('Multiplier =', mult)
sim.solve(loglevel, refine_grid)
# At this point, we should have the solution for the hydrogen/air problem.
sim.show_solution()
# Now switch the inlet to the methane/air composition.
sim.inlet.X = comp2
# set more stringent grid refinement criteria
sim.set_refine_criteria(100.0, 0.15, 0.2, 0.0)
# solve the problem for the final time
sim.solve(loglevel, refine_grid)
# show the solution
sim.show_solution()
# save the solution in XML format. The 'restore' method can be used to restart
# a simulation from a solution stored in this form.
sim.save("catcomb.xml", "soln1")
# save selected solution components in a CSV file for plotting in
# Excel or MATLAB.
sim.write_csv('catalytic_combustion.csv', quiet=False)
sim.show_stats(0)

54
interfaces/cython/cantera/examples/surface_chemistry/diamond_cvd.py Normal file
View file

@ -0,0 +1,54 @@
"""
A CVD example.
This example computes the growth rate of a diamond film according to a
simplified version of a particular published growth mechanism (see file
diamond.cti for details). Only the surface coverage equations are solved here;
the gas composition is fixed. (For an example of coupled gas- phase and
surface, see catalytic_combustion.py.) Atomic hydrogen plays an important
role in diamond CVD, and this example computes the growth rate and surface
coverages as a function of [H] at the surface for fixed temperature and [CH3].
"""
import csv
import cantera as ct
print('\n****** CVD Diamond Example ******\n')
# import the models for the gas and bulk diamond
g, dbulk = ct.import_phases('diamond.cti', ['gas', 'diamond'])
# import the model for the diamond (100) surface
d = ct.Interface('diamond.cti', 'diamond_100', [g, dbulk])
ns = d.n_species
mw = dbulk.molecular_weights[0]
t = 1200.0
x = g.X
p = 20.0 * ct.one_atm / 760.0 # 20 Torr
g.TP = t, p
ih = g.species_index('H')
xh0 = x[ih]
f = open('diamond.csv', 'w')
writer = csv.writer(f)
writer.writerow(['H mole Fraction', 'Growth Rate (microns/hour)'] +
d.species_names)
iC = d.kinetics_species_index(dbulk.species_index('C(d)'), 1)
for n in range(20):
x[ih] /= 1.4
g.TPX = t, p, x
d.advance_coverages(10.0) # integrate the coverages to steady state
carbon_dot = d.net_production_rates[iC]
mdot = mw * carbon_dot
rate = mdot / dbulk.density
writer.writerow([x[ih], rate * 1.0e6 * 3600.0] + list(d.coverages))
f.close()
print('H concentration, growth rate, and surface coverages '
'written to file diamond.csv')

357
interfaces/cython/cantera/examples/surface_chemistry/sofc.cti Normal file
View file

@ -0,0 +1,357 @@
#########################################################################
#
# This is a an example input file that defines models for phases and
# interfaces that could be used, for example, to simulate a solid
# oxide fuel cell. Note, however, that reaction rate coefficients and
# species thermochemistry ARE NOT REAL VALUES - they are chosen only
# for the purposes of this example.
#
#########################################################################
# since Cantera input files are actually executable Python scripts,
# we can put any valid Python statements in the input file. Here we
# import the value of R from Cantera.
from Cantera import GasConstant
# These units will be used by default for any quantities entered
# without units. Quantities with compound units (e.g. concentration)
# will be constructed from these - the units of concentration will be
# mol/cm^3, etc.
units(length = "cm", time = "s", quantity = "mol", act_energy = "kJ/mol")
# Turn on mechanism validation to detect unbalanced reactions, if any
validate()
#------------------------------------------------------------------
#
# parameters
#
#------------------------------------------------------------------
# a few numeric parameters are collected here to allow easy modification.
# this temperature is used to initialize objects. But since
# scripts/programs usually set the temperature, it is not really
# necessary.
tc = 800.0 # temperature in C
tt = tc + 273.15 # temperature in K
# these values are defined here only so they may be easily changed to
# assess the effects of the oxide thermochemistry. For work at a
# single temperature, all that we really need is g = h -
# Ts. Therefore, it is somewhat arbitrary to assign separately
# enthalpies and entropies (but this is what the input format
# requires).
hox = (-170.0, 'kJ/mol') # enthalpy of an oxygen ion
sox = (50.0, 'J/K/mol') # entropy of an oxygen ion
hhydrox = (-220.0, 'kJ/mol') # enthalpy of a surface hydroxyl group
shydrox = (87.0, 'J/mol/K') # entropy of a surface hydroxyl group
####################### BULK PHASES ####################################
# First we'll define the bulk (i.e. 3D) phases - a gas, a metal, and
# an oxide.
#------------------------------------------------------------------
#
# Gas phase.
#
#------------------------------------------------------------------
# The gas contains only the minimum number of species needed to model
# operation on hydrogen. The species definitions are imported from
# gri30.cti. The initial composition is set to hydrogen + 5% water, but
# usually this is reset in the program importing this definition.
#
ideal_gas(name = "gas",
elements = " H O N",
species = "gri30: H2 H2O N2 O2",
transport = "Mix",
initial_state = state( temperature = tt,
pressure = OneAtm,
mole_fractions = 'H2:0.95, H2O:0.05'))
#------------------------------------------------------------------
#
# Bulk solid metal phase.
#
#------------------------------------------------------------------
#
# This phase will be used for the electrodes. All we need is
# a source/sink for electrons, so we define this phase as only
# containing electrons. Note that the 'metal' entry type requires
# specifying a density, but it is not used in this simulation and
# therefore is arbitrary.
#
metal(name = "metal",
elements = "E",
species = "electron",
density = (9.0, 'kg/m3'),
initial_state = state( temperature =tt,
mole_fractions = 'electron:1.0'))
# The electron is set to have zero enthalpy and entropy. Therefore,
# the chemical potential of the electron is zero, and the
# electrochemical potential is simply -F * phi, where phi is the
# electric potential of the metal. Note that this simple model is
# adequate only because all we require is a reservoir for electrons;
# if we wanted to do anything more complex, like carry out energy or
# charge balances on the metal, then we would require a more complex
# model. Note that there is no work function for this metal.
species( name = "electron", atoms = "E:1",
thermo = const_cp(h0 = (0.0, 'kcal/mol')))
# Note: the "const_cp" species thermo model is used throughout this
# file (with the exception of the gaseous species, which use NASA
# polynomials imported from gri30.cti). The const_cp model assumes a
# constant specific heat, which by default is zero. Parameters that
# can be specified are cp0, t0, h0, and s0. If omitted, t0 = 300 K, h0
# = 0, and s0 = 0. The thermo properties are computed as follows: h =
# h0 + cp0*(t - t0), s = s0 + cp0*ln(t/t0). For work at a single
# temperature, it is sufficient to specify only h0.
#-------------------------------------------------------------------
#
# Bulk solid oxide electrolyte
#
#--------------------------------------------------------------------
# Here too, we create a very simple model for the bulk phase. We only
# consider the oxygen sublattice. The only species we define are a
# lattice oxygen, and an oxygen vacancy. Again, the density is a
# required input, but is not used here, so may be set arbitrarily.
incompressible_solid(name = "oxide_bulk",
elements = "O E",
species = "Ox VO**",
density = (0.7, 'g/cm3'),
initial_state = state( temperature = tt,
pressure = OneAtm,
mole_fractions = "Ox:0.95 VO**:0.05")
)
# The vacancy will be modeled as truly vacant - it contains no atoms,
# has no charge, and has zero enthalpy and entropy. This is different
# from the usual convention in which the vacancy properties are are
# expressed relative to the perfect crystal lattice. For example, in
# the usual convention, an oxygen vacancy has charge +2. But the
# convention we will use is that an oxygen ion has charge -2, and a
# vacancy has charge 0. It all works out the same, as long as we are
# consistent.
# A bulk lattice vacancy
species( name = "VO**", atoms = "",
thermo = const_cp(h0 = (0.0, 'kJ/mol')))
# A bulk lattice oxygen
species( name = "Ox", atoms = "O:1 E:2",
thermo = const_cp(h0 = hox, s0 = sox))
####################### SURFACE PHASES ####################################
#--------------------------------------------------
#
# Metal surface
#
#--------------------------------------------------
# The surface of a bulk phase must be treated like a separate phase, with its
# own set of species. Here we define the model for the metal surface.
# We allow the following species:
# (m) - an empty metal site
# H(m) - a chemisorbed H atom
# O(m) - a chemisorbed O atom
# OH(m) - a chemisorbed hydroxl
# H2O(m) - a physisorbed water molecule
# Notes:
# 1. The site density is in mol/cm2, since no units are specified and
# 'mol' and 'cm' were specified in the units directive above as the
# units for quantity and length, respectively.
# 2. The 'reactions' field specifies that all reaction entries in this file
# that have ID strings beginning with "metal-" are reactions belonging
# to this surface mechanism.
ideal_interface(name = "metal_surface",
elements = "H O",
species = " (m) H(m) O(m) OH(m) H2O(m) ",
site_density = 2.60e-9,
phases = 'gas',
reactions = ["metal-*"],
initial_state = state( temperature = 973.0,
coverages = '(m):0.5 H(m):0.5') )
species( name = "(m)", atoms = "",
thermo = const_cp(h0 = (0.0, 'kJ/mol'),
s0 = (0.0, 'J/mol/K')))
species( name = "H(m)", atoms = "H:1",
thermo = const_cp(h0 = (-35.0, 'kJ/mol'),
s0 = (37.0, 'J/mol/K')))
species( name = "O(m)", atoms = "O:1",
thermo = const_cp(h0 = (-220.0, 'kJ/mol'),
s0 = (37.0, 'J/mol/K')))
species( name = "OH(m)", atoms = "O:1, H:1",
thermo = const_cp(h0 = (-198.0, 'kJ/mol'),
s0 = (102.0, 'J/mol/K')))
species( name = "H2O(m)", atoms = "H:2, O:1",
thermo = const_cp(h0 = (-281.0, 'kJ/mol'),
s0 = (123.0, 'J/mol/K')))
# Surface reactions on the metal. We assume three dissociative
# adsorption reactions, and three reactions on the surface
# among adsorbates. All reactions are treated as reversible.
surface_reaction( "H2 + (m) + (m) <=> H(m) + H(m)",
stick(0.1, 0, 0), id = 'metal-rxn1')
surface_reaction( "O2 + (m) + (m) <=> O(m) + O(m)",
stick(0.1, 0, 0), id = 'metal-rxn2')
surface_reaction( "H2O + (m) <=> H2O(m)",
stick(1.0, 0, 0), id = 'metal-rxn3')
surface_reaction( "H(m) + O(m) <=> OH(m) + (m)",
[5.00000E+22, 0, 100.0], id = 'metal-rxn4')
surface_reaction( "H(m) + OH(m) <=> H2O(m) + (m)",
[5.00000E+20, 0, 40.0], id = 'metal-rxn5')
surface_reaction( "OH(m) + OH(m) <=> H2O(m) + O(m)",
[5.00000E+21, 0, 100.0], id = 'metal-rxn6')
#--------------------------------------------------------
#
# Oxide surface.
#
#--------------------------------------------------------
#H
# On the oxide surface, we consider four species:
# 1. (ox) - a surface vacancy
# 2. O''(ox) - a surface oxygen with charge -2
# 3. OH'(ox) - a surface hydroxyl with charge -1
# 4. H2O(ox) - physisorbed neutral water
ideal_interface(name = "oxide_surface",
elements = "O H E",
species = "(ox) O''(ox) OH'(ox) H2O(ox)",
site_density = 2.0e-9,
phases = 'gas oxide_bulk',
reactions = 'oxide-*',
initial_state = state( temperature = tt,
coverages = "O''(ox):2.0, (ox):0.0") )
# Note: hox, sox, hhydrox, and shydrox are defined near the top of
# this file.
# An oxygen ion at the surface, with charge = -2
species( name = "O''(ox)", atoms = "O:1 E:2",
thermo = const_cp(h0 = hox,
s0 = sox))
# An OH at the surface, with charge = -1
species( name = "OH'(ox)", atoms = "O:1 H:1 E:1",
thermo = const_cp(h0 = hhydrox,
s0 = shydrox))
# A surface vacancy in the oxygen sublattice
species( name = "(ox)", atoms = "",
thermo = const_cp(h0 = (0.0, 'kJ/mol'),
s0 = (0.0,'J/mol/K')))
species( name = "H2O(ox)", atoms = "H:2, O:1",
thermo = const_cp(h0 = (-265.0, 'kJ/mol'),
s0 = (98.0,'J/mol/K')))
# This reaction represents the exchange of a surface oxygen vacancy and
# a subsurface vacancy. The concentration of subsurface vacancies is
# fixed by the doping level. If this reaction is given a large rate,
# then the surface vacancies will stay in equilibrium with the bulk
# vacancies.
surface_reaction("(ox) + Ox <=> VO** + O''(ox)",
[5.0e8, 0.0, 0.0], id = "oxide-vac")
# Desorption of physisorbed water. This is made fast.
surface_reaction("H2O(ox) <=> H2O + (ox)",
[1.0e14, 0.0, (0.0, 'kJ/mol')], id = "oxide-water")
# chemisorption of water as surface hydroxyls. In reality, this
# reaction would surely be activated and have a lower pre-exponential
surface_reaction("H2O(ox) + O''(ox) <=> OH'(ox) + OH'(ox)",
[1.0e14, 0.0, (0.0, 'kJ/mol')], id = "oxide-oh")
####################### TRIPLE PHASE BOUNDARY #########################
# The triple phase boundary between the metal, oxide, and gas. A
# single species is specified, but it is not used, since all reactions
# only involve species on either side of the tpb. Note that the site
# density is in mol/cm. But since no reactions involve TPB species,
# this parameter is unused.
edge(name = "tpb",
elements = "H O",
species = "(tpb)",
site_density = 5.0e-17,
reactions = "edge-*",
phases = 'metal metal_surface oxide_surface',
initial_state = state( temperature = tt,
coverages = '(tpb):1.0 ') )
# dummy species
species( name = "(tpb)", atoms = "")
# Here we define two charge transfer reactions. Both reactions are
# reversible, and can be used to model either anodes or cathodes
# (although real anodes and cathodes would usually have different
# reaction mechanisms, except in a symmetric cell).
# in this reaction, a proton from the metal crosses the TPB to the
# oxide surface to make a hydroxyl and deliver an electron to the
# metal.
edge_reaction("H(m) + O''(ox) <=> (m) + electron + OH'(ox)",
[5.0e13, 0.0, 120.0], beta = 0.5, id="edge-f2")
# in this reaction, an oxygen on the metal surface plus 2 electrons
# from the metal bulk fill a surface vacancy in the oxide lattice.
edge_reaction("O(m) + (ox) + 2 electron <=> (m) + O''(ox)",
[5.0e13, 0.0, 120.0], beta = 0.5, id="edge-f3")
# this reaction is commented out, but you can explore its effects by
# uncommenting it. Be careful, if you are not solving for the OH'
# concentration that the system does not become overdetermined
# (i.e. impossible for all reactions to be simultaneously in
# equilibrium). If this happens, the wrong OCVs will result.
#edge_reaction("H(m) + OH'(ox) <=> H2O(ox) + (m) + electron",
# [5.0e13, 0.0, 120.0], beta = 0.5, id="edge-f")

261
interfaces/cython/cantera/examples/surface_chemistry/sofc.py Normal file
View file

@ -0,0 +1,261 @@
"""
SOFC
This script implements a simple model of a solid oxide fuel cell. Unlike most
SOFC models, however, it does not use semi-empirical Butler-Volmer kinetics
for the charge transfer reactions, but uses elementary, reversible reactions
obeying mass-action kinetics for all reactions, including charge transfer. As
this script will demonstrate, this approach allows computing the OCV (it does
not need to be separately specified), as well as polarization curves.
NOTE: The parameters here, and in the input file sofc.cti, are not to be
relied upon for a real SOFC simulation! They are meant to illustrate only how
to do such a calculation in Cantera. While some of the parameters may be close
to real values, others are simply set arbitratily to give reasonable-looking
results.
It is recommended that you read input file sofc.cti before reading or running
this script!
"""
import cantera as ct
import math
import csv
import inspect
import os
ct.add_module_directory()
# parameters
T = 1073.15 # T in K
P = ct.one_atm
# gas compositions. Change as desired.
anode_gas_X = 'H2:0.97, H2O:0.03'
cathode_gas_X = 'O2:1.0, H2O:0.001'
# time to integrate coverage eqs. to steady state in
# 'advanceCoverages'. This should be more than enough time.
tss = 50.0
sigma = 2.0 # electrolyte conductivity [Siemens / m]
ethick = 5.0e-5 # electrolyte thickness [m]
TPB_length_per_area = 1.0e7 # TPB length per unit area [1/m]
def show_coverages(s):
"""Print the coverages for surface s."""
print('\n{}\n'.format(s.name))
cov = s.coverages
names = s.species_names
for n in range(s.n_species):
print('{:16s} {:13.4g}'.format(names[n], cov[n]))
def equil_OCV(gas1, gas2):
return (-ct.gas_constant * gas1.T *
math.log(gas1['O2'].X / gas2['O2'].X) / (4.0*ct.faraday))
def NewtonSolver(f, xstart, C=0.0):
"""
Solve f(x) = C by Newton iteration.
- xstart starting point for Newton iteration
- C constant
"""
f0 = f(xstart) - C
x0 = xstart
dx = 1.0e-6
n = 0
while n < 200:
ff = f(x0 + dx) - C
dfdx = (ff - f0)/dx
step = - f0/dfdx
# avoid taking steps too large
if abs(step) > 0.1:
step = 0.1*step/abs(step)
x0 += step
emax = 0.00001 # 0.01 mV tolerance
if abs(f0) < emax and n > 8:
return x0
f0 = f(x0) - C
n += 1
raise Exception('no root!')
#####################################################################
# Anode-side phases
#####################################################################
# import the anode-side bulk phases
gas_a, anode_bulk, oxide_a = ct.import_phases('sofc.cti',
['gas', 'metal', 'oxide_bulk',])
# import the surfaces on the anode side
anode_surf = ct.Interface('sofc.cti', 'metal_surface', [gas_a])
oxide_surf_a = ct.Interface('sofc.cti', 'oxide_surface', [gas_a, oxide_a])
# import the anode-side triple phase boundary
tpb_a = ct.Interface('sofc.cti', 'tpb', [anode_bulk, anode_surf, oxide_surf_a])
anode_surf.name = 'anode surface'
oxide_surf_a.name = 'anode-side oxide surface'
# this function is defined to use with NewtonSolver to invert the current-
# voltage function. NewtonSolver requires a function of one variable, so the
# other objects are accessed through the global namespace.
def anode_curr(E):
"""
Current from the anode as a function of anode potential relative to
electrolyte.
"""
# the anode-side electrolyte potential is kept at zero. Therefore, the
# anode potential is just equal to E.
anode_bulk.electric_potential = E
# get the species net production rates due to the anode-side TPB reaction
# mechanism. The production rate array has the values for the neighbor
# species in the order listed in the .cti file, followed by the tpb phase.
# Since the first neighbor phase is the bulk metal, species 0 is the
# electron.
w = tpb_a.net_production_rates
# the sign convention is that the current is positive when
# electrons are being delivered to the anode - i.e. it is positive
# for fuel cell operation.
return ct.faraday * w[0] * TPB_length_per_area
#####################################################################
# Cathode-side phases
#####################################################################
# Here for simplicity we are using the same phase and interface models for the
# cathode as we used for the anode. In a more realistic simulation, separate
# models would be used for the cathode, with a different reaction mechanism.
# import the cathode-side bulk phases
gas_c, cathode_bulk, oxide_c = ct.import_phases('sofc.cti',
['gas', 'metal', 'oxide_bulk'])
# import the surfaces on the cathode side
cathode_surf = ct.Interface('sofc.cti', 'metal_surface', [gas_c])
oxide_surf_c = ct.Interface('sofc.cti', 'oxide_surface', [gas_c, oxide_c])
# import the cathode-side triple phase boundary
tpb_c = ct.Interface('sofc.cti', 'tpb', [cathode_bulk, cathode_surf,
oxide_surf_c])
cathode_surf.name = 'cathode surface'
oxide_surf_c.name = 'cathode-side oxide surface'
def cathode_curr(E):
"""Current to the cathode as a function of cathode
potential relative to electrolyte"""
# due to ohmic losses, the cathode-side electrolyte potential is non-zero.
# Therefore, we need to add this potential to E to get the cathode
# potential.
cathode_bulk.electric_potential = E + oxide_c.electric_potential
# get the species net production rates due to the cathode-side TPB
# reaction mechanism. The production rate array has the values for the
# neighbor species in the order listed in the .cti file, followed by the
# tpb phase. Since the first neighbor phase is the bulk metal, species 0
# is the electron.
w = tpb_c.net_production_rates
# the sign convention is that the current is positive when electrons are
# being drawn from the cathode (i.e, negative production rate).
return -ct.faraday * w[0] * TPB_length_per_area
# initialization
# set the gas compositions, and temperatures of all phases
gas_a.TPX = T, P, anode_gas_X
gas_a.equilibrate('TP') # needed to use equil_OCV
gas_c.TPX = T, P, cathode_gas_X
gas_c.equilibrate('TP') # needed to use equil_OCV
phases = [anode_bulk, anode_surf, oxide_surf_a, oxide_a, cathode_bulk,
cathode_surf, oxide_surf_c, oxide_c, tpb_a, tpb_c]
for p in phases:
p.TP = T, P
# now bring the surface coverages into steady state with these gas
# compositions. Note that the coverages are held fixed at these values - we do
# NOT consider the change in coverages due to TPB reactions. For that, a more
# complex model is required. But as long as the thermal chemistry is fast
# relative to charge transfer, this should be an OK approximation.
for s in [anode_surf, oxide_surf_a, cathode_surf, oxide_surf_c]:
s.advance_coverages(tss)
show_coverages(s)
# find open circuit potentials by solving for the E values that give
# zero current.
Ea0 = NewtonSolver(anode_curr, xstart=-0.51)
Ec0 = NewtonSolver(cathode_curr, xstart=0.51)
print('\nocv from zero current is: ', Ec0 - Ea0)
print('OCV from thermo equil is: ', equil_OCV(gas_a, gas_c))
print('Ea0 = ', Ea0)
print('Ec0 = ', Ec0)
print()
# do polarization curve for anode overpotentials from -250 mV
# (cathodic) to +250 mV (anodic)
Ea_min = Ea0 - 0.25
Ea_max = Ea0 + 0.25
csvfile = open('sofc.csv', 'w')
writer = csv.writer(csvfile)
writer.writerow(['i (mA/cm2)', 'eta_a', 'eta_c', 'eta_ohmic', 'Eload'])
# vary the anode overpotential, from cathodic to anodic polarization
for n in range(100):
Ea = Ea_min + 0.005*n
# set the electrode potential. Note that the anode-side electrolyte is
# held fixed at 0 V.
anode_bulk.electric_potential = Ea
# compute the anode current
curr = anode_curr(Ea)
# set potential of the oxide on the cathode side to reflect the ohmic drop
# through the electrolyte
delta_V = curr * ethick / sigma
# if the current is positive, negatively-charged ions are flowing from the
# cathode to the anode. Therefore, the cathode side must be more negative
# than the anode side.
phi_oxide_c = -delta_V
# note that both the bulk and the surface potentials must be set
oxide_c.electric_potential = phi_oxide_c
oxide_surf_c.electric_potential = phi_oxide_c
# Find the value of the cathode potential relative to the cathode-side
# electrolyte that yields the same current density as the anode current
# density
Ec = NewtonSolver(cathode_curr, xstart=Ec0+0.1, C=curr)
cathode_bulk.electric_potential = phi_oxide_c + Ec
# write the current density, anode and cathode overpotentials, ohmic
# overpotential, and load potential
writer.writerow([0.1*curr, Ea - Ea0, Ec - Ec0, delta_V,
cathode_bulk.electric_potential -
anode_bulk.electric_potential])
print('polarization curve data written to file sofc.csv')
csvfile.close()

2
interfaces/cython/setup.py.in
View file

@ -25,4 +25,4 @@ setup(name="Cantera",
ext_modules = exts,
package_data = {'cantera.data': ['*.*'],
'cantera.test.data': ['*.*'],
'cantera.examples': ['*/*.py']})
'cantera.examples': ['*/*.*']})