\documentclass[12pt]{elsarticle}

%\usepackage[nomarkers,figuresonly]{endfloat}
\usepackage{endfloat}

%\usepackage{lineno,hyperref}


\usepackage{subfig}
\usepackage{multirow}
\usepackage{graphicx}
%\usepackage[autostyle]{csquotes}
\usepackage{csquotes}
\usepackage{leftidx}
\usepackage{amsmath}

 
%%%%  HYPERREF  HYPERREF  HYPERREF  HYPERREF  HYPERREF 
\usepackage[pdfauthor={Boris Jeremic}, 
            colorlinks=true,
            linkcolor=blue,
            citecolor=blue,
            urlcolor=blue]
            {hyperref}




\journal{Comput. Methods Appl. Mech. Engrg.}

%\bibliographystyle{elsarticle-num}
%\usepackage[numbers,round,colon,authoryear]{natbib}
\usepackage{natbib}

%% Water mark
%\usepackage{type1cm}
%\usepackage{eso-pic}
%\usepackage{color}

% \makeatletter
% \AddToShipoutPicture{%
%  \setlength{\@tempdimb}{.05\paperwidth}% X coord placement
%  \setlength{\@tempdimc}{.25\paperheight}% Y coord placement
%  \setlength{\unitlength}{1pt}
%  \put(\strip@pt\@tempdimb,\strip@pt\@tempdimc){
%   \makebox(-5,5){\rotatebox{90}{\textcolor[gray]{0.7}
%     {\fontsize{0.6cm}{0.5cm}\selectfont{DRAFT of Han et al, work in progress!}}}}
%  }
% }

%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% spacing REMOVE FOR official version
%
%
%baselinestretch
%
\renewcommand{\baselinestretch}{1.40}
\small\normalsize  % trick from Kopka and Daly p47.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{document}

\begin{frontmatter}

\title{Energy Dissipation Analysis of Elastic-Plastic Materials }
%\title{Energy dissipation analysis of elastic-plastic materials\tnoteref{mytitlenote}}
%\tnotetext[mytitlenote]{This research is supported by the US-DOE}

%% or include affiliations in footnotes:
\author[mymainaddress]{Han Yang}
%\ead{hhhyang@ucdavis.edu}

\author[mymainaddress]{Sumeet Kumar Sinha}

\author[mymainaddress]{Yuan Feng}

\author[mysecondaddress]{David B McCallen}

