In this section, we will follow the main function, in the file sofcFoam.C, through the model execution. Like all OpenFOAM applications, the model begins by including the OpenFOAM src files setRootCase.H to check the case path, and createTime.H to read the system/controlDict file and instantiate the Time object runTime. This is followed by the creation of meshes, reading of properties, and creation of fields for the global cell mesh and the region meshes.
The first mesh to be created is the mesh for the entire cell. This is accomplished by the included file $FOAM_SRC/OpenFOAM/lnInclude/createMesh.H.. Mesh data is obtained from the constant/polyMesh directory. Global cell properties, or parameters, as listed in Table 1, above, are read from case file constant/cellProperties into model variables by appSrc file readCellProperties.H. Fields for the entire cell are created as IOobjects in appSrc file createCellFields.H. Some of the IOobjects have the MUST_READ attribute and must have a file in the starting time directory. The file name must be the same as the name specified in the IOobject. Others have the READ_IF_PRESENT attribute. The corresponding files will be read if they are present in the starting time directory. Still others have the NO_READ attribute and any file with their name will be ignored at field creation. Fields with IOobjects having the AUTO_WRITE attribute will write a file in the output time directories, whereas those with the NO_WRITE attribute aren’t written.
Similarly, meshes, properties, and fields are established for the regions interconnect0, air, electrolyte, fuel, and interconnect1. The meshes are created by the appSrc files create<Region>Mesh.H. These files specify the location of the region’s polyMesh directory and also create face-, cell-, and patch-maps to the global mesh. Constant properties are read by read<Region>Properties.H from case files constant/<region>/<region>Properties, and fields are read from case files 0/<region>/<fieldName> by appSrc file create<Region>Fields.H. Note that there are no fields on either of the interconnect meshes, and both interconnects are assumed to have the same properties.
A number of air and fuel fields are specific to the species that comprise the fluid. The key appSrc files are createAirSpecies.H and createFuelSpecies.H. The key case files are constant/air/sofcSpeciesProperties and constant/fuel/sofcSpeciesProperties. See these for details of input data format.
As just indicated, the air species names and associated data are read from the constant/air/sofcSpeciesProperties file at run time by appSrc file createAirSpecies.H. A speciesTable object airSpeciesNames is instantiated from the list of species following the key word speciesList in the sofcSpeciesProperties file.
A pointerList, airSpecies, of pointers to sofcSpecie objects is created. The i-th air species can then be referenced as airSpecies[i]. For each specie in the air mixture, the specie properties of name, molar weight, molar charge (for Faraday’s law), reaction sign (produced=1, non-reacting=0, consumed=-1), enthalpy of formation, and standard entropy are stored in an sofcSpecie object, and can be accessed by class functions name(), MW(), ne(), rSign(), hForm() and sForm(), respectively. The sofcSpecie class can be found in src/libSrc.
After reading the species data, one of the species is designated as airInertSpecie. It is found after the keyword inertSofcSpecie in the sofcSpeciesProperties file. Note that the inertSpecie may well be chemically inert (non-reacting), but need not be. Here “inert” means that the mass fraction of the specie will not be computed from a partial differential equation, but rather calculated by adding the mass fractions of the other components and subtracting that sum from 1. The “inert” nomenclature follows that of OpenFOAM’s thermophysicalModels.
The sofcSpeciesProperties file concludes with a toddYoung dictionary of Todd-Young (2002) polynomial coefficients for molar isobaric heat capacity. These are used to create a pointerList, molarCpAir, of pointers to polyToddYoung objects. The polyToddYoung class has functions to evaluate the polynomial and to also evaluate definite integrals which correspond to enthalpy and entropy in the case of isobaric heat capacity (see src/libSrc/polyToddYoung). Thus, eg, we can evaluate the isobaric heat capacity of the i-th air species at ambient temperature with the expression molarCpAir[i].polyVal(Tair.internalField()).
The specie names are used to create a pointerList, Yair, to mass fraction fields with names of the form Ysp, where “sp” is one of the specie names. There is one such file for each specie in the mixture. The mass fraction field IOobjects are MUST_READ, so must exist as files in the starting time directory.
Mole fraction fields are calculated from the mass fraction fields. Their names have the form Xsp. The read/write attributes are NO_READ, AUTO_WRITE. Again we have a pointerList, Xair, so that Xair[i] references the mole fraction field of the i-th air species.
Finally, createAirSpecies.H establishes a pointerlist, diffSpAir, to scalar fields for the diffusivities of the individual species in the mixture. These have READ_IF_PRESENT and AUTO_WRITE attributes, and are created with initial value 1.
A completely analogous discussion applies to the fuel side.
The chemical species involved in the reaction and their stoichiometric coefficients are listed after the rxnSpecies keyword in the constant/rxnProperties case file. The species “e” must always be in the list. Its coefficient is the number of moles of electrons transferred in the reaction. The list is read by appSrc file readRxnProperties.H, where a hash table, rxnSpCoef, is created. Thus a specie’s stoichiometric coefficient can be obtained via its name, e.g., rxnSpCoef["O2"], or rxnSpCoef[airSpecies[i].name()].
For convenience, global variables are established for the IDs (indices) of a number of patches that are frequently referenced. These patches were assigned to variable names in appSrc file readCellProperties.H, reading from case file constant/cellProperties. Their IDs are found and assigned in appSrc file setGlobalPatchIds.H.
Electrolyte thickness is used in the calculation of the electrolyte’s volumetric heat source term for the energy equation. It is calculated, in appSrc file electrolyteThickness.H, as the electrolyte volume divided by the average area of the anode and cathode interfaces. Average area is used to take into account cylindrical cells.
Electrochemistry is assumed to occur within the electrolyte. Mole fraction patch fields of oxidants on the air side (cathode) interface and of hydrogen and water on the fuel side (anode) must be therefore interpolated to fields in the electrolyte. The resulting current density is calculated on the fuel side of the anode interface and must be interpolated to the electrolyte side in order to calculate the electrolyte’s volumetric heat source terms. These interpolations are carried out by OpenFOAM patchToPatchInterpolation objects, created in appSrc file createPatchToPatchInterpolation.H.
The diffusivity models are created in appSrc file createDiffusivityModels.H. For both air and fuel, a pointerList of diffusivity models is declared. On the air side, this is airDiffModels. Creation of a new diffusivity model requires a scalar field for the calculated diffusivity values, a label list of cells for which the diffusivity is to be calculated (normally corresponding to a cell zone), and a dictionary specifying the model type along with the corresponding type-specific parameters.
On the air side, a scalarField, airDiff, is used by each airDiffusivityModel[m] to return its calculated values. Then, a diffusivity model is established for each air zone and a pointer to it is added to the pointerList. There is one model for the entire air zone. Its dictionary is located in the constant/air/airProperties case file. There is also one model for each porous zone within the air region, with dictionaries in the corresponding zone in the case file constant/air/porousZones. The cellZone labelLists are found in the polyMesh directory for the air region.
An analogous story holds on the fuel side. How the SOFC model uses the diffusivity models will be described below.
A number of calculations are repeated until convergence.
The fluid regions need local temperature fields in order to calculate local fluid density. The mapping is done in appSrc file mapFromCell.H, using the cell maps established during mesh creation.
Fluid densities ρ are calculated from pressure p, temperature T, constituent mass fractions yi and molar weights Mi as
$$\rho =\frac{p}{RT\mathop{\sum }^{}\left( {{y}_{i}}/{{M}_{i}} \right)}$$
where R is the universal gas constant. The calculations are coded in appSrc files rhoAir.H and rhoFuel.H.
Pressure and momentum are solved using the PISO iteration in appSrc files solveAir.H and solveFuel.H. Following the solution, the Reynolds numbers are calculated in appSrc file ReynoldsNumber.H. These are informative only, and require the input of hydraulic diameter of the gas channels, in case files constant/air/airProperties and constant/fuel/fuelProperties.
The diffusivities are calculated in appSrc files diffusivityAir.H and diffusivityFuel.H. The following discussion refers only to air, but applies equally to fuel. Of interest is the diffusivity, D(a), of specie a in the mixture. A diffusivity field is thus calculated for each specie, except for the background “inert” specie. A specie a can be given a fixed diffusivity value using the fixedDiffusivity model, or its diffusivity in the mixture, D(a), can be modelled from the pairwise binary diffusivities D(a,b) of specie a diffusing in specie b.
The sofcFoam model combines pairwise binary diffusivities following Wilke (1950):
$$D\left( a \right)=~\frac{1-~{{x}_{a}}}{{{S}_{a}}}$$
where
$${{S}_{a}}~=\underset{b\ne a}{\mathop \sum }\,\frac{{{x}_{b}}}{D\left( a,b \right)}$$
with mole fractions x and specie indices a,b.
The calculation proceeds by first fixing specie a, stepping through the species b and accumulating the results to obtain the sum Sa, which is then used to get D(a). The computed diffusivity is stored in volScalarField diffSpAir[a].
The diffusivity calculations have the following algorithmic structure:
forAll(airSpecies, a)
{
if(airSpecies[a].name() != airInertSpecie)
{
forAll(airDiffModels, m)
{
if(airDiffModels[m]->isFixed())
{
//obtain fixed diffusivity
}
else if(!airDiffModels[m]->isBinary())
{
Error: must have fixed or binary: exit
}
else
{
initialize sumA = 0
forAll(airSpecies,b)
{
if (b != a)
{
set specie b in airDiffModel[m]
calculate binary diffusivity
accumulate sumA
}
}
obtain diffusivity in mixture using sumA
} //isBinary
} //m
} //!inert
} //a
Mass fractions are computed in appSrc files YairEqn.H and YfuelEqn.H. For each specie other than the background “inert” specie, a partial differential equation is solved for Yair[i], where i is the specie index. The mass fraction of the “inert” specie is obtained by subtracting the sum of all the other species’ mass fractions from 1.
Electrochemistry is assumed assumed to occur on the anode interface with the electrolyte. It is calculated in appSrc file solveElectrochemistry.H. Global temperature is interpolated to the fuel mesh anode patch. Reactant oxidant species mole fractions are computed on the air mesh cathode patch and interpolated to the fuel mesh anode patch. Reactant and product fuel species are also calculated on the fuel mesh anode patch.
With all fields defined on the anode patch, the mole fractions, temperature, and specie properties are used to calculate the Nernst potential, E, via included appSrc file NernstEqn.H. Area specific resistance, R, is modelled by a function in included appSrc file ASRfunction.H. Current density i is then calculated as
$$i=\frac{E-V-\eta}{R}$$
where V and \(\eta\) are the present values of voltage and activation overpotential. The voltage is subsequently corrected using the present and prescribed mean current densities.
Electrochemical heating is calculated in included appSrc file electrochemicalHeating.H. Specie enthalpies are calculated and combined with enthalpy of formation and Joule heating in the electrolyte volume.
Species electrochemical mass fluxes are calculated and used to set Neumann boundary conditions on the cathode and anode patches for the calculation of air and fuel mass fractions, and Dirichlet conditions for the air and fuel velocities.
For the reaction
$$\text{rxn:}\ \ \sum\limits_{i}{{}}{{a}_{i~}}{{R}_{i~}}=~\underset{j}{\mathop \sum }\,{{b}_{j}}{{P}_{j}}$$
having reactants Ri with stoichiometric coefficients ai and products Pj with stoichiometric coefficients bj, the Nernst potential, E, is calculated as
$$E=~{{E}_{0}}-~\frac{RT}{zF}\text{ln}Q$$
where R is the universal gas constant, F is Faraday’s constant, z is the number [moles] of electrons transferred,
$$Q=~\frac{\mathop{\prod }_{j}{{\left[ {{P}_{j}} \right]}^{{{b}_{j}}}}}{\mathop{\prod }_{i}{{\left[ {{R}_{i}} \right]}^{{{a}_{i}}}}}$$
with [ denoting mole fraction, and
$${{E}_{0}}=~-\frac{\Delta {{G}_{\text{rxn }\!\!~\!\!\text{ }\!\!~\!\!\text{ }}}}{F}=~-\frac{\left( \Delta {{H}_{\text{rxn}}}-T\Delta {{S}_{\text{rxn}}} \right)}{F}$$
where
$$\Delta {{H}_{\text{rxn}}}~=~\underset{j}{\mathop \sum }\,{{b}_{j}}\Delta H\left( {{P}_{j}} \right)-\underset{i}{\mathop \sum }\,{{a}_{i}}\Delta H({{R}_{i}})$$
and
$$\Delta {{S}_{\text{rxn}}}~=~\underset{j}{\mathop \sum }\,{{b}_{j}}\Delta S\left( {{P}_{j}} \right)-\underset{i}{\mathop \sum }\,{{a}_{i}}\Delta S({{R}_{i}})$$
For species X,
$$\Delta H\left( X \right)=~\underset{{{T}_{0}}}{\overset{T}{\mathop \int }}\,{{C}_{p,X}}\left( T \right)dT$$
and,
$$\Delta S\left( X \right)=~\underset{{{T}_{0}}}{\overset{T}{\mathop \int }}\,\frac{{{C}_{p,X}}\left( T \right)dT}{T}$$
where T0 and T are reference and ambient temperatures, respectively. Todd-Young polynomials for molar heat capacity are used to evaluate the and integrals, using polyToddYoung class functions polyInt and polyIntS, respectively.
The cell activation overpotential is presumed to be composed of anodic and cathodic contributions: \[\eta ={{\eta }_{a}}+{{\eta }_{c}}\] Defining the electrode overpotential, \({{\eta }_{e}}={{\eta }_{a}}\) at the anode and \({{\eta }_{e}}={{\eta }_{c}}\) at the cathode. Values of \({{\eta }_{e}}\left( i \right)\) at both electrodes are obtained using a root finder based on Ridder's method coded in appSrc file, riddersRoot.C, to obtain numerical solutions of the Butler-Volmer equation \[{i}={{{i}}_{0}}\left[ \exp \left( A{{\eta }_{e}} \right)-\exp \left( B{{\eta }_{e}} \right) \right]\] where \[{A=2{{\alpha }_{e}}F}/{RT}\;\]\[{B=-2\left( 1-{{\alpha }_{e}} \right)F}/{RT}\;\] at each electrode. The exchange current densities are given by; \[{{{i}}_{0,c}}={{\gamma }_{c}}p_{{{\text{O}}_{2}}}^{m}\exp \left( -\frac{{{E}_{act,c}}}{RT} \right)\] \[{{{i}}_{0,a}}={{\gamma }_{a}}p_{{{\text{H}}_{2}}}^{a}p_{{{\text{H}}_{\text{2}}}\text{O}}^{b}\exp \left( -\frac{{{E}_{act,a}}}{RT} \right)\]
The reaction orders, a, b, m, the pre-constants \({{\gamma }_{a}}\), \({{\gamma }_{c}}\), activation energies Eact,a , Eact,c are prescribed together with the transfer coefficients, \({{\alpha }_{e}}\), \({{\alpha }_{c}}\) in the dictionary activationParameters, located in constant/electrolyte. The prescibed values are those recommended in Leonide at el. (2009). Note the partial pressures are in Atmospheres.
Heating sources for the energy solution are comprised of enthalpy of reactions (assumed to occur on the electrolyte/fuel interface), enthalpy changes of gaseous species from reference to ambient temperature, and resistive heating. It is assumed that there is no reaction on the air side. The volumetric source terms are ascribed to the electrolyte region.
Enthalpy of formation and enthalpy changes, normalized by number of contributing electrons, are calaculated by specie, accumulating products and reactants separately for each of air and fuel. The separate (normalized) enthalpy accumulations are then combined according to
\[{{h}_{src}}=\left( {{h}_{form}}+{{h}_{P,fuel}}-{{h}_{R,fuel}}+{{h}_{P,air}}-{{h}_{R,air}} \right)\frac{i}{F{{l}_{e}}}\]
where i is current density, F is Faraday’s constant, le is the electrolyte thickness, hform is accumulated enthalpy of formation, hP,fuel is accumulated enthalpy change of products on the fuel side, hR,fuel is accumulated enthalpy change of reactants on the fuel side, and similarly for the air side. The final source for the energy equation then becomes
$${{S}_{\text{energy}}}=-{{h}_{\text{src}}}-\frac{iV}{{{l}_{e}}}$$
An electrochemical mass flux is calculated for each specie taking part in the reaction, in accordance with Faraday’s law, and taking into account whether the specie is consumed or produced:
$${{\dot{m}}^{''}}=~\pm \frac{Mi}{\nu F}$$
where M is molar mass, i is current density, ν is valence, and F is Faraday’s constant. Air species are treated separately from fuel species, and separate sums of air and fuel fluxes are accumulated for later use. The plus sign is used for products and the negative sign for reactants. On the air side, the above calculations are coded at lines 164 to 177 of solveElectrochemistry.H.
Mass fractions Y of all species, except the background “inert” specie, are found by solving a partial differential equation. Mass is transferred through the electrode boundaries, so these boundary fields are cast to a fixedGradient type, to which a (generally non-uniform) gradient value is assigned. The mass fraction gradient of a specie i will be due to both mass flux of i and to the mass fluxes of the other species. Letting Yi be the mass fraction boundary field of specie i on the electrode boundary, we have
$$\frac{\partial {{Y}_{i}}}{\partial n}={1\over{\rho D}}\left(\dot{m}_{i}^{''}\left( 1-{{Y}_{i~}} \right)-{{Y}_{i}}\underset{j\ne i}{\mathop \sum }\,\dot{m}_{j}^{''}\right)$$
The mass flux sums are used to calculate Dirichlet velocity boundary conditions for air on the cathode and for fuel on the anode. Note that the calculations take place with fields defined on the anode interface, so patchToPatchInterpolation to the air side is required. Then, we have for U on the boundary
\[{{\left. \mathbf{U} \right|}_{\text{electrode}~\text{ }}}=-\frac{\sum{{{\dot{m}}''}}}{\rho }\frac{\mathbf{A}}{\left\| \mathbf{A} \right\|}~\]
where \(\mathop{\sum }^{}{{\dot{m}}^{''}}\) is the accumulated mass flux, ρ is the fluid density, and \({{\mathbf{A}}/{\left\| \mathbf{A} \right\|}\;}\) is the unit outward normal.
In order to solve the energy equation, regional values of thermal conductivity, heat capacity, density, velocity, etc., are required on the global mesh. Prior to any mapping, the internal values of the recipient fields on the global mesh are reset to zero by appSrc file mapToCell.H. Regional fields are then mapped to the global mesh by appSrc files map<Region>ToCell.H for each of the regions air, fuel, electrolyte, interconnect0 and interconnect1.
In the case of the solids, there is no convective heat transfer, so velocity and heat capacity for the corresponding space within the global mesh will simply be set to zero. Density and thermal conductivity are assumed uniform and are given by the values in their properties files. In the case of the electrolyte, the volumetric energy source terms computed there are also mapped to a global energy source field.
For the fluids, we have both convective and diffusive heat transfer. Mass based heat capacity of the mixture, cp, is calculated as a linear combination of the specie molar heat capacities, Cp,i, divided by their molar mass, Mi, using mass fraction, Yi, as the linear coefficients:
$${{c}_{p}}=\underset{i}{\mathop \sum }\,\frac{{{Y}_{i}}{{C}_{p,i}}}{{{M}_{i}}}$$
Thermal conductivity is assumed to be uniform within each fluid zone. In the fluid channels it takes the value read from the fluid properties file, whereas in the porous zones it is a linear combination of the channel value and the porous zone value weighted by porosity. Porosity and porous zone thermal conductivity are obtained from the porous zone dictionary.
Rather than mapping velocity directly, it is the fluxes \(\phi =\rho \mathbf{U}\centerdot \mathbf{dA}\) (computed during the pressure-velocity calculation such that mass conservation is satisfied) that are required for the finite volume equations. The fluxes form a surface field on the mesh faces and are mapped onto the global faces corresponding to the regional meshes interior faces and also to the regional mesh boundary patches. At this stage, we have fluxes on the fluid sides of the fluid/electrolyte interfaces, but not on the electrolyte side, the electrolyte being a solid. This can result in heat leaving or entering a fluid through its boundary with the electrolyte, without a corresponding gain or loss within the electrolyte. This situation is remedied by zeroing all fluxes that touch the electrolyte and introducing an additional compensating source term into the energy equation. The zeroing of fluxes is done in mapElectrolyteToCell.H. The source term is equivalent to the continuity error introduced by the zeroing of the fluxes on the electrolyte boundaries.
After the regional fields are mapped to the corresponding locations in the global mesh, the global energy is computed from a partial differential equation in appSrc file solveEnergy.H. The equation is essentially
\[\operatorname{div}\left( \rho \mathbf{U}{{c}_{p}}T \right)-{{S}_{\phi }}-\operatorname{div}\left( k\operatorname{grad}T \right)={{S}_{q}}\]
where \({{S}_{q}}\) is a volumetric heat source and \({{S}_{\phi }}\) is the source due to the zeroed fluxes on the electrolyte boundary faces. The term \(-{{S}_{\phi }}\) appears in the discrete equation as
fvm::SuSp(-fvc::div(rhoCpPhiCell), Tcell)
where rhoCpPhiCell is equal to the face interpolated cp multiplied by the surface scalar field phiCell. Note that phi already incorporates rho. The finite volume matrix operator SuSp linearizes the source term and adds the linear part to the matrix diagonal or to the source in such a way as to maximize diagonal dominance of the matrix.