Model-Based Analysis of PFSA Membrane Mechanical Response to Relative Humidity and Load Cycling in PEM Fuel Cells

The hydration/dehydration cycling of the membrane during the fuel cell operation and the resulting mechanical stress are in part responsible for the mechanical failures of a polymer electrolyte membrane. To thoroughly investigate the mechanical behaviors of the membrane under in-situ cyclic conditions, in this paper, we have interfaced a comprehensive two-dimensional transient fuel cell transport model with a viscoelastic-plastic membrane mechanical model. The transport model is used to produce the spatiotemporal proﬁles of membrane water content and temperature in an operating cell, which are then sent to the membrane mechanical model to calculate the mechanical parameters of interest. This providesthe extendedcapability ofstudying the membrane mechanical response under in-situ conditions, which was not possible with the original mechanical model. The effects of cycling relative humidity and voltage and current at different temperatures on the membrane stresses are studied using the coupled model. It is found that the location of maximum in-plane tensile stress can vary signiﬁcantly with the operating temperature during voltage cycling, whereas the highest in-plane compressive stress occurs under the land for all cases, particularly near the cathode. The simulation results conﬁrm the need for coupling the two models to capture comprehensive transport phenomena in studying the membrane mechanical behaviors, and represent an important step toward improved understanding of various synergistic mechanical failure mechanisms that affect the membrane in an operating fuel cell.

Mechanical failure of the polymer electrolyte membrane in the form of cracks, pin-holes and delamination has been identified as a limiting factor in the durability of the fuel cell.The hydration/dehydration cycling of the membrane during the fuel cell operation and the resulting mechanical stress are in part responsible for the mechanical failures. 1,2For example, the cyclic mechanical stresses have been shown to play a significant role in propagating a crack through the thickness of a membrane. 3Recent in-situ experiments have shown that membrane damage can be dominated by these cracks. 4Since operation of a fuel cell causes hydration/dehydration cycling along with spatiotemporal variations of temperature that can influence the membrane response as well, 5 it is critical to thoroughly understand the mechanical behavior of the membrane under in-situ conditions, so that alleviating strategies focusing on the fuel cell operating conditions can be developed.
][8][9][10] In this prior work, the hydration states (water volume fraction) and temperatures at the interface of gas diffusion media and the gas channels are the model inputs.The spatiotemporal water transport in the electrodes and membrane is governed by diffusion only, without the in-situ features including electro-osmotic drag (EOD), water generation by reaction, and reaction-generated spatiotemporal temperature variation.Therefore, even though the existing model 10 has provided the capability to study the membrane behavior in a cell subject to relative humidity (RH) cycling, its response to electrochemical load cycles (dynamic variations of the cell voltage or current density) has remained unexplored thus far.
The operating conditions for a stack, including current, voltage, system pressure, stack temperature, inlet relative humidity, and flow rate affect the spatiotemporal water distribution in the membrane.Modern automotive fuel cell systems (e.g., Toyota Mirai 11 ) operate without external humidification to reduce the system cost and volume.This self-humidification strategy typically results in even greater spatiotemporal variations of local current, species concentration, temperature, and hydration as observed by comprehensive modeling 12 and experiments. 13Furthermore, it is the operating parameters listed above that can most effectively be controlled in an automotive system and serve as the inputs to a fuel cell transport model.Thus, it is worthwhile to examine the functional relationship, i.e., the direct transfer function from the fuel cell operating conditions to the membrane mechanical response.To that end, a comprehensive fuel cell transport model is needed to convert the typical fuel cell operating parameters to the detailed membrane hydration profile.
In this work, we have interfaced a two-dimensional transient fuel cell transport/performance model in COMSOL Multiphysics 14 with the aforementioned membrane mechanical model in ABAQUS.The transport model is used to produce the spatiotemporal profiles of membrane water content along with the temporal evolution of the average membrane electrolyte assembly (MEA) temperature, which are then sent to the membrane mechanical model to calculate the mechanical parameters of interest.The transport model has comprehensive consideration of species transport and reaction kinetics in a fuel cell, which leads to increased fidelity in predicting the water distribution, and therefore the stress and strain, in the membrane.Using the combined model, we study the impacts of dynamically cycling reactant RH, current and voltage on the membrane to show the feasibility of this approach, and examine the effects of RH/voltage/current cycling on the membrane mechanical behavior with preliminary simulation results.
Chemical degradation of the membrane also plays an important role in fuel cell lifetime, along with various interactive effects between mechanical and chemical stressors that have been proposed or observed. 15Although the chemical degradation effect is beyond the scope of this paper, the model developed here to predict the mechanical behaviors of a membrane under in-situ conditions is a significant step forward.In the future, the modeling capability may be extended to capture the chemical stressors within the current mechanical modeling framework.
The rest of the article is organized as follows.First, both the transport model and the membrane mechanical model are explained in Mathematical Models section with details regarding the numerical implementations for each model.Relative Humidity and Load Cycles section describes the cycles considered in this work.Results and discussions are then presented followed by conclusions and a brief summary of the results.

