\documentclass[twocolumn,10pt]{asme2e}

\usepackage[T1]{fontenc}
\usepackage[latin9]{inputenc}
\usepackage{float}
\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{esint}

\bibliographystyle{asmems4}
%% The class has several options
%  onecolumn/twocolumn - format for one or two columns per page
%  10pt/11pt/12pt - use 10, 11, or 12 point font
%  oneside/twoside - format for oneside/twosided printing
%  final/draft - format for final/draft copy
%  cleanfoot - take out copyright info in footer leave page number
%  cleanhead - take out the conference banner on the title page
%  titlepage/notitlepage - put in titlepage or leave out titlepage
%
%% The default is oneside, onecolumn, 10pt, final

%% Replace here with information related to your conference
\confshortname{DSCC 2009} \conffullname{2009 ASME Dynamic Systems and Control Conference
and Bath/ASME Symposium on Fluid Power & Motion Control} \confdate{12-14}
\confmonth{October} \confyear{2009} \confcity{Hollywood,
California} \confcountry{USA}

\papernum{DSCC2009-XXXXX}


\title{Modeling electric currents in the electrodes of IPMC}

%%% first author
\author{Deivid Pugal
    \affiliation{
    Active Materials and Processing Laboratory\\
    Mechanical Engineering Department\\
    University of Nevada, Reno\\
    Reno, NV, 89557\\
    Email: pugald@unr.nevada.edu
    }
    \affiliation{
    IMS Lab\\
    Institute of Technology\\
    Tartu University\\
    Estonia\\
    }
}

%%% second author
%%% remove the following entry for single author papers
%%% add more entries for additional authors
\author{Alvo Aabloo
	\affiliation{
    IMS Lab\\
    Institute of Technology\\
    Tartu University\\
    Estonia\\
    Email: alvo.aabloo@ut.ee
    }
}

\author{Kwang J. Kim\thanks{Address all correspondence to this author.} \\
    \affiliation{Active Materials and Processing Laboratory\\
    Mechanical Engineering Department\\
    University of Nevada, Reno\\
    Reno, NV, 89557\\
    Email: kwangkim@unr.edu
    }
}

\begin{document}
%\pagestyle{empty}

\maketitle

%%% ABSTRACT %%%
\begin{abstract}
This paper presents a Finite Element Analysis (FEA) of an ionic 
polymer-metal composite (IPMC) material. The IPMC materials are 
known to bend when electric field is applied on the electrodes and
also produce a potential difference on the electrodes when bent.
Several authors have used FEA to describe that fenomenon and therefore
rather precise basic Finite Element (FE) models already exist.
However, the focus of the current study is to model the electric
currents in the electrodes. The basis of the calculations is the
Ramo-Shockley theorem, which has been used in other areas of physics
to describe the currents in a circuit due to a charge movement in a
media. We have used the theorem to calculate the current density
in the continuous electrodes of IPMC due to the ion migration
in the backbone polymer. Along the current densities we are able
to calculate exact voltage on the electrode at a given time moment.
In this stage, the electrodes are modeled as static domains, i.e.
no dynamic resistance changes have been considered. Although the
work is still in progress, we have got some physically reasonable
results to demonstrate and a section of the paper considers the
extendibility of the model. Those results have already been
previously published. However, in the paper we also consider
some possible applications of the developed model and also couple
of  concerns in regard to modeling the whole IPMC in three
dimensional domain with FE method.
\end{abstract}

\section{INTRODUCTION}

Electroactive polymer actuators have gained a lot of attention in
many fields such as robotics and micro electronics. The advantages
of EAP actuators are relatively simple mechanics and noiseless actuation.
Additionally some EAPs, such as IPMCs\cite{shahinpoor2001ipm}, are
able to function in aqueous environments. Those qualities make the
materials possible to use as so called artificial muscles. In this
paper we consider three dimensional time dependent simulations of
IPMC type materials with the Finite Element Method.

One of the most important qualities of IPMC materials is relatively
large amplitude bending in response to electrical stimulation. An
ion exchange polymer membrane, such as Nafion, is covered
with metal layers. The metal is typically platinum or gold. During
the fabrication process the polymer membrane is saturated with certain
solvent and ions. When voltage is applied to the metal electrodes,
the ions start migrating due to the applied electric field. Migrating
ions usually drag some solvent with them, causing expansion and contractions
respectively near the surface layers. That in turn causes bending
like actuation of IPMC sheet. 

