Introduction

Membranes play a prominent role in many chemical, biological, and physical processes, where they separate two environments with different chemical composition1,2. Their microscopic structure and features determine the type of chemical compounds that cannot diffuse across them. Particularly important is the capability of charged membranes to interact with charged species3,4. Cell membranes, which bear a negative fixed charge5, are at the basis of cell physiology6, and they are often considered as one of the most important component of cellular life7. Similarly, artificial ion-exchange membranes, wherein ionic groups are covalently bonded to the polymeric skeleton, are widely employed in fuel cells, hydrogen production, and energy storage technologies4,8.

Modeling of charged membranes is central to understand biological processes and designing new, sustainable energy devices. Typically, charged membranes have been studied through the lens of electrochemistry, wherein the motion of charges is modeled by means of homogenized continuum treatments9,10,11, with roots in the classical Poisson-Nernst-Planck system that describes electrolytes in solution12. Extensions of these models include complex, non-ideal behavior of electrolytes within charged membranes13,14.

In this picture, membrane mechanics has often been neglected and relegated to a handful of phenomena in which deformations play a major role, such as mechanosensitive ionic gates in cell membranes15. A more thorough analysis of membrane chemoelectromechanics has been proposed for actuators and sensors based on ionic polymers, which transduce mechanical to electrochemical energy and vice versa9,11,16,17,18,19. Despite these inquiries, the coupling between mechanics and electrochemistry has yet to be explored in many of the natural and artificial processes in which charged membranes play a key role.

Recent theoretical and experimental endeavors highlighted profound, unexpected interactions between the mechanics and electrochemistry of charged membranes, with broad implications on biological membranes and energy storage and production devices. On the one hand, the non-ideal behavior of the solution permeating charged membranes can significantly affect the membrane deformations20. Specifically, solvation, the interaction between solvent molecules and counterions21,22, can elicit macroscopic bending of charged membranes in electrolyte solutions. On the other hand, the mechanical deformations of charged membranes can have sizeable effects on membrane electrochemistry. For example, the migration of solvent within charged membranes, which normally occurs from high to low counterions’ concentrations, can be reversed by volumetric deformations of the membrane, contradicting our intuition based on osmosis23.

Understanding the role of solution non-idealities on the coupled mechanics and electrochemistry of charged membranes constitutes an untapped problem. We know that solution non-idealities are expected to bear important consequences on the electrochemistry of the solution permeating the charged membrane. In fact, several studies have previously highlighted the critical role of non-idealities in solutions. For example, solvation is responsible for the lower conductivity of small ions compared to ions of intermediate size, contrary to our intuition based on the radii of isolated ions, the so-called Stokes-Nernst-Einstein’s predictions24. Similarly, long-range electrostatic interactions between ions can explain the remarkable discrepancies from ideality in the evaluation of the activity of the ions in an electrolyte solution or a charged membrane, within the classical Debye-Hückel theory25. Other types of non-idealities can also contribute to experimentally measurable changes in the membrane electrochemistry13. Quantifying how these non-idealities work together or at odds to modulate the coupled mechanics and electrochemistry of charged membranes is the objective of our effort.

Clearly, there is a dire need for mathematical models that can describe the chemoelectromechanics of charged membranes, accounting for large deformations and non-ideal behavior of solutions in the membrane. Here, we bridge this gap in the literature by proposing a continuum-level, faithful theory of the coupled mechanics and electrochemistry of charged membranes. The theory includes the incompressibility constraint accounting for the finite volume of counterions, which is responsible for the inversion of solvent migration23, along with four different types of non-ideal solution behavior that are expected to play a substantial role in membrane response.

Following the nomenclature in refs. 13,14, the four modeled sources of non-idealities are (Fig. 1): solvation, long-range electrostatic interactions, short-range physical interactions, and steric effects due to finite pore size. Solvation indicates the interaction between solvent and solute molecules, which affects ideal mixing by causing the formation of solvation shells of solvent around counterions21,22. Long-range electrostatic interactions are well-known in the electrochemistry literature, since the seminal work of Debye and Hückel25 aimed at explaining departures from the ideal behavior of electrolyte solutions. These departures are associated with the rearrangement of ions in solution, whereby cations tend to be surrounded by anions and vice versa26,27. Short-range physical interactions are related to the interactions between pairs of ions in the non-diluted limit26,28. Such interactions occur between ions of opposite sign and depend on the specific characteristics of each ion in the couple28. Finally, steric effects are related to the finite volume of pores in the membrane13. Due to these effects, the concentrations of solvent and counterions cannot grow indefinitely, unless the volume of pores is increased due to volumetric deformation.

Fig. 1: Schematics of the four types of non-idealities considered in this manuscript.
figure 1

Here, counterions and fixed ions are represented as green and blue spheres, while solvent molecules are shown as red-gray molecules (focusing on water, oxygen, and hydrogen atoms are shown in red and gray, respectively). The insets show the four non-idealities: (1) solvation (the formation of solvation shells around ions); (2) long-range electrostatic interactions (the arrangement of ions such that charges of one sign tend to be surrounded by charges of the opposite sign); (3) short-range physical interactions (the interaction of closeby ions with opposite charges); and (4) steric effects due to finite pore size (the limitation to the maximum concentration of counterions and solvent due to the finite volume of the pores of the membrane). The orange lines highlight the interactions associated with each type of non-ideality.

Full size image

We highlight that the ‘ideal’ model we define within our modeling framework is one in which the free energy of mixing only contains the ideal mixing term. In classical models of electrolyte solutions, the definition of the ‘ideal’ model is different, whereby long-range electrostatic interactions are typically included in addition to the ideal mixing term. Our choice of not including this term in the ‘ideal’ model is motivated by the previous literature on the mechanics and electrochemistry of charged membranes, where long-range electrostatic interactions have typically been neglected9,10,11.

As a key application of our theory, we study the steady-state actuation of an ionic polymer actuator. Through a series of simplifying hypotheses, we reduce the three-dimensional chemoelectromechanical problem to a nonlinear system of integro-differential equations that can be solved numerically. We study membranes with small (lithium) and large (cesium) monovalent counterions, which are characterized by two different sets of parameters that have a dramatic effect on actuation. Our results unravel a significant interplay between non-ideal solution behavior and mechanics, which are of major interest for the study of charged membranes.

Results

Continuum theory of mechanics and electrochemistry of charged membranes

Without loss of generality, we consider a charged membrane with monovalent negative fixed ions and positive counterions, mirroring cell membranes and common membranes utilized in fuel cells and ionic polymer transducers29. We assume perfect permselectivity of the membranes, such that no ion with a charge of the same sign of fixed ions can cross the membrane. We assume that the membrane undergoes small deformations with respect to a reference configuration that is stress-free, electroneutral, and with zero electric field. This hypothesis allows us to confound the reference and deformed configurations, and therefore material and spatial variables. Thus, variables are functions of the position x in the reference configuration and time t. The extension to finite deformations is tackled in the Supplementary Information.

The mechanics of the membrane is defined by the displacement u(x, t) of the material particle x at time t. Membrane deformation is described by the infinitesimal strain tensor \({{{\boldsymbol{\varepsilon }}}}({{{\bf{x}}}},t)=\frac{1}{2}[{{{\boldsymbol{\nabla }}}}{{{\bf{u}}}}({{{\bf{x}}}},t)+{\left({{{\boldsymbol{\nabla }}}}{{{\bf{u}}}}({{{\bf{x}}}},t)\right)}^{{{{\rm{T}}}}}]\), where ∇(⋅) is the nabla operator. The trace of the small strain tensor \({{{\rm{tr}}}}\,({{{\boldsymbol{\varepsilon }}}})\) determines the difference in volume compared to the reference configuration. The electrochemistry is described by three scalar fields: the concentration of solvent C0(x, t), the concentration of counterions C+(x, t), and the electric potential φ(x, t) with respect to a common ground. We assume that the concentration of fixed anions C−(x, t) is uniform, \({C}_{-}({{{\bf{x}}}},t)={\bar{C}}_{-}\), such that \({C}_{+}={\bar{C}}_{-}\). In the reference configuration, we consider that the concentration of solvent is also uniform, \({C}_{0}={\bar{C}}_{0}\).

