Abstract
An intrinsic time-dependent one-dimensional (1D) model and a macro two-dimensional (2D) model for the anode of the direct methanol fuel cell (DMFC) are presented. The two models are based on the dual-site mechanism, which includes the coverage of intermediate species of methanol, OH, and CO (, , and ) on the surface of Pt and Ru. The intrinsic 1D model focused on the analysis of the effects of operating temperature, methanol concentration, and overpotential on the transient response. The macro 2D model emphasises the dimensionless distributions of methanol concentration, overpotential and current density in the catalyst layer which were affected by physical parameters such as thickness, specific area, and operating conditions such as temperature, bulk methanol concentration, and overpotential. The models were developed and solved in the PDEs module of COMSOL Multiphysics, giving good agreement with experimental data. The dimensionless distributions of methanol concentration, overpotential, and current density and the efficiency factor were calculated quantitatively. The models can be used to give accurate simulations for the polarisations of methanol fuel cell.
1. Introduction
Compared with hydrogen fuel cells, direct methanol fuel cells (DMFCs) have significant advantages such as higher efficiency, easier design and operation, simple storage, and refilling of liquid fuel. However, DMFCs show serious disadvantages such as lower current density, larger polarisation, and lower limiting currents resulting from more complex kinetics of methanol oxidation in the anode of DMFCs than that for hydrogen fuel cells [1–3].
A reason for the disadvantages of DMFCs was the adsorption coverage of intermediates, which accompany the kinetics of methanol oxidation, and play a crucial function in the behavior of DMFCs [4–8]. The micro structure of the anode, which results in different transport coefficients for various species, also determines the performance of DMFCs to a certain extent [9–11].
Meyers and Newman [12] developed a comprehensive model which describes the thermodynamics, transport phenomena, and electrode kinetics of the system in DMFCs. Kauranen and Skou [13] indicated that saturation of OH on surface of Pt leads to limiting current when oxidation of methanol was catalyzed by pure platinum. Nordlund and Lindbergh [14] adopted a porous model with agglomerates and kinetic equations based on surface coverage to study the influence of anode structure on the anode of DMFCs. Argyropoulos, Scott et al. [7, 8, 15] investigated the anode reaction with Pt-Ru catalyst layers, based on surface coverage, to produce simple current-voltage expressions, which were solved analytically, and they showed that saturation of methanol on the Pt surface and OH on the Ru surface were the reasons for limiting currents occurring for the aspect of intrinsic kinetics.
In general, describing electrode performance by Tafel or dual-site kinetics, the dynamic equations were always nonlinear because of the exponential functions of potential. Since Newman [16] developed BAND(J) in order to solve such nonlinear equations, software such as COMSOL, FLUENT has been adopted to model the DMFCs [17–19]. Thus, models became more complex and more detailed. Modeling the catalyst layer, as the core area of anode, was considered as the first priority since the electrochemical reaction determines its performance and thus that overall of the DMFCs [20].
The present model aims to combine the dual-site methanol oxidation mechanism with material and charge balances in order to simulate time-dependent current-voltage response and predict the anode behavior under different operating temperatures and concentrations, and to especially focus on the distributions of methanol concentration, overpotential, and current density inside of the catalyst layer.
2. Mathematical Model
Figure 1 presents a schematic view of the modeling domain of a porous anode with Pt-Ru catalysis for a typical DMFC unit. Methanol was fed into the anode channel, and it diffuses into catalyst layer through the diffusion layer. The specifications and relative parameters of the anode are shown in Table 1.
The assumptions adopted in the present model were as follows.(1)Methanol concentration was defined as constant at the interface of catalyst layer () along the -direction which was adjacent to diffusion layer. This meant that at the interface of catalyst layer, the methanol concentration was the same as the bulk concentration in the anode liquid channel. (2)Carbon dioxide bubbles were formed beyond the catalyst layer. It was possible for nucleation of carbon dioxide to take place in diffusion layer or limited to a partial region of the catalyst layer by choosing appropriate operating condition [21, 22].(3)The porous catalyst layer was assumed to be isothermal, isotropic, and homogeneous.(4)Steady state for mass and charge transport was assumed in the porous catalyst layer for the 2D model.
2.1. Intrinsic Kinetic Equations
According to dual-site kinetics which was widely accepted for methanol absorption and electrochemical oxidation on the surface of Pt-Ru catalyst, the methanol oxidation mechanism can be described as the following four elemental steps below:
It was assumed that steps (1), (2), and (4) occurred on platinum sites (Pt), step (3) occurred on ruthenium sites (Ru), and step (5) was catalyzed on the sites of both Pt and Ru. Generally, it was believed that the reaction of COads to CO2 occurs on Ru, and Pt serves as an active surface of adsorption and dehydrogenation of methanol [4, 5, 13, 23]. As a result, the rate controlling step was reaction (5), which in turn depended on elemental steps (1), (2), (3), and (4) for the formation of adsorbed intermediates. Thus, the rate expression of the overall reaction can be written as where, .
The rates of changes of surface coverage of different intermediates with respect to time were as follows:
Numerous papers indicate that OHads was preferentially formed on the surface of Ru, not on the surface of Pt [4, 23–25]. Our previous modeling [10, 11] confirmed that was almost zero when the overpotential was lower than 0.5 V (versus dynamic hydrogen reference electrode (DHE)). As a result, it was reasonable to assume adsorption of hydroxyl ions on Pt sites could be neglected (). Thus, if the water activity can be defined as unity (), (6) was simplified to where and can be obtained from the result of solving (7)–(9).
Intrinsic kinetic current density can be obtained by combining Faraday’s Law: where the exchange current density .
2.2. Macro Kinetic Equations
2.2.1. Mass Transport Balance
Mass transport of methanol in porous catalyst layer can be described by Fick’s first law Because the transport process of methanol in a differential volume of catalyst layer was only described by diffusion, the divergence of was written as:
According to a mass balance [9–11, 20] Combining (14) and (15), we have Equation (16) was used to describe the effect of concentration changes on polarisation of anode. The boundary conditions for the second-order differential equation above were
2.2.2. Charge Transport Balance
Charge transport of methanol in the porous catalyst layer can be described by Ohm’s law: The divergence of was written as According to a charge balance [9–11, 20] The overpotential was written as [9–11, 14]. Both and could be considered as constants. Substitute (19) into (20), the differential charge balance which described the potential field in the porous anode becomes a nonlinear Poisson equation: Equation (21) was used to describe the effect of ionic resistance on anode polarisation.
The boundary conditions for the second-order differential equation above were
2.2.3. Generalization of the Model Equations
Equations (16) and (21) could be generalized by dimensionless variables such as and . As a result, the dimensionless equations were presented as below: For a two-dimensional porous electrode [26], (23) become: with boundary conditions where , and .
The dimensionless modulus, , in (24) characterizes resistances of mass transport in the porous anode and the dimensionless modulus, , in (25) characterizes the relative resistance of charge transport when applying different overpotentials at the boundary .
2.2.4. Current Density and Effectiveness Factor
According to Ohm’s law, the local current density described by concentration and charge flux was
Thus, the dimensionless current density was defined as Equations (28) and (29) can be used to describe the current distribution in the anode catalyst layer. Then the relationship between and can be obtained by the expression of and , We can derive the total current density as This was the equation for predicting the macro current density of the anode. Here the relationship of and typically describe the apparent and dimensionless polarisation curve.
Effectiveness factor was introduced to evaluate the impact of physical parameters such as thickness and specific surface area of the catalyst layer, effective diffusion coefficient, and effective conductivity of the anode, and was defined as, which was calculated by solving (7)–(9), and was only the function of initial methanol concentration and overpotencial and .
Substituting (12) into (32), we obtain
Hence, the expressions of and were applied to calculate the effectiveness of the Pt-Ru catalyst layer. Assume that the thickness and width of the catalyst layer would be of the same order of magnitude. The modeling domain was chosen as a rectangle with the width five-times larger than the thickness.
2.3. Numerical Method
1D and 2D models were established and solved in the PDE module of COMSOL Multiphysics software in order to analysis firstly the change of surface coverage of different intermediates (, , and ) and time-dependant voltage-current density relationship which were described by (7)–(9). Secondly, the dimensionless distributions of methanol concentration, overpotential, and current density, described by (24)–(29), and thirdly, effectiveness factor, which was described by (33), were determined.
Equations (7)–(9) were solved by the time dependent solver of COMSOL in the 1D model, time stepping was from 0 to 2000 seconds with an interval of 1 second. In the 2D model, coupled equations (7)–(9), (24) and (25) were solved by the nonlinear stationary solver. The inlet boundary () was considered as Dirichlet boundary (the first kind of boundary condition) where the initial methanol concentration and overpotential were defined as one (), the outlet boundary () was defined as Neuman-boundaries (the second kind of boundary condition) where the first order derivative of methanol concentration and overpotential were defined as zero () and the other boundaries ( and ) were defined as Dirichlet boundary where the initial methanol concentration and overpotential were defined as zero () in the present 2D model.
It was worthwhile to note that the macro kinetic model was established by coupling (7)–(9), (24) and (25). These equations were included in the same domain of modeling and solved at the same time, which did not need to be decoupled. In order to make sure the accuracy of calculation, 1285 nodes in 1D model and 12544 grid nodes in 2D model were created in the finite element mesh of COMSOL, where the process of calculation was undertaken step by step along the nodes from to .
The parameters used in the present model were listed in Table 1 and Table 2.
3. Result and Discussion
3.1. Model Validation
The simulated polarisation curves of the DMFC anode at different operating temperature and methanol concentration are shown in Figure 2. The kinetic parameters were selected from the paper of Shivhare et al. [7]. As shown in Figure 2, the simulated polarisation curves were obtained by using two groups of thicknesses and specific areas of the catalyst layer selected from literatures [7, 10, 11, 15, 21–23, 27].
Figure 2 indicates that the simulated polarisation curves show good agreement with the experimental data at lower bulk concentration, lower operating temperature, and lower overpotential. Differences between the simulated values and experimental data become more apparent at higher bulk concentration, higher temperatures, and higher overpotentials, with the simulated values being higher than experimental values. The simulations over-predict the current densities of the fuel cell at the higher temperature, and methanol concentration in comparison to experimental data. This may be due to liberation of carbon dioxide bubbles, which were generated by the increased electrochemical reaction and temperature (lower solubility), limiting methanol transfer from the bulk into the porous catalyst layer.
3.2. Time-Dependant Performance
The influence of different operating temperature, methanol concentration, and overpotential on the transient surface coverage of methanol is shown in Figure 3. Elemental step (1) happens on the surface of Pt site and thus results in the typical change of from zero to saturation. In Figure 3, it is clear that took a longer period of time to approach saturation at lower operating temperatures, methanol concentrations, and overpotentials (see the dot lines). However, reached the saturation point faster when higher operating temperatures, methanol concentrations, and overpotentials were applied (see the solid lines). The comparison of Figures 3(a), 3(b), and 3(c) indicates that the operating temperature has a greater influence on compared to the other two parameters, methanol concentration, and overpotential.
(a)
(b)
(c)
Figure 4 shows the transient variation in current density at different operating conditions. Typically the steady state is approached within 200–400 seconds. The time for the transient current density to approach a steady state decreased when temperature, methanol concentration and overpotential increased. The operating temperature affects the current density more than the methanol concentration, and overpotential. The influence on was similar to its influence on . Moreover, the curves of and (not shown here) have great similarities in respect to the shape and the time variation, although was some four to eight orders of magnitude smaller than . This was because elemental step (3), which happens on the surface of Ru site, was extremely fast and resulted in the very fast and high saturation (very close to one) of on the surface of Ru site [11].
(a)
(b)
(c)
As stated in (12), current density was a function of the product of and . The current density would have no relationship to if could be considered as constant (the value is one). We may reasonably conclude that step (3) was not the rate determining step; thus, the whole speed of the electrochemical reaction was determined by step (2). In addition, step (1) might be the rate determining step because it provides the reactant to step (2).
3.3. Distributions in the Porous Catalyst Layer
3.3.1. The Distribution of Methanol Concentration
Figures 5 and 6 show the distributions of dimensionless methanol concentration with different thickness and specific area of the catalyst layer, operating temperature, methanol concentration, and overpotential. The dimensionless methanol concentration was defined at the value of 1.0 at the inlet boundary () and it decreased whilst the diffusion process occurred through the catalyst layer. The decrease was greater with larger thickness and specific area of the catalyst layer and with higher operating temperature and lower methanol concentration. Comparing Figure 5(a) with Figure 6(b) we see that the dimensionless methanol concentration decreased more rapidly at higher overpotential. This is because the electrochemical reaction was accelerated by higher temperature, methanol concentration, and overpotential in respect to expression (7). Larger specific areas benefited the reaction near the inlet boundary which led to faster consumption of methanol there. As a result, the concentration of methanol became lower in the direction inside the catalyst layer, especially for the thicker catalyst layer. This will result in lower efficiency deep inside the catalyst layer. It is worthwhile to note that the methanol concentration was extremely low, almost zero, at both edges of the modeling domain.
The arrows in Figure 5 and Figure 6 show the gradient of the dimensionless methanol concentration (). We can see that methanol diffuses towards the two edges of the modeling domain (boundaries of and ) when operated at smaller thickness, smaller specific area, lower temperature, and lower overpotential. With increasing of thickness, specific area, operating temperature, and overpotential, methanol diffuses towards the centre due to the decrease of dimensionless methanol concentration in the middle of the modeling domain.
3.3.2. The Distribution of Overpotential
The contours in Figures 7 and 8 show the distributions of dimensionless overpotential with different thickness and specific area of the catalyst layer, operating temperature, methanol concentration, and overpotential. In general, the higher density of the contours the more complicated the distribution of dimensionless overpotential. In each of Figure 7 and Figure 8, three contours closest to the boundary were collected. The values of these three were 0.98, 0.94, and 0.90 from the left to the right. It was clear that the contours of the same value appeared at the different position of the modeling domain due to the influence of thickness and specific area of the catalyst layer and operating conditions such as temperature, methanol concentration, and overpotential. This indicated a greater distribution of the dimensionless overpotential inside the catalyst layer with larger thickness, specific area, and higher operating temperature, methanol concentration, and overpotential.
Similar to the distribution of methanol concentration, the dimensionless overpotential was very low at the edge whilst it was very high close to the middle of the modeling domain. The distributions of dimensionless overpotential in Figures 7 and 8 were not as great as the distributions of dimensionless methanol concentration in Figures 5 and 6. The reason was that the conductivity of the catalyst layer (3.4 Ω−1) was approximately ten orders of magnitude higher than the efficient diffusion coefficient ( m2 s−1) of the catalyst layer. This led to the significant difference of the restriction between the diffusions of methanol and charge.
The lines in Figures 9, 10, and 11 show the dimensionless overpotential at particular coordinates of the catalyst layer (shown in Figure 7(a)). The variation in the distributions of dimensionless overpotential become greater with larger thickness and specific area of the catalyst layer, higher temperature and methanol concentration, and higher overpotential. Furthermore, from Figure 11 we can see that the distribution of dimensionless overpotential was greater at the edge of the surface than that close to the middle.
3.3.3. The Distribution of Current Density
Figures 12 and 13 show the dimensionless current density at particular coordinates on the surface of the domain chosen for modeling. The dimensionless current density decreased to 0 at the boundary of (including A′ points) in both Figures 12 and 13. As shown in Figures 12 and 13, the maximum dimensionless current density appeared at the two corners of the rectangular domain (gray color) close to the boundary of the diffusion layer. The maximum values in Figure 12 are (a) 9.254, (b) 1.481, (c) 0.370, (d) 92.542, and (e) 0.925, and those in Figure 13 are (a) 4.993, (b) 9.254, (c) 16.120, (d) 4.627, and (e) 18.509. Besides the corners, the dimensionless current density appears quite small in the other areas. This was because the gradients of both the dimensionless methanol concentration and overpotential were significantly different in the corners than in the other areas. The more detailed comparisons were shown in Figure 16.
Apparently, an increase in operating temperature, methanol concentration, and overpotential increases the dimensionless current density. However, higher dimensionless current density is observed in the catalyst layer with smaller thickness and specific area. It is because of the higher utilization rate of the interface of the catalyst layer. Moreover, the apparent current density () is observed in the catalyst layer with larger thickness and specific area because the apparent current density is determined not only by the dimensionless current density () but also by thickness () and specific area () as well as exchange current density () according to expression (31).
Figure 12 shows that an increase in thickness and specific area of catalyst layer increases the nonlinearity of the distribution of dimensionless current density. As shown in Figure 13, increasing the temperature and methanol concentration increased the nonlinearity of the distribution of dimensionless current density. However, the effect was not as great as that due to the increased thickness and specific area of the catalyst layer. Comparing Figure 12(a) and Figure 13(b) we see that the nonlinearity of the distribution of dimensionless current density was linked to the increased overpotential.
The effects of temperature and methanol concentration were not significant with a ten-micrometer thick layer because the distributions of methanol concentration and overpotential were not as great with such a thin catalyst layer. As a result, the polarisation curves simulated by macro kinetics were very close to those for intrinsic kinetics.
Briefly, the reason for the nonlinearity of the distribution of dimensionless current density was the speed of the electrochemical reaction which was accelerated by larger thickness and specific area and higher temperatures, methanol concentrations and overpotentials.
Figures 14, 15, and 16 show the detailed distribution of dimensionless current with different thickness and specific area of the catalyst layer, operating temperature, methanol concentration, and overpotential at particular coordinates of the catalyst layer (shown in Figure 12 (a)). The dimensionless current and its nonlinearity increased with smaller thickness and specific area of the catalyst layer, higher operating temperature, overpotential, and larger methanol concentration. From (25), the definition of the dimensionless current , we can conclude that the dimensionless modulus and , which were influenced by the change of the thickness and specific area of the catalyst layer, the operating temperature, the methanol concentration, and the overpotential, affect the dimensionless current. The distribution of dimensionless current was also more serious and nonlinear at the edges of the surface than at the middle if we compare the distribution of dimensionless overpotential in Figure 11.
3.3.4. The Effectiveness of the Porous Catalyst Layer
Figure 17 describes the relationship between effectiveness factor and dimensionless current density with different thickness and specific area of the catalyst layer. In Figure 17, the end points of each curve (the pentagram symbols), were the dimensionless limiting current density obtained under specific physical parameters such as thickness and specific area. The efficiency factor decreased with increasing dimensionless current density until it reached the limit at each end of the curve. The effectiveness factor was higher for a thinner catalyst layer, with lower specific area because of the decreasing influence of mass transport restrictions.
Figure 18 shows the relationship between effectiveness factor and dimensionless overpotential at different operating temperature and bulk methanol concentrations. The effectiveness factor was close to 1.0 at low dimensionless overpotentials, and it fell when the dimensionless overpotential increased. At higher temperature and methanol concentration, the effectiveness factor went through a minimum with increased overpotential. This may be due to a change in the rate determining step from (2) to step (1). This transition of the efficiency factor could have been associated with the additional consumption of adsorbed methanol which was generated by step (1).
4. Conclusions
The 1D model of the DMFC, based on the intrinsic kinetics, demonstrated that the change in coverage of OH () was very fast, with values close to 1.0 in the whole range of overpotential. As a result, the element step (3) was not the rate determining step in the mechanism. The surface coverage, , on the surface of Pt determined the current density at lower overpotential (lower than 0.5 V) whilst the element step (2) could be considered as the rate determining step. The rate determining step changed to element step (1) when the overpotential increased to above 0.6 V.
The 2D model indicated that the performance of the porous Pt-Ru anode depended on the kinetics parameters of the dual-site mechanism of methanol oxidation as well as physical parameters such as thickness, specific area of the catalyst layer. The distributions of methanol concentration, overpotential and current density became more nonlinear with increased operating temperature, methanol concentration, and overpotential, as well as increased thickness and specific area of the catalyst layer. The nonlinearity of all the distributions, including methanol concentration, overpotential, and current density, were higher near the edges of the modeling surface and lower in the middle. At the two corners of the modeling surface, the highest current density appeared due to the maximum difference of the derivative of methanol concentration and overpotential. The distribution of methanol concentration was more obvious than that of overpotential and current density, because the coefficient of mass transport was approximate ten orders of magnitude smaller than the coefficient of charge transport.
The effectiveness factor of the catalyst layer was higher with lower operating temperature, methanol concentration, overpotential, thickness, and specific area. The efficiency factor also went through a minimum value with increasing dimensionless overpotential, at higher temperatures and methanol concentrations. This was probably a result of the change of the rate determining step from element step (2) to (1). This change accelerated the consumption of absorbed methanol which was generated by element step (1).
Nomenclature
| : | Specific area (m−1) |
| : | Dimensionless concentration |
| : | Concentration (mol m−3) |
| : | Diffusion coefficient (m2 s−1) |
| : | Effective diffusion coefficient (m2 s−1) |
| : | Overpotential (V) |
| : | Faraday constant (C mol−1) |
| : | Dimensionless current density |
| : | Current density (mA cm−2) |
| : | Exchange current density (mA cm−2) |
| : | Surface reaction rate constant (mol m−2 s−1) |
| : | Thickness of porous catalyst layer (m) |
| : | Flux of methanol (mol·m−2·s−1) |
| : | Number of electrons involved in reaction |
| : | Apparent reaction rate (mol·m−2·s−1) |
| : | Universal gas constant (J mol−1 K−1) |
| : | Surface reaction rate (mol·m−2·s−1) |
| : | Mass diffusion modulus |
| : | Operating temperature (K) |
| : | Width of porous catalyst layer (m) |
| : | x-direction dimensionless distance |
| : | x-direction distance (m) |
| : | y-direction dimensionless distance |
| : | y-direction distance (m). |
| : | Site density (mol m−2) |
| : | Aspect ratio of the rectangular electrode =l/w |
| : | Electrochemical transfer coefficient |
| : | (V−1) |
| : | Porosity |
| : | Coverage ratio |
| : | Effective conductivity (Ω−1 m−1) |
| : | Charge transfer modulus |
| : | Effectiveness factor |
| : | Potential of open-circuit (V) |
| : | Potential of electrolyte (V) |
| : | Potential of matrix (V) |
| : | Dimensionless relative overpotential . |
| : | Adsorbed |
| : | Effective |
| : | Species j |
| : | Local state |
| : | Methanol |
| : | Total or apparent state. |
| 0: | At the boundary of the diffusion layer . |
Acknowledgments
The authors of this paper gratefully acknowledge the financial support from EPSRC Supergen Fuel Cell Consortium award. Mr. L. Xing would like to say thanks to Wenxiao Huang, who is the master student in Bristol University, for her great help for the figures processing in this paper.