An Improved Tank in Series Model for the Direct Methanol Fuel Cell

Abstract

This study presents an improved TSR model for Direct Methanol Fuel Cell (DMFC) with mass and energy balance equations taking into account gas evolution in the anode channels. The TSR model includes a modified charge balance equation suitable for potentiostatic fuel cell operation mode. Polarization curves calculated with the improved TSR model agree with experimental data from the literature. The developed TSR model is able to approximate two-dimensional concentration, temperature, and current density profiles in DMFC with parallel flow field.

1. Introduction

The direct methanol fuel cell is a potential energy source to replace batteries for portable electronics. The main attractive features of DMFC include high energy density and nearly zero recharge time. Numerous 1D and 3D models are available for DMFC in the literature. Kulikovsky presented a 1D isothermal model for DMFC [1]. Murgia et al. developed 1D multicomponent steady-state isothermal model for liquid-feed DMFC [2]. Hyun et al. predicted the DMFC performance for different anode flow field designs using computational simulations [3]. Danilov et al. developed the two-phase model with a submodel for interphase transfer. They used CFD-based DMFC model to visualize and analyze the gas evolution and flow patterns in the anode channels [4–6]. In contrast to the CFD technique with highly complex process description, lumped parameter models with a reactor network are computationally fast, approximating flow patterns and concentration profiles in fuel cells. Sundmacher et al. [7–10] developed CSTR and TSR models for DMFC operating in galvanostatic mode. They used a tank in series reactor model to investigate the influence of the anode flow fields on the steady-state and dynamic behavior of DMFC with parallel, spot, and rhomboidal designs. Our previous studies focused on TSR model development for solid oxide fuel cells [11, 12]. The objective of this research is to modify this TSR model for DMFC operating in potentiostatic mode.

2. Model Formulation

2.1. TSR Model for Cocurrent Flow DMFC

Fuel and oxidant flow inside the fuel cell with parallel channels can be approximated as a reactor network complying with a tank in series model. Figure 1 illustrates the application of a tank in series reactor (TSR) model for cocurrent flow mode in a planar fuel cell with parallel flow field design.

Figure 1 

TSR model for cocurrent flow DMFC.

The developed TSR model for DMFC is based on the following assumptions.(i)Each anode and cathode compartment is treated as a TSR.(ii)Electrochemical reactions occur at the electrode/membrane interfaces.(iii)Fuel cell operates under the potentiostatic operating mode and constant pressure.(iv)Negligible contact resistances between components.(v)Negligible evaporation and condensation in cathode channels.(vi)Uniform distribution of fuel and oxidant in channels.

Balance equations are composed for channel compartment and catalyst layer in each tank. For fuel cells with cocurrent flow mode, the mass balance equations for the anode and cathode channels are written as follows: 𝑉^{𝐴𝑖,𝑗}𝑑𝐶^{(𝑘)𝐴,𝑖,𝑗}𝑑𝑡=𝐶^{(𝑘)𝐴,𝑖,𝑗−1}𝐹𝑘=^{𝐴𝑖,𝑗−1}M−e𝐶OH^{(𝑘)𝑖,𝑗},𝐹CO^{𝐴𝐿,𝑖,𝑗}₂−𝑁^{(𝑘)𝐴,𝑖,𝑗},𝑉,H₂O,_{𝐶𝑖,𝑗}𝑑𝐶^{(𝑘)𝐶,𝑖,𝑗}𝑑𝑡=𝐶^{(𝑘)𝐶,𝑖,𝑗−1}𝐹^{𝐶𝑖,𝑗−1}𝑘−=𝐶O^{(𝑘)𝐶,𝑖,𝑗}₂𝐹^{𝐶𝑖,𝑗}−𝑁^{(𝑘)𝐶,𝑖,𝑗},,H₂O,N₂.(1) Gas evolution in the anode channels is defined by the following balance equation for the gas phase: 𝑉^{𝐴𝑖,𝑗}𝑑𝜑^{𝐴𝐺,𝑖,𝑗}𝑑𝑡=𝜑^{𝐴𝐺,𝑖,𝑗−1}𝐹^{𝐴𝑖,𝑗−1}−𝜑^{𝐴𝐺,𝑖,𝑗}𝐹^{𝐴𝑖,𝑗}−𝑟^{𝐴desorp,𝑖,𝑗}𝜌^𝐴_mol,𝐿,(2) In cocurrent flow mode, the fuel and oxidant outgoing streams from tank 𝑖,𝑗−1, are the inlet streams into tank 𝑖,𝑗. Therefore, the inlet flow rate for each tank depends on the outgoing flow rate of the previous tank: 𝐹^{𝐴𝑖,𝑗}=𝐹^{𝐴𝑖,𝑗−1}−𝑁^{MeOH𝐴,𝑖,𝑗}+𝑁^CO2_{𝐴,𝑖,𝑗}+𝑁^H2O_{𝐴,𝑖,𝑗}𝐹_{𝐶𝑖,𝑗}𝐶^{MeOH𝐴,𝑖,𝑗}+𝐶^CO2_{𝐴,𝑖,𝑗}+𝐶^H2O_{𝐴,𝑖,𝑗},=𝐹^{𝐶𝑖,𝑗−1}−𝑁^{O2𝐶,𝑖,𝑗}+𝑁^H2O_{𝐶,𝑖,𝑗}𝐶^{O2𝐶,𝑖,𝑗}+𝐶^H2O_{𝐶,𝑖,𝑗}+𝐶^N2_{𝐶,𝑖,𝑗},(3) Here 𝑉^{𝐴𝑖,𝑗}, 𝑉^{𝐶𝑖,𝑗} represent the volume of anode and cathode channels in tank (𝑖,𝑗), 𝑉^{𝐴𝑖,𝑗}=𝑉^𝐴/(𝑛_𝑖𝑛_𝑗), 𝑉^{𝐶𝑖,𝑗}=𝑉^𝐶/(𝑛_𝑖𝑛_𝑗); 𝑆_𝑖,𝑗 electrode area in tank (𝑖,𝑗), 𝑆_𝑖,𝑗=𝑆_cell/(𝑛_𝑖𝑛_𝑗).

