Stability and convergence of the numerical solution for the species equation of a model for PEMFCs

Computers and Mathematics with Applications xxx (xxxx) xxx

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage: www.elsevier.com/locate/camwa

Stability and convergence of the numerical solution for the species equation of a model for PEMFCs M.M. De Souza, A.L. De Bortoli

∗

Federal University of Rio Grande do Sul, Av. Bento Gonçalves 9500, PO Box 15080, Porto Alegre, RS 91509-900, Brazil

article

info

Article history: Received 8 June 2018 Received in revised form 5 October 2018 Accepted 2 November 2018 Available online xxxx Keywords: Numerical analysis Proton exchange membrane fuel cell Finite element method Mole fraction

a b s t r a c t Many mathematical models for fuel cells are presented in the literature, however, few studies of numerical analysis of these models are found. The novelty of this work is the proof of stability and convergence of the equation of chemical species, providing greater reliability to the mathematical model. In this work, a mathematical model is proposed that calculates the flow and concentration of species for proton exchange membrane fuel cells using the finite element method for spatial discretization and the Crank–Nicolson method for temporal discretization. The model takes into account the losses overpotentials at the anode and at the cathode. The numerical results confirm the proof of convergence which is demonstrated in this work. Finally, a numerical result is shown for the fuel cell power relative to current density, and the result has good relation with the experimental data. © 2018 Elsevier Ltd. All rights reserved.

1. Introduction Oil, coal and natural gas are the most widely used fuels in the world, around 80% of the world’s energy matrix. As they are not renewable, the burning of these fuels contributes directly to global warming, acid rain, among other problems. In addition, non-renewable energy sources are at risk of scarcity, as energy consumption has been increasing over time. Therefore, sustainable development through clean and efficient energy sources is increasingly desired. In this context, a device that has attracted great interest precisely for its efficiency and low emission of pollutants is the fuel cell. Fuel cells are generating increasing interest, because they have low emission of pollutants and transform the chemical energy of the fuel directly into electricity [1]. Fuel cells provide clean energy and high efficiency in a wide variety of applications. They are candidates to replace current batteries and have the potential to be sustainable energy sources in the future [2]. Among fuel cell types, the proton exchange membrane fuel cell (PEMFC) has been preferred, because it operates at low temperatures (up to 373 K). Mathematical models involving fuel cells have been developed to better understand the physical and chemical behavior within the cell, and the conversion of chemical energy into electricity. Generally, the performance of a fuel cell is given by means of the polarization curve, which registers the electrical potential difference between the two electrodes according to the electric current produced. Thus, some mathematical models have been developed. Siegel [3] made a review of computational heat and mass transfer in PEMFCs, whose articles began to be published in the early 1990s. Among some recent works, Wang et al. [4] developed a model of hybrid PEMFC for power management control. Goel and Basu [5] developed a one-dimensional steady-state model to analyze the performance of a direct ethanol fuel cell. Iranzo et al. [6] analyzed some mathematic models for proton ∗ Corresponding author. E-mail address: [email protected] (A.L. De Bortoli). https://doi.org/10.1016/j.camwa.2018.11.010 0898-1221/© 2018 Elsevier Ltd. All rights reserved.

Please cite this article as: M.M. De Souza and A.L. De Bortoli, Stability and convergence of the numerical solution for the species equation of a model for PEMFCs, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.11.010.

2

M.M. De Souza and A.L. De Bortoli / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 1. Sketch of a PEMFC.

exchange membrane fuel cells fed by hydrogen. Gong and Cai [7] develop a method to optimize PEMFC parameters. Carral and Mélé [8] developed a model to investigate the influence of the PEMFC stacking phase on the mechanical state of the active layer. Hu et al. [9] described a three-dimensional model to investigate fluid flow, species transport, and electrochemical reactions in a PEMFC and Wu et al. [10] presented a numerical analysis of dynamics of processes in fully humidified PEMFCs. In this work, a mathematical model is used to obtain the velocity field and the mole fractions of the main species involved in the oxidation of ethanol and the reduction of oxygen within a PEMFC. Numerical results are obtained using the finite element method, which is known to have a good theoretical mathematical basis [11,12]. The equations are considered in the transient form, and are discretized by the Crank–Nicolson method [13]. In addition, the model considers the overpotential losses, due to activation, ohmic resistance and concentration. These losses are calculated using the species concentration on the catalyst surface, in according to the current density. Mathematical models capable of simulating the operations of PEMFCs and investigating the fluid flow and species concentrations in these devices were developed, but there are no studies that show the stability and convergence of these models. Therefore, we propose to make a numerical analysis of the equations that model the molar fraction of the species involved in the chemical reactions for the oxidation and reduction processes that occur in fuel cells. The stability and convergence analysis of equations that model a PEMFC using the finite element method are unpublished. So, we believe that this will increase the reliability of these models. We present two results: the first confirms the convergence rate obtained in the demonstration of the convergence theorem and the second presents the behavior of the molar fraction of the species in a direct ethanol fuel cell, which is used to calculate the overpotential losses and finally compare the voltage and power density of the model with experimental data. 2. Proton exchange membrane fuel cell In a typical proton exchange membrane fuel cell, the fuel is inserted on the anode side, and reacts to form carbon dioxide, protons, and electrons. The anode is located on the left side and the cathode is on the right side, as shown in Fig. 1. δdl is the thickness of the diffusion layer, δcl the thickness of the catalyst layer, δfc the thickness of the flow channel and δm the thickness of the membrane. The protons pass to the cathode through the membrane and the electrons through an external circuit, providing a potential difference. On the cathode side, air reacts with the protons and electrons formed at the anode to produce water vapor [14]. For example, using a mixture of ethanol and water as fuel, the electrochemical reactions occurring inside the cell are given by [15]: C2 H5 OH + 3H2 O → 2CO2 + 12H+ + 12e− 3O2 + 12H+ + 12e− → 6H2 O

at the anode,

at the cathode,

(1) (2)

corresponding to the overall reaction C2 H5 OH + 3O2 → 2CO2 + 3H2 O.

(3)

The catalysts used in the PEMFCs are aimed at accelerating reactions and are generally made of platinum, because of the ease of separation of the hydrogen molecules. But, depending on the fuel used (methanol, ethanol), researchers opt to use bi-metallic catalysts, as the PtRu/C catalyst, because of the difficulty in break C–C bonds, to form intermediates that are absorbed at the surface of the catalyst. Please cite this article as: M.M. De Souza and A.L. De Bortoli, Stability and convergence of the numerical solution for the species equation of a model for PEMFCs, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.11.010.

M.M. De Souza and A.L. De Bortoli / Computers and Mathematics with Applications xxx (xxxx) xxx

3

The membrane most commonly used in PEMFCs is DuPont’s Nafion, a perflurosulfonic acid polymer membrane consisting of a Teflon-type backbone with side chains terminating with -SO3 H groups [16]. Nafion membrane is mechanically strong, acidic and chemically resistant, very absorbent to water and good proton-conducting (H+ ), if well hydrated. 3. Notation and preliminaries The L2 (Ω ) norm and inner product will be denoted by ∥·∥ and (·, ·), respectively. The Lp (Ω ) norms and the Sobolev Wpk (Ω ) norms are denoted by ∥ · ∥Lp and ∥ · ∥W k , respectively. For the semi-norm in Wpk (Ω ) we use |·|W k . The space of Sobolev W2k p

p

is represented by H k and ∥ · ∥k denote the norm in H k . { } Consider the space V = H01 (Ω ), Q = L20 (Ω ) and Z = H01 (Ω ), where H01 (Ω ) = v ∈ H 1 : v = 0 in ∂ Ω and L20 = {q ∈ ∫ L2 (Ω ) : Ω q dxdy = 0}. The dual space of V is V ∗ , with the norm ∥ · ∥∗ . Choose finite dimension spaces Vh ⊂ V , Qh ⊂ Q and Zh ⊂ Z . Let Vh and Qh be given by Taylor–Hood element, i.e., quadratic polynomials for velocity and linear polynomials for pressure [17]. For the equations of the species, we choose quadratic polynomials. We make use of the following properties [11]: inf ∥u − v∥ ≤ Chk+1 ∥u∥k+1 ,

v∈V h

inf ∥u − v∥1 ≤ Chk ∥u∥k+1 ,

v∈V h

inf ∥p − q∥ ≤ Chs+1 ∥p∥s+1 ,

q∈Q h

inf ∥X − w∥ ≤ Chr +1 ∥X ∥r +1 ,

w∈Z h

inf ∥X − w∥1 ≤ Chr ∥X ∥r +1 ,

(4)

w∈Z h

where u ∈ H k+1 (Ω ), p ∈ H s+1 (Ω ) and X ∈ H r +1 (Ω ). Discrete time error analysis requires norms that are analogues to the norms used in the case of continuous time [18,19]

|||v|||∞,k := max ∥vn ∥k ,

(5)

0≤n≤Nt

( |||v|||p,k :=

N ∑

)1/p p n k

∆t ∥v ∥ ,

.

(6)

n=0

Some results and definitions will be presented below, as they will be fundamental for the demonstrations of the next sections [17–20]. Definition 3.1. Define the skew-symmetric trilinear form b : V × V × V → R as b(u, v, w ) :=

