For the reaction

$$\text{rxn:}\ \ \sum\limits_{i}{{}}{{a}_{i~}}{{R}_{i~}}=~\underset{j}{\mathop \sum }\,{{b}_{j}}{{P}_{j}}$$

having reactants R_i with stoichiometric coefficients a_i and products P_j with stoichiometric coefficients b_j, 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 T₀ 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 E_act,a , E_act,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, l_e is the electrolyte thickness, h_form is accumulated enthalpy of formation, h_P,fuel is accumulated enthalpy change of products on the fuel side, h_R,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 Y_i 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.