Component balance equations for the anode and cathode catalyst layers reflect changes due to consumption and production of species via reaction and mass transfer: 𝑉^𝐴_{cat,𝑖,𝑗}𝑑𝐶^(𝑘)𝐴,_{cat,𝑖,𝑗}𝑑𝑡=𝑁^{(𝑘)𝐴,𝑖,𝑗}𝑉−𝑟^{(𝑘)𝐴,𝑖,𝑗},𝑘=MeOH,CO₂,H₂O,(4)_{𝐶cat,𝑖,𝑗}𝑑𝐶^(𝑘)𝐶,_{cat,𝑖,𝑗}𝑑𝑡=𝑁^{(𝑘)𝐶,𝑖,𝑗}−𝑟^{(𝑘)𝐶,𝑖,𝑗},𝑘=O₂,H₂O.(5) Here 𝑉^𝐶_{cat,𝑖,𝑗}, 𝑉^𝐴_{cat,𝑖,𝑗} correspond to the volume of anode and cathode catalyst layer in tank (𝑖,𝑗), 𝑉^𝐶_{cat,𝑖,𝑗}=𝑉^𝐶_cat/(𝑛_𝑖𝑛_𝑗), 𝑉^𝐴_{cat,𝑖,𝑗}=𝑉^𝐴_cat/(𝑛_𝑖𝑛_𝑗).

The energy balances for the channels and the MEA structure in tank 𝑖,𝑗 are written as follows: 𝑉^{𝐴𝑖,𝑗}𝜌^𝐴_mol,𝐿𝐶^{𝐴𝑝,𝐿}𝑑𝑇^{𝐴𝐿,𝑖,𝑗}=𝜌𝑑𝑡_𝐴mol,𝐿+𝑞_{𝐴𝐿,𝑖,𝑗}𝐹^{𝐴𝑖,𝑗−1}Δℎ^{𝐴𝐿,𝑖,𝑗−1}−𝜌^𝐴_mol,𝐿𝐹^{𝐴𝑖,𝑗}Δℎ^{𝐴𝐿,𝑖,𝑗}𝑉+𝑞^𝐴_{coll,𝑖,𝑗},(6)_{𝐶𝐺,𝑖,𝑗}𝜌^𝐶_mol,𝐺𝐶^{𝐶𝑝,𝐺}𝑑𝑇^{𝐶𝐺,𝑖,𝑗}=𝜌𝑑𝑡_𝐶mol,𝐺+𝑞_{𝐶𝐺,𝑖,𝑗}𝐹^{𝐶𝑖,𝑗−1}Δℎ^{𝐶𝐺,𝑖,𝑗−1}−𝜌^𝐶_mol,G𝐹^{𝐶𝑖,𝑗}Δℎ^{𝐶𝐺,𝑖,𝑗}𝑉+𝑞^𝐶_{coll,𝑖,𝑗},(7)_{𝑆𝑖,𝑗}𝜌𝐶^𝑆𝑝𝑑𝑇^{𝑆𝑖,𝑗}𝑑𝑡=𝑞^{𝑆𝑖,𝑗}−𝑞^{𝐴𝐿,𝑖,𝑗}−𝑞^{𝐶𝐺,𝑖,𝑗}.(8) Here 𝑞^{𝐶𝐺,𝑖,𝑗}, 𝑞^{𝐴𝐿,𝑖,𝑗} represent the convective heat transferred from the MEA to the channel in tank 𝑖,𝑗; 𝑞^𝐶_{coll,𝑖,𝑗}, 𝑞^𝐴_{coll,𝑖,𝑗} represent heat source due to heat transfer from channel to the collector; 𝑞^{𝑆𝑖,𝑗} represents the heat generation in MEA in tank 𝑖,𝑗.

The charge balance equations for anode and cathode electrode/membrane interfaces are given by 𝐶^𝐴𝑑𝑙𝑑𝜂^{𝐴𝑖,𝑗}𝑑𝑡=𝐼^{cell𝑖,𝑗}𝐶_𝐶𝑑𝑙−𝐼^{𝐴𝑖,𝑗},𝑑𝜂^{𝐶𝑖,𝑗}𝑑𝑡=−𝐼^{cell𝑖,𝑗}−𝐼^{𝐶𝑖,𝑗}+𝐼^par_𝑖,𝑗.(9)

