Cooling Design for PEM Fuel-Cell Stacks Employing Air and Metal Foam: Simulation and Experiment

: A new study investigating the cooling efﬁcacy of air ﬂow inside open-cell metal foam embedded in aluminum models of fuel-cell stacks is described. A model based on a commercial stack was simulated and tested experimentally. This stack has three proton exchange membrane (PEM) fuel cells, each having an active area of 100 cm 2 , with a total output power of 500 W. The state-of-the-art cooling of this stack employs water in serpentine ﬂow channels. The new design of the current investigation replaces these channels with metal foam and replaces the actual fuel cells with aluminum plates. The constant heat ﬂux on these plates is equivalent to the maximum heat dissipation of the stack. Forced air is employed as the coolant. The aluminum foam used had an open-pore size of 0.65 mm and an after-compression porosity of 60%. Local temperatures in the stack and pumping power were calculated for various air-ﬂow velocities in the range of 0.2–1.5 m/s by numerical simulation and were determined by experiments. This range of air speed corresponds to the Reynolds number based on the hydraulic diameter in the range of 87.6–700.4. Internal and external cells of the stack were investigated. In the simulations, and the thermal energy equations were solved invoking the local thermal non-equilibrium model—a more realistic treatment for airﬂow in a metal foam. Good agreement between the simulation and experiment was obtained for the local temperatures. As for the pumping power predicted by simulation and obtained experimentally, there was an average difference of about 18.3%. This difference has been attributed to the poor correlation used by the CFD package (ANSYS) for pressure drop in a metal foam. This study points to the viability of employing metal foam for cooling of fuel-cell systems.


Introduction
A fuel cell is a direct-energy-conversion system: It transforms the chemical energy of hydrogen to electricity. Theoretically, the thermodynamic efficiency of a fuel cell can be 80% or greater. There are no harmful emissions of fuel cells [1]. Hence, fuel cells are seriously investigated as a clean power source for automobiles [2]. Fuel cells are also being considered as a replacement for various kinds of power plants [3]. A good example of the use of fuel cells in buses is given in [4], where it is indicated that these buses had a 62% lower harmful emission compared to buses with diesel engines. Small-size fuel cells may be applied to power portable devices [5].
Proton exchange membrane fuel cells (PEMs) are most popular. PEMs operate at low temperatures, typically 60 to 80 • C, while solid-oxide fuel cells usually operate at about 700 • C [6]. PEMs have a quick startup, and they are less likely to suffer from corrosion [7]. Compared to other types of fuel cells, PEMs seem to be more robust with little or no sealing issues [8]. Not only do PEMs possess high power density, but they also have higher power output compared to direct-methanol fuel cells. The light weight and tininess of PEMs make them the prime contender for replacing the internal combustion engine [8].
Indeed, very few studies investigated metal foam as a potential cooling medium for fuel cells. Boyd and Hooman published a purely numerical study on a small metalfoam cooling system for fuel cells [3]. Odabaee et al. [41] studied metal-foam air-cooling potential for fuel cells. This was intended as a replacement for water cooling. Their experiment examined the enhancement of heat transfer for a thin layer of aluminum foam positioned between adjacent bipolar plates. The study showed that to extract the same amount of heat, water-based systems required twice the power compared to air-cooled metal-foam systems. Afshari et al. [42] simulated a heat-transfer model in a 15 cm × 15 cm square area of cooling channels with different designs: straight flow; serpentine flow with multipasses; and metal foam. They proved that metal foam is the best cooling medium. Santamaria et al. [43] investigated a nickel-foam microchannel for coolingfuel cells. They determined the pressure drop experimentally at different Reynolds numbers for nickel foams having different porosities. They also analyzed heat transfer and fluid flow and used available data to validate their model. A study by Toghyani et al. [44] investigated PEMs' performance using a distributor made from metal foam and compared it to three common flow fields. Current density, temperature, hydrogen mass fraction, and pressure drop were compared for the different models. 3D designs were modeled and analyzed numerically. Results showed that both pressure drop and stack performance are affected by metal-foam permeability. Vazifeshenas et al. [45] studied the effects of using metal foam in the cooling channel of a fuel cell. They designed a 3D model for three flow paths: parallel, serpentine, and serpentine with multichannels. The inclusion of metal foam enhanced heat transfer at the cost of increased pressure drop. They found that copper foam provided larger Nusselt numbers.
Some of the limitations of previous studies include the fact that some of these studies were purely numerical/simulation [3,42], while others were focused on liquid cooling [41,42,45]. Also, some of these studies were intended for a single fuel cell, and they employed a single metal foam with one porosity and one pore size [39]. For air cooling using metal foam, previous studies invoked the local thermal equilibrium assumption [3,43]. This assumption asserts that the temperatures of the air and the metal inside the foam are locally equal, and as such there is no convection heat transfer between them. It has been shown that his assumption is not valid for gas flow in metal foam [46,47]. None of the previous studies investigated an end cell in a stack. Such a cell is subjected to asymmetric heating since the heat comes from only one side of the cell. This study numerically and experimentally investigates an air-cooling system for a commercial PEMs stack employing aluminum foam. In the numerical part, the local thermal equilibrium assumption is relaxed and is replaced by the more physically meaningful assumption of local thermal nonequilibrium. Moreover, both internal and external cells of the stack are investigated.

