# Linear Finite-Difference Bond Graph Model of an Ionic Polymer Actuator

Mohamed Bentefrit, Sébastien Grondel, Caroline Soyer, Adelyne Fannir, Eric Cattan, John D.W. Madden, Giao T.M. Nguyen, Cedric Plesse, Frederic Vidal

► To cite this version:

Mohamed Bentefrit, Sébastien Grondel, Caroline Soyer, Adelyne Fannir, Eric Cattan, et al.. Linear Finite-Difference Bond Graph Model of an Ionic Polymer Actuator. Smart Materials and Structures, 2017, 26 (9), pp.095055. 10.1088/1361-665X/aa7f7f. hal-03258700

HAL Id: hal-03258700

https://hal.science/hal-03258700

Submitted on 14 Sep 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Linear finite-difference bond graph model of an ionic polymer actuator

M Bentefrit<fnref n="1" />, S Grondel<fnref n="1" />, C Soyer<fnref n="1" />, A Fannir<fnref n="2" />, E Cattan<fnref n="1" />,
J D Madden<fnref n="3" />, T M G Nguyen<fnref n="2" />, C Plesse<fnref n="2" /> and F Vidal<fnref n="2" />

E-mail: eric.cattan@univ-valenciennes.fr

## Abstract

With the recent growing interest for soft actuation, many new types of ionic polymers working in air have been developed. Due to the interrelated mechanical, electrical, and chemical properties which greatly influence the characteristics of such actuators, their behavior is complex and difficult to understand, predict and optimize. In light of this challenge, an original linear multiphysics finite difference bond graph model was derived to characterize this ionic actuation. This finite difference scheme was divided into two coupled subparts, each related to a specific physical, electrochemical or mechanical domain, and then converted into a bond graph model as this language is particularly suited for systems from multiple energy domains. Simulations were then conducted and a good agreement with the experimental results was obtained. Furthermore, an analysis of the power efficiency of such actuators as a function of space and time was proposed and allowed to evaluate their performance.

Keywords: ionic EAP, finite difference method, bond graph

## 1. Introduction

Over the last two decades, electroactive polymers have been extensively studied [1–3] and developed due to their huge potential in terms of actuators, soft sensors, and for applications in robotics [4], biomedical science [5–7], and artificial muscles [8]. Electroactive polymers can be divided into two families: dielectric polymers and ionic polymers. The first one based on electrostatic actuation is the most predominant in the literature thanks to suited properties in terms of large displacements and the ability to work in air, but requires high voltage for activation. The second family based on a redox process causes the insertion or expulsion of ions supplied by an electrolyte into or from an active layer of the conducting polymer which in turn creates a change in volume in the same layer, and therefore a global mechanical displacement.

With the recent development of ionic actuators that work in air and host their own ionic liquid [9–14], new opportunities [8] have been opened up. In particular, a new type of

ionic actuator working in air composed of a three-layer interpenetrated architecture has been recently developed (figure 1) at the LPPI together with IEMN [15]. This specific architecture improves mechanical behavior, since there is no delamination between each layer and increases ionic conductivity within the actuator thickness [16].

For such a configuration, the central layer acts as an ionic host matrix for the electrolyte. It is formed from an interpenetrated network of poly(ethylene oxide) (PEO) and nitrile butadiene rubber (NBR) and is filled with the ionic liquid 1-Ethyl-3-methylimidazolium bis-trifluoromethylsulfonyl-imide (EMITFSI). The two remaining layers behave as working electrodes and are made of a conducting polymer—poly(3,4-ethylenedioxythiophene) (PEDOT)—semi-interpenetrated with the central matrix. The actuation mechanism relies on the volume variation of the PEDOT electrodes. Indeed, when a voltage is applied between those two electrodes, one of them is reduced while the other one is oxidized. This results in a migration of ion, in or from the host matrix, which in turns

<footnote n="1">¹Univ. Valenciennes, CNRS, Univ. Lille, Yncréa, Centrale Lille, UMR 8520—IEMN, DOAE, F-59313 Valenciennes, France</footnote>

<footnote n="2">²LPPI, EA 2528, UCP, Neuville sur Oise, France</footnote>

<footnote n="3">³Molecular Mechatronics Laboratory, University of British Columbia, Vancouver, Canada</footnote>
<figure><img src="image_2.png" /><figcaption>**Figure 1.** Ionic actuator structure. The blue layers are the PEDOT polymer electrodes whereas the yellow layer corresponds to the ionic PEO-NBR host matrix, and the gold layers to the Au electrodes required for electrical contact.</figcaption></figure>

causes the expansion of the electrode being reduced and the contraction of the electrode being oxidized. Those combined phenomena produce the global bending actuation of the device.

