Critique and Improvement of a One-Dimensional Semianalytical Model of a Direct Methanol Fuel Cell

Direct methanol fuel cells (DMFC) have gained much interest as an emerging technology to support the design of power supply systems for portable electronic devices. Key advantages lie in the high energy density and easier storage facilitated by the use of a liquid fuel. However, the technology is affected by some acknowledged drawbacks, such as the relatively slower kinetics of the methanol oxidation process and the undesirable effect of methanol crossover [1], which remain as challenges for the commercialization process. The availability of fuel-cell models that make reliable predictions of operating performance can play an important role in overcoming these challenges.

Methanol crossover, a most intriguing topic in DMFC analysis, is a phenomenon where unreacted methanol fuel is transported from the anode to the cathode, passing through the membrane that separates the electrodes. The driving forces for the methanol crossover are mass transport mechanics involving electro-osmosis, diffusion, and differential pressure between the electrodes. Once at the cathode the crossover fuel is oxidized, causing an adverse thermodynamic effect by consuming oxygen and contributing to a reduction of electrical conversion efficiency by poor fuel utilization and the generation of the equivalent of a leakage current. The effects that multiple physical parameters have on methanol crossover have been explored experimentally in the literature [2–4]; however, practical impediments to making relevant measurements inside the DMFC have proven to be a limiting factor. The complexity of the fluid flow, heat transfer, and chemical reaction processes involved, coupled with physical accessibility constraints of the DMFC architecture, make it difficult to acquire multipoint measurements of temperature, pressure, and species concentrations along a path inside the system for the purpose of characterizing the spatial distribution of these properties. A number of DMFC models have been recently developed with the goal of elucidating details of the transport phenomena inside the cells [5–7] and overcome to a certain extent the physical impediments that inhibit more detailed experimental analysis.

Models for DMFC systems can be formulated with a high degree of complexity. Two-dimensional and three-dimensional models involving two-phase flow are proposed in the literature [8–13]. Coupling such complex flow models to detailed multidimensional models of the diffusional, chemical, and electrochemical processes appears to be a most promising approach for supporting predictive analysis. However, simplified models based on reasonable assumptions that lead to analytical, semianalytical, empirical, or one-dimensional descriptions of DMFC systems, have proven far superior for the purposes of system design and for providing physical insights. The reduced complexity of such models makes them more suitable for use in an engineering design context or in system-level optimization operations, taking advantage of the opportunity for fast numerical execution and simple implementation, though perhaps at a cost in reduction of accuracy and incomplete description of physical-phenomena details.

A particularly attractive simplified model of a DMFC developed by García et al. [14] presents an one-dimensional analysis of a single-membrane assembly configuration that yields a semianalytical description for the system. Simplicity is attained through the adoption of a number of assumptions, including isothermal and steady-state operation, and the restriction that all carbon dioxide is present in dissolved form in a liquid aqueous phase. Using numerical simulation software implemented in current desktop type of computational hardware, the execution time of a simulation is less than 1 min, which enables fast prediction of the cell’s polarization curve and makes the DMFC model easier to implement as part of a more elaborate system-level simulation environment that includes a stack of multiple membrane assemblies and the requisite balance-of-plant components.

The model evaluation carried out in [14] suggests that there is adequate agreement between model-predicted polarization curves and experimental curves at several bulk-concentration levels of methanol fuel. However, upon closer examination of the model it was determined that an implicit equation given in [14] had not been solved to a sufficient level of convergence. When allowed to converge, the implicit equation has the deleterious effect of eroding the accuracy of the match between the predicted and experimental polarization-curve data, particularly at higher current values, leading in some cases to erroneous predictions of limiting current values. In spite of this shortcoming, the model remains potentially successful for the prediction of polarization curve values over a wide range of currents, and hence, the approach holds significant promise as a modeling procedure after the high-current inaccuracies are mitigated.

Another issue of relevance regarding the formulation of the model given in [14] involves the suitability of the concentration boundary condition imposed at the end of the proton exchange membrane (MEM) that neighbors the cathode catalyst layer (CCL). Most of the one-dimensional literature assumes that the methanol concentration at this interface is identically zero. The justification argument typically advanced is that the chemical reaction taking place in the CCL is very fast [14–16]. However, our analysis shows that this boundary condition plays a critical role in the resulting limiting current density, and hence the suitability of the zero boundary condition warrants further analysis.

The objective of this work is to provide a critique of the model proposed by García et al. [14], identifying advantages and potential shortcomings, and to propose a number of improvements that lead to better predictions of experimental data. The contributions of this work are in the form of a critique that quantifies the extent of discrepancy with experimental data experienced by the model as proposed in [14], most notably at high current-density regimes, and include the proposition of a number of concrete modeling modifications that can improve the quality of the model’s predictions.

