Thermoelectric Generator as the Waste Heat Recovery Unit of Proton Exchange Membrane Fuel Cell: A Numerical Study

: The proton exchange membrane fuel cell (PEMFC) is a prominent environmentally friendly alternative candidate to internal combustion engines in automotive applications. The recovery of the waste heat of light-duty diesel engines has been investigated recently, which is similarly relevant for PEMFCs. Thermoelectric generators (TEG) applied on the stack’s walls have been already proposed and tested as a cooling method for small scale applications of the PEMFC. For the medium scale usages of the PEMFC stack, TEG technology may be further used to recover heat lost through the cooling water required for stack thermal management, which was the focus of the present study. Using an agglomerate model for the PEMFC and a computational fluid dynamic (CFD) thermal model for the TEG heat exchanger unit, the operation and performance of the PEMFC stack and heat recovery unit were simulated, respectively. After validation, results indicated that the transferred heat from the PEMFC to the cooling channel increased the temperature of the coolant from room temperature to 330.5 K at the current density of 0.8 A/cm 2 . CFD analysis revealed that 37.7 W of the heated wasted by the PEMFC stack could be recovered by the currently available TEG material and geometry.


Introduction
Concentrated efforts have been made to reduce the usage of fossil fuels after ascertaining their adverse contributions to climate change and pollution [1]. The automotive industry still largely relies on internal combustion engines fueled by diesel [2]. However, fuel cells and batteries are becoming viable alternatives to diminish and also control the localization of pollutant emissions [3]. Herein, proton exchange membrane fuel cells (PEMFCs) are a prominent technology for the automotive industry that electrochemically converts hydrogen and oxygen into water, electricity, and heat [4].
In different industries, there exist considerable amounts of heat that can be recovered [5]. Figure 1 illustrates the share of waste heat energy in United Kingdom. Approximately half of the hydrogen utilized by the low-temperature PEMFC changes to electricity, while the rest mainly converts to heat, which should be removed to keep the operational temperature range around 65 • C to 85 • C [6]. This temperature range stems from constraints related to the electrolyte membrane, usually made of Nafion [7]. Operation at temperatures above or below leads to drying of the membrane or flooding, respectively [8]. The heat produced by PEMFCs can be removed through natural convection, water evaporation, or extra reactants, but mainly through the cooling system [9]. Thermal and water management, that is, the combined control of the spatial distribution of temperature and humidity in the stack, is required to guarantee high performance and durability for PEMFCs [10]. In this regard, a precise and accurate model is needed to characterize the performance of the PEMFC. Both experimental and numerical models are used to evaluate the output power density and voltage of the cell, hence the corresponding amount of waste heat [12]. The electrochemical performance of the catalyst layer (CL) can be simulated with varying levels of detail and accuracy in modeled phenomena and discretization, from zero to threedimensional [13]. The corresponding physics encompasses electrochemistry coupled with multiphase liquid, gas, and solid phase transport in the reticulated porous and Pt-loaded carbon networks [14]. Among the proposed modeling approaches for the CL [15], the socalled agglomerate model proved to be precise, at the cost of higher complexity compared to the porous electrode, interface, and microscopic models.
The agglomerate model provides an improved representation of the CL by averaged microstructural parameters, in particular the electrocatalyst and carbon agglomerates covered by a thin layer of membrane. The reaction rate is further assumed to be homogeneous. Porous electrode models also consider the presence of porous agglomerates of ionomercovered carbon particles in the CL. They, however, represent the system at the scale of the catalyst layer, while the agglomerate model accounts explicitly for the diffusion of products and reactants throughout the electrolyte film, which surrounds the electrocatalyst and carbon supports [16]. Simpler microscopic and single-pore models consider cylindrical Teflon-coated pores for gas diffusion, whereas interface models treat the CL as interface conditions between the gas diffusion layer (GDL) and membrane. In addition to the selection of the right model for the PEMFC, the most suitable strategy should be selected to recover the wasted heat and improve the thermal management [17][18][19].
Based on the size of the PEMFC stack, different cooling methods, such as air, passive, liquid, and phase-change cooling, can be used. The studies by Kwan [20][21][22] and Shen [23,24] explored thermoelectric generators (TEGs) implemented on the walls of the bipolar plates (BPs) to recover and convert the waste heat directly into electricity in small-scale applications. The performance is currently very modest: in the design developed by Kwan et al. [21], each kg of TEG material is expected to produce 0.15 W of electricity under optimized conditions. TEGs generate electricity when subjected to a temperature gradient, and they are used in different applications such as temperature sensors, wearable energy harvesters, or for waste heat recovery [25]. Based on the trend described by He et al. [26] on the increase in the figure of merit for different materials, Shen et al. [24] evaluated the corresponding efficiency of different TEG units with different materials, i.e., considering different figures of merit, by changes in hot-side temperature. Figure 2 shows a similar trend for efficiency based on the results of Shen et al. [24] and Jouhara et al. [27] in 2016 concerning changes in the hot-side temperature and materials of TEG units in the constant cold-side temperature of 300 K. It should be noted that the data mentioned in Figure 2 were the highest possible values that could be obtained with recent TEG materials in 2016. For example, in a similar temperature range, the efficiency of Bi 2 Te 3 varied between 2% and 9% in 2010, based on the values reported by Shen et al. [24]. Fernández-Yañez et al. [28] and Zorbas et al. [29] considered conventional TEG modules to recover the waste heat from the output exhaust gas of light-duty diesel engines. Bi 2 Te 3 modules of 2.5 cm × 2.5 cm in size were installed at different positions on the exhaust pipe of a 96 kW engine car (1995 cm 3 ). The power and recovery efficiency for the TEG unit in the range of 300 W and 5.5% were achieved, respectively, with a cold-side temperature of 30 • C and at part-load engine conditions. Projections by Zorbas et al. [29] about potential progress in TEG technology suggest that a reduction in fuel consumption by approximately 20% may be reached in the future, compared to 5% with current technology [29]. From an economic standpoint, the cost of the technology considered by Zorbas et al. [29], with an estimated 5% reduction in fuel consumption using TEG technology, was in the range of 500 EUR for a 96 kW engine implementation with amortization expected within 2 to 3 years. Although there have been suggestions for using the thermoelectric generators for waste heat recovery at small scales, TEGs have not been used on a medium-scale level accounting for liquid cooling methods.
The main novelty of the current study was to investigate the performance of TEG modules for the recovery of waste heat from the PEMFC coolant channels using the liquid cooling method with medium-scale stacks of PEMFCs (see Figure 3). The expected recovery performance was assessed assuming standard TEG material performance in terms of temperature levels. An agglomerate PEMFC model and three-dimensional thermal modeling of the TEG-based recovery module were then used to simulate the performance of the current standard material technology and component geometries. The PEMFC model was validated at low and high current densities using the I-V characteristic curve. Furthermore, a 3-D CFD thermal model of the TEG unit was developed, and pressure drops compared with numerical data by Fernández-Yañez et al. [28]. Empirical relationships were used for calculations by post-processing of TEG performance corresponding to the thermal conditions. The waste heat from the PEMFC was computed by a modeling approach with agglomerate formalism for the CL. The recovery performance and corresponding spatial distribution of temperature along the TEG heat exchangers were simulated and discussed.

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. 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. 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.