The deformation of the membrane is governed by equilibrium, which, upon neglecting inertia and body forces, requires that the stress tensor is divergence-free9, see “Methods”. Migration of solvent and counterions is governed by mass conservation. We assume that no reaction occurs within the membrane, such that concentrations’ variations over time are associated to fluxes of solvent or counterions, see “Methods”. We neglect electrodynamics phenomena, such that Maxwell equations reduce to the Gauss law30 that relates the electric displacement D to the net charge density, see “Methods”.

Following the examples of the literature11,23,31,32,33, we put forward two constraints that allow us to significantly simplify the description of the motion of the solution in the membrane. Specifically, we assume that connected pores are fully saturated, and we impose that the individual constituents of the mixture (solid phase, counterions, and solvent) are incompressible, such that23

$${{{{\mathcal{V}}}}}_{0}({C}_{0}-{\bar{C}}_{0})+{{{{\mathcal{V}}}}}_{+}({C}_{+}-{\bar{C}}_{+})={{{\rm{tr}}}}\,({{{\boldsymbol{\varepsilon }}}}),$$
(1)

where \({{{{\mathcal{V}}}}}_{(\cdot )}\) is the molar volume of the corresponding species. The inclusion of the finite volume of ions in Eq. (1) is critical to capture the coupling between mechanical deformations and electrochemistry in charged membranes, including the inversion of solvent migration23.

Constitutive equations and excess energy

We employ a thermodynamically consistent approach to determine the constitutive equations for the charged membrane9,32,34. To this end, we define a Helmholtz free-energy density Ψ per unit volume of the mixture. We hypothesize that this quantity can be additively decomposed in four contributions9,13: the mechanical strain energy Ψmec(ε), the ideal energy of mixing Ψmix(C0, C+), the excess energy that encapsulates solution non-idealities within the membrane Ψexc(C0, C+), and the polarization energy \({{{\Psi }}}_{{{{\rm{pol}}}}}({{{\boldsymbol{\varepsilon }}}},{{{\bf{D}}}})\). We modify the free-energy density to impose the constraint in Eq. (1), obtaining the following functional:

$$\begin{array}{l}{{{\Psi }}}_{{{{\rm{f}}}}}({{{\mathbf{\varepsilon }}}},{C}_{0},{C}_{+},{{{\bf{D}}}},\pi )={{{\Psi }}}_{{{{\rm{mec}}}}}({{{\mathbf{\varepsilon }}}})+{{{\Psi }}}_{{{{\rm{mix}}}}}({C}_{0},{C}_{+})+{{{\Psi }}}_{{{{\rm{exc}}}}}({C}_{0},{C}_{+})\\\qquad\qquad\qquad\qquad\quad\;+\,{{{\Psi }}}_{{{{\rm{pol}}}}}({{{\mathbf{\varepsilon }}}},{{{\bf{D}}}})-\pi \left[\right.{{{\rm{tr}}}}\,({{{\mathbf{\varepsilon }}}})-{{{{\mathcal{V}}}}}_{0}\left({C}_{0}-{\bar{C}}_{0}\right)\\ \qquad\qquad\qquad\qquad\quad\;-\,{{{{\mathcal{V}}}}}_{+}({C}_{+}-{\bar{C}}_{+})\left.\right].\end{array}$$
(2)

Here, the hydraulic pressure π serves as a Lagrange multiplier to impose the constraint in Eq. (1), whereby taking the variation of Eq. (2) with respect to π and setting it to zero provides the incompressibility constraint. Constitutive equations that satisfy the second principle of thermodynamics can be derived upon differentiation of the Helmholtz free-energy density34, see “Methods”. The only quantities that are not defined by the functional are the solvent and counterions’ fluxes, for which we consider a Nernst-Planck form11, see “Methods”. Due to the presence of ε, C0, and C+ in the incompressibility constraint, the modification of the free-energy density with π affects the stress and electrochemical potentials of solvent and cations.

The mechanical, mixing, and polarization free-energy contributions are relatively standard in the study of charged membranes, and are defined in “Methods”. On the other hand, the definition of excess energy associated with non-idealities in the solution constitutes a fundamental advance in the study of coupled mechanics and electrochemistry of charged membranes. We decompose the excess energy into four contributions13, associated with solvation, long-range electrostatic interactions, short-range physical interactions, and steric effects due to finite pore size.

The term solvation identifies the interaction between solvent molecules and ions in electrolyte solutions21,22. Depending on the nature of the solvent, interactions may include dipole-dipole interactions, ion-induced dipole, and hydrogen bonds. Thus, solvent molecules dispose themselves around ions, creating so-called solvation shells21,22. The rearrangement of solvent molecules and ions due to solvation generates an excess energy compared to the ideal mixing case. In this case, both mobile counterions and fixed anions contribute to the energy. We follow the example of the literature13,28,35 toward defining the free energy of solvation. The model relies on the hypothesis of a stepwise solvation process, from which one can write the free energy as

$$\begin{array}{l}{{{\Psi }}}_{{{{\rm{slv}}}}}({C}_{0},{C}_{+})={{{\mathcal{R}}}}{{{\mathcal{T}}}}\left[{C}_{0}\log {\alpha }_{0}({C}_{0},{C}_{+})+{C}_{+}\log \left(\frac{{\alpha }_{+}({C}_{0},{C}_{+})}{{\alpha }_{+}^{\infty }}\right)\right.\\ \qquad\qquad\qquad\;\;\left.+\,{C}_{-}\log \left(\frac{{\alpha }_{-}({C}_{0},{C}_{+})}{{\alpha }_{-}^{\infty }}\right)\right].\end{array}$$
(3)

Here, \({{{\mathcal{R}}}}\) is the universal gas constant, \({{{\mathcal{T}}}}\) is the temperature, α0 is the ratio between the free solvent molar fraction in the solvated and unsolvated states, αi is the ratio between the molar fraction of unbound ions i = +, − in the solvated and unsolvated states, and \({\alpha }_{i}^{\infty }\) is the fraction of unbound ions i = +, − in the infinite dilute state.

To determine α0, α+, α−, \({\alpha }_{+}^{\infty }\), and \({\alpha }_{-}^{\infty }\), we ought to define the free molar fractions of solvent and ions,

$${\chi }_{0}=\frac{{C}_{0}}{{C}_{0}+{C}_{+}+{C}_{-}},$$
(4)
$${\chi }_{i}=\frac{{C}_{i}}{{C}_{0}+{C}_{+}+{C}_{-}},$$
(5)

for i = +, −. We also have to define the free molar fractions of solvent and ions in the solvated state. For the solvent, we need to reduce the concentration by the fraction of solvent molecules bonded to cations and to anions, \({C}_{+}{\bar{h}}_{+}\) and \({C}_{-}{\bar{h}}_{-}\), respectively, where \({\bar{h}}_{i}\) is the average solvation number of ion species i, which quantifies the average number of solvent molecules bonded to each ion. Thus, we obtain

$${\chi }_{0}^{{{{\rm{slv}}}}}=\frac{{C}_{0}-{C}_{+}{\bar{h}}_{+}-{C}_{-}{\bar{h}}_{-}}{{C}_{0}+{C}_{+}(1-{\bar{h}}_{+})+{C}_{-}(1-{\bar{h}}_{-})}.$$
(6)

For the ions, the free molar fractions in the solvated state are defined as

$${\chi }_{i}^{{{{\rm{slv}}}}}=\frac{\frac{{C}_{i}}{{(1+{k}_{i}{\alpha }_{0}{\chi }_{0})}^{{{{{\mathcal{N}}}}}_{i}}}}{{C}_{0}+{C}_{+}(1-{\bar{h}}_{+})+{C}_{-}(1-{\bar{h}}_{-})},$$
(7)

where ki is the binding constant of ion species i and \({{{{\mathcal{N}}}}}_{i}\) is the number of binding sites for the solvent on ion species i. The binding constant is related to the enthalpy of solvation \({{\Delta }}{H}_{i}^{{{{\rm{slv}}}}}\) through \({k}_{i}=\exp \left[-{{\Delta }}{H}_{i}^{{{{\rm{slv}}}}}/({{{\mathcal{R}}}}{{{\mathcal{T}}}})\right]\). Based on these definitions, we can finally write

