Form Follows Function: a brief survey on the morphology of generic solvated

polymer electrolyte membrane

Gillian Kwan1 and Lauren Clark-Johnson2
1)
Department of Chemistry, Mount Holyoke College
2)
Department of Physics, Mount Holyoke College
To improve the performance of proton exchange membrane fuel cells (PEMFCs), the proton conduction
mechanism needs to be thoroughly understood. The proton conduction mechanism is dictated by the structure
of the membrane, therefore the elucidation of operating form of the polymer electrolyte membrane (PEM) is
the crucial first step to advancing fuel cell technology. A series of generic solvated bead-spring polymer were
simulated to explore the general hydrated morphology of polymer electrolyte membrane used in fuel cells.
Simulations were run using classical molecular dynamics program LAMMPS and its micelle tool. Structural
and potential parameters were varied in generic PEM systems to see their impact on phase segregations upon
coarse grained amphiphilic polymer molecules interaction with generic hydrophilic solvent. It was found
that an intermediate to high concentration of surfactant stabilizes the generic polymer membrane, as well
as longer surfactant chain length. These two factors help explain the popularity of Nafion, a sulfonated
polytetrafluoroethylene based polymer, as the standard in PEMs.

I. INTRODUCTION known as polymer electrolyte membrane fuel cells be-

cause the electrolyte layer consists of a water-based,
A Lamborghini Aventador LP 700-4 can accelerate acidic polymer membrane. The polymer material is com-
from 0 to 100 km/h in just 2.9 seconds thanks to its posed of hydrophilic segments that conduct ions, and hy-
V12 petrol engine, which has a maximum power output drophobic segments that provide mechanical stability in
of 515 kW at 8250 rpm.1 The internal combustion en- solution. PEMs operate at relatively low temperatures
gine liberates the energy stored in fuel to do mechanical (below 100◦ C) because proton conduction occurs via sol-
work, thereby propelling the car to move. In spite of such vation by water.3 SOFCs on the other hand operate be-
impressive speed, the V12 engine is subjected to Carnot tween 800 and 1000◦ C for the ceramic solid electrolytes
limitation as heat requires a temperature gradient, thus do not become electrically and ionically active until such
limiting the maximum efficiency of the engine by the fac- high temperature range.4 These two types of fuel cells
tor T1T−T 2
, which translates to about 15-20% efficiency are of particular interest because even though they em-
1
in real life.2 In addition, the car emits carbon dioxide at ploy similar mechanism to generate power, they operate
a rate of 370 g/km.1 In wake of energy crisis and global at the opposite ends of the fuel cell operation tempera-
warming, there is not enough fossil fuels to spare or room ture scale. The optimal operating temperature obviously
for more green house gases. depends on the characteristics of the electrodes and elec-
Fuel cell technology is a promising development for a trolyte materials, which in turn are derived from their
safe, environmental-friendly sustainable energy source. molecular and atomic structures. The principle “form
Unlike a combustion engine, a fuel cell converts chem- follows function” was first coined by American architect
ical energy into electrical energy directly without pollu- Louis Sullivan and became the defining pillar of mod-
tants as by-products and its efficiency is not hindered by ern architecture, but it is so universal that it extends to
Carnot limitation. The fuel cell uses an external supply chemical engineering as well.
of chemical energy, usually hydrogen and oxygen, and can It is the pervading law of all things organic
operate indefinitely provided there is continuous supply and inorganic,
of hydrogen and its components are functional. All fuel of all things physical and metaphysical,
cells consists of a layer of electrolyte sandwiched between of all things human and all things superhu-
an anode and a cathode; the reaction taking place inside man,
fuel cell is summarized in Equation 1. of all true manifestations of the head, of the
heart, of the soul,
Anode − 4H+ + 4e−
2H2 (g) → that the life is recognizable in its expression,
Cathode O2 (g) + 4H+ + 4e− → − 2H2 O(`) that form ever follows function.5
Overall 2H2 (g) + O2 (g) →
− 2H2 O(`) (1) Following this principle, before investigating proton con-
duction mechanism in fuel cell electrolyte, we first need
Fuel cells can be engineered to sizes ranging from small to determine the “form” that allows such function. Such
handheld devices to power plants by varying the mate- task is relatively straightforward for SOFCs, as ceramic
rial used for the electrolyte. The focus of this paper solid electrolytes are normally mechanically robust and
are proton-conducting fuel cells, namely proton-exchange they retain their main framework in their operating en-
membrane fuel cells (PEMFCs) and solid oxide fuel cells vironment. However it is much more complex in PEM-
(SOFCs). Proton-exchange membrane fuel cells are also FCs because the membrane form that facilitates pro-

Typeset by REVTEXand AIP

 2

