Electrochemo-poromechanics of Ionic Polymer Metal Composites: Towards the Accurate Finite Element Modelling of Actuation and Sensing

Ionic polymer metal composites (IPMCs) consist of an electroactive polymeric membrane, which is plated with metal electrodes and includes a fluid phase of ions in a solvent, whose diffusion allows for actuation and sensing applications. We build on a previous finite-deformation theory of our group that accounts for the cross-diffusion of ions and solvent and couples the mass balances of these species with the stress balance and the Gauss law. Here, we abandon the assumption that the fluid phase is a dilute solution, with benefits on both modelling and computation. A reliable finite element (FE) implementation of electrochemomechanical theories for IPMCs is challenging because the IPMC behaviour is governed by boundary layers (BLs) occurring in tiny membrane regions adjacent to the electrodes, where steep gradients of species concentrations occur. We address this issue by adopting the generalized FE method to discretise the BLs. This allows unprecedented analyses of the IPMC behaviour since it becomes possible to explore it under external actions consistent with applications, beside obtaining accurate predictions with a reasonable computational cost. Hence, we provide novel results concerning the influence of the membrane permittivity on the species profiles at the BLs. Additionally, by leveraging on the mobility matrix, we establish that the initial peak deflection in actuation strongly depends on the constitutive equations for the species transport and discuss the predictions of some experimental results from the literature. Overall, we demonstrate the potential of the proposed model to be an effective tool for the thorough analysis and design of IPMCs.


Introduction
(Asaka et al. [3], Shahinpoor and Kim [74]). The ionomer, which may be for instance made by Nafion™ or less often by Flemion™, is soaked with a solution including ions that can move within the ionomer macromolecules. These mobile ions are referred to as counterions because they have opposite sign with respect to those fixed to the ionomer backbone chains, usually consisting of anions.
In actuation, homogeneous electroactive elastomers, including ionomers, requires a too high applied voltage to exhibit a suitably large displacement because of their low dielectric permittivity. As nowadays well-known, this shortcoming can be overcome by appropriately designing composite materials including soft dielectric phases [8,14,46,55,57,72,76]. In IPMCs, this is accomplished by plating the membrane between thin sheets of noble metal, acting as electrodes. This forces the occurrence, during actuation and short-circuit sensing, of boundary layers (BLs) of very large variation of counterion concentration in thin membrane regions at the interfaces with the electrodes. Hence, most of the voltage drop across the IPMC actually occurs in the BLs, thus leading to a large electric field therein and, in turn, to a high polarisation stress triggering complex electrochemomechanics (Kilic et al. [48], Porfiri [67], Cha and Porfiri [28], Boldini and Porfiri [16], Boldini et al. [15]). It should be noted that, although the IPMC performance is very good in actuation, weak electrostatic interactions are observed for relatively large mechanical loads in IPMC sensing. This lack of "symmetric performance" in the dual basic components of the IPMC electrochemomechanics has to be ascribed to the IPMC nonlinear behaviour.
This investigation focuses on the thermodynamically-consistent modelling of IPMCs, which is an active area of research since the pivotal efforts of Nemat-Nasser and Li [62], Nemat-Nasser and Zamani [63], Wallmersperger et al. [88], Nardinocchi et al. [61], and Cha and Porfiri [28]. The essential features to be coupled in a multiphysical theory for IPMCs are the electrostatics, the mechanics, and the transport of species. Although the importance of the solvent motion has been pointed out since earlier contributions (Shahinpoor and Kim [75], Schicker and Wallmersperger [73], Zhu et al. [90]), only recent efforts have included its explicit and consistent modelling through additional flux law and mass balance with respect to those describing the counterion transport (Boldini and Porfiri [17], Bardella and Leronni [54], Olsen and Kim [64]). These theoretical studies, along with the recent experimental campaign of Arnold at al. [1], in fact demonstrate the importance of modelling the solvent flux, whose nonlinear interplay with the counterion flux sets the IPMC dynamic response.
Since the electrodes are substantially impermeable to mobile ions, electric double layers of charge form in the membrane at the interfaces with the electrodes [67]. The thickness of the electric double layers, which turns out to be comparable to that of the BLs, is on the order of the Debye screening length [7,10,20] where F ≈ 96,485 C/mol is the Faraday constant, ε is the absolute permittivity of the membrane, R ≈ 8.3145 J/(mol K) is the gas constant, θ is the absolute temperature, and C 0 is the initial nominal molar concentration of counterions, which is assumed to be uniform and coincident with the concentration of the anions fixed to the membrane backbones. Under this circumstance, the system constituted by the ionomer saturated with the solution of counterions is electroneutral. At room temperature, D may amount to a few tens of nanometres only for conventional membranes [67], while the IPMC thickness is on the order of magnitude of tenths of mil-limetres. Therefore, the ratio with H denoting the membrane semi-thickness, is extremely large. Large gradients of electric potential and counterion concentration in thin BLs have allowed leveraging on the magnitude of H/ D to develop useful analytical solutions, typically through the method of matched asymptotic expansions [85], both in short-circuit sensing [53,69,86,87] and in actuation [65,70,71]. These efforts, mostly relying on the electrochemomechanical theory of Cha and Porfiri [28], are however limited to relatively small loads (e.g., low applied voltages in actuation), since these solutions are acceptable estimates for mild nonlinearities only. Properly studying the IPMC nonlinear behaviour requires numerical methods [15,16], where the large value of the ratio in Eq. (2) poses an extremely challenging issue, still unaddressed to the best of our knowledge.
This numerical issue ensuing from the presence of BLs has been pointed out in other engineering applications involving the transport of charged species (Bauer et al. [9]) and the difficulty of modelling the electric double layer at solid-state electrochemical interfaces has been recently discussed by Swift et al. [81]. We are unaware of efficient computational models to accurately and consistently resolve the BLs in IPMCs [15,16,54,60].
In order to deal with this numerical issue, Narayan et al. [60] propose to fictitiously augment the BL size by largely increasing the permittivity (by a factor 10 7 ) in the Poisson equation, which, in fact, can be conveniently written in a form clearly displaying the role of D in setting the BL thickness [67]. However, this holds only in the linearised models, where the electrochemistry, which is governed by the classical Poisson-Nernst-Planck (PNP) system of equations, can be decoupled from the mechanics. We infer that including nonlinearities, because of either the application of large external actions or the need to account for more complex physics (typically in the species fluxes [54], or to account for non-ideal interactions [19]), may make the response more difficult to foresee by relying on basic physical models. Moreover, Narayan et al. [60] propose to maintain the "uninflated" material value of the permittivity in the constitutive law defining the polarisation (Maxwell) stress as a function of the electric field, in spite of its large modification in the Poisson equation, while other authors consistently augment the permittivity value in all the equations where it enters the adopted theory (see, e.g., [15,54] and references therein). This has led to a debate in the literature [20,60] on the appropriate way to deal with the numerical issue ensuing from the extremely small BL thickness for actual values of the IPMC permittivity. Boldini and Porfiri [20] have argued that, in actuation, the permittivity should be consistently augmented in all the equations governing the theory because of a compensation between the increase in the Maxwell stress for a given electric field and the increase in the BL thickness (proportional to √ ε, see Eq. (1)), the latter leading to a decrease in the electric field for a fixed voltage drop in the BLs. This conclusion of Boldini and Porfiri [20] not only relies on physical models such as the conventional form of the Poisson equation and the constitutive law for the Maxwell stress, but also on an implementation of the multiphysics theory proposed by Narayan et al. [60] including some of its nonlinearities.
In addition, we note that this issue is made even more complicated by the uncertainty on the actual value of the permittivity [20] and the lack of consensus on the definition of Maxwell stress [23]. The arguments of [20] hold if the Maxwell stress in the current configuration (i.e., its contribution to the Cauchy stress) reads σ pol = ε e ⊗ e − 1 2 (e · e)I , which is the most common expression for the Maxwell stress in the IPMC literature and is adopted in this investigation as well. 1 Here and henceforth, e is the electric field in the current configuration, ⊗ indicates the tensor product, · denotes the inner product, and I is the second order identity tensor. Moreover, all the foregoing models, adopted in order to draw conclusions on the role of the permittivity on the IPMC behaviour, disregard an explicit description of the transport of the solvent phase, which makes the whole picture much more entangled.
In this investigation we offer a contribution to this debate by proposing a robust finite element (FE) implementation of an electrochemo-poromechanical IPMC theory obtained by extending the proposal of [54], whose main feature, with respect to previous endeavors, is the modelling of the transport of both counterions and solvent, including their crossdiffusion. For the sake of brevity, henceforth, the theory of Leronni and Bardella [54] is referred to as LB21. Here, we extend this theory (in Sect. 2) by removing the hypothesis that counterions and solvent constitute a dilute solution. Hence, we change, with respect to [54], the contribution to the Helmholtz free energy accounting for the mixing of species and, more importantly, we enrich the constraint linking the membrane volume ratio to the species motion by also including the counterion transport.
About the proposed FE implementation (reported in great detail in Sects. 3, 4, and 5), its crucial feature is the use of generalized FEs (GFEs), whose implementation is a quite simple extension of that required for standard isoparametric FEs [39]. Basically, it needs the introduction of non-polynomial shape functions that enrich the capability to approximate the solution in the BLs. In this investigation we propose the use of a single enriching function that is designed to control the gradients of the primal fields in the BLs. In order to demonstrate the benefits of this implementation technique, we apply it both to the original LB21 theory (whose algorithm is reported in Appendix C) and to its extension developed in this work. To the best of our knowledge, the numerical models developed in this study are the first ones in the literature that allow the analysis of the IPMC behaviour by linking it to a detailed descriptions of the BLs for applied external actions consistent with applications and experiments.
In order to demonstrate our findings, first, in Sect. 6, we validate the GFE algorithm against the numerical results presented in [54], where the theory was implemented in COM-SOL Multiphysics ® . Notably, our implementation allows us to obtain accurate numerical results for cases that could not be solved in [54]. More specifically, both in actuation (Sect. 7.1) and short-circuit sensing (Sect. 7.2), we can obtain results for much larger applied loads, whose values match those of real applications. For instance, with the model parameters of [54], in actuation the proposed GFE model converges for applied voltage drop across the electrodes ψ 0 = 2 V, against the maximum value of 0.25 V that could be applied with the implementation in [54].
Also, we can obtain results for a wide range of model parameters, which was not possible with the implementation in [54]. This allows us to establish some shortcomings of the LB21 theory that are fixed by the extension developed in this work, as discussed through the numerical results of Sect. 7.3.
Then, in Sect. 8 we elaborate on the influence of the permittivity ε on the actuation response, with focus on both the scaling of the IPMC response with ε and the capability of the proposed GFE models to provide results for low ε. We highlight that the scaling expected from inspection of Eqs. (1) and (3) holds only for suitably small applied voltages, thus for instance avoiding that at the BL subject to counterion depletion the counterion concentration vanishes, which would trigger a quite strong nonlinearity leading to different thicknesses of the two BLs, followed by back-relaxation [70]. We also discuss the form of the mobility matrix, which governs the constitutive laws for the solvent and counterion fluxes, showing that the form proposed in [54], including cross-diffusion and convection, delivers better predictions than the simplest diagonal form, the latter being however the sole one granting us to observe, not only at steady state, the scaling with ε argued by Boldini and Porfiri [20]. This comparison also allows us to demonstrate that the constitutive laws for the fluxes have an enormous impact on the prediction of the initial peak deflection in actuation.
Finally, in Sect. 9 we discuss the capability of the new theory implemented through GFEs to provide reasonable predictions of the experimental results of He et al. [42] on the actuation of Nafion™-Pt IPMCs subjected to applied voltage up to 3 V.
As detailed in the concluding remarks of Sect. 10, although several important open issues remain to be addressed, we expect the present proposal to constitute a solid basis towards the development of a reliable tool for the analysis and optimal design of IPMCs.

