**2. Problem Description and Governing Equations**

The functional layers in a PEMFC are the membrane, the gas diffusion backing layers (either a single GDL or an assembly of GDL and microporous layers), bipolar plates (BPs), and CL. Studies have further shown the beneficial effect of inserting a microporous layer (MPL) between the GDL and CL to improve the effective liquid and gas-phase transport properties of the multilayered system for better water management characteristics [30]. In the absence of contaminants, hydrogen diffuses through the PEMFC anode side GDL and MPL until the CL to produce protons and electrons by electrochemical reaction (Figure 4). The membrane selectively transports protons to the cathode, where electrochemical reaction with oxygen produces water.

**Figure 4.** Schematic of PEMFC working principles and components.

Nafion is the common membrane material, which has a thickness ranging from 0.01 to 0.1 mm. Carbon-supported catalyst and ionomer porous composites are used for the CL with thicknesses spanning from 100 nm to 0.05 mm. GDL and MPL are porous carbon structures coated by polytetrafluoroethylene (PTFE), with the porosity of the MPL usually being lower to improve water management. In the present study, only variations and fluxes along the x direction were considered, as the corresponding gradients in composition and potentials were expected to be the largest. This assumption was considered acceptable to predict the output electrical power of the system with sufficient accuracy for the objective of the present study.

The PEMFC model was one-dimensional along the z-direction. It assumed incompressible flows, ideal gases, a laminar regime and single-gas phase transport in the voids. Additionally, the concentration of protons in the membrane was considered constant, that is, the potential gradient was the driving force for the flux of protons. The input flow was assumed to be pure hydrogen, that is, potential contamination of the cell due to trace elements in practice was neglected.

The geometry of the heat recovery module modeled in the current study corresponds to that analyzed by Fernández-Yañez [28]. The present study was indeed a feasibility investigation using a configuration for which data for validation were available. The module consisted of a TEG unit implemented in a water/water co-flow heat exchanger: the flow of cooling water exiting the PEMFC was the input for the water exchanger in the TEG unit (straight internal channels), and the cooling water from the environment to TEG (serpentine channel) acted as the cold source. Figure 5 shows the location of the suggested TEG unit to be integrated in the cooling channel of the PEMFC as a cooling method on medium-scale stacks given by Figure 3.

**Figure 5.** The suggested location of the TEG unit in the PEMFC's balance of plant.

TEGs were mounted on steel sheets below the cooling exchanger and above the water exchanger. The temperature difference between the hot and cold sides of TEG was the driving force for the generation of electricity. Fins were implemented in the water compartment for better heat transfer between the water and the steel layer. The locations and the geometry of these fins were adapted from Fernández-Yañez [28]. Figure 6 shows the schematic of the TEG unit.

**Figure 6.** Schematic of the TEG unit for waste heat recovery: (**a**) assembled unit; (**b**) water channels on the hot side of the TEG; (**c**) arrangement of the 40 TEG modules in the assembly.

### *2.1. Membrane*

The modeling of water and proton transport follows standard governing equations (see Equation (1)):

$$\mathbf{N}\_{\mathbf{k}} = \mathbf{J}\_{\mathbf{k}} + \mathbf{c}\_{\mathbf{k}} \mathbf{u}^{\text{mix}} \tag{1}$$

where umix (m.s<sup>−</sup>1) is the mixture velocity in the z-direction in Figure 4, k refers to either water or protons, Nk (mol.m<sup>−</sup>2.s<sup>−</sup>1) is the molar flux due to electro-osmotic driving forces and convection, Jk (mol.m−.s<sup>−</sup>1) is the molar diffusion flux, and ck (mol.m−3) is the molar concentration. The molar flux of water accounts for back diffusion and electroosmotic drag:

$$\mathrm{J}\_{\mathrm{H\_2O}} = -\mathrm{D\_{H\_2O}} \frac{\partial \mathrm{c\_{H\_2O}}}{\partial \mathrm{x}} + \mathrm{n\_{drag}} \frac{\mathrm{i\_{proticic}}}{\mathrm{F}} \tag{2}$$