Another improvement involves a new method to predict the maximum current density based on further analysis of the underlying kinetic reaction-rate expression is included. The method can reduce modeling errors when during a simulation study the operating current density approaches the limiting values.

A second improvement involves the adoption of systematic least-squares regression technique that permits the extraction of parameters values using limited experimental data. This advancement is of relevance because DMFC models typically include a large number of physical parameters. Due to an incomplete literature base, or lack of experimental resources, it is not always feasible to gain access to all needed parameters. A common approach is to estimate parametric values based on known systems that bear a physical similarity to the DMFC environment being analyzed, or alternatively, resort to the extreme of assuming parametric values without further justification when no information is available. Our new approach may prove of particular utility in conditions where parametric values are scarce.

An advancement of relevance concerns the analysis of the impact of the specification of an appropriate concentration boundary condition at the interface between the MEM and the CCL for the purpose of improving the predictive accuracy of the simplified model.

For self-consistency, this section catalogs the modeling equations proposed by García et al. [14], followed by the closed-form solution of the system equations derived by those authors. The treatment is deliberately succinct, and the reader is referred to the original publication source for details. A description of the symbols utilized is given in the Nomenclature.

A schematic diagram of the DMFC considered consists of a sequence of four layers defining the membrane electrode assembly (MEA), configured as indicated in Fig. 1. The intermediate layers indicated in the figure are the anode catalyst layer (ACL) and the MEM. The thickness of each layer is indicated by the subscripted variable $δ$⁠. A horizontal coordinate $z$ spans the MEA, ranging from $z=0$ corresponding to the left-most boundary of the anode backing layer (ABL), to $z=zIV$ corresponding to the right-most boundary of the CCL. The dashed lines indicate the methanol concentration profile along the $z$ direction that is predicted in [14] using the modeling equations. The concentration discontinuities at coordinates $zI$ and $zII$ are due to classical partition-coefficient distributions associated with mass-transfer equilibrium across a boundary.

where $KI$ and $KII$ are the partition coefficients that, respectively, characterize the local equilibrium across the ABL/ACL and ACL/MEM interfaces.

Note that the cell current expression given in Eq. (20) can be calculated using the pointwise methanol-concentration profiles given in the analytical expressions (1) to (7). However, because these analytical concentration expressions are themselves functions of the cell current, we refer to Eq. (20) as a semianalytical equation since it represents an equation that is implicit in the cell current variable $Icell$⁠.

is the leakage current density due to methanol that crosses over from the anode side and is oxidized in the CCL layer

The steady-state DMFC model proposed by García et al. [14] has a number of significant merits. In particular, it is a highly compact model due to its ability to produce analytical methanol concentration expressions and semianalytical cell-current expressions. As a consequence it possesses the advantage of computational simplicity, allowing fast calculation of polarization curves, including a prediction of the effect that the bulk-fuel concentration in the anode channel has on the polarization curve. Our experience suggests that the model is particularly adept at predicting polarization curves that match the Ohmic regime of the experimental data reported in [14].

In spite of its numerous advantages, the model [14] has relevant shortcomings that need to be addressed for the purpose of strengthening its predictive capabilities and better match experimental data over a wider current-density range. In the ensuing subsections we identify specific modeling issues that need to be revisited, and in the next section we propose specific modeling modifications and practices that can potentially overcome inadequacies of the model.

The original publication of García et al. [14] contains errors in the reported values of two model parameters of the ABL, namely, the effective diffusion coefficient for methanol $DB$ and the partition coefficient $KI$⁠. The corrected values for these two parameters are given in Table 1 [19]. Note that the partition coefficient in question, which is taken from the work of Baxter et al. [20], is corrected to adopt a value equal to the inverse of the value documented in [14]. Note also that the corrected value is consistent with the physical expectation that the methanol concentration in the ABL side of the ABL/ACL interface should be higher than the value on the ACL side of the interface.

Once these two parameters are corrected, the numerical execution of the model shows that the quality of match to the experimental data presented in [14] is significantly eroded. For the purpose of recovering the model’s ability to match experimental polarization curve data, a corrected value for $λ$ the kinetic rate parameter appearing in the ACL Eq. (11) was adopted as reported in Table 1 [19]. The modeling results with the corrected parameters is shown in Fig. 2, where the polarization curve represented with a dashed line is obtained using the arbitrarily selected total number of iteration $ni=5$ suggested in [19] to solve the implicit Eq. (20), and where the circle markers represent experimental data. Hence, the polarization curve described by the dashed line, corresponding to an errata free model and an improved kinetic rate parameter, shows that under these conditions the corrected model [14] with fixed iteration parameter $ni=5$ appears to yield a reasonably good match to experimental data. Unfortunately, as argued in the next section, this result is not acceptable because the polarization curves given by the dashed lines are obtained when the procedure for solving Eq. (20) is arbitrarily interrupted after five iterations, when the process has not yet converged to a solution. Consequently, the polarization-curve result is not an accurate representation of the prediction that would be made by Eq. (20) when it is appropriately solved to a satisfactory level of convergence.

