Mireille Tixier
email: mireille.tixier@uvsq.fr

Joël Pouget
email: pouget@lmm.jussieu.fr

& J Pouget

Validation of a model for an ionic electro-active polymer in the static case

Keywords: Electro-active polymers, Multiphysics coupling, Polymer mechanics, Nafion, EAP modelling, Ionic polymer, EAP beam PACS numbers : PACS 47.10.ab, PACS 83.60.Np, PACS 82.47.Nj, PACS 77.65.-j

ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

The present work proposes a study of a thin strip of electro-active polymer (EAP). More precisely, the study focuses on Ionic Polymer-Metal Composite (IPMC) belonging to the ionic class. We investigate responses of the EAP subject to an applied difference of electric potential on the metallic electrodes and an applied punctual force at the tip of the strip.

In previous works we have constructed step by step a micro-mechanical model to establish the conservation laws for electro-active polymers (Tixier & Pouget 2014). Following this work, the constitutive equations were deduced from the hypothesis of local thermodynamics equilibrium and the Gibbs relation using the thermodynamics of irreversible processes (Tixier & Pouget 2016). The present study follows on from the spirit of the previous works and we want to characterize the behavior of a thin blade, especially the chemical, electrical and mechanical parameters under electro-mechanical loading. The behavior of electro-chemical-mechanical interactions of EAP is of great interest for research and engineering advanced technology. Simply, we can say that an EAP is a polymer exhibiting a mechanical response, such a stretching, contracting or bending for example, when subject to an electric field. Conversely the EAP can produce electric energy in response to a mechanical load [START_REF] Shahinpoor | [END_REF], Shahinpoor & Kim 2001, Bar-Cohen 2005, Pugal et al. 2010, Park et al. 2010). This particular property is highly attractive for applications. For instance, we can quote biomimetic devices (robotics, bio-inspired underwater robots such as fishes [START_REF] Chen | [END_REF], haptic actuators (artificial skin, tactile displays or artificial muscle (Deole et al. 2008[START_REF] Matysek | Proc. SPIE[END_REF], Bar-Cohen 2005) and this material is an excellent candidate for energy harvesting (Farinholt et al. 2009, Tiwari et al. 2008, Aureli et al. 2010[START_REF] Jean-Mistral | [END_REF], Cha et al. 2013). In addition promising applications to micro-mechanical systems (MEMS) at the sub-micron scale are also expected in medical engineering for accurate medical control or investigation, for instance Fang et al. (2007) and Chikhaoui et al. (2018).

Electro-active polymers can be categorized into two main groups depending on their activation mechanisms. The first one is electronic electro-active polymers and dielectrics which are subject to Coulomb force ; their volume change is due to the application of an electric field. For instance dielectric elastomers belongs to this category as well as piezoelectric and electrostrictive polymers. The second group of EAP is the ionic electro-active polymers. The latter are driven by the displacement of ions inside the material. The polyelectrolyte gels, ionic polymer-metal composites, conductive polymers are such media. One of the advantages of the these polymers is that they can be activated by very low difference of electric potential of about 1 -5 Volts. However, they can only be operated within electrolyte medium.

Most efforts have been devoted on electronic conducting polymers (E.C.P.) belonging to the ionic class in order to improve strain, output forces and response times. These polymers are used for multilayer polymer composite combined with a polymer which can be considered as ion reservoir to improve ion transfer. Trilayer actuators made of solid polymer electrolyte film sandwich between two electronic conducting polymers have been reported in Festin et al. (2014) and Nguyen et al. (2018). One of the advantage of these electro-active polymers is that they can operate in air and they are good candidates for biomimetic devices. Efforts have been conducted to increase the integrating biocompatible conducting polymers into continuum micro-robots. Developments of such biocompatible conducting polymers with the aim at designing accurate position control for the trajectory of the telescopic soft robots have been studied by Chikhaoui et al. (2018).

Modeling EAP is not an easy task, especially the model must include electro-mechanical and chemical-electric couplings of ion transport, electric field action and elastic deformation. Different kinds of approach have been proposed in the literature according to the underlying physics of chemical activation (Shahinpoor & Kim 2001, Brunetto et al. 2008, Deole et al. 2008, Bahramzadeh & Shahinpoor 2014). An instructive and comprehensive review paper devoted to IPMC has been proposed by Jo et al. (2013). The authors present the chemico-physical mechanisms involved in IPMCs. As a complement to the previous references, a set of interesting studies has been devoted to the modeling of EAP. Among them, Wallmersperger et al. (2009) develope a thermodynamically based mechanical model involving the chemico-electric transport phenomenon. Numerical simulations are proposed for a strip of EAP made of Nafion ® 117 Li + . The authors deduce the profiles, within the strip thickness, of the electric charge density, electric potential, electric field and strain. [START_REF] Nardinocchi | [END_REF] deduce a model based on the 3-D theory of linear elasticity. The thermodynamics allows them to introduce chemo-electro-mechanical coupling and they deduce the constitutive equations of the material, especially a Nernst-Planck like equation is deduced. Their study continues with numerical illustrations for a thin strip of EAP. Moreover, time evolution at low frequency is proposed as well. It is noteworthy that these papers are partly, in their spirit, rather close to the present model. Nevertheless, their results depend strongly on adjustable parameters to fit the experimental results. In their model, the Fourier law, Darcy law and the generalized rheological constitutive equation are not presented. Moreover, using their adjustable phenomenological parameters the authors deduce a dielectric permittivity greater than ours.

A continuum approach for multiphasic materials has been investigated by Bluhm et al. (2016). The authors write down the conservation laws and the entropy inequality based on the theory of porous media for an arbitrary number of individual constituents. In spite of the rather thorough investigation, the authors do not compare their results with those proposed in the literature. However, their constitutive laws are related to the material constituents and not to the macroscopic medium, and their coefficient are not numerically evaluated.

The present model is mainly based on an averaging method of the microscopic description of phases in order to deduce the macroscopic behavior of the material. Thanks to this approach, the constitutive coefficients for the whole material are computed with the help of the microscopic physico-chemical properties of the constituents (volume fractions, mass densities, chemical potentials, electric charges, etc.) and physical meanings are identified.

An important polymer property that we would like to address is the dependency of the dielectric permittivity with the cation concentration in the polymer. As matter of fact, the dielectric constant of the polymer is not absolutely homogeneous within the strip thickness. Accordingly, it seems that the ionic transport when a difference of electric potential is applied to the polymer strongly modifies the dielectric constant along the thickness direction of the blade. This question will be examined and discussed in details in the forthcoming sections.

The paper is organized as follows: the next section is devoted to the description of the IPMC and chemicophysical process of activation is briefly described. This section reports also the way of modeling the EAP, especially, the conservation laws and constitutive equations are summarized. Section 3 focuses on the application of the model to a slender beam made of thin layer of EAP. A part of the section highlights the influence of the cation concentration within the thickness on the dielectric permittivity. Three kinds of dielectric laws will be considered. Numerical simulations are presented in Section 4. The profiles of the relevant variables and the scaling laws are reported. Comparisons to the experimental data available in the literature are discussed according to the chosen law of the dielectric permittivity.

Modelling of the polymer

Description of the material and hypothesis

We study an IPMC (ionic polymer-metal composite): the system consists in an ionic electro-active polymer (EAP) blade of which the two faces are covered with thin metal layers acting as electrodes. The model we developed applies for example to Nafion ® , an ionic polymer well documented in the literature that we will use for the validation of our constitutive equations. When saturated with water, this polymer dissociates quasi completely, releasing small cations in water whereas anions remain bound to the polymer backbone (Chabe 2008). When a potential difference is applied between the two electrodes, the cations migrate towards the negative electrode, carrying the water away by osmosis. As a result, the blade contracts on the side of the positive electrode and swells on the opposite side, causing its bending (Figure 1). To model the electro-active polymer, we used a "continuous medium" approach. Negatively charged polymer chains are assimilated to a deformable, homogeneous and isotropic porous medium in which flows an ionic solution (water and cations). The system is therefore composed of two phases and three components which move relative to each other: the cations, the solvent (water) and the porous solid. The solid and liquid phases (water + cations) are separated by an interface without thickness. The three components are respectively identified by the subscripts 1, 2 and 3; 4 denotes the solution (1 + 2) and i the interfaces. Gravity and magnetic induction are supposed to be negligible; the only external force exerted is therefore the electric action. The different phases are supposed to be incompressible and the solution diluted. It is further recognized that solid deformations are small.

Basis of the model

We used a coarse-grained model developed for twocomponent mixtures [START_REF] Ishii | Thermo-fluid dynamics of twophase flow[END_REF]. We define two scales. The conservation equations are first written at the microscopic level for each phase and for the interfaces. At this scale (typically about 100 A • ), the elementary volume contains one phase only but is large enough so that the medium can be considered as continuous. The macroscopic equations of the material are deduced by averaging the microscopic ones using a presence function for each phase and interface. For each phase k and for the interfaces, a Heaviside-like function of presence is defined. The macroscale quantities are obtained by averaging the corresponding microscale quantities weighted by the functions of presence.

This volume average is assumed to be equivalent to a statistical average (ergodic hypothesis). At the macroscopic scale, the representative elementary volume (R.E.V.) must be large enough so that these averages are relevant, but small enough so that the average quantities could be considered as local. According to Gierke et al. (1981) and Chabe (2008), its characteristic length is about 1 µm.

To write the balance equations, it is necessary to calculate the variations of the extensive quantities for a closed system in the thermodynamic sense. At the microscopic scale, we use the particle derivative or derivative following the motion of a constituent or an interface.

At the macroscopic scale, the three constituents velocities are different; we introduce a "material derivative" or derivative following the movement of the three constituants D Dt , which is a weighted average of the particle derivatives related to each constituent [START_REF] Coussy | Mechanics of porous continua[END_REF], Biot 1977).

Conservation laws

We thus obtain balance equations of mass (1), momentum (2), total, kinetic, potential and internal energy densities (3), entropy, electric charge (4) and the Maxwell equations (5) for the complete material (Tixier & Pouget 2014)

∂ρ ∂t + div ρ -→ V = 0 , (1) ρ
D -→ V Dt = div σ + ρZ -→ E , (2)
ρ D Dt U ρ = 3,4 σ k : grad V k +   I - k=3,4,i ρ k Z k V k   • E -div Q , (3) div -→ I + ∂ρZ ∂t = 0 , (4) rot -→ E = -→ 0 , div -→ D = ρZ , -→ D = ε -→ E , (5)
where ρ k denotes the densities relative to the volume of the whole material, -→ V the velocity, σ the stress tensor, Z the electrical charge per unit of mass,

-→ E the electric field, U the internal energy density, -→ I the current density vector, -→ Q the heat flux, -→ D the electrical induction and ε the dielectric permittivity. Subscripts refer to a phase, interface or a constituent and quantities without subscript are relative to the whole material.

We verify that the stress tensor of the whole material is symmetrical. The second member of equation (3) highlights the source terms of the internal energy (viscous dissipation and Joule heating) and its flux -→ Q . Equations ( 4) and ( 5) show that an EAP behaves like an isotropic homogeneous linear dielectric. In the last Maxwell equation, the permittivity of the whole material is obtained by a local mixing law, and is therefore likely to vary over space and time.

Constitutive equations

We make the hypothesis of local thermodynamic equilibrium. The Gibbs relation of the whole material is deduced from the Gibbs relations introduced by de Groot et [START_REF] De Groot | Non-equilibrium thermodynamics[END_REF] for a deformable solid and for a two-constituent fluid

T d dt S ρ = d dt U ρ + p d dt 1 ρ - k=1,2,3 µ k d dt ρ k ρ - 1 ρ p1 + σ e : grad -→ V ,( 6
)
where T is the absolute temperature, S the entropy density, p the pressure, µ k the mass chemical potentials, 1 the identity tensor and σ e the equilibrium stress tensor. d dt denotes the derivative following the barycentric velocity -→ V . At equilibrium, it is assumed that the material satisfies Hooke's law and that the liquid phase is newtonian and stokesian. By combining the balance equations of internal energy and entropy with the Gibbs relation, we deduced the entropy production of the system. The thermodynamics of linear irreversible processes makes then possible to identify the generalized forces and fluxes and to deduce the constitutive equations of the electro-active polymer. According to the Curie symmetry principle, a coupling between a force and a flux of different tensorial ranks is impossible because of the isotropy of the medium.

We thus obtained a Kelvin-Voigt type stress-strain relation ( 7) and generalized Fourier's, Darcy's (8) and Fick's laws (9). Given the orders of magnitude of the different physico-chemical parameters of the polymer (in particular, we admit that the solution is diluted), these equations can be written in the isothermal case on a first approximation (Tixier & Pouget 2016)

σ = λ trǫ 1 + 2Gǫ + λ v trǫ 1 + 2µ v ǫ , (7)
-→ V 4 - -→ V 3 ≃ - K η 2 φ 4 grad p -CF -ρ 0 2 Z 3 -→ E , (8)
-→ V 1 = -→ V 2 (9) - D C grad C - Z 1 M 1 C RT -→ E + Cv 1 RT 1 - M 1 M 2 v 2 v 1 grad p ,
where ǫ denotes the strain tensor, λ the first Lamé constant, G the shear modulus, λ v and µ v the viscoelastic coefficients; η 2 is the dynamic viscosity of the solvent, φ k the volume fractions, K the intrinsic permeability of the solid, C the cations molar concentration, F = 96487 C mol -1 the Faraday's constant, ρ 0 2 the mass density of the solvent, D the mass diffusion coefficient of the cations in the liquid phase, M k the molar masses, v k the partial molar volumes and R = 8.31J K -1 the gaz constant.

The generalized Darcy's law models the motion of the solution compared with the solid phase. This movement is caused by the pressure forces, and also by the electric field, which reflects the electro-osmosis phenomenon. The third equation expresses the motion of cations by convection ( -→ V 2 ), by mass diffusion, under the actions of the electric field and the pressure field; it can be identified with a Nernst-Planck equation [START_REF] Lakshminarayanaiah | Transport Phenomena in Membranes[END_REF]). An estimation of these different terms in the case of Nafion ® Li + shows that the pressure one is negligible.

Application of the model to a static cantilevered beam

Beam model on large displacements

To validate this model, we apply it to an EAP strip clamped at its end O under the action of a permanent electric field (static case). The other end A is either free or subject to a shear force preventing its displacement (blocking force). Consider an EAP strip of length L, of width 2l and of thickness 2e; we define a coordinate system Oxyz such that the Ox axis is along the length, the Oy axis along its width and the Oz axis parallel to the imposed electric field (see figure 2). For all numerical applications, we choose the Nafion ® Li + , an EAP well documented; in the nominal case, the dimensions of the strip are L = 2 cm, l = 2.5 mm and e = 100 µm and it is subject to an electric potential difference ϕ 0 = 1 V . Considering these values and the exerted forces, this is a two-dimensional (x, z) problem. The strip being thin and given the high values of the deflection, we used a beam model in large displacements.

Z 0 L A X - → p p --→ M p A -→ F p Z O θ θ n R - → - → Figure 2.
Forces exerted on the beam

The force system applied to the beam can be modelled by a distributed lineic force -→ p p , a bending moment --→ M p A applied to the end A and in some cases a shear force -→ F p in A. The strip being very slender, we venture the hypothesis that the distributed force is independent of the x coordinate along the beam and is orthogonal to it. The internal electrostatic forces of the strip cancel each other. The electric force produced by the electrodes also vanishes due to the electroneutrality condition. We deduce that the distributed force is zero everywhere.

L x z O A p p F p + + + + + + + + + + + + + + _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ E Figure 3. EAP bending strip
When a potential difference ϕ 0 is applied between the two faces, the cations and the solvent move towards the negative electrode, causing a volume variation and the bending of the strip (figure 3). The bending moment is therefore exerted along the Oy axis and results from the pressure forces p = -σxx

We made the usual hypotheses of a beam model: we assumed that the straight sections of the strip remain flat and normal to the neutral fibers after deformation (Bernoulli hypothesis) and that the stress and strain distributions are independent of the application points of the external forces (Barré Saint Venant hypothesis). We define a local coordinate system where -→ t and -→ n are the tangent and normal vectors; s and s denote the curvilinear abscissas along the beam at the rest and deformed configurations respectively, n the coordinate in the normal direction and θ the angle of rotation of a cross-section (figure 4). Let choose the point O as the origin of the curvilinear abscissas. No normal effort is applied, so we assumed that there is no beam elongation (ds = ds). The bending moment in the current section M p and the radius of curvature R p are

Z O θ θ t n R - → e x - → e z
M p = F p (s -L) + M p A , 1 R p = dθ ds . (11)
Let -→ u the displacement vector; its gradient with respect to the reference configuration (beam at rest) is

Grad -→ u = 1 -n R p cos θ -1 -sin θ 1 -n R p sin θ cos θ -1 . (12)
One deduces the strain tensor

ǫ = 1 2 Grad -→ u + 1 T Grad -→ u + 1 -1 . ( 13
)
The beam is thin, so |n| << R p and

ǫ xx = - n R p 1 - n 2R p ≃ - n R p . ( 14
)
In the case of a pure bending beam, the strain is ǫ xx = M p EI p n using the constitutive law for the bending moment M p = EI p R p , where E is the Young's modulus and I p = 4le 3 3 the moment of inertia with respect to the Oy axis. The shear force has a negligible effect on the deflection, so we deduce

1 R p = dθ ds = F p EI p (L -s) - M p A EI p , ( 15
) θ = F p 2EI p s (2L -s) - M p A EI p s . ( 16
)
The deflection w is obtained by integrating the relation

dz ds = sin θ. If F p = 0, the beam is circle shaped (R p is constant) w = EI p M p A cos M p A EI p L -1 , θ = - M p A EI p L . (17)
It should be noted that this deflection calculation becomes incorrect if the angle of rotation exceeds 90 • . With the hypothesis of small displacements, we would have obtained

w s = - M p A 2EI p L 2 . ( 18
)
If a blocking force is exerted, the deflection is zero which provides on small displacements

F p s = 3M p A 2L . ( 19
)
On large displacements, F p verify

w = L 0 sin F p 2EI p x(x -2L) + M p A EI p x dx = 0 . ( 20
)
Let us compare this result with the calculation on small displacements. The above condition can be written

w = 0 = L 0 sin 1 x * 2 (x + x 0 ) 2 - x 2 0 x * 2 dx = x * cos x 2 0 x * 2 [S (x 1 /x * ) -S (x 0 /x * )] -x * sin x 2 0 x * 2 [C (x 1 /x * ) -C (x 0 /x * )] , (21)
where

x * = 2EI p F p , x 0 = M p A F p -L , x 1 = M p A F p , (22)
and where S and C denote the Fresnel functions

S(x) = x 0 sin t 2 dt = +∞ n=0 (-1) n x 4n+3 (2n+1)! (4n+3) , C(x) = x 0 cos s 2 dt = +∞ n=0 (-1) n x 4n+1 (2n)! (4n+1) . (23)
In the case of Nafion

≃ F p L 2 6EI p 3 M p A F p -2L ≃ 0 . ( 24
)
As a consequence, small displacements give a good approximation of the blocking force, that we will verify with our numerical simulations. On the contrary, the calculation of the deflection of the cantilevered beam requires a large displacements model.

Additional hypothesis and static equations for a bending strip

Let us evaluate the variations of the volume fraction of the solution φ 4 . Consider a small volume dV located at a distance z from the beam axis. According to Bernoulli's hypothesis, this volume takes the value

|R p |+z
|R p | dV when the beam bends with a radius of curvature R p . The volume of the solid phase does not change, and the variation of liquid phase volume fraction is about

dφ 4 = φ 3 z |R p | = (1 -φ 4 ) M p EI p z ≃ (1 -φ 4 ) 2|w| L 2 z . ( 25
)
φ 4 ≃ 0.38 (Cappadonia et al. 1994[START_REF] Chabé | Etude des interactions moléculaires polymère-eau lors de l'hydratation de la membrane Nafion[END_REF], Nemat-Nasser & Li 2000). In the nominal case, w ∼ 1 mm (Nemat-Nasser 2002, [START_REF] Newbury | [END_REF][START_REF] Newbury | Characterization, modeling and control of ionic-polymer transducers[END_REF]. φ 4 therefore varies less than 0.1% over the thickness of the beam. As a consequence, we assume that the volume fraction φ 4 is constant.

Considering the dimensions of the strip, we assume that the problem is two-dimensional in the Oxz plane. On a first approximation, we suppose that the electric field and induction are parallel to the Oz axis:

-

→ E ≃ E z -→ e z and -→ D ≃ D z -→ e z .
We further admit that C, E z , D z , p, ρZ and the electrical potential ϕ only depend on the variable z. Finally, we neglect the pressure term of the equation ( 9), an assumption that we will verify later. The equation system of our model then becomes

E z = -dϕ dz , dDz dz = ρZ , D z = εE z , ρZ = φ 4 F (C -C moy ) , dp dz = CF -ρ 0 2 Z 3 E z , dC dz = F C RT E z , (26)
where

C moy = - (1 -φ 4 ) ρ 0 3 Z 3 φ 4 F . ( 27
)
C moy denotes the cations average concentration. The anions being attached to the polymer chains, they are uniformly distributed within the material; their concentration is therefore constant and equal to C moy considering the electroneutrality condition. From Nemat-Nasser & Li (2000), the mass density of dry Nafion is ρ 0 3 = 2 10 3 kg m -3 , and its equivalent weight, that is the weight of polymer per mole of sulfonate groups, is M eq = 1.1 kg eq -1 [START_REF] Chabé | Etude des interactions moléculaires polymère-eau lors de l'hydratation de la membrane Nafion[END_REF], Colette 2008), which provides Z 3 = -F Meq = -9 10 -4 C kg -1 . We deduce C moy = 3080 mol m -3 . We also choose an absolute temperature T = 300 K.

The boundary conditions and the electroneutrality condition are The permittivity strongly depends on the conductivity, hence of the electric charge; it increases with the water uptake [START_REF] Deng | [END_REF], Nemat-Nasser 2002), therefore with the cations concentration. We assume that it satisfies a mixing law

ε = ε 0 + αC , with α = ε moy -ε 0 C moy , (29)
where ε moy denotes the average permittivity of the material.

We have considered three models of permittivity: constant, linear and affine. We choose our permittivity values in order that deflections and blocking forces are in agreement with the literature data, namely 0.5 < w < 1.5 mm (Nemat-Nasser 2002, Newbury 2002) and 0.6 < F p < 1.3 mN [START_REF] Newbury | Characterization, modeling and control of ionic-polymer transducers[END_REF][START_REF] Newbury | [END_REF], Newbury & Leo 2003) in the nominal case (table 1).

