%% LyX 1.4.4 created this file.  For more info, see http://www.lyx.org/.
%% Do not edit unless you really know what you are doing.
\documentclass[twocolumn,english]{revtex4}
\usepackage[T1]{fontenc}
\usepackage[latin1]{inputenc}
\usepackage{longtable}
\usepackage{amsmath}
\usepackage{graphicx}

\makeatletter

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% LyX specific LaTeX commands.
%% Because html converters don't know tabularnewline
\providecommand{\tabularnewline}{\\}

\usepackage{babel}
\makeatother
\begin{document}

\title{Finite Element Model of a Self-Oscillating IPMC}


\author{Deivid Pugal}


\affiliation{Active Materials and Processing Laboratory, Mechanical Engineering
Department, University of Nevada, Reno, U.S.A}


\affiliation{IMS Lab, Institute of Technology, Tartu University, Estonia}


\author{Kwang J. Kim}


\email{kwangkim@unr.edu}


\affiliation{Active Materials and Processing Laboratory, Mechanical Engineering
Department, University of Nevada, Reno, U.S.A}


\author{Andres Punning, Heiki Kasemägi and Alvo Aabloo}


\affiliation{IMS Lab, Institute of Technology, Tartu University, Estonia}

\begin{abstract}
This paper presents a electro-mechanical model of an ionic polymer-metal
composite (IPMC) material. The technique of the modelling is Finite
Element method (FEM). An applied electric field causes the drift of
counter-ions (e.g. Na+) which in turn drags water molecules along.
The mass and charge imbalance inside the polymer is the main cause
of the bending motion of the IPMC. All foregoing physical effects
have been considered as time dependent and modelled with FEM. The
model takes into account the mechanical properties of Nafion polymer
as well as thin coating of the platinum electrodes. (Furthermore,
the water drift and the counter ion drift inside the polymer backbone
is considered separately.) As self-oscillating behavior of an IPMC
was described in our previous paper, the modelling of the electrochemical
reactions is also considered. Reactions occuring on the surface of
the platinum electrode, which is immersed into formaldehyde (HCHO)
solution, are described by partial differential equations and also
modeled using FEM. By coupling the equations to the rest of the model,
we are able to simulate self-oscillatory behavior of an IPMC sheet.
\end{abstract}
\maketitle

\section{Introduction}

Electroactive polymer (EAP) based materials are valuable for many
applications from micro robotics to military and space applications.
Some of the advantages of EAP materials are light weight, noiseless
actuation, simple mechanics and large displacement. In addition some
EAPs, such as IPMCs, are able to function in aquous environments.
Those qualities make the materials possible to use as artificial muscles.
In this paper we consider simulations of IPMC type materials with
the Finite Element method.

IPMC materials are highly porous polymer materials such as Nafion$^{TM}$,
Teflon, filled with some kind of ionic conductive liquid. There are
water based IPMCs which operate in aquatic environment and current
is caused by ions such as Na+, K+ dissociated in water. Other kind
of IPMCs - ionic liquid based, do not need wet environment for operating.
A sheet of an ionic polymer is coated with a thin metal layer, usually
platinum or gold. All freely movable cations inside the polymer migrate
towards an electrode due to applied electric field, causing expansion
of the material at the one end of the sheet and contraction at the
other end, which results in bending of the sheet.

For simulating actuation of an IPMC sheet we need to solve coupled
problems due to the complex nature of bending of an IPMC. It includes
simulations in different domains such as mechanical, electrostatic
and mass transfer, and even electrochemical for more advanced models.
Some authors like Wallmersperger and Nemat-Nasser have already simulated
mass transfer and electrostatic effects. We have used similar approach
in our model. Toi{[}toi2005fea] has shown a Finite Element model including
viscocity terms in transportation processes. However, the described
model is good for a rectangular beam with 2 pairs of electrodes. 
Our approach for simulating mechanical bending is taking advantage
of the numerical nature of FEM problems - we use continuum mechanics
equations instead of analytical Euler beam theroy as used by many
authors{[}]. By coupling equations from different domains, we get
a model for an IPMC muscle sheet, which allows us to use it as a starting
point for solving more complex problems, thus we have introduced a
simulation of electrochemical reactions on the platinum electrode
of an IPMC sheet which leads to self-oscillating actuation. Spontaneous
oscillations are common phenomena in nature and have been studied
in many experiments, including electrochemical systems such as oxidation
of organic materials and metals{[}miller2006oid]. Under certain conditions
such system can generate oscillations{[}krischer2001sfs]. We have
conducted series of tests, where IPMC sheet have been immersed into
acidic formaldehyde (HCHO) solution and exposed to a constant potential.
The measurements however show current oscillations, which in turn
result in oscillating bending of the IPMC sheet{[}]. Hence we have
introduced a FEM model in this paper to describe time dependent bending
of a self-oscillating IPMC.