Numerical Determination of Metal Foam Porosity and Pore Density
A preliminary numerical study was conducted in order to elucidate the effect of porosity and pore density of aluminum foam on its cooling capability. The study was also used to ascertain the effect of air velocity through the foam on heat dissipation. A summary of this study will be given here, however, more details relevant to this study will be elaborated on in Section 2.2 in conjunction with the simulation presented therein. The preliminary study was 2-dimensional ( Figure 1). It simulated a channel between two parallel plates filled with aluminum foam and cooled by air. The two plates comprising the channel represented adjacent fuel cells in a stack. Each plate was subjected to a heat flux of 2.5 W/cm 2 . The dimensions of the metal foam sheet between them were 50 mm × 3.56 mm. Different pore densities of 20, 40, 80, and 100 pores per inch (ppi) with different porosities of 60, 70, 80, and 90% of the foam were investigated. Average inlet air speeds for each case ranged from 0.2 to 3 m/s. One-half of this channel was simulated taking advantage of the symmetry. The solution domain and boundary conditions are shown in Figure 1.  The following governing equations were solved using ANSYS-Fluent:

Continuity :
∂u ∂x  The following governing equations were solved using AN Continuity: where ɛ and k are the porosity and permeability of the foam, r Energy equation for fluid: Energy equation for solid: where and are the velocity components in the flow direc to the flow direction (y), and is the pressure. and µ are t the fluid, respectively. is the form drag coefficient for the f fective thermal conductivities of the fluid and the solid, respec rosity and permeability of the foam, respectively. ℎ is the in ficient, and is the interfacial surface area density (m 2 /m 3 ) ins the fluid and solid temperatures, respectively. The following boundary conditions were imposed: At the channel's inlet: = , = 0, an At the heated walls: ,, = − = − For symmetry: At the channel outlet: = 0 and = Results showed that the higher the ppi, the higher was t of the foam. This is because the surface-area density of the foam ppi. Even though 100 ppi can be achieved, commercial alumin ble at no more than 40 ppi. Therefore, 40-ppi aluminum foam investigated in terms of porosity. Figure 2 shows the heated w tion of distance from the inlet at various inlet speeds. Part (a 70%, and part (c) is for 60% porous foam. A few trends can figure. The wall temperature decreases with increasing inlet a fact that convection heat transfer in the foam is more efficien where  The following governing equations were solved using ANSYS-Fluent: Continuity: Momentum: where ɛ and k are the porosity and permeability of the foam, respectively.
Energy equation for fluid: Energy equation for solid: where and are the velocity components in the flow direction (x) and perpendicular to the flow direction (y), and is the pressure. and µ are the density and viscosity of the fluid, respectively. is the form drag coefficient for the foam. and are the effective thermal conductivities of the fluid and the solid, respectively. ɛ and are the porosity and permeability of the foam, respectively. ℎ is the interfacial heat transfer coefficient, and is the interfacial surface area density (m 2 /m 3 ) inside the foam. and are the fluid and solid temperatures, respectively.
The following boundary conditions were imposed: At the channel's inlet: = , = 0, and = At the heated walls: ,, = − = − , = = 0 For symmetry: At the channel outlet: = 0 and Results showed that the higher the ppi, the higher was the heat-removal capability of the foam. This is because the surface-area density of the foam increases with increasing ppi. Even though 100 ppi can be achieved, commercial aluminum foam is usually available at no more than 40 ppi. Therefore, 40-ppi aluminum foam was chosen to be further investigated in terms of porosity. Figure 2 shows the heated wall's temperature as a function of distance from the inlet at various inlet speeds. Part (a) is for 90%, part (b) is for 70%, and part (c) is for 60% porous foam. A few trends can be identified based on this figure. The wall temperature decreases with increasing inlet air speed. This is due to the fact that convection heat transfer in the foam is more efficient at a higher coolant speed. and k are the porosity and permeability of the foam, respectively.
Energy equation for solid : k s ∂ 2 T s ∂y 2 + where u and v are the velocity components in the flow direction (x) and perpendicular to the flow direction (y), and p is the pressure. ρ and µ are the density and viscosity of the fluid, respectively. f is the form drag coefficient for the foam. k f and k s are the effective thermal conductivities of the fluid and the solid, respectively. Results showed that the higher the pp of the foam. This is because the surface-area ppi. Even though 100 ppi can be achieved, ble at no more than 40 ppi. Therefore, 40-p and K are the porosity and permeability of the foam, respectively. h s f is the interfacial heat transfer coefficient, and σ is the interfacial surface area density (m 2 /m 3 ) inside the foam. T f and T s are the fluid and solid temperatures, respectively.
The following boundary conditions were imposed: At the channel s inlet : u = u in , v = 0, and T f = T in (5) At the heated walls : For symmetry : At the channel outlet : p = 0 and ∂T f ∂x Results showed that the higher the ppi, the higher was the heat-removal capability of the foam. This is because the surface-area density of the foam increases with increasing ppi. Even though 100 ppi can be achieved, commercial aluminum foam is usually available at no more than 40 ppi. Therefore, 40-ppi aluminum foam was chosen to be further investigated in terms of porosity. Figure 2 shows the heated wall's temperature as a function of distance from the inlet at various inlet speeds. Part (a) is for 90%, part (b) is for 70%, and part (c) is for 60% porous foam. A few trends can be identified based on this figure. The wall temperature decreases with increasing inlet air speed. This is due to the fact that convection heat transfer in the foam is more efficient at a higher coolant speed. The wall's temperature increases slightly in the flow direction as air absorbs heat while moving through the foam. The model with 60% porosity shows the lowest range of maximum wall temperatures. In the 60% porosity case, the maximum plate temperatures reached the safe operating region of PEMs at an inlet velocity of about 1.5 m/s or higher. The operationally safe temperature range for PEMs is below 90 • C or 363 K [24]. It should be noted here that these results were for geometry and heat flux that are different from the stack to be investigated in this paper. Nonetheless, the results mentioned above give guidance to the porosity and pore density of the foam (60% and 40 ppi) that is recommended for designing foam-based cooling systems for actual fuel-cell stacks.