To predict the actuation, a good understanding of underlying processes
must exist. Electrostatics, mass transfer and mechanical effects must
be taken account to get a minimal functional base model which could
predict actuation. Usually two dimensional time dependent model would
be enough to get reasonable results. However, in this paper we consider
three dimensional model of IPMC. This allows to take into account
surface resistance changes for whole area of the metallic layer. Some
authors\cite{wallmersperger2004cce,nematnasser2006mer,branco2006dcm,analysisofelec}
have already simulated mass transfer and electrostatic effects. We
used similar approach in our model. Toi\cite{toi2005fea} has shown
a Finite Element model including viscosity terms in transportation
processes explicitly. The simulation is performed as time dependent
and for three dimensions. However, the basis of the described model
is a rectangular beam with 2 pairs of electrodes. 

The current work is an extension of our already published paper\cite{pugalspie2009x}, where we
studied the surface currents and the voltage drops on the electrodes. 
The analytical surface resistance model is derived by Chen\cite{acontrolorient}
for 2D domain.The proposed FE simulations work both in 2D and 3D domain.
However, the current paper provides more insight of the challenges and some necessary
optimizations in regard to the simulations.


\section{THE SURFACE CURRENT MODEL}

We have used Nafion$^{TM}$ 117, coated with thin layer of platinum
in our experiments and therefore in theory. Mass transfer and electrostatic
simulations are done only for backbone polymer. Surface current calculations
are carried out for electrodes. All simulations are done for an IPMC
strip of $180\mu m$ thick polymer coated with $3\mu m$ thick platinum,
in a cantilever configuration - one end of the strip is fixed (see
Figure~\ref{fig:figure_IPMC}). %
\begin{figure}[t]
\begin{center}
\includegraphics[scale=0.23]{images/IPMC.eps}
\end{center}
\caption{The simulated IPMC strip - the latter voltage drop simulation
was conducted for an IPMC with the length of $40mm$. the image is out of
scale for illustraion purpose.
}
\label{fig:figure_IPMC}
\end{figure}

\subsection{The base model}

The cation migration in the polymer backbone is described by the Nernst-Planck
equation, which covers migration and diffusion part. 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$ universal gas constant,
$F$ Faraday constant, $\vec{u}$ velocity, $z$ charge number and
$\phi$ electric potential. The equations is solved only for cations
as anions are fixed in the polymer backbone. As voltage is applied
to the platinum electrodes, all free cations start migrating towards
cathode, causing current in the outer electric circuit. As ions cannot
move beyond the boundary of the polymer, local charge intensity starts
to increase near the surface of the platinum electrodes, resulting
in increase of electric field in the opposite direction to the applied
one. This effect could be described by by Gauss' Law:\begin{equation}
\nabla\cdot\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 and can be
also expressed as $\nabla\phi=-\vec{E}$. The formed steady state
of the cations and corresponding electric field distribution is shown
in Figure~\ref{fig:Electric-field-and}. The values of simulation
constants are shown in Table~\ref{table_sim_values}.

So far we have described the base model, which is usable in both two
dimensional and three dimensional modeling. Next we extend the model
to include the current flow in the electrodes and also to estimate
the currents and therefore the voltage distribution as well.%

%
\begin{figure}[t]
\begin{center}
\includegraphics[scale=0.45]{images/efield_conc.eps}
\end{center}
\caption{The cross section electric field and cation concentration
in IPMC. the shape of the graphs is generally same for 2D and 3D simulations.}
\label{fig:Electric-field-and}
\end{figure}


\begin{table}[t]
\caption{The simulation input values for the base model.}
\begin{center}
\label{table_sim_values}
\begin{tabular}{c c c l}
& & \\ % put some space after the caption
\hline 
Variable & Value & Dimension  & Comment\\
\hline 
$D_{cation}$ & $2\cdot10^{-11}$ & $\frac{m^{2}}{s}$ & Diffusion const.\cite{WALLMERSPERG08}\\
\hline 
$\varepsilon$ & $2.5\cdot10^{-2}$ & $\frac{F}{m}$ & Measured.\\
\hline 
$\mu$ & $8\cdot10^{-15}$ & $\frac{mol\cdot s}{kg}$ & $\mu=\frac{D}{R\cdot T}$, $T=293K$\\
\hline
\end{tabular}
\end{center}
\end{table}

