diff --git a/Cantera/python/examples/fuel_cells/sofc.cti b/Cantera/python/examples/fuel_cells/sofc.cti new file mode 100644 index 000000000..df4f3daad --- /dev/null +++ b/Cantera/python/examples/fuel_cells/sofc.cti @@ -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") + + + + + + diff --git a/Cantera/python/examples/fuel_cells/sofc.py b/Cantera/python/examples/fuel_cells/sofc.py new file mode 100644 index 000000000..d523b313b --- /dev/null +++ b/Cantera/python/examples/fuel_cells/sofc.py @@ -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()