- Email: [email protected]
Mass transfer formulation for polymer electrolyte membrane fuel cell cathode
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 1 1 6 4 1 e1 1 6 5 0
Available online at www.sciencedirect.com
ScienceDirect journal homepage: www.elsevier.com/locate/he
Mass transfer formulation for polymer electrolyte membrane fuel cell cathode Steven B. Beale Institute of Energy and Climate Research, IEK-3, Forschungszentrum Ju¨lich GmbH, 92425 Ju¨lich, Germany
article info
abstract
Article history:
A standard mechanical engineering mass transfer analysis is adapted to obtain perfor-
Received 3 February 2015
mance calculations for the cathode of a typical polymer electrolyte membrane fuel cell.
Received in revised form
Mathematical details are provided and sample calculations performed, based on a number
4 May 2015
of simplifying assumptions. The relative advantages of low and high-mass transfer for-
Accepted 8 May 2015
mulations are critically appraised and compared to a highly-simplified method.
Available online 22 July 2015
Copyright © 2015, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved.
Keywords: Fuel cells Transport phenomena Mass transfer Concentration polarisation Diffusion over-potential
Introduction It is popular within some elements of the fuel cell modelling community to consider accounting for diffusion losses with expressions with the following generic form, 00 m_ i ¼ g ye yb ;
(1) 00
where the rate-of-reaction of species i (oxygen), m_ i kg/(m ,s), is given by a Butler-Volmer equation or the like. Phenomenological models for mass transfer mechanisms are typically invented to enumerate the mass transfer coefficient or conductance, g, in Eq. (1), which may easily be written in a volumetric form, as opposed to per unit area, and is often caste in spherical polar coordinates for idealised particles, or agglomerations of such particles. The molar form of Eq. (1), 00
n_i ¼ kðxe xb Þ;
2
(2)
is also widely found in the electrochemical literature, and is to be considered equivalent [1]. The purpose of this paper is to consider the class of formulation exemplified by Eq. (1), investigate the range of circumstances under which it is realistic, and to propose some modifications and extensions where necessary. In a previous work [2], the present author identified limitations of Eq. (1) when applied to problems where gases are employed as fuel and oxidant. It is generally accepted that the range of validity is confined to two situations: (i) very dilute mixtures and (ii) chemical catalysis with no net mass transfer at the boundary. The reader is referred to [3] for further details. Consider the idealised case shown in Fig. 2 where it is assumed that the cathode is in the form of a porous media, and that the solid electrode material is surrounded by a liquid, or liquid-like, water/Nafion® film, through which oxygen is transported. Water also is transferred either (a) into the liquid film and/or as vapour in the gas phase. For mass transfer in porous media it is common to employ volume-averaging techniques, and Fig. 3 illustrates schematically the volume-
E-mail address: [email protected]. http://dx.doi.org/10.1016/j.ijhydene.2015.05.074 0360-3199/Copyright © 2015, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved.
11642
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 1 1 6 4 1 e1 1 6 5 0
Fig. 1 e Typical PEM Cathode Assembly (not to scale).
Fig. 3 e Unravelled cathode showing nomenclature. averaged values of the bulk, wall, and interfacial values of the mass fraction which may be considered as locally 1-D for the purpose of the analysis. Let ye denote the mass fraction of the dissolved oxygen in the liquid film near the electrode wall, yl the value in the film near to the liquidegas interface, yg the oxygen mass fraction in the gas-phase near to the liquidegas interface and yb the bulk value in the gas free stream. One further concept is introduced, namely the notion of the transferred-substance, or T-state, value denoted by yt. The unfamiliar reader is referred to [3e5] and below, for details. Fig. 3 represents the very simplest possible idealisation of a cathode of a polymer electrolyte membrane fuel cell (PEMFC) for the purpose of the analysis. The air layer transfers oxygen to the liquid film, and thence to the solid electrode interface. The overall reaction at the electrode assembly may be written, 2H2 þ O2 þ 4aH2 O ¼ 2ð1 þ 2aÞH2 O
(3)
where the drag coefficient, a, is the net number of water molecules per proton emerging at the cathode: In a PEMFC, electro-
osmotic forces drag nd water molecules through the membrane with each proton. Electro-osmotic drag is modified (reduced) by back diffusion from the cathode to the anode, and also by a convective flux in the event of there being a significant pressure gradient between the two electrodes, to yield a net value a, typically less than 1. Dai et al. [6] review the experimentally determined values, and the results of calculations from mathematical models for water balance in PEMFC membrane electrode assemblies, and present ranges of values for both np and a. A mass balance may be written as follows, 00
m_ ¼
00 00 m_ H2 þ m_ H2 O
anode
ð4 þ 4a 18Þkg
¼ ¼
00 00 m_ O2 þ m_ H2 O
cathode
32kg þ 2ð1 þ 2aÞ18kg
(4)
The T-state value for the oxygen phase is defined as; . 00 00 yt ¼ m_ O2 m_
(5)
so that yt¼8/(1þ18a), as shown in Fig. 4. The water T-state is just 1eyt. Oxygen mass transfer is generally considered as being an important and potentially rate-limiting factor in the
Fig. 2 e Schematic of porous media, illustrating volumeaveraging process.
Fig. 4 e Transferred-state values of oxygen and water as a function of the net water drag coefficient, a
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 1 1 6 4 1 e1 1 6 5 0
operation of PEMFCs. The total mass flux is computed from Faradays law; 00
00
m_ ¼
MO2 i ð1000Þ 4Fyt
(6)
where i00 is current density and F is Faraday's constant. Eq. (6) equates mass flux and current density.
11643
More general analysis Let it be assumed that both water and oxygen are transferred at the electrode into the liquid film. Furthermore consider two limiting cases where; (a) no water is transferred to the gaseous phase and (b) all water is transferred to the gas. (a) Water in liquid phase
Simplified analysis based on equivalent film theory For a dilute mixture [5]
1
00 m_ O2 ¼ gg yg yb ¼ gl ye yl
(7)
00 m_ O2 ¼ g ye yb H
(8)
g ¼
gg gl
gg þ gl H gg
(9)
(13)
00 00 00 m_ O2 ¼ m_ yt ¼ m_ yt yg þ gg yg yb
(14)
which may be written in terms of the liquid-side and gas-side driving forces, Bl, Bg, respectively, 00 m_ ¼ gg Bg yt ¼ gl Bl
(15)
ye yt Bl ¼ yl ye
(16)
. yg 1 Bg ¼ yb yg
(17)
gl
where and are gas-side and liquid-side conductances or mass transfer coefficients, and g is an overall conductance. The non-dimensional Henry number [3] is, H¼
00 00 m_ yt ¼ m_ ye þ gl ye yl
yg ð1000Þ kH rl ¼ yl pg Mg
(10)
where .
kH ¼ pO2 ;g cl
(b) Water in gas phase 00 00 m_ yt ¼ m_ ye þ gl ye yl
(18)
00 00 m_ yt ¼ m_ yg þ gg yg yb
(19)
which leads to, (11)
00
m_ ¼ gg Bg ¼ gl Bl
(20)
is Henry's constant, cl is molar concentration of oxygen in the liquid film, and Ml and Mg are the mean mixture molecular weights of the liquid and gas (in g/mol). A molar Henry constant, H0 ¼ xg =xl , may also be defined. Various correlations for kH for the solubility of oxygen in water/Nafion® may be found in the literature [7,8]. For the liquid film let it be assumed that,
ye yt Bl ¼ yl ye
(21)
. yg yt Bg ¼ yb yg
(22)
gl ¼ rl Dl =dl
At high rates of mass transfer, the quantity g is no longer constant, and the expression, g/g*¼ln(1þB)/B, or equivalently, 00 m_ ¼ g* lnð1 þ BÞ, is exact in 1-D, and has been shown to be approximately true in many 3-D situations. Under these circumstances, in the liquid film
(12)
where Dl and dl are the oxygen diffusion coefficient and film thickness of the liquid, as appropriate. For large values of gg H=gl , the process is liquid-controlled 00 with m_ O2 ¼ gl ðye yb =HÞ, for small values, the gas phase con00 trols with m_ O2 ¼ gg ðHye yb Þ. The inversion of Eq. (8) to obtain a value of 00 ye ¼ yb =H þ m_ O2 =g for any given i00 is straightforward, and as such is referred to as the ‘simplified’ approach below. The reader will note that the right side of Eq. (8) is similar, but not identical, to that in Ref. [9]. It might be argued that in a PEMFC the dissolved oxygen is sufficiently dilute that Eq. (8) can be applied, compared to the situation in a solid oxide fuel cell where the gases are not dilute and Beale [2] shows Eq. (1) is not necessarily appropriate, consistent with the discussion on p. 67 of ref. [3]. The matter is further explored below.
(I) High mass transfer rate
00 yl ye m_ ¼ gl ln 1 þ ye yt (a) If all liquid water remains in the liquid film 00 yb yg m_ yt ¼ gg ln 1 þ yg 1
The use of the asterisk, g*, denotes the value of g in the limit, 00 g ¼ g lim m_ /0: By analogy Bird et al. [1] denote a similar value of k by k . At low mass flow rates, these may be considered constant. At high mass flow rates, however, g and k are variable. 1
(24)
(b) If all the water is transferred to the gas,
00
(23)
m_ ¼
gg ln
yb yg 1þ yg yt
! (25)
The solution to the system of Eqs 23e25 are referred to as the ‘full’ formulation, below.
11644
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 1 1 6 4 1 e1 1 6 5 0
(II) Low mass transfer rate In the PEMFC example given in ref. [10], it was noted that the blowing parameter was sufficiently small to linearise the 00 problem, g¼g*, m_ ¼ g* B. For low mass flow rates Eqs. 23 and 24 simplify to the following linearised form; In the liquid film 00
m_ ¼
h¼
00 m_ O2 yl ye ¼ gl yt ye yt
(26)
(a) If all liquid water remains in the liquid film yb yg yg 1
00
m_ yt ¼ gg
! (27)
(b) If all the water is transferred to the gas Eq. (25) reduces to 00
m_ ¼ gg
yb yg yg yt
to obtain the mass fraction at the electrode wall, ye, as a function of this and the bulk value, yb, by inverting the above system of equations. In general the (I) full, (II) linearised, and (III) simplified approaches will generate different numerical results for ye. In this paper we shall compare the results in terms of a concentration overpotential, h, defined here as
! (28)
Solution to Eqs 26e28 are referred to as the ‘linearised’ form, below. Let the non-dimensional blowing parameter, 00 b ¼ m_ g
(29)
clet be defined as the ratio of convection to diffusion i.e., a Pe number for mass transfer with bl and bg defined as mass flux normalised with respect to liquid and gas conductances, respectively.
Problem statement and solution methodology In many mass transfer problems in the chemical engineering sector, the rate of mass transfer at the boundary is unknown, and the system of equations above is used to obtain a solution to the problem. In fuel cells however, the mass transfer rate, 00 m_ , corresponding to some current density, i00 , may be considered as known (at least superficially), and it is required
RT ye H RT xe H0 ln ln ¼ 4F yb 4F xb
(30)
The calculation procedure for h(i00 ) is detailed in Table 1: For both (I) full and (II) linearised formulations; Given the 00 current density, i00 , the mass flux, m_ and gas blowing parameter, bg, are computed from Faraday's law, Eq. (6), and the definition of b, Eq. (29), respectively. This allows the gas driving force, Bg, and hence the gas-side oxygen mass fraction at the interface, yg, to be computed from known gas bulk value, yb, and conductance, gg . Precise mathematical details vary depending according to whether the linearised or fullform is assumed, and also whether it is presumed all water is transferred to (a) liquid or (b) gas bulk state. Full details of the method for the sequence of calculations inverting the mass transfer equations to obtain ye from yb for each case are detailed in Table 1. After obtaining the oxygen mass fraction on the liquid-side of liquidegas interface from Henry's law, Eq. (10), the above procedure is repeated for the liquid film, and the value at the electrode wall finally computed from the blowing parameter/driving force identity for the liquid film. For the simplified case (III) it is not necessary to proceed from the gas bulk to the electrode wall in a sequential fashion, since Eq. (8) may easily be inverted to write ye in terms of yb and g . Another parameter of interest is the so-called limiting 00 current density, ilim , which for cases (I) and (II) is obtained by 00 iterating on values of the mass flux, m_ , i.e. i00 , until an electrode-value, ye¼0 is obtained. This is achieved using any simple root-finding algorithm. For case (III) the simplified 00 method, no such iteration is needed and ilim is computed directly, as detailed in Table 1. It is, of course, straightforward to write analytical expressions for ye(bl,bg), for cases (I) and (II), though these are non-linear expressions, and inversion is therefore most easily accomplished numerically.
00
Table 1 e Methods used to compute ye and ilim for all water transferred to (a) liquid state, (b) gas state. Method I. Full: Eqs. 23e25
(a) (b)
II. Linearised: Eqs. 26e28
(a) (b)
III. Simplified: Eq. (8)
00
ye
ilim
Bg ¼ exp(ytbg1) yg ¼ (ybþBg)/(1þBg) Bg ¼ exp(bg1) yg ¼ (ybþBgyt)/(1þBg) yl ¼ yg/H Bl ¼ exp(bl1) ye ¼ (ylþBlyt)/(1þBl) Bg ¼ ytbg yg ¼ (ybþBg)/(1þBg) Bg ¼ b g yg ¼ (ybþBgyt)/(1þBg) yl ¼ yg/H Bl ¼ b l ye ¼ (ylþBlyt)/(1þBl) 00 ye ¼ yb =H þ m_ =g
Iterate on values of i00 until ye ¼ 0
Iterate on values of i00 until ye ¼ 0
ilim ¼ 4000Fg yb =MO2 H 00
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 1 1 6 4 1 e1 1 6 5 0
11645
Table 2 e Prescribed values, baseline case. Temperature, T Pressure, p O2 gas mass fraction Net drag coefficient, a Gas conductance, gg Henry constant, kH Liquid density, rl Oxygen liquid diffusivity, Dl Liquid film thickness, dl
353 K 1 Atm 0.232(21%vol.)0.116(10.5%vol.) 0.2 0.1 kg(m2 s) 0.2 Atm m3/mol 9.718 kg/m3 108m2/s 109m
Results Fig. 6 e Limiting current density as a function of the conductance (a) no water transferred to gas.
Table 2 show values prescribed to perform sample mass transfer calculations. For the gas phase, yb is taken as the bulk value in the air channel, Fig. 1, and the gas mass transfer coefficient, gg , is estimated as the harmonic average of the combined channel and gas diffusion layer values from Ref. [11]. Two bulk values were considered for the air corresponding to 21% vol. oxygen (i.e., atmospheric air) and 10.5% vol. (i.e., 50% depleted air). Sample calculations were performed to investigate the sensitivity of the results to the choice of equations. Fig. 5 show concentration overpotential, defined by Eq. (30) vs. current density, i00 , for (a) no gaseous water transfer. The (I) high mass rate formulation is compared to (II) the low mass rate linearised formulation as well as (III) the simplified theory. It can be seen that a limiting current density is reached for yb ¼ 0.116 at (I) 1.39105 A/m2, (II) 1.32105 A/m2 and (III) 1.31105 A/m2. For yb ¼ 0.232 the values are (I) 2.97105 A/m2, (II) 2.65105 A/m2, and (III) 2.61105 A/m2. Thus, for the set of parameters chosen in Table 1, there is only a difference of 00 about 6e12% in ilim , with the full approach (I) generating a 00 maximum value of ilim and the simplified method (III) a minimum value. The Henry number has a value of 6.7103 corresponding to kH¼0.2Atm,m3/mol. The corresponding value of gg H=gl is 0.069. Fig. 6 shows limiting current density computed as a function of the non-dimensional parameter gg H=gl on the basis that (a) all produced water remains in a liquid state, with KH¼0.2 Atm , m3/mol. This figure was generated by varying
the ratio Dl/dl, i.e. gl , over a range of values reported in the literature. It can be seen that the limiting current increases with bulk mass fraction of oxygen, and decreases as gg H=gl increases. Fig. 7 shows results similar to Fig. 5, this time on the assumption that (b) all the produced water enters the gas phase. It can be seen that for yb¼0.116 a limiting current density is reached, at values of (I) 1.27105 A/m2, (II) 1.30105 A/m2, and (III) 1.31105 A/m2. For yb¼0.232 a limiting current density is reached, at values of (I) 2.46105 A/m2, (II) 2.60105 A/m2, and (III) 2.61105 A/m2, with the full approach (I) 00 generating the minimum ilim and the simplified theory (III) about 3e6% higher. Fig. 8 shows the limiting current density on the premise (b) that all the produced water enters the gas phase. It can be seen that the detailed mass transfer formulation predicts limiting current densities that are lower than the simple equivalent film theory in contrast to the results of Fig. 6, where the converse is true. Fig. 9 shows the limiting current density as a function of the net drag coefficient, a, for case (b) produced water transferred to gas. All other parameters are held to the baseline values given in Table 2. It can be seen that as a is increased from 0.5 to þ1.0, the full model (I) indicates a
Fig. 5 e Concentration overpotential for the baseline case (a) No water transferred to gas.
Fig. 7 e Concentration overpotential for the baseline case (b) All water transferred to gas.
11646
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 1 1 6 4 1 e1 1 6 5 0
Fig. 8 e Limiting current density as a function of the conductance (b) all water transferred to gas.
substantial reduction in limiting current from 3.0105 A/m2 to 2.1105 A/m2 compared to (II) the much smaller decline from 2.6105 A/m2 to 2.5105 A/m2 for the linearised model at yb¼0.232. The (III) simplified model does not account for changes in limiting current as a function of water content and therefore displays a horizontal slope. The three profiles coin00 cide precisely at a value of ilim ¼ 2:6 105 A=m2 , with a¼1/18 corresponding to yt¼±∞, Fig. 4. For case (a) not shown, there is 00 no significant change in ilim as a function of a. The matter is discussed further below.
Discussion Inspection of Figs. 5 and 7 suggests that the results of the simple equivalent film theory (III) agree well with the low mass formulation (II), and quite well with the high mass-rate formulation (I) when the oxygen mass fraction is dilute, yb¼0.116. This is consistent with ref. [2] which indicates that for yb¼0.1, the simplified theory (III) generates reasonable results for oxygen consumption. However, as the oxygen bulk mass fraction is increased, e.g. to yb¼0.232, the errors increase, and will do so further in the event that enriched air were employed as oxidant.
Fig. 9 e Limiting current density as a function of the net transfer coefficient, a. (b) All water transferred to gas.
Within the liquid film, oxygen mass fractions are always extremely small with H¼6.7103. This means the maximum ye is 1.7105 for yb¼0.116 and 2.4105 for yb¼0.232 for i00 ¼0. Clearly the condition jyej≪jytj is met for the liquid film. For the choice of parameters in Table 2, the gas controls mass transfer more than the liquid, as gg H=gl is 0.069 i.e., the gas mass transfer resistance is larger. In fact the value of gg H is only about 7% lower than the value of g computed by Eq. (9). With (a) no water being transferred into the gas phase, Fig. 5, at i00 ¼105A/m2, Bg and Bl are of the order of 0.08 and þ5107, respectively i.e., Bg is quite small and Bl is tiny. (The sign of B indicates that the net mass transfer is out of the gas but into the liquid.) The values are the same at both yb¼0.116 and 0.232. For i00 ¼2105A/m2, Bg and Bl are of the order of 0.15 and 106. Thus the simplified theory (III) approximates the low mass transfer linearised method (II) quite well, with some minor differences between these and the full formulation (I), due to the fact that even with bg small, values of Bg¼exp(bg)1 and Bg¼bg are slightly different, with the difference increasing 00 00 as i /ilim . In other words, yg is also quite small compared to unity in (though not as small as ye), Eq. (17). For this reason, there are relatively minor differences in the limiting current density computed by the three methods. Increasing the oxygen bulk value, the differences start to become more significant, as might be expected. With (b) all the water being transferred into the gas phase, Fig. 7, values of the liquid driving force, Bl, are similar to those for the results of Fig. 5, but the gas values differ both in magnitude and sign with Bg around þ0.05 and þ 0.1 for i00 ¼105A/m2 and 2105A/m2,respectively. This change of sign is due to the fact that the while direction of the net mass flux of water and oxygen is (a) out of the gas for the results of Fig. 5, it is into the gas for the latter case (b), due to the addition of the water vapour to the air. Comparison of Figs. 5 and 7 reveal that whereas the simple theory (III), predicts the same rate of mass transfer when either (a) all water remains in the film or (b) all water is transferred to the gas; both the high-mass (I), and linearised low-mass (II) formulations suggest that higher limiting current densities and lower concentration overpotential magnitudes occur when the water remains in the liquid rather than the vapour state. Thus a comparison of method (I) for limiting current densities of (a) 1.39105 A/m2 and 2.97105 A/m2 with (b) 1.27105 A/m2 and 2.46105 A/m2 at yb¼0.116 and 0.32 suggests a 10e21% difference between cases (a) and (b) for this performance measure. Inspection of Figs. 6 and 8 reveals the limiting current is reduced as values of the parameter gg H=gl increase. There is 00 some difference between values of ilim predicted by the highmass theory and the others, notably for small gg H=gl ; however the linearised and simplified agreements are by-and-large, in agreement. The oxygen concentration overpotential, defined by Eq. (30) is a measure of decrease in potential due to mass transfer effects, that is the difference in Nernst potential associated with the electrode mass fraction being ye rather than yb/H. Because of this, the ordinates in Figs. 5 and 7 are proportional to a logarithmic quantity, ln(yeH/yb). This has the effect of compressing the results, thereby minimising differences. If instead the change in current density as opposed to potential
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 1 1 6 4 1 e1 1 6 5 0
is considered; and if it is assumed cathode reaction occurs in the Tafel zone. Then at any given electrode potential, the current density will be reduced by the exponential of (the change in) the potential compared to the value if the electrode were at yb/H, i.e., the reduction in current density is proportional to yeH/yb and not ln(yeH/yb). The current density, and mass flow rate are nominally proportional to the mass transfer driving force, and it can be shown that for a given cell voltage, the current density may be written, . 00 00 (31) i ib ¼ expðb=2Þ ye H yb 00
00
where b is a charge transfer coefficient, and ib denotes i when ye¼yb/H'. Note that Eq. (31) is an idealisation, since the exchange current density will also be a function of oxygen concentration, further complicating the relationship. There will also be a small water concentration-overpotential; since the water mass fraction, will also differ. This should not be confused with the problem of liquid-water blocking oxygen transport in porous media [10]. The ratio of the parameters Dl and dl is important, and clearly affects the value of gl in Eq. (12). Various values may be found in the literature; dl-values ranging from 2 nm to 80 nm have been recorded [12,13]. Oxygen liquid-diffusivity correlations also vary significantly, with published values of Dl in the range of 1.21010m2/s to 2108m2/s [7,9,14]. Variations in Henry constants are smaller, but by no means insignificant. The parameter gl could be regarded simply as a fitted parameter, to be obtained from limiting current density data for a given cathode, assuming gg is known. As the parameter gg H=gl approaches zero, the mass transfer calculation can be performed based on the gas-side alone [3]. If all that is of interest is the oxygen concentration overpotential, it may be argued that one may adopt the Henry's law-modified simplified theory (III) as a reasonable first approximation to predict oxygen mass transfer losses in the cathode of a PEMFC. The oxygen concentration overpotential is readily obtained [2] from Eq. (8) in the well-known form, 00
RT i ln 1 00 h¼ 4F ilim
! (32)
where ilim ¼ 4000Fg yb =MO2 H as given in Table 1. The reader 00 will note that ilim , as predicted by method (III), is independent of yt and hence a. The same is not true for methods (I) and (II) 00 where it can be seen that ilim decreases as a is increased. Method (III) is a fair approximation under the present circumstance, because the driving forces are small, and the magnitude of yt is for the most-part substantially greater than that of yg, and thus Eq. (7) is therefore a reasonable approximation to the full equations. Simplicity in a mathematical model is a good thing, however, when it comes to the calculation for the rate of mass transfer of water, Eq. (7) is not suitable except when the water mass fractions are very small. Thus if Eq. (7) is used to compute or fit a value of the overall mass transfer coefficient, simply substituting the value of g into the conductance in Eq. (1) for mass transfer of water is illadvised. The matter is the subject of a future parametric study. Thus the simplified method does not provide a rationale foundation to obtain all the mass fractions of the component gases simultaneously, and a basic understanding 00
11647
of the full analysis is important. This is corroborated by the fact that many authors, for example [15], have noted problems with use of the term RT/4F in Eq. (32) and suggest instead the 00 00 use of a fitted-constant, C, with h ¼ Cln 1 i =ilim Þ. This has no theoretical basis in the field of mass transfer, and is merely a reflection that the logic underlying methodology (III) is too simplified to be universally applicable, except for dilute mixtures [2]. The results of Fig. 9 show clearly that the impact of net water drag coefficient, a, upon performance cannot be predicted by the simplified approach (III). Here too some significant differences between the (I) full and (II) linearised models are apparent. A diverse range of a values have been reported in the literature, with some authors claiming it to be a function of current density and others not [6]. The three methodologies predict identical limiting currents for a¼1/18 corresponding to yt¼±∞, when there is no net mass transfer with the mass flux of oxygen being precisely equal and opposite to that of the water. As the magnitude of a deviates from this value, so does agreement between the models for 00 ilim . Because the magnitude of Bl is small in comparison to Bg, the effect is important for case (b) all water transferred to gas, but not case (a) all water remains as liquid. In most circumstances, it is a requirement to obtain the mass fractions of all the species; oxygen, water, nitrogen, etc. as part of global or local calculations, where activation and ohmic terms will also affect the operating voltage/current 00 density which in turn affects the value of i00 and hence m_ . If the diffusion coefficients of the gases are all nominally equal, the full/linearised methodology provides a simple method for computing these values since BO2 ¼ BH2 O ¼ BN2
(33)
allows each gas-phase mass fraction to be computed from the mass flow rate, bulk and T-state values. This type of engineering approximation was frequently made for atmospheric air and other gases, prior to the advent of detailed computational models. When the binary diffusivities of the gas component-pairs differ considerably, as might be the case for mixtures of gases with hydrogen, there is some cause for concern over the applicability of this type of model, per se. In PEMFCs often only a binary mixture of humidified hydrogen is employed on the fuel side, though in other fuel cells and electrolysers, situations could arise where there is a need to employ higher-order numerical schemes [16] based on MaxwelleStefan equations or generalised Fick's law [17]. On the other hand, it is to be noted that many mass transfer Stanton number correlations (used to calculate the conductance g*) exhibit a relatively weak dependence upon the Schmidt number, the ratio of diffusivity and kinematic velocity. In any event, there is a role for both analytical [18] and detailed 3-D models in the development of fuel cell technology. For a lumped-capacity analysis as typically found in 0-D and 1-D analyses [19], the choice of yb, corresponding to the bulk value in the channel passage as reference, is appropriate. The associated length-scale is the sum of the size of the diffusion layer and half-hydraulic diameter of the channel. In a fuel cell operating at high air utilisation, there will be a significant decrease in oxygen mass fraction from inlet to outlet, which must be accounted-for: In a lumped capacity model, the bulk value should be a mean value of the inlet and
11648
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 1 1 6 4 1 e1 1 6 5 0
outlet mass fractions, the latter being computed from a mass balance for oxygen and water. This is a function of the total current and the partial densities, ri¼ryi, of oxygen, water (and nitrogen). Therefore calculations on the water mass fraction are important. The same is true for simple multi-dimensional ‘along the channel’ models, where the bulk mass fractions are adjusted cell-by-cell. For all these reasons, the (I) full or (II) linearised mass transfer formulations based on general theory, as outlined in this paper are recommended for consideration when performing calculations on mass transfer in all but the simplest PEMFC problems. In 3-D computational fluid dynamics (CFD) calculations, the fuel cell geometry is typically tessellated with a computational mesh, and the species equations solved-for in the channels, diffusion layers and electrodes together with the equations for mass, momentum, and energy [19]. The mass transfer formulation, described here, may be considered a sub-scale closure within the electrode and yb chosen as the local bulk gas value within the electrode pores, as opposed to within the gas channel, Fig. 2. This local bulk value is not constant but changes as a function of position, local current density, etc. The gas length scale is therefore of the order of the pore size, and the gas conductance greatly increased, accordingly. Thus gg H=gl [0, i.e., the liquid-film is the ‘controlling’ or rate-limiting factor, however the CFD code should at the very least account for solubility effects on the Nernst equation and current density and if possible be capable of calculating bulk-to-interface driving forces, although depending on the pore size these may be small. As discussed above, an estimate of the overall mass transfer coefficient may be made from the limiting current density, and vice versa. However it is only an estimate. Moreover Henry's law is to be considered as an approximation at finite rates-of mass transfer [3,20]. A major assumption of the mass transfer model presented here, is the presence of a simple liquid film of effective height, dl, since the height will almost certainly vary with current density and protonic drag. In this paper the limiting cases of all produced water, either (a) remaining in the liquid film or (b) entering the gas phase were considered. In practice this needs to be computed either locally in a CFD or presumed-flow type code, or globally in a lumped-parameter analysis. Although this is considered beyond the scope of this paper, some brief discussion is in order. The situation is more typical of conventional mass transfer problems, where the rate of mass transfer is not known, a priori, but rather is the required dependent variable, and may be classified as a ‘mass-from-heat transfer’ problem for which there have been many solved examples in chemical and nuclear engineering over the past four decades [21,22]. Let the mass fraction of water vapour be denoted by yv. Under the assumptions that (i) the latent heat of vapourisation, Hfg, is constant, (ii) the vapour specific-heat is much larger than that of the liquid, and (iii) the gas is perfect, the Clausius-Clapeyron relationship may be stated in the form, pv ðTÞ ¼ C exp Hfg RT
(34)
Empirical correlations such as the August-Roche-Magnus equation, are readily employed in place of Eq. (34) to formulate an explicit expression for pv as a function of T. If the vapour partial-pressure, pv, is known at the interface then the
temperature T(pv) is also known, and vice-versa. Once the interfacial temperature is fixed, as a Dirichlet condition, the unknown liquid-to-vapour mass flux can be obtained from the enthalpy equation, and local bulk gas and liquid temperatures. The combined flux of both dissolved oxygen, and water from liquid to vapour, is a function of the T-state value which is dependent on the individual mass fluxes (of oxygen and water). Thus, in addition to Henry's law, Eq. (34) gives the unknown water-steam mass flux, and hence the correct transferred substance states in the gas mixture 00
m_ ¼
hg Tb;g T þ hl ðTb;l TÞ Hfg
(35)
where hl and hg are liquid-side and gas-side heat transfer coefficients, Tb,l and Tb,g are bulk values in the two phases, and Hfg is evaluated at the interface temperature, T. This would normally be done using a sub-scale closure in a CFD-based calculation, e.g. in the gas diffusion layer and/or electrodes. An alternative approach [23] is to compute the mass fraction and mass transfer rate from the partial vapour pressure, pv(T), given the interface temperature via Eq. (34). The associated inter-phase heat flux may then be prescribed as a Neumann condition in the energy equation. In multi-phase CFD codes it is not generally necessary nor desirable to presume thermal equilibrium between phases, since the phase temperatures can differ, and moreover if there is no heat transfer, then there can be no mass transfer by evaporation/condensation and vice versa. Regarding equilibrium of species, it is generally true that the partial vapour pressure will only deviate by a small amount from the equilibrium value for finite rates of mass transfer [3]. In catalyst layers, where the liquid film is not saturated water, many authors use the closure due to Springer et al. [24] where the dependence of the liquid concentration on pv is prescribed empirically as a function of temperature and membrane composition (humidification). This is a natural extension of the approach above and may be considered as a similar, but more complex case than can be described by idealised solution theory (Raoult's law). The additional water mass flux must therefore be accounted for by modifying the mass flux and transferred substance-states from the two limiting cases prescribed in Eqs. (14) and (19) to additionally account for the fraction of water which is evaporating (or condensing). The mass transfer analysis described in this paper may also be used to consider heat transfer, by analogy; Driving force(s) and blowing parameter(s) may all be similarly defined in terms of temperature or enthalpy. The challenge facing the engineer in this context is to obtain reliable correlations for phase heat and mass transfer coefficients for the complex geometries and flow regimes frequently encountered in PEMFC designs.
Conclusions A standard mechanical-engineering type analysis was developed and applied to an idealised PEMFC cathode. Both high and low mass transfer models were derived and compared to a highly-simplified method. It was shown that if the oxidant is composed of atmospheric air with a maximum value of 21%
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 1 1 6 4 1 e1 1 6 5 0
vol., the simplified method gives a reasonable first estimate of the limiting current density. The simplified method agrees well with the (linearised) low mass transfer method, and has the advantage that the limiting current density may be computed directly and simply. If, however, it is required that the gas mass fractions of all species be computed simultaneously, these should be derived from the driving force approach-based detailed analysis presented in this manuscript. Moreover, if enriched oxygen is employed, the standard approach is a prerequisite requirement for a reliable solution. It was shown that the limiting current density is reduced by up to 21% if water is transferred to the gaseous state, rather than remaining in entirely in the liquid state; something which cannot be predicted by the simplified model. The ratio of the gas and liquid conductances, gg H=gl , governs the degree to which the gas or liquid controls mass transfer, with limiting current density decreasing as gg H=gl increase. This value was quite small for the particular baseline parameters prescribed in this study.
Acknowledgement Prof. Werner Lehnert and Dr. Uwe Reimer both read the manuscript, and contributed useful commentary. I would also like to thank Dr. Michael Malin for remarks on the phase change problem and especially Prof. Brian Spalding for personally sharing his unique insight into heat and mass transfer, which directly influenced the subject matter of this paper. The support of the Forschungszentrum Ju¨lich is gratefully acknowledged.
Nomenclature English Letters A area, m2 b fitted constant, V c concentration, mol/m3 D diffusion coefficient, m2/s F Faraday's constant, 96485.3 C/mol g mass transfer coefficient, kg/(m2,s) 00 mass transfer conductance in the lim m_ /0, kg/ g* 2 (m ,s) h heat transfer coefficient, W/(m2,K) latent heat of vapourisation, J/kg Hfg k molar conductance, mol/(m2,s) 00 molar conductance in the lim n_ /0, mol/(m2,s) k Henry's constant, (Pa,m3)/mol kH 00 mass flux, kg/(m2,s) m_ M molecular weight, g/mol current density, A/m2 i00 00 molar flux, mol/(m2,s) n_ electro-osmotic drag coefficient nd p pressure, Pa partial vapour pressure, Pa pv R gas constant, 8.3144 J/(K,mol) t time, s
T x y
11649
temperature, K molar fraction mass fraction
Greek Letters a net protonic drag coefficient b charge transfer coefficient dl liquid film thickness, m h concentration overpotential, V r mixture density, kg/m3 t wall shear stress, Pa Non-dimensional number b blowing parameter B driving force H Henry number molar Henry number H0 Subscripts b bulk value e electrode wall value g gas, gas-side interface value hydrogen H2 H2O water i species i l liquid, liquid-side interface value lim limiting value nitrogen N2 oxygen O2 t transferred substance state value Superscripts overall value
references
[1] Bird RB, Stewart WE, Lightfoot EN. Transport phenomena. New York: Wiley; 1960. [2] Beale SB. Calculation procedure for mass transfer in fuel cells. J Power Sources 2004;128(2):185e92. [3] Spalding DB. Convective mass transfer: an introduction. Edward Arnold; 1963. [4] Kays WM, Crawford ME, Weigand B. Convective heat and mass transfer. 4th ed. New York: McGraw-Hill; 2005. [5] Spalding DB. A standard formulation of the steady convective mass transfer problem. Int J Heat Mass Transf 1960;1:192e207. [6] Dai W, Wang H, Yuan X-Z, Martin JJ, Yang D, Qiao J, et al. A review on water balance in the membrane electrode assembly of proton exchange membrane fuel cells. Int J Hydrog Energy 2009;34(23):9461e78. [7] Bernardi DM, Verbrugge MW. Mathematical model of a gas diffusion electrode bonded to a polymer electrolyte. AIChE J 1991;37(8):1151e63. [8] Parthasarathy A, Srinivasan S, Appleby AJ, Martin CR. Pressure dependence of the oxygen reduction reaction at the platinum microelectrode/Nafion interface: electrode kinetics and mass transport. J Electrochem Soc 1992;139(10):2856e62. [9] Um S, Wang C-Y, Chen KS. Computational fluid dynamics modelling of proton exchange membrane fuel cells. J Electrochem Soc 2000;147:4485e93.
11650
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 1 1 6 4 1 e1 1 6 5 0
[10] Beale SB, Schwarz DH, Malin MR, Spalding DB. Two-phase flow and mass transfer within the diffusion layer of a polymer electrolyte membrane fuel cell. Comput Therm Sci 2009;1(2):105e20. [11] Beale SB. Conjugate mass transfer in gas channels and diffusion layers of fuel cells. ASME J Fuel Cell Sci Technol 2007;4:1e10. [12] Harvey D, Pharoah JG, Karan K. A comparison of different approaches to modelling the PEMFC catalyst layer. J Power Sources 2008;179(1):209e19. [13] Sadeghi E, Putz A, Eikerling M. Hierarchical model of reaction rate distributions and effectiveness factors in catalyst layers of polymer electrolyte fuel cells. J Electrochem Soc 2013;160(10):F1159e69. [14] Parathasarthy A, Srinivasan S, Appleby AJ, Martin CR. Temperature dependence of the electrode kinetics of oxygen reduction at the platinum/Nafion interface - a miroelectrode investigation. J Electrochem Soc 1992;139. 2530e2357. [15] Larminie J, Dicks A. Fuel cell systems explained. 2nd ed. New York: Wiley; 2003. nchez-Insa A, Fueyo N. Multimodal [16] Garcı´a-Camprubı´ M, Sa mass transfer in solid-oxide fuel-cells. Chem Eng Sci 2010;65(5):1668e77. [17] Taylor R, Krishna R. Multicomponent mass transfer. John Wiley & Sons; 1993.
[18] Kulikovsky AA. Analytical modelling of fuel cells. Elsevier; 2010. [19] Weber AZ, Borup RL, Darling RM, Das PK, Dursch TJ, Gu W, et al. A critical review of modeling transport phenomena in polymer-electrolyte fuel cells. J Electrochem Soc 2014;161(12):F1254e99. [20] Perry RH, Green DW, Maloney JO. Perry's chemical engineers' handbook. New York: McGraw-Hill; 1997. [21] Spalding DB. Numerical computation of multi-phase fluid flow and heat transfer. In: Taylor C, Morgan K, editors. Recent advances in numerical methods in fluids, vol. 1. Pineridge Press; 1980. p. 139e67. Ch. 5. [22] Spalding DB. Mathematical methods in nuclear-reactor thermal hydraulics. In: Lahey R, editor. Proceedings of American nuclear society, meeting on nuclear-reactor thermal hydraulics. Saratoga, New York: American Nuclear Society; 1980. p. 1979e2023. [23] Berning T, Djilali N. A 3D, multiphase, multicomponent model of the cathode and anode of a PEM fuel cell. J Electrochem Soc 2003;150(12):A1589e98. [24] Springer TE, Zawodzinski TA, Gottesfeld S. Polymer electrolyte fuel cell model. J Electrochem Soc 1991;138(8):2334e42.
Available online at www.sciencedirect.com
ScienceDirect journal homepage: www.elsevier.com/locate/he
Mass transfer formulation for polymer electrolyte membrane fuel cell cathode Steven B. Beale Institute of Energy and Climate Research, IEK-3, Forschungszentrum Ju¨lich GmbH, 92425 Ju¨lich, Germany
article info
abstract
Article history:
A standard mechanical engineering mass transfer analysis is adapted to obtain perfor-
Received 3 February 2015
mance calculations for the cathode of a typical polymer electrolyte membrane fuel cell.
Received in revised form
Mathematical details are provided and sample calculations performed, based on a number
4 May 2015
of simplifying assumptions. The relative advantages of low and high-mass transfer for-
Accepted 8 May 2015
mulations are critically appraised and compared to a highly-simplified method.
Available online 22 July 2015
Copyright © 2015, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved.
Keywords: Fuel cells Transport phenomena Mass transfer Concentration polarisation Diffusion over-potential
Introduction It is popular within some elements of the fuel cell modelling community to consider accounting for diffusion losses with expressions with the following generic form, 00 m_ i ¼ g ye yb ;
(1) 00
where the rate-of-reaction of species i (oxygen), m_ i kg/(m ,s), is given by a Butler-Volmer equation or the like. Phenomenological models for mass transfer mechanisms are typically invented to enumerate the mass transfer coefficient or conductance, g, in Eq. (1), which may easily be written in a volumetric form, as opposed to per unit area, and is often caste in spherical polar coordinates for idealised particles, or agglomerations of such particles. The molar form of Eq. (1), 00
n_i ¼ kðxe xb Þ;
2
(2)
is also widely found in the electrochemical literature, and is to be considered equivalent [1]. The purpose of this paper is to consider the class of formulation exemplified by Eq. (1), investigate the range of circumstances under which it is realistic, and to propose some modifications and extensions where necessary. In a previous work [2], the present author identified limitations of Eq. (1) when applied to problems where gases are employed as fuel and oxidant. It is generally accepted that the range of validity is confined to two situations: (i) very dilute mixtures and (ii) chemical catalysis with no net mass transfer at the boundary. The reader is referred to [3] for further details. Consider the idealised case shown in Fig. 2 where it is assumed that the cathode is in the form of a porous media, and that the solid electrode material is surrounded by a liquid, or liquid-like, water/Nafion® film, through which oxygen is transported. Water also is transferred either (a) into the liquid film and/or as vapour in the gas phase. For mass transfer in porous media it is common to employ volume-averaging techniques, and Fig. 3 illustrates schematically the volume-
E-mail address: [email protected]. http://dx.doi.org/10.1016/j.ijhydene.2015.05.074 0360-3199/Copyright © 2015, Hydrogen Energy Publications, LLC. Published by Elsevier Ltd. All rights reserved.
11642
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 1 1 6 4 1 e1 1 6 5 0
Fig. 1 e Typical PEM Cathode Assembly (not to scale).
Fig. 3 e Unravelled cathode showing nomenclature. averaged values of the bulk, wall, and interfacial values of the mass fraction which may be considered as locally 1-D for the purpose of the analysis. Let ye denote the mass fraction of the dissolved oxygen in the liquid film near the electrode wall, yl the value in the film near to the liquidegas interface, yg the oxygen mass fraction in the gas-phase near to the liquidegas interface and yb the bulk value in the gas free stream. One further concept is introduced, namely the notion of the transferred-substance, or T-state, value denoted by yt. The unfamiliar reader is referred to [3e5] and below, for details. Fig. 3 represents the very simplest possible idealisation of a cathode of a polymer electrolyte membrane fuel cell (PEMFC) for the purpose of the analysis. The air layer transfers oxygen to the liquid film, and thence to the solid electrode interface. The overall reaction at the electrode assembly may be written, 2H2 þ O2 þ 4aH2 O ¼ 2ð1 þ 2aÞH2 O
(3)
where the drag coefficient, a, is the net number of water molecules per proton emerging at the cathode: In a PEMFC, electro-
osmotic forces drag nd water molecules through the membrane with each proton. Electro-osmotic drag is modified (reduced) by back diffusion from the cathode to the anode, and also by a convective flux in the event of there being a significant pressure gradient between the two electrodes, to yield a net value a, typically less than 1. Dai et al. [6] review the experimentally determined values, and the results of calculations from mathematical models for water balance in PEMFC membrane electrode assemblies, and present ranges of values for both np and a. A mass balance may be written as follows, 00
m_ ¼
00 00 m_ H2 þ m_ H2 O
anode
ð4 þ 4a 18Þkg
¼ ¼
00 00 m_ O2 þ m_ H2 O
cathode
32kg þ 2ð1 þ 2aÞ18kg
(4)
The T-state value for the oxygen phase is defined as; . 00 00 yt ¼ m_ O2 m_
(5)
so that yt¼8/(1þ18a), as shown in Fig. 4. The water T-state is just 1eyt. Oxygen mass transfer is generally considered as being an important and potentially rate-limiting factor in the
Fig. 2 e Schematic of porous media, illustrating volumeaveraging process.
Fig. 4 e Transferred-state values of oxygen and water as a function of the net water drag coefficient, a
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 1 1 6 4 1 e1 1 6 5 0
operation of PEMFCs. The total mass flux is computed from Faradays law; 00
00
m_ ¼
MO2 i ð1000Þ 4Fyt
(6)
where i00 is current density and F is Faraday's constant. Eq. (6) equates mass flux and current density.
11643
More general analysis Let it be assumed that both water and oxygen are transferred at the electrode into the liquid film. Furthermore consider two limiting cases where; (a) no water is transferred to the gaseous phase and (b) all water is transferred to the gas. (a) Water in liquid phase
Simplified analysis based on equivalent film theory For a dilute mixture [5]
1
00 m_ O2 ¼ gg yg yb ¼ gl ye yl
(7)
00 m_ O2 ¼ g ye yb H
(8)
g ¼
gg gl
gg þ gl H gg
(9)
(13)
00 00 00 m_ O2 ¼ m_ yt ¼ m_ yt yg þ gg yg yb
(14)
which may be written in terms of the liquid-side and gas-side driving forces, Bl, Bg, respectively, 00 m_ ¼ gg Bg yt ¼ gl Bl
(15)
ye yt Bl ¼ yl ye
(16)
. yg 1 Bg ¼ yb yg
(17)
gl
where and are gas-side and liquid-side conductances or mass transfer coefficients, and g is an overall conductance. The non-dimensional Henry number [3] is, H¼
00 00 m_ yt ¼ m_ ye þ gl ye yl
yg ð1000Þ kH rl ¼ yl pg Mg
(10)
where .
kH ¼ pO2 ;g cl
(b) Water in gas phase 00 00 m_ yt ¼ m_ ye þ gl ye yl
(18)
00 00 m_ yt ¼ m_ yg þ gg yg yb
(19)
which leads to, (11)
00
m_ ¼ gg Bg ¼ gl Bl
(20)
is Henry's constant, cl is molar concentration of oxygen in the liquid film, and Ml and Mg are the mean mixture molecular weights of the liquid and gas (in g/mol). A molar Henry constant, H0 ¼ xg =xl , may also be defined. Various correlations for kH for the solubility of oxygen in water/Nafion® may be found in the literature [7,8]. For the liquid film let it be assumed that,
ye yt Bl ¼ yl ye
(21)
. yg yt Bg ¼ yb yg
(22)
gl ¼ rl Dl =dl
At high rates of mass transfer, the quantity g is no longer constant, and the expression, g/g*¼ln(1þB)/B, or equivalently, 00 m_ ¼ g* lnð1 þ BÞ, is exact in 1-D, and has been shown to be approximately true in many 3-D situations. Under these circumstances, in the liquid film
(12)
where Dl and dl are the oxygen diffusion coefficient and film thickness of the liquid, as appropriate. For large values of gg H=gl , the process is liquid-controlled 00 with m_ O2 ¼ gl ðye yb =HÞ, for small values, the gas phase con00 trols with m_ O2 ¼ gg ðHye yb Þ. The inversion of Eq. (8) to obtain a value of 00 ye ¼ yb =H þ m_ O2 =g for any given i00 is straightforward, and as such is referred to as the ‘simplified’ approach below. The reader will note that the right side of Eq. (8) is similar, but not identical, to that in Ref. [9]. It might be argued that in a PEMFC the dissolved oxygen is sufficiently dilute that Eq. (8) can be applied, compared to the situation in a solid oxide fuel cell where the gases are not dilute and Beale [2] shows Eq. (1) is not necessarily appropriate, consistent with the discussion on p. 67 of ref. [3]. The matter is further explored below.
(I) High mass transfer rate
00 yl ye m_ ¼ gl ln 1 þ ye yt (a) If all liquid water remains in the liquid film 00 yb yg m_ yt ¼ gg ln 1 þ yg 1
The use of the asterisk, g*, denotes the value of g in the limit, 00 g ¼ g lim m_ /0: By analogy Bird et al. [1] denote a similar value of k by k . At low mass flow rates, these may be considered constant. At high mass flow rates, however, g and k are variable. 1
(24)
(b) If all the water is transferred to the gas,
00
(23)
m_ ¼
gg ln
yb yg 1þ yg yt
! (25)
The solution to the system of Eqs 23e25 are referred to as the ‘full’ formulation, below.
11644
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 1 1 6 4 1 e1 1 6 5 0
(II) Low mass transfer rate In the PEMFC example given in ref. [10], it was noted that the blowing parameter was sufficiently small to linearise the 00 problem, g¼g*, m_ ¼ g* B. For low mass flow rates Eqs. 23 and 24 simplify to the following linearised form; In the liquid film 00
m_ ¼
h¼
00 m_ O2 yl ye ¼ gl yt ye yt
(26)
(a) If all liquid water remains in the liquid film yb yg yg 1
00
m_ yt ¼ gg
! (27)
(b) If all the water is transferred to the gas Eq. (25) reduces to 00
m_ ¼ gg
yb yg yg yt
to obtain the mass fraction at the electrode wall, ye, as a function of this and the bulk value, yb, by inverting the above system of equations. In general the (I) full, (II) linearised, and (III) simplified approaches will generate different numerical results for ye. In this paper we shall compare the results in terms of a concentration overpotential, h, defined here as
! (28)
Solution to Eqs 26e28 are referred to as the ‘linearised’ form, below. Let the non-dimensional blowing parameter, 00 b ¼ m_ g
(29)
clet be defined as the ratio of convection to diffusion i.e., a Pe number for mass transfer with bl and bg defined as mass flux normalised with respect to liquid and gas conductances, respectively.
Problem statement and solution methodology In many mass transfer problems in the chemical engineering sector, the rate of mass transfer at the boundary is unknown, and the system of equations above is used to obtain a solution to the problem. In fuel cells however, the mass transfer rate, 00 m_ , corresponding to some current density, i00 , may be considered as known (at least superficially), and it is required
RT ye H RT xe H0 ln ln ¼ 4F yb 4F xb
(30)
The calculation procedure for h(i00 ) is detailed in Table 1: For both (I) full and (II) linearised formulations; Given the 00 current density, i00 , the mass flux, m_ and gas blowing parameter, bg, are computed from Faraday's law, Eq. (6), and the definition of b, Eq. (29), respectively. This allows the gas driving force, Bg, and hence the gas-side oxygen mass fraction at the interface, yg, to be computed from known gas bulk value, yb, and conductance, gg . Precise mathematical details vary depending according to whether the linearised or fullform is assumed, and also whether it is presumed all water is transferred to (a) liquid or (b) gas bulk state. Full details of the method for the sequence of calculations inverting the mass transfer equations to obtain ye from yb for each case are detailed in Table 1. After obtaining the oxygen mass fraction on the liquid-side of liquidegas interface from Henry's law, Eq. (10), the above procedure is repeated for the liquid film, and the value at the electrode wall finally computed from the blowing parameter/driving force identity for the liquid film. For the simplified case (III) it is not necessary to proceed from the gas bulk to the electrode wall in a sequential fashion, since Eq. (8) may easily be inverted to write ye in terms of yb and g . Another parameter of interest is the so-called limiting 00 current density, ilim , which for cases (I) and (II) is obtained by 00 iterating on values of the mass flux, m_ , i.e. i00 , until an electrode-value, ye¼0 is obtained. This is achieved using any simple root-finding algorithm. For case (III) the simplified 00 method, no such iteration is needed and ilim is computed directly, as detailed in Table 1. It is, of course, straightforward to write analytical expressions for ye(bl,bg), for cases (I) and (II), though these are non-linear expressions, and inversion is therefore most easily accomplished numerically.
00
Table 1 e Methods used to compute ye and ilim for all water transferred to (a) liquid state, (b) gas state. Method I. Full: Eqs. 23e25
(a) (b)
II. Linearised: Eqs. 26e28
(a) (b)
III. Simplified: Eq. (8)
00
ye
ilim
Bg ¼ exp(ytbg1) yg ¼ (ybþBg)/(1þBg) Bg ¼ exp(bg1) yg ¼ (ybþBgyt)/(1þBg) yl ¼ yg/H Bl ¼ exp(bl1) ye ¼ (ylþBlyt)/(1þBl) Bg ¼ ytbg yg ¼ (ybþBg)/(1þBg) Bg ¼ b g yg ¼ (ybþBgyt)/(1þBg) yl ¼ yg/H Bl ¼ b l ye ¼ (ylþBlyt)/(1þBl) 00 ye ¼ yb =H þ m_ =g
Iterate on values of i00 until ye ¼ 0
Iterate on values of i00 until ye ¼ 0
ilim ¼ 4000Fg yb =MO2 H 00
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 1 1 6 4 1 e1 1 6 5 0
11645
Table 2 e Prescribed values, baseline case. Temperature, T Pressure, p O2 gas mass fraction Net drag coefficient, a Gas conductance, gg Henry constant, kH Liquid density, rl Oxygen liquid diffusivity, Dl Liquid film thickness, dl
353 K 1 Atm 0.232(21%vol.)0.116(10.5%vol.) 0.2 0.1 kg(m2 s) 0.2 Atm m3/mol 9.718 kg/m3 108m2/s 109m
Results Fig. 6 e Limiting current density as a function of the conductance (a) no water transferred to gas.
Table 2 show values prescribed to perform sample mass transfer calculations. For the gas phase, yb is taken as the bulk value in the air channel, Fig. 1, and the gas mass transfer coefficient, gg , is estimated as the harmonic average of the combined channel and gas diffusion layer values from Ref. [11]. Two bulk values were considered for the air corresponding to 21% vol. oxygen (i.e., atmospheric air) and 10.5% vol. (i.e., 50% depleted air). Sample calculations were performed to investigate the sensitivity of the results to the choice of equations. Fig. 5 show concentration overpotential, defined by Eq. (30) vs. current density, i00 , for (a) no gaseous water transfer. The (I) high mass rate formulation is compared to (II) the low mass rate linearised formulation as well as (III) the simplified theory. It can be seen that a limiting current density is reached for yb ¼ 0.116 at (I) 1.39105 A/m2, (II) 1.32105 A/m2 and (III) 1.31105 A/m2. For yb ¼ 0.232 the values are (I) 2.97105 A/m2, (II) 2.65105 A/m2, and (III) 2.61105 A/m2. Thus, for the set of parameters chosen in Table 1, there is only a difference of 00 about 6e12% in ilim , with the full approach (I) generating a 00 maximum value of ilim and the simplified method (III) a minimum value. The Henry number has a value of 6.7103 corresponding to kH¼0.2Atm,m3/mol. The corresponding value of gg H=gl is 0.069. Fig. 6 shows limiting current density computed as a function of the non-dimensional parameter gg H=gl on the basis that (a) all produced water remains in a liquid state, with KH¼0.2 Atm , m3/mol. This figure was generated by varying
the ratio Dl/dl, i.e. gl , over a range of values reported in the literature. It can be seen that the limiting current increases with bulk mass fraction of oxygen, and decreases as gg H=gl increases. Fig. 7 shows results similar to Fig. 5, this time on the assumption that (b) all the produced water enters the gas phase. It can be seen that for yb¼0.116 a limiting current density is reached, at values of (I) 1.27105 A/m2, (II) 1.30105 A/m2, and (III) 1.31105 A/m2. For yb¼0.232 a limiting current density is reached, at values of (I) 2.46105 A/m2, (II) 2.60105 A/m2, and (III) 2.61105 A/m2, with the full approach (I) 00 generating the minimum ilim and the simplified theory (III) about 3e6% higher. Fig. 8 shows the limiting current density on the premise (b) that all the produced water enters the gas phase. It can be seen that the detailed mass transfer formulation predicts limiting current densities that are lower than the simple equivalent film theory in contrast to the results of Fig. 6, where the converse is true. Fig. 9 shows the limiting current density as a function of the net drag coefficient, a, for case (b) produced water transferred to gas. All other parameters are held to the baseline values given in Table 2. It can be seen that as a is increased from 0.5 to þ1.0, the full model (I) indicates a
Fig. 5 e Concentration overpotential for the baseline case (a) No water transferred to gas.
Fig. 7 e Concentration overpotential for the baseline case (b) All water transferred to gas.
11646
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 1 1 6 4 1 e1 1 6 5 0
Fig. 8 e Limiting current density as a function of the conductance (b) all water transferred to gas.
substantial reduction in limiting current from 3.0105 A/m2 to 2.1105 A/m2 compared to (II) the much smaller decline from 2.6105 A/m2 to 2.5105 A/m2 for the linearised model at yb¼0.232. The (III) simplified model does not account for changes in limiting current as a function of water content and therefore displays a horizontal slope. The three profiles coin00 cide precisely at a value of ilim ¼ 2:6 105 A=m2 , with a¼1/18 corresponding to yt¼±∞, Fig. 4. For case (a) not shown, there is 00 no significant change in ilim as a function of a. The matter is discussed further below.
Discussion Inspection of Figs. 5 and 7 suggests that the results of the simple equivalent film theory (III) agree well with the low mass formulation (II), and quite well with the high mass-rate formulation (I) when the oxygen mass fraction is dilute, yb¼0.116. This is consistent with ref. [2] which indicates that for yb¼0.1, the simplified theory (III) generates reasonable results for oxygen consumption. However, as the oxygen bulk mass fraction is increased, e.g. to yb¼0.232, the errors increase, and will do so further in the event that enriched air were employed as oxidant.
Fig. 9 e Limiting current density as a function of the net transfer coefficient, a. (b) All water transferred to gas.
Within the liquid film, oxygen mass fractions are always extremely small with H¼6.7103. This means the maximum ye is 1.7105 for yb¼0.116 and 2.4105 for yb¼0.232 for i00 ¼0. Clearly the condition jyej≪jytj is met for the liquid film. For the choice of parameters in Table 2, the gas controls mass transfer more than the liquid, as gg H=gl is 0.069 i.e., the gas mass transfer resistance is larger. In fact the value of gg H is only about 7% lower than the value of g computed by Eq. (9). With (a) no water being transferred into the gas phase, Fig. 5, at i00 ¼105A/m2, Bg and Bl are of the order of 0.08 and þ5107, respectively i.e., Bg is quite small and Bl is tiny. (The sign of B indicates that the net mass transfer is out of the gas but into the liquid.) The values are the same at both yb¼0.116 and 0.232. For i00 ¼2105A/m2, Bg and Bl are of the order of 0.15 and 106. Thus the simplified theory (III) approximates the low mass transfer linearised method (II) quite well, with some minor differences between these and the full formulation (I), due to the fact that even with bg small, values of Bg¼exp(bg)1 and Bg¼bg are slightly different, with the difference increasing 00 00 as i /ilim . In other words, yg is also quite small compared to unity in (though not as small as ye), Eq. (17). For this reason, there are relatively minor differences in the limiting current density computed by the three methods. Increasing the oxygen bulk value, the differences start to become more significant, as might be expected. With (b) all the water being transferred into the gas phase, Fig. 7, values of the liquid driving force, Bl, are similar to those for the results of Fig. 5, but the gas values differ both in magnitude and sign with Bg around þ0.05 and þ 0.1 for i00 ¼105A/m2 and 2105A/m2,respectively. This change of sign is due to the fact that the while direction of the net mass flux of water and oxygen is (a) out of the gas for the results of Fig. 5, it is into the gas for the latter case (b), due to the addition of the water vapour to the air. Comparison of Figs. 5 and 7 reveal that whereas the simple theory (III), predicts the same rate of mass transfer when either (a) all water remains in the film or (b) all water is transferred to the gas; both the high-mass (I), and linearised low-mass (II) formulations suggest that higher limiting current densities and lower concentration overpotential magnitudes occur when the water remains in the liquid rather than the vapour state. Thus a comparison of method (I) for limiting current densities of (a) 1.39105 A/m2 and 2.97105 A/m2 with (b) 1.27105 A/m2 and 2.46105 A/m2 at yb¼0.116 and 0.32 suggests a 10e21% difference between cases (a) and (b) for this performance measure. Inspection of Figs. 6 and 8 reveals the limiting current is reduced as values of the parameter gg H=gl increase. There is 00 some difference between values of ilim predicted by the highmass theory and the others, notably for small gg H=gl ; however the linearised and simplified agreements are by-and-large, in agreement. The oxygen concentration overpotential, defined by Eq. (30) is a measure of decrease in potential due to mass transfer effects, that is the difference in Nernst potential associated with the electrode mass fraction being ye rather than yb/H. Because of this, the ordinates in Figs. 5 and 7 are proportional to a logarithmic quantity, ln(yeH/yb). This has the effect of compressing the results, thereby minimising differences. If instead the change in current density as opposed to potential
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 1 1 6 4 1 e1 1 6 5 0
is considered; and if it is assumed cathode reaction occurs in the Tafel zone. Then at any given electrode potential, the current density will be reduced by the exponential of (the change in) the potential compared to the value if the electrode were at yb/H, i.e., the reduction in current density is proportional to yeH/yb and not ln(yeH/yb). The current density, and mass flow rate are nominally proportional to the mass transfer driving force, and it can be shown that for a given cell voltage, the current density may be written, . 00 00 (31) i ib ¼ expðb=2Þ ye H yb 00
00
where b is a charge transfer coefficient, and ib denotes i when ye¼yb/H'. Note that Eq. (31) is an idealisation, since the exchange current density will also be a function of oxygen concentration, further complicating the relationship. There will also be a small water concentration-overpotential; since the water mass fraction, will also differ. This should not be confused with the problem of liquid-water blocking oxygen transport in porous media [10]. The ratio of the parameters Dl and dl is important, and clearly affects the value of gl in Eq. (12). Various values may be found in the literature; dl-values ranging from 2 nm to 80 nm have been recorded [12,13]. Oxygen liquid-diffusivity correlations also vary significantly, with published values of Dl in the range of 1.21010m2/s to 2108m2/s [7,9,14]. Variations in Henry constants are smaller, but by no means insignificant. The parameter gl could be regarded simply as a fitted parameter, to be obtained from limiting current density data for a given cathode, assuming gg is known. As the parameter gg H=gl approaches zero, the mass transfer calculation can be performed based on the gas-side alone [3]. If all that is of interest is the oxygen concentration overpotential, it may be argued that one may adopt the Henry's law-modified simplified theory (III) as a reasonable first approximation to predict oxygen mass transfer losses in the cathode of a PEMFC. The oxygen concentration overpotential is readily obtained [2] from Eq. (8) in the well-known form, 00
RT i ln 1 00 h¼ 4F ilim
! (32)
where ilim ¼ 4000Fg yb =MO2 H as given in Table 1. The reader 00 will note that ilim , as predicted by method (III), is independent of yt and hence a. The same is not true for methods (I) and (II) 00 where it can be seen that ilim decreases as a is increased. Method (III) is a fair approximation under the present circumstance, because the driving forces are small, and the magnitude of yt is for the most-part substantially greater than that of yg, and thus Eq. (7) is therefore a reasonable approximation to the full equations. Simplicity in a mathematical model is a good thing, however, when it comes to the calculation for the rate of mass transfer of water, Eq. (7) is not suitable except when the water mass fractions are very small. Thus if Eq. (7) is used to compute or fit a value of the overall mass transfer coefficient, simply substituting the value of g into the conductance in Eq. (1) for mass transfer of water is illadvised. The matter is the subject of a future parametric study. Thus the simplified method does not provide a rationale foundation to obtain all the mass fractions of the component gases simultaneously, and a basic understanding 00
11647
of the full analysis is important. This is corroborated by the fact that many authors, for example [15], have noted problems with use of the term RT/4F in Eq. (32) and suggest instead the 00 00 use of a fitted-constant, C, with h ¼ Cln 1 i =ilim Þ. This has no theoretical basis in the field of mass transfer, and is merely a reflection that the logic underlying methodology (III) is too simplified to be universally applicable, except for dilute mixtures [2]. The results of Fig. 9 show clearly that the impact of net water drag coefficient, a, upon performance cannot be predicted by the simplified approach (III). Here too some significant differences between the (I) full and (II) linearised models are apparent. A diverse range of a values have been reported in the literature, with some authors claiming it to be a function of current density and others not [6]. The three methodologies predict identical limiting currents for a¼1/18 corresponding to yt¼±∞, when there is no net mass transfer with the mass flux of oxygen being precisely equal and opposite to that of the water. As the magnitude of a deviates from this value, so does agreement between the models for 00 ilim . Because the magnitude of Bl is small in comparison to Bg, the effect is important for case (b) all water transferred to gas, but not case (a) all water remains as liquid. In most circumstances, it is a requirement to obtain the mass fractions of all the species; oxygen, water, nitrogen, etc. as part of global or local calculations, where activation and ohmic terms will also affect the operating voltage/current 00 density which in turn affects the value of i00 and hence m_ . If the diffusion coefficients of the gases are all nominally equal, the full/linearised methodology provides a simple method for computing these values since BO2 ¼ BH2 O ¼ BN2
(33)
allows each gas-phase mass fraction to be computed from the mass flow rate, bulk and T-state values. This type of engineering approximation was frequently made for atmospheric air and other gases, prior to the advent of detailed computational models. When the binary diffusivities of the gas component-pairs differ considerably, as might be the case for mixtures of gases with hydrogen, there is some cause for concern over the applicability of this type of model, per se. In PEMFCs often only a binary mixture of humidified hydrogen is employed on the fuel side, though in other fuel cells and electrolysers, situations could arise where there is a need to employ higher-order numerical schemes [16] based on MaxwelleStefan equations or generalised Fick's law [17]. On the other hand, it is to be noted that many mass transfer Stanton number correlations (used to calculate the conductance g*) exhibit a relatively weak dependence upon the Schmidt number, the ratio of diffusivity and kinematic velocity. In any event, there is a role for both analytical [18] and detailed 3-D models in the development of fuel cell technology. For a lumped-capacity analysis as typically found in 0-D and 1-D analyses [19], the choice of yb, corresponding to the bulk value in the channel passage as reference, is appropriate. The associated length-scale is the sum of the size of the diffusion layer and half-hydraulic diameter of the channel. In a fuel cell operating at high air utilisation, there will be a significant decrease in oxygen mass fraction from inlet to outlet, which must be accounted-for: In a lumped capacity model, the bulk value should be a mean value of the inlet and
11648
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 1 1 6 4 1 e1 1 6 5 0
outlet mass fractions, the latter being computed from a mass balance for oxygen and water. This is a function of the total current and the partial densities, ri¼ryi, of oxygen, water (and nitrogen). Therefore calculations on the water mass fraction are important. The same is true for simple multi-dimensional ‘along the channel’ models, where the bulk mass fractions are adjusted cell-by-cell. For all these reasons, the (I) full or (II) linearised mass transfer formulations based on general theory, as outlined in this paper are recommended for consideration when performing calculations on mass transfer in all but the simplest PEMFC problems. In 3-D computational fluid dynamics (CFD) calculations, the fuel cell geometry is typically tessellated with a computational mesh, and the species equations solved-for in the channels, diffusion layers and electrodes together with the equations for mass, momentum, and energy [19]. The mass transfer formulation, described here, may be considered a sub-scale closure within the electrode and yb chosen as the local bulk gas value within the electrode pores, as opposed to within the gas channel, Fig. 2. This local bulk value is not constant but changes as a function of position, local current density, etc. The gas length scale is therefore of the order of the pore size, and the gas conductance greatly increased, accordingly. Thus gg H=gl [0, i.e., the liquid-film is the ‘controlling’ or rate-limiting factor, however the CFD code should at the very least account for solubility effects on the Nernst equation and current density and if possible be capable of calculating bulk-to-interface driving forces, although depending on the pore size these may be small. As discussed above, an estimate of the overall mass transfer coefficient may be made from the limiting current density, and vice versa. However it is only an estimate. Moreover Henry's law is to be considered as an approximation at finite rates-of mass transfer [3,20]. A major assumption of the mass transfer model presented here, is the presence of a simple liquid film of effective height, dl, since the height will almost certainly vary with current density and protonic drag. In this paper the limiting cases of all produced water, either (a) remaining in the liquid film or (b) entering the gas phase were considered. In practice this needs to be computed either locally in a CFD or presumed-flow type code, or globally in a lumped-parameter analysis. Although this is considered beyond the scope of this paper, some brief discussion is in order. The situation is more typical of conventional mass transfer problems, where the rate of mass transfer is not known, a priori, but rather is the required dependent variable, and may be classified as a ‘mass-from-heat transfer’ problem for which there have been many solved examples in chemical and nuclear engineering over the past four decades [21,22]. Let the mass fraction of water vapour be denoted by yv. Under the assumptions that (i) the latent heat of vapourisation, Hfg, is constant, (ii) the vapour specific-heat is much larger than that of the liquid, and (iii) the gas is perfect, the Clausius-Clapeyron relationship may be stated in the form, pv ðTÞ ¼ C exp Hfg RT
(34)
Empirical correlations such as the August-Roche-Magnus equation, are readily employed in place of Eq. (34) to formulate an explicit expression for pv as a function of T. If the vapour partial-pressure, pv, is known at the interface then the
temperature T(pv) is also known, and vice-versa. Once the interfacial temperature is fixed, as a Dirichlet condition, the unknown liquid-to-vapour mass flux can be obtained from the enthalpy equation, and local bulk gas and liquid temperatures. The combined flux of both dissolved oxygen, and water from liquid to vapour, is a function of the T-state value which is dependent on the individual mass fluxes (of oxygen and water). Thus, in addition to Henry's law, Eq. (34) gives the unknown water-steam mass flux, and hence the correct transferred substance states in the gas mixture 00
m_ ¼
hg Tb;g T þ hl ðTb;l TÞ Hfg
(35)
where hl and hg are liquid-side and gas-side heat transfer coefficients, Tb,l and Tb,g are bulk values in the two phases, and Hfg is evaluated at the interface temperature, T. This would normally be done using a sub-scale closure in a CFD-based calculation, e.g. in the gas diffusion layer and/or electrodes. An alternative approach [23] is to compute the mass fraction and mass transfer rate from the partial vapour pressure, pv(T), given the interface temperature via Eq. (34). The associated inter-phase heat flux may then be prescribed as a Neumann condition in the energy equation. In multi-phase CFD codes it is not generally necessary nor desirable to presume thermal equilibrium between phases, since the phase temperatures can differ, and moreover if there is no heat transfer, then there can be no mass transfer by evaporation/condensation and vice versa. Regarding equilibrium of species, it is generally true that the partial vapour pressure will only deviate by a small amount from the equilibrium value for finite rates of mass transfer [3]. In catalyst layers, where the liquid film is not saturated water, many authors use the closure due to Springer et al. [24] where the dependence of the liquid concentration on pv is prescribed empirically as a function of temperature and membrane composition (humidification). This is a natural extension of the approach above and may be considered as a similar, but more complex case than can be described by idealised solution theory (Raoult's law). The additional water mass flux must therefore be accounted for by modifying the mass flux and transferred substance-states from the two limiting cases prescribed in Eqs. (14) and (19) to additionally account for the fraction of water which is evaporating (or condensing). The mass transfer analysis described in this paper may also be used to consider heat transfer, by analogy; Driving force(s) and blowing parameter(s) may all be similarly defined in terms of temperature or enthalpy. The challenge facing the engineer in this context is to obtain reliable correlations for phase heat and mass transfer coefficients for the complex geometries and flow regimes frequently encountered in PEMFC designs.
Conclusions A standard mechanical-engineering type analysis was developed and applied to an idealised PEMFC cathode. Both high and low mass transfer models were derived and compared to a highly-simplified method. It was shown that if the oxidant is composed of atmospheric air with a maximum value of 21%
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 1 1 6 4 1 e1 1 6 5 0
vol., the simplified method gives a reasonable first estimate of the limiting current density. The simplified method agrees well with the (linearised) low mass transfer method, and has the advantage that the limiting current density may be computed directly and simply. If, however, it is required that the gas mass fractions of all species be computed simultaneously, these should be derived from the driving force approach-based detailed analysis presented in this manuscript. Moreover, if enriched oxygen is employed, the standard approach is a prerequisite requirement for a reliable solution. It was shown that the limiting current density is reduced by up to 21% if water is transferred to the gaseous state, rather than remaining in entirely in the liquid state; something which cannot be predicted by the simplified model. The ratio of the gas and liquid conductances, gg H=gl , governs the degree to which the gas or liquid controls mass transfer, with limiting current density decreasing as gg H=gl increase. This value was quite small for the particular baseline parameters prescribed in this study.
Acknowledgement Prof. Werner Lehnert and Dr. Uwe Reimer both read the manuscript, and contributed useful commentary. I would also like to thank Dr. Michael Malin for remarks on the phase change problem and especially Prof. Brian Spalding for personally sharing his unique insight into heat and mass transfer, which directly influenced the subject matter of this paper. The support of the Forschungszentrum Ju¨lich is gratefully acknowledged.
Nomenclature English Letters A area, m2 b fitted constant, V c concentration, mol/m3 D diffusion coefficient, m2/s F Faraday's constant, 96485.3 C/mol g mass transfer coefficient, kg/(m2,s) 00 mass transfer conductance in the lim m_ /0, kg/ g* 2 (m ,s) h heat transfer coefficient, W/(m2,K) latent heat of vapourisation, J/kg Hfg k molar conductance, mol/(m2,s) 00 molar conductance in the lim n_ /0, mol/(m2,s) k Henry's constant, (Pa,m3)/mol kH 00 mass flux, kg/(m2,s) m_ M molecular weight, g/mol current density, A/m2 i00 00 molar flux, mol/(m2,s) n_ electro-osmotic drag coefficient nd p pressure, Pa partial vapour pressure, Pa pv R gas constant, 8.3144 J/(K,mol) t time, s
T x y
11649
temperature, K molar fraction mass fraction
Greek Letters a net protonic drag coefficient b charge transfer coefficient dl liquid film thickness, m h concentration overpotential, V r mixture density, kg/m3 t wall shear stress, Pa Non-dimensional number b blowing parameter B driving force H Henry number molar Henry number H0 Subscripts b bulk value e electrode wall value g gas, gas-side interface value hydrogen H2 H2O water i species i l liquid, liquid-side interface value lim limiting value nitrogen N2 oxygen O2 t transferred substance state value Superscripts overall value
references
[1] Bird RB, Stewart WE, Lightfoot EN. Transport phenomena. New York: Wiley; 1960. [2] Beale SB. Calculation procedure for mass transfer in fuel cells. J Power Sources 2004;128(2):185e92. [3] Spalding DB. Convective mass transfer: an introduction. Edward Arnold; 1963. [4] Kays WM, Crawford ME, Weigand B. Convective heat and mass transfer. 4th ed. New York: McGraw-Hill; 2005. [5] Spalding DB. A standard formulation of the steady convective mass transfer problem. Int J Heat Mass Transf 1960;1:192e207. [6] Dai W, Wang H, Yuan X-Z, Martin JJ, Yang D, Qiao J, et al. A review on water balance in the membrane electrode assembly of proton exchange membrane fuel cells. Int J Hydrog Energy 2009;34(23):9461e78. [7] Bernardi DM, Verbrugge MW. Mathematical model of a gas diffusion electrode bonded to a polymer electrolyte. AIChE J 1991;37(8):1151e63. [8] Parthasarathy A, Srinivasan S, Appleby AJ, Martin CR. Pressure dependence of the oxygen reduction reaction at the platinum microelectrode/Nafion interface: electrode kinetics and mass transport. J Electrochem Soc 1992;139(10):2856e62. [9] Um S, Wang C-Y, Chen KS. Computational fluid dynamics modelling of proton exchange membrane fuel cells. J Electrochem Soc 2000;147:4485e93.
11650
i n t e r n a t i o n a l j o u r n a l o f h y d r o g e n e n e r g y 4 0 ( 2 0 1 5 ) 1 1 6 4 1 e1 1 6 5 0
[10] Beale SB, Schwarz DH, Malin MR, Spalding DB. Two-phase flow and mass transfer within the diffusion layer of a polymer electrolyte membrane fuel cell. Comput Therm Sci 2009;1(2):105e20. [11] Beale SB. Conjugate mass transfer in gas channels and diffusion layers of fuel cells. ASME J Fuel Cell Sci Technol 2007;4:1e10. [12] Harvey D, Pharoah JG, Karan K. A comparison of different approaches to modelling the PEMFC catalyst layer. J Power Sources 2008;179(1):209e19. [13] Sadeghi E, Putz A, Eikerling M. Hierarchical model of reaction rate distributions and effectiveness factors in catalyst layers of polymer electrolyte fuel cells. J Electrochem Soc 2013;160(10):F1159e69. [14] Parathasarthy A, Srinivasan S, Appleby AJ, Martin CR. Temperature dependence of the electrode kinetics of oxygen reduction at the platinum/Nafion interface - a miroelectrode investigation. J Electrochem Soc 1992;139. 2530e2357. [15] Larminie J, Dicks A. Fuel cell systems explained. 2nd ed. New York: Wiley; 2003. nchez-Insa A, Fueyo N. Multimodal [16] Garcı´a-Camprubı´ M, Sa mass transfer in solid-oxide fuel-cells. Chem Eng Sci 2010;65(5):1668e77. [17] Taylor R, Krishna R. Multicomponent mass transfer. John Wiley & Sons; 1993.
[18] Kulikovsky AA. Analytical modelling of fuel cells. Elsevier; 2010. [19] Weber AZ, Borup RL, Darling RM, Das PK, Dursch TJ, Gu W, et al. A critical review of modeling transport phenomena in polymer-electrolyte fuel cells. J Electrochem Soc 2014;161(12):F1254e99. [20] Perry RH, Green DW, Maloney JO. Perry's chemical engineers' handbook. New York: McGraw-Hill; 1997. [21] Spalding DB. Numerical computation of multi-phase fluid flow and heat transfer. In: Taylor C, Morgan K, editors. Recent advances in numerical methods in fluids, vol. 1. Pineridge Press; 1980. p. 139e67. Ch. 5. [22] Spalding DB. Mathematical methods in nuclear-reactor thermal hydraulics. In: Lahey R, editor. Proceedings of American nuclear society, meeting on nuclear-reactor thermal hydraulics. Saratoga, New York: American Nuclear Society; 1980. p. 1979e2023. [23] Berning T, Djilali N. A 3D, multiphase, multicomponent model of the cathode and anode of a PEM fuel cell. J Electrochem Soc 2003;150(12):A1589e98. [24] Springer TE, Zawodzinski TA, Gottesfeld S. Polymer electrolyte fuel cell model. J Electrochem Soc 1991;138(8):2334e42.