Kinematics and Balance Equations
The theory assumes that the IPMC membrane consists of a solid phase identifying with a charged elastomer that is soaked with a fluid phase consisting of an uncharged solvent in which counterions are immersed. Analogously to mixture theory, all the phases coexist within each single material point (Ateshian [4]) and the fluid phase is assumed to always saturate the membrane. The mechanical deformation of each macroscopic material point is expressed in terms of its deformation gradient where u = x − X is displacement vector, ∇ denotes the material gradient, such that (∇u) iJ = ∂u i /∂X J , 2 and x and X are the current and the reference (material) position vectors. The reference (initial) configuration is undeformed and electroneutral. In [54] it is assumed that the solution of counterions and solvent is dilute, that is, C C s , with C and C s indicating the nominal molar concentrations of counterions and solvent, respectively. Importantly, here we disregard this restriction. As specified in detail in the following, this requires to modify a few key governing equations of the theory. However, the major impact on the predictive capability of the model undoubtedly lies in the kinematic constraint where J is the volume ratio, v is the hydrated counterion molar volume (if the solvent is water), v s is the solvent molar volume, and C 0 and C 0 s are the initial values of C and C s , respectively.
Equation (5) assumes that polymer chains, solvent molecules, and counterions are intrinsically incompressible, and, contrary to the simpler proposal in [54], establishes that the volume ratio depends also on the counterion motion. This is crucial for the purposes of this investigation, aiming at clarifying the role of the fluxes of counterions and solvent in the IPMC response, which is strictly related to the capability of these species to either deplete or accumulate at the BLs. To this purpose, Eq. (5) accounts for the volume fraction of the BL region occupied by the counterions, thus describing a physically more appropriate picture than that in [54], disregarding the term v(C − C 0 ). Note that v is an extra material parameter with respect to the LB21 theory. Moreover, since v is the molar volume of hydrated counterions, C s is the concentration of the water free to move independently of the counterions only, which does not represent the total molar concentration of the solvent.
The rate form of Eq. (5) is where( ) indicates partial time derivative, that is,( )(X, t) ≡ ∂( )(X, t)/∂t . As in [54], the membrane behaviour is governed by four balance equations written with respect to the reference configuration. They are the overall linear momentum balance in the absence of body forces Div P = 0 , the solvent mass balanceĊ the counterion mass balanceĊ and the Gauss law Here and henceforth, P is the nominal stress tensor and Div is the material divergence, such that (Div P) i = ∂P iJ /∂X J ; J s and J are the nominal molar fluxes of free solvent and hydrated counterions, respectively; D is the nominal electric displacement. Among the assumptions behind these balance equations, we mention the rapidity in attaining the mechanical equilibrium with respect to the time scale of the fluid phase transport (Hong et al. [45], Cha and Porfiri [28]) and the unit valency assigned to both fixed and mobile ions, which is usually the case for IPMCs.

Assumptions on the Electrode Behaviour
The balance of the electrodes relies on Eq. (7) only, because they are assumed to be perfect conductors whose interface with the membrane is impermeable to solvent and counterions. Moreover, they are assumed to undergo purely elastic deformations governed by the Fig. 1 Cantilever IPMC subjected to an applied voltage ψ 0 across the electrodes (non-zero in actuation) and to an imposed nominal surface load T 0 (non-zero in sensing): geometrical parameters, reference system, and physical constituents. Taken from [54], where there is no distinction between free solvent molecules and solvent molecules fixed to the counterions isotropic Saint-Venant-Kirchhoff strain energy density where E = (C − I)/2 is the Green-Lagrange strain tensor, C = F T F is the right Cauchy-Green deformation tensor, and λ e = E e ν e /[(1 + ν e )(1 − 2ν e )] and G e = E e /[2(1 + ν e )] are the Lamé constants, here expressed in terms of the Young modulus E e and the Poisson ratio ν e . Although IPMCs may undergo large displacement fields, the strains in the electrodes are expected to be relatively small, such that any isotropic elastic law involving both Lamé constants should provide about the same response. From the FE implementation perspective, we find it convenient to adopt Eq. (11).

Geometry and Initial and Boundary Conditions
With reference to Fig. 1, the membrane thickness is 2H , each electrode thickness is h, and the IPMC length is L. Usually, the IPMC is a slender sandwich structure with thin skins, i.e. L H h. The IPMC is clamped at the left end, such that u = 0 at X = 0, and is stress-free at the edges X = L and Y = H + h. The remaining side is either stress-free in the actuation problem or subjected to in the sensing problem, where T 0 (t) is the time-varying amplitude of the distributed nominal load. The electrode impermeability requires vanishing fluxes at the interfaces, that is Because of the IPMC slenderness, the edge effects at the IPMC ends are unimportant, thus allowing us to impose zero-flux at X = 0 and X = L as well.
The electrostatic boundary conditions are with ψ denoting the electric potential and ψ 0 (t) the time-varying applied voltage drop across the electrodes, which is zero in the short-circuit sensing problem. At the sides X = 0 and X = L it is convenient to prevent charge accumulation by enforcing D X = 0. Finally, the initial conditions read the latter ensuring electroneutrality at rest.

Internal Energy and Dissipation
Here we follow the Coleman-Noll procedure as in [54], with some differences due to adopting the constraint (5), instead of (80) as assumed in [54]. By introducing the nominal electric field the solvent chemical potential μ s , the counterion chemical potential μ, and the counterion electrochemical potentialμ the rate of the nominal internal energy in the membrane readṡ Next, we impose non-negative dissipation by postulating the existence of a nominal Helmholtz free-energy density W : where π , physically describing the osmotic pressure field, is a Lagrange multiplier [29,44,89] needed to account for the constraint (5). In [54], because of the dilute solution approximation, the analogous Lagrange multiplier (therein denoted as p s ) is the solvent pressure.
As in [54], W is a function of F, D, C s , and C. Hence, Eq. (19) becomes from which we can establish the general constitutive relations where, differently from [54], the counterion chemical potential depends also on the Lagrange multiplier. Then, from Eqs. (19) and (20), the second law of thermodynamics results into the constraint The standard way to fulfill this dissipation inequality consists of writing each flux as a linear combination of ∇μ s and ∇μ:

Specific Form of the Membrane Free-Energy Density
While it is postulated that the last two contributions toU in (18) have dissipative nature, P, E, μ s , and μ are obtained by deriving with respect to F, D, C s , and C the Helmholtz free energy which is assumed to consist of the three contributions W mec , W mix , and W pol separately accounting for the membrane deformation, the mixing of free solvent molecules and hydrated counterions, and the membrane polarisation, respectively. A more sophisticated free energy, considering non-idealities, has been recently put forward by Boldini and Porfiri [19]. As in [54], the membrane elasticity is governed by the isotropic compressible Neo-Hookean model proposed by Simo and Pister [79] where I 1 = tr C and λ and G are the Lamé parameters of the membrane, having Young modulus E and Poisson ratio ν. The coupling of the deviatoric and volumetric behaviours in Eq. (29) is a feature usually needed in order to capture the elastomeric large deformation response (Holzapfel [44], Boyce and Arruda [22]). As shown in [54], the use of a model involving both Lamé parameters allows a better control of the interaction between the solvent flux and the volumetric deformation, which is important in the adopted poromechanical context, where it plays a significant role in the counter-diffusion phenomena occurring in actuation and sensing. In the present theory, the availability of two elastic constants is even more important because constraint (5) also relates the volume ratio to the hydrated counterion motion. 3 We remark that the membrane elastic deformation could be constrained to be isochoric by adopting a richer kinematics in which F would involve inelastic contributions, such as the viscoplastic contribution proposed by Silberstein and Boyce [77] to describe the Nafion™ behaviour, or the swelling contribution adopted by Narayan et al. [60] to model IPMC actuation. In the terminology of Giantesio et al. [40], this swelling contribution to F is an "active strain" that accounts for osmosis, while in the present theory we follow Hong et al. [45] and Cha and Porfiri [28], thus adopting an "active stress" framework for both polarisation and osmosis.
By assuming that the fluid phase is an ideal solution of solvent and counterions (Fermi [37], Ateshian [4], Kondepudi and Prigogine [50]) that displays purely entropic behaviour, the free energy of mixing may be specified as which is consistent with the assumed incompressibility of all the membrane phases (Boldini and Porfiri [18]) and is different from the form employed in [54], where W mix (C, C s ) is simplified under the hypothesis that C C s . The membrane is assumed to behave as an ideal dielectric, whose polarisation free energy reads Next, we write the constitutive laws for the recoverable quantities, which are obtained from the foregoing equations.
The total nominal stress in the membrane reads with denoting the mechanical contribution, denoting the osmotic contribution proportional to the Lagrange multiplier π , and denoting the polarisation contribution, that is usually referred to as the Maxwell stress. It is convenient to write the total pressure p (positive if compressive) in terms of the Cauchy stress σ = PF T /J : where d = J −1 FD = ε e is the current electric displacement (Dorfmann and Ogden [33]). By inverting (20b) and using (16), the nominal electric displacement reads The solvent chemical potential is obtained by combining Eqs. (20c) and (30) into whereas the counterion chemical potential μ ensues from combination of Eqs. (20d) and (30), then leading to the counterion electrochemical potential (17) in the form

Non-dimensional Electrochemical Primal Variables
In view of the FE implementation, it is convenient to introduce non-dimensional primal variables, as customary when analytically solving linearised electrochemomechanical problems for IPMCs (see, e.g., [15,65,86] and references therein). Selecting the optimal set of them is a non-trivial task that we have accomplished by evaluating several numerical aspects, here unreported for the sake of brevity. First, the non-dimensional chemical potentialμ is defined as By adopting the symbol r for the ratio between counterion and solvent concentrations, we can re-write Eq. (36) in the simple and compact form Also, for later reference, it is convenient to combine Eqs. (37) and (5) to determine C s as a function of J and r: Second, the non-dimensional osmotic pressureπ is defined as Third, the non-dimensional electric potentialψ is defined through the thermal voltage Equations (40) and (41) allow the solvent and counterion chemical potentials to be written as Therefore, the counterion electrochemical potential (17) results On the basis of definitions (42) and (44), the gradients of the potentials that determine the fluxes (23) are Finally, we introduce the non-dimensional timē where D is the counterion diffusivity entering the coefficients (24) of the mobility matrix M , as explained next.
In this investigation, we consider where D s is the solvent diffusivity and the non-dimensional parameter γ can switch on and off the cross-diffusion and the counterion transport by convection. Setting γ = 1 recovers the form adopted in [54], while γ = 0 leads to the simplest possible diagonal form of M , suppressing both cross-diffusion and counterion transport by convection. From Eq. (25), in order to ensure non-negative dissipation, γ must be limited as γ (γ − 1) < DC s /(D s C), which is straightforwardly satisfied in the cases here of interest γ = 1 and γ = 0.
The form of (47) of M is convenient as we will resort to the case γ = 1 both to validate our numerical results against those in [54] and to demonstrate that our implementation allows overcoming the limits of the FE model in [54], while we will assume the simple diagonal form of (47) (that is, γ = 0) in Sect. 8, to discuss the outcomes of other IPMC models from the literature, with focus on the expected scaling of the BL size with the permittivity value, as suggested by Eq. (2).
Let us clarify that within the LB21 theory (γ = 1), consistently with the dilute solution approximation, D s is the solvent diffusivity in the polymer network and D is the counterion diffusivity in the solvent. Instead, in the present theory, where the solution is not necessarily dilute, both D s and D assume the role of diffusivity of a species in the medium constituted by the polymer network and the other species.
We observe that in the context of non-dilute solutions there is room for more sophisticated constitutive laws than those in Eq. (47) for the coefficients of the mobility matrix (Vanag and Epstein [84], Kondepudi and Prigogine [50]). However, the prescriptions in Eq. (47) are appropriate for the purposes of this investigation.
By combining Eqs. (23), (45), and (47), we obtain the fluxes: As in [54], setting γ = 1 implies that the solvent flux is independent of ∇μ (i.e., independent of both ∇C and ∇C s ), while it includes the electro-osmosis through its dependence on ∇ψ, which some investigators indicate as an important contribution to solvent transport in IPMC actuation (Asaka and Oguro [2], Shahinpoor and Kim [75], Zhu et al. [90]). The opposite is obtained by choosing γ = 0, whereby J s depends on ∇μ and there is no electro-osmotic contribution. Moreover, rearranging the contributions in Eq. (49) allows one to single out the convective contribution to the counterion flux with the solvent [54], which reads CJ s /C s and plays a role in sensing, where a mechanical load moves the solvent along its pressure gradient, transporting some counterions. The solvent motion sets a volume ratio gradient that leads to an additional counterion Fickian diffusion governed by the terms dependent on ∇C s in Eq. (49), which can be unveiled by substituting definition (36) in Eq. (49). Instead, the value of γ has much less impact on the functional form of the counterion flux, always displaying dependence on ∇μ , ∇π , and ∇ψ . Specifically, in the limit case of fixed solvent (D s = 0) and for γ = 1, in Eq. (49) it is possible to recognise the classical Nernst-Planck law, written with respect to the reference configuration, establishing that both ∇C (Fick effect) and ∇ψ (electrophoretic effect) concur to the counterion migration (Porfiri [67]).

Weak Form of the Balance Equations
For the nodal degrees of freedom (DOFs), we select the displacement vector u along with the non-dimensional electric potentialψ , chemical potentialμ , and osmotic pressureπ , as defined in Eqs. (41), (36), and (40), respectively. It should be noted that adopting a nondimensional displacement such as u/H does not change much at all the numerics. When discussing the results, we will comment on other possible choices as a function of the BLs size.
The weak form of the balance equations is obtained by integrating over the whole reference domain of the membrane the variation of its internal energy rate provided in Eq. (18). By further making use of the equations presented in Sect. 2, we obtain where δ denotes an admissible variation of the field it is applied to andJ and C s are obtained from Eqs. (6) and (39) in terms of the time and spatial derivatives of u (involving the deformation gradient). The FE formulation adopted in the electrodes is standard, and it is reported in Appendix A for the sake of completeness.

Isoparametric FE Formulation and Array Notation for the IPMC Two-Dimensional Problem
The conventional FE discretisation relies on the shape functions where N denotes the number of nodes of each FE and ξ ∈ [−1, 1] and η ∈ [−1, 1] are the intrinsic coordinates of the parent element. We use the same shape functions to interpolate all the nodal DOFs, by adopting isoparametric FEs such that whereX i andŶ i indicate the nodal values of the coordinates in the reference configuration and, as usual, the polynomial shape functions build a partition of unit, i.e.
The shape function vector is defined as with ( ) T denoting the transpose operation. Hence, the non-dimensional discretised potential μ and osmotic pressureπ read, respectively, whereμ 0 = ln(C 0 /C 0 s ) and the vectorsμ andπ π π, both having dimension N , contain the nodal DOFs. Note that, to avoid a too complex notation, for these two scalar quantities we adopt the same symbols for the exact and the discretised fields.
The displacement is described by the vectorial field whose values at the nodes are collected in the vector In the FE implementation, the second order tensors are stored in a single column array, whereas fourth order tensor components are stored in two-dimensional arrays. By considering the deformation gradient F under plane-strain conditions, the relevant components of the matrix [F] are converted into a single-column array of dimension 4 according to the rule Here and henceforth, the symbol → denotes the conversion from tensorial notation to array notation, the latter being exactly that followed in the implementation, as in de Souza Neto et al. [31]. The deformation gradient can then be computed as a function of the displacements at nodes as The volume ratio is written as (54) and the relevant components of the right Cauchy-Green deformation tensor are collected into The two non-vanishing components of the gradient of scalar quantities are computed through the intrinsic coordinates ξ , η by using the matrix operator such that, for example, the gradient of the non-dimensional electric potential reads where the components of the vectorψ are the nodal values ofψ .
The nominal stress, electric displacement, solvent flux, and counterion flux are, respectively, represented by the vectors By referring to the non-dimensional time defined in Eq. (46), the incremental initial value problem aims at evaluating all the unknowns at the end of a suitably small time step t n ,t n+1 , where all quantities are known at timet n . All quantities at timet n are indicated with a subscript n, while, to keep the notation as simple as possible, we drop the analogous subscript n + 1 for the quantities evaluated att n+1 .
We adopt a fully implicit time-integration scheme where the increment of a generic variable q is defined as Upon integration over the time step the discretised increment of internal energy (50) reads where, by inverting Eq. (46), In the present notation, the inverse of the tensor C can be computed as: It is convenient to further introduce the matrix-vector product between the generic plane tensor τ and 2D vector v in array notation: By using Eqs. (33) and (56), the nominal electric displacement can then be computed as By computing and the mechanical contribution to the total nominal stress P (32) is while the osmotic stress reads and the Maxwell stress can be written as where the second order tensor ω is defined in terms of the electric displacement D as and • denotes the matrix product between the two second order tensors F and ω, with the latter a dyadics as in (64), such that The solvent concentration C s is computed as in (39) and the fluxes can be written in array notation as The internal pseudo-forces read The consistent "stiffness matrix" K, of dimension 5N × 5N , required to complete the FE algorithm and ensure second order convergence of the global Newton-Raphson scheme is defined as the differentiation of the right-hand sides of Eqs. (67)-(69) with respect to the incremental nodal variables and reads whose submatrices are provided in Appendix B.