Stack Cooling Simulations
The investigation then moved to study the above metal foam (40 ppi and 60% porosity) for cooling of a common commercial stack, EFC-100-03-6-ST made by ElectroChem Inc., Woburn, MA, USA. This stack is intended for use in the automotive industry and power backup systems. It has three PEM fuel cells, each having an active area of 100 cm 2 . The total output power of the stack is 500 W. Currently, the stack is water-cooled to maintain the operating temperature of 60-90 • C. The water is pumped through serpentine flow fields in four channels; each channel is 1.78 mm × 0.89 mm × 73.66 mm (0.07 × 0.035 × 2.9 ). Table 1 lists the key features and the specifications of this stack. For simulation purposes, the geometry of the water-cooling field of the commercial stack, EFC-100-03-6-ST, was built in a 3D physical design in SolidWorks and imported into ANSYS-CFD. However, all shoulders between the serpentine small channels in the bipolar plate were removed, which made the cooling field one large channel having dimensions (L × W × H) 100 mm × 73.66 mm × 3.56 mm. The new design for air cooling stipulates that this large channel is filled with metal foam in order to dissipate heat more efficiently to air that will flow through it, as shown in Figure 3. In the simulation, the actual fuel cells including the bipolar plates were replaced by 100 mm × 100 mm × 7.12 mm aluminum plates. One large air-cooling channel was machined along with each plate. The depth of this channel was 3.65 mm.
The heat produced by each fuel cell was applied to each metal-foam cooling channel in the form of constant heat flux. Since the total heat removal requirement for the stack is 500 W, as indicated by the manufacturer, this amount was divided equally among the fuel cells, which resulted in a heat flux of 1.55 W/cm 2 applied to each side of an internal cell and on only one side of an end cell. Each cooling channel was filled with metal foam having the properties shown in Table 2. A drawing of the 3-dimensional heat transfer problem is shown in Figure 4.    The following governing equations were solved using Continuity: where ɛ and k are the porosity and permeability of the foa Energy equation for fluid: Energy equation for solid: where and are the velocity components in the flow d to the flow direction (y), and is the pressure. and µ a the fluid, respectively. is the form drag coefficient for t fective thermal conductivities of the fluid and the solid, re rosity and permeability of the foam, respectively. ℎ is th ficient, and is the interfacial surface area density (m 2 /m 3 ) the fluid and solid temperatures, respectively. The following boundary conditions were imposed: Results showed that the higher the ppi, the higher w of the foam. This is because the surface-area density of the ppi. Even though 100 ppi can be achieved, commercial alu ble at no more than 40 ppi. Therefore, 40-ppi aluminum investigated in terms of porosity. Figure 2 shows the heate tion of distance from the inlet at various inlet speeds. Pa 70%, and part (c) is for 60% porous foam. A few trends c figure. The wall temperature decreases with increasing in fact that convection heat transfer in the foam is more effic