\subsection{Extended model}

\begin{figure}[t]
\begin{centering}
\includegraphics[scale=0.45]{images/3dcurrents}
\caption{The base model with the surface electrode at $t=0.1s$.
The arrows indicate the cation flux and the whe white-gray area
is the voltage inside the polymer. The blue-red area is the 
voltage on the electrode. Blue equals $0$V.}
\label{fig:A-base-model}
\par\end{centering}

\end{figure}

The model described in the previous section is sufficient for most
cases and for two dimensional or three dimensional domains. However,
solving the base model in three dimensional domain does not give us
any kind of extra information - it is basically the same as an extended
two dimensional solution. At the same time there is a significant
increase in the complexity in terms of meshing and increase of degrees
of freedom of the system, which in turns leads to longer solving times.
The usefulness of the third dimension becomes eminent when the model
takes into account also the surface resistance of the electrodes and
when the surface resistance is not uniform in all three dimensions.
The surface resistance is an important characteristic of an IPMC strip.
Besides of being different for different IPMC sheets, it tends to
depend on the curvature of the IPMC strip\cite{punning2007ssi}. At
the same time, the surface resistance is the parameter which could
be rather easily altered. For instance, it is possible to make some
areas of the muscle sheet less conductive. That is the case where
three dimensional model would be more precise. However, even when
the surface resistance is considered constant, there still exists
a voltage drop along the surface and the extended model helps to take
this into account. 

The idea of the extended model is to put together part of the electrical
model, which includes active resistances of the electrodes and the
base model described in the previous section. Thus the cation transportation
and continuum mechanics is coupled to the currents/voltages in the
surface layers. Once having this model, the surface resistance could
be easily varied as well, however, it is not in the scope of the current
paper. The theory and simulation of electrical parameters of the electrodes
are considered in the following subsections.


\begin{figure}[t]
\begin{centering}
\includegraphics[scale=0.45]{images/comsol_mapped}
\caption{The mapped mesh generated with Comsol Multiphysics meshing tool.}
\label{fig:mapped_mesh}
\par\end{centering}

\end{figure}

\subsubsection{Theory}

An IPMC strip with an external power source and connections form a
closed electric circuit. There are roughly two types of conduction
mechanisms in the circuit: electron conduction in the outer part of
the circuit and the ion migration inside the IPMC. Even though the
ions move only inside the polymer backbone, there is a connection
between current in outer circuit and displacement of the ions. The
theorem which is more often used in plasma physics\cite{paris2005lic},
is called Ramo-Shockley theorem. It connects movement of charged particles
in a confined space to the currents in a connected electric circuit.
The theory, for instance, has been used to create a model of ion channels\cite{nonner2004rmc}.
The equation to describe the current in the external circuit due to
the relocation of the ions is:\begin{equation}
I=\frac{1}{V}\sum_{i}q_{i}\times\vec{W\left(r_{i}\right)}_{i}\vec{v}_{i},\label{eq:Ramo}\end{equation}
where $j$ is the index of a particle, $q$ is the charge, and $v$
is the velocity of a particle. $W$ corresponds to a electric field
which would exist without any charged particles present\cite{nonner2004rmc}.
By using Eq.~\eqref{eq:Ramo}, we can calculate the current flowing
in an electrode. This in turn is related to the voltage in the surface
layer. The surface layer could be depicted as a series of resistors
which are also connected to the polymer. For this kind of two electrode
system, the Eq.~\eqref{eq:Ramo} could be written as\begin{equation}
J=\frac{F}{d}\int_{0}^{d}\vec{f}\cdot\vec{dy},\label{eq:Ramo2}\end{equation}
where $d$ is the distance between electrodes and $\vec{f}$ is the
ion flux inside the polymer $\left[\frac{mol}{s\ m^{2}}\right]$,
$F$ is Faraday constant. The $J$ is the current density in the electrode.
By knowing the $J$, the current and also voltage in the electrode
could be calculated using numerical methods. Therefore the electrodes
are described as a conductive media by using the following equation:\begin{equation}
-\nabla\cdot\frac{\partial}{\partial t}\left(\varepsilon\nabla V\right)-\nabla\cdot\left(\sigma\nabla V\right)=0,\label{eq:currents}\end{equation}
where $V$ is voltage in the electrode, $\sigma$ is electric conductivity.
The boundary conditions of the electrodes are the applied voltage,
the ground for the contact points, and insulation for all other surfaces.