The approach for solving Eq. (20) requires that the total number of iterations $ni$ be selected so that the process converges to a solution. If fewer iterations are made, prematurely interrupting the process before full convergence, the numerical result does not produce a prediction of $ηA$ that is consistent with Eq. (20). In fact, the resulting $ηA$ fails to satisfy the implicit Eq. (20).

Paradoxically, in spite of this procedural shortcoming the dashed-line curves in Fig. 2 show a good fit to the experimental data over the entire range of currents. This result is due in part to the fact that the good fit was obtained by adjusting the values of other modeling parameters to compensate for the model-data mismatch introduced by the early interruption of the iterative solution procedure. The polarization results obtained in this fashion are therefore not acceptable because the modeling equations are not properly solved due to lack of convergence of the required iterative procedure.

A natural remedial measure for consideration is to increase the number of iterations until satisfactory convergence is obtained. Unfortunately, the models’ polarization curve predictions under full convergence conditions are still affected by error in the high-current regime, as shown by the solid lines in Fig. 2, leading to a significant level of mismatch with experimental data. This suggests that the proposed model suffers from a structural inadequacy to describe the polarization curve at high current densities.

The iterative procedure used for solving Eq. (20) requires that the user specify an initial guess for the current density. The iterative procedure fails to converge when the initial guess is improperly chosen to be higher than the maximum current density that can be realized by Eq. (20). Hence, the execution of the model using computer code requires further guidance for the user on how to produce feasible initial guesses. The scope of the analysis in [14] does not include a technique for estimating the maximum current density under given operating conditions. We propose an initial-guess bounding technique in Sec. 4.2.

It is a common practice for fuel-cell models to adopt anodic and cathodic transfer coefficients reported in the literature. However, due to its relatively complex nature, the reaction mechanism in an operating fuel cell typically involves a number of parameters that are difficult to measure experimentally. Equation (20) reveals that two parameters that have a significant impact on the model’s prediction, namely the anodic transfer coefficient $αA$ and the activation overpotential $ηA$ (which appear as the product $αAηA$⁠), do not have widely reported values. Because the parameters appear as a product, the practice of arbitrarily setting the value of one of the two parameters necessarily leads to an erroneous estimate of the other parameter given that the product must remain constant. The problem also involves with the cathodic analog of these two parameters (which also appear as the product $αCηC$⁠). One possible method for determining these parameters in a more robust fashion is to extract their values from the experimental data using a standard least-squares method (LSM). A detailed and systematic LSM procedure is presented in Sec. 4.4.

Several parameters in the DMFC model under consideration are not available in the literature. Parameter values may be available under idealized experimental conditions that do not adequately represent the multivariable conditions of an operational DMFC; hence, some of the experimentally measured parameters, such as diffusion coefficient, conductivity, water content of the membrane, etc., may be of low practical utility for modeling purposes. The challenge arises from the fact that the complicated nature of an experimentally operational DMFC system makes it difficult to measure some parameters in situ.

A possible approach to overcome this data deficiency consists of identifying suitable parameter values by comparing the model-generated polarization curve with experimental data. This venue, however, may prove onerous if carried out by a trial-and-error mechanism or by an ad hoc method. The usefulness of the model proposed in [14] would be significantly enhanced if supplemented by systematic methods for parameter estimation from experimental data. A method illustrating a systematic philosophy for guiding the extraction of physical parameters from experimental data using least-squares regression is given in Sec. 4.4.

The boundary condition at the MEM/CCL interface has a significant effect on the modeling results. Most of the literature commonly assumes that the methanol concentration at the MEM/CCL interface is equal to zero, based on the argument that the reaction rate is very fast relative to the diffusion rate controlling methanol crossover into the CCL. This assumption also enables the attainment of an analytical solution from Eqs. (1) to (7). However, the validity of the assumption depends on the operating conditions of the DMFC. As demonstrated in Fig. 5 of García et al. [14], where the operating current density used to generate the polarization curve is 0.015 A cm⁻², which is approximately one tenth of the limiting current density, and yet the model is still based on the zero-concentration boundary condition even though under such low current conditions one may expect a non-negligible methanol concentration at the boundary in question.