Numerical Solution
The governing equations were solved invoking the local thermal nonequilibrium in ANSYS-CFD. This assumption allows the local temperatures of the fluid and the solid in the foam to be different, and hence there is convection heat transfer from the solid to the fluid. The governing equations for this 3D case are: where The following governing equations were solved using ANSYS-Fluent: Continuity: where ɛ and k are the porosity and permeability of the foam, respectively.
Energy equation for fluid: Energy equation for solid: where and are the velocity components in the flow direction (x) and perpendicular to the flow direction (y), and is the pressure. and µ are the density and viscosity of the fluid, respectively. is the form drag coefficient for the foam. and are the effective thermal conductivities of the fluid and the solid, respectively. ɛ and are the porosity and permeability of the foam, respectively. ℎ is the interfacial heat transfer coefficient, and is the interfacial surface area density (m 2 /m 3 ) inside the foam. and are the fluid and solid temperatures, respectively.
The following boundary conditions were imposed: At the channel's inlet: = , = 0, and = At the heated walls: ,, = − = − , = = 0 For symmetry: At the channel outlet: = 0 and = = 0 Results showed that the higher the ppi, the higher was the heat-removal capability of the foam. This is because the surface-area density of the foam increases with increasing ppi. Even though 100 ppi can be achieved, commercial aluminum foam is usually available at no more than 40 ppi. Therefore, 40-ppi aluminum foam was chosen to be further investigated in terms of porosity. Figure 2 shows the heated wall's temperature as a function of distance from the inlet at various inlet speeds. Part (a) is for 90%, part (b) is for 70%, and part (c) is for 60% porous foam. A few trends can be identified based on this figure. The wall temperature decreases with increasing inlet air speed. This is due to the fact that convection heat transfer in the foam is more efficient at a higher coolant speed. and k are the porosity and the permeability of the foam, respectively, and → u is the velocity vector.
Energy equation for fluid : Energy equation for solid : The following boundary conditions were imposed: At the channel s inlet : u = u in , v = w = 0 and T f = T in (13) For an insulated wall : For a heated wall : For symmetry : At the channel outlet : p = 0 and ∂T f ∂x The equations were solved for inlet air velocities from 0.2 to 1.6 m/s. The viscous and inertial resistances (correspond to the last two terms in Equation (10)) needed to be entered for metal foam. They were calculated according to the following equations given by ANSYS-CFD: where D p is the pore diameter of the foam. Both interfacial surface area density and heat transfer coefficient were calculated and manually entered to solve the energy equations. Dukhan and Patel [48] calculated the interfacial area density for 40 PPI of metal foam. They used a correlation provided by ERG Materials and Aerospace given as Equation (20).  where and are the velocity components in the flow to the flow direction (y), and is the pressure. and µ the fluid, respectively. is the form drag coefficient for fective thermal conductivities of the fluid and the solid, rosity and permeability of the foam, respectively. ℎ is ficient, and is the interfacial surface area density (m 2 /m ) + 3580 (20) where σ is the surface area per unit volume of the foam. Kuwahara et al. [34] provided the following correlation for the coefficient of interfacial convective heat transfer in metal foam:    where Re p is the Reynolds number based on the pore diameter. This correlation was used to calculate the heat transfer coefficient.