ton conduction is very dynamic – so dynamic that there conduction mechanism can be investigated directly with-
exists a cornucopia of theories on the structure of hy- out prior structure determination. As a side comparison
drated perfluorosulfonic acid (PFSA) membrane, the pre- with PEM simulation, a few of the proton conduction
dominant class of polyelectrolyte material in PEMFCs, pathways of yttrium-doped BaZrO3 were inspected.
alone. These structural models arise from the complex
phase behavior of PFSA-water system, and they include
spherical water clusters, parallel cylindrical water chan- II. MOLECULAR MODELS AND METHODS
nel networks, layered structures, polymer-bundles, and
polymer-ribbon models among many others.6–8 In con- A. Simulation of PEMs
sidering a generic PEMFC, given the amphiphilic nature
of the polymer electrolyte, there is no surprise that phase To understand the proton conduction mechanism in
separation occurs when PEM is solvated. At high surfac- polymer electrolyte membrane, simulation study requires
tant concentration (as in PEM), the polymer molecules elucidation of the electrolyte structure in its operating
self-assemble into micelles, aggregates in which the hy- environment – the solvated membrane. Instead of mod-
drophobic regions are sequestered from solvent by the eling Nafion, the archetypal PSFA membrane manufac-
hydrophilic segments.9 This nanoscale phase segregation tured by DuPont and the most popular material in PEM-
resulting from surfactant interaction with solvent obvi- FCs, a generic surfactant is used so that an unbiased
ously alters the morphology of the membrane, which in exploration of PEM structure is possible. For faster
turn determines the transport properties.10 computational time, all simulation models were two-
Previous simulation studies have attempted to elu- dimensional; such simplification was acceptable for mod-
cidate the morphology of solvated PEMs on various eling a thin PEM in a stacked membrane electrode assem-
levels of resolution, ranging from microscale atomistic- bly. In order to function in a PEMFC, the membrane
quantum mechanical methods to mesoscale coarse- material must meet the following criteria: i) high in-
grained methods.11 While atomistic simulations provide trinsic proton conductivity; ii) long-term chemical, ther-
detailed information on structures, they only apply to mal and mechanical stability; iii) impermeable to gases,
local structures; extrapolating atomistic results to larger methanol, and contaminating ions.15,16 In this study, cri-
scale simulations is not sensible as each scale is governed terion i was assessed by the shape and connectivity of
by different laws of physics.12 Fortunately, the different water channels formed in the membrane. Criteria ii and
scales can be bridged by using a multi-scale approach, iii were not evaluated as they were out of the scope of
which normally entails starting from atomistic level and this study.
from that deriving a coarse grained model.13 An accu- The model The simulation box contained an initial
rate multi-scale modeling requires meticulous parameter- square lattice of solvents of the size 30 × 30Å, which was
ization of macroscopic properties, as degrees of freedom then divided into smaller grids of length rlattice in which
increase in scaling-up.14 Fastidiousness aside, the multi- molecules only interacted with other molecules that were
scale approach seems to be the most appropriate method in the same or neighboring grids. rlattice was the square
in investigating the morphology of hydrated PEMs. How- root of the inverse of the number density of the sys-
ever on the other hand, regardless of the precise morphol- tem (ρ∗ ). The actual box size was the product of ini-
ogy of the solvated dynamic membrane and the actual tial square lattice of solvents and rlattice . The surfac-
proton conduction mechanism, the fact that proton con- tants were introduced to the system by replacing a cer-
ductivity depends on membrane hydration in all types tain number of solvents specified in the definition file. In
of PEMs that operate in temperature range of 60-90◦ C this coarse grained model, the generic surfactant-solvent
establishes firmly that proton transport occurs within a system was presented as a sequence of two-dimensional
water network.2 In other words, atomistic scale details beads. There were four types of beads: solvent (1), head
notwithstanding, on the macroscopic level all PEMFCs (2), tail (3) and tail (4). The non-bonded interactions
contain a water channel network of some sort as means among them are summarized as follows: (a) solvent-
of proton conduction. Since the main goal of this pa- head, full-size and long-range; (b) tail-tail, size-averaged
per was to survey the structures that facilitate proton and long-range; (c) solvent-tail, full-size and repulsive;
conduction in PEMFCs in general, atomistic simulations (d) head-tail, size-averaged and repulsive. These sites in-
were both unnecessary and unfeasible for generic PEMs. teracted with one another via truncated Lennard-Jones
For these reasons, the investigation of proton conduc- potential. Lennard-Jones potential is chosen because it
tion in PEMFCs was limited to the morphological aspect is the standard pairwise interaction potential for generic
only, specifically mesoscale simulations of water channel systems; the potential is truncated because it is custom-
topology. Such limitation emerged due to generic na- ary to establish a cut-off distance to conserve computa-
ture of the system; quantitative characterizations such tional resources.17 The expression of Lennard-Jones 12-6
as conductivity within membrane and proton conduc- potential is given as
tion mechanism are impossible without explicit molecular   
species. Unlike PEMFCs, there is no contention about σ 12  σ 6
V (r) = 4ε − r < rc (2)
the structure of electrolytes in SOFCs, hence the proton r r
 3

