*5.1.1 Electrical conductivity model*

merit to highlight separately single events (such as the ohmic losses in the electrolyte), thus allowing to easily study the effects of the variation of single quantities/parameters on the FC behavior. The main drawback of circuit models is that complex interactions and nonlinearity are not simulated in detail. More

system dynamics [16].

Jacobian matrix.

**34**

**5.1 Proton exchange membrane fuel cells**

**5. Distributed parameter models**

*Thermodynamics and Energy Engineering*

importantly, they can be implemented into circuit simulation software, to study the electrical interface of the FC with the power management electronics and system supervisor that is devoted to provide FC control, in order to study the overall

The knowledge of how physical fields (electric field, current density, flow, velocity, temperature, species concentrations) are distributed within internal components constitutes a pivotal aspect in FC analysis and design, since gradients and irregularities hamper the achievement of optimal performance, but they can not be gripped by zero-dimensional, lumped models. On the other hand, distributed models have to cope with additional challenges deriving from the huge number of grid points needed for a complete tessellation of the multilayer 3-D domain, resulting in issues of "curse of dimensionality" which can hardly be faced without resorting to supercomputers. Parallel computing with domain decomposition can overcome this challenge if less powerful computers are used, by assigning one subdomain to each processor and implementing the few interactions between subdomains. The two electric potentials, for the electronic and electrolyte phases, are coupled by the surface overpotentials at the catalyst layers, where reaction kinetics is modeled by the Butler-Volmer equation. In the case of a one-dimensional formulation, Newton's method is an efficient algorithm to integrate the equations, with LU factorization used at each interaction. Nevertheless, in the case of 2-D and 3-D formulations, the sparse Jacobian matrix produced by Newton's method is too large to be efficiently handled. In this case, Gauss elimination with a generalized minimal residual subroutine (GMRES) preconditioned with a Gauss-Seidel block and a multigrid algorithm has proven to be more suitable to face the non-symmetric

Typical multiphysics coupled models include, among others, proton conduction,

water and fuel transport, joule dissipation, and thermal diffusion. The models, typically discretized with the finite element method (FEM), pose significant numerical challenges. Some commercial simulation tools like COMSOL®

Multiphysics allow the solution of general time-dependent systems of partial differential equations [25] and are therefore very useful tools for this class of problems. For the computation of the fluid-dynamic field, particularly in the case of turbulent motion at high Reynolds numbers, the finite volume method is also used. Ansys Fluent® is a commercial package based on this method particularly efficient in modeling fuel cells. PEM fuel cells, as the name implies, are based on protonconducting polymeric membranes. The most commonly used material for their realization is persulfonated polytetrafluoroethylene, commercialized as Nafion® by

Chemours. This material has a structure similar to the one of PTFE, but is

functionalized with sulfonic acid groups providing charge sites for proton conduction [7]. If the membrane is properly hydrated, protons can form hydronium

As briefly mentioned above, protonic conduction in Nafion® strongly depends on the temperature, since the mechanism is based on charged particles jumping from site to site, with a rate described by the diffusivity *D*. This statistical parameter depends on the activation barrier energy which exhibits an exponential dependence on the temperature *T* according to the law:

$$D = D\_o \ e^{-W\_{ai}/k \cdot T} \tag{21}$$

where *Do* is a diffusivity reference value, *Wai* is activation barrier energy, and *k* is Boltzmann's constant. The charged particle mobility *μ* is proportional to *D* according to the Einstein relation:

$$
\mu = \frac{|z|FD}{RT} \tag{22}
$$

with |*z*| the ion charge number, so that the proton conductivity *σ* = *ρ<sup>c</sup> μ* can be written as:

$$
\sigma = \frac{(zF)^2 c D\_o}{RT} \quad e^{-W\_{ii}/kT} = \sigma\_o \ e^{-W\_{ii}/kT} \tag{23}
$$

being *ρ<sup>c</sup>* = |*z*|, *Fc* the charge density, and *c* the molar concentration.

