The article information
 Ang Li, Zijing 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): 139146
 化学物理学报, 2017, 30(2): 139146
 http://dx.doi.org/10.1063/16740068/30/cjcp1610198

Article history
 Received on: October 23, 2016
 Accepted on: November 5, 2016
Solid oxide fuel cell (SOFC) has considerable potential in future clean energy industry due to its high efficiency and low level of pollutant emissions [15]. 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 triplephaseboundary, 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 indepth 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 inhouse developed software employing simplified SOFC geometries and coarse grids [614]. 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 [1552]. 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 [1739], even though some may be claimed as stack models [31, 33, 3537]. 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 [4052]. 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 5and 10cell 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 cm^{2} 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 indepth 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 modelUnderstanding 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 gaschannelinterconnectrib pitches consisting of the cathode (positive)electrolyteanode (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 O_{2} transport and electrochemistry in the interconnect rib covered cathode area, as schematically illustrated in Fig. 1. Due to the relatively low O_{2} diffusion coefficient and the thinness of cathode layer (50200 μm) for O_{2} flux, the O_{2} concentration decreases rapidly away from the gas channel, as shown in Fig. 1(c). To describe the O_{2} transport in the cathoderib area requires a very fine mesh setting, as illustrated in Fig. 1(b) and (d). Moreover, the strong coupling between the rapidly varying O_{2} concentration and the electrochemistry may cause instability in the numerical simulation.
To solve the O_{2} 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 10^{5} grid points are necessary. A 30cell 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 singlecell stack models or stacks with a low number of cells and small cell sizes.
B. Determination of cathode O_{2} transport and current productionAs discussed above, a fine grid setting is required due to the need to describe the O_{2} transport and electrochemistry in the cathoderib area. However, the contribution of the cathoderib area to the physical processes (O_{2} 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 cathoderib area with reasonable accuracy and eliminate the associated fine mesh requirement. Such a method is derived below.
Denoting the O_{2} concentration at the cathodeelectrolyte interface, as C_{O2}^{ce} and at the cathodechannel interface, as C_{O2}^{cc}. As the cathode is very thin, the O_{2} flux along the cathode thickness direction is mainly driven by O_{2} 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, D_{O2}^{eff} the effective O_{2} diffusion coefficient in the cathode, j the local current density, F the Faraday constant.
The equality of O_{2} flux and electrochemical consumption is also valid for the cathoderib region. As C_{O2}^{ce} decays rapidly inside the cathoderib 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 O_{2} molar flow rate and electrochemical consumption can be expressed as:
(2) 
Here w_{eff} is the effective width of the cathoderib area that is electrochemically active and j_{crib} 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 w_{eff}, we may write
(3) 
(4) 
It is natural to expect that j_{crib} increases with the increased current density in the cathodechannel region, j_{cc}. As a first order approximation, one may simply write,
(5) 
where h is a parameter independent of j_{cc}. Combining Eqs.(2)(5), we have
(6) 
Eq.(6) means that the contribution of the cathoderib region to the current generation per unit cell area is known once j_{cc} is given. That is, it is not necessary to solve j_{crib} directly, eliminating the need of the described above fine mesh in the cathoderib region. Naturally, the usefulness of Eq.(6) depends on whether a constant parameter h for different j_{cc} of practical importance can be found.
C. Implementation of multiphysics modelEq.(1) shows that the O_{2} distribution in the cathodechannel region is known once j_{cc} is given. Combining with Eq.(6), for the mass transport and electrochemistry processes, it is necessary only to solve j_{cc}. Therefore, the overall electrical potential balance of an SOFC stack may be expressed as:
(7) 
(8) 
Here V_{stack} is the total operating voltage of the stack, V_{i} the output voltage of the ith cell in the stack. E^{Nernst}, η^{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 E^{Nernst}, η^{Ohm}, η_{a}^{act}, η_{c}^{act}, η_{a}^{con}, and η_{c}^{con} 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 cathodechannel area of width w_{cc}. For the computation of activation and Ohmic polarizations, the relevant local current density, j_{cc}, is used, in accordance with their definitions. For computing η^{con} for the cathode, C_{O2}^{ce} is required. Notice that, as both cathoderib sides of a given cathodechannel contribute to the current generation, the overall cathoderib contribution to the unit area current generation is 2j_{crib}w_{eff}. Consequently, C_{O2}^{ce} is calculated with Eq.(1) using the effective current density, j=j_{cc}+2j_{crib}w_{eff}/w_{cc}, to account for the overall electrochemical consumption of O_{2}. Similarly, j is used for the computations of source/sink of H_{2}O/H_{2} 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 w_{cc} and w_{eff} regions, (ⅱ) the current density distribution and the associated ohmic and activation polarization losses in the w_{cc} 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 O_{2} 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 multicomponent 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 selfconsistent 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 modelTo validate the use of Eq.(1) and Eq.(6) to simply the multiphysics in SOFC, a singlecell stack model with and without the analytical simplification is constructed to simulate the cell IV 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 10^{5} 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.
Figure 3 compares the IV 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.
Moreover, it is important to note that the numerical instability induced by the coupled O_{2} transport and electrochemical consumption in the cathoderib 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.
B. Multiphysics simulation of production scale SOFC stackWith 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 30cell pSOFC stack with parallel flow design. The overall stack size is L_{x}×L_{y}×L_{z}=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×10^{7}. 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/(m^{2}·K), contacting a free stream at 300 K.
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 watergas 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.
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, CH_{4} 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).
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 selfsustainability, but is also expected to be helpful to the stack performance.
The numerical simulations of the above 30cell stack are carried out in a 2CPU/16core 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.
Ⅳ. ConclusionA 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 O_{2} 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 30cell 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.
Ⅴ. AcknowledgmentsThis 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, Energy 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/annurevchembioeng061010114148 
[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/S03787753(01)007923 
[7]  E. Achenbach, J. Power Sources 49 , 333 (1994). DOI:10.1016/03787753(93)018334 
[8]  J. R. Ferguson, J. M. Fiard, and R. Herbin, J. Power Sources 58 , 109 (1996). DOI:10.1016/03787753(95)022694 
[9]  H. Yakabe, M. Hishunuma, M. Uratani, Y. Matsuzaki, and I. Yasuda, J. Power Sources 86 , 423 (2000). DOI:10.1016/S03787753(99)004449 
[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/S01672738(00)006457 
[11]  P. V. Hendriksen, Model Studies of Internal Steam Reforming in SOFC Stacks. Solid Oxide Fuel Cells (SOFC5) Proceedings, U. Stimming, S. C. Singhal, H. Tagawa, and W. Lehnert Eds. , Pennington, NJ: Electrochemical 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. Hydrogen Energy 36 , 10288 (2011). DOI:10.1016/j.ijhydene.2010.10.023 
[24]  H. Iwai, Y. Yamamoto, M. Saito, and H. Yoshida, Energy 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/S03787753(02)004871 
[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. Energy 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. AlMasri, 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/S03787753(03)000673 
[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. MartinezBotas, 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/16740068/28/cjcp1503033 