# An Electrohydrodynamic Equilibrium Shape Problem for Polymer Electrolyte Membranes in Fuel Cells

Sven-Joachim Kimmerle, Peter Berg, Arian Novruzi

► To cite this version:

Sven-Joachim Kimmerle, Peter Berg, Arian Novruzi. An Electrohydrodynamic Equilibrium Shape Problem for Polymer Electrolyte Membranes in Fuel Cells. 25th System Modeling and Optimization (CSMO), Sep 2011, Berlin, Germany. pp.387-396, 10.1007/978-3-642-36062-6_39. hal-01347560

HAL Id: hal-01347560

https://inria.hal.science/hal-01347560

Submitted on 21 Jul 2016

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# An Electrohydrodynamic Equilibrium Shape Problem for Polymer Electrolyte Membranes in Fuel Cells

Sven-Joachim Kimmerle<fnref n="1" />*, Peter Berg<fnref n="2" />, and Arian Novruzi<fnref n="3" />

**Abstract.** We present a novel, thermodynamically consistent, model for the charged-fluid flow and the deformation of the morphology of polymer electrolyte membranes (PEM) in hydrogen fuel cells. The solid membrane is assumed to obey linear elasticity, while the pore is completely filled with protonated water, considered as a Stokes flow. The model comprises a system of partial differential equations and boundary conditions including a free boundary between liquid and solid. Our problem generalizes the well-known Nernst-Planck-Poisson-Stokes system by including mechanics. We solve the coupled non-linear equations numerically and examine the equilibrium pore shape. This computationally challenging problem is important in order to better understand material properties of PEM and, hence, the design of hydrogen fuel cells.

**Keywords:** Nernst-Planck-Poisson-Stokes system, Free boundary problem, Equilibrium shape, Fluid-structure interaction, Polymer electrolyte membrane, Proton exchange membrane fuel cell, Nafion, Mechanical deformation of pores, Ohmic interface resistance

## 1 Introduction

Fuel cells running at low temperature provide a possibility for the electrification of the power train in automotive devices. Proton exchange membrane fuel cells (PEMFC) are based on hydrogen as fuel and do not rely on the use of fossil combustible material, which is likely getting increasingly scarce and expensive in the future. PEMFC allow to produce electric current by only emitting water as a byproduct and no carbon dioxide. In a hydrogen fuel cell, hydrogen enters

<footnote n="1">¹ Universität der Bundeswehr München, Institut für Mathematik und Rechneranwendung, Werner-Heisenberg-Weg 39, 85577 Neubiberg, Germany
{sven-joachim.kimmerle@unibw.de}</footnote>

<footnote n="2">² Faculty of Science, University of Ontario Institute of Technology, 2000 Simcoe Street N, Oshawa, ON, L1H 7K4, Canada
{peter.berg@uoit.ca}</footnote>

<footnote n="3">³ Department of Mathematics and Statistics, University of Ottawa, 585 King Edward Avenue, Ottawa, ON, K1N 6N5, Canada
{novruzi@uottawa.ca}</footnote>

<footnote>* Corresponding author</footnote>
into the fuel cell at the anode (negative), while oxygen flows in at the cathode (positive). The design of the PEM allows for control of the potentially explosive chemical reaction between hydrogen and oxygen. The reaction product, just water, leaves the fuel cell predominantly through the cathode outlet. The electric load is applied between anode and cathode, which closes the electric circuit. Within the PEM, the protons migrate from anode to cathode.

The PEM consists of a polymer with a nanopore structure. Negatively charged sulfonic acid groups ($SO_3^−$) at the walls of the nanoscale pores allow the dissociation of protons ($H_3O^+$) in the presence of water. A typical material for a PEM is Nafion, a perfluorosulfonic acid ionomer. A hydrated Nafion membrane exhibits a hydrophobic elastic backbone and hydrophilic pores, that are filled with protons and water molecules [1]. The precise morphology of Nafion on nanoscales remains a controversially discussed issue. We follow the widely accepted approach by Schmidt-Rohr and Chen [2], that a Nafion membrane consists mainly of parallel cylindrical channels surrounded by hydrophilic side chains.