The Nafion relative permittivity found in the literature are very scattered, but compatible with those we have chosen. The average permittivity was measured by [START_REF] Deng | [END_REF] for hydrated perfluorosulfonate ionomer membranes of the Nafion ® family with different water contents. They obtained permittivity values between 10 -7 F m -1 and 10 -6 F m -1 by electrical impedance measurements. Nemat-Nasser (2002) deduced the permittivity of capacity measurements by assimilating the IPMC strip to a capacitor; for Nafion ® 117 Li + , he got ε ≃ 2.7 10 -3 F m -1 . Farinholt & Leo (2014) deduced a close value from their measurements but did not specify their method. Wang et al. (2014) measured the permittivity of Nafion ® 117 N a + by time domain dielectric spectroscopy. For samples obtained with different manufacturing methods and a water content about 22%, the permittivity ranges from 5 10 -5 to 5 10 -3 . (30)

Resolution with different permittivity models

The system of equations and the boundary conditions become

E = -dϕ dz , D = εE , dD dz = ρZ , ε = A 0 C + A 1 , dp dz = C + A 3 E , dC dz = A 2 CE , ρZ = C -1 , D(1) = D(-1) , ϕ(-1) = 1 , ϕ(1) = 0 , (31)
with

A 0 = ϕ0(εmoy -ε 0 ) e 2 φ4F Cmoy , A 1 = ϕ0ε 0 e 2 φ4F Cmoy , A 2 = F ϕ0 RT ∼ 38.7 , A 3 = - ρ 0 2 Z3 CmoyF ∼ 0.303 , (32)
in the nominal case. We deduce the following relations