For j = 1 tanks, representing the inlet sections of the anode and cathode channels, we obtain the following balance equations: 𝑉^𝐴𝑖,1𝑑𝐶^{(𝑘)𝐴,𝑖,1}𝑑𝑡=𝐶^{(𝑘)𝐴,𝑖,in}𝐹𝑘=^𝐴𝑖,_in−M𝐶eOH^{(𝑘)𝐴,𝑖,1},𝐹CO^𝐴𝑖,1₂−𝑁^{(𝑘)𝐴,𝑖,1},𝑉,H₂O,(10)_𝐶𝑖,1𝑑𝐶^{(𝑘)𝐶,𝑖,1}𝑑𝑡=𝐶^{(𝑘)𝐶,𝑖,in}𝐹^𝐶𝑖,_in𝑘−=𝐶O^{(𝑘)𝐶,𝑖,1}₂𝐹^𝐶𝑖,1−𝑁^{(𝑘)𝐶,𝑖,1},𝑉,H₂O,N₂,(11)_𝐴𝑖,1𝑑𝜑^{𝐴𝐺,𝑖,1}𝑑𝑡=𝜑^{𝐴𝐺,𝑖,in}𝐹^𝐴𝑖,_in−𝜑^{𝐴𝐺,𝑖,1}𝐹^𝐴𝑖,1−𝑟^{𝐴desorp,𝑖,1}𝜌^𝐴_mol,𝐿𝐹,(12)_𝐴𝑖,1=𝐹^𝐴𝑖,_in−𝑁^{MeOH𝐴,𝑖,1}+𝑁^CO2_𝐴,𝑖,1+𝑁^H2O_𝐴,𝑖,1𝐹𝐶^{MeOH𝐴,𝑖,1}+𝐶^CO2_𝐴,𝑖,1+𝐶^H2O_𝐴,𝑖,1,(13)_𝐶𝑖,1=𝐹^𝐶𝑖,_in−𝑁^{O2𝐶,𝑖,1}+𝑁^H2O_𝐶,𝑖,1𝐶^{O2𝐶,𝑖,1}+𝐶^H2O_𝐶,𝑖,1+𝐶^N2_𝐶,𝑖,1,(14) where 𝐹^𝐴𝑖,_in=𝐹^𝐴_in/𝑛_𝑖 is the inlet anode volumetric flow rate for (𝑖,𝑗) tank with j = 1; 𝐹^𝐶𝑖,_in=𝐹^𝐶_in/𝑛_𝑖 is the inlet cathode volumetric flow rate for (𝑖,𝑗) tank with j = 1.

For 𝑗=𝑛_𝑗 tanks, corresponding to the outlet section of fuel and oxidant gas channels, we define mean outlet variables as follows: 𝐹^𝐴_out=^𝑛𝑖_𝑖=1𝐹^{𝐴𝑖,𝑛}_𝑗,𝐶^MeOHout=∑^{𝑛𝑖𝑖=1}𝐶^MeOH_{𝐴,𝑖,𝑛𝑗}𝐹^{𝐴𝑖,𝑛}_𝑗𝐹^𝐴_out𝑇_𝐴𝐿,out,=∑^{𝑛𝑖𝑖=1}𝑇^{𝐴𝐿,𝑖,𝑛}_𝑗𝐹^{𝐴𝑖,𝑛}_𝑗𝐹^𝐴_out,𝐹^𝐶_out=^𝑛𝑖_𝑖=1𝐶_O2out𝐹^{𝐶𝑖,𝑛}_𝑗,=∑^{𝑛𝑖𝑖=1}𝐶^O2_{𝐶,𝑖,𝑛𝑗}𝐹^{𝐶𝑖,𝑛}_𝑗𝐹^𝐶_out,𝑇^𝐶𝐺,_out=∑^{𝑛𝑖𝑖=1}𝑇^{𝐶𝐺,𝑖,𝑛}_𝑗𝐹^{𝐶𝑖,𝑛}_𝑗𝐹^𝐶_out.(15) The developed mathematical model includes the following phenomena:(i)electrochemical oxidation of methanol at the anode electrode/membrane interface,(ii)electrochemical reduction of oxygen at the cathode electrode/membrane interface,(iii)charge balances at anode and cathode electrode/membrane interfaces,(iv)energy balances in gas channels and in MEA.

2.2. The Electrochemistry Submodel

2.2.1. Electrode Current

The rate of electrochemical reactions is defined by Butler-Volmer equation. The anode reaction rate is I^{𝐴𝑖,𝑗}=𝐼^𝐴0𝐶𝐶^{MeOHcat,𝑖,𝑗MeOHref}^0.5×exp𝛼^𝐴𝐴𝐹𝑅𝑇_𝑖,𝑗𝜂^{𝐴𝑖,𝑗}−exp−𝛼_𝐴𝐶−𝜂^𝐴_eq𝐹𝑅𝑇_𝑖,𝑗𝜂^{𝐴𝑖,𝑗}−𝜂^𝐴_eq,(16) where 𝐼^𝐴0 is the anode exchange current density.