\section{Bending Simulation Details}

An IPMC sheet consists of a backbone polymer and a metal coating.
We have used Nafion$^{TM}$ 117, coated with thin layer of platinum
in our experiments. Mass transfer and electrostatic simulations are
done in three mechanical domains - pure backbone polymer and mixture
of polymer and platinum - it comes from the fact that some platinum
is diffused into the polymer during coating process {[}]. Physical
properties and dimensions of pure, $2\mu m$ thick platinum coating
are considered only when calculating bending. That gives us five mechanical
domains as shown in Figure \ref{fig:Illustration-of-domins}. Most
simulations are carried out for an IPMC strip, 2 - 4 cm long, 200
$\mu m$ thick polymer, including 2 $\mu m$ thick Pt diffusion region
on each side, coated with 2 $\mu m$ thick platinum, in a cantilever
configuration - one end of the strip is fixed. Gravitational forces
are not considered in any of the simulations.

%
\begin{figure}
\includegraphics[scale=0.32]{IPMC}


\caption{\label{fig:Illustration-of-domins}Illustration of domins and dimensions
used in modeling}
\end{figure}


Nernst-Planck equation describes diffusion, convection and in presence
of electric field and charges, migration of particles. The general
form of the equation is \begin{equation}
\frac{\partial C}{\partial t}+\nabla\cdot(-D\nabla C-z\mu FC\nabla\phi)=-\vec{u}\cdot\nabla C,\label{eq:nernst planck}\end{equation}
where $C$ is concentration, $\mu$ mobility of species, $D$ diffusion
constant, $T$ absolute temperature, $R$ gas constant, $F$ Faraday
constant, $\vec{u}$ velocity, $z$ charge number and $\phi$ electric
potential. This equation is described for movable cations. As anions
are fixed, they maintain constant charge density throughout the polymer.
As voltage is applied to the electrodes of an IPMC, all free cations
start migrating towards cathode, causing current in the outer electric
circuit. As ions cannot move beyond the boundary of the polymer, charges
start to accumulate, resulting in increase of electric field, which
cancels out the applied one. It could be described by Gauss' Law:\begin{equation}
\nabla\vec{E}=-\Delta\phi=\frac{F\cdot\rho}{\varepsilon},\label{eq:gauss}\end{equation}
where $\rho$ is charge density, $\varepsilon$ is absolute dielectric
constant and $E$ is the strength of the electric field. A steady
state of the cations forms when electric field created by distribution
of cations cancels out applied electric field, i.e. the strength of
electric field inside the polymer is approximately zero as shown in
Figure \ref{fig:The-concentration-of}. The steady state cation concentration
of $1200\frac{mol}{m^{3}}$ is also brought out in the same figure.
%
\begin{figure}
\includegraphics[scale=0.32]{field_concentration}


\caption{\label{fig:The-concentration-of}The concentration of counter ions
and electric field strength inside the polymer according to simulations.}
\end{figure}
 It is interesting to notice that there are fluctuations in charge