$${\alpha }_{0}=\frac{{\chi }_{0}^{{{{\rm{slv}}}}}}{{\chi }_{0}}=\frac{\frac{{C}_{0}-{C}_{+}{\bar{h}}_{+}-{C}_{-}{\bar{h}}_{-}}{{C}_{0}+{C}_{+}(1-{\bar{h}}_{+})+{C}_{-}(1-{\bar{h}}_{-})}}{\frac{{C}_{0}}{{C}_{0}+{C}_{+}+{C}_{-}}},$$
(8)
$${\alpha }_{i}=\frac{{\chi }_{i}^{{{{\rm{slv}}}}}}{{\chi }_{i}}=\frac{1}{{\chi }_{0}{(1+{k}_{i}{\alpha }_{0}{\chi }_{0})}^{{{{{\mathcal{N}}}}}_{i}}\left(1+\frac{{C}_{+}}{{C}_{0}}(1-{\bar{h}}_{+})+\frac{{C}_{-}}{{C}_{0}}(1-{\bar{h}}_{-})\right)}.$$
(9)

In the infinitely dilute state, we have C+ ≪ C0 and C− ≪ C0, such that χ0 → 1 and α0 → 1. As such, from Eq. (9), we find that the fraction of unbound ions in the infinite dilute state is

$${\alpha }_{i}^{\infty }=\frac{1}{{(1+{k}_{i})}^{{{{{\mathcal{N}}}}}_{i}}}.$$
(10)

The average solvation number of ion species i, \({\bar{h}}_{i}\), can be expressed as \({\bar{h}}_{i}={{{{\mathcal{N}}}}}_{i}{k}_{i}{\alpha }_{0}{\chi }_{0}/(1+{k}_{i}{\alpha }_{0}{\chi }_{0})\). Since \({\bar{h}}_{i}\) depends on α0, the contribution of solvation to the excess energy is defined implicitly. By substituting the expression of \({\bar{h}}_{i}\) into Eq. (8), we obtain a nonlinear equation that can be solved for α0 for given values of C0 and C+.

Similar to solvation, long-range electrostatic interactions play a significant role in determining the activity coefficients of electrolytes in solution26. Due to these interactions, ions of one sign tend to be surrounded by ions of opposite sign, contrary to what one would expect from ideal mixing. To model this behavior, we consider a modified version of Debye-Hückel theory27 that accounts for the presence of fixed anions in the membrane13. We write

$$\begin{array}{ll}{{{\Psi }}}_{{{{\rm{els}}}}}({C}_{0},{C}_{+})=-\frac{4}{3}{{{\mathcal{R}}}}{{{\mathcal{T}}}}A{I}^{\frac{3}{2}}({C}_{0},{C}_{+}){{{{\mathcal{V}}}}}_{0}{C}_{0}\tau [aB{I}^{\frac{1}{2}}({C}_{0},{C}_{+})]\\ \qquad\qquad\qquad\;\;\,-\,{{{\mathcal{R}}}}{{{\mathcal{T}}}}{C}_{-}\displaystyle\frac{2A}{bB}\log (bB{I}^{\frac{1}{2}}),\end{array}$$
(11)

where A and B are the Debye-Hückel limiting slope and solvent parameter (known for each solvent), respectively, a is the distance of closest approach of counterions (that is, their average diameter), b is the spacing between fixed charges in the polymer, \(\tau (x)=3\left[\log (1+x)-x+{x}^{2}/2\right]/{x}^{3}\), and I is the ionic strength36, computed only for mobile counterions as \(I({C}_{0},{C}_{+})={C}_{+}/(2{C}_{0}{{{{\mathcal{V}}}}}_{0})\).

The third term in the excess energy represents the free energy associated with short-range physical interactions between cations and anions in solution13,28. Physical interactions between ions of the concordant sign are instead much more sporadic and are here neglected13. This term is of semi-empirical origin, based on simple physical considerations on the probability of encounter of cations and anions in solution (similar to kinetic theory) and fitting of solution data. Following the example of ref. 13, specialized for the case of a single species of positive counterions and fixed negative ions, we define

$${{{\Psi }}}_{{{{\rm{phy}}}}}({C}_{0},{C}_{+})=2\frac{{{{\mathcal{R}}}}{{{\mathcal{T}}}}}{{C}_{0}}{\beta }_{+-}{C}_{+}{C}_{-},$$
(12)

where β+− is the interaction coefficient between cations and anions.

Finally, we consider the free energy associated with steric effects due to finite pore size13. Due to the finite size of the pores, we need to consider the excluded volume of counterions in the saturating solution, which would affect the mixing of free energy. Following the approach in ref. 13, which considers the pores as the space between randomly oriented walls37, we define

$${{{\Psi }}}_{{{{\rm{stc}}}}}({C}_{0},{C}_{+})=\frac{{{{\mathcal{R}}}}{{{\mathcal{T}}}}}{d({C}_{0},{C}_{+})-{d}_{0}}{C}_{+}{a}_{+},$$
(13)

where d(C0, C+) is the actual spacing between hydrophilic domains, and d0 is the spacing between hydrophilic domains in the dry membrane. These two quantities are related experimentally by the volumetric deformation of the membrane13, which we express as a function of C0 and C+ through the incompressibility constraint in Eq. (1).

Static solution for ionic polymer actuators

As a key application of our theory to describe the complex interplay between mechanics and electrochemistry in charged membranes, we focus on ionic polymer actuators. We analyze the steady-state response of an ionic polymer actuator under the premises of small deformations and plane-strain, representing the master problem to be studied to describe the multiaxial response of these actuators38,39. Under these two assumptions (small-strain and plane-strain), the two-dimensional problem can be effectively decomposed into a (highly nonlinear and computationally challenging) one-dimensional problem through the thickness of the membrane and a one-dimensional problem along the width (with a trivial solution). Such a decomposition has been thoroughly validated through nonlinear, two-dimensional finite element simulations38,39. Ultimately, we obtain a ‘semi-analytical’ solution that allows us to solve the two-dimensional mechanics and electrochemistry through the numerical solution of a one-dimensional electrochemical problem. We distinguish this type of solution from the ‘numerical’ solution of the two-dimensional mechanics and electrochemistry.

We consider an ionic polymer actuator of length L and thickness 2H ≪ L, with zero-thickness electrodes (Fig. 2). We define a reference frame centered on the left end of the membrane mid-axis, with the x-axis along the actuator length, y-axis along the membrane thickness, and z-axis along the membrane width to form a right-handed coordinate frame. We assume that all boundaries are impermeable to both solvent and counterions, such that the total mass of solution in the membrane is conserved. For classical ionic polymer actuators, while displacements may be very large, the strains typically remain small (below 5%40), due to the small thickness of the actuators. We focus on plane-strain deformations, which are non-zero only within the x-y plane. The plane-strain hypothesis is verified for standard ionic polymer actuators, in which the width is several times their thickness. We study the static actuation of the ionic polymer for a voltage \(\bar{V}\) applied across its electrodes at y = H (anode) and y = −H (cathode). We consider simply supported boundary conditions for the actuator, such that the left end of the membrane mid-axis is constrained not to translate in the x-y plane, whereas the right end of the membrane mid-axis can only translate along the x-axis, see Fig. 2.

Fig. 2: Schematics of the plane-strain actuation problem for an ionic polymer with corresponding boundary conditions and reference frame.
figure 2

Here, positive mobile counterions are depicted as green spheres, fixed negative charges as blue cubes, and solvent molecules are shown as red-gray molecules (focusing on water, oxygen and hydrogen atoms are shown in red and gray, respectively).

Full size image

As we focus our attention on static actuation, we eliminate the dependence on time. Similar to ref. 38, we hypothesize that the variation of electrochemical quantities along the length and width of the membrane is negligible compared to the variation along the thickness. As a consequence, all electrochemical variables depend only on the y coordinate. We expect this hypothesis to be valid along the entire length of the ionic polymer actuator, apart from small regions in proximity of the supports, as extensively verified through finite element simulations in previous studies38,39. In the static regime, fluxes of solvent and counterions must be equal to zero, such that their (electro)chemical potentials must be constant. Analogously, Gauss law reduces to a second-order ordinary differential equation along the thickness of the membrane.

