MathJax.Hub.Config({tex2jax: {inlineMath: [['$','$'], ['\$','\$']]}}); Efficient Mass Transport and Electrochemistry Coupling Scheme for Reliable Multiphysics Modeling of Planar Solid Oxide Fuel Cell Stack
 Chinese Journal of Chemical Physics  2017, Vol. 30 Issue (2): 139-146

#### The article information

Ang Li, Zi-jing Lin

Efficient Mass Transport and Electrochemistry Coupling Scheme for Reliable Multiphysics Modeling of Planar Solid Oxide Fuel Cell Stack

Chinese Journal of Chemical Physics, 2017, 30(2): 139-146

http://dx.doi.org/10.1063/1674-0068/30/cjcp1610198

### Article history

Accepted on: November 5, 2016
Efficient Mass Transport and Electrochemistry Coupling Scheme for Reliable Multiphysics Modeling of Planar Solid Oxide Fuel Cell Stack
Ang Li, Zi-jing Lin
Dated: Received on October 23, 2016; Accepted on November 5, 2016
Hefei National Laboratory for Physical Sciences at the Microscales & CAS Key Laboratory of Strongly-Coupled Quantum Matter Physics, Department of Physics, University of Science and Technology of China, Hefei 230026, China
Author: Zi-jing Lin, E-mail:zjlin@ustc.edu.cn, Tel:+86-551-63606345
Abstract: A multiphysics model for a production scale planar solid oxide fuel cell (SOFC) stack is important for the SOFC technology, but usually requires an unpractical amount of computing resource. The major cause for the huge computing resource requirement is identified as the need to solve the cathode O2 transport and the associated electrochemistry. To overcome the technical obstacle, an analytical model for solving the O2 transport and its coupling with the electrochemistry is derived. The analytical model is used to greatly reduce the numerical mesh complexity of a multiphysics model. Numerical test shows that the analytical approximation is highly accurate and stable. A multiphysics numerical modeling tool taking advantage of the analytical solution is then developed through Fluent®. The numerical efficiency and stability of this modeling tool are further demonstrated by simulating a 30-cell stack with a production scale cell size. Detailed information about the stack performance is revealed and brie y discussed. The multiphysics modeling tool can be used to guide the stack design and select the operating parameters.
Key words: Simulation     Mesh setting     Analytical model     Computational efficiency     Numer-ical stability
Ⅰ. INTRODUCTION

Solid oxide fuel cell (SOFC) has considerable potential in future clean energy industry due to its high efficiency and low level of pollutant emissions [1-5]. The planar design is competitive due to its high volumetric power density and relatively low cost of manufacturing. For a widespread adoption of the technology, substantial improvement of SOFC stack performance to approach its theoretical limit is required. A theoretical modeling tool that can reliably analyze and predict the operating characteristics of production scale SOFC stacks is critically important for advancing the technology.

The performance of SOFC stack are governed by the interplay of a number of physical processes such as fuel and air transport in flow channels and porous electrodes, electrochemical reaction at the triple-phase-boundary, current, and heat conduction in the stack. These physical processes are strongly coupled, resulting in considerable complexity about the stack operating behavior. Reliable improvement of stack performance can only be made based on in-depth understanding of the coupled multiphysics phenomenon. Multiphysics fully coupled models of production scale SOFC stacks are necessary for providing credible predictions of the stack performance in practical scenarios.

The existing multiphysics numerical models of planar SOFC (pSOFC) stacks can be roughly classified into two categories. The first category or the early model was characterized by in-house developed software employing simplified SOFC geometries and coarse grids [6-14]. The second category or the recent model focuses on developing user defined functions (UDFs) of commercial computational fluid dynamics (CFD) or finite element codes to take advantage of the well developed commercial tool [15-52]. Even though the commercial code based models were reported more than a decade ago [31, 46], the progress made since then is quite limited.