Generalized Finite Element Extension
In order to obtain a numerical model able to accurately represent the IPMC response, we adopt the generalized FE (GFE) method (GFEM, see, e.g., Fries and Belytschko [39] and references therein) to enrich the model in thin membrane regions adjacent to the electrodes, where high gradients of most physical quantities occur. In fact, resolving these BLs, whose size is on the order of magnitude of the Debye screening length D (1), through classical FE requires too tiny elements in the BL regions. This causes fatal floating point errors due to the extremely large ratio between the IPMC thickness 2H and D (see Eq. (2)). Alternatively, the use of an appropriately regular mesh would require a number of standard FEs resulting in an unbearably large numerical cost. Following Fries and Belytschko [39], we adopt the GFE approximation of a generic function g(ξ, η), whereâ i are conventional nodal DOFs, N i (ξ, η) are standard polynomial functions as in Eq. (51),b i are additional nodal DOFs, and f (Y ) is referred to as the enrichment function, whose role is to bring into the FE discretisation some knowledge of the solution. In the case of IPMCs, f (Y ) accounts for the high gradients in the BLs and it is convenient to refer it to the real reference coordinate Y .
Since we consider the standard FE geometric discretisation (51), for which Eq. (52) holds in the whole domain , it can be demonstrated that Eq. (71) can exactly reproduce any enrichment function f (Y ) in any subset of discretised by GFEs (Babuška and Melenk [6,59]).
One of the main complications arising with the GFEM is that, obviously, the field value g at i-th node does not coincide withâ i nor withb i . For this reason, special care must be taken to ensure the continuity at the interfaces between standard and enriched elements, whereb i must be constrained to vanish to ensure a smooth variation of the field g. Such interfaces are located at Y = ±H (see Fig. 1) and Y = ±(H − ε ), where ε ≥ D is the thickness of the rectangular domains (of length L) that are discretised by GFEs.
Moreover, it is desirable that the Dirichlet boundary conditions (13) and (14) could be imposed by directly setting the valueḡ at Y = ±H .
The foregoing issues are addressed by resorting to the so-called "δ ij property", which consists in properly shifting f (Y ) such that, at each i-th node, the nodal DOFb i multiplies a vanishing enrichment function. This requires the modification of Eq. (71) into It should be noted that this modification, strongly advised by Fries and Belytschko [39], does not compromise the capability of reproducing f (Y ) in the GFE domain (Belytschko et al. [11]). Although more enrichment functions f i (Y ) (with i = 1, . . . , N) could be introduced in the discretisation, here we adopt a single enrichment function. The derivatives of g(ξ, η) can then be computed with respect to the intrinsic coordinates as Finally, the computation of the derivatives of g(ξ, η) with respect to X, Y is totally standard, requiring the inverse of the FE Jacobian matrix determined in terms of ∂N i (ξ, η)/∂ξ , ∂N i (ξ, η)/∂η, and the nodal coordinatesX i ,Ŷ i .
The enrichment function f (Y ) plays a fundamental role in the predictive capability of the GFE model. Here, f (Y ) is obtained by modifying the enrichment function proposed by Benvenuti et al. [12,13]. Specifically, we assume where α is a non-dimensional parameter to be tuned in order to represent the gradients of the relevant fields occurring in the BLs, over which f (Y ) behaves as shown in Fig. 2 for a specific choice of α, ε , and H . An important feature of this f (Y ) is that it tends to become spatially uniform for |Y | → H − ε , i.e., far enough from the high-gradient region for a suitable choice of ε > D . This strategy moves the problematic "linear dependencies" that may affect the GFEM away from the region whose behaviour is difficult to capture and governs the IPMC response. Such "linear dependencies" refer to redundancy in the spatial discretisation [39] and may lead to ill-conditioned global tangent "stiffness" matrix. This is an important issue in commercial codes, such as ABAQUS/Simulia [30] adopted in this investigation, where the user may not control the solver preconditioner nor the accuracy in solving linear systems.
It should be noted that the proposed f (Y ) in Eq. (73) has the unit of length. Although this leads to an annoying dimensional non-homogeneity of theâ i andb i DOFs, it grants a good control of the gradients in the BLs, because df/ dY is non-dimensional. Our numerical experiments have suggested Eq. (73) as the best option among several attempts we tried for f (Y ). Specifically, definition (73) for f (Y ) allows us to avoid as much as possible illconditioned matrices, which may occur for non-dimensional definitions of f (Y ), lacking a reliable control of the gradients.

Spatial Integration
Due to the high nonlinearity of the enrichment function (73), a large number of Gauss points must be adopted in the GFEs for the numerical integration of both the pseudo-forces and the stiffness matrices. However, the nature of the expected solution allows us to save some computational cost by adopting a non-isotropic distribution of the Gauss point positions.
Specifically, each GFE has 3 Gauss points along the X direction and 20 Gauss points along the Y directions. 5 Although this simple approach does not introduce any particular issue at the implementation level, it leads to a convoluted post-processing of the results. This is because contourplots need the field variables computed in the Gauss points to be mapped at nodes, and their accurate evaluation requires the use of the extended shape functions, usually unavailable in standard post-processing codes.

Enriched DOFs and GFE Operators
In the proposed FE model, the enriched discretisation (72) is adopted for the four scalar fields u y ,ψ ,μ , andπ . By denoting with the enriched shape functions referred to the i-th node and by defining the array of size 2Ñ the chemical potential and the osmotic pressure are discretised as where the 2N nodal unknowns are stored in the two arrays wherẽ Each node of the enriched FE has then 9 DOFs. The internal pseudo-forces and the consistent "stiffness matrix" (now of dimension 9N × 9N ) for the GFEs have the same forms as those in Eqs. (67)-(69), (70), and Appendix B, except that the operators N, G, and A must be simply substituted withÑ,G, and A, respectively.
The FE algorithms employed in this investigation have been implemented in user element subroutines (uel) for the commercial FE code ABAQUS/Simulia [30].

Spatial Discretisation
The simple geometry of IPMCs allows their FE mesh to be "structured", i.e., constituted by rectangular elements. For the standard FEs, we always adopt quadratic shape functions that, in conjunction with the structured mesh, minimise the integration errors of both the internal forces and the residuals, even when we have to resort to FEs with "bad" aspect ratio to reach the required fine discretisation along the membrane thickness without unnecessarily enriching the discretisation along the IPMC axis.
In order to have control on the numerical results, we have analysed the IPMC behaviour by using several meshes. Here, we mainly report the results concerned with two of them, which are referred to as "coarse" and "fine" meshes in the following. Figure 3 shows some details the "coarse" mesh. Each electrode is discretised along its thickness with four 8-noded quadratic isoparametric FEs with reduced integration (4 Gauss points per FE). Eight rows of GFEs (displayed in dark grey in the right image of Fig. 3) discretise the membrane regions, of size equal to ε along the Y direction, next to the electrodes. The remaining part of the membrane is discretised along Y , overall by fifty rows of elements, by using the following double bias strategy. In most analyses, we have set a ratio ≈ 1800 between the largest and smallest sizes of the FEs, where the maximum thickness occurs at the IPMC mid axis while the minimum thickness is attained next to the GFE regions; in other analyses, which Fig. 3 refers to, the minimum and maximum sizes of the FEs are H/100 (at |Y | = H − ε ) and H/10 (at Y = 0), respectively. We have thus established that these kind of details on the discretisation of the bulk membrane have little impact on the quality of the numerical results.
Along the X direction, we divide the IPMC in two regions characterised by two different element widths. The region from the clamped cross-section (X = 0, see Fig. 1) to X = 2H (that is, the length of the membrane thickness) is divided into 16 equally spaced FE columns. In the remaining part of the domain the element width in the X direction is equal to the membrane semi-thickness H . This discretisation along the IPMC axis is common to the all the meshes and, overall, in the case L/H = 200 (that is the case of all our simulations except those of Sect. 9) it leads to meshes constituted by 214 columns of elements.
In summary, the "coarse" mesh has 48, 085 nodes and consists of: 1, 712 standard isoparametric FEs (with 2 DOFs for each node) with quadratic shape functions and reduced integration for the electrodes; 10, 700 quadratic, fully integrated isoparametric FEs (5 DOFs for each node) to discretise the bulk membrane region; 3, 424 quadratic GFEs (with 9 DOFs for each node) to discretise the BL regions. The total number of DOFs is 269, 617.
The features of the "fine" mesh are displayed in Fig. 4. This mesh differs from the "coarse" one of Fig. 3 in two aspects. First, the BL regions are discretised with twelve rows of GFEs, instead of eight. Second, the bulk membrane is discretised by using seventysix rows of elements (instead of fifty), with minimum and maximum sizes of the FEs H/500 (at |Y | = H − ε ) and H/10 (at Y = 0), respectively. Overall, this "fine" mesh has 69, 981 nodes and consists of: 1, 712 quadratic isoparametric FEs with reduced integration for the electrodes; 16, 264 quadratic, fully integrated isoparametric FEs plus 5, 136 quadratic GFEs to discretise the membrane. The total number of DOFs is 399, 705.
In order to thoroughly discuss the results, we also consider a mesh with standard FEs only. In this model, the mesh geometry is identical to the "coarse" mesh, where the GFEs are substituted with the same quadratic isoparametric FEs (5 DOFs for each node) used in the membrane outside the BLs. The total number of DOFs is in this case equal to 224, 969.