As far as we know, few models have been developed with such a family of actuators. Several physical domains are involved during actuation: electrochemical and mechanical, which are coupled. From the literature review, methods of modeling ionic actuators can be divided non-exhaustively into three categories. The first category called white box, contains analytical methods that require perfect knowledge of all the phenomena occurring, but in return provide an exact solution to the problem. For instance, Martinez *et al* [8] propose a theoretical model which describes electrochemical and mechanical phenomena since it relates mechanical work to flowing current and elements involved in the chemical reaction. This model is well adapted for actuators working at low frequencies but requires many physical parameters to measure and identify. The second category, gray box, only involves partial knowledge of the physical phenomena that occur and provides an approximated solution. For example, Nishida Gou *et al* who work on IPMC actuators [17], use successfully distributed port-Hamiltonian systems on multiple spatial scales and are therefore able to analyze influence of non-uniform and local properties. In the same category, we can also quote the work of Moghadam *et al* [18] who takes into account the geometric nonlinear displacement with rigid finite elements, whereas Moeinkhah *et al* [19] propose a distributed equivalent electrical circuit in order to better represent the electrochemical part. The last category, black box, also gives an approximated solution but does not require any knowledge of the mechanisms occurring. This method is often used in control applications since it prioritizes the use of simple models adapted to a given sample [20, 21].

The main issue of this paper is to understand and predict the behavior of the new ionic polymer composed of three PEDOT/PEO-NBR layers, as well as to quantify its power efficiency. Given the partial knowledge we have of electrochemical phenomena occurring and in order to cope with the challenge of space description and power analysis, we decided to develop an original and simple gray box model based

on a linear multiphysics finite difference bond graph representation. The primary benefits of such a representation are first the balance between accuracy and complexity supplied by the finite difference scheme [22] as it gives a spatial solution for both electrochemical and mechanical domains and second the schematized and simple representation of energy exchanges and storage within the system when using the Bond Graph language.

This article is organized as follows. In section 2, physical equations are briefly recalled for each domain and finite difference scheme are then developed in Bond Graph language. In section 3, physical parameters are estimated experimentally and injected into the model. Finally, the built model is used to estimate the actuator efficiency.

## 2. Presentation of the model

As usual, Bond Graph modeling requires first to establish the system Word Bond Graph. As this representation is the first step toward describing the energy map and the composition of the system, a physical analysis of the system and its power interactions are sufficient to determine the Word Bond Graph. As shown in figure 2, it is divided into three parts in our case. The first part is dedicated to the electrochemical behavior representation, the last one describes the mechanical behavior, and the middle one enables coupling between the first two. Contrary to conventional block diagrams, the inputs and outputs of each subsystem here define power variables represented by a conjugated effort-flow pair (e, f) indicated by a half arrow.

In this case, the effort and flow variables are (F, v) and ($\sigma A, \dot{\epsilon}$) for the mechanical domain, and (E, I) and (u, i) for the electrochemical domain, respectively. Variables (F, v) and ($\sigma A, \dot{\epsilon}$) represent mechanical force and velocity or strain rate, respectively, while (E, I) and (u, i) correspond to electrical voltage and current, respectively.

Each Word Bond Graph can be developed using elements that can either produce or store and dissipate energy, as shown figure 2. Elements relating to storage are usually defined by *C* for potential energy and *I* for kinetic energy, while dissipative elements are described by *R*. Effort and flow sources for energy creation are represented by *Se* and *Sf* elements, respectively. Finally, elements such as ‘0’ and ‘1’ junctions or TF and GY elements can be used for transmitting energy without loss [23].

### 2.1. Electrochemical aspect

As related in [5], when driven at low frequencies by a voltage *u*(*x*, *t*), conducting polymers act as electrical capacitors. The three layers' actuator can be thus electrically represented by an RC circuit with resistances that translate the electronic and ionic transport. The resistance part will be provided by surface and thickness resistance of the PEDOT on one hand, and by the thickness resistance of the NBR-PEO hosting matrix on a second hand. We will also consider in some case a third resistance when leakage in current *i*(*x*, *t*) is observed, this
<figure><img src="image_3.png" /><figcaption>**Figure 2.** Word bond graph of the global system. Boxes represent physical domains. Parameters above and below arrows represent efforts and flows provided or received by those domains.</figcaption></figure>

<figure><img src="image_4.png" /><figcaption>**Figure 3.** Equivalent circuit of an elementary actuator of length dx.</figcaption></figure>

shunt resistance will describe thus local short circuits that can occur between the two PEDOT electrodes. The resistance per unit length of PEDOT is $R_e$, the resistance per unit length of NBR-PEO is $R_{EHM}$, the shunt resistance per unit length is $R_{cc}$, and the equivalent capacitance per unit length is $C_{eq}$. We assume that PEDOT is electronically isotropic so we can then express the same resistance per unit length along the surface or the thickness.

From the electrical diagram, figure 3, we can then determine the partial differential equations (PDEs) governing current and voltage behavior in the circuit. Using voltage divider law we can express $u(x + dx, t)$ as function of $u(x, t)$ and thereby we find:

$$ \frac{\partial i(x, t)}{\partial x} = -u(x, t) * \frac{R_{cc} + R_{EHM} + \frac{1}{jC_{eq}\omega}}{2R_e \left( \frac{1}{jC_{eq}\omega} + R_{EHM} + R_{cc} \right) + R_{cc} \left( R_{EHM} + \frac{1}{jC_{eq}\omega} \right)}. \quad (2.1) $$