The cathode reaction rate is 𝐼^{𝐶𝑖,𝑗}=𝐼^𝐶0𝐶^O2_{cat,𝑖,𝑗}𝐶^O2_refexp𝛼^𝐶𝐴𝐹𝑅𝑇_𝑖,𝑗𝜂^{𝐶𝑖,𝑗}−exp−𝛼_𝐶𝐶−𝜂^𝐶_eq𝐹𝑅𝑇_𝑖,𝑗𝜂^{𝐶𝑖,𝑗}−𝜂^𝐶_eq,(17) where 𝐼^𝐶0 is the cathode exchange current density.

Parasitic current is calculated as methanol molar flux from anode catalyst layer to cathode catalyst layer with a linear profile in the membrane: 𝐼^par_𝑖,𝑗=𝑛^𝐴𝑒𝐹𝐷^MeOH𝑚𝛿_𝑚𝐶^{MeOH𝐴,cat,𝑖,𝑗}−0.(18)

Source terms in the component balance equations reflect changes due to the consumption or production of species via reaction or mass transfer. For component balance equations in anode catalyst layer (4), the rate of electrochemical reaction is given by 𝑟^MeOH_{𝐴,𝑖,𝑗}=𝑆_𝑖,𝑗𝜈^MeOH𝑛^𝐴𝑒𝐹𝐼^{𝐴𝑖,𝑗}𝑟_{CO2𝐴,𝑖,𝑗},=𝑆_𝑖,𝑗𝜈^CO2𝑛^𝐴𝑒𝐹𝐼^{𝐴𝑖,𝑗}𝑟_{H2O𝐴,𝑖,𝑗},=𝑆_𝑖,𝑗𝜈^H2O𝑛^𝐴𝑒𝐹𝐼^{𝐴𝑖,𝑗}−𝑆_𝑖,𝑗𝑛_{𝑑,𝑖,𝑗}𝐼^{𝐴𝑖,𝑗}𝐹,(19) where 𝑛_{𝑑,𝑖,𝑗} represents the electro-osmotic drag coefficient.

For balance equations in cathode catalyst layer (5), the rate of electrochemical reaction is 𝑟^O2_{𝐶,𝑖,𝑗}=𝑆_𝑖,𝑗𝜈^O2𝑛^𝐶𝑒𝐹||𝐼^{𝐶𝑖,𝑗}||,(20) The source term for water on the cathode side includes electrochemical reaction and transfer from anode side 𝑟^H2O_{𝐶,𝑖,𝑗}=𝑆_𝑖,𝑗𝜈^H2O𝑛^𝐶𝑒𝐹||𝐼^{𝐶𝑖,𝑗}||+𝑆_𝑖,𝑗𝑛_{𝑑,𝑖,𝑗}||𝐼^C_𝑖,𝑗||𝐹,(21) An empirical equation for calculating electro-osmotic drag coefficient 𝑛_𝑑 is taken from Ren et al. [13]. For the inert nitrogen component, the corresponding source term is equal to zero (𝑟^N2_{𝐶,𝑖,𝑗}=0).

2.2.2. Electrolyte Current

For fuel cells operating under potentiostatic mode (𝐸_cell=const), the current density in PEM electrolyte media is estimated from the voltage equation: 𝐼^cell_𝑖,𝑗=𝐸^OCV𝑖,𝑗−𝐸_cell−𝜂^𝐴_{act,𝑖,𝑗}+𝜂^Cact_,𝑖,𝑗𝑅^Ohmic_𝑖,𝑗,(22) where 𝜂^𝐴_{act,𝑖,𝑗} is the anode activation overpotential, 𝜂^𝐴_{act,𝑖,𝑗}=𝜂^{𝐴𝑖,𝑗}−𝜂^𝐴_eq; 𝜂^𝐶_{act,𝑖,𝑗} is the cathode activation overpotential, 𝜂^𝐶_{act,𝑖,𝑗}=𝜂^{𝐶𝑖,𝑗}−𝜂^𝐶_eq; 𝜂^𝐴_eq, 𝜂^𝐶_eq is the anodic and cathodic equilibrium potential differences; 𝜂^{𝐴𝑖,𝑗}, 𝜂^{𝐶𝑖,𝑗} is the anode and cathode potential differences at the electrode/membrane interfaces; 𝑅^Ohmic_𝑖,𝑗 is the ohmic resistance.

As shown in Appendix B, voltage (22) is applicable for calculating the current in PEM under a linear membrane phase potential profile. Membrane conductivity is given by the following empirical relation [14]: 𝜎_𝑚=0.0005139𝐶_𝑚−0.000326exp12681303−𝑇1.(23)

2.3. Heat Transfer Submodel

