\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 HyperelastoPlasticity of Anisotropic Materials}
\title{On Large Deformation HyperelastoPlasticity 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}{EEC9701568, CMS0337811}
\noreceived{}
\norevised{}
\noaccepted{}
\keywords{anisotropy, hyperelasticplastic, 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
hyperelastoplastic material models. Of particular interest is the appropriate
treatment of cases with noncoaxiality of principal stress and its conjugate
hyperelastic strain. The noncoaxiality can result from hyperelastic
anisotropy, anisotropic flow directions or anisotropic yield function. In
addition to discussing issues related to noncoaxiality, a constitutive
integration algorithm is presented in some detail that is designed to be
applicable to anisotropic noncoaxial (as well as to isotropic)
hyperelasticplastic 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 hyperelastoplastic format mostly used today in large deformation
computations is based on work by Simo and Ortiz \cite{local12}, \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 largedeformation
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
elastoplasticity that is based on the early work of Green and Naghdi
\cite{Green1965}. Developments include provisions for noncollinearity of
principal stress and strain measures. The framework was tested using VonMises
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 noncoaxiality 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 noncoaxiality of principal stresses and
its conjugate hyperelastic strains and its influence on large deformation
hyperelasticplastic integration algorithms. This noncoaxiality can stem from
a number of features present in general hyperelasticplastic 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 noncoaxiality of principal stresses and its conjugate hyperelastic
strains can result if the hyperelastic strain energy function is not an
isotropic function of the elastic EulerGreen strain $\bar{b}^{e}_{ij}$ (or
the elastic LagrangeGreen 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 EulerGreen 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 hyperelasticplastic response precludes them from being applied
to anisotropic hyperelasticplastic models. Our presented algorithm removes
such limitation.
\item[Nonisotropic Flow Rules:]
The noncoaxiality 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 hyperelasticplastic 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 nonisotropic flow
% rule.
\item[Kinematic Hardening:]
The noncoaxiality 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 noncoaxial 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 noncoaxial 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 EulerianGreen 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, nonisotropic 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 EulerianGreen 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
hyperelasticplastic 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 backstress is used to illustrate our approach.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{General Anisotropic HyperelastoPlastic 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 PiolaKirchhoff 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 LagrangeGreen stiffness tensor in the
intermediate configuration, and $\bar{E}_{ij}^{e} := {1}/{2}(
\bar{C}_{ij}^{e}  \delta_{ij})$ is the elastic LagrangeGreen 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 LagrangeGreen 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 stresslike internal variable $q$,
and the kinematic stresslike 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 nonassociated (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 stresslike isotropic internal variables,
% $\beta$ is a parameter to give how fast $K$ changes, one obtain a nonlinear 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 stresslike internal variables are only depend on the corresponding conjugate strainlike 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 EulerianGreen 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 EulerianGreen strain $\bar{b}^{e}_{ij}$
%=== (or the elastic LagrangianGreen strain $\bar{C}^{e}_{ij}$). In general
%=== material modeling (for anisotropic hyperelasticity), no such assumption can be made.
%===
%=== \item[Nonisotropic 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 nonassociated 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 nonisotropic 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 VonMises 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 EulerianGreen 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, nonisotropic 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 EulerianGreen 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 KarushKuhnTucker
(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{KuhnTucker50})
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 nonsymmetric 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 nonsymmetric 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 predictorcorrector 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 predictorcorrector 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 rank2 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
PiolaKirchhoff 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 upperleft $(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^{p1}_{ki}F^{p1}_{lj} \bar{\mathcal L}^{AS}_{klmn}
% 
% 2 \bar{\mathcal L}^{P}_{ijkl} {\mathcal Q}_{klmn}
% \right)
% {}^{n} \! F_{mp}^{p1} \; {}^{n} \! F_{nq}^{p1}
% \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{FEMIV1} 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 PiolaKirchhoff stress tensor, and $t_i$ and $b_i$ are surface
% tractions and body forces respectively.
%
%%
%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section {Numerical Example}
A numerical example with noncoaxial 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 noncoaxiality 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, stresslike) 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 coaxial. However, if $x_{ij}$ is any nonzero 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 NeoHookean 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 LagrangeGreen 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 nonzero tensor $x_{ij}$ can be used to show noncoaxiality, 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 CauchyGreen 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 nonzero 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}{crrr}
%% \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 nonzero back stress (anisotropic case) that the
eigenvectors of $b_{ij}$ and $\sigma_{ij}$ are noncoaxial.
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}{crrr}
% \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
noncoaxial ($x = 5.0$~kPa) and coaxial ($x=0.0$~kPa) cases.}
\begin{center}
\begin{tabular}{crrrrrr}
\hline
\multicolumn{1}{c}{} &
\multicolumn{3}{c}{ noncoaxial, $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 noncoaxial.
% %
% %
% \begin{table}[!htbp]
% \caption{ \label{NCAeigenvectors_zero}
% Calculated eigenvectors at $\gamma = 5\%$ with x = 0.0.}
% \begin{center}
% \begin{tabular}{crrr}
% \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 noncoaxial 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
noncoaxiality
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) noncoaxial (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) noncoaxial (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 nonlinearly 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 hyperelasticplastic material models which develop noncoaxiality of principal
stress and its conjugate hyperelastic strain. In addition to discussion
on noncoaxiality, 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 backstress 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 NSFEEC9701568 (cognizant program director Dr. Joy
% Pauschke);
% %
% Civil and Mechanical System program, Directorate of Engineering of the National
% Science Foundation, under Award NSFCMS0337811 (cognizant program director Dr.
% Steve McCabe);
% %
% California Department of Transportation (Caltrans) under Award \# 59A0433,
% (cognizant program director Dr. Saad ElAzazy).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{thebibliography}{10}
\bibitem{Bathe85}
K.J. Bathe, R.~Slavkovi{\'c}, and M.~Koji{\'c}.
\newblock On large strain elastoplastic and creep analysis.
\newblock In Wunderlich Bergan, Bathe, editor, {\em Finite Element Methods for
Nonlinear Problems, EuropeUS Symposium, Trondheim, Norway}, pages 176190,
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 3544, Universit{\`e} de Bruxelles, 50.
Av{\`e}nue Franklin Roosevelt, 1957.
\bibitem{Crisfield1997}
M.~A. Crisfield.
\newblock {\em NonLinear 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
101580012, USA, 1997.
\bibitem{Cuitino92}
A.~Cuiti{\~n}o and M.~Ortiz.
\newblock A material  independent method for extending stress update
algorithms from smallstrain plasticity to finite plasticity with
multiplicative kinematics.
\newblock {\em Engineering Computations}, 9:437451, 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:10991114, 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):6775, 1968.
\bibitem{Green1965}
A.~Green and P~Naghdi.
\newblock A general theory of an elasticplastic continuum.
\newblock {\em Arch Rational Mechanical Analysis}, 18:251281, 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 elasticplastic pressure sensitive materials subjected
to large deformations.
\newblock {\em International Journal of Solids and Structures},
36(31/32):49014918, 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):273334,
1960.
\bibitem{KuhnTucker50}
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 Elasticplastic deformation at finite strains.
\newblock {\em Journal of Applied Mechanics}, 36(1):16, 1969.
\bibitem{Lee67}
E.~H. Lee and D.~T. Liu.
\newblock Finitestrain elasticplastic theory with application to
planewave analysis.
\newblock {\em Journal of Applied Physics}, 38(1):1927, 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:53835425, 2002.
\bibitem{Miehe2004}
Christian Miehe and Nikolas Apel.
\newblock Anisotropic elasticplastic 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
elastoplasticity.
\newblock {\em International Journal for Computer Methods in Applied Mechanics
and Engineering}, 159:118, 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 elastoplasticity based on logarithmic
strains: Computational issues.
\newblock {\em Computer Methods in Applied Mechanics and Engineering},
94:3561, 1992.
\bibitem{Prevost87a}
Jean~H. Prevost.
\newblock Modeling the behavior of geomaterial.
\newblock In Sayed~M. Sayed, editor, {\em Geotechnical Modeling and
Applications}, pages 875. 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:61112, 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:41104, 1992.
\bibitem{local12}
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:221245, 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:199219, 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:131, 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}