distribution only in really thin boundary layers. There is no charge
imbalance inside the polymer. General conclusion by many authors is
that locally generated charge imbalance nearby platinum electrodes
is directly connected to the bending of an IPMC{[}nematnasser2000eri].
As it worked out very well for our model, we followed the example
of Leo{[}leo:cmi] and defined longitudinal force per unit area in
each point in the polymer of an IPMC as follows: \begin{equation}
\vec{F}=(A\cdot\rho+sgn(\rho)\cdot B\cdot\rho^{2})\cdot\hat{x},\label{eq:force}\end{equation}
where $\rho$ is charge density and $A$ and $B$ are constants which
could be found from different experiments. $sgn(\rho)$ is to preserve
correct sign of the excess charge. The shape of the force inside the
polymer is also shown in Figure . To relate the force to the physical
bending of an IPMC sheet, we use continuum mechanics equations, described
in Comsol Multiphysics structural mechanics software package. Normal
and shear strain are defined as: \begin{equation}
\varepsilon_{i}=\frac{\partial u_{i}}{\partial x_{i}},\;\varepsilon_{ij}=\frac{1}{2}\left(\frac{\partial u_{i}}{\partial x_{j}}-\frac{\partial u_{j}}{\partial x_{i}}\right),\label{eq:normal shear strain}\end{equation}
where $u$ is the displacement vector, $x$ denotes a coordinate and
indices $i$ and $j$ are from 1 to 3 and denote components correspondingly
to x, y, or z direction. The stress-strain relationship is \begin{equation}
\sigma=D\varepsilon,\label{eq:stress strain}\end{equation}
where $D$ is $6x6$ elacticity matrix, consisting of components of
Young's modulus and Poisson's ratio. The system is in equilibrium,
if the following relation is satisfied:\begin{equation}
-\nabla\cdot\sigma=\vec{F},\label{eq:equilibrium,contiinum}\end{equation}
which is basically Navier's equation for displacement. The stress
generated in the polymer is shown in Figure .

As our simulations are rather dynamic than static, we have to introduce
an equation to describe the motion of an IPMC sheet. To do that, we
use Newton's Second law \begin{equation}
\rho\frac{\partial^{2}\vec{u}}{\partial t^{2}}-\nabla\cdot c\nabla\vec{u}=\vec{F},\label{eq:newton II}\end{equation}
where the second term is the static Navier's equation. Dynamic part
is introduced by the first term of the equation. To describe damping
of the system, we use Rayleigh damping model, which is for a system
of one degree of freedom as follows: \begin{equation}
m\frac{d^{2}u}{dt^{2}}+\xi\frac{du}{dt}+ku=f(t),\label{eq:single degree of freedom}\end{equation}
 where the damping parameter $\xi$ could be expressed as $\xi=\alpha m+\beta k$.
Parameter $m$ is a mass, $k$ is a stiffness and $\alpha$ and $\beta$
are correspondingly damping coefficients. The equation for multiple
degrees of freedom is \begin{equation}
\rho\frac{\partial^{2}\vec{u}}{\partial t^{2}}-\nabla\cdot\left[c\nabla\vec{u}+c\beta\nabla\frac{\partial\vec{u}}{\partial t}\right]+\alpha\rho\frac{\partial\vec{u}}{\partial t}=\vec{F}.\label{eq:damping,multifreedom}\end{equation}
By coupling Eq. \eqref{eq:damping,multifreedom} to the previously
described model, good basic model for IPMC actuation has been got.
Figure \ref{fig:Experimental-and-simulation} shows simulation result
in comparsion to an experiment. All values used in the simulations
have been brought out in Table \ref{tab:Parameter-values-used}.

%
\begin{figure}
\includegraphics[scale=0.32]{punn_d}


\caption{\label{fig:Experimental-and-simulation}Experimental and simulation
results of tip displacement.}
\end{figure}


%
\begin{table}

\caption{\label{tab:Parameter-values-used}Parameter values used in bending
simulations.}

\begin{longtable}{ccc}
\hline
\hline 
Parameter&
Value&
Unit\tabularnewline
\hline
\endhead
\hline
\hline 
&
&
\tabularnewline
\endfoot
&
&
\tabularnewline
&
&
\tabularnewline
&
&
\tabularnewline
&
&
\tabularnewline
\end{longtable}
\end{table}



\section{Self-oscillations}

We have conducted a series of tests with IPMCs in a constant electric
field in formaldehyde (HCHO) solution. Measurements show that current
oscillations begin from applied potential of ca. $0.75V$. More about
experiments and conclusions could be read from our previous letter
. However, studies show that there are sequential electrochemical
reaction taking place on the platinum cathode. The initial burst of
the current is caused by the reaction \begin{equation}
CO_{ads}+OH_{ads}\rightarrow CO_{2}+H^{+}+e^{-}+2*,\label{eq:inital current reaction}\end{equation}
where subscript $ads$ denotes species adsorbed to the platinum and
$*$ denotes an active platinum site. The result of the reaction \eqref{eq:inital current reaction}
is freeing up 2 platinum sites, whichs causes $CO$ adsorb again.
Chronopotentiomentry scans show that before reaction \eqref{eq:inital current reaction},
the following reactions occur:\begin{eqnarray}
HCHO+* & \rightarrow & Pt-CHO_{ads}+H^{+}+e^{-},\label{eq:HCHO adsorb}\\
Pt-CHO & \rightarrow & Pt-CO_{ads}+H^{+}+e^{-},\label{eq:CO formation}\\
H_{2}O+* & \leftrightarrow & Pt-OH_{ads}+H^{+}+e^{-}.\label{eq:OH adsorb}\end{eqnarray}
HCHO is dissociated on the electrode surface at lower anodic potential.
Higher anodic potential causes dehydrogenation of water which results
in water oxidation with intermediate Pt-OH formation. We believe that
these reaction lead to oscillating potentials, which in turn lead
to self-oscillating motion of the IPMC sheet.