Mathematical Models
Transport model.-Onemajor contribution of this work is the use of a comprehensive transport model to obtain water and temperature distributions in the MEA during various relative humidity and load cyclings.In particular, the model presented in Ref. 14 is employed in this work.The 2D modeling domain, illustrated in Fig. 1, includes the membrane and catalyst, microporous and gas diffusion layers on both the anode and cathode sides.Since this work focuses on the membrane mechanical response, the transport model equations are not provided for space consideration.However, the model draws from earlier works by Balliet et al. 16 and Zenyuk et al. 17 and is briefly discussed below.
In terms of the transport of different species, Fick's law is used for gas phase diffusion.Pressure drop in both the gas and liquid phases is modeled with Darcy's law, while modified Ohm's law is used for transport of charged species.Water transport across the MEA includes the effects of electro-osmotic drag, diffusion, and thermoosmosis 18 as well as the interfacial resistances. 19It is also important to note that, as described in the literature, 20,21 hydraulic permeation through the membrane is implicitly accounted for by an effective bi-modal water diffusion coefficient. 22,23Finally, the mixed wettability porous medium model 24 is employed to obtain the water retention curves as well as the average Knudsen radii, relative permeabilities, and interfacial liquid-vapor area in the porous media.The simulations in this paper use the porous media parameters reported in Ref. 16.
As for the reaction kinetics, hydrogen oxidation reaction (HOR) is described using the dual pathway kinetics model, 25 while the double trap model 26 is used for the oxygen reduction reaction (ORR).A homogeneous model is used for the cathode catalyst layer (CL), where the local transport resistance for oxygen is described as a summation of diffusional and interfacial resistances. 27inally, recent literature reports on incorporating the effects of ionomer relaxation in the transport models to obtain a more realistic view of the cell transient response. 14,28However, as noted below, the effects of the stresses on the transport phenomena are neglected in the current study and ionomer relaxation is not included in the transport model here.Therefore, a comprehensive investigation of how the cell's mechanical response impacts the transport and performance remains an important future research endeavor.