The source terms in energy balance equations for anode (6) and cathode channels (7) are defined as follows: 𝑞𝑞^{𝐴𝑖,𝑗}=𝛼^𝐴𝑆_𝑖,𝑗𝑇^{𝑆𝑖,𝑗}−𝑇^{𝐴𝑖,𝑗},(24)_{𝐶𝑖,𝑗}=𝛼^𝐶𝑆_𝑖,𝑗𝑇^{𝑆𝑖,𝑗}−𝑇^{𝐶𝑖,𝑗}.(25) The source term in energy balance equation (8) for MEA is 𝑞^{𝑆𝑖,𝑗}=𝛿_𝑚𝑆_𝑖,𝑗𝑅^Ohmic_𝑖,𝑗𝐼^cell_𝑖,𝑗²−Δ𝐻_𝑅𝑁^H2O_{𝐶,𝑖,𝑗}+𝑆_𝑖,𝑗𝐸_cell𝐼^cell_𝑖,𝑗+𝑆_𝑖,𝑗𝛼^𝐶+_𝑘𝐶^𝑝(𝑘)𝑁^{(𝑘)𝐶,𝑖,𝑗}𝑇^{𝐶𝑖,𝑗}+𝑆_𝑖,𝑗−𝑇^{𝑆𝑖,𝑗}𝛼^𝐴+_𝑘𝐶^𝑝(𝑘)𝑁^{(𝑘)𝐴,𝑖,𝑗}𝑇^{𝐴𝑖,𝑗}−𝑇^{𝑆𝑖,𝑗}.(26) Additional source terms in energy equations (6) and (7) take into account heat transfer from channels to the anode and cathode collectors 𝑞𝑞^𝐴_{coll,𝑖,𝑗}=𝛼^𝐴_coll𝑆_𝑖,𝑗𝑇_coll−𝑇^{𝐴𝑖,𝑗},(27)_{𝐶coll,𝑖,𝑗}=𝛼^𝐶_coll𝑆_𝑖,𝑗𝑇^coll−𝑇^{𝐶𝑖,𝑗}.(28) Convective heat transfer coefficients 𝛼^𝐶 and 𝛼^𝐴 are calculated from empirical correlations for laminar heat transfer in channels. Heat transfer coefficients 𝛼_coll are found by a technique proposed by Siegel et al. [15].

2.4. Mass Transfer Submodel

Component molar fluxes for species on the cathode side are defined by mass transfer in gas phase: 𝑁^{(𝑘)𝑖,𝑗}=𝑆_𝑖,𝑗𝜌_mol,𝐺𝛽^(𝑘)𝐺,_eff𝑦^{(𝑘)𝑖,𝑗}−𝑦^(𝑘)_{cat,𝑖,𝑗},𝑘=O²,H₂O,(29) where 𝛽^(𝑘)𝐺,_eff is the effective mass transfer coefficient in gas phase for k component; 𝑦^{(𝑘)𝑖,𝑗},𝑦^(𝑘)_{cat,𝑖,𝑗} are mole fraction of k component in (𝑖,𝑗) tank in channel and catalyst surface, respectively.

The multicomponent mixture in gas phase includes carbon dioxide, methanol and water in equilibrium with the liquid phase. Component molar fluxes for species on the anode side are defined by mass transfer in liquid phase: 𝑁^{(𝑘)𝑖,𝑗}=𝑆_𝑖,𝑗𝜌𝑘=^mol,𝐿M𝛽eO^{(𝑘)𝐿,eff}H,𝑥CO^{(𝑘)𝑖,𝑗}₂−𝑥^(𝑘)_{cat,𝑖,𝑗},,H₂O.(30)

2.5. Gas Evolution

Gas evolution results from interphase mass transfer of carbon dioxide in the anode channels. The source term in balance equation for gas phase (2) describes the rate of desorption and absorption processes in the anode channels. The conventional submodel for estimating the source term is based on mass transfer equation 𝑟_{desorp,𝑖,𝑗}=𝑆_𝑖,𝑗𝑘_𝑉𝜌_mol,𝐺𝑦^CO2𝑖,𝑗−𝑦^CO2_{sat,𝑖,𝑗},(31) where 𝑘_𝑉 is an empirical volumetric mass transfer coefficient.

It is generally assumed that the multicomponent gas phase is in equilibrium with the liquid phase. Using the equilibrium condition, Danilov et al. [4–6] proposed the following equation for estimating the source term in anode channels: 𝑟_{desorp,𝑖,𝑗}=𝛾𝜓+𝑟^CO2_{𝐴,𝑖,𝑗},(32) where 𝛾 is the splitting factor; 𝜓 is the coefficient. The derivation of the auxiliary equation for coefficient 𝜓 is given in [4–6]. The local splitting factor 𝛾 is found from solving the equilibrium flash equation [4–6]. It should be noted that new submodel (32) corresponds to an equilibrium model of multicomponent mass transfer between liquid and gas in the anode channels [4–6].

According to the physical meaning of the mass balance equations, molar component concentration is defined as mixture concentrations for gas-liquid flow in channels. Component mole fraction in liquid phase (𝑥^(𝑘)) can be expressed through mixture concentration (𝐶^(𝑘)) with known gas volume fraction (𝜑_𝐺). As shown in Appendix C, the condition of thermodynamic equilibrium in the multicomponent gas-liquid mixture gives the next relationship for the specie mole fraction in liquid phase 𝑥^(𝑘)=𝐶^(𝑘)𝐾^(𝑘)𝜑_𝐺𝜌_mol,𝐺+1−𝜑_𝐺𝜌_mol,𝐿,(33) where 𝐾^(𝑘) is the distribution coefficient for k-component, 𝐾^(𝑘)=𝑦^(𝑘)/𝑥^(𝑘).