C = B 1 exp (-A 2 ϕ) , (33)
D 2 2 = 1 A 2 A 0 2 C 2 + (A 1 -A 0 )C + A 1 ϕ + B 2 2 , (34)
p = C A 2 -A 3 ϕ + B 3 , (35)
where B 1 , B 2 and B 3 are three constants. The polymer strip behaves like a conductive material. The electric field, displacement and charge are then zero throughout the strip except near the sides. We can deduce

B 2 = 1 A 2 (A 0 -2A 1 -2A 1 ln B 1 ) . (36)
The positive constant B 1 satisfies the electroneutrality condition (31)

A 0 (1 + e -A2 )B 2 1 + 2(A 1 -A 0 )B 1 - 2A 1 A 2 1 -e -A2 = 0 .( 37
)
When ϕ 0 1 V , e -A2 << 1. In the case of a constant permittivity (A 0 = 0)

B 1 = A 2 1 -e -A2 ≃ A 2 . ( 38
)
In case of a linear permittivity (A 1 = 0)

B 1 = 2 1 + e -A2 ≃ 2 , ( 39
)
and in the general case of an affine permittivity B 1 is of the order of 2

B 1 = A0-A1 A0(1+e -A 2 ) 1 + 1 + 2A0A1A2 (A0-A1) 2 1+e -A 2 1-e -A 2 ∼ 2 , (40)
assuming A 1 << A 0 and 2A2A1