We follow the approach in ref. 38 to obtain a Saint-Venant solution for the membrane deformation, which can accurately describe multiaxial deformations of simply supported ionic polymer actuators. In this vein, we hypothesize that shear strains are negligible, such that εxx and εyy are the only non-zero strain components. We assume a linear form of the longitudinal strain εxx, an assumption that has been widely verified through finite element simulations38,39. As equilibrium requires the through-the-thickness stress to be zero, we can express the through-the-thickness strain εyy in terms of the longitudinal strain εxx and electrochemical quantities, see “Methods”. Since the actuator is unloaded, the longitudinal strain is written as a function of the electrochemical quantities by requiring the longitudinal force and bending moment per unit width of the membrane to be zero. Thus, the incompressibility constraint in Eq. (1) can be written in terms of electrochemical variables only.

The resulting system of equations, including the constant solvent and counterions’ electrochemical potentials, through-the-thickness Gauss law, and the incompressibility constraint, constitutes an integro-differential system in the variables C0(y), C+(y), π(y), and φ(y). We complement this system with two boundary conditions for φ, \(\varphi (\pm H)=\pm \bar{V}/2\). The electrochemical potentials of solvent and counterions, defined up to an additive constant, are obtained by imposing the mass conservation of solvent and counterions, which is not automatically satisfied in the static case, as well-known from the classical Poisson-Boltzmann equations12.

We solve this system by adapting the numerical method for two-point boundary value problems in ref. 41, see “Methods”. We solve the resulting nonlinear algebraic system numerically. We utilize parameters for ionic polymer actuators based on a typical commercial Nafion™ membrane with two different types of counterions, lithium and cesium, and water as the solvent, see “Methods”. We apply a voltage of 1 Vth, where \({V}_{{{{\rm{th}}}}}={{{\mathcal{R}}}}{{{\mathcal{T}}}}/{{{\mathcal{F}}}}\) is the thermal voltage (~25 mV at room temperature), being \({{{\mathcal{F}}}}\) the Faraday constant.

Electrochemistry of ionic polymer actuators

In Figs. 3 and 4, we show the steady-state profiles of the electrochemical quantities through the thickness of membranes with lithium and cesium counterions, respectively. Specifically, we consider the concentration of solvent C0, the concentration of cations C+, the pressure π that imposes incompressibility and generates a stress in the membrane, and the voltage ϕ within the membrane. These variables are scaled by the reference concentration of solvent \({\bar{C}}_{0}\), the reference concentration of cations \({\bar{C}}_{+}\), a reference pressure \({{{\mathcal{R}}}}{{{\mathcal{T}}}}{\bar{C}}_{+}\), and the thermal voltage \({{{\mathcal{R}}}}{{{\mathcal{T}}}}/{{{\mathcal{F}}}}\), respectively.

Fig. 3: Electrochemistry of ionic polymer actuators with lithium counterions.
figure 3

Steady-state profiles of solvent concentration C0 (a, b), cations' concentration C+ (c, d), pressure π (e, f), and voltage φ (g, h) for an ionic polymer actuator with lithium counterions, close to the cathode (a, c, e, g), and close to the anode (b, d, f, h). The blue solid line indicates the results with solution non-idealities, whereas the red dashed line represents the results for an ideal solution.

Full size image

Similar to the ideal case, in the bulk of the membrane, concentrations of solvent and counterions remain the same as in the reference configuration, for both lithium and cesium counterions. Boundary layers of solvent (Figs. 3a, b, 4a, b) and counterions (Figs. 3c, d, 4c, d) are formed close to the cathode and anode. Non-idealities widen the boundary layers and considerably decrease the pile-up and depletion of counterions in the boundary layers compared to the ideal case. For example, for the pile-up at the cathode, we register −54% and −66% of the maximum concentration for lithium and cesium counterions, respectively.

With lithium counterions, solvent migrates in the same direction of counterions, as one would expect from osmosis, regardless of the presence of non-idealities (Fig. 3a, b). Contrary to counterions, non-idealities increase the pile-up of solvent at the cathode (33%) and the depletion at the anode (37%), by increasing both the change of solvent concentration with respect to the reference value and the thickness of boundary layers. With cesium counterions, for the ideal solution, we observe the solvent migration inversion predicted in ref. 23, opposing osmosis (Fig. 4a, b). We find that non-idealities completely hinder this phenomenon, both at the cathode and anode (Fig. 4a, b), such that we observe solvent migration in the direction expected from osmosis. The suppression of inversion of solvent migration is an important qualitative difference introduced by non-ideal solution behavior.

Fig. 4: Electrochemistry of ionic polymer actuators with cesium counterions.
figure 4

Steady-state profiles of solvent concentration C0 (a, b), cations' concentration C+ (c, d), pressure π (e, f), and voltage φ (g, h) for an ionic polymer actuator with cesium counterions, close to the cathode (a, c, e, g), and close to the anode (b, d, f, h). The blue solid line indicates the results with solution non-idealities, whereas the red dashed line represents the results for an ideal solution.

Full size image

For both counterions, the pressure in the bulk is zero, similar to the reference configuration. The pressure is higher than this reference value at the cathode and lower at the anode, regardless of the presence of non-idealities (Figs. 3e, f, 4e, f). For lithium membranes, we observe an increase in the difference of pressure with respect to the reference configuration when considering non-idealities (21% at the cathode and 40% at the anode), which also widen the boundary layers (Fig. 3e, f). When considering cesium membranes with non-idealities, we register a smaller peak pressure difference in the boundary layers compared to the ideal case (Fig. 4e, f), with a 39% reduction at the cathode and 33% at the anode. Such a decrease is accompanied by a smaller rate of change in the profile and thicker boundary layers compared to ideality.

The electric potential is the only variable that is not zero in the membrane bulk (Figs. 3g, h, 4g, h), where it achieves a negative value that is slightly reduced by non-idealities. Non-idealities cause a larger voltage drop at the cathode, while decreasing the one at the anode (Figs. 3g–h, 4g–h). Further, non-idealities widen the boundary layers at both electrodes.

Free-energy contributions in ionic polymer actuators

We now focus on the differences in free-energy contributions with respect to the reference configuration \({{{\Psi }}}_{(\cdot )}-{\bar{{{\Psi }}}}_{(\cdot )}\) (Fig. 5). For both counterions, all the changes are concentrated in the boundary layers.

Fig. 5: Free-energy changes in ionic polymer actuators.
figure 5

Steady-state profiles of the differences in free-energy contributions with respect to the reference configuraton for an ionic polymer actuator with lithium (a, b) and cesium (c, d) counterions, close to the cathode (a, c), and close to the anode (b, d). The black solid line is the mixing free energy, the red dash-dotted line is the solvation free energy, the green loosely dotted line is the free energy associated with electrostatic interactions, the blue finely dotted line is the free energy related to physical interactions, and the brown dashed line is the free-energy associated with steric effects due to finite pore size.

Full size image

The contribution of the mixing free energy is negative, and is the dominant one in absolute value for both counterions at both electrodes. For cesium counterions, its value is almost half that for lithium ones. The second most important contribution is the one associated with electrostatic interactions, which is negative and reaches more than 50% of the mixing free-energy values. This contribution is slightly smaller in cesium membranes compared to lithium ones. The solvation term is positive and of similar magnitude for both counterions. The one associated with physical interactions is almost identical to solvation in the case of lithium membranes, whereas it has the opposite sign for cesium membranes. Finally, the contribution related to steric effects is negligible at this voltage level for both counterions.

Mechanics of ionic polymer actuators

Finally, we look at the effect of non-idealities on the mechanics of ionic polymer actuators. First, we investigate the longitudinal and through-the-thickness strains (Fig. 6). As expected from our hypotheses, the longitudinal strain is linear throughout the membrane thickness and does not present boundary layers in the proximity of the electrodes (Fig. 6a, d). The slope of the strain indicates the actuator curvature, while the value at y = 0 is the mid-axis strain38. The inverse of the curvature is the radius of curvature, which for our simulations is around 5 m, in line with what we would expect at these levels of applied voltage from experimental results42. While for cesium counterions the longitudinal strain profile is analogous to the ideal case (Fig. 6d), non-idealities for lithium counterions cause a larger steady-state curvature and smaller mid-axis strain compared to the ideal scenario (Fig. 6a). The effect of non-idealities on the final curvature is remarkable, whereby the value of the curvature almost doubles compared to the ideal case.

Fig. 6: Strains in ionic polymer actuators.
figure 6

Steady-state profiles of the longitudinal strain εxx (a, d), and through-the-thickness strain εyy close to the cathode (b, e) and close to the anode (c, f), for an ionic polymer actuator with lithium counterions (a–c) and cesium counterions (d–f). The blue solid line indicates the results with solution non-idealities, whereas the red dashed line represents the results for an ideal solution.