1 2

(u · ∇ V , w ) −

1 2

(u · ∇w, v ).

Lemma 3.1. For u, v, w ∈ V the trilinear term b(u, v, w ) can be bounded in the following ways: b(u, v, w ) ≤ C0 (Ω )∥∇ u∥∥∇v∥∥∇w∥, b(u, v, w ) ≤ C0 (Ω )∥u∥

1/2

∥∇ u∥

1/2

∥∇v∥∥∇w∥.

(7) (8)

Lemma 3.2. Let ∆t be a fixed positive real number. If u is smooth enough, then 2

∥un+1/2 − u(tn+1/2 )∥ ≤

∆t 3 48

tn + 1

∫

∥utt ∥2 dt ,

(9)

tn

 2 ∫ tn+1 3  un+1 − un    ≤ ∆t − u (t ) ∥uttt ∥2 dt , t n + 1 / 2  ∆t  1280 tn ∫ ∆t 3 tn+1 2 ∥∇ (un+1/2 − u(tn+1/2 ))∥ ≤ ∥∇ utt ∥2 dt . 48

(10) (11)

tn

Definition 3.2. For f ∈ L2 (Ω ), the V ∗ norm of f is

∥f ∥∗ := sup v∈V

(f , v )

∥∇v∥

.

Please cite this article as: M.M. De Souza and A.L. De Bortoli, Stability and convergence of the numerical solution for the species equation of a model for PEMFCs, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.11.010.

4

M.M. De Souza and A.L. De Bortoli / Computers and Mathematics with Applications xxx (xxxx) xxx Table 1 Source term f .

f

Flow channel

Diffusion layer

0

−εdl

µ u κ

Catalyst layer

−εcl

µ u κ

Thus, the following inequality is valid, (f , v ) ≤ ∥f ∥∗ ∥∇v∥,

(12)

for all v ∈ V . 4. Governing equations In the mathematical model used in this work we consider that the membrane remains well hydrated, the flow is incompressible, laminar and isothermal. On the side of the anode and the cathode the equations for the flow, continuity and momentum, and the concentration of species have the same form. The equations of momentum, continuity and concentration of the species used in this article were adapted from the works of Liu and Wang [21], Le et al. [22] and Wang and Wang [23]. In addition, constant viscosity is considered. In this way, the equations are given by [21,22]

∂u 1 − ν∇ 2 u + u · ∇ u + ∇ p = f in (0, T ] × Ω , ∂t ρ ∇ · u = 0 in [0, T ] × Ω , u(0; x, y) = u0 in Ω ,

(13) (14) (15)

where ρ is the density, ν the viscosity and Ω the domain. Table 1 shows the source term f , where εdl and εcl are the porosities of the diffusion layer and of the catalyst layer, respectively, and κ is the permeability. The velocity of the fluid in the diffusion and catalyst layers is described by Darcy’s law, which is an equation describing the flow of a fluid through a porous medium [21]. The use of Darcy’s law is valid because the flow is laminar and viscous with very low velocity, the local Reynolds number is less than 1 in the diffusion layer and in the catalyst layer, since the flow rate at the anode is 1 ml/min and at the cathode is 120 ml/min. The general conservation equation of a species can be written as [23]

ρ

∂ Xi 2 + ρ u · ∇ Xi − ρ Deff i ∇ Xi = Si , with ∂t

∑

Xi = 1,

(16)

i

Xi (0; x, y) = Xi,0 in Ω ,

(17)

Deff i

where Xi is the mole fraction of the species i, is the effective diffusion coefficient and Si is the source term. Si is zero in the channels and in the diffusion layer, and in the catalyst layer is given by

γi Mi

j, n F where j is the volumetric current density in the fuel cell, Mi is the molecular weights of species i, γi stoichiometric coefficient of i, n total number of electrons produced in reaction and F is the Faraday’s constant. The effective diffusivity is given by [21] Si =

3/2 Deff , i = Di ε

where ε is the media porosity and Di is the diffusion of the ith component. The variational formulation of (13)–(14) for velocity u : [0, T ] → V and pressure p : [0, T ] → Q is given by

(

) ∂u 1 , v + ν (∇ u, ∇v ) + b(u, u, v ) − (p, ∇ · v ) = (f , v ), ∂t ρ (∇ · u, q) = 0,

for all v ∈ V and q ∈ Q . In the same way, assuming w ∈ Z , we have the variational formulation for species Xi : [0, T ] → Z given by

ρ

(

) ∂ Xi , w + ρ Deff i (∇ Xi , ∇w ) + ρ b(u, Xi , w ) = (Si , w ). ∂t

The domain discretization is obtained by triangulation, subdividing Ω into a set of non-overlapping triangles. We consider the finite dimensional spaces Vh ⊂ V , Qh ⊂ Q and Zh ⊂ Z for velocity, pressure and species, respectively. Please cite this article as: M.M. De Souza and A.L. De Bortoli, Stability and convergence of the numerical solution for the species equation of a model for PEMFCs, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.11.010.

M.M. De Souza and A.L. De Bortoli / Computers and Mathematics with Applications xxx (xxxx) xxx

5

Thus, the solution by finite elements of momentum and continuity equations are given by, respectively

(

) 1 ∂ uh h , v + ν (∇ uh , ∇v h ) + b(uh , uh , v h ) − (ph , ∇ · v h ) = (f , v h ), ∂t ρ (∇ · uh , qh ) = 0,

(18) (19)

where v ∈ Vh and q ∈ Qh . The finite element solution of species equations is given by h

ρ

(

h

) ∂ Xih h h h h h h h , w + ρ Deff i (∇ Xi , ∇w ) + ρ b(u , Xi , w ) = (Si , w ), ∂t

(20)

∀wh ∈ Zh . The equations considered in this work are in the transient form, since the time discretization method is essential in the process of obtaining the approximate solution. The Navier–Stokes equations have characteristics that influence the choice of the methods used for time discretization, such as the stiffness of the problem in diffusion-dominated flow regions, the general sensitivity of the physical problem, and the lack of regularity of its solution [17]. Therefore, numerical simulations were performed using the Crank–Nicolson method (CN), which is semi-implicit, second order accurate and A-stable [13]. The momentum and continuity equations in the discretized form are given by, respectively [24,25]

(

uhn+1 − uhn

∆t

) ,v

+ ν (∇ uhn+ 1 , ∇v h ) + b(uhn+ 1 , uhn+ 1 , v h )

h

2

2

2

1

− (phn+1 , ∇ · v h ) = (fn+ 1 , v h ), 2 ρ (∇ · uhn+1 , qh ) = 0, λ

(21) (22)

+λ

where v h ∈ Vh , qh ∈ Qh and λn+ 1 = n+12 n . 2 The Crank–Nicolson method leads at each discrete time step to a nonlinear system of equations, which is solved by fixedpoint iteration. The set of species equations in the discretized form is given by

( ρ

Xih,n+1 − Xih,n

∆t

+ρ

b(uhn+1

,

) ,w

h

Xh 1 i,n+ 2

h h + ρ Deff i (∇ Xi,n+ 1 , ∇w ) 2

, w ) = (Si,n+ 1 , wh ), h

(23)

2

where w h ∈ Zh . To facilitate notating, we will use Di instead of Deff i , X in place of Xi and S rather than Si . 5. Numerical analysis In this section, we will show that the solution of the discrete form the species Eq. (23) is unconditionally stable and convergent for the species Eq. (16). First, we show the stability of the discrete Navier–Stokes equations (21)–(22), whose results will be used later. Lemma 5.1 (Stability). Consider the approach of Crank–Nicolson FEM for (21)–(22). If a solution uhi , i = 1, . . . , N, exists at every time-step and the scheme is unconditionally stable, then it satisfies the following a priori condition: 2

∥uhN ∥ + ∆t ν

N −1 ∑

2

2

∥∇ uhn+1/2 ∥ ≤ ∥uh0 ∥ +

n=0

N −1 ∆t ∑ ∥fn+1/2 ∥2∗ . ν

(24)

n=0

Proof. The proof of Lemma 5.1 follows the same idea as that of Lemma 5.2 and a notion about it can be found in Layton’s work [17]. We show below that Eq. (23) is unconditionally stable if solved using the Crank–Nicolson method. Lemma 5.2 (Stability). Consider the approximation scheme given by Eq. (23). A solution Xih , i = 1, . . . , N, exists at each time-step and the scheme is unconditionally stable. Then, it satisfies the following a priori condition: 2

∥XNh ∥ + ∆tDi

N −1 ∑ n=0

2

2

∥∇ Xnh+1/2 ∥ ≤ ∥X0h ∥ +

N −1 ∆t ∑ ∥Sn+1/2 ∥2∗ . ρ 2 Di

(25)

n=0

Please cite this article as: M.M. De Souza and A.L. De Bortoli, Stability and convergence of the numerical solution for the species equation of a model for PEMFCs, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.11.010.

6

M.M. De Souza and A.L. De Bortoli / Computers and Mathematics with Applications xxx (xxxx) xxx

Proof. In Eq. (23) let w h = Xnh+1/2 , then Xnh+1 , Xnh+1/2 − Xnh , Xnh+1/2 + ∆tDi (∇ Xnh+1/2 , ∇ Xnh+1/2 )

)

