Numerical Simulation And Optimization Of A Direct Methanol Fuel Cell

Computers and Chemical Engineering xxx (2005) xxx–xxx

Numerical simulation and optimization of a direct methanol fuel cell Cong Xu, Peter M. Follmann, Lorenz T. Biegler ∗ , Myung S. Jhon Department of Chemical Engineering, Carnegie Mellon University, Pittsburgh, PA 15213, USA Received 17 July 2004; received in revised form 15 March 2005; accepted 15 March 2005

Abstract We examine a direct methanol fuel cell (DMFC) model that captures the essence of electrode kinetics and methanol crossover through the polymer electrolyte membrane. Model parameters and key factors for the DMFC model including methanol crossover are identified. Moreover, we establish a relationship between the methanol feed concentration and the power density at a given current density. To gain insight into the effect of the methanol feed concentration, we also performed sensitivity analysis between the cell voltage and the methanol feed concentration. From this sensitivity curve, we are able to identify the optimal feed concentration, which provides the highest power density output for a set value of the current density. Finally, the optimal response of the cell voltage to the change of the feed concentration is studied via dynamic optimization to the differential algebraic equation system. These dynamic optimization results provide a constant feeding strategy that achieves the highest power density at given operating conditions specified by a set current density. © 2005 Elsevier Ltd. All rights reserved. Keywords: Direct methanol fuel cell; Differential algebraic equation; Sensitivity; Methanol crossover; Power density; Dynamic optimization

1. Introduction Fuel cells have attracted considerable attention recently as a possible replacement for power generation systems. For optimal performance and design, accurate system modeling is needed. A fundamental understanding in electrochemistry, materials, manufacturing, and heat and mass transfer are, therefore, critical for developing an optimum model, which can satisfy the goals of cost competitiveness, high performance, and reliability for the fuel cells. There are two types of polymer electrolyte membrane (PEM) fuel cells, i.e., hydrogen fuel cells and direct methanol fuel cells (DMFC), both of which utilize PEM to transfer protons. The DMFC feeds a solution of methanol and water (typically 0.5–2 M; Mench, Wang, & Thynell, 2001) to the anode, where it is internally reformed by the catalyst. Compared to hydrogen fuel cells, DMFC is advantageous for its ease of fuel delivery and storage, lack of humidification requirement, and ∗

Corresponding author. Tel.: +1 412 268 2232; fax: +1 412 268 7139. E-mail address: [email protected] (L.T. Biegler).

0098-1354/$ – see front matter © 2005 Elsevier Ltd. All rights reserved. doi:10.1016/j.compchemeng.2005.03.007

reduced design complexity. Due to the absence of ancillary equipment, the DMFC is ideally suited for portable electronics. Despite these compelling advantages, the DMFC’s low electro-activity of methanol oxidation at the anode and large amount of undesired methanol transported through the PEM from the anode to cathode still remain as important technological concerns. Here, the methanol crossover is defined as the loss of methanol from the anode compartment to the cathode compartment through the membrane and reacted in the cathode compartment. Methanol crossover in DMFC has been extensively studied both experimentally and theoretically (Baxter, Battaglia, & White, 1999; Dohle, Divisek, & Jung, 2000; Narayanan et al., 1995; Ren, Springer, & Zawodzinski, 2000; Ren, Zawodzinski, Uribe, Dai, & Gottesfeld, 1995). Narayanan et al. (1995) and Ren et al. (1995) showed that the methanol crossover rate is inversely proportional to the membrane thickness, which suggests that the diffusion process dominates the methanol transport through the membrane. In addition, Ren et al. (2000) compared the diffusion process with methanol electro-osmotic drag processes and demonstrated

2

C. Xu et al. / Computers and Chemical Engineering xxx (2005) xxx–xxx

Nomenclature AS

(m2 )

cross-sectional electrode area concentration (mol/m3 ) concentration of component i in anode catalyst layer (mol/m3 ) ciF feed concentration of component i (mol/m3 ) F∗ cCH3 OH optimal feed concentration of CH3 OH (mol/m3 ) M cH proton concentration in membrane, 1200 mol/m3 C double layer capacity (F/m2 ) M d PEM thickness (m) M DCH diffusion coefficient of methanol in membrane 3 OH (m2 /s) M DH diffusion coefficient of protons in membrane, 5.4 × 10−9 m2 /s F Faraday constant, 96485 C/mol icell current density (A/m2 ) rate constant of reaction j, (mol/(m2 s)) kj LS k mass transfer coefficient (m/s) Kj equilibrium constant of reaction (j) 2 nM CH3 OH methanol crossover (mol/(m s)) M ni molar flux density of component i in membrane (mol/(m2 s)) p total pressure (Pa) q model parameters r reaction rate (mol/(m2 s)) R universal gas constant, 8.314 J/(mol K) t time (s) T cell temperature (K) U electrode potential (V) Ucell total cell voltage (V) anode feed flow rate, 1.36 ml/min VF V anode compartment volume (m3 ) c ciCL