Apart from the temperature, the proton conductivity also depends on the water content in the membrane. A common modeling approach to represent such dependence is based on the hydration *λ*, that is, the ratio between the number of water molecules and the number of charge sites available for proton conduction. In the specific case of Nafion®, such ratio can be rewritten in a modified form using the water concentration *cw* and the sulfonic acid concentration *cas*, that is, *λ* = *cw*/*cas*. In absence of more sophisticated models, a linear dependence of conductivity on hydration can be assumed:

$$
\sigma = a \; \lambda = 0.5139 \; \lambda \quad \text{S/m} \tag{24}
$$

where *λ* is derived from a correlation empirically derived for Nafion® [7]. Combining the above, the factorized expression of *σ*(*λ*,*T*) is obtained:

$$
\sigma(\lambda, T) = a\lambda \text{ } e^{\frac{W\_{\text{ai}}}{k} \left(\frac{1}{30} + \frac{1}{7}\right)} \tag{25}
$$

with *Wai*/*k* = 1268 K for ions hopping. This is a more sophisticated model than the one expressed by Eq. (9). The conductivity *σ* influences the scalar potential *φ* according to the charge conservation equation in quasi-static conditions:

$$\nabla \cdot \sigma(\lambda, T) \nabla \varphi = 0 \tag{26}$$

which shows that the distribution of *φ* depends on *λ* and *T*.

#### *5.1.2 Hydration model*

A critical issue in modeling PEMFCs consists in providing an accurate description of the hydration effects, which rules proton conductivity [26]. The distribution of *λ* in the membrane can be computed resorting to specific equations at the surfaces and in the bulk. For the membrane bulk, two mechanisms are taken into account, namely electro-osmotic drag and back-diffusion, giving rise to the following expression of the water molar flow:

$$\mathbf{N}\_w = \mathbf{N}\_{w\epsilon} + \mathbf{N}\_{wd} = \xi \lambda \frac{J}{|z|F} - D\_w c\_{as} \nabla \lambda \tag{27}$$

*<sup>ρ</sup>cp <sup>∂</sup><sup>t</sup> <sup>T</sup>* � <sup>∇</sup> � *<sup>k</sup>*ð Þ*<sup>λ</sup>* <sup>∇</sup>*<sup>T</sup>* � *σ λ*ð Þ *, T* j j <sup>∇</sup>*<sup>φ</sup>*