Full size image

The through-the-thickness strain is negligible in the bulk of the material and strain concentration occurs close to the electrodes, in correspondence to the electrochemical boundary layers (Fig. 6b, c, e, f), with positive strains at the cathode (Fig. 6b, e) and negative strains at the anode (Fig. 6c, f). Regardless of non-idealities and counterions’ type, through-the-thickness strain at the electrodes is three orders of magnitude larger than longitudinal strain. For lithium membranes, the presence of non-idealities causes an increase in the absolute values of the through-the-thickness strain and in the thickness of its boundary layers, further exacerbating strain concentration at the electrodes (Fig. 6b, c). In the case of cesium, non-idealities decrease the peak absolute values of the strain in the boundary layers, while increasing the thickness of the boundary layers and reducing the slope of their profile (Fig. 6e, f).

In Figs. 7 and 8, we display the steady-state profiles for lithium and cesium counterions, respectively, of the three contributions to the axial stress, which are related to the mechanical stress \({\sigma }_{xx}^{{{{\rm{mec}}}}}\), Maxwell stress \({\sigma }_{xx}^{{{{\rm{pol}}}}}\), and incompressibility \({\sigma }_{xx}^{{{{\rm{inc}}}}}\). Regardless of the counterions’ type, all the three contributions to the stress are zero in the bulk. Boundary layers of stress occur for all contributions, even to the mechanical one, due to Poisson effects associated with εyy. The mechanical stress contribution always shows positive values in the proximity of the cathode (Figs. 7a and 8a) and negative values near the anode (Figs. 7b and 8b). The effect of non-idealities on the stress contribution associated with mechanical stress is different depending on the type of counterion. Specifically, for lithium counterions, non-idealities cause an increase in the thickness of the boundary layer of mechanical stress and of its absolute value at both electrodes (27% at the cathode and 32% at the anode), compared to the ideal case (Fig. 7a, b). For cesium counterions, we still register a larger thickness of the boundary layers when considering non-idealities, but we find an overall smaller peak mechanical stress (Fig. 8a, b).

Fig. 7: Axial stresses in ionic polymer actuators with lithium counterions.
figure 7

Steady-state profiles of the three contributions to the axial stress, associated with mechanical stress (a, b), Maxwell stress (c, d), and incompressibility (e, f), for an ionic polymer actuator with lithium counterions, close to the cathode (a, c, e), and close to the anode (b, d, f). The blue solid line indicates the results with solution non-idealities, whereas the red dashed line represents the results for an ideal solution.

Full size image
Fig. 8: Axial stresses in ionic polymer actuators with cesium counterions.
figure 8

Steady-state profiles of the three contributions to the axial stress, associated with mechanical stress (a, b), Maxwell stress (c, d), and incompressibility (e, f), for an ionic polymer actuator with cesium counterions, close to the cathode (a, c, e), and close to the anode (b, d, f). The blue solid line indicates the results with solution non-idealities, whereas the red dashed line represents the results for an ideal solution.

Full size image

For both counterions’ type, non-idealities reduce the contribution of Maxwell stress (Figs. 7c, d and 8c, d). This contribution is always positive, and it is the smallest in absolute value among the three, likely due to the moderate levels of applied voltage. The two boundary layers are almost symmetric, since the nonlinearity is modest. As non-idealities elicit a widening of the boundary layers of voltage (Figs. 3g, h and 4g, h), the electric field is reduced compared to the ideal case, and so is Maxwell stress (Figs. 7c, d and 8c, d). The intensity of the reduction is higher for cesium cations, around 63% (Fig. 8c, d), compared to lithium cations, that experience a 50% reduction (Fig. 7c, d).

Finally, we consider the largest stress contribution, which is associated with incompressibility (Figs. 7e, f and 8e, f). This term is exactly the opposite of the hydraulic pressure shown in Figs. 3e, f and 4e, f. Thus, we always record a negative contribution at the cathode (Figs. 7e and 8e) and a positive contribution at the anode (Figs. 7f and 8f). Similar to the other quantities, non-idealities widen the boundary layers, regardless of the counterions’ type. Analogous to the pressure, for lithium counterions (Fig. 7e, f), non-ideal behavior increases in absolute value the stress contribution associated with incompressibility, compared to the ideal scenario. On the other hand, non-idealities in membranes with cesium counterions (Fig. 8e, f) cause a reduction of the peak stress contribution in absolute value compared to an ideal solution, along with a less sharp change of the stress profile from the electrode to the bulk.

Discussion

The coupling between mechanics and electrochemistry has so far been mostly neglected by the literature on charged membranes. Only recently, a few efforts highlighted a significant effect of non-ideal behavior of electrolyte solutions on membrane mechanics20, as well as counterintuitive consequences of including membrane deformations on the electrochemical response23. In this vein, it becomes critical to build a unifying, high-fidelity model of the interplay between mechanics and electrochemistry in charged membranes to inform their applications from cell physiology to energy storage and production devices.

In this manuscript, we established a continuum theory for charged membranes accounting for potentially large mechanical deformations and four non-idealities that can play a significant role in the mechanics and electrochemistry of charged membranes, namely: solvation, long-range electrostatic interactions, short-range physical interactions, and steric effects due to finite pore size. Toward obtaining insight into the importance of these four non-idealities, we focused on a specific application of our theory, analyzing the plane-strain, steady-state actuation of ionic polymer actuators. The resulting system of integro-differential equations was solved numerically to study the mechanics and electrochemistry of membranes with different counterions.

We found that non-idealities significantly affect the membrane electrochemistry, mediated by the size of the counterions. For small counterions, we observed an increase in the absolute value of the concentration of solvent, accompanied by a rise in the absolute value of the pressure. Larger counterions displayed an inversion in solvent migration in the ideal case, opposite to the one expected from osmosis23. Non-idealities contrasted the inversion of solvent migration, restoring migration in the direction predicted by osmosis. Hence, non-idealities increase the molar volume of counterions required for solvent migration inversion, compared to the values from the ideal case23. Such an increase could explain why inversion of solvent migration has not been observed yet in charged membranes, but only for large macromolecules in nanochannels43. In this vein, the model with non-idealities seems to describe better the physics of the inversion of solvent migration than the ideal one.

The breakdown of the various contributions to the membrane free-energy showed that non-idealities were not negligible compared to ideal mixing free energy. In particular, we discovered that long-range electrostatic interactions cannot be disregarded, especially for large counterions. Solvation and short-range interactions constitute other two important contributions, which should be accounted for in any high-fidelity charged membrane theory. Only steric effects due to finite pore size could be neglected, although we anticipate that they could play a more significant role at higher values of the applied voltage, similar to steric effects due to finite ion size44,45.

The effect of non-idealities on the electrochemistry was mirrored by membrane mechanics. Non-ideal behavior with small counterions increased the absolute values of membrane curvature and through-the-thickness strain, with no change in the curvature and a decrease in through-the-thickness strain for large counterions. Thus, non-idealities magnify the difference between actuators with small and large counterions, which have been experimentally observed in previous efforts46,47. Non-ideal behavior also affected the contributions to the axial stress, in a way that depends on the type of counterion.

Our continuum theory constitutes the first unifying framework for the study of the complex interplay between mechanics and electrochemistry in charged membranes, paving the way for a series of new efforts to quantify the effects of these couplings on charged membranes used in a variety of applications. Future studies should analyze dynamic phenomena in counterions’ and solvent migration, neglected in this fundamental work. Further, extensive simulations should be performed to unravel the nonlinear interplay between deformations, ideal mixing, polarization, and non-idealities at various applied voltages. These simulations should span the variability range of material parameters, to obtain the sensitivity of the results on each parameter. The high dimensionality of the resulting space suggests the adoption of dimensionality reduction techniques for the modeling of charged membranes48,49, toward a high-fidelity description of their mechanics and electrochemistry at a limited computational cost.

Methods

Governing equations

By neglecting inertia and body forces, mechanical equilibrium reads9

$${{{\rm{Div}}}}\,{{{\boldsymbol{\sigma }}}}={{{\bf{0}}}},$$
(14)

where σ is the Cauchy stress tensor34 and \({{{\rm{Div}}}}\,(\cdot )\) is the divergence operator.

Mass conservation of solvent and counterions are11

