投稿論文
投稿
論文

Paper

Analysis of Water Cluster Structures and Mass Transport
in Polymer Electrolyte Membranes Containing Cerium
Ions Using Molecular Dynamics Simulations
Kyohei Ishikawa＊１，２, Takuya Mabuchi＊２，３，Takashi Tokumasu＊２＊１

Graduate School of Engineering, Tohoku University,＊２＊３

Institute of Fluid Science, Tohoku University,

The Frontier Research Institute for Interdiscipliary Sciences, Tohoku University

Abstract：We analyzed the effect of cerium ions, which are added to improve the chemical durability of a polymer
electrolyte membrane（PEM）, on the membrane structure and proton transport properties using molecular
dynamics simulations. We found that at low water content the water cluster connectivity is enhanced by theaddition of a small amount of cerium ions（0 . 6 - 1 . 8 wt%）that attract the surrounding water molecules. Hence,
better proton diffusion pathways are formed, leading to an increase in the proton diffusion coefficient. At high watercontent, however, because the water clusters are well connected, regardless of the cerium ion content in the PEM,
the presence of cerium ions hinders proton transport, resulting in a decrease in the proton diffusion coefficient.

Key Words ：Polymer electrolyte fuel cell, Molecular dynamics, Proton transport, Cerium ion, Durability.

１．Introduction

by the molecular motion of hydronium ions. In the Grotthuss mechanism, protons move by the connecting andPolymer electrolyte fuel cells（PEFCs）are highly re-disconnecting of the O-H bond between the hydronium

garded as a promising technology for stationary power

ion and the water molecule. Protons diffuse much fastergeneration and automobiles. To expand the wider use

by the Grotthuss mechanism than by the vehicle mecha-of PEFCs, improvements in the power density of unit

nism, but the effect of Grotthuss mechanism on proton１）

cells are necessary . Proton transport in the membrane

conductivity is reduced without an abundance of waterelectrode assembly（MEA）is closely associated withmolecules around the hydronium ions. The proton trans-achieving improved power density. Hence, clarifying theport properties in a PEM are greatly influenced by thecorrelation between the structure of the membrane andstructure of the water clusters３）−６）. One of the factors

proton transport at the molecular level is important. In

affecting the structure of a water cluster is the waterparticular, when considering proton conductivity, it is thecontent of the PEM. Experimental and computational

polymer electrolyte membranes（PEMs）that require op-

studies have reported that an increase in the water con-２）

timization . A perfluorosulfonic acid membrane having atent leads to improvement in the connectivity of clusters

sulfonic acid group（-SO３H）, such as Nafion, Flemion, orand an increase in the diffusivity of protons in a PEM

Aciplex, is mainly used as a PEM. In PEMs, high proton

due to the increased diffusivity by both transport mecha-

conductivity is required to improve the power genera-

nisms.

tion efficiency. Protons are transported in a PEM by two

An improvement in the durability of PEFCs is clearly

different mechanisms: i.e., vehicle and Grotthuss mecha-

required, therefore, it is important to clarify the causes

nisms. A proton typically exists as a hydronium ion,of deterioration and take the necessary remedial action７）.

which is the smallest stable hydrated form of an excess

Chemical degradation of PEMs is considered to be one ofproton in water, upon combination with a water mol-the causes of deterioration8）. An hydroxyl radical（OH·）

ecule. In the vehicle mechanism, protons are transported

is generated by the decomposition of hydrogen perox-

80

燃料電池 Vol.20 No.2 2020

ide, generated by an oxygen reduction reaction during

es in water content and cerium ion addition on the water

PEFC operation. This OH has high reactivity and it will

structure, ion distribution, and mass transport properties

degrade the polymer membrane９）−11）. Based on this deg-in the membrane.

radation mechanism, the method of adding a substancethat inactivates OHs before they react with the polymer（i.e., a radical scavenger）has been put into practical use.
One of the most useful radical scavengers is the cerium
ion

12）
−14）

２．Simulation methods
We constructed a Nafion bulk membrane system. However, previous studies have reported that

containing cerium ions and analyzed the effects of thecerium ions in a PEM migrate during operation and thewater content and cerium ion addition rate（CAR）on