where *ρ* is the hydrated Nafion® density, *cp* its specific heat, and *k*(*λ*) the thermal conductivity. According to data given in [28], *k* = 0.12 + 0.81*λ* (W m�<sup>1</sup> K�<sup>1</sup>

The complete model to be solved is assembled from the above equations together with proper boundary (time-dependent Dirichlet and homogeneous Neumann) and initial conditions. The specific characteristics of the set of partial differential equations together with the high aspect ratio of the geometry (i.e., the very small thickness of the membrane compared to its extension in the plane) lead, after discretization with FEM, to a badly conditioned system of linear equations. The system tends to be quite large so that direct solvers may become inapplicable and it is also difficult to precondition so that iterative solvers tend to converge slowly or to fail altogether. Due to the strong nonlinearity of the complete coupled problem, a further numerical challenge concerns the nonlinear solver, typically the Newton-Raphson (NR) algorithm, which usually converges only, if at all, with strong under-

The numerical solution of the final system may require substantial effort in correctly setting linear and nonlinear solver parameters to achieve convergence. Extreme care is needed in the choice of the drop tolerance if GMRES coupled to an ILU pre-conditioner with threshold is applied to solve the linear system arising in Newton's method. Numerical experiments have shown that the selection of the drop tolerance needed to obtain convergence is highly problem-dependent so that the only way to reliably obtain a solution was to apply an efficient direct solver (PARDISO®). This however limited the maximum admissible problem size on the available hardware. Furthermore, the fully coupled nonlinear system of equations can typically not be solved in its complete assembled form by the standard NR technique, even with strong under-relaxation. The problem was therefore solved with a further iterative loop, that is, a so-called "segregated solution," in which one or more blocks of equations are fed into the next one in a simple iteration until convergence. In our specific case, the electrical model formed the first block, and its solution was inserted into the thermal model; then, the solution of both was used in the diffusion model. An important technological problem in the construction of PEM fuel cells is the reproducibility of the production process in the case of large membranes. In particular, variations in the thickness can result in hot spots which can cause aging effects and may lead to membrane failure. In order to investigate such effects, a lens-shaped compression of 2 mm radius and 50 μm depth (1/4 of the total PEM thickness of 200 μm) was studied. Results show an increase in current density in the order of 100%, leading to strong localized overheating (**Figure 5**).

Spatially resolved analyses of DMFCs allow studying their limitations, like sluggish kinetics and methanol crossover. These models resort to a detailed description of the real device geometry and materials, and of local nonlinear physics phenomena such as heat, momentum, multicomponent mass transport, and electrochemical processes. Early distributed parameter models were multi-domain ones, in the sense that the problem variables were defined in separated domains by introducing appropriate internal boundary conditions. As an example, in [12], the equations

can be assumed. The last term at the left-hand side represents Joule's losses.

*5.1.4 Coupled multiphysics model*

*Distributed and Lumped Parameter Models for Fuel Cells*

*DOI: http://dx.doi.org/10.5772/intechopen.89048*

**5.2 Direct methanol fuel cells**

**37**

relaxation.

<sup>2</sup> <sup>¼</sup> <sup>0</sup> (32)

)

where *N<sup>i</sup>* is the ionic molar flow vector and *J* the current density vector, and *Dw* is the water diffusivity in the membrane. This equation is nonlinear because *Dw* itself depends on *λ* and also on *T.* Such dependences can be expressed by factorizing the statistical mechanics exponential dependence on *T* with a polynomial regression obtained from experimental data:

$$D\_w(\lambda, T) = \sum\_{i=0}^{\overline{3}} d\_i \lambda^i e^{\frac{W\_{av}}{k} \left(\frac{1}{305} - \frac{1}{T}\right)}\tag{28}$$