Computational Cost
The aim of this subsection is to give information on the efficiency of the proposed GFE algorithm. The numerical analyses have been performed by using a workstation equipped with 256 GB of RAM and an AMD Ryzen Threadripper 3970X having 32 cores, 64 threads, and a CPU clock speed of 3900 MHz. Fig. 3 Details of the adopted "coarse" mesh: the left image shows the discretisation over the whole IPMC thickness, whereas the right image shows a detail about the bottom electrode (whose domain is displayed by the bottom four rows in light grey) with the adjacent membrane region including the GFE discretisation, the latter being displayed in dark grey Fig. 4 Details of the adopted "fine" mesh: the left image shows the discretisation over the whole IPMC thickness, whereas the right image shows a detail about the bottom electrode (whose domain is displayed by the bottom four rows in light grey) with the adjacent membrane region including the GFE discretisation, the latter being displayed in dark grey The most critical analysis reported in this paper is that in Sect. 8 for the IPMC actuation at applied voltage ψ 0 = 3 V with permittivity ε = 10 −13 C/(mV mm). In this case, the required computational time to run the analysis untilt = 5 amounts to 16 min and 36 sec to complete 154 loading steps and 339 iterations.
The results of the most expensive analysis, still reported in Sect. 8, are concerned with the IPMC actuation at ψ 0 = 1 mV with ε = 10 −12 C/(mV mm) and counterion molar volume v = 500 mm 3 /mol, run up tot = 100. In this case, the required computational time amounts to 1 h, 7 min and 39 sec to complete 1011 loading steps and 1226 iterations. In this section we validate the proposed numerical model by comparison with the results obtained in [54], where the theory was implemented in the commercial code COMSOL Multiphysics ® and the analyses, both in actuation and short-circuit sensing, were run for relatively small applied loads and limited ranges of material parameters. For the details of the FE model adopted in [54] the reader is referred to that paper; here, we only report its total number of DOFs equal to 329, 133, which is comparable to that of the models developed in the present study.
Our new FE algorithm employed to implement the LB21 theory, concerned with the results presented here and in Sect. 7, is slightly different from that presented in previous Sects. 3 and 4, as it implements a theory with a different constraint governing the volumetric deformation, as reported in Eq. (80), and a different expression for the free energy of mixing (specifically, Eq. (81), which is the dilute solution approximation of Eq. (30)). For the sake of completeness, in Appendix C we also provide our GFE algorithm of the LB21 theory.

Geometrical and Material Parameters
The geometrical and material parameters are the same as those selected in [54]. With reference to Fig. 1, the IPMC thickness is 2H = 0.2 mm, the IPMC length L is equal to 20 mm, and the electrode thickness h is equal to 0.001 mm. The material parameters are reported in Table 1, which includes also the GFE parameters required by the proposed numerical model. In all the analyses the temperature is set to θ = 300. K.
With the parameters of Table 1, application of Eq. (1) provides a Debye screening length D ≈ 1.5 · 10 −4 mm = ε /6.6. In Sect. 8, we will discuss the impact of changing ε and we will present results for lower permittivities, up to three orders of magnitude lower than that of Table 1. 6

Actuation
In actuation, [54] presents results for an applied step voltage drop across the electrodes ψ 0 = 0.25 V, which is about the maximum that could be applied with the FE implementation 6 For Nafion™-based IPMCs, Wallmersperger et al. [88] have identified a relative permittivity ε/ε 0 = 120, which results in ε ≈ 10. −15 C/(mV mm) and, from Eq. (1) and with the parameters of Table 1, D ≈ 4.7 · 10 −7 mm. However, in IPMCs, ε should be increased to account for its non-homogeneity. Specifically, in membrane regions adjacent to the membrane-electrode, referred to as "intermediate layers" [82] or "composite layers" [26,27] in the literature, metal particles from the electrode manufacturing partly fill the interstices between the ionomer macromolecules, thus resulting in a much larger ε and a much lower diffusivity than in the bulk membrane [71,86]. Additionally, given the large through-the-thickness direct strain component in the BLs, the deformation-dependent permittivity observed in some electroactive polymers could play a role [52,56].   [54]. In all the analyses this voltage drop is instantaneously applied. Figures 5, 6, 7, 8, and 9 compare the results obtained from our numerical model with that of [54] in terms of the time-evolution of the deflection and, for the three different instants t = 0.01 s, 0.1 s, and 10 s, the through-the-thickness variation of the electric potential (including the detail of the BL at the anode), the BLs of counterion concentration, the BLs of solvent concentration, and the BLs of solvent pressure, respectively. Let us remind that, because of the dilute solution approximation, the solvent pressure, which is denoted as p s , is, in the LB21 theory, the  Lagrange multiplier analogous to the osmotic pressure π in the extended theory developed in this investigation.
In this actuation problem, t = 0.1 s is close to the instant in which the maximum deflection is attained. Then, back-relaxation occurs and steady state is reached at about t = 50 s. This can be seen in Fig. 5, which, specifically, compares the LB21 numerical results (solid line) with those obtained in this investigation by using the models described in Sect. 5, i.e., the FE model constituted only by standard isoparametric FEs (represented with filled squares in the graphs) and the two models enriched with GFEs, with either eight or twelve rows of such elements to discretise the BLs (represented with either filled triangles or void circles in the graphs, respectively). Figure 5 shows that our richest model gives the results closest to those of [54], thus validating the results of [54] as well. Here and henceforth, when representing the IPMC deflection, the symbol u y is adopted for the time-dependent maximum transversal displacement of the IPMC mid-axis, u y (X = L, Y = 0, t).
In Figs. 6, 7, 8, and 9 and throughout the paper, the through-the-thickness variation of the relevant fields is always computed in the IPMC cross-section at X = 2H ≡ L/100, as in [54]. 7 Overall, this comparison demonstrates that the FE model of [54] is accurate enough when it converges (at least, for the low applied voltage of 0.25 V) and that, for what concerns the present numerical models, the GFE extension is needed in order to predict the correct results. This is quite clear from the BLs of counterion concentration, solvent concentration, and solvent pressure at the cathode (see the left graphs in Figs. 7, 8, 9, respectively), mostly by looking at the peak values at Y/H = −1, which are largely underestimated in the absence of GFEs.
The difficulty in accurately describing the BLs where counterions pile up should be due to the fact that the LB21 theory disregards steric effects that would hamper a too large concentration of counterions at the cathode [19,28,71]. In Sect. 7, we will show that the extended theory proposed in this investigation addresses this issue by accounting for the counterion molar volume v, which, most of all, allows a more appropriate description of the IPMC response, thus overcoming some physical problems of the LB21 theory that we have faced in actuation at large applied voltages.
As reported in the left graphs of Figs. 8 and 9, the larger discrepancy between our GFEenriched numerical models and that of [54] concerns the solvent concentration at t = 10 s and the pressure at t = 50 s, the former being underestimated and the latter being overestimated by the FE model in [54]. Note that we assume as reference our most refined numerical model, given that it can provide results for large applied voltage, where the LB21 implementation in COMSOL Multiphysics ® fails, and it has been tested to be close to convergence through suitable mesh refinements.

Short-Circuit Sensing
In short-circuit sensing, in [54] a uniformly distributed nominal load of magnitude T 0 = 0.05 kPa could be imposed, linearly increasing with time up to t = 0.1 s and then maintained until steady state. Here, to be consistent with the COMSOL simulations in [54], we maintain the same load application modality, although our algorithm would also allow an instantaneous application of the load. Figures 10 and 11 compare the results obtained from our numerical models with those in [54] in terms of the time-evolutions of the deflection u y , the charge accumulated at the electrodes, defined as For these integrated quantities all the numerical models provide similar results, although the LB21 model is softer (Fig. 10) and our plain FE model underestimates the accumulated charge (right graph in Fig. 11). Let us also observe, from Fig. 11, that, before going to zero at steady state, the current I becomes negative because of the counter-diffusion of mobile ions, in qualitative agreement with experimental observations [35]. Figures 12, 13, 14, and 15 instead focus, at the instants t = 0.1 s, 0.5 s, 5. s, 50. s, on the through-the-thickness variation of the electric potential, on the BLs of counterion concentration, on the BLs of solvent concentration, and on the through-the-thickness variation of the solvent pressure.  While the agreement on the time-evolutions of deflection, electric current, and accumulated charge is remarkable, we observe a discrepancy between our results and those of [54] in the through-the-thickness fields at steady state. This is clearly shown in Fig. 12 for the electric potential, which displays a special feature, which we have carefully checked by further post-processing of our results and those of [54]: the prediction of its slope (i.e., the electric field) outside the BLs is identical in our models and that of [54], although in the bulk membrane our profiles of ψ are significantly shifted from that of [54]. This is due to different predictions of the BLs of counterion concentration at steady state (see Fig. 13). Interestingly, this discrepancy is substantially absent in the comparison of Fig. 14 on the BLs of solvent concentration. Moreover, Fig. 15, on the through-the-thickness variation of the solvent pressure, shows a discrepancy in the BLs in the initial part of the sensing response, at t = 0.1 s; such a discrepancy emerges only with respect to the more accurate GFE models, such that we infer that the poorer discretisation adopting only standard FEs leads to an inaccuracy. Let us now provide more details on the observed discrepancy. Our predictions of ψ in the bulk membrane start to drift away from that of [54] much after the peak in the electrical response. In our analyses, at steady state, ψ changes sign in the bulk closer to the IPMC axis (Y ≈ 0), while the computational model of [54] predicts that ψ changes sign at Y ≈ −H/2. In investigating on this discrepancy, we have discovered that these results are very sensitive to the integration time step. Also, the fact that ψ at steady state is not anti-symmetric (i.e., ψ does not change sign in Y = 0) is due to the complex interplay between various nonlinearities. A way to obtain an anti-symmetric ψ at steady state is to specialise the theory to the framework of small strains and rotations. As a minor indication, the lack of antisymmetry of ψ is largely smoothed out by setting the membrane Poisson ratio ν = 0, such as λ = 0, which eliminates the nonlinear contribution in the membrane elasticity differently weighing increase and decrease of volume (see Eq. (29)). This plays a role in determining the IPMC response because the BLs experience opposite volumetric deformations that may be relatively large, due to the localisation of the through-the-thickness direct strain component [15]. In fact, the nonlinear volumetric deformation in the BLs certainly influences the profile of the solvent pressure, in turn affecting the counterion transport through both the cross-diffusion and the convective flux of counterions with the solvent, which are features of the LB21 theory.