$$\frac{\partial {C}_{0}}{\partial t}+{{{\rm{Div}}}}\,{{{{\bf{J}}}}}_{0}=0,$$
(15)
$$\frac{\partial {C}_{+}}{\partial t}+{{{\rm{Div}}}}\,{{{{\bf{J}}}}}_{+}=0.$$
(16)

Here \(\frac{\partial }{\partial t}(\cdot )\) indicates the partial derivative with respect to time, while J0 and J+ are the molar fluxes of solvent and counterions per unit area.

The Gauss law for the charged membrane is30

$${{{\rm{Div}}}}\,{{{\bf{D}}}}=Q,$$
(17)

where Q is the net charge density,

$$Q={{{\mathcal{F}}}}({C}_{+}-{C}_{-}).$$
(18)

Free-energy contributions

Constitutive equations that satisfy the second principle of thermodynamics can be derived upon differentiation of the Helmholtz free-energy density34. For the strain energy, we consider a linear isotropic constitutive behavior, with9

$${{{\Psi }}}_{{{{\rm{mec}}}}}({{{\boldsymbol{\varepsilon }}}})=\frac{{\lambda }_{{{{\rm{L}}}}}}{2}{[{{{\rm{tr}}}}({{{\boldsymbol{\varepsilon }}}})]}^{2}+{\mu }_{{{{\rm{L}}}}}{{{\rm{tr}}}}\,({{{{\boldsymbol{\varepsilon }}}}}^{2}),$$
(19)

where λL and μL are the Lamé parameters50. The Lamé parameters are related to the Young modulus E and Poisson ratio ν through λL = Eν/[(1 + ν)(1 − 2ν)] and \({\mu }_{{{{\rm{L}}}}}=E/\left[2(1+\nu )\right]\)50.

Ideal mixing of free energy is associated with mixing entropy, describing the ideal mixing between distinguishable particles51. We write

$${{{\Psi }}}_{{{{\rm{mix}}}}}({C}_{0},{C}_{+})={{{\mathcal{R}}}}{{{\mathcal{T}}}}\left[{C}_{0}\log \left(\frac{{C}_{0}}{{C}_{0}+{C}_{+}}\right)+{C}_{+}\log \left(\frac{{C}_{+}}{{C}_{0}+{C}_{+}}\right)\right].$$
(20)

Contrary to ref. 11 and mirroring refs. 23,31, we do not assume a dilute solution in the ionic polymer, whereby this hypothesis may not be verified in boundary layers.

Finally, we assume that the ionic polymer behaves as a linear dielectric with dielectric constant ϵ. Thus, the polarization free-energy density is written as9

$${{{\Psi }}}_{{{{\rm{pol}}}}}({{{\boldsymbol{\varepsilon }}}},{{{\bf{D}}}})=\frac{1}{2\epsilon [1+{{{\rm{tr}}}}\,({{{\boldsymbol{\varepsilon }}}})]}\left[({{{\bf{I}}}}+2{{{\boldsymbol{\varepsilon }}}})\cdot (\tilde{{{{\bf{D}}}}}\otimes \tilde{{{{\bf{D}}}}})\right].$$
(21)

Here, ‘⋅’ is the inner product, whereas ‘⊗’ is the outer product.

Utilizing well-known procedures to identify relevant derivatives32,34, we find

$${{{\boldsymbol{\sigma }}}}=\frac{\partial {{\Psi }}}{\partial {{{\boldsymbol{\varepsilon }}}}}-\pi {{{\bf{I}}}},$$
(22)
$${\mu }_{0}=\frac{\partial {{\Psi }}}{\partial {C}_{0}}+\pi {{{{\mathcal{V}}}}}_{0},$$
(23)
$${\mu }_{+}=\frac{\partial {{\Psi }}}{\partial {C}_{+}}+{{{\mathcal{F}}}}\varphi +\pi {{{{\mathcal{V}}}}}_{+},$$
(24)
$${{{\bf{E}}}}=\frac{\partial {{\Psi }}}{\partial {{{\bf{D}}}}}.$$
(25)

Here E(x, t) = −∇φ(x, t) is the electric field, while μ0 and μ+ are the (electro)chemical potentials of solvent and counterions, respectively.

Solvent and counterions’ fluxes are not determined by the functional, as the second principle of thermodynamics only provides an inequality constraint −J0 ⋅ ∇μ0 − J+ ⋅ ∇μ+ ≥ 0. We assume a Nernst-Planck relationship between fluxes and electrochemical potentials11, \({\left[\begin{array}{cc}{{{{\bf{J}}}}}_{0}^{{{{\rm{T}}}}},&{{{{\bf{J}}}}}_{+}^{{{{\rm{T}}}}}\end{array}\right]}^{{{{\rm{T}}}}}=-{{{\bf{M}}}}{\left[\begin{array}{cc}{{{\boldsymbol{\nabla }}}}{\mu }_{0}^{{{{\rm{T}}}}},&{{{\boldsymbol{\nabla }}}}{\mu }_{+}^{{{{\rm{T}}}}}\end{array}\right]}^{{{{\rm{T}}}}}\), where M is a symmetric, positive semidefinite mobility tensor. Note that cross-diffusion terms are symmetrical, according to Onsager theory36. The choice of the mobility tensor can differ between different formulations. A potential choice is the one proposed in ref. 11,

$${{{{\bf{J}}}}}_{0}=-\frac{{{{{\mathcal{D}}}}}_{0}}{{{{\mathcal{R}}}}{{{\mathcal{T}}}}}\left({C}_{0}{{{\boldsymbol{\nabla }}}}{\mu }_{0}+{C}_{+}{{{\boldsymbol{\nabla }}}}{\mu }_{+}\right),$$
(26)
$${{{{\bf{J}}}}}_{+}=\frac{{C}_{+}}{{C}_{0}}{{{{\bf{J}}}}}_{0}-\frac{{{{{\mathcal{D}}}}}_{+}}{{{{\mathcal{R}}}}{{{\mathcal{T}}}}}{C}_{+}{{{\boldsymbol{\nabla }}}}{\mu }_{+},$$
(27)

where \({{{{\mathcal{D}}}}}_{0}\) and \({{{{\mathcal{D}}}}}_{+}\) are the diffusivities of solvent and counterions, respectively. Since we will only study the steady-state problem in which we neglect the variation of solvent and counterions’ concentration over time, this specific choice of fluxes does not affect our results.

Steric effects due to finite ion size

The actual spacing between hydrophilic domains d and the spacing between hydrophilic domains in the dry membrane d0 are related experimentally by ref. 13 as \(d({{{\boldsymbol{\varepsilon }}}})={d}_{0}{\phi }_{{{{\rm{solid}}}}}^{-m}({{{\boldsymbol{\varepsilon }}}})\), where m is a parameter that depends on the size of the hydrophilic domains and ϕsolid is the volumetric fraction of the solid phase in the deformed configuration. This quantity is related to the volumetric fraction of the solid phase in the undeformed configuration, Φsolid, through \({{{\Phi }}}_{{{{\rm{solid}}}}}=[1+{{{\rm{tr}}}}\,({{{\boldsymbol{\varepsilon }}}})]{\phi }_{{{{\rm{solid}}}}}\). Due to the incompressibility constraint in Eq. (1), \({{{\Phi }}}_{{{{\rm{solid}}}}}=1-{{{{\mathcal{V}}}}}_{0}{\bar{C}}_{0}-{{{{\mathcal{V}}}}}_{+}{\bar{C}}_{+}\) is a constant. Upon substituting the expression of the volumetric deformation through Eq. (1), we obtain \(d({C}_{0},{C}_{+})={d}_{0}{{{\Phi }}}_{{{{\rm{solid}}}}}^{-m}{\left[1+{{{{\mathcal{V}}}}}_{0}({C}_{0}-{\bar{C}}_{0})+{{{{\mathcal{V}}}}}_{+}({C}_{+}-{\bar{C}}_{+})\right]}^{m}\).

Equations for the static actuation of ionic polymer actuators

