Electrochemo-poromechanics of ionic polymer metal composites: identification of the model parameters

We propose a procedure to identify the parameters of a model for the multiphysics response of ionic polymer metal composites (IPMCs). Aiming at computational efficiency and accuracy, the procedure combines analytical structural mechanics and fully-coupled electrochemo-poromechanics, additionally resorting to an evolutionary algorithm. Specifically, we consider the finite-deformation electrochemo-poromechanical theory recently developed by our group, which couples the linear momentum balance, the mass balances of solvent and mobile ions, and the Gauss law. Remarkably, the theory constitutively accounts for the cross-diffusion of solvent and mobile ions. This, in conjunction with a generalised finite element implementation that we have recently proposed, allows us to accurately capture the boundary layers of ions and solvent concentrations occurring at the membrane–electrode interfaces, which govern the IPMC behaviour in actuation and short-circuit sensing. Thus, we can explore the IPMC behaviour under external actions consistent with applications and obtain accurate predictions with a reasonable computational cost for wide ranges of model parameters. We focus on experimental data from the literature that are concerned with Nafion™-Pt IPMCs of variable membrane thickness and subjected to peak voltage drop across the electrodes ranging from 2 to 3.5 V (under alternating current). Importantly, the considered tests deal with both the tip displacement of cantilever IPMCs and the blocking force of propped-cantilever IPMCs. Overall, the adopted theory and the proposed procedure allow unprecedented agreement between predictions and experimental data, thus marking a step forward in the IPMC characterisation.