where *Waw*/*<sup>k</sup>* = 2416 K for water in Nafion® and *<sup>d</sup>*<sup>0</sup> = 2.563671 � <sup>10</sup>�<sup>6</sup> *, <sup>d</sup>*<sup>1</sup> <sup>=</sup> �0.33671 � <sup>10</sup>�<sup>6</sup> *, d*<sup>2</sup> = 0.0264 � <sup>10</sup>�<sup>6</sup> *,* and *<sup>d</sup>*<sup>3</sup> = 0.000671 � <sup>10</sup>�<sup>6</sup> for *Dw* (in cm2 s �1 ). The dynamics of *N<sup>w</sup>* is ruled by Fick's second law, which can be written in terms of *λ*:

$$\nabla \cdot \mathbf{N}\_w + \mathfrak{c}\_{ds} \partial\_t \lambda = \mathbf{0} \tag{29}$$

Letting Eq. (27) in Eq. (29) and assuming |*z*| = 1 for protons provide the following diffusion equation:

$$
\nabla \cdot D\_w \nabla \lambda - \partial\_t \lambda - \nabla \cdot \frac{\xi \lambda \mathbf{J}}{c\_{as} F} = \mathbf{0} \tag{30}
$$

According to Maxwell's equations, the current density *J* can be expressed in terms of *φ*, so that *Eq.* (30) becomes:

$$\nabla \cdot D\_w \left( \lambda, T \right) \nabla \lambda - \partial\_t \lambda - \nabla \cdot \frac{\xi \lambda}{c\_{ds} F} \sigma(\lambda, T) \nabla \varphi = \mathbf{0} \tag{31}$$

Imposing interfaces conditions at the electro-catalyst layers is troublesome since the precise distribution of *λ* at such surfaces cannot be determined with sufficient accuracy. However, hydration can be expressed in terms of a more easily representable quantity, that is, the water vapor activity (that is, the relative humidity). The relationship between these two physical quantities can be modeled by resorting to an empirical relationship [27], which in turn can be expressed as a function of the water vapor pressure by another empirical model.

#### *5.1.3 Thermal model*

Since most of quantities appearing in the conductivity and hydration models depend on temperature, it is also necessary to take into account the transient heat equation:

*Distributed and Lumped Parameter Models for Fuel Cells DOI: http://dx.doi.org/10.5772/intechopen.89048*

$$\rho c\_p \partial\_t T - \nabla \cdot k(\lambda) \nabla T - \sigma(\lambda, T) |\nabla \rho|^2 = \mathbf{0} \tag{32}$$

where *ρ* is the hydrated Nafion® density, *cp* its specific heat, and *k*(*λ*) the thermal conductivity. According to data given in [28], *k* = 0.12 + 0.81*λ* (W m�<sup>1</sup> K�<sup>1</sup> ) can be assumed. The last term at the left-hand side represents Joule's losses.

## *5.1.4 Coupled multiphysics model*

*5.1.2 Hydration model*

ing expression of the water molar flow:

*Thermodynamics and Energy Engineering*

obtained from experimental data:

*<sup>d</sup>*<sup>1</sup> <sup>=</sup> �0.33671 � <sup>10</sup>�<sup>6</sup>

following diffusion equation:

*5.1.3 Thermal model*

equation:

**36**

terms of *φ*, so that *Eq.* (30) becomes:

water vapor pressure by another empirical model.

cm2 s �1

terms of *λ*:

A critical issue in modeling PEMFCs consists in providing an accurate description of the hydration effects, which rules proton conductivity [26]. The distribution of *λ* in the membrane can be computed resorting to specific equations at the surfaces and in the bulk. For the membrane bulk, two mechanisms are taken into account, namely electro-osmotic drag and back-diffusion, giving rise to the follow-

where *N<sup>i</sup>* is the ionic molar flow vector and *J* the current density vector, and *Dw* is the water diffusivity in the membrane. This equation is nonlinear because *Dw* itself depends on *λ* and also on *T.* Such dependences can be expressed by factorizing the statistical mechanics exponential dependence on *T* with a polynomial regression

3

*i*¼0 *di λ<sup>i</sup> e Waw k* 1

). The dynamics of *N<sup>w</sup>* is ruled by Fick's second law, which can be written in

*ξ λ J*

where *Waw*/*<sup>k</sup>* = 2416 K for water in Nafion® and *<sup>d</sup>*<sup>0</sup> = 2.563671 � <sup>10</sup>�<sup>6</sup>

Letting Eq. (27) in Eq. (29) and assuming |*z*| = 1 for protons provide the

According to Maxwell's equations, the current density *J* can be expressed in

Imposing interfaces conditions at the electro-catalyst layers is troublesome since the precise distribution of *λ* at such surfaces cannot be determined with sufficient accuracy. However, hydration can be expressed in terms of a more easily representable quantity, that is, the water vapor activity (that is, the relative humidity). The relationship between these two physical quantities can be modeled by resorting to an empirical relationship [27], which in turn can be expressed as a function of the

Since most of quantities appearing in the conductivity and hydration models depend on temperature, it is also necessary to take into account the transient heat

<sup>∇</sup> � *Dw*∇*<sup>λ</sup>* � *<sup>∂</sup><sup>t</sup> <sup>λ</sup>* � <sup>∇</sup> �

<sup>∇</sup> � *Dw* ð Þ *<sup>λ</sup>, T* <sup>∇</sup>*<sup>λ</sup>* � *<sup>∂</sup><sup>t</sup> <sup>λ</sup>* � <sup>∇</sup> � *ξ λ*

j j *<sup>z</sup> <sup>F</sup>* � *Dw cas*∇*<sup>λ</sup>* (27)