Membrane mechanical model.
[8][9][10] The material properties and behaviors used for each component of this FE model are as follows: r The membrane swells in response to changes in water content and temperature.The total strain tensor is assumed to be the summation of swelling, thermal and mechanical strain where ε H ij , ε T ij , and ε M ij are swelling, thermal, and mechanical strain components, respectively.Swelling strain is given by the empirical Here θ and θ o are the current and reference temperature and φ P is the polymer volume fraction.This swelling behavior was determined experimentally and the strain is implemented via a user-defined subroutine.
r The thermal strain is given by The stress is derived from the mechanical strain, via a constitutive viscoelastic-plastic model involving two parallel load paths: an elasticviscous load path and an elastic-plastic one 30,31 (described in detail in previous work 10 ), such that Here, ε EP ij and ε EV ij are the elastic-plastic and elastic-viscous strain components of the constitutive model, respectively.All the membrane material parameters used in the constitutive model are explicit functions of temperature, water content and loading rate.These functional relationships were determined by fitting experimental data from a battery of uniaxial tensile and stress relaxation tests conducted on PFSA Nafion 211 membrane to this viscoelastic-plastic constitutive model.Then, this model was incorporated in the FE analysis using the commercial software ABAQUS.The constitutive relations have been validated with experimental data in previous publications. 10,32The membrane is the only component subjected to swelling due to water absorption and desorption.The swelling is assumed to be isotropic.
r The material properties of the electrodes are also viscoelastic- plastic with temperature and water content dependent properties.These properties were determined in experiments similar to those used to characterize the membranes. 33The diffusion media (i.e.GDL and MPL) and bipolar plates are assumed to be linear elastic, where the diffusion medium is transversely isotropic with in-plane and out-of-plane elastic moduli of 1500 MPa and 9 MPa, respectively 34 and the bipolar plate is isotropic with elastic modulus and Poisson's ratio of 10 GPa and 0.25, respectively.Thermal expansion for both of these components has been neglected.
The model is loaded by the imposition of a continuous change in water volume fraction (φ) distribution and average temperature throughout the MEA.In response to this change, components of the model undergo different degrees of hygrothermal expansion or contraction resulting in the development of stresses.
Model coupling.-Inthe earlier implementation of the mechanical model the transport and stresses were calculated simultaneously, providing a mechanism for continuous coupling between the phenomena.As input, RH was imposed on the surfaces of the GDL exposed to the gas channels, and a temperature field was imposed throughout the domain (Fig. 2a).Water transport was predicted by simple Fickian diffusion, and mechanical stresses were calculated from the strain using the constitutive relations.However, as described in the previous section, water and heat transport during fuel cell operation can be far more involved than what is captured by a simple diffusion model.Therefore, to obtain a more accurate picture of the water content and temperature distributions, the aforementioned transport model is used (Fig. 2b), the results of which are imposed directly on the nodes of the mechanical finite element model at each time step.To accommodate this new coupling from separate models, water transport across the mechanical mesh is suppressed.The temperature variation across the membrane is small and has been shown in preliminary studies to only have secondary effects on the stress distributions.Therefore, only the temporal average of the MEA temperature along with spatiotemporal water distribution are imposed in the current mechanical model.
Coupling between the separate transport and mechanical models is achieved as follows.The transport model is used to generate the spatiotemporal distribution of the membrane water content and the temporal average MEA temperature.The membrane water content values are then translated into water volume fraction in the membrane.The spatiotemporal distribution of the water volume fraction along with the node locations for the transport model's mesh grid is used to generate linear interpolation functions in both space and time.These interpolation functions are then queried by the mechanical model at any given time step (determined adaptively by the solver in the mechanical model) to evaluate the water volume fractions at the corner nodes of the mechanical model's mesh grid.Meanwhile, the average MEA temperature is imposed on all the nodes in the mechanical model.This process is carried out through in-house scripts written in MATLAB and Python.
To verify the accuracy of the coupling algorithm, a series of preliminary studies were conducted to select an optimum time step and suitable locations of nodal boundary conditions.In these studies, spatiotemporal water distributions produced by the earlier mechanical model were imposed directly onto the current mechanical model.Once the proper numerical parameters were selected, the resultant stresses were found to be identical.
At this point, it is worth pointing out that both the previously published mechanical model and the current model assume one-way coupling between the transport and mechanics, i.e., the mechanical response is affected by the transport phenomena through the water distribution in the membrane, while the effect of the mechanical response on the transport is assumed to be negligible.Therefore, the key difference between this work and previously published ones [6][7][8][9][10] lies in the linking with the comprehensive transport model to enable direct simulation of electrochemical load cycles.
Since the material behavior in the MEA is time-dependent, this process must be started from a point in time with known stress and water/temperature distributions.The way this is addressed in the current study is to start both the transport and the stress simulations from an equilibrium configuration, such as a zero-stress state under uniform water and temperature distributions.Furthermore, to reduce transient effects seen only at the beginning of loading and ensure a repeatable cyclic behavior, several cycles are simulated using the transport model.However, due to computational time demands, only two cycles are simulated in the mechanical model, since the transport results illustrate a repeatable pattern starting from the second cycle.
Numerical implementation.-Thetransport model is implemented in COMSOL Multiphysics 5.2a.A mapped mesh consisting of 3060 quadrilateral elements is used throughout the domain with increased mesh density in the MEA region.Furthermore, the mesh density is exponentially increased close to the boundaries between adjacent layers to accommodate the different material properties.The generalized-alpha method 35 is used for time stepping with an amplification coefficient of 0.3 and the maximum time step size is limited to 200 milliseconds.The resulting linear system is solved using the MUMPS direct solver provided in COMSOL.To improve the computational efficiency of the model, an under-relaxation scheme is employed, where the value of liquid saturation at the previous time step is used to calculate effective properties such as the diffusion coefficients at the current time step.This was achieved using the Previous Solution operator in COMSOL 5.2a.In a preliminary study, it was found that the under-relaxation scheme can result in up to five times faster solutions in the two-phase regime.Overall, the results for the cycles presented in Results and Discussion section were obtained in 9 to 13 hours on a personal computer utilizing two cores at 3.60 GHz.
The mechanical model is constructed in the software package ABAQUS using eight-node, biquadratic, coupled, temperaturedisplacement, generalized plane strain elements (CPEG8T).To simu-late the clamping of fuel cell stack, a fixed clamping pressure (1MPa) is applied at the top of the bipolar plate.The left end (middle of flow channel), right end (middle of the land) and the bottom edge of the bipolar plate are subjected to zero displacement boundary condition (to approximate symmetry boundary conditions 34 ).A convergence study resulted in a mesh of 30805 elements and 95512 nodes with finer meshing in the MEA and at the corners between the land and channel.The minimum time step is 50 milliseconds, and overall it took 3.5 to 4 hours to complete each simulation on a personal computer using 2 processors at 3.16 GHz.