TABLE I: S0 Lennard-Jones potential parameters

Interaction Sites ε σ Cut-off

1-1 1.0 1.0 2.5
solvent-head 2-2 1.0 1.0 2.5
1-2 1.0 1.0 2.5
3-3 1.0 0.75 2.5
tail-tail 4-4 1.0 0.50 2.5
3-4 1.0 0.67 2.5
1-3 1.0 1.0 1.12246
solvent-tail
1-4 1.0 1.0 1.12246
2-3 1.0 0.88 1.12246
head-tail
2-4 1.0 0.75 1.12246

where r is the distance between sites, rc is the cut-off dis-

tance, ε is the well depth, and σ is the zero-crossing dis-
tance for the potential. Physically ε and σ are constants
characteristic of each Lennard-Jones system, in which ε
denotes the magnitude of attraction between two sites
and σ represents the maximum distance at which non-
bonded interaction remains zero between two sites. The FIG. 1: A simplified scheme illustrating interactions in
values of ε and σ were varied to see their impact on the MD simulation
geometry of micelles. Although the Lennard-Jones po-
tentials defined for solvent molecules do not physically
correspond to water, for the purpose of this study the tems were generated from running the definition file in
solvent is regarded as interchangeable with water. The the Fortran-based program. The data containing the co-
system that served as the starting point was referred to as ordinates were then read by a LAMMPS input script
System 0 (S0). The Lennard-Jones potential parameters for atom definition; but because the coordinates referred
for S0 are listed in Table I. S0 had a number density of to the positions of polymer beads rather than those of
0.7, 150 surfactants, and surfactant bond length of 0.75 atoms, FENE (finite extensible nonlinear elastic) poten-
in the units of σ (where 1 σ is equivalent to 1Å as per tial was used to account for weighted pairwise energy and
program definition). Given the polymeric nature of sur- force contributions from bonded atoms within each bead
factants, concentration is expressed in terms of weight molecule.20 The equation for FENE is
throughout this paper, with each bead weighing 1 mass "
unit.  2 #   
2 r σ 12  σ 6
Molecular Dynamics theory and computational details E = −0.5KR0 ln 1 − + 4ε − +ε
R0 r r
Simulations were run using classical molecular dynamics
program LAMMPS (Large-scale Atomic/Molecular Mas- (4)
sively Parallel Simulator). Molecular Dynamics (MD)
is a statistical mechanics simulation technique in which where K is constant coefficient with the unit of en-
atomic trajectories within a system are calculated by ergy per squared distance, and R0 is the maximum ex-
integrating Newton’s equation of motion at each time tent of the FENE bond. Each bead was connected to
step.18 Newton’s second law is given as follows: one another via springs like in harmonic oscillators. A
schematic of the model is shown in Figure 1. Interac-
d2 ri (t) ∂V tions among beads were first computed using soft poten-
mi = Fi = − i = 1, ..., N (3) tial push-off at fixed NVE condition, whose formula is
dt2 ∂ri
 
where ri is the position of particle i. From Equation 3, πr
E = A 1 + cos r < rc (5)
it is clear that a function V (ri , ...rN ) representing the rc
potential energy of the system needs to be chosen prior
to simulation, as it is from this function that forces are This soft pairwise potential is also known as cosine po-
derived as gradients of potential with respect to rela- tential, in which A is a pre-factor. It is to clean up the
tive coordinates of particles.19 As discussed previously, molecular topology, such as separate overlapping beads,
Lennard-Jones 12-6 potential was the potential model of before the actual simulation.20 The soft potential was
choice. computed for 1000 time steps, then the main computa-
Simulation systems were constructed using LAMMPS tion with Lennard-Jones pairwise potential was run for
tool micelle2d.f. Coordinates of the two-dimensional sys- 60000 timesteps.
 4

B. Simulation of SOFCs TABLE II: SE0 Lennard-Jones potential parameters