Greek letters α charge transfer coefficient ηa anode overpotential (V) ηc cathode overpotential (V) κM conductivity of membrane, 17 −1 m−1 θ surface coverage τ mean residence time in anode compartment (s) Superscript CL catalyst layer F feed LS from liquid bulk to surface M membrane * optimal 0 at standard conditions i mass component j reaction steps in the electrode

Subscript a anode c cathode CH3 OH methanol CO2 carbon dioxide H+ proton O2 oxygen Pt free active sites of platinum catalyst Pt3 -COH three sites of platinum catalyst for methanol

their importance in methanol transport. In their analysis, the electro-osmotic drag is considered as a convection effect and the diluted methanol moves with electro-osmotically dragged water molecules. Additionally, the cathode potential is influenced by the mixed potential phenomenon due to simultaneous methanol oxidation and oxygen reduction as well as poisoning of Pt catalysts by methanol oxidation products. Baxter et al. (1999) developed a one-dimensional model for a liquid-feed DMFC, focusing primarily on the anode catalyst layer. A major assumption in their study is that carbon dioxide remains dissolved in the liquid and their model for transport and electrochemical processes in the anode catalyst layer is, therefore, based on a single phase. On the other hand, Dohle et al. (2000) presented a one-dimensional model for the vapor-feed DMFC and the crossover phenomenon. The effects of methanol concentration on the cell performance were examined. Sundmacher et al. (2001) developed a lumped parameter DMFC model, and compared their results with the experimental data and examined the steady state and dynamic behavior of their model. This study considers the analysis and optimization of DMFC models. In the next section we will present a differential algebraic equation (DAE) model based on the work by Sundmacher et al. (2001) and will perform parameter estimation for this model. Section 3 deals with a number of parametric cases along with steady state optimization. This is followed by dynamic optimization in Section 4 and conclusions in Section 5. 2. The differential algebraic equation model for DMFC 2.1. The model equations Based on the work of Sundmacher et al. (2001), we summarize the DMFC DAE model in Table 1. The mathematical model considers the electrode kinetics, and mass and charge balances. The key assumptions for this model are summarized as: - the anode compartment is treated as a continuous stirred tank reactor (CSTR). The ohmic drops in current collectors and electric connections are negligible,

C. Xu et al. / Computers and Chemical Engineering xxx (2005) xxx–xxx

3

Table 1 The DAE model for DMFC in various components (taken from Sundmacher et al., 2001) Electrode kinetics Anode Cathode

r1 = k1 exp r5 = k5 exp


α F  1 RT

ηa

3 CL θPt cCH3 OH −


α F  5 RT

ηc




1 − exp −






1 F exp − ηa θPt3 -COH K1 RT

F ηc RT




p 3/2  O2

(1) (2)

p0

Mass balance Anode compartment

Anode catalyst layer