Relative Humidity and Load Cycles
The US Department of Energy (DOE) provides guidelines for accelerated stress tests designed to induce mechanical degradation of membranes without electrochemical load. 36We simulate these tests by cycling between a completely dry (RH = 0%) and a supersaturated (RH = 148% -corresponding to a dew point of 90 • C) condition at an operating temperature of 80 • C. The RH at both the wet and dry states is maintained for 19 seconds with a transition (hydration/dehydration) time of 1 second between the two.
In addition, electrochemical load cycling simulations are carried out to study the stress evolution in the membrane in an operating fuel cell.Voltage cycles between an upper potential limit (UPL) of 0.7 V and a lower potential limit (LPL) of 0.4 V are simulated at three different operating temperatures (60, 70, and 80 • C).The hold time at each potential level is 29 seconds with a transition time of 1 second.An RH value of 50% is used for the voltage cycles.This value is chosen to highlight the dynamics of membrane hydration and the corresponding stress evolution induced by load cycles with relatively dry gas feeds, which is of interest in automotive applications.Finally, a current density cycle is simulated to determine if one mode of operation is preferable in terms of mechanical stress development in the membrane.For comparison between voltage and current density cycling, the limits for the current density cycle were chosen according to the steady state power output of the voltage cycling case.In this set of simulations, the hold time at each load was 28 seconds with a transition time of 2 seconds and the RH was held at 60%.The transition time and the RH were increased relative to prior voltage cycles to minimize voltage reversals during current cycling due to severe anode dehydration.The details of the simulated electrochemical load cycles are summarized in Table I.
All the cycles presented in this paper (RH and load cycles) use a pressure of 2 bars in both anode and cathode channels.A mixture of hydrogen and nitrogen is assumed to exist in the anode channel (with a hydrogen molar fraction of 0.75), while air is used on the cathode side.The anode gas composition was chosen to give a realistic representation of the conditions seen in automotive stacks with hydrogen recirculation.