Within the production of PEMs, a possibility is to press together several layers of thinner membranes. An increased ohmic resistance between interfaces of two joint PEM is observed in experiments. Our objective is to establish a mathematical model, describing the charged fluid flow and the change of the morphology of the pore due to mechanical deformations, that allows to understand better material properties of PEM. In particular we are interested in finding an equilibrium pore shape and, finally, in explaining this ohmic resistance.

The behaviour of water in a nanochannel may be very different compared to that of bulk water. The physics of water in a confined small channel depends significantly on the type of surface, i.e. whether the interface is hydrophilic or hydrophobic, and on the presence of surface charges [3]. There are several approaches, adapted to different scales, in order to study the charged fluid within PEM pores. On a microscopic level these are mainly molecular dynamics and Brownian motion, on a mesoscopic level there are continuum models, e.g. the Nernst-Planck-Poisson-Stokes (NPPS) system. We follow a continuum approach that is applicable in our situation since the Debye length is small, see [4].

We generalize the well-known NPPS model [5] by fluid-structure interaction between the charged fluid flow and the elastic wall of the channel. Previous electrohydrodynamic models, e.g. [4–6], do not incorporate the coupling to the mechanical displacement field. The equations for the flow stated by Castellanos [5] and analysed by Schmuck [7] are more general. They deal with the full non-stationary Navier-Stokes flow and allow also negative charge carriers, but the liquid domain is fixed. The equations in [4, 6] represent the first-order approximation for a stationary version of our model without mechanical deformations. The significant effects of the radial variation of system parameters in this situation has been underlined in [4]. However, we work with a varying viscosity instead of a no-slip surface as in [4], modelling essentially the Stern layer in the pore. Furthermore, we consider slightly different boundary conditions (b.c.), e.g. we work with homogeneous Neumann b.c. (13), (14) for the chemical potential on in-/outlet and on the interface instead of Neumann b.c. for the proton concen-
<figure><img src="image_2.png" /><figcaption>**Fig. 1.** Domains: *L*, the liquid, *S* the solid, *I* the interface, *O*₁ inlet (anode), *O*₂ outlet (cathode). Unknowns: in *L*: **u**, *p*, *φ*, *c*; in *S*: **U**. Dashed part of the interface: initial pore shape, part of the interface with continuous lines: equilibrium pore shape.</figcaption></figure>

tration, which are only an approximation of our thermodynamically consistent b.c. A free boundary problem for an interface between elastic solids and fluids is described in [8, 9]. These models are similar to our model, when any coupling of the flow to electric field and concentration is neglected. A similar setting as ours has been considered by [10], but the authors focus on possible equilibrium configurations depending on the hydration of the pore and they make the assumption, arguably too strong in nature, of a constant proton concentration. Secondly, we study a stationary fluid flow in a fully saturated PEM.

To the best knowledge of the authors, a full problem involving both charged fluid flow and fluid-structure interaction has not been considered yet. We derive partial differential equations that we state in Sect. 2. We focus on a single cylindrical channel that joins another cylindrical channel from the adjacent PEM layer. Both channels are assumed to be completely filled with protonated water. In §2.1 we define our geometry and introduce the crucial physical quantities. The model is solved numerically by a commercial finite element software in Sect. 3. We discuss the impact of our results in the last part, Sect. 4, of this short paper.

## 2 An Elasto-Electrohydrodynamical Model for Polymer Electrolyte Membranes

### 2.1 Geometry and Relevant Physical Quantities

