diff --git a/interfaces/cython/cantera/examples/dusty_gas.py b/interfaces/cython/cantera/examples/dusty_gas.py new file mode 100644 index 000000000..52115a64c --- /dev/null +++ b/interfaces/cython/cantera/examples/dusty_gas.py @@ -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)) diff --git a/interfaces/cython/cantera/examples/multiphase/plasma_equilibrium.py b/interfaces/cython/cantera/examples/multiphase/plasma_equilibrium.py new file mode 100644 index 000000000..17560a44e --- /dev/null +++ b/interfaces/cython/cantera/examples/multiphase/plasma_equilibrium.py @@ -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() diff --git a/interfaces/cython/cantera/examples/onedim/flame_fixed_T.py b/interfaces/cython/cantera/examples/onedim/flame_fixed_T.py new file mode 100644 index 000000000..91f802a28 --- /dev/null +++ b/interfaces/cython/cantera/examples/onedim/flame_fixed_T.py @@ -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() diff --git a/interfaces/cython/cantera/examples/onedim/tdata.dat b/interfaces/cython/cantera/examples/onedim/tdata.dat new file mode 100644 index 000000000..aec20fa7b --- /dev/null +++ b/interfaces/cython/cantera/examples/onedim/tdata.dat @@ -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 + diff --git a/interfaces/cython/cantera/examples/surface_chemistry/catalytic_combustion.py b/interfaces/cython/cantera/examples/surface_chemistry/catalytic_combustion.py new file mode 100644 index 000000000..cc6434399 --- /dev/null +++ b/interfaces/cython/cantera/examples/surface_chemistry/catalytic_combustion.py @@ -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) diff --git a/interfaces/cython/cantera/examples/surface_chemistry/diamond_cvd.py b/interfaces/cython/cantera/examples/surface_chemistry/diamond_cvd.py new file mode 100644 index 000000000..4a0ecf30d --- /dev/null +++ b/interfaces/cython/cantera/examples/surface_chemistry/diamond_cvd.py @@ -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') diff --git a/interfaces/cython/cantera/examples/surface_chemistry/sofc.cti b/interfaces/cython/cantera/examples/surface_chemistry/sofc.cti new file mode 100644 index 000000000..b605815d7 --- /dev/null +++ b/interfaces/cython/cantera/examples/surface_chemistry/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 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") + + + + + + diff --git a/interfaces/cython/cantera/examples/surface_chemistry/sofc.py b/interfaces/cython/cantera/examples/surface_chemistry/sofc.py new file mode 100644 index 000000000..e516920ac --- /dev/null +++ b/interfaces/cython/cantera/examples/surface_chemistry/sofc.py @@ -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() diff --git a/interfaces/cython/setup.py.in b/interfaces/cython/setup.py.in index c9842037f..451e9821a 100644 --- a/interfaces/cython/setup.py.in +++ b/interfaces/cython/setup.py.in @@ -25,4 +25,4 @@ setup(name="Cantera", ext_modules = exts, package_data = {'cantera.data': ['*.*'], 'cantera.test.data': ['*.*'], - 'cantera.examples': ['*/*.py']}) + 'cantera.examples': ['*/*.*']})