Membrane
The modeling of water and proton transport follows standard governing equations (see Equation (1)): where u mix (m.s −1 ) is the mixture velocity in the z-direction in Figure 4, k refers to either water or protons, N k (mol.m −2 .s −1 ) is the molar flux due to electro-osmotic driving forces and convection, J k (mol.m − .s −1 ) is the molar diffusion flux, and c k (mol.m −3 ) is the molar concentration. The molar flux of water accounts for back diffusion and electroosmotic drag: where i protonic (A.m −2 ) is the protonic current density in the z-direction (Figure 4), n drag is the dimensionless drag coefficient, 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]: 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 m (kg.mol −1 ) is the membrane molecular mass, ρ m dry (kg.m −3 ) is the dry membrane density, and D (m 2 .s −1 ) is the diffusion coefficient at the reference temperature. The total molar flux of water is consequently: Proton transport is computed according to Equations (7) and (8): where φ m (V) is the electrostatic potential within the membrane, provided by Equation (9) with membrane conductivity Equation (10), is the universal gas constant, and D H + is the proton diffusivity, assumed constant (4.5 × 10 −9 m 2 .s −1 ) Here, σ m 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].

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). The volumetric current density i c , (A.m −3 ) is calculated using Equation (11), considering uniform distribution of identical agglomerates in the cathode CL [32]: where, P O 2 (Pa) is the partial pressure of oxygen, H O 2 ,agg (Pa.m 3 .mol −1 ) is Henry's constant for oxygen in the agglomerate, k (1.s −1 ) is the reaction rate parameter, a agg,w (1.m −1 ) is the liquid water effective agglomerate surface area, and a agg, i (1.m −1 ) is the ionomer effective agglomerate surface area. D O 2 ,i (m 2 .s −1 ) and D O 2 ,w (m 2 .s −1 ) are the oxygen diffusivity in the ionomer and liquid water, respectively, while δ i (m) and δ w (m) are the corresponding values for the thickness of ionomer film and the thickness of the liquid water film. r agg (m) is the agglomerate radius, while ε CL is the volume fraction of void in the catalyst layer: Here, L i and L Pt/C are the corresponding values of volume fractions of the ionomer phase and the Pt/C particles, as follows [32]: where t CL (m) is the thickness of the CL, ρ pt (kg.m −3 ) is the density of platinum, and ρ c (kg.m −3 ) is that of carbon. m pt (mg.cm −2 ) is the platinum loading, while L i,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): Afterwards, the thickness of the ionomer film (δ i ) and the thickness of the liquid water film (δ w ) can be calculated [32]: Here, s is the liquid water saturation, while N v (1.m −3 ) is the number of agglomerate particles per CL volume [32]: 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.
Herein, ϕ is the Thiele modulus: is the effective oxygen diffusion coefficient in the catalyst layer [31].
Here, i 0,ORR (A.m −2 ) is the reference exchange current density at cathode, c ref is the reference molar concentration of oxygen, η ORR (V) is the cathode overpotential, α c is the symmetry cathodic charge transfer coefficient, and a agg (1.m −1 ) is the effective agglomerate surface area, as follows: 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 O 2 diffusion in the agglomerate compared to the ORR rate yielding E ∼ = 1 ϕ , 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 −2 (here referred to as the reference current density I ref ), Equations (23) and (24) provide the mass flow (kg.s −1 ) of reactants at anode and cathode, respectively: where ζ a = 1.5 is the stoichiometric ratio at the anode, while that of cathode is ζ c = 2.0. M H 2 (kg/mol) and M O 2 (kg.mol −1 ) are the molecular weights of hydrogen and oxygen, respectively, while Y k indicates the mass fraction of component k, and A m (m 2 ) 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 (V cell ) at the operating temperature of 353.15 K. The overpotential terms to calculate the potential of the cell V cell 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]: Here, V 0 (V) is the Nernst potential, while V act (V), V ohmic (V), and V conc (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.

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. 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][36][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 Bi 2 Te 3 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 −2 . It should be noted that the specific heat capacity of this TEG module was 156.1 (J.kg −1 .K −1 ), while the other materials utilized in the TEG unit were aluminum and steel with specific heat capacities of 434 (J.kg −1 .K −1 ) and 903 (J.kg −1 .K −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 V TEG (V) in the system, assuming a constant cold-side solid temperature of 30 • C.
Here, T h (K) is the hot-side solid temperature. Similarly, the TEG internal resistance R TEG (Ω) is approximated by the following empirical relationship: The TEG output electrical power P TEG (W) is then: For rapid screening analyses, the efficiency of a TEG module (η TEG ) can be expressed as [39]: The first right-hand side term corresponds to the Carnot efficiency, and T c is the cold source temperature. ZT m is a dimensionless parameter that characterizes the TEG material performance at an average temperature T m (K): where α (V.K −1 ) is the Seebeck coefficient, σ (1.m −1 .Ω −1 ) is the electrical conductivity, and κ (W.m −1 .K −1 ) is the thermal conductivity of the TEG material. Table 1 lists the input parameters to compute the output thermal power and mass flow rates of the current PEMFC model. The latter are then used as input data for TEG module simulations, to investigate the recovery performance along with the flow behavior and spatial temperature distribution in the TEG unit.

Model Validation
The aim of the current developed PEM fuel cell model was to predict the output temperature and the mass flow rate of the water in the cooling channels, which enters the TEG heat exchanger. In this regard, to validate the output results of the PEMFC model, the I-V characteristic curves of the model should be compared with the verified results in the literature at both low and high current densities. Once the validation is performed, the model is suitable to predict the needed output values. The PEMFC model predictions were tested with both experimental and simulation results by Ticianelli et al. [42] and Mohammedi et al. [43], respectively. It is noteworthy to mention that the input model parameter values for the present verification were not those listed in Table 1, but correspond to the values reported in Mohammedi et al. [43] to allow for direct comparison.   [42] and simulations of Mohammedi et al. [43]; (b) model comparison with the 500 W BCS stack using the experimental data of Xue et al. [44], Correa et al. [45], and simulation data of Sharifi-Asl et al. [46]. Furthermore, the current model was also validated with three different sets of experimental and numerical data at low current densities. Figure 9b compares the output of the current PEM fuel cell model with the 500 W BCS stack manufactured by BCS technologies [44]. The number of cells for the 500 W BCS stack [44] was 32, with an active area of 64 cm 2 operating at 333 K. It should be noted that in this stack, the resistance to the electron flow was 0.0003 Ω, while the limiting current density was 0.469 A.cm −2 . Accounting for the high current densities, Figure 10 illustrates the I-V characteristic curve of the current model with the output performance of the single-cell Ballard Mark V PEM fuel cell [44] at the operating temperature of 343 K. In this case, the limiting current density was 1.5 A.cm −2 , the active area of the cell was 50.6 cm 2 , and the resistance to the electron flow was 0.0003 Ω. The heat transfer verification of the TEG model consisted of the comparison of pressure drops computed by Fernández-Yañez et al. [28] (see Table 2), which was a follow-up of the investigation in ref. [47]. The published dataset included nine stationary modes of a light-duty diesel engine at different torques. In that article, simulations were performed to obtain the corresponding temperature and pressure drop of the exhaust gases of the diesel engine. The present TEG model verification used a subset of four of the modes reported in [26] (A, I, G and D), which means four different input temperatures and mass flow rates. The agreement shown in Table 2 indicates that the present model pressure drop predictions were in line with the mixed numerical and experimental study by Fernández-Yañez et al. [28]. In every CFD analysis, the values of the output results should be independent from the size of the cells in the mesh structure. In this regard, a grid independency analysis is needed to prove the suitability of the simulation model for the TEG unit. As the validation of this unit (see Table 2) was conducted with the given data by Fernández-Yañez et al. [28], the grid independency of the TEG thermal model was also performed using the temperature corresponding to engine mode D. Figure 11 shows that grid independency of the temperature simulated at 475 K was reached starting around 4.2 × 10 6 cells. In this condition, the average surface area of the cells was 1.57 × 10 −3 m 2 , the minimum edge length was 1 × 10 −3 m, and the maximum cell edge size was 0.1132 m, with highest skewness of 0.9. Figure 12 shows the mesh structure of the developed model, only when the number of cells is lowest (1 × 10 6 ), for illustration reasons. Figure 11. Grid independency study of the three-dimensional CFD model of the TEG unit considering the temperature in the engine mode D given by Table 2 as the verification value.

PEMFC Modeling
In this study, the mass flow rates and temperature of the exhaust water from the PEMFC's coolant channel were used as inputs for the TEG unit simulations. Figure 13a,b shows the computed characteristic curves of the PEM fuel cell with the input parameters listed in Table 1. The open circuit voltage at low current density decreased because of fuel crossover. Activation losses caused the pronounced non-linear decrease in voltage at low current density. Ohmic losses resulted in a close to linear dependence until the onset of mass transport limitations became dominant at high current densities and fuel conversion.  Table 1 considering fuel utilization of 78% and 800 mA/cm 2 ; (a) I-V, (b) P-I, (c) simulated dependence of the water content and effective conductivity on the thickness of the membrane.
The effects of selected model parameters relevant to thermal and water management were further evaluated using published data, as a further verification after the I-V curve comparison in Section 3.1. The one-dimensional numerical study by Falcão et al. [48], which was compared to experimental data [49,50] and a three-dimensional model [51], evaluated differences in water content and local transport properties, depending upon the membrane thickness using a one-dimensional model. Their evaluation considering four thicknesses (0.0051 cm, 0.0127 cm, 0.0178 cm, 0.003 cm) concluded that low thicknesses result in lower ionic resistance and voltage loss, hence better cell performance. However, the fabrication of MEA with thin membranes is a challenge, and other factors such as mechanical integrity of the membrane and manufacturing cost should be considered. Figure 13c illustrates the variations in water content and effective conductivity computed for different membrane thicknesses, showing qualitative agreement.
The purpose of the present agglomerate PEMFC model was to calculate the output thermal power. It should be noted that all types of heat transport and transfer by conduction, radiation, and convection were at play in the electrodes. The sources of heat in the polarized CLs were ohmic losses induced by charge transport, entropy changes, and activation overpotentials due to electrochemical reactions. The effectiveness factor (Equation (19)) is a metric for the assessment of the cell performance. An effectiveness factor equal to one corresponds to the limiting case where all the catalyst surface area is active under polarization. It can be determined for both cathode and anode sides after calculation of the Thiele modulus (Equation (20)), which is the ratio of the reaction rate on the catalyst surface to the reaction rate without any transport losses. Figure 14 illustrates the changes in the effectiveness factor and the flux density of hydrogen upon variations in current density. Figure 14b shows that transport losses increase with increasing current, leading to an effectiveness factor approaching zero. Figure 14a indicates that hydrogen crossover increases drastically at current densities higher than 0.8 A/cm 2 , which is detrimental for the cell's performance and is known to increase the degradation. A current density of 0.8 A/cm 2 was, therefore, the selected operation point for the coupling with the TEG module. The corresponding simulated electrical power and efficiency were 3.34 kW and 39%.

TEG Simulation
As mentioned, the TEG unit will not have a direct impact on the performance of the PEMFC unit, although the integrated system of the PEMFC-TEG will have higher output power in comparison to the standalone PEMFC. The simulated output flow of water from the PEMFC's coolant channel (3 kg/s at 330.5 K, at 0.8 A/cm 2 ) is used as the hot-side water exchanger inlet condition for the analysis of the TEG unit for waste heat recovery. Prescribed boundary conditions are the pressure outlet of the TEG unit and a temperature of 30 • C at all exterior walls of the PEMFC and TEG coolant domains, except those in contact with the TEG modules. A uniform convection heat transfer coefficient of 15 W m 2 K is assumed. The cold side coolant exchanger (see Figure 6) is fed with 1.37 kg/s of water at 30 • C. The reason for the mentioned flow rates of water in the TEG heat exchange module is to validate the results with the given data by Fernández-Yañez et al. [28]. Ahn and Choe [52] had already suggested the minimum mass flow rate of 0.8 kg/s for the PEM fuel cell in the coolant channel. The inlet temperature at the hot side of the water exchanger has been calculated using the agglomerate model for the PEMFC while the inlet temperature at the cold side of the water exchanger is based on the study developed by Liso et al. [53]. The high resolution Rhie-Chow [54,55] scheme with a convergence criterion of 10 −5 is applied for the coupling between pressure and velocity. Additionally, the 40 elements forming the TEG unit are assumed to be thermally and electrically insulated from each other.
Recovering the output flow of the PEMFC by the TEG unit, Figure 15 shows the corresponding contours of temperature, voltage, and power density in the TEG unit both on the water and coolant exchanger sides. The cooling effect of cold water on the flow of water from the PEMFC entering the water exchanger of the TEG unit is as expected in Figure 15. The regions with higher temperature in Figure 15b also result from back flow pressure on the edges, leading the water to flow around a circle. In the water exchanger, the presence of a boundary layer on the sides of the heat exchanger causes lower temperature compared to the other parts of the heat exchanger. Consequently, the temperature on the side of the coolant exchanger is comparatively lower as well. The voltage of the TEG units computed by Equation (27) is dependent upon the solid hot water exchanger side temperature, leading as expected to a potential distribution such as shown in Figure 15c. The maximum voltage is 1.36 V, which is much lower than values around 8 V in ref. [28]. The reason is the higher hot source temperature in the light-duty diesel engines due to higher exhaust heat. The cold-side temperature varies by approximately 2 • C, showing that the assumption of constant cold-side temperature in Equations (27)- (29) was acceptable for the present study, even though it does not allow the analysis of the effects of flow patterns in Figure 15b on the output electrical power. Figure 15d shows the contour of produced electrical power in footprint mesh cell of the TEG unit. It indicates that the maximum local power in the module is 3.61 ×10 −3 W/cm 2 , compared to the total power of 37.7 W and the averaged power density of 0.00806 W/cm 2 .
The corresponding value of the TEG surface-specific power density for the whole unit was also 3.41 W/cm 2 .
Other parameters to be considered were the PEMFC's current density, which was set at 0.8 A/cm 2 in this study, and the PEMFC water cooling flux. Under the corresponding thermal conditions, the visual inspection of the temperature profile indicated that the TEG modules produce up to roughly twice less electrical power than the central ones. Moreover, the temperature difference between the cold and hot side is almost constant along the central lines in Figure 15a. This suggests that higher efficiency recovery may be achieved with a larger, yet more expensive heat exchanger module and/or improved heat exchanger designs. PEMFC operation at higher current density produces more excess heat. In the case where the stack temperature is kept at 353 K (80 • C) and the coolant water flow is kept constant, the temperature of the latter will increase, compared to 330 K (57 • C) in Figure 15a at the inlet. Therefore, the TEG output electrical power is expected to increase. Herein, the beneficial effect on TEG performance may well compensate for the heat exchanger limitations and lead to higher recovery efficiencies. Alternatively, the recovery efficiency in the present case of 0.8 A/cm 2 can be increased by decreasing the PEMFC cooling water flow, which was not the case in the current study, as the coolant flow of PEMFC was in a separate water loop of the PEMFC's stack, and the water mass flow rate was assumed to be constant. As expected, the TEG module and PEMFC must be designed and operated in synergy for viable recovery efficiencies over a sufficiently large window of operation conditions for automotive applications.

Conclusions
This study investigated the relevance of TEG technology to recover the waste heat from the PEMFC. It was inspired by the idea of producing electricity from the exhaust heat of light-duty diesel engines by TEG. Herein, an agglomerate one-dimensional PEMFC model was developed to estimate the maximum output heat from the fuel cell stack. The validation of the I-V characteristic curves of the current model with six different series of experimental and simulation data proved the suitability of the developed model. Under the considered PEMFC operating conditions and assumptions, the simulated output water temperature and flow rate from the PEMFC's coolant channel were 330.5 K and 3 kg/s, respectively. The CFD thermal simulations with cooling water at 30 • C combined with post-processing based on empirical relationships for the local TEG electrical performance resulted in a modest total output electrical power of 37.7 W. The voltage was in the range of 1.25 V and showed non-uniformity (0.77-1.36 V), suggesting potential for improvement of the TEG heat exchanger design, as well as of the PEMFC operating conditions for waste heat recovery. The figures were, therefore, modest. Although this study presented the integration of the TEG unit as a cooling method for a medium-scale size of the PEMFC, further research can be performed for future studies:

•
This study only evaluated the possibility of using a TEG unit as a cooling method in a medium-scale PEMFC stack by calculating the output power and the voltage of the TEG unit. However, exergy and economic aspects of the system could be assessed through exergo-economic and thermo-economic analyses to evaluate the suitability of this integration.

•
The design of the current TEG unit was inspired from the geometry given in ref. [28], which was mainly for diesel engine applications. Due to the novelty of this topic, there is a lack of optimized geometry for the application of PEMFC. In this regard, further studies on the geometry of the TEG unit, the needed number of TEG modules in the TEG unit, and the required mechanical/physical properties of the unit for maximum performance are needed.

•
The current modeling was performed in the steady-state condition for both the PEMFC numerical model and the TEG simulation model. However, performing a dynamic and real-time analysis would be of value to better investigate the possibility of the current integration.

•
The materials in the TEG modules and their figures of merit play a crucial role in the performance of the TEG unit and the waste heat recovery. This study assumed Bi 2 Te 3 as the base material of the TEG modules, while many other alternatives can be analyzed and evaluated to reach the highest waste heat recovery.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.