{Large scale model predictions on the effect of GDL thermal conductivity and porosity on PEM fuel cell performance} doi:10.5599/jese.413 223 J. Electrochem. Sci. Eng. 7(4) (2017) 223-235; DOI: http://dx.doi.org/10.5599/jese.413 Open Access : : ISSN 1847-9286 www.jESE-online.org Original scientific paper Large scale model predictions on the effect of GDL thermal conductivity and porosity on PEM fuel cell performance Obaid Ur Rehman1,, Amber Fishan Zafar2 1Chemical Engineering Department, University of Karachi, Karachi, Pakistan 2Automotive and Marine Engineering Department, NED University of Engineering & Technology, Karachi, Pakistan Corresponding authors E-mail: orehman@uok.edu.pk; Tel. +92-321-2536798 Received: July 25, 2017; Revised: September 24, 2017; Accepted: September 25, 2017 Abstract The performance of proton exchange membrane (PEM) fuel cell majorly relies on properties of gas diffusion layer (GDL) which supports heat and mass transfer across the membrane electrode assembly. A novel approach is adopted in this work to analyze the activity of GDL during fuel cell operation on a large-scale model. The model with mesh size of 1.3 million computational cells for 50 cm2 active area was simulated by parallel computing technique via computer cluster. Grid independence study showed less than 5% deviation in criterion parameter as mesh size was increased to 1.8 million cells. Good approximation was achieved as model was validated with the experimental data for Pt loading of 1 mg cm-2. The results showed that GDL with higher thermal conductivity prevented PEM from drying and led to improved protonic conduction. GDL with higher porosity enhanced the reaction but resulted in low output voltage which demonstrated the effect of contact resistance. In addition, reduced porosity under the rib regions was significant which resulted in lower gas diffusion and heat and water accumulation. Keywords Computational fluid dynamics; Parallel computing; Experimental validation Introduction The impact of immense utilization of fossil fuels in the automotive and power generation sector has drawn enormous damage on environment. The development of renewable and clean energy resources has been a key topic of researchers in past few decades. Among other energy conversion devices fuel cells have been widely recognized as a promising clean energy resource due to their high energy density and conversion efficiency [1,2]. PEMFC has emerged as a substitute for internal combustion engines in automotive sector as it requires low operating temperatures and offers quick http://dx.doi.org/10.5599/jese.413 http://www.jese-online.org/ mailto:orehman@uok.edu.pk J. Electrochem. Sci. Eng. 7(4) (2017) 223-235 THERMAL CONDUCTIVITY AND POROSITY ON PEM FUEL CELL 224 startup. Though yet the fuel cells have not been fully commercialized due to high manufacturing cost and lack of infrastructure for hydrogen storage and its utilization, but still it has a great potential to replace conventional energy conversion systems. The main components of PEMFC are represented in Figure 1. Air is injected from cathode side of the cell while hydrogen gas from anode. At the surface of anode side catalyst layer (CL), hydrogen oxidation reaction (HOR) occurs as shown in the reactions below. Simultaneously, the oxygen reduction reaction (ORR) occurs at the cathode side CL. The hydrogen ions (protons) from anode travel through the PEM and reach the cathode catalyst where they react with oxygen to produce water. Product heat is released as the reaction is exothermic. Figure 1. Schematic of PEM fuel cell The performance of PEMFC is highly dependent upon the thermal and transport characteristics of the porous CL and GDL [3]. GDLs are employed for distribution of reactant gases to the reaction sites of the catalysts which increase the diffusion capacity and enhance reaction rate. They also provide the pathway for removal of water and heat from CL to gas flow channels which in turn limit the concentration overpotential. To achieve high current densities, the GDL must be porous and allow for the flow of both water and reactant gases. It must also be thermally and electrically conductive for the flow of product heat and electric current in both in-plane and through plane directions[4]. Water droplets form at low operating temperature in fuel cell which block pores of GDL and reduce gas diffusivity and number of reaction sites at CL [5]. On the other hand, high water content (H2O/SO-3 ratio or ) in PEM promotes proton conductivity [6]. This trade-off between water content and reaction rate should be taken into consideration for an efficient operation of PEMFC. Favorable characteristics of GDL and CL are required for improved gas diffusivity and thermal management in the cell. Various studies have been done experimentally to analyze the physical and thermal effects of GDLs on PEMFC performance [7-10]. When compared with the experimental procedures numerical simulation provides relatively agile methods to design and analyze complex systems and it also offers access to diversified results. Currently various numerical models of PEMFC are available in the published literature covering the O. Rehman et al. J. Electrochem. Sci. Eng. 7(4) (2017) 223-235 doi:10.5599/jese.413 225 transport phenomena and electrochemical kinetics [11-21]. To analyze the performance of PEMFC at different operating conditions and design configurations various simulations have been produced but very few provide analysis of full GDL on large-scale model. Zhang et al. [22] investigated the effect of porosity of cathode side GDL on catalyst potential distribution and pressure drop along flow channels. They noticed that at porosity of 0.6 the potential was maximal but as the value of porosity was increased potential dropped, which indicated that the contact resistance impeded the transfer of current through GDL. Inamuddin et al. [23] studied a single channel, three dimensional model of PEMFC to evaluate its performance at different GDL porosities and thickness. They noticed a gradual increase in current density with GDL porosity. However, performance at porosity greater than 0.7 was not estimated which may indicate its limiting value. Khazaee et al. [24] developed a three dimensional model to investigate the effect GDL and membrane characteristics on performance of annular PEMFC. Their results suggested that high GDL porosity was not favorable at high current densities as it led to increased contact resistance. They also suggested high thermal conductivity of GDL to prevent PEM from drying. Alhazmi et al. [25] developed a 11-channel, three dimensional model to estimate the performance of PEMFC at different in-plane and through plane GDL thermal conductivities. Improved power density was observed for high GDL thermal conductivities which corresponded to low PEM temperature. Low temperature operations were favorable for low electrical and protonic resistance in GDL and PEM respectively. Fadzillah, Nee and Rosli [26] simulated a two dimensional model to investigate the distribution of oxygen on cathode side GDLs with different porosities and thicknesses. It was observed that porosity of GDL played a key role to facilitate the reactant to reach more reaction sites which resulted in improved PEMFC performance. Maslan et al. [27] developed a three dimensional single channel model to predict the performance of PEMFC with respect to GDL properties. Effect of porosity and PTFE content of GDL was analyzed. Their results showed that at low porosity the concentration overpotential dominated the PEMFC performance because the water droplets were trapped in GDL pores which resulted in reduced reaction sites. In this paper a simulation work is presented for a full cell model using commercial code of ANSYS FLUENT® with parallel computing technique. Effects of GDL thermal conductivity and porosity on PEMFC performance were investigated. Cell performance was analyzed with relation to water content and temperature across PEM. Oxygen concentration and reaction heat production at GDL- catalyst interface was also examined to observe the impact of GDL porosity. Model description The governing transport phenomena and reaction kinetics of PEMFC has been modelled in numerous works. Some good reviews on the model of PEMFC can be found in [28] and [29]. The geometrical and mathematical models are described in subsequent paragraphs. Geometric model The full cell geometry with 50 cm2 of active area and 45-channel serpentine flow design was developed using GAMBIT pre-processor. The model consists of seven layers as shown by the schematic in Figure 1. The geometry was meshed by hexahedral scheme as shown in Figure 2. The whole geometry comprised approximately 1.3 million computational cells. The photograph of the flow channels is presented in Figure 3. The geometrical dimensions are presented in Table 1. J. Electrochem. Sci. Eng. 7(4) (2017) 223-235 THERMAL CONDUCTIVITY AND POROSITY ON PEM FUEL CELL 226 Figure 2. Meshed geometry of PEMFC Figure 3. Engraved serpentine gas flow channels in bipolar plate Table 1. Model physical properties Dimension Length, mm Channel height 1 Channel length 70 Channel width 0.75 Rib width 0.82 GDL thickness 0.19 Catalyst layer thickness 0.01 Membrane thickness 0.0508 Current collector thickness 0.25 Mathematical model - governing equations A multiphase mathematical model was employed in present work. The main reaction at cathode side takes place at triple phase boundary. Hydrogen ions formed at anode CL travel to cathode CL through PEM where they are combined with oxygen gas to form water. Conservation equations In a finite volume method the basic equation for a conservation of a general property φ over a control volume in a steady state problem is given as [30]. VSA)n(ΓAu)n(ρ CV V AA ddd    (1) where n represents the vector normal to a differential surface dA and Γ is a diffusion coefficient. The left-hand side gives the convective flux and right-hand side gives the diffusive flux plus the generation or consumption of the property. The continuity equation for the model is given as     mSuερ t ρε     (2) where ρ is the density, ε is the porosity, u  is the fluid velocity vector and Sm denotes the mass production term. O. Rehman et al. J. Electrochem. Sci. Eng. 7(4) (2017) 223-235 doi:10.5599/jese.413 227 The momentum conservation equation is given as follows       uSuεμpεuuερ t uρε      (3) where p represents the pressure, and Su is the force per unit volume. For the species conservation following equation is employed       kk eff kk k SCDCuερ t εC     (4) where Ck is the concentration of species and eff kD represents the diffusivity, which is given for porous material by Bruggeman equation as follows k 1.5 m eff k DεD  (5) The species production term Sk is related to electrochemical reaction in the fuel cell which is given in the following form j kw, k R nF M S  (6) where k shows the specie, n is the electron transfer number, and subscript j used for anode or cathode. Energy conservation equation is given as follows       h eff STkhuερ t ρεh     (7) keff is the effective thermal conductivity of porous material. The energy source (Sh) is a sum of different source terms such as heat of reaction, ohmic loss and electric work and latent heat of evaporation for water, and can be given as Lohm 2 catan,catan,reacth hRIηR-h=S  (8) The equation for conservation of charge is given as   0e eff e  ΦSΦσ (9) where eΦ is the charge potential for membrane or solid phase, eff eσ is the ionic conductivity, ΦS is the source term which depends on exchange current density (A/cm2). The dependence of exchange current density on reactant concentration can be expressed by Butler-Volmer’s [13] equation. Simplified form of which known as Tafel equation is given below                                       RTFηα γ RTFηα γ Φ eII eII S / ref2 2ref catcat / ref2 2ref anan catcat cat anan an O O H H (10,11) where Iref is the reference exchange current density, γ is the concentration dependence factor, α is the transfer coefficient, F is Faraday constant, and η is the overpotential. PEM properties PEM proton conductivity ( memσ ) is related to membrane water content () and temperature by following equation           Teλσ 1 303 1 1268 mem 0.3260.514 (12) J. Electrochem. Sci. Eng. 7(4) (2017) 223-235 THERMAL CONDUCTIVITY AND POROSITY ON PEM FUEL CELL 228 The membrane water transport by osmotic drag is given by the following relation [11] )(2d drag IαJW  (13) The osmotic drag coefficient ( dα ) for proton conductivity of membrane also depends upon λ and is related by following equation 22 2.5 λ αd  (14) The back diffusion flux of water molecules through membrane can be related to λ as λDM M ρ JW  lOH m mdiff 2 (15) where mρ and mM are density and equivalent weight of dry membrane. The diffusion coefficient, lD , is represented by the following equation 632 1 303 1 2416 l 10)0.0006710.02640.332.563          λλλ(eD T (16) Water content of membrane was estimated using equation [11] as follows       1for1)1.4(14 1for3639.8517.180.043 32 aa aaaa λ (17) where a is water activity and is defined mathematically as: sat WV P P a  (18) The vapor pressure can be related to molar fraction and total pressure as follows: PP OHWV 2x (19) Model parameters and boundary conditions The simulations were carried out at steady state and non-isothermal conditions. The reaction parameters were set at 1 atm and 353 K (80°C). Output voltages were calculated at a fixed current density of 0.6 A/cm2. Flow rate for hydrogen gas was set at 0.4 SLPM (9.477×10-7 kg/s) and for air at 1.26 SLPM (2.82×10-5 kg/s). Both streams were entered with 100 % relative humidity. The model was calibrated in order to generate comparable solutions with the experimental results by tuning the reference exchange current density. The model parameters are presented in Table 2. Table 2. Model parameters Quantity Value Anode reference current density [31] 100 A m-2 Anode reference molar concentration [31] 0.04 kmol m-3 Anode concentration exponent [17] 0.5 Anode exchange coefficient [32] 0.5 Cathode reference exchange current density 0.00035 A m-2 Cathode reference molar concentration [31] 0.04 kmol m-3 Cathode concentration exponent [17] 1 Cathode exchange coefficient 0.6 Open circuit voltage 1.05 V Hydrogen reference diffusivity [31] 1.1028×10-4 m2 s-1 Oxygen reference diffusivity [31] 3.2348×10-5 m2 s-1 Water reference diffusivity [31] 7.65×10-5 m2 /s s-1 Nitrogen reference diffusivity [31] 3×10-5 m2 s-1 O. Rehman et al. J. Electrochem. Sci. Eng. 7(4) (2017) 223-235 doi:10.5599/jese.413 229 Material properties In order to solely recognize the effect of thermal conductivity and porosity of GDL its other properties were taken identical in all solutions. Isotropic model was adopted for specifying the values of viscous and electrical resistances. The properties of materials used in this model are given in Table 3. Table 3. Material properties Material Property Value Nafion Density 1968 kg m-3 Specific heat capacity 4188 J kg-1 K-1 Thermal conductivity, dry at 65°C 0.12 W m-1 K-1 Equivalent weight 1100 kg kmol-1 Toray Carbon Bulk density 440 kg m-3 Specific heat capacity 685 J kg-1 K-1 Electrical conductivity 1250 -1 m-1 Viscous resistance 1.02×1011 m-2 Platinum, Pt Thermal conductivity at 60°C 73 W m-1 K-1 Carbon support (Vulcan XC 72) Bulk density 264 kg m-3 Specific heat capacity 685 J kg-1 K-1 Thermal conductivity 7.63 W m-1 K-1 Electrical conductivity 400 -1m-1 Graphite Plate Density 1990 kg m-3 Specific heat capacity 710 J kg-1 K-1 Thermal conductivity 117 W m-1 K-1 Electrical conductivity 92600  -1m-1 Solver specification A finite volume based FLUENT® solver was implemented to solve the governing equations. The large scale computational domain was handled by parallel processing technique. The URFs are employed to control the solution of highly coupled equations. To achieve convergence the URFs were tuned to an optimum value in order to lower the residuals for each variable without large oscillations. About 700 iterations were performed to achieve converged solutions. Stopping criteria was set at a residual value of 1×10-6 for the equation of continuity and 1×10-5 for other variables (i.e. potential fields, water content, species outlet mole fractions etc.) which took more time to converge than the scaled residuals. Results and discussion Grid independence study To make sure that the large-scale model is independent of the meshing criteria the model was meshed with three different sizes as sown in Table 4. The current densities were calculated for each mesh size at a fixed potential of 0.67 V. About 8.5 % deviation in the criterion parameter was found when mesh size was increased from 1 to 1.3 million cells. The deviation reduced to about 4.4% when mesh size was further increased to 1.8 million cells. The study showed that the deviation diminished as the mesh became finer. Due to high computational load the grid size with 1.3 million computational cells was selected for all the calculations. J. Electrochem. Sci. Eng. 7(4) (2017) 223-235 THERMAL CONDUCTIVITY AND POROSITY ON PEM FUEL CELL 230 Table 4. Average current densities for different mesh sizes at 0.67 V Mesh Cells Average CD, A cm-2 Fine 1858725 0.425 Medium 1352475 0.408 Coarse 1099350 0.375 Experimental validation The simulation results were validated with experimental data by comparing the polarization curves, as shown in Figure 4. In-house experiments were carried out to generate a polarization curve. A standard test PEM fuel cell with 50 cm2 active area, developed by Electro Chem Inc. was used for experiments at 50°C temperature and atmospheric pressure. Toray carbon (TP 60) material was used as GDL and the catalyst layer was made with Pt nanoparticles on carbon substrate (0.2 mg Pt/mg Vulcan XC 72) with loading of 1 mg cm-2. The measurement of current density was carried out by galvanostatic control. Fuel cell testing system (FCTS) was employed for data acquisition. Constant gas flow rates were chosen throughout the experiments to achieve a minimum stoichiometry of 1.5 and 2 for hydrogen and oxygen at 0.8 A cm-2. Figure 4. Polarization and power density curves at 50 °C and 1 atm A good agreement was achieved as the numerical results followed the experimental curve as shown in Figure 4. However, at high current densities above 0.4 A cm-2 the curve for numerical solution deviated due to inadequacy of the model to reproduce the actual behavior of PEMFC. At high current density, the concentration overpotential dominate the cell performance as water droplets blocks the diffusion pathways for reactants to reach reaction sites as depicted by the experimental curve. However, the discrepancies at high current density could be minimized by decreasing the URFs for saturation source term and membrane water content but the number of iterations would be very high and would result in stalling of convergence. O. Rehman et al. J. Electrochem. Sci. Eng. 7(4) (2017) 223-235 doi:10.5599/jese.413 231 Effect of GDL thermal conductivity Commercially available data [32] for four types of carbon fiber paper used as GDL was incorporated namely Toray carbon, E-Tek, Spectracarb and Sigarcet. The effect of thermal conductivity of GDL on the water content and temperature of PEM is illustrated by contours shown in Figure 5 and 6 which represent the iso-surface at the middle of PEM. Figure 5 shows the overall decrease in water content () of PEM with thermal conductivity of GDL. The gradual increase of water content (Figure 5) along the channel clearly shows the saturation of the reactant gas with product water generated by ORR which is removed through GDL. PEM drying is evident in Figure 5 as GDL thermal conductivity is decreased from 1.7 to 0.16 W / m K. The drying phenomenon is related to high temperature operation as shown in Figure 6. The effect of GDL thermal conductivity on PEM temperature is distinctively revealed in Figure 6. PEM temperature increased as GDL thermal conductivity is reduced from 1.7 to 0.16 W / m K. Moreover, a low temperature profile is eminent in both Figures making an interdigitated pattern in the regions under the rib where water accumulates due to low gas flow. The temperature of membrane plays a significant role in the PEMFC performance. High temperature causes drying and consequently lowers the proton conductivity ( memσ ). Moreover, enhanced proton conductivity of PEM effectively reduced the ohmic overpotential which resulted in high output voltages as depicted in Figurer 7 which shows the calculated output voltages for each GDL at fixed current density of 0.6 A cm-2. Figure 5. Contours of PEM water content () at GDL thermal conductivity of (a) 1.7 W/m K; (b) 1.4 W / m K; (c) 0.6 W / m K; (d) 0.16 W / m K H20/SO-3 x / mm x / mm x / mm x / mm x / mm x / mm J. Electrochem. Sci. Eng. 7(4) (2017) 223-235 THERMAL CONDUCTIVITY AND POROSITY ON PEM FUEL CELL 232 Figure 6. Temperature, K contours of PEM at GDL thermal conductivity of (a) 1.7 W / m K; (b) 1.4 W / m K; (c) 0.6 W / m K; (d) 0.16 W / m K Figure 7. Output voltages at different GDL thermal conductivities Figure 8. Output voltages at different GDL porosities Effect of GDL porosity Among other functions of GDL the distribution of reactant gases to the active surface area of CL is very critical. The diffusivity of gases highly depends upon the porosity of material as shown in Equation 5. A highly porous GDL will lead to better transport of reactant gases. On the other hand, high porosity promotes contact resistance between GDL and bipolar plate which reduces electrical x / mm x / mm x / mm x / mm O. Rehman et al. J. Electrochem. Sci. Eng. 7(4) (2017) 223-235 doi:10.5599/jese.413 233 conductivity [33]. Therefore, an optimum porosity is desired for efficient process. Figure 8 shows the output voltages that are generated by the simulations at different GDL porosities obtained from commercially available data. The results depict that the output voltage was increased as the porosity of GDL was lowered from 0.88 to 0.63 due to improved electrical conductivity which indicate reduced contact resistance. The effect of GDL porosity on reaction rate at GDL-catalyst interface on cathode side of PEMFC is shown in Figures 9 and 10. The GDL was treated as a single domain with uniform porosity. A subtle change is noticeable in both contours of Figure 9 which represents mole fraction of oxygen that increased from 0.1 to 0.22 along the channel. The results showed that oxygen was diffused at a higher rate in GDL with 0.88 porosity (Figure 9(a)) than with 0.63 porosity (Figure 9(b)). Furthermore, pathways of gas to CL are obstructed in the regions under the rib which exhibit the accumulation of water. Figure 9. Oxygen mole fraction at cathode side GDL-catalyst interface at GDL porosity of (a) 0.88; (b) 0.63 The heat generated in the ORR is illustrated in Figure 10. A considerable change in both diagrams 10(a) and 10(b) is noticed which signifies the dependence of ORR rate on GDL porosity. Product heat was increased along the channel corresponding to the same phenomenon represented in Figure 9. Heat accumulated under rib while in the channel it was swept away by the gas stream. This suggests that an optimum porosity of GDL is highly decisive in efficient performance of PEMFC. x / mm x / mm J. Electrochem. Sci. Eng. 7(4) (2017) 223-235 THERMAL CONDUCTIVITY AND POROSITY ON PEM FUEL CELL 234 Figure 10. Reaction heat source (W/cm3) at GDL-catalyst interface at GDL porosity of (a) 0.88; (b) 0.63 Conclusions This study provides a numerical investigation of the effect of thermal and transport properties of GDL on the performance of PEMFC. The mathematical model was validated with experimental results by the comparison of the polarization and power density curves. Contours of different iso- surfaces were presented to show the impact of thermal conductivity and porosity of GDL on the electrochemical behavior of PEMFC. The results provide a substantial basis in understanding the actual phenomenon occurred inside the complex system of fuel cell. It is reported that higher thermal conductivity of GDL led to improved proton conductivity of PEM by maintaining low temperatures. It is also found that higher porosity of GDL promoted the reactant gas transport to CL but at the same time raised contact resistance which resulted in lower electrical conductivity. The compressive force in the cell also affected the performance by impeding the gas flow through porous GDL. These results suggest that the optimum GDL characteristics and compressive force is required for maximum efficiency of PEMFC. Further study can be done by incorporating other GDL properties like gas permeability, thickness, electrical conductivity and interface resistance between the layers adjacent to GDL. x / mm x / mm O. Rehman et al. J. Electrochem. Sci. Eng. 7(4) (2017) 223-235 doi:10.5599/jese.413 235 References [1] M. F. B. Sunden, Transport Phenomena in Fuel Cells (Developments in Heat Transfer). WIT Press, Southampton, UK, 2005. [2] N. Sammes, Fuel Cell Technology Reaching Towards Commercialization. Springer, London, UK, 2006. [3] X. L. David P. Wilkinson, Jiujun Zhang, Rob Hui, Jeffrey Fergus, Proton Exchange Membrane Fuel Cells: Materials Properties and Performance. CRC Press, Florida, USA, 2009. [4] M. M. Mench, Fuel Cell Engines. John Wiley & Sons, Hoboken, USA, 2008. [5] U. Pasaogullari, C. Y. Wang, J. Electrochem. Soc. 151(3) (2004) A399–A406. [6] P. C. Rieke, N. E. Vanderborgh, J. Memb. Sci. 32(2) (1987) 313–328. [7] N. Parikh, J. S. Allen, R. S. Yassar, Fuel Cells 12(3) (2012) 382–390. [8] R. Banerjee, J. Hinebaugh, H. Liu, R. Yip, N. Ge, A. Bazylak, Int. J. Hydrogen Energ. 41(33) ( 2016) 14885–14896. [9] Y. Zhang, A. Verma,R. Pitchumani, Int. J. Hydrogen Energy 41(20) (2016) 8412–8426. [10] Y. Gao, A. Montana, F. Chen, J. Power Sources 342 (2017) 252–265. [11] T. E. Springer, T. A. Zawodzinski, S. Gottesfeld, J. Electrochem. Soc. 138(8) (1991) 2334–2342. [12] T. V NguyenR. E. White, J. Electrochem. Soc. 140(8) (1993) 2178–2186. [13] S. Um, C.‐Y. Wang,K. S. Chen, J. Electrochem. Soc. 147(12) (2000) 4485–4493. [14] H. Ju, H. Meng, C. Y. Wang, Int. J. Heat Mass Transf. 48 (2005) 1303–1315. [15] K. W. LumJ. J. McGuirk, J. Power Sources 143(1) (2005) 103–124. [16] Y. Wang and C.-Y. Wang, J. Electrochem. Soc. 153(6) (2006) A1193–A1200. [17] Y. M. FerngA. Su, Int. J. Hydrogen Energ. 32(17) (2007) 4466–4476. [18] A. A. Shah, G.-S. Kim, P. C. Sui, and D. Harvey, J. Power Sources 163(2) (2007) 793–806. [19] A. D. LeB. Zhou, J. Power Sources 182(1) (2008) 197–222. [20] Y. Wang, S. Basu, C.-Y. Wang, J. Power Sources 179(2) (2008) 603–617. [21] J. J. BaschukX. Li, Appl. Energy 86(2) (2009) 181–193. [22] Z. Zhang, X. Wang, X. Zhang, L. Jia, J. Fuel Cell Sci. Technol. 5(3) (2008) 31007–31009. [23] Inamuddin, T. A. Cheema, S. M. J. Zaidi, S. U. Rahman, Renew. Energ. 36(2) (2011) 529–535. [24] I. Khazaee, M. Ghazikhani, M. N. Esfahani, Appl. Surf. Sci. 258(6) (2012) 2141–2148. [25] N. Alhazmi, D. B. Ingham, M. S. Ismail, K. J. Hughes, L. Ma, M. Pourkashanian, Int. J. Hydrogen Energ. 38(1) (2013) 603–611. [26] D. M. Fadzillah, C. L. Nee, M. I. Rosli, Int. J. Mechan. Mechatron. Eng. 15(5) (2015) 83–89. [27] N. H. Maslan, M. M. Gau, M. S. Masdar, M. I. Rosli, J. Eng. Sci. Technol. 11(1) (2016) 85–95. [28] N. Djilali, Energy 32(4) (2007) 269–280. [29] C. Siegel, Energy 33(9) (2008) 1331–1352. [30] W. M. H. Versteeg, An Introduction to Computational Fluid Dynamics: The Finite Volume Method, Second. Pearson Education Limited, London, UK, 2007. [31] S. Kamarajugadda S. Mazumder, Comput. Chem. Eng. 32(7) (2008) 1650–1660. [32] A. Pfrang, D. Veyret, G. Tsotridis, Convection and Conduction Heat Transfer, Intech, Rijeka, Croatia, 2011, p. 215. [33] Y.-X. Huang, C.-H. Cheng, X.-D. Wang, J.-Y. Jang, Energy 35(12) (2010) 4786–4794. ©2017 by the authors; licensee IAPC, Zagreb, Croatia. This article is an open-access article distributed under the terms and conditions of the Creative Commons Attribution license (http://creativecommons.org/licenses/by/4.0/) http://creativecommons.org/licenses/by/4.0/)