One may argue that the methanol concentration can be non-negligible at the MEM/CCL interface, especially at current densities significantly lower than the limiting current density. The assumption of zero concentration is equivalent to assuming an infinite Damköhler number [21]. This assumption may be reasonable in many simple surface-reaction-driven catalysis problems [22]. However, the present problem is rather complex, as the oxygen and methanol reactants diffuse from opposite directions towards the MEM/CCL interface, and the reaction takes place not only at the interface but also within the bulk of the CCL. Furthermore, consumption of methanol within the CCL can be inferred under certain conditions where residual crossover methanol has been measured in the cathode-side exit stream [23]. Hence, there is no experimental evidence that the reaction must occur exclusively or primarily at the interface, which would be necessary for establishing an infinite Damköhler number.

In Sec. 4.5 we demonstrate that using nonzero methanol concentration values at the MEM/CCL interface in the DMFC model under consideration sensitively affects the maximum current density predicted by the model as well as the methanol concentration distribution profile throughout the MEA. Hence the adoption of an accurate boundary condition at this interface is critical.

Five approaches are pursued in this section with the goal of improving the ability of the model of García et al. [14] to produce predictions that better match available experimental data. The first approach consists of advancing the iteration algorithm to guarantee the modeling results converge to a solution. Second, the maximum current density under a given condition is predicted numerically to avoid failures in iteration convergence. Third, modifying the modeling equations with the objective of improving the polarization-curve quality at high currents, where significant deviations from experimental measurements are observed. Fourth, a systematic method for identifying relevant modeling parameters from available experimental data is formulated for the goal of improving parametric consistency. And finally, the effect of the methanol concentration value used to define a boundary condition at the MEM/CCL interface is investigated in terms of its effects on the methanol concentration profile across the membrane assembly, as well as in terms of its effect on the polarization curve characteristics.

where the $abs$ operation designates the absolute value. If the conversion criterion is not satisfied in step (v), then the current estimate for $ηA$ is updated and steps (iii) to (v) are repeated. The convergence criterion (23) calls for termination when the relative error in current density obtained with the current $ηA$ update is less than 1%. The update at step (v) can be done using standard numerical root-finding routines for single-variable problems [24].

Step (i) of the $ηA$ iteration scheme described in Sec. 4.1 involves an initial guess of current density that satisfies the constraint $Icell≤Imax$⁠, where $Imax$ is the limiting current of the DMFC. The convergence scheme fails when this constraint is violated because the difference calculated in step (iv) may never evolve to reach an absolute value sufficiently small to meet the specified convergence criterion. This lack-of-convergence problem can be avoided if the user can first obtain a reliable estimate of $Imax$⁠, and hence ensure the satisfaction of the constraint in step (i) of the algorithm.

Note that the methanol concentration in the anode catalyst layer (appearing in the integrand) is itself a function of the current density (in this case $Icell=Imax$⁠). Hence an iterative algorithm analogous to that described in Sec. 4.1 is used to solve for $Imax$ from Eq. (25). Note that every iteration involves the update of the variables $C1$ and $C2$ in Eqs. (4) and (5).

Further insight into the behavior of the model near the maximum current-density point is gained by recognizing that the voltage should be expected to drop dramatically in the high current-density regime. This is a consequence of the fact that the anodic overpotential $ηA$ becomes infinitely large at such large current densities, and hence from the term $-ηA$ appearing on the right-hand side of the voltage Eq. (19) it then follows that the cell voltage must experience a theoretical singularity as the maximum current-density point is approached. Hence, in agreement with other theoretical models, this framework also predicts a singular voltage behavior at the maximum current-density point, even though the assumptions adopted for the model’s formulation do not include the logarithmic singularity induced by the concentration overpotential model that is commonly invoked in more classical modeling approaches [25]. The singularity at the maximum current density is depicted in Fig. 2 by the near vertical slopes of the polarization curves in the high-current limit for each concentration considered. Hence, the modeling predictions obtained by carrying to full convergence the iterative procedure described in Sec. 3.2 provide polarization-curve profiles that are consistent with the theoretical singularity inherent in the model’s formulation.

Note that we deliberately introduce the terminology “maximum current density” instead of “limiting current density” to highlight the fact that the former is found by solving a modeling equation and is hence only an estimate of the latter. Equation (25) also shows that the maximum current density is related to the specific surface area of the anode catalyst, the exchange current density, and the methanol concentration profile in the catalyst layer, among others. Those are indeed the design parameters available to increase the maximum current density in a DMFC.

