Energy solution

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.