where iprotonic (A.m<sup>−</sup>2) is the protonic current density in the z-direction (Figure 4), ndrag is the dimensionless drag coefficient, cH2O (mol.m<sup>−</sup>3) is the molar concentration of water, JH2O (mol.m<sup>−</sup>2.s<sup>−</sup>1) is the molar diffusion flux of water, <sup>F</sup> (C.mol<sup>−</sup>1) is Faraday's constant, and DH2O (m2.s<sup>−</sup>1) is the diffusion coefficient. The electrolyte's resistance to transport is dependent upon the dimensionless water content in the membrane λm, and a semiempirical relationship is used for the effective diffusion coefficient [31]:

$$\mathbf{n}\_{\text{drag}} = 2.5 \frac{\lambda\_{\text{m}}}{22} \tag{3}$$

$$
\lambda\_{\rm m} = \frac{\text{c}\_{\rm H\_2O}}{\frac{\rho\_{\rm dry}^m}{\text{M}^m} - \text{b} \text{c}\_{\rm H\_2O}} \tag{4}
$$

$$\mathbf{D}\_{\rm H\_2O} = \mathbf{D}' \left[ \exp \left( 2416 \cdot \left( \frac{1}{303} - \frac{1}{\Gamma} \right) \right) \right] \cdot \frac{1}{\mathbf{a} \left( 17.81 - 78.9 \mathbf{a} + 108 \mathbf{a}^2 \right)} \cdot \lambda\_{\rm m} \tag{5}$$

where b is the dimensionless membrane extension coefficient in the z-direction (based on experiments, b = 0.0126 [31]), a is a dimensionless parameter called water activity, T (K) is temperature, M<sup>m</sup> (kg.mol<sup>−</sup>1) is the membrane molecular mass, ρ<sup>m</sup> dry (kg.m<sup>−</sup>3) is the dry membrane density, and D (m2.s<sup>−</sup>1) is the diffusion coefficient at the reference temperature. The total molar flux of water is consequently:

$$\text{N}\_{\text{H}\_2\text{O}} = \text{J}\_{\text{H}\_2\text{O}} + \text{c}\_{\text{H}\_2\text{O}}\text{u}^{\text{mix}} \tag{6}$$

Proton transport is computed according to Equations (7) and (8):

$$\mathbf{J}\_{\rm H^{+}} = -\frac{\mathbf{F}}{\mathbf{R\_{u}T}} \mathbf{D\_{H^{+}}c\_{\rm H^{+}}} \frac{\partial \phi\_{\rm m}}{\partial \mathbf{x}} \tag{7}$$

$$\mathbf{N}\_{\rm H^{+}} = \mathbf{J}\_{\rm H^{+}} + \mathbf{c}\_{\rm H^{+}} \mathbf{u}^{\rm mix} \tag{8}$$

where φ<sup>m</sup> (V) is the electrostatic potential within the membrane, provided by Equation (9) with membrane conductivity Equation (10), cH<sup>+</sup> (mol.m<sup>−</sup>3) is the molar concentration of protons, J <sup>H</sup><sup>+</sup> (mol.m<sup>−</sup>2.s<sup>−</sup>1) is the molar diffusion flux of protons, Ru (J.mol<sup>−</sup>1.K<sup>−</sup>1) is the universal gas constant, and DH<sup>+</sup> is the proton diffusivity, assumed constant (4.5 <sup>×</sup> <sup>10</sup>−<sup>9</sup> m2.s<sup>−</sup>1)

$$\frac{\partial \phi\_{\rm m}}{\partial \mathbf{x}} = -\frac{\mathbf{i}\_{\rm protmic}}{\sigma\_{\rm m}} + \frac{\mathbf{F}}{\sigma\_{\rm m}} \mathbf{c}\_{\rm H^{+}} \mathbf{u}^{\rm mix} \tag{9}$$