A0

<< 1, assumption checked since we choose ε 0 << ε moy . We moreover verify that B 1 tends toward 1 when ϕ 0 tends toward 0.

The system of equations (31) enable us to calculate the values of the various parameters at the center and at the boundaries of the strip (table 2). Note that p is defined to within an additive constant.

-1 center 1 C B 1 exp (-A 2 ) 1 B 1 ρZ B 1 exp (-A 2 ) -1 0 B 1 -1 ϕ 1 ln B 1 A 2 0 D D(1) 0 D(1) ε A 1 A 0 + A 1 A 0 B 1 + A 1 p -B 3 B 1 exp(-A 2 ) A 2 -A 3 1 A 2 -A 3 ln B 1 A 2 B 1 A 2 where D(1) = 1 A 2 [A 0 + 2A 1 (A 2 -1 -ln B 1 )]
The bending moment is given by

M p A = A 5 M = A 5 1 -1 ( C A2 -A 3 ϕ)z dz , with M = 1 -1 p z dz , A 5 = 6le 2 F ϕ 0 C moy , 1 -1 Cz dz = 2D(1) -A0 A2 B 1 -A 1 .
(41)

I ϕ =
1 -1 ϕ z dz must be numerically evaluated. We deduce w, w s , F p , F p s and θ using ( 17) to ( 20). The resolution of the equation system is tricky from a numerical point of view because of the steepness of the functions near the boundaries. We used different methods according to the permittivity model to obtain the best precision.