\subsubsection{Simulations}

The Comsol Multiphysics Finite Element software package is used for
simulations. That allowed to use predefined set of equations for different
domains. The simulation of the metallic electrode coupled with the
base model is rather complex, because the currents in the surface
cause the changes of voltage acting on the polymer, i.e this is a
feedback system. So building the model involved roughly the following
steps:
\begin{enumerate}
\item A domain for base model was created and the equations~\eqref{eq:nernst planck}
and~\eqref{eq:gauss} were applied - see Figure~\ref{fig:A-base-model}
\item Two surface layer domains were created. The surface layer characteristics
were set close to the values of pure platinum. The conductivity, however,
was altered a bit as the surface is not uniform for the case of IPMC.
At the same time the conductivity was held constant during the simulation
(see Figure~\ref{fig:A-base-model}). The electrode
thickness was assumed to be around $3\mu m.$ It is suggested though,
that platinum layer is $1\mu m,$ but as the model does not consider
the platinum diffusion zone, the layer was chosen to be a bit thicker
instead\cite{nematnasser2002mai}.
\item The voltage of the electrode boundaries which coincided with polymer,
were set to be equal to the electric potential $\phi$ (see Eq.~\eqref{eq:gauss}).
\end{enumerate}
Meshing in a three dimensional domain is not as straightforward as
it is for two dimensions. There are couple of things which should
be taken into account. First of all, two dimensions of an IPMC sheet
are relatively large (width and length are in range of centimeters)
but the thickness is really small, much less than a millimeter. In
addition, the thickness consists of three layers - a polymer backbone
and two layers of metal coating which are considered as separate domains.
For this kind of system, the mapped mesh (as shown in Figure~\ref{fig:mapped_mesh})
would be the most reasonable.
However, there are some restriction due to the solver (Comsol Multiphysics).
The solver does not allow value projection for the mapped mesh. Therefore
the tetrahedral mesh is used instead. The downside of that mesh is
increased solution time.

\begin{figure}[t]
\begin{centering}
\includegraphics[clip,scale=0.45]{images/voltage_surface}
\caption{The voltage drop on the electrodes at $t=0.1s$.}
\label{fig:The-voltage-drop}
\par\end{centering}

\end{figure}
%

\begin{figure}[t]
\begin{centering}
\includegraphics[clip,scale=0.45]{images/experiment_simulation}
\caption{The experimentla voltage drop at $(t=0.2s).$ The simulated voltage drop and
the experimental voltage drop show quite good correspondence specially
at the ends of the sheet. However, in the middle of the sheet, the
current example seems to be more non-linear than the simulation.}
\label{fig:Experiment}
\par\end{centering}

\end{figure}

The given model is able to calculate surface currents and voltages
of the IPMC. Figure~\ref{fig:The-voltage-drop}
the shows voltage drop along the surface quantitatively. As it could
be seen, the voltage drop for ideal surface is rather negligible.
This result is also supported by the rough measurements for some water
based, platinum coated IPMCs produced in the lab. The small voltage
drop for enhanced surface electrodes is also supported by Shahinpoor
and Kim\cite{Shahinpoo2000}.

However, for some similar materials, rather large voltage drop at
the beginning of actuation is observed\cite{punning2007sre,punningandres}.
This could be caused by the bigger resistance or better ionic conductivity,
which allows to draw more current near the electrode and therefore
causing relatively bigger voltage drop along the electrode. To simulate
the described situation, the surface resistance is altered in the
model. The following relation was used for the cathode conductivity:\begin{equation}
\sigma=\frac{\sigma_{Pt}}{10^{5}\times x},\label{eq:sigma}\end{equation}
where $x$ is the distance from the electrode and $\sigma_{Pt}=8.9\cdot10^{6}\, S/m.$
The value of $x$ is kept always positive in the simulations and it
should be noted that it stays considerably less than one. Figure~\ref{fig:Experiment}
shows the simulated voltages and the voltages measured by Punning\cite{punningandres}
at different points on the electrode. As the voltage was measured
only at three different locations on cathode, the experimental voltage
curve could somewhat differ from the one interpolated on the graph.
%