In the same way, by applying Kirchhoff's voltage law, we end up with a similar relation with voltage:

$$ \frac{\partial u(x, t)}{\partial x} = -i(x, t) * R_e. \quad (2.2) $$

By combining (2.1) with (2.2) and assuming that $dx \to 0$, we can then express the two PDEs that govern current and voltage in an elementary equivalent circuit:

$$ \frac{\partial^2 i(x, t)}{\partial x^2} - i(x, t) * \frac{Z_2}{Z_1} = 0, \quad (2.3) $$

$$ \frac{\partial^2 u(x, t)}{\partial x^2} - u(x, t) * \frac{Z_2}{Z_1} = 0. \quad (2.4) $$

<figure><img src="image_5.png" /><figcaption>**Figure 4.** Bond graph of an electrochemical cell.</figcaption></figure>

By approximating (2.3) and (2.4) using the forward finite difference method, the Bond Graph scheme of an electrochemical actuator part can be rewritten.

As illustrated in figure 4, the schematic established, provides a micro-scale visual representation of the power transfer between the resistive and capacitive parts of the electrochemical element. Replication and serial linking of the same element enable the rapid and easy visualization of the same power transfer on a macro-scale.

## 2.2. Mechanical aspect

From a mechanical point of view, the polymer actuator is assimilated to a beam with a rectangular cross section.

In order to simplify the study and with regards to the electrochemical phenomena involved, it is assumed that the actuator only generates bending deformations. We will also put forth the hypothesis that in our field of study the material remains viscoelastic and linear, i.e. that no creeping occurs, and that it can be seen as a homogenous, isotropic material with an equivalent Young's modulus $E$ and an area moment of inertia $I$ along the $x$ axis. Let us name $F$, $\rho$, $A$, $v$, $\xi$, $M$, $\psi$, $k$, and $G$, the external applied force, the beam density, the cross-sectional area, the transverse displacement, the damping, the external applied momentum, the angular displacement, the numerical factor related to the cross section shape, and the shear modulus, respectively. The equations of
<figure><img src="image_6.png" /><figcaption>Figure 5. Bond graph of a mechanical cell.</figcaption></figure>

The Bond graph schematic of the mechanical element can be then developed from the previous approximations, as shown figure 5.

## 2.3. Assembly of the global model

The expulsion of ions from the ionic liquid and insertion into the two PEDOT conducting polymer layers causes one layer to expand and the other to contract, resulting in mechanical-momentum and the displacement of the actuator. It then appears, as reported by Tina Shoa et al [5], that the displacement of the actuator (or the momentum generated) is directly related to the charges flowing inside the conducting polymer. Equation (2.8) can be rewritten as follows:

$$M(x, t) = EI \frac{\partial\psi(x, t)}{\partial x} + EA\Delta x\epsilon(x, t). \quad (2.13)$$

Equation (2.13) is called the electro-chemo-mechanical coupling equation, where $\epsilon(x, t)$ is the strain caused by the insertion (expulsion) of ions. It can be expressed as follows:

$$\epsilon(x, t) = \alpha\rho_c(x, t), \quad (2.14)$$

where $\rho_c$ is the charge transfer per unit volume and $\alpha$ is an empirically determined ratio that expresses the strain generated per volumetric charge transfer.

From equations (2.13) and (2.14) we can then transform the global model to a Bond Graph for a given elementary actuator of length $dx$, see figure 6.

motion can then be written as follows:

$$\frac{\partial F(x, t)}{\partial x} = \rho A \frac{\partial^2 v(x, t)}{\partial t^2} + \xi \frac{\partial v(x, t)}{\partial t}, \quad (2.5)$$

$$F(x, t) + \frac{\partial M(x, t)}{\partial x} = \rho I \frac{\partial^2 \psi(x, t)}{\partial t^2}, \quad (2.6)$$

$$F(x, t) = kGA \left( \frac{\partial v(x, t)}{\partial x} - \psi(x, t) \right), \quad (2.7)$$

$$M(x, t) = EI \frac{\partial \psi(x, t)}{\partial x}. \quad (2.8)$$

With the additional assumptions that there is no inertial rotation of the cross sections and no shear force (Euler-Bernoulli assumptions), equation (2.7) disappears and equation (2.6) is simplified as follows:

$$F(x, t) + \frac{\partial M(x, t)}{\partial x} = 0. \quad (2.9)$$

Obtained PDEs can be approximated with the forward finite difference method:

$$\frac{F(x + \Delta x, t) - F(x, t)}{\Delta x} = \rho A \frac{\partial^2 v(x, t)}{\partial t^2} + \xi \frac{\partial v(x, t)}{\partial t}, \quad (2.10)$$

$$F(x + \Delta x, t) + \frac{M(x + \Delta x, t) - M(x, t)}{\Delta x} = 0, \quad (2.11)$$