Actuation and Short-Circuit Sensing: Results for Large External Actions and Limits of the LB21 Model
In this section we start by demonstrating that the GFE-enriched implementation proposed in this work allows one to obtain results for much larger applied voltages in actuation (Sect. 7.1) and much larger applied loads in short-circuit sensing (Sect. 7.2) within the LB21 theory. In these two sections the adopted algorithm is still that reported in Appendix C and the numerical model is the "fine" one described in Sect. 5. Then, in Sect. 7.3, by focusing on actuation, we demonstrate that the novel extended theory presented in Sect. 2 overcomes some limits of the LB21 theory.

Actuation: Behaviour Predicted by the LB21 Theory at Large Applied Voltage
While the COMSOL implementation in [54] allowed a maximum applied voltage drop across the electrodes of about 0.25 V only, in this section we show that the proposed GFE  model can accurately predict the IPMC actuation response up to ψ 0 =2 V, which is in line with experiments and applications. In particular, Fig. 16 reports the maximum displacement u y as a function of time. We observe two drawbacks of the response predicted by the LB21 theory: for large ψ 0 the backrelaxation becomes anomalously large, such as the steady state displacement has magnitude much larger than that of the initial peak (having opposite sign), and, as shown by the inset in Fig. 16, the initial displacement peak is almost independent of ψ 0 .
As shown in Table 2, the instant t p in which the initial peak response is attained strongly depends on ψ 0 . This has to be accounted for in interpreting the results presented in the following of this subsection when in the same graph different curves related to different ψ 0 are all referred to t p . Instead, the instant t ss in which steady state occurs is conventionally assumed to be equal to 30 s independent of ψ 0 .
The drawbacks emerging from inspection of Fig. 16 can be explained on the basis of the BLs of counterion concentration that are reported in Fig. 17. Even if increasing ψ 0 , the BL at the cathode, where counterions pile up, maintains almost the same BL thickness because in the LB21 theory there is no limit on the counterion peak value at the interface between membrane and cathode (see left graph in Fig. 17). On the contrary, counterion depletion at the anode is limited by C = 0, such that this BL must enlarge to compensate for large quantities of counterions migrating towards the cathode. This asymmetry in the distribution of the counterions at the two BLs leads to asymmetric electric fields at the BLs, which,  however, reach the same maximum value at the Y = H and Y = −H because the quantity of counterions leaving the anode is about the same as that accumulating at the cathode (note that in the bulk membrane C ≈ C 0 and the sole significant derivative in the Gauss law (10) is that with respect to Y ). The Maxwell stress profile mirrors this behaviour of the electric field, such that the much larger thickness of the BL at the anode is responsible for a polarization moment (to be balanced by the mechanical moment) that contributes to back-relaxation [70]. In the actuation predicted by the LB21 theory at large ψ 0 this back-relaxation effect becomes significant quite soon, hampering large values of the initial displacement peak, and becomes questionably dominant at steady state.
At large ψ 0 , this difference in size of the BLs is so large that it leads to the anomalous electric potential profile of Fig. 18, which shows that for ψ 0 =2 V the voltage drop at the cathode BL is almost negligible in comparison with that the anode BL (experiencing counterion depletion), where, in fact, the same electric field is integrated over a much larger size to establish the voltage drop. We remark two aspects of this result: first, it makes the numerical convergence to an accurate steady state response a very challenging task; second,  most likely it is an unreliable response with respect to the actual IPMC behaviour, where the counterion accumulation at the cathode must be somehow bounded. Figure 19 reports the BLs of solvent concentration, which follow those of the counterions although the peak at the cathode may not increase too much because of the constraint (80) between the solvent volume fraction v s C s and the volume ratio.
Finally, Fig. 20 illustrates the BLs of solvent pressure, showing that the maximum negative value of p s , reached at the interface between the membrane and the anode (Y = H ), decreases with time and is independent of ψ 0 at stead state. This is due to the combined effects of C ≈ 0 therein and of the solvent counter-diffusion.

Short-Circuit Sensing: Behaviour Predicted by the LB21 Theory at Large Applied Loads
While in [54] the maximum applied load for which the FE model converges is of about T 0 = 0.05 kPa, here we can apply far larger loads. In Fig. 21, we report the time-evolution of the   Fig. 21 are made nondimensional by scaling the displacement through that predicted for the smallest load (T 0 = 0.05 kPa) by the linear elastic Euler-Bernoulli theory, which is denoted as u EB y . From Fig. 21 it is clear that increasing the load leads to a less-than-linear relationship between u y and T 0 . This behaviour is ascribed to the chosen hyperelastic strain energy for the membrane (29), which provides a stiffening response for strains far beyond the linear regime, such large strains occurring at the BLs [15,54]. Figure 22 reports the electric current I (left) and the accumulated charge Q (right), defined as in Eqs. (78) and (77), respectively. In order to understand these two graphs in Fig. 22, it is important to observe that their abscissas span different time scales. Of course, both I and Q increase with the applied load T 0 . When counterion counter-diffusion begins, at t ≈ 0.5 s independent of T 0 , I becomes negative and Q starts diminishing. The amplitude of this negative value of I , again, increases with T 0 . Also the steady state value of Q increases with T 0 .  Figures 23 and 24 report the time-evolution of the electric potential and the counterion concentration at the BLs, respectively. These figures show that the electrochemistry remains linear even for large applied loads in sensing. This means that, as well known in the IPMC literature, the IPMC sensing performance is quite poor, in contrast with actuation, where IPMCs can undergo large deformations under relatively small applied voltages. Figures 25 and 26, finally, illustrate the time-evolution of the solvent concentration at the BLs and the solvent pressure, respectively. The BLs of C s increase with time, contrary to the gradient of C s , corresponding to a decrease of the magnitude of p s towards zero. Also, the amplitude of the BLs of C s increases with the applied load T 0 .
To briefly summarise, Sects. 7.1 and 7.2 illustrate the superior capabilities of the GFE numerical model proposed in this investigation, along with some issues of the LB21 theory that are addressed, as shown next, by the extended electrochemo-poromechanical theory presented in Sect. 2.

Comparison Between the LB21 Electrochemo-poromechanical Theory and Its New Extension
From now on, we focus on actuation only. The constraint (5) drastically changes the behaviour of migrating species at the boundary layers with respect to that predicted by the LB21 theory, because the latter neglects the contribution of the counterion motion to the volume ratio. With the extended theory, on the one hand, the solvent must leave room for the counterions at the cathode. On the other hand, the constraint (5) hampers the possibility of counterions to pile up in a very large peak at the cathode, thus leading therein to a boundary layer occupying a larger membrane region, which is characterised by a smaller gradient of counterion concentration. As a comment on the IPMC models in the literature, the constraint (5) includes the counterion "steric effect", which has been differently introduced by Cha and Porfiri [28] through a sophisticated free energy of mixing directly delivering the nominal osmotic stress by differentiation with respect to the unconstrained deformation gradient.
In this section we illustrate the significant improvement in the predictions of the IPMC response due to the extension of the LB21 theory presented in Sect. 2. Henceforth, all the results are obtained by using the new GFE implementation, specifically adopting the "coarse" model described in Sect. 5. For the LB21 theory the GFE algorithm is that reported in Appendix C, while the numerical algorithm for the extended theory developed in this investigation is that provided in detail is Sects. 3 and 4. In this section the IPMC geometrical parameters are still the same as those reported in Sect. 6.1.