With the hypotheses of small deformations and electrochemical variables varying only along the y coordinate, Eqs. (15) and (16) become \(\frac{{{{\rm{d}}}}{J}_{{0}_{y}}}{{{{\rm{d}}}}y}=0\) and \(\frac{{{{\rm{d}}}}{J}_{{+}_{y}}}{{{{\rm{d}}}}y}=0\), such that \({J}_{{0}_{y}}\) and \({J}_{{+}_{y}}\) are constants that do not depend on y. Since the normal fluxes of solvent and counterions at the electrodes are zero, \({J}_{{0}_{y}}\) and \({J}_{{+}_{y}}\) must be zero as well throughout the thickness. For this condition to be identically fulfilled, we need the (electro)chemical potentials of the solvent μ0 and counterions μ+ to be constant,

$${{{\mathcal{R}}}}{{{\mathcal{T}}}}\log \left(\frac{{C}_{0}}{{C}_{0}+{C}_{+}}\right)+{\mu }_{{0}_{{{{\rm{slv}}}}}}+{\mu }_{{0}_{{{{\rm{els}}}}}}+{\mu }_{{0}_{{{{\rm{phy}}}}}}+{\mu }_{{0}_{{{{\rm{stc}}}}}}+{{{{\mathcal{V}}}}}_{0}\pi \equiv {\underline{\mu }}_{0},$$
(28)
$${{{\mathcal{R}}}}{{{\mathcal{T}}}}\log \left(\frac{{C}_{+}}{{C}_{0}+{C}_{+}}\right)+{\mu }_{{+}_{{{{\rm{slv}}}}}}+{\mu }_{{+}_{{{{\rm{els}}}}}}+{\mu }_{{+}_{{{{\rm{phy}}}}}}+{\mu }_{{+}_{{{{\rm{stc}}}}}}+{{{{\mathcal{V}}}}}_{+}\pi +{{{\mathcal{F}}}}\varphi \equiv {\underline{\mu }}_{+}.$$
(29)

Here, \({\mu }_{{(\cdot )}_{{{{\rm{slv}}}}}}=\frac{\partial {{{\Psi }}}_{{{{\rm{slv}}}}}}{\partial {C}_{(\cdot )}}\), \({\mu }_{{(\cdot )}_{{{{\rm{els}}}}}}=\frac{\partial {{{\Psi }}}_{{{{\rm{els}}}}}}{\partial {C}_{(\cdot )}}\), \({\mu }_{{(\cdot )}_{{{{\rm{phy}}}}}}=\frac{\partial {{{\Psi }}}_{{{{\rm{phy}}}}}}{\partial {C}_{(\cdot )}}\), and \({\mu }_{{(\cdot )}_{{{{\rm{stc}}}}}}=\frac{\partial {{{\Psi }}}_{{{{\rm{stc}}}}}}{\partial {C}_{(\cdot )}}\).

By neglecting the effect of the deformation on the dielectric displacement from Eq. (25)52 and substituting in Eq. (17), we find

$$-\epsilon \frac{{{{{\rm{d}}}}}^{2}\varphi }{{{{\rm{d}}}}{y}^{2}}={{{\mathcal{F}}}}\left({C}_{+}-{C}_{-}\right).$$
(30)

We now focus on satisfying equilibrium from Eq. (14). We consider small deformations, such that we can confound variables in the undeformed and deformed configurations. The stress tensor can be additively decomposed into three contributions: the mechanical stress σmec, Maxwell stress σpol, and the stress due to the incompressibility constraint σinc; σpol and σinc are so-called eigenstresses.

Due to our hypothesis on the shear strain and the form of eigenstresses, equilibrium along the y-axis requires that σyy does not depend on y. Since σyy must be zero at the electrode due to boundary conditions, we have that σyy = 0 everywhere. Thus, from the expression of σyy, we recover the through-the-thickness strain εyy as a function of the longitudinal strain εxx and electrochemical quantities,

$${\varepsilon }_{yy}=\frac{1}{{\lambda }_{{{{\rm{L}}}}}+2{\mu }_{{{{\rm{L}}}}}}\left[-{\lambda }_{{{{\rm{L}}}}}{\varepsilon }_{xx}+\pi -\frac{\epsilon }{2}{\left(\frac{{{{\rm{d}}}}\varphi }{{{{\rm{d}}}}y}\right)}^{2}\right].$$
(31)

We consider a Saint-Venant form of the solution by assuming that38 εxx(y) = −ky + ε0, where k and ε0 are the curvature and the mid-axis strain of the membrane, respectively. These quantities are obtained by requiring the longitudinal force and bending moment per unit width of the membrane to be zero, that is, \(N=\int\nolimits_{-H}^{H}{\sigma }_{xx}\,{{{\rm{d}}}}y=0\) and \(M=-\int\nolimits_{-H}^{H}{\sigma }_{xx}y\,{{{\rm{d}}}}y=0\), respectively, so that

$${\varepsilon }_{0}=\frac{\epsilon }{8{\mu }_{{{{\rm{L}}}}}H}\int\nolimits_{-H}^{H}{\left(\frac{{{{\rm{d}}}}\varphi }{{{{\rm{d}}}}y}\right)}^{2}\,{{{\rm{d}}}}y,$$
(32)
$$k=-\frac{3}{8{\mu }_{{{{\rm{L}}}}}({\lambda }_{{{{\rm{L}}}}}+{\mu }_{{{{\rm{L}}}}}){H}^{3}}\left[2{\mu }_{{{{\rm{L}}}}}\int\nolimits_{-H}^{H}\pi y\,{{{\rm{d}}}}y+({\lambda }_{{{{\rm{L}}}}}+{\mu }_{{{{\rm{L}}}}})\epsilon \int\nolimits_{-H}^{H}{\left(\frac{{{{\rm{d}}}}\varphi }{{{{\rm{d}}}}y}\right)}^{2}y\,{{{\rm{d}}}}y\right].$$
(33)

Substituting the expressions of εxx and εyy in the incompressibility in Eq. (1), we find

$${{{{\mathcal{V}}}}}_{0}({C}_{0}-{\bar{C}}_{0})+{{{{\mathcal{V}}}}}_{+}({C}_{+}-{\bar{C}}_{+})=\frac{1}{{\lambda }_{{{{\rm{L}}}}}+2{\mu }_{{{{\rm{L}}}}}}\left[2{\mu }_{{{{\rm{L}}}}}(-ky+{\varepsilon }_{0})+\pi -\frac{\epsilon }{2}{\left(\frac{{{{\rm{d}}}}\varphi }{{{{\rm{d}}}}y}\right)}^{2}\right].$$
(34)

The system of equations in Eqs. (28), (29), (30), and (34) constitutes an integro-differential system in the variables C0(y), C+(y), π(y), and φ(y). We complement this system with two boundary conditions for φ, \(\varphi (\!\pm\! H)=\pm \bar{V}/2\). Finally, to obtain the constant electrochemical potentials \({\underline{\mu }}_{0}\) and \({\underline{\mu }}_{+}\), we impose mass conservation of solvent and counterions, such that \(\int\nolimits_{-H}^{H}{C}_{0}\,{{{\rm{d}}}}y=2H{\bar{C}}_{0}\) and \(\int\nolimits_{-H}^{H}{C}_{+}\,{{{\rm{d}}}}y=2H{\bar{C}}_{+}\).

Numerical method

We solve the integro-differential system in Eqs. (28), (29), (30), and (34) with boundary and integral conditions by adapting the numerical method for two-point boundary value problems in ref. 41. To this end, we define a grid with N nodes yi, i = 1, …, N and N − 1 elements, with centroid \({y}_{j+\frac{1}{2}}=({y}_{j+1}+{y}_{j})/2\) and size \({{\Delta }}{y}_{j+\frac{1}{2}}={y}_{j+1}-{y}_{j}\). Note that the element size is non-uniform, whereby we consider a finer mesh close to the electrodes, where we anticipate the formation of boundary layers.

We need to solve for the 5N variables \({C}_{{0}_{i}}={C}_{0}({y}_{i})\), \({C}_{{+}_{i}}={C}_{+}({y}_{i})\), πi = π(yi), φi = φ(yi), and \({\xi }_{i}=\frac{{{{\rm{d}}}}\varphi }{{{{\rm{d}}}}y}({y}_{i})\), along with the unknown constants \({\underline{\mu }}_{0}\) and \({\underline{\mu }}_{+}\), for a total of 5N + 2 variables. We obtain 3N equations by evaluating the algebraic equations in Eqs. (28), (29), and (34) at each node yi. k and ε0 in Eq. (34) are evaluated from Eqs. (32) and (33), upon discretizing the integrals with the trapezoidal rule, which approximates the integral of any function f(y) on the [−H, H] interval as53\(\int\nolimits_{-H}^{H}f(y)\,{{{\rm{d}}}}y\approx \frac{1}{2}\mathop{\sum }\nolimits_{j = 1}^{N-1}\left[f({y}_{j+1})+f({y}_{j})\right]{{\Delta }}{y}_{j+\frac{1}{2}}\).