The tortuosity and porosity of PEMs attests to how Interaction Sites ε σ Cut-off
forms follow functions and vice versa. As a side com- 1-1 2.0 1.0 2.5
parison, a brief survey of proton conduction pathway in solvent-head 2-2 2.0 1.0 2.5
yttrium-doped barium metazirconate (YBaZrO3 ), a solid 1-2 2.0 1.0 2.5
ceramic electrolyte whose functional form has already 3-3 2.0 0.75 2.5
been elucidated, was conducted. tail-tail 4-4 2.0 0.50 2.5
Kinetic Monte Carlo method and computational de-
3-4 2.0 0.67 2.5
tails Although molecular dynamics can be used to com-
pute the dynamic evolution of the proton conduction in 1-3 0.75 1.0 1.12246
solvent-tail
YBaZrO3 , the time steps required for accurate simulation 1-4 0.75 1.0 1.12246
are usually of the order of femtoseconds (1×10−15 s), and head-tail
2-3 0.75 0.88 1.12246
because proton conduction proceeds on a much longer 2-4 0.75 0.75 1.12246
time scale, MD simulation of this scenario would require
infelicitously long computational time. Kinetic Monte
Carlo (KMC) is an algorithm that bypasses integration inefficient. S1 has the same potential and box size as
of motion along all time steps by considering transition does S0 but twice the number of surfactants, resulting
between states in a system directly.21 in a system with 60% surfactants and 40% solvent. As
The YBaZrO3 system contains 96 proton binding seen in Figure 2b, doubling the number of surfactants
sites and 48 redundant sites. The proton can move from 150 to 300 has lead to narrower water channels,
by (1) rotation; (2) intra-octahedral transfer; (3) and but the aggregates are still appearing to be distributed
inter-octahedral transfer. KMC computation was run at unevenly throughout the membrane. Based on this ob-
1000K starting at binding site 35 for 100 ps. Long range servation, a series of investigation in varying ε and σ was
pathways in trajectory were located using find program. conducted to see their impact on clusters formation. Fig-
ure 3a shows a modified S0 in which the attractive and
repulsive interactions have been made more extreme (the
III. RESULTS parameters are listed in Table II). This modified poten-
tial system is referred to as SE0, where E stands for
A. PEMs
epsilon. The VMD snapshot in Figure 3a shows that the
more extreme interactions among sites result in voids.
SE1 and SE2 are defined by the same set of potential
parameters as SE0, except the number of surfactants has
been changed from 150 to 300 and 500 respectively. From
Figure 3, it appears that increased surfactant concentra-
tion stabilizes the surface voids. However the surfactant
concentration of 78.9% in SE2 is too high that the wa-
ter channels are blocked. From these preliminary runs, it
could be gathered that the desirable number of surfactant
present in a functional PEM should be high enough to
stabilize strong site interactions, but not so high that the
surfactant clusters interrupt water clusters connectivity.
(a) S0 (b) S1

FIG. 2: Snapshots from VMD of the MD simulations of

S0 and S1. Solvent is represented by cyan beads, head
pink, tail 1 purple, and tail 2 lime for all systems in this
paper.

A snapshot from VMD of S0 is shown in Figure 2a.

The actual box size has expanded to 35.8569 × 35.8569Å (a) SE0 (b) SE1 (c) SE2
after simulation and contains 1200 molecules in total,
of which 62.5% by weight is solvent. At this hydra- FIG. 3: Comparison
tion level, the surfactant concentration (37.5%) is too
low to delineate any water channels for more efficient Varying ε in SEs did not yield noticeable change in mi-
proton conduction. The broad and wide water networks celles organization, so another series of simulation tests
mean that protons could travel via myriad of pathways to were run (hence systems T). The T systems are derived
the cathode, which could be time-consuming and energy- from SE1: they all contain 1500 molecules with an initial
 5

TABLE III: T Systems Lennard-Jones potential

parameters

Interaction Sites ε σ
T1 T2 T3 T4 T1 T2 T3 T4
1-1
solvent-head 2-2 3.0 3.0 3.5 5.0 1.0 0.5 1.0 1.0
1-2
3-3 0.75
tail-tail 4-4 3.0 3.0 3.5 5.0 0.50
3-4 0.67
1-3
solvent-tail 0.5 0.5 0.5 0.25 1.0 0.5 1.2 1.2
1-4 (a) Solvent-head interactions
2-3 0.88
head-tail 0.25 0.25 0.25 0.2
2-4 0.75
The cut-off values of T systems are virtually the same
as the previous systems, except rc were changed from
1.12246 to 2.5 for solvent-tail interactions in T3 and T4
to accommodate for a shift in equilibrium distance due
to a larger σ.

(b) Solvent-tail interactions

FIG. 5: Lennard-Jones potential plots