\author[mymainaddress,mysecondaddress]{Boris Jeremi{\'c}\corref{mycorrespondingauthor}}
\ead{jeremic@ucdavis.edu}

\cortext[mycorrespondingauthor]{Corresponding author}

\address[mymainaddress]{Department of Civil and Environmental Engineering, University of California, Davis, CA, USA}
\address[mysecondaddress]{Earth Science Devision, Lawrence Berkeley National Laboratory, Berkeley, CA, USA}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%\setcounter{tocdepth}{4}
%\tableofcontents
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
\begin{abstract}
Presented    is    an      energy  dissipation  analysis  framework for granular
material that is based on thermodynamics.
%
%Through  many  papers  on  seismic energy analysis of soil-structure-interaction
%(SSI) systems, the confusion between plastic work and plastic dissipation caused
%apparent  violation  of  the  second  law  of  thermodynamics.  A seismic energy
%analysis  framework  based  on  thermodynamics and computational geomechanics is
%presented.   
%
Theoretical    formulations    are    derived    from    the    second   law  of
thermodynamics,  in  conjunction  with  a  few  plausible  assumptions on energy
transformation and dissipation.
%
The  role  of  plastic  free  energy is emphasized by a  conceptual  experiment 
showing its physical nature.
%
Theoretical  formulation  is  adapted  in order to be applied in elastic-plastic
finite element method (FEM) simulations.
%
Developed  methodology  is  verified  through   comparison of input work, stored
energy, and energy dissipation of the system.
%
Separation  of  plastic  work  into  plastic  free energy and energy dissipation
removes  a  common  mistake,  made  in  a  number  of publications, where energy
dissipation can attain negative values (energy production) which is impossible.
% %
% Examples  presented show that, in most cases, more that 50\% of the input energy
% (coming  from, for example, seismic excitations) is dissipated due to plasticity
% on material level.
% %For materials with linear kinematic hardening,
% %linear  increase  of  plastic  dissipation  is  observed.  For  nonlinear
% %hardening  materials,  energy  is dissipated nonlinearly. Large-scale simulation
% %results on realistic SSI systems is being processed and will presented in future
% %papers.  The  proposed  method  is  capable  of  energy analysis for large-scale
% %realistic  SSI systems and hence can be used to improve structure and foundation
%designs on the aspect of proper seismic energy dissipation.
\end{abstract}

\begin{keyword}
Seismic  energy  dissipation  \sep  FEM  \sep  Computational  geomechanics  \sep
Thermodynamics \sep Elastic-plastic materials \sep Plastic free energy
\end{keyword}

\end{frontmatter}

%\linenumbers


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}

Energy  dissipation in elastic plastic solids and structures is  the result of an
irreversible dissipative process in which energy is transformed from one form to
another and entropy is produced.
%
The transformation and dissipation of energy is related to permanent deformation
and damage within an elastic-plastic material.
%
Of particular    interest  here  is the dissipation of mechanical energy that is input into
elastic-plastic solids by static or dynamic excitations.
% 
% Upon   deformation, elastic-plastic materials will absorb energy. 
% %
% Absorbed  energy  is    composed    of  the  elastic strain energy, plastic free
% energy and plastic dissipation.



Early work on plastic dissipation was done by \citet{Farren1925} and \citet{Taylor1934}.
 They performed experiments on
metals and proved that a large part, but not all, of the input mechanical energy
is  converted  into heat. The remaining part of the non-recoverable plastic work
is  known as the stored energy of cold work. The ratio of plastic work converted
into heating (Quinney--Taylor coefficient), usually denoted as $\beta$, has been
used   in  most later work on this topic.  Based on large number of experiments,
this   ratio   was   determined   to   be   a   constant   between  0.6  to  1.0
\citep{Clifton1984,Belytschko1991,Zhou1996,Dolinski2010,Ren2010,Osovski2013}.
%It  has  been  realized that this assumption is not valid in all cases, but it's
%simply  too  complicated to involve the evolution of Quinney--Taylor coefficient
%in thermomechanical constitutive models.


More recently  Rittel
\citep{Rittel2000a,Rittel2000b,Rittel2003}  published  several insightful papers
on  the  energy dissipation (heat generation) of polymers during cyclic loading,
presenting both experimental and theoretical works. 
%
Rosakis et al. \cite{Rosakis2000} presented
a   constitutive  model  for metals based  on  thermoplasticity that is able to 
calculate the evolution of energy dissipation.
%
Follow up papers
\citep{Hodowany2000,Ravichandran2002} present assumptions to
simply the problem. 
%
%One widely used assumption is the adiabatic condition, since
%air conducts heat much slower than metal. This assumption is reasonable in rapid
%monotonic  or cyclic loading (impact, vibration, earthquake). 
%
One direct application of
plastic dissipation  to  geotechnical   engineering   is   presented  by
Veveakis et al. \citep{Veveakis2007,Veveakis2012}, using thermoporomechanics to model
the heating and pore pressure increase in large landslides, like the 1963 Vajont
slide in Italy.


%OVDE


In  the  past  few  decades,  extensive  studies   have been conducted on energy
dissipation in structures and foundations. 
%
Work    by     \citet{Uang1990} has been considered a source and a reference for
many recent publications dealing with  energy  as  a  measure of structural demand.
% by many
%researchers
%\citep{Leger1992,Cosenza1993,Kalkan2007a,Kalkan2008,Symans2008,Gajan2011,
%Moustafa2011,Moustafa2014,Mezgebo2017,Deniz2017}.        
%
\citet{Uang1990} developed an  energy   analysis
methodology based
on  absolute  input  energy  (or  energy  demand).
Numerical  analysis  results  were  compared  with  experiments on a multi-story
building.  
%
In their work, Uang and Bertero \cite{Uang1990}, 
calculated hysteretic energy indirectly by taking
the   difference  of  absorbed  energy  and  elastic  strain  energy.  
%
The  term
{absorbed  energy}  of  each  time step is simply defined as restoring force times
incremental  displacement.  
%
It     is  also    stated    that   {hysteretic energy} is irrecoverable,  which
indicates   that   this   parameter  was  considered  the  same  as  {hysteretic
dissipation} or {plastic dissipation}.
%
An equation for
energy balance, is given by  (\citet{Uang1990}) as:
% 
\begin{equation}
E_i = E_k + E_{\xi} + E_a = E_k + E_{\xi} + E_s + E_h
\label{Eq1}
\end{equation}
%
where  $E_i$  is  the  (absolute)  input energy, $E_k$ is the (absolute) kinetic
energy,  $E_{\xi}$  is the viscous damping energy, $E_a$ is the absorbed energy,
which is composed of elastic strain energy $E_s$ and hysteretic energy $E_h$.

The  problem  with  this approach  is  the absence of plastic free energy, which is
necessary  to correctly evaluate energy dissipation of elastic-plastic materials
and  to  uphold  the  second  law of thermodynamics. While there is no direct plot of
plastic  dissipation (hysteretic energy) in \cite{Uang1990}, since it was not defined
directly, there  are  plots of other energy components. Plastic dissipation
can  be  easily  calculated from  these  plots. 
%
After doing this,
indications of negative incremental energy dissipation, which violates the basic
principles of thermodynamics, were found in various sections of the paper.


This   misconception  could  be  clarified  by  renaming  hysteretic  energy as
{plastic  work}, a sum  of plastic dissipation and
plastic  free  energy.  Both  plastic  work  and  plastic  free  energy  can  be
incrementally  negative,  but  plastic dissipation (defined as the difference of
plastic  work and plastic free energy) must be incrementally non-negative during
any  time  period.  Unfortunately, this misconception has been inherited (if not
magnified)  by many following  studies on energy analysis of earthquake
soils and structures (hundreds of papers).



%   %  % 
%   %   In the past few decades, extensive
%   %  studies  have been conducted on energy dissipation in structures and foundations
%   %  \cite{Uang1990},   \cite{Leger1992},  \cite{Symans1998},
%   %  \cite{Soong2002},    \cite{Symans2008},    \cite{Wong2008},
%   %  \cite{Nehdi2010}.
%   %
%   A numbe
%   Despite  different methodology used,  the calculations of energy  dissipation due to
%   hysteretic damping (material plasticity) are all performed using  stress-strain   and/or
%   force-displacement  response  curves. 
%   %
%   Traditionally in elastic-plastic computations, 
%   it  is  assumed that all plastic deformation will lead to energy dissipation.
%   %
%   In  other  words,  plastic  energy  dissipation is calculated by integrating the
%   product  of  generalized  plastic  strain and generalized stress. 
%   %
%   In  some  cases,  researchers    simply   take the area of hysteresis loop as an
%   evaluation of energy dissipation.  


%   %
%   However, close examination of results (in a number  of  plots)  from the previously
%   mentioned  papers reveals  that  accumulated  (total)  energy  dissipation was
%   observed to decrease during dynamic loading.
%   %
%   This is a clear violation of the second law of thermodynamics.
%   %
%   It  also  means that the approach used to compute energy dissipation is
%   inappropriate.
%   




%OVDE





The basic principles of
thermodynamics  are frequently used to derive new constitutive models, for example by
\citet{Dafalias1975},                   \citet{Ziegler1987},
\citet{Collins97},                 
\citet{Houlsby2000}, \citet{Collins2002a},
\citet{Collins2002b}, \citet{Collins2003b} and
\citet{Feigenbaum2007}.
%
%In these studies, the definition of plastic energy dissipation is clarified.
%
The  concept  of  plastic free energy is introduced to enforce the second law of
thermodynamics for developed constitutive models. 
%
It is important to distinguish  between  energy  dissipation  due to plasticity
and plastic work, which is often a source of a confusion.
%
The    physical   nature  of  plastic  free energy is illustrated  later in this
 paper through a conceptual example that is analyzed on particle scale.
%
Essentially,   development   of  plastic  free  energy  is  caused  by  particle
rearrangement in granular assembly under external loading.

Specific formulation of free energy depends on whether the elastic and plastic behavior
of the material is coupled.
%
According      to      Collins     et     al.     \cite{Collins97},
\cite{Collins2002a},
\cite{Collins2003b},  material  coupling
behavior  can be divided into {modulus coupling}, where the instantaneous
elastic  stiffness  (or  compliance)  moduli  depend  on the plastic strain, and
{dissipative  coupling},  where  the rate of dissipation function depends
not  only  on  the  plastic  strains  and  their rates of change but also on the
stresses (or equivalently the elastic strains).
%
The modulus coupling  describes the degradation of
stiffness  as in for  rock and concrete, and is usually modeled by employing a coupled
elastic-plastic  constitutive  model  or  by  introducing  damage variables. 
%
The  dissipative  coupling    is  considered  to  be one of the main reasons for
non-associative    behavior   in   geomaterials
\cite{Collins97}, \cite{Ziegler1981}.

%
%In  this  paper,  various elastic-plastic material
%models  are  applied  and  tested,  involving  both coupled stiffness tensor and
%non-associative plasticity.
%

A  number  of  stability  postulates are commonly used to prevent
violation of principles of thermodynamics.
%
Stability postulates include Drucker's
stability                     condition            
\cite{Drucker1956}, \cite{Drucker1957},   Hill's   stability   condition
\cite{Bishop1951}, \cite{Hill1958},  and  Il'Iushin's  stability  postulate
\cite{Il1961}, \cite{local-63}.
%
As  summarized  in  a  paper  by Lade
\cite{Lade2002},
theoretical         considerations        by        
\citet{Nemat-Nasser83} and \citet{local-25} have suggested that they are
sufficient  but not necessary conditions for stability. 
%
These  {stability
postulates}  can  indeed  ensure the admissibility of the constitutive models by
assuming  certain  restrictions  on  incremental  plastic work. As demonstrated by
\citet{Collins2002a}, if the plastic strain rate is replaced by
the  irreversible  stain  rate  in  Drucker's  postulate,  then all the standard
interpretations  of  the  classical  theory  still  apply for coupled materials.
\citet{Dafalias77b}  also  modified  Il'Iushin's  postulate in a
similar way and applied it to both coupled and uncoupled materials.





%Other  than  theoretical  considerations,  deficiency in the current practice of
%seismic energy analysis also exists on numerical modeling. A number of numerical
%studies  on  SSI  systems  have  been  conducted  with  significant  results and
%inspirations \cite{Chen1977},
%\cite{Makris1994},  \cite{McKenna97}, \cite{Bielak2001},
%\cite{Sett2005a},   \cite{Sett2005b},
%\cite{Fenves1998},    \cite{Elgamal2008},   \cite{Jeremic2008d},
%\cite{Jeremic2013},                          \cite{Jeremic2004g},
%\cite{Boulanger1999},  \cite{Yang2003}.
%
% HAN:
%
% What are you trying to say with the above references?
%
% Han's comment:
% This part is not useful, should be deleted...
%




%
%
%However,  among all
%the papers on SSI problems, most of them have documented over-simplified or even
%no  energy  analysis in their cases. 
%

It is important to note that development of inelastic deformation in geomaterials
involves large changes in entropy, and significant energy dissipation.
%
It is thus useful to perform energy dissipation (balance) analysis for all 
models with inelastic deformation.
%
%Thus,  it  is  necessary to perform energy dissipation analysis in all numerical
%studies  of SSI systems. Unfortunately, little effort has been done regarding to
%this topic. 
In this paper we focus on  energy dissipation
on  material  level. 
%More results of large-scale modeling is being
%processed and will be discussed in our following papers.
%
Focus is on proper modeling that follows thermodynamics. 
% 
Comparison  is  made between accumulated  plastic  dissipation and accumulated 
plastic work, since these quantities can
be  quite different  in most cases. 
%Depending on the material constitutive model
%used  in simulation, the energy dissipation accumulates with unique patterns. 
As
a  way of verification, the input work, which is introduced by applying external
forces, is compared with the stored energy and dissipation in the entire system.
Finally,  conclusions  on plastic energy dissipation are drawn from the verified
results.




%  Han's comment:
%  I saw many papers have a paragraph like this to explain the contents of the paper
%
%  
%  In  this  study,  a  seismic  energy  analysis framework is established based on
%  thermodynamics  and  computational geomechanics. Different mechanisms of seismic
%  energy  transformation and dissipation are identified and analyzed. Derived from
%  the  second  law  of  thermodynamics and reasonable assumptions, formulations of
%  energy  parameters  in  multiple  forms are presented. In order to apply them to
%  numerical  simulations,  these  theoretical  equations  are  modified to satisfy
%  computational  requirements  of  FEM.  To  show  the  effectiveness of our newly
%  developed  methodology,  a  series  of  numerical  simulations  under  different
%  boundary conditions are conducted with emphasize on energy dissipation analysis.
%  In  particular,  values  of  accumulated  plastic  dissipation are compared with
%  accumulated  plastic  work. It can be observed later that the two quantities can
%  be  very  different  in most cases. Depending on the material constitutive model
%  used  in simulation, the energy dissipation accumulates with unique patterns. As
%  a  way of verification, the input work, which is introduced by applying external
%  forces, is compared with the stored energy and dissipation of the entire system.
%  Finally,  conclusions  on plastic energy dissipation are drawn from the verified
%  results.
%  




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Theoretical and Computational Formulations}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
\subsection{Thermo-Mechanical Theory}

For  energy  analysis  of elastic-plastic materials undergoing {isothermal
process},  it  is beneficial to start from the statement of the first and second
laws of thermodynamics:
%
\begin{equation}
\hat{W} = \dot{\Psi} + \Phi  
%, \quad \mathrm{where} \; \Phi \ge 0
\end{equation}
%
where $\Phi \ge 0$ and   $\hat{W}    \equiv   \sigma   :   \dot{\epsilon}   =   \sigma_{ij}
\dot{\epsilon}_{ij}$  is  the  rate  of  work  per unit volume. 
%
The function $\Psi$ is the {Helmholtz free
energy},  and  $\Phi$  is  the rate of dissipation; both defined per volume. The
free  energy $\Psi$ is a function of the {state variables} (also known as
 {internal  variables}),  but  $\Phi$  and  $\hat{W}$  are  not  the  time