3. Results and Discussion

Modeling and simulations are valuable tools for improving our understanding of the heat and mass transfer processes in fuel cells. To validate the improved TSR model, we used experimental data reported by Murgia et al. [2] for a 25 cm² DMFC operating with 1.5 M and 1 M aqueous methanol feeds at 90°C. Operating conditions for this DMFC are listed in Table 1.

Table 1 

Operating conditions for the DMFC.

The number of tanks is set as 𝑛_𝑖=4 and 𝑛_𝑗=4, corresponding to a TSR model with 272 nonlinear coupled first-order ordinary differential equations. The accepted number of tanks corresponds to the laminar flow regime in the channels with parallel flow field design. The developed TSR model was implemented in MATLAB, and it was initialized with the feed composition and temperature. The geometry and electrochemical parameters of DMFC are listed in Table 2.

Table 2 

Geometry and electrochemical parameters for DMFC.

Figure 2 compares the experimental and calculated polarization curves for DMFC with operating conditions in Table 1. The solid line represents a DMFC model prediction with the improved TSR model.

Figure 2 

Comparison of the polarization curves predicted by TSR model with experimental data [2].

Figure 3 displays histograms with steady-state simulation results predicted by the TSR model for DMFC with 1.0 M feed of fuel at the anode and pure oxygen feed at the cathode. It should be noted that each tank is characterized by complete mixing, and distribution of concentration, temperature, and current density is presented by step change for tanks in series. Concentration of methanol in channel is decreased due to the anode electrochemical reaction. For the anode channels, carbon dioxide is the product of the anode electrochemical reaction. Gas content in the anode channels is increased from inlet to the outlet following the trend in molar CO₂ concentration. The two-dimensional histograms indicate a strong coupling between the reactants and the current density distribution. The predicted current density profile is following the trend in methanol concentration in anode channels. Murgia et al. [2] obtained experimental polarization data by circulating aqueous methanol from a reservoir where CO₂ was released to the atmosphere. The inlet CO₂ concentration in the liquid fuel in Table 1 corresponds to the gas-liquid equilibrium condition in the fuel reservoir. The predicted temperature profile in MEA is following the main trend in the current density distribution.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 3 

Histograms with steady-state simulation results predicted by TSR model for DMFC with coflow mode of fuel and oxidant (𝐸cell=0.2V). Numerical modelling of fuel cell with anode 1.0 M feed of fuel and pure oxygen feed at the cathode. Temperature variation 𝛿𝑇𝑆𝑖,𝑗=𝑇𝑆𝑖,𝑗−𝑇𝐶in. Methanol molar concentration variation 𝛿𝐶MeOH𝑖,𝑗=𝐶MeOH𝑖,𝑗−𝐶MeOHin. Carbon dioxide molar concentration variation 𝛿𝐶CO2𝑖,𝑗=𝐶CO2𝑖,𝑗−𝐶CO2in.

Figure 4 presents histograms with steady-state simulation results predicted by the TSR model for DMFC with 1.5 M feed of fuel at the anode and pure oxygen feed at the cathode.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 4 

Histograms with steady-state simulation results predicted by TSR model for DMFC with coflow mode of fuel and oxidant (𝐸cell=0.2V). Numerical modelling of fuel cell with anode 1.5 M feed of fuel and pure oxygen feed at the cathode. Temperature variation 𝛿𝑇𝑆𝑖,𝑗=𝑇𝑆𝑖,𝑗−𝑇𝐶in. Methanol molar concentration variation 𝛿𝐶MeOH𝑖,𝑗=𝐶MeOH𝑖,𝑗−𝐶MeOHin. Carbon dioxide molar concentration variation 𝛿𝐶CO2𝑖,𝑗=𝐶CO2𝑖,𝑗−𝐶CO2in.

Gas management greatly influences the performance of the fuel cell. On the anode side, carbon dioxide is produced by electrochemical oxidation of the methanol. Inefficient removal of CO₂ bubbles may block anode channels and decrease efficiency of the fuel cells due to limited mass transport and maldistribution of reactants. According to the conventional submodel, the interface flux is completely dependent on empirical mass transfer coefficient (𝑘_𝑉). The predicted gas content in the anode channels can vary from 0 to 90% depending on the value of the mass transfer coefficient. Simulation results reveal that high values of the mass transfer coefficient (𝑘_𝑉>6000s⁻¹) predict gas evolution in the anode channels with high gas content (𝜑>40%). For low mass transfer coefficients (𝑘_𝑉<6000s⁻¹), on the other hand, 𝜑<40%. The volumetric mass transfer coefficient 𝑘_𝑉=6000s⁻¹ used in gas evolution submodel (31) provides the condition 𝜑≈40% in the anode channels, which is close to the equilibrium condition given by the flash equation.