$$
\sigma\_{\rm m} = (0.005139 \lambda\_{\rm m} - 0.00326) \exp\left[1268 \left(\frac{1}{T\_{\rm ref}} - \frac{1}{T}\right)\right] \tag{10}
$$

Here, σ<sup>m</sup> is the dimensionless conductivity of the membrane.

The applied solution procedure and boundary conditions to solve the system of equations in the membrane follows the detailed description in Spiegel [31].

### *2.2. Catalyst Layer*

The agglomerate CL model assumes the uniform distribution of the electrochemical reaction rate and aims at capturing the effect of the main microstructural feature of the CL, that is, an agglomerate type structure where the electrocatalyst is supported by a carbon agglomerate (see Figure 7).

**Figure 7.** A schematic of the catalyst layer, membrane, and microporous layer indicating the agglomerates inside the catalyst layer.

The volumetric current density ic, (A.m<sup>−</sup>3) is calculated using Equation (11), considering uniform distribution of identical agglomerates in the cathode CL [32]:

$$\mathbf{i}\_{\rm c} = 4 \mathbf{F} \frac{\mathbf{P}\_{\rm O\_2}}{\mathbf{H}\_{\rm O\_2, \rm agg}} \left( \frac{1}{\mathbf{Ek}^{\prime} (1 - \varepsilon\_{\rm CL})} + \frac{\mathbf{r\_{agg}} + \delta\_{\rm i} + \delta\_{\rm w}}{\mathbf{r\_{agg}}} \left( \frac{\delta\_{\rm i}}{\mathbf{a\_{agg}} \mathbf{i\_{D\_2,i}}} + \frac{\delta\_{\rm w}}{\mathbf{a\_{agg}} \mathbf{D\_{O\_2,w}}} \right) \right)^{-1} \tag{11}$$

where, PO2 (Pa) is the partial pressure of oxygen, HO2,agg (Pa.m3.mol<sup>−</sup>1) is Henry's constant for oxygen in the agglomerate, k (1.s<sup>−</sup>1) is the reaction rate parameter, aagg,w (1.m<sup>−</sup>1) is the liquid water effective agglomerate surface area, and aagg, i (1.m<sup>−</sup>1) is the ionomer effective agglomerate surface area. DO2,i (m2.s<sup>−</sup>1) and DO2,w (m2.s<sup>−</sup>1) are the oxygen diffusivity in the ionomer and liquid water, respectively, while δ<sup>i</sup> (m) and δ<sup>w</sup> (m) are the corresponding values for the thickness of ionomer film and the thickness of the liquid water film. ragg (m) is the agglomerate radius, while εCL is the volume fraction of void in the catalyst layer:

$$
\varepsilon\_{\rm CL} = \mathbf{1} - \mathbf{L}\_{\rm \frac{\rm D}{\rm C}} - \mathbf{L}\_{\rm i} \tag{12}
$$

Here, Li and LPt/C are the corresponding values of volume fractions of the ionomer phase and the Pt/C particles, as follows [32]:

$$\mathcal{L}\_{\text{Pt}} = \frac{\mathbf{m}\_{\text{Pt}}}{\mathbf{t}\_{\text{CL}}} \left( \frac{1}{\rho\_{\text{Pt}}} + \frac{1 - \mathbf{f}}{\mathbf{f}} \frac{1}{\rho\_{\text{c}}} \right) \tag{13}$$

$$\mathbf{L}\_{\rm i} = \frac{\mathbf{L}\_{\rm g}}{\mathbf{r}\_{\rm agg}^{3} \left(1 - \mathbf{L}\_{\rm i, agg}\right)} \left[ \mathbf{r}\_{\rm agg}^{3} \mathbf{L}\_{\rm i, agg} + \left( \left( \mathbf{r}\_{\rm agg + \delta\_{i}} \right)^{3} - \mathbf{r}\_{\rm agg}^{3} \right) \right] \tag{14}$$