derivatives of the state functions. 
%
The choice of state variables depends on the
complexity  of  constitutive  model  that  is being used, as cyclic loading with
certain  hardening behaviors usually requires more state variables. This will be
elaborated  in  the  following  sections  as  we discuss specific elastic-plastic
material  models.  Note that in this paper all stresses are defined as effective
stresses. In order to avoid confusion, the common  notation ($\sigma_{ij}^{'}$) will not be used. 
Standard definition of stress from mechanics of materials, i.e. positive in tension, is used.



For  general  elastic-plastic  materials,  the  free  energy  depends on both the
elastic  and  plastic  strains. In most material models, it can be assumed
that the free energy $\Psi$ can be decomposed into elastic and plastic parts:
%
%
% HAN: is this Phi or Psi in the equation below, I am skipping editing 
%
% Han's comment:
% Should be Psi, this was a typo. Corrected.
%
\begin{equation}
\Psi = \Psi_{el} + \Psi_{pl}
\label{Eq2}
\end{equation}

The  total  rate  of work associated with the effective stress can be written as
the sum of an elastic and plastic component:
%
% HAN: here you use effective stress, while it should be just stress, right?
%
% Han's comment:
% Plastic energy dissipation comes from inter-particle frictions, which should be associated 
% with effective stresses. So, as stated in line 370, all stresses are defined as effective
% stresses. In order to avoid confusion, the common  notation ($\sigma_{ij}^{'}$) will not be used. 
%
\begin{equation}
\hat{W}^{el} \equiv \sigma_{ij} \dot{\epsilon}^{el}_{ij} = \dot\Psi_{el}
\end{equation}
%
and
% 
\begin{equation}
\hat{W}^{pl} \equiv \sigma_{ij} \dot{\epsilon}^{pl}_{ij} = \dot\Psi_{pl} + \Phi
\label{Eq4}
\end{equation}


%
%
%
In  the  case  of a  {decoupled material},  the elastic free energy $\Psi_{el}$ depends only on the
elastic  strains,  and the plastic free energy $\Psi_{pl}$ depends only on the plastic strains, as shown
by \citet{Collins97}:
%
\begin{equation}
\Psi = \Psi_{el}(\epsilon^{el}_{ij}) + \Psi_{pl}(\epsilon^{pl}_{ij})
\label{Eq5}
\end{equation}

The effective stress can also be decomposed into two parts:
\begin{equation}
\sigma_{ij} = \alpha_{ij} + \chi_{ij}, \quad \text{where} \; \alpha_{ij} \equiv \frac{\partial \Psi_{pl}}{\partial \epsilon^{pl}_{ij}} \quad \text{and} \quad \chi_{ij} \equiv \frac{\partial \Phi}{\partial \dot{\epsilon}^{pl}_{ij}}
\label{Eq6}
\end{equation}

The  stress  tensors  $\alpha_{ij}$  and $\chi_{ij}$ are termed the shift (drag,
back or quasi-conservative) stress and dissipative stress respectively.

Ziegler's orthogonal postulate \cite{Ziegler1987} ensures the validity
of  Equation~\ref{Eq6}.  It  is  equivalent  to  the maximum entropy production
criterion, which is necessary to obtain unique formulation. Also, this is a weak
assumption  so  that  all  the  major  continuum  models of thermo-mechanics are
included.  Equation~\ref{Eq4}  of  plastic  work  rate  can hence be
rewritten as:
%
\begin{equation}
\hat{W}^{pl} \equiv \sigma_{ij} \dot\epsilon^{pl}_{ij} = \dot\Psi_{pl} + \Phi = \alpha_{ij} \dot\epsilon^{pl}_{ij} + \chi_{ij} \dot\epsilon^{pl}_{ij}
\end{equation}



%
% CONTINUING editing from here, until we correct Psi or Phi above
%


The   {plastic work} is the product of the  {true stress} with
the  plastic  strain rate, while the  {dissipation rate} is the product of
the   {dissipative  stress}  with  the  plastic strain rate. They are only
equal  if  the  back stress is zero, or equivalently, if the free energy depends
only on the elastic strains.

In  kinematic  hardening models, where the back stress describes the translation (or rotation)
 of  the yield surface, the decomposition of the true stress (sum of back
stress  and  dissipative  stress) is a default assumption. Although such a shift
stress      is      important      for      anisotropic     material     models,
\citet{Collins2002b} have pointed out that it is also necessary in
isotropic  models  of  geomaterials  with  different  strength  in  tension  and
compression.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
\subsection{Plastic Free Energy}

A  popular  conceptual model, which focused on particulate materials and demonstrated the
physical     occurrence     of     shift     stresses,    was    described    by
\citet{Besseling1994}        and \citet{Collins2002b}.  On  macro (continuum)  scale, every
point  in  a  given element is at yield state and deforms plastically. But on
meso-scale,  only  part  of this element is undergoing plastic deformations, the
remaining  part  is  still  within  yield  surface  and respond elastically. The
elastic  strain  energy  stored  in  the  elastic part of a plastically deformed
macro-continuum  element  is considered to be locked into the macro-deformation,
giving  rise  to the plastic free energy function $\Psi_{pl}$ and its associated
back  stress  $\alpha_{ij}$.  This  energy can be released only when the plastic
strains are reversed.

For  better  explanation,  the  nature  of  plastic  free energy in particulate materials 
is  illustrated  through  a  finite  element simulation combined with
considerations    of    particle  rearrangement on mesoscopic scale. 
Figure~\ref{example_DP_AF_cyclic_response} shows stress-strain response of 
Drucker-Prager with nonlinear Armstrong-Frederick kinematic hardening, a typical 
elastic-plastic model for metals and geomaterials.
% As shown in
% Figure~\ref{example_DP_AF_cyclic_response},  the  stress-strain  curves  match the
% typical
% response  of  an  elastic-plastic material under cyclic loading. Numerically, the
% constitutive  model  is  Drucker-Prager  yield function with Armstrong-Frederick
% nonlinear  kinematic  hardening.
% which is usually used to simulate behavior of
% geotechnical  materials.  
Six  states  during  shear  are  chosen  to  represent
evolution of micro fabric of the numerical sample.
%
Correspondingly,  Figure~\ref{2D_assembly_cyclic}  shows  the process of particle
rearrangement of the 2D granular assembly under cyclic shearing from microscopic
level.  The  square  window can be roughly considered as a representative volume (a constitutive level or a finite element)
in FEM. 
%Note
%that  this  example  only  shows  one  arbitrary  result  of  infinite  particle
%arrangements, but it is representative nonetheless.

\begin{figure}[!htbp]
 \centering
 \subfloat[]{\includegraphics[width=0.49\columnwidth]{Figures/Theo_Stress-Strain.pdf}}
 \subfloat[]{\includegraphics[width=0.49\columnwidth]{Figures/Theo_Stress-Plastic_Strain.pdf}}
 \caption{
 \label{example_DP_AF_cyclic_response}
 Elastic-plastic material modeled with Drucker-Prager   yield  function  and  
 Armstrong-Frederick  kinematic hardening  under  cyclic  shear loading: (a) Stress-strain curve; 
 (b) stress and plastic strain versus time.}
\end{figure}

\begin{figure}[!htbp]
\centering
\subfloat[]{\includegraphics[width=0.4\columnwidth]{Figures/Plastic_Free_Energy_1.pdf}}
\quad
\subfloat[]{\includegraphics[width=0.4\columnwidth]{Figures/Plastic_Free_Energy_2.pdf}}
\\
\subfloat[]{\includegraphics[width=0.4\columnwidth]{Figures/Plastic_Free_Energy_3.pdf}}
\quad
\subfloat[]{\includegraphics[width=0.4\columnwidth]{Figures/Plastic_Free_Energy_4.pdf}}
\\
\subfloat[]{\includegraphics[width=0.4\columnwidth]{Figures/Plastic_Free_Energy_5.pdf}}
\quad
\subfloat[]{\includegraphics[width=0.4\columnwidth]{Figures/Plastic_Free_Energy_6.pdf}}
\caption{
\label{2D_assembly_cyclic}
Particle rearrangement of a 2D granular assembly under cyclic shearing:
(a)  Initial  state;  (b) Loading (accumulating plastic free energy); (c) End of
loading  (maximum  plastic  free  energy);  (d)  Unloading  (plastic free energy
unchanged);  (e)  Reverse  loading  (releasing  plastic free energy); (f) End of
reverse loading (plastic free energy released).}
\end{figure}

By        discussing     movement    and    energy    of    particle    A   in  
Figure~\ref{2D_assembly_cyclic}, the physical nature of plastic free energy is
illustrated.  At state (a), which is the beginning of deformation, particle A does not
bear    any    load    other  than  its  self  weight. State (b) is in middle of
loading, when particle B pushes downwards to particle A until it makes contact
with  particle  D  and  E.  Load reaches peak at state (c), and there's no
space  for  particle  A  to  move. Then the sample is unloaded to state (d).
 Particle A is now  {stuck} between particles C, D,
and F, which means that certain amount of elastic energy is stored due to particle elastic
deformation.  Compared  with  state  (a),  this  part  of  elastic energy is not
released  when  the  sample  is  unloaded,  which indicates that it's  not classic
strain energy. {This part of elastic energy on particle level which can't
be  released  by  unloading  is  defined  as the plastic free energy in granular
materials.}  Reverse loading starts at state (e), where particle  D pushes
particle A upwards, making it  {squeeze through} particle C and F.
Elastic energy on particle level, which is now defined as plastic free energy, is
released during reverse loading.


By  analyzing  this  example,  an explanation on particle scale is provided for the origin of
plastic  free  energy in granular materials. It is important
to  note   that  the  concept of plastic free energy also exists in metals and
other     materials,     as    studied    by    \citet{Dafalias2002}    and
\citet{Feigenbaum2007}. The physical nature of plastic free energy
in  these  materials  can  be  different and probably should be studied on
molecular and/or crystalline  level.

Collins \cite{Collins2002b}, \cite{Collins2003b} suggested
that  in  the  case  of  granular  materials,  the  particle-level plastic energy
dissipation  during  normal  compaction,  arises  from  the  plastic deformations
occurring  at  the  inter-granular contacts on the strong force chains, that  are
bearing  the bulk of the applied loads. Collins also suggested that the  {locked-in} elastic
energy  is produced in the weak force networks, where the local stresses are not
large  enough  to produce plastic deformation at the grain contacts. The plastic
strains  can be associated with the irreversible rearrangement of the particles,
whilst  the  elastic  energy arises from the elastic compression of the particle
contacts.  Part of this elastic strain energy will be released during unloading,
however other  part  of this energy will be trapped as a result of the irreversible
changes  in  the  particle  configuration.  
% This idea is best illustrated using
% Discrete  Element  Method  (DEM), where granular assemblies are modeled
% with  discrete  elements  and contacts with elastic-plastic contact models.
%More  detailed  explanations  of  DEM  and  its applications can be found in the
%papers        by       \cite{Cundall1979},       \cite{Radjai1996},
%\cite{Thornton2000}, \cite{Yang2015}, and \cite{Yang2016}.






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Plastic Dissipation}


As  pointed  out,  plastic work and energy dissipation are not the same physical
quantity. The confusion of these two
concepts  often  leads to incorrect results and conclusions, especially in seismic
energy  dissipation analysis.  Of major concern in this paper is the computation
of plastic dissipation, as elaborated in this section.


With   the  decoupling  assumption  (Equation~\ref{Eq5}),  the  second  law  of
thermodynamics  (positive entropy production) directly leads to the dissipation
inequality, which states that the energy dissipated due to the difference of the
plastic  work  rate  and the rate of the plastic part of the free energy must be
non-negative:
%
\begin{equation}
\Phi = 
\sigma_{ij} \dot{\epsilon}^{pl}_{ij} - \dot{\Psi}_{pl} = 
\sigma_{ij} \dot{\epsilon}^{pl}_{ij} - \rho \dot{\psi}_{pl} \ge 0
\label{Eq8}
\end{equation}
%
where $\dot{\psi}_{pl}$ is the rate of plastic free energy, per unit volume, and
$\rho$  is  the  density. In addition, $\psi_{pl}$ denotes plastic free energy density,
which is  generally not constant at different locations in a body. This expression is closer
to physics and makes it convenient for further derivations.



Now  we proceed to consider how to calculate plastic free energy, which can then
be   used   to  calculate  dissipation.  According  to  
\citet{Feigenbaum2007},  plastic  free  energy  density $\psi_{pl}$ is
assumed  to  be  additively  decomposed  into  parts  which  correspond  to  the
isotropic, kinematic and distortional hardening mechanisms as follows:
%
\begin{equation}
\psi_{pl} = 
\psi_{pl}^{iso} + 
\psi_{pl}^{ani}; 
\quad \psi_{pl}^{ani} = 
\psi_{pl}^{kin} - \psi_{pl}^{dis}
\end{equation}
%
where     $\psi_{pl}^{iso}$,     $\psi_{pl}^{ani}$,     $\psi_{pl}^{kin}$,    and
$\psi_{pl}^{dis}$  are  the  isotropic,  anisotropic, kinematic, and distortional
parts  of  the plastic free energy, respectively. The anisotropic part is assumed to decompose
into  kinematic  and distortional parts, which correspond to different hardening
models.  The  subtraction,  instead  of addition, of $\psi_{pl}^{dis}$ from
$\psi_{pl}^{kin}$, to obtain the overall anisotropic part $\psi_{pl}^{ani}$ of the
plastic  free  energy,  is  a  new  concept  proposed  by 
\citet{Feigenbaum2007}.  This  expression  can better fit experimental
data,  as  well  as  satisfy  the  plausible  expectations  for  a limitation of
anisotropy development.



As pointed out by \citet{Dafalias2002},  the  thermodynamic
conjugates  to  each  of  the internal variables exist and each part of the
plastic  free  energy  can be assumed to be  only  a  function  of  these  
conjugates. The explicit expressions for the isotropic and kinematic components 
of the plastic free energy are:
%
\begin{equation}
\psi_{pl}^{iso} = 
\psi_{pl}^{iso}(\bar{k}) = 
\frac{\kappa_1}{2 \rho} \bar{k^2}; 
\quad \psi_{pl}^{kin} = 
\psi_{pl}^{kin}(\bar{\alpha}_{ij}) = 
\frac{a_1}{2 \rho} \bar{\alpha}_{ij} \bar{\alpha}_{ij}
\label{Eq10}
\end{equation}
%
where $\bar{k}$ and $\bar{\alpha}_{ij}$ are the
thermodynamic  conjugates  to  $k$  (size  of  the yield surface) and $\alpha_{ij}$
(deviatoric back stress tensor representing the center of the yield surface),  respectively.
Material constants $\kappa_1$ and $a_1$ are  non-negative  material  constants whose 
values depend on the choice of elastic-plastic material models.   

According to  definition, the thermodynamic conjugates are related to the
corresponding internal variables by:
%
\begin{equation}
k = 
\rho\frac{\partial \psi_{pl}^{iso}}{\partial \bar{k}} = 
\kappa_1 \bar{k}; 
%
\quad \alpha_{ij} = 
\rho\frac{\partial \psi_{pl}^{kin}}{\partial \bar{\alpha}_{ij}} = 
a_1 \bar{\alpha}_{ij}
\label{Eq11}
\end{equation}
%
By substituting Equation~\ref{Eq11} back into Equation~\ref{Eq10}, the plastic free
energy can be expressed in terms of the internal variables:
%
\begin{equation}
\psi_{pl}^{iso} = 
\frac{1}{2 \rho \kappa_1} k^2; 
%
\quad \psi_{pl}^{kin} = 
\frac{1}{2 \rho a_1} \alpha_{ij} \alpha_{ij}
\label{Eq12}
\end{equation}

With Equation~\ref{Eq12}, the components of plastic free energy can be computed,
as  long as the internal variables are provided. 
Combining Equation~\ref{Eq8} with \ref{Eq12}, the plastic dissipation in a given 
elastic-plastic material can be accurately obtained at any location, at any time.
This  approach  allows  engineers and designers  to correctly identify energy 
dissipation in time and space and make appropriate conclusions on material behavior.

% In the above presented methodology, the important role
% of  plastic free energy is recognized and applied to numerical modeling studies.
% %Compared  to earlier energy analysis approaches, one major advantage is that the
% Using this approach, plastic  dissipation  of an SSI model can be obtained at any
% location, at any time during  simulation.
% 
%This  allows  us  to  identify  components with high energy
%dissipation,  which  indicates  significant  deformation or damage, from a large
%numerical  model.  Based  on this information, structures can be better designed
%against seismic or other types of dynamic loading.





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Energy Computation in Finite Elements}




Formulations  from  the previous section are applied to FEM analysis in order
to follow energy dissipation. 
%
%Since  the  intent  is  to  analyze the  location of energy dissipation in large
%Earthquake Soil Structure Interaction (ESSI) models,
Energy density is chosen as the physical parameter for energy analysis.
%
Energy  density  in  this study is  defined  as the amount of energy stored in a
given region of space per unit volume.


For  FEM  simulations, both external forces and displacements can be prescribed.
The   finite   element program accepts either (or both) forces and/or displacements as input and
solves for the other.
Either  way,  the   {rate  of  input  work}  can  be  calculated by simply
multiplying  force  and displacement  within a time step. Therefor  {input
work} of a finite element model is:
%
\begin{equation}
W_{Input}(t) = 
\int_{0}^{t} \dot W_{Input}(T) dT = 
\int_{0}^{t} \sum_{i} F_{i}^{ex}(\boldsymbol{x},T) \dot{u}_{i}(\boldsymbol{x},T) dT
\label{Eq_input_work}
\end{equation}
%
where $F_i^{ex}$ is the external force and $u_i$ is the displacement computed at
the location of the applied load, at given time step, for a load controlled analysis.
The  external  load can have many forms, including nodal loads, surface loads,
and  body  loads.  All  of  them  are  ultimately transformed into nodal forces.  As  shown
in Equation~\ref{Eq_input_work}, input
work  is  computed  incrementally  at  each  time  step,  in order to obtain the
evolution of total input work at certain time.



As      shown     in    Figure~\ref{energy_classification},  when  loads  and/or
displacements are introduced
into   a finite element model, the input energy will be converted in a number of
different forms as it propagates
through  the  system.  
%
Input energy will be converted into {kinetic
energy},   {free  energy},  and  {dissipation}. As mentioned before,
free  energy  can be further separated into elastic part, which is traditionally
defined as strain energy, and plastic part, which is defined as the plastic free
energy.  Kinetic  energy  and strain energy can be considered as the recoverable
portion  of  the total energy since they are transforming from one to
another. Plastic free energy is more complicated in
the  sense  that  it is conditionally recoverable during reverse loading, as has
been  discussed  in  detail  in previous sections. Other than kinetic energy and
free    energy,    the  rest  of  the  input  energy  is dissipated, transformed
into heat or other forms of energy that are irrecoverable.


\begin{figure}[!htbp]
\centering
\includegraphics[width=0.7\columnwidth]{Figures/Energy_Transformation.pdf}
\caption{\label{energy_classification}  
% FIGURE  TO  BE  CHANGED  SO  THAT  IS A REGULAR FINITE ELEMENT, WITH INPUT LOADS
% AND/OR DISPLACEMENTS, AND THEN THE REST (DISSIPATION) IS THE SAME...
%
Different  forms  of energy in a dynamic soil-structure system.}
\end{figure}



Calculation of  kinetic energy and strain energy is rather straight
forward:
%

\begin{equation}
U_K(\boldsymbol{x},t) = 
\frac{1}{2} \rho \dot{u}_{ij}(\boldsymbol{x},t) \dot{u}_{ij}(\boldsymbol{x},t)
\label{Eq_kinetic_energy_density}
\end{equation}

\begin{equation}
U_S(\boldsymbol{x},t) = 
\int_{0}^{t} \dot {U_S}(\boldsymbol{x},T) dT = 
\int_{0}^{t} \sigma_{ij}(\boldsymbol{x},T) \dot \epsilon_{ij}^{el}(\boldsymbol{x},T) dT
\label{Eq_strain_energy_density}
\end{equation}
%
where  $U_K$  and $U_S$ are the  {kinetic energy density} and  {strain energy   
density}, respectively.
%

%
Similar  to  the input energy, strain energy density and plastic free energy are
also  computed  incrementally. Integrating energy density over the entire model,
corresponding energy quantities are expressed as:
%
\begin{eqnarray}
E_K(t) = \int_V U_K(\boldsymbol{x},t) dV
\label{EK_kinetic_energy}
\end{eqnarray}

\begin{eqnarray}
E_S(t) = \int_V U_S(\boldsymbol{x},t) dV 
\label{ES_strain_energy}
\end{eqnarray}


\begin{eqnarray}
E_{P}(t) = \int_V \Psi_{pl}(\boldsymbol{x}, t) dV
\label{EP_plastic_free_energy}
\end{eqnarray}

where $E_K$, $E_S$, and $E_P$ are the kinetic energy, strain energy, and plastic
free  energy  of  the  entire  model, respectively. Energy densities, defined in
Equations~\ref{Eq_kinetic_energy_density} and \ref{Eq_strain_energy_density}  are functions of
both    time  and space, while energy components, defined in the above equations
(Equation~\ref{EK_kinetic_energy}, \ref{ES_strain_energy}, and \ref{EP_plastic_free_energy}) are only functions of time, since
they  are  integrated  over  the  whole  model.  
%
% It is more intuitive to compare
% results  between  different  models  by the total energy, since it is one scalar
% quantity.  However,  in  order to understand the distribution of energy inside a
% model, it is necessary to analyze energy density.

Although  the  plastic  free  energy  is  conditionally recoverable, it is still
considered  to  be  stored in the system, rather than dissipated. Summing up all
the stored energy $E_{Stored}$, one obtains:
%
\begin{equation}
E_{Stored} = E_K + E_S + E_P
\end{equation}

Rate  of plastic dissipation, given by Equation~\ref{Eq8}, can be integrated over
time and space:
%
\begin{equation}
D_P(t) = \int_V \int_{0}^{t} \Phi(\boldsymbol{x},T) dT dV
\end{equation}
%
where  $D_P$ is the dissipation due to plasticity of the entire model at certain
time.

Finally the energy balance of a finite element model is given by:
%
\begin{equation}
W_{Input} = E_{Stored} + D_P = E_K + E_S + E_P + D_P
\end{equation}

%Another important issue is how to choose the location for energy density calculation, which is originated from FEM. As shown in Figure~\ref{energy_density_calculation}, in this case we choose to use the Gauss points. Two aspects are considered in this decision. Firstly, we have the solution for stress and strain on the Gauss points, which makes it straight forward to calculate the strain energy density of this location. Secondly, in our 1D beam case with 8-node-brick elements, displacements at Gauss points can be linearly interpolated from nodal displacements, which we already calculated and stored in output files. Then we just need to take its derivative to get velocity, thus the kinetic energy density.

%\begin{figure}
% \centering
% \includegraphics[width=0.5\columnwidth]{energy_density_calculation.pdf}
% \caption{\label{energy_density_calculation}
% Position to calculate energy density components in 1D case with 8-node-brick elements}
%\end{figure}









%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Numerical Studies}



Numerical simulation results presented in this paper are performed using
the      Real      ESSI     (Real     Earthquake Soil Structure     Interaction)
Simulator \cite{Real_ESSI_Simulator}. 
% Real  ESSI  Simulator  is  a software,
% hardware  and  documentation  system  for  high fidelity, high performance, time
% domain,   nonlinear,   3D,   finite   element   modeling   and   simulation   of
% earthquake-soil/rock-structure  interaction  of  Infrastructure Objects (Nuclear
% Facilities,  Nuclear Power Plants, Dams, Bridges, Buildings). Real ESSI contains
% options  for  various  materials,  elements,  constrains,  loads,  and computing
% algorithms,  especially newly developed constitutive models for geotechnical and
% structural components. These features make Real ESSI highly effective on seismic
% energy dissipation studies, which is important yet rare to exist in many classic
% seismic numerical simulation frameworks. The proposed methodology is implemented
% and verified in the Real ESSI Simulator.
% 
Examples  in  this  paper  focus  on  constitutive  behavior  of elastic-plastic
material from the perspective of energy dissipation. All cases are modeled with
solid brick elements, using static, load control analysis.
% Static algorithm is employed
%to  eliminate the influence of inertial forces and dynamic reactions, which will
%be the focus of our following papers.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Elastic Material}