As discussed in Sec. 3.2, there is a significant level of polarization-curve matching error at high current densities realized by the model of García et al. [14] after the model is executed so that the iterative solution of Eq. (20) meets the convergence criterion (23) and the corrected values of modeling parameters listed in Table 1 are adopted. The model/data mismatch is illustrated in Fig. 2, where the solid lines are the modeling results after convergence is reached. It can be observed that, for all concentrations of anodic methanol solution reported in the figure, there is a marked discrepancy between the model-predicted polarization curve and the experimental points in the high current-density region. This consistent mismatch error in the transport-limited region suggests that the model formulation needs to be fine tuned to improve its ability to describe the physics of the transport-limited regime.

We propose revisiting the ion conductivity assumption adopted in García et al. [14] as one possible venue to reduce the extent of the documented mismatch. In [14] the ion conductivity is assumed to adopt the constant value $κ=3.6 s m-1$⁠. However, the ion conductivity may be expected to vary with increasing current density, depending on the operating conditions. In operational fuel cells the chemical reactions are most likely nonuniformly distributed across the catalyst layer, resulting in local hot spots that would be more pronounced as the cell operates at higher current densities. In turn, a local hot spot may further promote dry-out or material cracking inside the cell assembly. As is well known in the case of proton-exchange-membrane fuel cells, the dry-out phenomena results in an increase of internal resistance [26,27], leading to an increase in the downward slope of the polarization curve in the regime where such hot spots are more likely, namely in the transport-limited regime characteristic of higher current densities.

where $Imax$ is the maximum current density of the cell, and where the non-negative parameters $β0$ and $β1$ are the degrees of freedom of the model. In particular, the parameter $β1$ is intended to characterize the phenomenon commonly reported in the literature where there is a rate of reduction in conductivity as the current increases [26,27]. We propose that the estimation of these two parameters by carried out by least-squares [24] linear regression on experimental polarization data (details are omitted for brevity of exposition).

The anodic and cathodic overpotentials $ηA$ and $ηC$ and the transfer coefficients $αA$ and $αC$ in the current-modeling Eqs. (20) to (21) must be specified by the user. In principle, all four parameters can be estimated via laboratory experiments carried out under closely controlled conditions. In practical applications, however, literature sources are often unavailable and the parameters must be estimated by numerical fit to experimental data obtained from a working fuel cell.

The technique adopted in García et al. [14] for addressing this parameter estimation problem involves assuming a priori the values of $αA$ and $αC$⁠, and then proceeding to identify values of $ηA$ and $ηC$ that when incorporated into the DMFC model produce an acceptable match to experimental polarization-curve data. While this procedure may succeed in yielding a good model-based polarization-curve prediction, it may nevertheless fail to identify reasonable values for the individual overpotential parameters. The problem stems from the fact that in the modeling Eq. (20), the anodic overpotential and mass-transfer coefficient always appear in the form of the product $αAηA$⁠. As a consequence, a given numerical value $γA=αAηA$ for the product that leads to good current-voltage predictions can unfortunately be realized by an infinite number of pairs of $(αA,ηA)$ values. For example, if a pair of values $(αA,ηA)$ yields a unique product-value $γA$⁠, all pairs of parametric estimates $(pαA,ηA/p)$⁠, where $p$ is an arbitrary nonzero number, yield the same product $γA$ because the factor $p$ cancels out when both parameters are multiplied to obtain $γA$⁠. The same argument applies to the product $αCηC$ appearing in Eq. (21). In summary, if the a priori estimates of $αA$ and $αC$ are affected by error, the resulting estimates for the overpotentials $ηA$ and $ηC$ may be highly uncertain due to the error propagation induced by the parameter-product constraint discussed.

We propose a more robust and systematic method for identifying both transfer coefficients. To that end, consider the parameter-product variables $γA=αAηA$ and $γC=αCηC$⁠. Next, consider the experimentally determined quantities $Yi=-(Vcell-UO2+UMeOH+δMIcell/κ)i$⁠, where $i$ runs from $1$ to the total number of experimental data points $n$⁠. Then, $αA$ and $αC$ can be obtained using standard least-squares regression techniques [24] after recognizing that the relationship $Yi=αA-1γA,i+αC-1γC,i$ is linear with respect to the unknown inverse parameters $αA-1$ and $αC-1$⁠.

Hence, the methanol concentration gradient in Eq. (14) is replaced by Eq. (31). Note that the assumed boundary value $CIIIM$ adopted in Eq. (27) effectively influences the concentrations $CIA$ and $CIIA$⁠, as indicated in Eqs. (29) and (30). Note also that when one sets $CIIIM=0$⁠, Eqs. (28), (29), and (30), respectively, reduce to the original expressions (3), (6) and (7) derived by García et al. [14].

