A 3D analytical ion transport model for ionic polymer metal composite actuators in large bending deformations

Ionic polymer metal composites (IPMCs) are a kind of soft electroactive polymer composites. An IPMC strip commonly has a thin polymer membrane coated with a noble metal as electrodes on both sides. Whenever an electric voltage is applied to the IPMC, it bends and whenever it is deformed, a low voltage is measurable between its electrodes, hence IPMC is an actuator as well as a sensor. They are well known for their promising features like low density, lightness, high toughness and remarkable stimulus strain, also, they have the potential for low-voltage operation while exhibiting acceptable large bending deformation. In this paper, a three-dimensional (3D), dynamic and physics-based model is presented analytically and experimentally for IPMC actuators. The model combines the ion transport dynamics within the IPMC and the bending dynamics of it as a beam under an electrical stimulation. In particular, we present an analytical model to create a relation between the input voltage and the output tip displacement of an IPMC actuator for large bending deformations. Experimental results show that the proposed model captures well the tip displacement.

www.nature.com/scientificreports/ I Moment of inertia ω(z, t) Shape function of the IPMC κ Math (z, t) Mathematical definition for the curvature of the shape function κ Mech (z, t) Mechanical definition for the curvature of the shape function δ Tip (t) Tip displacement of the IPMC Ionic polymer metal composites (IPMCs) are one of the main ionic electro-active polymers that have potential applications especially in the field of medicine [1][2][3][4][5][6] . As depicted in Fig. 1, an IPMC strip has a thin polymer membrane (e.g. Nafion) coated with a noble metal as electrodes (e.g. Pt) on both sides (anode and cathode). A relatively low voltage applied through the thickness of an IPMC causes the internal hydrated cations to transport toward the cathode. This transportation leads to the Nafion inflation in the cathode side and consequently produces a bending response toward the anode side [7][8][9][10][11] . Reversely, we can measure a low voltage signal between its electrodes in response to an input mechanical deformation, hence IPMC is an actuator and a sensor as well.
IPMCs have received great interest in recent years for their potential applications both as sensor and actuator in a large variety of engineering areas. They are bio-compatible 5 and require low voltages (less than 5 V), produce large bending deformations, can work properly in air and in aqueous environments using appropriate encapsulation layers 6 , and have very large stimulus strain and low density 1 . As mentioned in the advantages of an IPMC actuator, it can produce large bending deformations, however, in this situation, its behavior is nonlinear [12][13][14][15] . IPMC's behavior is time-dependent and except its steady state analysis, it should be modeled dynamically which means that static models cannot be used for tracking the time-dependent performance of an IPMC. Additionally, IPMC shows nonlinear behavior in large deformation situations. Hence, it is required to use dynamic nonlinear models to reach a time-dependent model in large deformation situations of IPMC. In recent years, there have been extensive efforts in the dynamic modeling and understanding of IPMC actuators, however, they have studied small bending deformations (tip displacement of less than 1 mm) in their experiments and have not discussed large bending deformations. Although, due to the widespread use of IPMC actuators, an accurate and practical dynamic model is desirable. In general, approaches for the modeling of IPMC actuators can be divided into two main groups. The first group includes analytical modeling and the second one consists of predictive identification methods. Various analytical models have been studied for IPMC actuators. For example, distributed resistor-capacitor (RC) equivalent circuits [15][16][17][18] , physical and multi-physical approaches using partial differential equations (PDEs), finite element methods, COMSOL modeling [19][20][21] and also several miscellaneous uncategorized models 22,23 are some of the presented analytical models. Moreover, there are two approaches in the predictive identification group, classical 24,25 and intelligent methods 9,10,26 . Nevertheless, the important point is that all mentioned methods in the two main groups have not addressed the issue of large bending deformations of an IPMC actuator except 12 . On the other hand, reference 12 has proposed a different predictive identification method but this method is completely a black-box model and cannot describe the dynamic chemo-electric response of an IPMC actuator.
The contribution of the current paper is an analytical, three-dimensional (3D), dynamic and physics-based model which can describe the true electromechanical behavior of IPMC actuators even in large bending deformations. Various physics-based models have also been studied for IPMCs which can be classified into three categories, thermodynamics of irreversible process (TIP) models [25][26][27] , frictional models (FR) 20,21 and Nernst-Planck based (NP) models [28][29][30][31][32][33] . Among these three categories, NP models propose the most straightforward way to explain ion transport in IPMC sensors and actuators.
In the current paper, the starting point of the model development is the same governing PDE as in NP-based models presented until now. However, this work improves previous studies significantly in three aspects. First, it considers that the ionic content of an IPMC strip migrates through all directions and needs to solve the problem in a three-dimensional space. All the above mentioned NP-based methods are one dimensional except 28 which is a two dimensional, non-analytical and semi-numerical method. Second, the proposed model is the first dynamic and full analytical 3D physics-based ion transport model that covers the nonlinear behavior of IPMC actuators even in large bending deformations. Third, it is a promising result that for the first time we get a mathematical well-defined relationship between the input voltages applied to the IPMC and the curvature of the IPMC in large bending deformations.
Experiments were conducted to validate the proposed dynamic model for an IPMC actuator in a cantilever configuration. A good match was achieved between the measured bending response and the model prediction.
The rest of this paper consists of three main sections. In Section II we describe our proposed 3D physics-based ion transport model and its details. The estimation of the model's unknown parameters and assessment of the model accuracy are presented in Section III. The discussion along with a conclusion is provided in Section IV.