Material and GFE Parameters
Here and henceforth, we adopt the material parameters as identified in next Sect. 9 to model the experimental data of He et al. [42] through the extended theory of Sect. 2. These parameters are reported in Table 3 for the membrane only, given that the Lamé constants of the electrodes are as in Table 1. We note that the parameters that change with respect to those in [54] are D, v s , and C 0 . Clearly, the extended theory additionally includes the hydrated counterion molar volume v. The values of these parameters will be discussed in Sect. 9. With the parameters of Table 3, the Debye screening length (1) results D ≈3.66×10 −4 mm. Therefore, given that ε = 0.001 mm > D , the GFE parameters of Table 1 would seem to be still appropriate for the present analyses. However, in the following, we will show the benefit to increase ε and decrease α at the anode BL in order to capture its nonlinear behaviour under large applied voltage ψ 0 . Figure 27 displays the time-evolution of the maximum displacement u y for varying ψ 0 . The extended theory suffers much less back-relaxation and provides significantly larger influence of ψ 0 on the initial peaks of u y , although this relation is still less than linear. Let us remark that the recent experimental results of Arnold et al. [1] show that back-relaxation may lead to a final steady state deflection of opposite sign with respect to initial peak deflection, such phenomenon being still predicted by the new theory at large ψ 0 . All the benefits of the extended theory come from the larger size of the counterion BL it predicts at the cathode, due to accounting, through the constraint (5), for the relation between the size of the hydrated counterions and the volume ratio.  Consistently with this picture and the discussion of the results of Sect. 7.1, note that if in the LB21 model C 0 is increased, as for instance in Table 1, the initial peaks of displacement turn out to be even closer to each other.
The huge difference, between the two theories, in the counterion BLs is shown in Fig. 28, where the dashed curves represent the results of the extended theory. Note that these curves, at the anode and at steady state, are not as smooth as one would expect. This is a numerical issue that can be fixed by enlarging the size discretised with GFEs at the anode, as done to obtain the results of Fig. 29.
The results of Fig. 29 have been obtained by setting the thickness ε of the anode region discretised with GFEs to 0.004 mm, instead of 0.001 mm, the latter still being used for the cathode. Moreover, in these asymmetric meshes, at the anode, the GFE parameter α is diminished, with respect to that of Table 1, to α = 250.   Figures 30 and 31 represent, by using the symmetric and asymmetric meshes, respectively, the BLs of solvent concentration. Also these profiles highlight substantial differences in the predictions between the two theories considered. A common feature of the two theories is that C s at steady state diminishes at the cathode and increases at the anode, with respect to the levels reached when the initial displacement peak is attained. This is due to the solvent counter-diffusion. However, in the extended theory C s is very much limited at the cathode, with respect to the LB21 theory, because of the volume occupied by the counterions piling up therein. This results in C s even attaining a minimum at Y = −H , being much lower than its initial value C 0 s in that region. Under the constraint of conservation of the total solvent mass imposed by the first relation in (13) and its analogous condition at the IPMC ends X = 0 and X = L, this behaviour at the cathode implies that C s must be larger than C 0 s somewhere else. Given that in the bulk membrane C s ≈ C 0 s , as illustrated in Fig. 32, it turns out that C s > C 0 s at the anode. This indicates that, with the chosen material parameters, in the extended theory the solvent flux is more influenced by the need of hydrated counterions to occupy volume in accumulating at the cathode (thus pushing away the free solvent) than by transport through electro-osmosis, which is provided by the dependence of the solvent flux on the electric field (see Eq. (48)). Figure 33 finally shows the remarkable difference between the two theories in the predicted profiles of the electric potential ψ . The most important quantitative discrepancy is in the voltage drop at the anode at steady state, which is much larger in the LB21 theory. The The results in Fig. 33 also suggests that, in order to facilitate the IPMC characterisation, it would be very useful if in actuation experiments the voltage in the membrane mid-plane were measured, up to steady state.

Influence of the Permittivity on the Actuation Response
With this section, by applying the theory proposed in this work, we aim at contributing to the recent discussion about the role of the permittivity ε on the actuation response of IPMCs [20,60]. Boldini and Porfiri [20] argue that a compensation between the thickness of the BLs, which is governed by Eq. (1), and the constitutive law for the Maxwell stress (3) makes the magnitude of the Maxwell stress independent of ε. Moreover, given that increasing ε leads to larger BLs with the same Maxwell stress, for a given applied voltage ψ 0 , the bending displacement scales with √ ε, i.e. it is proportional to the Debye screening length D . However, in order to precisely observe this scaling, one should focus on the steady state behaviour in the absence of nonlinearities, such as one should, first of all, apply an appropriately low voltage ψ 0 . The same discussed scaling would hold for the whole history if, additionally, the species fluxes were governed by linear equations. In order to meet this requirement as close as possible, in our extended theory, we need to set both γ = 0 in Eq. (47), in such a way as to resort to the simplest diagonal mobility matrix M , and very low molar volumes v and v s .
In this section, we initially apply ψ 0 = 1 mV, we test for M both the diagonal form and the form proposed in [54] (γ = 1 in Eq. (47)), and we consider three values of permittivity, namely ε = 10 −10 , 10 −12 , and 10 −13 C/(mV mm), along with three values of the hydrated counterion molar volume v, namely v = 500,000, 20,000, and 500 mm 3 /mol. The GFE parameter α and the Lamé constants for the electrodes are as those listed in Table 1, while the other material parameters for the membrane are as those in Table 3. Consistently with Eq. (1) and the scaling law we want to investigate, the regions discretised by GFEs have size ε dependent on ε. Specifically, for ε = 10 −10 C/(mV mm) that size is 10 −3 mm (as in the foregoing analyses), whereas for ε = 10 −12 C/(mV mm) and ε = 10 −13 C/(mV mm) ε is set to 10 −4 mm and 3.16228×10 −5 mm, respectively. The IPMC geometry is that provided in Sect. 6.1.
The results are illustrated in Figs. 34 and 35, reporting the time-evolution of the maximum transversal displacement, that is normalised by D , which, in turn, depends on √ ε. Therefore, when the curves coincide, the scaling discussed in [20] holds. Note that also the non-dimensional time in the abscissas, defined as in Eq. (46), scales with √ ε; specifically, t ∝ 1/ √ ε. First, we note that independent of the mobility matrix and of the value of v, all the steady state responses scale with √ ε, as established by Boldini and Porfiri [20]. This is because the fluxes nonlinearities do not affect the steady state, which is governed only by the Poisson equation (that is, the Gauss law (10) combined with relation (33) between the electric displacement and the electric potential) and the Navier-like equation (that is, the mechanical equilibrium (7) written in terms of the displacement field through the constitutive and kinematic laws), the latter being substantially linear for small ψ 0 . This steady state response can certainly be computed analytically by first solving the electrochemistry uncoupled from mechanics (Porfiri [67]) and, then, using the electrochemical "active stresses" (or eigenstresses) as known external actions to solve the mechanical problem and obtain the IPMC deformation (see, e.g., Boldini et al. [15] and references therein).
Second, as expected, the sole case in which the displacement is close to be proportional to √ ε all along the time history is that in which γ = 0 with small v (right graph in Fig. 34). We

Fig. 36
Actuation under large applied voltage ψ 0 : modelling capability of the extended formulation with mobility matrix proposed in [54] (γ = 1) for variable permittivity and counterion volume v = 500,000 mm 3 /mol observe that this scaling is more accurate for v = 20,000 mm 3 /mol (blue curves in Fig. 34) than for 500 mm 3 /mol (green curves in Fig. 34), thus suggesting that the effect of some nonlinearity is still present (most likely, related to the constraint (5) on the volume ratio). Third, as shown in Fig. 35, when M as in [54] is adopted, the influence of v only slightly changes the initial peak values of the displacement, whereas they strongly depend on the permittivity. With the adopted values of ε, increasing it of two orders of magnitude leads to an initial displacement peak only about five times larger.
Also, the comparison between Figs. 34 and 35 clearly shows that changing M has an enormous impact on the initial displacement peaks. To the best of our knowledge, this important result has never been pointed out in the IPMC literature, in which, so far, little attention has been paid on the interaction between the counterion and solvent transports. In the specific case here of concern, M as in [54] allows peaks more than ten times larger than those obtained with the simple diagonal form of M (γ = 0). This different behaviour should be ascribed to the electro-osmosis, which is accounted for in [54], where the solvent transport depends also on the electric field. In order to eliminate this contribution, one should set a very small value of the solvent molar volume v s .
The present analyses show that the proposed GFE implementation of the IPMC multiphysical response allows one to obtain results also for small values of ε. In the benchmark problem of this section, further decreasing ε below 10 −13 C/(mV mm) hampers the convergence of our code because of floating point errors due to the fact that the nodal coordinates of the GFEs turn out to be too much close to each other. It should be possible to address this issue by adopting also for the displacement nodal DOFs an appropriate non-dimensional variable.
As a final remark on the numerical efficiency of the proposed GFE algorithms, Fig. 36 displays results demonstrating that our code does converge in the case of large applied voltage and relatively small ε.