(

)

(

1

+∆t b(uhn+1 , Xnh+1/2 , Xnh+1/2 ) = ∆t (Sn+1/2 , Xnh+1/2 ). ρ

(26)

Note that Xnh+1 , Xnh+1/2 − Xnh , Xnh+1/2 =

)

(

)

(

) 1 ( h h) 1( h Xn+1 , Xnh+1 − X ,X , 2 2 n n

(27)

and b(uhn+1/2 , Xnh+1/2 , Xnh+1/2 ) = 0, as Temam [13]. Then, Eq. (26) is given by

) 1 ( h h) 1( h − X , Xh X , X + ∆tDi (∇ Xnh+1/2 , ∇ Xnh+1/2 ) 2 n+1 n+1 2 n n 1 = ∆t (Sn+1/2 , Xnh+1/2 ). ρ

Using the norm corresponds, 1

1

1

∥X h ∥2 − ∥Xnh ∥2 + ∆tDi ∥∇ Xnh+1/2 ∥2 = ∆t (Sn+1/2 , Xnh+1/2 ) 2 n+1 2 ρ and using the inequality (12), results 1 2

1

1

∥Xnh+1 ∥2 − ∥Xnh ∥2 + ∆tDi ∥∇ Xnh+1/2 ∥2 ≤ ∆t ∥Sn+1/2 ∥∗ ∥∇ Xnh+1/2 ∥. 2 ρ

(28)

Applying the Young’s inequality with ϵ = Di ρ in the right side of Eq. (28), results in 1 2

1

∆t

2

2ρ 2 Di

∥Xnh+1 ∥2 − ∥Xnh ∥2 + ∆tDi ∥∇ Xnh+1/2 ∥2 ≤ ∆t

∥Sn+1/2 ∥2∗

Di ∥∇ Xnh+1/2 ∥2 . 2 Grouping similar terms, multiplying by 2, and adding terms from n = 0 to n = N − 1, results

+

2

∥XNh ∥ + ∆tDi

N −1 ∑

2

2

∥∇ Xnh+1/2 ∥ ≤ ∥X0h ∥ +

n=0

N −1 ∆t ∑ ∥Sn+1/2 ∥2∗ ρ 2 Di

□

n=0

Below is the convergence theorem of the set of Eqs. (21)–(22). Theorem 5.1 (Convergence). Let (u(t), p(t)) be a sufficiently smooth, strong solution of Eqs. (13)–(14). Suppose that (uh0 , ph0 ) are approximations of (u(0), p(0)) based on the conditions (4). Then, there is a constant C = C (u, p) such that

|||u − uh |||∞,0 ≤ B(∆t , h) + Chk+1 |||u|||∞,k+1 , ( )1/2 N −1 ∑ h 2 ν ∆t ∥∇ (un+1/2 − un+1/2 )∥ ≤ B(∆t , h) + C ν 1/2 ∆t 2 ∥∇ utt ∥2,0

(29)

n=0

+C ν 1/2 hk |||u|||2,k+1 .

(30)

where B(∆t , h) = (∆t , h ): 2

k

B(∆t , h) := C ν −1/2 hk+1/2 (|||u|||24,k+1 +|||∇ u|||24,0 )

+C ν −1/2 hk (|||u|||24,k+1 +ν −1/2 (∥uh0 ∥ + ν −1/2 |||f |||2,∗ ))

+C ν 1/2 hk |||u|||2,k+1 +C ν −1/2 hs+1 |||p1/2 |||2,s+1 +C ∆t 2 (∥uttt ∥2,0 + ν −1/2 ∥ptt ∥2,0 + ∥ftt ∥2,0 + ν 1/2 ∥∇ utt ∥2,0 +ν −1/2 ∥∇ utt ∥24,0 + ν −1/2 |||∇ u|||24,0 +ν −1/2 |||∇ u1/2 |||24,0 ).

(31)

Proof. Available in Layton [17] on page 169. Next, the convergence theorem of Eq. (23) is stated and demonstrated, which is the most important result of this paper. Theorem 5.2 (Convergence). Let X (t) be a sufficiently smooth, strong solution of Eq. (16) satisfying either no-slip or periodic with zero-mean boundary conditions. Suppose that X0h is an approximation of χ (0) accurately according to the conditions (4). So, there Please cite this article as: M.M. De Souza and A.L. De Bortoli, Stability and convergence of the numerical solution for the species equation of a model for PEMFCs, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.11.010.

M.M. De Souza and A.L. De Bortoli / Computers and Mathematics with Applications xxx (xxxx) xxx

7

is a constant C = C (X , u, p) such that

|||X − X h |||∞,0 ≤ A(∆t , h) + Chr +1 |||X |||∞,r +1 , )1/2 ( N −1 ∑ 1/2 2 h ≤ A(∆t , h) + CDi ∆t 2 ∥∇ Xtt ∥2,0 Di ∆t ∥∇ (Xn+1/2 − Xn+1/2 )∥ n=0 1/2

+CDi hr |||X |||2,r +1 .

(32)

(33)

Proof. At time tn+1/2 the true solution of species X satisfies,

(

Xn+1 − Xn

∆t

=

1

ρ

) , wh + Di (∇ Xn+1/2 , ∇wh ) + b(un+1/2 , Xn+1/2 , wh )

(Sn+1/2 , w h ) + Ω (Xn , un ; w h ),

(34)

for all w h ∈ Zh , where Ω (Xn , un ; w h ) is the error of consistency given by

) ∂X (tn+1/2 ), w h ∆t ∂t +Di (∇ Xn+1/2 − ∇ X (tn+1/2 ), ∇wh ) + b(un+1/2 , Xn+1/2 , wh )

Ω (Xn , un ; w h ) =