In the case of a constant permittivity (A 0 = 0), we use the variable y = ln C which verifies

d 2 y dz 2 = A2
A1 (e y -1) ,

dy dz 1 ≃ dy dz -1 = 2 A2 A1 (A 2 -lnA 2 -1) , I ϕ = 1 -1 ϕ z dz = -1 A2 1 -1 yz dz . (42)
This differential equation can be numerically integrated near the boundaries, but not over the entire thickness.

In the case of a linear permittivity (A 1 = 0), the cation concentration can be calculated using the following differential equation

d 2 C dz 2 = 1 z 0 2 C -1 , ( 43
)
where z 0 = A0 A2 << 1, which provides

C = 1 + (B 1 -1) sinh (z/z 0 ) sinh (1/z 0 ) . ( 44
)
The functions thus obtained are especially steep near the boundaries in this case, which makes the calculation of I ϕ very delicate; we calculate the piecewise integral using Taylor expansions on some intervals.

In the case of an affine permittivity, the following equation can be numerically integrated over the entire interval

d 2 D dz 2 = A 2 1 + dD dz A 0 1 + dD dz + A 1 D .
(45)

Simulation results

Hypothesis validation

The dimensionless equation ( 9) provides

dC dz = A 2 C E 1 + A 4 C + A 3 , (46)
where (Tixier & Pouget 2016)