As shown in Figure 5, empirical model (31) and new submodel (32) for gas evolution predict similar trends in the gas content distribution. In contrast to nonequilibrium submodel (31) with empirical mass transfer, new rate expression (32) determines the gas content under equilibrium conditions without empirical coefficients.

(a)

(b)

Figure 5 

Histograms with gas content in anode channels predicted by TSR model for DMFC with (a) empirical mass transfer (31) and (b) new gas evolution submodel (32) (𝐸cell=0.2V).

Comparison of DMFC performance with different fuel composition is given in Table 3. Changing the fuel composition from 1 M to 1.5 M gives a small improvement in fuel cell performance. The mass transfer coefficient is presented in the form of limiting current density. Average current density approaches the value of limiting current densities under cell voltage 𝐸_cell=0.2V. The possible way of DMFC performance improvement is to intensify mass transfer from the channels to the catalyst layer.

Table 3 

Single DMFC cell simulation results (𝐸_{𝑐𝑒𝑙𝑙}=0.2𝑉).

4. Conclusions

This study presents a TSR model for DMFC operating in potentiostatic mode. The developed TSR model includes mass and energy balance equations in anode and cathode compartments together with gas evolution in anode channels. Modified charge balance equations are defined as interface boundary conditions. The possibility of evaluating two-dimensional profiles in the DMFC is one of the main advantages of the developed TSR model. The simulation results indicate the strong coupling between concentration, temperature and current density distribution in coflow DMFC. Taking into account the assumption of uniform distribution of fuel and oxidant in channels, the TSR model predicts the limiting performance of the DMFC under the given flow regime (number of tanks), mass transfer and catalyst activity. The improved TSR model allows studying the influence of different parameters on the DMFC performance. Results can be used to better understand and investigate the effects of various parameters and operating conditions on DMFC performance.

Appendices

A. Derivation of Charge Balance Equation

The electric potential fields are governed by the charge conservation equations. The charge balance at the interface between electron-conducting and ion-conducting media is given by [16] 𝜕𝑄𝜕𝑡+∇⋅𝑖^s=𝐼₁−𝐼₂,(A.1) where 𝐼₁ is the current in electron-conducting media normal to the boundary, 𝐼₂ is the current in ion-conducting media normal to the boundary, 𝑖_𝑠 is the superficial current density, 𝑄is the charge. For the TSR model, the net charge flux is zero (∇⋅𝑖_𝑠=0). The interfaces between ionic and electronic media behave like a capacitor in which the charge density is a function of potential difference across the double layer. Charge or discharge rate at the electrode-electrolyte double layer can be represented as 𝜕𝑄𝜕𝑡=𝐶^dl𝜕𝜂𝜕𝑡.(A.2)

For the difference of potential of the electron conducting media 𝜙 and the potential of the electrolyte phase 𝜙_𝑚, the following charge conservation equation is valid: 𝜕𝜂𝜕𝑡=𝐶1_dl𝐼₁−𝐼₂,(A.3) where 𝜂 is potential differences or overpotential, 𝜂=𝜙−𝜙_𝑚.

The improved fuel cell models with new charge balance equation provide a better understanding of main phenomena governing electrochemical reactions in fuel cells [4–6, 11, 12].

B. Electrolyte Current

By definition, the electrolyte current is written as follows: 𝐼_𝑚=−𝜎_𝑚𝜕𝜙_𝑚𝜕𝑛.(B.1) Here 𝜙_𝑚 represents the potential in the ionic conducting media; 𝑛 is the normal to the interface; 𝜎_𝑚 is the membrane conductivity. Using linear approximation of the potential profile in the ionic conducting media, the electrolyte current is defined as follows: 𝐼_𝑚≈−𝜎_𝑚𝜙^𝐶𝑚−𝜙^𝐴𝑚𝛿_𝑚,(B.2) where 𝛿_𝑚 is the membrane thickness; 𝜙^𝐶𝑚 is the membrane phase potential at cathode electrode/membrane interface; 𝜙^𝐴𝑚 is the membrane phase potential at anode electrode/membrane interface. The values of potential in ionic conducting media at electrode/membrane interfaces can be expressed through potential difference 𝜂 given by the definition of the activation overpotential 𝜂𝜂^𝐴_act𝐶act=𝜙^𝐴−𝜙^𝐴𝑚−𝜂^𝐴_eq,=𝜙^𝐶−𝜙^𝐶𝑚−𝜂^𝐶_eq,(B.3) where 𝜙^𝐶 is the potential at cathode electrode/membrane interface in the electron conducting media; 𝜙^𝐴 is the potential at anode electrode/membrane interface in the electron conducting media.

