Energy analysis of a proton exchange membrane fuel cell (PEMFC) with an open-ended anode using agglomerate model: A CFD study

Energy analysis of a proton exchange membrane fuel cell (PEMFC) with an open-ended anode using agglomerate model: A CFD study

Energy 188 (2019) 116090 Contents lists available at ScienceDirect Energy journal homepage: www.elsevier.com/locate/energy Energy analysis of a pro...

2MB Sizes 0 Downloads 6 Views

Energy 188 (2019) 116090

Contents lists available at ScienceDirect

Energy journal homepage: www.elsevier.com/locate/energy

Energy analysis of a proton exchange membrane fuel cell (PEMFC) with an open-ended anode using agglomerate model: A CFD study Mirollah Hosseini a, Hamid Hassanzadeh Afrouzi b, Hossein Arasteh c, Davood Toghraie d, * a

Department of Mechanical Engineering, Islamic Azad University, Qaemshahr Branch, Qaemshahr, Mazandaran, Iran Babol Noshirvani University of Technology, Babol, Iran c Department of Mechanical Engineering, Isfahan University of Technology, Isfahan, Iran d Department of Mechanical Engineering, Khomeinishahr Branch, Islamic Azad University, Khomeinishahr, Iran b

a r t i c l e i n f o

a b s t r a c t

Article history: Received 5 June 2019 Received in revised form 4 September 2019 Accepted 7 September 2019 Available online 11 September 2019

In this study, a single two-dimensional cell and open-ended anode proton exchange membrane fuel cell (PEMFC) is numerically studied using agglomerate model. The working fluids are considered water and air. The flow has been simulated using the two-phase model to consider the effects of bubble generations in the fuel cell. The numerical results show a better agreement using the agglomeration model with experimental data compared to the other methods. The effects of various parameters including, the stoichiometric coefficient, the amount of saturated water in the cathode gas diffusion layer, operating temperature and pressure, and relative humidity on the fuel cell performance, have been examined. The obtained results revealed that by increasing this coefficient from 1.5 to 2 and 2 to 2.3, the fuel cell output power enhances by 1.68% and 0.53%, respectively. It was also found that increasing the operating pressure has enhanced the mass fraction consumptions of both hydrogen and oxygen. In addition, it was deduced that the maximum local temperature occurs in the middle of the polymer fuel cell. Finally, the numerical results showed that increasing the relative humidity enhances the water formation from cathode to the anode side. © 2019 Elsevier Ltd. All rights reserved.

Keywords: Energy analysis Proton exchange membrane fuel cell Agglomerate model CFD study

1. Introduction The energy crisis and environmental pollution caused by fossil fuels have inflicted irreversible damages on the global ecosystem. Nowadays, conventional technology is one of the major sources of pollution and energy consumption due to its low efficiency and high fuel consumption [1e4]. Fuel cells are so helpful and ecofriendly according to their high efficiency and the absence of atmospheric and noise pollutants. In fact, the fuel cells convert the chemical energy into electrical energy during a reaction between hydrogen and oxygen across the electrodes and the polymer membrane with almost no sign of any pollution. Among the known fuel cells, the proton exchange membrane fuel cell (PEMFC) is more applicable due to its advantages in the mobile and transportation sections in comparison to the other low-temperature fuel cells [5e12].

* Corresponding author. Department of Mechanical Engineering, Islamic Azad University, Khomeinishahr Branch, Khomeinishahr, 84175-119, Iran. E-mail address: [email protected] (D. Toghraie). https://doi.org/10.1016/j.energy.2019.116090 0360-5442/© 2019 Elsevier Ltd. All rights reserved.

The PEMFCs with open-ended anode are studied by many researchers. There are different methods to model the catalyst layer in fuel cells. In many studies, no detail of catalyst layer is considered in the solution such as Springer et al. [13] and Bernardi et al. [14,15] in which the catalyst layer is considered as a very thin layer to cover the gas distributor layer (GDL) and a source in boundary with no mass transfer through it. Meanwhile, the presence of ions with carbon and platinum particles on the catalyst layer is highly influential on the amount of gas distribution and consequently, the current density obtained using this model is overestimated compared to the experimental data. Another method of catalyst layer modeling is the macroscopic homogeneous model [16,17]. Although in this model, the volume of the catalyst layer is considered, details of the multicomponent structure of the catalyst layer are not considered. In the cathode catalyst layer, the oxygen reduction reaction is influenced by the carbon and Nafion fibers and the platinum particles on them and also the geometric structure. These carbons are affected by this reaction and the available moisture as well as the released positive ions are mass accumulations. The size of these mass accumulations is affected by parameters like catalyst layer material type, charging rate of platinum

2

M. Hosseini et al. / Energy 188 (2019) 116090

m r

Nomenclature ðgÞ

cp

ðgÞ Di ðmÞ DH2 O

Ecell F h k L m_ H2 O Mi M(m) nd ðgÞ

specific heat capacity of gas mixture, J/kgK 2

diffusivity of species i, m /s diffusivity of water in the membrane, m2/s cell voltage (V) Faraday's constant, As/mole height, m thermal conductivity, W/mK length, m interphase mass transfer due to condensation or evaporation of water, kg/m3s molecular mass of species i, kg/mole equivalent weight of the dry membrane, kg/mole electroosmotic drag coefficient

ni

mass flux of species i, kg/m2s

pðcÞ ; pðgÞ

capillary and gas pressure, Pa liquid saturation source term temperature, K velocities, m/s coordinates, m

s S T u, u, v x, y, z

Greek symbols water content ε porosity

l

particles on carbon fibers, fibers effective area, the amount of water in the catalyst, and the penetration coefficient of oxygen in the water. Furthermore, researchers were able to model more details of the phenomena in the cathode catalyst layer using the agglomerate model [18e20]. This model was proposed for coupling of oxygen penetration and reduction reaction in the cathode catalyst layer in polymer fuel cells and considering the effects of various details such as platinum and carbon particles and ions in the catalyst layer. They considered the cylindrical geometry for these mass accumulations and assumed that the potential drop only changes in the direction of the axis of this cylinder, and also the dissolved gases penetrate the radial direction of these masses. They examined the effect of the radius of these cylinders on the current of fuel cells. One of the basic problems in their model is the absence of any complexity for the gas phase transfer in the masses. Celiker et al. [21] studied the spherical model of floating masses instead of the cylindrical model. Their method was used to model the alkaline fuel cells. Yin [22] investigated a two-dimensional and steady-state model based on the agglomeration model in the catalyst layer. They considered the mass diameters as a fixed mean value, and their results indicated that the fuel cell performance is completely affected by its structure. In most implemented models, the fixed mean value is considered for the spherical carbon masses [23]. A reasonable assumption in this model is that the mass diameter size is different from the membrane of the GDL with a constant slope as has been performed by Yoon et al. [24]. In this model, the mass diameters along the channel and at equal depths are yet considered to be equivalent. Wang et al. [25] investigated the cinematic of transportation and reaction in the cathode catalyst layer with spherical masses. They considered two types of spherical masses; the first consisting of a mixture of carbon and platinum, and proton ions, and the other type, including carbon, catalyst particles, and water in the cavities. In fact, one of the main reasons for not considering the catalyst layer structure in the polymer fuel cell

