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

Ang Li,Zi-jing Lin

Hefei National Lab oratory for Physical Sciences at the Microscales&CAS Key Lab oratory of Strongly-Coupled Quantum Matter Physics,Department of Physics,University of Science and Technology of China,Hefei 230026,China

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

Ang Li,Zi-jing Lin?

Hefei National Lab oratory for Physical Sciences at the Microscales&CAS Key Lab oratory of Strongly-Coupled Quantum Matter Physics,Department of Physics,University of Science and Technology of China,Hefei 230026,China

(Dated:Received on October 23,2016;Accepted on November 05,2016)

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 O2transport and the associated electrochemistry.To overcome the technical obstacle,an analytical model for solving the O2transport 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 briefly discussed.The multiphysics modeling tool can be used to guide the stack design and select the operating parameters.

Simulation,Mesh setting,Analytical model,Computational efficiency,Numerical stability

I.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-phaseboundary,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 onstacks 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 cm2that 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.

II.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 gaschannel-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 O2transport and electrochemistry in the interconnect rib covered cathode area,as schematically illustrated in Fig.1.Due to the relatively low O2diffusion coefficient and the thinness of cathode layer(50-200μm)for O2flux,the O2concentration decreases rapidly away from the gas channel, as shown in Fig.1(c).To describe the O2transport 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 O2concentration and the electrochemistry may cause instability in the numerical simulation.

To solve the O2transport 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 105grid 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 O2transport and current production

As discussed above,a fine grid setting is required due to the need to describe the O2transport and electrochemistry in the cathode-rib area.However,the contribution of the cathode-rib area to the physical processes (O2distribution,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.

FIG.1 Representative structural domain of SOFC stack and the associated O2distribution 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)O2distribution in the cathode-rib region,and(d)mesh setting for the cathode-rib domain.

Here weffis the effective width of the cathode-rib area that is electrochemically active and jcribis 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

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

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

Eq.(6)means that the contribution of the cathoderib region to the current generation per unit cell area is known once jccis given.That is,it is not necessary to solve jcribdirectly,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 jccof practical importance can be found.

C.Implementation of multiphysics model

Eq.(1)shows that the O2distribution in the cathodechannel region is known once jccis 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:

Here Vstackis the total operating voltage of the stack,Vithe output voltage of the i th cell in the stack.ENernst, ηOhm,ηactand ηconare the Nernst potential,the Ohmic polarization,activation polarization,and the concentration polarization,respectively.The subscripts i,a,c denote the i th cell,the anode and the cathode,respectively.

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.

As there is no need to explicitly solve the O2transport 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.

III.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 105grids 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 themesh settings are shown in Fig.2.

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.

Moreover,it is important to note that the numerical instability induced by the coupled O2transport 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-2shown in the model presented in Fluent?[53].The result shows that the new model is numerically highly efficient and stable.

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

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 I.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.4 Examples of the residual convergent behaviors obtained by multiphysics simulation of an SOFC stack with the analytical relationship between the O2transport and current production.

FIG.5(a)Geometric model of a 30-cell planar SOFC stack, (b)fuel or air flow path within the stack.

Figure 6 shows the temperature distribution of thestack.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.

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

TABLE I Basic operating parameters and boundary conditions for the stack simulations.

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.

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

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

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

[2]R.M.Ormerod,Chem.Soc.Rev.32,17(2003).

[3]A.Hawkes,I.Staffell,D.Brett,and N.Brandon,Energy Environ.Sci.2,729(2009).

[4]R.J.Gorte and J.M.Vohs,Ann.Rev.Chem.Biomol. Eng.2,9(2011).

[5]I.Dincer and C.Acar,Int.J.Energy Res.39,585 (2015).

[6]H.Yakabe,T.Ogiwara,M.Hishinuma,and I.Yasuda, J.Power Sources 102,144(2001).

[7]E.Achenbach,J.Power Sources 49,333(1994).

[8]J.R.Ferguson,J.M.Fiard,and R.Herbin,J.Power Sources 58,109(1996).

[9]H.Yakabe,M.Hishunuma,M.Uratani,Y.Matsuzaki, and I.Yasuda,J.Power Sources 86,423(2000).

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

[11]P.V.Hendriksen,M odel S tudies of I nternal S team Ref or ming in S OF C S tack s.S olid Oxide F uel C ells (S OF C-5)P roceeding s,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).

[14]D.H.Jeon,Electrochim.Acta 54,2727(2009).

[15]P.-W.Li and M.K.Chyu,J.Power Sources 124,487 (2003).