Meshing and Mesh Independence
The computational domain was meshed using three different meshes: coarse, medium, and fine, as shown in Figure 5a-c, respectively. All meshes employed tetrahedron elements. A mesh-independence study was carried out on these meshes. The results are shown in Table 3. As can be seen in this table, the difference between the fine and medium meshes was very small: the difference between the obtained pressure drop for the fine and medium meshes is only 0.11%. Therefore, the medium mesh was used in subsequent analysis, in order to save time. The meshed geometry using the medium mesh is shown in Figure 6.   Medium-smoothing-quality meshes were applied to all designed geometries with a de-feature size of 5 × 10 −3 mm. Mesh quality and smooth transition inflation were checked with 0.272 of the transition ratio. The smooth transition option was recommended by ANSYS based on the geometry of the model. This option calculates the local initial height of tetrahedral elements according to the total height of the elements, so that the rate of volume change, and the transition between each layer, are smooth. The transition ratio is the volume-based growth rate of the last tetrahedral elements. It works best for planar surfaces. ANSYS automatically assigns a transition ratio of 0.272 for a smooth transition. Double precision and serial processing options were selected. For pressure-velocity, the coupling simple scheme solution method was selected. Least squares cell-based gradient spatial discretization, PRESTO solution method for pressure, and second-order upwind for momentum and energy were applied. For all monitors for checking convergence, absolute criteria were chosen to be equal to 10 −6 . Figure 7 shows 2-dimensional contour plots of the temperature inside the foam at different air speeds. As the air speed increases, the temperature along the foam decreases. Also, the temperature at the bottom is maximum due to proximity to the heated plate. As expected, the temperature increases along the plate as air absorbs more heat.

Experimental Stack Model
An aluminum model was constructed in order to experimentally verify the simulation results. For that reason, this model was identical to the geometry used in the simulation. As in simulations, all shoulders between air channels in the bipolar plates were removed forming one large channel in each bipolar plate, Figure 8a. Commercial aluminum foam was mechanically compressed into each channel, Figure 8b.

Metal Foam and Test Section
Aluminum foam having 40 ppi with a porosity of 90.3% was obtained commercially [49]. The foam needed to be mechanically compressed in order to achieve 60% porosity (same porosity as in the simulation). The compression was done according to Boomsma et al. [50]. To house the experimental model, a test section was constructed of Plexiglas pieces glued together as shown in Figure 9. The inlet of the test section had dimensions 3.98 cm × 7.62 cm (1.57 in × 3 in). In order to get a uniform velocity profile entering the experimental model (as in the simulations), a piece of metal foam was inserted into the inlet of the test section at 3.81 cm (1.5 in) from the experimental stack. The outlet of the test section was enlarged to a square cross-section of 12.50 cm in order to be connected to the inlet of a suction machine by latches. Twelve holes were drilled on the left side of the test section, through the Plexiglas, for inserting thermocouples for measuring plate temperatures at various locations, Figure 9. One other hole was drilled on the right side of the test section for heater wires.

Experimental Setup
The larger opening of the test section was connected to a suction unit (SuperFlow SF-600) to form an open-loop wind tunnel, as shown in Figure 10. A connecting duct between the inlet of the suction unit and the outlet of the test section was made of Plexiglas and measured 60 cm in length. The suction unit could force ambient air to flow through the stack model and across the tunnel. A speed controller on the suction unit was used to control air speed through the tunnel.
Beads of ready-made insulated thermocouples (type K) having a diameter of 0.08 mm were inserted in each hole from one slide of the test section. The other ends of the thermocouples were connected to a data acquisition module (Omega OMB-DAQ-56), which was connected to a computer. The module had multiple channels and could receive voltage signals from thermocouples and other transducers. A Honeywell industrial vacuum pressure transducer model PX2EG1XX001BAAAX measured absolute vacuum pressure in the range 0 to 1 bar downstream from the test section. The pressure drop through the stack model was then calculated as the difference between the vacuum pressure after the model and the atmospheric pressure at the inlet, which was practically constant during experimental runs at 98.69 kPa, approximately. A hot-wire anemometer manufactured by Kanomax was used to measure the average air speed in the tunnel. The measuring range of this device was 0.01 to 30 m/s with a resolution of 0.01 m/s and an accuracy of 0.015 m/s. This device also could measure the test section's outlet temperature, which was within its range (−20 to 70 • C). For each run, the suction unit was turned on, then the power input to each heater was adjusted to produce the same heat flux that was applied in the simulations (1.55 W/cm 2 ). The average velocity through the tunnel was adjusted to the desired value by controls on the suction unit. Before recording any measurements, steady-state conditions were reached, which took 25 min, approximately. At steady state, all thermocouple readings, as well as pressures before and after the test section, were recorded. The air speed and the temperature after the test section were also recorded. The measurements were taken three times and the average was considered in the data analysis. This was repeated for various average air speeds between 0.2 and 1.6 m/s.