We assume that a single Nafion pore consists of several nanochannels and that the length of nanochannels is larger than their diameter *d*, typically a few nanometre. Instead of solving the full problem for many channels, we examine the situation around the region where one channel from one PEM layer meets one other channel from another PEM layer. We consider a nanochannel segment of length *l* including the interface (see Fig. 1). The joint channel extends further on in both directions and is assumed to be connected to other pores at both
ends. Let $C = L \cup S \cup I$ denote the whole domain with the open domains $L$ and $S$, modelling the liquid channel and the solid Nafion backbone (of both PEM). $I = \partial L \cap \partial S$ denotes the interior free boundaries (or interface) between $S$ and $L$ and $\partial C$ are all outer boundaries of $C$. We denote the inlet by $O_1$ and the outlet by $O_2$. Here the outer boundaries are considered as fixed boundaries, since in our model the channel extends further out of the considered domain $C$. Free boundaries may move with a normal speed $\omega$. We make the following convention: $\nu$ denotes the outer normal on $I \cup \partial C$, pointing on $\partial C$ outside $C$ and pointing on $I \setminus \partial C$ always from the liquid into the solid. Consequently $\partial_{\nu} f = \nabla f \cdot \nu$ denotes the derivative in direction of the outer normal for a differentiable function $f$. The mean curvature (times a factor 2) of a surface is the tangential divergence of the outer normal, i.e. $\kappa = \nabla_{\tau} \cdot \nu$, being positive on $I$, if $I$ is convex (seen from $L$). $\tau_1$ and $\tau_2$ denote the two orthonormal tangential vectors on $I$.

In this study, we suppose that the channel is completely filled with protonated water. Counter-ions, i.e. protons, are considered to be the only charge carriers in our model. On the interface $I$ between liquid and solid we have negatively charged sulfonic acid groups, that are modelled by the negative surface charge density $\sigma_C$, see [6]. We like to solve for the velocity field $u$, the pressure $p$, the proton concentration $c$, and the electric potential $\phi$ in the liquid $L$, and we are looking for the mechanical displacement field $U$ in the solid $S$.

## 2.2 Modelling of Strains and Stresses

Mechanical strains and stresses are defined w.r.t. a reference configuration, where the system is free of strains and stresses or the stresses are at least known *a-priori*, either from experiments or from calculations (as it is possible e.g. in the case of a symmetrical geometry). In order to define a mechanical reference configuration, we consider at first another case where we have a straight circular channel with a fixed diameter cut out of a box of solid Nafion. An outer pressure $p^*$ is exerted on the liquid channel without flow. At the interface the well-known Young-Laplace law $p^* + \gamma\kappa^R = p^R$ holds. Here $\gamma$ is the surface tension, $\kappa^R = 2/d$ is the mean curvature of the straight channel, and $p^R$ is the reference pressure in the solid for a straight cylinder.

As a reference configuration we consider the case of two joint straight cylindrical channel segments with the same radius and parallel axes. These channels are shifted by a fixed offset $s$ between the axes. The offset is the distance between the circle centres at the interface plane. The reference pressure $\bar{p}$ in the solid, corresponding to the situation of two channels with an offset, is varying in space. Namely we have $\bar{p} = p^R = p^* - 2\gamma/d$ on the part of $I$ that belongs to the cylinder barrels, while $\bar{p} = p^*$ on the part of $I$ near the offset that belongs to the cylinder covers. Stresses and strains are to be formulated w.r.t. this reference configuration. The geometry within the reference configuration will be denoted by $S_0$, $L_0$ and $I_0$, while in the current (actual or deformed) configuration we write $S$, $L$ and $I$. We assume that we may neglect here inelastic deformations that are due to changes of the chemical composition.
We consider mechanical deformations and displacements $U_0$ in the reference configuration on $S_0$ and then as $U$ in the current (actual) configuration on $S$ w.r.t. the reference configuration. We consider a material point $X$ in the reference configuration, whose location at time $t$ is given by $x$ in the current configuration. $\chi(t, X) = x$ is called the deformation of material points. The displacement of a material point is defined by $U_0(t, X) := \chi(t, X) - X$. We assume that we may invert the deformation $\chi$ at any time w.r.t. material points $X$ and, thus, we express the displacement w.r.t. the current configuration, $U(t, x) := U_0(t, \chi^{-1}(t, x))$. It is convenient to consider fluid flow in the current configuration, i.e. in Eulerian coordinates, while mechanical displacements in a solid are considered in the reference configuration, i.e. in Lagrangian coordinates. We may work in good approximation with linear elasticity instead of nonlinear elasticity, which would be needed for a full description of large deformations of polymers, see e.g. [11]. Therefore solid stresses are represented by the Cauchy stress tensor $\sigma_S(\nabla U) = (-\bar{p} + \lambda_S \text{tr}(\nabla U))\mathbf{1} + \mu_S e(\nabla U)$, where $\mathbf{1}$ is the unit tensor, $\text{tr}$ denotes the trace, and $e(t) = t + t^\top$ is the symmetrization of a tensor $t$ times 2. The Lamé constants $\lambda_S, \mu_S$ are given material parameters. In the liquid we deal with a Newtonian fluid and the stress tensor reads $\sigma_L = -(p + \frac{2}{3}\mu_d \text{tr}(\nabla u))\mathbf{1} + \mu_d e(\nabla u)$, with $\mu_d$ being the dynamic viscosity.

