1. Introduction
Polymer electrolyte membrane (PEM) fuel cell, being more efficient, environmentally friendly, and modular, is one of the most outstanding candidates to replace the internal combustion engine [
1,
2,
3,
4]. Fluorinated or hydrocarbon-based proton exchange polymeric ionomers serve as the central component in PEM fuel cells owing to their crucial roles. In order to use them in viable technological applications, these materials need to have superior properties like high proton conductivity, thermal and chemical stability, impermeability to gases, and easier compatibility with electrodes [
5,
6,
7]. Morphological peculiarities of these materials mostly stem from the micro-phase separation of hydrophobic and hydrophilic parts of the macromolecules [
8]. Many studies have displayed that the conduction of protons takes place via hydrophilic domains that are resulted from micro-phase separation [
9]. Morphological attributes of these membrane materials are very much affected by water and the backbone’s chemistry.
High-temperature (100–200 °C) PEM fuel cells have been of great interest to researchers because of their outstanding benefits in terms of higher catalytic activity, CO poisoning elimination, low cost, device cooling, and water management. Several authors have studied novel polyelectrolyte systems with high conductivity at temperatures from 100 °C to 200 °C as an alternative solution to the commercial perfluorinated polyelectrolytes (e.g., Nafion), which show high conductivities only through water-assisted proton conduction [
10,
11,
12] and operate at temperatures below 100 °C [
13]. To date, building strong acid-base complexes among functionalities connected to the polymeric backbones has been one of the most well-known methods that enable high proton conductivity under low humidity or even anhydrous conditions [
14,
15,
16,
17,
18,
19,
20]. In these materials, a homogeneous polymeric blend is formed by strong multiple acid-base interactions and hydrogen bonding networks, allowing proton conductivity by Bronsted acid-base pairs.
Understanding the process of micro-phase separation in various ionomers is important and leads researchers to employ different modeling methods. It is surmised that molecular features and processing conditions result in mesoscale alterations at different length scales in morphology. Atomic simulations by molecular dynamics (MD) and Monte Carlo (MC) simulation not only help to understand the molecular structures of small systems, but they are also helpful in analyzing the aggression of the sulfonic acid groups and local transport of protons in perfluorinated sulfonic acid (PFSA) ionomers. On the other hand, these simulations are not useful to analyze the nanoscale morphology in complicated inhomogeneous phases because they require more atoms and excessive time for equilibration. For this reason, many modeling methods are used on a comparable scale with experimental SAXS (small-angle X-ray scattering) data in these systems [
21,
22,
23]. The coarse-grained and continuum are mainly two methods for convenient modeling of intermolecular interactions. Individual conceptual sites with appropriate non-bonded interaction parameters can be put in place of chemical groups of a complex polymer in particle-based models. However, local concentration fields with collective variables are used to figure out self-organizing structures for field-based methods. The length is scaled in 10–103 nm, and time is scaled up to milliseconds (even seconds), enabling the simulation of larger size systems and getting results comparable to experimental data. The morphology evolution and eventual phase separation of several ionomer systems have been investigated via mesoscale modeling using dissipative particle dynamics (DPD) simulations [
24,
25,
26,
27]. Moreover, self-consistent mean-field (SCMF) simulations are also employed in order to study how temperature and water content affect phase separation and morphology in PFSA membranes [
28].
There are various studies on PFSA-based proton conducting blend membranes as an alternative to costly fluorinated membranes. Although mesoscale morphologies of hydrated PFSA-based membranes were investigated both experimentally and computationally, very few studies have investigated the mesoscale phase behavior of blend membranes for polymer electrolyte fuel cells [
26]. The present study is complementary to our previous experimental studies on proton-conducting PFSA-based polymer blends. In these studies, Nafion/poly(vinyl phosphonic acid) (Nafion-PVPA) [
16] and Nafion/poly(1-vinyl-1,2,4-triazole) (Nafion-PVTri) [
29] blend membranes with significant anhydrous proton conductivity were fabricated. However, it is not possible to experimentally analyze the mesoscale morphologies of these blends. In this work, mesoscale phase morphologies of Nafion-PVTri and Nafion-PVPA blend membranes with different compositions were investigated using the DPD simulation method. Simulation results show that both blend membranes can form a highly separated microstructure due to the hydrophobic and hydrophilic character of different polymer chains forming the blend of different segments in the same polymer chain. This study uses dissipative particle dynamics to construct the mesoscale phase morphologies of Nafion-PVPA and Nafion-PVTri blend membranes. Herein, we establish a connection between the simulated morphologies and experimental properties.
2. Theoretical Background
In order to understand complex fluid dynamics, mesoscopic simulation methods are significant and have been widely used. Hoogerbrugge and Koelman were the first to create the DPD method, which is later improved by Groot and Warren [
30,
31,
32]. In DPD, Fluid material’s mesoscopic regions with similar chemical features are acted as fundamental particles called “beads. It is unlike MD that these particles are assumed to have no atomic properties. This is because all degrees of freedom that are greater than a bead radius in the DPD unit are neglected. Thus, coarse-grained interactions among beads are computed, and all atomic details are gone.
Simulations are done on a combination of particles that interacts by Newton’s equations of motion.
The masses are adjusted to 1; thus, the acceleration equals the force acting on a particle. The sum of three pairwise additive components gives the force
:
The summation runs over all other particles j in a given cutoff radius for particle i. It is this short-range cutoff that makes the interactions local. The cutoff radius is assumed to be unity for simplicity as this is the only length–scale in the system. So, forces beyond the bead radius are not taken into account for all beads. Soft repulsion interactions are taken into consideration, whereas the excluded-volume effect is not functional.
is the conservative force, a soft repulsive interaction between the particles’ centers i and j. is the drag force (or dissipative force), and it depends on the relative velocities of the beads i and j. It normalizes velocities to decrease the relative momentum of the particles. denotes a random force stabilizing the temperature.
The repulsive conservative force
exists between the centers of the
ith and
jth particles, and it is given as:
stands for the maximum repulsion force between the
ith and
jth. particles; and
is given as
,
,
.
, the dissipative force acting to reduce the relative momentum of two beads depends on the relative velocity of these beads, and it is given as:
is a friction coefficient,
is a short-range weight function, and
. Due to the dissipative force structure chosen, the total momenta for any pair of particles and hence the system overall is conserved.
For the random force given below, a different distance-dependent
is employed. This force is between all pairs of beads with a similar short-range cutoff:
is a random variable with Gaussian statistics satisfying
, and
.
It has been proved that one can define only one of the two weight functions
or
, and the other one is defined accordingly [
33].
In DPD simulations, beads move freely. Some molecules like polymers and surfactants are represented by two or more beads. Therefore, to assure the connectivity of the chain, an additional force between successive beads must be used:
In the DPD method, the particles refer to beads rather than molecules or atoms, and beads contain coarse-grained groups of atoms. Bead by bead interactions take the places of atomic-level interactions. The engagement of springs between two successive beads along the chain further enables the use of this system for polymers. The random and dissipative forces perform in unison. Their total effect is a thermostat. This fine conservative repulsive force reveals the system’s leading chemistry. The parameters represented by are named bead–bead repulsion parameters. It is because these parameters are dependent on underlying atomistic interactions that they are also referred to as DPD interaction parameters.
The simulations of liquid-liquid and liquid-solid interfaces can be done by employing the DPD method. Therefore, it is very similar to the mean-field theory of Flory–Huggins for polymer chains and can be considered a continuous version of the lattice model established therein [
34,
35,
36]. The mean-field theory is used to describe polymer’s miscibility in a given solvent, and this is done via comparing the free energy of the system before and after mixing. Description of thermodynamics of polymer blends is done with the help of a similar theory, and block copolymers and their blends with homopolymers can also be described with some modifications. The mixing energy is directly related to the dimensionless Flory–Huggins interaction parameter
χ, which represents the change in energy in units of
kT when a bead of
A is drawn from an environment of bulk
A and surrounded with a bead
B in an environment of bulk
B.
A and
B would be polymer or molecules.
The only computational parameter that the beads’ chemistry is involved in the conservative force as noise couples with dissipation. The relationship between the repulsion parameter
acting between neighboring beads and the Flory–Huggins interaction parameter
is described by Groot and Warren [
31]. The repulsion parameter and the energy of mixing are proportional to each other.
The Flory–Huggins model in binary systems is the most well-known theory employed for the thermodynamics of phase separation and mixing. The general equation for the system’s free energy of mixing (
) with bead types A and B can be expressed as:
and
stand for the volume fraction for type A and B beads, respectively. In Equation (7), the first two quantities refer to the entropy of mixing, and these quantities vary with respect to the number of distinct ways in which the chains can come together on a lattice, and the last term is the free energy of interaction. The positive interaction parameter
leads to various phase diagrams.
The Flory–Huggins interaction parameters can be defined as:
is the energy of mixing describing the change in free energy because of the interaction between the pure and mixed state. It can be calculated using the traditional Flory–Huggins model:
denotes coordination number and
is the binding energy for the chosen
ith and
jth components.
, the Flory–Huggins parameter can also be calculated from Hildebrand solubility parameters, which is calculated from cohesive energy density values:
where
is the beads’ average molar volume,
and
denote the parameters for solubility for beads
A and
B, respectively, and
is the molar enthalpy for vaporization. Parameters for solubility vary depending on the species’ chemical nature in question.
3. Computational Methodology
Chemical structures for Nafion, PVPA, and PVTri are shown in
Figure 1. To construct systems for the DPD, copolymer Nafion is divided into three beads (A, B, and S), and homopolymers of PVPA and PVTri are represented by a single bead (P for PVPA and T for PVTri). Since the fluorine-containing backbone of Nafion is chemically different from the ether and the sulfonic acid-containing side chain, the backbone is constructed as a single bead (A) while the side chain is built as two different beads (B and S). By attaching extra fluorine and hydrogen atoms to the connection points of Nafion and homopolymer beads, respectively, consistency is obtained.
The DPD modules of the Material Studio suite of programs were used in all DPD simulations carried out. The optimization of the particles’ structures was done via molecular mechanics along with COMPASS force field parameters by employing the Forcite module in Materials Studio. The number of particles per unit, namely, the simulated density value, was set at ρ = 3. Calculation of interaction radius (
rc = 4.5 A) was done using the beads’ average size. The Blends and Synthia modules in Materials Studio are used in order to calculate the Flory–Huggins parameters corresponding to bead–bead interactions. Interaction parameters of DPD are computed by using the parameter
. Groot established a relationship between the parameter
and the repulsion parameter
aij between the beads
i and
j [
31]. The density affects this relationship, and the density is set at 3 in our computations:
In a simulation box whose edges are 36 nm, about 192,000 DPD beads were positioned in arbitrary locations at the beginning. It is equilibrated for a period of 6 ns. The examination of the blend morphologies is done via homopolymer density contour plots, which are formed from the equilibrated mesostructures.
4. Results and Discussion
4.1. Flory–Huggins Parameters and Interaction Energies
Conventional PFSA-based membranes’ proton conductivity highly depends on the water content in their protonated form. This dependency limits the fuel cells’ working temperature to below 100 °C due to the membranes’ dehydration from hydrophilic cavities at elevated temperatures. These membranes typically show phase separation where the network of hydrophilic nanopores and nanochannels are embedded in a hydrophobic phase domain. The hydrophilic nanoporous phase contains water and acidic sites, facilitating proton conductivity by transporting free protons. The hydrophobic phase serves as a mechanical force to stabilize the morphology of the membrane. Due to the drawbacks of low-temperature PEM fuel cells, developing high-temperature proton conducting membranes is needed. Recently, having thermal and chemical stability, acid-base polymer blends have been examined as alternative membranes for PEM fuel cells operating above 100 °C. Among them, Nafion-PVPA and Nafion-PVTri show promising properties and considerable proton conductivities.
Flory–Huggins interaction parameters
for chosen components
i and
j are calculated through Equation (8), using the binding energies and coordination number. The
parameters (
Table 1) give pairwise information about how the beads will interact with each other. The backbone segment of Nafion is immiscible with all other beads except for the bead B. The molar energy of vaporization for the nonpolar materials should be less because of their weak intermolecular energies. The polar parts result in cohesive energy density, which is what is calculated in this study for the fluorine-containing groups. A highly positive
indicates the contacts of the bead A with the beads S, P, and T. The contact of the bead B with the beads S and P are less favorable when compared with the contact of bead A with bead B and bead B with bead T. A negative value of
indicates the interactions of the bead S with P and T is highly preferred due to complexation tendency between corresponding segments.
The various trends investigating the solubility parameters of different beads can be elucidated by partitioning short-range interaction energy into hydrogen-bonding, London dispersion (or van der Waals type), and dipole-dipole types in the Hansen solubility parameters. There is only the former one for the case of nonpolar molecules. Therefore, It is due to due to weaker intermolecular interactions that cohesive energy density or molar energy of vaporization for the nonpolar beads is less than those for polar beads. Hence, parameters of the solubility for the nonpolar fluorinated segments and polar triazole and oxygen moieties are placed on distinct solubility spectrum sites. For the molecules’ consistency, attaching hydrogen atoms to the head and tail points might cause little variation in the computation of solubility parameters, and this might be discarded in this study.
In theory, liquids having similar cohesive energy densities are expected to carry similar characteristics for solubility. The parameter indicates that segments containing fluorine are immiscible with other segments. The smaller solubility parameter of fluorine-containing bead has a smaller solubility than other beads, which shows that by increasing the structure’s fluorine content, the immiscibility of fluorine-containing segments can be increased. The other beads, on the other hand, have larger solubility parameters due to their chemical structure.
4.2. Morphology
There have been many computational studies to understand the phase morphologies of hydrated polymer electrolytes. The control of phase behavior is significant to improve membranes’ proton conductivity. In this work, the effects of the molar composition of Nafion/PVPA and Nafion/PVTri on the mesoscopic morphologies of blend membranes were investigated. As the concentration of PVPA increases, PVPA clusters become continuous after a 1:2 molar ratio (
Figure 2). Proton conductivity increases with PVPA content. On the other hand, Nafion-PVPA blend membranes are not stable with high PVPA content. In our experimental work, maximum proton conductivity was obtained with 1:3 composition [
16]. In that composition, it was clear that the PVPA phase was continuous enough to promote proton conductivity through interaction with sulfonic acid groups of Nafion.
Herein, Nafion, having an amphiphilic character because of incompatible polar and nonpolar segments, was considered. These segments tend to reside away from each other due to having distinct interaction parameters, which results in phase separation of chains. However, they are joined together due to a covalent bond that increases their involuntary compatibility, enabling only microscopic phase separation. While this transition from disorder to order causes high order continuous phases in the case of Nafion-PVPA blends, there is a comparatively low-order phase of cylinders in Nafion-PVTri blends (
Figure 3), which is, in fact, consistent with the previous observations done for Nafion-containing blends.
4.3. End-to-End Distance of Blend Membranes
A single polymer chain can appear in many possible conformations ranging from a tight coil to a straight chain. The probability of a single polymer chain having a specific end-to-end distance rises as the number of possible conformations that can have this end-to-end distance rises. A straight chain can be produced by only one conformation; however, more possible conformations exist as the molecule coils up more. Therefore, a polymer chain is virtually inclined to coil up. A model that assumes that molecules are composed of many segments can predict the end-to-end distance of a polymer chain. All the segments in this model are rigid, but they are freely jointed at both ends and hence capable of making any angle with the next segment. Therefore, a molecule model can be constructed by connecting each consecutive segment at an arbitrary angle. The name of this process is a random walk. In the present study, the end-to-end distance for Nafion chains, an important structural property of polymeric materials, is computed according to PVPA and PVTri by employing DPD simulations. Since there is a high interaction parameter between sulfonic and phosphonic acid groups, root means square (RMS) of end-to-end distances of polymer chains in Nafion-PVPA are smaller than those in Nafion-PVTri (
Figure 4).
RMS end-to-end distance is highly dependent on the type and composition of the blends. The smallest RMS end-to-end distance is obtained in 1:1 composition in both systems, which may be attributed to the fact that polymer chains in the exact one-to-one ratio prefer tight coil conformation. The length of individual chains in the system does not change in any environment, while their relative positioning in the bulk of the system changes to lead to diverse morphologies. As PVPA content increases RMS end-to-end distance increases, this is because the continuous PVPA phase starts. However, RMS end-to-end distance does not vary significantly through changing the composition of the Nafion-PVTri blend system. It is because DPD is a coarse-grained method, and the only parameter that carries the atomistic information is the soft repulsion parameter that these results should be confirmed using other conformational search methods (e.g., high-temperature Molecular Dynamics).
4.4. Diffusion Rates of Blend Membranes
Diffusivity rates of beads can be seen in
Figure 5. Diffusivity of the bead C is the highest among other bead diffusivities. Since C is located on the polymer’s side chains, it is relatively more mobile than the other beads like A on the chain, making its diffusivity the highest.
4.5. Density Profiles of Blend Membranes
In the x-direction, the density profiles for the beads A and P in the Nafion-PVPA system and the density profiles for the beads A and T in the Nafion-PVTri system are presented in
Figure 6a,b, respectively. The results are in line with the previously conducted studies demonstrating morphology. Both P and T are highly immiscible with A that densities of them will decrease as the density of A increases and vice versa. However, the degree of repulsive interaction parameter directly affects the density profiles of the blend membranes. The density variation in Nafion-PVPA is higher than in Nafion-PVTri due to higher repulsion between the beads A and P.