Measurements Error
For the thermocouples measuring the temperatures in the stack model, the error was 0.4%. The test section's outlet air temperature was measured by the hot-wire anemometer, which had a resolution of 0.1 • C and an accuracy of ±0.5 • C. As for the pressure measurement, the pressure transducer had an error of ±2%. The root-sum-squares, Equation (22), was used to calculate the uncertainty in the pressure drop measurement [51]. A fixed error of 2% and a reading error of 0.25% were stated by the pressure transducer's manufacturer.
δ p = ± e f 2 + e r 2 (22) where δ p is the pressure uncertainty, e f is the fixed error, and e r is the reading error. The inlet velocity was measured by a digital nanometer. The uncertainty of 3% as a fixed error and 0.015 m/s as a reading error were reported by the device's manufacturer. The pressure uncertainty was obtained as 2.06%, and air velocity uncertainty was calculated to be 3%.

Local Temperatures
Figure 11a-c shows a comparison between the temperatures obtained experimentally and numerically. These plots are for an internal plate (plate 2), representing an internal cell in the stack. The temperatures are plotted against the inlet velocity at different axial distances from the air inlet. Part (a) is for location 2.5 cm, part (b) is for 5 cm, and part (c) is for 7.5 cm. The agreement between the numerical and experimental values for all locations is excellent. As expected, the temperature decreases as the velocity increases. This is due to the fact that the heat-transfer coefficient increases with velocity, and more heat is swept by the air from the heated plate. This is clearly shown by Equation (21) for the interstitial heat-transfer coefficient that is proportional to the Reynolds number (or the velocity) to the power 0.6. One can notice that the agreement improves as the velocity increases for all locations: the numerical and experimental temperatures for a velocity of 1.5 m/s coincide. This could be attributed to possibly better prediction of the volume-averaged porous-media model (Equations (18)- (21)) within Fluent at higher velocities, or the improved validity of the various input quantities calculated by Equations (18), (19) and (21). The temperatures at the three axial locations are about the same because of the efficient cooling of metal foam, and because of the relatively short length of the model in the flow direction (10 cm).  Figure 12 shows that the general behavior of the local temperature for an end plate is similar to that of an internal plate, meaning that the local temperature increases in the flow direction since the air heats up as it travels through the foam. In addition, the agreement between the experimental and numerical temperatures is very good, especially for the lower velocities 0.2 and 0.6 m/s (Figure 12a,b). For the higher velocity of 1 m/s, part (c) of the figure, there is a little difference between the experimental and numerical tempera-  Figure 12 shows that the general behavior of the local temperature for an end plate is similar to that of an internal plate, meaning that the local temperature increases in the flow direction since the air heats up as it travels through the foam. In addition, the agreement between the experimental and numerical temperatures is very good, especially for the lower velocities 0.2 and 0.6 m/s (Figure 12a,b). For the higher velocity of 1 m/s, part (c) of the figure, there is a little difference between the experimental and numerical temperatures, with the numerical analysis over-predicting the local plate temperatures. As expected, the local temperature decreases as the velocity increases due to increased convection rates at higher velocities. The temperature value for an end plate is lower than that of an internal plate because an end plate is heated from one side only. A useful comparison from a practical point of view is between internal and external plate temperatures. Such a comparison is shown in Figure 13; part (a) is for the low air velocity of 0.2 m/s, while part (b) is for 1.5 m/s. For the low velocity of 0.2 m/s, the internal and external plates have practically identical temperatures at the three locations along the flow direction. The temperature increases in the flow direction, as the air absorbs heat from the heated plates. For the high velocity of 1.5 m/s, the internal and external plate temperatures are close to each other, but there is a noticeable difference between them, as shown in part b of Figure 13. The internal plate temperature is slightly higher. This means that no special design is needed for an external cell in terms of cooling.  A useful comparison from a practical point of view is between internal and external plate temperatures. Such a comparison is shown in Figure 13; part (a) is for the low air velocity of 0.2 m/s, while part (b) is for 1.5 m/s. For the low velocity of 0.2 m/s, the internal and external plates have practically identical temperatures at the three locations along the flow direction. The temperature increases in the flow direction, as the air absorbs heat from the heated plates. For the high velocity of 1.5 m/s, the internal and external plate temperatures are close to each other, but there is a noticeable difference between them, as shown in part b of Figure 13. The internal plate temperature is slightly higher. This means that no special design is needed for an external cell in terms of cooling.