\begin{figure}[!htbp]
 \centering
 \subfloat[]{\includegraphics[width=0.3\columnwidth]{Figures/Model_Single_Block.pdf}}
 \subfloat[]{\includegraphics[width=0.6\columnwidth]{Figures/Model_Cantilever.pdf}}
 \caption{\label{figure_models} 
 Numerical models used in this paper: (a) Single brick element; (b) Cantilever with 10 brick elements.}
\end{figure}


Initial investigation of energy dissipation is focused on linear elastic material.
%
It  is  noted  that  linear  elastic  material  does not dissipate energy. 
%
However,  use  of  linear  elastic  material  model is suitable for  preliminary
verification of the newly developed energy analysis methodology. In this section, 
energy balance in a single brick element and a cantilever beam is studied,
as shown in Figure~\ref{figure_models}.



The simplest case is a single element model under uniform shear load. The model is
constrained  appropriately  to  simulate simple shear test. In order to show the
influence  of  different material parameters and loads, a set of simulations are
performed      and      the      results      are     presented     in     
Table~\ref{table_linear_elastic_single}  and Figure~\ref{figure_linear_elastic}. 



\begin{table}[htb]
 \centering
 \caption{Energy analysis results for linear elastic materials (single element).}
 \label{table_linear_elastic_single}
 \resizebox{\textwidth}{!}{
  \begin{tabular}{c|cc|cc|cccc|c}
   \hline
   \multirow{2}{*}{Case} & \multicolumn{2}{c|}{Material Property} & \multicolumn{7}{c}{Simulation Results}                                                                                                                                        \\ \cline{2-10}
   & $E$ (GPa) & $\nu$ & $u$ (m) & \boldmath$W_{Input}$ (J) & $E_{K}$ (J) & $E_{S}$ (J) & $E_{P}$ (J) & \boldmath$E_{Stored}$ (J) & \boldmath$D_P$ (J) \\
   \hline
   1a & 100 & 0.30 & 2.60E-5 & 13.00 & 0.00 & 13.00 & 0.00 & 13.00 & 0.00 \\
   1b & 150 & 0.30 & 1.73E-5 & 8.67 & 0.00 & 8.67 & 0.00 & 8.67 & 0.00 \\
   1 & 200 & 0.30 & 1.30E-5 & 6.50 & 0.00 & 6.50 & 0.00 & 6.50 & 0.00 \\
   1c & 250 & 0.30 & 1.04E-5 & 5.20 & 0.00 & 5.20 & 0.00 & 5.20 & 0.00 \\
   1d & 300 & 0.30 & 8.67E-6 & 4.33 & 0.00 & 4.33 & 0.00 & 4.33 & 0.00 \\
   \hline
   1e & 200 & 0.20 & 1.20E-5 & 6.00 & 0.00 & 6.00 & 0.00 & 6.00 & 0.00 \\
   1f & 200 & 0.25 & 1.25E-5 & 6.25 & 0.00 & 6.25 & 0.00 & 6.25 & 0.00 \\
   1 & 200 & 0.30 & 1.30E-5 & 6.50 & 0.00 & 6.50 & 0.00 & 6.50 & 0.00 \\
   1g & 200 & 0.35 & 1.35E-5 & 6.75 & 0.00 & 6.75 & 0.00 & 6.75 & 0.00 \\
   1h & 200 & 0.40 & 1.40E-5 & 7.00 & 0.00 & 7.00 & 0.00 & 7.00 & 0.00 \\
   \hline
  \end{tabular}
 }