Results and Discussion
Transport modeling results.-Weprovide a brief overview of the transport modeling results in this subsection to highlight features of the average dynamic response of the cell.However, discussions about the local transport phenomena are provided in the subsections on the membrane mechanical response.
RH Cycling.-TheRH cycling results are presented in Fig. 3.In particular, the RH cycle is shown along with the average water volume fraction (φ) in the ionomer phase and the average MEA temperature.A few observations are in order.First, we see that even after the channels reach supersaturated conditions at 20 seconds (Fig. 3a), the ionomer phase continues to absorb water before reaching its maximum water content at about 26 seconds, after which it reaches steady state (Fig. 3b).During this time, however, water continues to condense in the porous media due to the supersaturated conditions.The water absorption into the ionomer phase along with the condensation in the porous media releases heat, which results in an increase in the average MEA temperature (Fig. 3c).
Second, as the cycle goes through the dehydration phase, a somewhat counterintuitive increase in the ionomer water content is observed initially at about 35 seconds (Fig. 3b).To understand this behavior, note that as dehydration begins, the liquid water in the catalyst layer starts to evaporate, leaving the gas adjacent to the ionomer phase in saturated conditions.However, the evaporation results in a decrease in the MEA temperature (Fig. 3c, at around 35 seconds), which in turn increases the ionomer water storage capacity.The combination of these effects results in the slight increase in the ionomer water content at the start of the dehydration process.
Finally, we note that the relatively fast diffusion through the MEA thickness results in negligible water content gradient in the throughplane direction.Nevertheless, a slight difference is observed between the water volume fraction values in the anode CL and the rest of the MEA (the difference between the lines in Fig. 3b).This can be attributed to the fact that water vapor diffusion is faster on the anode side due to the anode gas composition.
Voltage and current cycling.-Tostudy the membrane mechanical response to hydration and dehydration processes induced by electrochemical load cycling in an operating cell, we ran simulations for the cycles described in Relative Humidity and Load Cycles section.The voltage cycle, along with the resulting current density and average MEA temperature dynamics are shown in Fig. 4 for the different operating temperatures.As seen in Fig. 4a, the performance (in terms of current density) deteriorates with increasing temperature.This can be attributed to the fact that the simulations are run under relatively dry conditions (RH = 50%).Consequently, increasing the operating temperature results in an observable increase in the ohmic losses due to membrane dehydration.In addition, we see that the current density dynamics show faster response with increased temperature due to the fact that the ionomer water content reaches the equilibrium value more rapidly.In particular, note that at 80 • C the current density reaches the steady state value in about 10 seconds after each step in the cell voltage (e.g., see the dotted green line from around the 25 seconds point to the 35   around 55 seconds).This rather slow increase in the current density at lower temperatures can be attributed to ionomer humidification due to the electrochemical water generation.The figure also shows the average MEA temperature in each case (Fig. 4b), the dynamics of which follow the current density dynamics closely.
Fig. 5 depicts the resulting water content dynamics in the ionomer phase across the MEA for the different operating temperatures.The changes in the ionomer water content due to voltage cycling are more significant at 60 • C (compare Fig. 5a, Fig. 5b, and Fig. 5c).As mentioned above, this is due to continued humidification of the ionomer phase with electrochemical water generation at lower temperatures.As shown in Membrane Mechanical Response subsection below, such a significant hydration and dehydration of the ionomer induced by electrochemical load cycling results in considerable mechanical stress development in the membrane.One final point to notice from the figure is the considerable impact that EOD has on the water dynamics in the ionomer.In particular, we see that the water content in the anode CL demonstrates a response to the voltage cycle that is characterized as "nonminimum phase" in the system dynamics literature. 37With a step down in cell voltage, current density increases, which results in a higher water removal rate from the anode side due to the EOD.This in turn causes dehydration of the ionomer in the anode CL (e.g., as seen by the downward spikes in the dashed red line at around 30, 90 and 150 seconds in Fig. 5).As time goes on, increased water generation due to the electrochemical reaction results in ionomer hydration on the cathode side (the gradual rises in the dotted green lines in Fig. 5 starting at 30, 90 and 150 seconds).The water generated on the cathode diffuses through the MEA, consequently replenishing the anode side of the membrane.The reverse process occurs during the voltage step up, where initial decrease in EOD results in an increase of water content on the anode side due to the increased effect of back diffusion (e.g., the upward spikes in the dashed red line at around 60, 120 and 180 seconds in Fig. 5).This excess water, however, is eventually removed by desorption into the anode CL pore space.
As a final case study, the voltage and current density cyclings are compared to investigate the effects of different load changing schemes.The conditions for this case study are presented in Table I.The resulting power dynamics are shown in Fig. 6 along with the average water volume fraction across the MEA.The figure shows that the power dynamics demonstrate faster response during current cycling (dashed red vs. solid blue lines in Fig. 6).Nevertheless, significant oscillations are observed during current density step up, which, in some cases, result in voltage reversal due to severe dehydration of the anode side (e.g., the downward spike in the dashed red line at 30 seconds in Fig. 6a).This is the result of the step change in EOD corresponding to the change in current density.In particular, note that the current density changes rather smoothly during voltage cycling (as can be seen from the solid blue line in Fig. 6a).This enables the anode CL to hydrate  with back diffusion of water.However, during the current cycling, there is a step change in EOD, which does not allow enough time for the water back diffusion to replenish the dehydrated anode CL.As such, a significant voltage undershoot is observed during current step ups.This phenomenon can also be observed in Fig. 6b, where the water content in the anode CL goes through more severe changes during current cycling compared to voltage cycles.The voltage undershoot is more significant during the first cycle, which results in voltage reversal.This type of behavior is also reported elsewhere in the literature. 38s for the membrane and cathode CL, voltage and current cycling seem to induce very similar changes in the water dynamics (Fig. 6b,  6c, and 6d).The effects of such behavior on the mechanical response of the membrane are explored in the next subsection.
Membrane mechanical response.-Inthe following, we import water content and temperature data, calculated by the transport model, into the mechanical model to calculate the corresponding stresses.The three types of cycles considered earlier, are revisited here: RH cycling, voltage, and current load cycling.Previous work has shown that inplane stresses are the largest and most critical components of stress in the membrane during fuel cell operation. 6,7,39,40Therefore, contour plots of the 2D in-plane stress distribution in the membrane are examined to determine the locations of the maximum stress magnitudes (tension/compression).From these contour plots, "observation points" are chosen for temporal stress plots for the different load cases.
RH cycling.-For the RH cycle described earlier, the spatiotemporal water volume fraction distribution and average temperature in the MEA are imposed at the nodes of the mechanical model to evaluate the stresses in the membrane.Results are shown for various times in Fig. 7a.
The maximum in-plane compressive and tensile stresses in the membrane occur at 24 (at the start of holding at wet state during the first RH cycle) and 97 (at the end of holding at dry state during the second RH cycle) seconds, respectively.From the contour plot snapshots at those times (Fig. 7b), the location of the maximum compressive in-plane stress is at the right end of the domain under the land and the maximum residual tensile in-plane stress is at the left end under the flow channel.The locations of maximum stress are in general agree-ment with previous studies. 7,10,39,40This shows consistency between the previous model and the current model for RH cycling, where water generation, EOD and reaction kinetics are not present.However, some of the details of the water content, and therefore, stress distribution are quite different than in the previous model.
To study the stress evolution, one node in each of these regions is taken as an observation point (Fig. 7) and the temporal in-plane water content (Fig. 8a) and stress profiles (Fig. 8b) are generated.Fig. 8b shows the maximum compression that occurs under the land, and the maximum tension that occurs under the channel.Even though the observation points are both chosen on the cathode side of the membrane, similar trends are observed on the anode side as well due to the rather symmetric distribution of water across the MEA thickness.
Looking at the water volume fraction dynamics (Fig. 8a), a longer delay in dehydration is noticeable for the observation point with maximum compressive stress.This can be attributed to the fact that this point is under the land, where the water has a longer path to get to the channels.At the start of the second hydration phase (60 seconds), the membrane water content starts to increase during the hydration process (Fig. 8a), causing the membrane to swell (Fig. 8b).However, the mechanical constraint from the GDLs and bipolar plates prevents expansion in the through-plane direction, generating a comparatively large compressive stress, on the order of 7-9 MPa.Then, as the water volume fraction levels off, the compressive stress starts to relax.Therefore, the in-plane stresses in this region alternate between tension and compression.The kink at the beginning of the dehydration process (about 80 seconds in Fig. 8b), corresponds to the slight increase in water volume fraction at that time (explained in RH cycling subsection) resulting in the increase in the compressive stresses.The dehydration and hold at low humidity (82 -100 seconds) leads to shrinkage of the membrane and residual tensile stress, on the order of 5-7 MPa (about 60 and 97 seconds in Fig. 8).The magnitude of the tensile stress is slightly higher under the channel suggesting that this region is marginally more prone to failure during RH cycling.While these results are qualitatively similar to previous simulations, 41 the previous membrane model did not consider the above-mentioned comprehensive transport phenomena, which limits its application to RH cycling without the electrochemical processes.Voltage cycling.-Thetransport model is used to simulate the voltage cycles described earlier and the resulting water volume fractions and average temperatures are imposed on the MEA in the mechanical model to capture the corresponding stress profile.
As before, 2D contour plots of in-plane stresses are generated as illustrated in Fig. 9.The figure shows that for the simulated conditions, the maximum compressive stress consistently develops under the land and on the cathode side of the membrane.The location for maximum tensile stress, on the other hand, changes for different operating temperatures.The maximum tensile stress is under the flow channel on the cathode side of the membrane at 60 • C, under the land on the cathode side of the membrane at 70 • C, and under the flow channel on the anode side of the membrane at 80 • C.
To help explain the above observation, Fig. 10 demonstrates the temporal evolution of water volume fractions at six different nodes in the membrane.The nodes are selected at the top, middle and bottom of the membrane on both the land and channel sides.The figure shows that the operating temperature has a profound impact on the water distribution in the MEA.To further explain this, it is helpful to describe the notion of "ignition" in PEM fuel cells, first proposed in the works of Benziger et al. 42,43 and further explored by others. 44In these references, they note that during low humidity operation, the membrane is auto-humidified by taking up the electrochemically generated water during ORR.This increases the membrane ionic conductivity, which in turn results in higher current generation and the cell is said to be in an "ignited" state.On the other hand, if the membrane is too dry, the current density will be retarded to the point that the amount of water generated is not enough to hydrate the membrane leading to an "extinguished" state of the cell.For the cycling conditions reported here, the notions of ignited and extinguished states of the cell are not directly applicable, since the gas feeds are humidified.Nevertheless, the underlying concept remains critically important.
The results in Fig. 10 can be explained in the context of ignition.At 60 • C (the first column in Fig. 10), the relatively low temperature results in higher RH values throughout the domain.Nevertheless, the RH remains slightly higher under the land for reasons explained earlier.Therefore, with a step down in voltage, the current density (and therefore the water volume fraction in the membrane) responds quickly under the land.Meanwhile, the lower RH values under the channel mean that the membrane has to undergo an auto-humidification process, which slows down the response.However, the initial membrane hydration level is enough to "ignite" the auto-humidification process.This means that the membrane under the channel will eventually reach a hydration level on par with that under the land, albeit at a slower rate.Therefore, a rather uniform in-plane water distribution is obtained after reaching steady state at 60 • C. When temperature is increased to 70 • C (the second column in Fig. 10), the RH decreases throughout the domain, which disturbs the uniform distributions observed at the lower temperature.The conditions under the land still remain favorable for the auto-humidification process.However, this is no longer the case under the channel.The result is a rather significant in-plane gradient in the water volume fraction in the MEA.In particular, the water content under the channel does not seem to vary significantly, while that under the land still goes through noticeable transients.This explains the fact that the highest in-plane tensile stress at 70 • C is observed under the land.Finally, at 80 • C (the last column in Fig. 10) the RH remains low throughout the domain, "extinguishing" auto-humidification.Therefore, this condition again promotes a rather uniform water distribution in the in-plane direction and the location for the maximum in-plane tensile stress is observed to be under the channel.The fact that this location happens to be on the anode side of the membrane can be attributed to minor dynamic variations in the water volume fractions during the cycles as well as the effects of EOD.Since the membrane under the land is more mechanically constrained than under the channel, the location of maximum in-plane compressive stress is consistently under the land, which is also where the water content is typically highest.Nevertheless, it is conceivable that different operating conditions, particularly those that result in reactant starvation (and therefore very small current generation) could alter this.
To investigate further, Fig. 11 illustrates the temporal stress evolutions for specific observation points taken in the maximum in-plane stress (compressive/tensile) regions shown in Fig. 9.The figure shows that the stress induced during voltage cycling at 60 • C is significantly higher than at the other two temperatures, especially in the maximum tensile stress region (the bottom row of Fig. 11).In particular, the in-plane maximum tensile stress at 60 • C is 2 and 8 times higher than that at 70 and 80 • C, respectively.This can be attributed to two factors.First, the transport phenomena lead to a significant difference between the maximum and minimum water content of the membrane during hydration and dehydration processes, resulting in increased strain rates.Second, the lower temperature makes the membrane stiffer.This observation is consistent with previous experimental results. 32o study the changes in temporal stress evolution at the core of the membrane due to changes in operating conditions, stresses are evaluated at observation points for all three temperatures as shown in Fig. 12.One observation point is chosen under the channel, since different studies 32,39,40,45 suggest that membranes are more prone to mechanical failure there.The other point is chosen under the land, where the in-plane stresses are mostly compressive, which may tend to influence the chemical decomposition of ionomers. 15,46Figs.12c and  12d show that the in-plane stresses are higher at lower temperatures for the simulated conditions (solid vs. broken lines).As mentioned before, this can be attributed to the fact that these simulations are carried out under relatively dry conditions and the operating temperature has a profound impact on membrane water uptake.Furthermore, we observe that when the water volume fraction dynamics are the same between two conditions, the resulting stress profiles are similar, as expected.This can be seen, for instance, in the case of in-plane stresses under the channel at 70 and 80 • C (broken lines in Figs.12a and 12c).
In summary, the results of this subsection illustrate the critical role of transport phenomena in determining the location of maximum in-plane stresses, both tensile and compressive.The contrast that is observed between the location of maximum in-plane tensile stress during load cycling and that during RH cycling (see Figs. 7 and  9) highlights the importance of accurately taking transport into account while designing accelerated stress tests.8][49][50] Therefore, while the RH cycling results suggest that membrane cracks and/or Under the Channel Under the Land pinholes might form under the channel first (Fig. 7b), the load cycling results point to the possibility of initial crack and/or pinhole formation under the land depending on the operating conditions (Fig. 9 70 • C).In contrast, the region under the land, where the maximum inplane compressive stress is always found, might be considered more prone to chemical failures, since it has been reported in the literature that chemical degradation tends to increase under compression. 46herefore, a comprehensive framework to investigate the impact of transport phenomena on mechanical behavior of an operating fuel cell is needed and the one-way coupling between the transport model and mechanical model is a first step toward developing such a framework.
Voltage vs. current cycling.-Inthis subsection, the water content profiles and stresses are predicted in response to i) voltage cycling and ii) current cycling under the same operating conditions (explained in Voltage and current cycling section).
Fig. 13 shows the stress contours at selected times.The maximum compressive and tensile stresses during voltage cycling occur under the land and channel, respectively, both on the cathode side (Fig. 13a).In the case of current cycling, the maximum in-plane compressive stress is under the land near the cathode and the maximum tensile stress develops under the channel and close to the anode (Fig. 13b).
The temporal profiles of the in-plane stresses in the regions of maximum stress (indicated in Fig. 13) are illustrated in Fig. 14.Fig. 14a shows significantly different dynamic profiles in the tension region for the two loading schemes.This difference can be attributed to the fact that maximum tension develops on the cathode side of the membrane during voltage cycling, whereas it appears on the anode side of the membrane during current cycling, and that the transients on the anode side of the membrane are severely affected by the EOD.In addition, the anode and cathode sides of the membrane typically experience opposite trends during load cycles due to the combined effects of diffusion and EOD.In particular, water depletion on the anode side usually results in better hydration of the cathode (through EOD) and vice versa (through back diffusion).The dynamic profiles in the compressive region, on the other hand, are very similar for the two loading schemes due to the fact that maximum compression is in the location for of these load cycles (Figs.14b and 14d).The small kinks in both the tension and compression regions at about 60 and 120 seconds can be explained as follows.As the electrical load decreases at these times, so does the EOD and water Therefore, the cathode side of the membrane gets depleted of water quickly by back diffusion.However, some liquid water exists in cathode CL under the land, which starts to evaporate slightly later resulting in a slight increase in membrane water at that location before its eventual dehydration.The notably changes in EOD during current cycling may provide an avenue for faster crack initiation on the anode side of the membrane.Therefore, it might seem tempting to suggest that voltage cycling would be preferable from a purely mechanical point of view.However, the presented results are far from being quantitatively conclusive given that only one operating condition is tested here.Furthermore, the synergistic effects of mechanical and chemical degradation have to be taken into account provide a comprehensive picture of the mechano-chemical stressors before a conclusion can be made in this regard.Studies like this should help pave the toward better understanding of such synergistic effects and provide insight about the optimal operating conditions to enhance the cell's durability.

