Linear finite-difference bond graph model of an ionic polymer actuator

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.


Introduction
Over the last two decades, electroactive polymers have been extensively studied [1][2][3] and developed due to their huge potential in terms of actuators, soft sensors, and for applications in robotics [4], biomedical science [5][6][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][10][11][12][13][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,4ethylenedioxythiophene) (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 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 nonuniform 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.

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 (sA,  ) for the mechanical domain, and (E, I) and (u, i) for the electrochemical domain, respectively. Variables (F,  v) and (sA,  ) 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].

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 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 x t d , as function of ( ) u x t , and thereby we find: In the same way, by applying Kirchhoff's voltage law, we end up with a similar relation with voltage: 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.

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 r x F A v , , , , , y M 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   motion can then be written as follows: , , , .

. 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: , 0.

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

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

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 mechanicalmomentum 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: 13 is the strain caused by the insertion (expulsion) of ions. It can be expressed as follows: where r c is the charge transfer per unit volume and a 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.

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 spincoated 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 (trifluoromethylesulfonyl) 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.

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).
According with that equivalent circuit, the current during charging ( ) i x t , 1 can be related to voltage as follows: where t is the time constant of the circuit and ( ) u x t , is the applied voltage. + R R e 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 W and for the capacitive part a value of 2.09 mF.

Coupling coefficient
The coupling coefficient a, also called charge-to-strain ratio, is expressed as the ratio of the peak-to-peak magnitude of the strain DS and the peak-to-peak magnitude of the charge  The bending strain  b is expressed in per cent and calculated from the displacement as follows [25]: 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 × 10 −10 m 3 C −1 .

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 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: 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  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  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  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  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: where C and f depend on the initial conditions, and g and w a are the damping ratio and the pseudo pulsation, respectively.
The damping ratio can then be expressed as follows: 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 −1 s −1 for a beam density of 1600 kg m −3 .
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 temperature in [26,27]. In the next part, the simulation will be performed using measured parameters on the first actuator.

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 t 0, 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 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 prediction 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 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 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.

Power and efficiency study
The results exposed in figures 11 and 12 validate the multiphysics model for small deflections of the ionic actuator. This section is now devoted to a power analysis both in electrochemical 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 electrochemically stored inside the equivalent capacitance C eq near location m = x 0 m and time = t 0 s, the stored electrochemical 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 μ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 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.

Conclusion
In this paper, a simple one-dimensional electro-chemomechanical 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 comparing 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.