\end{table}

\begin{figure}[!htbp]
 \centering
 \subfloat[]{\includegraphics[width=0.45\columnwidth]{Figures/Elastic_E.pdf}}
 \subfloat[]{\includegraphics[width=0.45\columnwidth]{Figures/Elastic_Poisson.pdf}}
 \caption{\label{figure_linear_elastic} 
 Relationships  between  energy  storage  and  different  simulation  parameters
 (single element model): (a) Young's modulus; (b) Poisson's ratio.}
\end{figure}





Since  linear  elastic  material is used with  static  algorithm,
energy  components  related to dynamics (kinetic energy) and plasticity (plastic
free energy and plastic dissipation) are equal to zero. This means that all input
work is stored in the system, as observed in all cases.
%The match
%between  input work and the summation of stored energy and dissipation, which is
%zero  in  this  case,  verifies the reliability of the presented energy analysis
%methodology in single element cases.


%  HAN: is is inversely proportional to Poisson or directly proportional (as per
%       equation
%
% Han's comment:
% Should be proportional to Poisson's ratio. Corrected.

Figure~\ref{figure_linear_elastic} shows that   energy stored in the system is
inversely  proportional  to  Young's  moduli  and  proportional to  Poisson's ratio. This is
expected  because  of  the  following  equations for strain energy under
static shear loading:
%
\begin{equation}
E_S = \frac{1}{2} \tau \gamma = \frac{1}{2G} \tau^2 = \frac{1+\nu}{E} \tau^2
\label{equation_20}
\end{equation}

Note  that these relationships are only valid at constitutive level. For models with
more  finite elements,  stress  and strain are generally not uniform.  
The computation of energy depends on the distribution of energy density, and nonuniform stress/strain
distribution  will  result in nonuniform energy density distribution.


%  HAN:  need  to make a sketch/drawing of this model and of the simple block of
%        material from the previous example,  just so that readers can see what
%        is going on
% Han's comment:
% Added two figures at the beginning of the section



In  order  to  study  the  influence  of simulation parameters in larger models,
another  set  of simulations with cantilever model (Figure~\ref{figure_models}b) are performed.
Vertical loads
are  applied  to  the nodes of the free end. In this case, both shearing
and    bending  occurs, which means that in general a full 3D state of
stress and strain is present. The results are
presented      in    
Table~\ref{table_linear_elastic_cantilever} and Figure~\ref{figure_linear_elastic}. As
 expected,    energy    behavior   of cantilever is different than the
 single-element/constitutive example.


\begin{table}[htb]
 \centering
 \caption{Energy analysis results for linear elastic materials (cantilever model).}
 \label{table_linear_elastic_cantilever}
 \resizebox{\textwidth}{!}{
  \begin{tabular}{c|cc|cc|cccc|c}
   \hline
   \multirow{2}{*}{Case} & \multicolumn{2}{c|}{Material Property} & \multicolumn{7}{c}{Simulation Results}
   \\
   \cline{2-10}
   & $E$ (GPa) & $\nu$ & $u$ (m) & \boldmath$W_{Input}$ (J) & $E_{K}$ (J) & $E_{S}$ (J) & $E_{P}$ (J) & \boldmath$E_{Stored}$ (J) & \boldmath$D_P$ (J) \\
   \hline
   2a & 100 & 0.30 & 2.33E-3 & 116.57 & 0.00 & 116.57 & 0.00 & 116.57 & 0.00 \\
   2b & 150 & 0.30 & 1.55E-3 & 77.71 & 0.00 & 77.71 & 0.00 & 77.71 & 0.00 \\
   2 & 200 & 0.30 & 1.17E-3 & 58.28 & 0.00 & 58.28 & 0.00 & 58.28 & 0.00 \\
   2c & 250 & 0.30 & 9.33E-4 & 46.63 & 0.00 & 46.63 & 0.00 & 46.63 & 0.00 \\
   2d & 300 & 0.30 & 7.77E-4 & 38.86 & 0.00 & 38.86 & 0.00 & 38.86 & 0.00 \\
   \hline
   2e & 200 & 0.20 & 1.20E-5 & 65.89 & 0.00 & 65.89 & 0.00 & 65.89 & 0.00 \\
   2f & 200 & 0.25 & 1.26E-3 & 62.97 & 0.00 & 62.97 & 0.00 & 62.97 & 0.00 \\
   2 & 200 & 0.30 & 1.17E-3 & 58.28 & 0.00 & 58.28 & 0.00 & 58.28 & 0.00 \\
   2g & 200 & 0.35 & 1.02E-3 & 51.17 & 0.00 & 51.17 & 0.00 & 51.17 & 0.00 \\
   2h & 200 & 0.40 & 8.12E-4 & 40.60 & 0.00 & 40.60 & 0.00 & 40.60 & 0.00 \\
   \hline
  \end{tabular}
 }
\end{table}




For all cases, the energy balance between input and storage is maintained, which
gives us confidence on the energy calculation methodology for elastic material.
According to results in Figure~\ref{figure_linear_elastic},
energy  stored in the system is still inversely proportional to Young's modulus.
This is because the general equation for elastic strain energy density is:
%
\begin{equation}
E_S = 
\frac{1}{2E}
\left( \sigma^2_{xx} + \sigma^2_{yy} + \sigma^2_{zz} + 
2(1+\nu)(\sigma^2_{xy} + \sigma^2_{yz} + \sigma^2_{zx}) \right)
\label{equation_21}
\end{equation}

\noindent
So  as  long as all the elements have the same Young's modulus, the relationship
between  stored  energy  and  Young's  modulus  will  remain  valid.  
% 
%  Due to the
%  existence  of  bending  (normal  stresses),  the  curve of stored energy against
%  Poisson's  ratio  is  generally  nonlinear  for realistic models. Even more, the
%  increase  of  Poisson's  ratio  can  lead to a rise or fall of energy storage in
%  different  cases.  The  influence  of loading is more complicated in models with
%  multiple  elements.  Large  amount of tests have shown that in a cantilever-like
%  model,  the  plot of elastic strain energy stored in the system against vertical
%  force on the free end is usually close to a quadratic curve.
%
%    
%  
%  
%  From  two  shown examples, is is apparent that even for
%  linear  elastic  material, calculated energy can vary non-linearly.
%  




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Elastic-Perfectly Plastic Material}