On the Capability of the Extended Theory to Predict Experimental Results
Although experimental results on the IPMCs response are still affected by very large scatter [1], in this section we focus on part of the data provided by He et al. [42] to give a hint on the predictive capability of the electrochemo-poromechanical model developed in this investigation.
Note that He et al. [42] characterise the mechanical properties of the employed Nafion™ membranes by resorting to Berkovich nano-indentation, which results in measuring an elastic modulus that turns out to increase with the membrane thickness, from ≈ 0.2 GPa for 2H = 0.22 mm to ≈ 0.8 GPa for 2H = 0.8 mm. He et al. [42] ascribe this behaviour to the larger annealing time required in the manufacturing of thicker Nafion™ membranes, possibly leading to larger degree of polymer cross-linking and crystallinity.
For the purposes of this investigation, we have followed a very simple calibration of the material parameters, resulting in the constants reported in Table 3 for the membrane (and in Table 1 for the electrode and GFE constants). First, we have neglected the dependence of the membrane moduli on the membrane thickness. Second, we have measured the maximum displacement under the application of a single voltage cycle. Third, we have calibrated the material parameters on a single datum, that is the maximum displacement obtained for 2H = 0.42 mm at peak applied voltage equal to 2 V; in doing so, we have taken the liberty to decrease C 0 as needed, while increasing v such that C 0 v is suitably smaller than 1. It turns out that C 0 v is about five times smaller than that expected for the Nafion™ and, correspondingly, v is about five times larger than the molar volume of hydrated counterions. This is surely a shortcoming, not only of the present calibration procedure, but also of the adopted theory, which should be further improved as explained in the concluding remarks of Sect. 10.1. We feel however that a rigorous identification procedure would be worth in the presence of a much more precise experimental characterisation of IPMCs.
However, on the basis of the results obtained in this investigation and the foregoing discussion, we infer that there is room for improvement in the comparison discussed next, which is reported in Fig. 37.
The match is quite satisfactory for ψ 0 = 2 and 2.5 V. The discrepancy increases for ψ 0 = 3 V, where, by the way, we expect that the influence of phenomena unaccounted for by the present theory becomes more important. Among them, we mention the water electrolysis, the roughness of the membrane-electrode interfaces, the spatial heterogeneity of the actual membrane permittivity and diffusivities, and the mechanical behaviour beyond hyperelasticity of the membrane, which is expected to involve some relaxation. Our numerical model fails to converge for ψ 0 = 3.5 V, for which we should further refine the GFE discretisation.
Finally, we note that, in order to identify the material parameters, it would be useful that the frequency of the applied sinusoidal voltage scaled with the membrane thickness, as suggested by Eq. (46), which is not the case in the experimental protocol followed in [42].

Concluding Remarks
This work has focused on two advances in modelling actuation and sensing of ionic polymer metal composites (IPMCs). First of all, we have presented in detail generalized finite  [42] element (GFE) algorithms that can accurately compute the IPMC response up to load levels consistent with experiments, which has been so far considered an extremely challenging task because of the highly nonlinear IPMC behaviour. This difficulty is due to very thin boundary layers (BLs) of concentration of mobile ions (counterions) occurring at the membraneelectrode interfaces in both actuation and short-circuit sensing. The availability of the GFE codes has allowed us to discover some shortcomings of the electrochemo-poromechanical theory developed by Leronni and Bardella [54]. Therefore, as a second contribution, we have proposed a theory that slightly extends that in [54], providing evidence and the physical reasons for this need.
The novel theory extends that in [54] since it removes the hypothesis of dilute solution of counterions and solvent, which has an impact on two governing equations. One is the definition of the contribution to the Helmholtz free energy due to the mixing of counterions and solvent. The other, much more important, is the constraint relating the volume ratio of the membrane with the species motion. In particular, here we account also for the counterion motion in that constraint, thus hampering a too large peak of counterion concentration at the electrode-membrane interface where counterions pile up in actuation and short-circuit sensing. This, for a given voltage-driven amount of counterions leaving the anode region to accumulate at the cathode region of the membrane in actuation, avoids that the thicknesses of the two BLs grow too much asymmetrically, which would lead to a questionably large back-relaxation [54,70].
The novel computational models adopt GFEs to enrich the discretisation of the membrane regions adjacent to the electrodes where BLs occur. This allows us to capture, therein, the very large gradients of counterion concentration, solvent concentration, and electric potential, along with localisation of the through-the-thickness direct strain component [15]. This technique constitutes an effective way to address the most important computational issue on the IPMC modelling, which has required, so far, to artificially play with model parameters such as the membrane permittivity [20,60], or to limit the analysis to too small applied voltages/loads. The proposed implementation strategy along with the comparison between different theoretical models have allowed us to demonstrate, for the first time to the best of our knowledge, the key role of the specific constitutive prescriptions affecting the equations governing the fluxes of counterions and solvent. They in fact have an enormous impact on the model prediction of the initial peak deflection in actuation. This has further consolidated the knowledge on the crucial importance of being able to precisely model the capability of counterions and solvent to either leave or accumulate at the boundary layers.
Overall, the proposed model gives the basis for an effective and reliable tool for the analysis and design of IPMCs. However, it surely requires further efforts, as detailed next.

Open Issues
As very recently demonstrated by the thorough experimental campaign of Arnold et al. [1] on Nafion™-Pt IPMC actuators under controlled environmental conditions, there still is at least an important aspect to include in the IPMC modelling, that is the actual humidity of the membrane. The membrane is in fact assumed to be completely saturated in the present investigation and in most theories for IPMCs. However, even if desirable, this may be not the case and the results in [1], obtained by carefully maintaining different levels of membrane humidity during the tests, show that the actual solvent content significantly impacts the IPMC behaviour. This, of course, agrees with the important role of the solvent as highlighted by our present modelling.
Also, Arnold et al. [1] have shown that the scattering in the experimental results on IPMCs, as nowadays produced, can be extremely large, thus challenging the reliability of the material parameters identification of any theory. Novel additive manufacturing techniques might help in improving the consistency of the IPMC performance, towards large-scale industrial applications [24,43]. Other experimentalists (see, e.g. [83]) are instead investigating on the possibility of stabilising the IPMC performance at different relative humidities by tuning the concentration of ionic liquids.
A drawback of the present model is that it seems to require to set a too low initial value of the counterion concentration C 0 and a too large value of the hydrated counterion molar volume v in order to predict, in actuation, the dependence on the applied voltage of the initial peak deflection u Y (t p ). In the light of the present results on the key role of the fluxes in establishing u Y (t p ), we infer that there is room for improvement of this modelling capability. One option is to leverage on the mobility matrix governing the solvent and counterion fluxes (Vanag and Epstein [84], Kondepudi and Prigogine [50]), whose modification with respect to that proposed in [54] might be consistent with the removal of the hypothesis of dilute solution of solvent and counterions.
Additionally, another possibility would be the enrichment of the free energy of mixing (30), where one could also account for the enthalpic contribution governed by the nondimensional Flory-Huggins interaction parameter χ [32,38,47], thus extending Eq. (30) to About the numerics, we underline that we have employed a single enrichment function f (Y ) for all the fields experiencing large gradients at the BLs. It would be worth to explore different possibilities in order to further improve the computational efficiency. Another option might consist in the use of suitable solver preconditioners [39], whose investigation is hampered by the use of a commercial FE code such as ABAQUS.
The IPMC modelling should also account for the interface roughness (Porfiri [68]) and for the modified electrochemo-poromechanical properties of the membrane regions next to the electrodes, due to a couple of factors. The first one is process-dependent and is related to the counterion steric effect where counterions pile up (Borukhov et al. [21], Kilic et al. [49]). The second one is concerned with "composite layers" of much larger permittivity and much lower diffusivity in the membrane where metal particles can fill ionomer interstices (Tiwari and Kim [82], Aureli and Porfiri [5], Cha and Porfiri [27]). Changing material properties in membrane layers would be straightforward with the present FE model if data on this were available. To this purpose, a thorough and vast experimental campaign, along with an appropriate identification procedure, would be needed.
The overall behaviour of the BLs could be significantly modified also by a constitutive law for the membrane accounting for the behaviour beyond the elastic regime, especially regarding relaxation ensuing from viscoelastoplasticity (e.g., Silberstein and Boyce [77], Silberstein et al. [78]). Also the membrane small-scale behaviour should be investigated, being possibly important in the BLs (Feng et al. [36]).
The non-ideal behaviours in the species interactions in the membrane, namely long-range electrostatic interactions, short-range physical interactions, and steric effects due to finite pore size, as discussed and modelled by Boldini and Porfiri [19], would deserve to be further studied at experimental boundary conditions, through numerical models as that proposed here.
Finally, an important test which would require little effort to be simulated with the present 2D GFE model is that of a propped-cantilever IPMC subjected to a voltage drop across the electrodes, where experimentalists measure the "blocking force", that is the reaction force developed by the support constraining the otherwise free IPMC end.

Appendix A: FE Isoparametric Implementation for the Electrodes
The electrodes are modelled using standard isoparametric elements, under the assumptions or arbitrary large displacements and deformations. The internal virtual energy rate for the balance equations in the electrodes is where E is the domain occupied by the electrodes in the reference configuration. Upon time integration over the time step t , the discretised internal energy increment reads The deformation gradient and the right Cauchy-Green deformation tensor are as in Eqs. (53) and (55), whereas the vectorial form of the Green-Lagrange strain tensor is The second Piola-Kirchhoff stress S can be computed by differentiating Eq. (11)  and I is the fourth order identity tensor, coincident with a 4x4 identity matrix in array notation.

Appendix B: Consistent "Stiffness Matrix" for the Membrane
The submatrices constituting the FE "stiffness matrix" K in (70) have dimension N × N , except when otherwise indicated below. They are obtained in such as way as to be consistent

B.1 Common Derivatives to Compute the Stiffness Matrix
where I is the 4x4 identity matrix and which, in contrast with (5), neglects the counterion motion. Moreover, for the free energy of mixing, [54] adopts the approximation of (30) obtained by assuming C/C s 1. Therefore, in this formulation, one has C s = (J + v s C 0 s − 1)/v s and, from the definition ofμ, it results r = C/C s = exp(μ). Moreover, the nodal DOFs are u,ψ ,μ, andp s .
Equation (32) still holds for the total nominal stress in the membrane P, but the osmotic contribution P osm now reads The internal pseudo-forces read where D = D − D n , μ =μ −μ n , and J = J − J n . The consistent "stiffness matrix" K, of dimension 5N × 5N , required to complete the FE algorithm and ensure second order convergence of the global Newton-Raphson scheme is defined as the differentiation of the right-hand sides of Eqs.

Competing interests
The authors declare no competing interests.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.