Series of chronopotentiometry scans were conducted to characterize
oscillations for different HCHO concentration and current densities.
As was brought out in our previous paper {[}], the oscillations start
at approximately $7\ mA/cm^{2}.$ The experiments were conducted up
to current density values of $40\ mA/cm^{2}$. Tests with HCHO concentrations
of 1 M, 2 M, and 3 M show that oscillations frequency remains constant
up to current density of $14\ mA/cm^{2}$, but after further increasing
current, in 1 M HCHO, the frequency decreases, in 2 M HCHO, the frequency
remains constant, and in 3 M HCHO solution, the frequency starts to
increase, as also shown in Fig . Our goal is to develope a model
of the frequency behavior depending on HCHO concentration and current
density. The basic model and concepts are taken from PhD thesis of
D. Kim {[}] and work of P. Strasser {[}]. To describe oscillations,
the four dynamic parameters, therefore 4 differential equations, must
be observed: concentration of adsorbed $OH$, $CO$, change of double
layer potential due to electrochemical reactions, and change of concentration
of $HCHO$ near the surface of platinum. First two variables could
be expressed for a certain current density and $HCHO$ concentration
as {[}]:\begin{eqnarray}
\dot{\theta}_{CO} & = & k_{2}M-k_{4}\theta_{CO}\theta_{OH},\label{eq:osc_CO}\\
\dot{\theta}_{OH} & = & k_{3}\theta_{CO}M-k_{-3}\theta_{OH}-k_{4}\theta_{CO}\theta_{OH},\label{eq:osc_OH}\end{eqnarray}
where $\theta_{CO}$ and $\theta_{OH}$ are normed adsorbtion coverages
of CO and OH. Variables $k_{i}$ and $M$ are described by equations:\begin{eqnarray}
k_{i}\left(\phi\right) & = & exp\left[s_{i}\left(\phi-\phi_{i}\right)\right],\label{eq:osc_k}\\
M & = & \left(1-\theta_{CO}-\theta_{OH}\right),\label{eq:osc_D}\end{eqnarray}
where $s_{i}$ is a modeling coefficient and $\phi_{i}$ are potentials
of reactions. As our model is highly dynamic, the double layer with
thickness $\delta$ near the platinum electrode is introduced. At
the far end of the layer, the concentration of formic acid is considered
constant, and due to adsorbtion of $HCHO$ on Pt, the concentration
of the solution is changing in time near the electrode. There are
two components responsible of decrease of the concentration. The first
one is direct oxidation of formic acid to $CO_{2}$ and $2H^{+}$,
the second one is adsorbtion of $CO$ on platinum surface due to electrochemical
reaction.{[}] The mechanism restoring the $HCHO$ concentration near
the surface is diffusion. So the amount of formic acid is decreasing
significantly while the adsorption rate is high and increasing due
to the diffusion during the low adsorbtion period. The equation describing
the diffusion process is\begin{equation}
\frac{\partial c_{FA}}{\partial t}=\nabla\cdot\left(D_{FA}\nabla c_{FA}\right)\label{eq:hcho_D}\end{equation}
with constant concentration at the far end of the double layer and
flux\begin{equation}
f=k_{2}M\left(1+kc_{r}\right)\cdot S_{tot},\label{eq:hcho_F}\end{equation}
as boundary condition on the electrode. Here $S_{tot}$ denotes the
total number of platinum sites per surface area, $c_{r}$ is normed
concentration near boundary layer and equals $c_{r}=c_{FA}/c_{0}$,
where $c_{0}$ is initial concentration. Variable $k$ is a simulation
constant. The second term of EQ \eqref{eq:hcho_F} represents simplified
version of the direct oxidation path, shown by Strasser. Considering
those equations and interesting nature of frequency characteristic
for different amounts of HCHO (Figure {[}]), we can now describe
somewhat empirical equation for the last dynamic variable - double
layer potential:\begin{eqnarray}
\dot{\phi} & = & \frac{1}{C_{dl}}\left[j_{th}-j_{d}+A\times B\times j^{2}\times\left(j-j_{th}\right)\times c_{r}-\right.\nonumber \\
 &  & \left.-S_{tot}F\left(k_{1}M+k_{4}\theta_{CO}\theta_{OH}\right)\right],\label{eq:osc_phi}\end{eqnarray}