A 3D physics-based ion transport model for IPMC
Preliminary equations and obtaining the governing partial differential equations. A geometric definition of an IPMC beam is illustrated in Fig. 2. The beam is clamped at one end (z = 0) and is subject to an electrical voltage. In the following equations, the electric displacement, the electric field, the electric potential and the electric charge density are denoted, respectively, by D x, y, z, t , E(x, y, z, t), φ(x, y, z, t) and ρ x, y, z, t . The following equations hold: , κ e is the effective dielectric constant of the membrane and ρ x, y, z, t is also defined as follows [29][30][31] : In (4), F is the Faraday constant and C + x, y, z, t and C − are cation and anion concentrations in the membrane, respectively.
The continuity equation for cations is given by: that J x, y, z, t is the cation flux vector. The ion flux consists of diffusion, migration and convection terms in the Nernst-Planck PDE as follows: where d is the ionic diffusivity, R is the gas constant, T is the absolute temperature, p is the fluid pressure, v is the free solvent velocity field and V represents the volumetric change of the membrane 31,32 . Note As considered until now, in this model all variables are assumed three-dimensional time-variant. However, for ease of notation, (x, y, z, t) is omitted in the following equations, except in equations that may lead to confusion. Equation (7) is a time-variant 3D PDE that describes the electric potential of IPMC,φ x, y, z, t . For more details, we refer the reader to Appendix A.
(2) E x, y, z, t = −∇φ x, y, z, t , (3) ∇. D x, y, z, t = ρ x, y, z, t , (4) ρ x, y, z, t = F C + x, y, z, t − C − .  www.nature.com/scientificreports/ The objective of the next part is to solve Eq. (7). (7)  Finding a solution for φ(x, y, z, t). The solution of Eq. (7),φ x, y, z, t , represents the electric potential for the IPMC. It is obvious that the electric potential difference between the surface electrodes at the clamped end of the IPMC is similar to the applied input voltage, V I (t) . Thus, the governing boundary conditions for this boundary problem are expressed as follows.