$$M(x + \Delta x, t) = EI \frac{\psi(x + \Delta x, t) - \psi(x, t)}{\Delta x}. \quad (2.12)$$

# 3. Estimation of the physical parameters

## 3.1. Beams synthesis and preparation

In order to simulate the ionic actuator, all the physical parameters of the finite difference Bond graph model must be fully determined. In this way, a beam measuring 9375 µm × 2000 µm × 200 µm was initially synthesized [24] and shaped from a bulk PEDOT-PEO/NBR material as follows:

- The bulk material was as described in [24] and was spin-coated onto a rigid silicon wafer with a polyvinyl alcohol (PVA) acting as an adhesive and a sacrificial layer.

- The beam was then cut directly on the silicon wafer using an infrared nano-pulsed laser (E-Series, OXFORD Laser Ltd). With this cutting system, the beam can be shaped to the required dimensions with a resolution of around one micrometer. The luminous power of the laser beam was adjusted so that the etching did not go through the silicon wafer.

- The beam was then released from the wafer by dissolving the PVA layer in deionized water for 10 min.

- Once released, the beam was dried by blowing nitrogen over it and swollen with ionic liquid 1-Ethyl-3-methylimidazolium bis (trifluoromethylethylsulfonyl) imide (SOL-VIONIC) at ambient temperature for over two weeks.

The aim is now to extract the physical characteristics i.e. both mechanical and electrochemical from this sample.
<figure><img src="image_7.png" /><figcaption>**Figure 6.** The electro-chemo-mechanical element of length dx. Mechanical and electrochemical elements are coupled through a modulated source element which is the translation of equation (2.13).</figcaption></figure>

## 3.2. Electrochemical parameters

All the electrochemical parameters can be determined from a simple current reading with an applied square voltage of ±1.5 V swept at a frequency of 0.1 Hz. Notice here that it was decided to clamp the entire sample between two glass wafers covered with a 200 nm thin gold layer in order to eliminate the surface resistance effect.

The readings, figure 7(a), show a charging and discharging phase, during positive and negative square voltage. For each phase, the current goes back to zero, which means that there is no current leakage due to the large thickness of the material and low voltage used. As main assumption, we can now neglect the $R_{cc}$ element for this experiment and therefore the equivalent circuit can be simplified as in figure 7(b).

<figure><img src="image_8.png" /></figure>

## b) Equivalent circuit

<figure><img src="image_9.png" /><figcaption>**Figure 7.** (a) Current flowing through the thickness of the three layers actuator PEDOT/PEO-NBR under ±1.5 V squared voltage driven at 0.1 Hz, (b) simplified equivalent circuit.</figcaption></figure>

According with that equivalent circuit, the current during charging $i_1(x, t)$ can be related to voltage as follows:

$$i_1(x, t) = \frac{C_{eq}u(x, t)}{\tau} \exp\left(-\frac{t}{\tau}\right), \quad (3.1)$$

where $τ$ is the time constant of the circuit and $u(x, t)$ is the applied voltage. $R_e + R_{EHM}$ can then be estimated from the maximal current magnitude during the charging phase whereas $C_{eq}$ can be determined from the slope of the current curve during the same phase. It must then be verified that the current reading matches an RC charging or discharging circuit.

Figure 7(a) shows that the current signal follows an exponential behavior and hence can be approximated by equation (3.1). We can then extract $R_e$, $R_{EHM}$, and $C_{eq}$, we find for the resistive part a total value of 645 Ω and for the capacitive part a value of 2.09 mF.

## 3.3. Coupling coefficient

The coupling coefficient $α$, also called charge-to-strain ratio, is expressed as the ratio of the peak-to-peak magnitude of the strain $ΔS$ and the peak-to-peak magnitude of the charge
<figure><img src="image_10.png" /><figcaption>Figure 8. (a) Transferred charge to the three layers of the PEDOT/ PEO-NBR actuator with a square voltage of ±1.5 V driven at 0.1 Hz, (b) produced strain under ±1.5 V squared voltage swept at 0.1 Hz.</figcaption></figure>

transferred $\Delta Q$.

$$ \alpha = \frac{\Delta S}{\Delta Q}. \qquad (3.2) $$

The bending strain $\epsilon_b$ is expressed in per cent and calculated from the displacement as follows [25]:

$$ \epsilon_b = 100 \frac{hv}{l^2 + v^2}. \qquad (3.3) $$

With $v$ being the displacement, $w$ the width and $l$ the length.

From the previous charge curve and displacement curve under ±1.5 V, see figures 8(a) and (b), we find for the strain to charge ratio a value of $4.56 \times 10^{-10} \text{ m}^3 \text{ C}^{-1}$.

### 3.4. Mechanical characterization

In the case of mechanical characterization, Young's modulus can be measured in several ways, each introducing a different level of accuracy and information. The traditional method is

<figure><img src="image_11.png" /><figcaption>Figure 9. (a) Stiffness-measurement bench, (b) beam stiffness.</figcaption></figure>