We note that the reversible Nernst potential at the electrode/membrane interface is linked with anodic and cathodic equilibrium potential differences: 𝜂^𝐶_eq−𝜂^𝐴_eq=𝐸^OCV.(B.4) It is known that the difference of the anodic and cathodic electric potential is equal to the cell voltage 𝐸_Cell=𝜙^𝐶−𝜙^𝐴.(B.5) The values of potential in ionic conducting media at electrode/electrolyte interfaces can be expressed from (B.3) as follows: 𝜙𝜙^𝐴𝑚_𝐶𝑚=𝜙^𝐴−𝜂^𝐴_act−𝜂^𝐴_eq,=𝜙^𝐶−𝜂^𝐶_act−𝜂^𝐶_eq.(B.6) Taking into account equations (B.2)–(B.6), we obtain the next expression for the electrolyte current density: 𝐼_𝑚=𝐸^OCV−𝐸_cell−𝜂^𝐴_act+𝜂^𝐶_act𝑅^Ohmic,(B.7) where 𝑅^Ohmic is the ohmic resistance, 𝑅^Ohmic=𝛿_𝑚/𝜎_𝑚.

Application of new expression for electrolyte current with fuel cell models is shown in our papers [4–6, 11, 12].

C. Definition of Mixture Concentration

Mass balance equations for channels and catalyst layer are defined with molar component concentration. In view of gas-liquid flow in channels, the physical meaning of component concentration corresponds to mixture concentration. By definition, the mixture concentration is 𝐶^∗=𝐶^∗𝐺𝜑_𝐺+𝐶^∗𝐿1−𝜑_𝐺,(C.1) where 𝐶^∗ is mixture concentration, kg/m³; 𝐶^∗𝐺 is concentration in gas phase, kg/m³; 𝐶^∗𝐿 is concentration in liquid phase, kg/m³. Dividing both parts by molecular weight, we transfer to molar mixture concentration 𝐶^∗𝑀=𝐶^∗𝐺𝜑_𝐺𝑀+𝐶^∗𝐿1−𝜑_𝐺𝑀(C.2) or 𝐶^(𝑘)=𝑦^(𝑘)𝜑_𝐺𝜌_mol,𝐺+𝑥^(𝑘)1−𝜑_𝐺𝜌_mol,𝐿,(C.3) where 𝐶 is mixture concentration, mol/m³; 𝐶_𝐺 is concentration in gas phase, mol/m³; 𝐶_𝐿 is concentration in liquid phase, mol/m³. The condition of thermodynamic equilibrium in multicomponent gas-liquid mixture gives the relationship among species concentrations in phases 𝑥^(𝑘)=𝐶^(𝑘)𝐾^(𝑘)1−𝜑_𝐺𝜌_mol,𝐺+1−𝜑_𝐺𝜌_mol,𝐿.(C.4) Equation (C.4) is used for calculating the driving mass transfer force in mass transfer equation.

Abbreviations

Latin Letters

C:	Molar concentration/mol m−3
𝐶dL:	Double layer capacitance
Cp:	Specific heat/J mol−1 K−1
Cm:	Water content in membrane/mol m−3
𝐸cell:	Cell voltage/V
𝐸ocv:	Open circuit voltage/V
FA,FC:	Anode and cathode volumetric flow rate/m3 s−1
F:	Faraday’s constant/96485 C mol−1
K:	Distribution coefficient
kV:	Volumetric mass transfer coefficient/s−1
h:	Enthalpy/J mole−1
Δ𝐻𝑅:	Reaction enthalpy/J mol−1
I:	Current density/A m−2
𝐼0:	Exchange current density/A m−2
MEA:	Membrane electrode assembly
N:	Component molar flux/mol s−1
n:	Number of tanks
nd:	Net drug coefficient
ne:	Number of electrons
P:	Pressure/Pa
r:	Mass source term/mol s−1
R:	Ideal gas constant/J mol−1 K−1
𝑅ohmic:	Ohmic resistance/Ω m2
q:	Energy source term/J s−1
S:	Electrode area/m2
t:	Time/s
T:	Temperature/K
TSR:	Tank in series reactor
V:	Volume/m3
x:	mole fraction in liquid phase
y:	mole fraction in gas phase.

Greek Symbols

𝛼:	Heat transfer coefficient/W m−2 K−1
𝛼𝐴𝐴:	Anodic charge transfer coefficients for anode
𝛼𝐴𝐶:	Cathodic charge transfer coefficients for anode
𝛼𝐶𝐴:	Anodic charge transfer coefficients for cathode
𝛼𝐶𝐶:	Cathodic charge transfer coefficients for cathode
𝛽eff:	Effective mass transfer coefficient/m s−1
𝜈:	Stoichiometry coefficient
𝜂:	Potential difference/V
𝜌mol:	Molar density/mol m−3
𝜎:	Conductivity/Ω−1m−1
𝛿𝑚:	Membrane thickness/m
𝛿cat:	Catalyst thickness/m
𝜑𝐺:	Gas content/m3m−3.

Subscripts/Superscripts

in:	Inlet
out:	Outlet
eff:	Effective
i,j:	i,j tank
A:	Anode
C:	Cathode
cell:	Fuel cell
G:	Gas phase
lim:	Limiting
mol:	Molar
m:	Membrane
par:	Parasitic
ref:	Reference
S:	Interface; solid
sat:	Saturated
MeOH:	Methanol
O2:	Oxygen
H2O:	Water
N2:	Nitrogen
CO2:	Carbon dioxide.

Acknowledgment

The authors are grateful to IWT Vlaanderen, for financial support for OCPEC “Optimised Chemical Production with Electricity Cogeneration”.