For example, most models are in fact a single cell or even a partial single cell model [17-39], even though some may be claimed as stack models [31, 33, 35-37]. The main progress is reflected by using refined grids to capture the true structure of an SOFC cell. However, it is well recognized that a cell model is far from representative of an operating stack [54], but to extend similar grids to a production scale pSOFC stack demands a huge computing resource and is unrealistic. Efforts of developing stack models all used unrealistic stack geometries as well as partially coupled multiphysics descriptions [40-52]. Simulations are often performed on stacks with no gas manifolds and interconnect ribs, or with very small cell sizes [47, 52]. Moreover, the solution for the species transport in the electrodes is omitted by assuming a direct relationship between the concentration polarizations and the local current [47, 52]. For example, the 5-and 10-cell stack models recently developed by Sudaprasert et al. [38] do not solve the gas transport in the electrodes and use an active cell area of only 11.6 cm2 that is one order of magnitude smaller than the production cell size. Stack models with realistic geometries have been developed, but only the physical processes that are easy to couple are actually solved. Such stack models often avoid the electrochemistry and charge transport modeling by using the average current density and the average heat source [41, 42, 44, 48, 49]. Due to the strong coupling among the flow, electrochemistry and thermal conduction, however, the results of the physical fields obtained through simplified stack geometries or partially coupled scheme are often quite different from the reality. To summarize, none of the models reported so far possess the true features of a multiphysics pSOFC stack model that solves the full couplings of mass flows, chemical and electrochemical reactions, thermal and current conductions in a realistic pSOFC stack geometry.

To fully realize the power of numerical models for the in-depth understanding and development of operating SOFC stacks, the development of realistic geometry based multiphysics stack model is mandatory. In this report, a major technical obstacle impeding the development of a full multiphysics pSOFC stack model is analyzed. The obstacle is overcome by designing an algorithm for the electrochemistry and mass transport coupling. The algorithm is capable of drastically reducing the grid complexity of stack model and tested to be highly accurate for practical requirements. The algorithm is implemented as a UDF interacting with the commercial CFD software of Fluent®, enabling a creation of an authentic multiphysics model for pSOFC stack. The fully coupled multiphysics stack model is shown to be numerically highly efficient.

Ⅱ. METHOD A. Analysis of technical problem in developing multiphysics SOFC stack model

Understanding the difficulties of creating an authentic multiphysics model for production scale SOFC stacks is a necessary step towards solving the problem. Therefore, a brief analysis about the major difficulty is presented first. The analysis is based on realistic stack geometry parameters.

A production scale SOFC stack may consist of, 30 cell units. Each cell unit may include about 30 gas-channel-interconnect-rib pitches consisting of the cathode (positive)-electrolyte-anode (negative) (PEN) assembly, the air and fuel channels and interconnects. One major difficulty of developing a multiphysics SOFC stack model is about the description of O2 transport and electrochemistry in the interconnect rib covered cathode area, as schematically illustrated in Fig. 1. Due to the relatively low O2 diffusion coefficient and the thinness of cathode layer (50-200 μm) for O2 flux, the O2 concentration decreases rapidly away from the gas channel, as shown in Fig. 1(c). To describe the O2 transport in the cathode-rib area requires a very fine mesh setting, as illustrated in Fig. 1(b) and (d). Moreover, the strong coupling between the rapidly varying O2 concentration and the electrochemistry may cause instability in the numerical simulation.

 FIG. 1 Representative structural domain of SOFC stack and the associated O2 distribution and mesh setting. (a) 2D cross section of a repeating structural unit in a SOFC planar stack, (b) mesh setting in a repeating structural unit, (c) O2 distribution in the cathode-rib region, and (d) mesh setting for the cathode-rib domain.

To solve the O2 transport and electrochemistry process with reasonable accuracy, tests show that a minimum of 2000 grid points are required for a 2D repeating pitch unit model (Fig. 1) with a pitch width of 4 mm and a cell height of 4.1 mm. For a 3D repeating pitch unit with a channel length of 120 mm, about 105 grid points are necessary. A 30-cell stack with 30 repeating pitch units in each cell asks for 90 million grid points. Considering the grids for the stack flow channels, frame, walls, etc., and the requirement of grid continuity, the total number of grid points is about 150 millions, unacceptably large for a multiphysics fully coupled simulation. Consequently, multiphysics simulations are only reported for single-cell stack models or stacks with a low number of cells and small cell sizes.

B. Determination of cathode O2 transport and current production

As discussed above, a fine grid setting is required due to the need to describe the O2 transport and electrochemistry in the cathode-rib area. However, the contribution of the cathode-rib area to the physical processes (O2 distribution, current generation, heat production, etc.) is relatively small [18]. That is, major resources are devoted to solutions of minor influences. Such an undesirable situation points to the need of finding a method to estimate the contribution of the cathode-rib area with reasonable accuracy and eliminate the associated fine mesh requirement. Such a method is derived below.