tensile testing, which also determines the strain hardening property of the test material. Young's modulus can also be estimated by vibrating the beam at its resonance frequency or acoustically by measuring the celerity of a known acoustic wave travelling inside the material. As our material is highly damped, the latter is not appropriate. The tensile test could not be considered either as it is difficult to adapt small-geometry samples for the tensile test bench.

The most suitable method we have found so far is the flexure test. Indeed, when a stress is applied to the beam that is below the yield stress of the material, according to Hooke's law, the force applied is proportional to displacement as follows:

$$ F(x, t) = kv(x, t), \qquad (3.4) $$

where $k$ is the stiffness of the material and $v(x, t)$ the resulting displacement under the applied force $F$. We can then relate Young's modulus $E$ to stiffness $k$ through the width $w$ and length $l$ of the beam, through equation (3.5) which stays valid only for cantilever beams.

$$ E = \frac{12kl^3}{3wh^3}. \qquad (3.5) $$

Measurements were performed with a FemtoTools FT-S microforce probe that can detect forces down to 1 µN, see figure 9. The probe is mounted on an x, y, z piezo-actuated platform and driven by software that allows us to run
<figure><img src="image_12.png" /><figcaption>Figure 10. (a) Displacement measurement bench, (b) free tip motion after sinusoidal excitation. The zoom-in-box shows a pseudo period $T_a$ of the signal.</figcaption></figure>

automated measurements. The probe, which can be positioned with extremely high precision, is initially placed in contact with the beam and moved 100 µm along $\vec{y}$ in 10 µm increments. Every 10 µm, the probe reads the force and the values are then used to plot the force-displacement curve, see figure 9(b). The curve shows that the data are perfectly spaced along $\vec{y}$ and that the points are well aligned. This implies that the resulting stress was way below the yield stress.

In order to obtain the damping ratio of the material, the beam was clamped to a mechanical shaker (TIR-Vibration Test System TV 500018) and excited by an oscillating force. Once the excitation stopped, the tip displacement was determined using a laser system Keyence LKG32, see figure 10(a). From the logarithmic decrement, see figure 10(b), the damping ratio can be determined.

Indeed, in a spring-damper system the displacement $v(t)$ can be described as shown below:

$$v(t) = C \cdot \exp(-\gamma \cdot t) \cdot \cos(\omega_a t + \phi), \quad (3.6)$$

where C and $\phi$ depend on the initial conditions, and $\gamma$ and $\omega_a$ are the damping ratio and the pseudo pulsation, respectively.

The damping ratio can then be expressed as follows:

$$\gamma = - \frac{\ln\left(\frac{A_2}{A_1}\right)}{T_a} \quad (3.7)$$

With $T_a$ being the pseudo period, and $A_1$ and $A_2$ two displacement points away from $T_a$. For a beam with dimensions of 9375 µm × 2000 µm × 200 µm, applying the previously described techniques we find a young modulus of 200 MPa ± 20 MPa, and a damping equal to 0.38 N m⁻¹ s⁻¹ for a beam density of 1600 kg m⁻³.

All physical parameters for the specimen with dimensions of 9375 µm × 2000 µm × 200 µm are recalled in table 1. In order to focus on the uncertainties that may occur between two specimens from a same batch, another actuator with a different width was also characterized whereas its physical parameters were deduced from the first actuator by taking into account this width difference. Results show a slight deviation on the electrochemical part between measurements and estimated values. This means that the electrochemical parameters contrarily to the mechanical parameters seem sensitive to ambient conditions. This phenomenon was already observed by the authors who analyzed physical parameters variations as function of
<figure><img src="image_13.png" /><figcaption>**Figure 11.** (a) Simulated current *i*(*x*, *t*) along the ionic beam location, (b) simulated current versus experimental current at clamping point *i*(*x* = 0, *t*), (c) simulated transverse beam displacement *v*(*x*, *t*) as function of time, (d) simulated displacement versus experimental displacement at free end tip *v*(*x* = 9375 µm, *t*). Entry voltage is set at a ±1.5 V swept at 0.1 Hz.</figcaption></figure>

<table><caption>**Table 1.** Measured parameters on 9375 × 2000 × 200 beam, measured versus estimated parameters on 10 000 × 3000 × 200 beam.</caption><thead><tr><th></th><th colspan="3">Measured</th><th colspan="3">Estimated</th></tr></thead><tbody><tr><td>Geometry (µm³)</td><td>9375</td><td>2000</td><td>200</td><td>10 000</td><td>3000</td><td>200</td></tr><tr><td>Resistance (Ω)</td><td colspan="3">645</td><td colspan="3">400</td></tr><tr><td>Capacitance (mF)</td><td colspan="3">2.09</td><td colspan="3">5.53</td></tr><tr><td>Strain to charge ratio (m³C⁻¹)</td><td colspan="3">4.56 × 10⁻¹⁰</td><td colspan="3">1.58 × 10⁻¹⁰</td></tr><tr><td>Young modulus (MPa)</td><td colspan="3">200</td><td colspan="3">200</td></tr><tr><td>Damping (N m⁻¹ s⁻¹)</td><td colspan="3">0.38</td><td colspan="3">0.38</td></tr><tr><td>Density (kg m⁻³)</td><td colspan="3">1600</td><td colspan="3">1600</td></tr></tbody></table>