A 4 = C moy v 1 1 - M 1 M 2 v 2 v 1 ≃ C 2 moy M 1 M 2 ρ 02 2 ≃ 1.2 10 -3 ,
(M 1 = 6.9 g mol -1 and M 2 = 18 g mol -1 for Nafion ® Li + ). A 4 C + A 3 corresponds to the pressure term; its maximum value (about 0.04) is reached near the negative electrode in the case of a constant permittivity, and it is of the order of A 4 << 1 in the center of the strip. We can therefore neglect the pressure term for the calculation of the bending moment.

For each of our simulations, we calculated the different parameters in small and large displacements. The difference between blocking forces calculated in small and large displacements is less than 1% in all simulations and largely below 0.1% as soon as the thickness exceeds 50 µm. Small displacements are thus sufficiently accurate for the blocking force calculation. Concerning the deflection, the difference between both models is less than 1% when the deflection is less than 3.5 mm; it otherwise becomes very large (over 100% sometimes, especially for low thicknesses). The large displacements model is thereby essential for the deflection evaluation. Moreover, we can verify that the small displacements calculation is an upper bound on deflection.

Profiles of different quantities

We have plotted the profiles of the various quantities in the thickness of the strip in the nominal case described in section 3 (Nafion ® Li + with L = 2 cm, l = 2.5 mm, e = 100 µm and ϕ 0 = 1 V ). The profiles of dimensionless concentration, electric potential and displacement and pressure have similar aspects: the curves are almost constant in the central zone and very steep near the boundaries (figures 5 to 8).

The cation concentration profiles displays near the positive electrode a zone almost free of cations whose size depends on the permittivity model: 0.08 µm in the constant case and 0.03 µm in the affine one. The linear model predicts a permittivity tending towards 0 and is therefore incorrect in this range. Near the negative electrode, there is an accumulation of cations over a characteristic length depending on the chosen permittivity model: close to 0.1 µm for a constant permittivity and 1 µm in linear and affine cases. The concentration on the negative electrode is twenty times higher with the constant permittivity model than in the two other cases. Nemat-Nasser (2002), Wallmersperger et al. (2009) and [START_REF] Nardinocchi | [END_REF] obtained a similar profile, although less steep.

The electric potential profiles look similar to those obtained by Wallmersperger et al. (2009) and [START_REF] Nardinocchi | [END_REF], although they are less steep in the vicinity of the electrodes. The linear and affine models give almost identical results for the other profiles. The constant model distinguishes by its steepness near the boundaries, with a characteristic length close to 0.02 µm near the negative electrode, twenty times smaller than the other models; near the positive electrode, its characteristic length is 0.1 µm, five times larger than the other models for electric potential and pressure profiles. All three models provide nearly identical electrical displacement profiles with close boundary values. The pressure profiles are also very similar, except near the negative electrode where the constant model provides a boundary twenty times higher than the other two.

We furthermore verify that the electrical displacement is null in the central part of the strip with the three models, which confirms that the material behaves like a conductor.

Expected scaling laws

Let us determine the dependence of the mechanical quantities with respect to l, L, e and ϕ 0 . These quantities are written

M p A = 6F C moy le 2 ϕ 0 M , θ = -9 2 F Cmoy E L e ϕ 0 M , F p = 9F C moy le 2 L ϕ 0 M , w s = -9 4 F Cmoy E L 2 e ϕ 0 M , w = 2 9 E F Cmoy e ϕ0 1 M cos 9 2 F Cmoy E L e ϕ 0 M -1 , ( 47
)
where

F is a constant. M = 1 -1 ( C A2 -A 3 ϕ)zdz is a function of the constants A 0 , A 1 , A 2 , A 3 et B 1 . Let A 0 = a 0 ϕ 0 e 2 , A 1 = a 1 ϕ 0 e 2 , A 2 = a 2 ϕ 0 . (48)
a 0 , a 1 , a 2 and A 3 like C moy and E only depend on the chosen material.

Let us first look at the case where ϕ 0 1; B 1 checks the approximate equation (see eq. ( 37)) two values -z 1 and z 2 with z 1 , z 2 0.99. Assuming for example 0

a 0 B 1 (B 1 -2) = 2a 1 (a 2 ϕ 0 -B 1 ) , (49)
< z 1 < z 2 M = -z1 -1 p z dz + z1 -z1 p z dz + z2 z1 p z dz + 1 z2 p z dz . (50)
Given the low values of δ 1 = 1z 1 and δ 2 = 1 -z 2 , we can expand in Taylor series up to the second order the first and last integrals

-z1 -1 p z dz = -p(-1)δ 1 + o δ 2 1 , z1 -z1 p z dz = 0 , z2 z1 p z dz = p(0) 2 z 2 2 -z 1 2 , 1 z2 p z dz = p(1)δ 2 + o δ 2 2 . (51)
Hence, in all cases and at the first order in δ

M = A 3 δ 1 + B 1 A 2 δ 2 + 1 A 2 (1 -A 3 ln B 1 ) (δ 1 -δ 2 ) . (52)
δ 1 and δ 2 can be roughly evaluated using the following formulas (thanks to eq. ( 31))

dp dz -1 ∼ p (-z 1 ) -p(-1) δ 1 , dp dz 1 ∼ p(1) -p (z 2 ) δ 2 . (53)
M is therefore independent of L and l and approximately inversely proportional to e. Its dependence on ϕ 0 is more complex and linked with the chosen permittivity model. Given the values of a 0 , a 1 , a 2 , A 3 and B 1 in the nominal case, we obtain with a precision of about 20%

Constant case M ≃ a1 2 a 3 √ ϕ0 e , Linear case M ≃ a0 a 3 2 1 eϕ0 , Affine case M ≃ a0 a 3 2 f (ϕ0) eϕ0 , with f (ϕ 0 ) = a 1 a 0 a 2 2 a3ϕ 2 0 +2 a 1 a 0 a2ϕ0+1 1+2 a 1 a 0 a2ϕ0 . (54)
We deduce the scaling laws presented in table 3.