Pumping Power
The numerical and experimental results of the pressure drop were used to calculate the pumping power. Equation (23) was used to calculate the pumping power: where .
W represents the pumping power in watts needed to pump the air through the cooling system, Q is the volumetric flow rate across the foam (m 3 /s), and ∆p is the pressure drop between inlet and outlet (Pa). Figure 14 has a comparison between the numerical and experimental pumping power needed to force air to flow through the foam channels. As shown in the figure, the pumping power is low at low flow rates. Increasing the flow rate of the air substantially increases the pumping power. While experimental and numerical curves have the same behavior, the numerical result shows a lower trend. The considerable difference warrants some explanation. Considering the inputs for the reciprocal of the permeability and the Iner.Resis that were calculated by Equations (18) and (19) then entered into ANSYS-CFD, one realizes that these two terms are based on Ergun's correlation [52]. It is likely that the Ergun equation is not applicable to metal foam, especially because of the vast geometrical difference between the foam's structure and that of packed spheres, for which the correlation was developed.
The disagreement between the pressure drop in the foam obtained experimentally, and calculated numerically by applying Ergun's equation, is expected to be larger in the form drag term. Thus, the difference will increase with increasing air velocity through the foam, since the form drag is a strong function of velocity. These trends are reflected in the pumping power as shown in Figure 14. The maximum difference in pumping power as predicted by ANSYS-CFD and as obtained experimentally reached 22.4% at the highest velocity of 1.6 m/s. The end-user of ANSYS-CFD does not seem to have provisions to alter the way this package calculates pressure drop for metal foam.

Conclusions
A new design for cooling fuel-cell stacks was examined numerically and experimentally. The design employs metal foam in flow channels for air cooling. The purpose of the foam was to increase the surface area available for heat transfer. The commercial aluminum foam that was used had 40 pores per linear inch and a porosity of 60%. This low porosity was achieved by mechanical compression. Air speeds were in the range of 0.1-1.6 m/s. In the simulation, the local thermal nonequilibrium model was utilized. The cooling design mimics that of commercial PEM fuel-cell stacks having three cells. The heat dissipated by the cells was represented by heaters that provided constant heat flux. The local temperature in the flow direction for both internal and end plates was obtained. The required pumping power across the metal foam was also determined. The new cooling system could remove the 500 watts of waste heat generated by the stack, keeping the stack within the safe range of operating temperatures of 60-90 • C. The numerical and experimental temperatures agreed very well at low velocities, with the small difference between them becoming greater for high velocities. As expected, the temperature decreased as the velocity increased. The temperature value for an end plate was lower than that of an internal plate because an end plate was heated from one side. However, the internal and external plates' temperatures were generally close to each other. This means that no special design is needed for an external cell in terms of cooling. ANSYS-CFD employs the Ergun equation for entering the viscous and form drags and calculating the pressure drop through the foam. This correlation is not applicable to metal foam, and thus ANSYS-CFD pressure drop (and thus pumping power) predictions can significantly be different from experimental values. This calls for considering more appropriate correlations for pressure drop in metal foam to be used in ANSYS-CFD. In general, the results suggest that the design can be used for cooling fuel-cell stacks for a reasonable pumping power. The thickness of the metal-foam-filled cooling channels can possibly be reduced, which should be investigated in the future. Moreover, including the actual electro-chemical performance of the fuel-cell stack in the modeling would be a logical next step for verifying the efficacy of this proposed cooling design.

Data Availability Statement:
The data presented in this study are available in this article.