Multiphysics Simulator for the IPMC Actuator: Mathematical Model, Finite Difference Scheme, Fast Numerical Algorithm, and Verification

The article is devoted to the development and creation of a multiphysics simulator that can, on the one hand, simulate the most significant physical processes in the IPMC actuator, and on the other hand, unlike commercial products such as COMSOL, can use computing resources economically. The developed mathematical model is an adjoint differential equation describing the transport of charged particles and water molecules in the ion-exchange membrane, the electrostatic field inside, and the mechanical deformation of the actuator. The distribution of the electrostatic potential in the interelectrode space is located by means of the solution of the Poisson equation with the Dirichlet boundary conditions, where the charge density is a function of the concentration of cations inside the membrane. The cation distribution was obtained by means of the solution of the equation system, in which the fluxes of ions and water molecules are described by the modified Nernst-Planck equations with boundary conditions of the third kind (the Robin problem). The cantilever beam forced oscillation equation in the presence of resistance (allowing for dissipative processes) with assumptions of elasticity theory was used to describe the actuator motion. A combination of the following computational methods was used as a numerical algorithm for the solution: the Poisson equation was solved by a direct method, the modified Nernst-Planck equations were solved by the Newton-Raphson method, and the mechanical oscillation equation was solved using an explicit scheme. For this model, a difference scheme has been created and an algorithm has been described, which can be implemented in any programming language and allows for fast computational experiments. On the basis of the created algorithm and with the help of the obtained experimental data, a program has been created and the verification of the difference scheme and the algorithm has been performed. Model parameters have been determined, and recommendations on the ranges of applicability of the algorithm and the program have been given.


Introduction
The development of propulsors for industrial and medical microrobots is an urgent task. Due to the requirements for such propulsors (workability, energy efficiency, the possibility of creating large forces and displacements in space), ionic polymer-metal composites (IPMCs) attract special attention. IPMCs are a sandwich structure consisting of a porous ion-exchange membrane impregnated with electrolyte and coated on both sides with metal electrodes [1,2].
The main inherent properties of IPMCs are flexibility, light weight, ease of manufacture, and processing. In addition, IPMCs are capable of large deformation in response to a voltage applied to the metallization of several volts. Due to this, they can be used as soft robotic actuators, artificial muscles, and dynamic sensors in the field of bionic engineering [3].
The operating principle of an IPMC actuator is based on transport processes in a polymer membrane. During water saturation, the dry polymer is structured in such a way that the hydrophilic ends of the polymer chains are oriented towards the water-filled membrane pores. Under the action of an electric field caused by the voltage applied to the electrodes, the charged liquid component in the membrane moves along the through pore system. The resulting electroosmotic water flow causes an increase in the liquid pressure at one electrode and a decrease in the vicinity of the other electrode. The pressure difference leads to bending of the IPMC actuator [1,4].
There are numerous attempts to describe the actuation mechanism of IPMCs. Since the middle of the last century, many researchers have been trying to simulate actuators and sensors based on IPMC.
To study the behavior and predict the motion of an IPMC actuator, it is necessary to have reliable models. Currently, in order to describe the physical processes in the IPMC actuator, a large number of models have been developed, which are conventionally divided into three groups: multiphysics models (white box models), which describe the behavior of the problem parameter fields-such as ion concentration, potential, and deformation [5,6]-models based on the method of analogies or models with lumped parameters (grey box models) [7,8], and empirical models (black box models), which are approximation curves of experimental data [9]. In [6], a continuum electromechanical model is proposed to improve the production process and the composition of the IPMC. In [10], the authors presented the IPMC electromechanical reaction model, based on the electrostatic forces of attraction/repulsion in the IPMC, to predict the behavior of actuators, in particular, to explain the relaxation of the IPMC. These models can predict the bending behavior of the IPMC with a sufficient degree of accuracy, but they require many predetermined physicochemical properties that are measured in experiments. Black box models, also called empirical and phenomenological models-presented in [11]-are applicable only to specific shapes and boundary conditions. A more reasonable black box model was introduced in [12]. The new model was validated using the experimental results from some materials [13]. Grey box models proposed at [7,8] are based on physical laws and parameters determined by experiments.
Despite the long research time and a large number of different approaches, none of the models has proven its applicability to IPMCs with different thicknesses and different compositions of the ionic polymer membrane [14].
Further development led to the creation of a model that explicitly took into account the solvent (water) dynamics [15]. This model has a significant advantage of being able to describe the dynamic response that is highly dependent on the IPMC moisture saturation. However, it is difficult to perform analysis and simulation [16], since the mathematical model presented in [15] consists of complex nonlinear partial differential equations. We have made an attempt to solve this problem by creating a calculation scheme and a program code. A combination of the following computational methods was used as a numerical algorithm for the solution: the Poisson equation was solved by a direct method, the modified Nernst-Planck equations were solved by the Newton-Raphson method, and the mechanical oscillation equation was solved using an explicit scheme.