## 2.3 Governing Equations

Now we may state our mathematical problem. For details of its derivation from first principles and the choice of thermodynamically consistent constitutive relations, see [12]. We remark that we assume that the dynamic viscosity $\mu_d$, the electric permittivity $\epsilon_r$, and the diffusion coefficient of protons, $D$, may vary in space [4]. Furthermore, we have as given data the outer pressure at in-/outlet $p_0^{(1)}/p_0^{(2)}$, a given displacement $g_0$ on the boundary $\partial S_0 \setminus I_0$, and a constant external electric field $E_{ext}$.

Our problem consists of Stokes equations for $u$ and $p$,

$$ [\text{Momentum balance}] \qquad -\nabla \cdot (\mu_d e(\nabla u)) + \nabla p = -Fc\nabla\phi \quad \text{in } L, \quad (1) $$

$$ [\text{Incompressibility}] \qquad \nabla \cdot \mathbf{u} = \mathbf{0} \qquad \text{in } L, \quad (2) $$

where $Fc\nabla\phi$ accounts for electro-osmotic pressure, $F$ being Faraday's constant. This is complemented by the boundary conditions

$$ [\text{Normal pressure bal.}] \qquad \mu_d e(\nabla u)\nu - p\nu = -p_0^{(i)}\nu \qquad \text{on } O_i, i = 1, 2, \quad (3) $$

$$ [\text{Tangential moment. bal.}] \qquad u \cdot \tau_j = 0 \qquad \text{on } O_i, i, j = 1, 2, \quad (4) $$

$$ [\text{Momentum balance}] \qquad u = \partial_t U \qquad \text{on } I. \quad (5) $$

The mechanical displacement field $U_0$ is determined by the following problem of linear elasticity, formulated in the reference configuration,

$$ -\nabla \cdot \sigma_S(\nabla U_0) = \mathbf{0} \qquad \text{in } S_0, \quad (6) $$

$$ U_0 = g_0 \qquad \text{on } \partial S_0 \setminus I_0, \quad (7) $$

$$ -\sigma_S(\nabla U_0)\nu_0 = -\mu_d e(\nabla u_0)\nu_0 + (p_0 - \gamma\kappa_0)\nu_0 \qquad \text{on } I_0, \quad (8) $$
where the last line is the Young-Laplace equation, a pressure balance. $u_0$, $p_0$ and $\kappa_0$ denote $u$, $p$ and $\kappa$ expressed in the reference configuration. For the electric potential $\phi$, we solve the Poisson equation with Neumann b.c.

$$ [\text{Electrostatics}] \qquad -\varepsilon_0 \nabla \cdot (\varepsilon_r \nabla \phi) = Fc \qquad \text{in } L, \quad (9) $$

$$ [\text{External electric field}] \qquad -\partial_{\nu}\phi = (-1)^{i}E_{\text{ext}} \quad \text{on } O_{i}, i = 1, 2, \quad (10) $$

$$ [\text{Interface charges}] \qquad -\varepsilon_0 \varepsilon_r \partial_{\nu} \phi = -\sigma_c \qquad \text{on } I, \quad (11) $$

where $\varepsilon_0$ is the vacuum permittivity. For the proton concentration $c$, we solve the Nernst-Planck equation

$$ -\nabla \cdot \left( D \nabla c + \frac{F}{RT} Dc \nabla \phi \right) + u \cdot \nabla c = 0 \qquad \text{in } L, \quad (12) $$

$$ \partial_{\nu}c + (-1)^{i+1} \frac{F}{RT} E_{\text{ext}} c = 0 \quad \text{on } O_i, i = 1, 2, \quad (13) $$

