Validation of a model for an ionic electro-active polymer in the static case
VValidation of a model for an ionic electro-activepolymer in the static case
M. Tixier & J. Pouget Laboratoire de Math´ematiques de Versailles (LMV), UMR 8100, Universit´e deVersailles Saint Quentin, 45, avenue des Etats-Unis, F-78035 Versailles, France Sorbonne Universit´e, CNRS, Institut Jean le Rond d’Alembert, UMR 7190,F-75005 Paris, FranceE-mail: [email protected], [email protected]
April 2020
Abstract.
IPMCs consist of a Nafion ® ionic polymer film coated on both sideswith a thin layer of metallic electrodes. The polymer completely dissociates whenit is saturated with water, releasing small cations while anions remain boundto the polymer chains. When this strip is subject to an orthogonal electricfield, the cations migrate towards the negative electrode, carrying water awayby osmosis. This leads to the bending of the strip. We have previously publisheda modelling of this system based on the thermodynamics of irreversible processes.In this paper, we use this model to simulate numerically the bending of a strip.Since the amplitude of the deflection is large, we use a beam model in largedisplacements. In addition, the material permittivity may increase with ionconcentration. We therefore test three permittivity models. We plot the profiles ofthe cations concentration, pressure, electric potential and induction, and we studythe influence of the strip geometry on the tip displacement and on the blockingforce. The results we obtain are in good agreement with the experimental datapublished in the literature. The variation of these quantities with the imposedelectric potential allow us to discriminate between the three models. 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.-jpublished in
Smart Materials and Structures (085019), 2020https://doi.org/10.1088/1361-665X/ab8fca a r X i v : . [ phy s i c s . c l a ss - ph ] J a n alidation of a model for an ionic electro-active polymer in the static case
1. Introduction
The present work proposes a study of a thin stripof electro-active polymer (EAP). More precisely, thestudy focuses on Ionic Polymer-Metal Composite(IPMC) belonging to the ionic class. We investigateresponses of the EAP subject to an applied differenceof electric potential on the metallic electrodes andan applied punctual force at the tip of the strip.In previous works we have constructed step bystep a micro-mechanical model to establish theconservation laws for electro-active polymers (Tixier& Pouget 2014). Following this work, the constitutiveequations were deduced from the hypothesis of localthermodynamics equilibrium and the Gibbs relationusing the thermodynamics of irreversible processes(Tixier & Pouget 2016). The present study follows onfrom the spirit of the previous works and we want tocharacterize the behavior of a thin blade, especially thechemical, electrical and mechanical parameters underelectro-mechanical loading.The behavior of electro-chemical-mechanical interac-tions of EAP is of great interest for research and engi-neering advanced technology. Simply, we can say thatan EAP is a polymer exhibiting a mechanical response,such a stretching, contracting or bending for example,when subject to an electric field. Conversely the EAPcan produce electric energy in response to a mechani-cal load (Shahinpoor et al. 1998, Shahinpoor & Kim2001, Bar-Cohen 2005, Pugal et al. 2010, Park et al.2010). This particular property is highly attractive forapplications. For instance, we can quote biomimeticdevices (robotics, bio-inspired underwater robots suchas fishes (Chen 2017), haptic actuators (artificial skin,tactile displays or artificial muscle (Deole et al. 2008,Matysek et al. 2009, Bar-Cohen 2005) and this ma-terial is an excellent candidate for energy harvesting(Farinholt et al. 2009, Tiwari et al. 2008, Aureli etal. 2010, Jean-Mistral et al. 2010, Cha et al. 2013).In addition promising applications to micro-mechanicalsystems (MEMS) at the sub-micron scale are also ex-pected in medical engineering for accurate medical con-trol or investigation, for instance Fang et al. (2007) andChikhaoui et al. (2018).Electro-active polymers can be categorized intotwo main groups depending on their activationmechanisms. The first one is electronic electro-activepolymers and dielectrics which are subject to Coulombforce ; their volume change is due to the applicationof an electric field. For instance dielectric elastomersbelongs to this category as well as piezoelectric andelectrostrictive polymers. The second group of EAP isthe ionic electro-active polymers. The latter are drivenby the displacement of ions inside the material. Thepolyelectrolyte gels, ionic polymer-metal composites, conductive polymers are such media. One of theadvantages of the these polymers is that they can beactivated by very low difference of electric potential ofabout 1 − ® Li + . The authors deduce the profiles,within the strip thickness, of the electric charge density,electric potential, electric field and strain. Nardinoc-chi et al. (2011) deduce a model based on the 3-Dtheory of linear elasticity. The thermodynamics allowsthem to introduce chemo-electro-mechanical couplingand they deduce the constitutive equations of the ma-terial, especially a Nernst-Planck like equation is de-duced. Their study continues with numerical illustra-tions for a thin strip of EAP. Moreover, time evolutionat low frequency is proposed as well. It is notewor-thy that these papers are partly, in their spirit, ratherclose to the present model. Nevertheless, their results alidation of a model for an ionic electro-active polymer in the static case
2. Modelling of the polymer
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 withthin metal layers acting as electrodes. The modelwe developed applies for example to Nafion ® , anionic polymer well documented in the literature thatwe will use for the validation of our constitutiveequations. When saturated with water, this polymerdissociates quasi completely, releasing small cations inwater whereas anions remain bound to the polymerbackbone (Chabe 2008). When a potential differenceis applied between the two electrodes, the cationsmigrate towards the negative electrode, carrying thewater away by osmosis. As a result, the blade contractson the side of the positive electrode and swells on theopposite side, causing its bending (Figure 1). - + - + ++ + + + ++++ + ++ ____ ____ _____ __ __ + ++ +++++ - --- - - - + +++ + ++ +++ - -- -- -- -- -- Figure 1.
Deformable porous medium: (a) Undeformed strip(b) Strip bending under an applied electric field
To model the electro-active polymer, we used a”continuous medium” approach. Negatively chargedpolymer chains are assimilated to a deformable,homogeneous and isotropic porous medium in whichflows an ionic solution (water and cations). Thesystem is therefore composed of two phases and three alidation of a model for an ionic electro-active polymer in the static case i the interfaces. Gravity and magnetic inductionare supposed to be negligible; the only external forceexerted is therefore the electric action. The differentphases are supposed to be incompressible and thesolution diluted. It is further recognized that soliddeformations are small. We used a coarse-grained model developed for two-component mixtures (Ishii & Hibiki 2006). We definetwo scales. The conservation equations are first writtenat the microscopic level for each phase and for theinterfaces. At this scale (typically about 100 A ◦ ),the elementary volume contains one phase only but islarge enough so that the medium can be consideredas continuous. The macroscopic equations of thematerial are deduced by averaging the microscopicones using a presence function for each phase andinterface. For each phase k and for the interfaces,a Heaviside-like function of presence is defined. Themacroscale quantities are obtained by averaging thecorresponding microscale quantities weighted by thefunctions of presence. This volume average isassumed to be equivalent to a statistical average(ergodic hypothesis). At the macroscopic scale, therepresentative elementary volume (R.E.V.) must belarge enough so that these averages are relevant, butsmall enough so that the average quantities could beconsidered as local. According to Gierke et al. (1981)and Chabe (2008), its characteristic length is about1 µm .To write the balance equations, it is necessaryto calculate the variations of the extensive quantitiesfor a closed system in the thermodynamic sense. Atthe microscopic scale, we use the particle derivativeor derivative following the motion of a constituentor an interface. At the macroscopic scale, thethree constituents velocities are different; we introducea ”material derivative” or derivative following themovement of the three constituants DDt , which is aweighted average of the particle derivatives related toeach constituent (Coussy 1995, Biot 1977).
We thus obtain balance equations of mass (1),momentum (2), total, kinetic, potential and internalenergy densities (3), entropy, electric charge (4) and the Maxwell equations (5) for the complete material(Tixier & Pouget 2014) ∂ρ∂t + div (cid:16) ρ −→ V (cid:17) = 0 , (1) ρ D −→ VDt = div σ (cid:101) + ρZ −→ E , (2) ρ DDt (cid:18) Uρ (cid:19) = (cid:88) , (cid:16) σ k (cid:102) : grad (cid:126)V k (cid:17) + (cid:126)I − (cid:88) k =3 , ,i (cid:16) ρ k Z k (cid:126)V k (cid:17) · (cid:126)E − div (cid:126)Q , (3)div −→ I + ∂ρZ∂t = 0 , (4)rot −→ E = −→ , div −→ D = ρZ , −→ D = ε −→ E , (5)where ρ k denotes the densities relative to the volumeof the whole material, −→ V the velocity, σ (cid:101) the stresstensor, 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 theelectrical induction and ε the dielectric permittivity.Subscripts refer to a phase, interface or a constituentand quantities without subscript are relative to thewhole material.We verify that the stress tensor of the wholematerial is symmetrical. The second member ofequation (3) highlights the source terms of the internalenergy (viscous dissipation and Joule heating) and itsflux −→ Q . Equations (4) and (5) show that an EAPbehaves like an isotropic homogeneous linear dielectric.In the last Maxwell equation, the permittivity of thewhole material is obtained by a local mixing law, andis therefore likely to vary over space and time. We make the hypothesis of local thermodynamicequilibrium. The Gibbs relation of the whole materialis deduced from the Gibbs relations introduced by deGroot et Mazur (1962) for a deformable solid and fora two-constituent fluid
T ddt (cid:18) Sρ (cid:19) = ddt (cid:18) Uρ (cid:19) + p ddt (cid:18) ρ (cid:19) − (cid:88) k =1 , , µ k ddt (cid:18) ρ k ρ (cid:19) − ρ (cid:16) p (cid:101) + σ e (cid:102)(cid:17) : grad −→ V , (6)where T is the absolute temperature, S the entropydensity, p the pressure, µ k the mass chemicalpotentials, 1 (cid:101) the identity tensor and σ e (cid:102) the equilibriumstress tensor. ddt denotes the derivative following thebarycentric velocity −→ V . At equilibrium, it is assumed alidation of a model for an ionic electro-active polymer in the static case σ (cid:101) = λ (cid:0) tr(cid:15) (cid:101)(cid:1) (cid:101) + 2 G(cid:15) (cid:101) + λ v (cid:0) tr ˙ (cid:15) (cid:101)(cid:1) (cid:101) + 2 µ v ˙ (cid:15) (cid:101) , (7) −→ V − −→ V (cid:39) − Kη φ (cid:104) grad p − (cid:0) CF − ρ Z (cid:1) −→ E (cid:105) , (8) −→ V = −→ V (9) − D C (cid:20) grad C − Z M CRT −→ E + Cv RT (cid:18) − M M v v (cid:19) grad p (cid:21) , where (cid:15) (cid:101) denotes the strain tensor, λ the first Lam´econstant, G the shear modulus, λ v and µ v theviscoelastic coefficients; η is the dynamic viscosityof the solvent, φ k the volume fractions, K theintrinsic permeability of the solid, C the cations molarconcentration, F = 96487 C mol − the Faraday’sconstant, ρ the mass density of the solvent, D the massdiffusion coefficient of the cations in the liquid phase, M k the molar masses, v k the partial molar volumes and R = 8 . J K − the gaz constant.The generalized Darcy’s law models the motionof the solution compared with the solid phase. Thismovement is caused by the pressure forces, and alsoby the electric field, which reflects the electro-osmosisphenomenon. The third equation expresses the motionof cations by convection ( −→ V ), by mass diffusion, underthe actions of the electric field and the pressure field;it can be identified with a Nernst-Planck equation(Lakshminarayanaiah 1969). An estimation of thesedifferent terms in the case of Nafion ® Li + shows thatthe pressure one is negligible.
3. Application of the model to a staticcantilevered beam
To validate this model, we apply it to an EAP stripclamped at its end O under the action of a permanent electric field (static case). The other end A is either freeor subject to a shear force preventing its displacement(blocking force). Consider an EAP strip of length L ,of width 2 l and of thickness 2 e ; we define a coordinatesystem Oxyz such that the Ox axis is along the length,the Oy axis along its width and the Oz axis parallelto the imposed electric field (see figure 2). For allnumerical applications, we choose the Nafion ® Li + ,an EAP well documented; in the nominal case, thedimensions of the strip are L = 2 cm , l = 2 . mm and e = 100 µm and it is subject to an electricpotential difference ϕ = 1 V . Considering these valuesand the exerted forces, this is a two-dimensional ( x, z )problem. The strip being thin and given the highvalues of the deflection, we used a beam model in largedisplacements. + + + + + + + + + + + + + + _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ (cid:31)E Z0 L A X −→ p p −−→ M pA −→ F p Z O θ θ (cid:31)t(cid:31)n R −→ e x −→ e z Figure 2.
Forces exerted on the beam
The force system applied to the beam can bemodelled by a distributed lineic force −→ p p , a bendingmoment −−→ M pA applied to the end A and in some casesa shear force −→ F p in A . The strip being very slender,we venture the hypothesis that the distributed force isindependent of the x coordinate along the beam and isorthogonal to it. The internal electrostatic forces of thestrip cancel each other. The electric force produced bythe electrodes also vanishes due to the electroneutralitycondition. We deduce that the distributed force is zeroeverywhere. L x z O A p p M pA F p + + + + + + + + + + + + + + _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ E z R O t n e x e z Figure 3.
EAP bending strip
When a potential difference ϕ is applied betweenthe two faces, the cations and the solvent move towardsthe negative electrode, causing a volume variation andthe bending of the strip (figure 3). The bendingmoment is therefore exerted along the Oy axis andresults from the pressure forces p = − σ xx M pA = − (cid:90) l − l (cid:90) e − e σ xx z dz dy = 6 l (cid:90) e − e p z dz . (10) alidation of a model for an ionic electro-active polymer in the static case −→ t and −→ n are the tangent and normalvectors; s and s denote the curvilinear abscissas alongthe beam at the rest and deformed configurationsrespectively, n the coordinate in the normal directionand θ the angle of rotation of a cross-section (figure 4).Let choose the point O as the origin of the curvilinearabscissas. No normal effort is applied, so we assumedthat there is no beam elongation ( ds = ds ). + + + + + + + + + + + + + + _ _ _ _ _ _ _ _ _ _ _ _ _ _ _ (cid:31)E Z0 L A X −→ p p −−→ M pA −→ F p Z O θ θ (cid:31)t(cid:31)n R −→ e x −→ e z Figure 4.
Beam on large displacements: coordinate system
The bending moment in the current section M p and the radius of curvature R p are M p = F p ( s − L ) + M pA , R p = dθds . (11)Let −→ u the displacement vector; its gradient withrespect to the reference configuration (beam at rest)isGrad −→ u = (cid:18) (cid:0) − nR p (cid:1) cos θ − − sin θ (cid:0) − nR p (cid:1) sin θ cos θ − (cid:19) . (12)One deduces the strain tensor (cid:15) (cid:101) = 12 (cid:104)(cid:0) Grad −→ u + 1 (cid:101)(cid:1) T (cid:0) Grad −→ u + 1 (cid:101)(cid:1) − (cid:101)(cid:105) . (13)The beam is thin, so | n | << R p and (cid:15) xx = − nR p (cid:18) − n R p (cid:19) (cid:39) − nR p . (14)In the case of a pure bending beam, the strain is (cid:15) xx = M p EI p n using the constitutive law for the bendingmoment M p = EI p R p , where E is the Young’s modulusand I p = le the moment of inertia with respect tothe Oy axis. The shear force has a negligible effect onthe deflection, so we deduce1 R p = dθds = F p EI p ( L − s ) − M pA EI p , (15) θ = F p EI p s (2 L − s ) − M pA EI p s . (16)The deflection w is obtained by integrating the relation dzds = sin θ .If F p = 0, the beam is circle shaped ( R p isconstant) w = EI p M pA (cid:104) cos (cid:16) M pA EI p L (cid:17) − (cid:105) , θ = − M pA EI p L . (17)It should be noted that this deflection calculationbecomes incorrect if the angle of rotation exceeds 90 ◦ .With the hypothesis of small displacements, we wouldhave obtained w s = − M pA EI p L . (18)If a blocking force is exerted, the deflection is zerowhich provides on small displacements F ps = 3 M pA L . (19)On large displacements, F p verify w = (cid:90) L sin (cid:20) F p EI p x ( x − L ) + M pA EI p x (cid:21) dx = 0 . (20)Let us compare this result with the calculationon small displacements. The above condition can bewritten w = 0 = (cid:82) L sin (cid:104) x ∗ ( x + x ) − x x ∗ (cid:105) dx = x ∗ cos (cid:16) x x ∗ (cid:17) [ S ( x /x ∗ ) − S ( x /x ∗ )] − x ∗ sin (cid:16) x x ∗ (cid:17) [ C ( x /x ∗ ) − C ( x /x ∗ )] , (21)where x ∗ = (cid:113) EI p F p , x = M pA F p − L , x = M pA F p , (22)and where S and C denote the Fresnel functions S ( x ) = (cid:82) x sin t dt = + ∞ (cid:80) n =0 ( − n x n +3 (2 n +1)! (4 n +3) ,C ( x ) = (cid:82) x cos s dt = + ∞ (cid:80) n =0 ( − n x n +1 (2 n )! (4 n +1) . (23)In the case of Nafion ® Li + , E = 1 , P a (Baueret al. 2005, Barclay Satterfield & Benziger 2009,Silberstein & Boyce 2010). Using the blocking forcevalues provided by the literature in the nominal case(Newbury 2002, Newbury & Leo 2002, Newbury & Leo2003), x x ∗ < x x ∗ ∼ .
45. We deduce, with a relativeerror around 2% w (cid:39) F p L EI p (cid:18) M pA F p − L (cid:19) (cid:39) . (24)As a consequence, small displacements give a goodapproximation of the blocking force, that we will verifywith our numerical simulations. On the contrary, thecalculation of the deflection of the cantilevered beamrequires a large displacements model. alidation of a model for an ionic electro-active polymer in the static case Let us evaluate the variations of the volume fraction ofthe solution φ . Consider a small volume dV locatedat a distance z from the beam axis. According toBernoulli’s hypothesis, this volume takes the value | R p | + z | R p | dV when the beam bends with a radius ofcurvature R p . The volume of the solid phase doesnot change, and the variation of liquid phase volumefraction is about dφ = φ z | R p | = (1 − φ ) M p EI p z (cid:39) (1 − φ ) 2 | w | L z . (25) φ (cid:39) .
38 (Cappadonia et al. 1994, Chab´e 2008,Nemat-Nasser & Li 2000). In the nominal case, w ∼ mm (Nemat-Nasser 2002, Newbury & Leo 2002,Newbury 2002). φ therefore varies less than 0 . φ is constant.Considering the dimensions of the strip, weassume that the problem is two-dimensional in the Oxz plane. On a first approximation, we suppose thatthe electric field and induction are parallel to the Oz axis: −→ E (cid:39) E z −→ e z and −→ D (cid:39) D z −→ e z . We further admitthat C , E z , D z , p , ρZ and the electrical potential ϕ only depend on the variable z . Finally, we neglect thepressure term of the equation (9), an assumption thatwe will verify later. The equation system of our modelthen becomes E z = − dϕdz , dD z dz = ρZ ,D z = εE z , ρZ = φ F ( C − C moy ) , dpdz = (cid:0) CF − ρ Z (cid:1) E z , dCdz = F CRT E z , (26)where C moy = − (1 − φ ) ρ Z φ F . (27) C moy denotes the cations average concentration. Theanions being attached to the polymer chains, theyare uniformly distributed within the material; theirconcentration is therefore constant and equal to C moy considering the electroneutrality condition. FromNemat-Nasser & Li (2000), the mass density of dryNafion is ρ = 2 10 kg m − , and its equivalent weight,that is the weight of polymer per mole of sulfonategroups, is M eq = 1 . kg eq − (Chab´e 2008, Colette2008), which provides Z = − FM eq = − − C kg − .We deduce C moy = 3080 mol m − . We also choose anabsolute temperature T = 300 K .The boundary conditions and the electroneutralitycondition are ϕ ( − e ) = ϕ , ϕ ( e ) = 0 , (cid:82) e − e ρZ d z = 0 . (28)According to (26), this last condition is equivalent to D z ( e ) = D z ( − e ). Table 1.
Dielectric permittivity values ε ( F m − ) ε moy ( F m − )constant 5 10 − − linear 0 10 − affine 5 10 − − The permittivity strongly depends on the conduc-tivity, hence of the electric charge; it increases with thewater uptake (Deng & Mauritz 1992, Nemat-Nasser2002), therefore with the cations concentration. Weassume that it satisfies a mixing law ε = ε + αC , with α = ε moy − ε C moy , (29)where ε moy denotes the average permittivity of thematerial. We have considered three models ofpermittivity: constant, linear and affine. We chooseour permittivity values in order that deflections andblocking forces are in agreement with the literaturedata, namely 0 . < w < . mm (Nemat-Nasser 2002,Newbury 2002) and 0 . < F p < . mN (Newbury2002, Newbury & Leo 2002, Newbury & Leo 2003) inthe nominal case (table 1).The Nafion relative permittivity found in theliterature are very scattered, but compatible withthose we have chosen. The average permittivitywas measured by Deng & Mauritz (1992) forhydrated perfluorosulfonate ionomer membranes of theNafion ® family with different water contents. Theyobtained permittivity values between 10 − F m − and10 − F m − by electrical impedance measurements.Nemat-Nasser (2002) deduced the permittivity ofcapacity measurements by assimilating the IPMC stripto a capacitor; for Nafion ® Li + , he got ε (cid:39) . − F m − . Farinholt & Leo (2014) deduceda close value from their measurements but did notspecify their method. Wang et al. (2014) measuredthe permittivity of Nafion ® N a + by time domaindielectric spectroscopy. For samples obtained withdifferent manufacturing methods and a water contentabout 22%, the permittivity ranges from 5 10 − to5 10 − . Let us introduce dimensionless variables E = E z eϕ , ϕ = ϕϕ , D = D z eφ F C moy ,C = CC moy , ρZ = ρZφ F C moy , p = pF ϕ C moy ,z = ze , ε = εϕ e φ F C moy . (30)The system of equations and the boundary conditions alidation of a model for an ionic electro-active polymer in the static case E = − dϕdz , D = εE , dDdz = ρZ , ε = A C + A , dpdz = (cid:0) C + A (cid:1) E , dCdz = A CE ,ρZ = C − , D (1) = D ( − ,ϕ ( −
1) = 1 , ϕ (1) = 0 , (31)with A = ϕ ( ε moy − ε ) e φ F C moy , A = ϕ ε e φ F C moy ,A = F ϕ RT ∼ . , A = − ρ Z C moy F ∼ . , (32)in the nominal case. We deduce the following relations C = B exp ( − A ϕ ) , (33) D A (cid:18) A C + ( A − A ) C (cid:19) + A ϕ + B , (34) p = CA − A ϕ + B , (35)where B , B and B are three constants. Thepolymer strip behaves like a conductive material. Theelectric field, displacement and charge are then zerothroughout the strip except near the sides. We candeduce B = 1 A ( A − A − A ln B ) . (36)The positive constant B satisfies the electroneu-trality condition (31) A (1 + e − A ) B + 2( A − A ) B − A A − e − A = 0 . (37)When ϕ (cid:38) V , e − A <<
1. In the case of a constantpermittivity ( A = 0) B = A − e − A (cid:39) A . (38)In case of a linear permittivity ( A = 0) B = 21 + e − A (cid:39) , (39)and in the general case of an affine permittivity B isof the order of 2 B = A − A A (1+ e − A ) (cid:104) (cid:113) A A A ( A − A ) e − A − e − A (cid:105) ∼ , (40)assuming A << A and A A A <<
1, assumptionchecked since we choose ε << ε moy . We moreoververify that B tends toward 1 when ϕ tends toward 0.The system of equations (31) enable us to calculatethe values of the various parameters at the center andat the boundaries of the strip (table 2). Note that p isdefined to within an additive constant. Table 2.
Center and boundary values of the different quantities − C B exp ( − A ) 1 B ρZ B exp ( − A ) − B − ϕ ln B A D D (1) 0 D (1) ε A A + A A B + A p − B B exp( − A ) A − A A − A ln B A B A where D (1) = (cid:113) A [ A + 2 A ( A − − ln B )] The bending moment is given by M pA = A M = A (cid:82) − ( CA − A ϕ ) z dz , with M = (cid:82) − p z dz ,A = 6 le F ϕ C moy , (cid:82) − Cz dz = 2 D (1) − A A B − A . (41) I ϕ = (cid:82) − ϕ z dz must be numerically evaluated. Wededuce w , w s , F p , F ps and θ using (17) to (20).The resolution of the equation system is trickyfrom a numerical point of view because of the steepnessof the functions near the boundaries. We used differentmethods according to the permittivity model to obtainthe best precision.In the case of a constant permittivity ( A = 0),we use the variable y = ln C which verifies d ydz = A A ( e y − , dydz (cid:12)(cid:12)(cid:12) (cid:39) dydz (cid:12)(cid:12)(cid:12) − = (cid:113) A A ( A − lnA − ,I ϕ = (cid:82) − ϕ z dz = − A (cid:82) − yz dz . (42)This differential equation can be numerically inte-grated near the boundaries, but not over the entirethickness.In the case of a linear permittivity ( A = 0),the cation concentration can be calculated using thefollowing differential equation d Cdz = 1 z (cid:0) C − (cid:1) , (43)where z = (cid:113) A A <<
1, which provides C = 1 + ( B −
1) sinh ( z/z )sinh (1 /z ) . (44)The functions thus obtained are especially steepnear the boundaries in this case, which makes thecalculation of I ϕ very delicate; we calculate thepiecewise integral using Taylor expansions on someintervals. alidation of a model for an ionic electro-active polymer in the static case d Ddz = A (cid:16) dDdz (cid:17) A (cid:16) dDdz (cid:17) + A D . (45)
4. Simulation results
The dimensionless equation (9) provides dCdz = A C E (cid:2) A (cid:0) C + A (cid:1)(cid:3) , (46)where (Tixier & Pouget 2016) A = C moy v (cid:18) − M M v v (cid:19) (cid:39) C moy M M ρ (cid:39) . − , ( M = 6 . g mol − and M = 18 g mol − for Nafion ® Li + ). A (cid:0) C + A (cid:1) corresponds to thepressure term; its maximum value (about 0 .
04) isreached near the negative electrode in the case of aconstant permittivity, and it is of the order of A << .
1% as soon as thethickness exceeds 50 µm . Small displacements are thussufficiently accurate for the blocking force calculation.Concerning the deflection, the difference between bothmodels is less than 1% when the deflection is lessthan 3 . mm ; it otherwise becomes very large (over100% sometimes, especially for low thicknesses). Thelarge displacements model is thereby essential for thedeflection evaluation. Moreover, we can verify that thesmall displacements calculation is an upper bound ondeflection. We have plotted the profiles of the various quantitiesin the thickness of the strip in the nominal casedescribed in section 3 (Nafion ® Li + with L = 2 cm , l = 2 . mm , e = 100 µm and ϕ = 1 V ). The profilesof dimensionless concentration, electric potential anddisplacement and pressure have similar aspects: thecurves are almost constant in the central zone and verysteep near the boundaries (figures 5 to 8).The cation concentration profiles displays near thepositive electrode a zone almost free of cations whosesize depends on the permittivity model: 0 . µm in the constant case and 0 . µm in the affine one. Thelinear model predicts a permittivity tending towards0 and is therefore incorrect in this range. Near thenegative electrode, there is an accumulation of cationsover a characteristic length depending on the chosenpermittivity model: close to 0 . µm for a constantpermittivity and 1 µm in linear and affine cases.The concentration on the negative electrode is twentytimes higher with the constant permittivity modelthan in the two other cases. Nemat-Nasser (2002),Wallmersperger et al. (2009) and Nardinocchi et al.(2011) obtained a similar profile, although less steep.The electric potential profiles look similar tothose obtained by Wallmersperger et al. (2009) andNardinocchi et al. (2011), although they are lesssteep in the vicinity of the electrodes. The linearand affine models give almost identical results for theother profiles. The constant model distinguishes byits steepness near the boundaries, with a characteristiclength close to 0 . µm near the negative electrode,twenty times smaller than the other models; near thepositive electrode, its characteristic length is 0 . µm ,five times larger than the other models for electricpotential and pressure profiles. All three modelsprovide nearly identical electrical displacement profileswith close boundary values. The pressure profiles arealso very similar, except near the negative electrodewhere the constant model provides a boundary twentytimes higher than the other two.We furthermore verify that the electrical displace-ment is null in the central part of the strip with thethree models, which confirms that the material behaveslike a conductor. Let us determine the dependence of the mechanicalquantities with respect to l , L , e and ϕ . Thesequantities are written M pA = 6 F C moy le ϕ M ,θ = − F C moy
E Le ϕ M ,F p = 9 F C moy le L ϕ M ,w s = − F C moy
E L e ϕ M ,w = EF C moy eϕ M (cid:104) cos (cid:16) F C moy
E Le ϕ M (cid:17) − (cid:105) , (47)where F is a constant. M = (cid:82) − ( CA − A ϕ ) zdz is afunction of the constants A , A , A , A et B . Let A = a ϕ e , A = a ϕ e , A = a ϕ . (48) a , a , a and A like C moy and E only depend on thechosen material.Let us first look at the case where ϕ (cid:38) B checks the approximate equation (see eq. (37)) a B ( B −
2) = 2 a ( a ϕ − B ) , (49) alidation of a model for an ionic electro-active polymer in the static case Figure 5.
Variation of the dimensionless cation concentration in the thickness of the strip; the distribution close to the boundariesare detailled in insets. The constant permittivity model is in blue, the linear one in red and the affine one in green. D i m e n s i o n l e ss e l ec t r i c p o t e n t i a l Figure 6.
Variation of the dimensionless electric potential in the thickness of the strip; the distribution close to the boundaries aredetailled in insets (same colours as in figure 5). D i m e n s i o n l e ss e l ec t r i c d i s p l a ce m e n t Figure 7.
Variation of the dimensionless electric displacement in the thickness of the strip; the distribution close to the boundariesare detailled in insets (same colours as in figure 5). and therefore depends only on ϕ and the material.The pressure profile is almost constant throughout the center of the strip and is very steep near theboundaries. It can be modelled by a constant between alidation of a model for an ionic electro-active polymer in the static case -0.4-0.200.20.40.60.81 -1 -0.8 -0.6 -0.4 -0.2 0 0.2 0.4 0.6 0.8 1 -0.10.10.30.995 0.9975 1 -0.4-0.20 -1 -0.9975 -0.995 Figure 8.
Variation of the dimensionless pressure in the thickness of the strip; the distribution close to the boundaries are detailledin insets (same colours as in figure 5). two values − z and z with z , z (cid:46) .
99. Assumingfor example 0 < z < z M = (cid:82) − z − p z dz + (cid:82) z − z p z dz + (cid:82) z z p z dz + (cid:82) z p z dz . (50)Given the low values of δ = 1 − z and δ = 1 − z ,we can expand in Taylor series up to the second orderthe first and last integrals (cid:82) − z − p z dz = − p ( − δ + o (cid:0) δ (cid:1) , (cid:82) z − z p z dz = 0 , (cid:82) z z p z dz = p (0)2 (cid:0) z − z (cid:1) , (cid:82) z p z dz = p (1) δ + o (cid:0) δ (cid:1) . (51)Hence, in all cases and at the first order in δM = A δ + B A δ + 1 A (1 − A ln B ) ( δ − δ ) . (52) δ and δ can be roughly evaluated using the followingformulas (thanks to eq. (31)) dpdz (cid:12)(cid:12)(cid:12)(cid:12) − ∼ p ( − z ) − p ( − δ , dpdz (cid:12)(cid:12)(cid:12)(cid:12) ∼ p (1) − p ( z ) δ . (53) M is therefore independent of L and l and approxi-mately inversely proportional to e . Its dependence on ϕ is more complex and linked with the chosen permit-tivity model. Given the values of a , a , a , A and B in the nominal case, we obtain with a precision ofabout 20%Constant case M (cid:39) (cid:112) a a √ ϕ e , Linear case M (cid:39) (cid:113) a a eϕ , Affine case M (cid:39) (cid:113) a a f ( ϕ ) eϕ , with f ( ϕ ) = a a a a ϕ +2 a a a ϕ +1 (cid:113) a a a ϕ . (54)We deduce the scaling laws presented in table 3. Table 3.
Expected scaling laws
L l e ϕ ϕ ϕ ε constant linear affine M pA - l e ϕ / - f ( ϕ ) θ L - e − ϕ / - f ( ϕ ) F p L − l e ϕ / - f ( ϕ ) w s L - e − ϕ / - f ( ϕ ) When the imposed electric potential is very small, B tends to 1 for all permittivity models, C tends to C moy and ϕ to ϕ over the entire thickness of the strip.We check that all the mechanical quantities becomenull, which is in agreement with the experimentalresults. Our numerical simulations ascertain the results of theprevious section with an excellent correlation: for thethree permittivity models, the bending moment varieslinearly with the width and thickness of the strip andis independent of its length.The blocking force is proportional to the widthand inversely proportional to the length, which is ingood agreement with the results of Newbury & Leo(2003). It is also proportional to the thickness.The predictions of the previous section are alsowell fitted by the deflection in large displacements; itis independent of the width and approximately pro-portional to L : according to the chosen permittivitymodel, the power law that best approximates our sim-ulations has an exponent between 1 .
90 and 1 .
96, andthe variations of the ratio w/L are less than 15% inall cases (figure 9). This is consistent with the resultsof Shahinpoor (1999). It varies almost like e − : we alidation of a model for an ionic electro-active polymer in the static case Length (mm) D e fl ec t i o n ( mm ) B l o c k i n g f o r ce ( m N ) Length (mm) ( b ) ( c ) Length (mm) A n g l e o f r o t a t i o n ( d e g ) Figure 9.
Influence of the length: (a) on the deflection; (b) on the blocking force; (c) on the angle of rotation. The constantpermittivity model is in blue disks, the linear one in red diamonds and the affine one in green triangles. Fitting by power laws (solidcurves) or linear law (dashed curves). D e fl ec t i o n ( mm ) Thickness ( µ m) ( a ) B l o c k i n g f o r ce ( m N ) Thickness ( µ m) ( b ) Thickness ( µ m) A n g l e o f r o t a t i o n ( d e g ) ( c ) Figure 10.
Influence of the thickness: (a) on the deflection; (b) on the blocking force; (c) on the angle of rotation (same colours,marks and curves as in figure 9). find an exponent between − .
90 and − .
96 accordingto the chosen permittivity model and the variations ofthe product we are less than 14% (figure 10). Thisresult is corroborated by the measurements of He etal. (2011) as well as by the simulations of Vokoun etal. (2015).We also observe that the charge of the negativeelectrode F C ( e ) = F B C moy is independent of thethickness e with the three permittivity models, whichagrees with the results obtained by Lin et al. (2012)for a close material. Unlike scaling laws, the relation between the differentmechanical quantities and ϕ depends on the chosenpermittivity model. According to equations (47), theangle of rotation, the blocking force, the deflection insmall displacements and the bending moment vary overthe imposed potential in the same way. First, we checkthat the bending moment tends to 0 when ϕ tendsto 0 in all three cases using Taylor expansions. Forimposed potentials close to 1 V , we have seen that M pA is proportional to ϕ / in the case of a constantpermittivity, is almost constant if the permittivity islinear and is a complex function f ( ϕ ) in the case ofan affine permittivity (table 3, figure 11).Bakhtiarpour et al. (2016) observed experimen- B e nd i n g m o m e n t( µ N m ) Imposed electric potential (V)
Figure 11.
Influence of the electric potential on the bendingmoment: fitting by power law (solid curve), affine laws (dashedcurves) and with equations (54) (dotted curves); same coloursand marks as in figure 9. tally an approximately linear relation between the de-flection and the potential for actuators, and Shahin-poor et al. (1998) and Mojarrad & Shahinpoor (1997)did the same for the sensor effect. Our simulationsshow a quite good correlation with a linear law ifthe permittivity is constant (the variation of the ra-tio w/ϕ is about 30%); on the contrary, the resultsobtained with the other two models of permittivity donot agree with these experimental results. More pre- alidation of a model for an ionic electro-active polymer in the static case Imposed electric potential (V) D e fl ec t i o n ( mm ) ( a ) Imposed electric potential (V) B l o c k i n g f o r ce ( m N ) ( b ) Imposed electric potential (V) A n g l e o f r o t a t i o n ( d e g ) ( c ) Figure 12.
Influence of the electric potential: (a) on the deflection; (b) on the blocking force; (c) on the angle of rotation (samecolours marks and curves as in figure 9). cisely, in the constant case, the curve that best fits ournumerical results is a power law of exponent 1 .
16 (fig-ure 12); this curve can hardly be distinguished from anexperimental linear curve (Tixier & Pouget 2018). Theblocking force and the angle of rotation follow approx-imately the same trend (good correlation with a linearlaw and with a power law of exponent 1 . .
5. In the linear case, the momentis independent of the imposed potential for ϕ (cid:38) e andsurface Ll subject to a potential difference ϕ . Wesaw that the cations accumulate near the negativeelectrode over a thickness of e (cid:48) << e , and that therewas a similar thickness area without cations near thepositive electrode. The electric charge density is about ρZ + = F C moy ( B −
1) near the positive electrode and ρZ − = − F C moy near the negative one. The capacitivecharge Q can therefore be estimated by Q = ( ρZ + − ρZ − ) φ lLe (cid:48) = F C moy φ lLe (cid:48) B . (55)The permittivity is then worth ε = QeLlϕ = F C moy φ ee (cid:48) ϕ B . (56)In the case of constant permittivity, B = A hence ε ∼ − − − . In the case of a linear or affinepermittivity, B (cid:39) ε ∼ − . These valuesare consistent with the permittivity values derivedfrom capacity measurements published in the literature(Nemat-Nasser 2002, Farinholt & Leo 2004). The hydrated IPMC can be considered as a five-layer capacitor: the electrodes of thicknesses e and e ,the central area of thickness e very large compared tothe previous thicknesses and of zero electric charge,and two very thin zones of the polymer of respectivethicknesses e and e in the vicinity of the positive andnegative electrodes. Each element has a permittivity ε i and a capacity C i = ε i Le i , and the five elementsare in series. The permittivity of a material beingclosely related to its conductivity, so here to the cationconcentration, ε = ε >> ε >> ε (cid:38) ε . In addition e is much greater than e , e , e and e . We deduce ε = e e + e ε + e ε + e ε + e ε (cid:39) ε , where e denotes the total thickness. The overallpermittivity of the strip subject to an electric field istherefore very close to the permittivity of the centralpart, which is that of the hydrated polymer withoutelectric field. Our average permittivity values rangefrom 5 10 − F m − to 10 − F m − depending onthe permittivity model and are thus compatible withthe permittivity values measured by Deng & Mauritz(1992) and Wang et al. (2014), which range from10 − F m − to 5 10 − F m − .
5. Conclusion
We studied the flexion of an ionic polymer stripusing a model based on the thermomechanical ofcontinuous media that we had previously developed.The strip, clamped at one of its ends, is subject toa continuous potential difference applied between itstwo faces (static case). The other end is either free orblocked by a shear force. The mechanical quantities(deflection, blocking force and angle of rotation) weredetermined using a beam model in large displacements.The material chosen to perform the simulations isNafion ® Li + . Three models of permittivity (constant,linear and affine functions of cation concentration)are examined. The permittivity values we used arein good agreement with dielectric spectroscopy and alidation of a model for an ionic electro-active polymer in the static case ® and to studyother configurations such as a strip clamped at its twoends. Notations k = 1 , , , , i subscripts respectively representcations, solvent, solid, solution (water and cations)and interface; quantities without subscript refer to thewhole material. Superscript denotes a local quantity;the lack of superscript indicates average quantity atthe macroscopic scale. Superscript T indicates thetranspose of a second-rank tensor. Overlined lettersdenote dimensionless quantities. A i , B i : dimensionless constants; C , C moy : cations molar concentrations (relative tothe liquid phase); D : mass diffusion coefficient of the cations in theliquid phase; (cid:126)D : electric displacement field; e : half-thickness of the strip; E , G , λ : Young’s and shear modulus, first Lam´econstant; (cid:126)E : electric field; F = 96487 C mol − : Faraday’s constant ; (cid:126)F p ( (cid:126)F ps ): blocking force on large (small) displace-ments; (cid:126)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 moleof sulfonate groups); (cid:126)M p ( (cid:126)M pA ): bending moment; (cid:126)n ( n ): normal vector (coordinate) to the beam; p : pressure; −→ p p : distributed electric force (cid:126)Q : heat flux; R = 8 , J K − : 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; (cid:126)u : displacement vector; U : internal energy density; v k : partial molar volume of component k (relative tothe liquid phase); (cid:126)V ( (cid:126)V k ): velocity; w ( w s ): deflection of the beam on large (small)displacements; Z ( Z k ): total electric charge per unit of mass; ε ( ε , ε moy ): permittivity (average permittivity); (cid:15) (cid:101) : strain tensor; η : 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 thewhole material; ρ k : mass density relative to the volume of the phase; σ (cid:101) ( σ k (cid:102) ), σ e (cid:102) : stress tensor, equilibrium stress tensor; φ k : volume fraction of phase k ; ϕ ( ϕ ): electric potential (imposed electric potential); References
Aureli, M. and Prince, C. and Porfiri, M. and Peterson, S.D.2010
Smart Materials and Structures , 015003.Bahramzadeh, Y. and Shahinpoor, M. 2014 Soft Robotics (1),38–52.Bakhtiarpour, P. and Parvizi, A. and M¨uller, M. and Shahin-poor, M. and Marti, O. and Amirkhani, M. 2016 SmartMaterials and Structures , 015008 doi:10.1088/0964-1726/25/1/015008.Bar-Cohen, Y. 2005 WIT Transactions on State of the Art inScience and Engineering , 66–81.Barclay Satterfield, M. and Benziger, J. B. 2009 J. Polym. Sci.Pol. Phys. (1), 11–24. alidation of a model for an ionic electro-active polymer in the static case Bauer, F. and Denneler, S. and Willert-Porada, M. 2005
Journalof Polymer Science part B : Polymer Physics (7), 786–795.Biot, M.A. 1977 International Journal of Solids and Structures , 579–597.Bluhm, J. and Serdas, S. and Schr¨oder, J. 2016 Archive ofApplied Mechanics , 3–19.Brunetto, P. and Fortuna, L. and Graziani, L. and Strazzeri, S.2008 Smart Materials and Structures , 025029.Cappadonia, M. and Erning, J. and Stimming, U. 1994 Journalof Electroanalytical Chemistry (1), 189–193.Cha, Y. and Shen, L. and Porfiri M. 2013
Smart Materials andStructures , 055027.Chab´e, J. 2008 Etude des interactions mol´eculaires polym`ere-eaulors de l’hydratation de la membrane Nafion, ´electrolytede r´ef´erence de la pile `a combustible. PhD thesisUniversit´e Joseph Fourier Grenoble I http://tel.archives-ouvertes.fr/docs/00/28/59/99/PDF/THESE JCS.pdf .Chen, Z. 2017 Robotics and Biomimetics (24),https://doi.org/10.1186/s40638-017-0081-3Chikhaoui, M.T. and Benouhiba, A. and Rougeot, P. andRabenorosoa, K. and Ouisse, M. and Andreff, N. 2018 Annals of Biomedical Engineering (10), 1511–1521.https://doi.org/10.1007/s10439-018-2038-2.Collette, F. 2008 Vieillissement hygrothermique duNafion. PhD Thesis Ecole Nationale Sup´erieuredes Arts et M´etiers http : //tel.archives − ouvertes.fr/docs/ / / / /P DF/T hese F loraineCOLLET T E .pdf .Coussy, O. 1995 Mechanics of porous continua.
WileyChichester.de Groot, S. R. and Mazur, P. 1962
Non-equilibrium thermody-namics.
North-Holland publishing company AmsterdamDeng, Z.D. and Mauritz, K.A. 1992
Macromolecules (10),2739–2745.Deole, U. and Lumia, R. and Shahinpoor, M. and Bermudez, M.2008 Journal of Micro-Nano Mechatronics , 95–102.Fang, B.K. and Ju, M.S. and Lin, C.C.K. 2007 Sensors andActuators A (2), 321–329.Farinholt, K. and Leo, D.J. 2004
Mechanics of Materials (5),421–433.Farinholt, K. and Pedrazas, N.A. and Schluneker, D.M. andBurt, D.W. and Farrar, C.R. 2009 Journal of IntelligentMaterial Systems and Structures (5), 633–642.Festin, N. and Plesse, C. and Pirim, P. and Chevrot, C. andVidal, F. 2014 Sensors and Actuators B: Chemical ,82–88.Gierke, T.D. and Munn, G.E. and Wilson, F.C. 1981
Journal ofPolymer Science : Polymer Physics Edition (11), 1687–1704.Hasani, M. and Alaei, A. and Mousavi, M. S. S. and Esmaeili, E.and Kolahdouz, M. and Naeini, V. F. and Masnadi-Shirazi,M. 2019 Journal of Micromechanics and Microengineering , 085008 https://doi.org/10.1088/1361-6439/ab272c.He, Qingsong and Yu, Min and Song, Linlin and Ding, Haitaoand Zhang, Xiaoqing and Dai, Zhendong 2011 Journal ofBionic Engineering , 77–85.Ishii, M. and Hibiki, T. 2006 Thermo-fluid dynamics of two-phase flow.
Springer New-York.Jean-Mistral, C. and Basrour, S. and Chaillout, J.J. 2010
SmartMaterials and Structures , 085012.Jo, C. and Pugal, D. and Oh, I. and Kim, K.J. and Asaka, K.2013 Progress in Polymer Science , 1037–1066Lakshminarayanaiah, N. 1969 Transport Phenomena in Mem-branes.
Academic Press New-York.Lin, J. and Liu, Y. and Zhang, Q.M. 2012
Macromolecules (4),2050–2056.Matysek, M. and Lotz, P. and Schlaak H.F. 2009 Proc.SPIE 7287, Electroactive Polymer Actuators and Devices(EAPAD) , 72871D doi: 10.1117/12.819217. Mojarrad, M. and Shahinpoor, M. 1997
Proc. SPIE , 52–60.Nardinocchi, P. and Pezzulla, M. and Placidi, L. 2011
Journal ofIntelligent Material Systems and Structures (16), 1887–1897.Nemat-Nasser, S. 2002 Journal of Applied Physics (5,) 2899–2915.Nemat-Nasser, S. and Li, J. 2000 Journal of Applied Physics (7), 3321–3331.Newbury, K.M. 2002 Characterization, modeling and controlof ionic-polymer transducers. PhD thesis Faculty ofthe Virginia Polytechnic Institute and State UniversityBlacksburg, Virginia.Newbury, K.M. and Leo, D.J. 2002 Journal of IntelligentMaterial Systems and Structures (1), 51–60.Newbury, K.M. and Leo, D.J. 2003 Journal of IntelligentMaterial Systems and Structures (6), 343–357.Nguyen, N.T. and Dobashi, Y. and Soyer, C. and Plesse, C.and Nguyen, G.T.M. and Vidal, F. and Cattan, E. andGrondel, S. and Madden, J.D.W. 2018 Smart Materials andStructures , 115032.Park, I. and Kim, S.M. and Pugal, D. and Huang, L. and Tam-Chang, S.W. and Kim K.J. 2010 Applied Physics Letters (4), 043301.Pugal, D. and Jung, K. and Aabloob, A. and Kim, K.J. 2010 Polymer International , 279–289.Shahinpoor, M. 1999 Proc. of SPIE , 109–121.Shahinpoor, M. and Bar-Cohen, Y. and Simpson, J.O. andSmith, J. 1998
Smart Materials and Structures (6), R15–R30.Shahinpoor, M. and Kim, K.J. 2001 Smart Materials andStructures (4), 819–833.Silberstein, M. N. and Boyce, M. C. 2010 J. Power Sources (17), 5692–5706.Tiwari, R. and Kim, K. J. and Kim, S. M. 2008
Smart Structuresand Systems (5), 549. DOI: 10.12989/sss.2008.4.5.549.Tixier, M. and Pouget, J. 2014 Continuum Mechanics andThermodynamics (4), 465–481.Tixier, M. and Pouget, J. 2016 Continuum Mechanics andThermodynamics (4), 1071–1091.Tixier, M. and Pouget, J. 2018 Chapter 39: Modelling of anIonic Electroactive Polymer by the Thermodynamics ofLinear Irreversible Processes in Generalized Models andNon Classical Mechanical Approaches in Complex MaterialsSpringer Berlin.Vokoun, D. and Qingsong, He and Heller, L. and Min, Y. andZhen, D.D. 2015
Journal of Bionic Engineering (1), 142–151.Wallmersperger, T. and Horstmann, A. and Kr¨oplin, B. andLeo, D.J. 2009 Journal of Intelligent Material Systems andStructures , 741–750.Wang, Y. and Zhu, Z. and Chen, H. and Luo, B. andChang, L. and Wang, Y. and Li, D. 2014 SmartMaterials and Structures23