Denoting the O2 concentration at the cathode-electrolyte interface, as CO2ce and at the cathode-channel interface, as CO2cc. As the cathode is very thin, the O2 flux along the cathode thickness direction is mainly driven by O2 concentration gradient and should be equal to the amount consumed by the local electrochemistry at j/4F. Therefore, we have,

 (1)

where d is the cathode thickness, DO2eff the effective O2 diffusion coefficient in the cathode, j the local current density, F the Faraday constant.

The equality of O2 flux and electrochemical consumption is also valid for the cathode-rib region. As CO2ce decays rapidly inside the cathode-rib region, however, the air pressure gradient in the region may be substantial and the convective flux may be not negligible and should be considered. For a unit length along the channel flow direction (perpendicular to the 2D cross section shown in Fig. 1), the relationship between the O2 molar flow rate and electrochemical consumption can be expressed as:

 (2)

Here weff is the effective width of the cathode-rib area that is electrochemically active and jcrib is the corresponding average local current density. k is the effective permeability coefficient and μ the viscosity coefficient, P the air pressure. Based on the effective parameter weff, we may write

 (3)
 (4)

It is natural to expect that jcrib increases with the increased current density in the cathode-channel region, jcc. As a first order approximation, one may simply write,

 (5)

where h is a parameter independent of jcc. Combining Eqs.(2)-(5), we have

 (6)

Eq.(6) means that the contribution of the cathode-rib region to the current generation per unit cell area is known once jcc is given. That is, it is not necessary to solve jcrib directly, eliminating the need of the described above fine mesh in the cathode-rib region. Naturally, the usefulness of Eq.(6) depends on whether a constant parameter h for different jcc of practical importance can be found.

C. Implementation of multiphysics model

Eq.(1) shows that the O2 distribution in the cathode-channel region is known once jcc is given. Combining with Eq.(6), for the mass transport and electrochemistry processes, it is necessary only to solve jcc. Therefore, the overall electrical potential balance of an SOFC stack may be expressed as:

 (7)
 (8)

Here Vstack is the total operating voltage of the stack, Vi the output voltage of the ith cell in the stack. ENernst, ηOhm, ηact and ηcon are the Nernst potential, the Ohmic polarization, activation polarization, and the concentration polarization, respectively. The subscripts i, a, c denote the ith cell, the anode and the cathode, respectively.

The expressions for computing ENernst, ηOhm, ηaact, ηcact, ηacon, and ηccon can be found in Ref.[25]. All the quantities are dependent on the local temperature and, except for ηOhm, on the gas species' concentrations. Therefore, the electrical potential balance equation, Eq.(7) and Eq.(8) are naturally coupled with the thermal conduction and the transport of gas species. The thermal conduction is in turn dependent on the heat sources of chemical and electrochemical reactions and ohmic heating as well as thermal boundary. The gas transport is also dependent on the source/sink term due to the chemical and electrochemical reactions. That is, all the physical processes are strongly coupled to each other. A complete set of multiphysics governing equations can be found in Ref.[30] and will be not repeated for brevity. Here we focus on the technique of applying Eq.(1) and (6) in computing the polarization terms of Eq.(8).

Due to Eq.(1) and (6), Eq.(8) may be confined to the block concerning the cathode-channel area of width wcc. For the computation of activation and Ohmic polarizations, the relevant local current density, jcc, is used, in accordance with their definitions. For computing ηcon for the cathode, CO2ce is required. Notice that, as both cathode-rib sides of a given cathode-channel contribute to the current generation, the overall cathode-rib contribution to the unit area current generation is 2jcribweff. Consequently, CO2ce is calculated with Eq.(1) using the effective current density, j=jcc+2jcribweff/wcc, to account for the overall electrochemical consumption of O2. Similarly, j is used for the computations of source/sink of H2O/H2 in the anode as well as the total cell operating current. The principle governing the implementation is to maintain the following physical principle of the stack operations: (ⅰ) the mass balance and the gas species distributions due to the electrochemical conversions occurring in the wcc and weff regions, (ⅱ) the current density distribution and the associated ohmic and activation polarization losses in the wcc region, (ⅲ) the total current generation and the overall electrical potential balance of each cell in the stack.