dynamic viscosity, kg/m.s density, kg/m3 phase potential, V mass fraction

f

u

Superscripts (c) capillary (g) gas phase (l) liquid phase (m) membrane ref reference (s) solid Subscripts a ave c cl eff ff gdl i in m mass pot sp

anode average cathode catalyst layer effective flow-fields gas diffusion layer species i inlet membrane mass potential separator plate

model was the absence of the proper and accurate imaging of nanoscale dimensions which was done in order to find out how the ions were transferred into the catalyst. Later, with access to advanced microscopes such as the scan electron microscope (SEM) and transport electron microscope (TEM), researchers were able to examine the complicated nanostructures of fuel cell catalyst layers. Middleman [26] examined the membrane structure with this advanced tool. He showed that the catalyst layer is a random mass of holes and particles and also showed that the particles in the catalyst layer are like masses of much smaller particles that are covered with a layer of Nafion. Sun et al. [27] studied a twodimensional model of spherical masses in the ion-impregnated environment and examined the ion volume fraction and the amount of platinum loading on carbon. The effect of liquid water inside the catalyst layer was investigated by Rao et al. [28]. The geometric properties of the masses were not considered in their work. Moreover, it was assumed that there is a steady and uniform layer of liquid water on the masses. Sasmito et al. [29] performed a two-phase model of polymer electrolyte fuel cell with a dead-end anode. They examined efficient purging strategies to manage gas and water. Their obtained results revealed that the purge duration and purging frequency play important roles in affecting the fuel cell performance [30, 31] and [32e39]. The aforementioned literature review shows that, although many researchers have investigated the performance of PEMFCs experimentally and numerically, performing a study of a twodimensional single cell and open ended-anode PEMFC with the inputs of water and air benefiting from the agglomerate model by means of computational fluid dynamics numerically is a science gap which this work aims to bridge. 2. Mathematical modeling The computational domain for a single-channel polymer

M. Hosseini et al. / Energy 188 (2019) 116090

electrolyte fuel cell consists of 9 regions for the anode and cathode which are: electricity collector, gas channel, gas diffusion layer, catalyst layer, and also the membrane which the cooling channels are not considered reducing the computational expenses. The physical model comprises a membrane (m) sandwiched by two catalyst layers (cl), two gas diffusion layers (gdl), two porous flow fields (ff), and two separator plates (sp) shown in Fig. 1 (derived from Ref. [28]). The current simulation model consists of three transport phenomena: firstly, the transfer of species, mass and momentum which two-phase species, mass and momentum conservation. Since the solubility of gasses is low, the liquid phase is assumed to contain just liquid water. Secondly, the heat transfer which considers the convection, conduction, evaporation and condensation, ohmic heating, entropy generation and irreversibilities according to the electrochemical reactions in the fuel cell. Thirdly, the charge and Ohm's law conservations equations are considered during the simulation. Tables 1 and 2 represent the geometrical and physical properties of the open-ended fuel cell, respectively.

3

Table 1 Geometrical properties of present fuel cell. Fuel cell section

Length (mm)

Channel length Channel height Gas diffusion layer height Catalyst layer height Membrane height Current collector height

90 0.5 0.3 0.01 0.05 0.5

  v  ðgÞ  _ H2 O εr þ V: rðgÞ uðgÞ ¼ Smass  m vt

(1)

 Mass conservation equation in liquid phase:

  v  ðlÞ  _ H2 O εsr þ V: rðlÞ uðlÞ ¼ m vt

(2)

 Momentum equation: 2.1. Governing equations The polymer electrolyte fuel cell characteristics are mathematically modeled by the mass conservation, momentum, energy, chemical species transfer, and electric charge transfer equations. The liquid water is also calculated by the volume fraction obtained from the saturation model. The subscripts (g), (l), (s), and (m) pertain properties are associated with the gas, liquid, solid, and membrane, respectively. Also, (c) defines any quantity related to the capillary pressure.  Mass conservation equation in gas phase:

   v  ðgÞ mðgÞ ðgÞ εr þ uðgÞ þ V: rðgÞ uðgÞ uðgÞ ¼ V:s  u vt k

(3)

 Energy equation:



rcp



v

eff vt

T þ V:



rcp

 eff

   uðgÞ T ¼ V: keff VT þ Stemp

 Species equation:

Fig. 1. A schematic diagram of present study computational domain (derived from Ref. [28]).

(4)

4

M. Hosseini et al. / Energy 188 (2019) 116090 Table 2 Physical and geometrical properties of different components of present fuel cell. Parameter

Value  A 9  109 m3  kmol 0:0564 m3 0.5 1  A 275 m2  kmol 3:39  103 m3 1 1.5 1:23ðVÞ  2 m 11:03  105 s  2 m 3:23  105 s  2 m 7:356  105 s  2 m 5 3  10 s 0.4  1 5:68  1012 m2 0.4  1 5:68  1012 m2  1 2  105 m  kg 1100 kmol

Anode reference current density Anode reference concentration Anode concentration power Anode change coefficient Cathode reference current density Cathode reference concentration Cathode concentration power Cathode change coefficient Open circuit voltage Hydrogen reference diffusion Oxygen reference diffusion Water reference diffusion Other reference diffusion Diffusion layer porosity Viscous resistance of diffusion layer Catalyst layer porosity Viscous resistance of catalyst layer Catalyst layer area to volume ratio Membrane equivalent weight Membrane surface area

0:09ðm2 Þ s 1370000 s m 500   sm 17 m W 16:3 m:K W 1:5 m:K W 0:95 m:K