(a) T1 (b) T2
is observed that even though the values of ε become more
extreme in successive systems, T3 and T4 still retain the
same semblance as that of T1. While there appears to
be slight dimpling of the membrane with more extreme
parameterization of ε, there is no significant change in
the number of lone polymer chains (monomers) in these
three systems. Table III suggests that it is the sub-unity
values of σ in T2 that engender the void-riddled struc-
ture. To better understand the peculiarity of T2, con-
sider the Lennard-Jones potential plots of solvent-head
(c) T3 (d) T4 interactions and solvent-tail interactions in T systems
in Figure 5. (The Lennard-Jones potentials in tail-tail
FIG. 4: Series T simulation snapshots interactions and head-tail interactions are omitted be-
cause the σ values do not vary across series.) Recalling
that rc = 2.5σ in all solvent-head interactions, Figure 5a
solvent square lattice of 30×30Å and a surfactant concen- shows that the repulsion in T2 occurs at shorter range
tration of 60%. The Lennard-Jones potential parameters due to a smaller value of σ. Consequently, despite that a
are cataloged in Table III, and the configurations of the diminished magnitude of ε in T2, the long-range poten-
systems after MD simulations are included in Figure 4. tial is more prominent, which explains the denser pack-
T1, T3 and T4 all exhibit leopard spots-like patterns in ing seen in Figure 4b. Similarly, the long-range potential
their micelle arrangements, where as T2 shows a radi- dominates in T2 in solvent-tail interactions as demon-
cally different topology. It is interesting to note that the strated in Figure 5b. The repulsive force between hy-
only potential parameter that is different in T1 and T2 drophobic tail sites and hydrophilic solvent sites in T2
is the σ corresponding to solvent-head interaction and occurs over the shortest range out of all, followed by that
solvent-tail interaction. Comparing across the series, it in T1, then T3 and T4. In T2, the cut-off distance
 6

is 21/6 , which is on the tail of the attractive force ap- surfactants were reduced from 300 to 200 in the next two
proximating London dispersion force. Again the negative trials to see if fewer polymer chains would prompt fewer
force contributes to larger micelle aggregates and fewer blocked water channels. T2c and T2d both have a sur-
monomers observed in this system. Even though the factant concentration of 46%, which is 14% less than the
irregularly shaped, holes-riddled membrane structure is previous two; the former has the same number density as
highly impractical in a fuel cell system, the denser pack- T2a and T2b, while ρ∗ in the latter has been increased to
ing and the more aligned, elongated micelle islands are 0.90. T2c In Figures 6c and 6d, the water channels seem
closer to the traits sought after in an ideal proton con- to be less congested, which is the desired result; how-
ducting polyelectrolyte membrane, as they may delineate ever the micelle clusters are oriented in a more disorga-
solvent channels from surfactant chains for more direct nized fashion than in T2a and T2b, and the membranes
proton pathway. appear more fragmented as well. Once again fewer sur-
factants point to destabilization of the membrane, and
this condition combined with high number density may
have disturbed further the balance between surfactant
effect and long range interactions in these particular T2
potential systems. So far T2a is the best variant of T2,
and since improvement of T2a is not achieved by varying
number density or surfactant concentration, that leaves
bond length of surfactant as the remaining parameter to
be varied (initial solvent square lattice and number of
tails need to be kept constant for direct comparison with
(a) T2a (b) T2b its parent configuration T2). The bond length of surfac-
tants in T2e and T2f has been extended from 0.75Å to
1.0Å the number of surfactants is 300 in the former and
250 in the latter. Overall T2e in Figure 6e is not an im-
provement of T2a; in fact there are more discontinuities
of the membrane than its progenitor. A 6.4% decrease of
surfactant concentration lets T2f fare marginally better
out of the two, but its fragmentary character still proves
that it is not an improvement over T2a.

(c) T2c (d) T2d

(a) TT0 (b) TT1 (c) TT2

FIG. 7: TT systems

(e) T2e (f) T2f The results from T and T2-related systems suggest
that a smaller σ would cause ellipsoidal micelles to adopt
FIG. 6: Simulation snapshots of systems based on T2 a more prolate geometry, which may form alternating wa-
ter channels à la lamellar filters given an auspicious pa-
The unusual feature of T2 led to another series of sim- rameterization. Based on the aforementioned Lennard-
ulations based on the T2 potentials. The solvent density Jones potential plots, the value of σ = 0.5 in T2-type
of the systems were altered to see if porosity would be systems may be too small that the membranes become
lessened, and the series results are displayed in Figure 6. over constricted. In testing this hypothesis, another sys-
The number density of solvent, ρ∗ , indicates that there is tems, called TT1 and TT2 are constructed and simu-
a certain number of solvents per width of a solvent grid lated. The Lennard-Jones potential parameters for TT
squared, and that ratio is 0.7 for all systems considered systems (name stemming from the fact they are tests of
thus far. ρ∗ has been increased to 0.85 in T2a and 0.87 the previous trials) are identical to those defined for T2
in T2b. The configurations in Figure 6a and Figure 6b systems, except σ in solvent-head and solvent-tail inter-
show greater continuous surface area as well as longer actions have been changed from 0.5 and 0.5 to 0.8 and
micelle chains than T2, but there are slightly more in- 1 respectively. The difference among TT0, TT1 and
terstitial voids in these two systems. The numbers of TT2 is that the surfactant bond length is 0.75Å 1.0Å
 7