%  By  performing  energy analysis for linear elastic materials, the input work and
%  storage  computations  are  verified.  The  next  step  is to apply the analysis
%  framework to elastic-plastic materials in order to analyze energy dissipation due
%  to plasticity.

In  this section, elastic-perfectly plastic material is used.
Equations~\ref{Eq8} and \ref{Eq12} indicate that in the case of no hardening the
rate of plastic free energy is zero. Then
the incremental plastic work is equal to incremental plastic dissipation. Note
that    this    is    one  of  the rare cases where plastic dissipation equals to plastic work.

Figure~\ref{figure_elastic_perfectly_plastic}  shows 
stress--strain  curve (left) and energy calculated for elastic-perfectly plastic
constitutive model (right) used here.


\begin{figure}[!htbp]
 \centering
 \subfloat[]{\includegraphics[width=0.49\columnwidth]{Figures/Perfectly_Plastic_Stress-Strain.pdf}}
 \subfloat[]{\includegraphics[width=0.49\columnwidth]{Figures/Perfectly_Plastic_Energy.pdf}}
 \caption{\label{figure_elastic_perfectly_plastic}
 Energy analysis of elastic perfectly plastic material: (a) Stress--strain curve;
 (b)  Input  work, plastic  dissipation, strain energy and plastic work.}
\end{figure}



In this case,  the plastic dissipation is equal to the
plastic  work. This means that the plastic free energy does not
develop  at  all  during  loading  and unloading. Zero plastic free energy
points   out   the absence of fabric evolution 
of  a particulate, elastic-plastic material, as all the input work is dissipated
through particle to particle friction.
% 
Since   there  is  no  plastic free energy $E_P$ in this case, the stored energy
equals to
mechanical  energy,  which  is  the combination of 
%
% HAN: always put variable names before their symbol!
%
% NAME THAM ALL!
%
strain energy $E_S$ and kinetic energy $E_K$. Total stored energy $E_{Stored}$
develops  nonlinearly  and  always  has the same value at the beginning of every
loop  after  the  first  one. Plastic dissipation $D_P$ increases linearly when the material yields.
This can be explained by rewriting
Equation~\ref{Eq8} with $\Psi_{pl}=0$:

\begin{equation}
\Phi = \sigma_{ij} \dot\epsilon^{pl}_{ij}
\end{equation}
%
where stress $\sigma_{ij}$ is constant after elastic perfectly plastic material yields,
and rate plastic deformation $\dot\epsilon^{pl}_{ij}$ is also constant. Then the rate
of plastic dissipation is constant which makes the plastic dissipation $D_P$ increase linearly.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Elastic-Plastic Material with Isotropic Hardening}



Next material model used is elastic-plastic with
linear isotropic hardening.  First  used to model monotonic behavior of elastic-plastic
materials,  isotropic  hardening assumes that the yield surface maintains shape, while isotropically 
(proportionally) changing its size. 
%
%  Failed to
%  model  the  Bauschinger effect by itself, isotropic hardening is usually used to
%  derive  more  sophisticated  elastic-plastic  materials or applied in combination
%  with  other  hardening rules. Either way, it is valuable to be analyzed from the
%  perspective of energy dissipation.
%  
%  
Figure~\ref{figure_isotropic_hardening}      illustrates    the    stress-strain
response  as well as energy balance for elastic-plastic material with
isotropic hardening.