Lemma 1 Equation
where the function a L (V I (t)) is called distributed surface attenuation. Ideally, if IPMC electrodes are assumed perfectly conducting surface and the fractal penetration of the electrodes into the membrane is ignored, the voltage of the clamped region will be constant along the length of the electrodes without any attenuation in its longitudinal direction 23 . Since IPMC electrodes are not perfectly conducting surfaces due to the electrode deposition process and also fractal penetration of the electrodes into the membrane, the applied voltage is attenuated along the z-axis ( Fig. 2) proportional to the amplitude, frequency and other components of the input voltage.a L (V I (t)) is a continuous function that produces a real number in the (0, 1] interval and a L (V I (t)) = 1 means electrodes are perfectly conducting surfaces.
That U 1 x, y, z, t and U 2 x, y, z are defined as Eq. (12) and Eq. (13), respectively. k 1 and are constant coefficients and α 1 (t),α 2 and α 3 (t) are defined as follows: www.nature.com/scientificreports/ In the above equations, θ is 1 W and µ is a constant coefficient that should be bigger than 1. Also α 1 (t),α 2 and α 3 (t) have to always satisfy the following condition.
where and ϒ 1 to ϒ 8 and γ 1 to γ 3 are constant coefficients. The electro-mechano-chemical coupling. As mentioned in Section I, when a voltage is applied to the IPMC, hydrated sodium cations move to the cathode side in response to the applied input voltage across the membrane. This migration leads to increasing concentration of water molecules in the cathode side and gradually decreases them in the anode side. When the cathode negative charge equals the positive charge of sodium ions, this motion is stopped. Increasing water molecules in the cathode side leads to the Nafion inflation around the cathode and consequently, actuation toward the anode is observed 23,33 . Hence, to create an electro-mechanochemical coupling, the first step is to calculate the concentration of cations (electrochemical component) and then the next step is to find a charge-stress coupling (electromechanical component). Nemat-Naser and Li 31 assumed that the charge density at the surface of the polymer (ρ) is proportional to the induced stress (σ) in the form of σ(x,t) = αρ(x,t), where α is the charge-stress coupling constant. This assumption was utilized in the most papers 16,29,30,32,33 in the field of IPMC modeling. Therefore, one obtains the following equation at the surface of the polymer (x = ± h): However, we consider a three-dimensional time-variant form of this assumption, which will be shown to produce more accurate predictions in experiments, especially for large bending deformations.
As shown in Fig. 3, based on the governing principles in mechanics for a cantilever beam (i.e. an IPMC), bending in the direction of the x-axis (the main direction of the IPMC bending) needs a component of an induced www.nature.com/scientificreports/ stress in the direction of the z-axis. On the other hand, the most of papers which presented a physical model for IPMC actuators assumed that cations move in the direction of the x-axis between the anode and cathode. It means they assumed that the concentration of cations is a one-dimensional variable 16,19,[35][36][37][38][39][40][41] . Nevertheless, in this paper, we modified this assumption and extended it to a three-dimensional time-variant problem. Hence, it is assumed that ions move in the direction of vector A. Components of vector A depend on the electric resistance of the IPMC, which is defined as follows: In (23), � x (r M , t) , � y (r ew , t) and � z (r el , t) are functions of the membrane resistance (r M ), the transverse resistance (r ew ) and the longitudinal resistance (r el ) of electrodes, respectively.
Indeed, migration and concentration changes of cations, which lead to mechanical stress and bending of the IPMC, are done in the direction of vector A. Mathematically, to calculate the changes of a multivariable function in the direction of a vector, the directional derivative of the function should be calculated in the direction of the desired vector. For example, to find the changes of ion's concentration in the direction of the movement's vector, the directional derivative of C + x, y, z, t in the direction of vector A should be calculated. It is shown by ∇ A C + x, y, z, t and defines as follows: a is the unit vector of A and defined by a = A |A| . Also,∇C + is a gradient of C + x, y, z, t . Now, if it is considered that the bending displacement of a stimulated IPMC is attenuated over the time due to the resistance, hence it is sensible to choose time-ascending functions for modeling effects of the resistance (functions like a ramp, sigmoid, tangent hyperbolic, etc.). In the following, for Eqs. (25)(26)(27), we choose a ramp function for � x (r M , t) , � y (r ew , t) and � z (r el , t) and we define the mentioned resistances as functions of dimensions (h, W, and L).
where u(t) is the Heaviside step function. Consequently, ∇ A C + is obtained as follows: and r(h, W, L) is: ∇ A C + determines changes in the cations concentration in the direction of the optimum movement's vector. These changes induce the mechanical stress. Hence, there is a linear relationship between ∇ A C + and the mechanical stress: That σ 0 is a constant coupling coefficient that should be estimated in the model estimation stage.
(30) σ x, y, z, t = σ 0 ∇ A C + x, y, z, t . www.nature.com/scientificreports/ The mechanical stress. To find the induced stress using Eq. (30) we should calculate ∇ A C + and for calculation of ∇ A C + we need to know C + . It is straightforward to drive Eq. (31) for C + by combining Eqs. (1-4).
Then, from Eq. (28), ∇ A C + can be determined and finally σ is obtained as: where U 1 and U 2 already were introduced in Eqs. (12) and (13).
The bending moment. The global bending moment of a cantilevered IPMC beam, where the stimulus is applied at the clamped end of the beam (z = 0), is described by a vector like M as follows: ans r is the position vector and is defined by r = x i + y j . Equation (37) can be expressed in terms of: Figure 4 shows our geometric definitions of an IPMC beam in the cantilevered configuration. As illustrated in Fig. 4, M zx is exactly zero, M zz is almost zero and only M zy dominates the IPMC bending in the direction of the x-axis. Indeed, the only effective bending moment, M zy , can be expressed as: σ zz x, y, z, t is the renamed version of σ x, y, z, t defined in Eq. (32). Then, by substituting Eqs. (32)(33)(34)(35)(36) into Eq. (42) and solving integration terms, M zy is determined as Eq. (43).
A i x, y, z ,  www.nature.com/scientificreports/ Here, "csch(.)" is the hyperbolic cosecant operation and b mnp is a constant coefficient that is related to the numbers of harmonics in calculating U 2 and Υ 2 , also,Υ 5 and Υ 6 are constant coefficients that are produced through the combination of ϒ 1 to ϒ 8 and γ 1 to γ 3 and we should predict them in the model prediction stage. Consequently, M zy can be rewrite as Eq. (46).    16,30,32,33,42 . Based on LEB cantilever beam theory, the relation between the shape function of the IPMC,ω(z, t) , and its bending moment, M zy (z, t) , is defined as follows: Y is Young's modulus of the Nafion membrane, and I is the moment of inertia and calculated by the following equation: However, LEB cantilever beam theory is valid just for infinitesimal strains which means this theory can only model small deformation (vibration) of a beam. Hence, we cannot use it for modeling large deformation situations of an IPMC beam.
Mathematically the curvature of a function like ω(z, t) can be calculated by κ Math (z, t) which is defined as follows: On the other hand, the curvature of a beam under an applied bending moment,M zy (z, t), is mechanically defined by: Consequently, our aim is to find a function like ω(z, t) that its curvature matches with the curvature of a cantilever beam (IPMC) with an applied bending moment M zy (z, t) to it. To find this function, κ Math (z, t) is set equal to κ Mech (z, t) and the resulted PDE is solved. Then, the following nonlinear PDE is obtained. Eq. (54), if we ignore the square of the beam slope (i.e. term ∂ω(z,t) ∂z 2 ), we reach to the PDE resulted by LEB beam theory (Eq. 50) for small deformations. Equation (54) is a nonlinear PDE but it is solvable analytically. So, ω(z, t) will be obtained as follows:

Parameter estimation and model validation.
In this section, we aim to validate our claims and find an accurate well-defined relationship between the input applied voltage ( V I (t) ) and the output tip displacement of the IPMC ( δ Tip (t) ). Totally, to produce V I (t) and acquire δ Tip (t) , we used a set of hardware apparatus. As shown in this Fig. 5, the setup was composed of several components: a computer with MATLAB 2012b, a data acquisition board, a differential electronic amplifier, an accurate mechanical camera handler, a tin-plated copper clips, and a camera. To record the IPMC actuation data, the desired input signals ( V I (t) ) were produced in MATLAB and written serial to send to the data acquisition board (DAQ board). The output voltage of one of the digital to analog converter (DAC) outputs of the DAQ was amplified and shifted by a designed differential amplifier and then applied to the IPMC using a clip. When the signal was applied to the IPMC, a video was taken of the IPMC's actuation by a video camera, simultaneously. Then, by employing appropriate video and image processing techniques (we refer the reader to Appendix D), the required features were extracted from these video files 9 . The sampling rate of the data acquisition was 29.7 S/s which was equal to the frame per second of the recorded video (29.7 FPS). Table 1 summarizes specifications of the IPMC used in this research. In all physics-based models there are some specified physical parameters and some unknown parameters. The specified physical parameters were reported in the references and are listed in Table 2 but unknown parameters should be estimated in the parameter estimation stage.
For parameter estimation and validation of the model we used the GA toolbox of MATLAB software. Three different signals were applied as the input voltage to the IPMC and measured its tip displacement in response to each stimulus. Applied stimulus signals were chirp, Pseudo Random Binary Sequence (PRBS) and sinusoidal signals with a peak voltage of 2.3 V. Then, we estimated the unknown parameters of the proposed model using the (60) www.nature.com/scientificreports/ half of each input-output data pairs (voltage as the input and the tip displacement of the IPMC in the x-direction as the output). For model validation the defined model was tested by the rest of the data pairs. There are a variety of classic and intelligent methods to estimate unknown parameters and we chose a genetic algorithm-based (GA) optimization method to find unknown parameters. The estimated parameters are reported in Table 3. Figures 6 and 7 illustrate the experimental tests of the model for Chirp and PRBS applied voltages. Moreover, the displacement vs. voltage state-trajectory for chirp dataset is shown in Fig. 8. Results indicate that the proposed model follows the actual output precisely. Moreover, to quantify the accuracy of the model by a numerical criterion, the normalized mean-square error (NMSE) was selected. The NMSE of the model in response to sinusoidal, PRBS, and chirp input signals are 0.07, 0.025 and 0.0047, respectively. Accordingly, the proposed model is accurate and can work properly in large deformation situations and practical applications. Furthermore, we compared the result of the proposed model with two other methods, Laguerre expansion method (LEM) and ANFIS-NARX paradigm 10 . The reason for choosing these two methods is that we used only the information of the input to predict the output in the proposed model which means the proposed method is    www.nature.com/scientificreports/ a non-autoregressive model, hence, we should compare our model with other non-autoregressive methods. One of the best non-autoregressive models is Laguerre Expansion Technique (LET) 43 that in Table 4 we compared the NMSE of the proposed model with LET. The mean NMSE of the proposed model is about 24 times smaller than the mean NMSE of LEM which means that, the proposed model is 24 times more accurate than LEM. As the second method, we compared the proposed model with ANFIS-NARX paradigm. ANFIS-NARX paradigm is an autoregressive method that is valid for IPMC large deformation situation 10 . As reported in Table 4,  www.nature.com/scientificreports/ ANFIS-NARX paradigm is more accurate than the proposed model, but it doesn't mean it is more appropriate because ANFIS-NARX is an autoregressive model which uses output feedback to achieve a precise prediction model. Using output data for an accurate modeling eliminates the causal relationship between the applied voltage and the IPMC bending. Moreover, in most of practical applications acquiring output feedback of the IPMC is not feasible. Hence, from the viewpoint of practical applications, a non-autoregressive model even with less accuracy is preferred to a more accurate autoregressive model.

Conclusion
It this paper a fully analytical and physics-based 3D ion transport model was presented for large deformable IPMC actuators. This model exhibits a compact and explicit relationship between the input voltage and the output tip displacement of the IPMC actuator which is uniquely valid for large deformation situation. During the various stages of modeling, four Lemmas were presented which support that the proposed model passed a true way as well as the previous models for small deformations and moreover, is an explanation of them for large deformations. Furthermore, several experimental results were presented to demonstrate its applicability to arbitrary electrical inputs. The agreement between model predictions and experimental results also provides insight into the underlying actuating mechanisms of IPMC materials due to the physical based of the model.