dcCH3 OH 1 F kLS AS CL − cCH3 OH ) − = (cCH (cCH3 OH − cCH ) 3 OH 3 OH dt τ Va

(3)

dcCO2 1 F kLS AS CL = (cCO − cCO2 ) − (cCO2 − cCO ) 2 2 dt τ Va

(4)

dcCL kLS AS AS AS CL = (cCH3 OH − cCH ) − CL nM CH3 OH − CL r1 3 OH dt VaCL Va Va

(5)

kLS AS dcCL AS CL = (cCO2 − cCO ) + CL r1 2 dt VaCL Va

(6)

Charge balance Anode Cathode

dηa 1 (−icell − 6Fr1 ) = dt Ca dηc 1 (−icell − 6F (r5 + nM = CH3 OH )) dt Cc

(7) (8)

- isothermal operation, - oxygen is fed in excess, i.e., oxygen conversion in the cathode compartment is negligible, and therefore oxygen mass balance is not required, - oxygen and carbon dioxide do not diffuse into the PEM, - water concentration is assumed to be constant (excess component in liquid mixture), - mass transport resistances in the catalyst layers are negligible as these are thin (10 ␮m) in comparison to the diffusion layers (100 ␮m) and PEM (200 ␮m), - mass transport coefficients of methanol and carbon dioxide in the anode diffusion layer are equal, and - the mixtures in the anode compartment and in the anode diffusion layer are treated as a pure liquid phase mixture, i.e., gas-phase formation, by release of carbon dioxide bubbles, is not taken into account.

LS means from liquid bulk to surface, VaCL the volume of the anode catalyst layer, nM CH3 OH the methanol crossover, Va and τ are the volume and mean residence time in the anode compartment, respectively. Equations (7) and (8) are the charge balance equations, where Ca and Cc the double layer capacity in anode and cathode areas, and icell is the current density. The eight equations in Table 1 form the DAE system to solve for the eight system variables, which can be used to calculate the overall cell voltage Ucell as follows:

In Table 1, Equations (1) and (2) give the rate expressions in both anode and cathode with rj as the rate of reaction j (j = 1, . . ., 5) at the electrode area, αj the charge transfer coefficient, and kj is the rate constant of reaction j. (Note: reactions (2)–(4) were implicitly eliminated in the model of Sundmacher et al., 2001.) Kj is the equilibrium constant of reaction j, θ Pt and θPt3 -COH are the surface coverages for platinum and adsorbed methanol (which occupies three sites), respectively. Also, ciCL is the concentration of component i (i = CH3 OH and CO2 ) in anode catalyst layer, p0 the standard pressure, pO2 the partial pressure of O2 , R the gas constant, T the absolute temperature, F the Faraday constant, and ηa and ηc are the anode and the cathode overpotential, respectively. Equations (3)–(6) are the mass balance equations, where AS is the cross-sectional electrode area, ci the component i concentration in the anode compartment, ciF the feed concentration of component i, kLS the mass transfer coefficient,

2.2. Parameter estimation

0 Ucell = Ucell − ηa + η c −

dM icell κM

(9)

0 is the open-circuit cell voltage, κM the PEM conwhere Ucell ductivity, and dM is the PEM thickness.

We found the given model parameters do not fit the reported data. Therefore, we estimated the model parameters. We describe the parameter estimation problem in the following manner (Arora & Biegler, 2004): min

{q,y}

s.t.

f (q, y) c(q, y) = 0 qL

≤q≤

(10)

qU ,

Here, f(q, y) is an objective function derived from a statistical basis, y the set of state variables, q the set of parameters to be estimated, and L and U denote the lower and upper bounds, respectively. Here, we assume:f : n → , c : n → m , q ∈ (n−m) , and n > m. The equality constraint c(q, y) = 0 represents the steady state equations of the DMFC

4

C. Xu et al. / Computers and Chemical Engineering xxx (2005) xxx–xxx

Table 2 Values for kLS (taken from Sundmacher et al., 2001) F cCH (mol/m3 ) 3 OH

(A/m2 )

ilim kLS (10−5 m/s)

125

250

500

2000

200 0.31

480 0.38

1150 0.47

4000 0.44

Table 3 Estimated parameters and their standard deviations with four sets of data Parameter (mol/(m2 s))

Mean value

Standard deviation

k1 k5

4.9284 × 10−7 1.6145 × 10−7

2.9856 × 10−7 0.8013 × 10−7

model in Section 2.1. The functions f(q, y) and c(q, y) are assumed to be at least twice differentiable. Equation (10) can be solved by successive quadratic programming (SQP). In our parameter estimation, the degrees of freedom, (n − m), are small in comparison to the total number of variables. Details on the trust region SQP algorithm can be found in Arora and Biegler (2004). We adopt a least square objective made Cal. ) up of differences between calculated voltage values (Ucell Exp. and the experimental ones (Ucell ) for 14 different current densities: f =

14 

2

Cal. k {Ucell (icell ) − Ucell (ikcell )} Exp.

(11)

k=1

Due to the complexity of the reactions involved in the anode catalyst layer, parameters such as reaction rate and equilibrium constants are very difficult to determine experimentally. Thus, parameter fitting method has been largely adopted in developing numerical models for the DMFC (Baxter et al., 1999; Bernardi & Verbrugge, 1992; Pisani, Murgia, Valentini, & D’Aguanno, 2002; Sundmacher et al., 2001). By analyzing the model equations and parameters, we found that two parameters are closely related to the calculation of anode overpotential, Ucell ; these are the rate constants k1 and k5 . Also, the calculated values for kLS are given in Table 2. To estimate parameters for k1 and k5 , we used the four sets of steady state polarization data (Sundmacher et al., 2001) for methanol feed concentrations of 125, 250, 500, and 2000 mol/m3 . By applying our in-house code (Arora & Biegler, 2004), we first obtained initial parameter estimates using the better behaved data sets for 500 and 2000 mol/m3 . With these initial parameter estimates, we re-ran the parameter estimation and obtained the parameter values listed in Table 3. The regions for k1 and k5 with 90, 95 and 99% confidence are plotted in Fig. 1.

Fig. 1. The elliptical regions for k1 and k5 with 90, 95, 99% confidence levels.

voltage versus current density, (see Fig. 2) for four different feed concentrations. Our curves generally agree well with the results provided in Sundmacher et al. (2001). With this model we now examine the effect of two design parameters through an optimization study. 3.1. Analysis of the DMFC performance We now examine several parametric studies for DMFC performances. For a performance measure, power density (Pcell ) is defined as: Pcell ≡ Ucell icell

(12)

Fuel efficiency (ηCH3 OH ) and energy efficiency (ηenergy ) are also defined as: ηCH3 OH ≡

reacted NCH 3 OH total NCH 3 OH

(13)

and ηenergy ≡

Pcell Ucell icell ≡ Pcombustion NCH3 OH !H

(14)

3. Steady state analysis Given the model and parameter estimates described in Section 2 we obtain the steady state polarization curves, i.e., cell

Fig. 2. Simulated steady state polarization curve for four methanol feed concentrations.

C. Xu et al. / Computers and Chemical Engineering xxx (2005) xxx–xxx

Here, NCH3 OH is the molar flux of methanol consumed by the cell, which can be calculated as: F NCH3 OH = V F (cCH − cCH3 OH ) 3 OH

(15)

where VF is the methanol feed volumetric rate and !H is the enthalpy of the following chemical reaction: CH3 OH + 23 O2 = CO2 + 2H2 O + !H

(16)

Our study is based upon the DAE model in Table 1 with F a focus on the two operational parameters, dM and cCH , 3 OH while maintaining all the other operational parameters fixed. 3.1.1. PEM Thickness, dM The effect of the parameter dM is described by the current density criteria as follows: (1) At a low current density, less than 1000 A/m2 , increasing dM is preferable as it leads to reduced methanol transport through the PEM. This allows us to achieve higher cell potentials, benefiting overall cell performance and energy efficiency. (2) At a current density greater than 1500 A/m2 , the ohmic loss becomes significant as can be observed in Equation (9), which outweighs the effect of reduced crossover. These competing effects are illustrated in Fig. 3 for methanol concentrations of 2000 and 3000 mol/m3 . At lower

Fig. 3. Polarization curve with two different PEM thicknesses (0.2 and 0.1 mm) for methanol feed concentrations of (a) 2000 mol/m3 and (b) 3000 mol/m3 .

5

current densities, the cell voltage for dM = 0.2 mm is higher than that for dM = 0.1 mm. However, at higher current densities, the result is reversed, which is in agreement with the above statement. Methanol crossover is also related to the methanol feed concentration. We note that the higher feed concentration (Fig. 3b) leads to greater methanol crossover, which ultimately leads to reduced energy efficiencies. As a result, the optimization of PEM thickness should always be considered in context of other parameters including the operational current density and the feed concentration. Therefore, an optimal dM not only improves energy efficiency but also satisfies the power density requirements as well. F 3.1.2. Methanol feed concentration, cCH 3 OH F If cCH3 OH is too low, as shown in Fig. 2, the DMFC cannot satisfy the high current density requirements during operaF tion. If cCH is too high, then the methanol crossover due 3 OH to diffusion will increase significantly and inversely affect the cell performance by decreasing the overall cell voltage F and energy efficiency. Therefore, an optimal cCH exists 3 OH between these two extremes.

3.2. Optimization strategy By using the sensitivity approach with Differential Algebraic Solver—Preconditioned Krylov (DASPK, numerical software in FORTRAN which can solve DAE systems and obtain the sensitivity simultaneously), we are able to study the nonlinearities as well as assess the dynamic properties of this system. Moreover, this approach then becomes a continuation method applied to the optimality conditions. While the maximization approach is approximate because of a graphical interpolation, it also provides us a much better picture of the surface and the location of both maximum and minimum stationary points. This would be more difficult with a typical optimization method. While an optimal thickness can be determined for DMFC design, it is more convenient to control the methanol concentration. We simulate the DMFC for two different PEM F thicknesses. For a given dM , we seek the optimal cCH that 3 OH achieves the highest energy efficiency at a given power density. If the maximum power density of the same value can be F achieved for a certain range of cCH , then increasing energy 3 OH efficiency should be our objective. That is, at a certain power density, the feed concentration that yields the highest energy efficiency is the desired optimal methanol feed concentration. A similar analysis can be found in Meyers and Newman (2002). Fig. 4 shows power density versus current density for the methanol feed concentrations of 2000 and 3000 mol/m3 . The maximum power densities achieved in both cases are almost identical. In Fig. 5, we plotted energy efficiency versus power density curves for six different methanol feed concentrations. As the methanol feed concentration increases, the maximum

6

C. Xu et al. / Computers and Chemical Engineering xxx (2005) xxx–xxx

3.3. DAE solver and sensitivity evaluation In the initial stage of our work, we formulated steady state optimization using standard nonlinear programming solvers and models in GAMS (Brooke, Kendrick, & Meeraus, 1992). However, as seen in Figs. 2 and 5, DMFC models exhibit extreme nonlinearity and even loss of continuity in the optimization space. As a result, it was very difficult to initialize the problem and obtain solutions this way. Instead, we opted to use a continuation approach for the steady state optimization that is based on the use of dynamic models and solvers, coupled with sensitivity analysis; these solutions were later validated with an NLP solver. Here, DASPK 3.0 (Li, Petzold, & Zhu, 2000) was used as the DAE solver and sensitivity package. DASPK can solve initial value problems of the form f (t, y, y˙ , q) = 0 (f, y, and y˙ are N-dimensional vectors, and q is the decision variable vector) via a combination of backward differentiation formulae (BDF) methods and different linear system solution methods. It also produces a sensitivity analysis of the parameters dy/dqi (t) concurrently within the DAE system. By integrating the DAE system up to steady state, we obtain steady state sensitivity values, dy/dqi , for the DMFC model. 3.4. Sensitivity of Ucell versus methanol feed concentration at various current densities

Fig. 4. Cell power curves with two different PEM thicknesses (0.2 mm (lower) and 0.1 mm (upper)) for methanol feed concentrations of (a) 2000 mol/m3 and (b) 3000 mol/m3 .

value of power density also increases. However, beyond the methanol feed concentration of 2000 mol/m3 , further increases in the feed concentration fail to increase the maximum power density. A number of feed concentrations can be found to satisfy a given power density requirement. Different feed concentrations, however, will result in different energy efficiencies. We have one specific feed concentration that provides the highest energy efficiency and, thus, the optimal feed concentration from the energy utilization point of view.

We performed a sensitivity analysis for the DAE model in Table 1 using DASPK 3.0. The sensitivity was calculated under different current densities. Typical results are shown in Fig. 6. From these curves, we can see that the sensitivity curve F of dUcell /dcCH crosses the abscissa, beyond which Ucell 3 OH F decreases as cCH3 OH increases. In Fig. 6, we note that there F are two stationary points (where dUcell /d(cCH ) = 0). By 3 OH examining the slopes at these points, we notice that the left stationary points are local maxima while the right stationary points are local minima. The optimal feed concentrations, therefore, correspond to the left stationary points derived from the curves. Note that stationary points close to 2000 mol/m3 should be avoided as they are local minima where methanol crossover dominates with a negative effect on the cell voltage. 3.5. Optimal methanol feed concentrations versus current density

Fig. 5. Energy efficiency vs. power density for various methanol feed concentrations.

We obtained the stationary points for the methanol feed concentrations from the sensitivity curves in Fig. 6 and plotted the relationship between the optimal methanol feed conF∗ centration (cCH ) and current density in Fig. 7. These opti3 OH mization results were also validated using the CONOPT NLP solver in GAMS (Brooke et al., 1992). From Fig. 7, we can see that the corresponding optimal methanol feed concentration increases with an increase in the current density. However, the curve levels off as the current density increases further, suggesting an optimal feed

C. Xu et al. / Computers and Chemical Engineering xxx (2005) xxx–xxx

7

Fig. 8. 3D polarization curves for different optimal methanol feed concentrations.

  max  F cCH

0

tf

F Ucell (cCH , y˙ , y, t)dt 3 OH

3 OH

s.t.

Fig. 6. Sensitivity analysis of cell voltage vs. methanol feed concentration under various current densities (from 10 to 1000 A/m2 ) for (a) low current densities and (b) high current densities.

F f (cCH , y˙ , y, t) = 0 3 OH

(17)

F,U F 0 ≤ cCH ≤ cCH . 3 OH 3 OH

F Here, cCH is our control variable, the integration of Ucell 3 OH over the time domain (0 to tf ) is the objective function, and the constraints are composed of the DAE model we presented F in Section 2.1 (denoted as f (cCH , y˙ , y, t) = 0 in Equation 3 OH (17)) along with the lower and upper bounds for the control variables. Note that for a fixed current density, this objective function corresponds to maximization of power density.

4.1. Methods and tools

Fig. 7. Optimal methanol feed concentrations vs. current density.

concentration of about 1200 mol/m3 at high current densities. We calculated the polarization data to verify the optimal methanol feed concentrations obtained from our sensitivity analysis. The 3D polarization plots are provided in Fig. 8. We observe that the cell voltage obtained becomes the maximum at the optimal feed concentration for a given current density.

4. Dynamic simulation and optimization Here, we consider an optimal dynamic profile for methanol feed concentration. The problem can be stated as:

In general, the solution of large-scale nonlinear optimal control problems can be obtained by discretizing either the control variables (Vassiliadis, Sargent, & Pantelides, 1994) or both the control and state variables (Cervantes, Wachter, Tutuncu, & Biegler, 2000). Both approaches transform the continuous optimal control problem into a nonlinear programming (NLP) problem. In our study, both control and state profiles are discretized. In order to formulate the continuous dynamic optimization as an NLP problem, the state variables are approximated as polynomials on finite elements. The type of polynomials and the distribution of the discrete time steps are closely related to implicit Runge–Kutta methods. In our study, we applied three point collocation at Radau points with piecewise constant control variables for each finite time element. A dynamic optimization tool called OptControl Centre (OCC; Jockenhoevel, Biegler, & Wachter, 2003) was used for the analysis. In order to apply the discretization to general optimal control problems in a user friendly way, the software package Optimal Control Code Generator for MAPLE (OCOMA) has been used. The basic structure of OCOMA is illustrated in Fig. 9. With OCOMA we only define the continuous dynamic optimization problem. This can be done easily by coding the

8

C. Xu et al. / Computers and Chemical Engineering xxx (2005) xxx–xxx

Fig. 9. Basic structure of the software package OCOMA for dynamic optimization.

Fig. 10. Simulated results for a constant methanol feed concentration of 1500 mol/m3 using OCC from 0 to 1000 time steps (Ac, anode catalyst layer; Fe, feed; cCo2, CO2 concentration; cMe, methanol concentration; nMe, methanol mass flux; Mem, membrane; An, anode; Me, methanol; ovpo, overpotential; Cc, cathode catalyst layer; c, concentration; r, reaction rate).

C. Xu et al. / Computers and Chemical Engineering xxx (2005) xxx–xxx

9

system consisting of DAEs, the inequality constraints, the bounds for variables, and the objective function in a simple syntax in MAPLE. The main tasks are conducted by MAPLE to formulate the collocation equations and to derive the first and second derivative information. In the second step, the continuous system is discretized using a three-point Radau implicit Runge–Kutta scheme. In the third step, a FORTRAN code consisting of the objective function, the constraints, and the symbolic gradients is generated. The resulting files are compiled and linked with the optimization code IPOPT to form an executable file. After the optimization, the results are stored in a file that can be read and visualized by MATLAB. 4.2. Dynamic simulation and optimization results To make use of dynamic optimization tools, we rearrange the anode model equations. From Sundmacher et al. (2001), the following relationship holds between θPt3 -COH and θPt : θPt3 -COH =

CL cCO 2

˜ 3 K3 K4 K 2

where ˜ 2 ≡ K2 exp K



3 θPt

F ηa RT

(18)

(19)

K2 , K3 , and K4 are the equilibrium constants of the anode reaction steps. Substituting Equations (18) and (19) into Equation (1), we obtain

α1 F r1 = k1 exp ηa RT 

CL cCO 4F CL 3 2 × cCH3 OH − θPt exp − ηa ˜ 3 K3 K4 RT K1 K 2 (20) The dynamic simulation and optimization are carried out for a fixed value of current density, e.g., 300 A/m2 . Similar to the steady state case, we designate the methanol feed concentration as the control variable and the cell voltage as the objective function. For each study, at a particular current density, we assume θ Pt is independent of time. With the above assumption, we end up with a form that is consistent with the general form of the Butler–Volmer equation, which is extensively used in studying the electro-kinetic behavior at the electrodes. Also, because the dynamics of the surface concentration appear to be slower than the liquid concentrations, assuming a steady state value for θ Pt seems justified. Some evidence for this assumption can be observed from the pulsed experiments by Zhou, Schultz, Peglow, & Sundmacher (2001). Here, we chose θ Pt as the value at the optimal steady state, corresponding to any given value of current density (note that θ Pt depends on current density), which

Fig. 11. Selected variables from the OCC simulation results for 1500 mol/m3 methanol feed concentration for comparison with the data in Zhou et al. (2001): (a) anode and cathode overpotential and (b) methanol concentration and methanol mass flux in the membrane in the anode catalyst layer.

is calculated from the steady state model. By making the above assumption, we observed the following form for the anode reaction kinetics expression:

α1 F ˜ r1 = k1 exp ηa RT 

CL cCO 4F CL 2 × cCH3 OH − exp − ηa ˜ 3 K 3 K4 RT K1 K 2 (21) 3 . Parameters for the simulation are modified where k˜ 1 ≡ k1 θPt accordingly in the dynamic case due to a change of catalyst (Sundmacher et al., 2001; Zhou et al., 2001). The cathode reaction rate constant k5 and mass transfer coefficient kLS are taken from Zhou et al. (2001). For the anode equilibrium constants, as they are insensitive to the cell voltage, we thus use the same values as for the steady state. The only fitted parameter in the dynamic case is the anode reaction rate constant k˜ 1 . Fig. 10 shows the dynamic simulation results when the methanol feed concentration jumps from 0 to 1500 mol/m3 at time zero and then remains unchanged in the entire time domain; this is shown as “cMeFe”. Some important variable profiles we obtained are those of anode overpotential (ovpoAn), cathode overpotential (ovpoCc), methanol mass flux in the membrane (nMeMem), and methanol concentration at the anode catalyst layer (cMeAc), which are plotted separately in Fig. 11 for comparison with Zhou et al.

10

C. Xu et al. / Computers and Chemical Engineering xxx (2005) xxx–xxx

Table 4 Model parameters used in the dynamic simulation and optimization Parameter

Units

Dynamic case

kLS

m/s mol/(m2 s) mol/(m2 s) mol/(m2 s) mol/(m2 s) mol/(m2 s) – – – –

8 × 10−6 4.9 × 10−7 1 1 × 10−2 1 1 × 10−8 1 0.5 1 0.7

k˜ 1 k2 k3 k4 k5 K1 K2 K3 K4

mal control profile is either at its bounds or exhibits singular arcs. This makes the numerical solution of this problem difficult. In most numerical approaches (Cervantes et al., 2000; Vassiliadis et al., 1994) singular problems lead to illconditioned second derivative matrices which lead to slow convergence or even nonunique profiles with arbitrary oscillations; these were observed in our preliminary study. To stabilize these solutions we modify the objective function by adding a quadratic regularization term 2 F λ(cCH − c1 ) ; which leads to a nonsingular control prob3 OH lem: 

(2001). We verified that the value of k˜ 1 used in our simulation can provide good agreement with the results from Zhou et al. (2001). Table 4 provides the parameters used in our dynamic simulation. Note that the methanol feed control profile in problem (17) appears linearly in Equation (3). Consequently, the opti-

 max  F cCH

0

tf

F Ucell (cCH , y˙ , y, t) 3 OH

3 OH

2

F −λ(cCH − c1 ) dt 3 OH

s.t.

(22)

F f (cCH , y˙ , y, t) = 0 3 OH F,U F 0 ≤ cCH ≤ cCH . 3 OH 3 OH

Fig. 12. Dynamic simulation and optimization results for icell = 300 A/m2 with λ = 1.0 × 10−6 . Similar results were obtained for different values of λ.

C. Xu et al. / Computers and Chemical Engineering xxx (2005) xxx–xxx

11

Fig. 13. Methanol feed concentration in the dynamic optimization for: (a) λ = 1.0 × 10−3 , (b) λ = 1.0 × 10−6 and (c) λ = 1.0 × 10−8 .

Here, we run three cases for λ = 1.0 × 10−3 , 1.0 × 10−6 , and 1.0 × 10−8 . Note that the regularization term is very small compared to Ucell and thus can be ignored upon the convergence of the program. c1 is a constant value of methanol feed concentration (we set it to be 1000 mol/m3 ). We also set the upper bound of the methanol feed concentration to be 3000 and 2000 mol/m3 for testing purposes. We obtained the same optimization results in both cases. This agrees with the findings of other research groups (Dohle et al., 2000; Mench et al., 2001) which conclude that it is unfavorable for the methanol feed concentration to be increased to any value greater than 2000 mol/m3 . The results shown in Fig. 12 are for 2000 mol/m3 upper bound case only. In Fig. 12 Ucell integral is our objective function. It almost increases linearly with time as Ucell remains constant most of the time. From these results, we found the optimal feed concentration is approximately 1000 mol/m3 at three different λ values. Also, we displayed our optimal objective values corresponding to different values of λ in Table 5. Note that the invariant values of the objective function at different λ confirm that our program has converged and the choice of regularization term is able to deal with this type of singular problem. To provide a closer look at methanol feeding, we plot them in Fig. 13 for these cases. At λ = 1.0 × 10−3 or 1.0 × 10−6 , the change of methanol concentration is insignificant, and at λ = 1.0 × 10−8 , the optimal methanol concentration remains in the range from 970 to 1015 mol/m3 . Note that the optimal methanol feed profile does not show a pulsed feeding strategy that was described in Sundmacher Table 5 Objective function value vs. λ λ

tf 0

Ucell dt

1.0 × 10−3 0.4105

1.0 × 10−6 0.4105

1.0 × 10−8 0.4105

et al. (2001) and Zhou et al. (2001) because of the correction in the singular control problem. This may be due to our assumption in simplifying the model and setting θ Pt to be constant. On the other hand, in their studies none of the cases were optimized either. As a result, we conclude that dynamic pulsed feeding may be unnecessary to achieve the optimal cell voltage for the cases we studied.

5. Conclusions In this paper, we have examined the DMFC model for both steady and dynamic state with a consideration of electrode kinetics and transport phenomena. The model is adopted from Sundmacher et al. (2001), which provides an excellent description of the experimental results. However, certain reaction constants are extremely difficult to quantify experimentally. Typically, they are fitted or estimated in the model to obtain desirable simulation results. A trust region SQP algorithm (Arora & Biegler, 2004) is applied to the DMFC model for parameter estimation, and the simulation results are in excellent agreement with the published data. We identified the key factors affecting the DMFC performance and energy utilization efficiency. From our analysis, we found that PEM thickness and methanol feed concentration are determining factors for the methanol crossover. We also performed a sensitivity analysis for the DMFC. Our preliminary sensitivity results for cell voltage versus methanol feed concentration at a series of current densities are encouraging. From the analysis, we were able to obtain the optimal methanol feed concentrations in steady state. Using these optimal feed concentrations, we were able to obtain the maximum cell voltages for a broad range of current densities. Generally, at current densities below 1000 A/m2 , the

12

C. Xu et al. / Computers and Chemical Engineering xxx (2005) xxx–xxx

optimal methanol feed concentration increases correspondingly with current density. However, this trend gets less and less significant as current density increases until the optimal feed concentration becomes approximately 1200 mol/m3 at high current density above 1000 A/m2 . This provides insight into the improvement of performance and energy efficiency by treating the methanol feed concentration as the decision variable. Finally, a dynamic simulation and optimization for the DAE system has also been carried out, which yield promising results favoring a constant feeding strategy of the methanol solution at the optimal feed concentration in order to improve the DMFC performance. Our simulation and optimization were based on a simple DMFC model. As more complex physical phenomena (e.g., two-phase transport in the anode, methanol distribution in the PEM) and intensive numerical schemes (e.g., PDE models) are incorporated into the model, the optimization problem will become more challenging. Efforts for simulation and optimization of complex and micro-scale models are ongoing. We intend to incorporate the fuel cell stack models for real-world applications, e.g., automobiles or power generation systems. We believe that the framework demonstrated in this paper will give advancement towards the optimization of real world fuel cell systems. Acknowledgement Pennsylvania Infrastructure Technology Alliance (administered via the Institute for Complex Engineered Systems at Carnegie Mellon University) is gratefully acknowledged for the financial support of this project. References Arora, N., & Biegler, L. T. (2004). A trust region SQP algorithm for equality constrained parameter estimation with simple parametric bounds. Computational Optimization and Applications, 28, 51–86. Baxter, S. F., Battaglia, V. S., & White, R. E. (1999). Methanol fuel cell model: Anode. Journal of the Electrochemical Society, 146(2), 437–447.

Bernardi, D. M., & Verbrugge, M. W. (1992). A mathematical model of the solid-polymer-electrolyte fuel-cell. Journal of the Electrochemical Society, 139(9), 2477–2491. Brooke, A., Kendrick, D., & Meeraus, A. (1992). GAMS: A user’s guide, release 2. 25. South San Francisco (CA): The Scientific Press. Cervantes, A. M., W¨achter, A., Tutuncu, R., & Biegler, L. T. (2000). A reduced space interior point strategy for optimization of differential algebraic systems. Computers & Chemical Engineering, 24(1), 39– 51. Dohle, H., Divisek, J., & Jung, R. (2000). Process engineering of the direct methanol fuel cell. Journal of Power Sources, 86(12), 469– 477. Jockenhoevel, T. B., Biegler, L. T., & Wachter, A. (2003). Dynamic optimization of the tennessee eastman process using the optcontrolcentre. Computers & Chemical Engineering, 27(11), 1513–1531. Li, S. T., Petzold, L., & Zhu, W. J. (2000). Sensitivity analysis of differential-algebraic equations: A comparison of methods on a special problem. Applied Numerical Mathematics, 32(2), 161–174. Mench, M. M., Wang, C. Y., & Thynell, S. (2001). An introduction to fuel cells and related transport phenomena. The International Journal of Transport Phenomena, 3(3), 1. Meyers, J. P., & Newman, J. (2002). Simulation of the direct methanol fuel cell—iii. Design and optimization. Journal of the Electrochemical Society, 149(6), A729–A735. Narayanan, S. R., Frank, H., Jeffries-Nakamura, B., Smart, M., Chun, W., Halpert, G., Kosek, J., & Cropley, C. (1995). Proton conducting membrane fuel cells i. The electrochemical society proceedings series. Pisani, L., Murgia, G., Valentini, M., & D’Aguanno, B. (2002). A working model of polymer electrolyte fuel cells—Comparisons between theory and experiments. Journal of the Electrochemical Society, 149(7), A898–A904. Ren, X. M., Springer, T. E., & Zawodzinski, T. A. (2000). Methanol transport through nafion membranes—Electro-osmotic drag effects on potential step measurements. Journal of the Electrochemical Society, 147(2), 466–474. Ren, X. Z., Zawodzinski, T. A., Uribe, F., Dai, H., & Gottesfeld, S. (1995). Proton conducting membrane fuel cells i. The electrochemical society proceedings series. Sundmacher, K., Schultz, T., Zhou, S., Scott, K., Ginkel, M., & Gilles, E. D. (2001). Dynamics of the direct methanol fuel cell (DMFC): Experiments and model-based analysis. Chemical Engineering Science, 56(2), 333–341. Vassiliadis, V. S., Sargent, R. W. H., & Pantelides, C. C. (1994). Solution of a class of multistage dynamic optimization problems. 2. Problems with path constraints. Industrial & Engineering Chemistry Research, 33(9), 2123–2133. Zhou, S., Schultz, T., Peglow, M., & Sundmacher, K. (2001). Analysis of the nonlinear dynamics of a direct methanol fuel cell. Physical Chemistry Chemical Physics, 3(3), 347–355.