Constrained port Hamiltonian formulation of multiscale distributed parameter IPMC systems ⁎

8th IFAC Symposium on Mechatronic Systems 8th IFACAustria, Symposium Systems Vienna, Sept. on 4-6,Mechatronic 2019 8th IFAC IFAC Symposium Symposium on Mechatronic Systems 8th Systems Available online at www.sciencedirect.com Vienna, Austria, Sept. on 4-6,Mechatronic 2019 Vienna, Austria, Sept. on 4-6,Mechatronic 2019 8th IFACAustria, Symposium Systems Vienna, Sept. 4-6, 2019 Vienna, Austria, Sept. 4-6, 2019

ScienceDirect

IFAC PapersOnLine 52-15 (2019) 495–500

Constrained port Hamiltonian formulation Constrained port Hamiltonian formulation Constrained port Hamiltonian formulation of multiscale distributed parameter IPMC Constrained port Hamiltonian formulation of multiscale distributed parameter of multiscale distributed parameter IPMC IPMC systems parameter IPMC of multiscale distributed systems systems  systems ∗ Ning Liu, ∗ Yongxin Wu, ∗∗ Yann Le Gorrec ∗∗

Ning Liu, ∗∗ Yongxin Wu, ∗∗ Yann Le Gorrec ∗∗ Ning Liu, Yongxin Wu, Yann Le Gorrec ∗ Ning Liu, ∗ Yongxin Wu, ∗ Yann Le Gorrec ∗ e Bourgogne ∗ FEMTO-ST CNRS UMR 6174, Universit´ CNRS UMR 6174, e Bourgogne ∗ ∗ FEMTO-ST FEMTO-ST CNRS de UMR 6174, Universit´ Universit´ Bourgogne Franche-Comt´ e, 26 chemin l’´epitaphe, F-25030eBesan¸ con, France. ∗ Franche-Comt´ e , 26 chemin de l’´ e pitaphe, F-25030 ccon, FEMTO-ST CNRS de UMR 6174, Universit´ Bourgogne Franche-Comt´ e, [email protected]; 26 chemin l’´epitaphe, F-25030eBesan¸ Besan¸ on, France. France. (e-mail: [email protected]; (e-mail: [email protected]; [email protected]; Franche-Comt´ e, [email protected]; 26 chemin de l’´epitaphe, F-25030 Besan¸con, France. [email protected]). (e-mail: [email protected]; [email protected]). (e-mail: [email protected]; [email protected]; [email protected]). [email protected]). Abstract: In this paper, a constrained distributed parameter port-Hamiltonian model of Abstract: In this paper, a constrained distributed parameter port-Hamiltonian model of the ionic polymer composite actuator is proposed. This model describes the multiscale Abstract: In thismetal paper, a constrained distributed parameter port-Hamiltonian model of multiscale the ionic metal composite actuator is This model describes the Abstract: paper, a constrained distributed parameter port-Hamiltonian model of the ionic polymer polymer metal composite actuator is proposed. proposed. Thismulti-scale model describes theIn multiscale structure of In thethis system. Submodels are coupled by boundary elements. order to structure of system. Submodels are coupled by boundary elements. In order to the ionic the polymer metal composite actuator is proposed. Thismulti-scale model describes multiscale structure of the the system. areLagrangian coupled bymultipliers boundary multi-scale elements. In orderthe to preserve causality ofSubmodels the system, are introduced to the deal with preserve the causality ofSubmodels the system, Lagrangian multipliers are introduced to deal with structure of the system. are coupled boundary elements. Instructure orderthe to coupling between the electro-stress diffusion in by the polymer multi-scale and the flexible beam preserve the causality of the system, Lagrangian multipliers are introduced deal with the coupling between the electro-stress diffusion in the polymer and the flexibleto beam structure preserve the causality of the system, Lagrangian multipliers are introduced to deal with the of the actuator. Finally, a structure-preserving discretization scheme and some appropriate coupling between the electro-stress diffusion in the polymer and the flexible beam structure of Finally, aa structure-preserving discretization scheme and some appropriate coupling between thetoelectro-stress diffusion in the polymer and the flexible beam structure of the the actuator. actuator. Finally, structure-preserving discretization scheme and appropriate projections are used derive an explicit model suitable for simulation. Thesome accuracy of the projections are used to derive an explicit model suitable for simulation. The accuracy of of the isactuator. Finally, a structure-preserving scheme and appropriate projections are used toexperimental derive an explicit suitable for simulation. Thesome accuracy of the the model verified using data. model discretization model is verified usingtoexperimental data. model suitable for simulation. The accuracy of the projections are used derive an explicit model is verified using experimental data. © 2019,isIFAC (International Federation ofdata. Automatic Control) Hosting by Elsevier Ltd. All rights reserved. model verified using experimental Keywords: Constrained port Hamiltonian system, infinite dimensional system, multi-scale Keywords: Constrained port Hamiltonian system, Keywords: model Constrained portIPMC Hamiltonian system, infinite infinite dimensional dimensional system, system, multi-scale multi-scale modeling, reduction, actuator modeling, model reduction, IPMC actuator Keywords: Constrained port Hamiltonian system, infinite dimensional system, multi-scale modeling, model reduction, IPMC actuator modeling, model reduction, IPMC actuator 1. INTRODUCTION tems interconnected each other using boundary multi1. tems interconnected each other using boundary 1. INTRODUCTION INTRODUCTION tems (BMS) interconnected other By using boundarythemultimultiscale couplingeach elements. considering outscale (BMS) couplingeach elements. By considering themultiout1. INTRODUCTION tems interconnected other using boundary domain variables as uniform (Nishida et al., 2011), the (BMS) coupling elements. By considering the outIonic polymer metal composites (IPMCs) are electro- scale domain variables as uniform (Nishida et al., 2011), the Ionic polymer metal composites (IPMCs) are electroscale (BMS) coupling elements. By considering the outBMS works as a differential gyrator, which lets the outIonic metalcancomposites (IPMCs) electrovariables as uniform (Nishida et al., 2011), the active polymer systems that be used either as an are actuator or domain BMS works as a differential gyrator, which lets the outactive systems that can be used either as an actuator or variables as uniform (Nishida et al.,lets 2011), the BMS works as a differential gyrator, which thefuncoutpolymer metal composites (IPMCs) electrodomain variables be multiplied by a characteristic active systems that be usedofeither as an are actuator or domain aIonic sensor. Among thecan diversity electro-active materials domain variables be multiplied by a characteristic funca sensor. Among the diversity of electro-active materials BMS works as a differential gyrator, which lets the outactive systems that can be used either as an actuator or tion, meanwhile, makes the in-domain variables be inteasuch sensor. Among the diversity of electro-active materials domain variables be multiplied by a characteristic funcas piezoelectric materials, magnetostrictive meanwhile, makes the in-domain variables be funcintesuch as magnetostrictive materials tion, variablesHowever, be multiplied by a characteristic a sensor. Among thematerials, diversity ofused electro-active gratedmeanwhile, spatially. there exists a variables conflict ofbe causalsuch as piezoelectric piezoelectric materials, tion, makes the in-domain inteetc., IPMCs are more and moremagnetostrictive in differentmaterials applica- domain grated spatially. However, there exists a conflict of causaletc., IPMCs are more and used different applica- tion, meanwhile, makes the in-domain integrated spatially. However, exists a variables conflict ofbe causalsuch fields, as piezoelectric magnetostrictive ity due to the coupling ofthere the mechanical properties of etc., IPMCs arebiomedical morematerials, and more more used in in bio-manipulation differentmaterials tion e.g applications, applicaity due to the coupling of the mechanical properties of tion fields, e.g biomedical applications, bio-manipulation grated spatially. However, there exists a conflict of causality due to the coupling of the mechanical properties of etc., IPMCs are more and more used in different applicathe gel and the mechanical structure (passive moment tion fields, e.g biomedical applications, systems bio-manipulation and microor macro-electromechanical (Shahin- the gel and the mechanical structure (passive moment and microor macro-electromechanical systems (Shahinity due to the coupling of the mechanical properties of tion fields, e.g biomedical applications, bio-manipulation coupling of equation (54) in (Nishida et al., 2011)). To and microor macro-electromechanical systems (Shahinthe gel and the mechanical structure (passive moment poor, 2016) due to their low-cost voltage, large deforma- coupling of equation (54) in (Nishida et al., 2011)). To poor, 2016) due to their voltage, large deformaand theconflict, mechanical structure (passive moment and microor (Shahin- the deal gel with this weinconsider a et multiscale model poor, 2016)working duemacro-electromechanical to frequency their low-cost low-cost voltage, largecapability deformacoupling of equation (54) (Nishida al., 2011)). To tion, wide ranges andsystems their deal with this conflict, we aa et multiscale model tion, wide ranges and their capability equation (54) (Nishida al., To deal withofLagrange this conflict, weinconsider consider multiscale model poor, 2016)working dueaqueous to frequency theirenvironments. low-cost voltage, largeconsist deformaincluding multiplier to account for 2011)). these metion, wide working frequency ranges and their capability of working in IPMCs of coupling including Lagrange multiplier to account for these meof working in aqueous environments. IPMCs consist of deal with this conflict, we consider a multiscale model including Lagrange multiplier to account for these metion, wide working frequency ranges and their capability chanical constraints, and numerically simulate the model of working in aqueous IPMCs consistgel. of chanical constraints, and numerically simulate the model a double electrode layerenvironments. filled with a polyelectrolyte a electrode layer filled a gel. Lagrange to coupling account for the these meof double working insolvent aqueous environments. IPMCs of including more precisely, whichmultiplier includes all relations. The aCations double electrode layer filled with with a polyelectrolyte polyelectrolyte gel. chanical constraints, numerically simulate model and molecules migrate toward consist the cathmore precisely, whichand includes all coupling relations. The Cations and solvent molecules migrate toward cathconstraints, and numerically simulate the model a double electrode with apotential polyelectrolyte gel. chanical resulting systemwhich of differential algebraic equation (DAE) is Cations solventlayer molecules migrate towardisthe the cathmore precisely, includes all coupling relations. The ode whenand a difference in filled the electric imposed resulting system of algebraic equation (DAE) is ode when a difference in electric potential imposed precisely, which includes all coupling relations. The resulting system of differential differential algebraic equation (DAE) is Cations and solvent molecules migrate towardis the cathreduced to an ordinary differential equation (ODE) using ode when atwo difference in ofthe the potential is imposed across the terminals itselectric double electrode layer. As a more reduced to an ordinary differential equation (ODE) using across the two terminals of its double electrode layer. As a resulting system of differential algebraic equation (DAE) is reduced to an ordinary differential equation (ODE) using ode when a difference in the electric potential is imposed coordinates projection. across the twothe terminals its double electrode layer. As a coordinates projection. consequence, cathodeofside swells while the anode side consequence, cathode swells while the side to an ordinary differential using across the twothe terminals ofside its double electrode layer. As a reduced The present paper is organized asequation follows. (ODE) In Section 2 consequence, the cathode side swells while the anode anode side coordinates projection. shrinks, entailing a bending effect to the anode side (Park The presentprojection. paper is organized as follows. In Section 2 shrinks, entailing a effect to the anode side (Park consequence, the cathode swells while theand anode side coordinates is given the constrained port Hamiltonian model of the shrinks, entailing a bending bending effect tostructure the anode side (Park The present paper is organized as follows. In Section 2 et al., 2010). Based on its side physical working given the constrained port model of et al., Based on physical and working The present paperInisSection organized follows. In Section 2 is given the constrained port3,Hamiltonian Hamiltonian model of the the shrinks, entailing amodels bending effect tostructure the anode side (Park is multiscale IPMC. a as finite difference method et al., 2010). 2010). Based on its its for physical structure and working principle, various IPMCs have been proposed multiscale IPMC. In Section 3, aa finite difference method principle, various models for IPMCs have been proposed is given the constrained port Hamiltonian model of the multiscale IPMC. In Section 3, finite difference method et al., 2010). Based on its physical structure and working on staggered grids is applied to discretize the system and principle, various models for IPMCs have been proposed in the literature, going from the black box model (Xiao on staggered grids is applied to discretize the system and in the literature, from the box model (Xiao multiscale IPMC. In 3, finite difference method principle, various going models for IPMCs have proposed the final model is reduced by using coordinates projection. in the literature, going the black black boxbeen model (Xiao on staggered grids is Section applied to adiscretize the system and and Bhattacharya, 2001)from to models using more physical final model is reduced by using coordinates projection. and Bhattacharya, 2001) to models using more physical on staggered grids is applied to discretize the system and in the literature, going from the black model (Xiao the Numerical simulation and by conclusions are given in Section and Bhattacharya, 2001) to models more physical the final model is reduced using coordinates projection. insight (Shahinpoor, 2016; Branco et using al.,box 2012). Numerical simulation and conclusions are given in Section insight (Shahinpoor, 2016; Branco et al., 2012). the final model is reduced by using coordinates projection. Numerical simulation and conclusions are given in Section and Bhattacharya, 2001) to models using more physical 5, respectively. insight (Shahinpoor, 2016; Branco and et al., 2012).of complex 44 and A powerful tool for the modeling control A powerful tool for modeling and control simulation and conclusions are given in Section 4 and and 5, 5, respectively. respectively. insight (Shahinpoor, 2016; Brancocalled et al., 2012).of A powerful toolnonlinear for the the modeling and control of complex complex Numerical multi-physical systems, port-Hamiltonian multi-physical nonlinear systems, called port-Hamiltonian 4 and 5, respectively. 2. MODELING OF IPMC A powerfulhas toolnonlinear for introduced the modeling and control of multi-physical systems, called port-Hamiltonian approach, been and developed in complex the last 2. approach, has been introduced and developed in 2. MODELING MODELING OF OF IPMC IPMC multi-physical called port-Hamiltonian approach, has nonlinear been and developed in the the last last decade (Maschke andintroduced vansystems, der Schaft, 1992). decade (Maschke and van der Schaft, 1992). 2. MODELING OF IPMC is of length L, The IPMC under investigation (cf. Fig. approach, has been introduced and developed in the last decade (Maschke and van der modeling Schaft, 1992). The first port-Hamiltonian of IPMC actua- The IPMC under investigation (cf. Fig. 1) 1) is of length L, The first port-Hamiltonian modeling of IPMC actuawidth b and thickness h. It consists three The IPMC under investigation (cf. Fig. 1) is sub-systems of length L, decade (Maschke and van der Schaft, 1992). The first port-Hamiltonian modeling of IPMC actuators has been proposed in (Nishida et al., 2011). This width b and thickness h. It consists of of three tors been proposed in al., 2011). This IPMC investigation (cf. 1. Fig. 1) is sub-systems of length L, width b andunder thickness h. It in consists of three sub-systems at different scales as shown Fig. The has first port-Hamiltonian modelinget tors has been proposed in (Nishida (Nishida et of al.,IPMC 2011). This The model consists in three sub-components which are actuamultiat different scales as shown in Fig. 1. model consists in three sub-components which are multiwidth b and thickness h. It consists of three sub-systems at different scales as shown in Fig. 1. First, an electrical model, which is at a scale of tors has been proposed in by (Nishida et which al., 2011). This model consists in three sub-components are multiscale, and are all described distributed parameter sys- First, an electrical model, which is at a scale of nanometer, nanometer, scale, and are all described by distributed parameter sysdifferent scales as shown in Fig. 1.a scale of is used represent the fractal-like structure ofnanometer, the double First, anto electrical model, which is at model consists in three sub-components which are multiscale, and are all described by distributed parameter sys- at is used to represent the fractal-like structure of the double First, an electrical model, which is at a scale of nanometer,  electrical layers. The dynamics of structure the polyelectolyte gel, is used to represent the fractal-like of the double scale, and are all described by distributed parameter sysThis work is supported by the INFIDHEM project and the  This work is supported by the INFIDHEM project and the electrical layers. The dynamics of the polyelectolyte gel, is to layers. represent thedynamics fractal-like the double  electrical of structure thebypolyelectolyte gel, at used a scale of 100The µm, is described an of electro-stress  Bourgogne-Franche-comt´ e Region ANER project under the reference This work work is is supported supported by the the INFIDHEM project and the the This by INFIDHEM project and at a scale of 100 µm, is described by an electro-stress Bourgogne-Franche-comt´ e Region ANER project under the reference  electrical layers. The dynamics of the polyelectolyte gel, at a scale of 100 µm, is described by an electro-stress diffusion coupling model. At last, the global mechanical code ANR-16-CE92-0028 and 2018Y-06145, respectively. This work is supported by the INFIDHEM project and the Bourgogne-Franche-comt´ e Region ANER project under the reference Bourgogne-Franche-comt´ Region ANER project under the reference diffusion last, thebyglobal mechanical code ANR-16-CE92-0028eand 2018Y-06145, respectively. at a scalecoupling of 100 model. µm, is At described an electro-stress Bourgogne-Franche-comt´ Region ANER project under the reference code ANR-16-CE92-0028 ANR-16-CE92-0028eand and 2018Y-06145, respectively. diffusion coupling model. At last, the global mechanical code 2018Y-06145, respectively. coupling Atreserved. last, the global mechanical code ANR-16-CE92-0028 and 2018Y-06145, respectively. 2405-8963 © 2019, IFAC (International Federation of Automatic Control) diffusion Hosting by Elsevier Ltd.model. All rights

Copyright © 2019 IFAC 1242 Copyright © under 2019 IFAC 1242Control. Peer review responsibility of International Federation of Automatic Copyright © 2019 2019 IFAC IFAC 1242 Copyright © 1242 10.1016/j.ifacol.2019.11.724 Copyright © 2019 IFAC 1242

2019 IFAC MECHATRONICS 496 Vienna, Austria, Sept. 4-6, 2019

Ning Liu et al. / IFAC PapersOnLine 52-15 (2019) 495–500

and combining equation (2), equation (1) can be written 1 on  the  form    f1 0 ∂ξ e1 = , fr1 ∂ξ 0 er1       T T (V + Vc e1 (Lξ )) f∂ξ (e1 (0) e1 (Lξ )) = . = T T e∂ξ (−er1 (0) er1 (Lξ )) (−I er1 (Lξ )) (3) Assuming that the impedances are infinite, the current at the endpoint of each fractal structure is zero, namely i(Lξ , t) = er1 (Lξ ) = 0 (Nishida et al., 2011).

electrodes polymer

diffusion

fixed anion

water molecule

hydrated cation

mobile cation

Fig. 1. IPMC structure and shape (Shahinpoor, 2016). deformation of IPMC is described by the Timoshenko beam model, whose scale is centimeter. Both the electrical system and the electro-stress diffusion system are modeled locally, whereas the mechanical beam system is modeled globally. In this section, we present each sub-system and their coupling through boundary or in domain multiscale elements. 2.1 Electrical system The two electrodes of Fig. 1 are modeled by the distributed RC circuit. The voltage V is supposed to be uniformly distributed on the double layers electrodes.

2.2 Electro-stress diffusion system An electro-stress diffusion coupling model is considered to describe the swelling and shrinking dynamics in the gel (Nishida et al., 2011). Compared to the diffusion in the liquid phase of the gel, the deformation of the solid phase is so fast that it is considered as quasi-static. Consequently, the mechanical dynamics of the gel is not represented explicitly, and the radius of curvature of the gel is derived from the rotational angle of the beam, leading to an algebraic constraint. Deformation of the solid phase This deformation is assumed to be symmetric. The schematic is shown in Fig. 3. The radius of curvature R(x, t) fluctuates along the x-axis,

Fig. 2. Infinite dimension electrical system.

-

With the idea of (Nishida et al., 2008), (Nishida et al., 2011), each fractal-like structure (see the red circle in Fig. 1) is referenced as ξ in a virtual coordinate. Lξ denotes the length of each fractal-like structure, R1 is the resistance between two adjacent branches of fractal-like structure, and R2 and C2 are the resistive and capacitive impedance of each branch, respectively (Bao et al., 2002). Each fractal-like structure is connected to the electroactive gel through its boundary at ξ = 0. The continuity equation and the Kirchhoff’s current law (KCL) yield: ∂i(ξ, t) ∂Q(ξ, t) =− , (1) ∂t ∂ξ where Q(ξ, t) is the charge density of each capacitor. By applying the Kirchoff’s voltage law (KVL), one gets : ∂v(ξ, t) + R1 (ξ)i(ξ, t) = 0. (2) ∂ξ Let v(0, t) = V + Vc , and i(0, t) = I, where Vc corresponds to the voltage coming from the gel. By defining (Nishida et al., 2008): Q(ξ, t) ∂Q(ξ, t) e1 (ξ, t) = v(ξ, t) = + R2 (ξ) , C2 (ξ, t) ∂t fr1 (ξ, t) = ∂/∂ξ (vc2 (ξ, t) + R2 (ξ)∂Q(ξ, t)/∂t) , ∂Q(ξ, t) f1 (ξ, t) = − , ∂t er1 (ξ, t) = −fr1 (ξ, t)/R1 (ξ),

0

Fig. 3. Deformation of gel in one dimension. but is assumed to be locally homogeneous (Nishida et al., 2008), i.e. ∂R(x, t)/∂x = 0 always holds in each z domain. The displacement of each volume point projected in Cartesian coordinate is given by: z z uz = uz (z, t), ux = x, uy = y. (4) R(x, t) R(x, t) The swelling ratio fs1 (z, x, t) of the solid part is defined as the divergence of the displacement tensor:   T ∂ ∂ ∂ fs1 (z, x, t) = ∇ · u = ∂x ∂y ∂z [ux uy uz ] (5) ∂uz 2z + . = R(x, t) ∂z According to the hypothesis of symmetric deformation, the linear formulation about the stress tensor and the displacement is expressed as:  ∂uk ∂ui ∂uj 2  ∂uk σij = K δij + G( + − δij ), (6) ∂xk ∂xj ∂xi 3 ∂xk k

k

where K and G are the bulk modulus and the shear modulus of the gel, respectively, and δij is the Dirichlet function. The stresses are average variables in the IPMC. Further informations on different stresses of IPMC model

1

For the purpose of simplicity, ∂/∂ξ is denoted as ∂ξ and the symbol t is omitted in the following context.

1243

2019 IFAC MECHATRONICS Vienna, Austria, Sept. 4-6, 2019

Ning Liu et al. / IFAC PapersOnLine 52-15 (2019) 495–500

are available in (Zhu et al., 2012). As a result, equation (4) and (6) yield the following expression:   2 2G σxx (z, x, t) = K − G fs1 (z, x, t) + z, 3 R(x, t)   (7) 4 4G σzz (z, x, t) = K + G fs1 (z, x, t) − z. 3 R(x, t) Dynamics of the liquid phase Two coupled phenomena can be distinguished in the liquid phase: electro-osmosis, and water transport. These coupled phenomena were formulated in the work of (De Gennes et al., 2000), covering the transport of ions and the solvent: je = −σe ∇ψ − λ∇p, js = −k∇p − λ∇ψ,

(8)

where je and js represent the electric current density and the water flux density, respectively. σe is the conductance, λ stands for the Onsager’s coupling constant, k denotes the Darcy’s permeability, ψ is the electric field, and p represents the water pressure in the network (De Gennes et al., 2000). It is supposed that the liquid goes only in the z direction, so ∇p = ∂p/∂z is the mechanical force. The bulk region of the gel satisfies the charge neutrality condition, namely ∇je = 0. The incompressibility of poly-electrolyte gel is also assumed in this work (Yamaue et al., 2005), i.e.:   4 4G z. (9) p = σzz = K + G fs1 (z, x, t) − 3 R(x, t) Thus, the gradient of pressure can be calculated as:   ∂p 4 ∂fs1 (z, x, t) 4G ∇p = = K+ G − . (10) ∂z 3 ∂z R(x, t) So equation (8) can be rewritten as:  2  λ λ js (z, x, t) = je + − k ∇p σe σe ∂f (z, x, t) λ s1 + 1Z je (t) + 1Z Φ(x, t), = −D ∂z σe (11) with      4 λ2 λ2 4G D = k − K + G , Φ(x, t) = k − , σe 3 σe R(x, t) (12) where 1Z stands for the characteristics function of domain z. It distributes the boundary values λ/σe je (t) and Φ(x, t) as uniform constants into z domain. In the liquid phase, a swelling ratio fs2 is also introduced. It follows the conservation law that: ∂js (z, x, t) ∂fs2 (z, x, t) =− . (13) ∂t ∂z This equation can then be reformulated in the PHS framework as:    ∂fs2     fs2 − ∂t f2 0 ∂z =  ∂fs2  = . (14) s1 fr2 ∂z 0 −D ∂f ∂z ∂z The effort and the boundary variables are:

497



 h  s1 −D ∂f −2 ∂z h         ∂f  s1   fs2 f∂z e2 2   −D ∂z   , = =  . s1 er2 e∂z −D ∂f  −fs2 − h2  ∂z     fs2 h2

(15) Since the solid and liquid phases are strongly mixed with each other, we have fs1 = fs2 = fs . As hinted by equation (15), 1Z λ/σe je (t) and 1Z Φ(x, t) do not appear explicitly in the dynamics, while they play a role of input in order to match the impermeable assumption js (±h/2, t) = 0. Bending moments generated in the gel According to (Nishida et al., 2008), the stress σxx can be divided into two parts: the active one σa = (K −2/3G)fs (z, x, t) related to the active swelling of the gel, and its passive counterpart σp = 2G/R(x, t)z corresponding to the storing energy. The active stress can generate an active moment Ma :  h2  h2 σa (z, x, t)bzdz = Ba fs (z, x, t)dz, Ma (x, t) = −h 2

−h 2

(16) with Ba (z) = (K − 2/3G)bz. Besides, the passive moment Mp comes from the passive stress σp :  h2 Gbh3 . (17) σp (z, x, t)bzdz = Mp (x, t) = 6R(x, t) −h 2

Regarding to the mechanical model along x-axis, the curvature 1/R is related to the angular strain ∂θ/∂x of the IPMC via the geometric relationship: ∂θ 1 + = 0. (18) R(x, t) ∂x At the initial phase of actuation, the active moment Ma is much larger than the passive moment Mp , as evident from the phenomenon of quick bending of IPMC. As the curvature increases gradually, Mp is getting larger than Ma , which makes the IPMC to bend back slowly.

Coupling with the electrical system In light of equation (11), the interconnection between the electro-stress diffusion system and the electrical system is through boundary variables as 1Z λ/σe je (t), ∂fs /∂z and I, Vc . je can be related to I by: 1 I(t). (19) je (t) = Lb Given that the two pairs of energy variables 1Z λ/σe je (t), ∂fs /∂z and I, Vc are of different scales and are defined in domains z and ξ respectively, a coupling element BMS is proposed to realize the interconnection (Nishida et al., 2011), as represented in Fig. 4. By crossing the BMS, λ/σe je (t) is multiplied by the characteristic function 1Z , which signifies an integration in domain z. fs |∂z denotes the space integration of ∂fs /∂z with fs |∂z = fs (h/2, t) − fs (−h/2, t). Based on the power conservation law, fs |∂z is transformed into voltage Vc (t) via the gyrator GY :   λ λ h h Vc (t) = − fs (t)|∂z = − fs ( , t) − fs (− , t) . σe b σe b 2 2 (20)

1244

2019 IFAC MECHATRONICS 498 Vienna, Austria, Sept. 4-6, 2019

-

Ning Liu et al. / IFAC PapersOnLine 52-15 (2019) 495–500

-

Fig. 4. Bond graph of the coupling between the domain z and ξ. Coupling with the mechanical system At the macroscale, the electro-stress diffusion model connects with the mechanical model through two bending moments, Ma and Mp , and the angular velocity ∂θ(x, t)/∂t. 2.2.5.1. Coupling through the active moment In view of the power conservation, an additional term is added into equation (13) to match the output of active moment: ∂js (z, t) ∂fs (z, x, t) ∂θ(x, t) =− + 1Z Ba . (21) ∂t ∂z ∂t The latter term can be regarded as a diffusion term of the mass conservation, since the gel consists of multiple molecules. The bond graph of this interconnection is similar to the one in Fig. 4. 2.2.5.2. Coupling through the passive moment As for the coupling via the passive moment Mp , it is supposed to make the connection with Φ(t), since both the gel model and the beam model have the same curvature, 1/R(t) and ∂θ/∂x. Because Φ is a flow source for this electro-stress diffusion system and Mp is the output of this system, a Lagrangian multiplier λL is proposed here to deal with the causality, as shown in Fig. 5. where A = Bp with

Mp (x, t) = Φ(x, t)Bp ,      (22) ∂θ(x, t) h h −1 = Bp fs , x, t − fs − , x, t . ∂t 2 2 With the Lagrangian λL = Φ, equation (22) is rewritten as:       1 Φ λL , λL = = Bp Mp B p λL   (23) fs |∂z ∂θ (1 Bp ) −∂θ = fs |∂z − Bp = 0, ∂t ∂t which reveals that the arrow of the Lagrangian multiplier in the bond graph Fig. 5 is an effort source with zero flow. This ensures the power conservation. Accordingly, equation (15) changes to:        0 ∂z e2 f2 −Ba ∂θ(x,t) ∂t = + , fr2 ∂z 0 er2 0 (24)     h  h T  s s f∂z −D ∂f − 2 −D ∂f ∂z ∂z 2 . =     T e∂z −fs − h2 fs h2

2.3 Mechanical system

The mechanical deformation of IPMC can be represented by a classic Timoshenko beam with x ∈ [0, L]. The dynamics equation is reformulated under port Hamiltonian framework as (Villegas, 2007): x (x, t)  0 ∂ 0 −1 e (x, t) 0 x

3

∂ ∂t

3

x4 (x, t) = ∂x 0 0 0  e4 (x, t) + 0 Mext , x5 (x, t) x6 (x, t)

0 1

0 0 ∂x 0 ∂x 0

0 1

e5 (x, t) e6 (x, t)

(25)

where x3 (x, t) = ∂x ω(x, t) − θ(x, t), x4 (x, t) = ρA(x)∂t ω(x, t), x5 (x, t) = ∂x θ(x, t),x6 (x, t) = 1 x4 (x, t), ρI(x)∂t θ(x, t), e3 = GA(x)x3 (x, t), e4 = ρA(x) 1 e5 = EI(x)x5 (x, t), e6 = ρI(x) x6 (x, t). The distributed bending moment comes from the electrostress diffusion system, and reads: Mext = Ma + Mp . (26) According to (Le Gorrec et al., 2005), the boundary variables are calculated as:     T f∂x (e4 (0) e3 (L) e6 (0) e5 (L)) . (27) = T e∂x (−e3 (0) e4 (L) −e5 (0) e6 (L)) 2.4 Global system The above three subsystems can be connected to a global system, which is expressed as:      f1

Fig. 5. Bond graph of the coupling through Ma , Mp and ∂θ/∂t.  −1 λ2 bh3 Bp = k− . 24 σe The relations between Φ and Mp , and fs |∂z and ∂θ/∂t are expressed as follows: 1245

0 ∂ξ

fr1   ∂ξ     f2   0    fr2   0  =  f3   0     f4   0    f   0  5  f6

   f

0



0

0

0

0

0

0

0

0

0

0

0

0

∂z

0

0

0

0

0

∂z

0

0

0

0

0 ∂x 0

0

0

0

0

0

0 ∂x 0

0

0

0

0

0

0 int 0

0

0

1

0 ∂x

 J

0

e1

 er1  0   −1Z Ba   e2     er2  0    +A λ ,   1 L −1    e3     0   e4     ∂x   e5  0

e6

    e

(28)

2019 IFAC MECHATRONICS Vienna, Austria, Sept. 4-6, 2019

Ning Liu et al. / IFAC PapersOnLine 52-15 (2019) 495–500



T

Ba (·)dz, A1 = (0 0 −(∂z 1Z ) 0 0 0 0 Bp ) ,     h h ∗ − e2 − − Bp e6 = 0. A1 e = e2 (29) 2 2 Note that the operator ∂z 1Z equals to zeros when multiplied by a variable outside the z domain (e.g. λL ), while equals to the difference at the boundary when multiplied by a variable belongs to its domain (e.g. ∂z 1Z e2 = −e2 (h/2) + e2 (−h/2)). Suppose the parameters C2 , K, A, ρ, E, and I are homogeneous, the Hamiltonian of this system is:     1 Q2 1 dξdx + f 2 dzdx H(t) = 2 x ξ C2 2 x z s    x4 (x)2 x6 (x)2 1 2 2 GAx3 (x) + + EIx5 (x) + dx. + 2 x ρA ρI (30) where int = with

Z

3. DISCRETIZATION AND MODEL REDUCTION In order to preserve the geometric structure of the overall PHS, the finite difference method on staggered grids (Trenchant et al., 2018) is employed to discretize each scale of the distributed parameter model of the system. The principle of this method is to approximate effort and flow variables on different grids in order to preserve the power balances. A particular care has to be paid on the boundary conditions and to interconnexions between two different scales. In a second instance a projection is used in order to get rid of the Lagrangian multipliers. 3.1 Multiscale discretization Considering the fact that we have a multi-scale model, ξ and z are local coordinates while x is the global coordinate, which leads to the assumption that for each point in x, there is one corresponding ξ and z. Hence, there will be Ne (= Nξ × Nb ) elements for the electrical system, Ng (= Nz × Nb ) ones for the electro-stress diffusion system, and Nb elements for the mechanical system. For a sake of conciseness the discretization method is not detailed here. It yields the final dimensional model below: 









x˙d



 Jr



−M2













0

    e d

B

diag(Bp )

 gc



T T  b b , x6d = x16 · · · xN , e d = L d x d , Ld = x15 · · · xN 5 6 blockdiag (1/C2 , 1, GA, 1/(ρA), EI, 1/(ρI)) of proper dimension, Lr1 = diag(−1/R1 ), Lr2 = diag(−R2 ), Lr3 = diag(−D ), M 1 = σeλLb (g21 + g22 )g T Lr1 , M 2 = (I + D1 Lr1 D1T Lr2 )−1 D1 Lr1 , P1 = σeλLb M2 g2T , P2 = D2T Lr3 D2 − σeλLb M1 (D1T Lr2 M2 + I)g(g21 + g22 )T , P3 = M1 D1T (I − Lr2 M2 D1T ). Di with i ∈ {1, 26, 6, 63} and gk with k ∈ {21, 22} are matrices depending on the discretization steps and systems parameters. 3.2 Elimination of Lagrangian multipliers In this section the DAE (31) together with (32) will be reduced to an ODE, in order to perform the simulation and apply control strategies afterwards. The proposed method is based on the idea of coordinate projection as in (Wu et al., 2014). Given   S , M= (gcT gc )−1 gcT where S satisfies: S · gc = 0. T Now define X1 = (x1d x2d x3d x4d x5d ) , X2 = x6d ,  T ˜ 2 = M (X1 X2 )T , ˜1 X X T    B1 = −M2T M1 (D1T Lr2 M2 − I) T 0 0 0 and B2 = 0. Equation (31) is multiplied by the matrix M , and becomes:    ˙     ˜ ˜ X1 ˜˙ X 2

= M Jr M T M −T Ld M −1

    J˜r



˜d L

X1

 X˜ 2

+

M B1 0 V + λLd I 0

(33)

Meanwhile, equation (32) can be rewritten as:   ˜d gcT M T L

  

˜1 X ˜2 X

=0

(34)

g ˜c

˜ 1 does not depend It is implied from equation (33) that X on λLd . Hence the second line of equation (33) is substituted by equation (34), reforming a descriptor system as follows:          ˜˙ ˜ ˜ ˜ ˜ ˜ I 0 0 0

X1 ˜˙ 2 X

=

J11 J12 g˜c

Ld11 Ld12 ˜ d21 L ˜ d22 L

X1 ˜2 X

+

M B1 V 0

4. SIMULATION

M1 (D1T Lr2 M2 − I) g21 + g22     0  0     +  V +  0  λLd , 0     0 0 0



(35)

x1d ˙ M2 D1T P1 0 0 0 0 e1d P2 0 0 0 D26  e2d  x˙ 2d   P3 x˙ 3d   0 T  e  0 0 −D6T 0 −D63   3d    = x˙ 4d   0 0 D6 0 0 0  e4d   e  x˙   0 0 0 0 0 −D6T 5d 5d T x˙ 6d e6d 0 −D26 D63 0 D6 0

  

499

(31)

In this section, the proposed model will be verified by comparing the simulation with the experimental results. The simulation carries out with an IPMC actuator of dimension 45mm (in length)×5mm (in width)×0.2mm (in height), which contains the tetraethyl-ammonium ion TEA+ , and whose total resistance and capacitance equal to 23.6178 Ω and 0.0635 F, respectively. Mechanical parameters of this IPMC are illustrated in table 1. With a voltage 1V applied, the deformation of the endpoint of the IPMC strip and the output current are shown in Fig. 6. From Fig. 6, one can observe the simulation reproduces the same behaviors of the experimental results.

and the constraint equation 5. CONCLUSION (32) gcT ed = 0.    1,1  1,1 Nξ ,Nb T Nz ,Nb T , x2d = fs · · · fs , In this article, a detailed IPMC model is characterized where x1d = Q · · · Q  1   1  Nb T Nb T , x4d = x4 · · · x4 , x5d = under the port-Hamiltonian framework, and its dynamic x3d = x3 · · · x3 1246

2019 IFAC MECHATRONICS 500 Vienna, Austria, Sept. 4-6, 2019

Ning Liu et al. / IFAC PapersOnLine 52-15 (2019) 495–500

Table 1. Parameters values Par. ρ1 E2 v3 k 1 2 3

Value 1.633 × 103 9 × 107 0.3 8.53 × 10−14

Unit kg/m3 Pa 1 m3 s/kg

Par. λ σe D

Value 16.6 × 10−9 3.274 × 10−3 1.375 × 10−11

Unit m2 /(V s) 1/(Ωm) m2 /s

Material density. Young’s modulus. Poisson ratio.

Tip displacement

2

10

12

-3

Current

10

Experimental

1.5

Simulation

Amplitude

Displacement (mm)

8 1

0.5

6

4

2 0

-0.5

Experimental Simulation

0

5

10

15

20

25

0

-2

0

1

Time (seconds)

2

3

4 Time (seconds)

5

6

7

8

Fig. 6. Displacement of the endpoint of the IPMC and current. performance is investigated numerically. The Lagrangian multiplier method is used to model the geometric constraints between the gel and the beam. The global system forms a stokes-Dirac structure, that guarantees the energy preservation. This system is later discretized by means of the finite difference method on staggered grids. Thus, the model can be reduced into a descriptor port-Hamiltonian form with the elimination of the Lagrangian multiplier. Finally, the proposed model has been validated by the experimental measurements. The ongoing work is to deal with the modeling of the 2-D tubular IPMC actuator. Also, the passivity based control design for the IPMC actuator would be investigated under the port-Hamiltonian framework in the future. REFERENCES

Nishida, G., Takagi, K., Maschke, B., and Osada, T. (2011). Multi-scale distributed parameter modeling of ionic polymer-metal composite soft actuator. Control Engineering Practice, 19(4), 321–334. Park, K., Yoon, M.K., Lee, S., Choi, J., and Thubrikar, M. (2010). Effects of electrode degradation and solvent evaporation on the performance of ionic-polymer–metal composite sensors. Smart Materials and Structures, 19(7), 075002. Shahinpoor, M. (ed.) (2016). Ionic Polymer Metal Composites (IPMCs), volume 1 of Smart Materials Series. The Royal Society of Chemistry. doi: 10.1039/9781782622581. Trenchant, V., Ramirez, H., Le Gorrec, Y., and Kotyczka, P. (2018). Finite differences on staggered grids preserving the port-hamiltonian structure with application to an acoustic duct. Journal of Computational Physics, 373, 673–697. Villegas, J.A. (2007). A port-hamiltonian approach to distributed parameter systems. Wu, Y., Hamroun, B., Le Gorrec, Y., and Maschke, B. (2014). Port hamiltonian system in descriptor form for balanced reduction: Application to a nanotweezer. IFAC Proceedings Volumes, 47(3), 11404–11409. Xiao, Y. and Bhattacharya, K. (2001). Modeling electromechanical properties of ionic polymers. In Smart Structures and Materials 2001: Electroactive Polymer Actuators and Devices, volume 4329, 292–301. International Society for Optics and Photonics. Yamaue, T., Mukai, H., Asaka, K., and Doi, M. (2005). Electrostress diffusion coupling model for polyelectrolyte gels. Macromolecules, 38(4), 1349–1356. Zhu, Z., Chen, H., Wang, Y., and Li, B. (2012). Multiphysical modeling for electro-transport and deformation of ionic polymer metal composites. In Electroactive Polymer Actuators and Devices (EAPAD), volume 8340. Int. Soc. for Optics and Photonics.

Bao, X., Bar-Cohen, Y., and Lih, S.S. (2002). Measurements and macro models of ionomeric polymer-metal composites (ipmc). In Smart Structures and Materials 2002: Electroactive Polymer Actuators and Devices (EAPAD), volume 4695, 220–228. International Society for Optics and Photonics. Branco, P.C., Lopes, B., and Dente, J. (2012). Nonuniformly charged ionic polymer–metal composite actuators: Electromechanical modeling and experimental validation. IEEE Trans. on Ind. Elec., 59(2), 1105–1113. De Gennes, P., Okumura, K., Shahinpoor, M., and Kim, K.J. (2000). Mechanoelectric effects in ionic gels. EPL (Europhysics Letters), 50(4), 513. Le Gorrec, Y., Zwart, H., and Maschke, B. (2005). Dirac structures and boundary control systems associated with skew-symmetric differential operators. SIAM journal on control and optimization, 44(5), 1864–1892. Maschke, B. and van der Schaft, A. (1992). Port-controlled hamiltonian systems: Modelling origins and systemtheoretic properties. IFAC Proceedings Volumes, 25(13), 359 – 365. 2nd IFAC Symposium on Nonlinear Control Systems Design 1992, Bordeaux, France, 24-26 June. Nishida, G., Takagi, K., Maschke, B., and Luo, Z.w. (2008). Multi-scale distributed port-hamiltonian representation of ionic polymer-metal composite. In Proceedings of the 17th IFAC World Congress, volume 17. 1247