where tCL (m) is the thickness of the CL, ρpt (kg.m<sup>−</sup>3) is the density of platinum, and ρ<sup>c</sup> (kg.m<sup>−</sup>3) is that of carbon. mpt (mg.cm<sup>−</sup>2) is the platinum loading, while Li,agg is the volume fraction of the ionomer phase inside the agglomerate. Equation (15) also indicates the platinum loading divided by the sum of platinum loading and carbon loading (f):

$$\mathbf{f} = \frac{\mathbf{m\_{pt}}}{\mathbf{m\_{pt}} + \mathbf{m\_c}} \tag{15}$$

Afterwards, the thickness of the ionomer film (δi) and the thickness of the liquid water film (δw) can be calculated [32]:

$$\delta\_{\rm i} = \text{r}\_{\rm agg} \left[ \sqrt[3]{\frac{\text{L}\_{\rm i} \left(1 - \text{L}\_{\rm i,agg} \right)}{\text{L}\_{\rm f}} - \text{L}\_{\rm i,agg} + 1} - 1 \right] \tag{16}$$

$$\delta\_{\rm w} = \sqrt[3]{\left(\mathbf{r\_{agg}} + \delta\_{\rm i}\right)^{3} + \frac{3\mathbf{s}\boldsymbol{\varepsilon}}{4\pi\mathbf{N}\_{\rm v}}} - \left(\mathbf{r\_{agg}} + \delta\_{\rm i}\right) \tag{17}$$

Here, s is the liquid water saturation, while Nv (1.m<sup>−</sup>3) is the number of agglomerate particles per CL volume [32]:

$$\text{N}\_{\text{v}} = \frac{\text{3L}\_{\text{ft}}}{4\pi \text{r}\_{\text{agg}}^3 \left(1 - \text{L}\_{\text{i,avg}}\right)}\tag{18}$$

Equivalent governing equations are used for the cathode and anode, but only the former case is detailed hereafter. Assuming a first-order Tafel kinetic reaction yields, by solving the mass conservation equation in a spherical agglomerate, a relation for the

effectiveness factor (Equation (19)) [33]. The effectiveness factor is the ratio of the average reaction rate in the agglomerate core to the reaction at the surface of the core.

$$\mathcal{E} = \frac{1}{3\varrho^2} (3\varrho \coth(3\varrho) - 1) \tag{19}$$

Herein, ϕ is the Thiele modulus:

$$\mathbf{q}^{\prime} = \frac{\mathbf{r\_{\text{agg}}}}{3} \sqrt{\frac{\mathbf{k}^{\prime}}{\mathbf{L}\_{\text{i,agg}}^{1.5} \mathbf{D}\_{\text{O}\_2^{\text{eff}}}^{\text{eff}}}} \tag{20}$$

Deff O2 (m2.s<sup>−</sup>1) is the effective oxygen diffusion coefficient in the catalyst layer [31].

$$\mathbf{k}' = \frac{\mathbf{\dot{i}\_{0, \text{ORR}}} \mathbf{\dot{a}\_{\text{agg}}}}{4 \mathbf{F} (1 - \varepsilon\_{\text{CL}}) \mathbf{C}\_{\text{O}\_2}^{\text{ref}}} \exp \left( - \frac{\alpha\_{\text{c}} \mathbf{F}}{\mathbf{R}\_{\text{u}} \mathbf{T}} (\mathfrak{m}\_{\text{ORR}}) \right) \tag{21}$$

Here, i0,ORR (A.m<sup>−</sup>2) is the reference exchange current density at cathode, cref O2 (mol.m<sup>−</sup>3) is the reference molar concentration of oxygen, ηORR (V) is the cathode overpotential, α<sup>c</sup> is the symmetry cathodic charge transfer coefficient, and aagg (1.m<sup>−</sup>1) is the effective agglomerate surface area, as follows:

$$\mathbf{a}\_{\text{agg}} = \frac{\text{3L}\_{\text{\textdegree F}} \varepsilon\_{\text{CL}}}{\text{r}\_{\text{agg}}^3 \left(1 - \text{L}\_{\text{i,agg}}\right)} \left(\mathbf{r}\_{\text{agg}} + \boldsymbol{\delta}\_{\text{i}}\right)^2 \tag{22}$$

The surface concentration can be assumed to be equal to the bulk concentration based on Fick's law, considering the low solubility of oxygen, steady-state conditions, and thinness of the film. The details of the derivation for this equation are available in refs. [31,33] for the two limiting cases ϕ >> 1 and ϕ << 1. The former indicates slow O2 diffusion in the agglomerate compared to the ORR rate yielding E ∼= <sup>1</sup> <sup>ϕ</sup> , while the latter opposite situation results in E ∼= 1.

Equation (11) holds for the cathode if the external mass transfer limitations can be neglected. These mass transfer limitations occur when the diffusion rate of the molecules is lower than the reaction rate.

The boundary conditions at the inlet of the gas flow channel are the temperature (353.15 K), mass flow rate, and the species mass fraction. At a current density of 1.0 A.cm−<sup>2</sup> (here referred to as the reference current density Iref), Equations (23) and (24) provide the mass flow (kg.s<sup>−</sup>1) of reactants at anode and cathode, respectively:

$$\mathbf{Q\_{a}} = \frac{\zeta\_{\rm a} \mathbf{M\_{H\_{2}}}}{2 \mathbf{F} \mathbf{Y\_{H\_{2}}}} \mathbf{I\_{ref}} \mathbf{A\_{m}} \tag{23}$$

$$\mathbf{Q}\_{\rm c} = \frac{\zeta\_{\rm c} \mathbf{M}\_{\rm O\_2}}{2 \mathbf{F} \mathbf{Y}\_{\rm O\_2}} \mathbf{I}\_{\rm ref} \mathbf{A}\_{\rm m} \tag{24}$$

where ζ<sup>a</sup> = 1.5 is the stoichiometric ratio at the anode, while that of cathode is ζ<sup>c</sup> = 2.0. MH2 (kg/mol) and MO2 (kg.mol<sup>−</sup>1) are the molecular weights of hydrogen and oxygen, respectively, while Yk indicates the mass fraction of component k, and Am (m2) is the representative of the membrane's area. The inlet liquid water saturation is assumed to be zero, like the electrostatic potential at the anode terminal. The electrostatic potential at the cathode terminal is also equal to the cell's voltage (Vcell) at the operating temperature of 353.15 K. The overpotential terms to calculate the potential of the cell Vcell under polarization, which is based upon the Nernst equation (see Equation (25)), are computed from the above governing equations, as discussed in detail by Siegel et al. [34]:

$$\mathbf{V\_{cell}} = \mathbf{V\_0} - \mathbf{V\_{act}} - \mathbf{V\_{chmic}} - \mathbf{V\_{corr}} \tag{25}$$

Here, V0 (V) is the Nernst potential, while Vact (V), Vohmic (V), and Vconc (V) are the activation, ohmic, and concentration overpotentials, respectively.

The generated present study considers an upper bound limiting case, that is, all the output thermal power of the PEMFC system (Equation (26)) is assumed transferred to the coolant loop as an increase in fluid enthalpy, since the flux is constant.

$$P\_{\text{thermal}} = -(\mathbf{V\_{act}} + \mathbf{V\_{ohm}} + \mathbf{V\_{conc}}) . \mathbf{N\_{cell}} \,\mathrm{i} \tag{26}$$

#### *2.3. Thermoelectric Generator (TEG)*

ANSYS CFX module 19.2 has been used to solve the three-dimensional continuity, momentum, and energy equations for viscous, laminar, and incompressible flow of water in the hot- and cold-side heat exchangers, as well as heat transport in the solid domains (Figure 6). The CFD model, therefore, corresponds to a water–water heat exchanger. Figure 8 shows the details of the simulation procedure to analyze the performance of the TEG unit.