\section{DISCUSSION}

The comparison of the developed model and experiments show that taking
account the voltage drop in the electrode makes the model more precise.
That could benefit even more when the developed continuum model is
linked to a external circuit model to describe the system as one complete
circuit. The developed model could be extended even further. As it
is known, the platinum layer is not uniform on the surface of the
polymer\cite{shahinpoor2001ipm}. Therefore the resistance of the
platinum is also somewhat increasing at the bottom layers, so by suggesting
different resistance profiles, the model could be tuned to be more
precise\@. Another benefit what could have when considering the surface
electrodes, is the estimation the of heat generation on the electrodes.

\subsection{Required optimizations}

However, as it is brought out in the earlier sections of this paper,
the model is rather complex and the calculations are time consuming.
This limits the usage of the model. However, there 
are some things which could be done to optimize the model and therefore
reduce the solution time.

For the three-dimensional model, the number of meshed elements is very 
large and therefore the solver must deal with more than a hunderd 
thousand degrees of freedom. As brought out before, the current solver
does not allow using the mapped meshing due to the projection restrictions.
So the goal is to generate as efficient tetrahedral mesh as possible.
One way to achieve that is to test other preprocessors besides Comsol's
integrated one. For instance, we have already tested Gmsh\cite{geuzaine2008gtd}
meshing for a small piece of IPMC with electrodes (see Figure ~\ref{fig:Gmsh_example}).

Another possible way to reduce the calculation time is to use a scaled
model of IPMC - at least for the preliminary results. Final calculations
could be done with a full scale model, however, the reduced scale helps 
significantly save the model development time. In some extent, we
have used this for the current simulations already. To further
simplify the calculations, timesteps could be used. It means that
the static model is used for electrodes and by using a scripting
language, we could solve the model for different timesteps. This 
would result close to the surface dynamic model.

Yet not explored area is using the weak form for the electrodes. As the coupling
is complicated and nonlinear, the weak form (see Eq.~\eqref{eq:General}) could give us
way better convergence. For the general equation
\begin{equation}\nabla\cdot\Gamma=F,\label{eq:General}\end{equation} in
the domain $\Omega$, the weak form 
\begin{equation}\int_\Omega\nu\nabla \cdot\Gamma dA =\int_\Omega\nu F dA\label{eq:weak}\end{equation}
is obtained by multiplying the Eq~\eqref{eq:General} by a test function $\nu$.

So as it could be seen, altough the proposed model gives some physically
reasonable results, it still requires more development.
Though not suitable for real-time systems, the model could still be used 
in material development process to quantitatively validate the material properties.
Furthermore, if only the purpose of a model is to describe bending of the IPMC in
the linear region, then considering the voltage drop on the electrode
would not contribute to the outcome significantly and could be omitted. 
%

\begin{figure}[t]
\begin{centering}
\includegraphics[clip,scale=0.34]{images/gmsh_example}
\caption{A tetrahedral mesh for an IPMC with electrodes. Generated using Gmsh.}
\label{fig:Gmsh_example}
\par\end{centering}

\end{figure}

\section{CONCLUSIONS}

We have developed a base model to model the simple physical processes
such as ion migration and electric field change in an IPMC. In this
paper we have extended the base model to three dimensions, added surface
electrodes and modeled the currents and voltage drop on them. Some
improvements in meshing techniques have been necessary to be able
to solve the three dimensional model within reasonable time, though
even more improvements are necessary. To get the full use of the three
dimensional model, the variable surface resistance part could be added
later without significantly increasing the complexity of the existing
model. The method used to calculate the currents in external circuit
was derived from Ramo-Shockley theorem, mostly known in other fields
of physics. Altough the developed model gives physically meaningful 
results, more optimizations are needed to improve the model in terms
of solution time and convergence. After the improvements, the
model will be applicable for variety of IPMCs to estimate the
current consumption and the voltage drop on the electrodes.


\begin{acknowledgment}
The authors from University of Nevada, Reno (UNR) acknowledge the
financial support from the U.S. National Science Foundation (Award\#:
0713075). We also acknowledge the financial support from the European
Science Foundation DoRa programme and CRDF-ETF joint grant no. ESE1-2902-TR-07.

\end{acknowledgment}

\bibliography{paper-pugal-2009-dscc}

\end{document}