Summary Conclusions
This study creates a tool improved understanding of the PFSA membrane behavior by coupling a comprehensive transport model with viscoelastic-plastic membrane mechanical model.Simulations are carried out for an RH cycle as well as several load cycles under different operating conditions.The RH cycling results indicate that a cycle based on US DOE guidelines can result in significant compressive and tensile stress development in the membrane.The load cycling cases, which could not be simulated by the original membrane mechanical model, showcase the important role of various transport phenomena on the stress distribution in the membrane.
The simulations focus on the locations of maximum in-plane stresses under RH and load cycling, which correspond to potential damage initiation sites.particular, we found that the location of maximum in-plane stress can vary significantly with the operating temperature.Under voltage cycling at 70 • C, the maximum in-plane tensile stresses are found under the land, which could make that site more vulnerable to fatigue failure.In contrast, at other tested temperatures the maximum in-plane stress during voltage cycling occurs under the channel.In all cases, the highest in-plane compression stress is found under land, particularly near the cathode.Since stresses and types (compressive/tensile) at key locations are observed differ due to in the temperature under otherwise similar operating temperature can effectively influence the mechanical durability by controlling the spatiotemporal water distribution besides changing the material properties.The comparison of voltage and current cycling has shown significantly different dynamic profiles in the region of maximum tensile stress for the two loading schemes, as the maximum tension develops on the cathode side during voltage cycling and it appears on the anode side during current cycling.From a purely mechanical point of view and based on the chosen simulation conditions, voltage control might seem preferable due to the transport phenomena and the corresponding stress development associated with immediate changes in EOD in response to current steps.
The simulation results confirm the need to account for complex and coupled transport phenomena in studying the mechanical behavior and failure mechanisms of the membrane.These results are first steps toward improved understanding of various synergistic mechanical failure mechanisms that affect the membrane in an operating fuel cell.Such understanding can be utilized to devise degradationconscious control algorithms that aim to enhance the lifetime while meeting the performance requirements.One of the future directions is to add the chemical degradation capability into the mechanical framework of this work.