[16]M.Lockett,M.J.H.Simmons,and K.Kendall,J. Power Sources 131,243(2004).

[17]G.L.Wang,Y.Z.Yang,H.O.Zhang,and W.S.Xia, J.Power Sources 167,398(2007).

[18]S.X.Liu,C.Song,and Z.J.Lin,J.Power Sources 183, 214(2008).

[19]F.Arpino and N.Massarotti,Energy 34,2033(2009).

[20]S.X.Liu,W.Kong,and Z.J.Lin,J.Power Sources 194,854(2009).

[21]W.S.Xia,Y.Z.Yang,and Q.S.Wang,J.Power Sources 194,886(2009).

[22]T.X.Ho,P.Kosinski,A.C.Hoffmann,and A.Vik,J. Power Sources 195,6764(2010).

[23]A.Mauro,F.Arpino,and N.Massarotti,Int.J.Hydrogen Energy 36,10288(2011).

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

[26]W.Kong,H.Y.Zhu,Z.Y.Fei,and Z.J.Lin,J.Power Sources 206,171(2012).

[27]M.Ni,Energy Convers.Manag.70,116(2013).

[28]X.Han,D.Zheng,and B.F.Bai,Energy 67,575 (2014).

[29]H.R.Amedi,B.Bazooyar,and M.R.Pishvaie,Energy 90,605(2015).

[30]B.X.Wang,J.Zhu,and Z.Lin,Appl.Energy 176,1 (2016).

[31]K.P.Recknagle,R.E.Williford,L.A.Chick,and M. A.Khaleel,J.Power Sources 113,109(2003).

[32]N.Autissier,D.Larrain,J.Van herle,and D.Favrat, J.Power Sources 131,313(2004).

[33]R.T.Leah,N.P.Brandon,and P.Aguiar,J.Power Sources 145,336(2005).

[34]C.M.Huang,S.S.Shy,and C.H.Lee,J.Power Sources 183,205(2008).

[35]A.A.Kulikovsky,J.Fuel Cell Sci.Technol.7,011015 (2010).

[36]S.Hosseini,K.Ahmed,and M.O.Tad,J.Power Sources 234,180(2013).

[37]B.Lin,Y.X.Shi,M.Ni,and N.S.Cai,Int.J.Hydrogen Energy 40,3035(2015).

[38]M.Fardadi,D.F.McLarty,and F.Jabbari,Appl.Energy 178,43(2016).

[39]W.X.Bi,D.F.Chen,and Z.J.Lin,Int.J.Hydrogen Energy 34,3873(2009).

[40]M.Peksen,Int.J.Hydrogen Energy 36,11914(2011).

[41]M.Peksen,Int.J.Hydrogen Energy 39,5137(2014).

[42]S.S.Wei,T.H.Wang,and J.S.Wu,Energy 69,553 (2014).

[43]A.Al-Masri,M.Peksen,L.Blum,and D.Stolten,Appl. Energy 135,539(2014).

[44]L.Petruzzi,S.Cocchi,and F.Fineschi,J.Power Sources 118,96(2003).

[45]A.C.Burt,I.B.Celik,R.S.Gemmen,and A.V. Smirnov,J.Power Sources 126,76(2004).

[46]M.A.Khaleel,Z.Lin,P.Singh,W.Surdoval,and D. Collin,J.Power Sources 130,136(2004).

[47]B.A.Haberman and J.B.Young,J.Fuel Cell Sci. Technol.5,011006(2008).

[48]H.Mounir,A.El Gharad,M.Belaiche,and M Boukalouch,Energy Convers.Manag.50,2685(2009).

[49]C.K.Lin,L.H.Huang,L.K.Chiang,and Y.P.Chyou, J.Power Sources 192,515(2009).

[50]S.F.Lee and C.W.Hong,Int.J.Hydrogen Energy 35, 1330(2010).

[51]K.Sudaprasert,R.P.Travis,and R.F.Martinez-Botas, J.Fuel Cell Sci.Technol.7,011002(2010).

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

[53]ANSYS,AN S Y S F LU E N T 14.5 in AN S Y S W ork bench U sers Guide,Canonsburg,PA:ANSYS, Inc.,(2012).

[54]A.Li,X.Fang,and Z.Lin,ECS Trans.68,3025(2015).

[55]B.X.Wang,J.Zhu,and Z.J.Lin,Chin.J.Chem.Phys. 28,299(2015).

?Author to whom correspondence should be addressed.E-mail: zjlin@ustc.edu.cn,Tel:+86-551-63606345