distribution becomes nonuniform. At points where there

the structure, molecular distribution, and mass transport

are only a few cerium ions in the membrane, the effect asproperties of the membrane using MD simulation. Fig.

15）
，16）

.

1 shows the chemical structure of Nafion model with

Therefore, it is important to elucidate the transport

the equivalent weight of 1146 used in the present study.

mechanism of the cerium ion in a PEM. Furthermore,

The all-atom DREIDING force field was adopted forthe cerium ion is a trivalent ion, hence its large charge

the molecular interaction of Nafion, 19） and the particle

probably affects proton transport and water distribution

mesh Ewald method was used to calculate the Coulomb

17）

in the membrane . It is very difficult to investigate thestructure in a PEM experimentally and the mechanismthat determines the relationship between structure andforce. 20）
7 -CF-CF 2 ] 10 -[(CF 2 -CF 2）―

a radical scavenger is lost and degradation occurs

transport in a PEM has not yet been elucidated. There-

O-CF 2 -CF-CF 3―

fore, analysis by numerical simulation is considered effective.
Mabuchi and Tokumasu used molecular dynamics（MD）simulations to analyze proton transport properties

O-CF 2 -CF 2 -SO 3 -H+
Fig. １

Chemical structure of the Nafion polymer.

in a Nafion bulk membrane and elucidated the relation-

The cutoff distance of the Lennard-Jones and Coulomb

ship between the cluster structure and the contributions

potentials was set at r c = 12 Å. The anharmonic two-of the vehicle mechanism and the Grotthuss mechanism

state empirical valence bond（aTS-EVB）model was. Kawai et al.

used for water molecules and hydronium ions 21）. Using

used MD simulations to analyze the structure and proton

this model, we can perform simulations incorporatingtransport in a Nafion bulk membrane containing ferrous

the Grotthuss mechanism（one of the proton conduction

ions and elucidated that clusters will bind and proton

mechanisms, mentioned earlier）. The positive charge ofdiffusivity probably improve when a certain amount of

the hydronium ion and the cerium ion was balanced within electro-osmosis and proton transport

４）
，
６）

17）

ferrous ions is added at low water contents . Endoh et

the negative charge of the sulfonic acid group in the sys-al. experimentally clarified the increase in durability oftem. The number of water molecules was determined bya membrane upon the addition of cerium ions under low

the water content, λ, as shown in equation（1）,temperature and low humidity conditions

14）. Wong etal. developed a comprehensive chemical decomposition

λ =（N H 2 O+N H 3 O+）/N SO 3 -,

（1）

model for the ceria-added MEA. They indicated that

whereN H 2 O, N H 3 O+, and N SO 3 - represent the numbers of

the cerium ions will be transported to the cathode cata-

each of the molecules. The water content was set at λlyst layer side at low cell voltages and the degradation= 3 , 6 , 9 , and 12 , which correspond to the water content

18）

rate of the membrane probably increases . However,in the bulk Nafion membrane at approximately RH =molecular analysis of water cluster structure and mass

20 %, 60 %, 80 %, and 85 %, respectively 22）. The CAR wastransport of PEMs containing cerium ions is insufficient.

determined to correspond to 0 - 15 % of the total number

In this study, we modeled PEMs containing cerium ionsof sulfonic acid groups in the membrane. Three-dimen-

using MD simulations and analyzed the effects of chang-sional periodic boundary conditions were applied to the

燃料電池 Vol.20 No.2 2020

81

投稿論文
投稿
論文

Paper

simulation system. This simulation system is shown in

was judged to have reached the most stable state.

Fig. 2 The construction procedure of the simulation sys-[５] NPT ensemble at 350 K and 0 . 1 MPa for 100 ps.

tem is as follows. First, Nafion chains, water molecules,[６] NVT ensemble at 350 K for 1500 ps.

hydronium ions, and cerium ions were randomly placed

After annealing, the production run was performed forin a simulation system（x × y × z = 200 × 200 × 200

24 ns using the NVT ensemble. The time step was set

3

Å ）.

at 1 fs and the sampling interval was 10 , 000 steps. The
simulation was performed once under each condition,and two times additional simulations were performed forthe λ = 3 with different diffusion coefficient trends. The
temperature and density applied during the annealing
process are shown in Fig. 3 In the fourth and fifth NPT
ensembles of Step 4 , to reach equilibrium state the system was judged because there was almost no change indensity.

Fig.2

 napshot of the simulation system:λ = 12 and CAR.
S
= 10%.（White beads: Nafion, red beads:sulfur, green
beads:cerium ions, blue beads:water, orange beads:
hydronium ions.）

Table 1 shows the number of molecules in the system.
Next, annealing was performed. Temperature was conFig.3

trolled by the Nosé-Hoover method 23） and pressure was

Relaxation process for annealing.

24）

controlled by the Andersen method .
Table1

Number of molecules in the calculation system.

Water
content λNafionCe 3 +H2 O

H 3 O+

3

4

0-6

80 - 98

22 - 406

40-6

200 - 218

22 - 4094

0-6

320 - 338

22 - 40

12

4

0-6

440 - 458

22 - 40

３．Result and discussion
３．
１

Diffusion coefficient

The self-diffusion coefficient was calculated from themean square displacement（MSD）using Einstein’s equation（2）,,The annealing procedure was as follows（where thefollowing abbreviations apply: N is the number of mole-where〈（2）〉is the time average, t is the diffusion time,cules, V is the volume of the system, T is the tempera-r is the atomic position vector, and D is the self-diffusion

ture of the system, and P is the pressure of the system）:

coefficient. The self-diffusion coefficient was calculated[１]　Simulation system was compressed by the NPT en-

from the liner regime of MSD curves from 0 . 2 to 0 . 3

semble at 350 K and 1 MPa for 190 ps.

ns. The self-diffusion coefficient was converted to the[２] NVT ensemble at 800 K for 40 ps.

membrane resistance and compared with the experi-[３] NPT ensemble at 350 K and 0 . 1 MPa for 40 ps.

mental value 15）. The conductivity was obtained using the[４] Step [ 2 ] and step [ 3 ] were repeated five times.

Nernst-Einstein equation（3）and converted into a mem-

Since the densities of the system were almost thebrane resistance using equation（4）,same after the fourth and fifth steps, the system

82

燃料電池 Vol.20 No.2 2020

lar, a decrease in the diffusion coefficient is dominant,


（3）

,

compared to the decrease in the proton concentration,
at low water content. However, at λ = 3 and CAR = 5 ,



（4）

,

10 %, the increase in the membrane resistance is reduced
because the diffusion coefficient of proton increases,

where e is the elementary charge, n is the number

compared with when no cerium is added. Therefore, indensity, kB is the Boltzmann constant, and T is the tem-the next section（3 . 2）
, we will compare the structures

perature. Fig. 4 and Fig. 5 show the proton diffusion coef-

when λ = 3 and λ = 12 , and consider the causes of theficient and membrane resistance, respectively. The diffu-

increase in diffusivity.

sion coefficient value at CAR = 0 % and the decrease in

Fig.6 shows the self-diffusion coefficient of cerium

proton conductivity with increasing CAR at high waterion. The diffusion coefficient increases with increasing

content are well consistent with previous experimental

water content. This is in agreement with experimental

researches 22），25）. Therefore, we believe that the simulated

results 26）. There was no clear relationship between thesystem reproduces reasonable physical properties. The

CAR and the diffusion coefficient. The mass transport

diffusion coefficient increases with increasing water con-characteristics described above are probably greatly af-

tent regardless of the CAR. The membrane resistancefected by a change in the membrane structure, upon antends to increase as the CAR increases due to a decrease

increase in water content; therefore, in the next sectionin the diffusion coefficient of protons or a decrease in（3 . 2）, we will discuss these reasons, based on structural

15）

the proton concentration in the membrane . In particu-

analysis.

Fig.6
Fig.4

Self-diffusion coefficient of cerium ions.

Self-diffusion coefficient of protons.

３．
２Structural analysis

The dependence of the molecular structure in a PEM
on the water content and the CAR was analyzed using
the radial distribution function（RDF）. Because the membrane resistance, depending on the presence or absenceof cerium ions, differs between at low water content andat high water content, the analysis was performed focusing on the structural change inside the PEM with an increase in the water content. Fig. 7 shows the RDF gS-Ce（r ）between the sulfonic acid group and the cerium ionwith an increase in the water content when the CAR =10 %. The first peak emerges at r = 4 . 1 Å , the height of
Fig.5

Membrane resistance at various water contents.

which decreases as the water content increases. There-

燃料電池 Vol.20 No.2 2020

83

投稿論文

Fig.7

Paper

RDF of sulfur-cerium ion（CAR 10%）
.

Fig.9

RDF of sulfur-hydronium oxygen（λ = 3）.

（a）Fig.10

RDF of cerium ion-water oxygen（CAR 10%）.

（b）

sulfonic acid group and the hydronium ions when λ =3 . The first peak is observed at r = 4 . 1 Å in both thewater molecule and the hydronium ion, indicating thatthese are concentrated around the sulfonic acid groupas well as the distribution of cerium ion. In addition, thesecond peaks appear at r = 6 . 2 Å in the 8 （a）,（b）due tothe addition of cerium ions. This indicates that the ceriFig.8

RDF of sulfur-water oxygen:（a） when λ = 3 and（b）when λ = 12.

um ions that are around the sulfonic acid groups attract
water molecules, hence the water distribution changes.
The second peak is more apparent because the watercontent is lower. From the above, we observed that ce-

fore, the indications are that the distribution of cerium

rium ions, water molecules, and hydronium ions exist

ions in the membrane is concentrated around the sul-around the sulfonic acid group, and at low water content,

fonic acid group and spread out into the solvent, away

more water molecules are attracted by the addition of

from the sulfonic acid group, as the number of water

the cerium ions. This is probably related to the decrease

molecules increases. This also occurs, even when there isin membrane resistance at low water content due to thea change in the CAR.

small amount of cerium addition（described in the previ-

Fig. 8 （a）,（b）shows the RDF gS-Ow（r ） between the sul-

ous section）because the cerium ions would bind water

fonic acid group and the water molecules when λ = 3clusters. This point will be examined in the later clusterand λ = 12 . Fig. 9 shows the RDF gS-Oh（r ） between theanalysis（see Section 3 . 3）.

84

燃料電池 Vol.20 No.2 2020

Fig. 10 shows the RDF gCe-Ow（r ） between the cerium

are distributed throughout the whole membrane.   This

ion and the water molecules with an increase in the wa-change in the distribution of cerium ions is thought to be

ter content when CAR = 10 %. There is a peak at r = 2 . 7

related to the increase in the diffusion coefficient of ce-Å , and the position of the peak does not change withrium ions with increasing water content.

the water content. Figs 11 and 12 show the coordination
numbers of cerium ions to sulfonic acid groups and wa-

３．
３

Cluster analysis

ter molecules to cerium ions, respectively. As the water

Next, we analyzed the number of clusters and thecontent increases, the coordination number of the cerium

cluster size. In this calculation, a water cluster is defined

ion around the sulfonic acid group decreases, and the co-

as a water and hydronium ion aggregate in which theordination number of water molecules around the cerium

intermolecular distance of the oxygen atom of the waterion increases. Furthermore, the coordination number ofmolecule is within 3 . 4 Å. This distance of 3 . 4 Å is thewater molecules around the cerium ion is sufficiently

terminal value of the first peak in the RDF between oxy-large（5 . 5 - 8 . 2）, regardless of the water content. This

gen atoms of water molecules in the Nafion bulk mem-

implies that there are normally many water molecules

brane（see Fig. 13）. This value is consistent with thearound the cerium ion, indicating that the cerium ion

value used for cluster analysis in the previous study.６）In

exists inside the water cluster. From the above, we con-this study, the number of water molecules contained in asidered that the cerium ions are concentrated aroundcluster is defined as the cluster size.

the sulfonic acid group at low water content but as theFig. 14 shows the average number of clusters of eachwater content increases the cerium ions move away

CAR value at λ = 3 and 12 and Fig. 15 shows the aver-from the sulfonic acid group into the water cluster andage cluster size of each CAR value at λ = 3 and 12 . At
λ = 3 , the number of clusters decreases and the average cluster size increases with an increase in CAR. At λ= 12, however, no significant change was observed in thenumber of clusters and the average cluster size at each
CAR. Furthermore, upon comparing the results according to water content, we observed that the number ofclusters decreases and the average cluster size increases
as the water content increases.
From the above, we deduced that the cerium ionsgather around the surrounding water molecules by electric charges and bind clusters at low water content, as

Fig.11

 oordination number of cerium ions around sulfonicCgroups.

Fig.12

 oordination number of water oxygens around ceriumC
ions.

shown in Fig. 16 （a）. At low water content, CAR = 15 %,Fig.13

燃料電池 Vol.20 No.2 2020

RDF of water oxygen-water oxygen.

85

投稿論文

Paper
（a）

（b）

Fig.14

Number of clusters at λ = 3 and 12.

Fig.16 chematic of the structure of the sulfonic groups andSwater clusters:（a）when λ = 3 and （b）when λ = 12.

rounding water molecules and the clusters are connected. Hence, the proton diffusion pathways were improved.
At the water content of 3 and CAR = 5 , 10 %, although
the proton number density decreases and the membrane
Fig. 15 Average cluster size at λ = 3 and 12.

resistance increases with the increase CAR, the change
is significantly small because the diffusion coefficient ofproton increased.  

the connectivity of clusters is enhanced, but the proton

At high water content, the diffusion coefficient ofdiffusivity is considered to have decreased because theprotons decreased due to the addition of cerium ions be-increase of Ce ion concentration increases the hindrance

cause water clusters are already well connected, and ce-of proton transport by the electric charge of Ce ion.

rium ions have negligible effect on the cluster structure.

Then, because the clusters have already connected suffi-

In this case, the positive charge of the cerium ion repels

ciently at high water content（Fig. 16 （b））
, we considered

the protons and prevents proton diffusion, resulting thethat the cerium ion has almost no effect on the cluster

increase in the membrane resistance.
The diffusion coefficient of cerium ions increases withstructure.

increasing water content because the water molecules

4．Conclusion

gathering around the sulfonic acid group shield thecharge as the water content increases and the cerium

In this study, we analyzed in the effects of the cerium

ions diffuse into the water cluster and move freely.

ion addition rate and the water content on the protonand cerium ion transport properties and the structuralproperties in a PEM using MD simulations that takes

Acknowledgments
This research was supported by the New Energy
and Industrial Technology Development Organization

into account the Grotthuss mechanism.
At low water content, upon the addition of a small（NEDO）of Japan. We used the integrated super compu-

amount of cerium ions in the membrane, an improve-

tation system at the Institute of Fluid Science of Tohoku

ment in the proton diffusion coefficient was found. The

University.

cerium ions around the sulfonic acid group attract sur-86

燃料電池 Vol.20 No.2 2020

12）E. Endoh: Development of highly durable PFSA

References
１）NEDO: Fuel cell / hydrogen technology develop-

membrane and MEA for PEMFC under high tem-ment roadmap Detailed version（Fuel cell field）,

perature and low humidity conditions, ECS Trans .,（2017）.

16, 1229-1240,（2008）.

２）D. D. Borges, G. Gebel, A. A. Franco, et al.: Mor-

13）P. Trogadas, J. Parrondo and V. Ramani: Degrada-

phology of supported polymer electrolyte ultrathintion Mitigation in Polymer Electrolyte Membranes

films: A numerical study, J. Phys. Chem. C , 119,

Using Cerium Oxide as a Regenerative Free-

1201-1216,（2015）.

Radical Scavenger, Electrochem. Solid-State Lett .,３）T. Mabuchi and T. Tokumasu: Effects of water

11, B113,（2008）.

nanochannel diameter on proton transport in14）E. Endoh, N. Onoda, Y. Kaneko, et al.: Membrane

proton-exchange membranes, J. Polym. Sci. Part B

Degradation Mitigation of PEFC during Cold-

Polym. Phys ., 57, 867-878,（2019）.

Start Application of the Radical Quencher Ce3+-,４）T. Mabuchi and T. Tokumasu: Relationship between Proton Transport and Morphology of Perfluorosulfonic Acid Membranes: A Reactive Molecular Dynamics Approach, J. Phys. Chem. B , 122,

ECS Electrochem. Lett ., 2, F73-F75,（2013）.
15）Current status and issues of FCV fuel cells, NEDO

Issue Shar. Forum .
16）A. M. Baker, D. Torraco, E. J. Judge, et al.: Cerium
migration during PEM fuel cell assembly and op-

5922-5932,（2018）.
５）H. L. Yeager and A. Steck: Cation and water diffu-

eration, ECS Trans ., 69, 1009-1015,（2015）.

sion in nafion ion exchange membranes: Influence17）K. Kawai, T. Mabuchi and T. Tokumasu: Molecu-

of polymer structure, J. Electrochem. Soc ., 128,

lar Dynamics Analysis of Proton Diffusivity in Hy-

1880-1884,（1981）.

drated Nafion Membranes Contaminated with Fer-

６）T. Mabuchi and T. Tokumasu: Dependence ofelectroosmosis on polymer structure in proton exchange membranes, Mech. Eng. J ., 4, 17-00054-17-

rous Ions, Macromol. Theory Simulations , 1900047,
1-8,（2019）.
18）K. H. Wong and E. Kjeang: In-Situ Modeling ofChemical Membrane Degradation and Mitigation

00054,（2017）.
７）M. P. Rodgers, L. J. Bonville, H. R. Kunz, et al.:Fuel cell perfluorinated sulfonic acid membrane
degradation correlating accelerated stress testingin Ceria-Supported Fuel Cells, J. Electrochem.

Soc ., 164, F1179-F1186,（2017）.
19）T. Mabuchi and T. Tokumasu: Effect of boundand lifetime, Chem. Rev ., 112, 6075-6103,（2012）.

state of water on hydronium ion mobility in hy-８）C. Lim, L. Ghassemzadeh, F. Van Hove, et al.:drated Nafion using molecular dynamics simula-

Membrane degradation during combined chemicaltions, J. Chem. Phys ., 141,（2014）.

and mechanical accelerated stress testing of poly-20）U. Essmann, L. Perera, M. L. Berkowitz, et al.: A

mer electrolyte fuel cells, J. Power Sources , 257,

smooth particle mesh Ewald method, J. Chem.

102-110,（2014）.

Phys ., 103, 8577-8593,（1995）.

９）M. Inaba, T. Kinumoto, M. Kiriake, et al.: Gas21）T. Mabuchi, A. Fukushima and T. Tokumasu: A

crossover and membrane degradation in polymer

modified two-state empirical valence bond model

electrolyte fuel cells, Electrochim. Acta , 51, 5746-

for proton transport in aqueous solutions, J. Chem.

5753,（2006）.

Phys ., 143,（2015）.

10）J. Healy, C. Hayden, T. Xie, et al.: Aspects of the22）S. Ochi, O. Kamishima, J. Mizusaki, et al.: Investi-

chemical degradation of PFSA ionomers used ingation of proton diffusion in Nafion®117 membrane

PEM fuel cellsx, Fuel Cells , 5, 302-308,（2005）.
11）L. Gubler, S. M. Dockheer and W. H. Koppenol:
Radical（HȮ, Ḣ and HOȮ）Formation and ionomer degradation in polymer electrolyte fuel cells, J.

Electrochem. Soc ., 158,（2011）.

by electrical conductivity and NMR, Solid State

Ionics , 180, 580-584,（2009）.
23）W. G. Hoover: Constant-pressure equations of motion, Phys. Rev. A , 34, 2499-2500,（1986）.
24）H. C. Andersen: Molecular dynamics simulations

燃料電池 Vol.20 No.2 2020

87

投稿論文

Paper

at constant pressure and/or temperature, J. Chem.

Phys ., 72, 2384-2393,（1980）.

F136,（2019）.
26）S. Kumaraguru: Durable High Power Membrane25）K. H. Wong and E. Kjeang: Simulation of Perfor-

Electrode Assembly with Low Pt Loading, Energy.

mance Tradeoffs in Ceria Supported Polymer Elec-gov Hydrog. Fuel Cells Progr. Annu. Merit Rev ., 1,

trolyte Fuel Cells, J. Electrochem. Soc ., 166, F128-

1-9,（2018）.

88

燃料電池 Vol.20 No.2 2020