Introduction
Ionic polymer metal composites (IPMCs) are functional materials that leverage on electrical stimuli to work as actuators, where they exhibit potential in a wide range of applications [1][2][3][4].Specifically, IPMCs feature a sandwich structure with the purpose of improving the transduction capability of the ionic electroactive polymer (EAP) constituting their soft core, in most cases made by Nafion™.In fact, plating such an ionomeric membrane between metal electrodes allows overcoming the major drawback of soft and homogeneous EAPs, that is the need of high applied voltage to display appropriately large displacement in actuation.In particular, the membrane is characterised by backbone chains with fixed ions of equal charge, usually anions, and is made globally electroneutral by soaking it in a solution including counterions that may move within the membrane.The benefit achieved with the IPMC sandwich structure ensues from the formation, under both actuation and short-circuit sensing, of boundary layers (BLs) of counterion concentration in tiny membrane regions at the interfaces with the electrodes.
Under ideal conditions [5,6], the BL thickness is comparable to the Debye screening length ℓ D , which is the size of the electric double layers of charge forming at the interfaces between membrane and electrodes, the latter being substantially impermeable to counterions [7].The value of ℓ D depends on the absolute permittivity of the membrane ε, the absolute temperature θ, and the initial nominal molar concentration of counterions C 0 , reading [5,8] in which F ≈ 96 485 C mol −1 and R ≈ 8.3145 J (mol K) −1 are the Faraday and gas constants, respectively.Note that C 0 is assumed to be spatially uniform and equal to the concentration of fixed (an)ions in the membrane.The very large variation of counterion concentration occurring in the BLs adjusts most of the voltage drop across the IPMC, thus corresponding to large Maxwell and osmotic stresses in the BLs, where the mechanical stress required by equilibrium leads to a complex deformation field that governs the IPMC electrochemomechanics (Kilic et al [9], Cha and Porfiri [10], Boldini et al [11], Panteghini and Bardella [6]).One of the strongest nonlinearities is related to the tendency of the counterion concentration to vanish at the BL subject to ion depletion, leading both to very large magnitude of the chemical potential therein and to different size of the two BLs, this latter phenomenon being a source of back-relaxation in actuation [12,13].
Although IPMCs work both in actuation and sensing (Shahinpoor and Kim [14]), the performance of these dual behaviours is very much different because of the highly nonlinear IPMC electrochemistry, which is governed by a Poisson-Nernst-Planck system of equations significantly modified by the coupling with finite-deformation (poro)mechanics [10,15].In fact, IPMC sensing exhibits weak electrostatics, even under relatively large mechanical input.In this investigation, we propose a procedure for the identification of the IPMC model parameters that focuses on the IPMC most prominent feature in applications, that is their soft actuation, consisting in the capability of displaying large and fast deflection under relatively low applied voltage (≈2-3 V) and working under both AC and DC conditions (as documented since the pivotal contribution of Asaka et al [16] to the most recent experimental efforts of Arnold et al [17,18]).
We refer to the electrochemo-poromechanical theory developed by Leronni and Bardella [15] and Panteghini and Bardella [6], which, for the sake of brevity, henceforth, is referred to as 'LBP theory'.The essential feature of this theory, with respect to other contributions [10,[19][20][21][22][23], is the thermodynamically-consistent modelling of the transport of both counterions and solvent, including their cross-diffusion.In particular, in [6], we have proposed a finite element (FE) implementation of the LBP theory adopting specific generalised FEs (GFEs) to discretise the BL regions.The accurate modelling of the BLs has been a crucial numerical issue to address, due to the fact that in usual IPMCs the ratio H/ℓ D , with H denoting the membrane semi-thickness, is extremely large.
The GFE method [24] allows us to accurately solve actuation and sensing problems under relatively large applied actions (consistent with applications and experiments) and, therefore, to demonstrate that the nonlinear interplay between the counterion and solvent fluxes strongly affects the IPMC dynamic response.Additionally, the GFE implementation permits us to obtain results for a wide range of model parameters, which is clearly crucial for the parameters identification.The LBP theory and its FE implementation are summarised in section 2.
With this computational model at hand, this investigation gives a contribution to the IPMC characterisation by showing that the use of an adequate multiphysics theory allows the identification of a single set of model parameters to reasonably predict several different experimental results on IPMCs constituted by the same materials.Ideally, the considered experiments should encompass both actuation and sensing.However, this is very much complicated by the large scatter of the IPMC performance depending on details of the IPMC manufacturing and environmental conditions that are often left unsaid in the literature.This issue has led us to look for IPMCs consistently produced and tested, eventually resulting in the use of a single paper reporting the largest possible number of different experimental results.
Therefore, here we propose a procedure to identify the parameters of the LBP theory by focusing on the two most common experiments in which the IPMC is subjected to a voltage drop across the electrodes, denoted as ψ 0 .The first experiment measures the displacement at the free end, here denoted as u Y (L), of a cantilever IPMC of length L.Here and henceforth, Y is the coordinate normal to the IPMC axis at rest, the latter being denoted as X.The second experiment is referred to as the 'blocking force' experiment in the literature.As depicted in figure 1, it refers to a propped-cantilever IPMC where the reaction force developed by the support, F Y , which is called the blocking force (BF), is measured.The BF experiment is an interesting benchmark because, contrary to the actuation of a cantilever IPMC, it features a non-uniform curvature.In principle, for not-too-slender IPMCs, this might even promote non-negligible fluxes along the IPMC axis because of the coupling between electrochemistry and poromechaincs.
In particular, in section 3, we propose a simple yet accurate model for the parameters identification.This model relies on decoupling the structural mechanics, which admits straightforward analytical solutions, thus enabling a strong reduction of the computational cost.This allows us to employ an evolutionary algorithm in the identification procedure.After the identification, the results are anyway checked and provided by running fully-coupled GFE analyses for the whole IPMC.
Although there are several experimental data in the literature on the actuation problems modelled here (for what concerns the BF, see, e.g. the review [25]), in this investigation (in section 4) we select the experimental data of He et al [26] to test the proposed identification procedure.In fact, He et al [26] tested Nafion™-Pt IPMCs of variable membrane thickness to measure both u Y (L) and F Y under sinusoidal voltage signal ψ 0 of peak magnitude ranging from 2.0 to 3.5 V.By assuming that all these IPMCs are manufactured such as they have the same electrodes and, most of all, about the same membraneelectrode interfaces, the work of He et al [26], differently from others, gives a large enough set of self-consistent data to test the proposed procedure.
On the basis of the obtained results, we infer that the proposed model is very much promising for the effective analysis and design of IPMCs, although several important issues remain to be addressed, as reported in the concluding remarks of section 5.

Summary of the LBP electrochemo-poromechanical theory for IPMCs and its FE implementation
The mechanical deformation of each macroscopic material point is expressed in terms of the deformation gradient where x is the current position vector, ∇ denotes the material gradient, such that (∇x) iJ = ∂x i /∂X J , in which X is the reference (material) position vector and, when index notation is adopted, small case and capital case subscripts indicate the spatial and material coordinates, respectively, by referring to a fixed orthonormal rectangual Cartesian system.The reference (initial) configuration is undeformed and electroneutral.Additionally, the fluid phase is assumed to always saturate the membrane.By neglecting the microscopic volumetric deformation of both the fluid phase and polymer chains, the macroscopic volume ratio J ≡ det F is subject to the kinematic constraint where C and C s are the nominal molar concentrations of counterions and free solvent, respectively, v is the molar volume of hydrated counterions (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.Here and henceforth, the adjective nominal refers to the reference total volume.
In the LBP theory four balance equations govern the membrane behaviour and read Div P = 0, (4a) Ċs + Div J s = 0, (4b) The laws in equation ( 4) are written in the reference configuration, such as they all involve the material divergence Div.In the overall linear momentum balance disregarding externally applied body forces, equation (4a), P is the nominal stress tensor, so that (Div P) i = ∂P iJ /∂X J .In the solvent and counterion mass balances, equations (4b) and (4c), J s and J are the nominal molar fluxes of free solvent and hydrated counterions, respectively, and ( ) indicates partial time derivative, that is, ( )(X, t) ≡ ∂( )(X, t)/∂t.In the Gauss law, equation (4d), D is the nominal electric displacement.For all the details on the assumptions behind balance equations ( 4), the reader is referred to [15].
An assumption that is crucial for the predicted behaviour is that the interfaces between membrane and electrodes are impermeable to both solvent and counterions.Under the further hypothesis that the electrodes are perfect conductors, they turn out to be subject to balance (4a) only.
The IPMC relevant geometrical parameters, as displayed in figure 1, are the membrane thickness, 2H, the electrode thickness, h ≪ H, and the IPMC length, L ≫ H. Contrary to the actuation and sensing problems on the cantilever IPMC modelled as a two-dimensional (2D) continuum, the IPMC width B also enters the BF benchmark, as F Y scales linearly with B. Usually, B ≫ H such that 2D models for IPMCs assume plane-strain conditions.
About the initial and boundary conditions, we begin with those concerned with the electrochemistry and holding both in sensing and actuation (i.e. for any IPMC), under the assumptions of the present model.Electrode impermeability implies vanishing fluxes at the membrane-electrode interfaces: which may be adequate for IPMCs operating in a highly humid environment [17,18].The IPMC slenderness allows us to disregard edge effects, such as at the IPMC ends, X = 0 and X = L, we conveniently impose both absence of charge accumulation, that is D X = 0, and J s X = J X = 0.The electrochemistry also requires initial conditions, which read where the latter relation denotes electroneutrality at rest.Actuation is characterised by the application of a possibly time-varying voltage drop across the electrodes ψ 0 (t) such that where ψ denotes the electric potential.By now referring to mechanical boundary conditions, those shared by the two actuation boundary value problems (BVPs) here considered require fully-clamped IPMC left-end side, such that, by denoting with u = x − X the displacement field, u (t) = 0 at X = 0, and traction-free sides of outwards normal ±Y, which, in our 2D setting, require We also always assume P xX (t) = 0 at X = L, whereby distinction between the two BVPs depend on the tangential mechanical boundary condition at the right IPMC end.In the cantilever IPMC where u Y (L) is measured one has Note that with this last condition we spread the support all along the right-end side, thus avoiding as much as possible numerical issues and stress concentrations that might trigger electrochemomechanics difficult to interpret.Given that the IPMC is a sandwich structure with soft core, we however note that the actual way the support is realised might have an impact on the mechanical behaviour [27].
For all the details about the derivation of the constitutive laws through the Coleman-Noll procedure the reader is referred to [6,15].Here, we just provide the rate of the nominal internal energy density in the membrane, which is also the starting point to write the weak form of the balance equations to be implemented in a FE code: (7) where • denotes the inner product, is the nominal electric field, µ s is the solvent chemical potential, µ is the counterion chemical potential, and is the counterion electrochemical potential.It is assumed that the last two contributions to U in (7) have dissipative nature, while P, E, µ s , and µ are obtained in terms of the Helmholtz free energy W as where π is the Lagrange multiplier needed to impose the constraint (3) [28] and W admits the additive decomposition W mec , W mix , and W pol being the contributions due to the membrane deformation, the mixing of solvent molecules and counterions, and the membrane polarisation, respectively.The membrane elasticity is governed by the isotropic compressible Neo-Hookean model proposed in [29]: where are the Lamé constants, here expressed in terms of the Young modulus E and the Poisson ratio ν, and C = F T F is the right Cauchy-Green deformation tensor.W mec in equation ( 12) displays two features that are not necessarily common in the literature when resorting to nonlinear elasticity in a complex theory as the present one.First, the deviatoric and volumetric behaviours are coupled, aiming at more reliable predictions of the membrane elastomeric behaviour when strains are large [28,30], which may be the case in the IPMC BLs [11,15].Second, the membrane elasticity depends on two elastic constants, aiming at an optimal control of the interplay between the flux of the fluid phase and the volumetric deformation, which is surely important in the LBP theory coupling poromechanics with counter-diffusion phenomena 1 .The free energy of mixing is assumed in the form which is consistent with the hypotheses that the fluid phase is an ideal solution of solvent and counterions [32,33] displaying purely entropic behaviour and all the membrane phases are incompressible [34].
The membrane is assumed to behave as an ideal dielectric, whose electrostatic energy density reads From the equations above, the constitutive laws for the recoverable fields follow.The total nominal stress in the membrane turns out to be the sum of mechanical, P mec , osmotic, P osm , and Maxwell, P pol , contributions: where ⊗ denotes the tensor product 2 .By combining equations (10b) and ( 8), we rewrite the nominal electric displacement as whereas combination of equations (10c) and ( 13) leads to the solvent chemical potential The counterion electrochemical potential follows from of equations (10d), (13), and ( 9), as As observed by Boldini and Porfiri [34], the use of the constraint (3) and the chemical potentials ( 17) and ( 18) allows one to account for the steric effects of the migrating species in a way similar to that resulting form approaches relying on statistical mechanics [9,36].

Non-dimensionalisation of the problem and resulting forms of the fluxes
An important step towards an efficient FE implementation consists in defining appropriate non-dimensional primal variables.Note that, because of the crucial role played by the BLs, this step is also crucial if one aims at analytically solving linearised electrochemomechanical problems for IPMCs (see, e.g.[11,[37][38][39] and references therein).
In [6] we define the non-dimensional chemical potential the non-dimensional osmotic pressure 2 Hence, by introducing the left Cauchy-Green deformation tensor b = FF T and the electric displacement in the current configuration d = J −1 FD [35], the total pressure p (positive if compressive) can be conveniently written as the trace of the Cauchy stress σ = PF T /J: the non-dimensional electric potential and the non-dimensional time where Rθ/F is the thermal voltage and D is the diffusivity of the hydrated counterions in the medium constituted by the polymer network and the free solvent molecules.This diffusivity enters the expressions for species fluxes, which are obtained through an appropriate definition of the mobility matrix [15,32,40] that allows for fulfilling the dissipation inequality J s • ∇µ s + J • ∇μ < 0 ∀(J s , J) ̸ = (0, 0) and for describing the counterion and solvent cross-diffusion [15,41].They read where D s is the diffusivity of the free solvent in the medium constituted by the polymer network and the hydrated counterions, the non-dimensional parameter γ can switch on and off the cross-diffusion and the counterion transport by convection, and Introducing the variable r reduces equation ( 19) to r = exp (μ ∆ ) and implies that, in our solution algorithm, we determine C s in terms of J and r, as Setting γ = 1 in equations ( 23) and ( 24) recovers the form proposed in [15], while γ = 0 leads to the simplest possible diagonal form of the mobility matrix, suppressing both cross-diffusion and counterion transport by convection.If one chooses γ = 1 it turns out that the solvent flux is independent of both ∇C and ∇C s , whereas electro-osmosis is included through the dependence on ∇ ψ [15].On the contrary, by setting γ = 0, J s depends on ∇μ ∆ and electro-osmosis is neglected.Moreover, inspection of equation ( 24) allows singling out the convective flux of counterions with the solvent, namely CJ s /C s , which gives an important contribution in sensing, where the application of a mechanical load triggers the solvent motion, which is accompanied with some counterions.Within its motion, the solvent establishes a volumetric deformation gradient leading to a further counterion Fickian diffusion, again included in equation ( 24) [15].Note that, independent of γ, the counterion flux is always a function of ∇μ ∆ , ∇π, and ∇ ψ.In particular, the classical Nernst-Planck law, written in the reference configuration, is recovered by setting γ = 1 in the limit of immobile solvent, that is D s = 0.In this case, both ∇C and ∇ψ contribute to the counterion migration, respectively governing the classical Fick and electrophoretic effects [7].
In [6] we have demonstrated that the value of γ has an enormous impact on the peak response and, in particular, accounting for cross-diffusion and convection leads to much better results than those obtained with a diagonal mobility matrix (i.e.γ = 0).Therefore, in this work all the calculations are performed by setting γ = 1.More details on the fluxes are provided in appendix.

Constitutive behaviour of the electrodes
The metal electrodes are known to undergo relatively small strains [6,11], albeit they may experience large displacements.This motivates the assumption that they deform within a purely elastic regime governed by the computationally-convenient De Saint-Venant-Kirchhoff strain-energy density, which is isotropic and quadratic in the Green-Lagrange strain tensor E = (C − I)/2, with I denoting the second-order identity tensor: Here, λ e and G e are the Lamé parameters of the electrodes, having Young modulus E e and Poisson ratio ν e .

FE implementation
The nodal degrees of freedom (DOFs) in the membrane are the displacement vector u along with the non-dimensional electric potential ψ, chemical potential μ∆ , and osmotic pressure π, as defined in equations ( 21), (19), and (20), respectively.All the details of the FE implementation are provided in [6].
Here, we just recall the crucial GFE discretisation of the membrane regions of thickness ℓ ε that are adjacent to the electrodes, where ℓ ε is a parameter of the numerical model to be suitably calibrated.The GFE extension enriches the standard polynomial shape functions and, as in [6], it relies on the use of the single enrichment function in which α is a non-dimensional parameter of the numerical model to be properly selected to capture the gradients of the relevant fields in the BLs.These fields are u Y , ψ, μ∆ , and π, since the standard isoparametric discretisation is appropriate for u X .Note that f (Y) has dimension of a length as its purpose is to control the gradients of the primal fields.Although one might expect that ℓ ε should be suitably larger than ℓ D , as defined in equation ( 1), our numerical experiments, in this investigation and in [6], have shown that this is not necessary in order to obtain an accurate and efficient numerical model.

Spatial discretisation and integration.
In [6] we have provided several details about different meshes adopted to model the IPMC behaviour.Here, we employ the mesh referred to as 'coarse' in [6], which is represented in figure 2 and has been shown to provide accurate results at a reasonable computational cost.
We describe this FE model by first focusing on its most important feature, which is the discretisation along the Y direction.The electrodes are discretised through four rows of eight-noded standard isoparametric FEs (each having 2 DOFs per node: u X and u Y ) with reduced integration.The membrane regions of size equal to ℓ ε next to the electrodes are discretised with eight rows (displayed in dark grey in figure 2, right) of GFEs (each having nine DOFs per node: u X plus 2 for each of u Y , ψ, μ∆ , and π).The rest of the membrane is discretised by 50 rows of eight-noded fully-integrated isoparametric FEs (each having five DOFs per node: u X , u Y , ψ, μ∆ , and π) through a double bias strategy leading to a ratio ≈1828 between the largest and smallest FEs, located at Y = 0 and |Y| = H − ℓ ε , respectively.About the discretisation along the X direction, the region from the clamped cross-section (X = 0, see figure 1) to X = 2H is divided into 16 identical FE columns, while from X = 2H to X = L the FE columns have width equal to H. Given that we always consider L/H = 200, the mesh is constituted by 214 columns of FEs.Overall, the FE model has 48 085 nodes and a total number of DOFs equal to 269 617.
About the spatial integration, let us point out that the enrichment function (25) requires a large number of Gauss points for an accurate integration in the GFEs.Given the nature of the expected solution, we adopt a non-isotropic distribution of the Gauss points, in which each GFE has 20 Gauss points along the Y directions and three of them along the X direction.

Summary of the main assumptions of the model and of its parameters
We believe that the strongest assumptions of the adopted model consist of using (i) perfectly flat electrodes that are impermeable to the free solvent, (ii) spatially uniform properties of the membrane, and (iii) a completely saturated membrane.Among other assumptions, we recall the we adopt a 2D model, isotropic elasticity for both membrane and electrodes, electrodes that are perfect conductors and are impermeable to hydrated counterions, ideal dielectric behaviour of the membrane, incompressibility of all the phases in the membrane (such as its macroscopic volumetric deformation strictly depends on the redistribution of free solvent and hydrated counterions), and ideal mixing behaviour of the phases disregarding the enthalpic contribution.
Here we summarise the parameters of the adopted computational model.
-Geometrical parameter: Of course, some of these parameters are obviously known (as θ), or kind of trivial to measure (e.g.L and B), while others are more difficult to be determined than one might think of (e.g.λ e , G e , and h because of the irregularity and porosity of the electrodes).A crucial difficulty is due to the heterogeneity in the membrane, both in terms of fields under operating conditions and in terms of composition after electrode deposition, which should require the use of spatially variable permittivity and diffusivities within the membrane.Note that ℓ ε may assume different values at anode and cathode, because of the asymmetry of the two BLs.The influence of these parameters on the model predictions has been discussed in our previous studies [6,15].

A simple yet accurate model for the material parameters identification relying on actuation
Plainly adopting the model summarised in section 2 for the identification of the material parameters easily leads to a computationally impossible task.To address this issue, we propose a relatively simple approach, where the IPMC subjected to a voltage drop across the electrodes ψ 0 is modelled as a Euler-Bernoulli beam (or, more precisely, given that in IPMCs usually B ≫ 2(H + h), as a Kirchhoff plate under cylindrical bending) subjected to an electrochemical bending moment M 0 (t) that is uniform, i.e. independent of the position along the IPMC axis X.
Within this approach, M 0 (t), which is due to the Maxwell and osmotic eigenstresses, is evaluated, as presented in section 3.1.2,by disregarding the structural mechanics, thus solving the electrochemo-poromechanical problem by considering only a 'through-the-thickness slice' of IPMC, i.e. a single column of FEs, with appropriate boundary conditions, such as no discretisation is required along X.This is possible under the reliable assumption that the electrochemical fields depend on Y only.
Then, once M 0 (t) is known, the solution of the structural problem, which admits the simple analytical expressions obtained in section 3.1.1,delivers the mechanical fields of interest.
Note that the proposed approach is significantly different from those in the literature that fully rely on the one-way coupling between electrochemistry and mechanics, where, in actuation, electrochemistry is solved first, usually analytically, by completely disregarding mechanics [11].

Structural model for the IPMC subjected to a known
electrochemical bending moment M 0 .In the proposed structural model the electrochemical moment M 0 is applied to the membrane only, whereby the electrodes hamper the free growth of the elastic curvature ensuing from the mechanical bending moment M m developed by the membrane in order to fulfill the equilibrium.Because of this and as shown below, the longitudinal elastic modulus of the electrodes affects the actuation displacement, while it has a mild influence on the BF, thus confirming previous findings by Vokoun et al [42]  3 .
Hence, we assume that each electrode develops an axial forces of magnitude N e , contributing to the equilibrium where, by restricting attention to linear elasticity, M m and N e are linked, respectively, to the elastic curvature in the membrane χ E and to the longitudinal strain in the electrodes ε through the membrane bending stiffness and the electrode axial stiffness, as ) Here and henceforth, Ẽ = E/(1 − ν 2 ) and Ẽe = E e /(1 − ν 2 e ) are the longitudinal moduli of the membrane and the electrodes, respectively, under plane strain and the usual structural assumption that the through-the-thickness normal stress component is negligible.
Compatibility requires that the longitudinal strain at the interface is the same if determined in the membrane or in the electrodes, that is x . ( Combination of equations ( 26)-( 29) allows the computation of the total curvature χ E experienced by the membrane Under the above assumptions, the IPMC tip deflection under actuation is which, from equation (30), is clearly dependent on Ẽe , with u Y (L, t) decreasing as Ẽe increases (e.g.lim Ẽe→∞ u Y (L, t) = 0). 3This is in stark contrast with the results obtained from oversimplified models in the literature, in which an electrochemical (eigen)curvature is assumed to be applied over the whole IPMC cross-section.Under this wrong hypothesis, the electrode elastic modulus turns out to affect the BF instead of uY(L, t).We further remark that the active strain approach in which one applies an (eigen)curvature is dual to the active stress approach followed here, where we apply an (eigen)moment, and, under the assumptions of this section, the two approaches lead exactly to the same conclusions if the electrochemical action is applied to the membrane only.Now, by employing the foregoing results, we can compute the blocking force F Y by resorting to conventional beam analysis of the propped-cantilever IPMC, thus obtaining where is the bending stiffness of the whole IPMC cross-section.By combining equations ( 30), (32), and (33), we finally have the BF as a function of the applied electrochemical moment and the geometrical and mechanical parameters of the linear elastic IPMC: From inspection of equation ( 34), we deduce that Ẽe should have a very low influence on F Y .In the very common case in which the second contribution at the right-hand side of equation ( 33) is negligible with respect to the sum of the first and third ones and H ≫ h/2, F Y is even totally independent of the elastic moduli, just reading Finally, by combining equations ( 30), ( 31), (33), and (34), we obtain which establishes that a given linear elastic IPMC subjected to a given ψ 0 (t) has ratio between the tip displacement in actuation and the BF that is independent of both time and ψ 0 (t).
Obviously, the main limit of the foregoing analysis is that it disregards the quite complex strain field in the membrane, for instance encompassing localisation of the throughthe-thickness direct strain component in the BLs [11,15].However, the computation of M 0 will instead account for the actual profiles of the electrochemical (eigen)stresses, along with the ensuing strain field.Moreover, the present analysis, by assuming that M 0 (t) is uniform along the IPMC axis, neglects the fact that the extra constraint needed to measure F Y (t) introduces, with respect to the cantilever IPMC where u Y (L) is measured, an additional elastic curvature χ F that is nonuniform along the IPMC axis.By assuming again absence of warping in the regime of small strains and rotations, χ F = F Y (L − X)/D.Such elastic curvature potentially affects the species fluxes, adding non-vanishing flux components along the IPMC axis and, thus, leading to an electrochemical bending moment dependent on X as well.This effect is however expected to play a role only in relatively short IPMCs, which are seldom adopted in practical applications.Boundary conditions applied to the reduced FE model employed in the identification of the model parameters: through the control point CP1 we impose vanishing rotation φ Z of the otherwise unconstrained right side of the modelled IPMC region and, consequently, we obtain the mechanical bending moment as a reaction.Note that, in order to have optimal numerics, we also impose, at the membrane sides A i and B i , periodic boundary conditions on ψ, μ∆, and π, although, in an exact solution, these conditions would ensue from D X = J s X = J X = 0, which are applied therein.
In section 4.1, we will further discuss equation (36) in the light of experimental results to motivate its use as an extremely convenient constraint in the challenging identification of the IPMC material parameters.

FE model for the determination of the electrochemical moment M 0 (t).
With reference to figure 2, the electrochemical problem can be solved by a single column of FEs meshing the IPMC along its thickness.By adopting the same spatial discretisation described in section 2.3.1, this FE model consists of 74 FEs only.In this FE model, there are not constraints to the direct deformation component along the through-thethickness direction Y and the two membrane sides of normal ±X are subjected to D X = J s X = J X = 0.About the other mechanical boundary conditions, two strategies can be followed.
The first one consists of fixing the longitudinal displacements of all the nodes at the two sides of normal ±X to compute the mechanical bending moment developed by the reaction forces, M m + N e (2H + h), which is equal and opposite to the electrochemical bending moment M 0 (t).As illustrated in figure 3, in ABAQUS we have conveniently imposed this condition through the use of a 'control point' [43] that, actually, constrains only the rotation of the right side of this reduced model.
The second strategy consists of fixing the longitudinal displacements of all the nodes of one side only, say that of normal −X, while imposing that the nodes at the other side must follow a reasonable kinematics, e.g.either fully remaining on a straight line or obeying a zigzag warping [37], where the rotation of the membrane is independent from that of the electrodes.In this case, one can directly estimate the elastic curvature of the membrane χ E .Although assuming a zigzag warping provides a more precise parallel with the simplified theory of section 3.1.1,and is in general more suitable for sandwich panels [27], it has a mild influence in IPMCs and, most of all, it leads to numerical issues in the FE model at the membrane-electrode interfaces.
In our numerical attempts the two strategies have always provided negligible differences, such that in this paper we present only results obtained by following the first one.

Evolutionary algorithm
The model parameters are those minimising a target function T , which is referred to as the fitness function in the terminology of evolutionary algorithms, that sums up the squares of the relative errors with respect to N considered experimental data.Hence the fitness function is defined as where D exp i and D num i are, respectively, the experimental and numerical data of the ith test.In our case concerned with actuation, each datum is either the value of the tip displacement, u Y (L, t), or the value of the BF, F Y (t), at selected instants t.
Since the analytical form of the gradient of T is unknown, we resort to a specific evolutionary algorithm, which is the non-gradient global optimiser referred to as Coliny in the software Dakota [44].In computational intelligence, these algorithms inspired by biology optimise on the basis of a population of individuals.Here, an individual is a single selection of the set of model parameters to optimise plus the related value of T .As per advise of the Dakota manual, our population consists of 40 sets of model parameters that, initially, are randomly generated, although within certain bounds set by us.
In order to numerically obtain T , we have implemented a code coupled with Dakota that runs the required FE analyses by using the model of section 3.1.2and, then, computes u Y (L, t) and F Y (t) by resorting to the simplified relations of section 3.1.1.The FE analyses required depend, of course, on the experimental data that one wants to use, among those available, in the parameter optimisation.
As a result of each iteration, Dakota rebuilds the population of 40 individuals by changing 30 of them, while keeping the best ten [44].In section 4.2 we provide the algorithms 1 and 2 specific of the identification performed in this investigation.
It is important to emphasise that, after the model parameters are determined through the foregoing identification procedure, the final results presented in the following are obtained by running simulations with the full 2D GFE model.

Results and discussion
Although we are aware of the very large scatter in the experimental data on IPMC actuation and on the influence of the environmental conditions on them [17,18,45], here we test our model against the data obtained by He et al [26], who, unfortunately, do not provide error bars on the recorded data.These data are concerned with Nafion™-Pt IPMCs saturated with a solution of water and sodium counterions.We have Table 1.Experimental results of He et al [26] for the maxima of tip displacement and blocking force.by applying variable peak voltage drops ψ p 0 = 2, 2.5, 3, 3.5 V, through a sinusoidal voltage signal of 1 Hz frequency.Unfortunately, He et al [26] do not provide the time variation of the tip displacement and the BF, nor the times at which their maxima are attained.This notwithstanding, the data in [26] are rich enough to constitute a good test for our proposal.
Although He et al [26], through Berkovich indentation on plain membranes, found a dependence of the membrane elastic moduli on the membrane thickness, we have neglected such dependence and identified a single couple of elastic constants for the membrane of all the samples.The main reason for this choice is that if we considered such dependence, to be consistent, we would have to do it also for other material parameters of the membrane, such as the permittivity, whose value should probably be even assumed to be a function of Y, as discussed later.Although this would probably allow an excellent prediction of all the experimental results in [26], it would make our contribution much less significant because of the relatively small number of data that we have available.
The IPMCs tested by He et al [26] have length, electrode thickness, and width, respectively, equal to L = 20 mm, h = 0.01 mm, B = 5 mm.The slenderness of these IPMCs suggests that their electrochemistry should be very slightly affected by the non-uniform curvature χ F in the BF tests.Even though it is very unlikely that the measure h = 0.01 mm is accurate and uniform for all the samples, we keep this value fixed, as provided in [26].
We begin by discussing the simplified structural model presented in section 3.1.1,where the tip displacement u max Y and the BF F max Y are provided in terms of a spatially uniform electrochemical moment M 0 (t).

Preliminary validation of the simplified model relying on decoupling structural mechanics
Table 1 reports the experimental results of He et al [26].Of course, in reality, u max Y and F max Y are not necessarily attained at the same time, even for the same IPMC under the same ψ 0 (t), because of experimental uncertainties.
Table 2. Ratio between the maxima of tip displacement and blocking force: experimental results of He et al [26] against equation (36) with elastic moduli of electrodes and membrane as obtained from the identification procedure presented in section 4.2.Table 2 shows that the simplified analysis of section 3.1.1,which leads to the result in equation (36) suggesting that the dependence on both ψ 0 (t) and the electrochemical parameters of the ratio u Y /F Y may be unimportant, is closer to the actual behaviour for thicker IPMCs.This is inferred by evaluating the uniformity of the results in each row of the second column in table 2.

2H(mm) u max
The results in the last column of table 2 are obtained by plugging in equation ( 36) the elastic moduli λ = 2028.95MPa, G = 20 MPa (such as Ẽ ≈ 79 MPa), and Ẽe = 6568.41MPa that are obtained from the identification procedure as presented in section 4.2.These values of the elastic constants are a bit smaller than those adopted in other investigations on IPMCs and, in particular, those of the membrane are significantly smaller than those measured by He et al [26] through Berkovich nano-indentation, who obtained an elastic modulus value increasing with the membrane thickness, from ≈0.2 GPa for 2H = 0.22 mm to ≈0.8 GPa for 2H = 0.8 mm.However, the values obtained by the present identification procedures seem to be reliable because they enter the bending stiffness D in equation (36), whose application matches quite well the experimental results.The large discrepancy for the thinnest IPMC at highest ψ p 0 might be due to larger influence of the BL regions, where highly nonlinear phenomena occur.In next section 4.2 we will show that also the full nonlinear GFE model cannot capture u Y (L) for 2H = 0.22 mm at ψ p 0 = 3.5 V, while anyway providing a good estimate for the BF.

Identification of the model parameters
For the identification procedure described in section 3.2, we have selected the six model parameters to be optimised by the evolutionary algorithm, although in the electrodes only the longitudinal modulus Ẽe plays a significant role, as verified later.
About the other parameters, from the data provided in [26] we estimate that C 0 v = 0.1 and introduce this as a constraint; moreover, from [26] we know that the total initial concentration of solvent (i.e. the free solvent plus the solvent that is bonded to the counterions), C tot s , is equal to 1.55 • 10 −5 mol mm −3 , and that three molecules of solvent are bonded to each counterion.Hence, once we have an estimate of C 0 , we obtain v = 0.1/C 0 and C 0 s = C tot s − 3C 0 .For v s , D, and D s we maintain the values obtained in [6].The diffusivities basically affect the instant in which u max Y and F max Y are attained, such instants being unknown to us as they are not provided in [26].The molar volume of the water is set to v s = 20 000 mm 3 mol −1 (some studies adopt the slightly different value v s = 18 000 mm 3 mol −1 ).
Also the GFE parameters ℓ ε and α are kept as selected in [6], and their suitability has been checked by numerical experiments.
We identify the model parameters by considering only the three intermediate values of the membrane thickness, 2H = 0.32, 0.42, 0.64 mm and only the two peak applied voltages ψ p 0 = 2 and 3 V.Hence, the fitness function T in equation ( 37) to minimise is given by N = 12 addends.
Note that we have computed u max Y and F max Y under the application of a single voltage cycle, because we have checked that, in our simulations, the amplitude of the response diminishes, although slightly, with increasing the cycle number.
In all the analyses the absolute temperature is set to θ = 300 K.
By adopting the foregoing assumptions, algorithms 1 and 2 further illustrate and summarise how this specific identification of the model parameters works with reference to the general procedure described in section 3.2.
4 for each H i ∈ {0.16 mm, 0.21 mm, 0.32 mm} do 5 compute 6 for each ψ p 0j ∈ {2V, 3V} do 7 write ABAQUS input file for H i , ψ p 0j and the active set of parameters; 8 run ABAQUS standard with the GFE UELsubroutine [6] untilM 0 (t) reaches the value M p 0ij , which is the peak value for the current parameters H i and ψ p 0j ; 9 read M p 0ij from ABAQUS output file; 14 write T into 'results.out'file.

Computational cost of the identification procedure.
Overall, the identification required 15 217 evaluations of T in equation (37), which corresponds to 91 302 FE analyses (with the reduced model described in section 3.1.2).The total number of iterations, in which the procedure evaluates all the set of parameters constituting the population, was 503.The total computational cost resulted equal to 368 019 s by employing 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.Note that the above data account for the fact that we had to stop and restart Dakota three times to adjust the bounds of the parameters to optimise (this happens when a parameter gets too close to one of its imposed bounds).Details on the computational cost of the full 2D FE model are provided in [6].However, to give an idea of the computational convenience of the simplified model adopted in the identification procedure, a preliminary attempt of identification by employing the full 2D FE model, where we used only the BF data, required a computational time of 250 738 s for 91 iterations and 2766 evaluations of T .by heuristically optimising them on a single datum, that was u max Y for 2H = 0.42 mm at ψ p 0 = 2.0 V, and we completely disregarded the BF data.Note that, as already explained, in the present identification we adopted the same diffusivities and GFE parameters as in [6].
The comparison between the experimental data and our predictions is shown in figures 4 and 5, in terms of u max Y and F max Y , respectively.
As reported in table 3, the sole modulus of the electrodes that influences the IPMC behaviour is the longitudinal modulus Ẽe , which results equal to 6568.41 MPa, thus well agreeing with the value of the longitudinal modulus adopted by Vokoun et al [42] to model the same experimental data, equal to 5 GPa4 .All the results presented here refer to ν e = 0.110 94, corresponding to λ e = 832.61MPa, G e = 2919.85MPa, E e = 6487.57MPa 5 .
By plugging the parameters of table 3 in equation ( 1) we obtain a Debye screening length ℓ D ≈ 5.1 • 10 −3 mm ≈ 5ℓ ε , meaning that the thickness of the BLs may be larger than the size of the regions discretised with GFEs.We have however addressed this issue by verifying that FE analyses carried out with ℓ ε increased by one order of magnitude lead to the same response in terms displacement and BF as functions of time.
Given that in Nafion™ C 0 ≈ 10 −6 mol mm −3 , ℓ D is basically determined by ε, whose value impacts the amplitude of the actuation response, although it does not much affect the relative contributions of Maxwell and osmotic stress [5,6].There has been a long debate in the literature about the value that should be assigned to the Nafion™ permittivity in IPMCs.For instance, Wallmersperger et al [20] identified a relative permittivity ε/ε 0 = 120 (where ε 0 ≈ 8.8542 • 10 −12 C (V m) −1 is the vacuum permittivity), which results in ε ≈ 10 −15 C (mV mm) −1 .Instead, Arnold et al [18] have measured permittivity values about four orders of magnitude larger than that reported in [20], and Barramba et al [46], still on Nafion™-Pt IPMCs, even measured ε ≈ 3.0 • 10 −8 C (mV mm) −1 , which agrees well with our findings.Hence, the literature values of the Nafion™ permittivity span at least seven orders of magnitude.
The permittivity, in fact, is very much affected by the IPMC manufacturing and material adopted for the electrodes.In our view, one of the main issue is that in IPMC membranes ε may be largely non-homogeneous depending on the species redistribution, which in turn is a complex function the boundary conditions and the IPMC dynamics.Moreover, the dependence of the permittivity on the deformation, as observed in some EAPs [47], might need to be accounted for in the BLs, where the through-the-thickness direct strain component may be quite large [11,15].Additionally, even IPMCs at rest display important sources of heterogeneity in the membrane.In fact, its regions adjacent to the electrodes, referred to as 'composite layers' [45,48] or 'intermediate layers' [49] in the literature, are partly filled with metal particles from the electrode manufacturing, hence leading to a much larger permittivity, and much lower diffusivity, with respect to the bulk membrane [13,38].Some efforts in the literature [50] even adopt a fractal approach to model the rugosity of membrane-electrode interfaces.Summing up, we believe that in IPMCs the permittivity should not be regarded as a property of the bulk membrane (e.g.Nafion™), but as a property of the whole system.As a more sophisticated alternative approach, our FE modelling would allow one to straightforwardly assign a spatially heterogeneous permittivity ε(Y); however, this would require to be supported by an ad hoc experimental characterisation.Back to the results listed in table 3, the value identified for C 0 , equal to ≈7.8 • 10 −7 mol mm −3 , is quite close to that usually adopted in the literature, C 0 ≈ 1.0 • 10 −6 mol mm −3 , up to 1.2 • 10 −6 mol mm −3 .In particular, the obtained value for C 0 is much closer to its reference value than that calibrated in [6], which resulted C 0 = 2.0 • 10 −7 mol mm −3 .Given that the product vC 0 should be approximately constant, the increase of C 0 corresponds to a decrease of v, from v = 500 000 mm 3 mol −1 in [6] to the much more reliable value v = 128 015 mm 3 mol −1 obtained in this work.
The model parameters identified in the present investigation allow us to obtain results from the full 2D FE model for all the applied voltages, which was not possible with the parameters as calibrated in [6], that hampered convergence of the FE analyses for ψ p 0 = 3.5 V, for which a more refined GFE was needed, albeit not pursued on that occasion.
As displayed in figures 4 and 5, the match is very satisfactory overall.In the case ψ p 0 = 3.5 V, although the match on the BF is good, there is a large disagreement on u max Y for the thinner IPMCs, increasing with decreasing H. Since we have checked that our full FE model is consistent with the simplified analysis of section 3.1.1,adopting linear elasticity, we deduce that the problem is not in the adopted nonlinear elastic laws for membrane and electrodes.As already mentioned, there might be several reasons for this discrepancy, beside the scatter in the experimental results that is unavailable to us.Most of them are concerned with the increased relative importance of BLs on thinner IPMCs.In fact, BLs govern the IPMC behaviour being subject to large concentration gradients and, incidentally, to large through-the-thickness direct strain [11,15].They therefore surely constitute the difficult part of the modelling, so far neglecting water electrolysis, the roughness of the membrane-electrode interfaces, the presence of metal particles in the 'composite layers', the membrane mechanics beyond hyperelasticity, and so on.
For the sake of completeness, in figure 6 we provide a graphical output of the full 2D FE analysis of the BF problem at ψ p 0 = 3.5 V.At the instant in which F max Y is recorded, it represents the IPMC deformed shape and the contour of the through-the-thickness component of the counterion flux, J Y .
Figures 7-10 compare the results obtained with the simplified model adopted in the identification procedure and the full 2D FE model.In particular, figures 7 and 8 focus on the timevarying responses, for ψ p 0 = 2 V and 3 V, respectively.The responses are provided for both F Y (t) (left) and u Y (L, t) (right) and for three thickness values, 2H = 0.32, 0.42, 0.64 mm.The discrepancy between the two models is relatively small; however, we observe that, as expected on the basis of the foregoing discussion, the discrepancy is larger for u Y (L, t) and it increases with the applied voltage and with diminishing the membrane thickness.
Figures 9 and 10 are instead concerned with the contours of the free solvent concentration and the electric potential at the peak response in the case ψ p 0 = 2 V and 2H = 0.64 mm.Noticeably, there is almost no difference at all between the reduced model (left images) and the full 2D model (right images), the latter being referred to the BF benchmark and displayed only in a representative region by limiting it along the IPMC axis X.Note also that the size along X of the reduced model is irrelevant given its boundary conditions and its purpose (see section 3.1.2).As a final note, we remark that we obtained poor results, surely unable to predict the experimental data of He et al [26], by neglecting the cross-diffusion of solvent and counterions, i.e. by setting γ = 0 in equations ( 23) and (24).

Concluding remarks
We have proposed an identification procedure to determine the material parameters entering the electrochemoporomechanical theory for IPMCs developed in [6,15] and here referred to as 'LBP theory'.A crucial prerequisite to apply an identification procedure is the capability of the model to deliver solutions for external actions consistent with experiments and applications, and for a wide range of material parameters.This is the case of our computational model, which adopts the GFE implementation developed in [6].The GFE technology constitutes a numerical tool able to overcome a fundamental issue in the numerics concerned with IPMCs, that is the crucial role of the BLs of counterion concentration forming in the membrane at the interfaces with the electrodes [5,7,9,36].
Although the theory works both in actuation and sensing, here the focus is on actuation, which, by the way, is the most promising field of applications for IPMCs.The proposed identification procedure has been developed in order to be computationally affordable.It consists of decoupling the structural mechanics problem, thus limiting the electrochemoporomechanical GFE analysis to a model that discretises the IPMC only along its thickness.The main assumptions are that the electrochemical fields depend on the IPMC through-thethickness coordinate only and that the IPMC structural behaviour can be obtained by applying to the IPMC membrane a spatially uniform electrochemical bending moment which depends on the Maxwell and osmotic stresses.This approach is applied to two experiments typical of IPMC actuation, which aim at measuring the tip displacement in a cantilever IPMC and the BF in a propped-cantilever IPMC.In the second case, the BF denotes the reaction force developed by the support constraining the otherwise free IPMC end.
The simplified model has been plugged into the evolutionary algorithm Coliny in the software Dakota [44] to identify the material parameters against the experimental results of He et al [26], who provide the maximum tip displacement and maximum BF of IPMCs of variable membrane thickness subjected to sinusoidal applied voltages of variable peak level ψ p 0 .After the material parameters are identified, the final results are provided by running GFE analyses on the full 2D electrochemo-poromechanical plane-strain model for IPMCs.
In view of both the large scatter of the experimental data on IPMCs [17,18,45] and of the strong uncertainties on the IPMC characterisation, the obtained results are very satisfactory to us.The match between predictions and experimental data is always good in terms of BF, while significant discrepancies on the maximum tip displacement arise for the thinner samples at the largest ψ p 0 , which is equal to 3.5 V. We have provided arguments to justify the observed discrepancy, ascribing it to phenomena and IPMC features unaccounted for by our model, the latter mainly concerned with the precise description of the interfaces between membrane and electrodes.A conceptually simple option that might be able to fix the observed discrepancy consists of assuming spatiallyheterogeneous material parameters in the membrane.
Although, for the sake of brevity, we have not presented the results of failed attempts to model the experimental results of He et al [26], we can state that our study has confirmed the crucial role of the specific constitutive prescriptions for the fluxes of counterions and solvent.They in fact hugely impact the IPMC peak response and accounting for the cross-diffusion of counterions and solvent, as proposed in [15], seems to be almost mandatory to be able to identify material parameters that may be physically reasonable.

Open issues
About the modelling, the most important issues that still should be addressed are concerned with the accurate description of the interface regions between membrane and electrodes [36,48,49,51,52] and with accounting for the possibility of the solvent to exit the IPMC, as a function of the environmental conditions [17,18].In the present modelling, the interfaces are perfectly flat and impermeable to both hydrated counterions and free solvent, whereby the membrane is completely saturated.
Several other details could improve the modelling capability, including (i) a further improvement of the mobility matrix governing the species fluxes [40], (ii) a mechanical law for the membrane that describes its behaviour beyond hyperelasticity, thus enabling stresses in the BLs to relax [53], (iii) accounting for a so far missing enthalpic contribution in the free energy of mixing [54], (iv) the analysis of the membrane small-scale behaviour [55], (v) the implementation of nonideal behaviours in the species interactions [56].
Before all of this, in our opinion, a much stronger interaction between experimental and theoretical investigations is essential to make further progresses.Given the complexity of the IPMC multiphysics and the huge scatter in the experimental data on IPMCs [17,18], ad hoc experimental campaigns should be designed in order to more efficiently develop the actual state of the art on the modelling side.For instance, in order to test an electrochemo-poromechnical theory that aims at being comprehensive as that adopted here, it would be extremely important to have both actuation and sensing data on samples consistently produced.In this investigation we could focus on actuation only, but the dual tests on short-circuit sensing would provide data on the electric current that in the model could be predicted as is the charge accumulated at the electrodes [15,37,38,57].Moreover, it would be crucial to produce IPMCs whose behaviour displays far less scatter.Novel techniques are expected to improve and stabilise the IPMC performance.Among them, we mention additive manufacturing [4,58], the membrane hot-pressing [59], and the use of additives such as ionic liquids [60] or the dielectric gel coating [46] for IPMCs operating in open air or at variable humidity.

Figure 1 .
Figure 1.Schematic of the blocking force problem: geometrical parameters, reference system, and physical constituents.The inset at the right shows the composition of a macroscopic material point of the membrane, where counterions are hydrated with a shell of solvent.

Figure 2 .
Figure 2.Images of the FE mesh, from[6]: discretisation over the whole IPMC thickness at the fully-clamped cross-section region (left) and detail about the bottom electrode (right), whose domain is displayed by the bottom four rows in light grey, with the adjacent membrane region including the GFE discretisation, which is displayed in dark grey.
thickness of each electrode • B IPMC width.-Mechanical material constants: • λ first Lamé constant of the membrane • G shear modulus of the membrane • λ e first Lamé constant of the electrodes • G e shear modulus of the electrodes.-Electrochemical material constants (of the membrane): • ε dielectric permittivity • v molar volume of hydrated counterions • v s molar volume of the solvent • C 0 initial nominal molar concentration of counterions • C 0 s initial nominal molar concentration of free solvent • D diffusivity of the hydrated counterions • D s diffusivity of the free solvent.-Temperature:• θ. -Parameters of the GFE numerical model:• ℓ ε thickness of the membrane regions discretised with GFEs • α parameter controlling the gradients in the BLs.

Figure 3 .
Figure3.Boundary conditions applied to the reduced FE model employed in the identification of the model parameters: through the control point CP1 we impose vanishing rotation φ Z of the otherwise unconstrained right side of the modelled IPMC region and, consequently, we obtain the mechanical bending moment as a reaction.Note that, in order to have optimal numerics, we also impose, at the membrane sides A i and B i , periodic boundary conditions on ψ, μ∆, and π, although, in an exact solution, these conditions would ensue from D X = J s X = J X = 0, which are applied therein.

Algorithm 1 .end 7 endAlgorithm 2 .
Calibration procedure with Dakota.Result: minimise the fitness function T in equation (37) by using genetic algorithms Output: set of optimised parameters C 0 , ε, G, λ, Ge, λe 1 for each Dakota iteration do / * Generate a population of 40 individuals (set of parameters and fitness function): / * 2 for each individual do 3 compute C 0 , ε, G, λ, Ge, λe by using the Coliny-ea evolutionary algorithm and write their current values into 'params.in'file; 4 run algorithm 2 to compute the fitness T of the current set of parameters; 5 read the fitness function T of the current individual from 'results.out'file; 6 Compute the fitness function T of a single individual.Input : 'params.in'file from Dakota Data : Geometry: L = 20 mm, B = 5 mm, h = 0.01 mm, 2H = {0.32,0.42, 0.64} mm [26] Data : Experimental maximum values of the tip displacement, u max,exp Yij , and of the BF, F max,exp Yij , for specific values of membrane thickness H i and peak applied voltage ψ p 0j [26] Data : fixed material parameters: Ds = 10 −3 mm 2 s −1 , D = 10 −6 mm 2 s −1 , vs = 20 000 mm 3 mol −1 , C 0 v = 0.1, C tot s = 1.55 • 10 −5 mol mm −3 Data : GFE parameters: ℓε = 0.001 mm, α = 1000; Output : 'results.out'file containing the value of the fitness function T

Figure 4 .
Figure 4. Displacement tip in actuation: comparison between predictions and experimental results of He et al [26].

Figure 5 .
Figure 5. Blocking force in actuation: comparison between predictions and experimental results of He et al [26].

Figure 6 .
Figure 6.Displacement field (amplified 16.45 times) in the BF problem along with the contour of the through-the-thickness counterion flux J Y for the case 2H = 0.8 mm with ψ p 0 = 3.5 V at the instant in which the BF attains its maximum value.

Figure 7 .
Figure 7. Comparison, in terms of blocking force and tip displacement, between the reduced FE model used in the identification procedure (1D, by referring to the electrochemistry) and full (2D) FE model for ψ p 0 = 2.0 V.

Figure 8 .
Figure 8. Comparison, in terms of blocking force and tip displacement, between the reduced FE model used in the identification procedure (1D, by referring to the electrochemistry) and full (2D) FE model for ψ p 0 = 3.0 V.

Figure 9 .
Figure 9. Peak response for ψ p 0 = 2 V and 2H = 0.64 mm in terms of Cs: comparison between the full 2D FE model (right) and the reduced FE model used in the identification procedure (left).

Figure 10 .
Figure 10.Peak response for ψ p 0 = 2 V and 2H = 0.64 mm in terms of ψ: comparison between the full 2D FE model (right) and the reduced FE model used in the identification procedure (left).