TABLE IV: MT Lennard-Jones potential parameters the presence of an extra tail segment; the potential pa-
rameters for the MT systems are tabulated in Table IV.
Interaction Sites ε σ rc MT1 sees an increase of surfactants by 50 compared to
1-1 1.0 1.0 2.5 MT0; the initial solvent square lattice has also been en-
solvent-head 2-2 1.0 1.0 2.5 larged by 5Å to accommodate simulation error. MT2
1-2 1.0 1.0 2.5
was meant to be the system with further surfactant con-
centration increase than MT1, but due to programming
3-3 1.0 0.75 2.5
constraint, the number density was lowered from 0.7 to
4-4 1.0 0.50 2.5 0.65 instead. The simulation snapshots in Figure 8 show
tail-tail 5-5 1.0 0.25 2.5 a generally more promising class of micelle systems. A
3-4 1.0 0.67 2.5 direct comparison between MT0 in Figure 8a and its
3-5 1.0 0.34 2.5 2-tailed surfactant counterpart S1 in Figure 2b reveals
4-5 1.0 0.67 2.5 that the aggregates have not only grown in size due to
1-3 1.0 1.0 1.12246 extra tail segment in the former, but there also appears
solvent-tail 1-4 1.0 1.0 1.12246
to be some micelles amalgamated to form a large contin-
uous bulk. However, there also seems to be a moderate
1-5 1.0 1.0 1.12246
increase of monomers in MT0 relative to S1. Figure 8b
2-3 1.0 0.88 1.12246 suggests that while MT1 may have fewer numbers of
head-tail 2-4 1.0 0.75 1.12246 monomers, there seems to be less large micelle amalga-
2-5 1.0 0.62 1.12246 mation and more smaller micelles dispersed in the mem-
brane. This is most likely due to simulation box dilation
accompanying the minute increase in surfactant concen-
and 1.1Å respectively. 0.75Å is the typical bond length tration; even though there are more surfactants available
between two hydrogens, and extending the bond length for micellization, the expanded surface area discourages
in TT1 and TT2 is effectively moving away from the global agglomeration, hence there are more localized mi-
simplest diatomic bond towards a more realistic bond celles. MT2 mimics having higher surfactant concentra-
length for a coarse-grained polymeric system. Neverthe- tion by lowering its number density, and the aftermath of
less the new bond lengths are still considerably shorter lowered solvent density may be manifest in the dimpling
than an average carbon-carbon bond of 1.54Å.22 The sim- of the membrane seen in Figure 8c. The micelle clusters
ulation results seen in Figure 7 shows a vast diminution are connected to one another like a network, which is
of porosity in the TT systems. The micelle clusters are an arrangement sought after in forming water channels.
not aligned neatly, but there are canal-like expanses in However upon closer inspection, this surfactant network,
between micelles that could serve as proton conduction resembling branched nerves, bisects in a manner that
channels. obstructs direct proton transportation across the mem-
brane. Nevertheless, the MT systems are closest to the
theoretical configuration of a polymeric proton conduct-
ing membrane out of all trials.

(a) MT0 (b) MT1 (c) MT2 B. SOFCs

FIG. 8: VMD snapshots of MT systems. The

additional tail segment is represented by the color The average time in KMC simulation of proton conduc-
mauve. tion in YBaZrO3 was 0.21ps. Out of 205 possibilities, 2
of the proton conduction pathways found are pictured in
Figure 9. The pathway in Figure 9a shows proton mov-
The last variable to be explored is surfactant chain ing through the perovskite oxide in a loop starting at
length. In MD simulations, this is achieved by insert- site H63 in a 10-step TRRTRTRRTR pathway, where T
ing more tail segments. For reason unbeknownst to the and R stand for transfer and rotation respectively. Fig-
author, simulations of systems with more than 3 tail seg- ure 9b shows proton migrating from site H43, moving in
ments were unsuccessful even after lengthy debugging a 20-step RRTTRTRRTRTRRTRTRTRT pathway that
sessions. Because of this, only systems with 3 tail seg- eventually leads to the periphery of the unit cell towards
ments were simulated; these systems were later referred another image. The brief computation of the same sys-
to as MT, short for “micelle tails”. The basis system tem was re-run at temperature 1200K, but since no sig-
MT0 has the same specifications as S1, plus additional nificant different pathways were found, the results are not
Lennard-Jones parameters that characterize the effect of presented here.
 8

brane in other computational investigations to see if the