Other 2(N − 1) equations are found from discretization of the first-order system of differential equations equivalent to Eq. (30), \(\frac{{{{\rm{d}}}}\varphi }{{{{\rm{d}}}}y}=\xi\) and \(\frac{{{{\rm{d}}}}}{{{{\rm{d}}}}y}\left(\xi \right)=-{{{\mathcal{F}}}}\left({C}_{+}-{C}_{-}\right)/\epsilon\). We collocate these equations at each midpoint of each element, \({y}_{j+\frac{1}{2}}\). We approximate first derivatives in y with central differences over the element, and the value of functions at \({y}_{j+\frac{1}{2}}\) as the average of the function at the nodes of the element, according to \(\frac{{{{\rm{d}}}}f}{{{{\rm{d}}}}y}({y}_{j+\frac{1}{2}})\approx [f({y}_{j+1})-f({y}_{j})]/{{\Delta }}{y}_{j+\frac{1}{2}}\) and \(f({y}_{j+\frac{1}{2}})\approx [f({y}_{j+1})+f({y}_{j})]/2\).

The last four equations are obtained from the boundary and integral conditions. The boundary conditions are directly applied, such that \({\varphi }_{1}=-\bar{V}/2\) and \({\varphi }_{N}=\bar{V}/2\). In the integral conditions, we utilize the trapezoidal rule for discretizing the integrals, such that \(\mathop{\sum }\nolimits_{j = 1}^{N-1}\left({C}_{{0}_{j+1}}+{C}_{{0}_{j}}\right){{\Delta }}{y}_{j+\frac{1}{2}}/2=2H{\bar{C}}_{0}\) and \(\mathop{\sum }\nolimits_{j = 1}^{N-1}\left({C}_{{+}_{j+1}}+{C}_{{+}_{j}}\right){{\Delta }}{y}_{j+\frac{1}{2}}/2=2H{\bar{C}}_{+}\).

These discretized equations constitute a nonlinear system of 5N + 2 algebraic equations in 5N + 2 unknowns. We solve this problem numerically, by utilizing a nonlinear solver. The output of the solver are the values of \({C}_{{0}_{i}}\), \({C}_{{+}_{i}}\), πi, φi, ξi, \({\underline{\mu }}_{0}\), and \({\underline{\mu }}_{+}\), for i = 1, …, N. From these quantities, we can evaluate the curvature and the mid-axis strain of the membrane in Eqs. (32) and (33) and the contributions to the free energy in Eq. (2).

We solve the static problem for membranes with lithium and cesium counterions through two independent simulations. We utilize a mesh with N = 502 nodes and 501 elements, corresponding to a total of 2512 variables. Such a mesh includes a refinement in the proximity of the electrodes with a resolution of λ/20 for a region of thickness 10 λ, where \(\lambda =\frac{1}{{{{\mathcal{F}}}}}\sqrt{\frac{\epsilon {{{\mathcal{R}}}}{{{\mathcal{T}}}}}{{\bar{C}}_{+}}}\) is the Debye screening length54. We consider the static solution for a voltage of \(\bar{V}=1\,{V}_{{{{\rm{th}}}}}\), where \({V}_{{{{\rm{th}}}}}={{{\mathcal{R}}}}{{{\mathcal{T}}}}/{{{\mathcal{F}}}}\approx 25\,{{{\rm{mV}}}}\) is the thermal voltage.

We solve the nonlinear systems of equations in Matlab®, by utilizing the nonlinear solver fsolve. We use a nested structure, with an outer solver that finds the values of \({C}_{{0}_{i}}\), \({C}_{{+}_{i}}\), πi, φi, ξi, \({\underline{\mu }}_{0}\), and \({\underline{\mu }}_{+}\), and an inner solver that, given the tentative values of these variables, finds \({\mu }_{{0}_{{{{{\rm{slv}}}}}_{i}}}\) and \({\mu }_{{+}_{{{{{\rm{slv}}}}}_{i}}}\) through the implicit definition of α0. As an initial guess for the outer solver, we utilize the values of \({C}_{{0}_{i}}\), \({C}_{{+}_{i}}\), πi, φi, ξi, \({\underline{\mu }}_{0}\), and \({\underline{\mu }}_{+}\) from the ideal case. Once the values of \({C}_{{0}_{i}}\), \({C}_{{+}_{i}}\), πi, φi, ξi, \({\underline{\mu }}_{0}\), and \({\underline{\mu }}_{+}\) are known, we evaluate the changes in energy contributions in Eq. (2), as well as the curvature and mid-axis strain in Eqs. (32) and (33).

Material parameters

The semi-thickness of membranes is assumed to be 10−4 m, which is typical of commercial ionic membranes9. We set the temperature to \({{{\mathcal{T}}}}=300\,{{{\rm{K}}}}\).

Molar volume of water can be obtained from water’s molar mass and density, and equals 1.802 × 10−5 m3 mol−1. Molar volume of cations in solution is retrieved by previous studies, and corresponds to 1.5 × 10−6 m3 mol−1 for lithium and 2.2 × 10−5 m3 mol−1 for cesium counterions22.

Concentrations in the reference configuration are obtained from experiments with fully saturated membranes. For counterions’ concentration in the reference configuration, we use the data from ref. 55. Since that work does not consider actuators with cesium counterions, we use instead counterions’ concentration for actuators with sodium counterions. Thus, we select 1100 m3 mol−1 for lithium and 1200 m3 mol−1 for cesium. Water concentration in the reference configuration is obtained from absorption experiments with fully saturated membranes. Based on the ratio between water and counterions’ concentrations from ref. 56, we estimate water concentration to be 17,600 mol m−3 and 12,000 mol m−3 for membranes with lithium and cesium counterions, respectively. The reference volume fraction of solid is found from the values of \({\bar{C}}_{0}\) and \({\bar{C}}_{+}\) as Φsolid = 0.681 for lithium and Φsolid = 0.757 for cesium counterions.

The dielectric constant of the membranes is determined to maintain the Debye screening length equal to 10−7 m. This value, which is larger than the actual λ in Nafion™ membranes, is utilized to model high surface electrodes that significantly increase the capacity of the electric double layers57,58,59 and, indirectly, to enhance numerical convergence. For this value of the Debye screening length, we expect the continuum hypothesis to hold60. Due to the differences in \({\bar{C}}_{-}\) for membranes with lithium and cesium counterions, we obtain two slightly different dielectric constants of 4.106 × 10−5 F m−1 and 4.480 × 10−5 F m−1, respectively.

The mechanical properties of ionic membranes have been previously determined experimentally. We utilize the Young modulus of Nafion™ membranes in water saturated form measured in ref. 46. Since membranes with lithium counterions were not considered in that study, we utilize instead the Young modulus of membranes with sodium counterions. Thus, we use 80 × 106 Pa and 165 × 106 Pa for membranes with lithium and cesium counterions, respectively. Similar to our previous studies38,39, we set the Poisson ratio to 0.45.

The parameters \({{{{\mathcal{N}}}}}_{i}\) and ki for solvation are taken from ref. 13, such that \({{{{\mathcal{N}}}}}_{+}=5\), \({{{{\mathcal{N}}}}}_{-}=0\), k+ = 3.59 for lithium and k+ = 2.74 for cesium counterions, and k− = 0. We utilize the values from the same reference and its supplementary material for the parameters A = 1.177 m3/2 mol−1/2 and B = 3.291 m3/2 mol−1/2 nm−1 of Debye-Hückel theory, spacing d0 = 2.7 × 10−9 m between hydrophilic domains in the dry membrane, counterions’ diameter a = 3 × 10−10 m, fixed ions’ distance b = 4.7 × 10−10 m, physical interaction parameter β+− = 7.83 for lithium and β+− = −4.44 for cesium counterions (scaled by the molar mass of water), and the scaling parameter m = 1.33 for the pore deformation. Note that the expression of β+− is different from that of βij in ref. 13 since we incorporated M0, the molar mass of solvent, in the β+− coefficient to make it non-dimensional.