**Figure 8.** The flowchart of simulation procedure to perform the CFD analysis of the TEG unit.

Empirical relationships were derived to calculate during post-processing the TEG electrical performance under the thermal conditions simulated by the CFD exchanger model. In this regard, the possible additional resistance of the TEG unit, which may result in more energy consumption [35–37], was considered. The CFD simulations and TEG calculations were performed in an uncoupled manner. The local relationship between temperature gradient, potential, and current was computed using the freely available module TEG1-4199-5.3 [38], which shows close to linear current–voltage characteristics. For instance, the tabulated output voltage for Bi2Te3 was 6.7 V at a current, power, hot-side temperature, and cold-side temperature of 1.12 A, 7.5 W, 300 ◦C, and 30 ◦C, respectively. The corresponding open-circuit voltage of the module was 13.4 V under a heat flow flux of 9.5 W.cm<sup>−</sup>2. It should be noted that the specific heat capacity of this TEG module was 156.1 (J.kg<sup>−</sup>1.K<sup>−</sup>1), while the other materials utilized in the TEG unit were aluminum and steel with specific heat capacities of 434 (J.kg<sup>−</sup>1.K<sup>−</sup>1) and 903 (J.kg<sup>−</sup>1.K<sup>−</sup>1), respectively. As the hot-side temperature was different in this study, a correlation was extrapolated based on the data in ref. [38]. Equation (27) is the fitted empirical expression used to calculate the TEG open voltage VTEG (V) in the system, assuming a constant cold-side solid temperature of 30 ◦C.

$$\mathbf{V\_{TEG}} = -0.00004 \times \left(\mathbf{T\_h} - 273.15\right)^2 + 0.0638 \times \left(\mathbf{T\_h} - 273.15\right) - 1.9012 \tag{27}$$

Here, Th (K) is the hot-side solid temperature. Similarly, the TEG internal resistance RTEG (Ω) is approximated by the following empirical relationship:

$$R\_{\rm TEG} = -0.00001087 \times \left(\text{T}\_{\rm h} - 273.15\right)^2 + 0.008726 \times \left(\text{T}\_{\rm h} - 273.15\right) + 3.922 \tag{28}$$

The TEG output electrical power PTEG (W) is then:

$$P\_{\rm TEG} = \frac{\left(-0.00004 \times \left(\text{T}\_{\text{h}} - 273.15\right)^2 + 0.0638 \times \left(\text{T}\_{\text{h}} - 273.15\right) - 1.9012\right)^2}{-0.00001087 \times \left(\text{T}\_{\text{h}} - 273.15\right)^2 + 0.008726 \times \left(\text{T}\_{\text{h}} - 273.15\right) + 3.922} \tag{29}$$

For rapid screening analyses, the efficiency of a TEG module (ηTEG) can be expressed as [39]:

$$\eta\_{\rm TEG} = \frac{\mathbf{T\_h} - \mathbf{T\_c}}{\mathbf{T\_h}} \frac{\sqrt{1 + Z\mathbf{T\_m}} - 1}{\sqrt{1 + Z\mathbf{T\_m}} + \frac{\mathbf{T\_c}}{\mathbf{T\_h}}} \tag{30}$$

The first right-hand side term corresponds to the Carnot efficiency, and Tc is the cold source temperature. ZTm is a dimensionless parameter that characterizes the TEG material performance at an average temperature Tm (K):

$$\mathbf{T\_m} = \frac{\mathbf{T\_h} + \mathbf{T\_c}}{2} \tag{31}$$

$$Z = \frac{\sigma \alpha^2}{\kappa} \tag{32}$$

where α (V.K<sup>−</sup>1) is the Seebeck coefficient, σ (1.m<sup>−</sup>1.Ω<sup>−</sup>1) is the electrical conductivity, and κ (W.m<sup>−</sup>1.K<sup>−</sup>1) is the thermal conductivity of the TEG material.