temperature in [26, 27]. In the next part, the simulation will be performed using measured parameters on the first actuator.

# 4. Results

## 4.1. Time domain responses

The aim of this section is to compare the results of the simulated time domain responses with those of the experiments in order to assess the performance of the proposed modeling approach. To simulate the actuator behavior, it has been chosen to use ten finite difference elements in the bond

graph model to achieve a balance between accuracy and ease of calculation. Moreover, the simulated actuator was driven in the same conditions than in the experiments one, i.e. by a ±1.5 V square voltage at a frequency of 0.1 Hz in air and at ambient temperature.

Both electrochemical parts and mechanical parts are analyzed in this study. Hence, figure 11(a) represents the simulated electrical current *i*(*x*, *t*) along the ionic actuator length whereas figure 11(b) shows a comparison between the simulated and experimental current *i*(*0*, *t*) at the clamping point, i.e. where the gold electrodes are located. In a same way, figure 11(c) illustrates the simulated transverse beam displacement *v*(*x*, *t*) at different locations along the ionic
<figure><img src="image_14.png" /><figcaption>**Figure 12.** Top: simulated current versus experimental current at clamping point *i*(*x* = 0, *t*) under various voltages. Bottom: maximum simulated displacement versus maximum experimental displacement at free end tip *v*(*x* = 6400 µm, *t*), under various voltages. Entry voltages are squared and set at ±0.56 V, ±1.22 V, ±1.47 V and ±2.06 V and swept at 0.05 Hz.</figcaption></figure>

actuator length whereas figure 11(d) shows a comparison
between the simulated and experimental transverse beam
displacement *v*(*x*, *t*) at the tip position, i.e. for *x* = 9375 µm.

First of all, results show a close matching between the
proposed model and the experimental data. Hence, it can be
seen from figure 11(d) that the obtained simulated and
experimental displacements at the ionic beam extremity for a
±1.5 V square voltage excitation and a frequency of 0.1 Hz
are in good agreement. Furthermore, it can also be noticed in
figure 11(b) that the simulated current slope fits well with the
experimental one, despite a discard of 1.8 mA noticed on the
maximum value. This discrepancy in the current pic values is
mainly due to the time discretization in our numerical model,

but is not sufficient to affect the deflection results of the ionic
actuator beam.

To estimate the validity of the proposed model, another
simulation is proposed on a similar sample but with different
dimensions (length of 6400 µm), voltages amplitude and
frequency i.e. with an entry square voltage of ±0.56 V,
±1.22 V, ±1.47 V and ±2.06 V and a frequency of 0.05 Hz.
Correlation between the displacements and the currents on
figure 12 shows that the model is able to make a good pre-
diction for samples with different geometries at different
voltages and frequencies.

Additionally, one of the advantages of our modeling
approach using finite difference element is to allow the direct
<figure><img src="image_15.png" /></figure>

<figure><img src="image_16.png" /></figure>

<figure><img src="image_17.png" /></figure>

<figure><img src="image_18.png" /></figure>

<figure><img src="image_19.png" /></figure>

<figure><img src="image_20.png" /><figcaption>Figure 13. (a) Injected electrical power $P_i(x, t) = u(x, t) \cdot i(x, t)$. (b) Injected mechanical power. (c) Stored electrochemical power. (d) Stored mechanical power. (e) Dissipated electrochemical power. (f) Dissipated mechanical power.</figcaption></figure>

simulation of the ionic actuator behavior as function of the space. Therefore, it can be noticed from figure 11(a) and as expected that the simulated electrical current was not absorbed in the same way along the ionic actuator length due to the

surface resistance. Furthermore, figure 11(c) illustrates the linear deflection of the ionic actuator. It must be mentioned here that despite the fact that our model is not developed for the case of large deflection, it could be adapted for such a case
<figure><img src="image_21.png" /><figcaption>**Figure 14.** Evolution of total injected power and energy (in blue), total stored power and energy (in green) and total dissipated power and energy (in red), for the electrochemical domain (a), (c) and mechanical domain (b), (d).</figcaption></figure>

using similar modeling than in [18] but this issue requires at the same time to experimentally characterize the nonlinearity of the polymer behavior in case of high voltage. As this part is out of the scope of this paper since we concentrate us on the power efficiency and exchanges inside the ionic actuator, it will not be treated here.

## 4.2. Power and efficiency study

The results exposed in figures 11 and 12 validate the multi-physics model for small deflections of the ionic actuator. This section is now devoted to a power analysis both in electro-chemical and mechanical parts, in order to determine the global actuator efficiency.

As illustrated in figure 13(a), the power electrically conducted through PEDOT electrodes is non-uniformly transmitted along the ionic actuator length due to the electrical surface resistance. A large decrease is observed which means that for a better efficiency of the actuator, the electrical