Figure 3 shows the experimental data along with modeled polarization curves obtained after the adoption of all the improvements proposed in this work. The circle markers denote experimental data reported in [14] at four bulk methanol compositions. The solid lines represent polarization curves produced by an improved model that adopts the modifications described in Sections 4.1–4.4. In particular, the model includes full convergence of the implicit-equation algorithm, and the implementation of the proposed ion conductivity model given in Eq. (26) with the parameter values $β0=7.69$ and $β1=6.55$ for all concentrations. The concentration-dependent transfer coefficients $αA$ and $αB$⁠, obtained as discussed in Sec. 4.4, are reported in Table 3. The model incorporates the errata corrections reported in Table 1, and the standard MEM/CCL boundary condition $CIIIM=0$ is adopted.

Figure 3 demonstrates that the incorporation of the proposed modeling refinements leads to a significant improvement in the model’s ability to better match experimental polarization-curve data over the range of reported concentrations. The extent of improvement can be readily assessed in a qualitative fashion when comparing the plots of this figure with those given in Fig. 2 (where the errata-free model is implemented without the proposed refinements). The modeling curves in Fig. 3 are obviously much closer to the experimental data over a very wide range of current densities.

A quantitative analysis shows that the errata-free model proposed in [14] and used to generate the polarization curves reported in Fig. 2 leads to a worst-case absolute error in cell voltage in excess of $100%$ with respect to the experimental data for the case of a bulk methanol concentration of 0.5M. In contrast, the worst-case prediction error associated with the improved-model predictions of Fig. 3 for the same concentration does not exceed $1%$⁠. Another quantitative assessment of the extent of modeling improvement can be made from Fig. 3 by focusing on the predictions of maximum current densities for each case of bulk methanol concentration considered. Note that the maximum current density predicted by the model using Eq. (25) can also be identified from Fig. 3 as the point of intersection of the horizontal axis and the solid-line polarization-curve for each bulk concentration of methanol. The results are reported in Table 2, where the first column lists the bulk methanol concentration (⁠$Cb$⁠) and the second column the corresponding maximum current (⁠$Imax$⁠) extracted from Fig. 3. It is reasonable to accept that the best experimental estimate of the maximum current density afforded at each bulk methanol concentration value is given by the experimental data point corresponding to the largest current density reported (i.e., the right-most circle marker in each polarization curve of Fig. 3). Using those estimates as a reference basis, the results in Table 2 can be used to readily establish that the improved-model predictions for $Imax$ match the experimental estimates with a discrepancy no greater than a $1%$ absolute error at all bulk methanol concentrations. This analysis lends quantitative support to the claim of that the improved model leads to more accurate predictions.

The effect of the proposed current-dependent ion conductivity model given in Eq. (26) is captured in Fig. 3 through the adoption of the least-squares estimated parameters $β0=7.69$ and $β1=6.55$ obtained using only data at a bulk methanol concentration of 0.5M, and assumed to remain fixed at those values for all other concentrations. As discussed in the previous subsection, Fig. 2 indicates that the errata-free model proposed in García et al. [14] realizes a worst-case absolute error in excess of 100% with respect to the experimental data for a bulk methanol concentration of 0.5M. Figure 3 shows that at the same methanol concentration the incorporation of the linear conductivity model reduces the discrepancy to a level that is lower than $1%$⁠. Hence, for this methanol concentration, the improved model leads to a significant improvement in polarization-curve fit.

Further comparison of the plots presented in Figs. 2 and 3 reveals that the improved model also leads to partially enhanced predictions at a methanol concentration of 0.2M, but to only marginal enhancements or no visible enhancements at the two lowest concentrations of 0.1M and 0.05M. The large accuracy improvement in the fit to the 0.5M concentration is an expected consequence of the choice made to extract ion-conductivity parameters $β0$ and $β1$ using least-squares regression using only the 0.5M data. This implies that an improved overall fit would be obtained if new pairs of ion-conductivity parameters $(β0,β1)$ are regressed at each composition, and then used to build a model that reflects the dependency of these parameters on the bulk composition.

The formulation of the ion conductivity model given in Eq. (26) is designed to induce an increase of the internal resistance to ion transport as the operating current density increases. The corresponding increase of the overpotential contributes to a reduction of discrepancy with the experimental data. This internal resistance increase is consistent with phenomena observed when a DMFC operates in the high current-density region, especially under an excess of hot dry air and at a high cell temperature. As mentioned in Sec. 4.3, under these conditions a nonuniform electrochemical reaction could create local hot spots inside the DMFC, and the available internal heat generated at the hot spots could evaporate adsorbed liquid water, causing local dry-out regions [28–30]. These conditions in turn eventually cause local material fractures or cracks, leading to an effective resistance increase. Hence, the effect of our proposed conductivity model in the form of an increased internal resistance is effectively established as consistent with the proposition that under such conditions the overpotential increases as a function of the operating current density.