As there is no need to explicitly solve the O2 transport in the cathode, the electrolyte and cathode layers may be combined into a single material layer with properly adjusted electrical and thermal properties. This observation can be used to further simplify the mesh setting and is implemented in our stack model.

The multiphysics model is implemented through commercial CFD software of Fluent® [53] with a number of user defined functions (UDFs) to account for the specific nature of the SOFC physics. The UDFs account for the species source terms and the effective current generation reflected by Eq.(1) and Eq.(6) as well as the corresponding electrochemical reaction and heat source terms. The UDFs also include a module for the multi-component gas diffusion in the porous anode as described by an equivalent of the Dusty gas model [26] a dynamic reaction rate equation for the methane steam reforming [55]. A self-consistent iteration scheme is also coded as a UDF to determine the cell and stack outputs of the current and voltages.

Ⅲ. RESULTS AND DISCUSSION A. Validation of the stack model

To validate the use of Eq.(1) and Eq.(6) to simply the multiphysics in SOFC, a single-cell stack model with and without the analytical simplification is constructed to simulate the cell I-V curve. The material properties and geometric parameters described in Ref.[25] are used. The grid dependence test shows that the original model without the analytical approximation requires at least 105 grids for each repeating pitch unit, while only 8400 grids are required by the model using the analytical expressions of Eq.(1) and (6). The details of the mesh settings are shown in Fig. 2.

 FIG. 2 Comparison of the numerical grids for repeating pitch unit of SOFC required by a fully coupled multiphysics model without (a) and with (c) the analytical approximation developed in this work. (b) and (d) are the close-up view of (a) and (c), respectively

Figure 3 compares the I-V curves of the two models. The parameter h in Eq.(6) is set at 0.91. As shown in Fig. 3, the results obtained with the analytical approximation are in excellent agreement with the results obtained by the fully coupled numerical model for all practical operating current outputs. The agreement demonstrates that the method developed here can be used to simulate SOFC stack and obtain reliable results.

 FIG. 3 Comparison of I-V curves obtained with (red) and without (blue) the analytical approximation.

Moreover, it is important to note that the numerical instability induced by the coupled O2 transport and electrochemical consumption in the cathode-rib region is removed in the new model. Figure 4 shows the numerical convergence behavior of the new multiphysics model. As shown in Fig. 4, the residual curves are rather smooth for a multiphysics simulation of SOFC. The worst convergence of the continuity equations is about 10-4, in comparison with a typical convergence level of 10-2 shown in the model presented in Fluent® [53]. The result shows that the new model is numerically highly efficient and stable.

 FIG. 4 Examples of the residual convergent behaviors obtained by multiphysics simulation of an SOFC stack with the analytical relationship between the O2 transport and current production.
B. Multiphysics simulation of production scale SOFC stack

With the help of Eq.(1) and Eq.(6), a truly multiphysics fully coupled numerical model can be constructed for a production scale pSOFC stack and simulated with high numerical efficiency. To demonstrate the success of the new method, numerical examples are provided for the simulations of the complex multiphysics behaviors of production scale SOFC stack with a heat convection boundary condition.

Figure 5 shows a geometric model for a production scale 30-cell pSOFC stack with parallel flow design. The overall stack size is Lx×Ly×Lz=153 mm×131.5 mm×134 mm, where x is along the direction of flow in the gas channel indicated in Fig. 1 and z is along the stack thickness direction. The PEN structures and relevant mesh setting were the same as that presented in Fig. 2 (c) and (d). The computational domain includes all cells, gas channels, manifolds, interconnects, seals and frames. The total number of grids for the model is about 1.2×107. Some operating parameters and boundary conditions for the simulations are indicated in Table Ⅰ. In addition, the thermal boundary for the stack walls is of heat convection type. The film coefficient of the heat exchange is set to 15 W/(m2·K), contacting a free stream at 300 K.

 FIG. 5 (a) Geometric model of a 30-cell planar SOFC stack, (b) fuel or air flow path within the stack.
Table Ⅰ Basic operating parameters and boundary conditions for the stack simulations.

