Fluid dynamic study of fuel cell devices--simulation and experimental validation

ELSEVIER Journal of Power Sources 52 (1994) 243-249 Fluid dynamic study of fuel cell devices: simulation and experimental validation P. Costamagna a,*, E. Arato a, E. Achenbach b, U. Reus b a Istitito di Scienze e Tecnotogie dell’lngegneria Chimica, Facoltd di Ingegneria dell’lJniversitci di Geneva, V?a Opera Pia 15, 16145 Genoa, Itab b Institut fiir Energieverfahrenstechnik (IEV), Forschungszentmm Jiilich (KFA), 52425 Jiilich, Gemany Received 8 July 1994; accepted 20 August 1994 Abstract The paper is concerned with the mass flow distribution in fuel cell stacks. In particular, the flow through the manifold system connected to the parallel arrangement of the cell channels is modelled and numerically treated. The numerical results are recognized to be more realistic than those obtained by means of an approximate analytical solution since more detailed effects could be accounted for. This evidence is confirmed by experiments carried out at a stack model device consisting of 100 cells. Pressure and velocity distributions were measured for various Reynolds numbers and geometrical shapes of the manifolds. The agreement between the experimental and numerical results is good. Keywordsr Fuel cell stacks; Fluid dynamic study 1. Introduction Steady-state and dynamic simulations of fuel cell stacks have recently been accomplished. Details about these subjects can be found in Refs. [l-18]. Particular problems which may lead to undesirable or dangerous operating conditions arise from the flow distribution in such devices [19-211. In the present contribution, a fluid dynamic computational model is presented, and the reliability of the theoretical pre- dictions is demonstrated by means of experimental data. Fig. 1 shows a sketch of a planar fuel cell stack [22]. The monolithic electrochemical reactor is composed of a number of single cells mounted one upon the other I t I ” inlet flow outlet llow Fig. 1. Drawing of a fuel cell stack with inlet and outlet manifolds * Corresponding author. and electrically connected in series. The reactants are fed to and discharged from each of the cells through feeding and exhaust manifolds. The separator or bipolar plates contain the cathodic and anodic channels which carry the oxidant (oxygen) and fuel (hydrogen), re- spectively. In such fuel cell stacks differences in the voltage output have been observed between the top and the bottom of the device. They seem to be caused by a non-uniform distribution of the feeding gas along the cell stack. In addition, a bad distribution across the channels of a cell may occur. As a consequence of these two phenomena, a lack of reactants in some portion of the cell can arise, which usually causes confined but irreversible damage. In phosphoric acid fuel cells (PAFC), for instance, electrochemical reac- tions may destroy the electrode graphite. Moreover, other severe problems may result. Since the cells are connected in series the maximum electrical current is dominated by the minimum flow rate supplied to any cell. In particular under strong operating con- ditions, recirculation arises, i.e., some of the cells get their feedstock from the outlet manifold and discharge gases into the inlet duct. In this case, the stack does not supply any net electric current. Therefore it is necessary to predict and then avoid those operating conditions where considerable non- uniformities or even recirculation occur. Thus, simu- 0378-7753/94/$07.00 0 1994 Elsevier Science S.A. All rights resewed SSDI 0378-7753(94)02014-T 244 P. Costamagna et ai. / Journal of Power Sources 52 (1994) 243-249 lation is an important tool for the design of a fuel cell stack. A simulation model for the prediction of the pressure and flow distribution in a fuel cell system has been described in the above-mentioned papers [19-211. The aim of the present work is to verify the calculations made for a special case in order to demonstrate the reliability with regard to a general application. For this purpose, an experimental facility simulating a fuel cell stack has been provided. Since the subject of this study is related to fluid dynamic aspects rather than to electrochemical ones, the apparatus was designed for operation with air under ambient conditions and no reaction takes place in it. 2. The fluid dynamic simulation model 2.1. General remarks The mass flow rate in a channel of a stack depends on the difference between the values of static pressure prevailing in the particular positions of the inlet and outlet manifold. Since the static pressure in both mani- folds varies along the streaming length due to gravi- tational, viscous and inertial effects, the pressure dif- ference acting on each of the channels will be locally different causing a non-uniform flow distribution. 2.2. Equation of motion for the channels For the flow through the channels it is assumed that the inertial effects are negligibly low compared to the friction forces. Therefore, the pressure drop of the rectangular channels under laminar flow conditions can be determined by the well-known relationship [23]: &friction =K (1) where K is a shape parameter dependent on the ratio a/b of both sides of the rectangular cross section of channels. As this equation only holds for the fully established channel flow, the inlet and outlet losses are accounted for by 1.0 and 0.5 times the dynamic pressure [23], respectively. Thus, the flow rate can be expressed as a function of the pressure difference across the cell: where: (2) c=” p ’ 4 (abn)2 2.3. Equation of motion for the manifolds If we consider a duct with lateral inlet and outlet flow rates, as shown in Fig. 2, the projection in the z- direction of the momentum balance applied to the control volume formed by a portion of the pipe between sections 1 and 2, as shown in Fig. 2, has the following form: &PSV; + ~Qout~ut +pzS + ~2s =N,pSv:+pQi,V,,+plS+pgz,S+~F (3) where we have assumed the density of the gas and the area of the cross section to be constant; CF is the sum of the components of the external forces exerted by the duct walls on the fluid, while N is a coefficient which accounts for the non-uniform velocity distribution across the duct cross section. The terms pQout~Uout and pQi,,uin stand for the momentum components of lateral effluxes and influxes, respectively. In some reports [24,2.5], it is suggested on the basis of experimental tests that lateral inlet flow rates have velocities exactly perpendicular to the manifold one, so their contribution to the momentum balance equation is zero. On the other hand, when considering small lateral effluxes, the same authors assume that thevelocity at the outlet of the control volume has a component along the main duct direction which is equal to the velocity of the main flow. Considering outflows and inflows as continuously distributed and expressing the XF term - which rep- resents the friction forces - as a function of the friction factor A, the equation can be written in the following differential form (the z-axis is assumed vertical and upward): dP - =_-2Npl:g -N;vq,,,,-pg- ;;P”.~ dz (4) h C,=K& Fig. 2. Sketch of a duct with inlet and outlet lateral flows: a volume for the momentum balance is indicated. control I? Costamagna et al. 1 Journal of Power Sources 52 (1994) 243-249 245 Substituting the continuity equation & _ qin-40ul _p dz s (5) in Eq. (4) yields: (6) The friction factor can be expressed by the formulae for the established flow through ducts. The Reynolds number varies in a wide range from laminar to the turbulent flow regime. However, the problems associated with the transition phenomena have not been taken into account. The following equations have been applied: (i) laminar flow: 21 z (7) (ii) turbulent flow: A = 0.0868Re;Ji4 (8) where Re = 4pvRhIp; Remod =fiRe and R, is the hydraulic radius defined as the ratio of the cross section of the pipe to its wetted perimeter; fi, as well as K’, depend on the geometrical dimension of the duct [23]. Eq. (4) has been used to model the exhaust manifold, where qoUt = 0, while Eq. (6) has been used to model the feed manifold, where qin = 0. These two expressions become identical except for the inertial term, which, in the exhaust duct, is twice that in the feeding one. The gravitational terms in the balance equations can be neglected under the operating conditions of the actual experimental device. An a priori estimate of the inertial and viscous effects on the pressure distribution in the manifolds shows that the most important is the first one. In addition to the above treatment, an alternative approach to the problem can be made by using handbook correlations [26] for the various contributions of the pressure losses. The basic equation is the integrated form of the Euler’s equation: (9) where AffH,,,,,, means the sum of both the viscous and the momentum loss terms. The friction losses can be expressed by means of the friction factor A defined in Eqs. (7) and (8). The losses caused by the suction and bleeding phe- nomena are modelled by applying the empirical equa- tions of branching tubes. For this purpose, the inter- sections consisting of the manifold and each of the channels have been considered as to be a series of T- crossing rectangular ducts with extreme ratios of the main to the secondary cross section area [26]. Partic- ularly the following relationships were applied: (i) for the feeding manifold (numbers 1 and 2 refer to the sections before and behind the branching, re- spectively): (ii) for the outlet manifold: A&,,,,,=[l.55(1- z)-(l- ::)3&: (11) The calculations performed on the basis of both methods came out with the same computational results. 2.4. Approximated analytical solution Since the present fluid dynamic model should be a subroutine of the main code predicting the operating conditions ‘of a stack, it could be useful to have an analytical simplified expression for the solution of the flow problem. For this purpose, the equation of motion was integrated neglecting the friction terms in the manifolds and the inlet and outlet losses in the channels. In this way, the problem can be expressed introducing the dimensionless variables cp= v/v*, ~=p -p*&w*“, {=z/h, and the parameter A=N/ShC,, by the set of the following dimensionless differential equations: dr’ d2 ’ cp - =_ - d5 c&Z2 d+ -= * d[ -2 d12 W - =A(& _ #) dJ !!!!f =A(& _ +‘) dC p’ = (+y The boundary conditions are: (12) l=I cp’=rp”=I r=O The solution of Eq. (12) is: (13) rr’ = 1 -[B tan(ABQ12 +‘=l-B2-2[B tan(ABn]’ cp’ = cp” = B tan(ABn (14) where B is given by: B tan(AB)=l. The results are evaluated and discussed below in con- nection with the experimental data. 246 I? Costamagna et al. I Journal of Power Sources 52 (1994) 243-249 3. Experimental Fig. 3 shows a sketch of the experimental apparatus. It is a plastic device having the overall dimensions 276 mm x 393 mm X 1100 mm. The cell is simulated by a 276 mm x 200 mm X 11 mm plate crossed by 10 channels, each of them 200 mm long and 20 mmX4 mm by cross section. The stack consists of 100 of these plates to which the inlet and outlet manifolds are connected. Their cross section has 245 mm x 80 mm internal di- mensions. A blower passes air through the system. The flow enters at the top of the inlet manifold and then moves downward while simultaneously being distributed to the cells. The air is collected in the outlet manifold streaming upwards and leaving it at the top. The mass flow rate is measured by means of an orifice mounted in the feeding duct. The pressure distributions along the inlet and outlet manifolds are used to carry out the comparison between the theoretical and experimental results. For this pur- pose, three rows of 23 holes each are drilled into the external walls of the manifolds. Each pressure tape is connected via a scanivalve system to a pressure trans- ducer adapted with respect to its range to the pressure INLET FLOW r-7 OUTLET FLOW \ I !-: Fig. 3. Drawing of the experimental test facility. / / level expected (10 mbar). The accuracy is given by 0.5% of the total range. Each of the pressure data results from the mean value of five subsequent measurements to account for the effect of fluctuating pressure. The experimental error is about 2%. Usually the ambient pressure was used as reference. Thus, negative numbers will occur in the diagrams, particularly when the air was sucked through the model. In preliminary tests, the pressure distribution mea- sured along the central row of taps was compared with that taken from the lateral rows to confirm the as- sumption of constant pressure across the cross section of the manifold. Having ascertained that this is true the subsequent measurements were performed using only the central pressure taps. The channels of the top fourteen cells were closed and a grid placed at the entrance of the feed manifold. Thus, effects of entrance disturbances on the mea- surement should be excluded. 4. Results In Fig. 4 the experimental data of the pressure along the manifolds are presented togetherwith the theoretical results. The figure ‘0’ denotes the top, ‘100’ the bottom plate of the stack. The agreement between theory and experiment is quite good everywhere except at the top of the stack where flow separation occurs near the leading edge. As mentioned above, the top fourteen cells were closed therefore forming an entrance section. Since this entrance effect is not modelled, the agreement between theory and experiment cannot be satisfactory. The increase of pressure down along the inlet manifold caused by the deceleration of the flow demonstrates that the inertial forces play a predominant role. This is also the reason for the strong pressure drop close to the exit of the outlet manifold. Additional experiments have been performed varying the flow rate in the range Fig. 4. Pressure distribution along feeding and exhaust manifolds: (*) experimental data, and (-) numerically simulated results. P. Coskzmagna et al. / Journal of Power Sources 52 (1994) 243-249 247 from 0.045 to 0.153 kg/s. The same agreement of the results was observed. Fig. 5 shows the same experimental results compared with those from the approximated analytical integration (Eq. (14)). As far as the above-mentioned simplifying assumptions are acceptable for a particular situation, the result is satisfactory in the sense of a first ap- proximation. Higher precision can only be achieved by the numerical solution of the differential equation sys- tem. With the purpose of obtaining a better understanding of the flow mechanisms in the manifolds these have been studied separately. Therefore, experiments have been conducted setting the boundary conditionsp = con- stant at the inlet/outlet of channels by removing one of the manifolds. Thus the interaction of both manifolds could be avoided. Fig. 6 exhibits the experimental and theoretical results for the inlet manifold keeping the pressure constant at the outlet: the experimental data show a certain scatter due to the instabilities of the decelerated flow. The simulated results exhibit an acceptable agreement Measurement Point Fig. 5. Comparison behveen pressure distributions obtained from (*) the experimental data, and (-) the analytically simulated results. -0 05 Lo.___ I 10 20 30 40 50 60 70 80 90 LOO Measurement Point Fig. 6. Pressure distribution along the inlet manifold: ( * ) experimental data; (-) simulation results, and (- - -) analytical solution. with experimental data: the error, which is calculated with reference to the total pressure difference along the manifold, is of about 4%. Thus, it can be concluded that the flow in the inlet manifold is sufficiently well simulated. The approximated solution for the inlet manifold only (the treatment is analogous to that one reported in Section 2.4) shows a pressure varying with the hyperbolic tangent of the streaming length. This result does not agree as well with the experimental data. The departures increase from 8 to 15% with increasing flow rate, since the contribution of the neglected terms becomes gradually relevant. In addition, a multiplicity of steady-state solutions has been analytically predicted. It was noticed that the solutions were not unique: there exist two different mass flow rates for a given pressure drop across the whole stack. The reason is that - due to the neglected quantities of the viscous terms in the manifold and of the entrance and outlet losses in the channels - the inertial terms of the manifold can compensate for the viscous terms in channels with increasing mass flow. The numerical model does not predict such an effect, so it can be concluded that the hypothesis of disregarding inlet and outlet channel losses and viscous terms in manifolds in this case is a quite approximated one, although both these terms are less significant than the other ones. The outlet manifold was tested by removing the inlet system. Furthermore, it was necessary to suck the air through the device. As shown in Fig. 7, the numerical and experimental results collapse within a scatter of only 1%. Again, the approximated solution, which in this case depends on the trigonometric tangent of the streamwise coordinate, yields a minor agreement ex- pressed by departures of 8%. Finally, some results of the flow distribution in the channels along the manifolds are shown in Fig. 8. The curves represent the experimental and computed data. 1 10 20 30 40 50 60 70 80 90 100 Measurement Point Fig. 7. Pressure distribution along the outlet manifold: (*) exper- imental data; (-) simulated results, and (- - -) analytical solution. 248 P. Costamagna et al. J Journal of Power Sources 52 (1994) 243-249 0.5 t 1 measurement po,nt Fig. 8. Flow rate distribution along the stack: (* ) experimental data, and (-) simulation results. The measurement technique applied is the hot wire anemometry. The mean velocity of each channel was determined by placing the probe immediately behind the channel exit. With increasing approach to the manifold exit the hot wire signal is affected by the continuously increasing velocity in the outlet manifold. There was made an attempt to correct the results for this effect. By this way a satisfactory agreement between experiment and theory can be achieved. 5. Conclusions A mathematical model is presented for the simulation of the flow through fuel cell stacks. The momentum equations established for the manifold system and the cell channels are numerically solved. The theoretical results are verified by experimental data obtained from a special experimental device consisting of 100 bipolar plates and adapted to fluid dynamic investigations. The agreement between theory and experiment is good. It is obvious that in the Reynolds number range covered by the actual experiments the inertial terms play a considerable role with respect to the pressure distri- bution along the manifold system. An alternative method of solving the flow problem was investigated. Particularly, the model of subsequent branching tubes was applied. The agreement between the results of both methods was very good. Finally, the suitability of an analytic solution which required some simplifications was checked. It was seen that remarkable departures from the experimental find- ings occurred. 6. List of symbols A a, b dimensionless coefficient (-) long and short sides of the channel section, respectively (m) parame