The experimental polarization data reported in [14] is used to deploy the parameter estimation algorithm described in Sec. 4.4 for the purpose of extracting values for the anodic and cathodic transfer coefficients $αA$ and $αC$⁠. Only data lying in the Ohmic region are considered, given that this regime is of the most practical interest. The results are reported in Table 3, where the first column shows the methanol bulk concentrations considered, and the remaining columns show the values obtained for the two transfer coefficients, numerically rounded to three decimals.

Figure 4 shows the effects that the new MEM/CCL methanol concentration boundary condition proposed in Eq. (27) has on the methanol concentration distribution along the MEA when the bulk methanol concentration is $Cb=0.5M$⁠. The figure contains concentration profiles corresponding to three boundary-condition values for the methanol concentration $CIIIM$⁠, namely 0.0M (the assumption most commonly made in the literature), 0.1M, and 0.2M. The horizontal coordinate $z$ is as defined in Fig. 1, with origin $z=0$ at the beginning of the ABL, and with coordinates $zI=0.015 cm$⁠, $zII=0.0173 cm$⁠, and $zIII=0.0353 cm$⁠. The concentration profiles are calculated under a current density of 0.015 A cm⁻² along the MEA.

A first conclusion derived from these plots is that the concentration gradients of methanol at the ABL, ACL, and MEM regions change to a noticeable extent when varying the methanol concentration at the MEM/CCL interface from 0.0M to 0.2M. Note that in the ABL and MEM regions the gradients become smaller as the boundary concentration increases. On the other hand, the gradient increases in the ACL region. Given the linear form of the profiles in the ABL and MEM regions, the decrease in slope corresponds to a net decrease of concentration change across each one of these regions. As a consequence of the lower slope in the ABL region, the concentration of methanol at location $zI$ increases, and in turn this defines a higher boundary concentration at the beginning of the ACL layer, where the reaction takes place. Since the reaction rate of methanol increases with increases of concentration, the resulting concentration slope observed in the ACL region necessarily increases.

This analysis supports the conclusion that the adoption of different boundary concentrations in (27) can significantly affect the methanol concentration profile across all the relevant MEA sublayers. Furthermore, a departure from the standard zero-concentration boundary condition leads to the realization of higher methanol concentrations in each layer. It follows that the appropriate specification of this boundary condition is of high physical significance, given that it may play a dominant role in the model’s ability to predict cell performance in a number of scenarios of interest, such as the quantification of methanol crossover effects.

Figure 5 shows the effects of the methanol-concentration boundary condition defined in (27) on the polarization curve predicted by the model at a bulk methanol concentration of 0.5M. The MEA/CCL boundary concentration values are the same as those reported in Fig. 4. The circle markers on the figure correspond to the experimental data reported in [14] at $Cb=0.5M$⁠. To avoid the potentially confounding effects of introducing all the modeling improvements proposed in this work, the modeled polarization curves included in the figure do not incorporate any of the improvements advocated in Secs. 4.1 to 4.4, but include the boundary-condition modification given in Eq. (27). Note that as a result of these modeling choices there is an anticipated model-data mismatch that is particularly salient at high current densities. An errata-free model is adopted in all cases. The maximum current density is read from these plots as the point on the horizontal axis that corresponds to zero voltage.

From the plots shown in Fig. 5 it is concluded that an increase in the boundary concentration in question causes an effective increase in the maximum current density predicted by the model. In fact, the maximum current density rises from 0.150 to 0.179 A cm⁻² in the boundary concentration range considered. Furthermore, as the boundary concentration increases, the polarization curve is translated to the right in the sense that a given voltage is obtained at higher current densities; thus leading to the prediction of improved cell performance. These observed trends are consistent with the fact that the polarization curve and the maximum current density are governed by the methanol concentration profile. Indeed, as shown in Eqs. (20) and (25), the methanol concentration profile in the ACL region affects the polarization curve and the maximum current density. Given that as argued in the context of Fig. 4, an increase of the methanol concentration at the MEM/CCL interface causes an increase in the methanol concentration in the ACL region, it follows that the model consistently predicts increases in cell performance and in maximum current density.

Errata in the model published in [14] obtained per communication with the authors of that reference is documented and incorporated as needed to correct the model. It is shown that, in spite of its significant potential as a useful tool, the errata-free model [14] still experiences discrepancies with experimental results, pointing out the need for modeling refinements. A particular improvement proposed and successfully tested is in the form of a method to predict the maximum current density for given operating conditions, providing a venue to eliminate a convergence failure of the iterative algorithm used to solve the implicit modeling equation that predicts the current density. A side benefit of this result is that the dominant parameters that affect the maximum current density are made evident; hence providing DMFC designers with a useful tool for engineering the system’s performance.

