297 lines
9.5 KiB
Python
297 lines
9.5 KiB
Python
# 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()
|