(

Xn−1 − Xn

−

−b(u(tn+1/2 ), X (tn+1/2 ), wh ) +

1

ρ

(S(tn+1/2 ) − Sn+1/2 , w h ).

Subtracting Eq. (23) from Eq. (34) and writing eXn = Xn − Xnh and eun = un − uhn we have

(

eXn−1 − eXn

∆t

) ,w

h

+ Di (∇ eXn+1/2 , wh ) + b(un+1/2 , Xn+1/2 , wh )

−b(uhn+1/2 , Xnh+1/2 , wh ) = Ω (Xn , un ; wh ),

(35)

and after adding and subtracting b(uhn+1/2 , Xn+1/2 , w h ), Eq. (35) results in eXn−1 − eXn , w h + ∆tDi (∇ eXn+1/2 , w h ) + ∆t b(eun+1/2 , Xn+1/2 , w h )

(

)

+∆t b(uhn+1/2 , eXn+1/2 , wh ) = ∆t Ω (Xn , un ; wh ).

(36)

Decompose the errors as eXn = (Xn − χn ) − (Xnh − χn ) = ηn − φnh and eun = (un − Un ) − (uhn − Un ) = τn − ϕnh , where χ is the L2 projection of X in Zh and U is the L2 projection of u in Vh . Let w h = φnh+1/2 , then (36) can be rewritten as

( h ) ( ) φn+1 − φnh , φnh+1/2 + ∆tDi (∇φnh+1/2 , ∇φnh+1/2 ) = ηn+1 − ηn , φnh+1/2 +∆tDi (∇ηn+1/2 , ∇φnh+1/2 ) + ∆t b(ϕnh+1/2 , Xn+1/2 , φnh+1/2 ) −∆t b(τn , Xn+1/2 , φnh+1/2 ) + ∆t b(uhn+1/2 , φnh+1/2 , φnh+1/2 ) −∆t b(uhn+1/2 , ηn+1/2 , φnh+1/2 ) + ∆t Ω (Xn , un ; φnh+1/2 ),

(37)

The first term of left-hand side of Eq. (37) can be decomposed in the same way from Eq. (27), the term b(uhn+1/2 , φnh+1/2 ,

( ) φnh+1/2 ) = 0 as explained in Temam [13], the term ηn+1 − ηn , φnh+1/2 = 0 since χ is the L2 projection of X in Zh , and after applying the norm we get 1 2

1

∥φnh+1 ∥2 − ∥φnh ∥2 + ∆tDi ∥∇φnh+1/2 ∥2 = ∆tDi (∇ηn+1/2 , ∇φnh+1/2 ) 2

+∆t b(ϕnh+1/2 , Xn+1/2 , φnh+1/2 ) − ∆t b(τn , Xn+1/2 , φnh+1/2 ) −∆t b(uhn+1/2 , ηn+1/2 , φnh+1/2 ) + ∆t Ω (Xn , un ; φnh+1/2 ),

(38)

Now we treat the right side of Eq. (38) term by term. For the first term on the right side we use the Cauchy–Schwarz inequality and then the inequality of Young with ϵ = 1/6, resulting in

∆tDi (∇ηn+1/2 , ∇φnh+1/2 ) ≤ ∆tDi ∥∇ηn+1/2 ∥∥∇φnh+1/2 ∥ ∆tDi ≤ ∥∇φnh+1/2 ∥2 + C ∆tDi ∥∇ηn+1/2 ∥2 . 12

Please cite this article as: M.M. De Souza and A.L. De Bortoli, Stability and convergence of the numerical solution for the species equation of a model for PEMFCs, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.11.010.

8

M.M. De Souza and A.L. De Bortoli / Computers and Mathematics with Applications xxx (xxxx) xxx

In the second, third and fourth terms on the right side of (38), we use lemma 3.1 and inequality of Young with ϵ = Di /6, getting

∆t b(ϕnh+1/2 , Xn+1/2 , φnh+1/2 )

≤ C ∆t ∥∇ϕnh+1/2 ∥∥∇ Xn+1/2 ∥∥∇φnh+1/2 ∥ ∆tDi C ∆t ≤ ∥∇φnh+1/2 ∥2 + ∥∇ϕnh+1/2 ∥2 ∥∇ Xn+1/2 ∥2 , 12

Di

∆t b(τn , Xn+1/2 , φnh+1/2 )

≤ C ∆t ∥τn ∥1/2 ∥∇τn ∥1/2 ∥∇ Xn+1/2 ∥∥∇φnh+1/2 ∥ ∆tDi C ∆t ≤ ∥∇φnh+1/2 ∥2 + ∥τn ∥∥∇τn ∥∥∇ Xn+1/2 ∥2 , 12

Di

∆t b(uhn+1/2 , ηn+1/2 , φnh+1/2 )

≤ C ∆t ∥uhn+1/2 ∥1/2 ∥∇ uhn+1/2 ∥1/2 ∥∇ηn+1/2 ∥∥∇φnh+1/2 ∥ ∆tDi C ∆t h ≤ ∥∇φnh+1/2 ∥2 + ∥un+1/2 ∥∥∇ uhn+1/2 ∥∥∇ηn+1/2 ∥2 . 12

Di

Combining the last estimates with (38), multiplying by 2 and adding from n = 0 to n = N − 1, results in

∥φNh ∥2 − ∥φ0h ∥2 + ∆tDi +

+

N −1 C ∆t ∑

Di

3

N −1 ∑

∥∇φnh+1/2 ∥2 ≤ C ∆tDi

n=0

∥ηn+1/2 ∥2

n=0

∥∇ϕnh+1/2 ∥2 ∥∇ Xn+1/2 ∥2 +

n=0

N −1 C ∆t ∑

Di

N −1 4∑

N −1 C ∆t ∑

Di

∥τn ∥∥∇τn ∥∥∇ Xn+1/2 ∥2

n=0

∥uhn+1/2 ∥∥∇ uhn+1/2 ∥∥∇ηn+1/2 ∥2 + ∆t

n=0

N −1 ∑

|Ω (Xn , un ; φnh+1/2 )|.

(39)

n=0

Using the properties of approximation (4) and definition (6), the first term on the right-hand side of Eq. (39) is given by N −1 ∑

C ∆tDi

∥ηn+1/2 ∥2 ≤ C ∆tDi

n=0

N ∑

∥ηn ∥2 ≤ C ∆tDi

n=0

N ∑

h2r |Xn |2r +1

n=0

≤ CDi h2r |||X |||22,r +1 . To the third term of the right side we use the triangular inequality, the properties (4), inequality of Young and definition (6), which gives N −1 C ∆t ∑

Di

∥τn ∥∥∇τn ∥∥∇ Xn+1/2 ∥2 ≤

n=0

N −1 C ∆t ∑ (∥τn ∥∥∇τn ∥ + ∥τn ∥∥∇τn+1 ∥ Di n=0

+∥τn+1 ∥∥∇τn ∥ + ∥τn+1 ∥∥∇τn+1 ∥)∥∇ Xn+1/2 ∥2 ≤

N −1 C ∆t ∑ k+1 (h |un+1 |k+1 hk |un+1 |k+1 + hk+1 |un |k+1 hk |un |k+1 Di n=0

+2hk+1 |un |k+1 hk |un+1 |k+1 )∥∇ Xn+1/2 ∥2 ≤

N −1 C ∆t ∑

Di

h2k+1 (|un+1 |2k+1 ∥∇ Xn+1/2 ∥2 + |un |2k+1 ∥∇ Xn+1/2 ∥2

n=0

+2|un |k+1 |un+1 |k+1 ∥∇ Xn+1/2 ∥2 ) ≤

≤ ≤

C Di C Di C Di

h2k+1 ∆t

N ∑

|un |2k+1 ∥∇ Xn ∥2

n=0

h2k+1 (∆t

N ∑ n=0

|un |4k+1 + ∆t

N ∑

∥∇ Xn ∥4 )

n=0

h2k+1 (|||u|||44,k+1 +|||∇ X |||44,0 ),

Please cite this article as: M.M. De Souza and A.L. De Bortoli, Stability and convergence of the numerical solution for the species equation of a model for PEMFCs, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.11.010.

M.M. De Souza and A.L. De Bortoli / Computers and Mathematics with Applications xxx (xxxx) xxx

9

and for the fourth term, we use the triangular inequality, the approximation properties (4), Young’s inequality, definition (6) and a priori estimation (Lemma 5.1), resulting in N −1 C ∆t ∑

Di

∥uhn+1/2 ∥∥∇ uhn+1/2 ∥∥∇ηn+1/2 ∥2

n=0

≤

≤

≤

≤

N −1 C ∆t ∑

Di

∥∇ uhn+1/2 ∥∥∇ηn+1/2 ∥2

n=0

N −1 C ∆t ∑

Di

) ( ∥∇ uhn+1/2 ∥ ∥∇ηn+1 ∥2 + ∥∇ηn ∥2

n=0

N −1 C ∆t ∑

Di C Di

) ( ∥∇ uhn+1/2 ∥ h2r |Xn+1 |2r +1 + h2r |Xn |2r +1

n=0

h2r ∆t

N ∑

∥∇ uhn+1/2 ∥|Xn |2r +1 ≤

n=0

(

C

∆t

h2r

Di

N ∑

2

∥∇ uhn+1/2 ∥ + ∆t

n=0

N ∑

) |Xn |4r +1

n=0

[ (

) ] N ∆t ∑ 2 4 ∥ ∥ + ∥fn+1/2 ∥∗ + |||X |||4,r +1 ≤ h Di ν ν n=0 [ ( ) ] C 1 1 2 ≤ h2r ∥uh0 ∥ + |||f |||22,∗ + |||X |||44,r +1 . Di ν ν C

1

2r

2 uh0

Now, let us analyze the term Ω (Xn , un ; φnh+1/2 ), which consists of the error of consistency. In the first and second terms, we use the inequalities of Cauchy–Schwartz and Young, and the lemma 3.2, which gives N −1 ( ∑ Xn−1 − Xn

) ∂X (tn+1/2 ), φnh+1/2 ∆t ∂t n=0   N −1 ∑  h  Xn−1 − Xn ∂X  − (tn+1/2 ) ≤ ∆t  ∥φn+1/2 ∥  ∆t ∂t n=0  2 N −1 N −1 ∑ ∑  1  Xn−1 − Xn ∂X 1 h   + ∆t ≤ ∆t − (t ) ∥φ ∥2 n + 1 / 2   2 ∆t ∂t 2 n+1/2 n=0 n=0 ∫ tn + 1 N −1 N −1 4 ∑ ∑ ∆t 1 h ≤ ∥Xttt ∥2 dt + ∆t ∥φn+1/2 ∥2

∆t

2560

n=0

≤

∆t

∆t 4 2560

−

tn

n=0

|||Xttt |||22,0 +∆t

N −1 ∑ 1 n=0

2

2

∥φnh+1/2 ∥2 .

N −1 ∑

(∇ Xn+1/2 − ∇ X (tn+1/2 ), ∇φnh+1/2 )

n=0

≤ ∆t

N −1 ∑

∥∇ Xn+1/2 − ∇ X (tn+1/2 )∥∥∇φnh+1/2 ∥

n=0

≤ ∆t

N −1 ∑ C n=0

≤

Di 48

C ∆t Di

4

N −1 ∑ Di n=0

∫ N −1 ∑ C ∆t 4 n=0

≤

Di

∥∇ Xn+1/2 − ∇ X (tn+1/2 )∥2 + ∆t tn+1

∥Xtt ∥2 dt + ∆t

tn

|||Xtt |||22,0 +∆t

N −1 ∑ Di n=0

N −1 ∑ Di n=0

6

6

6

∥∇φnh+1/2 ∥2

∥∇φnh+1/2 ∥2

∥∇φnh+1/2 ∥2 .

Please cite this article as: M.M. De Souza and A.L. De Bortoli, Stability and convergence of the numerical solution for the species equation of a model for PEMFCs, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.11.010.

10

M.M. De Souza and A.L. De Bortoli / Computers and Mathematics with Applications xxx (xxxx) xxx

After analyzing the subtraction of the non-linear terms, we add and subtract a term and then apply the lemma 3.1, the Young’s inequality with ϵ = Di /3 and the lemma 3.2, to obtain

∆t

N −1 ∑

b(un+1/2 , Xn+1/2 , φnh+1/2 ) − ∆t

N −1 ∑

n=0

b(u(tn+1/2 ), X (tn+1/2 ), φnh+1/2 )

n=0 N −1

≤ ∆t

∑

b(un+1/2 − u(tn+1/2 ), Xn+1/2 , φnh+1/2 )

n=0

+∆t

N −1 ∑

b(u(tn+1/2 ), Xn+1/2 − X (tn+1/2 ), φnh+1/2 )

n=0

≤ ∆t

N −1 ∑ (

C ∥∇ (un+1/2 − u(tn+1/2 ))∥∥∇ Xn+1/2 ∥

n=0

) +C ∥∇ u(tn+1/2 )∥∥∇ (X (tn+1/2 ) − Xn+1/2 )∥ ∥∇φnh+1/2 ∥ N −1 ( ∑ C ∥∇ (un+1/2 − u(tn+1/2 ))∥2 ∥∇ Xn+1/2 ∥2 ≤ ∆t Di

n=0

+

C Di

≤ ∆t

) N −1 ∑ Di ∥∇ u(tn+1/2 )∥2 ∥∇ (X (tn+1/2 ) − Xn+1/2 )∥2 + ∆t ∥∇φnh+1/2 ∥2 n=0

( ∫ N −1 ∑ C ∆t 3 ∥∇ Xn+1/2 ∥2 Di

n=0

48

+∥∇ u(tn+1/2 )∥

≤

( N −1 ∑ C ∆t 4 Di 96

n=0 tn+1

∫ +

tn

≤

2 ∆t

3

48

tn+1

∫

∥∇ utt ∥2 dt

tn

)

∥∇ Xtt ∥ dt + ∆t 2

tn

N −1 ∑ Di n=0

∥∇ Xn+1/2 ∥4 +

∫

tn+1

6

∥∇φnh+1/2 ∥2

∥∇ utt ∥4 dt + ∥∇ u(tn+1/2 )∥4

tn

) N −1 ∑ Di ∥∇ Xtt ∥4 dt + ∆t ∥∇φnh+1/2 ∥2 n=0

C ∆t 4 ( Di 96

+∆t

tn+1

6

6

∥∇ Xn+1/2 ∥4 + ∥∇ utt ∥44,0 + ∥∇ u(tn+1/2 )∥4 + |||∇ Xtt |||44,0

N −1 ∑ Di n=0

)

∥∇φnh+1/2 ∥2 .

6

In the source term, we use the inequalities of Cauchy–Schwartz and Young, and the lemma 3.2, which gives

∆t

N −1 ∑

(S(tn+1/2 ) − Sn+1/2 , φnh+1/2 )

n=0

≤ ∆t

N −1 ∑

∥S(tn+1/2 ) − Sn+1/2 ∥∥φnh+1/2 ∥

n=0

≤ ∆t

N −1 ∑ 1 n=0

≤

∫ N −1 ∑ ∆t 4 n=0

≤

2

∥S(tn+1/2 ) − Sn+1/2 ∥2 + ∆t

∆t

4

96

96

N −1 ∑ 1 n=0

tn + 1

∥Stt ∥2 dt + ∆t

tn

|||Stt |||22,0 +∆t

N −1 ∑ 1 n=0

N −1 ∑ 1 n=0

2

2

2

∥φnh+1/2 ∥2

∥φnh+1/2 ∥2

∥φnh+1/2 ∥2 .

Please cite this article as: M.M. De Souza and A.L. De Bortoli, Stability and convergence of the numerical solution for the species equation of a model for PEMFCs, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.11.010.

M.M. De Souza and A.L. De Bortoli / Computers and Mathematics with Applications xxx (xxxx) xxx

11

Combining the estimates, the consistency error in Eq. (39) is given by

∆t

N −1 ∑

|Ω (Xn , un ; φnh+1/2 )| ≤

n=0

N −1 ∑

∆t ∥φnh+1/2 ∥ +

n=0

N −1 ∑

∆t

Di

n=0

3

∥φnh+1/2 ∥2

) ( ∆t 4 ( +C ∆t |||Xttt |||22,0 +|||Stt |||22,0 + C |||Xtt |||22,0 +∥∇ Xn+1/2 ∥4 Di ) +|||∇ utt |||44,0 +∥∇ u(tn+1/2 )∥4 + |||∇ Xtt |||44,0 . 4

(40)

h 0

Substituting (40) in (39) and assuming that ∥φ ∥ = 0, results

∥φNh ∥2 + ∆tDi

N −1 ∑

∥∇φnh+1/2 ∥2 ≤

n=0

+ +

N −1 C ∆t ∑

Di C

∆t ∥φnh+1/2 ∥ + CDi h2r |||X |||22,r +1

n=0

C

∥∇ϕnh+1/2 ∥2 ∥∇ Xn+1/2 ∥2 +

Di

n=0

h

Di

N −1 ∑

[

2r

1

ν

uh0 2

(∥

∥ +

1

ν

2 2,∗ )

|||f |||

h2k+1 (|||u|||44,k+1 +|||∇ X |||44,0 )

4 4,r +1

]

+ |||X |||

∆t 4

(Di |||Xttt |||22,0 +Di |||Stt |||22,0 +|||Xtt |||22,0 Di +∥∇ Xn+1/2 ∥4 + |||∇ utt |||44,0 +∥∇ u(tn+1/2 )∥4 + |||∇ Xtt |||44,0 ).

+C

Knowing that ∆t is sufficiently small, i.e. ∆t < C , we use the discrete lemma of Gronwall, to obtain

∥φNh ∥2 + ∆tDi

N −1 ∑

∥∇φnh+1/2 ∥2 ≤ CDi h2r |||X |||22,r +1

n=0

+ +

N −1 C ∆t ∑

Di C

n=0

h

Di

∥∇ϕnh+1/2 ∥2 ∥∇ Xn+1/2 ∥2 +

2r

[

1

ν

(∥uh0 ∥2 +

1

ν

C Di

h2k+1 (|||u|||44,k+1 +|||∇ X |||44,0 )

|||f |||22,∗ ) + |||X |||44,r +1

]

∆t 4

(Di |||Xttt |||22,0 +Di |||Stt |||22,0 +|||Xtt |||22,0 Di +∥∇ Xn+1/2 ∥4 + |||∇ utt |||44,0 +∥∇ u(tn+1/2 )∥4 + |||∇ Xtt |||44,0 ).

+C

(41)

Note that C ∆t

N −1 ∑

∥∇ϕnh+1/2 ∥2 ≤ C ∆t

n=0

N −1 ∑

∥∇τn+1/2 ∥2 + C ∆t

n=0

≤ C ∆t

N ∑

h2k |∇ un |2k+1 + C ∆t

n=0

N −1 ∑

∥∇ eun+1/2 ∥2

n=0 N −1 ∑

∥∇ eun+1/2 ∥2

n=0

≤ Ch2k |||∇ u|||22,k+1 +C ∆t

N −1 ∑

∥∇ eun+1/2 ∥2 ,

(42)

n=0

The second term on the right-hand side of Eq. (42) can be described by inequality (30), i.e. C ∆t

N −1 ∑

2

∥∇ϕnh+1/2 ∥ ≤ Ch2k |||∇ u|||22,k+1

n=0

+

C

ν

1/2

(

B (∆t , h) + C ν 1/2 ∆t 2 ∥utt ∥2,0 + C ν 1/2 hk |||u|||2,k+1

)2

.

Thus, Eq. (41) results in

∥φNh ∥2 + ∆tDi

N −1 ∑

∥∇φnh+1/2 ∥2 ≤ CDi h2r |||X |||22,r +1

n=0

+

C Di

[

h2k |||∇ u|||22,k+1 +

(

C

ν 1/2

B(∆t , h) + C ∆t 2 ∥utt ∥2,0

Please cite this article as: M.M. De Souza and A.L. De Bortoli, Stability and convergence of the numerical solution for the species equation of a model for PEMFCs, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.11.010.

12

M.M. De Souza and A.L. De Bortoli / Computers and Mathematics with Applications xxx (xxxx) xxx

)2 ] C +Chk |||u|||2,k+1 ∥∇ Xn+1/2 ∥2 + h2k+1 (|||u|||44,k+1 +|||∇ X |||44,0 ) Di ] [ C 2r 1 1 h 2 2 + h (∥u0 ∥ + |||f |||2,∗ ) + |||X |||44,r +1 Di ν ν ∆t 4 ( +C Di |||Xttt |||22,0 +Di |||Stt |||22,0 +|||Xtt |||22,0 +∥∇ Xn+1/2 ∥4 Di ) +|||∇ utt |||44,0 +∥∇ u(tn+1/2 )∥4 + |||∇ Xtt |||44,0 .

(43) h M

eXM

To obtain the first inequality (32) we apply the triangular inequality, i.e. ∥φ ∥ = ∥ use (4),

− ηM ∥ ≥ ∥

eXM

∥ − ∥ηM ∥, then we

C

|||eX |||2∞,0 ≤ CDi h2r |||X |||22,r +1 + h2k+1 (|||u|||44,k+1 +|||∇ X |||44,0 ) D [ ( i C C 2k 2 h |||∇ u|||2,k+1 + B(∆t , h) + C ∆t 2 ∥utt ∥2,0 + Di ν 1/2 )2 ] +Chk |||u|||2,k+1 ∥∇ Xn+1/2 ∥2 [ ] C 2r 1 1 h 2 2 4 + h (∥u0 ∥ + |||f |||2,∗ ) + |||X |||4,r +1 Di ν ν 4 ( ∆t Di |||Xttt |||22,0 +Di |||Stt |||22,0 +|||Xtt |||22,0 +∥∇ Xn+1/2 ∥4 +C Di ) +|||∇ utt |||44,0 +∥∇ u(tn+1/2 )∥4 + |||∇ Xtt |||44,0 + |||η|||2∞,0 . Taking the square root, we obtain C

1/2

|||X − X h |||∞,0 ≤ CDi hr |||X |||2,r +1 + +|||∇ X |||24,0 ) +