\begin{figure}[!htbp]
 \centering
 \subfloat[]{\includegraphics[width=0.49\columnwidth]{Figures/IH_Stress-Strain.pdf}}
 \subfloat[]{\includegraphics[width=0.49\columnwidth]{Figures/IH_Energy.pdf}}
 \caption{\label{figure_isotropic_hardening}
 Energy  analysis  of  elastic-plastic  material  with  isotropic  hardening: (a)
 Stress-strain  curve;  (b) Input work, plastic dissipation, strain energy, and 
 plastic work.}
\end{figure}



As   can be observed from Figure~\ref{figure_isotropic_hardening}, plastic  free
energy  is equal to the plastic work, which means that the plastic dissipation is
zero  during  cycles of loading. Even though this might sound surprising, it can
be explained     using    basic  thermodynamics. Linear isotropic hardening,
used  in  this  case,  can be described  through a rate of the internal variable
(size of the yield surface) $\dot{k}$ as:
%
\begin{equation}
\dot k = \kappa_1 | \dot \epsilon^{pl}_{ij} |
\end{equation}
%
where  $|  \dot  \epsilon^{pl}_{ij}  |$  is the magnitude of the rate of plastic
strain  while  $\kappa_{1}$  is  a  hardening  constant.  Substituting  previous
equation into Equation~\ref{Eq12} yields:
%
% HAN: put names for variables, name \psi_{pl) and so on...
%
%
\begin{equation}
\psi_{pl} 
= 
\psi_{pl}^{iso} 
= 
\frac{\kappa_1}{2 \rho} \epsilon^{pl}_{ij} \epsilon^{pl}_{ij}
\end{equation}

Take the time derivative of the above equation:
%
\begin{equation}
\dot{\psi}_{pl} 
= 
\frac{\kappa_1}{\rho} \epsilon^{pl}_{ij} \dot{\epsilon}^{pl}_{ij}
\end{equation}

Then the rate of dissipation due to plasticity can be expressed as:
%
\begin{equation}
\Phi 
= 
\sigma_{ij} \epsilon^{pl}_{ij} - \rho \dot{\psi}_{pl}
=
(\sigma_{ij} - \kappa_1 \epsilon^{pl}_{ij}) \dot \epsilon^{pl}_{ij} 
=
(\sigma_{ij} - k m_{ij}) \dot \epsilon^{pl}_{ij}
\end{equation}
%
where   $m_{ij}$  is  the  plastic  flow direction. The plastic flow direction
defines the direction of incremental plastic strain, which can be different from the 
direction of total plastic strain. But in the case of von Mises type elastic-plastic material
with only isotropic hardening, whose yield surface is always a circle with center at the origin,
the plastic flow direction $m_{ij}$ is the same as the direction of the total plastic strain
$\epsilon^{pl}_{ij}$. Thus we have $\kappa_1 \epsilon^{pl}_{ij} = k m_{ij}$ in the above equation.

If we assume, for simplicity
sake, that plastic flow direction is associated with the yield function, that is
there is only deviatoric plastic flow, as yield function is of von Mises type,
the gradient of the yield surface $n_{ij} (=\partial F/ \partial \sigma_{ij})$ 
is equal to the
plastic    flow    direction    $m_{ij} (=n_{ij})$.    Noting   that   $\sigma_{ij}   \dot
\epsilon^{pl}_{ij}  =  s_{ij}  \dot  \epsilon^{pl}_{ij}$,  where $s_{ij} (=\sigma_{ij} - 1/3 \sigma_{kk} )$ is the
deviatoric  part  of  the  stress tensor, the rate of plastic dissipation can be
rewritten as:
%
\begin{equation}
\Phi = (s_{ij} - k n_{ij}) \dot \epsilon^{pl}_{ij} = \alpha_{ij} \dot \epsilon^{pl}_{ij}
\end{equation}
%

Realizing  that  the back stress $\alpha_{ij}$ is always zero since we assume no
kinematic  hardening,  then  the rate of plastic dissipation becomes zero, which
means  there  is no energy dissipation during cycles of loading for isotropically hardening material. 
Obviously, the observed  response  is  not physical from the perspective of energy dissipation.
%
%
%  Also,  this  happens not only during cyclic but monotonic loading as well, which
%  is  commonly  considered suitable for isotropic hardening. Therefore, one should
%  pay  extra  attention  to  energy  computation when applying isotropic hardening
%  rules.
% 
% Therefore,  isotropically hardening material models, can mimic/replicate material
%   stress-strain  behavior  for  monotonic  loading.  
%
% It  is  known that
% isotropic  hardening  material  models  cannot  be used to properly model cyclic
% loading response, as they are unable to model 
%
Therefore,  isotropic  hardening  material  models  cannot  properly  model energy
dissipation, even for monotonic loading.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Elastic-Plastic Material with Kinematic Hardening}



Compared  with  isotropic hardening, kinematic hardening can better describe the
constitutive,  stress-strain  behavior of elastic-plastic materials, particularly
for cyclic loading. 
%In order to model the Bauschinger effect, kinematic
%hardening  assumes  that  the  yield surface remains the same size and shape but
%translates  or rotates in stress space. 
Elastic-plastic material that relies on  kinematic
hardening is used to analyze  energy dissipation. 
Both linear and nonlinear kinematic hardening rules are
investigated in relation to energy dissipation.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Prager Linear Kinematic Hardening}


Prager's  linear  kinematic  hardening rule is given as:
%
\begin{equation}
\dot \alpha_{ij} = a_1 \dot{\epsilon}^{pl}_{ij}
\label{Eq16}
\end{equation}
% 
where $a_1$ is a hardening constant.
%
If  only  linear  kinematic hardening (Equation~\ref{Eq16}) is assumed, the back
stress  $\alpha_{ij}$  is  expressed  explicitly,  and  can  be substituted into
Equation~\ref{Eq12} yielding:
%
\begin{equation}
\psi_{pl} = \psi_{pl}^{kin} = \frac{a_1}{2 \rho} \epsilon^{pl}_{ij} \epsilon^{pl}_{ij}
\end{equation}

Take the time derivative of the above equation:
\begin{equation}
\dot{\psi}_{pl} = \frac{a_1}{\rho} \epsilon^{pl}_{ij} \dot{\epsilon}^{pl}_{ij}
\label{Eq32}
\end{equation}

If we again assume that the gradient of the yield surface $n_{ij}$ 
is equal to the plastic    flow    direction    $m_{ij}$,  as in the case of linear isotropic hardening,
then the rate of dissipation due to plasticity can be rewritten as:
%
\begin{equation}
\Phi
=
\sigma_{ij} \dot{\epsilon}^{pl}_{ij} - \rho \dot{\psi}_{pl}
=
(s_{ij} - \alpha_{ij}) \dot\epsilon^{pl}_{ij}
=
k m_{ij} \dot\epsilon^{pl}_{ij}
\label{Eq33}
\end{equation}
%

Notice that the term $m_{ij} \dot\epsilon^{pl}_{ij}$ denotes the
magnitude of the rate of plastic strain. Since only linear kinematic hardening is assumed,
the internal variable $k$ will remain constant.
%
% HAN, what do you mean if load are applied in a certain way
% 
% It was bad expression, I see you changed this already... This is what I intended to say...
%
So  if  loads are applied in such a  way
that  the rate of plastic strain is constant, then the
rate  of  dissipation  will  also  remain  constant. In other words, the accumulated
dissipation will be linearly increasing under the assumption of linear kinematic
hardening.




Figure~\ref{figure_linear_kinematic_hardening}  shows  stress--strain response (left) and   energy  computation
results (right) of  an  elastic-plastic  material  with  linear  kinematic hardening. 

\begin{figure}[!htbp]
 \centering
 \subfloat[]{\includegraphics[width=0.49\columnwidth]{Figures/Linear_KH_Stress-Strain.pdf}}
 \subfloat[]{\includegraphics[width=0.49\columnwidth]{Figures/Linear_KH_Energy.pdf}}
 \caption{\label{figure_linear_kinematic_hardening}
 Energy  analysis  of  elastic-plastic  material  with  linear kinematic
 hardening: (a) Stress--strain curve; (b) Input work, plastic dissipation 
 strain energy and plastic work.}
\end{figure}


%
% HAN: below, this 50% is somewhat ad hoc, that is there is nothing significant 
% about 50%, or any other number, as they are influenced by material parameters...


% Agree. Changing 50% into other expressions
%
As  expected,    the    plastic dissipation increases linearly once the material
yields.  In contrast to the isotropic hardening case, a significant amount of the input
work  is  dissipated due to material plasticity. The ratio of dissipated energy to input work is largely influenced by the material parameters. However, in general, energy dissipation will be
observed   if   kinematic   hardening  model is  used.  
% From  the  perspective  of
% thermodynamics,  this  is  physically  reasonable, which provides support to the
% point that kinematic hardening is more realistic than isotropic hardening.

Another  important observation is that the plastic work decreases during certain
phases of  reverse  loading,  while  the  actual rate of energy dissipation is
always  nonnegative.  It is important not to confuses the definitions of plastic work to plastic
dissipation,  as  plastic  work  can  increase  and decrease which can lead to a
(impossible) conclusion, violating the second law of thermodynamics.
% And more importantly, all the energy analysis will be invalid in such
%circumstances.  Unfortunately,  this  is  a  fairly  common oversight among many
%papers on seismic energy dissipation.





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Armstrong-Frederick Kinematic Hardening}