M p A - l e ϕ 3/2 0 - f (ϕ 0 ) θ L -e -2 ϕ 3/2 0 - f (ϕ 0 ) F p L -1 l e ϕ 3/2 0 - f (ϕ 0 ) ws L 2 -e -2 ϕ 3/2 0 - f (ϕ 0 )
When the imposed electric potential is very small, B 1 tends to 1 for all permittivity models, C tends to C moy and ϕ to ϕ 0 over the entire thickness of the strip. We check that all the mechanical quantities become null, which is in agreement with the experimental results.

Influence of the strip geometry

Our numerical simulations ascertain the results of the previous section with an excellent correlation: for the three permittivity models, the bending moment varies linearly with the width and thickness of the strip and is independent of its length.

The blocking force is proportional to the width and inversely proportional to the length, which is in good agreement with the results of Newbury & Leo (2003). It is also proportional to the thickness.

The predictions of the previous section are also well fitted by the deflection in large displacements; it is independent of the width and approximately proportional to L 2 : according to the chosen permittivity model, the power law that best approximates our simulations has an exponent between 1.90 and 1.96, and the variations of the ratio w/L 2 are less than 15% in all cases (figure 9). This is consistent with the results of [START_REF] Shahinpoor | Proc. of SPIE[END_REF]. It varies almost like e -2 : we find an exponent between -1.90 and -1.96 according to the chosen permittivity model and the variations of the product we 2 are less than 14% (figure 10). This result is corroborated by the measurements of He et al. (2011) as well as by the simulations of [START_REF] Vokoun | [END_REF].

We also observe that the charge of the negative electrode F C(e) = F B 1 C moy is independent of the thickness e with the three permittivity models, which agrees with the results obtained by [START_REF] Lin | [END_REF] for a close material.

Influence of the imposed electric potential

Unlike scaling laws, the relation between the different mechanical quantities and ϕ 0 depends on the chosen permittivity model. According to equations (47), the angle of rotation, the blocking force, the deflection in small displacements and the bending moment vary over the imposed potential in the same way. First, we check that the bending moment tends to 0 when ϕ 0 tends to 0 in all three cases using Taylor expansions. For imposed potentials close to 1 V , we have seen that M p

A is proportional to ϕ 3/2 0 in the case of a constant permittivity, is almost constant if the permittivity is linear and is a complex function f (ϕ 0 ) in the case of an affine permittivity (table 3, figure 11). Bakhtiarpour et al. (2016) observed experimen- Imposed electric potential (V)

Figure 11. Influence of the electric potential on the bending moment: fitting by power law (solid curve), affine laws (dashed curves) and with equations (54) (dotted curves); same colours and marks as in figure 9.

tally an approximately linear relation between the deflection and the potential for actuators, and Shahinpoor et al. (1998) and [START_REF] Mojarrad | Proc. SPIE[END_REF] did the same for the sensor effect. Our simulations show a quite good correlation with a linear law if the permittivity is constant (the variation of the ratio w/ϕ 0 is about 30%); on the contrary, the results obtained with the other two models of permittivity do not agree with these experimental results. More pre- cisely, in the constant case, the curve that best fits our numerical results is a power law of exponent 1.16 (figure 12); this curve can hardly be distinguished from an experimental linear curve [START_REF] Tixier | Chapter 39: Modelling of an Ionic Electroactive Polymer by the Thermodynamics of Linear Irreversible Processes in Generalized Models and Non Classical Mechanical Approaches in Complex Materials[END_REF]. The blocking force and the angle of rotation follow approximately the same trend (good correlation with a linear law and with a power law of exponent 1.26). This result is in good agreement with Hasani et al. (2019) experimental data, which can be well fitted by a power law of exponent 1.5. In the linear case, the moment is independent of the imposed potential for ϕ 0 1, as well as the deflection, the blocking force and the angle of rotation, which does not correspond to the experimental observations. In the affine case, the correlations with a linear law and with a power law are wrong for all the quantities. The variation of the different quantities with the imposed electric potential is thus discriminating for the permittivity model: only a constant permittivity gives results compatible with the experimental studies.

We can evaluate the average permittivity of the strip by analogy with a capacitor of thickness e and surface Ll subject to a potential difference ϕ 0 . We saw that the cations accumulate near the negative electrode over a thickness of e ′ << e, and that there was a similar thickness area without cations near the positive electrode. The electric charge density is about ρZ + = F C moy (B 1 -1) near the positive electrode and ρZ -= -F C moy near the negative one. The capacitive charge Q can therefore be estimated by

Q = (ρZ + -ρZ -)φ 4 lLe ′ = F C moy φ 4 lLe ′ B 1 . ( 55
)
The permittivity is then worth

ε = Qe Llϕ 0 = F C moy φ 4 ee ′ ϕ 0 B 1 . (56)
In the case of constant permittivity, B 1 = A 2 hence ε ∼ 4 10 -3 -4 10 -2 . In the case of a linear or affine permittivity, B 1 ≃ 2 and ε ∼ 10 -2 . These values are consistent with the permittivity values derived from capacity measurements published in the literature (Nemat-Nasser 2002, Farinholt & Leo 2004).

The hydrated IPMC can be considered as a fivelayer capacitor: the electrodes of thicknesses e 1 and e 5 , the central area of thickness e 3 very large compared to the previous thicknesses and of zero electric charge, and two very thin zones of the polymer of respective thicknesses e 2 and e 4 in the vicinity of the positive and negative electrodes. Each element has a permittivity ε i and a capacity C i = ε i L ei , and the five elements are in series. The permittivity of a material being closely related to its conductivity, so here to the cation concentration, ε 1 = ε 5 >> ε 2 >> ε 3 ε 4 . In addition e 3 is much greater than e 1 , e 2 , e 4 and e 5 . We deduce