Figure 6 shows the temperature distribution of the stack. Temperature varies strongly due to the endothermic steam reforming reaction coupled with the exothermic electrochemical reaction, the water-gas shift reaction, and the Ohmic heating. The heat transport within the stack is coupled to the thermal exchange with the surrounding environment. The boundary heat convection affected the stack temperature profile in a profound way. Temperature in the cell area near the stack end plate is markedly lower than that in the middle, with a temperature difference of about 50 K. The maximum temperature difference with the whole stack is about 200 K. The maximum temperature difference within a cell is the highest for the middle cell that is least influenced by the heat exchange with the stack surrounding.

 FIG. 6 Temperature distribution in (a) the whole stack, (b) three representative cells and flow channles. The three cells are the top, middle, and bottom cells of the stack.

The reduced temperature in the cells near the stack end plates due to the boundary heat convection slowered down the endothermic steam reforming reaction therein. As a result, CH4 is reformed faster in the middle cells than that in the end cells, as can be seen in Fig. 7(a) and (b). The difference in the steam reforming reaction rate induced by the uneven temperature distribution causes different distributions of fuel gas species in different stack cells, as shown in Fig. 7(c). The uneven temperature distribution also causes different current distributions for different cells in the stack, as shown in Fig. 7(d).

 FIG. 7 Distributions of (a) CH4, (b) H2, (c) H2O in kmol/m3, and (d) current densities in 10−4A/cm2.

As the ionic conductivity of the electrolyte material is strongly temperature dependent, the different cell temperatures also result in different cell output volatges due to the constraint that every cell in the stack should generate the same amount of total current. Moreover, as the fuel flow viscosity is dependent on its gas composition and temperature, the flow distributions among different stack cells are necessarily different. Both the nonuniform flow distribution and cell output voltages are detrimental to the stack electrical efficiency and should be avoided as much as possible. Therefore, a reduced heat exchange between the stack and the environment is favorable. That is, a good insulation is not only beneficial for the stack thermal self-sustainability, but is also expected to be helpful to the stack performance.

The numerical simulations of the above 30-cell stack are carried out in a 2-CPU/16-core PC. The fully coupled multiphysics model of production scale SOFC stack developed here is highly efficient and stable numerically. Simulations with this numerical model can provide reliable information about the stack operations and has great potential for guiding the stack design and the selection of operating parameters.

Ⅳ. Conclusion

A crucial technical improvement for realizing a multiphysics fully coupled numerical model for production scale SOFC stack has been achieved. The improvement is made by developing an analytical model for solving the O2 transport in the cathode and its coupling with the electrochemistry. The model significantly reduces the grid complexity and numerical stability. The analytical model is validated for all practical working conditions by comparing with the results obtained with rigorous simulations. The new algorithm makes it feasible to construct a multiphysics model for production scale SOFC stack. Numerical examples are shown by simulating a 30-cell stack, revealing a wealthy of valuable information about the stack operation and demonstrating the numerical efficiency and stability of the model. The multiphysics model developed here can be used to speed up the development of the SOFC technology by selecting the stack design and operating parameters.

Ⅴ. Acknowledgments