Current collector electrical conductivity Electrical conductivity of gas and catalyst diffusion layer Membrane electrical conductivity Current collector thermal conductivity Thermal conductivity of gas and catalyst diffusion layer Membrane thermal conductivity

v  ðgÞ ðgÞ  ðgÞ εr ui þ V:ni ¼ Si vt

(5)

 Membrane water content:

.  v  ðmÞ ðmÞ r MH2 O l MðmÞ þ V:nH ¼ Sl 2O vt

gas, liquid and solid phases. Equation (6) only solves in the catalyst layer and membrane. To couple the water flux from equation (5), the source term of Sl is applied in the catalyst layer. The following equations represent the species of mass fluxes, current densities, liquid water velocity, and total stress tensor. ðgÞ

(6)

ni

ðmÞ

nH2 O ¼

 Membrane potential:

V:iðmÞ ¼ Spot

(7)

ðgÞ

ðgÞ

 rðgÞ Di;eff Vui

; ði ¼ H2 ; O2 ; H2 O; N2 Þ

nd MH2 O ðmÞ rðmÞ ðmÞ i  ðmÞ MH2 O DH O;eff Vl 2 F M ðmÞ

iðmÞ ¼  seff VfðmÞ ðsÞ

(8)

In above equations, the phases of liquid and gas are coupled by _ H2 O . Due to considering the the phase change source term of m porous flow field in this study, the Darcy term is also applied in the flow field in the last term of equation (3). The effective parameters are considered in energy equation in which is shared among the