electrodes need to be deposited all along the ionic actuator. Then, as seen figure 13(c), a part of power is electro-chemically stored inside the equivalent capacitance $C_{eq}$ near location $x = 0 \mu m$ and time $t = 0 s$, the stored electro-chemical power reaching a minimum of -1.11 mW. At that location, this means that only a part of the injected power takes places in the redox process, the other part being electrochemically dissipated through electrical and ionic resistance as seen figure 13(e). This dissipated power achieves a maximum value of 3.27 mW at $x = 0 \mu m$ and $t = 0 s$, which demonstrates that the electrochemical system is highly dissipative. Then a part of the stored electrochemical power is converted through the coupling (2.13) and used to generate mechanical displacement. Again, a part of the injected mechanical power figure 13(b) is stored figure 13(d) whereas the other part is dissipated figure 13(f). Notice that negative electrochemical power represents here the release of energy previously stored in the electrochemical parts, whereas the negative mechanical power is related to sign reversal of
<figure><img src="image_22.png" /><figcaption>**Figure 15.** Instantaneous efficiencies in (a) electrochemical system, (b) mechanical system, (c) global system.</figcaption></figure>

inertial forces and energy storage in elastic part of actuator for a short period.

Figures 14(a) and (b) illustrate the total powers, i.e. the input, dissipated and stored powers as function of the time of the electrochemical and mechanical parts of the ionic actuator respectively. From these curves, it was also possible to determine the input, dissipated and stored energies as function of time for both the electrochemical (figure 14(c)) and mechanical (figure 14(d)) domains.

The obtained responses reflect the behavior to a transient excitation of a multiple RC circuit and a multiple RLC circuit for the electrochemical domain and for the mechanical domain respectively. Indeed figure 14(a) shows that the stored electrochemical power increases first between 0 and 5 s (half a period) which means that electronics charge are getting accumulated inside the PEDOT layers, while input and dissipated powers are getting lower because of the decrease in current due to the charging capacitor phenomena. When the input voltage switches sign at $t = 5$ s, the stored power drops

below zero at -2.8 mW and therefore the stored energy is released to the system, the capacitors being discharging. Furthermore, we notice that the overcurrent is totally dissipated until the stored power goes above zero, at $t = 6.3$ s. Figure 14(b) shows that mechanical powers behave the same way except for the dissipative mechanical power which is really low since we are exciting the actuator at a frequency far from the mechanical eigenfrequency: the damping has here only a minimal effect at that frequency. From electrochemical and mechanical energies in figures 14(c) and (d), we come to the conclusions that mean stored energies are constant. This is because the energies are released back to the system in time with the difference that in mechanical domain the stored energy is greater due to the small dissipative part.

Figure 15(c) illustrates global efficiency i.e. the ratio between stored mechanical energies and input electrochemical energies, whereas figures 15(a) and (b) show electrochemical and mechanical efficiencies respectively, i.e. the ratio between stored and input energies. From the study of
these figures, it can be concluded that the overall maximal
efficiency is at about 0.01% and that energy dissipation
mainly occurs in the electrochemical part. In fact, maximal
electrochemical efficiency is equal to 7% (this is much lower
than a classical RC circuit at 50%) whereas maximal
mechanical efficiency is much greater at about 75% due to the
minimal effect of mechanical damping at this frequency.

**5. Conclusion**

In this paper, a simple one-dimensional electro-chemo-
mechanical model using a bond graph finite difference
method, was developed in order to simulate the behavior of a
new kind of ionic actuator. It is distinguished from existing
modeling work of ionic polymer actuators in that it is
amenable to time responses and energy exchange simulations
while capturing all the multi-physics domains. The simulated
times responses of the linear model were verified by com-
paring them with experimental results.

In summary, the built model enables us thus to predict
behavior and understand phenomena occurring in the ionic
actuator while visualizing the power and energy efficiencies
along the actuator length and over time. It is our aim that the
model presented in this paper may be used and developed
further for the use and design of ionic polymer actuators.

Future work will be focused on two aspects. First, the
proposed actuation model will be extended to incorporate
electrochemical and mechanical nonlinearities that become
pronounced at large actuations levels. The nonlinearities
include geometric nonlinearity, nonlinear elasticity and the
dependence of electrochemical parameters on the curvature
output. The second direction of future work is the application
of the modeling approach for the sensing case. There the
physics of the model has to be extended to account for force
interaction measurements in the presence of external objects.

Acknowledgments

This project received funding from the French Government
through the National Research Agency ANR-15-CE08-0032-
01 MICRO-TIP, under program PIA EQUIPEX LEAF
(ANR-11-EQPX-0025) and RENATECH.

ORCID iDs

M Bentefrit https://orcid.org/0000-0003-0888-8182

S Grondel https://orcid.org/0000-0003-1542-9494

E Cattan https://orcid.org/0000-0002-6576-5800

J D Madden https://orcid.org/0000-0002-0014-4712

References

[1] Smela E, Inganas O and Lundstrom I 1993 Conducting polymers as artificial muscles: challenges and possibilities
*J. Micromech. Microeng.* **3** 203