C

[

1/2

Di

1/2

Di

hk+1/2 (|||u|||24,k+1

hk |||∇ u|||2,k+1 +

C

ν 1/2

B(∆t , h) + C ∆t 2 ∥utt ∥2,0

] +Chk |||u|||2,k+1 ∥∇ Xn+1/2 ∥ ] [ 1 C r 1 h 2 (∥u0 ∥ + 1/2 |||f |||2,∗ ) + |||X |||4,r +1 + 1/2 h ν 1/2 ν Di +C

∆t 2 1/2

Di

1/2

(Di

1/2

|||Xttt |||2,0 +Di |||Stt |||2,0 +|||Xtt |||2,0 +∥∇ Xn+1/2 ∥2

+|||∇ utt |||24,0 +∥∇ u(tn+1/2 )∥2 + |||∇ Xtt |||24,0 ) + Chr +1 |||X |||∞,r +1 .

(44)

Define A(∆t , h) such that C

1/2 r

A(∆t , h) := CDi

+

C

h |||X |||2,r +1 +

[

hk |||∇ u|||2,k+1 +

1/2 Di k

1/2

Di

C

ν 1/2

hk+1/2 (|||u|||24,k+1 +|||∇ X |||24,0 )

B(∆t , h) + C ∆t 2 ∥utt ∥2,0

] +Ch |||u|||2,k+1 ∥∇ Xn+1/2 ∥2 [ ] 1 1 C r h 2 + 1/2 h (∥u0 ∥ + 1/2 |||f |||2,∗ ) + |||X |||4,r +1 ν 1/2 ν Di +C

∆t 2 1/2 Di

1/2

(Di

1/2

|||Xttt |||2,0 +Di |||Stt |||2,0 +|||Xtt |||2,0 +∥∇ Xn+1/2 ∥2

+|||∇ utt |||24,0 +∥∇ u(tn+1/2 )∥2 + |||∇ Xtt |||24,0 ). Thus, we have the result (32). For the second inequality of (33) we use the following,

∥∇ (X (tn+1/2 ) − Xnh+1/2 )∥2 ≤ ∥∇ X (tn+1/2 ) − Xn+1/2 ∥2 + ∥∇ηn+1/2 ∥2 + ∥φnh+1/2 ∥2 ∫ ∆ t 3 tn + 1 ≤ ∥∇ Xtt ∥2 dt + Ch2r |Xn+1 |2r +1 + Ch2r |Xn |2r +1 + ∥φnh+1/2 ∥2 . 48

tn

Please cite this article as: M.M. De Souza and A.L. De Bortoli, Stability and convergence of the numerical solution for the species equation of a model for PEMFCs, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.11.010.

M.M. De Souza and A.L. De Bortoli / Computers and Mathematics with Applications xxx (xxxx) xxx

13

Therefore, Eq. (43) can be written as

∆tDi

N −1 ∑

∥∇ (Xn+1/2 − Xnh+1/2 )∥2 ≤ CDi h2r |||X |||22,r +1 +

n=0

+|||∇ X |||44,0 ) + +Chk |||u|||2,k+1