Another successfully verified refinement is a systematic least-squares method for extracting key modeling parameters from experimental data, avoiding the need for making arbitrary assumptions on the values of such parameters when there is no reliable literature reference available. The models’ ability to fit experimental polarization-curve data increases dramatically when transfer coefficient values extracted via the least-squares method are adopted. The least-squares approach for extracting parameters is a more rigorous and systematic approach than the alternative technique based on adjusting parameters via a trial-and-error procedure. The linear regression approach could be used to characterize the parameters of a new DMFC design or an already manufactured DMFC, enabling the scientific and engineering team to better identify parametric values and thus build a more reliable mode, better capable of describing the system’s behavior in the areas of interest.

An ion conductivity relationship of a linear form, conceived to reduce discrepancies with experimental data at high current densities, succeeds in attaining the objective of yielding better modeling predictions. Our proposed ion conductivity model is based on a linear-effect conjecture whose plausibility is confirmed through improved polarization curve predictions at high current density. This is indeed an encouraging observation; however, the conjecture needs to be confirmed, rejected, or refined through careful follow-up experimental investigations. We plan to produce experimental studies in the future to address this conjecture, and we encourage the research and development community to undertake analytical and experimental work leading to the elucidation of the mechanistic changes that affect ion conductivity at higher current densities.

Finally, it is concluded that the methanol-concentration value adopted at the MEM/CCL interface has a large impact on the overall polarization curve, the maximum current density, and on the methanol concentration distribution throughout the MEA. This boundary interface condition is arguably a key assumption in the production of more accurate model predictions. It follows that the standard assumption of zero concentration at this boundary may not be an appropriate choice for all operating conditions. We encourage the research and development community to also embark in an effort to produce experimental evidence useful for quantifying the appropriate value of methanol concentration at this boundary as a function of operating conditions. This topic is within the scope of our near-future experimental plans.

The authors gratefully acknowledge financial support from the Communications-Electronics Research, Development, and Engineering Center (CERDEC)/U.S. Army REDCOM.

Nomenclature

a =

specific surface area of the anode, m⁻¹

$i$ =

index for experimental data point, dimensionless

$k$ =

constant in the rate expression, dimensionless

$n$ =

the total number of experimental data points

$z$ =

coordinate direction normal to the anode

$Cb$ =

bulk concentration of methanol in the flow field, mol m⁻³

$CI$ =

concentration of methanol at the ABL/ACL interface, mol m⁻³

$CII$ =

concentration of methanol at the ACL/MEM interface, mol m⁻³

$CIII$ =

concentration of methanol at the MEM/CCL interface, mol m⁻³

$CO2$ =

concentration of oxygen, mol m⁻³

$DA$ =

effective diffusion coefficient of methanol in the ACL, m² s⁻¹

$DB$ =

effective diffusion coefficient of methanol in the ABL, m² s⁻¹

$DM$ =

effective diffusion coefficient of methanol in the MEM, m² s⁻¹

$F$ =

Faraday’s constant, $96487 C equiv-1$

$Icell$ =

cell current density, A m⁻²

$Imax$ =

maximum current density, A m⁻²

$Ileak$ =

leak current density due to methanol crossover, A m⁻²

$I0,refMeOH$ =

exchange current density of methanol, A m⁻²

$I0,refO2$ =

exchange current density of oxygen, A m⁻²

$K$ =

partition coefficient

LSM =

least-squares method

$R$ =

universal gas constant, $8.314 J(mol K)-1$

$T$ =

temperature, K

$UMeOH$ =

thermodynamic equilibrium potential of methanol oxidation, V

$UO2$ =

thermodynamic equilibrium potential of oxygen reduction, V

$Vcell$ =

cell voltage, V

Greek

$αA$ =

anodic transfer coefficient

$αC$ =

cathodic transfer coefficient

$γA$ =

product of $αA$ and $ηA$

$γC$ =

product of $αC$ and $ηC$

$δA$ =

ACL thickness, m

$δB$ =

ABL thickness, m

$δC$ =

CCL thickness, m

$δM$ =

membrane thickness, m

$ηA$ =

anode overpotential, V

$ηC$ =

cathode overpotential, V

$κ$ =

ion conductivity of the membrane, S m⁻¹

$λ$ =

constant of the rate expression, mol m⁻³

$ξMeOH$ =

electro-osmotic drag coefficient of methanol

Subscripts

$b$ =

bulk

$A$ =

ACL

$B$ =

ABL

$C$ =

CCL

I =

ABL/ACL interface

II =

ACL/MEM interface

III =

MEM/CCL interface

IV =

end of CCL

$M$ =

membrane

$MeOH$ =

methanol

Superscripts

$M$ =

membrane

$MeOH$ =

methanol

$O2$ =

oxygen