$$ \partial_{\nu}c + \frac{F}{RT\varepsilon_{0}\varepsilon_{r}}\sigma_{c}c = 0 \qquad \text{on } I, \quad (14) $$

the last two lines representing homogeneous Neumann b.c. for the chemical potential $RT \ln(c/\bar{c}) + F\phi$, $R$ denoting the universal gas constant and $T$ the temperature. In (12), the first term represents the diffusion of protons, the second term the migration of protons within the electric field and the third term is the advection due to the moving fluid. The normal velocity $\omega$ of the free boundary $I$ is determined by

$$ [\text{Normal momentum balance}] \qquad \omega = \partial_t U \cdot \nu \quad \text{on } I. \quad (15) $$

In this study, we are looking for stationary solutions of (1) – (15). Consequently, we neglect the time-derivatives in (5) and (15). Otherwise, we would have to prescribe an initial condition for $I$, and consider time-dependent domains and boundaries. The equilibrium pore shape corresponds to a displacement on the interface, $U_0|_I$, that is constant in time. It is designated by (8).

We determine uniquely the electric potential $\phi$ by imposing $\phi = \bar{\phi}$ at the centre point $P$ of the inlet $O_1$. By solving a nonlinear eigenvalue problem, a typical value $\bar{c}$ for the proton concentration is derived in [6] that corresponds to setting $\phi = \bar{\phi}$ at the centre of a straight cylindrical channel in the case of constant permittivity. The arbitrariness of $\bar{\phi}$ or $\bar{c}$ corresponds to the gauge invariance of the electric potential.

We emphasize that the b.c. (10) is consistent with global electroneutrality, meaning that we rule out that the electric field extends into the solid and, hence, ions would cross the interface. This requires that the hydronium charges ($H_3O^+$) in the liquid and the negative ions ($SO_3^-$) on the boundary balance $\int_L Fc + \int_I \sigma_c = 0$. Together with (9), (11), and Gauss’ theorem we find $\int_{O_1 \cup O_2} \varepsilon_0 \varepsilon_r \partial_{\nu} \phi = 0$. The last equation is guaranteed by (10), a choice among many others.
# 3 Numerical Solution

For any solution approach, numerical or analytical, it is crucial to non-dimensionalize the system, and to split up the equations and, consequently, the quantities according to their different scales. We remark that by non-dimensionalization of the system we find the relevant dimensionless coupling parameters, e.g. a parameter resulting from non-dimensionalization of (9) is the Debye length of the nanochannel. The dimensionless parameter Pe arising in (12) is the Peclet number describing the relation between advection and diffusion for mass transfer. Since $Pe \ll 1$ we could be tempted to neglect the term $\mathbf{u} \cdot \nabla c$ in the Nernst-Planck equation. This observation motivates the first-order approximation that is considered and compared with the full model in [12]. However, for the flux of protons, $j_+ = -D(\nabla c + \frac{F}{RT}c\nabla\phi) + cu$, that enters significantly into the ohmic resistance within the nanochannel, we would obtain zero within this first-order approximation. Hence, it is crucial to consider the full equation (12). Sorting the r.h.s. of the PDE by their scales suggests to introduce an internal electric potential $\phi_0$ (with a non-zero b.c. only on the interface), an external electric potential $\phi_1$ (with a non-zero b.c. only at in-/outlet) and a remainder potential $\phi_2$. Corresponding to $\phi_0$ we introduce a concentration $c_0 = \bar{c}\exp(-F/(RT)\phi_0)$, yielding a remainder concentration $c_1$ with a well scaled PDE. Finally for the pressure it is convenient to replace $p$ by $q := p - RTc_0$, since $q = 0$ implies $u = 0$.

Since the analytical solution of our non-linear coupled system is quite ambitious we focus on a numerical solution. Schmuck's analytical results [7], for a model without coupling to linear elasticity, rely on the parabolic structure of the non-stationary equations for $\mathbf{u}$ and $c$ and cannot be transferred to our stationary case directly. We remark that we have a two-sided fluid-structure interaction. The charged fluid flow influences by means of (8) the deformation of the elastomer membrane, while a narrowing (e.g. a complete closing of the channel) or widening of $L$ has large influence on velocity or pressure.

## 3.1 Description of the Numerical Algorithm