where $j$ is applied current density, $j_{th}$ is threshold current
density with approximate value of $10mA/cm^{2}$, $j_{d}$ is direct
current density and is proportional to the second term of EQ. \eqref{eq:hcho_F}.
Variable $B$ could be explicitly written as $B=c_{0}-c_{neutral}$,
where $c_{netural}$ corresponds to concentration of 2 M. This is
denoted as {}``netural'' concentration, because as could be seen
in Figure , the oscillation frequency for the case $c_{0}=c_{neutral}=2M$
does not depend on the applied current density. $F$ is Faraday constant.
The numeric data could be found in Table \ref{tab:Variables-and-values}.
The third term in EQ. \eqref{eq:osc_phi} is empirical and reflects
the interesting behavior of the oscillations frequency for different
formic acid concentrations. Other terms are similar to the ones described
by Strasser{[}]. However, some values are adjusted to get realistic
simulation results. Measured voltage oscillation comparison to simulation
data could be seen in Figure \ref{fig:Potential-oscillations.}.

%
\begin{table}

\caption{\label{tab:Variables-and-values}Variables and values used in the
simulation of electrochemical oscillations.}

\begin{longtable}{ccc}
\hline
\hline 
Parameter&
Value&
Unit\tabularnewline
\hline
\endhead
\hline
\hline 
\multicolumn{2}{l}{Doyeon}&
\tabularnewline
\endfoot
$S_{tot}$&
$0.5\times10^{-6}$&
$mol/cm^{2}$\tabularnewline
$F$&
$96\times10^{6}$&
$mC/mol$\tabularnewline
$C_{dl}$&
$1$&
$mF/cm^{2}$\tabularnewline
$A$&
$1.2$&
$cm^{7}/\left(mA^{2}\times mol\right)$\tabularnewline
$k$&
$0.5\times10^{2}$ &
-\tabularnewline
$\phi_{1,2,3,-3,4}^{a}$&
$\left[0.2,\ 0.3,\ 0.01,\ 0.512,\ 0.77\right]$&
$V$\tabularnewline
$s_{1,2,3,-3,4}^{a}$&
$\left[10,\ -11,\ 9,\ -9,\ 20\right]$&
$V^{-1}$\tabularnewline
\end{longtable}

\begin{centering}

\par\end{centering}
\end{table}


By coupling described equations to the base Finite Element Model,
described in previous section, we can simulate oscillationg deflection
of an IPMC muscle. Two sample results for different HCHO concentrations
and current densities are shown in Figures \ref{fig:Oscillating-tip-displacement.1M}
and \ref{fig:Oscillating-tip-displacement.3M}.

%
\begin{figure}
\includegraphics[scale=0.32]{frequency}


\caption{\label{fig:Experimental-and-simulated frequency}Experimental and
simulated frequency dependence on concentration of HCHO and applied
current density. Simulations for 1 M HCHO concentration does not go
past $30\ mA/cm^{2}$, because given equation system did not give
reasonable results beyond that current density.}
\end{figure}


%
\begin{figure}
\includegraphics[scale=0.32]{potential_osc}


\caption{\label{fig:Potential-oscillations.}Potential oscillations. Measured
data and simulated data for 3M HCHO solution, with applied current
density of $20\ mA/cm^{2}.$}
\end{figure}


%
\begin{figure}
\includegraphics[scale=0.32]{disp_osc1m20ma}


\caption{\label{fig:Oscillating-tip-displacement.1M}Oscillating tip displacement.
Experimental data and simulation data for 1M solution, applied current
density of $20\ mA/cm^{2}.$}
\end{figure}


%
\begin{figure}
\includegraphics[scale=0.32]{disp_osc3m10ma}


\caption{\label{fig:Oscillating-tip-displacement.3M}Oscillating tip displacement.
Experimental data and simulation data for 3M solution, applied current
density of $10\ mA/cm^{2}.$}
\end{figure}

\end{document}