Mathematical Model
Taking into account the complex nature of non-stationary processes in the considered IPMC actuator, the developed mathematical model of the actuator is based on the thermodynamic theory of irreversible processes and includes the modified Nernst-Planck equations [15], the Poisson equation [17] and the mechanical oscillation equation of a cantilever beam with a fixed end [18,19].
The modified Nernst-Planck equations considered in [15] describe the transfer processes of ions and water molecules in the IPMC membrane with the following assumptions:

1.
Ionic polymer-metal composite is considered to be two-phase and includes the solid phase, which is a polymer porous structure, fixed negative charge and metal electrodes, and the liquid phase, which includes cations and water molecules, redistributed under an electric field and/or a mechanical load.

2.
The liquid phase flux consists of two components: diffusion (including electromigration) and convective. The diffusion fluxes of ions and water molecules are determined by the potential gradient, the concentration gradients of ions and water molecules, and the hydrostatic pressure gradient created by redistribution of ions and water molecules in polymer nanopores. The solid phase influences the diffusion fluxes through the nanopore structure and the electric field of fixed negative ions in the membrane. The convective fluxes are determined by the elastic force of the solid phase. 3.
In a short time interval, the hydraulic pressure and the inherent mechanical stress are balanced with the elastic stress of the composite solid phase.
As applied to the problem under consideration, the modified Nernst-Planck equations can be written as [15] where J I is the flux density of cations in the polymer volume; J W is the flux density of water molecules in the polymer volume; C I is the concentration of cations in the polymer volume; C W is the concentration of water molecules in the polymer volume; p is the hydraulic pressure in the polymer; V I is the molar volume of ions; V W is the molar volume of water; n dW is the number of water molecules associated with one cation; D II is the self-diffusion coefficient of cations; D WW is the self-diffusion coefficient of water; K is the filtration coefficient (Darcy's law); ϕ is the electrostatic potential; Z I is the relative charge of ion; T is the absolute temperature; F is the Faraday constant; R is the gas constant; y is the coordinate.
The summands −D II and −D II n dW ∂y reflect the interaction of the diffusion fluxes of cations and water molecules in the IPMC. Since a cation carries n dW water molecules by hydration, the diffusion coefficients of the interaction between the fluxes of cations and water D WI and, accordingly, D IW are expressed in these summands through the self-diffusion coefficient of ions as D WI = D II n dW ; The summands −C I K ∂p ∂y and −C W K ∂p ∂y describe, using Darcy's law [15], the convective components of the flux density, which determine the relationship between the hydraulic pressure of the liquid phase and the elastic mechanical stress of the solid phase of IPMC.
In order to study the dynamics of changes in the spatial distributions of the concentrations of cations and water molecules in time, taking into account that the rates of changes in the concentrations of ions and water molecules in time are determined by the divergence of the corresponding flux densities where t is time, the Nernst-Planck Equations (1) and (2) can be written as The concentration of cations at the initial time C I (y, t min ) (the initial condition for Equation (7)) is determined by the initial concentration of ions in the polymer C + in accordance with the expression where ρ SPN is the layer density of the dry membrane at normal ambient humidity; h is the thickness of the dry membrane (excluding metal electrodes) at normal ambient humidity; α is the expansion coefficient of the membrane at maximum humidification; P WN is the mass fraction of water in the dry polymer at normal ambient humidity. The boundary conditions for Equation (7) determine the ion flux densities through the outer boundaries of the polymer, caused by the evaporation processes where y min , y max are the coordinates of the polymer membrane boundaries; H is the thickness of metal electrodes; γ I is the dimensionless empirical coefficient determining the evaporation rate of cations into the external environment. The concentration of water molecules at the initial time C W (y, t min ) (the initial condition for Equation (8)) is determined by the degree of membrane humidification K w (the ratio of the concentration of water molecules in the polymer to the maximum possible concentration) in accordance with the expression where M W is the molar mass of water; P WS is the mass fraction of water in the maximum humidified polymer. The boundary conditions for Equation (8) determine the fluxes of water molecules evaporated through the outer boundaries of the polymer where γ W the dimensionless empirical coefficient determining the evaporation rate of water molecules into the external environment. The introduction of the empirical coefficients γ I , γ W into boundary conditions (10), (11), (13), (14) is due to the complexity of theoretical evaluation of the evaporation rates of water molecules and ions. This is due, in particular, to the considerable variation in the parameters of the granular structure of metal electrodes.
The Poisson equation connects the spatial distributions of the potential and the electric field in the IPMC membrane with the voltage at the electrodes and the spatial distribution of ions [17] where C − is the concentration of anions; q is the elementary charge; ε is the relative permittivity of water; ε 0 is the permittivity of vacuum. A uniform stationary spatial distribution of the concentration of anions is assumed The boundary conditions for the Poisson Equation (15) are determined by the potentials at the electrodes ϕ(y min , t) = 0; where U(t) is the time dependence of the voltage applied to the electrodes. The mechanical oscillation equation of a cantilever beam with a distributed mass, taking into account attenuation, describes the tip displacement of the beam relative to the equilibrium position under the action of an exciting force created by the transport of ions and water molecules and their interaction with the composite solid phase [18,19] where s(t) is the tip displacement of the beam relative to the equilibrium position; m L is the linear density of the cantilever beam; F L (t) is the exciting force per unit of beam length; ω 0 is the natural oscillation frequency of the cantilever beam; β is the coefficient characterizing dissipative processes. For the cantilever beam stationary bending, corresponding to the tip displacement s 0 of the beam relative to the equilibrium position, the right-hand side of Equation (19) can be expressed through the natural oscillation frequency ω 0 of the cantilever beam as follows where L S = L(1 + α) is the length of the humidified beam; k is the stiffness coefficient of the cantilever beam. Taking into account that the bending moment is created by the transport of ions and water molecules and their interaction with the composite solid phase is approximately uniform along the beam length, the tip displacement of the beam relative to the equilibrium position can be expressed through the bending moment as [19] where M is the bending moment; E eq is the equivalent Young's modulus of a three-layer cantilever beam (polymer membrane with electrodes); J y is the moment of inertia of the beam section. Then, substituting (20) and (21) into the right-hand side of Equation (19), we obtain The bending moment M(t) is determined by the spatial distribution of the hydraulic pressure p(y, t) in accordance with the expression where (23) determines the average value of the hydraulic pressure in the beam at the current time.
The spatial distribution of the hydraulic pressure at the current time can be expressed through the spatial distributions of the concentrations of water molecules and ions using the linear relation [1,20] as where t min is the initial time; η I , η W is the empirical coefficients having pressure units. The product of the equivalent Young's modulus of the cantilever beam and the moment of inertia of the section for the problem under consideration is expressed by the integral In expression (25), the dependence of the Young's modulus on the coordinate E(y, t) at the current time is considered not only in the polymer, but also in metal electrodes, as evidenced by the integration limits ± h S 2 + H . In this case, the dependence of the Young's modulus on the concentration of water molecules in the polymer is taken into account, which is determined experimentally in accordance with the procedure described in [21].
After substituting (23), (24), and (25) into (22), we obtain the cantilever beam mechanical oscillation equation in the form The initial conditions for Equation (26) determine zero values of the tip displacement and the tip velocity of the beam at the initial time s(t min ) = 0; The resonant frequency ω 0 of mechanical oscillations of the beam with a distributed mass in Equation (26) is determined as [19] where λ is the coefficient depending on the mode of beam bending oscillations. The linear density of the beam m L (t) is determined as the sum of the linear densities of the polymer, electrodes and water in the polymer and, taking into account the water evaporation process, is a function of time where m D is the mass of the beam dry polymer; ρ M is the density of the electrode material. Substituting expressions (25) and (30) into (29), we obtain In this study, the system of Equations (7), (8), (15), and (26) with initial conditions (9), (12), (27), (28); boundary conditions (10), (11), (13), (14), (17) and (18); and parameters (16) and (31) was solved numerically using the finite difference method without additional simplifications.

Discretization of the Mathematical Model, Numerical Simulation Technique
The discretization of the model described above was carried out within the framework of the finite difference method on uniform time G t and coordinate G y grids G t = t m = (m − 1)∆t|m = 1, . . . , M ; (32) where ∆t is the grid step in time; ∆y is the coordinate grid step; m is the point index t m of the time grid; j is the point index y j of the coordinate grid; M is the number of points in the time grid; J is the number of points in the coordinate grid. Moreover, the coordinate grid (33) covers only the polymer part of the three-layer beam between the boundaries with metal electrodes. In order to optimize the simulation time on a personal computer with insignificant RAM resources (8 GB), a self-consistent numerical solution of three subsystems of the proposed model at each time slice was carried out sequentially using a combination of the following methods: • Subsystem 1, which includes Poisson Equation (15) with boundary conditions (17), (18), was solved by a direct method • Subsystem 2, which includes modified Nernst-Planck Equations (7) and (8) with initial conditions (9) and (12) and boundary conditions (10), (11), (13), and (14), was solved using the Newton-Raphson method • Subsystem 3, which includes cantilever beam mechanical oscillation Equation (26) with initial conditions (27) and (28), was solved using an explicit scheme Taking into account the listed methods for the numerical solution of the proposed model equations, discretization schemes for subsystems 1-3 on time and coordinate grids (32) and (33) are obtained in the following form: • Subsystem 1 (15), (17), (18) where ϕ m j is the grid function of the potential; C I m j is the grid function of the concentration of cations; U m is the voltage applied to the electrodes at the time point t m .
• Subsystem 3 (26)-(28) was discretized on the time grid (32) and the extended non-uniform coordinate grid G yH including the coordinate grid G y (33), supplemented by the first and last nodes, the coordinates of which correspond to the outer boundaries of metal electrodes where s m is the grid function of the tip displacement of the cantilever beam; E m j is the grid function of the Young's modulus; ω 0m is the grid function of the resonant frequency; m Lm is the grid function of the linear density of the beam.  Due to the nonlinearity of the system of equations (34)-(50), the problem is solved iteratively. The variable k in the block diagram in Figure 1 reflects the iteration index, σ (k) is the residual at k-th iteration determined by the values of the concentration vector at k-th and (k + 1)-th iterations This technique is implemented as a specialized software created in the MATLAB programming environment.

Model Verification, Results, and Discussion
The proposed technique and software tools for the numerical implementation of model (34)-(50) allow for a detailed analysis of transients in the IPMC actuator polymer for an arbitrary mode of the control voltage oscillations on time, calculation the amplitude-frequency characteristics and dependences of the oscillation amplitude of the cantilever beam on the applied voltage amplitude at different degree of the polymer membrane humidification, taking into account the dependence of the dimensions and mass of the membrane on its humidification, and also taking into account the evaporation rates of ions and water molecules from the surface of metal electrodes.
In order to verify the considered model and the proposed technique for the numerical simulation of the IPMC actuator, a number of experimental studies have been performed.

Experimental Setup
The bench block diagram for the IPMC actuator investigation is presented in Figure 2. The investigated IPMC actuator was fixed with probes, through which the voltage was supplied from an Agilent 33500B Series waveform generator. The tip displacements of the IPMC actuator were recorded by an L-GAGE LG5A65PUQ laser displacement controller, from which information was transmitted to an Agilent DSO-X 3014A oscilloscope and then to a PC.

Numerical Simulation Results and Their Discussion-Comparison with Experimental Data
In order to verify the developed model (34)-(50), technique and software tools, the numerical simulation results of the IPMC actuator in a fairly wide range of amplitudes (peak-to-peak value U A = 0-5 V) and frequencies (f = 0.5-50 Hz) of the control voltage were obtained. The values of the physical constants and parameters of the IPMC actuator model used in the calculations are given in Table 1. The numerical simulation results of transients in the polymer membrane of the IPMC actuator, obtained using the developed software tools for implementation the model (34)-(50), are presented in Figures 3-20. These results were obtained on a coordinate grid containing 300 steps at a time grid step ∆t = 3 ms for a peak-to-peak control voltage U A = 5 V, harmonically changing in time with a frequency f = 1 Hz (Figures 3-11) and f = 10 Hz (Figures 12-20).                       The analysis of the developed model (34)-(50) and the presented graphs show that the time behavior in the spatial distributions of the concentrations of ions and water molecules in the polymer membrane volume of the IPMC actuator depends in a complicated way on a sufficiently large number of parameters given in Table 1, but it is most determined by the ratio between the diffusion coefficients of ions and water molecules in the polymer and the control voltage frequency. It can be seen from Figures 3-11 that at low control voltage frequencies f ≤ 1 Hz, wavelike changes in the concentrations have a significant amplitude throughout the entire volume of the polymer membrane, determining the noticeable on Figures 21 and 22 deviations in the form of changes in the acting force F(t) and, accordingly, the form of the beam tip mechanical oscillations s(t) from the harmonic form.
According to Figures 12-20, when the control voltage frequency increases, the amplitude of changes in the concentrations of ions and water molecules in the polymer membrane volume decreases markedly, which is mainly due to the inertia of diffusion processes. Moreover, a particularly strong decrease in the amplitude of concentration changes is observed in the central area of the membrane, where the spatial distributions of the concentrations during the entire transient remain almost uniform and even at a control voltage frequency f = 10 Hz change relative to equilibrium values by no more than 0.7% for water molecules and not more than 11% for cations (Figures 12-20). Figures 23-25, the change of the acting force in time F(t) with an increase in the control voltage frequency is more and more harmonic, and significant harmonic distortions of the form of the beam tip mechanical oscillations s(t) (see Figure 23b) reflect the interaction of the liquid and solid phases of the humidified actuator, taking into account the parameters that determine their inertia (diffusion coefficients D II , D WW for the liquid phase; elastic moduli E S , E M , densities and geometric dimensions of the beam layers for the solid phase; filtration coefficient K).

As a result, in accordance with
The result of the above-mentioned interaction of the liquid and solid phases of the humidified actuator is also a change in the phase shift between the control voltage oscillations and the beam tip mechanical oscillations, observed when increasing frequency in Figures 23-25. If at a frequency f = 10 Hz and at lower frequencies the phase shift between the control voltage and the mechanical oscillations is approximately π/2 (Figures 21-23), then at a frequency f = 36 Hz, which is close to the resonant frequency, the phase shift, in accordance with Figure 24, increases to π, and with a further increase in frequency to f = 40 Hz, in accordance with Figure 25, it reaches a value of 3π/5. It is important to note that the phase shift between the oscillations of the control voltage U(t) and the acting force F(t) remains unchanged and equal to approximately π/2 at any frequencies (Figures 21-25).
To obtain experimental dependences of the displacement amplitude on the peak-to-peak control voltage, IPMC actuators with dimensions of 20 × 5 mm based on the Nafion 117 ion-exchange membrane with a thickness of 175 µm were used. The deposition of Pt electrodes on the membrane surface was carried out according to the technology described in [22,23]. The manufactured IPMC actuators were kept in deionized water for a day, after which they were investigated on the bench described in Section 4.1. Figures 26-28 show the calculated and experimental dependences of the beam tip displacement amplitude on the peak-to-peak control voltage with a frequency of 1 Hz (Figure 26), the beam tip displacement from the DC control voltage (Figure 27), as well as the calculated and experimental amplitude-frequency characteristics (AFC) of the IPMC actuator ( Figure 28).   The analysis of the graphs presented in Figures 26-28 indicates a good agreement between the calculated and experimental dependences of the beam tip displacement amplitude on the peak-to-peak control voltage ( Figure 26) and the calculated and experimental AFC of the investigated IPMC actuator both in the beam mechanical oscillation amplitude and in the resonant frequency (f R ≈ 36 Hz), which corresponds to the characteristic maximum in Figure 28. In this case, the calculated dependence of the beam tip displacement on the DC control voltage level ( Figure 27) gives a good agreement with the experiment only at control voltage U A ≤ 1.5 V, which is probably due to the influence of factors not considered in the model (34)-(50).

Conclusions
This article focuses on the development and numerical implementation of a mathematical model of an actuator based on ionic polymer-metal composite (IPMC). To simulate the most significant physical processes in the IPMC actuator such as transients of the spatial distributions of the ion and water molecule concentrations under the changing in time the electric field and the process of cantilever beam mechanical oscillations using comparatively small computing resources, multiphysics model, corresponding technique, and software tools for numerical simulation were developed. The developed IPMC actuator model is based on the thermodynamic theory of irreversible processes and includes the modified Nernst-Planck equations for both ions and water molecules, the Poisson equation, and the mechanical oscillation equation of a cantilever beam with a fixed end.
Using the developed numerical simulation software, the spatial distributions of the concentrations of ions and water molecules in the IPMC polymer membrane, transients of the force and the beam tip displacement at different values of amplitude and frequency of the control voltage, as well as amplitude-frequency characteristics of the actuator were calculated and investigated.
For verification of the developed numerical model and software, experimental characteristics of the IPMC actuator were measured. According to the comparative analysis, the numerical simulation results demonstrate good agreement with the experimental data.
The proposed numerical model of the IPMC actuator can be implemented in any programming language and allows for fast computational experiments.