Our numerical algorithm solves for the state variables as well as for the free interface. It should be emphasized that it is a non-trivial matter in which order the coupled equations are solved so as to obtain fast convergence and numerically stable results. We make use of the ALE (Arbitrary Lagrangian Eulerian) method, i.e. we discretize the original geometry in Lagrangian coordinates. Our numerical algorithm is built up in the following way:

a) Considering the reference configuration, we initialize the mesh for $L_0$, $S_0$ and $I_0$. Formally, we set $U_0 = 0$, $u_0 = 0$, and $q_0 = 0$, where $q_0$ is $q$ transformed into the reference configuration.

b) We compute smoothed normal vectors on the boundaries by extending the normal vector field $\nu_0$ by means of $\Delta\nu_0 = 0$ into $L_0 \cup S_0$, using the definition of $\nu_0$ on $\partial C_0 \cup I_0$ as boundary condition. This enables us to compute the mean curvature $\kappa_0$, even at corners that would become smooth instantly due to surface tension anyways.
c) We store $U_0^{(old)} = U_0$, then we solve the linear elasticity problem (6) – (8) for $U_0$ in $S_0$, using the present values for $u_0$ and $q_0$.

d) We update the geometry. The free boundary is moved by the mechanical displacement $U_0 - U_0^{(old)}$ on $I_0$, yielding the update of $I$, and hence of $S$ and $L$.

e) We solve the electrohydrodynamical system (1) – (5), (9) –(14) as follows.

(i) First we solve for $\phi_1, \phi_0$, and $c_0$, that do not depend on other variables.

(ii) We solve the remaining equations iteratively. We start solving for $u$ and $q$ simultaneously, and then for $\phi_2$ and $c_1$ simultaneously.

(iii) Then we reiterate e)(i) and e)(ii) until we have a suitable residual error.

f) If $\max_{I_0} \|U_0 - U_0^{(old)}\| < \text{err}$, err a prescribed error tolerance, or if a topological event has occured (e.g. closing of the pore), or a specified maximal number of iterations has been reached, we terminate our algorithm with the obtained numerical solution. Else we take the updated geometry as new reference configuration. Should the situation of a geometrically unsuitable mesh arise after the deformation, we remesh. We restart with step b).

For more details of the algorithm see [13]. Our algorithm has been implemented in the commercial finite element software COMSOL 3.4. The use of a parallel solver, like PARDISO, is crucial for efficient calculations.

## 3.2 Numerical Results

We show two plots, corresponding to two different choices for external parameters, namely, (i) a situation where the initial interface is close to an equilibrium (due to large $p_0^{(1)} - p_0^{(2)}$ and $p^*$ the surface tension term is negligible) (see Fig. 2, top), and (ii) a situation where the charged fluid flow is close to zero ($p_0^{(1)} - p_0^{(2)}$ small, $E_{ext}$ negligible) (see Fig. 2, bottom). All remaining data for our simulations is discussed and summarized in [12]. We emphasize that the database, e.g. for surface tension of protonated water and typical pressures, is thin.

# 4 Conclusions and Open Questions

We have stated a continuum model describing the charged fluid flow within nanochannels of PEM. The numerically accurate simulation of the differential equations allows further investigation of system characteristics. In our study [12], we focus on the ohmic interface resistance between two circular cylindrical pores, and analyze its dependence on system parameters. By splitting our model into two models, we consider the main effects separately. In this context, our model suggests that this interface resistance depends mainly on two factors: (a) the offset value (the distance between the pore centres at the intersection plane) and its stochastic distribution and (b) on the deformed pore shape due to the balance of elastic and electrohydrodynamic forces. Furthermore, for a straight channel the electro-osmotic drag and the specific pore conductivity match results
<figure><img src="image_3.png" /><figcaption>**Fig. 2.** Top: $u_3$, velocity in $x_3$ direction, zoom towards the channel intersection, plotted on the cross section $L \cap \{x_1 = 0\}$ (s.t. the channels are sliced in half). Bottom: Cross section $S \cap \{x_1 = 0\}$. $U_2$, mechanical displacement field in $x_2$ direction, deformed configuration. A channel may close completely. In both figures the cylinders have diameter $d = 2$ nm in the reference configuration, length 10 nm and offset $s = 0.5$ nm; $g_0 = 0$.</figcaption></figure>