Figure 1 .
Figure 1.Modeling domain used in this work.

Figure 2 .
Figure 2. Schematic of the model coupling framework: (a) prior implementation of the mechanical model with a simple transport model (i.e.single-phase, no water and heat generation, no electro-osmosis or thermo-osmosis in the membrane), and (b) current implementation using the comprehensive transport model.

Figure 4 .
Figure 4. Voltage cycling: (a) Voltage cycle and current density response at different operating temperatures, (b) average MEA temperature.

Figure 6 .
Figure 6.Voltage and current cycling comparison: (a) Power density dynamics, and average water volume fraction in (b) anode CL, (c) membrane, (d) cathode CL.

Figure 7 .
Figure 7. RH cycling results: (a) spatial in-plane stress distribution of membrane at different phases of the cycle, (b) in-plane stress distributions of membrane at the time of the maximum compressive and maximum tensile stresses.

Figure 8 .Figure 9 .Figure 10 .Figure 11 .
Figure 8. RH cycling results: (a) temporal water volume fraction evolution at observation points in tension and compression regions, (b) temporal in-plane stress evolution at observation points in tension and compression regions.

Figure 12 .
Figure 12.Voltage cycling results: temporal water volume fraction evolution (a) under the channel and (b) under the land, and temporal in-plane stress evolution (c) under the channel and (d) under the land at observation points (taken at membrane's core) at different operating temperatures.

Figure 13 .
Figure 13.Spatial in-plane stress distributions of the membrane at the time of maximum compressive and tensile stress occurrence during (a) voltage and (b) current cycling.

Figure 14 .
Figure 14.Temporal water volume fraction evolution at observation points (shown in Fig. 13) in (a) tension and (b) compression and in-plane stress evolution at observation points in (c) tension and (d) compression regions during voltage and current cycling.