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

added fuel cell demo

Browse source
This commit is contained in:
Dave Goodwin 2006-11-09 13:04:39 +00:00
parent dde09cbd38
commit 371a3c4198
2 changed files with 654 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

357
Cantera/python/examples/fuel_cells/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 reservior 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 dissociave
# 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.
#
#--------------------------------------------------------
#
# 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, andd 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 exhange 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")

297
Cantera/python/examples/fuel_cells/sofc.py Normal file
View file

@ -0,0 +1,297 @@
# 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!
#---------------------------------------------------------------------
from Cantera import *
import math
#--------------------------------------------------
#
# parameters
#
#--------------------------------------------------
tc = 800.0 # T in C
temp = tc + 273.15
pres = OneAtm
# 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
# electrolyte conductivity
sigma = 2.0 # Siemens / m
# electrolyte thickness
ethick = 5.0e-5 # m
# TPB length per unit area
TPB_length_per_area = 1.0e7 # per meter
#----------------------------------------------------
#
# utility functions
#
#----------------------------------------------------
def showCoverages(s):
"""Print the coverages for surface s."""
print '\n '+s.name() + '\n'
cov = s.coverages()
names = s.speciesNames()
nsp = len(names)
for n in range(nsp):
print '%16s %13.4g ' % (names[n], cov[n])
def equil_OCV(gas1, gas2):
return -GasConstant*gas1.temperature()*math.log(gas1.moleFraction('O2')/
gas2.moleFraction('O2'))/(4.0*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
xlast = 999.0
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
xlast = x0
f0 = f(x0) - C
n += 1
raise 'no root!'
#####################################################################
#
# Anode-side phases
#
#####################################################################
# import the anode-side bulk phases
gas_a, anode_bulk, oxide_a = importPhases('sofc.cti',
['gas', 'metal', 'oxide_bulk',])
# import the surfaces on the anode side
anode_surf = importInterface('sofc.cti','metal_surface',[gas_a])
oxide_surf_a = importInterface('sofc.cti','oxide_surface',[gas_a, oxide_a])
# import the anode-side triple phase boundary
tpb_a = importEdge('sofc.cti', 'tpb', [anode_bulk, anode_surf, oxide_surf_a])
anode_surf.setName('anode surface')
oxide_surf_a.setName('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.setElectricPotential(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.netProductionRates()
# 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 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 = importPhases('sofc.cti',
['gas',
'metal',
'oxide_bulk',])
# import the surfaces on the cathode side
cathode_surf = importInterface('sofc.cti','metal_surface',[gas_c])
oxide_surf_c = importInterface('sofc.cti','oxide_surface',[gas_c, oxide_c])
# import the cathode-side triple phase boundary
tpb_c = importEdge('sofc.cti', 'tpb', [cathode_bulk,
cathode_surf, oxide_surf_c])
cathode_surf.setName('cathode surface')
oxide_surf_c.setName('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.
ee = E + oxide_c.electricPotential()
cathode_bulk.setElectricPotential(ee)
# 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.netProductionRates()
# the sign convention is that the current is positive when electrons
# are being drawn from the cathode (i.e, negative production rate).
return -Faraday * w[0] * TPB_length_per_area
# initialization
# set the gas compositions, and temperatures of all phases
gas_a.set(T = temp, P = pres, X = anode_gas_X)
gas_a.equilibrate('TP') # needed to use equil_OCV
gas_c.set(T = temp, P = pres, X = 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.setTemperature(temp)
# 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.advanceCoverages(tss)
showCoverages(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
file = open('sofc.csv','w')
writeCSV(file,['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.setElectricPotential(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.setElectricPotential(phi_oxide_c)
oxide_surf_c.setElectricPotential(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.setElectricPotential(phi_oxide_c + Ec);
# write the current density, anode and cathode overpotentials,
# ohmic overpotential, and load potential
writeCSV(file,[0.1*curr, Ea - Ea0, Ec - Ec0, delta_V,
cathode_bulk.electricPotential()
- anode_bulk.electricPotential()])
print 'polarization curve data written to file sofc.csv'
file.close()