( uðlÞ ¼

(10)

(12)

uðgÞ s  DðcÞ Vs ðffÞ ðgdl; clÞ DðcÞ Vs 

(9)

(11)

iðsÞ ¼  seff VfðsÞ

 Solid potential:

V:iðsÞ ¼  Spot

ðgÞ

¼ rðgÞ uðgÞ ui



s ¼  pðgÞ I þ mðgÞ VuðgÞ þ VuðgÞ

(13) T 

  2  mðgÞ VuðgÞ I 3

(14)

M. Hosseini et al. / Energy 188 (2019) 116090

In this model, all equations of components including hydrogen, oxygen, water, and nitrogen are solved in the whole domain. Equation (9), the water flux in the membrane, is represented with a phenomenological mode in terms of the membrane water content term, l which is calculated from electroosmotic drag and diffusion. Equations (10) and (11) state the electrical current from electron and hydrogen ion, respectively. According to the employed porous flow in this study, the liquid velocity in the gas phase is considered to be equal to the gas velocity with additional capillarity effect. The total stress tensor in equation (14) is due to the compressible gas flow and its pressure, viscous and volume dilation. The source terms are calculated as follows.

8 M J M H2 O J c ðmÞ O2 c > > >  4F þ 2F  V:nH2 O ðcathode clÞ > > < MH2 Ja Smass ¼ ðmÞ >  V:nH2 O ðanode clÞ  > > 2F > > : 0 ðelsewhereÞ

Si ¼

8 > > > > > > > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > > > > > > > :



MO2 Jc 4F

(15)

ðmÞ

ðsÞ vf

seff

¼ iset ;

MH2 Ja  2F

> > > :

ðH2 ; anode clÞ

ðclÞ

ðgÞ

ðgÞ

vfðsÞ ¼0; T¼Tin ; s¼0 vx (22)

 At the anode outlet:

pðgÞ ¼ pref ;

ðgÞ

vui vT vfðsÞ vs ¼ ¼ ¼ ¼0 vx vx vx vx

(23)

(17)

ðelsewhereÞ

0

(21)

in in pa ¼pref þpin a ; uH2 ¼ uH2 ;a ; uH2 O ¼ uH2 O;a ;

ðelsewhereÞ

ðmÞ

T ¼ Tcool

 At the anode inlet:

(16)

Sl ¼

vy

ðsÞ

ðgÞ

V:nH2 O

(20)

ðH2 O; gdl; ffÞ

0

8 > > > <

T ¼ Tcool

 At the lower wall:

ðH2 O; anode clÞ

_ H2 O m

2.2. Boundary conditions

fðsÞ ¼ 0;

ðH2 O; cathode clÞ

_ H2 O  V:nH2 O  m

The mass conservation source term of Smass consists of both mass production and consumption according to water transfer into the membrane and electrochemical reactions. On the other hand, the species conservation source terms of Si accounts for species production and consumption according to water transfer into the membrane and electrochemical reactions and also water interphase mass transfer. In the first row of equation (19), the first and second terms represent the reversible and irreversible entropic heat generated by the electrochemical reaction; the third and fourth terms state the Ohmic heating, and the last term indicate the interphase mass transfer.

 At the upper wall:

ðO2 ; cathode clÞ

MH2 O Jc ðmÞ _ H2 O  V:nH2 O  m þ 2F

5

 At the cathode inlet:

8 < Jc ðcathode clÞ Spot ¼ J ðanode clÞ : a 0 ðelsewhereÞ

ðgÞ ðgÞ in in _ in _ ðgÞ m a ¼ ma ; uH2 ¼ uH2 ;a ; uH2 O ¼ uH2 O;a ;

(18)

vfðsÞ ¼0; T¼Tin ; s¼0 vx (24)

 At the cathode outlet:

Stemp ¼

8 Jc > > > > > > > > > > > > > > > > > > > > < > > > > > > > > > > > > > > > > > > > > :

 T

 2  2 vErev ðmÞ ðsÞ _ H2 O Hvap ðcathode clÞ þ jhc j þ seff VfðmÞ þ seff VfðsÞ þ m vT ðmÞ

Ja ha þ seff



VfðmÞ

2

 2 ðsÞ _ H2 O Hvap ðanode clÞ þ seff VfðsÞ þ m 

ðmÞ VfðmÞ seff



VfðsÞ sðsÞ eff

2



2

_ H2 O Hvap þm

VfðsÞ sðsÞ eff

2

(19)

ðmÞ

ðff ; spÞ

ðgdlÞ

6

M. Hosseini et al. / Energy 188 (2019) 116090

Table 3 Present open ended anode fuel cell boundary conditions.

2

Parameter

Anode inlet

Cathode inlet

Mass flow rate Operating pressure Operating temperature Oxygen mass fraction Hydrogen mass fraction Water mass fraction Nitrogen mass fraction

e 111,325 Pa

0:0004 kg/s 101,325 Pa

333 K

333 K

e

0.20,218

0.3,380,353

e

0.661,964

0.1,323,057

e

0.6,655,143

kðgÞ ¼ (25)

u

¼ 0;

ðgÞ X xðgÞ a ka P ðgÞ a b xb Fab

(34)

keff ¼ εkf þ ð1  εÞkðsÞ ðsÞ

ðgÞ

vui vT vfðsÞ vfðmÞ vs ¼ ¼ ¼ ¼ ¼0 vx vx vx vx vx

(33)

  1 2 a 8 1þM Mb

The effective thermal conductivity is given by:

 At the left side and right side wall: ðgÞ

Fab ¼

The thermal conductivity of multicomponent gas mixture is calculated from:

ðgÞ

vu vfðsÞ vT vs ¼ ¼ ¼0 pðgÞ ¼ pref ; i ¼ vx vx vx vx

32  1  1 2 47 6 Mb 61 þ m a 7 mb 4 5 Ma

(26)

ðsÞ

(35) ðsÞ

ðsÞ which kðsÞ ¼ ðkff ; kgdl ; kcl ; kðsÞ m ; ksp Þ are the solid phases thermal

conductivity and the fluid thermal conductivity is defined as:

kf ¼ skðlÞ þ ð1  sÞkðgÞ

 At the interface of sp and ff:

(36)

The specific heat capacity of mixture gas is calculated from: ðgÞ vui

yðgÞ ¼ 0;

vy

¼

vs ¼0 vy

(27)

ðgÞ

cp ¼

X

ðgÞ uðgÞ cp;i i

(37)

i

The values for the boundary conditions of present open-ended anode fuel cell are shown in Table 3.

ðgÞ

ðgÞ

ðgÞ

ðgÞ

ðgÞ

which cp‫و‬i ¼ ðcp;H2 ; cp;o2 ; cp;H2 o ; cp;N2 Þ and the effective phase mixture properties is written as:



rcp

2.3. Constitutive relations



pðgÞ M RT

¼ εrf cp;f þ ð1  εÞrðsÞ cp

(38)

Each species mass diffusion coefficient based on the local pressure and temperature is determined as:

The gas density based on the ideal gas law is as follows:

rðgÞ ¼

ðsÞ

eff

ðgÞ

ðgÞ

(28)

Di



!   T 32 pðgÞ   ðgÞ ðgÞ 0 T; pðgÞ ¼ Di:0 T0 ; p0 ðgÞ T0 p

(39)

which MðgÞ accounts for the mixture molecular weight calculated as follows:

In the porous medium, a Brugeman correction is applied and according to the liquid water presence the pore blockage is considered.

. . . . 1  ðgÞ ðgÞ ðgÞ ðgÞ MðgÞ ¼ uO2 MO2 þ uH2 MH2 þ uH2 O MH2 O þ uN2 MN2

Di;eff ¼ ð1  sÞε2 Di

ðgÞ

3

ðgÞ

(40)

(29)

The relative humidity which shows the water content at the inlet of anode and cathode is defined as:

(30)



The nitrogen mass fraction is obtained from:

ðgÞ

uNðgÞ2

ðgÞ ¼ 1  uH2

ðgÞ  uO2



uðgÞ H2 O

The molecular concentrations are given by: ðgÞ

ci

¼

psat H2 O

 100

(41)

ðgÞ

which pH2 O is the water vapor partial pressure determined as:

uðgÞ rðgÞ i

(31)

Mi

The viscosity of gas mixture is calculated by: ðgÞ X xðgÞ a ma mðgÞ ¼ P ðgÞ a b xb Fab ðgÞ

pH 2 O

ðgÞ

; a; b ¼ H2 ; O2 ; H2 O; N2

ðgÞ

ðgÞ

pH2 O ¼ xH2 O pðgÞ

(42)

and psat H2 O is the water saturation pressure defined as below. In this

(32)

which xa and xb are the a and b species mole fractions, andFab is a dimensionless parameter which is defined as:

equation c1, c2, c3, and c4 are equal to 2.1794  102, 2.953  102 K1, 9.1837  105 K2, and 1.4454  107 K3, respectively as the constant magnitudes for saturation water. ref psat  10c1 þc2 ðTT0 Þþc3 ðTT0 Þ H2 O ¼ p

2

þc4 ðTT0 Þ3

The average current density is calculated from:

(43)

M. Hosseini et al. / Energy 188 (2019) 116090

iave ¼

1 L

ðL

7

Table 4 The agglomerations model data.

iðsÞ :ey dx

(44)

(45)

ðpÞ DO2

where cr is a constant accounting for the condensation or evaporation rate. The capillary diffusion for two-phase flow is defined as:

DO2

3:032  109

jref a

1  109

jref c;0

275

MðmÞ

1.1

rd aox a ; aa rd aox ; a c c bðmÞ

1,1

Units  atm:m3 mol  atm:m3 mol  atm:m3 mol V  A:s mol ðmÞ  2 m s  2 m s  A m3  A m3  kg mol e

1.5,1.5

e

0.7

e

εgdl ; εcl kgdl ; kcl

0.4,0.4

e

7:3  1013 ; 7:3  1013

ragg

107 0 0:01 ; 0:0003

ðm2 Þ ðmÞ

0

0 _ H2 O ¼ cr [email protected]ð1  sÞ m

D

Correlations and magnitudes  666 1:33 exp Tcell  498 5:08 exp Tcell  170 0:0258 exp Tcell

ðpÞ HO2

which “L” is the length of fuel cell. Concerning the condensation or evaporation of water, the interphase mass transfer is obtained as follows.

ðcÞ

Parameters

ðgÞ

pH2 O  psat H2 O RT

1 MH2 O ;  srðlÞ A

ðgdl; cl; ff Þ

ks3 dpðcÞ ¼ ðlÞ ds m

(46)

where p(c) is the capillary pressure and is given by:

 ε 12

pðcÞ ¼ t cos q

k

J

(47)

where t and q are surface tension and wetting angle, respectively and J is the Leverett function defined as:

J ¼ 1:417ð1  sÞ  2:12ð1  sÞ2 þ 1:263ð1  sÞ3

(48)

ðlÞ

HO2 ðpÞ

HH2

Erev;0 F

1.23 96487

hcl

0:01  103

 23848 4:38  106 exp RTcell

ðlÞ

q L

ðpÞ

;L

ðptÞ

uðptÞ rðmÞ ; rðptÞ

0.4

2000 ; 21450 ; 1800 ; 983

rðcÞ ; rðlÞ

ð Þ  kg m2 e  kg m2  kg m2

2.4. Phenomenological membrane model The water amount in the membrane in terms of water content of membrane is defined as:



2 3 l ¼ 0:043 þ 17:81a  39:85a þ 36:0a

14 þ 1:4ða  1Þ

a1 1
(49)

(54)

which “a” is the water activity determined as: ðgÞ pH O a ¼ sat2 pH2 O

 2s

 8 2436 > 7 > for l  3  l ½expð0:28 l Þ  1exp  3:1  10 > < T ðmÞ DH2 o ¼  > 2436 > > : 4:17  108  l½1 þ 161expðlÞexp  for l >3 T

Furthermore, the electro osmotic drag is stated as:

(50)

nd ¼ 2:5

l

(55)

22

Due to the use of a GORE membrane (is a microscopically reinforced composite membrane), for protonic conductivity, sðmÞ , and the water diffusivity of membrane, bðmÞ , a correction factor of are considered. The effective protonic conductivity is defined as: ðmÞ seff ¼ bðmÞ sðmÞ

(51)

which:



sðmÞ ¼ ð0:5193l  0:326Þexp 1268



 1 1  303:15 T

(52)

Then the effective water diffusivity of membrane is calculated by: ðmÞ 2 o;eff

DH

where:

¼b

ðmÞ

ðmÞ

D H2 o

(53)

2.5. Agglomerate model Table 4 presents the data correlations and magnitudes in agglomeration model. In Butlere-Volmer equation, in order to cover the agglomerates on both sides of anode and cathode:

0

11

2  ox  ox  ðgÞ c @ H2 A exp aa Fha  exp ac Fha Ja ¼ jref ð1  sÞ a ðgÞ RT RT cH ;ref 2

(56)

8

M. Hosseini et al. / Energy 188 (2019) 116090

0

1" ! !# ard ard ref a F c F @ A Jc ¼ jc ð1  sÞ ðgÞ h þ exp h  exp RT c RT c cO ;ref 2 ! gðpÞ RT 1 1  ðaggÞ x  ð1  gcl Þ ðpÞ 1 1 þ x þ x g 2 3 H ðgÞ

cO2

O2

uðPÞ ¼

LðCÞ ¼

(57) In above equations,

jref a=c

are the volumetric reference exchange

current densities. In addition,

aox=rd a=c

ðPtCÞ

g

are the transfer coefficients for

anode/cathode oxidation/reduction reactions, which it has been assumed that the transfer coefficient at oxidation and reduction reactions are equal for both cathode and anode. Moreover, ðgÞ are the hydrogen and oxygen reference concentrations. 2 =O2 ; ref ðpÞ Also, HO2 represents the Henry's constant for the air-polymer

cH

interface. The correction factors of x1 , x2 , and x3 are due to the resistances of the agglomerate, films of polymer and liquid water, respectively. Further, gcl , gðaggÞ , and gðpÞ are the pores volume fraction in layer of catalyst, polymer volume fraction, and agglomerate volume fraction, respectively. A correction for temperature of jref c by Arrhenius-type and also the volume concentrations of components are presented as:

   Ea 1 1 ref  jref ¼ j exp  c c;0 R T T1

LðpÞ LðPtÞ þ LðCÞ þ LðpÞ LðPtÞ

uðPtÞ

¼

pref

x1 ¼

1



1 1  F tanhð3FÞ 3F

(60)

hc ¼ fðsÞ  fðmÞ  Erev

(61)

where Erev shows the reversible potential and is defined as:

Erev ¼ Erev;0 þ [1 ðT  T2 Þ þ

RT ðgÞ lnxO2 4F

(62)

where the subscript of “000 denotes the standard conditions as [1 and T2 are constants. The components volume fractions and the used parameters are presented and calculated below by using the data in Table 4:

gcl ¼

Vvoid ¼ 1  gðaggÞ ¼ 0:475 Vtot

 (71)

(72)

! jref c

!

ðpÞ 1  ggðaggÞ

ðgÞ 2 ;ref

4FcO

ðpÞ ¼ DO2

exp 

gðpÞ gðaggÞ

!

ard a F RT

hc þ exp

ard c F RT hc

(73)

ðgÞ 2 ;ref

4FcO

!1:5 (74) ðaggÞ , 2 ;eff

In above equations, F, rðaggÞ , kc , DO

ðpÞ

and DO2 are Thiele

modulus, agglomerate radius, reaction rate, oxygen effective diffusion coefficient in polymer inside the agglomerate (using Bruggeman correlation), and oxygen diffusion coefficient in polymer film, respectively. Also, the correction factor of x2 and the used parameters are given by:

x2 ¼

ðpÞ

(63)

(70)

O2 ;eff

ðaggÞ DO ;eff 2

ha ¼ fðsÞ  fðmÞ

(69)

vffiffiffiffiffiffiffiffiffiffiffiffiffiffi u rðaggÞ u kc t F¼ ðaggÞ 3 D

(58)

The anode and cathode over potentials are determined as:

 LðPtÞ ¼ 0:45  103 kg m2

In above equations, uðpÞ , uðPtÞ , LðPtÞ , LðPÞ , LðcÞ , rðPtÞ and Vtot are the mass fraction of polymer loading in the layer of catalyst and platinum loading on carbon, platinum, polymer, and carbon loading, total volume in the layer of catalyst, and platinum density, respectively. The correction factor of x1 and the used parameters are determined as:

(59)

ðpÞ

Hi

(68)

" # VðPtCÞ 1 1  uðPtÞ LðPtÞ ¼ ¼ ðPtÞ þ ðCÞ ðPtÞ ¼ 0:0264 Vtot hcl r r u

kc ¼ ðgÞ ci;ref

¼ 0:93

d

dðpÞ x1 ðpÞ ðpÞ DO2 a

kc

(75)

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi ! u u 3 gðpÞ 3  t ðaggÞ ¼ r 1 þ ðPtCÞ  rðaggÞ ¼ 1:7082344  107 m

g

(76)

gðpÞ ¼

ðpÞ

uðpÞ

ðPtÞ

V 1 L ¼ ¼ 0:498 Vtot 1  uðpÞ rðmÞ hcl uðPtÞ

(64)

  ðpÞ 2 ¼ 5808950:658 m2 aðpÞ ¼ 4pnðaggÞ rðaggÞ þ d

(77)

VðaggÞ ¼ gðpÞ þ gðPtCÞ ¼ 0:5244 Vtot

(65)

nðaggÞ ¼

VðaggÞ ¼ VðPtCÞ þ VðpÞ ; Vvoid ¼ VðgÞ þ VðlÞ

(66)

In above equations, dðpÞ , aðpÞ , and nðaggÞ indicate the polymer thickness, catalyst layer agglomerate surface area per unit volume, and number of agglomerates per unit volume, respectively. Finally, the correction factor of x3 and the used parameters are presented below:

gðaggÞ ¼

Vtot ¼ VðaggÞ þ Vvoid ¼ VðPtCÞ þ VðpÞ þ VðgÞ þ VðlÞ

(67)

3gðaggÞ 1 ¼ 6:302535827  1018 3  ðpÞ 3 ðaggÞ m 4p r þd

(78)

M. Hosseini et al. / Energy 188 (2019) 116090

x3 ¼

ðlÞ

d

dðlÞ x1 ðlÞ

DO2

ðlÞ

HO kc ðpÞ2 ðlÞ a H

(79)

O2

vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi !ffi u   ðlÞ u  g 3  3 ðpÞ ¼ t rðaggÞ þ d 1 þ ðaggÞ  rðaggÞ þ dðpÞ

g

  ðpÞ ðlÞ 2 aðlÞ ¼ 4pnðaggÞ rðaggÞ þ d þ d

gðlÞ ¼

9

(80)

(81)

VðlÞ ¼ sgcl Vtot

(82) ðlÞ

ðlÞ

In above equations, dðlÞ , DO2 , aðlÞ , HO2 , gðlÞ , and s are liquid water film thickness, oxygen diffusion coefficient in the liquid water, the agglomerate surface area comprising liquid water per unit volume, Henry's constant for the air-water interface, liquid water volume fraction, and liquid saturation function, respectively.

Fig. 3. Grid independency analysis.

3. Numerical procedure 3.1. Grid study A structural and non-uniform grid is employed to mesh the computational domain of the present open-ended anode polymer membrane fuel cell displayed in Fig. 2. To evaluate the independency of grid size from the solution, the effect of different grid sizes (30,000e50,000- 75,000e104,000- 158,000e300,000) on the fuel cell current at a voltage of 0.61 (Volt) have been examined and indicated in Fig. 3. This figure shows that the relative error for the fuel cell current after the number of meshes equal to 158,000 becomes very low which proves that this number of meshes is appropriate to be employed for further numerical simulations. 3.2. Validation To validate the present two-phase numerical study of openended anode fuel cell, the obtained results for the voltage-current figure at both agglomerate and non-agglomerate models have been compared to the experimental data of Rao et al. [28]. This comparison has been displayed in Fig. 4 which demonstrates a good agreement between the present numerical study and experimental data. This figure shows a little difference between the results at low voltages and high currents. In fact, in this region, due to generating the bubbles in the catalyst and GDL, reactions can pass it hardly and

Fig. 4. Comparison of present numerical data for agglomerate and non-agglomerate models of fuel cell with experimental data.

Fig. 5. Variations of voltage with current at temperature of 333 K and different stoichiometric coefficients. Fig. 2. The present study generated mesh of open ended anode polymer membrane fuel cell.

10

M. Hosseini et al. / Energy 188 (2019) 116090

Fig. 6. Saturated water contour in the cathode gas diffusion layer. Fig. 8. Mass fraction contour of oxygen in the present fuel cell computational domain.

causes a reduction in the fuel cell performance. Since the present study is performed using the two-phase model and considering the effects of bubbles generation in this region, lower performances have been resulted.

4. Results and discussion 4.1. Stoichiometric coefficient effect on the fuel cell performance Fig. 5 indicates the variations of voltage with current at temperature of 333 K and different stoichiometric coefficients. It reveals that, by increasing this coefficient from 1.5 to 2 and 2 to 2.3, the fuel cell output power enhances by 1.68% and 0.53%, respectively. The increase of the stoichiometric coefficient is associated with a higher waste of fuel and reduces the fuel cell overall efficiency. In order to prevent fuel waste, the use of a circulating pump is necessary to conduct the unused fuel to the channel inlet, which requires the use of a part of the energy generated by the fuel cell and it may not be worth it.

4.2. Effect of saturated water in the cathode gas diffusion layer on the fuel cell performance Fig. 6 displays the contour of saturated water in the cathode gas diffusion layer of open-ended anode polymer membrane fuel cell. At this stage, the other parameters including stoichiometric coefficient and relative humidity of hydrogen and nitrogen are fixed to 2.3 and 100%, respectively. In fact, the role of saturated water is to block the gaseous penetration layers and reduce the ability of oxygen transfer from the channel to the catalyst level, while reducing moisture can cause the membrane to dry and severely inhibit ion penetration. Furthermore, at higher pressures, the amount of water leaving the gas diffusion layer surface, especially near the channel entrance, increases which leads to a reduction in the saturated water in this area. Although changing the operating pressure of the fuel cell results in a difference in the input mass fractions of oxygen and hydrogen due to fixing other parameters, the fact that how oxygen and hydrogen can be consumed along the channel, as well as the slope of the changes, can be important. Therefore, the amount of oxygen and hydrogen mass fraction consumption on the outer surface of the gas diffusion layer at different operating pressures is depicted in Fig. 7. This figure demonstrates that the

Fig. 7. Mass fraction consumption of (a) hydrogen and (b) oxygen along the channel at different operating pressures.

M. Hosseini et al. / Energy 188 (2019) 116090

mass fraction consumption of oxygen is higher than that of the hydrogen at both pressures of 1 and 2 bar. Also, the slope of mass fraction consumption for hydrogen and oxygen at a pressure of 1 bar is almost the same, but it differ at a higher pressure of 2.1 bar. In addition, increasing the operating pressure has enhanced both mass fraction consumptions which means at higher pressures, components penetration into the gas diffusion layer, as well as the consumption rate at these regions, increase significantly. Fig. 8 presents the mass fraction contour of oxygen at different fuel cell computational domains. It shows that on the cathode channel path, the oxygen gradient is uniform. However, in order to obtain the maximum efficiency of the fuel cell, strategies should be considered to cause a uniform distribution of species throughout the fuel cell, in a way that the percentage of the total surface of membrane participation in the fuel cell performance be equivalent. One of the measures that can be taken for this purpose is to change the cross-sectional area of the channel along the flow path so that the cross-sectional area of the channel decreases by the flow direction to alleviate the effects of species concentration reduction, to some extent.

Fig. 10. Hydrogen magnitudes.

11

consumption

distribution

for

different

relative

humidity

4.3. Operating temperature effect on the fuel cell performance Increasing the temperature in the fuel cell reduces the density of the existing water vapor of the penetration layers, which allows the components to penetrate more easily into these layers and thereby increase the performance of the fuel cell. However, the excessive increase in the operating temperature causes the membrane to dry and consequently reduces the ion transfer coefficient. In fact, studies have shown that the operating temperature range for polymer fuel cells is 30e80 Celsius. On the other hand, it should be taken into account that raising the temperature of its input components requires energy consumption, which requires a balance between these two parameters so that the overall efficiency of the fuel cell reaches its maximum. Fig. 9 indicates the temperature contour of the fuel cell and shows that the maximum local temperature occurs in the middle of the polymer fuel cell, which is due to the formation of saturated water at the end of the gas diffusion layer.

Fig. 11. Water formation distribution for different relative humidity magnitudes.

Fig. 12. Nitrogen consumption distribution for different relative humidity magnitudes.

Fig. 9. Temperature contour in the present fuel cell computational domain.

12

M. Hosseini et al. / Energy 188 (2019) 116090

4.4. Relative humidity effect on the fuel cell performance Evaluating the relative humidity for hydrogen concentration is one of the key parameters in water management of fuel cells that shown in Fig. 10. In this regard, with increasing the relative humidity, the hydrogen concentration reduces and water formation on the cathode side decreases as well which causes an enhancement in concentration gradient water on both sides of the membrane which leads to water penetration from cathode side to the anode channel. Fig. 11 depicts the water formation distribution in the openended anode fuel cell at various relative humidity magnitudes. It shows that increasing the relative humidity enhances the water formation from cathode to the anode side. However, by augmenting the relative humidity, the hydrogen concentration has been reduced due to the water concentration enhancement displayed in Fig. 12. 5. Conclusions In this study, a two-dimensional single cell and open endedanode PEMFC with the inputs of water and air employing the agglomerate model was numerically analyzed. The flow was simulated with the two-phase model to consider the effects of bubbles generations in the fuel cell. The numerical results showed a better agreement with experimental data upon using the agglomeration model. The effects of various parameters including the stoichiometric coefficient, the amount of saturated water in the cathode gas diffusion layer, operating temperature and pressure, and relative humidity on the fuel cell performance have been examined. The obtained results are summarized as follows. 1. By increasing stoichiometric coefficient from 1.5 to 2 and 2 to 2.3, the fuel cell output power enhances by 1.68% and 0.53%, respectively. The increase of the stoichiometric coefficient is associated with a higher waste of fuel and reduces the fuel cell overall efficiency 2. Increasing the operating pressure has enhanced both mass fraction consumptions of hydrogen and oxygen which means at higher pressures, components penetration into the gas diffusion layer increases significantly. 3. Increasing the temperature in the fuel cell reduces the density of the existing water vapor of the penetration layers, which allows the components to penetrate more easily into these layers and thereby increase the performance of the fuel cell. However, the excessive increase in the operating temperature causes the membrane to dry and consequently reduces the ion transfer coefficient 4. The maximum local temperature occurs in the middle of the polymer fuel cell, which is due to the formation of saturated water at the end of the gas diffusion layer. 5. By increasing the relative humidity, the hydrogen concentration reduces and water formation on the cathode side decreases as well which causes an enhancement in concentration gradient water on both sides of the membrane which leads to water penetration from the cathode side to the anode channel. 6. Increasing the relative humidity enhances the water formation from cathode to the anode side. However, by augmenting the relative humidity, the hydrogen concentration has been reduced due to the water concentration enhancement.

Conflicts of interest statement The authors whose names are listed immediately below certify

that they have NO affiliations with or involvement in any organization or entity with any financial interest (such as honoraria; educational grants; participation in speakers’ bureaus; membership, employment, consultancies, stock ownership, or other equity interest; and expert testimony or patent-licensing arrangements), or non-financial interest (such as personal or professional relationships, affiliations, knowledge or beliefs) in the subject matter or materials discussed in this manuscript. Also, the authors whose names are listed immediately below report the following details of affiliation or involvement in an organization or entity with a financial or non-financial interest in the subject matter or materials discussed in this manuscript. Please specify the nature of the conflict on a separate sheet of paper if the space below is inadequate. References [1] Ijaodola OS, et al. Energy efficiency improvements by investigating the water flooding management on proton exchange membrane fuel cell (PEMFC). Energy 2019;179:246e67. [2] Lim BH, Majlan EH, Daud WRW, Rosli MI, Husaini T. Three-dimensional study of stack on the performance of the proton exchange membrane fuel cell. Energy 2019;169:338e43. [3] Arasteh H, Mashayekhi R, Toghraie D, Karimipour A, Bahiraei M, Rahbari A. Optimal arrangements of a heat sink partially filled with multilayered porous media employing hybrid nanofluid. J Therm Anal Calorim 2019:1e14. [4] Toghraie D, Mashayekhi R, Arasteh H, Sheykhi S, Niknejadi M, Chamkha AJ. Two-phase investigation of water-Al2O3 nanofluid in a micro concentric annulus under non-uniform heat flux boundary conditions. Int J Numer Methods Heat Fluid Flow 2019. https://doi.org/10.1108/HFF-11-2018-0628. [5] Afra M, Nazari M, Kayhani MH, Sharifpur M, Meyer JP. 3D experimental visualization of water flooding in proton exchange membrane fuel cells. Energy 2019;175:967e77. [6] Dashti I, Asghari S, Goudarzi M, Meyer Q, Mehrabani-Zeinabad A, Brett DJL. Optimization of the performance, operation conditions and purge rate for a dead-ended anode proton exchange membrane fuel cell using an analytical model. Energy 2019;179:173e85. [7] Lai X, Long R, Liu Z, Liu W. A hybrid system using direct contact membrane distillation for water production to harvest waste heat from the proton exchange membrane fuel cell. Energy 2018;147:578e86. [8] Laoun B, Kasat HA, Ahmad R, Kannan AM. Gas diffusion layer development using design of experiments for the optimization of a proton exchange membrane fuel cell performance. Energy 2018;151:689e95. [9] Kong IM, Jung A, Kim YS, Kim MS. Numerical investigation on double gas diffusion backing layer functionalized on water removal in a proton exchange membrane fuel cell. Energy 2017;120:478e87. n-Gaste lum M, et al. Pt-Au nanoparticles on graphene for oxygen [10] Beltra reduction reaction: stability and performance on proton exchange membrane fuel cell. Energy 2019;181:1225e34. [11] Arasteh H, Mashayekhi R, Goodarzi M, Motaharpour SH, Dahari M, Toghraie D. Heat and fluid flow analysis of metal foam embedded in a double-layered sinusoidal heat sink under local thermal non-equilibrium condition using nanofluid. J Therm Anal Calorim 2019:1e16. [12] Arasteh H, Salimpour MR, Tavakoli MR. Optimal distribution of metal foam inserts in a double-pipe heat exchanger. Int J Numer Methods Heat Fluid Flow 2019;29(4):1322e3142. [13] Springer TE, Zawodzinski TA, Gottesfeld S. Polymer electrolyte fuel cell model. J Electrochem Soc 1991;138(8):2334e42. [14] Bernardi DM, Verbrugge MW. “A mathematical model of the solid-polymerelectrolyte fuel cell. J Electrochem Soc 1992;139(9):2477e91. [15] Bernardi DM, Verbrugge MW. Mathematical model of a gas diffusion electrode bonded to a polymer electrolyte. AIChE J 1991;37(8):1151e63. [16] Dhar HP. Medium-term stability testing of proton exchange membrane fuel cell stacks as independent power units. J Power Sources 2005;143(1e2): 185e90. [17] Song D, et al. A method for optimizing distributions of Nafion and Pt in cathode catalyst layers of PEM fuel cells. Electrochim Acta 2005;50(16e17): 3347e58. [18] Machado BS, Mamlouk M, Chakraborty N. Three-dimensional agglomerate model of an anion exchange membrane fuel cell using air at the cathodeeA parametric study. J Power Sources 2019;412:105e17. [19] Machado BS, Chakraborty N, Mamlouk M, Das PK. A three-dimensional agglomerate model of an anion exchange membrane fuel cell. J. Electrochem. Energy Convers. Storage 2018;15(1):11004. [20] Xing L. An agglomerate model for PEM fuel cells operated with non-precious carbon-based ORR catalysts. Chem Eng Sci 2018;179:198e213. [21] Celiker H, Al-Saleh MA, Gultekin S, Al-Zakri AS. A mathematical model for the performance of Raney metal gas diffusion electrodes. J Electrochem Soc 1991;138(6):1671e81. [22] Yin K-M. Parametric study of proton-exchange-membrane fuel cell cathode using an agglomerate model. J Electrochem Soc 2005;152(3):A583e93.

M. Hosseini et al. / Energy 188 (2019) 116090 [23] Xia Z, Wang Q, Eikerling M, Liu Z. Effectiveness factor of Pt utilization in cathode catalyst layer of polymer electrolyte fuel cells. Can J Chem 2008;86(7):657e67. [24] Yoon W, Weber AZ. Modeling low-platinum-loading effects in fuel-cell catalyst layers. J Electrochem Soc 2011;158(8):B1007e18. [25] Wang Q, Eikerling M, Song D, Liu Z. Structure and performance of different types of agglomerates in cathode catalyst layers of PEM fuel cells. J Electroanal Chem 2004;573(1):61e9. [26] Middelman E. Improved PEM fuel cell electrodes by controlled self-assembly. Fuel Cells Bull 2002;2002(11):9e12. [27] Sun W, Peppley BA, Karan K. An improved two-dimensional agglomerate cathode model to study the influence of catalyst layer structural parameters. Electrochim Acta 2005;50(16e17):3359e74. [28] Rao RM, Bhattacharyya D, Rengaswamy R, Choudhury SR. A two-dimensional steady state model including the effect of liquid water for a PEM fuel cell cathode. J Power Sources 2007;173(1):375e93. [29] Sasmito AP, Mujumdar AS. Performance evaluation of a polymer electrolyte fuel cell with a dead-end anode: a computational fluid dynamic study. Int J Hydrogen Energy 2011;36(17):10917e33. [30] Boushehri R, Hasanpour Estahbanati S, Ghasemi-Fare O. Controlling frost heaving in ballast railway tracks using low enthalpy geothermal energy. 2019. [31] Ramezanianpour AA, Zolfagharnasab A, Zadeh FB, Estahbanati SH, Boushehri R, Pourebrahimi MR, Ramezanianpour AM. Effect of supplementary cementing materials on concrete resistance against sulfuric acid attack. In: InHigh tech concrete: where technology and engineering meet. Cham: Springer; 2018. p. 2290e8. [32] Afshari A, Akbari M, Toghraie D, Eftekhari Yazdi M. Experimental investigation

[33]

[34]

[35]

[36]

[37]

[38]

[39]

13

of rheological behavior of the hybrid nanofluid of MWCNTealumina/water (80%)eethylene-glycol (20%). J Therm Anal Calorim 2018;132:1001e15. Saeedi AH, Akbari M, Toghraie D. An experimental study on rheological behavior of a nanofluid containing oxide nanoparticle and proposing a new correlation. Phys E 2018;99:285e93. Esfahani NN, Toghraie D, Afrand M. A new correlation for predicting the thermal conductivity of ZnOeAg (50%e50%)/water hybrid nanofluid: An experimental study. Powder Tech 2018;323:367e73. Arabpour A, Karimipour A, Toghraie D. The study of heat transfer and laminar flow of kerosene/multi-walled carbon nanotubes (MWCNTs) nanofluid in the microchannel heat sink with slip boundary condition. J Therm Anal Calorim 2018;131:1553e66. Toghraie D, Abdollah MMD, Pourfattah F, Akbari OA. B Ruhani Numerical investigation of flow and heat transfer characteristics in smooth, sinusoidal and zigzag-shaped microchannel with and without nanofluid. J Therm Anal Calorim 2018;131:1757e66. Ruhani B, Barnoon P, Toghraie D. Statistical investigation for developing a new model for rheological behavior of Silicaeethylene glycol/Water hybrid Newtonian nanofluid using experimental data. Phys Stat Mech Appl 2019;525: 616e27. Barnoon P, Toghraie D, Dehkordi RB, Abed H. MHD mixed convection and entropy generation in a lid-driven cavity with rotating cylinders filled by a nanofluid using two phase mixture model. J Magn Magn Mater 2019;483: 224e48. Sarlak R, Yousefzadeh S, Akbari OA, Toghraie D, Sarlak S. The investigation of simultaneous heat transfer of water/Al2O3 nanofluid in a close enclosure by applying homogeneous magnetic field,. Int J Mech Sci 2017;133:674e88.