This work is supported the National Natural Science Foundation of China (No.11374272 and No.11574284), the National Basic Research Program of China (No.2012CB215405) and Collaborative Innovation Center of Suzhou Nano Science and Technology are gratefully acknowledged.

 [1] B.C.H. Steele, and A. Heinzel, Nature 414 , 345 (2001). DOI:10.1038/35104620 [2] R. M. Ormerod, Chem. Soc. Rev. 32 , 17 (2003). DOI:10.1039/b105764m [3] A. Hawkes, I. Staffell, D. Brett, and N. Brandon, En-ergy Environ. Sci. 2 , 729 (2009). DOI:10.1039/b902222h [4] R.J. Gorte, and J.M. Vohs, Ann. Rev. Chem. Biomol. Eng. 2 , 9 (2011). DOI:10.1146/annurev-chembioeng-061010-114148 [5] I. Dincer, and C. Acar, Int. J. Energy Res. 39 , 585 (2015). DOI:10.1002/er.v39.5 [6] H. Yakabe, T. Ogiwara, M. Hishinuma, and I. Yasuda, J. Power Sources 102 , 144 (2001). DOI:10.1016/S0378-7753(01)00792-3 [7] E. Achenbach, J. Power Sources 49 , 333 (1994). DOI:10.1016/0378-7753(93)01833-4 [8] J. R. Ferguson, J. M. Fiard, and R. Herbin, J. Power Sources 58 , 109 (1996). DOI:10.1016/0378-7753(95)02269-4 [9] H. Yakabe, M. Hishunuma, M. Uratani, Y. Matsuzaki, and I. Yasuda, J. Power Sources 86 , 423 (2000). DOI:10.1016/S0378-7753(99)00444-9 [10] M. Iwata, T. Hikosaka, M. Morita M, T. Iwanari, K. Ito, K. Onda, Y. Esaki, Y. Sakaki, and S. Nagata, Solid State Ionics 132 , 297 (2000). DOI:10.1016/S0167-2738(00)00645-7 [11] P. V. Hendriksen, Model Studies of Internal Steam Reforming in SOFC Stacks. Solid Oxide Fuel Cells (SOFC-5) Proceedings, U. Stimming, S. C. Singhal, H. Tagawa, and W. Lehnert Eds. , Pennington, NJ: Elec-trochemical Society, Incorporated, 1319 (1997). [12] Z. J. Lin, Y. Gu, and X. H. Zhang, J. Electrochem. 8 , 445 (2002). [13] D. H. Jeon, J. H. Nam, and C. J. Kim, J. Electrochem. Soc. 153 , A406 (2006). DOI:10.1149/1.2139954 [14] D. H. Jeon, Electrochim. Acta 54 , 2727 (2009). DOI:10.1016/j.electacta.2008.11.048 [15] P.-W. Li, and M.K. Chyu, J. Power Sources 124 , 487 (2003). DOI:10.1016/j.jpowsour.2003.06.001 [16] M. Lockett, M.J. H. Simmons, and K. Kendall, J. Power Sources 131 , 243 (2004). DOI:10.1016/j.jpowsour.2003.11.082 [17] G. L. Wang, Y. Z. Yang, H. O. Zhang, and W. S. Xia, J. Power Sources 167 , 398 (2007). DOI:10.1016/j.jpowsour.2007.02.019 [18] S. X. Liu, C. Song, and Z. J. Lin, J. Power Sources 183 , 214 (2008). DOI:10.1016/j.jpowsour.2008.04.054 [19] F. Arpino, and N. Massarotti, Energy 34 , 2033 (2009). DOI:10.1016/j.energy.2008.08.025 [20] S. X. Liu, W. Kong, and Z. J. Lin, J. Power Sources 194 , 854 (2009). DOI:10.1016/j.jpowsour.2009.06.056 [21] W. S. Xia, Y. Z. Yang, and Q. S. Wang, J. Power Sources 194 , 886 (2009). DOI:10.1016/j.jpowsour.2009.06.009 [22] T. X. Ho, P. Kosinski, A. C. Hoffmann, and A. Vik, J. Power Sources 195 , 6764 (2010). DOI:10.1016/j.jpowsour.2010.03.090 [23] A. Mauro, F. Arpino, and N. Massarotti, Int. J. Hydro-gen Energy 36 , 10288 (2011). DOI:10.1016/j.ijhydene.2010.10.023 [24] H. Iwai, Y. Yamamoto, M. Saito, and H. Yoshida, En-ergy 36 , 2225 (2011). [25] W. Kong, J. Y. Li, S. X. Liu, and Z. J. Lin, J. Power Sources 204 , 106 (2012). DOI:10.1016/j.jpowsour.2012.01.041 [26] W. Kong, H. Y. Zhu, Z. Y. Fei, and Z. J. Lin, J. Power Sources 206 , 171 (2012). DOI:10.1016/j.jpowsour.2012.01.107 [27] M. Ni, Energy Convers. Manag. 70 , 116 (2013). DOI:10.1016/j.enconman.2013.02.008 [28] X. Han, D. Zheng, and B. F. Bai, Energy 67 , 575 (2014). DOI:10.1016/j.energy.2014.02.021 [29] H. R. Amedi, B. Bazooyar, and M. R. Pishvaie, Energy 90 , 605 (2015). DOI:10.1016/j.energy.2015.07.095 [30] B. X. Wang, J. Zhu, and Z. Lin, Appl. Energy 176 , 1 (2016). DOI:10.1016/j.apenergy.2016.05.049 [31] K. P. Recknagle, R. E. Williford, L. A. Chick, and M. A. Khaleel, J. Power Sources 113 , 109 (2003). DOI:10.1016/S0378-7753(02)00487-1 [32] N. Autissier, D. Larrain, J. Van herle, and D. Favrat, J. Power Sources 131 , 313 (2004). DOI:10.1016/j.jpowsour.2003.11.089 [33] R. T. Leah, N. P. Brandon, and P. Aguiar, J. Power Sources 145 , 336 (2005). DOI:10.1016/j.jpowsour.2004.12.067 [34] C. M. Huang, S. S. Shy, and C. H. Lee, J. Power Sources 183 , 205 (2008). DOI:10.1016/j.jpowsour.2008.04.059 [35] A. A. Kulikovsky, J. Fuel Cell Sci. Technol. 7 , 011015 (2010). DOI:10.1115/1.3119060 [36] S. Hosseini, K. Ahmed, and M. O. Tad, J. Power Sources 234 , 180 (2013). DOI:10.1016/j.jpowsour.2012.12.123 [37] B. Lin, Y. X. Shi, M. Ni, and N. S. Cai, Int. J. Hydrogen Energy 40 , 3035 (2015). DOI:10.1016/j.ijhydene.2014.12.088 [38] M. Fardadi, D.F. McLarty, and F. Jabbari, Appl. En-ergy 178 , 43 (2016). DOI:10.1016/j.apenergy.2016.06.015 [39] W. X. Bi, D. F. Chen, and Z. J. Lin, Int. J. Hydrogen Energy 34 , 3873 (2009). DOI:10.1016/j.ijhydene.2009.02.071 [40] M. Peksen, Int. J. Hydrogen Energy 36 , 11914 (2011). DOI:10.1016/j.ijhydene.2011.06.045 [41] M. Peksen, Int. J. Hydrogen Energy 39 , 5137 (2014). DOI:10.1016/j.ijhydene.2014.01.063 [42] S. S. Wei, T. H. Wang, and J. S. Wu, Energy 69 , 553 (2014). DOI:10.1016/j.energy.2014.03.052 [43] A. Al-Masri, M. Peksen, L. Blum, and D. Stolten, Appl. Energy 135 , 539 (2014). DOI:10.1016/j.apenergy.2014.08.052 [44] L. Petruzzi, S. Cocchi, and F. Fineschi, J. Power Sources 118 , 96 (2003). DOI:10.1016/S0378-7753(03)00067-3 [45] A. C. Burt, I. B. Celik, R. S. Gemmen, and A. V. Smirnov, J. Power Sources 126 , 76 (2004). DOI:10.1016/j.jpowsour.2003.08.034 [46] M. A. Khaleel, Z. Lin, P. Singh, W. Surdoval, and D. Collin, J. Power Sources 130 , 136 (2004). DOI:10.1016/j.jpowsour.2003.11.074 [47] B.A. Haberman, and J.B. Young, J. Fuel Cell Sci. Technol. 5 , 011006 (2008). DOI:10.1115/1.2786468 [48] H. Mounir, A. El Gharad, M. Belaiche, and M Boukalouch, Energy Convers. Manag. 50 , 2685 (2009). DOI:10.1016/j.enconman.2009.06.023 [49] C. K. Lin, L. H. Huang, L. K. Chiang, and Y. P. Chyou, J. Power Sources 192 , 515 (2009). DOI:10.1016/j.jpowsour.2009.03.010 [50] S.F. Lee, and C.W. Hong, Int. J. Hydrogen Energy 35 , 1330 (2010). DOI:10.1016/j.ijhydene.2009.11.095 [51] K. Sudaprasert, R. P. Travis, and R.F. Martinez-Botas, J. Fuel Cell Sci. Technol. 7 , 011002 (2010). DOI:10.1115/1.3080810 [52] K. Lai, B. J. Koeppel, K. S. Choi, K. P. Recknagle, X. Sun, L. A. Chick, V. Korolev, and M. Khaleel, J. Power Sources 196 , 3204 (2011). DOI:10.1016/j.jpowsour.2010.11.123 [53] ANSYS, ANSYS FLUENT 14. 5 in ANSY S Workbench Users Guide, Canonsburg, PA: ANSYS, Inc. , (2012). [54] A. Li, X. Fang, and Z. Lin, ECS Trans. 68 , 3025 (2015). DOI:10.1149/06801.3025ecst [55] B.X. Wang, J. Zhu, and Z. J. Lin, Chin. J. Chem. Phys. 28 , 299 (2015). DOI:10.1063/1674-0068/28/cjcp1503033