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.