Armstrong-Frederick  kinematic hardening model \cite{Armstrong66} is often used
to simulate elastic-plastic  material behavior   under cyclic loading. 
%Its nonlinear feature makes the
%constitutive  models  better  describe the behaviors of realistic materials. 
%The
%material  parameters  are  usually  defined differently among various models. 
Material  parameters  of  the Armstrong Frederic kinematic hardening rule can be
derived from basic thermodynamics.
%model To
%clarify  these  parameters,  it  is beneficial to derive the Armstrong-Frederick
%kinematic  hardening  rule  from thermodynamics, which has been the base for all
%formulations in this paper.

\citet{Feigenbaum2007}   gave   the  sufficient  (but  not  necessary)
conditions  of  the  inequality in Equation~\ref{Eq8}. These conditions are more
restrictive  but satisfy the frameworks of most classical elastic-plastic models.
One   of  the  sufficient  conditions  corresponds  to  kinematic  hardening  is
expressed as:
%
%
% HAN: perhaps expand a bit this part as that equation (AF) pops out of blue skies!
%  In fact, you need to show von Mises yield surface F, all of then, isotropic and kinematic hardening...).
% this one will become F=(s_ij-alpha_ij)(s_ij-alpha_ij) - k ...

% I rewrote this part. The previous discussion focuses on how to derive hardening 
% equations from thermodynamics, but that's not our purpose. Now I just cite the equation
% from Feigenbaum and Dafalias (2007) to give the expression of AF hardening. Then
% add a little further derivation to point out that the energy dissipation in
% this case is non-negative and nonlinear.
%
\begin{equation}
\alpha_{ij} (\dot{\epsilon}^{pl}_{ij} - \frac{1}{a_1} \dot{\alpha}_{ij}) \ge 0
\label{Eq14}
\end{equation}
%
%
The following relation is a sufficient condition to satisfy Equation~\ref{Eq14}:
%
%
\begin{equation}
\dot{\alpha}_{ij} = a_1 \dot{\epsilon}^{pl}_{ij} - a_2 \dot{\lambda} \alpha_{ij}
\label{Eq15}
\end{equation}
%
where $\dot{\lambda}$  is a non-negative scalar plastic multiplier and 
$a_2$ is a non-negative material hardening constant. It can be proven that
$a_1/a_2$  is  related  to the limit of back stress magnitude $\left| \alpha_{ij
}\right|$. 
%Note that during the derivation of Equation~\ref{Eq14}, it is assumed
%to  be  associated  plasticity.  
Equation~\ref{Eq15} should be recognized as the
classical  Armstrong-Frederick  nonlinear kinematic hardening.

Taking the time derivative of the kinematic part of plastic free energy
(Equation~\ref{Eq32}), and substituting the expression of back stress $\alpha_{ij}$
(Equation~\ref{Eq15}) gives:
% 
\begin{equation}
\dot{\psi}_{pl}^{kin}
=
\frac{1}{\rho a_1} \alpha_{ij} \dot{\alpha}_{ij}
=
\frac{1}{\rho} \alpha_{ij} (\dot{\epsilon}^{pl}_{ij} - \frac{a_2}{a_1} \dot{\lambda} \alpha_{ij}) 
\end{equation}

If the gradient of the yield surface $n_{ij}$ is assumed to be
equal to the plastic    flow    direction    $m_{ij}$,  as was done in previous sections, then the rate of plastic energy dissipation of an Armstrong-Frederick kinematic hardening 
elastic-plastic material is given by:
% 
\begin{equation}
\Phi
=
\sigma_{ij} \dot{\epsilon}^{pl}_{ij} - \rho \dot{\psi}_{pl}
=
s_{ij} \dot\epsilon^{pl}_{ij} - \alpha_{ij} \dot\epsilon^{pl}_{ij} + \frac{a_2}{a_1} \dot{\lambda} \alpha_{ij} \alpha_{ij}
=
k m_{ij} \dot\epsilon^{pl}_{ij} + \frac{a_2}{a_1} \dot{\lambda} \alpha_{ij} \alpha_{ij}
\end{equation}

Compared with Equation~\ref{Eq33}, the above expression has an additional term 
which makes the rate of plastic dissipation non-constant even if the rate of 
plastic strain is constant. As the back stress $\alpha_{ij}$ becomes larger when
load increases, the rate of plastic dissipation also increases. This indicates a 
nonlinear result of total plastic dissipation, which is exactly what we have 
observed in our computations.

% In previously described general developments, until Equation~\ref{Eq12},  no plastic flow
% direction were assumed.  
%so the above
%expression  are  considered  general.  
% Assuming associated elastic-plasticity, plastic flow rule can be expressed as:
% plasticity, plastic  will be
%%considered,  which is typical with metal. The plastic flow rule going to be used
%is expressed as:
%
% \begin{equation}
% \dot{\epsilon}^{pl}_{ij} = \dot{\lambda} \frac{\partial f}{\partial \sigma_{ij}} = \dot{\lambda} m_{ij}
% \label{Eq13}
% \end{equation}
%
% where  $\dot{\lambda}$  is a non-negative scalar plastic multiplier and $m_{ij}$
% denotes plastic flow direction.


% Substituting Equation~\ref{Eq12} and \ref{Eq13} into \ref{Eq8}, one rewrites the
% inequality  in  Equation~\ref{Eq8}  with  decoupled  components  correspond  to
% different       hardening      mechanisms.      
% of the evanescence
%memory type.


\begin{figure}[!htbp]
 \centering
 \subfloat[]{\includegraphics[width=0.49\columnwidth]{Figures/AF_KH_Stress-Strain.pdf}}
 \subfloat[]{\includegraphics[width=0.49\columnwidth]{Figures/AF_KH_Energy.pdf}}
 \caption{\label{figure_AF_kinematic_hardening}
 Energy  analysis  of elastic-plastic material with Armstrong-Frederick kinematic
 hardening:  (a)  Stress-strain  diagram;  (b)  Input work, plastic dissipation,
 strain energy, and plastic work.}
\end{figure}



Figure~\ref{figure_AF_kinematic_hardening} shows the energy computation results
of  an  elastic-plastic  material  with  Armstrong-Frederick kinematic hardening.
Compared    to   all previous cases, the material response of this model is more
sophisticated and  more  realistic.  
%Up  to  60\%  of  the input energy is
%dissipated  due  to plasticity in each cycle. 
Decrease of plastic work is
observed, again,  while  the  plastic  dissipation is always nonnegative during the entire
simulation. For both linear and nonlinear kinematic hardening cases, the plastic
free  energy  is relatively small compared to the plastic dissipation. 
%It can be
%easily  ignored  in  SSI  simulations under seismic loading, since the number of
%load cycles  is  large  and  the  response is wrong. Nevertheless, the
%plastic  free  energy  should never be neglected, in the sense of satisfying the
%basic laws of thermodynamics.





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusions}


%OVDE


Presented  was  a methodology for (correct) computation of energy dissipation in
elastic-plastic materials based on the second law of thermodynamics.  
% 
A  very important role of plastic free energy was analyzed, with highlights  on 
its  physical  nature and theoretical formulations. The proposed methodology has
been illustrated using a number of elasto-plastic material models.

An analysis of a common misconception that equates plastic work and dissipation,
which leads to the
violation  of  the  basic  principles  of  thermodynamics,  was addressed. 
%The
%introduction  of  plastic  free  energy  has  been  proved  to  be  effective on
%eliminating  such  violation. 
%
A  conceptual  example,  for    granular  materials,  was used to  explain  the 
physical  meaning  of  plastic  free energy. It was also shown that plastic free
energy is responsible  for the evolution of internal variables.

It  was shown that energy  balance  is ensured  by  taking into consideration all
 energy  components,  including    kinetic  and  strain  energy.  Input work was
 balanced  with  the stored and dissipated energy, expressed as the summation of
 all possible components. 


Presented  approach  was illustrated and  tested  using several
elastic-plastic  constitutive  models  with  various  hardening rules.
%
Elastic  materials  showed no  energy  dissipation (as expected), leading to the
input work being equal to the stored energy. 
%
Elastic-perfectly plastic materials had no change in
plastic  free  energy,  which  led  to the equality of plastic work and plastic
dissipation  and  indicated  no  evolution  of  particle arrangements. The
plastic  dissipation,  in  that  case,  was  observed  to  be increasing linearly.
Isotropic hardening materials experienced zero dissipation even after
yielding.  This observation was surprising, but verified by further derivation
of  energy  equations.  This  observation  also  serves  as  a reminder that the
isotropic hardening rules  can   be  used, but only with observed lack of energy
dissipation.   Prager's   linear  and  Armstrong-Frederick  nonlinear  kinematic
hardening  materials  both gave significant dissipations, with large fluctuation
of  plastic  free  energy  as well. In the case with linear kinematic hardening,
linear  increase  of  dissipation  was  derived  and  observed, while energy was
dissipated   nonlinearly  in the case of nonlinear kinematic hardening. Although
the plastic free energy was 
not  significant  for  some  materials,  it is noted that it  should  always  be
recognized and
considered   during   energy   analysis,   so   that  the  basic  principles  of
thermodynamics are maintained.







%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Acknowledgments}



This    work  was  supported  by  the  US-DOE.  We would like to thank Professor
Yannis  Dafalias  (UCD  and  NTUA)  for  inspiring  discussions.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


\bibliography{refmech}

\bibliographystyle{abbrvnat}


\end{document}