C

[

Di

)2 ]

h2k |||∇ u|||22,k+1 +

∥∇ Xn+1/2 ∥2 +

(

C Di

C

ν 1/2 [

h2r

C Di

h2k+1 (|||u|||44,k+1

B(∆t , h) + C ∆t 2 ∥utt ∥2,0 1

ν

(∥uh0 ∥2 +

1

ν

|||f |||22,∗ )

] ∆t 4 ( +|||X |||44,r +1 + C Di |||Xttt |||22,0 +Di |||Stt |||22,0 +|||Xtt |||22,0 +∥∇ Xn+1/2 ∥4 Di

N −1 ∑ ) ∆t 4 +|||∇ utt |||44,0 +∥∇ u(tn+1/2 )∥4 + |||∇ Xtt |||44,0 + Di n=0

+∆tDi

N −1 ∑

Ch2r |Xn+1 |2r +1 + ∆tDi

n=0

N −1 ∑

48

∫

tn+1

∥∇ Xtt ∥2 dt

tn

Ch2r |Xn |2r +1 ,

n=0

i.e.

( ∆tDi

N −1 ∑

)1/2 ∥∇ (Xn+1/2 −

Xnh+1/2 ) 2

∥

1/2

≤ A(∆t , h) + CDi ∆t 2 |||∇ Xtt |||2,0

n=0 1/2

+CDi hr |||X |||2,r +1 . □

(45)

What has just been proved is important because it guarantees the convergence rate of the numerical solution of the species equation using the Crank–Nicolson method. The smaller the mesh, the smaller the error. Consequently, there will be a better approximation; on the other hand, the computational cost will be higher. The inequalities (44) and (45) of the convergence theorem are in function of ∆t and h. As expected, with the use of the Crank–Nicolson method, the time dependence must be of order 2, ∆t 2 . In addition, the species equation depends on the velocity, so it also depends on the pressure. Thus, in space the dependence is hr for each chemical species, hk for speed, and hs+1 for pressure. In addition, the constant C depends of the velocity, pressure and species, i.e., C (X , u, p). In the following corollary, we assume degree 2 for velocity and species, and degree 1 for pressure. Corollary 5.1. Consider the hypothesis of the theorems 5.1 and 5.2, suppose that V h and Q h are composed of Taylor–Hood elements, Z h of quadratic elements and the use of the Crank–Nicolson method for time. The error is of the order

( |||X − X h |||∞,0 + ∆tDi

N −1 ∑

)1/2 ∥∇ (Xn+1/2 − Xnh+1/2 )∥2

= O(h2 + ∆t 2 ).

n=0

6. Numerical results A code was developed in Fortran 90 for the solution of Eqs. (21)–(23). First, we validate the computational code with exact solution whose convergence rate was proved in Theorem 5.2; after that we show some results for a fuel cell. 6.1. Problem with exact solution In the same way as John et al. [25] do in their article, to validate the computational code for the species equation, we will consider an example with a dominating space error that the exact solutions of Eqs. (13), (14) and (16) are given by u1 = t 3 y2 , u2 = t 2 x , Xi = t 3 y2 , in domain Ω = (0, 1)2 , a unit square. The source term Si , the initial condition and the boundary conditions are obtained from Xi , which is the solution of Eq. (16) for ρ = 1 and Deff i = 1. We show the reliability of the computational code through convergence rate. Note that we are working with quadratic elements for the equation of species and Taylor–Hood for equation of velocity and pressure. In this way, the expected species convergence rate is two, as proved in the theorem of the previous section, Eq. (33). In Table 2, the numerical error and convergence rate are shown. The more refined the mesh, the less numerical error, i.e., adding degrees of freedom, the numerical solution tends to the exact solution. Therefore, expected convergence rate was obtained for species. Please cite this article as: M.M. De Souza and A.L. De Bortoli, Stability and convergence of the numerical solution for the species equation of a model for PEMFCs, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.11.010.

14

M.M. De Souza and A.L. De Bortoli / Computers and Mathematics with Applications xxx (xxxx) xxx Table 2 Numerical error and convergence rate. h

∥∇ (Xi − Xih )∥0

Rate

1/4 1/16 1/36 1/64

1.7933e−3 1.1423e−4 2.2567e−5 7.1405e−6

– 1.986304 1.999822 1.999977

Fig. 2. Mesh for the fuel cell.

6.2. Direct ethanol fuel cell (DEFC) We consider a PEMFC fueled with ethanol, that is, a direct ethanol fuel cell (DEFC). The electrochemical reactions that occur in this type of cell are presented in (1) and (2). A DEFC works well and is considered a promising energy transformation equipment, due to the higher energy density of ethanol. In addition, ethanol is less toxic than methanol and does not have the storage limitations of hydrogen. The DEFCs belong to the family of proton exchange membrane fuel cells (PEMFCs), in which ethanol is used directly as fuel [26]. The DEFC functionality follows the same form of PEMFC, i.e., a fuel is oxidized electrochemically at the anode, and generates protons and electrons, while oxygen is reduced at the cathode [1]. Consider a catalyst PtSn/C at the anode and Pt/C at the cathode, and Nafion 117 membrane. Fig. 2 shows the regular mesh of the fuel cell with a grid of h = 0.2 cm, i.e., 2358 degrees of freedom for each velocity component and for each species, and 488 degrees of freedom for pressure. The mesh of Fig. 2 was obtained after a series of mesh tests, and gave us a good solution with small simulation time. The total flow time is T = 30 min and the time-step is ∆t = 0.1 s. Please cite this article as: M.M. De Souza and A.L. De Bortoli, Stability and convergence of the numerical solution for the species equation of a model for PEMFCs, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.11.010.

M.M. De Souza and A.L. De Bortoli / Computers and Mathematics with Applications xxx (xxxx) xxx

15

The initial conditions (t = 0 s) on the anode side are defined as follows: u1 = 0.02 cm s−1 ; u2 = 0; p = 1 bar; XEtOH = 0.25, and on the cathode side are u1 = 2.55 cm s−1 ; u2 = 0; p = 1 bar; XO2 = 0.21. The boundary conditions in the fuel cell inlet channel, on the anode side, are u1 = 0.02 cm s−1 ; u2 = 0;

∂p = 0; XEtOH = 0.25, ∂x

and on the exit channel are

∂ u2 ∂ XEtOH ∂ u1 = 0; = 0; p = 1 bar; = 0. ∂x ∂x ∂x The boundary conditions on the input channel on the cathode side are u1 = 2.55 cm s−1 ; u2 = 0;

∂p = 0; XO2 = 0.21, ∂x

and on the exit channel are

∂ XO2 ∂ u1 ∂ u2 = 0; = 0; p = 1 bar; = 0, ∂x ∂x ∂x and for all other boundaries of the fuel cell, the speed is zero. Pressure and species have the Neumann condition equal to zero. Fig. 3 shows the mole fraction of ethanol on the anode and oxygen on the cathode, for 120 mA/cm2 current density. The mixture of ethanol and water is inserted into the anode side with a temperature of 348 K and ethanol concentration of 1 M. The ethanol has flow rate of 1 ml/min and pressure of 1 bar on the anode side, and the air flow rate is 120 ml/min and pressure of 1 bar on the cathode side [23,27]. The mole fraction of the ethanol decreases near the catalyst, where oxidation reactions occur. The mole fraction of oxygen decreases near the catalyst, where the reduction reactions occur. The isolines on the anode side and on the cathode side are different due to the physical condition of the fluid (liquid or gaseous). The results of the molar fraction are similar to the expected values. Often, we are interested in the effective voltage that a fuel cell can deliver with a specific current density. The overall voltage of the cell can be obtained using the following relation: 0 Vcell = Ecell − (ηact + ηohm + ηcon ),

(46)

where ηact are the losses due to the activation, ηohm are the losses due to the ohmic resistance, ηcon are the losses due to 0 concentration, Vcell is the cell voltage, and Ecell is the theoretical potential of the DEFC, which is given by [27] 0 Ecell =−

∆G0

∼ 1.14 V,

nF

where ∆G0 is Gibbs free energy. The overpotential losses indicate that an irreversible chemical process occurs inside fuel cells, being responsible for the decrease of the electric potential. The equation of the model for concentration overpotential at the anode is given by [28]:

ηcon,a =

(

⎛

)

RT

ln ⎝

αa nF

⎞ J

(

u

)⎠ ,

(47)

J0 e kd + B

where J0 is the exchanged current density, αa is the transfer coefficient in the anode, n is the number of electrons transferred, F is the Faraday constant, R is the universal gas constant, T is the temperature and the constant B is given by B=

J nFuCEtOH

(

eu/kd (u

) , /kf ) + 1 − 1

in which CEtOH is the concentration of ethanol in the catalyst layer and u is the velocity calculated by Eq. (13), kd is the mass transfer coefficient of ethanol and kf is the mass transfer of the feed stream to the diffusion layer. The expression for the concentration overpotential at the cathode is given by [29]:

ηcon,c =

(

RT

αc nF

)

( ln

Jc J0 (1 − Jc N)

)

,

(48)

Please cite this article as: M.M. De Souza and A.L. De Bortoli, Stability and convergence of the numerical solution for the species equation of a model for PEMFCs, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.11.010.

16

M.M. De Souza and A.L. De Bortoli / Computers and Mathematics with Applications xxx (xxxx) xxx

Fig. 3. Contours of mole fractions of ethanol and oxygen in the fuel cell for 120 mA/cm2 .

where Jc is the current density of the cathode, αc is the cathodic transfer coefficient and N is a limiting current density, defined as δdl,c CO2 , N = nFDeff O 2

in which CO2 is the concentration of oxygen in the catalyst layer and DO2 the oxygen diffusivity in porous media. The relationship between activation overpotential and the current density in the anode was reported in the work of Pramanik and Basu [28] and Gomes and De Bortoli [29],

ηact ,a =

RT

αa nF

[ ln

J(CEtOH CH2 O )−0.25 KEtOH

]

,

(49)

where CH2 O is the concentration of water in the catalyst layer and KEtOH is a global parameter for the oxidation of ethanol. Activation losses at the cathode can be obtained from the Tafel equation [15], resulting in:

ηact ,c =

RT

αc nF

( ) ln

Jc

J0

,

(50)

where αc is the transfer coefficient of the cathode. The ohmic losses are due to the contact of the components (ηcont ) of the fuel cell and resistance to the transport of protons across the membrane (ηmemb ). Ohmic losses can be written as [27]

ηohm = ηcont + ηmemb , where

ηcont =

δm + δcl,a + δcl,c + δdl,a + δdl,c Kseff

J

Please cite this article as: M.M. De Souza and A.L. De Bortoli, Stability and convergence of the numerical solution for the species equation of a model for PEMFCs, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.11.010.

M.M. De Souza and A.L. De Bortoli / Computers and Mathematics with Applications xxx (xxxx) xxx

17

Table 3 Parameters used in the model of the DEFC. Parameter

Unit

Present model

Reference

κO2 κH2 O µO2 µ H2 O ρO2 ρH2 O

cm2 cm2 g cm−1 s−1 g cm−1 s−1 g cm−3 g cm−3 Coulomb mol−1 cm cm cm cm cm cm cm dimensionless dimensionless dimensionless dimensionless dimensionless g mol−1 g mol−1 g cm−3 cm2 s−1 cm2 s−1 mAcm−2 Coulomb (mol cm)−0.5 s−1 S cm−1 S cm−1

1.76 × 10−7 1.0 × 10−7 2.05 × 10−5 0.04061 0.0012 1 96,487 0.00145 0.001 0.001 1.0 1.0 1.0 1.0 0.65 0.4 0.1 0.0475 12 46.06844 32 1 1.19 × 10−6 0.29 3.6 1 8.13 × 106 0.1416

[30] [30] [23] [23] [30] [30] [28] [28] [31] [31] assumed assumed assumed assumed [31] [32] [28] [33] Eq. (3) [29] [29] [29] [31] [28] [33] [34] [35] [36]

F

δm δcl,a δcl,c δdl,a δdl,c δfc ,a δfc ,c εdl εcl αc αa n MEtOH MO 2

ρH2 O

DEtOH DO 2 j0 KEtOH Ks Km

Fig. 4. Current density versus cell voltage and power density for temperature of 348 K and ethanol concentration of 1 M.

and

ηmemb =

δm Kmeff

J,

where Kseff and Kmeff are effective conductivity of the solid phase and effective conductivity of protons in the catalyst layer, respectively. The parameters used to solve the equations of the model for DEFC are listed in Table 3. Fig. 4 shows the overpotential losses of a fuel cell with the catalyst PtSn/C on the anode side and Pt/C on the cathode side [37], for feed of 1 M ethanol concentration and 348 K temperature. The numerical results obtained were compared with the experimental values of Zhao et al. [38]. Fig. 4 shows the cell voltage (calculated by Eq. (46)) and cell power for a given current density. It is observed that the results of the model are consistent with the experimental values. The higher the current density, the lower the voltage of the cell, for the conditions investigated in this work. This loss voltage occurs due to activation overpotential, ohmic losses and concentration overpotentials. In addition, the peak power density achieved is approximately 40 mW/cm2 and occurs between current density values of 110 and 120 mA/cm2 . For low current densities, Please cite this article as: M.M. De Souza and A.L. De Bortoli, Stability and convergence of the numerical solution for the species equation of a model for PEMFCs, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.11.010.

18

M.M. De Souza and A.L. De Bortoli / Computers and Mathematics with Applications xxx (xxxx) xxx

even if the difference between the experimental and the numerical result for voltage is high, the power density will not be high, since it corresponds to the multiplication of the voltage by the current density. 7. Conclusions In this work, we developed the numerical analysis of the equation of the species applied in the modeling of a PEMFC. The demonstrations of stability and convergence presented in this work are unpublished. Through the convergence theorem applied to the species equation, it is observed that the rate of convergence of the species is related to velocity and pressure, since it depends on the degree of the polynomials used for velocity and pressure discretizations. Therefore, to improve the performance of the algorithm it is appropriate to choose polynomials of the same degree for species and velocity, and a lower degree for pressure. Finally, this work contributes to the research of fuel cell models and the development of this technology. The demonstrations presented in this paper give greater credibility to the results that are presented in this area of study. In addition, the model presents results of voltage and power density of a DEFC that are in accordance with experimental data found in the literature, reinforcing the validity of the mathematical model. Nomenclature D diffusion coefficient (cm2 s−1 ) Vcell cell voltage (V) 0 Ecell reversible voltage (V) Xi mole fraction F Faraday’s constant (C mol−1 ) R gas constant (J mol−1 K−1 ) u velocity (cm s−1 ) J transfer current density (mA cm−2 ) j volumetric current density (mA cm−3 ) J0 exchange current density at anode and cathode (mA cm−2 ) p pressure (bar) M molecular weight (g mol−1 ) n number of electrons transferred

δ α ε ρ ν η κ

Greek symbols: thickness (cm) transfer coefficient porosity density (g cm−3 ) viscosity (cm2 s−1 ) overpotential (V) permeability (cm2 )

k a c m dl cl fc act ohm con EtOH O2

Subscripts and superscripts: species anode cathode membrane diffusion layer catalyst layer flow channel activation ohmic concentration ethanol oxygen

Acknowledgments This research is being developed at the Federal University of Rio Grande do Sul - UFRGS. De Souza thanks the financial support of the Instituto Federal de Educação, Ciência e Tecnologia de Santa Catarina - IFSC - Brazil [grant number 23292.020757/2016-52]. Professor De Bortoli thanks the financial support of CNPQ - Conselho Nacional de Desenvolvimento Científico e Tecnológico [grant number 303816/2015-5]. Please cite this article as: M.M. De Souza and A.L. De Bortoli, Stability and convergence of the numerical solution for the species equation of a model for PEMFCs, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.11.010.

M.M. De Souza and A.L. De Bortoli / Computers and Mathematics with Applications xxx (xxxx) xxx

19

References [1] C. Lamy, C. Coutanceau, J.M. Leger, The direct ethanol fuel cell: a challenge to convert bioethanol cleanly into electric energy, in: Catalysis for Sustainable Energy Production, Wiley-VCH Verlag GmbH & Co. KGaA, 2009, pp. 1–46. [2] Z. Zakaria, S.K. Kamarudin, S.N. Timmiati, Membranes for direct ethanol fuel cells: An overview, Appl. Energy 163 (2016) 334–342. [3] C. Siegel, Review of computational heat and mass transfer modeling in polymer-electrolyte-membrane (PEM) fuel cells, Energy 33 (9) (2008) 1331–1352. [4] Y.-X. Wang, K. Ou, Y.-B. Kim, Modeling and experimental validation of hybrid proton exchange membrane fuel cell/battery system for power management control, Int. J. Hydrogen Energy 40 (35) (2015) 11713–11721. [5] J. Goel, S. Basu, Mathematical modeling and experimental validation of direct ethanol fuel cell, Int. J. Hydrogen Energy 40 (41) (2015) 14405–14415. [6] A. Iranzo, M. Muoz, F. Rosa, J. Pino, Numerical model for the performance prediction of a PEM fuel cell. model results and experimental validation, Int. J. Hydrogen Energy 35 (20) (2010) 11533–11550. [7] W. Gong, Z. Cai, Parameter optimization of PEMFC model with improved multi-strategy adaptive differential evolution, Eng. Appl. Artif. Intell. 27 (2014) 28–40. [8] C. Carral, P. Mélé, A numerical analysis of PEMFC stack assembly through a 3D finite element model, Int. J. Hydrogen Energy 39 (9) (2014) 4516–4530. [9] G. Hu, J. Fan, S. Chen, Y. Liu, K. Cen, Three-dimensional numerical analysis of proton exchange membrane fuel cells (PEMFCs) with conventional and interdigitated flow fields, J. Power Sources 136 (1) (2004) 1–9. [10] H. Wu, X. Li, P. Berg, Numerical analysis of dynamic processes in fully humidified PEM fuel cells, Int. J. Hydrogen Energy 32 (12) (2007) 2022–2031. [11] S. Brenner, R. Scott, The Mathematical Theory of Finite Element Methods, vol. 15, third ed., Texts in Applied Mathematics, Springer-Verlag, New York, 2008. [12] C. Johnson, Numerical Solution of Partial Differential Equations by the Finite Element Method, Dover Books on Mathematics Series, Dover Publications, 1987. [13] R. Temam, The Navier-Stokes Equations: Theory and Numerical Analysis, CHEL Series, American Math. Soc., 2001. [14] M.A.R.S. Al-Baghdadi, Modelling of proton exchange membrane fuel cell performance based on semi-empirical equations, Renew. Energy 30 (10) (2005) 1587–1599. [15] S. Abdullah, S.K. Kamarudin, U.A. Hasran, M.S. Masdar, W.R.W. Daud, Development of a conceptual design model of a direct ethanol fuel cell (DEFC), Int. J. Hydrogen Energy 40 (35) (2015) 11943–11948. [16] S. Kundu, L.C. Simon, M. Fowler, S. Grot, Mechanical properties of Nafion electrolyte membranes under hydrated conditions, Polymer (ISSN: 00323861) 46 (25) (2005) 11707–11715. [17] W. Layton, Introduction to the Numerical Analysis of Incompressible Viscous Flows, in: Computational Science and Engineering Series, Society for Industrial and Applied Mathematics, 2008. [18] W. Layton, C.C. Manica, M. Neda, L.G. Rebholz, Numerical analysis and computational testing of a high accuracy Leray-deconvolution model of turbulence, Numer. Methods Partial Differential Equations 24 (2) (2008) 555–582. [19] A. Labovsky, W.J. Layton, C.C. Manica, M. Neda, L.G. Rebholz, The stabilized extrapolated trapezoidal finite-element method for the Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. 198 (9) (2009) 958–974. [20] W. Layton, C.C. Manica, M. Neda, L.G. Rebholz, Numerical analysis and computational comparisons of the NS-alpha and NS-omega regularizations, Comput. Methods Appl. Mech. Engrg. 199 (13) (2010) 916–931. [21] W. Liu, C.-Y. Wang, Three-dimensional simulations of liquid feed direct methanol fuel cells, J. Electrochem. Soc. 154 (3) (2007) B352 – B361. [22] A.D. Le, B. Zhou, A general model of proton exchange membrane fuel cell, J. Power Sources 182 (1) (2008) 197–222. [23] Z.H. Wang, C.Y. Wang, Mathematical modeling of liquid-feed direct methanol fuel cells, J. Electrochem. Soc. 150 (4) (2003) A508–A519. [24] V. John, J. Rang, Adaptive time step control for the incompressible Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. 199 (9) (2010) 514–524. [25] V. John, G. Matthies, J. Rang, A comparison of time-discretization/linearization approaches for the incompressible Navier-Stokes equations, Comput. Methods Appl. Mech. Engrg. 195 (44) (2006) 5995–6010. [26] S. Song, G. Wang, W. Zhou, X. Zhao, G. Sun, Q. Xin, S. Kontou, P. Tsiakaras, The effect of the MEA preparation procedure on both ethanol crossover and DEFC performance, J. Power Sources 140 (1) (2005) 103–110. [27] G.M. Andreadis, A.K.M. Podias, P.E. Tsiakaras, The effect of the parasitic current on the direct ethanol PEM fuel cell operation, J. Power Sources 181 (2) (2008) 214–227. [28] H. Pramanik, S. Basu, Modeling and experimental validation of overpotentials of a direct ethanol fuel cell, Chem. Eng. Process. 49 (7) (2010) 635–642. [29] R.S. Gomes, A.L. De Bortoli, A three-dimensional mathematical model for the anode of a direct ethanol fuel cell, Appl. Energy 183 (2016) 1292–1301. [30] J. Ge, H. Liu, A three-dimensional mathematical model for liquid-fed direct methanol fuel cells, J. Power Sources 160 (1) (2006) 413–421. [31] N.S. Suresh, S. Jayanti, Cross-over and performance modeling of liquid-feed polymer electrolyte membrane direct ethanol fuel cells, Int. J. Hydrogen Energy 36 (22) (2011) 14648–14658. [32] D. Kareemulla, S. Jayanti, Comprehensive one-dimensional, semi-analytical, mathematical model for liquid-feed polymer electrolyte membrane direct methanol fuel cells, J. Power Sources 188 (2) (2009) 367–378. [33] G. Andreadis, S. Song, P. Tsiakaras, Direct ethanol fuel cell anode simulation model, J. Power Sources 157 (2) (2006) 657–665. [34] K. Scott, P. Argyropoulos, K. Sundmacher, A model for the liquid feed direct methanol fuel cell, J. Electrochem. Soc. 477 (2) (1999) 97–110. [35] S.F. Baxter, V.S. Battaglia, R.E. White, Methanol fuel cell model: Anode, J. Electrochem. Soc. 146 (2) (1999) 437–447. [36] K.T. Jeng, C.W. Chen, Modeling and simulation of a direct methanol fuel cell anode, J. Power Sources 112 (2) (2002) 367–375. [37] W.J. Zhou, Z. Zhou, S.Q. Song, W.Z. Li, G.Q. Sun, Q. Xin, Pt-based anode catalysts for direct ethanol fuel cells, Appl. Catal. B: Environ. 46 (2003) 273–285. [38] X.S. Zhao, L.H. Jiang, G.Q. Sun, S.H. Yang, B.L. Yi, Q. Xin, Electrocatalytic property of Pt-Sn anode catalyst for electro-oxidation of ethanol, Chin. J. Catal. 25 (12) (2004) 983–988.

Please cite this article as: M.M. De Souza and A.L. De Bortoli, Stability and convergence of the numerical solution for the species equation of a model for PEMFCs, Computers and Mathematics with Applications (2018), https://doi.org/10.1016/j.camwa.2018.11.010.