archetypal fuel cell membrane maps onto our generic
results. The simulations of generic surfactant-solvent
PEMs were unsuccessful in finding a proton conducting
structure with distinct water channels, but the endeavor
was not unfruitful in that much useful information for
preliminary exploration of solvated PEM structures was
gleaned from the trials. The trials with S0, S1 and
SE systems demonstrate that there needs to be an ade-
quate amount of surfactants in a PEM fuel cell to provide
structural guidance for water channel formation as well
as membrane mechanical stability. The ideal percent-
age of surfactant by weight depends on the surfactant
size and the values of ε and σ Lennard-Jones parameters
that characterize the material, but based on the results
of this study it should be no less than 60% generally
speaking. This is a very loose estimate when compar-
ing with literature studies of solvated Nafion, in which
simulation systems contain water content no more than
20% by weight.23 However direct comparison between a
two-dimensional generic surfactant-solvent system and a
(a) Proton pathway at displacement fraction 0.488585
three-dimensional Nafion-water system is improbable, es-
pecially considering the particles in the generic system
all have mass of 1 molecular weight unit, and their inter-
atomic/intermolecular potentials do not correspond to
any actual physical substance that may be suitable for a
PEMFC. Nevertheless, this study does prove the concept
that there ought to be a high surfactant-to-solvent ratio.
The venture into exploring various ε and σ in the
Lennard-Jones potential yielded some bizarre morpholo-
gies of hydrated membranes seen in Figure 4b and Fig-
ure 6. While there exist too much interstitial holes, the
micelles are arranged in a manner akin to parallel chan-
nels. As discussed in the Results section, these interest-
ing features of these systems most likely arose from the
decreased value of σ in the solvent-head solvent-tail in-
teractions. The reduced σ values shift the potential curve
to the left, thereby making the long range attraction
much more dominant within the cut-off intermolecular
distance, which physically equate to very tightly packed
particles in the membranes. In Nafion, the Lennard-
Jones parameter σ for each atom is at least 3Å whereas
ε corresponding to each atom is on average fraction of
those in the T2 systems.24 Based on these values, the re-
(b) Proton pathway at displacement fraction 0.754597 duced σ in T2 systems is highly unlikely to correspond to
any physical occurring element. These Nafion Lennard-
FIG. 9: Proton conduction pathways of Y-doped Jones parameters would predict a membrane with more
BaZrO3 . The light pink, red, green and blue atoms dispersed micelles if the polymer chain were comprised
represent hydrogen, oxygen, zirconium and yttrium of only 3 segments. Just as a higher ratio of surfactant-
atoms respectively. The barium atoms are made to-solvent can stabilize a membrane surface, so can the
invisible for easier visualization. extended length of surfactant tails, as demonstrated in
the MT systems in Figure 8. Since Nafion is normally
modeled as oligomers composed of at least ten monomer
units in simulation studies,24 interstitial gaps are rarely
IV. DISCUSSION observed.
The polymer electrolyte membrane is the heart of
In parts of this section, the findings from this study membrane electrode assembly; it only conducts protons,
are considered in the context of hydrated Nafion mem- which forces electrons at the anode to travel through an
 9

external circuit towards the cathode, generating electrical to exhibit the most favorable proton/water mobility
power in the process.25 This selective permeability to pro- due to its optimal polymer side chain length and
tons prevents reactant crossover, which would result in high water channel network connectivity. It is the
decreased fuel efficiency. Prevention of reactant crossover superior structural form that allows Nafion to perform
is problematic in thinner PEMs, which are developed as the function of proton conduction most efficiently in
response to demand for mobile high power density fuel the polytetrafluoroethylene-based membrane family.29
cells. In slimming down PEM, the resistance through Conversely, it is the same structure that limits the PEM-
the membrane is decreased and power output is thereby FCs to function at temperature below 100◦ C. On the
amplified; but reduced resistance also leads to compro- other end of operating temperature spectrum, the solid
mise of mechanical durability and increased vulnerability crystallite structures of solid ceramic electrolytes require
to reactant crossover.26 Furthermore, the water content operating temperature in the range of 800-1000◦ C to
dependence of Nafion-type PEMFCs confines these fuel reach a state that enables proton mobility. The low
cells to operate at conditions below 100◦ C which is in- operating temperature condition for PEMFCs makes
convenient for the automobile industries because the heat them ideal for mobile applications, where as the con-
management component would have to be cumbersomely tinuous heating implicit to high temperature activation
big to ensure uniform temperature distribution in the of conductivity for SOFCs relegates them to stationary
fuel cell stack; in fact, the lower the operating temper- applications. From these two types of proton conducting
ature, the larger the radiator needs to be.27 For these electrolyte fuel cells, we see concrete manifestation of
reasons, alternatives to Nafion-type PEMs that operate how forms and functions are interlinked.
at elevated temperature are being developed. PEMFCs
that operate above 130◦ C have been reported to have
advantages of: (a) increased catalyst tolerance of impu-
rities; (b) simplification of water-management; (c) and 1 “Aventador LP 700-4 Technical Specifications,” .
2 A. Macdonald and M. Berry, “Energy Through Hydrogen: Re-
enhanced rejection of heat (thus reduction of radiator
size in automobiles).28 The proton conductivity in novel search Notes,” Heliocentris (2000).
3 K. D. Kreuer, Journal of Membrane Science (2001).
PEMs is independent of water content, and the generic 4 L. Malavasi, C. Fisher, and M. S. Islam, Chemical Society Re-
system of this study should still be applicable to these views (2010).
PEMs, it is just that the topology desired would differ 5 L. Sullivan, Lippincott’s Magazine (1896).
6 X. Kong and K. Schmidt-Rohr, Polymer 52, 1971 (2011).
from that we look for in a PEMFC that utilizes water- 7 K. A. Mauritz and R. B. Moore, Chemical Reviews 104, 4535
dependent proton conduction (vehicular mechanism).
(2004).
The issue of temperature is one factor that was not 8 K. Schmidt-Rohr and Q. Chen, Nature Materials 7, 75 (2007).
addressed in our study of generic PEM systems. The 9 R. Pool and P. G. Bolhuis, Physical Chemistry Chemical Physics

