- Email: [email protected]
Pore-scale investigation of catalyst layer ingredient and structure effect in proton exchange membrane fuel cell
Applied Energy 253 (2019) 113561
Contents lists available at ScienceDirect
Applied Energy journal homepage: www.elsevier.com/locate/apenergy
Pore-scale investigation of catalyst layer ingredient and structure effect in proton exchange membrane fuel cell ⁎
Yuze Houa, Hao Denga, Fengwen Panb, Wenmiao Chenb, Qing Dua, , Kui Jiaoa, a b
T
⁎
State Key Laboratory of Engines, Tianjin University, 135 Yaguan Road, Tianjin 300350, China Weichai Power Co. Ltd., 197A Fushou St. E., Weifang 261000, China
H I GH L IG H T S
layers under various ingredient contents and structures are reconstructed. • Catalyst pore-scale model is developed to simulate the reactive transport processes. • AHigh platinum/catalyst ratio can enhance the transport and improve the performance. • The active catalyst area should be given precedence during the electrode fabrication. • A novel catalyst layer design is proposed and can improve the performance by 50%. •
A R T I C LE I N FO
A B S T R A C T
Keywords: Proton exchange membrane fuel cell (PEMFC) Catalyst layer Lattice Boltzmann method (LBM) Platinum/catalyst (Pt/C) ratio Ionomer/catalyst (I/C) ratio Electrode structure design
A pore-scale model based on the lattice Boltzmann method (LBM) is developed to simulate the reactive transport processes in the cathode catalyst layer of a proton exchange membrane fuel cell (PEMFC). The porous structures of the cathode catalyst layers are reconstructed in the process-based method with the consideration of carbon supporter, platinum, ionomer and pores. Its characteristics are analyzed including pore size distribution, phase connectivity and active catalyst area. The effects of two critical parameters, platinum/catalyst (Pt/C) and ionomer/catalyst (I/C) ratios, and structure design are investigated in terms of oxygen concentration distribution, reactive area, and reaction rate. The results indicate that, under the constant platinum loading (0.3 mg cm−2), a higher Pt/C ratio yields a thinner catalyst layer, which significantly enhances the oxygen transport and improves the performance. For the same Pt/C ratio, although a higher I/C ratio brings more mass transport loss, it increases the active catalyst area and ultimately yields better performance. Therefore, the active catalyst area should be given precedence during catalyst layer fabrication. To realize a large active catalyst area on the premise of low transport loss, an ideal catalyst layer structure design is proposed and capable of improving the performance by 50%.
1. Introduction Proton exchange membrane fuel cell (PEMFC) is considered as one of the most promising energy conversion devices for its attractive features such as high efficiency and fast response [1,2]. Among the components of a PEMFC, the catalyst layer (CL) sandwiched between the porous transport layer (PTL) and the PEM is the most complex and important since it is the place where energy conversion occurs [3]. Specifically, the CL is required to be heterogeneous for the transport of protons, electrons, reactants and products, and also provide an active catalyst area for electrochemical reactions [4]. To meet these needs, the state-of-the-art CL is porous in structure and commonly composed of
⁎
carbon-supported platinum-based catalyst (platinum particles supported on carbon, as shown in Fig. 1(a)), ionomer and pores [5]. Additionally, to achieve more efficient transport, larger active catalyst area, and ultimately better cell performance, the content of each ingredient and the structure should be optimized delicately [6]. Experimentally, Lee and Hwang [7] studied the effects of ionomer loading and distribution and found distribution schemes can influence the cell performance. Ozden et al. [8] investigated the effects of ionomer type and observed that both the loading and type of ionomer affect the degradation pattern and performance of PEMFC. Shahgaldi et al. [9] studied the electrode manufacturing processes and found the manufacturing strategy should vary with the Pt loading to reach a
Corresponding authors. E-mail addresses: [email protected] (Q. Du), [email protected] (K. Jiao).
https://doi.org/10.1016/j.apenergy.2019.113561 Received 22 April 2019; Received in revised form 14 July 2019; Accepted 15 July 2019 0306-2619/ © 2019 Published by Elsevier Ltd.
Applied Energy 253 (2019) 113561
Y. Hou, et al.
Nomenclature A C D d e F f H i i0 l n P R S T t x y z
Greek letters
area of the cross-section, m2 oxygen concentration, mol m−3 diffusivity, m2 s−1 local pore size, m discrete lattice velocity Faraday’s constant, C mol−1 discrete distribution function Henry’s constant local current density, A m−2; discrete direction exchange current density, A m−2 thickness of catalyst layer, nm neighboring active Pt node number pressure, atm universal gas constant; radius of hole structure, nm reaction source term temperature, K lattice time direction; lattice site direction direction
γPt θPt θI ε εPt εI εcarbon ρ α η τ Δ
platinum loading, mg cm−2 platinum/catalyst ratio ionomer/catalyst ratio porosity volume fraction of platinum volume fraction of ionomer volume fraction of carbon density, mg cm−3 cathode transfer coefficient overpotential, V relaxation time lattice resolution, m
Subscripts and superscripts b eq I kn p ref
bulk equilibrium ionomer Knudsen pore reference
Fig. 1. (a) TEM image of the catalyst layer; (b) Reconstructed microstructures of the catalyst layer (2670 × 450 × 450 nm); (c) Catalyst layer reconstruction processes; (d) Direction instruction for reconstruction and D3Q19 model.
capable of modeling the whole cell and predicting its performance. However, to achieve this, resolving the CL microstructure is neglected as it requires too much computation and empirical formulations are required to determine some critical transport parameters. Therefore, the macroscale model is not able to consider the effect of local structure in the CL. As for the microscopic model, the molecular dynamics (MD) model [17,18] is capable of considering the swelling of ionomer and diffusivity at the molecular scale. While aiming for the balance of ingredient and structure optimization, the computation domain should be in micrometer size, which is far beyond the scope of MD. In contrast, the mesoscopic pore-scale model, which can model the detailed processes in the CL with acceptable computational cost, is the most appropriate method for the mentioned targets and also the trend of CL modeling [19]. Wang et al. [20,21] simulated the transport of oxygen and proton in simplified CLs with the finite volume method. Lange et al. [22,23] further considered the local Knudsen diffusion effects and
better performance. They also designed a cathode catalyst layer (CCL) with two sub-layers filled by different ionomer loadings and concluded that an appropriate loading gradient could be beneficial for both transport and performance [10]. Zhao et al. [11] found that under constant Pt loading, the catalyst with low platinum content results in a thicker electrode. In their further study [12], the gas permeability in the catalyst electrode was tested, and it was found that the presence of a CL significantly reduced the gas permeability. However, it is too time consuming for an experiment to be conducted on every combination of variables. Due to the nanometer scale, it is still impractical for advanced equipment to realize in-situ observation of transport processes inside the CL and also difficult to control its structure and variables accurately at will [4]. Therefore, numerical simulation is considered a complementary approach to obtain a fundamental understanding of the underlying mechanisms in the CL [13]. Through the existing studies, the macroscale model [14–16] is 2
Applied Energy 253 (2019) 113561
Y. Hou, et al.
conducted to check the effects on the reconstructed structures, such as pore size distribution (PSD), ionomer connectivity, and active catalyst area. The specific implementation algorithms of reconstruction and analyses are shown below.
Table 1 Parameter values in this work. Parameters
Symbol
Value
Units
platinum loading platinum density ionomer density carbon density overpotential Temperature Pressure Henry’s constant exchange current density reference oxygen concentration Cathode transfer coefficient Oxygen concentration at TPL/CL interface Lattice resolution
γPt ρPt ρI ρcarbon η T P H i0 Cref α Cinlet Δ
0.3 21,450 2000 1800 0.3 353 1.5 6.955 0.015 40.96 0.616 10.875 3
mg cm−2 mg cm−3 mg cm−3 mg cm−3 V K atm – A m−2 mol m−3 – mol m−3 nm
2.1. Calculate the ingredient contents In this work, the constant Pt loading (γPt ) is 0.3 mg cm−2. According to Ref. [11], a lower Pt/C ratio (θPt ) results in a thicker electrode. And it is reported in Ref. [10] that the I/C ratio (θ I ) and porosity are in inverse proportion since the ionomer occupies the space of the pore and the I/C ratio between 1:3 and 1:2.3 is expected to yield better performance and durability. Therefore, it is rational to make the following assumptions: 1. Under the reference I/C ratio (1:3), the porosity (ε ) of CL is constant (35%), and the effect of Pt/C ratio on the structure is considered by changing CL thickness. 2. Under the same Pt/C ratio, the thickness (l) of CL is not affected by I/C ratio.
reconstructed the CL with a resolution of 2 nm. The CL was represented as carbon spheres covered by thin ionomer with a specified thickness. Chen et al. [24] firstly modeled the reactive transport processes with the lattice Boltzmann method (LBM) and used quartet structure generation set (QSGS) to reconstruct the four-constituent CL microstructure, which was considered a great advancement in model development. And in their later work, the dissolution resistance is coupled into the model, and the schemes for reducing transport resistance are proposed [19]. Salomov et al. [25] conducted the pore scale simulations to predict the permeability in the reconstructed GDL and CL. Through a brief review, it can be demonstrated that a pore-scale model is a powerful tool for gaining valuable insights into the CL. With such a potential model, many unexplored but meaningful works can be realized. For CL fabrication, the contents of ingredients are described by two parameters, platinum/catalyst ratio (Pt/C, note that in this paper C stands for catalyst whose weight is the sum of platinum and carbon) and ionomer/catalyst (I/C) ratio. And these two critical parameters are usually set by experience or repeated tests, which can be inappropriate and time-consuming [26]. Besides, the structured CL is promising with the capacity to improve the cell performance significantly, but difficult to fabricate [6]. To the best of our knowledge, there has been no pore-scale based study to work on the above subjects. In this work, CLs under various parameters and structure designs are reconstructed elaborately with the consideration of carbon, Pt, ionomer, and pores, as shown in Fig. 1(b). And the LBM based pore-scale model is developed to simulate the reactive transport processes in the reconstructed CLs to offer valuable advice for fabricating CL and propose novel structure designs. Due to the relatively low oxygen reduction reaction (ORR) rate, this work focuses on the cathode side. The large-scale simulations are realized by in-house code computed in parallel using graphics processing unit (GPU) workstation. In the following, the algorithms of reconstruction and characteristic analyses are presented in Section 2, and the LBM based pore-scale model is developed in Section 3. The results and discussion are given in Section 4, while the conclusions are drawn in Section 5.
Regarding the reconstructed CL size, both the width and height are 450 nm (y and z-direction in Fig. 1(b)), and the thickness is determined by the Pt/C ratio. The volume fraction of each ingredient (εPt , εI and εcarbon ) and the CL thickness can be figured out by:
εPt = (1 − ε )
εI = (1 − ε )
θPt / ρPt θ I/ ρI + (1 − θPt )/ ρcarbon + θPt / ρPt
θ I / ρI θ I/ ρI + (1 − θPt )/ ρcarbon + θPt / ρPt
εcarbon = (1 − ε )
l=
(1 − θPt )/ ρcarbon θ I/ ρI + (1 − θPt )/ ρcarbon + θPt / ρPt
γPt εPt ρPt
(1)
(2)
(3)
(4)
where ρPt , ρI and ρcarbon are the density of Pt, ionomer, and carbon, respectively. In this paper, the specific values of the parameters are given in Table 1. 2.2. Generate carbon support Various kinds of carbon, such as Ketjen Black and Vulcan XC-72, have distinctly different structural properties. Since such a factor is not focused on the present study, the morphology of carbon is simplified as a sphere, which is widely adopted [19,22,23]. Carbon spheres with random radii from 8 to 22 nm [10] are generated one by one until the required fraction εcarbon is reached. The specific algorithm inspired by Lange [22] is briefly explained below. 1. A group of random functions gives the center coordinate and radius of the first sphere. Then the first carbon sphere is generated based on the given coordinate and radius, and the corresponding pore nodes are converted to carbon nodes. 2. Another combination of radius and center coordinate is given randomly. And a function is set to check if the given combination meets the following two requirements: a. The newly generated carbon sphere should overlap with at least one existing carbon sphere. b. And the overlap rate cannot exceed the specified value (15% in this work).
2. CL reconstruction and analyses algorithms Microstructure reconstruction is the most basic element of the porescale CL model and many researchers have put much effort to build the microstructure more realistically with higher resolution [20,27,28]. Among the various reconstruction approaches, the process-based method proposed in Ref. [27] is widely adopted for its feasibility and rationality. Following the fabrication processes, the reconstruction algorithms can be divided into three steps: 1. generate carbon supporter, 2. generate Pt particles and 3. generate ionomer, as shown in Fig. 1(c). Weighing the computation domain size against the computation resource, the grid resolution is set as 3 nm per lattice cell, which is enough to consider the presence of Pt. The characteristic analyses are also
If yes, the next carbon sphere is generated based on the examined combination. Otherwise, the random functions keep giving new combinations until the requirements are met. 3
Applied Energy 253 (2019) 113561
Y. Hou, et al.
reached.
3. Repeat process 2 until the target volume fraction of the carbon is reached.
By comparing the obtained microstructure with the transmission electron microscope (TEM) image of CL shown in Fig. 1(a), it can be seen that the present reconstruction algorithm can well describe the morphology of CL, which proves its reliability.
2.3. Generate Pt particle In this work, it is assumed that all the Pt particles are cubes with a side length of 3 nm length, which just fit a boundary cell. To control the variables, the agglomerated Pt particle is not considered. Therefore, only the pore nodes adjacent to a carbon node can be converted to Pt nodes. To be specific, as shown in Fig. 1(d), if node 0 is carbon node, only the pore nodes from 1 to 6 can be potentially changed. The algorithm works as:
2.5. CL characteristics analyses After the reconstruction, some important characteristics of the reconstructed CLs need to be extracted and analyzed. In the nanoscale pore, the effect of Knudsen diffusivity should be counted, which requires the pore diameter data. To consider the Knudsen diffusivity more realistically, instead of using the mean pore diameter, the effective diameter at each pore node is obtained by the 13-length method [22]. The connectivity is another important characteristic which affects both transport and reactive area. In this work, the carbon-generating algorithm itself can ensure all the carbon spheres are connected. Here, a connectivity test function is developed for ionomer nodes and works as.
1. Each potential Pt node is identified and assigned a random value (1–1,000,000 in this work). 2. If the assigned value is below the specific value (10 in this work), the pore node is converted to a Pt node. Otherwise, it remains. 3. Repeat such processes until the target volume fraction of Pt is reached. Note that the conversion probability should be sufficiently low to ensure the relatively uniform distribution of Pt particles.
1. As shown in Fig. 2, in the thickness (x) direction, ionomer nodes on the rightmost layer (adjacent to PEM) are set as active nodes. 2. The active ionomer nodes (for example, node 0 in Fig. 1(d)) can further activate the neighboring ionomer nodes (nodes 1–6). 3. Repeat process 2 until the number of active ionomer nodes does not increase.
2.4. Generate ionomer It should be mentioned that under the same Pt/C ratio, the carbon supported Pt frame is only reconstructed once and applied to various I/ C conditions to reduce the error caused by reconstruction. And it is reported in Ref. [23] that the effective diffusivity can be overstated for uniformly distributed ionomer. Therefore, the ionomer phase is generated in the same method as Pt, which can ensure the nonuniform distribution.
The test results are listed in Table 2, indicating that the present ionomer generating algorithm can generally ensure the connectivity of ionomer. And the obtained connective rates are also close to the results in Ref. [24]. After checking the connectivity, the data of the active catalyst area are also required. Since the electrochemical reaction takes place at the ionomer/Pt interface, a function is developed to go through the entire domain and identify the interface between Pt and connected ionomer. The total active catalyst areas under various conditions are given in Table 2. It should be mentioned that the effective diffusivity of the reconstructed CL is 6.0346154 × 10−7 m2 s - 1 and in the same order of magnitude with Ref. [24,29], which further validates the reconstruction
1. Each potential ionomer node (pore node adjacent to carbon, Pt, or ionomer) is identified and assigned a random value (1–10,000 in this work). 2. If the assigned value is below the specific value (10 in this work), the pore node is converted to an ionomer node. Otherwise, it remains. 3. Repeat such processes until the target volume fraction of ionomer is
Fig. 2. Cross sections of the reconstructed catalyst layers under various Pt/C ratios and structure designs. 4
Applied Energy 253 (2019) 113561
Y. Hou, et al.
Table 2 Detailed information of the reconstructed CCLs. Case Number
Porosity
θPt
Thickness(μm)
θI
Connected rate
Reactive area (nm2)
R(nm)
1 2 3 4 5 6 7 8 9 10 11
0.350 0.397 0.284 0.350 0.415 0.268 0.350 0.436 0.245 0.271 0.279
0.4 0.4 0.4 0.6 0.6 0.6 0.8 0.8 0.8 0.6 0.6
4.986 4.986 4.986 2.670 2.670 2.670 1.518 1.518 1.518 2.670 2.670
1:3 1:4 1:2.3 1:3 1:4 1:2.3 1:3 1:4 1:2.3 1:2.3 1:2.3
99.9% 99.7% 99.9% 99.9% 99.6% 99.9% 99.9% 99.9% 99.9% 99.9% 99.9%
15,034,707 11,630,232 18,744,255 15,171,057 11,964,402 18,815,148 17,856,180 14,683,032 21,172,023 18,422,622 18,145,449
0 0 0 0 0 0 0 0 0 15 30
method.
Table 3 Comparison of simulation results with other works. Data from
Effective diffusivity
3. Numerical model development
This work Bruggeman [34] Lange [20] Chen [22]
0.3312 0.3288 0.3374 0.3317
3.1. Model description Restricted by computation resources and time scale gap, it is still impractical for a single model to consider all the physical and electrochemical phenomena in CCLs. Specifically, it involves the transports of gas, electron, proton, water and heat, electrochemical reaction and even the swelling of the ionomer, etc. Therefore, this work neglects the effects of liquid and vapor and assumes the CCLs operate under isothermal and steady condition. And since negligible overpotential (η ) drops across the CLs are observed in Ref. [24,30], a uniform distribution of overpotential is assumed in this work. Based on the above assumptions, the present model mainly focuses on the reactive transport processes, including the oxygen diffusion in pores and ionomer, oxygen diffusion from pores into ionomer, and the electrochemical reactions on reactive sites. All of which is believed to be sufficient to study the effects of ingredient contents and structures. In the CCL, the oxygen diffusion is dominated by concentration difference, and the governing equation is shown as (5)
∇ ·(D∇C ) = S −3
2
−1
where C (mol m ) is oxygen concentration, D (m s ) is its diffusivity and S represents the reaction source term. As aforementioned, the diffusivity in nanoscale pores should be coupled by bulk diffusivity and Knudsen diffusivity as
Fig. 3. Comparison of the current density under various conditions.
−1
1 1 ⎞ Dp = ⎛ + D kn ⎠ ⎝ Db ⎜
⎟
D b = 0.22 × 10−4
D kn = 48.5d p
(T /293.2)1.5 P /1 atm
T 32
(6)
(7)
(8)
where d p (m) is the mentioned local pore size. And the oxygen diffusivity in ionomer is obtained by
D I = 10−10 × (0.1543(T − 273) − 1.65)
(9)
Henry’s law solves the equilibrium oxygen concentration in ionomer.
CI =
Fig. 4. Comparison of the oxygen distribution along x direction under various Pt/C ratios.
Cp H
(10)
where H is Henry's constant [31]. Note that the discontinuous concentration across the pore/ionomer interface may cause model instability. The hypothetical concentration method proposed by Chen [24] is adopted, which considers the effect of Henry’s law in the diffusivity term as 5
Applied Energy 253 (2019) 113561
Y. Hou, et al.
(c)
(d)
Fig. 5. Data analyses under Pt/C = 40% condition (a) distribution of active catalyst area along x direction; (b) pore size distribution; (c) distribution of oxygen concentration in ionomer along x direction; (d) distribution of oxygen reaction rate along x direction.
d ⎡ dCI ⎤ d ⎡ D I dCp ⎤ DI → ⎥ ⎢ H dx ⎦ dx ⎣ dx ⎦ dx ⎣
(12)
(13)
αF η⎞ ds ∫ i0 CCrefI exp ⎛⎝ RT ⎠
(16)
3.2. Numerical method Owing to its simplicity of formulation and versatility, the lattice Boltzmann method is widely applied in multiscale flows and interfacial dynamics, especially in irregular complex geometries [32–34]. More details and applications of this promising numerical approach can be found in these representative reviews [35–37]. In this work, the LBM is employed to simulate the reactive transport processes in the CCL. With a view to computation accuracy and efficiency, the D3Q19 lattice model with a single relaxation time collision operator is adopted, and the evolution equation for the oxygen concentration distribution can be expressed as:
where n is the number of the neighboring active Pt sites and Δ is the lattice resolution. The total current density of the computation domain can be obtained by
i total A =
∂C =0 ∂x
Periodic boundary conditions are applied to the other two directions. In the boundary array, the pore and ionomer nodes are numbered 0 and 1, respectively, which represents the fluid zone for the computational execution. Carbon and Pt nodes are numbered 2 and 3, respectively, which belong to the solid zone and implement the non-slip boundary condition [25].
where i 0 is the exchange current density, Cref is the reference oxygen concentration, α is the cathode transfer coefficient, F represents the Faraday’s constant and R is the universal gas constant. Then, the local ORR rate can be related to the current as
ni S= Δ4F
(15)
Right boundary:
Regarding the reaction source term, the local current density i is calculated by the Butler-Volmer equation as
C αF η⎞ i = i 0 I exp ⎛ Cref ⎝ RT ⎠
Left boundary:C = Cinlet (11)
(14)
where A is the area of the cross-section (perpendicular to x), and the right side is the integration of the local current density in the entire CCL. Since the overpotential is uniformly distributed, the total current density is therefore only related to oxygen concentration and active catalyst area. Regarding the boundary conditions, along the x-direction shown in Fig. 2, the left side is assumed to be attached to the TPL and right side to the PEM. Therefore, a constant concentration is applied to the left boundary and a non-flux boundary on the right side as:
1 fi (x + ei Δt , t + Δt) − fi (x, t) = − (fi (x, t) − f eq i (x , t)) + α i SΔt τ
(17)
where fi (x, t ) represents the discrete distribution function at lattice site x and lattice time t along i direction, fieq (x, t) denotes its equilibrium distribution and S stands for the oxygen depletion source term. The direction distribution in D3Q19 model is shown in Fig. 1(d)) and the 6
Applied Energy 253 (2019) 113561
Y. Hou, et al.
Fig. 6. Data analyses under Pt/C = 60% condition (a) distribution of active catalyst area along x direction; (b) distribution of oxygen concentration in ionomer along x direction; (c) distribution of oxygen reaction rate along x direction.
discrete lattice velocity ei can be written as
3.3. Model validation
i=0 ⎧ (0, 0, 0) ei = (± 1, 0, 0), (0, ± 1, 0), (0, 0, ± 1) i=1−6 ⎨ ⎩ (± 1, ± 1, 0), (0, ± 1, ± 1), (± 1, 0, ± 1) i = 7 − 18
Following the methods in Refs. [22,24]; the code is validated by conducting the Bruggeman correlation benchmark [38]. To be specific, simulations run on a reconstructed 300 × 300 × 300 nm cubic domain with a resolution of 3 nm, and the structure of the domain is identical to the experiment. The predicted effective diffusivities by our code is compared with the results from Lange [22], Chen [24] and the experiment [38], as shown in Table 3, which verifies the model accuracy.
(18)
Δt = 1 and Δx =1 are the nondimensional time step and lattice spacing, respectively. Since the convection is neglected, the equilibrium distribution function fieq (x, t) is given by i=0 ⎧1/3 fieq = αi C , αi = 1/18 i = 1 − 6 ⎨ ⎩1/36 i = 7 − 18
4. Results and discussion (19)
In this work, a total of 11 cases are conducted with various combinations of I/C ratios (1:2.3, 1:3, 1:4), Pt/C ratios (0.4, 0.6, 0.8) and structure designs (R = 0, R = 15, R = 30 nm ). Table 2 gives more details of the reconstructed CCLs like thickness, porosity, etc. The output current densities under these conditions are shown in Fig. 3 and its structural characteristics, concentration distribution, and reaction rate distribution are also analyzed by the algorithms above and models.
The relaxation time τ is determined by the diffusivity as
D=
1 Δx 2 (τ − 0.5) 3 Δt
(20)
And since the local diffusivity varies with the pore size, the relaxation time should also be calculated separately. The concentration can be obtained by
C=
∑ fi
4.1. Effects of Pt/C ratio In this part, Cases 1, 4, 7 with the reference I/C but different Pt/C ratios are compared. As shown in Fig. 2, under the constant Pt loading, a lower Pt/C ratio yields a thicker CCL, which is consistent with the experiment [11]. Due to the disorganized low-porosity structure, mass transport loss is unavoidable, and various degrees of loss can be observed in Fig. 4. Besides, it has been found that the transport loss increases sharply with the decline in Pt/C ratio due to the longer tortuous transport path. In addition to the severe transport loss, Table 2 shows
(21)
Note that to ensure the computation accuracy and efficiency, the conversion factors which connect the lattice and physical space should be chosen properly. In this work, the length and time scale are chosen as 3 × 10−9 m and 10−11 s , respectively. All the mentioned models are realized by an in-house code, which is based on CUDA and parallelized on a GPU workstation. 7
Applied Energy 253 (2019) 113561
Y. Hou, et al.
Fig. 7. Data analyses under Pt/C = 80% condition (a) distribution of active catalyst area along x direction; (b) distribution of oxygen concentration in ionomer along x direction; (c) distribution of oxygen reaction rate along x direction.
To comprehensively evaluate the effects of I/C ratio, the integrations of the reaction rate in cross-sections are calculated, and its distributions are presented in Figs. 5(d), 6(c), 7(c). The integration of reaction rate is both affected by local oxygen concentration, and active catalyst area, and so is its variation trend. Under the low Pt/C ratio conditions (Figs. 5 and 6) where the transport loss is relatively serious, the variation trend of reaction rate is generally similar to the trend of oxygen concentration. While under the low transport loss condition (Fig. 7), it is dominated by the variation trend of active catalyst area. However, the fundamental rules are the same. The reaction rate benefits more from a larger active catalyst area than oxygen concentration, which explains why the reaction rate under high I/C condition is mostly higher in these plots. And Fig. 3 can further prove it as the highest I/C ratio yields the best performance. So, the active catalyst area should be given precedence when choosing the I/C ratio.
that under the same I/C ratio, the total reactive area also increases with the rise in Pt/C ratio because it is easier for ionomer to coat the densely distributed Pt particles as shown in Figs. 5(a), 6(a) and 7(a). Therefore, it is shown in Fig. 3 that the cell performance can be greatly improved by using high Pt/C ratio catalyst on account of the significantly higher oxygen concentration and larger active catalyst area. 4.2. Effects of I/C ratio In this part, for each Pt/C ratio, cases with various I/C ratios (1:2.3, 1:3 and 1:4) are compared. In Fig. 5(b), it can be found that the pore size distribution is sensitively affected by I/C ratio and local pore size tends to decrease with the rise in I/C ratio, which is also consistent with the experiment [9]. The variation trend of PSD in other Pt/C conditions are similar and therefore not presented. The oxygen diffusivity in the pore is related to the local pore size as explained in Eqs. (6) and (8), so a smaller mean pore size results in the decrease of oxygen diffusivity. Besides, the diffusivity in the pore is orders of magnitude higher than that in the ionomer. Therefore, Figs. 5(c), 6(b) and 7(b) all show the rise in I/C ratio causes more transport loss, which is against the reaction rate. On the other hand, ionomer serves as proton conductor, and only the Pt particles coated by connected ionomer can catalyze the reaction. Here, the ionomer connection rate and the active catalyst area are analyzed. The results presented in Table 2 indicate the chosen I/C ratios can all ensure a high connection rate. And it also shows I/C ratio has a huge influence on the active catalyst area as more ionomer is capable of activating more Pt particles. The same result is also observed in Figs. 5(a), 6(a) and 7(a). In this respect, high I/C ratio can also improve the performance by providing a larger active catalyst area.
4.3. Effects of structure design Based on the above discussion, it can be concluded that CCL is required to ensure a large active catalyst area on the premise of low transport loss to achieve better performance. However, to activate more Pt sites, more ionomer is indispensable, and the corresponding transport resistance is also inevitable for the traditional CCL. Therefore, we wonder if it is possible to develop a novel structured CCL to optimize the reactive area and transport simultaneously. The proposed structure design is shown in Fig. 2, where a nanoscale round hole is straightly punched through the CCL. Except for the hole, the carbon frame, ionomer, and Pt distributions are all identical to Case 6 (Pt/C = 0.6, I/ C = 2.3). As shown in Table 2 and Fig. 8(a), due to the presence of the hole, the active catalyst area is just slightly smaller since some Pt and 8
Applied Energy 253 (2019) 113561
Y. Hou, et al.
Fig. 8. Effects of the proposed structure design on (a) distribution of active catalyst area along x direction; (b) distribution of oxygen concentration in ionomer along x direction; (c) distribution of oxygen reaction rate along x direction.
proposed, which can simultaneously ensure large active catalyst area and extremely low transport loss and enhance the performance by 50%. It proves that the structured electrode can dramatically improve the performance, and there is still huge potential for proton exchange membrane fuel cell.
ionomer nodes are converted to pore nodes. However, the benefit it brings is tremendous. Due to the disappearance of a part of the porous structure, the local mass transport resistance is significantly decreased. Besides, local diffusivity is also increased due to the rise in local pore size. Therefore, Fig. 8(b) shows the proposed structure can greatly reduce the transport loss, and a relatively larger radius can better enhance the transport. Therefore, the performance is greatly improved, as presented in Figs. 8(c) and 3 shows the current density increases by up to 50%, which also indicates that there is still huge potential for CCL to develop with the help of novel structural design.
Acknowledgments This work is supported by the National Key Research and Development Program of China (2018YFB0105505), the China-UK International Cooperation and Exchange Project (Newton Advanced Fellowship) jointly supported by the National Natural Science Foundation of China (Grant No. 51861130359) and the UK Royal Society (Grant No. NAF/R1/180146), and the National Natural Science Foundation of Tianjin (China) for Distinguished Young Scholars (Grant No. 18JCJQJC46700).
5. Conclusion In this work, the effects of Pt/C ratio, I/C ratio, and structure design are studied in terms of oxygen concentration distribution, reactive area, and reaction rate. The reactive transport processes in the cathode catalyst layer are simulated with LBM based pore-scale model. The cathode catalyst layers under various combinations of parameters and structure designs are elaborately reconstructed with the consideration of carbon, Pt, ionomer, and pores. Developed algorithms are employed to ensure the conductivity and relatively uniform distributions of these ingredients. The results indicate that a higher Pt/C ratio can significantly enhance mass transport and therefore improve the performance. Ionomer plays an ambivalent role as it’s essential for activating Pt sites while against the mass transport. However, it’s found that the performance benefits more from a larger active catalyst area as high I/C ratio yields better performance. Therefore, the active catalyst area should be prioritized, and the optimal I/C ratio should be just enough to activate most of the catalyst area. Finally, an ideal structure design is
References [1] Wang Y, Chen KS, Mishler J, Cho SC, Adroher XC. A review of polymer electrolyte membrane fuel cells: technology, applications, and needs on fundamental research. Appl Energy 2011;88(4):981–1007. [2] Huo S, Jiao K, Park JW. On the water transport behavior and phase transition mechanisms in cold start operation of PEM fuel cell. Appl Energy 2019;233:776–88. [3] Chang Y, Qin Y, Yin Y, Zhang J, Li X. Humidification strategy for polymer electrolyte membrane fuel cells–a review. Appl Energy 2018;230:643–62. [4] Wu HW. A review of recent development: transport and performance modeling of PEM fuel cells. Appl Energy 2016;165:81–106. [5] Perng SW, Wu HW. Effect of the prominent catalyst layer surface on reactant gas transport and cell performance at the cathodic side of a PEMFC. Appl Energy 2010;87(4):1386–99. [6] Zamel N. The catalyst layer and its dimensionality – a look into its ingredients and
9
Applied Energy 253 (2019) 113561
Y. Hou, et al.
[23] Lange KJ, Sui PC, Djilali N. Pore scale modeling of a proton exchange membrane fuel cell catalyst layer: effects of water vapor and temperature. J Power Sources 2011;196(6):3195–203. [24] Chen L, Wu G, Holby EF, Zelenay P, Tao WQ, Kang Q. Lattice Boltzmann pore-scale investigation of coupled physical-electrochemical processes in C/Pt and non-precious metal cathode catalyst layers in proton exchange membrane fuel cells. Electrochim Acta 2015;158:175–86. [25] Salomov UR, Chiavazzo E, Asinari P. Pore-scale modeling of fluid flow through gas diffusion and catalyst layers for high temperature proton exchange membrane (HTPEM) fuel cells. Comput Math Appl 2014;67(2):393–411. [26] Kim KH, Lee KY, Kim HJ, Cho E, Lee SY, Lim TH, et al. The effects of Nafion® ionomer content in PEMFC MEAs prepared by a catalyst-coated membrane (CCM) spraying method. Int J Hydrogen Energy 2010;35(5):2119–26. [27] Siddique NA, Liu F. Process based reconstruction and simulation of a three-dimensional fuel cell catalyst layer. Electrochim Acta 2010;55(19):5357–66. [28] Ishikawa H, Sugawara Y, Inoue G, Kawase M. Effects of Pt and ionomer ratios on the structure of catalyst layer: a theoretical model for polymer electrolyte fuel cells. J Power Sources 2018;374:196–204. [29] Yu Z, Carter RN. Measurement of effective oxygen diffusivity in electrodes for proton exchange membrane fuel cells. J Power Sources 2010;195(4):1079–84. [30] Zhang J, Yang W, Xu L, Wang Y. Simulation of the catalyst layer in PEMFC based on a novel two-phase lattice model. Electrochim Acta 2011;56(20):6912–8. [31] Bernardi DM, Verbrugge MW. A mathematical model of the solid-polymer-electrolyte fuel cell. J Electrochem Soc 1992;139(9):2477–91. [32] Hou Y, Deng H, Du Q, Jiao K. Multi-component multi-phase lattice Boltzmann modeling of droplet coalescence in flow channel of fuel cell. J Power Sources 2018;393:83–91. [33] Jiang ZY, Qu ZG, Zhou L, Tao WQ. A microscopic investigation of ion and electron transport in lithium-ion battery porous electrodes using the lattice Boltzmann method. Appl Energy 2017;194:530–9. [34] Li X, Ma T, Liu J, Zhang H, Wang Q. Pore-scale investigation of gravity effects on phase change heat transfer characteristics using lattice Boltzmann method. Appl Energy 2018;222:92–103. [35] Chen L, Kang Q, Mu Y, He YL, Tao WQ. A critical review of the pseudopotential multiphase lattice Boltzmann model: methods and applications. Int J Heat Mass Transf 2014;76:210–36. [36] Aidun CK, Clausen JR. Lattice-Boltzmann method for complex flows. Annu Rev Fluid Mech 2010;42:439–72. [37] Li Q, Luo KH, Kang QJ, Chen Q, Liu Q. Lattice Boltzmann methods for multiphase flow and phase-change heat transfer. Prog Energy Combust Sci 2016;52:62–105. [38] Bruggeman VDAG. Berechnung verschiedener physikalischer Konstanten von heterogenen Substanzen. I. Dielektrizitätskonstanten und Leitfähigkeiten der Mischkörper aus isotropen Substanzen. Ann Phys 1935;416(7):636–64.
how to characterize their effects. J Power Sources 2016;309:141–59. [7] Lee D, Hwang S. Effect of loading and distributions of Nafion ionomer in the catalyst layer for PEMFCs. Int J Hydrogen Energy 2008;33(11):2790–4. [8] Ozden A, Shahgaldi S, Li X, Hamdullahpur F. The impact of ionomer type on the morphological and microstructural degradations of proton exchange membrane fuel cell electrodes under freeze-thaw cycles. Appl Energy 2019;238:1048–59. [9] Shahgaldi S, Alaefour I, Li X. Impact of manufacturing processes on proton exchange membrane fuel cell performance. Appl Energy 2018;225:1022–32. [10] Shahgaldi S, Ozden A, Li X, Hamdullahpur F. Cathode catalyst layer design with gradients of ionomer distribution for proton exchange membrane fuel cells. Energy Convers Manage 2018;171:1476–86. [11] Zhao J, Ozden A, Shahgaldi S, Alaefour IE, Li X. Effect of Pt loading and catalyst type on the pore structure of porous electrodes in polymer electrolyte membrane (PEM) fuel cells. Energy 2018;150:69–76. [12] Zhao J, Shahgaldi S, Alaefour I, Xu Q, Li X. Gas permeability of catalyzed electrodes in polymer electrolyte membrane fuel cells. Appl Energy 2018;209:203–10. [13] Jiao K, Zhou Y, Du Q, Yin Y, Yu S, Li X. Numerical simulations of carbon monoxide poisoning in high temperature proton exchange membrane fuel cells with various flow channel designs. Appl Energy 2013;104:21–41. [14] Jiao K, Bachman J, Zhou Y, Park JW. Effect of induced cross flow on flow pattern and performance of proton exchange membrane fuel cell. Appl Energy 2014;115:75–82. [15] Li W, Zhang Q, Wang C, Yan X, Shen S, Xia G, et al. Experimental and numerical analysis of a three-dimensional flow field for PEMFCs. Appl Energy 2017;195:278–88. [16] Ferreira RB, Falcão DS, Oliveira VB, Pinto A. 1D+ 3D two-phase flow numerical model of a proton exchange membrane fuel cell. Appl Energy 2017;203:474–95. [17] Rao Z, Wang S, Peng F. Self diffusion of the nano-encapsulated phase change materials: a molecular dynamics study. Appl Energy 2012;100:303–8. [18] Andersson M, Yuan J, Sundén B. Review on modeling development for multiscale chemical reactions coupled transport phenomena in solid oxide fuel cells. Appl Energy 2010;87(5):1461–76. [19] Chen L, Zhang R, He P, Kang Q, He YL, Tao WQ. Nanoscale simulation of local gas transport in catalyst layers of proton exchange membrane fuel cells. J Power Sources 2018;400:114–25. [20] Wang G, Mukherjee PP, Wang CY. Direct numerical simulation (DNS) modeling of PEFC electrodes: Part I. Regular microstructure. Electrochim Acta 2006;51(15):3139–50. [21] Wang G, Mukherjee PP, Wang CY. Optimization of polymer electrolyte fuel cell cathode catalyst layers via direct numerical simulation modeling. Electrochim Acta 2007;52(22):6367–77. [22] Lange KJ, Sui PC, Djilali N. Pore scale simulation of transport and electrochemical reactions in reconstructed PEMFC catalyst layers. J Electrochem Soc 2010;157(10):B1434–42.
10
Contents lists available at ScienceDirect
Applied Energy journal homepage: www.elsevier.com/locate/apenergy
Pore-scale investigation of catalyst layer ingredient and structure effect in proton exchange membrane fuel cell ⁎
Yuze Houa, Hao Denga, Fengwen Panb, Wenmiao Chenb, Qing Dua, , Kui Jiaoa, a b
T
⁎
State Key Laboratory of Engines, Tianjin University, 135 Yaguan Road, Tianjin 300350, China Weichai Power Co. Ltd., 197A Fushou St. E., Weifang 261000, China
H I GH L IG H T S
layers under various ingredient contents and structures are reconstructed. • Catalyst pore-scale model is developed to simulate the reactive transport processes. • AHigh platinum/catalyst ratio can enhance the transport and improve the performance. • The active catalyst area should be given precedence during the electrode fabrication. • A novel catalyst layer design is proposed and can improve the performance by 50%. •
A R T I C LE I N FO
A B S T R A C T
Keywords: Proton exchange membrane fuel cell (PEMFC) Catalyst layer Lattice Boltzmann method (LBM) Platinum/catalyst (Pt/C) ratio Ionomer/catalyst (I/C) ratio Electrode structure design
A pore-scale model based on the lattice Boltzmann method (LBM) is developed to simulate the reactive transport processes in the cathode catalyst layer of a proton exchange membrane fuel cell (PEMFC). The porous structures of the cathode catalyst layers are reconstructed in the process-based method with the consideration of carbon supporter, platinum, ionomer and pores. Its characteristics are analyzed including pore size distribution, phase connectivity and active catalyst area. The effects of two critical parameters, platinum/catalyst (Pt/C) and ionomer/catalyst (I/C) ratios, and structure design are investigated in terms of oxygen concentration distribution, reactive area, and reaction rate. The results indicate that, under the constant platinum loading (0.3 mg cm−2), a higher Pt/C ratio yields a thinner catalyst layer, which significantly enhances the oxygen transport and improves the performance. For the same Pt/C ratio, although a higher I/C ratio brings more mass transport loss, it increases the active catalyst area and ultimately yields better performance. Therefore, the active catalyst area should be given precedence during catalyst layer fabrication. To realize a large active catalyst area on the premise of low transport loss, an ideal catalyst layer structure design is proposed and capable of improving the performance by 50%.
1. Introduction Proton exchange membrane fuel cell (PEMFC) is considered as one of the most promising energy conversion devices for its attractive features such as high efficiency and fast response [1,2]. Among the components of a PEMFC, the catalyst layer (CL) sandwiched between the porous transport layer (PTL) and the PEM is the most complex and important since it is the place where energy conversion occurs [3]. Specifically, the CL is required to be heterogeneous for the transport of protons, electrons, reactants and products, and also provide an active catalyst area for electrochemical reactions [4]. To meet these needs, the state-of-the-art CL is porous in structure and commonly composed of
⁎
carbon-supported platinum-based catalyst (platinum particles supported on carbon, as shown in Fig. 1(a)), ionomer and pores [5]. Additionally, to achieve more efficient transport, larger active catalyst area, and ultimately better cell performance, the content of each ingredient and the structure should be optimized delicately [6]. Experimentally, Lee and Hwang [7] studied the effects of ionomer loading and distribution and found distribution schemes can influence the cell performance. Ozden et al. [8] investigated the effects of ionomer type and observed that both the loading and type of ionomer affect the degradation pattern and performance of PEMFC. Shahgaldi et al. [9] studied the electrode manufacturing processes and found the manufacturing strategy should vary with the Pt loading to reach a
Corresponding authors. E-mail addresses: [email protected] (Q. Du), [email protected] (K. Jiao).
https://doi.org/10.1016/j.apenergy.2019.113561 Received 22 April 2019; Received in revised form 14 July 2019; Accepted 15 July 2019 0306-2619/ © 2019 Published by Elsevier Ltd.
Applied Energy 253 (2019) 113561
Y. Hou, et al.
Nomenclature A C D d e F f H i i0 l n P R S T t x y z
Greek letters
area of the cross-section, m2 oxygen concentration, mol m−3 diffusivity, m2 s−1 local pore size, m discrete lattice velocity Faraday’s constant, C mol−1 discrete distribution function Henry’s constant local current density, A m−2; discrete direction exchange current density, A m−2 thickness of catalyst layer, nm neighboring active Pt node number pressure, atm universal gas constant; radius of hole structure, nm reaction source term temperature, K lattice time direction; lattice site direction direction
γPt θPt θI ε εPt εI εcarbon ρ α η τ Δ
platinum loading, mg cm−2 platinum/catalyst ratio ionomer/catalyst ratio porosity volume fraction of platinum volume fraction of ionomer volume fraction of carbon density, mg cm−3 cathode transfer coefficient overpotential, V relaxation time lattice resolution, m
Subscripts and superscripts b eq I kn p ref
bulk equilibrium ionomer Knudsen pore reference
Fig. 1. (a) TEM image of the catalyst layer; (b) Reconstructed microstructures of the catalyst layer (2670 × 450 × 450 nm); (c) Catalyst layer reconstruction processes; (d) Direction instruction for reconstruction and D3Q19 model.
capable of modeling the whole cell and predicting its performance. However, to achieve this, resolving the CL microstructure is neglected as it requires too much computation and empirical formulations are required to determine some critical transport parameters. Therefore, the macroscale model is not able to consider the effect of local structure in the CL. As for the microscopic model, the molecular dynamics (MD) model [17,18] is capable of considering the swelling of ionomer and diffusivity at the molecular scale. While aiming for the balance of ingredient and structure optimization, the computation domain should be in micrometer size, which is far beyond the scope of MD. In contrast, the mesoscopic pore-scale model, which can model the detailed processes in the CL with acceptable computational cost, is the most appropriate method for the mentioned targets and also the trend of CL modeling [19]. Wang et al. [20,21] simulated the transport of oxygen and proton in simplified CLs with the finite volume method. Lange et al. [22,23] further considered the local Knudsen diffusion effects and
better performance. They also designed a cathode catalyst layer (CCL) with two sub-layers filled by different ionomer loadings and concluded that an appropriate loading gradient could be beneficial for both transport and performance [10]. Zhao et al. [11] found that under constant Pt loading, the catalyst with low platinum content results in a thicker electrode. In their further study [12], the gas permeability in the catalyst electrode was tested, and it was found that the presence of a CL significantly reduced the gas permeability. However, it is too time consuming for an experiment to be conducted on every combination of variables. Due to the nanometer scale, it is still impractical for advanced equipment to realize in-situ observation of transport processes inside the CL and also difficult to control its structure and variables accurately at will [4]. Therefore, numerical simulation is considered a complementary approach to obtain a fundamental understanding of the underlying mechanisms in the CL [13]. Through the existing studies, the macroscale model [14–16] is 2
Applied Energy 253 (2019) 113561
Y. Hou, et al.
conducted to check the effects on the reconstructed structures, such as pore size distribution (PSD), ionomer connectivity, and active catalyst area. The specific implementation algorithms of reconstruction and analyses are shown below.
Table 1 Parameter values in this work. Parameters
Symbol
Value
Units
platinum loading platinum density ionomer density carbon density overpotential Temperature Pressure Henry’s constant exchange current density reference oxygen concentration Cathode transfer coefficient Oxygen concentration at TPL/CL interface Lattice resolution
γPt ρPt ρI ρcarbon η T P H i0 Cref α Cinlet Δ
0.3 21,450 2000 1800 0.3 353 1.5 6.955 0.015 40.96 0.616 10.875 3
mg cm−2 mg cm−3 mg cm−3 mg cm−3 V K atm – A m−2 mol m−3 – mol m−3 nm
2.1. Calculate the ingredient contents In this work, the constant Pt loading (γPt ) is 0.3 mg cm−2. According to Ref. [11], a lower Pt/C ratio (θPt ) results in a thicker electrode. And it is reported in Ref. [10] that the I/C ratio (θ I ) and porosity are in inverse proportion since the ionomer occupies the space of the pore and the I/C ratio between 1:3 and 1:2.3 is expected to yield better performance and durability. Therefore, it is rational to make the following assumptions: 1. Under the reference I/C ratio (1:3), the porosity (ε ) of CL is constant (35%), and the effect of Pt/C ratio on the structure is considered by changing CL thickness. 2. Under the same Pt/C ratio, the thickness (l) of CL is not affected by I/C ratio.
reconstructed the CL with a resolution of 2 nm. The CL was represented as carbon spheres covered by thin ionomer with a specified thickness. Chen et al. [24] firstly modeled the reactive transport processes with the lattice Boltzmann method (LBM) and used quartet structure generation set (QSGS) to reconstruct the four-constituent CL microstructure, which was considered a great advancement in model development. And in their later work, the dissolution resistance is coupled into the model, and the schemes for reducing transport resistance are proposed [19]. Salomov et al. [25] conducted the pore scale simulations to predict the permeability in the reconstructed GDL and CL. Through a brief review, it can be demonstrated that a pore-scale model is a powerful tool for gaining valuable insights into the CL. With such a potential model, many unexplored but meaningful works can be realized. For CL fabrication, the contents of ingredients are described by two parameters, platinum/catalyst ratio (Pt/C, note that in this paper C stands for catalyst whose weight is the sum of platinum and carbon) and ionomer/catalyst (I/C) ratio. And these two critical parameters are usually set by experience or repeated tests, which can be inappropriate and time-consuming [26]. Besides, the structured CL is promising with the capacity to improve the cell performance significantly, but difficult to fabricate [6]. To the best of our knowledge, there has been no pore-scale based study to work on the above subjects. In this work, CLs under various parameters and structure designs are reconstructed elaborately with the consideration of carbon, Pt, ionomer, and pores, as shown in Fig. 1(b). And the LBM based pore-scale model is developed to simulate the reactive transport processes in the reconstructed CLs to offer valuable advice for fabricating CL and propose novel structure designs. Due to the relatively low oxygen reduction reaction (ORR) rate, this work focuses on the cathode side. The large-scale simulations are realized by in-house code computed in parallel using graphics processing unit (GPU) workstation. In the following, the algorithms of reconstruction and characteristic analyses are presented in Section 2, and the LBM based pore-scale model is developed in Section 3. The results and discussion are given in Section 4, while the conclusions are drawn in Section 5.
Regarding the reconstructed CL size, both the width and height are 450 nm (y and z-direction in Fig. 1(b)), and the thickness is determined by the Pt/C ratio. The volume fraction of each ingredient (εPt , εI and εcarbon ) and the CL thickness can be figured out by:
εPt = (1 − ε )
εI = (1 − ε )
θPt / ρPt θ I/ ρI + (1 − θPt )/ ρcarbon + θPt / ρPt
θ I / ρI θ I/ ρI + (1 − θPt )/ ρcarbon + θPt / ρPt
εcarbon = (1 − ε )
l=
(1 − θPt )/ ρcarbon θ I/ ρI + (1 − θPt )/ ρcarbon + θPt / ρPt
γPt εPt ρPt
(1)
(2)
(3)
(4)
where ρPt , ρI and ρcarbon are the density of Pt, ionomer, and carbon, respectively. In this paper, the specific values of the parameters are given in Table 1. 2.2. Generate carbon support Various kinds of carbon, such as Ketjen Black and Vulcan XC-72, have distinctly different structural properties. Since such a factor is not focused on the present study, the morphology of carbon is simplified as a sphere, which is widely adopted [19,22,23]. Carbon spheres with random radii from 8 to 22 nm [10] are generated one by one until the required fraction εcarbon is reached. The specific algorithm inspired by Lange [22] is briefly explained below. 1. A group of random functions gives the center coordinate and radius of the first sphere. Then the first carbon sphere is generated based on the given coordinate and radius, and the corresponding pore nodes are converted to carbon nodes. 2. Another combination of radius and center coordinate is given randomly. And a function is set to check if the given combination meets the following two requirements: a. The newly generated carbon sphere should overlap with at least one existing carbon sphere. b. And the overlap rate cannot exceed the specified value (15% in this work).
2. CL reconstruction and analyses algorithms Microstructure reconstruction is the most basic element of the porescale CL model and many researchers have put much effort to build the microstructure more realistically with higher resolution [20,27,28]. Among the various reconstruction approaches, the process-based method proposed in Ref. [27] is widely adopted for its feasibility and rationality. Following the fabrication processes, the reconstruction algorithms can be divided into three steps: 1. generate carbon supporter, 2. generate Pt particles and 3. generate ionomer, as shown in Fig. 1(c). Weighing the computation domain size against the computation resource, the grid resolution is set as 3 nm per lattice cell, which is enough to consider the presence of Pt. The characteristic analyses are also
If yes, the next carbon sphere is generated based on the examined combination. Otherwise, the random functions keep giving new combinations until the requirements are met. 3
Applied Energy 253 (2019) 113561
Y. Hou, et al.
reached.
3. Repeat process 2 until the target volume fraction of the carbon is reached.
By comparing the obtained microstructure with the transmission electron microscope (TEM) image of CL shown in Fig. 1(a), it can be seen that the present reconstruction algorithm can well describe the morphology of CL, which proves its reliability.
2.3. Generate Pt particle In this work, it is assumed that all the Pt particles are cubes with a side length of 3 nm length, which just fit a boundary cell. To control the variables, the agglomerated Pt particle is not considered. Therefore, only the pore nodes adjacent to a carbon node can be converted to Pt nodes. To be specific, as shown in Fig. 1(d), if node 0 is carbon node, only the pore nodes from 1 to 6 can be potentially changed. The algorithm works as:
2.5. CL characteristics analyses After the reconstruction, some important characteristics of the reconstructed CLs need to be extracted and analyzed. In the nanoscale pore, the effect of Knudsen diffusivity should be counted, which requires the pore diameter data. To consider the Knudsen diffusivity more realistically, instead of using the mean pore diameter, the effective diameter at each pore node is obtained by the 13-length method [22]. The connectivity is another important characteristic which affects both transport and reactive area. In this work, the carbon-generating algorithm itself can ensure all the carbon spheres are connected. Here, a connectivity test function is developed for ionomer nodes and works as.
1. Each potential Pt node is identified and assigned a random value (1–1,000,000 in this work). 2. If the assigned value is below the specific value (10 in this work), the pore node is converted to a Pt node. Otherwise, it remains. 3. Repeat such processes until the target volume fraction of Pt is reached. Note that the conversion probability should be sufficiently low to ensure the relatively uniform distribution of Pt particles.
1. As shown in Fig. 2, in the thickness (x) direction, ionomer nodes on the rightmost layer (adjacent to PEM) are set as active nodes. 2. The active ionomer nodes (for example, node 0 in Fig. 1(d)) can further activate the neighboring ionomer nodes (nodes 1–6). 3. Repeat process 2 until the number of active ionomer nodes does not increase.
2.4. Generate ionomer It should be mentioned that under the same Pt/C ratio, the carbon supported Pt frame is only reconstructed once and applied to various I/ C conditions to reduce the error caused by reconstruction. And it is reported in Ref. [23] that the effective diffusivity can be overstated for uniformly distributed ionomer. Therefore, the ionomer phase is generated in the same method as Pt, which can ensure the nonuniform distribution.
The test results are listed in Table 2, indicating that the present ionomer generating algorithm can generally ensure the connectivity of ionomer. And the obtained connective rates are also close to the results in Ref. [24]. After checking the connectivity, the data of the active catalyst area are also required. Since the electrochemical reaction takes place at the ionomer/Pt interface, a function is developed to go through the entire domain and identify the interface between Pt and connected ionomer. The total active catalyst areas under various conditions are given in Table 2. It should be mentioned that the effective diffusivity of the reconstructed CL is 6.0346154 × 10−7 m2 s - 1 and in the same order of magnitude with Ref. [24,29], which further validates the reconstruction
1. Each potential ionomer node (pore node adjacent to carbon, Pt, or ionomer) is identified and assigned a random value (1–10,000 in this work). 2. If the assigned value is below the specific value (10 in this work), the pore node is converted to an ionomer node. Otherwise, it remains. 3. Repeat such processes until the target volume fraction of ionomer is
Fig. 2. Cross sections of the reconstructed catalyst layers under various Pt/C ratios and structure designs. 4
Applied Energy 253 (2019) 113561
Y. Hou, et al.
Table 2 Detailed information of the reconstructed CCLs. Case Number
Porosity
θPt
Thickness(μm)
θI
Connected rate
Reactive area (nm2)
R(nm)
1 2 3 4 5 6 7 8 9 10 11
0.350 0.397 0.284 0.350 0.415 0.268 0.350 0.436 0.245 0.271 0.279
0.4 0.4 0.4 0.6 0.6 0.6 0.8 0.8 0.8 0.6 0.6
4.986 4.986 4.986 2.670 2.670 2.670 1.518 1.518 1.518 2.670 2.670
1:3 1:4 1:2.3 1:3 1:4 1:2.3 1:3 1:4 1:2.3 1:2.3 1:2.3
99.9% 99.7% 99.9% 99.9% 99.6% 99.9% 99.9% 99.9% 99.9% 99.9% 99.9%
15,034,707 11,630,232 18,744,255 15,171,057 11,964,402 18,815,148 17,856,180 14,683,032 21,172,023 18,422,622 18,145,449
0 0 0 0 0 0 0 0 0 15 30
method.
Table 3 Comparison of simulation results with other works. Data from
Effective diffusivity
3. Numerical model development
This work Bruggeman [34] Lange [20] Chen [22]
0.3312 0.3288 0.3374 0.3317
3.1. Model description Restricted by computation resources and time scale gap, it is still impractical for a single model to consider all the physical and electrochemical phenomena in CCLs. Specifically, it involves the transports of gas, electron, proton, water and heat, electrochemical reaction and even the swelling of the ionomer, etc. Therefore, this work neglects the effects of liquid and vapor and assumes the CCLs operate under isothermal and steady condition. And since negligible overpotential (η ) drops across the CLs are observed in Ref. [24,30], a uniform distribution of overpotential is assumed in this work. Based on the above assumptions, the present model mainly focuses on the reactive transport processes, including the oxygen diffusion in pores and ionomer, oxygen diffusion from pores into ionomer, and the electrochemical reactions on reactive sites. All of which is believed to be sufficient to study the effects of ingredient contents and structures. In the CCL, the oxygen diffusion is dominated by concentration difference, and the governing equation is shown as (5)
∇ ·(D∇C ) = S −3
2
−1
where C (mol m ) is oxygen concentration, D (m s ) is its diffusivity and S represents the reaction source term. As aforementioned, the diffusivity in nanoscale pores should be coupled by bulk diffusivity and Knudsen diffusivity as
Fig. 3. Comparison of the current density under various conditions.
−1
1 1 ⎞ Dp = ⎛ + D kn ⎠ ⎝ Db ⎜
⎟
D b = 0.22 × 10−4
D kn = 48.5d p
(T /293.2)1.5 P /1 atm
T 32
(6)
(7)
(8)
where d p (m) is the mentioned local pore size. And the oxygen diffusivity in ionomer is obtained by
D I = 10−10 × (0.1543(T − 273) − 1.65)
(9)
Henry’s law solves the equilibrium oxygen concentration in ionomer.
CI =
Fig. 4. Comparison of the oxygen distribution along x direction under various Pt/C ratios.
Cp H
(10)
where H is Henry's constant [31]. Note that the discontinuous concentration across the pore/ionomer interface may cause model instability. The hypothetical concentration method proposed by Chen [24] is adopted, which considers the effect of Henry’s law in the diffusivity term as 5
Applied Energy 253 (2019) 113561
Y. Hou, et al.
(c)
(d)
Fig. 5. Data analyses under Pt/C = 40% condition (a) distribution of active catalyst area along x direction; (b) pore size distribution; (c) distribution of oxygen concentration in ionomer along x direction; (d) distribution of oxygen reaction rate along x direction.
d ⎡ dCI ⎤ d ⎡ D I dCp ⎤ DI → ⎥ ⎢ H dx ⎦ dx ⎣ dx ⎦ dx ⎣
(12)
(13)
αF η⎞ ds ∫ i0 CCrefI exp ⎛⎝ RT ⎠
(16)
3.2. Numerical method Owing to its simplicity of formulation and versatility, the lattice Boltzmann method is widely applied in multiscale flows and interfacial dynamics, especially in irregular complex geometries [32–34]. More details and applications of this promising numerical approach can be found in these representative reviews [35–37]. In this work, the LBM is employed to simulate the reactive transport processes in the CCL. With a view to computation accuracy and efficiency, the D3Q19 lattice model with a single relaxation time collision operator is adopted, and the evolution equation for the oxygen concentration distribution can be expressed as:
where n is the number of the neighboring active Pt sites and Δ is the lattice resolution. The total current density of the computation domain can be obtained by
i total A =
∂C =0 ∂x
Periodic boundary conditions are applied to the other two directions. In the boundary array, the pore and ionomer nodes are numbered 0 and 1, respectively, which represents the fluid zone for the computational execution. Carbon and Pt nodes are numbered 2 and 3, respectively, which belong to the solid zone and implement the non-slip boundary condition [25].
where i 0 is the exchange current density, Cref is the reference oxygen concentration, α is the cathode transfer coefficient, F represents the Faraday’s constant and R is the universal gas constant. Then, the local ORR rate can be related to the current as
ni S= Δ4F
(15)
Right boundary:
Regarding the reaction source term, the local current density i is calculated by the Butler-Volmer equation as
C αF η⎞ i = i 0 I exp ⎛ Cref ⎝ RT ⎠
Left boundary:C = Cinlet (11)
(14)
where A is the area of the cross-section (perpendicular to x), and the right side is the integration of the local current density in the entire CCL. Since the overpotential is uniformly distributed, the total current density is therefore only related to oxygen concentration and active catalyst area. Regarding the boundary conditions, along the x-direction shown in Fig. 2, the left side is assumed to be attached to the TPL and right side to the PEM. Therefore, a constant concentration is applied to the left boundary and a non-flux boundary on the right side as:
1 fi (x + ei Δt , t + Δt) − fi (x, t) = − (fi (x, t) − f eq i (x , t)) + α i SΔt τ
(17)
where fi (x, t ) represents the discrete distribution function at lattice site x and lattice time t along i direction, fieq (x, t) denotes its equilibrium distribution and S stands for the oxygen depletion source term. The direction distribution in D3Q19 model is shown in Fig. 1(d)) and the 6
Applied Energy 253 (2019) 113561
Y. Hou, et al.
Fig. 6. Data analyses under Pt/C = 60% condition (a) distribution of active catalyst area along x direction; (b) distribution of oxygen concentration in ionomer along x direction; (c) distribution of oxygen reaction rate along x direction.
discrete lattice velocity ei can be written as
3.3. Model validation
i=0 ⎧ (0, 0, 0) ei = (± 1, 0, 0), (0, ± 1, 0), (0, 0, ± 1) i=1−6 ⎨ ⎩ (± 1, ± 1, 0), (0, ± 1, ± 1), (± 1, 0, ± 1) i = 7 − 18
Following the methods in Refs. [22,24]; the code is validated by conducting the Bruggeman correlation benchmark [38]. To be specific, simulations run on a reconstructed 300 × 300 × 300 nm cubic domain with a resolution of 3 nm, and the structure of the domain is identical to the experiment. The predicted effective diffusivities by our code is compared with the results from Lange [22], Chen [24] and the experiment [38], as shown in Table 3, which verifies the model accuracy.
(18)
Δt = 1 and Δx =1 are the nondimensional time step and lattice spacing, respectively. Since the convection is neglected, the equilibrium distribution function fieq (x, t) is given by i=0 ⎧1/3 fieq = αi C , αi = 1/18 i = 1 − 6 ⎨ ⎩1/36 i = 7 − 18
4. Results and discussion (19)
In this work, a total of 11 cases are conducted with various combinations of I/C ratios (1:2.3, 1:3, 1:4), Pt/C ratios (0.4, 0.6, 0.8) and structure designs (R = 0, R = 15, R = 30 nm ). Table 2 gives more details of the reconstructed CCLs like thickness, porosity, etc. The output current densities under these conditions are shown in Fig. 3 and its structural characteristics, concentration distribution, and reaction rate distribution are also analyzed by the algorithms above and models.
The relaxation time τ is determined by the diffusivity as
D=
1 Δx 2 (τ − 0.5) 3 Δt
(20)
And since the local diffusivity varies with the pore size, the relaxation time should also be calculated separately. The concentration can be obtained by
C=
∑ fi
4.1. Effects of Pt/C ratio In this part, Cases 1, 4, 7 with the reference I/C but different Pt/C ratios are compared. As shown in Fig. 2, under the constant Pt loading, a lower Pt/C ratio yields a thicker CCL, which is consistent with the experiment [11]. Due to the disorganized low-porosity structure, mass transport loss is unavoidable, and various degrees of loss can be observed in Fig. 4. Besides, it has been found that the transport loss increases sharply with the decline in Pt/C ratio due to the longer tortuous transport path. In addition to the severe transport loss, Table 2 shows
(21)
Note that to ensure the computation accuracy and efficiency, the conversion factors which connect the lattice and physical space should be chosen properly. In this work, the length and time scale are chosen as 3 × 10−9 m and 10−11 s , respectively. All the mentioned models are realized by an in-house code, which is based on CUDA and parallelized on a GPU workstation. 7
Applied Energy 253 (2019) 113561
Y. Hou, et al.
Fig. 7. Data analyses under Pt/C = 80% condition (a) distribution of active catalyst area along x direction; (b) distribution of oxygen concentration in ionomer along x direction; (c) distribution of oxygen reaction rate along x direction.
To comprehensively evaluate the effects of I/C ratio, the integrations of the reaction rate in cross-sections are calculated, and its distributions are presented in Figs. 5(d), 6(c), 7(c). The integration of reaction rate is both affected by local oxygen concentration, and active catalyst area, and so is its variation trend. Under the low Pt/C ratio conditions (Figs. 5 and 6) where the transport loss is relatively serious, the variation trend of reaction rate is generally similar to the trend of oxygen concentration. While under the low transport loss condition (Fig. 7), it is dominated by the variation trend of active catalyst area. However, the fundamental rules are the same. The reaction rate benefits more from a larger active catalyst area than oxygen concentration, which explains why the reaction rate under high I/C condition is mostly higher in these plots. And Fig. 3 can further prove it as the highest I/C ratio yields the best performance. So, the active catalyst area should be given precedence when choosing the I/C ratio.
that under the same I/C ratio, the total reactive area also increases with the rise in Pt/C ratio because it is easier for ionomer to coat the densely distributed Pt particles as shown in Figs. 5(a), 6(a) and 7(a). Therefore, it is shown in Fig. 3 that the cell performance can be greatly improved by using high Pt/C ratio catalyst on account of the significantly higher oxygen concentration and larger active catalyst area. 4.2. Effects of I/C ratio In this part, for each Pt/C ratio, cases with various I/C ratios (1:2.3, 1:3 and 1:4) are compared. In Fig. 5(b), it can be found that the pore size distribution is sensitively affected by I/C ratio and local pore size tends to decrease with the rise in I/C ratio, which is also consistent with the experiment [9]. The variation trend of PSD in other Pt/C conditions are similar and therefore not presented. The oxygen diffusivity in the pore is related to the local pore size as explained in Eqs. (6) and (8), so a smaller mean pore size results in the decrease of oxygen diffusivity. Besides, the diffusivity in the pore is orders of magnitude higher than that in the ionomer. Therefore, Figs. 5(c), 6(b) and 7(b) all show the rise in I/C ratio causes more transport loss, which is against the reaction rate. On the other hand, ionomer serves as proton conductor, and only the Pt particles coated by connected ionomer can catalyze the reaction. Here, the ionomer connection rate and the active catalyst area are analyzed. The results presented in Table 2 indicate the chosen I/C ratios can all ensure a high connection rate. And it also shows I/C ratio has a huge influence on the active catalyst area as more ionomer is capable of activating more Pt particles. The same result is also observed in Figs. 5(a), 6(a) and 7(a). In this respect, high I/C ratio can also improve the performance by providing a larger active catalyst area.
4.3. Effects of structure design Based on the above discussion, it can be concluded that CCL is required to ensure a large active catalyst area on the premise of low transport loss to achieve better performance. However, to activate more Pt sites, more ionomer is indispensable, and the corresponding transport resistance is also inevitable for the traditional CCL. Therefore, we wonder if it is possible to develop a novel structured CCL to optimize the reactive area and transport simultaneously. The proposed structure design is shown in Fig. 2, where a nanoscale round hole is straightly punched through the CCL. Except for the hole, the carbon frame, ionomer, and Pt distributions are all identical to Case 6 (Pt/C = 0.6, I/ C = 2.3). As shown in Table 2 and Fig. 8(a), due to the presence of the hole, the active catalyst area is just slightly smaller since some Pt and 8
Applied Energy 253 (2019) 113561
Y. Hou, et al.
Fig. 8. Effects of the proposed structure design on (a) distribution of active catalyst area along x direction; (b) distribution of oxygen concentration in ionomer along x direction; (c) distribution of oxygen reaction rate along x direction.
proposed, which can simultaneously ensure large active catalyst area and extremely low transport loss and enhance the performance by 50%. It proves that the structured electrode can dramatically improve the performance, and there is still huge potential for proton exchange membrane fuel cell.
ionomer nodes are converted to pore nodes. However, the benefit it brings is tremendous. Due to the disappearance of a part of the porous structure, the local mass transport resistance is significantly decreased. Besides, local diffusivity is also increased due to the rise in local pore size. Therefore, Fig. 8(b) shows the proposed structure can greatly reduce the transport loss, and a relatively larger radius can better enhance the transport. Therefore, the performance is greatly improved, as presented in Figs. 8(c) and 3 shows the current density increases by up to 50%, which also indicates that there is still huge potential for CCL to develop with the help of novel structural design.
Acknowledgments This work is supported by the National Key Research and Development Program of China (2018YFB0105505), the China-UK International Cooperation and Exchange Project (Newton Advanced Fellowship) jointly supported by the National Natural Science Foundation of China (Grant No. 51861130359) and the UK Royal Society (Grant No. NAF/R1/180146), and the National Natural Science Foundation of Tianjin (China) for Distinguished Young Scholars (Grant No. 18JCJQJC46700).
5. Conclusion In this work, the effects of Pt/C ratio, I/C ratio, and structure design are studied in terms of oxygen concentration distribution, reactive area, and reaction rate. The reactive transport processes in the cathode catalyst layer are simulated with LBM based pore-scale model. The cathode catalyst layers under various combinations of parameters and structure designs are elaborately reconstructed with the consideration of carbon, Pt, ionomer, and pores. Developed algorithms are employed to ensure the conductivity and relatively uniform distributions of these ingredients. The results indicate that a higher Pt/C ratio can significantly enhance mass transport and therefore improve the performance. Ionomer plays an ambivalent role as it’s essential for activating Pt sites while against the mass transport. However, it’s found that the performance benefits more from a larger active catalyst area as high I/C ratio yields better performance. Therefore, the active catalyst area should be prioritized, and the optimal I/C ratio should be just enough to activate most of the catalyst area. Finally, an ideal structure design is
References [1] Wang Y, Chen KS, Mishler J, Cho SC, Adroher XC. A review of polymer electrolyte membrane fuel cells: technology, applications, and needs on fundamental research. Appl Energy 2011;88(4):981–1007. [2] Huo S, Jiao K, Park JW. On the water transport behavior and phase transition mechanisms in cold start operation of PEM fuel cell. Appl Energy 2019;233:776–88. [3] Chang Y, Qin Y, Yin Y, Zhang J, Li X. Humidification strategy for polymer electrolyte membrane fuel cells–a review. Appl Energy 2018;230:643–62. [4] Wu HW. A review of recent development: transport and performance modeling of PEM fuel cells. Appl Energy 2016;165:81–106. [5] Perng SW, Wu HW. Effect of the prominent catalyst layer surface on reactant gas transport and cell performance at the cathodic side of a PEMFC. Appl Energy 2010;87(4):1386–99. [6] Zamel N. The catalyst layer and its dimensionality – a look into its ingredients and
9
Applied Energy 253 (2019) 113561
Y. Hou, et al.
[23] Lange KJ, Sui PC, Djilali N. Pore scale modeling of a proton exchange membrane fuel cell catalyst layer: effects of water vapor and temperature. J Power Sources 2011;196(6):3195–203. [24] Chen L, Wu G, Holby EF, Zelenay P, Tao WQ, Kang Q. Lattice Boltzmann pore-scale investigation of coupled physical-electrochemical processes in C/Pt and non-precious metal cathode catalyst layers in proton exchange membrane fuel cells. Electrochim Acta 2015;158:175–86. [25] Salomov UR, Chiavazzo E, Asinari P. Pore-scale modeling of fluid flow through gas diffusion and catalyst layers for high temperature proton exchange membrane (HTPEM) fuel cells. Comput Math Appl 2014;67(2):393–411. [26] Kim KH, Lee KY, Kim HJ, Cho E, Lee SY, Lim TH, et al. The effects of Nafion® ionomer content in PEMFC MEAs prepared by a catalyst-coated membrane (CCM) spraying method. Int J Hydrogen Energy 2010;35(5):2119–26. [27] Siddique NA, Liu F. Process based reconstruction and simulation of a three-dimensional fuel cell catalyst layer. Electrochim Acta 2010;55(19):5357–66. [28] Ishikawa H, Sugawara Y, Inoue G, Kawase M. Effects of Pt and ionomer ratios on the structure of catalyst layer: a theoretical model for polymer electrolyte fuel cells. J Power Sources 2018;374:196–204. [29] Yu Z, Carter RN. Measurement of effective oxygen diffusivity in electrodes for proton exchange membrane fuel cells. J Power Sources 2010;195(4):1079–84. [30] Zhang J, Yang W, Xu L, Wang Y. Simulation of the catalyst layer in PEMFC based on a novel two-phase lattice model. Electrochim Acta 2011;56(20):6912–8. [31] Bernardi DM, Verbrugge MW. A mathematical model of the solid-polymer-electrolyte fuel cell. J Electrochem Soc 1992;139(9):2477–91. [32] Hou Y, Deng H, Du Q, Jiao K. Multi-component multi-phase lattice Boltzmann modeling of droplet coalescence in flow channel of fuel cell. J Power Sources 2018;393:83–91. [33] Jiang ZY, Qu ZG, Zhou L, Tao WQ. A microscopic investigation of ion and electron transport in lithium-ion battery porous electrodes using the lattice Boltzmann method. Appl Energy 2017;194:530–9. [34] Li X, Ma T, Liu J, Zhang H, Wang Q. Pore-scale investigation of gravity effects on phase change heat transfer characteristics using lattice Boltzmann method. Appl Energy 2018;222:92–103. [35] Chen L, Kang Q, Mu Y, He YL, Tao WQ. A critical review of the pseudopotential multiphase lattice Boltzmann model: methods and applications. Int J Heat Mass Transf 2014;76:210–36. [36] Aidun CK, Clausen JR. Lattice-Boltzmann method for complex flows. Annu Rev Fluid Mech 2010;42:439–72. [37] Li Q, Luo KH, Kang QJ, Chen Q, Liu Q. Lattice Boltzmann methods for multiphase flow and phase-change heat transfer. Prog Energy Combust Sci 2016;52:62–105. [38] Bruggeman VDAG. Berechnung verschiedener physikalischer Konstanten von heterogenen Substanzen. I. Dielektrizitätskonstanten und Leitfähigkeiten der Mischkörper aus isotropen Substanzen. Ann Phys 1935;416(7):636–64.
how to characterize their effects. J Power Sources 2016;309:141–59. [7] Lee D, Hwang S. Effect of loading and distributions of Nafion ionomer in the catalyst layer for PEMFCs. Int J Hydrogen Energy 2008;33(11):2790–4. [8] Ozden A, Shahgaldi S, Li X, Hamdullahpur F. The impact of ionomer type on the morphological and microstructural degradations of proton exchange membrane fuel cell electrodes under freeze-thaw cycles. Appl Energy 2019;238:1048–59. [9] Shahgaldi S, Alaefour I, Li X. Impact of manufacturing processes on proton exchange membrane fuel cell performance. Appl Energy 2018;225:1022–32. [10] Shahgaldi S, Ozden A, Li X, Hamdullahpur F. Cathode catalyst layer design with gradients of ionomer distribution for proton exchange membrane fuel cells. Energy Convers Manage 2018;171:1476–86. [11] Zhao J, Ozden A, Shahgaldi S, Alaefour IE, Li X. Effect of Pt loading and catalyst type on the pore structure of porous electrodes in polymer electrolyte membrane (PEM) fuel cells. Energy 2018;150:69–76. [12] Zhao J, Shahgaldi S, Alaefour I, Xu Q, Li X. Gas permeability of catalyzed electrodes in polymer electrolyte membrane fuel cells. Appl Energy 2018;209:203–10. [13] Jiao K, Zhou Y, Du Q, Yin Y, Yu S, Li X. Numerical simulations of carbon monoxide poisoning in high temperature proton exchange membrane fuel cells with various flow channel designs. Appl Energy 2013;104:21–41. [14] Jiao K, Bachman J, Zhou Y, Park JW. Effect of induced cross flow on flow pattern and performance of proton exchange membrane fuel cell. Appl Energy 2014;115:75–82. [15] Li W, Zhang Q, Wang C, Yan X, Shen S, Xia G, et al. Experimental and numerical analysis of a three-dimensional flow field for PEMFCs. Appl Energy 2017;195:278–88. [16] Ferreira RB, Falcão DS, Oliveira VB, Pinto A. 1D+ 3D two-phase flow numerical model of a proton exchange membrane fuel cell. Appl Energy 2017;203:474–95. [17] Rao Z, Wang S, Peng F. Self diffusion of the nano-encapsulated phase change materials: a molecular dynamics study. Appl Energy 2012;100:303–8. [18] Andersson M, Yuan J, Sundén B. Review on modeling development for multiscale chemical reactions coupled transport phenomena in solid oxide fuel cells. Appl Energy 2010;87(5):1461–76. [19] Chen L, Zhang R, He P, Kang Q, He YL, Tao WQ. Nanoscale simulation of local gas transport in catalyst layers of proton exchange membrane fuel cells. J Power Sources 2018;400:114–25. [20] Wang G, Mukherjee PP, Wang CY. Direct numerical simulation (DNS) modeling of PEFC electrodes: Part I. Regular microstructure. Electrochim Acta 2006;51(15):3139–50. [21] Wang G, Mukherjee PP, Wang CY. Optimization of polymer electrolyte fuel cell cathode catalyst layers via direct numerical simulation modeling. Electrochim Acta 2007;52(22):6367–77. [22] Lange KJ, Sui PC, Djilali N. Pore scale simulation of transport and electrochemical reactions in reconstructed PEMFC catalyst layers. J Electrochem Soc 2010;157(10):B1434–42.
10