[2] Pei Q and Inganäs O 1993 Electrochemical muscles: bending strips built from conjugated polymers *Synth. Met.* **57** 3718–23

[3] Pei Q, Inganas O and Lundstrom I 1993 Bending bilayer strips built from polyaniline for artificial electrochemical muscles *Smart Mater. Struct.* **2** 1

[4] Jager E W H, Inganäs O and Lundström I 2000 Microrobots for micrometer-size objects in aqueous media: potential tools for single-cell manipulation *Science* **288** 2335–8

[5] Shoa T, Madden J D W, Munce N R and Yang V 2009 Analytical modeling of a conducting polymer-driven catheter *Polym. Int.* **59** 343–51

[6] Jager E W, Immerstrand C, Peterson K H, Magnusson K-E, Lundström I and Inganäs O 2002 The cell clinic: closable microvials for single cell studies *Biomed. Microdevices* **4** 177–87

[7] Low L M, Seetharaman S, He K Q and Madou M J 2000 Microactuators toward microvalves for responsive controlled drug delivery *Sensors Actuators* B **67** 149–60

[8] Martinez J G and Otero T F 2014 Mechanical awareness from sensing artificial muscles: experiments and modeling *Sensors Actuators* B **195** 365–72

[9] Plesse C, Vidal F, Randriamahazaka H, Teyssié D and Chevrot C 2005 Synthesis and characterization of conducting interpenetrating polymer networks for new actuators *Polymer* **46** 7771–8

[10] Alici G, Devaud V, Renaud P and Spinks G 2009 Conducting polymer microactuators operating in air *J. Micromech. Microeng.* **19** 025017

[11] Gaihre B, Alici G, Spinks G M and Cairney J M 2011 Synthesis and performance evaluation of thin film PPy-PVDF multilayer electroactive polymer actuators *Sensors Actuators* A **165** 321–8

[12] Jeon J-H, Oh I-K, Kee C-D and Kim S-J 2010 Bacterial cellulose actuator with electrically driven bending deformation in hydrated condition *Sensors Actuators* B **146** 307–13

[13] Panwar V, Kang B-S, Park J-O and Park S-H 2011 New ionic polymer-metal composite actuators based on PVDF/PSSA/PVP polymer blend membrane *Polym. Eng. Sci.* **51** 1730–41

[14] Kim S S, Jeon J H, Kee C D and Oh I K 2013 Electro-active hybrid actuators based on freeze-dried bacterial cellulose and PEDOT: PSS *Smart Mater. Struct.* **22** 085026

[15] Goujon L J, Khaldi A, Maziz A, Plesse C, Nguyen G T M, Aubert P-H, Vidal F, Chevrot C and Teyssié D 2011 Flexible solid polymer electrolytes based on nitrile butadiene rubber/poly(ethylene oxide) interpenetrating polymer networks containing either LiTFSI or EMITFSI *Macromolecules* **44** 9683–91

[16] Khaldi A 2012 *PhD Thesis* University of Valenciennes

[17] Nishida G, Takagi K, Maschke B and Osada T 2011 Multi-scale distributed parameter modeling of ionic polymer-metal composite soft actuator *Control Eng. Pract.* **19** 321–34

[18] Moghadam A A A, Hong W, Kouzani A, Kaynak A, Zamani R and Montazami R 2014 Nonlinear dynamic modeling of ionic polymer conductive network composite actuators using rigid finite element method *Sensors Actuators* A **217** 168–82

[19] Amiri-M A-A, Moavenian M and Torabiz K 2010 Takagi-Sugeno fuzzy modelling and parallel distributed compensation control of conducting polymer actuators *Syst. Control Eng.* **224** 41–51

[20] John S W 2008 *PhD Thesis* University of Wollongong

[21] John S W, Alici G and Cook C D 2010 Inversion-based feedforward control of polypyrrole trilayer bender actuators *IEEE/ASME Trans. Mech.* **15** 149–56

[22] Spiteri P 2002 *Méthode des Différences Finies Pour les EDP Stationnaires* (Paris: Techniques de l'ingénieur)

[23] Borutzky W 2011 *Bond Graph Modelling of Engineering Systems* (New York, NY: Springer)
[24] Plesse C, Khaldi A, Wang Q, Cattan E, Teyssié D, Chevrot C and Vidal F 2011 Polyethylene oxide-polytetrahydrofurane-PEDOT conducting interpenetrating polymer networks for high speed actuators *Smart Mater. Struct.* **20** 124002

[25] Sugino T, Kiyohara K, Takeuchi I, Mukai K and Asaka K 2009 Actuator properties of the complexes composed by carbon

nanotube and ionic liquid: the effects of additives *Sensors Actuators* B **141** 179–86

[26] Maziz A, Plesse C, Soyer C, Chevrot C, Teyssié D, Cattan E and Vidal F 2014 Demonstrating kHz frequency actuation for conducting polymer microactuators *Adv. Funct. Mater.* **24** 4851–9

[27] Fannir A 2017 *PhD Thesis* University of Cergy Pontoise p 63