from experiments, see the discussion in [12], supporting the validity of our model. The effect (b) raises the question whether it is possible to find an equilibrium shape for the nanochannel pore, balancing the solid pressure with the liquid pressure and the *interfacial* pressure, due to surface tension and mean curvature of the interface. This question yields a mathematically challenging problem and requires refined numerical techniques. For suitable outer pressures and external electric field, the above described algorithm suggests fast convergence, but a mathematical justification thereof is missing. However, for a simplified version of our model we can prove existence and uniqueness of a solution and the shape differentiability of the solution. Using a variational energy formulation, we show that this numerically determined optimal shape minimizes the free energy. The latter results are the subject of an upcoming paper [13].
In summary, our full model, being derived from thermodynamical first prin-
ciples, presents a generalisation of previous models [4–6] by including surface
tension and bulk stresses in the solid Nafion, or compared to [10, 8, 9], by includ-
ing charged fluid flow. We emphasize that, contrary to our last contributions [4,
6], we have solved additionally, in the electrohydrodynamical part of our model,
the full Nernst-Planck-Poisson-Stokes system without neglecting higher order
terms. In particular, our novel approach, i.e. a continuum model of the electro-
chemical fluid flow within a nanochannel of a PEM including surface tension and
bulk stresses, allows to explain ohmic resistance and electro-osmotic drag. Our
mathematical results might turn out to be important for further advances in the
design of hydrogen fuel cells.

**Acknowledgements.** The authors would like to thank Toyota Motor Engi-
neering and Manufacturing North America (TEMA) for financial support of this
research. S.-J. K. thanks the University of Ottawa for its hospitality.

**References**

1. Promislow, K., Wetton, B.: PEM Fuel cells: A mathematical Overview. SIAM J. Appl. Math. 70, 369–409 (2009)

2. Schmidt-Rohr, K., Chen, Q.: Parallel cylindrical water nanochannels in Nafion fuel-cell membranes. Nat. Mater. 7, 75–83 (2008)

3. Karniadakis, G., Beskok, A., Aluru, N.: Microflows and nanoflows. Fundamentals and simulation. Springer, New York (2005)

4. Ladipo K., Berg, P., Kimmerle, S.-J., Novruzi, A.: Effects of radially depen-
dent parameters on proton transport in polymer electrolyte membrane nanopores.
J. Chem. Phys. 134, 074103-1-074103-12 (2011)

5. Castellanos, A.: Electrohydrodynamics. Springer, Wien/New York (1998)

6. Berg, P., Ladipo, K.: Exact solution of an electro-osmotic flow problem in a cylin-
drical channel of polymer electrolyte membranes. Proc. R. Soc. A 465, 2663–2679
(2009)

7. Schmuck, M.: Analysis of the Navier-Stokes-Nernst-Planck-Poisson system. M3AS 9, 993–1014 (2009)

8. Fortin, M.: Problèmes de surfaces libres en mécanique des fluides. In: Shape op-
timization and free boundaries, pp. 143–172. Kluwer Academic Publishers, Dor-
drecht/Boston/London (1992)

9. Discacciati, M., Fourestey, G., Quarteroni, A., Deparis, S.: Fluid-structure algo-
rithms based on Steklov-Poincaré operators. Comput. Methods Appl. Mech. Engrg.
195, 5797–5812 (2006)

10. Elfring, G.J., Struchtrup, H.: Thermodynamics of pore wetting and swelling in Nafion, J. Membr. Sci. 315, 125–132 (2008)

11. Müller, I., Strehlow, P.: Rubber and rubber balloons. Paradigms of thermodynam-
ics. Springer, Berlin/Heidelberg (2004)

12. Kimmerle, S.-J., Novruzi, A., Berg, P., Ladipo, K.: Ohmic resistance of charged fluid flow in deformable nanochannels connecting polymer electrolyte membranes. Preprint submitted

13. Berg, P., Kimmerle, S.-J., Novruzi, A.: Modeling, shape analysis and computation of equilibrium pore shape near a PEM-PEM intersection. Preprint