<sup>303</sup>�<sup>1</sup> ð Þ*<sup>T</sup>* (28)

*cas <sup>F</sup>* <sup>¼</sup> <sup>0</sup> (30)

*cas <sup>F</sup> σ λ*ð Þ *, T* <sup>∇</sup>*<sup>φ</sup>* <sup>¼</sup> <sup>0</sup> (31)

*,* and *<sup>d</sup>*<sup>3</sup> = 0.000671 � <sup>10</sup>�<sup>6</sup> for *Dw* (in

<sup>∇</sup> � *<sup>N</sup><sup>w</sup>* <sup>þ</sup> *cas <sup>∂</sup><sup>t</sup> <sup>λ</sup>* <sup>¼</sup> <sup>0</sup> (29)

*,*

*<sup>N</sup><sup>w</sup>* <sup>¼</sup> *<sup>N</sup>we* <sup>þ</sup> *<sup>N</sup>wd* <sup>¼</sup> *ξλ <sup>J</sup>*

*Dw*ð Þ¼ *<sup>λ</sup>, T* <sup>X</sup>

*, d*<sup>2</sup> = 0.0264 � <sup>10</sup>�<sup>6</sup>

The complete model to be solved is assembled from the above equations together with proper boundary (time-dependent Dirichlet and homogeneous Neumann) and initial conditions. The specific characteristics of the set of partial differential equations together with the high aspect ratio of the geometry (i.e., the very small thickness of the membrane compared to its extension in the plane) lead, after discretization with FEM, to a badly conditioned system of linear equations. The system tends to be quite large so that direct solvers may become inapplicable and it is also difficult to precondition so that iterative solvers tend to converge slowly or to fail altogether. Due to the strong nonlinearity of the complete coupled problem, a further numerical challenge concerns the nonlinear solver, typically the Newton-Raphson (NR) algorithm, which usually converges only, if at all, with strong underrelaxation.

The numerical solution of the final system may require substantial effort in correctly setting linear and nonlinear solver parameters to achieve convergence. Extreme care is needed in the choice of the drop tolerance if GMRES coupled to an ILU pre-conditioner with threshold is applied to solve the linear system arising in Newton's method. Numerical experiments have shown that the selection of the drop tolerance needed to obtain convergence is highly problem-dependent so that the only way to reliably obtain a solution was to apply an efficient direct solver (PARDISO®). This however limited the maximum admissible problem size on the available hardware. Furthermore, the fully coupled nonlinear system of equations can typically not be solved in its complete assembled form by the standard NR technique, even with strong under-relaxation. The problem was therefore solved with a further iterative loop, that is, a so-called "segregated solution," in which one or more blocks of equations are fed into the next one in a simple iteration until convergence. In our specific case, the electrical model formed the first block, and its solution was inserted into the thermal model; then, the solution of both was used in the diffusion model. An important technological problem in the construction of PEM fuel cells is the reproducibility of the production process in the case of large membranes. In particular, variations in the thickness can result in hot spots which can cause aging effects and may lead to membrane failure. In order to investigate such effects, a lens-shaped compression of 2 mm radius and 50 μm depth (1/4 of the total PEM thickness of 200 μm) was studied. Results show an increase in current density in the order of 100%, leading to strong localized overheating (**Figure 5**).

### **5.2 Direct methanol fuel cells**

Spatially resolved analyses of DMFCs allow studying their limitations, like sluggish kinetics and methanol crossover. These models resort to a detailed description of the real device geometry and materials, and of local nonlinear physics phenomena such as heat, momentum, multicomponent mass transport, and electrochemical processes. Early distributed parameter models were multi-domain ones, in the sense that the problem variables were defined in separated domains by introducing appropriate internal boundary conditions. As an example, in [12], the equations

accounted for in [31]. The two-phase mass transport in the anode and cathode porous regions is formulated based on the classical multiphase flow in porous media. A micro-agglomerate model, that is able to reflect the effect of the microstructure of the catalyst layer on cell performance, is proposed. The resulting polarization curves and methanol crossover rates at different concentrations are in very good agreement with experimental data. In [32], a realistic passive liquid-feed DMFC in transient charge/discharge conditions is simulated. The main contribution of this work is that effects of feed methanol concentration in the reservoir and current density on both mass transport and performance are investigated. Analyses show that when the initial feed concentration in the reservoir decreases, methanol crossover is minimized, but the fuel cell runtime is shortened. Recent works, for example [33], provide a detailed description of the whole DMFC system, including reservoirs. By developing a transient multiphase model of a passive cell, the effects of operating current density, voltage, micro-porous layer, and methanol feeding condition are comprehensively investigated for the whole operating processes, that is, with the fuel tank evolving from full to empty. Results highlight that for all operating conditions, it is necessary to operate at moderate current density or voltage to limit the methanol crossover and ensure the energy conversion efficiency. A 2-D multiphase non-isothermal mass transfer model for the DMFC is presented in [34]. The model includes the reaction of methanol and oxygen at the anode and cathode and the diffusion of every component involved in the DMFC, such as water, oxygen, and methanol at the diffusion layer and methanol crossover. It is shown that a maximum output power can be achieved for optimal temperature

*Distributed and Lumped Parameter Models for Fuel Cells*

*DOI: http://dx.doi.org/10.5772/intechopen.89048*

A few papers are reported on 3-D two-phase DMFC models, which can capture the species distributions and the transport limitations along any direction inside the DMFC. Ref. [35] proposes a 3-D, two-phase, multicomponent model. Catalyst layers are incorporated in the computational domain instead of being modeled as zero-thickness interfaces. This model includes the effects of the second phase on the reduction of active catalyst surface areas and the mixed potential effects due to methanol crossover. The amount of carbon dioxide obtained from 3-D models indicates that the porosity of the anode diffusion layer plays a very important role in the DMFC performance. With a low porosity, the produced carbon dioxide cannot be removed effectively from the catalyst layer, thus reducing the active catalyst surface area as well as blocking methanol from reaching the reaction area. Ref. [31] extends the 2-D two-phase mass transport model for liquid-feed DMFCs to a fully 3-D model. The two-phase mass transport in the anode and cathode porous regions is formulated based on the multiphase flow theory in porous media without defining the mixture pressure of gas and liquid and assuming a constant gas pressure in the porous regions. The interaction between the phases in this 3-D model is captured by taking into account the effect of non-equilibrium evaporation/ condensation at the phase interface, as opposed to the assumption of other models of thermodynamic equilibrium condition between the phases. From 3-D analysis, it can be observed that methanol concentration in the diffusion layer is higher in the channel than under the ribs, demonstrating that the flow-field channels cause

methanol to be distributed unevenly over the reaction area (**Figure 6**).

In [36], a commercial flow solver (i.e., Ansys Fluent®) is used to solve at the same time flow, species, and charge transport equations. 3-D simulations are carried

and concentration values.

**39**

*5.2.2 DMFC three-dimensional models*

**Figure 5.** *Local increase of current density due to a manufacturing defect (courtesy of IEEE Trans. on Magnetics).*

describing the concentration and potential distribution within the electrode were solved numerically using the finite difference method (FDM) and Newman's BAND algorithm for the resulting simultaneous nonlinear equations. After the introduction of computational fluid dynamics (CFD) in the simulation of fuel cells, mostly single-domain models have been developed. The main advantages over the multidomain approach is that internal BCs and continuity conditions at each domain interface are not required, thus simplifying the model geometry construction and speeding up problem set-up into a commercial CFD code. Single-domain CFD models can be classified into two-dimensional (2-D) and three-dimensional (3-D) models, depending on the simplifying assumptions on geometry and on boundary conditions. 2-D models generally provide a strong reduction in terms of computational cost, but their solution is less accurate compared to 3-D models.
