\documentclass{cnmauth}

%\usepackage[dvips,
%            colorlinks,
%            bookmarksopen,
%            bookmarksnumbered,
%            citecolor=blue,
%            urlcolor=blue]{hyperref}
%\usepackage{endfloat}
\usepackage{graphics}
\usepackage{graphicx}
%\usepackage{fancyheadings,psfig}
\usepackage{amsmath}
\usepackage{mathrsfs}
\usepackage{amsfonts}
%\usepackage{natbib}
%\bibpunct{(}{)}{;}{a}{,}{,}

%% For derivatives , must have {amsmath}
\newcommand{\ud}{\mathrm{d}}

\usepackage{remark}
\newremark{remark}{Remark}[section]



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

\def \baselinestretch{2}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{document}
\CNM{1}{6}{00}{28}{00}

\runningheads{Jeremi{\' c} and Cheng}
{On Large Deformation Hyperelasto--Plasticity of Anisotropic Materials}

\title{On Large Deformation Hyperelasto--Plasticity of Anisotropic Materials}



\author{ 
Boris Jeremi{\'c} 
(\raisebox{-1.0mm}{\includegraphics[width=2.7cm]{/home/jeremic/tex/works/IMENAuOriginaly/BorisJeremic.eps}})
\affilnum{1}
\comma\corrauth \\
%Boris Jeremi{\'c} ({\cyr Boris Jeremi\cj{}})$^{1}$ \\
%%%%%%%%
Zhao Cheng
(\raisebox{-1.2mm}{\includegraphics[width=1.1cm]{/home/jeremic/tex/works/IMENAuOriginaly/Zhao.eps}})
\affilnum{2} \\
%%%%%%%%
%Yannis Dafalias$^{4}$ \\
}

%\author{ 
%Boris Jeremi{\'c}\affilnum{1}\comma\corrauth
%%
%\;\;
%Zhao Cheng\affilnum{2} 
%}

\address{
\affilnum{1}\ {Associate Professor, 
Department  of  Civil  and  Environmental  Engineering,  University  of
California, Davis, CA 95616. \texttt{jeremic@ucdavis.edu}} \\%
\affilnum{2}\ {Staff Engineer, Earth Mechanics, Inc., Oakland, CA 94621}\\
%
}


%
%\author{Boris Jeremi{\'c} and Zhao Cheng}
%
%\address{Department  of  Civil  and  Environmental  Engineering,  University  of
%California, Davis, CA, 95616, U.S.A.}

\corraddr{Boris Jeremi{\'c}, Department    of  Civil  and  Environmental  Engineering,  One Shields
Ave., University of California, Davis, CA, 95616, U.S.A. \texttt{jeremic@ucdavis.edu}}


\cgsn{NSF}{EEC--9701568, CMS--0337811}

\noreceived{}
\norevised{}
\noaccepted{}

\keywords{anisotropy, hyperelastic-plastic, large deformations}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 
%  
%  \newpage
%  \tableofcontents
%  \newpage
%  
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{abstract}
% Text of abstract

%% Tentative Abstrcat (Fr MIT3Paper), need modification
The    main    goal  of  this  short paper is to bring forward and explore the  
issue  of  constitutive integrations in large deformation regime for anisotropic
hyperelasto-plastic  material  models. Of particular interest is the appropriate
treatment  of  cases  with non--coaxiality of principal stress and its conjugate
hyperelastic strain.   The   non--coaxiality   can  result  from  hyperelastic
anisotropy,  anisotropic  flow  directions  or  anisotropic  yield  function. In
addition  to  discussing  issues  related  to  non--coaxiality,  a  constitutive
integration  algorithm  is  presented  in  some  detail  that  is designed to be
applicable   to   anisotropic   non--coaxial   (as   well   as   to   isotropic)
hyperelastic--plastic solids.
A numerical example is used to illustrate the problem.

\end{abstract}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   %  
%   %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%   \begin{keyword}
%   % keywords here, in the form: keyword \sep keyword
%   Hyperelastoplasticity \sep Large Deformation \sep Return Mapping \sep Anisotropy 
%   % PACS codes here, in the form: \PACS code \sep code
%   
%   \end{keyword}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
\label{Introduction_HEPLD}

The    hyperelasto--plastic    format  mostly  used  today  in large deformation
computations  is based on work by Simo and Ortiz \cite{local-12}, \cite{Simo92},
Bathe  et  al. \cite{Bathe85}, Simo \cite{Simo88a}, \cite{Simo88b}, Eterovic and
Bathe  \cite{Eterovic90},  Peri{\'c} et al. \cite{Peric92} and Cuitino and Ortiz
\cite{Cuitino92}.    Earlier  work  on      multiplicative      split   of  the 
deformation gradient (Hill \cite{Hill50}, Bilby et al. \cite{Bilby57} Kr{\"o}ner
\cite{Kroner60},   Lee   and   Liu   \cite{Lee67},   Fox  \cite{Fox68}  and  Lee
\cite{Lee69}) is incorporated in all the previous developments.


However,   few   researchers  have  addressed  the  issue  of  large-deformation
hyperelastoplastic  computational  formulations  for  anisotropic  materials. We
mention  an  algorithm by Eterovic and Bathe \cite{Eterovic90} which is based on
additive   split   of  logarithmic  stress  and  strain  measures  (elastic  and
hyperelastoplastic).  They  have  also  explored  the  use  of  a  more  general
approximation  of  deformation  tensors  that  can  support  both  isotropic and
anisotropic    material    models.    More   recently,   Papadopoulos   and   Lu
\cite{Papadopoulos98}  developed  a  general  framework  for  finite deformation
elasto-plasticity  that  is  based  on  the  early  work  of  Green  and  Naghdi
\cite{Green1965}.   Developments  include  provisions  for  non-collinearity  of
principal  stress  and strain measures. The framework was tested using Von-Mises
type  yield  criteria  with  translational  kinematic  hardening,  which  retain
isotropy of yield and plastic potential functions. Somewhat similar developments
were  reported by Miehe et al. \cite{Miehe2002} and \cite{Miehe2004} . They used
initial anisotropy in hyperelastic models and initially anisotropic criterion by
Hill  \cite{Hill50}  to  successfully simulate various problems. However, it was
not  clear if, and how, the non--coaxiality of principal stress and its conjugate
hyperelastic strain measures evolve and how much it influences the results.



%\paragraph{Motivation.}
Addressed  here  is    the  issue  of  non--coaxiality of principal stresses and
its conjugate hyperelastic  strains    and    its    influence    on    large    deformation
hyperelastic--plastic integration algorithms. This non--coaxiality can stem from
a  number  of features present in general hyperelastic--plastic material models.
We list some of those features below:
%
%
%%\paragraph{Generality of Material Models}
%
%Aspects of  {\it  generality}  of  a described set of models, as oposed to isotropic ones
%(Simo  \cite{Simo92},  Simo and Miehe \cite{Simo92b}), are listed below:
%to explain more:
%  
\begin{description}
  
  \item[Anisotropic Hyperelasticity:] 
  The non--coaxiality  of  principal  stresses  and  its conjugate hyperelastic
  strains  can  result  if  the  hyperelastic  strain  energy function is not an
  isotropic  function  of  the elastic Euler-Green strain $\bar{b}^{e}_{ij}$ (or
  the  elastic  Lagrange-Green  strain $\bar{C}^{e}_{ij}$). In the case that the
  yield  function  is an isotropic one as well, the assumption can be made that the
  elastic  Euler--Green  strain  $\bar{b}^{e}_{ij}$  and  the  Kirchhoff  stress
  $\tau_{ij}$  commute.  This isotropic assumption makes basis for many previous
  developments  (for example Simo \cite{Simo92}, Simo and Miehe \cite{Simo92b}).
  However,  use of this assumption in developing consistent algorithms for large
  deformation  hyperelastic--plastic  response precludes them from being applied
  to  anisotropic  hyperelastic--plastic models. Our presented algorithm removes
  such limitation.
  
  
  \item[Non-isotropic Flow Rules:] 
  The  non--coaxiality  can  result  from  the  plastic  flow  rule that is not 
  isotropic function of the Kirchhoff stress tensor $\tau_{ij}$.
  
%  
%  In  general,  the flow rule of hyperelastic--plastic material models mostly presented
%  in literature is associated with the yield
%  criterion    (associated    plasticity).  The flow rule in those models can be
%  extended  to  be  associated  with  a  different (then yield function) plastic
%  potential function
%  flow rule but only on the
%  condition  that the flow rule is an isotropic function of the Kirchhoff stress
%  tensor  $\tau_{ij}$. Simo's model can not be extended to an non-isotropic flow
%  rule.

  \item[Kinematic Hardening:]
  The non--coaxiality can develop if kinematic  hardening  is  present.
  In that
  case, the induced anisotropy may destroy the isotropy of the yield function.
  One  such  simple  example  is  that  of von Mises model (yield function) with
  an initial back stress that  is non--coaxial with the developing conjugate strain. 
  Prevost \cite{Prevost87a} noted this case in which the existence of kinematic hardening 
  destroys the isotropy of the yield function.
%
  On the other hand, algorithms that assume that the yield function is an isotropic function of 
  the  Kirchhoff  stress  tensor  $\tau_{ij}$  in  order  to  satisfy  the frame
  indifference cannot deal with this type of anisotropy.
 
 \item[Reversal of Loading, Cyclic loading:] 
 The        non--coaxial      may  also  occur  during reversal of loading. That
 reversal of loading might
 change    the   ordering  of  three  principal  directions  for  stress and its
 conjugate hyperelastic  strain
 independently, thus destroying coaxiality.
%
% If the hyperelasticity, yielding, and flowing are isotropic, 
% the principal directions of the elastic Eulerian-Green strain $\bar{b}^{e}_{ij}$, 
% the Kirchhoff stress $\tau_{ij}$ and the principal flow directions change with the same manner. 
% Simo's model may be used for unloading and reloading if anisotropic hyperelasticity, non-isotropic flow rule, 
% and kinematic hardening are circumvented, but it should be careful to figure out 
% the orders the principal directions during implementation.                    
%  
  \end{description}
  

%   The commutes of 
%   the elastic Eulerian-Green strain $\bar{b}^{e}_{ij}$, 
%   the Kirchhoff stress $\tau_{ij}$ and 
%   the normal tensor of the yield function $\dfrac{\partial \mathscr{F}}{\partial \tau_{ij}}$
%   in Simo's model (associated flow rules are assumed)  
%   were used as an important assumption in the implementation of its implicit algorithms. 
%   Here the model based on the state variables in the intermediate configuration has no restriction to 
%   hyperelastic and yielding isotropy restriction. 
%   The implicit algorithm for this model is different from that of Simo's model.       
%        


In   what   follows,   a   brief   description   of   a   generic version    of
hyperelastic--plastic  material model is presented. Following that, the implicit
algorithm  is  developed  that  resolves  all  the  above  mentioned issues with
anisotropy.  In  addition to the fully implicit algorithm, presented is also the
consistent  tangent  operator  that  is  used in conjunction with global (finite
element)  level  Newton  iterations.  A simple  example which features a von Mises
material model with initial back--stress is used to illustrate our approach.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{General Anisotropic Hyperelasto--Plastic Models}
\label{GHPM}

In  this  section,  an  overview is made of a generic material model formulation
that incorporates anisotropic  feature described in the introduction.
%
Material      models    incorporating   described   features  are  mostly  found
in  dealing with continuum representation of particular materials (soils,
rock, concrete, powders, bone material, foams...). 
%
A  quite  general  description  of any incremental material model can be done by
separating the model into its components, namely (a) hyperelastic relations, (b)
yield function, (c) flow rules and (d) hardening/softening laws.

% 
%   It  is worth noting that although we rely on the multiplicative decomposition of the
%   deformation  gradiente  ($F_{ij}=F^e_{ik}F^p_{kj}$),  deformation  rate  will be
%   decomposed additively (in intermediate configuration) as
%   %%%
%   \begin{equation}
%   \bar{D}_{ij}
%   = \bar{D}^e_{ij} + \bar{D}^p_{ij}
%   = \dot{\bar{E}}^e_{ij} + \bar{D}^p_{ij}
%   \label{VelocityGradientDecom}
%   \end{equation}
%   %%%
%   which   is   resembles   the   additive  decomposition  from  the  infinitesimal
%   counterpart and will be used in subsequent developments.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Hyperelastic Relations}

The  hyperelastic  material  derives from a general function $\Psi$ which can be
anisotropic, and relates the second Piola-Kirchhoff stress $\bar{S}_{ij}$ 
in the intermediate configuration and its  conjugate
elastic Green strain $\bar{E}^{e}_{ij}$:
%%
\begin{equation}
\dot{\bar{S}}_{ij}
= \bar{\mathcal L}_{ijkl} \dot{\bar{E}}^{e}_{kl}
%\label{hyperDot}
%\end{equation}
%
\;\;\;\;\;  \mbox{;} \;\;\;\;\;
%
%\begin{equation}
\bar{\mathcal L}_{ijkl} = \dfrac{\partial \bar{S}_{ij}} { \partial \bar{E}^{e}_{ij}}
= \dfrac{\partial^{2} \Psi} { \partial \bar{E}^{e}_{ij} \partial \bar{E}^{e}_{kl}}
\label{hyperTangent}
\end{equation}
%
where  $\bar{\mathcal  L}_{ijkl}$  is the Lagrange-Green stiffness tensor in the
intermediate   configuration,   and   $\bar{E}_{ij}^{e}   :=  {1}/{2}(
\bar{C}_{ij}^{e}  -  \delta_{ij})$ is the elastic Lagrange--Green strain tensor,
while $\bar{C}_{ij}^{e}  := F_{ki}^{e} F_{kj}^{e}$ is the elastic Lagrange deformation
tensor.

It proves  convenient to define this relationship (Eq. \ref{hyperTangent}) using
Mandel   stress      $\bar{T}_{ij}   (=   \bar{C}_{ik}^e   \bar{S}_{kj}$)
and the elastic Green strain $\bar{E}_{ij}$,  both in the intermediate configuration
% %%
% %\begin{equation}
% %\bar{T}_{ij} = \bar{C}_{ik}^{e} \bar{S}_{kj} 
% %\label{MStress}
% %\end{equation}
% %%
% 
% \begin{equation}
$
\dot{\bar{T}}_{ij} = \bar{\mathcal L}_{ijmn}^{M}\dot{\bar{E}}^{e}_{mn}
$
% \label{MandelTime}
% \end{equation}
%
where
$
\bar{\mathcal L}_{ijmn}^{M} = \delta_{im}\bar{S}_{jn} + \delta_{in}\bar{S}_{jm} + \bar{C}_{ik}^{e}\bar{\mathcal L}_{kjmn} 
$  is  the  stiffness tensor connecting Mandel stress and Lagrange--Green strain
tensors.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Yield Conditions}

The yield function can be expressed in terms of 
the Mandel stress in the intermediate configuration $\bar{T}_{ij}$, 
the scalar stress-like internal variable $q$,  
and the kinematic stress-like internal variable $a_{ij}$ as
%%
\begin{equation}
\mathscr{F} = \mathscr{F}(\bar{T}_{ij}, q, a_{ij}) = 0
\label{yieldFun}
\end{equation}
%%
%
%\begin {remark}
Since  $\bar{T}_{ij}$  is  invariant  under  any  rigid  translations or rotations,  this  yield
function can represent any type of anisotropy.
%
%are not a restricted to isotropy . The anisotropy can be incorporated into this model.
%\end {remark}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Flow Rules}

% The principle of maximum dissipation can result in an associated flow rules. 
% Similar to those in the small deformation regime, 
% the associated flow rules may be extended a small step to general flow rules, 
% but make the associated flow rules as the special cases.   

The  non--associated  (or associated), anisotropic flow rule is defined in terms
of the plastic deformation gradient as
%
\begin{equation}
\bar{L}_{ij}^p 
= \dot{F}^{p}_{ik}(F^{p})_{kj}^{-1} 
= \dot{\lambda}\bar{M}_{ij} (\bar{T}_{ij}, q, a_{ij})
\label{LFlowRule}                       
\end{equation}
%%
where $\bar{M}_{ij}$ is the direction of plastic deformation.
%
By using the definition for the plastic deformation rate  
($\bar{l}_{ij}^{p} = \bar{C}_{ik}^{e}\bar{L}_{kj}^{p}$)
and  the symmetric part of the plastic deformation rate 
($\bar{d}_{ij}^{p}  =  Sym( \bar{C}_{ik}^{e}\bar{L}_{kj}^{p} ) $),
we obtain
%
\begin{equation}
\bar{d}^p_{ij} 
= \dot{\lambda} \bar{M}^{C}_{ij} 
%= %= \dfrac{1}{2} \dot{\lambda} \left( \bar{C}^{e}_{ik}\bar{M}_{kj} + \bar{M}_{jk}\bar{C}^{e}_{ki} \right)
%= \label{DFlowRule}
%= \end{equation}
%= %%
%
\;\;\;\;\;  \mbox{;} \;\;\;\;\;
%
%= where
%= %%
%= \begin{equation}
\bar{M}^{C}_{ij} 
= Sym ( \bar{C}^{e}_{ik}\bar{M}_{kj} )
\label{barMC}
\end{equation}
%%
%
For scalar ($\xi$) and kinematic ($\eta_{ij}$)
internal variables, flow rules can be written as,
%%
\begin{equation}
\dot{\xi} = -\dot{\lambda} \; n^{scalar} (\bar{T}_{ij}, q, a_{ij})
%
\;\;\;\;\;  \mbox{;} \;\;\;\;\;
%
\dot{\eta}_{ij} = -\dot{\lambda} \; n^{kinematic}_{ij} (\bar{T}_{ij}, q, a_{ij})
\label{FlowRuleIK}
\end{equation}

%%
%   
%   \begin {remark}
%    For associated flow rules, 
%    %%
%    \begin{equation}
%    \bar{M}_{ij} = \dfrac{\partial \mathscr{F} (\bar{T}_{ij}, q, a_{ij}) } {\partial \bar{T}_{ij}}
%    %%\label{AssociatedM}
%    %%\end{equation}
%    %%
%    ,\;\;\;
%    %%\begin{equation}
%    n^{i} = \dfrac{\partial \mathscr{F} (\bar{T}_{ij}, q, a_{ij}) } {\partial q}
%    %%\label{AssociatedNI}
%    %%\end{equation}
%    %%
%    ,\;\;\;
%    %%\begin{equation}
%    n^{k}_{ij} = \dfrac{\partial \mathscr{F} (\bar{T}_{ij}, q, a_{ij}) } {\partial a_{ij}}
%    \label{AssociatedNK} 
%   \end{equation}
%   %%
%   \end {remark}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\paragraph{Hardening Laws}

The scalar and kinematic hardening laws are defined as
%%
\begin{equation}
\dot{q} = K \dot{\xi}
%
\;\;\;\;\;  \mbox{;} \;\;\;\;\;
%
\dot{a}_{ij} = {\mathcal H}_{ijkl} \dot{{\eta}}_{kl}
\label{HdLawIK}
\end{equation}
%%
where hardening/softening function ${\mathcal H}_{ijkl}$  can be any anisotropic
function while $K$ remains isotropic function.
%
% When the isotropic hardening stiffness $K$ 
% and 
% the kinematic hardening stiffness ${\mathcal H}_{ijkl}$ are constants, 
% one obtains linear isotropic or kinematic hardening law.
% The hardening stiffnesses $K$ and ${\mathcal H}_{ijkl}$ should not be considered only as constants, for example, when
% \begin{equation}
% K(\xi) = q_{sat} \beta \exp (- \beta \xi)
% \label{SaturatedIHd}
% \end{equation}
% where $q_{sat}$ and $q_0$ is the saturated stress-like isotropic internal variables,
% $\beta$ is a parameter to give how fast $K$ changes, one obtain a non-linear isotropic hardening 
% law with a saturated value. 
% %%
% The second in Equation \ref{HdLawIK} may be simplified as 
% $\dot{a}_{ij} = H \dot{{\eta}}_{ij}$,    
% where $H$ is a scalar constant, 
% for this case ${\mathcal H}_{ijkl}$ in Equation \ref{HdLawIK} actually equals to $H\delta_{ik}\delta_{jl}$.
%
%  For simplicity, the is is assumed that the internal variables $q$  and
%  $a_ij}$ are uncoupled, that is
%  
%    is here assumed between the internal variables, 
%  or the stress-like internal variables are only depend on the corresponding conjugate strain-like internal variables.   
%  It is possible that these internal variables have coupling interaction for more complicated models. 
%  

%===%===========================================================================================================
%===\paragraph{Generality of Material Models}
%===
%===Aspects of  {\it  generality}  of  a described set of models, as oposed to isotropic ones
%===(Simo  \cite{Simo92},  Simo and Miehe \cite{Simo92b}), are listed below:
%===to explain more:
%===  
%===\begin{description}
%===  
%===  \item[Anisotropic Hyperelasticity:] 
%===  Material models restricted to isotropy make an assumption that 
%===  the elastic Eulerian-Green strain $\bar{b}^{e}_{ij}$ and 
%===  description    the     Kirchhoff  stress  $\tau_{ij}$  commute.  This  is valid 
%===  only if the   hyperelastic strain energy function is an isotropic function of
%===  the elastic Eulerian-Green strain $\bar{b}^{e}_{ij}$ 
%===  (or    the    elastic  Lagrangian-Green strain $\bar{C}^{e}_{ij}$). In general
%===  material modeling (for anisotropic hyperelasticity), no   such assumption can be made.
%===  
%===  \item[Non-isotropic Flow Rules:] 
%===  The  restriction  of isotropic models is also present if the plastic flow rule
%===  is an isotropic function of the Kirchhoff stress tensor $\tau_{ij}$. 
%===  Although  generally the flow rule of Simo's model is associated with the yield
%===  criterion,  it  may be extended to an non-associated flow rule but only on the
%===  condition  that the flow rule is an isotropic function of the Kirchhoff stress
%===  tensor  $\tau_{ij}$. Simo's model can not be extended to an non-isotropic flow
%===  rule.
%===  
%===  \item[Kinematic Hardening:]
%===  Simo's model also assumes that the yield function is an isotropic function of 
%===  the Kirchhoff stress tensor $\tau_{ij}$ to make sure the frame indifference. 
%===  If there exists kinematic hardening, the induced anisotropy may destroy the isotropy of the yield function. 
%===  One simple such example is that Von-Mises model with kinematic hardening, 
%===  but the initial back stress is non-
%===  coaxial to the developing conjugate strain. 
%===  Prevost \cite{Prevost87a} also commented on how the existence of kinematic hardening 
%===  may violates the isotropy of the yield function.
%===  Simo's model should therefore not be applied to all cases of the kinematic hardening.
%=== 
%=== %%make it out, since it is somewhat confusing. 
%=== %out \item[Unloading and Reloading, Cyclic loading:] 
%=== %out Basically, unloading and reloading will change the order of the three principal directions. 
%=== %out If the hyperelasticity, yielding, and flowing are isotropic, 
%=== %out the principal directions of the elastic Eulerian-Green strain $\bar{b}^{e}_{ij}$, 
%=== %out the Kirchhoff stress $\tau_{ij}$ and the principal flow directions change with the same manner. 
%=== %out Simo's model may be used for unloading and reloading if anisotropic hyperelasticity, non-isotropic flow rule, 
%=== %out and kinematic hardening are circumvented, but it should be careful to figure out 
%=== %out the orders the principal directions during implementation.                    
%===  
%===  \end{description}
%===  
%===  The commutes of 
%===  the elastic Eulerian-Green strain $\bar{b}^{e}_{ij}$, 
%===  the Kirchhoff stress $\tau_{ij}$ and 
%===  the normal tensor of the yield function $\dfrac{\partial \mathscr{F}}{\partial \tau_{ij}}$
%===  in Simo's model (associated flow rules are assumed)  
%===  were used as an important assumption in the implementation of its implicit algorithms. 
%===  Here the model based on the state variables in the intermediate configuration has no restriction to 
%===  hyperelastic and yielding isotropy restriction. 
%===  The implicit algorithm for this model is different from that of Simo's model.       
%===       



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Algorithmic Constitutive Formulation}




% \paragraph{Constitutive Implicit Integration Algorithm}

In this section, a constitutive integration algorithm is developed that does not
suffer  from  restrictions  to  isotropic  material  models (as described in the
introductory section).
%
It  should  be noted that in developments that follow, the  Karush--Kuhn--Tucker
(KKT) 
conditions  $\dot{\lambda}  \geq 0$, $\mathscr{F}(\bar{T}_{ij}, q, a_{ij}) \leq
0$ and $\dot{\lambda} \mathscr{F} (\bar{T}_{ij}, q, a_{ij}) = 0$ 
(Karush  \cite{Karush1939},  Kuhn  and  Tucker \cite{Kuhn-Tucker50})
which should be satisfied simultaneously, are used.
%%   %%
%%   \begin{equation}
%%   \dot{\lambda} \geq 0 
%%   %
%%   \;\;\;\;\;  \mbox{;} \;\;\;\;\;
%%   %
%%   \mathscr{F}(\bar{T}_{ij}, q, a_{ij}) \leq 0 
%%   %
%%   \;\;\;\;\;  \mbox{;} \;\;\;\;\;
%%   %
%%   \dot{\lambda} \mathscr{F} (\bar{T}_{ij}, q, a_{ij}) = 0
%%   \end{equation}
%%
In     addition     to     KKT     condition,          consistency     condition
$\dot{\mathscr{F}}(\bar{T}_{ij}, q, a_{ij}) = 0 $ is also enforced.
%%  %%
%%  \begin{equation}
%%  \dot{\mathscr{F}}(\bar{T}_{ij}, q, a_{ij}) = 0
%%  \label{consistCond}
%%  \end{equation}



We start by integrating the flow rule (Equation (\ref{LFlowRule})) from 
time step $t$ to $t+1$ as
%
\begin{equation}
{}^{n+1} \! F_{ij}^{p}
=
\exp{\left( \Delta \lambda {}^{n+1} \! \bar{M}_{ik}   \right)}
{}^{n} \! F_{kj}^{p}
\label{IIA040}
\end{equation}
%
By using the multiplicative decomposition of the deformation gradient 
and Equation (\ref{IIA040}) we obtain
%
\begin{equation}
{}^{n+1} \! F_{ij}^{e}
=
{}^{n+1} \! F_{im} \; \left( {}^{n} \! F_{mk}^{p} \right)^{-1}
\exp{\left( - \Delta \lambda {}^{n+1} \! \bar{M}_{kj}   \right)}
%\nonumber \\
=
{}^{n+1} \! \bar{F}_{ik}^{e,trial}
\exp{\left( - \Delta \lambda {}^{n+1} \! \bar{M}_{kj}   \right)}
\label{IIA060}
\end{equation}
%
where it was used that 
% %
% \begin{equation}
$
{}^{n+1} \! F_{ik}^{e,trial}
=
{}^{n+1} \! F_{im} \; \left( {}^{n} \! F_{km}^{p} \right)^{-1}
$.
% \label{IIA070}
% \end{equation}
%
%
\noindent
The elastic deformation is then
%
\begin{eqnarray}
{}^{n+1} \! \bar{C}_{ij}^{e}
%
&=&
%
\left( {}^{n+1} \! F_{im}^{e} \right)^{T} \;
{}^{n+1} \! \bar{F}_{mj}^{e}
\nonumber \\
%out&=&
%out\exp{\left( - \Delta \lambda {}^{n+1} \! \bar{M}_{ir}^{T} \right)}  \;
%out\left( {}^{n+1} \! F_{rk}^{e,trial} \right)^{T} \;
%out%
%out{}^{n+1} \! F_{kl}^{e,trial} \;
%out\exp{\left( - \Delta \lambda {}^{n+1} \! \bar{M}_{lj}   \right)}
%out%
%out\nonumber \\
&=&
\exp{\left( -\Delta \lambda {}^{n+1} \! \bar{M}_{ir}^{T} \right)} \;
%
{}^{n+1} \! \bar{C}_{rl}^{e,trial} \;
%
\exp{\left( -\Delta \lambda {}^{n+1} \! \bar{M}_{lj}   \right)}
%
\label{IIA080}
\end{eqnarray}
%
By  recognizing  that the exponent  of  a  tensor  can  be  expanded  in Taylor
series (Pearson \cite{Pearson74}),
%
\begin{equation}
\exp{\left( -\Delta \lambda {}^{n+1} \! \bar{M}_{lj} \right)}
=
\delta_{lj}
-
\Delta \lambda {}^{n+1} \! \bar{M}_{lj}
+
Sym (  \Delta \lambda {}^{n+1} \! \bar{M}_{ls} )
+
\cdots
\label{IIA090}
\end{equation}
%
and  by applying first order expansion to Equation (\ref{IIA080}), 
and    by  neglecting    the      higher    order    terms  (Jeremi{\'c}  et al.
\cite{Jeremic97f}) the
solution  for  the right elastic deformation tensor ${}^{n+1}\!\bar{C}_{ij}^{e}$
can be written as
%
\begin{equation}
{}^{n+1} \! \bar{C}_{ij}^{e}
%
=
{}^{n+1} \! \bar{C}_{ij}^{e,trial}
-
\Delta \lambda \; \left(
{}^{n+1} \! \bar{C}_{ik}^{e,trial} \; {}^{n+1} \! \bar{M}_{kj}  \;
+
{}^{n+1} \! \bar{M}_{ik} \; {}^{n+1} \! \bar{C}_{kj}^{e,trial}
\right)
\label{IIA130}
\end{equation}

% 
% one obtains
% %
% \begin{eqnarray}
% {}^{n+1} \! \bar{C}_{ij}^{e}
% %
% &=&
% \left(
% \delta_{ir}
% -
% \Delta \lambda {}^{n+1} \! \bar{M}_{ir}
% \right)
%
% {}^{n+1} \! \bar{C}_{rl}^{e,trial} \;
% %
% \left(
% \delta_{lj}
% -
% \Delta \lambda {}^{n+1} \! \bar{M}_{lj}
% \right)
% %
% \nonumber \\
% %out&=&
% %out\left(
% %out{}^{n+1} \! \bar{C}_{il}^{e,trial}
% %out-
% %out\Delta \lambda {}^{n+1} \! \bar{M}_{ir} \; {}^{n+1} \! \bar{C}_{rl}^{e,trial} \;
% %out\right)
% %out%
% %out\left(
% %out\delta_{lj}
% %out-
% %out\Delta \lambda {}^{n+1} \! \bar{M}_{lj}
% %out\right)
% %out%
% %out\nonumber \\
% &=&
% {}^{n+1} \! \bar{C}_{ij}^{e,trial}
% -
% \Delta \lambda {}^{n+1} \! \bar{M}_{ir} \; {}^{n+1} \! \bar{C}_{rj}^{e,trial} \;
% -
% \Delta \lambda \;  {}^{n+1} \! \bar{C}_{il}^{e,trial} \; {}^{n+1} \! \bar{M}_{lj}
% \nonumber \\
% & & +
% \Delta \lambda^{2} \;
% {}^{n+1} \! \bar{M}_{ir} \;
% {}^{n+1} \! \bar{C}_{rl}^{e,trial} \;
% {}^{n+1} \! \bar{M}_{lj}
% \label{IIA100}
% \end{eqnarray}
% 

It  is important to note  that  the    Taylor's  series  expansion  from  Equation
\ref{IIA090}  is  a  proper  approximation  for the general non-symmetric tensor
$\bar{M}_{lj}$. That is, the approximate solution given by Equation \ref{IIA130}
is valid for a general anisotropic solid.
 
%\begin{remark}
%\label{remark010}
%  The  Taylor's series expansion from Equation \ref{IIA090} is a proper
%  approximation   for  the general  non-symmetric  tensor  $\bar{M}_{lj}$.  That  is,
%  the approximate  solution  given by Equation \ref{IIA100} is valid for a general
%  anisotropic   solid.   This   contrasts with the  spectral   decomposition  family  of
%  solutions (Simo \cite{Simo92}, Simo and Miehe \cite{Simo92b}) which are  restricted to isotropic solids.
%\end{remark}


The  predictor--corrector relation for  elastic Green strain can then be written
as
%%  
\begin{equation}
{}^{n+1} \! \bar{E}_{ij}^{e}
%
=
{}^{n+1} \! \bar{E}_{ij}^{e,trial}
-
\Delta \lambda \; 
\bar{M}^{C, trial}_{ij}
%
\;\;\;\;\;  \mbox{;} \;\;\;\;\;
%
\bar{M}^{C, trial}_{ij}
=
Sym ( {}^{n+1} \! \bar{C}_{ik}^{e,trial} \; {}^{n+1} \! \bar{M}_{kj} )
\label{MCtrial}
\end{equation}
%%

Similar predictor--corrector relations can be obtain for internal variables,
\begin{equation}
{}^{n+1} \! \xi
%
=
{}^{n+1} \! \xi^{trial}
-
\Delta \lambda \; 
{}^{n+1} \! n^{scalar}
\label{IIA132}
\end{equation}
%%
\begin{equation}
{}^{n+1} \! \eta_{ij}
%
=
{}^{n+1} \! \eta_{ij}^{trial}
-
\Delta \lambda \; 
{}^{n+1} \! n_{ij}^{kinematic}
\label{IIA133}
\end{equation}
%%


The    iterative    algorithm    is    then    developed   using   a   residual
${}^{n+1}\!\tilde{\bf     r}_{\bf     A}$    (using    idea    from    Crisfield
\cite{Crisfield1997}) at step $n+1$
%%
\begin{equation}
{}^{n+1}\!\tilde{\bf r}_{\bf A} 
= {}^{n+1}\!\tilde{\bf s}_{\bf A} 
- \tilde{\bf r}^{trial}_{\bf A} 
+ \Delta \lambda \; {}^{n+1}\!\tilde{\bf m}_{\bf A}
\label{residual}
\end{equation}
%%
where   $\tilde{\bf   r}_A$,   $\tilde{\bf  s}_A$,  and  $\tilde{\bf  m}_A$  are
generalized  vectors  or  matrices  in which the element components are tensors 
and scalar: 
%The indices with bold capital letter $\textbf{A}$, may be,
%for example, (ij) for rank-2 tensors or nothing for scalars.
%%
\begin{equation}
\tilde{\bf r}_{\bf A}
=
\left\{ \begin{array}{c}
\bar{R}_{ij} \\  r^{scalar}  \\  r^{kinematic}_{ij}
\end{array} \right\}
%%\label{vectorRT}
%%\end{equation}
%%
,\;\;\;
%%
%%\begin{equation}
\tilde{\bf s}_{\bf A}
=
\left\{ \begin{array}{c}
\bar{E}^e_{ij} \\ \xi \\ \eta_{ij}
\end{array}\right\}
%%\label{vectorsT}
%%\end{equation}
%%
,\;\;\;
%%
%%\begin{equation}
\tilde{\bf m}_{\bf A}
=
\left\{ \begin{array}{c}
\bar{M}^{C}_{ij} \\ n^i \\ n^k_{ij}
\end{array}\right\}
\label{vectorMT}
\end{equation}
%%%%
%These notations will used in other formulations but will not be explicitly explained any more.


Linearization of residual relation (Eq. (\ref{residual})) 
while noting that
trial  variables  are  fixed  during  the  iteration  process  and  omitting  the
superscript $n+1$, yields
%%
\begin{equation}
\tilde{\bf r}^{old}_{\bf A}
+
\Delta^2 \lambda \; \tilde{\bf m}_{\bf A}
+
{\bf \mathbb C}^{-1}_{\bf AB} \Delta \tilde{\bf s}_{\bf B} 
= 
{\bf 0}
\label{vectorEeAlpha}
\end{equation}
%%
where
\begin{equation}
{\bf \mathbb C}_{\bf AB}
=
\left[ \begin{array}{ccc}
I_{ijkl} + \Delta\lambda \dfrac{ \partial \bar{M}^{C}_{ij}}{ \partial \bar{T}_{mn}} 
\bar{\mathcal L}_{mnkl}^{M}
& \Delta\lambda \dfrac{ \partial \bar{M}^{C}_{ij}}{ \partial q} K
& \Delta\lambda \dfrac{ \partial \bar{M}^{C}_{ij}}{ \partial a_{mn}} {\mathcal H}_{mnkl} \\
\Delta\lambda \dfrac{ \partial n^{scalar}}{ \partial \bar{T}_{ij}} \bar{\mathcal L}_{ijkl}^{M}
& 1 + \Delta\lambda \dfrac{ \partial n^{scalar}}{ \partial q} K
& \Delta\lambda \dfrac{ \partial n^{scalar}}{ \partial a_{ij}} {\mathcal H}_{ijkl} \\
\Delta\lambda \dfrac{ \partial n^{kinematic}_{ij}}{ \partial \bar{T}_{mn}} \bar{\mathcal L}_{mnkl}^{M}
& \Delta\lambda \dfrac{ \partial n^{kinematic}_{ij}}{ \partial q} K
& I_{ijkl} + \Delta\lambda \dfrac{ \partial n^{kinematic}_{ij}}{ \partial a_{mn}} {\mathcal H}_{mnkl}
\end{array} \right]^{-1}
\label{matrixC}
\end{equation}
and $I_{ijkl}$ is the unit fourth order tensor.


%%
The generalized matrix ${\bf \mathbb C}$ has the following structure
%%
\begin{equation}
{\bf \mathbb C}
=
\left[ \begin{array}{ccc}
{\bf \mathbb C}^{\left\langle 11 \right\rangle} & 
{\bf \mathbb C}^{\left\langle 12 \right\rangle} & 
{\bf \mathbb C}^{\left\langle 13 \right\rangle} \\
{\bf \mathbb C}^{\left\langle 21 \right\rangle} & 
{\bf \mathbb C}^{\left\langle 22 \right\rangle} & 
{\bf \mathbb C}^{\left\langle 23 \right\rangle} \\
{\bf \mathbb C}^{\left\langle 31 \right\rangle} & 
{\bf \mathbb C}^{\left\langle 32 \right\rangle} & 
{\bf \mathbb C}^{\left\langle 33 \right\rangle}  
\end{array} \right]
\label{matrixC20}
\end{equation}
%%
where 
${\bf \mathbb C}^{\left\langle 11 \right\rangle}$, 
${\bf \mathbb C}^{\left\langle 13 \right\rangle}$, 
${\bf \mathbb C}^{\left\langle 31 \right\rangle}$, 
${\bf \mathbb C}^{\left\langle 33 \right\rangle}$ are fourth order tensors,
${\bf \mathbb C}^{\left\langle 12 \right\rangle}$, 
${\bf \mathbb C}^{\left\langle 21 \right\rangle}$, 
${\bf \mathbb C}^{\left\langle 23 \right\rangle}$, 
${\bf \mathbb C}^{\left\langle 32 \right\rangle}$ are second order  tensors, and  
${\bf \mathbb C}^{\left\langle 22 \right\rangle}$ is a scalar. 
%%



Linearized form of the consistency condition 
$
{}^{n+1}\!\mathscr{F} \left( {}^{n+1}\!\bar{T}_{ij}, {}^{n+1}\!q, {}^{n+1}\!a_{ij}\right)
= 0
$
at time $t+1$ yields 
%
\begin{equation}
\tilde{\bf f}^{T}_{\bf A} \Delta \tilde{\bf s}_{\bf A} + \mathscr{F}^{old} = {\bf 0}
\label{vectorF}
\end{equation}
%%
where    $\Delta    \tilde{\bf    s}_{\bf  A}$  is  the  differential  form of a
generalized vector $\tilde{\bf
s}_{\bf A}$ from Equation (\ref{vectorMT}) and 
(generalized) vector $\tilde{\bf f}_{\bf A}$ is defined as
%
\begin{equation}
\tilde{\bf f}_{\bf A}
=
\left\{ \begin{array}{ccc}
\dfrac{ \partial \mathscr{F}}{ \partial \bar{T}_{ij}} \bar{\mathcal L}_{ijkl}^{M}
& 
\dfrac{ \partial \mathscr{F}}{ \partial q} K
&  
\dfrac{ \partial \mathscr{F}}{ \partial a_{ij}} {\mathcal H}_{ijkl}
\end{array} \right\}^T
\label{vectorFT}
\end{equation}
%%
%%\begin{equation}
%%\left\{ \begin{array}{ccc}
%%\dfrac{ \partial \mathscr{F}}{ \partial \bar{T}_{ij}} \bar{\mathcal L}_{ijkl}^{M}
%%&
%%\dfrac{ \partial \mathscr{F}}{ \partial q} K
%%& 
%%\dfrac{ \partial \mathscr{F}}{ \partial a_{ij}} {\mathcal H}_{ijkl}
%%\end{array} \right\}
%%\left\{ \begin{array}{l} 
%%\Delta\bar{E}^{e}_{kl} \\ \Delta \xi \\ \Delta \eta_{kl}
%%\end{array}\right\}
%%+ \mathscr{F}^{old} = 0
%%\label{vectorF}
%%\end{equation}



By  using  Equations (\ref{vectorEeAlpha}) and (\ref{vectorF}) one can solve for
$\Delta \lambda$
\begin{equation}
\Delta^{2} \lambda
= \dfrac
{ \mathscr{F}^{old} - \tilde{\bf f}^T_{\bf A} {\bf \mathbb C}_{\bf AB} \tilde{\bf r}_{\bf B} }
{ \tilde{\bf f}^T_{\bf C} {\bf \mathbb C}_{\bf CD} \tilde{\bf m}_{\bf D} }
\label{delta2gamma}
\end{equation}
%% 
%%%
which can be used in conjunction with Equation \ref{vectorEeAlpha} to solve for 
an increment in (quasi) vector  $\tilde{\bf s}_{\bf A}$ 
%
\begin{equation}
\Delta \tilde{\bf s}_{\bf A}
=
\left\{ \begin{array}{l} 
\Delta\bar{E}^{e}_{kl} \\ \Delta \xi \\ \Delta \eta_{kl}
\end{array}\right\}
= 
-
{\bf \mathbb C}_{\bf AB} 
\left( \tilde{\bf r}_{\bf B} 
+ \Delta^{2} \lambda \; \tilde{\bf m}_{\bf B} \right)
\label{vectorEeandAlpha}
\end{equation}
%
that contains
increments  of elastic Green strain tensor ($\Delta \bar{E}^e_{ij}$, used to find Second
Piola--Kirchhoff  stress),  scalar ($\Delta \xi$)  and kinematic ($\Delta \eta_{ij}$) internal
variables. 
Iterations  continue  until both consistency condition 
($
{}^{n+1}\!\mathscr{F} \left( {}^{n+1}\!\bar{T}_{ij}, {}^{n+1}\!q, {}^{n+1}\!a_{ij}\right)
= 0
$)
and minimization of the norm of the residual ($\tilde{\bf r}_{\bf A} = 0$) are satisfied  
(within some tolerance).










In addition to the fully implicit algorithm, described above, the 
algorithmic, consistent stiffness  tensor  in intermediate 
configuration $\bar{\mathcal L}^{AS}_{ijmn}$  is derived similarly as
%
\begin{equation}
\bar{\mathcal L}^{AS}_{ijmn}
= \bar{\mathcal L}_{ijkl} \hat{\bf \mathbb C}^{\left\langle 11 \right\rangle}_{klmn}
\label{ATSL}
\;\;\;\;\;\;\; \mbox{from} \;\;\;\;\;\;\;
d\bar{S}_{ij}
= \bar{\mathcal L}_{ijkl} d\bar{E}^{e}_{kl}
= \bar{\mathcal L}_{ijkl} \hat{\bf \mathbb C}^{\left\langle 11 \right\rangle}_{klmn} d\bar{E}^{e,trial}_{mn}
\end{equation}
%
where
$\hat{\bf \mathbb C}^{\left\langle 11 \right\rangle}$ 
is the upper-left $(1,1)$ block of the  generalized matrix of $\hat{\bf \mathbb C}$, 
\begin{equation}
\hat {\bf \mathbb C}_{\bf AB}
:=
%\left( 
{\bf \mathbb C}_{\bf AB} 
- \dfrac
{( {\bf \mathbb C}_{\bf AM} \tilde{\bf m}_{\bf M} ) ( \tilde{\bf f}^{T}_{\bf N} {\bf \mathbb C}_{\bf NB} )}
{\tilde{\bf f}^{T}_{\bf C} {\bf \mathbb C}_{\bf CD} \tilde{\bf m}_{\bf D}}
%  \right)
\label{ATShatC}
\end{equation}
%
and    $\tilde{\bf  m}_{\bf  M}$, ${\bf  \mathbb  C}_{\bf  AM}$  and  $\tilde{\bf
f}^{T}_{\bf N}$ were defined by equations   (\ref{vectorMT}), (\ref{matrixC}) and 
(\ref{vectorFT}).






% 
%  
%  %
%  %
%  The algorithmic stiffness tensor in reference configuration can be 
%  also obtained as
%  %
%  \begin{equation}
%  {\mathcal L}^{AS}_{ijpq}
%  :=
%  \dfrac{dS_{ij}}{dE^{trial}_{pq}}
%  = 
%  \left(
%  F^{p-1}_{ki}F^{p-1}_{lj} \bar{\mathcal L}^{AS}_{klmn} 
%  - 
%  2 \bar{\mathcal L}^{P}_{ijkl} {\mathcal Q}_{klmn}
%  \right)
%  {}^{n} \! F_{mp}^{p-1} \; {}^{n} \! F_{nq}^{p-1}  
%  \label{ATSLref}
%  \end{equation}
%  %
%  where
%  %
%  \begin{equation}
%  {\mathcal Q}_{klmn} := {\bf D}_{klmn}^{\left\langle 1 \right\rangle}
%  \label{CalQ}
%  \end{equation}
%  %
%  
%  {\Large \sc Zhao: what is this ${\bf D}_{klmn}^{\left\langle 1 \right\rangle}$, can you
%  tell me where it comes from (reference to your dissertation...)}
%  

The  algorithmic stiffness tensor is used on the finite element level to provide
for fast (quadratic once within Newton trust region) convergence.

%  The incremental finite element equations, are given in the initial configuration
%  as
%  %
%  \begin {equation} 
%  K_{AaBb} \; \Delta \bar{u}_{Bb} \
%  = \Delta f^{ex}_{Aa}
%  \label{FEFmInitial_L1}
%  \end{equation}
%  %%
%  where  $K_{AaBb}$  is the stiffness matrix/tensor\footnote{In this notation, the
%  stiffness  is  a tensor (see Zienkiewicz and Taylor \cite{FEM-IV-1} chapter 6.)},
%  $\Delta  \bar{u}_{Bb}$  is  an  increment  of  unknown nodal displacements while
%  $\Delta  f^{ex}_{Aa}$  is  the increment in nodal loads. The stiffness tensor is
%  given as 
%  %
%  \begin{eqnarray}
%  K_{AaBb} 
%  =
%  \int_{\Omega_0} 
%  H_{A,J} \left( 
%  F_{aI} F_{bK} \mathcal{L}^{AS}_{IJKL} + \delta_{ab} S_{JL}
%  \right) 
%  H_{B,L} d V
%  \label{FEK01}
%  \end{eqnarray}
%  %
%  while the load tensor/vector is given as
%  %
%  \begin{eqnarray}
%  \Delta f^{ex}_{Aa}
%  =
%    \int_{\partial \Omega_0}  H_A \Delta t_i  d A
%  + \int_{\Omega_0} \rho_0 H_A \Delta b_i d V
%  \label{FE02}
%  \end{eqnarray}
%  %
%  In  the  above  equations,  $H_A$  are  interpolating functions, $F_{aI}$ is the
%  deformation  gradient, $\mathcal{L}^{AS}_{IJKL}$ is the stiffness tensor, $S_{JL}$ is
%  the  second  Piola--Kirchhoff  stress  tensor,  and  $t_i$ and $b_i$ are surface
%  tractions and body forces respectively.
%  
%%


%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section {Numerical Example}

A  numerical example with  non-coaxial anisotropic model is used to demonstrated
described  developments.  Material  model is an extension of von Mises yield and
potential  surfaces  with  kinematic hardening. The kinematic hardening features
initial  back  stress (which might have developed from previous (cyclic) loading
stage)  that  is  not coaxial with the developing conjugate hyperelastic strain.
This  non--coaxiality  will  in  this  case  prevent  use  of  any  constitutive
integration algorithms that is relying on assumption of coaxiality of stress and
its conjugate hyperelastic strain.

The material model yield function is defined as
%%%%
\begin{equation}
\mathscr{F}(\tau_{ij}, a_{ij})
=
\dfrac{3}{2} (Dev(\tau_{ij}) - x_{ij} - a_{ij}) (Dev(\tau_{ij}) - x_{ij} - a_{ij})  - Y_0^2 = 0
\label{vMYield_Simo}
\end{equation}
%%%%
where $Dev(\tau_{ij})$
is the deviatoric part of Kirchhoff stress in the current configuration. 
The material constant $Y_0$ represents the yield strength, 
while  $x_{ij}$  represents  constant  back  stress  tensor, and $a_{ij}$ is the
kinematic hardening (internal variable, stress-like) tensor.
It  is easy to obtain the first derivative of the yield function with respect to
Kirchhoff stress as
%%%%
\begin{equation}
\dfrac{\partial \mathscr{F}}{\partial \tau_{ij}}
= 3 [Dev(\tau_{ij}) - x_{ij} - a_{ij}]
\label{vMYield_dev_Simo}
\end{equation}
%%%%
If $x_{ij}$ is zero tensor, 
Kirchhoff stress $\tau_{ij}$ and first derivatives of yield function $\partial \mathscr{F} / \partial \tau_{ij}$ 
are  co-axial.  However,  if  $x_{ij}$  is any non--zero tensor  (except for pure
hydrostatic  tensor,  like  Kronecker  delta  $\delta_{ij}$),  then the Kirchhoff
stress
$\tau_{ij}$  and  and  derivatives  of  yield  function  $\partial \mathscr{F} /
\partial  \tau_{ij}$ are not coaxial. It should be noted that the yield function
given   in   Equation  (\ref{vMYield_Simo})  is  defined  in  terms  of  Kirchhoff
stress,  while  it  will  be redefined (in the same form) in terms of the Mandel
stress  $\bar{T}_{ij}$  in  the intermediate configuration in order to use it in
developed algorithms.




The  hyperelasticity  is of compressible Neo-Hookean type with the strain energy
function given as
%%%%
\begin{equation} 
\Psi 
= \dfrac{1}{2}\lambda_0\left( \ln J \right)^2 
- G_0(\ln J) 
+ \dfrac{1}{2}G_0\left( tr(\bar{C}^{e}) -3 \right)
\label{cnHstrainEnergy}
\end{equation}
%%
where  $G_0$ is shear modulus and $\lambda_0=K_0-(2/3)G_0$ is Lam{\' e} constant
while $K_0$ is bulk modulus. 
Volume   ration $J$ is defined as  $J^2 = \det(\bar{C}_{ij}^{e})$ while trace of
the elastic  Lagrange-Green strain  is $tr(\bar{C}_{ij}^{e})= \bar{C}^{e}_{ii}$.
%
The  material  constants  are  $K_0  = 1971.67$~kPa, $G_0 = 4225.50$~kPa, $Y_0 =
10$~kPa, while linear kinematic hardening law is adopted with  hardening modulus
$H = 80$~kPa.

A  random  non--zero  tensor  $x_{ij}$  can be used to show non--coaxiality, but
here, without losing generality, $x_{ij}$ is given as 
%%%%
\begin{equation}
x_{ij} = 
\left[ 
\begin{array}{ccc}
0.0 &  0.0   &  0.0 \\
0.0 &  0.0   &  x \\
0.0 &  x     &  0.0
\end{array}
\right]
\label{x_ij}
\end{equation}
%%%%
where  $x$  is  a  given  scalar controlling the size of the back stress and has
units of stress.
%
%
The  deformation  is  defined  through  a simple deformation gradient for simple
shear given as 
%%%%
\begin{equation}
F_{ij} = 
\left[ 
\begin{array}{ccc}
1.0 &  0.0 &  \gamma   \\
0.0 &  1.0 &  0.0      \\
0.0 &  0.0 &  1.0
\end{array}
\right]
\label{f_ij}
\end{equation}
%%%%
where $\gamma$ is the shear ratio. 




In this example, a total of 100 increments  $\Delta \gamma = 0.05\%$ 
are used. 
The eigenvectors of the left Cauchy--Green deformation $b_{ij}$ and the 
Cauchy stress $\sigma_{ij}$  ($\sigma_{ij} = J \tau_{ij}$, so that $\sigma_{ij}$
and $\tau_{ij}$ are always coaxial) are compared for collinearity. 
Table \ref{eigenvectors} list those eigenvectors at $\gamma = 5\%$ 
with initial non--zero and zero back stress as expressed in Equation (\ref{x_ij}). 
%%     %%
%%     \begin{table}[!htbp]
%%     \caption{ \label{NCAeigenvectors}      
%%     Calculated eigenvectors at $\gamma = 5\%$ with x = 5.0.}
%%     \begin{center}
%%     \begin{tabular}{|c||r|r|r|} 
%%     \hline
%%     {}           & eigenvector 1   & eigenvector 2   & eigenvector 3 \\ 
%%     \hline \hline
%%     {}           & 0.71589 & 0.00000 & -0.69822 \\
%%     $b_{ij}$     & 0.00000 & 1.00000 &  0.00000 \\
%%     {}           & 0.69822 & 0.00000 &  0.71589 \\
%%     \hline 
%%     {}             & 0.64279 & -0.43400 & -0.63124 \\
%%     $\sigma_{ij}$  & 0.31581 &  0.90088 & -0.29780 \\
%%     {}             & 0.69792 & -0.00794 &  0.71613 \\
%%     \hline
%%     \end{tabular}
%%     \end{center}
%%     \end{table}
It  is  obvious  that  for the non--zero back stress (anisotropic case) that the
eigenvectors of  $b_{ij}$ and $\sigma_{ij}$ are non-coaxial.
On  the  other  hand, for zero back stress (isotropic case) the eigenvectors are
equal, that is, the stress and strain tensors are coaxial.
%
%  %
%  %
%  On the other hand, Table \ref{NCAeigenvectors_zero} shows those eigenvector at $\gamma = 5\%$ 
%  with initial zero back stress as expressed in \ref{x_ij}. 
%  %
%  \begin{table}[!htbp]
%  \caption{ \label{eigenvectors}      
%  Calculated eigenvectors at $\gamma = 5\%$ with x = 0.0.}
%  \begin{center}
%  \begin{tabular}{|c||r|r|r|} 
%  \hline
%  {}           & eigenvector 1   & eigenvector 2   & eigenvector 3 \\ 
%  \hline \hline
%  {}           & 0.71589 & 0.00000 & -0.69822  \\
%  $b_{ij}$     & 0.00000 & 1.00000 &  0.00000  \\
%  {}           & 0.69822 & 0.00000 &  0.71589  \\
%  \hline 
%  {}             & 0.71589 & 0.00000 & -0.69822  \\
%  $\sigma_{ij}$  & 0.00000 & 1.00000 &  0.00000  \\
%  {}             & 0.69822 & 0.00000 &  0.71589  \\ 
%  \hline
%  \end{tabular}
%  \end{center}
%  \end{table}
%  
%  
%  
%  
%  
%  
%  
%  
%  
%%
\begin{table}[!htbp]
\caption{ \label{eigenvectors}      
Calculated    eigenvectors    (EV)   at  shear deformation of $\gamma = 5\%$ for
non--coaxial ($x = 5.0$~kPa) and coaxial ($x=0.0$~kPa) cases.}
\begin{center}
\begin{tabular}{|c||r|r|r||r|r|r|} 
\hline
\multicolumn{1}{|c||}{} & 
\multicolumn{3}{|c||}{ non--coaxial, $x = 5.0$~kPa     } &
\multicolumn{3}{|c|}{ coaxial, $x = 0.0$~kPa    }
  \\ \hline
\hline
           & EV 1   & EV 2   & EV 3 & EV 1   & EV 2   & EV 3 \\ 
\hline
{}           & 0.71589 & 0.00000 & -0.69822 & 0.71589 & 0.00000 & -0.69822  \\
$b_{ij}$     & 0.00000 & 1.00000 &  0.00000 & 0.00000 & 1.00000 &  0.00000  \\
{}           & 0.69822 & 0.00000 &  0.71589 & 0.69822 & 0.00000 &  0.71589  \\
\hline 
{}             & 0.64279 & -0.43400 & -0.63124 & 0.71589 & 0.00000 & -0.69822 \\
$\sigma_{ij}$  & 0.31581 &  0.90088 & -0.29780 & 0.00000 & 1.00000 &  0.00000 \\
{}             & 0.69792 & -0.00794 &  0.71613 & 0.69822 & 0.00000 &  0.71589 \\
\hline
\end{tabular}
\end{center}
\end{table}
%  Indeed, $b_{ij}$ and $\sigma_{ij}$ are found to be non-coaxial. 
%  %
%  %
%  \begin{table}[!htbp]
%  \caption{ \label{NCAeigenvectors_zero}      
%  Calculated eigenvectors at $\gamma = 5\%$ with x = 0.0.}
%  \begin{center}
%  \begin{tabular}{|c||r|r|r|} 
%  \hline
%  {}           & EV 1   & EV 2   & EV 3 \\ 
%  \hline \hline
%  {}           & 0.71589 & 0.00000 & -0.69822  \\
%  $b_{ij}$     & 0.00000 & 1.00000 &  0.00000  \\
%  {}           & 0.69822 & 0.00000 &  0.71589  \\
%  \hline 
%  {}             & 0.71589 & 0.00000 & -0.69822  \\
%  $\sigma_{ij}$  & 0.00000 & 1.00000 &  0.00000  \\
%  {}             & 0.69822 & 0.00000 &  0.71589  \\ 
%  \hline
%  \end{tabular}
%  \end{center}
%  \end{table}
%  %
%  
%  %
%  This time $b_{ij}$ and $\sigma_{ij}$ are coaxial. 
%
It  follows  that  with even simple non--coaxial material model the constitutive
integrations  cannot  use  any of the algorithms that make coaxiality (of stress
and  its conjugate hyperelastic strain,  and/or  stress and the gradient of yield function in
stress   space)   assumption.   Rather,   integration   algorithms   that   take
non--coaxiality
into account, like the one described above, need to be used.



It    is   also  interesting to compare  results for shearing of material models
with back stress that is (a) coaxial  (using $x = 0.0$) and
(b) non--coaxial  (using  $x  =  5.0$). Such results, ($\sigma_{13}$
versus $\gamma$) are presented in Figure \ref{noncoaxial}(A).
%%%%%%%%%%%%%%
\begin{figure}[!htbp]
\begin{center}
\includegraphics[width=0.65\textwidth]{/home/jeremic/tex/works/Thesis/ZhaoCheng/thesis/file/LDEP/noncoaxial/noncoaxial_new.ps.eps}
%\includegraphics[width = 0.7\textwidth]{noncoaxial.eps}
\caption{\label{noncoaxial}  
Simulation results ($\gamma$ versus $\sigma_{13}$)  
for material model with 
(a) coaxial  (using $x = 0.0$) and
(b) non--coaxial  (using  $x  =  5.0$) back stress.} 
\end{center}
\end{figure}
%%%%%%%%%%%%%%
The  response  of the coaxial linear kinematic hardening material model (with $x=0.0$~kPa)
is  represented  by  a  bilinear  response (as expected). On the other hand, the
response  of the similar material  model with initial back stress ($x=5.0$~kPa),
is  nonlinear  from  the  moment it yields. This is also expected as the initial
back stress is applied to different component ($x_{23} = x_{32} = 5.0$~kPa) than
the  one  which  is  loaded  ($F_{13}$). The yielding also happens sooner as the
yield  surface  is  shifted  toward stress origin by the presence of the initial
back  stress.  The nonlinearity of response follows from the fact that the yield
and  potential  surfaces harden according to Prager's rule, following the increment
of  plastic  deformation,  which  now changes non--linearly as there was a shift (initial
back  stress)  in  position  of  those  two surfaces. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section {Summary}


This  short paper brings forward the issue of constitutive integration for large
deformation hyperelastic--plastic material models which develop  non--coaxiality     of    principal
stress       and    its conjugate hyperelastic strain. In addition to discussion
on    non--coaxiality,    briefly  described  was  a  constitutive  integration 
algorithm   that  resolve the limitations to hyperelastic isotropy, plastic flow
direction isotropy and yield function isotropy.
A simple example was used to
illustrate discussed issues through  an anisotropic model with initial back--stress and its constitutive integration. 




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%  \section*{Acknowledgment}
%  
%  The  work presented in this paper was supported in part by a number of Agencies listed below:
%  %
%  Earthquake    Engineering   Research  Centers  Program  of  the National Science
%  Foundation  under  Award  Number NSF--EEC--9701568 (cognizant program director Dr. Joy
%  Pauschke);
%  %
%  Civil  and Mechanical System program, Directorate of Engineering of the National
%  Science  Foundation, under Award NSF--CMS--0337811 (cognizant program director Dr.
%  Steve  McCabe);
%  %
%  California  Department  of  Transportation  (Caltrans)  under  Award \# 59A0433,
%  (cognizant program director Dr. Saad El-Azazy).





                    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
 
\begin{thebibliography}{10}

\bibitem{Bathe85}
K.J. Bathe, R.~Slavkovi{\'c}, and M.~Koji{\'c}.
\newblock On large strain elasto--plastic and creep analysis.
\newblock In Wunderlich Bergan, Bathe, editor, {\em Finite Element Methods for
  Nonlinear Problems, Europe--US Symposium, Trondheim, Norway}, pages 176--190,
  1985.

\bibitem{Bilby57}
B.~A. Bilby, L.~R.~T. Gardner, and A.~N. Stroh.
\newblock Continuous distributions of dislocations and the theory of
  plasticity.
\newblock In {\em {$IX^{e}$} Congr{\`e}s International de M{\`e}chanique
  Appliqu{\`e}e}, volume {VIII}, pages 35--44, Universit{\`e} de Bruxelles, 50.
  Av{\`e}nue Franklin Roosevelt, 1957.

\bibitem{Crisfield1997}
M.~A. Crisfield.
\newblock {\em Non--Linear Finite Element Analysis of Solids and Structures
  Volume~1: Advanced Topics}.
\newblock John Wiley and Sons, Inc. New York, 605 Third Avenue, New York, NY
  10158--0012, USA, 1997.

\bibitem{Cuitino92}
A.~Cuiti{\~n}o and M.~Ortiz.
\newblock A material -- independent method for extending stress update
  algorithms from small--strain plasticity to finite plasticity with
  multiplicative kinematics.
\newblock {\em Engineering Computations}, 9:437--451, 1992.

\bibitem{Eterovic90}
Adrian~Luis Eterovic and {Klaus}-{J{\"u}rgen} Bathe.
\newblock A hyperelastic -- based large strain elasto -- plastic constitutive
  formulation with combined isotropic -- kinematic hardening using the
  logarithmic stress and strain measures.
\newblock {\em International Journal for Numerical Methods in Engineering},
  30:1099--1114, 1990.

\bibitem{Fox68}
N.~Fox.
\newblock On the continuum theories of dislocations and plasticity.
\newblock {\em Quarterly Journal of Mechanics and Applied Mathematics},
  XXI(1):67--75, 1968.

\bibitem{Green1965}
A.~Green and P~Naghdi.
\newblock A general theory of an elastic-plastic continuum.
\newblock {\em Arch Rational Mechanical Analysis}, 18:251--281, 1965.

\bibitem{Hill50}
R.~Hill.
\newblock {\em The Mathematical Theory of Plasticity}.
\newblock The Oxford Engineering Science Series. Oxford at the Clarendon Press,
  1st. edition, 1950.

\bibitem{Jeremic97f}
Boris Jeremi{\' c}, Kenneth Runesson, and Stein Sture.
\newblock A model for elastic--plastic pressure sensitive materials subjected
  to large deformations.
\newblock {\em International Journal of Solids and Structures},
  36(31/32):4901--4918, 1999.

\bibitem{Karush1939}
W.~Karush.
\newblock Minima of functions of several variables with inequalities as side
  constraints.
\newblock Master's thesis, University of Chicago, Chicago, IL., 1939.

\bibitem{Kroner60}
Ekkehart Kr{\"o}ner.
\newblock Algemeine kontinuumstheorie der versetzungen und eigenspanungen.
\newblock {\em Archive for Rational Mechanics and Analysis}, 4(4):273--334,
  1960.

\bibitem{Kuhn-Tucker50}
H.~W. Kuhn and A.~W. Tucker.
\newblock Nonlinear programming.
\newblock In Jerzy Neyman, editor, {\em Proceedings of the Second Berkeley
  Symposium on Mathematical Statistics and Probability}, pages 481 -- 492.
  University of California Press, July 31 -- August 12 1950 1951.

\bibitem{Lee69}
E.~H. Lee.
\newblock Elastic--plastic deformation at finite strains.
\newblock {\em Journal of Applied Mechanics}, 36(1):1--6, 1969.

\bibitem{Lee67}
E.~H. Lee and D.~T. Liu.
\newblock Finite--strain elastic--plastic theory with application to
  plane--wave analysis.
\newblock {\em Journal of Applied Physics}, 38(1):19--27, January 1967.

\bibitem{Miehe2002}
C.~Miehe, N.~Apel, and M.~Lambrecht.
\newblock Anisotropic additive plasticity in the logarithmic strain space:
  modular kinematic formulation and implementation based on incremental
  minimization principles for standard materials.
\newblock {\em Computer Methods in Applied Mechanics and Engineering},
  191:5383--5425, 2002.

\bibitem{Miehe2004}
Christian Miehe and Nikolas Apel.
\newblock Anisotropic elastic-plastic analysis of shells at large strains. a
  comparison of multiplicative and additive approaches to enhanced finite
  element design and constitutive modelling.
\newblock {\em International Journal for Numerical Methods in Engineering},
  61(12):2067 -- 2113, 2004.

\bibitem{Papadopoulos98}
Panayaiotis Papadopoulos and Jia Lu.
\newblock A general framework for the numerical solution of problems in finite
  elasto-plasticity.
\newblock {\em International Journal for Computer Methods in Applied Mechanics
  and Engineering}, 159:1--18, 1998.

\bibitem{Pearson74}
Carl~E. Pearson, editor.
\newblock {\em Handbook of Applied Mathematics}.
\newblock Van Nostrand Reinhold Company, 1974.

\bibitem{Peric92}
Djordje Peri{\'c}, D.~R.~J. Owen, and M.~E. Honnor.
\newblock A model for finite strain elasto--plasticity based on logarithmic
  strains: Computational issues.
\newblock {\em Computer Methods in Applied Mechanics and Engineering},
  94:35--61, 1992.

\bibitem{Prevost87a}
Jean~H. Prevost.
\newblock Modeling the behavior of geomaterial.
\newblock In Sayed~M. Sayed, editor, {\em Geotechnical Modeling and
  Applications}, pages 8--75. Gulf Publishing Company, Houston, London, Paris,
  Tokyo, 1987.

\bibitem{Simo92}
J.~C. Simo.
\newblock Algorithms for static and dynamic multiplicative plasticity that
  preserve the classical return mapping schemes of the infinitesimal theory.
\newblock {\em Computer Methods in Applied Mechanics and Engineering},
  99:61--112, 1992.

\bibitem{Simo92b}
J.~C. Simo and C.~Miehe.
\newblock Associative coupled thermoplasticity at finite strains: Formulation,
  numerical analysis and implementation.
\newblock {\em Computer Methods in Applied Mechanics and Engineering},
  98:41--104, 1992.

\bibitem{local-12}
J.~C. Simo and M.~Ortiz.
\newblock A unified approach to finite deformation elastoplastic analysis based
  on the use of hyperelastic constitutive equations.
\newblock {\em Computer Methods in Applied Mechanics and Engineering},
  49:221--245, 1985.

\bibitem{Simo88a}
Juan~C. Simo.
\newblock A framework for finite strain elastoplasticity based on maximum
  plastic dissipation and the multiplicative decomposition: Part i. continuum
  formulation.
\newblock {\em Computer Methods in Applied Mechanics and Engineering},
  66:199--219, 1988.
\newblock TA345. C6425.

\bibitem{Simo88b}
Juan~C. Simo.
\newblock A framework for finite strain elastoplasticity based on maximum
  plastic dissipation and the multiplicative decomposition: Part ii.
  computational aspects.
\newblock {\em Computer Methods in Applied Mechanics and Engineering},
  68:1--31, 1988.
\newblock TA345. C6425.

\end{thebibliography}

          
%\bibliography{refmech,refcomp,zhao}
%\bibliographystyle{plain}
%%\bibliographystyle{unsrt}
%%\bibliographystyle{alpha}
%%\bibliographystyle{ascelike}
%%\bibliographystyle{acm}
%%\bibliographystyle{agsm}
%%\bibliographystyle{jmb}
 
\end{document}