ε = e e1+e5 ε1 + e2 ε2 + e3 ε3 + e4 ε4 ≃ ε 3 ,
where e denotes the total thickness. The overall permittivity of the strip subject to an electric field is therefore very close to the permittivity of the central part, which is that of the hydrated polymer without electric field. Our average permittivity values range from 5 10 -7 F m -1 to 10 -4 F m -1 depending on the permittivity model and are thus compatible with the permittivity values measured by [START_REF] Deng | [END_REF] and Wang et al. (2014), which range from 10 -7 F m -1 to 5 10 -3 F m -1 .

Conclusion

We studied the flexion of an ionic polymer strip using a model based on the thermomechanical of continuous media that we had previously developed. The strip, clamped at one of its ends, is subject to a continuous potential difference applied between its two faces (static case). The other end is either free or blocked by a shear force. The mechanical quantities (deflection, blocking force and angle of rotation) were determined using a beam model in large displacements.

The material chosen to perform the simulations is Nafion ® Li + . Three models of permittivity (constant, linear and affine functions of cation concentration) are examined. The permittivity values we used are in good agreement with dielectric spectroscopy and electrical impedance measurements, but significantly lower than those deduced from capacity measurements. The resolution of the equations of our model enabled us to plot the cation concentration, pressure, electric potential and displacement profiles over the thickness of the strip. These quantities are almost constant in the central part, but vary drastically in the vicinity of the electrodes, which is characteristic of a conductive material. The scaling laws obtained for the deflection and the blocking force are in good agreement with the experimental data published in the literature: in particular, the deflection varies as the square of the strip length and is inversely proportional to the square of its thickness; the blocking force is proportional to the width and the thickness and it is inversely proportional to the length. The variation of the mechanical quantities with the imposed electric potential depends on the chosen permittivity model; only the constant permittivity model provides results compatible with the experimental data and will therefore be retained for further works. We now plan to apply our model to other materials close to the Nafion ® and to study other configurations such as a strip clamped at its two ends.

Notations k = 1, 2, 3, 4, i subscripts respectively represent cations, solvent, solid, solution (water and cations) and interface; quantities without subscript refer to the whole material. Superscript 0 denotes a local quantity; the lack of superscript indicates average quantity at the macroscopic scale. Superscript T indicates the transpose of a second-rank tensor. Overlined letters denote dimensionless quantities.

Figure 1 .

 1 Figure 1. Deformable porous medium: (a) Undeformed strip (b) Strip bending under an applied electric field

σ

  xx z dz dy = 6l e -e p z dz .

Figure 4 .

 4 Figure 4. Beam on large displacements: coordinate system

  ϕ(-e) = ϕ 0 , ϕ(e) = 0 , e -e ρZ dz = 0 . (28) According to (26), this last condition is equivalent to D z (e) = D z (-e).

Figure 5 .Figure 6 .Figure 7 .Figure 8 .

 5678 Figure 5. Variation of the dimensionless cation concentration in the thickness of the strip; the distribution close to the boundaries are detailled in insets. The constant permittivity model is in blue, the linear one in red and the affine one in green.

Figure 9 .Figure 10 .

 910 Figure 9. Influence of the length: (a) on the deflection; (b) on the blocking force; (c) on the angle of rotation. The constant permittivity model is in blue disks, the linear one in red diamonds and the affine one in green triangles. Fitting by power laws (solid curves) or linear law (dashed curves).

Figure 12 .

 12 Figure12. Influence of the electric potential: (a) on the deflection; (b) on the blocking force; (c) on the angle of rotation (same colours marks and curves as in figure9).

A

  i , B i : dimensionless constants; C, C moy : cations molar concentrations (relative to the liquid phase); D: mass diffusion coefficient of the cations in the liquid phase; D: electric displacement field; e: half-thickness of the strip; E, G, λ: Young's and shear modulus, first Lamé constant; E: electric field; F = 96487 C mol -1 : Faraday's constant ; F p ( F p s ): blocking force on large (small) displacements; I: current density vector; I p : moment of inertia; K: intrinsic permeability of the solid phase; l: half-width of the strip; L: length of the strip; M k : molar mass of component k; M eq : equivalent weight (weight of polymer per mole of sulfonate groups); M p ( M p A ): bending moment; n (n): normal vector (coordinate) to the beam; p: pressure; -→ p p : distributed electric force Q: heat flux; R = 8, 314 J K -1 : gaz constant; R p : radius of curvature of the beam; s (s): curvilinear abscissa along the beam at rest (deformed) S: entropy density; T : absolute temperature; u: displacement vector; U : internal energy density; v k : partial molar volume of component k (relative to the liquid phase); V ( V k ): velocity; w (w s ): deflection of the beam on large (small) displacements; Z (Z k ): total electric charge per unit of mass; ε (ε 0 , ε moy ): permittivity (average permittivity); ǫ : strain tensor; η 2 : dynamic viscosity of water; θ: angle of rotation of a beam cross section; λ v , µ v : viscoelastic coefficients; µ k : mass chemical potential; ρ (ρ k ): mass density relative to the volume of the whole material; ρ 0 k : mass density relative to the volume of the phase; σ (σ k ), σ e : stress tensor, equilibrium stress tensor; φ k : volume fraction of phase k; ϕ (ϕ 0 ): electric potential (imposed electric potential);

Table 1 .

 1 Dielectric permittivity valuesε 0 (F m -1 ) εmoy (F m -1 )

	constant	5 10 -7	5 10 -7
	linear	0	10 -4
	affine	5 10 -7	10 -4

Table 2 .

 2 Center and boundary values of the different quantities

Table 3 .

 3 Expected scaling laws

	L	l	e	ϕ 0	ϕ 0	ϕ 0
	ε			constant	linear affine