generic system was not made to emulate any particular (2006).

10 M. E. Selvan, E. Calvo-Muñoz, and D. J. Keffer, The Journal of
type of PEMFCs (even though our interpretation of re-
Physical Chemistry B 115, 3052 (2011).
sults have been leaning towards Naifon-type PEMFCS), 11 R. Jorn and G. A. Voth, The Journal of Physical Chemistry C
and the reduced Lennard-Jones units used throughout 116, 10476 (2012).
simulations most likely would not correspond to oper- 12 E. Weinan and B. Engquist, Notices of the AMS 50, 1063 (2003).
13 D. Raabe, M. Scheffler, K. Kremer, W. Thiel, J. Neugebauer,
ational temperature range of typical PEMs, integrating
temperature as a variable seems incompatible with the and M. Jansen, “Multi-Scale Modeling in Materials Science and
Engineering,” (2010).
goal of exploring polymer membrane without any spe- 14 J. A. Elliott and S. J. Paddison, Physical Chemistry Chemical
cific chemical characteristics. Physics (2007).
15 J. Karo, A. Aabloo, J. O. Thomas, and D. Brandell, The Journal

of Physical Chemistry B 114, 6056 (2010).

16 E. Spohr, P. Commer, and A. A. Kornyshev, The Journal of
V. CONCLUSION Physical Chemistry (2002).
17 F. Ercolessi, “A Molecular Dynamics Primer,” .
18 J. Li, in Handbook of Materials Modeling, edited by S. Yip
The morphology of solvated PEMs is indeed an elusive (Springer Science+Business Media, 2005) Chap. 2.8.
subject; the difficulties encountered in this rudimentary 19 A. Hinchliffe, Molecular Modelling for Beginners, 2nd ed. (Wiley,
survey alone are already quite daunting, therefore it 2008).
20 S. Plimpton, A. Thompson, and P. Crozier, LAMMPS Docu-
is perfectly logical that there exists a multitude of
conflicting theories on hydrated Nafion membrane in mentation (2013).
21 A. F. Voter, Radiation Effects in Solids , 1 (2007).
computational literature. Without obtaining the mem- 22 T. R. Gilbert, R. V. Kirss, N. Foster, and G. Davies, Chemistry:
brane structure that facilitates proton conduction, there The Science in Context, 2nd ed. (W. W. Norton and Company,
is no way to study the proton conduction mechanism. 2009).
23 S. Cui, J. Liu, M. E. Selvan, D. J. Keffer, B. J. Edwards, and
For a PEMFC that relies on water transportation of
protons, the hydrophobic phases in the membrane W. V. Steele, The Journal of Physical Chemistry B 111, 2208
(2007).
should arrange themselves in a fashion that would de- 24 A. Vishnyakov and A. V. Neimark, The Journal of Physical
marcate water channels, hence maximizing water/proton Chemistry B 104, 4471 (2000).
mobility for greater fuel cell efficiency. Nafion is known 25 N. Djilali, Energy 32, 269 (2007).
 10

26 J. W. Weidner, V. A. Sethuraman, and J. W. Van Zee, The 28 Z. Yang, D. Coutinho, F. Feng, J. P. Ferraris, and K. Balkus,
Electrochemical Society Interface 12, 40 (2003). Jr, Prepr. Pap.-Am. Chem. Soc., Div. Fuel Chem 49, 599 (2004).
27 Arizona State University, “Innovative fuel cell project aimed at 29 D. Brandell, J. Karo, A. Liivat, and J. O. Thomas, Journal of

meeting large power needs,” (2007). molecular modeling 13, 1039 (2007).