%%%%%%%%%%%%%%%%%%%%%%%% file template.tex %%%%%%%%%%%%%%%%%%%%%%%%%
%
% This is a general template file for the LaTeX package SVJour3
% for Springer journals.          Springer Heidelberg 2010/09/16
%
% Copy it to a new file with a new name and use it as the basis
% for your article. Delete % signs as needed.
%
% This template includes a few options for different layouts and
% content for various journals. Please consult a previous issue of
% your journal as needed.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
% First comes an example EPS file -- just ignore it and
% proceed on the \documentclass line
% your LaTeX will extract the file if required
\begin{filecontents*}{example.eps}
%!PS-Adobe-3.0 EPSF-3.0
%%BoundingBox: 19 19 221 221
%%CreationDate: Mon Sep 29 1997
%%Creator: programmed by hand (JK)
%%EndComments
gsave
newpath
  20 20 moveto
  20 220 lineto
  220 220 lineto
  220 20 lineto
closepath
2 setlinewidth
gsave
  .4 setgray fill
grestore
stroke
grestore
\end{filecontents*}
%
\RequirePackage{fix-cm}
%
\documentclass{svjour3}                     % onecolumn (standard format)
%\documentclass[smallcondensed]{svjour3}     % onecolumn (ditto)
%\documentclass[smallextended]{svjour3}       % onecolumn (second format)
% \documentclass[twocolumn]{svjour3}          % twocolumn
%
\smartqed  % flush right qed marks, e.g. at end of proof
%
% \usepackage{mathptmx}      % use Times fonts if available on your TeX system
%
% insert here the call for the packages your document requires
%==========================


\usepackage{amsmath}
\usepackage{graphicx}
\usepackage{hyperref}
\usepackage{multirow}
\usepackage{bm}
% ==========================
% For subfigure
\usepackage{lastpage}
\usepackage{fancyhdr}
\usepackage{float}
% \usepackage{caption}
\usepackage{subcaption}
\captionsetup{compatibility=false}
% add by Kaveh
\usepackage{color}
\usepackage{soul}
% ==========================
% % For envolope sign before the corresponding authors.
% \usepackage[misc,geometry]{ifsym} 
%
% please place your own definitions here and don't use \def but
% \newcommand{}{}

\usepackage{tcolorbox}


% Bibliography package for abbrvnat
\usepackage[numbers,round,colon,authoryear]{natbib}



% ----------------------------------------------------------------------
%
% TIME OF DAY
%
%
\newcount\hh
\newcount\mm
\mm=\time
\hh=\time
\divide\hh by 60
\divide\mm by 60
\multiply\mm by 60
\mm=-\mm
\advance\mm by \time
\def\hhmm{\number\hh:\ifnum\mm<10{}0\fi\number\mm}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% DATE
\def\today
{\number\day.\space \ifcase\month\or
January\or
February\or
March\or
April\or
May\or
June\or
July\or
August\or
September\or
October\or
November\or
December\fi,\space
\number\year}


% =======================================================================
% Table setting.
\newcommand{\tabincell}[2]{\begin{tabular}{@{}#1@{}}#2\end{tabular}}
% =======================================================================
% For thick lines in the tables.
\usepackage{tabu}
% Insert the name of "your journal" with
\journalname{Engineering with Computers.}
%
\begin{document}

\title{Procedures to Build Trust in Nonlinear Elastoplastic Integration Algorithm:
Solution and Code Verification
%\thanks{Grants or other notes
%about the article that should go on the front page should be
%placed here. General acknowledgments should be placed at the end of the article.}
}
% \subtitle{Do you have a subtitle?\\ If so, write it here}

%\titlerunning{Short form of title}        % if too long for running head

\author{
 Yuan Feng$^1$ \and
 Kaveh Zamani$^2$ \and
 Han Yang$^1$ \and
 Hexiang Wang$^1$ \and
 % Dragan Kovacevic$^1$ \and
 Fangbo Wang$^1$ \and 
 Boris Jeremi{\'c}$^{1,3}$ 
}

%\authorrunning{Short form of author list} % if too long for running head

\institute{ 
          \ Professor Boris Jeremi\'c
          \email{jeremic@ucdavis.edu}  
          \\
          1 \at
          Department of Civil and Environmental Engineering,
          University of California, Davis, California, U.S.A. 
          \and
          2 \at
          Senior Flood Modeler, 
          Wood Rodgers, Sacramento, California, U.S.A. 
          \and
          3 \at
          Earth and Environmental Sciences Area,
          Lawrence Berkeley National Laboratory, 
          Berkeley, California, U.S.A.
}

\date{Received: date / Accepted: date \\ (The correct dates will be entered by the editor)}



\maketitle


version: \today{} \hhmm{}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\setcounter{secnumdepth}{4}
\setcounter{tocdepth}{4}
\tableofcontents


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\newpage
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\begin{abstract} 

In the last decades, with the development of a number of nonlinear elastic-plastic
integration algorithms,
the correctness or accuracy of the underlying solution  becomes the main concern,
in  both  academia  and  industry. Correctness or accuracy can be estimated (and
improved) using verification.
Verification is one of the main procedures to build trust 
in the numerical modeling of any phenomena. 
A full verification process comprises (a) solution verification (calculation verification) 
and (b) code verification. 

Presented   in   this   paper  are  verification  procedures  for  constitutive,
elastic-plastic integration algorithms, as used in computational nonlinear solid
mechanics. Both explicit and implicit integration algorithms for elastic-plastic
constitutive   equations   are   verified  using  existing  and  developed  new
verification technique. Verification techniques used include prescribed solution
forcing  (PSF)  and Richardson extrapolation (RE). In addition, grid convergence
index  (GCI)  is  applied  to  estimate  the  algorithmic uncertainty during the
integration process.




Verification of elastic-plastic integration algorithms is applied to a number of 
material  models:  from  simple  von-Mises  perfectly plastic  to sophisticated
hyperbolic   Drucker-Prager   with   nonlinear   Armstrong-Frederick  rotational
kinematic  hardening.  Besides,  algorithmic  uncertainty is estimated with both
associative and non-associative material model. In addition caveats and pitfalls
to consider in the code/solution verification processes are deeply discussed.

\keywords{ 
  Elastoplastic Algorithms 
  \and Prescribed Solution Forcing 
  \and Richardson Extrapolation 
  \and Grid Convergence Index (GCI)
}

\end{abstract}

% ====================================================================
% ====================================================================
% ====================================================================
% ====================================================================
% ====================================================================


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


\label{intro}

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


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Modern    computational  modeling systems are very often used in engineering
community for design and assessment.
%\cite{Gousseau2013,Brodtkorb2012,Dinwoodie2015}
Given  all  the  positive  aspects  of  using  computational modeling systems, a
fundamental question is this:
%\cite{Gollan2013}: 
to what extent can we trust  model results \citep{Oberkampf2002,Roy2011}. 
There are a number of engineering failures that can be tied to problems in the numerical codes 
\citep{Jazequel1997,grottke2007} or to misuse of models. 
The  set  of  activities  and  procedures  that  can  be used to  build trust in
numerical modeling results are called verification and validation \citep{Wise2013,Roy2011}. 
According to Roache's definition \citep{Roache1998} verification is 
the response to the following question: "Is the model solving mathematical equations 
correctly?". 
%\cite{Vassiliou2015,Vassiliou2017},
On  the  other  hand,  validation  can be defined through this question: "Is the
model  solving  the  correct  equations?" \citep{Oberkampf2002,Roy2011}. Focus of
this   paper   is   on   verification   of  numerical  algorithms  to  integrate
elastic-plastic  constitutive  equations. Validation  of elastic-plastic
models is not within current scope of this paper and is left for future publications.


A    general    elastic-plastic  constitutive  model  usually  consist  of  four
components:

%
\begin{enumerate}

 \item  {\it  Elastic  law},  that  controls  elastic response before solid
 plastifies. Elastic model/law response can be linear or nonlinear.

 \item   {\it   Yield   criterion   (function)}   that  separates  elastic  from
 elastic-plastic    region.  Yield function is a function of the state of stress
 and  a  set  of  internal  variables. The common yield criteria can be pressure
 independent
 (such   as  von  Mises  and  Tresca)  or  pressure  dependent  (Mohr-Coulomb  or
 Drucker-Prager)  \citep{Drucker1952}. In addition to the common yield criterion,
 a  hyperbolic  yield surface is usually used to smooth the apex
 for higher order derivation.  \citet{Abbo1995}  used this approach for Mohr-Coulomb.

 \item  {\it  Plastic  flow  directions},  that  prescribe  direction  of plastic
 deformation  once material plastifies. 
 Plastic flow directions are a function of the state of stress
 and  a  set  of  internal  variables.
 Plastic flow direction can be associated
 with  the  yield  function (plastic flow is parallel to the normal to the yield
 function/surface  (usually  for  metals)  or  non-associated, when plastic flow
 directions are not normal to yield function/surface (usually for geomaterials).


 \item {\it Hardening laws} that controls evolution of internal variables.
    Hardening laws can control evolution of scalar internal variables (isotropic
    hardening/softening \citep{local-91}), second order tensors (kinematic hardening \citep{Armstrong66})
    and fourth order tensor (distortional hardening \citep{local-99}).

\end{enumerate}


The four components, using different elastic laws, yield functions, plastic flow
directions   and   hardening   laws,   can  combined  to  develop  a  number  of
elastic-plastic  material  models.  In  this  paper  two  models  are  used  for
verification,  namely,  von  Mises and hyperbolic Drucker-Prager elastic-plastic
material  model,  both  with  linear  and  nonlinear  hardening laws. Hyperbolic
Drucker-Prager used here is developed based on the idea from \citet{Abbo1995}.



General numerical integration algorithms for elastic-plastic constitutes equations are 
state-updating procedure.
Commonly used are explicit (forward Euler) and implicit (backward Euler) elastic-plastic integrations
algorithms.
%
%%
%Since the elastic-plastic material is path dependent, a simultaneous update of 
%the internal variables is critical in the iteration process. 


The paper is structured as follows.  
Section~\ref{sec_Background} introduces the history of verification from computational
fluid dynamics to nonlinear finite element analysis of solid mechanics.
Section~\ref{sec_intro_elastic-plastic} shortly introduces the elastic-plastic algorithms. 
Section~\ref{sec_Verification_theory} presents the theoretical basis of the applied 
verification techniques.
Section~\ref{sec_elastic-plastic_material_model} briefly reviewed the verified 
elastic-plastic material model. 
Section~\ref{sec_verification_example} shows the verification examples of different
elastic-plastic material models.
% Section~\ref{sec_pitfall_verification} discusses the pitfalls and concerns 
% in the verification process.



% ====================================================================
% ====================================================================
% ====================================================================
% ====================================================================
% ====================================================================
% ====================================================================



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Background}

\label{sec_Background}


Computer  speed  is  increasing at a steady pace, following Moore's observation,
that is sometimes called Moore's law \citep{Moore1965}
that states that 
''The  number of transistors in a dense integrated circuit doubles approximately
every two years'', leading to ever faster computers.
In addition, software 
developments,  in  last number of decades, have made it possible to develop very
sophisticated  mechanics designs  for  a  large  number  (almost  all)  of  recent  civil
engineering  objects. Use of computational modeling in the area of mechanics in
other fields of engineering is also noted, particularly in mechanical, aerospace
and biomedical engineering.
It   is  also  worth  noting  Wirth's  law \citep{Wirth1995},  that states  that
''Software  is getting slower more rapidly than hardware is getting faster'', is
almost  balancing  hardware  speed  advancements. Nevertheless, sophistication in
modeling   capabilities   in   computational  mechanics  (and  other  fields  of
computational  modeling) is bringing forward significant number of advantages, with
economy and safety being the main beneficiaries.
% 
% In the ending years of the 20th century, as computers' speed increased with 
% astonishing pace and advanced numerical algorithms were developed, 
% computational modeling turned out to be a popular tool for analysis, 
% design,  and  research  (noted are few example references from a large number of
% references  describing  use  of  computational  tools  in  design and assessment
% \citep{Maufroy2015,Burkhart2013,Bao96}).
% 



Despite  significant advances in modeling and simulations, there exists a number
of  cases  where  (problematic)  results  from nonlinear analysis contributed to
engineering incidents or even failures \citep{Hatton1994,Hatton97,Selby1997}. 
%
% However, the "computational modeling" results are typically looked upon with 
% skepticism compared to its expensive counterpart: "physical modeling". 
% This lack of confidence in computational results are not surprising as 
% many disasters happened as direct result of bugs in the computational codes 
% or incorrect use of codes, for example: Statoil's offshore oil rig, 
% Sleipner A, was sank as the consequence of an incorrect stress analysis 
% that had been done with NASTRAN, a well-established finite element method 
% computer package \citep{selby1997}. 
% In another catastrophe in academia, a software bug destroyed the shining 
% career of an outstanding and well-decorated researcher \citep{miller2006}. 
% Those issues emerged a very important question: How can we rely on the 
% results of computational modeling? 

The  main  question  that  is frequently asked, particularly for nonlinear
analysis  results,  is  this:  ''how  much  can  we trust results of the
numerical analysis for use in design and assessment of infrastructure objects?''.
This question can be split in two questions:
%      
%
\begin{enumerate}

\item How much can (should) we trust  model implementations?

\item How much can (should) we trust mathematical models?

\end{enumerate}
%
Answer  to  both  questions  can  be found out through a process of verification
(question \#1) and validation (question \#2).
Field of the Verification and Validation (V\&V) has received significant interest in
view of ever growing use of numerical modeling in design and assessment
\citep{Roach1998,Oberkampf2002,Roy2011,Oden2004,Babuska2004,Oden2010a,Oden2010b}.



V\&V provide a rigorous framework to assessing 
the accuracy computational simulations. 
Verification has two parts: code verification and solution verification.
Code verification (solver verification) is a procedure to ensure 
that a code is bug free. 
Code verification is a one-time task and it should be only repeated 
when the source code is modified. 
Solution verification is a procedure to check 
if the simulation results are accurate.

%
% 
% 
% It is task of the modeler and it should have performed on each 
% simulation \citep{Roache1998,Oberkampf2010}. 



The simulation verification techniques were developed primarily by the 
U.S. Department of Energy (DOE) after adaptation of Comprehensive 
Nuclear-Test-Ban Treaty (CTBT) in 1996 \citep{Oberkampf2010}. 
Both code verification and calculation verification studies are based on mesh 
convergence test. 
%
%
In this test, behavior of numerical error is studied as mesh size shrinks. 
There are a number of methods to perform mesh convergence test: 
and Richardson Extrapolation \citep{Roache1993},
Method of Manufactured Solution (MMS) \citep{Roache2002} and
Method of Exact solutions (MES) \citep{Zamani2014}. 
Other methods of code/solution verification were suggested 
in the literature in order to overcome practical difficulties or shortcomings 
of the above mentioned methods: 
Prescribed Solution Forcing Method \citep{dee1991},  
Method of Nearby Solution \citep{roy2007}, 
and External Verification Analysis \citep{ingraham2013}. 
%
It is worth mentioning that there is no consensus on the naming of verification 
methods  in  the  verification literature, so different names will be used, for,
potentially  similar procedures, and full reference will be provided to original
development.



Application  of numerical verification techniques started in computational fluid
mechanics and eventually founds its way into
computational solid mechanics \citep{Kamojjala2013}. 
Although the fundamental concepts are similar, there are inherent differences 
in implementation of the verification methods. 
To the best of the authors' knowledge, 
full implementation of numerical verification test in computational 
geomechanics is scarce \citep{Sjogreen2011}. 
%In this paper, we introduce and implement methods for numerical verification 
%and uncertainty quantification under the arch of computational mechanics 
%for the first time. 
Described in this paper are code and solution verification techniques 
for elastic-plastic constitutive integration algorithms.
%
Additionally, implemented and described is Roache Grid Convergence Index \citep{Roache1997} test
to uncover numerical uncertainty in the elastic-plastic simulations.   






% ====================================================================
% ====================================================================
% ====================================================================
% ====================================================================
% ====================================================================
% ====================================================================

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Elastoplastic Constitutive Integration Algorithm}

\label{sec_intro_elastic-plastic}


Elastoplastic integration algorithm is used to solve for the  stress increment
from a strain increment.
The elastic-plastic integration algorithm is a two-step algorithm including the elastic predictor
and plastic corrector \citep{local-64}. 
The elastic predictor step is similar for different algorithms. 
Depending  on  the  direction  of  the  plastic  corrector,  the elastic-plastic
algorithms  are  generally  divided  into  explicit, forward Euler and implicit,
backward Euler algorithms.
Both  explicit  and implicit algorithms are special cases of a general mid-point
algorithm   \citep{local-18,local-64,Crisfield1997}.   These   two  constitutive
algorithms represent the most common constitutive integration algorithms.
Explicit algorithm uses the starting elastic-plastic point to 
define the direction of the plastic corrector.
Implicit algorithm uses the  final stress state to 
define the direction of the plastic corrector. 
Iterations  are  required  for the implicit algorithm in order to find the final
stress state.


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

Small        deformation,       incremental       elasto-plastic       equations
(\ref{eq:03_constitutive}), given below:
%
\begin{equation}
\begin{aligned}
  d \sigma_{ij}     & = E_{ijkl} \ (d \epsilon_{kl} - d \epsilon_{kl}^{pl}) \\
%  \nonumber \\
  d \epsilon_{ij}^{pl} & =  d \lambda \; m_{ij}(\sigma, q_*) \\
%  \nonumber \\      
  d q_*             & =  d \lambda  \; h_{*}(\sigma, q_*) \\
%  \nonumber \\
  F(\sigma, q_*) & = 0 \; \wedge \;  d F(\sigma, q_*) = 0
\end{aligned}
\label{eq:03_constitutive}
\end{equation}
%
%
%
are used to
solve  for an increment  of  stress ($d \sigma_{ij} $), an increment of internal variables ($d q_{*}$)
and for elastic-plastic stiffness tensor 
($E_{ijkl}^{el-pl}$)  at
a material point
(integration point, Gauss point) for a given increment  in  strain 
($d \epsilon_{kl}$). 
Small deformaton assumption is used so that  an increment of strain ($d \epsilon_{ij}$)
can be defined from increments in displacements ($d u_{i}$) as
$d \epsilon_{ij} = 1/2 ( d u_{i,j} + d u_{j,i})$.
Small  deformation assumption allows for additive split of strain increment into
elastic  and  plastic  parts  
($d  \epsilon_{kl}  =  d  \epsilon_{kl}^{el}+d\epsilon_{kl}^{pl}$).
%
%
Internal variables ($q_{*}$) can be (a) scalars, for isotropic hardening/softening, 
(b) second-order tensors, for kinematic hardening/softening or (c) fourth order
tensors, for distortional hardening/softening.
%Increment in strain  ($d {\epsilon}_{ij}$) is
%additively  decomposed small deformation assumption into elastic strain  
%($d \epsilon_{kl}^{el}$) and plastic strain ($d \epsilon_{kl}^{pl}$).
%
The scalar $d \lambda $ is the plastic multiplier. Plastic flow directions
are  defined  by the tensor $m_{ij}$, and hardening/softening rule is defined by
the tensor function $h_*(\sigma, q_*)$.
The yield function $F(\sigma, q_*)$ separates elastic response of material  for
($F(\sigma, q_*)<0$) from elastic-plastic response of material ($F(\sigma, q_*)=0$.
It is important to note that yield function $F$ cannot be larger than zero 
$F(\sigma,  q_*)  \not>  0$,  hence  $d  F(\sigma,  q_*) = 0$. Numerically, this
requirement  of  non-positive yield function ($d F(\sigma, q_*) \not> 0$) cannot
be  enforced  for  explicit  integration  algorithm, as there are no equilibrium
iterations,  and  the  stress  solution  in  general  drifts away from the yield
surface into plastic region.


%Above  equation can be solved using a general set of mid-point algorithms, where
%explicit  (forward  Euler)  and  implicit  (backward  Euler)  are  the  two most
%prominent and used examples.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Explicit, Forward Euler Algorithm}

\label{subsec:Forward_euler}


An    explicit,    forward  Euler  algorithm  uses  plastic  flow directions and
hardening/softening  function  at  the initial elastic-plastic state, that is on
the  yield  surface,  to calculate the increment in stress ($d \sigma_{ij}$) and
continuum, tangent elastic-plastic stiffness tensor ${}^{cont} E_{ijkl}^{ep}$ 
\citep{local-64,Crisfield1997,Chen94}.
%
The plastic-multiplier 
$d \lambda$ is solved for as
%
\begin{equation}
    d \lambda = \frac{ {}^{cross} n_{ij} E_{ijkl} d\epsilon_{kl} }
                {{}^{cross} n_{ab} E_{abcd} {}^{cross} m_{cd}  
                - {}^{cross} \xi_{*} h_{*}}
\end{equation}
%
%
%
where the tensor ${}^{cross} n_{ij} ={\partial F}/{\partial \sigma_{ij}} $ is the 
normal to the yield surface at the starting point, while 
the parameter ${}^{cross} \xi_* = {\partial F}/{\partial q_{*}}  $ is the derivative
with respect to the corresponding internal variables, 
both  for a stress --  internal variable state  on the yield surface,
%Also, the notation \textit{cross} represents that the starting pointing is 
%the cross point between the elastic predictor and the current yield surface. 
Continuum, tangent stiffness tensor is obtained as 
%
%
\begin{equation}
    {}^{cont} E_{ijkl}^{ep} =
    E_{ijkl} -
    \frac{  E_{ijab} {}^{cross} m_{ab} {}^{cross} n_{cd} E_{cdkl}  }
    {
    {}^{cross} n_{ot} E_{otrs}  {}^{cross}  m_{rs}
     -    \xi_{*}   h_{*}
    }
\end{equation}
%
%




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Implicit, Backward Euler Algorithm}

\label{subsec:Backward_euler}


An  implicit,  backward  Euler  algorithm  uses  state of  stress and  internal 
variables  at  the  final step (solution) of particular increment to develop the
consistent, algorithmic stiffness tensor ${}^{cons} E_{ijkl}^{ep}$ \citep{local-19,local-4}.
%
Since  the state of stress $\sigma_{ij}$ and internal variables $q_{*}$ at the final,
solution state is not known before the calculation, 
the implicit algorithm requires iterations to
converge to the final results.

In the implicit algorithm, the plastic-multiplier 
$d \lambda$ is iteratively solved for \citep{Jeremic94} as
%
%
\begin{equation}
    d\lambda = \frac{
                {}^{next}F 
                - 
                {}^{next} n_{kl} r_{ij} {}^{next} T^{-1}_{ijkl} 
                }
                {
                {}^{next} n_{ot} E_{abcd} {}^{next} H_{cd} {}^{next} T^{-1}_{abot} 
                -
                \xi_* h_*
                }
\end{equation}
%
%
where $r_{ij} = \sigma_{ij} - ( {}^{pred} \sigma_{ij} - \lambda E_{ijkl} m_{kl} ) $
is the stress residual from the last iteration,
and the tensors $T_{ijkl}$ and $H_{ij}$ are  given as
%
%
\begin{equation}
    \begin{aligned}
        T_{ijkl} &= \delta_{ik}\delta_{jl} + 
                    \lambda E_{ijab} \frac{\partial m_{ab}}{\partial \sigma_{kl}} \\
        H_{ij}   &= m_{ij} + 
                    \lambda \frac{\partial m_{ij}}{\partial q_*} h_* \\
    \end{aligned}
\end{equation}
%
%
Iterations continue until predetermined tolerance for the stress residual ($r_{ij}$) and the yield
function ($F$) are reached.

The  consistent\footnote{Consistent  with  the  Newton  iterative  algorithm  on
global,  finite element level.}, algorithmic elastic-plastic stiffness tensor is
then obtained as
%
%
\begin{equation}
    {}^{cons} E_{ijkl}^{ep} =
    E_{ijkl} -
    \frac{  R_{ijab} \  {}^{next} H_{ab} \  {}^{next} n_{cd} R_{cdkl}  }
    {
    {}^{next} n_{ot} R_{otrs}  \  {}^{next}  H_{rs}
     -    \xi_{*}   h_{*}
    }
\end{equation}
%
%
where the tensors $R_{ijab}$, $T_{ijmn}$, and ${}^{n+1} \! H_{kl}$  are defined as
%
%
\begin{samepage}
\begin{eqnarray*}
R_{ijab} &=& T_{ijkl}^{-1} E_{klab} 
\\
T_{ijmn} &=&
\delta_{im} \delta_{nj} +
\left. %just because of \right|
d \lambda E_{ijkl} \frac{ \partial m_{kl} }{\partial \sigma_{mn}}
\right|_{n+1}
\\
{}^{n+1} \! H_{kl} &=&
{}^{n+1} \!m_{kl} 
+
\left.
\frac{ \partial m_{kl} }{\partial q_{*}}
\right|_{n+1} h_{*}
\end{eqnarray*}
\end{samepage}
%







% ====================================================================
% ====================================================================
% ====================================================================
% ====================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Theoretical Basis for Verification}

\label{sec_Verification_theory}




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{The Asymptotic Regime of Convergence}


For  numerical solutions  that rely on discretization in space 
%(mesh size) 
or in
time,  
%(step size), 
mesh size or time step size choice is crucial for accuracy of
the solution.
%
When the mesh size or step size $\triangle x$ is large, 
the error oscillates since the discretization 
size is too large to represent the actual continua (space or time) 
in the discretized domain. 
As  mesh size --  step size ($\triangle x$) is reduced to the middle size region, 
the error starts to decrease asymptotically.
This region is called the asymptotic regime of convergence. 
However, as the mesh size --  step size ($\triangle x$) is reduced even further,
the  round-off  error contributes to a rapid increase in solution error. This is
illustrated in Figure~\ref{fig_asympt_err_zone}.

\begin{figure}[!htbp]
 \centering
 \captionsetup{justification=centering}
 \includegraphics[width = 10cm]{../Figure-files/assympt_err_zoons_constant.PNG}
 \caption{
 \label{fig_asympt_err_zone}
  Behavior of solution error in mesh size -- step size convergence study: 
    I) over resolved, II) appropriate size, and III) under resolved.}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Richardson Extrapolation Technique as a Verification Tool}


Richardson extrapolation technique was developed by
\citet{Richardson1911}. 
Recently,  Richardson  extrapolation  was     used     for     code/solution    
verification   of PDEs solvers \citep{Karatekin1997,Xing2010} when the Method of
Exact solutions (MES) is not applicable \citep{Roache1993}.





The accuracy of the elastic-plastic algorithm is controlled to a large extent by
the size of the strain
increment.  This  is particularly true for the explicit algorithm, while for the
implicit algorithm, size of tolerance also plays important role.



Consider for example the following experiments.
Assume an initial stress at $\sigma_{0}$. 
Assume also that the exact stress solution is $\sigma^{*}$. 
%
To achieve this exact stress solution, a large number of small steps 
is applied to perform the integration. The accuracy of the final stress
$\sigma$, thus achieved, is controlled by the size of strain increment $d\epsilon$. 
Ideally, the final stress $\sigma$ is accurate, that is $\sigma \to \sigma^{*}$, 
when $d\epsilon \to 0 $ \citep{Lax1956}. 
The Taylor series expansion of $\sigma$ with respect to $d\epsilon$ is 
%
\begin{equation}
  \sigma(d\epsilon)  = \sigma^{*} + C d\epsilon^{\beta} + O(d\epsilon^{\beta+1})
  \label{eq_stress_Taylor_wrt_strain}
\end{equation}
%
where $C$ is constant and $O$ is the big-O notation, which represents the higher order
errors. 

In Richardson Extrapolation technique, a new function is defined based on the stress
integration algorithm. 
%
\begin{equation}
  R(d\epsilon, k) =  
   \frac{k^{\beta} \sigma(d \epsilon) - \sigma (k d \epsilon ) }
   { k^{\beta} -1 }
  \label{eq_richardson_original}
\end{equation}
%
Here $k$ is a constant that defines a number of increments that are used to 
create Richardson Extrapolation. For example, in the above equation, 
$R(d\epsilon,    k)$  is  a  difference  between a solution ($\sigma(d \epsilon)$)
using single increment
$d\epsilon$, multiplied by $k^{\beta}$ and a solution  ($\sigma(k d \epsilon)$) 
obtained using $k d \epsilon$ increment, and divided by $k^{\beta} -1$. 
In addition, parameter $\beta$   is   the  accuracy order that  is  solved
for from Richardson Extrapolation tests.


Then,      by     substituting     Eq~\ref{eq_stress_Taylor_wrt_strain}     into
Eq~\ref{eq_richardson_original}, we obtain
%
\begin{equation}
  \begin{aligned}
     R(d\epsilon, k) 
     = & \frac{k^{\beta} ( \sigma^{*} + C d\epsilon^{\beta} + O(d\epsilon^{\beta+1})) }{ k^{\beta} -1} 
    - \frac{ \sigma^{*} + C k^{\beta} d\epsilon^{\beta} + O(d\epsilon^{\beta+1})  }{ k^{\beta} -1} \\
     =& \sigma^{*} + O(d\epsilon^{n+1})
  \end{aligned}
  \label{eq_richardson_expansion}
\end{equation}
%

Function R($d\epsilon$, k) is the Richardson Extrapolation of $\sigma (d\epsilon) $. 
By     comparing      the     final     Richardson    equation    in
Eq~\ref{eq_richardson_expansion} to the original
Eq~\ref{eq_stress_Taylor_wrt_strain},  it is noted that the discretization error
$C d\epsilon^{\beta} $ is cancelled out.
It  is  also  noted  that by repeating  Richardson equation one more time, one is
able to cancel out the next higher
order  discretization error. This can be continued any number of times to cancel
higher order discretization errors, however there is a tradeoff with efficiency.

The     remaining  problem  is  to  determine  the accuracy order  parameter 
$\beta$        for        the        elastic-plastic        algorithms        in
Eq~\ref{eq_stress_Taylor_wrt_strain}.


% =================================================================================
% =================================================================================
% =================================================================================
% =================================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Uncertainties in Numerical Simulation}



The    uncertainty    of    any   computational   simulation   is   defined   as
\citep{ASME_Vand_V_2009}:
%
\begin{equation}
\delta_{simulation} = (\delta_{model}+\delta_{numerical}+\delta_{input}) - \delta_{experimental}
\label{eq_error_in_modeling}
\end{equation}
%
%
where  $\delta_{simulation} $  is  the  total  uncertainty  in  the  simulation. The terms in 
parenthesis  in  the right hand side Eq.~\ref{eq_error_in_modeling} are the best
effort  of  a numerical modeler and include: 

\begin{itemize} 

  \item Modeling uncertainty, $\delta_{model}$,
  
  \item  Errors/uncertainties  due  to  numerical solution
  schemes (iterative error, round off error, discretization error, etc)  $\delta_{numerical}$,
  
  
  \item Parametric uncertainty, $\delta_{input}$   is  the  errors  that  propagate into the results of
  numerical    simulation    from  uncertainty  in  input   parameters,   including   geometry,
  initial/boundary conditions and material properties, 
  
\end{itemize} 

Experimental  uncertainty,   $\delta_{experimental}$   is   the  uncertainty  in
the experimental measurement of the phenomena.



Reduction of each of the above-mentioned uncertainties will increase the quality
of  the  simulation  results. 

Experimental uncertainty  $\delta_{experimental}$ is the problem of
interest      in      experimental   mechanics  and  is  related  to  validation
\citep{Oberkampf2002}, modeling uncertainty $\delta_{model}$ is a problem of model sophistication,
and  $\delta_{input}$  is  the  problem  of  uncertain  material  parameters and
uncertain loads, and their propagation through the system \cite{Sett2010a}.
Thus here we only
discuss  the  procedure to develop an upper and lower bound on the "numerical"
uncertainty  $\delta_{numerical}$.



The  numerical  uncertainty has its  roots  in the truncation error, round-off error and
iterative error.
%
First,  truncation  error  represents  the  error  due  to  solving  discretized
equations instead of continuum counterpart.
%
Second,   round-off  error  is an accumulation of errors due to finite precision
calculations, machine epsilon.
This error is not  that important  in modern computational mechanics simulations as it is 
usually overshadowed by other numerical errors.
% %% %
% %This  is  the  consequence  of using a fixed number of digits on the computer to
% %represent the decimal number.
%
Finally,  the  iterative  error  is  due to the limited number of iterations for
solving either an algebraic system of equations or a system of discretized PDEs.
%
The  round-off  errors  and  the  iterative  errors are usually considered to be
smaller by orders of magnitude compared to truncation error \citep{Oberkampf2002c,Saad2003}.
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Grid Convergence Index}

\label{sec_estimate_numerical_uncertainty}


Grid Convergence Index (GCI) is an uncertainty quantification technique 
firstly used in the computational fluid dynamics (CFD) \citep{ASME_Vand_V_2009}. 
At least three groups of incremental sizes are required to compute the GCI. 
The procedures of the uncertainty estimation are as follows. 

\begin{itemize}

\item
Choose three significantly different incremental sizes $h$ and run the simulation. 
The increment refinement factor, $r = h_{coarse}/h_{fine} $, is usually greater than 1.3
based on experience \citep{Roache1998}. 

\item
Calculate the observed order $\beta$. 
Let $h_1 < h_2 < h_3 $ and $ r_{21} = h_2 / h_1 $, $r_{32}= h_3 / h_2 $. 
%
\begin{equation} 
  \begin{aligned}
    & \beta  =  \frac{\ln |\frac{d_{32}}{d_{21}}| + \varphi(\beta) }{\ln(r_{21})}  \\
    & \varphi(\beta)  = \ln( \frac{ r_{21}^{\beta} - t }{  r_{32}^{\beta} - t }  )  \\
    & t  = sign( \frac{d_{32}}{d_{21}} )  \\
  \end{aligned}
\end{equation}
%

where $d_{32} = \omega_3 - \omega_2 $,  $d_{21} = \omega_2 - \omega_1 $ and 
$\omega_n $ denotes the component of the stress result in the $n^{th}$ incremental size.





The relative error $e$, as a dimensionless form, is 
%
\begin{equation}
  e_{21} = \left| \frac{ \omega_{1} - \omega_{2} }{\omega_{1}} \right|. 
\end{equation}
%

\item 
The GCI is 
%
\begin{equation}
  GCI_{21} = \frac{Fs \cdot e_{21} }{r_{21}^{\beta} - 1}
\end{equation}
%

where $Fs$ is the factor of safety,
and $GCI_{21} $ represents the numerical uncertainty at the 
strain incremental size $h_2$ with reference to 
the more accurate results at the strain incremental size $h_1$.
Based  on experience from computational fluid dynamics,
a factor of safety $Fs=1.25$ \citep{ASME_Vand_V_2009}.

\end{itemize}


It follows that the uncertainty caused by the numerical error $U_{num} $ is 
equivalent to  $GCI_{21}$. 
The range of uncertainty $\pm U_{num} $ is bounding  
the exact mathematical solution with a 95\% confidence or, in other words, it is
within
two  standard  deviations  ($\pm  \sigma$), It is noted that notation for 
two standard deviations
$\pm 2\sigma$ represent the probability level (95.4\%), and not the stress state.


% % additional for  CHECK meaning of uncertainty!
% % Explanation:
% \begin{tcolorbox}
% The form of the GCI is based on theory of Richardson extrapolation, but the use
% of absolute values for estimated errors and the factor Fs are based on 
% empiricism involving the examination of several hundred (CFD) case studies. 
% %
% 
% The empirical tests involved the determination of conservatism in 95\% of
% the cases, corresponding to (dimensional) GCI = $U_{num}$ at 95\% confidence. 
% %
% 
% No assumptions on the form of the error distributions were made 
% nor were necessary for these empirical studies, 
% since actual data was examined with a simple pass/fail criterion. 
% %
% Specifically, the common statistical assumption of a Gaussian distribution
% was not used. 
% %
% 
% To agree with the new international standard use of one standard deviation $\sigma$.
% 
% \end{tcolorbox}
% 
 


% =================================================================================
% =================================================================================
% =================================================================================
% =================================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Accuracy order of Elastoplastic Algorithm}

\label{sec_accuracy_order_of_elastic-plastic}

Before    the  implementation  of Richardson Extrapolation, the accuracy   order
($\beta$)  of  the  elastic-plastic  algorithm  needs to be determined. 
Both the analytical determination and the numerical estimation of accuracy order  for 
elastic-plastic algorithms  are discussed in this
section.  Analytical  determination  of  accuracy  order  is possible for simple
material  models.  On  the  other  hand,  for more sophisticated material models,
it is only possible to develop numerical estimation of accuracy order.
% =================================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Analytical Determination of Accuracy Order}
\label{subsec_analysis_accuracy_order}


An analytical study of accuracy order for
elastic perfectly plastic constitutive integration algorithms was developed by 
\citet{local-18}.   
In      their      study,   \citet{local-18}  proved  that   all elastic plastic
algorithms  within  generalized  trapezoidal  and  generalized  midpoint  return
mapping are at least first-order accurate.
In    the  case  of  Crank-Nicolson method, where plastic corrector direction is
chosen  as an average of explicit (forward Euler) and implicit (backward Euler),
elastic  perfectly  plastic  constitutive  integration algorithm is second order
accurate.
%
%
%When the center direction is selected in the trapezoidal return mapping or 
%the center stress is selected in the midpoint return mapping, 
%the elastic-plastic algorithm meets the second-order accuracy. 

\noindent
{\bf Remark:} {Conclusions   on  algorithm  accuracy order 
made by \citet{local-18} do apply to elastic perfectly plastic material models.
%
However,        it    is  important to note that when material model features
hardening  and/or  softening  response  (isotropic  or  kinematic  hardening and
softening) the accuracy order is unclear.
% 
In  those  cases,  numerical  estimation  of  the accuracy order ($\beta$) is
needed.}

% =================================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Numerical Estimation of Accuracy Order}

Another approach to determine the accuracy order is 
through an estimation using numerical experiments on various 
sizes of increments. 
The idea is to apply the Richardson Extrapolation on 
strain increment sizes. 
At least three different strain increment sizes are required 
to estimate the observed order. 

Specifically,    using    three  different    strain    increment    sizes   (a)
$d\epsilon$,  (b)  $d\epsilon/a$,  and (c) $d\epsilon/b$, where the scale factor
$b>a>1.3$,  the  expression  of  approximation  of  exact solution $\sigma^*$ is
written as:
%
\begin{equation}
  \begin{aligned}
    \sigma^{*} 
    & = \frac{a^{\beta} \sigma(\displaystyle{\frac{d\epsilon}{a}}) - \sigma(d\epsilon) }{a^{\beta} -1} 
         + O(d\epsilon^{n+1}) \\
    & = \frac{b^{\beta} \sigma(\displaystyle{\frac{d\epsilon}{b}}) - \sigma(d\epsilon)  }{b^{\beta} -1} 
      + O(d\epsilon^{n+1}) \\
  \end{aligned}
  \label{eq_estimate_accuracy_order}
\end{equation}
%

Assuming the higher order errors [$O(d\epsilon^{n+1})$] are negligible 
we can rewrite Eq.~\ref{eq_estimate_accuracy_order} as
%
\begin{equation}
%  \begin{aligned}
%    & \sigma(\frac{d\epsilon}{a}) + \frac{ \sigma(\displaystyle{\frac{d\epsilon}{a}}) - \sigma(d\epsilon)  }{a^{\beta} -1} \\ 
%  = & \sigma(\frac{d\epsilon}{b}) + \frac{ \sigma(\displaystyle{\frac{d\epsilon}{b}}) - \sigma(d\epsilon)  }{b^{\beta} -1} \\ 
%  \end{aligned}
    \sigma(\frac{d\epsilon}{a}) + \frac{ \sigma(\displaystyle{\frac{d\epsilon}{a}}) - \sigma(d\epsilon)  }{a^{\beta} -1} \\ 
  = \sigma(\frac{d\epsilon}{b}) + \frac{ \sigma(\displaystyle{\frac{d\epsilon}{b}}) - \sigma(d\epsilon)  }{b^{\beta} -1} \\ 
  \label{eq_richardson_estimate_final}
\end{equation}
%

It is noted that in  Eq.~\ref{eq_richardson_estimate_final}, stress solutions
$\sigma(d \epsilon)$, 
$\sigma(d\epsilon/a) $, 
$\sigma(d\epsilon/b)$
are measured, observed from the three numerical experiments. 
The  only  unknown  in  Eq.~\ref{eq_richardson_estimate_final}    is  thus  the 
accuracy order $\beta$, which can be easily solved for.









% =================================================================================
% =================================================================================
% ====================================================================
% ====================================================================
% ====================================================================
% ====================================================================
% ====================================================================

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Elastoplastic Material Models}
\label{sec_elastic-plastic_material_model}


Elastoplastic material models used for numerical experiments are based on von
Mises and hyperbolic Drucker-Prager yield surfaces.
%
Material  models of von-Mises type are used to simulate the pressure-independent
material behavior.
%
On the other hand, hyperbolic Drucker-Prager family of models are used to simulate
pressure-dependent material behavior.


For modeling pressure independent behavior of material, plastic volume change is
usually    negligible,    which    leads  to  plastic  flow  directions that are
perpendicular/normal to
the von-Mises yield surface.
%
This is the so called  associated plasticity behavior.
%
%
For  modeling pressure dependent material response, plastic volume change can be
significant,    however    plastic    flow    directions    are    usually   not
perpendicular/normal to the yield surface.
%
This is the so called  non-associated plasticity behavior.
%
%
In  addition,  different  types  of hardening and softening of yield surface and
plastic  flow  directions exist \citep{Chen94}. 
%
For example perfectly plastic material will not
feature any evolution of yield surface.
%
On  the  other  hand  isotropic hardening and/or softening will change the yield
surface    size.
%
For    kinematic  hardening,  yield  surface  will move (translate or rotate) in
stress space.
%
Distortional     hardening \citep{local-99} where  yield  surface changes shape,
is also possible , but it is rarely used.
%
%
Plastic  flow  directions  are  usually developed either based on associated or
non-associated,  or,  for  more  sophisticated  model,  directly as functions of
stress and internal variable space \citep{Dafalias82,local-63,Manzari97}.



Presented below are common examples of yield surface functions for von-Mises and
Drucker-Prager  models, and examples of non-associated plastic flow as well as a
nonlinear kinematic hardening rule.





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{The von-Mises Yield Surface}

The  von  Mises  yield  surface \citep{local-63}, a cylinder in stress space, is
given as
%
\begin{equation}
F (\sigma_{ij}, k, \alpha_{ij})
= \sqrt{\frac{3}{2}(s_{ij} - \alpha_{ij})(s_{ij} - \alpha_{ij}) } - k = 0
\label{VMYF2}
\end{equation}
%
%%%%
where  $s_{ij}$ is the deviatoric stress ($s_{ij} = \sigma_{ij} - 1/3\sigma_{kk}
\delta_{ij}$).
%
The yield surface radius $k$ is the scalar internal variable describing the size
of yield surface.
%
The back stress $\alpha_{ij}$ represents location of yield surface in the stress space.
%
Both    yield   surface  radius  $k$ and the back stress $\alpha_{ij}$ represent
internal  variables  (generally  annotated as $q_{\ast}$) that can harden and/or
soften.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Hyperbolic Drucker-Prager Yield Surface}


The Drucker-Prager yield criterion is a pressure-dependent material model, where
yielding  of  materials  is controlled by the deviatoric stress, and also by the
confining pressure. 
%
The original Drucker-Prager yield surface \citep{local-63} 
has  a  singularity  point at the location of zero deviatoric stress. Hyperbolic
yield surface \citep{Abbo1995} modifies
the Drucker-Prager yield surface by removing the singular point, 
% where the original yield line becomes the asymptotic line of 
% the hyperbolic function,
% as shown in Figure~\ref{fig_hyperbolic_function_diagram1}.
% A typical hyperbolic function is 
% %
% \begin{equation}
%   \left(\frac{x-d}{a}\right)^2 - \left(\frac{y}{b}\right)^2 = 1
% \end{equation}
% %
% Extend the yield surface to include isotropic and kinematic hardening.
The hyperbolic Drucker-Prager yield function is given as
%
\begin{equation}
  F (\sigma, \eta, \alpha_{ij}) =
         \sqrt{\frac{1}{2}(s_{ij} - p \alpha_{ij})(s_{ij} - p \alpha_{ij})  
         + a^2 \eta^2 } - \eta p - \xi
\end{equation}
%
%
where  $\eta$  is  controls the friction angle, and $\xi$ controls the
cohesion,  while   $p \alpha_{ij}$  controls rotation of the yield surface. Back
stress  $\alpha_{ij}$ is multiplied by mean stress $p=\sigma_{ii}/3$ in order to
force the surface to rotate.
 Variable $a$ represents the rounded length between the cut-off
on    $p-axis$    and    the    original    apex    points,    as    shown    in
Figure~\ref{fig_hyperbolic_function_diagram1}.
%
%
%
%Figure~\ref{fig_hyperbolic_function_diagram1} shows trace of Drucker-Prager yield surface in stress space.
%
\begin{figure}[!htbp]
  \centering
  \includegraphics[width = 7cm]{../Figure-files/hyperbolic_dp.pdf}
  \caption{The hyperbolic Drucker-Prager yield function.}
  \label{fig_hyperbolic_function_diagram1}
\end{figure}



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Non-associative Plastic Flow}


In pressure-dependent materials, the dilative or contractive behavior 
is usually not associated
with the normal to the yield surface $n_{ij}=\partial F/ \partial \sigma_{ij}$. 
The non-associative plastic flow direction $m_{ij} $ is expressed as 
%
\begin{equation}
m_{ij}=
n_{ij}
- 
\frac{1}{3}D\delta_{ij}
\end{equation}
%
%
where dilatancy parameter D is defined as
%
\begin{equation}
D = \zeta \left( \sqrt{\dfrac{2}{3}}k_d - \frac{\sqrt{s_{mn}s_{mn}}}{p} \right)
\label{eq_dilatancy_D}
\end{equation}
%
%
The volumetric plastic strain rate $\zeta$ governs the amplitude of plastic volume changes. 
The term ${\sqrt{s_{mn}s_{mn}}}/{p}$ represents the stress-ratio. 
The material constants $k_d$ governs the critical stress-ratio, 
which controls the direction of the plastic flow.
The dilative or contractive behavior is controlled by the stress-ratio 
(${\sqrt{s_{mn}s_{mn}}}/{p}$)
\citep{local-91}.
When  dilatancy  parameter  $D<0$  plastic  compaction takes place, whereas when
dilatancy parameter $D>0$ plastic
dilatancy takes place. 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Armstrong-Frederick Hardening Law}


Nonlinear    kinematic    hardening    law   developed  by  \citet{Armstrong66}
simulates  nonlinear  behavior  of  materials through control of the back stress
($\alpha_{ij}$)
%
\begin{equation}
  d {\alpha }_{ij} = 
  \frac{2}{3}\;  h_a \; (d \epsilon_{ij}^p)^{dev} 
  -
  c_r \alpha_{ij} \sqrt{\frac{2}{3}\;  (d \epsilon_{rs}^p)^{dev} \; \; (d \epsilon_{rs}^p)^{dev}} 
\end{equation}
%
where 
$d {\alpha_{ij}}$ is increment of the back-stress, 
$h_a$ is the initial rate of change of back-stress,
$h_a/c_r$ controls the limit, assymptote of the back-stress, 
and $(d \epsilon_{ij}^p)^{dev}$   is the deviatoric component of increment of 
plastic strain.







% =================================================================================
% =================================================================================
% =================================================================================
% =================================================================================
% =================================================================================
% =================================================================================
% =================================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Error Maps for Elastoplastic Constitutive Integration}

\label{sec_verification_example}

Verification examples in Sections~\ref{subsec_example_vonMises_yf},
~\ref{subsec_example_DruckerPrager_linearHarding}, and 
\ref{subsec_example_DruckerPrager_nonlinear_hardening} present  
error  maps  in  octahedral  stress  plane,  for different stress step size, 
different  constitutive  integration  algorithms, for von-Mises elastic perfectly
plastic material model.
%
Next two examples, Sections~\ref{subsec_example_associative},
and \ref{subsec_example_non_associative} quantitatively evaluated the accuracy of
forward and backward Euler algorithms for linear kinematic hardening, associated and non-associated
hyperbolic Drucker-Prager material model.
%

% Newly Added by Yuan
In the error maps, the relative stress norm $\delta_{norm}$ is defined as 
\begin{equation}
  \delta_{norm} = \dfrac{\| \sigma_{ij} - \sigma_{ij}^{exact} \| }{ \| \sigma_{ij}^{exact} \| }
\end{equation}

where stress $ \sigma_{ij} $ is the calculated stress result, and 
stress $\sigma_{ij}^{exact}$ is the exact stress 
from either analytical solutions or Richardson extrapolation.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Error Map for von-Mises Perfectly Plastic Material Model} 

\label{subsec_example_vonMises_yf}

In  this  section, the von-Mises perfectly plastic material model is verified by
the prescribed solutions forcing (PSF) method.
%
% The perfectly plastic material assumes that 
% the elastic modulus becomes zero outside the yield surface. 
This type of material uses the radial return algorithm, \citep{local-78,Lubarda2002a}.
%
%which can be prescribed for error quantification.

% Both forward and backward Euler integration algorithms are verified 
% with von-Mises perfectly plastic material model. 
At the beginning of the verification test, 
the current stress state is set on the yield surface 
at the Lode angle $\theta=-30^{\circ}$. 
% The radius of the yield surface is plotted as 1 unit. 
Different elastic predictors are applied to this initial stress state.
%
The explicit, forward Euler constitutive  integration  algorithm  solve  the problem exactly for radial
stress predictors as solution consists of an exact radial return.
%
This is observed in zero error  Figure~\ref{fig_error_map1}, 
%
% when the normal to the yield surface of the starting stress 
% is on the right boundary at the Lode angle $\theta=-30^{\circ} $, 
% the numerical error is zero along the right boundary 
% at the Lode angle $\theta=-30^{\circ} $. 
However,  when  the  direction  of  the elastic predictors deviates far from the
normal   to  the  yield  surface,  radial  direction,  explicit,  forward  Euler
constitutive  integration  algorithm  creates  an  error.
%
That  error  increases  as  the  predictor  stress  increases,  as  observed  in
Figure~\ref{fig_error_map1}.
% gets greater 
%when the elastic predictor size increases.
\begin{figure}[!htbp]
  \centering
  \includegraphics[width=9cm]{../Figure-files/ch-error-map/error_relative_norm_1vm_perfect_plastic.pdf}
  \caption{
  \label{fig_error_map1}
  Relative stress norm of von-Mises perfectly plastic material model 
  with the forward Euler algorithm.
%   The radius of the yield surface is 1 unit. 
%   The stress starts on the yield surface at Lode angle $-30^{\circ}$. 
%   All possible stress predictors are tested and compared to the exact solution
%   on the yield surface with the corresponding Lode angles. 
%   When the predictors are perpendicular to the yield surface, the relative 
%   stress errors are zeroes. 
%   When the predictors deviate from the normal line, the relative stress
%   errors increase.
  }
\end{figure}






The   explicit, forward Euler constitutive  integration  algorithm together with
sub-incrementation  technique,  using  100 sub-incrementation, is able to
reduce  the numerical errors and return the stress states more accurately to the
yield surface, as shown in Figure~\ref{fig_error_map2}.
%
However, sub-incrementation does require much more time, so benefits of accuracy
increase have to be balanced with more time required for simulation.
%
% However, the algorithm still has some deficiencies when the elastic predictors
% deviates far from the normal to the yield surface.
% The two error maps in 
% Figure~\ref{fig_error_map1} and Figure~\ref{fig_error_map2} 
% are from the forward Euler algorithm. 
% the forward Euler algorithm do not have accuracy checking.
%
\begin{figure}[!htbp]
  \centering
  \includegraphics[width=9cm]{../Figure-files/ch-error-map/error_relative_norm_1vm_perfect_plastic_multiStep.pdf}
  \caption{
  \label{fig_error_map2}
  Relative stress norm of von-Mises perfectly plastic material model 
  with \textit{sub-increment technique},
  which reduces the relative stress error. 
  }
\end{figure}



For  implicit,  backward  Euler constitutive  integration algorithm, accuracy is
prescribed by the analyst as a yield function and/or residual stress tolerance.
%
%guarantees desired accuracy.
The  numerical  errors  for  the  of    implicit,  backward  Euler constitutive 
integration algorithm are shown in Figure~\ref{fig_error_map3}.
%
\begin{figure}[!htbp]
  \centering
  \includegraphics[width=9cm]{../Figure-files/ch-error-map/error_relative_norm_1vm_perfect_plastic_backward.pdf}
  \caption{
  \label{fig_error_map3}
  Relative stress norm of von-Mises perfectly plastic material model 
  with \textit{backward Euler} algorithm,
  which achieves accurate stress results. 
  }
\end{figure}



By comparing error maps shown in Figs.~\ref{fig_error_map1}, 
~\ref{fig_error_map2}, and ~\ref{fig_error_map3},
it is obvious that backward Euler algorithm achieves the highest  accuracy, by far.
%
However,   backward   Euler   constitutive  integration  algorithm  incurs  high
computational  costs, as higher order derivatives of plastic flow directions and
equilibrium iterations are needed.
%
Forward  Euler  algorithm  is  rather  simple  and  fast  algorithm  that can be
alternatively  used,  perhaps  with  the sub-incrementation option to solve 
problems.

% 
% presents a relatively low accuracy
% because forward Euler algorithm is designed to be a simple return algorithm 
% for speed with sacrificing accuracy. 





% =================================================================================
% =================================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

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

\clearpage
\subsection{Error  Map  for Hyperbolic Drucker-Prager Linear Kinematic Hardening
Material Model}
\label{subsec_example_DruckerPrager_linearHarding}

For  linear  hardening,  modified hyperbolic Drucker-Prager material model, exact
solution does not exist.
%
Linear kinematic hardening simplifies solution, however the nonlinear part of the modified
Drucker-Prager yield surface adds additional nonlinearities that introduce error
in constitutive integrations.
%
In  this  case  Richardson  extrapolation  is  used  to  verify the constitutive
integration algorithms.
%
Richardson extrapolation provides means to obtain accurate results that are then
used for error calculations.
%
Accurate  results  are developed using Richardson extrapolation and three groups
of tests, in order to determine accuracy order.
\begin{itemize}

  \item {Test 1.} The elastic predictor is applied in 1 increment. 

  \item {Test 2.} The elastic predictor is applied in 10 sub-incrementation, 
                  where the total magnitude of  10 sub-incrementation is 
                  the same as the 1 elastic predictor in Test 1.

  \item {Test 3.} The elastic predictor is applied in 100 sub-incrementation, 
                  where the total magnitude of 100 sub-incrementation is 
                  the same as the 1 elastic predictor in Test 1.
\end{itemize}

Results from these numerical experiments are then used to calculate the accuracy
order using Equation~(\ref{eq_richardson_estimate_final}).
% 
With  calculated  accuracy  order,  accurate solution can be developed, and then
used in calculating errors for different algorithms, as shown below.




Figure~(\ref{fig_nonlinear_ys_explicit})  shows  results  for  explicit, forward
Euler constitutive integration algorithm, applied to linear  hardening,  modified
hyperbolic Drucker-Prager material model.
% %


% %******************************************************************************************
% \begin{figure}[!htbp]
% \centering
% \captionsetup{justification=centering}
% 
% \noindent
% \includegraphics[width = 9.6cm]{../Figure-files/Richardson/explicitErrMap.pdf}
% \caption{
%   \label{fig_iso_compare_LK_E}
%   Linear kinematic hardening model, explicit algorithm.  
% }
% \end{figure}
% %******************************************************************************************


\begin{figure}[!htbp]
\centering
\captionsetup{justification=centering}
\includegraphics[width = 9cm]{../Figure-files/Richardson/explicitErrMap.pdf}
\caption{
  \label{fig_nonlinear_ys_explicit}
  Error map of explicit algorithms for 
  hyperbolic Drucker-Prager linear kinematic hardening material model.
%  %
%    The horizontal axis represents the $p-axis$ (mean pressure),
%    while the vertical axis represents the $q-axis$ (deviatoric stress).
%    The yellow line is the hyperbolic Drucker-Prager yield surface
%    and the red dot in the middle is the starting point of the stress.
%    The elastic predictors are specified to different locations across 
%    the $p-q$ stress space. 
%    To each specified location, three groups of tests are conducted.
}
\end{figure}

It       is      noted      that      error      maps,      as      show      in
Figure~(\ref{fig_nonlinear_ys_explicit}),  are  presented  in  deviatoric stress
plane,  where  horizontal  axes represents mean stress ($p=-(1/3) \sigma_{ii}$),
and   vertical   axes   represents   shear   stress   invariant  
($q=  \sqrt{(3/2) s_{ij}s_{ji}}$) and stress deviator $s_{ij}$ is defined as 
$s_{ij} = \sigma_{ij} - ({1}/{3}) \delta_{ij} \sigma_{kk}$, 
and $\delta_{ij}$ is the Kronecker delta.
%
It is noted that starting point is at the location of $p=8$kPa and $q=8$kPa, and
is marked with a red dot.
%
Error   map   for   explicit   constitutive   integration  algorithm,  shown  in
Figure~\ref{fig_nonlinear_ys_explicit}  reveals a potentially problematic stress
region on the tensile side of the yield surface.
%
However, it is noted that regular Drucker-Prager yield surface would not be even
able  to  produce  results  for  this  side  of  stress  space, as yield surface
derivatives are not even defined in that region.


Figure~\ref{fig_nonlinear_ys_implicit} show results for the  implicit, backward
Euler constitutive integration algorithm, applied to linear  hardening,  modified
hyperbolic Drucker-Prager material model.
%
\begin{figure}[!htbp]
\centering
\captionsetup{justification=centering}
\includegraphics[width = 9cm]{../Figure-files/Richardson/implicitErrMap.pdf}
\caption{
  \label{fig_nonlinear_ys_implicit}
  Error map of implicit algorithms for 
  hyperbolic Drucker-Prager linear kinematic hardening material model.
}
\end{figure}
%
It is noted that errors are much smaller than for explicit algorithm. 
%
It  is  also  noted  that  tensile  region  still tends to introduce more error,
however errors are still much smaller than with the explicit algorithm.



% 
% \begin{figure}[!htbp]
% \centering
% \captionsetup{justification=centering}
% 
% \noindent
% \includegraphics[width = 6.0cm]{../Figure-files/Richardson/explicitErrMap.pdf}
% % \hspace{1mm}
% \includegraphics[width = 6.0cm]{../Figure-files/Richardson/implicitErrMap.pdf}
% 
% 
% \caption{
%   \label{fig_DP_linearKinematicHardening}
%   (a) Left: linear kinematic hardening using explicit algorithm.  \\
%   (b) Right: linear kinematic hardening using implicit algorithm. \\
% }
% \end{figure}
% 


% =================================================================================
% =================================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\clearpage
\subsection{Error Map for Hyperbolic Drucker-Prager Armstrong-Frederick Nonlinear Kinematic Hardening Material Model}
\label{subsec_example_DruckerPrager_nonlinear_hardening}

More    advanced    material  models  used  in  elastoplastic modeling of cyclic
behavior feature nonlinear kinematic hardening.
%
It  is thus important to verify constitutive integration for the nonlinear kinematic hardening
material models.
% is critical to the
%reliability of the nonlinear elastic-plastic integration algorithm codes.
%
In this section,  Richardson Extrapolation technique is used to verify 
constitutive integration of nonlinear kinematic hardening material model.
% are verified with 
%Richardson Extrapolation technique. 
%
Material model used here is based on the Armstrong-Frederick
nonlinear kinematic hardening function \citep{Armstrong66}.
%
Using    similar    approach       for   estimating   accuracy   order   as   in
Section~\ref{subsec_example_DruckerPrager_linearHarding},   and  the  Richardson
Extrapolation algorithm, the error map for constitutive integration of elastoplastic 
material with Armstrong-Frederick nonlinear kinematic hardening law is 
shown in Figures.~\ref{fig_DP_nonlinearKinematicHardening_Explicit} and 
\ref{fig_DP_nonlinearKinematicHardening_Implicit}
%  %
%  \begin{figure}[!htbp]
%  \centering
%  \captionsetup{justification=centering}
%  \includegraphics[width = 9cm]{../Figure-files/Richardson/af_explicitErrMap.pdf}
%  \caption{  Error  map  for  explicit  constitutive integration algorithm using a
%  hyperbolic Drucker-Prager Armstrong-Frederick material model.}
%  \label{fig_armstrong_frederick_explicit}
%  \end{figure}
%  %
%  \begin{figure}[!htbp]
%  \centering
%  \captionsetup{justification=centering}
%  \includegraphics[width = 9cm]{../Figure-files/Richardson/af_implicitErrMap.pdf}
%  \caption{  Error  map  for  implicit  constitutive integration algorithm using a
%  hyperbolic Drucker-Prager Armstrong-Frederick material model.}
%  \label{fig_armstrong_frederick_implicit}
%  \end{figure}
%  


%
\begin{figure}[!htbp]
\begin{center}
\includegraphics[width = 9.0cm]{../Figure-files/Richardson/af_explicitErrMap.pdf}
\caption{
  \label{fig_DP_nonlinearKinematicHardening_Explicit}
  Error map of explicit algorithms for 
  hyperbolic Drucker-Prager nonlinear kinematic hardening material model.
}
\end{center}
\end{figure}



%
\begin{figure}[!htbp]
\begin{center}
\includegraphics[width = 9.0cm]{../Figure-files/Richardson/af_implicitErrMap.pdf}
\caption{
  \label{fig_DP_nonlinearKinematicHardening_Implicit}
  Error map of implicit algorithms for 
  hyperbolic Drucker-Prager nonlinear kinematic hardening material model.
}
\end{center}
\end{figure}
%
%
%%{\large \bf TO ADD/EDIT more description} 
%
A number of observations can be made.
%


\begin{itemize}

\item 
The  errors  using  implicit  constitutive integration is much smaller than that
using  explicit constitutive integration. It is noted that the error in implicit
algorithm is at least two orders of magnitude smaller than the error in explicit
algorithm, as seen in Figures~\ref{fig_DP_nonlinearKinematicHardening_Explicit} and 
\ref{fig_DP_nonlinearKinematicHardening_Implicit}.
%

\item 
Comparison of errors for explicit constitutive algorithm for linear and nonlinear kinematic
hardening     models,     Figures~\ref{fig_nonlinear_ys_explicit} and 
\ref{fig_DP_nonlinearKinematicHardening_Explicit} ,
reveals  that  nonlinearity  in  hardening  did
introduce additional error, as expected.
%

\item 
However, comparison of errors for implicit constitutive algorithm for linear and
nonlinear                kinematic               hardening               models,
Figures~\ref{fig_nonlinear_ys_implicit}  and \ref{fig_DP_nonlinearKinematicHardening_Implicit}
reveals  very  similar  error  distribution.    This is also expected as for the
implicit  algorithm  convergence  tolerance  is  prescribed and it holds for any
material model used.


% \item Implicit, backward Euler algorithm has much less error than 
%       explicit, forward algorithm in both the linear and nonlinear 
%       kinematic hardening material models.
% 

\item Constitutive integration error  for both implicit and explicit algorithms 
increases  as  the  mean pressure is decreases, that is, as the stress states is
closer to the nonlinear part of the yield surface.



\end{itemize}




% %******************************************************************************************


% %******************************************************************************************
% \begin{figure}[!htbp]
% \centering
% \captionsetup{justification=centering}
% \includegraphics[width = 9.8cm]{../Figure-files/Richardson/af_explicitErrMap.pdf}
% \caption{
%   \label{fig_iso_compare_NL_E}
%   Nonlinear kinematic hardening model, explicit algorithm.  
% }
% \end{figure}
% %******************************************************************************************



% %******************************************************************************************
% \begin{figure}[!htbp]
% \centering
% \captionsetup{justification=centering}
% \includegraphics[width = 9.8cm]{../Figure-files/Richardson/af_implicitErrMap.pdf}
% \caption{
%   \label{fig_iso_compare_NL_I}
%     Nonlinear kinematic hardening model, implicit algorithm. 
% }
% \end{figure}
% %******************************************************************************************


% % changed all four figs to have same error iso maps
% %******************************************************************************************
% %*********************** TEMP For Comparison                  *************************
% %******************************************************************************************
% \begin{figure}[!htbp]
% \centering
% \captionsetup{justification=centering}
% 
% \noindent
% \includegraphics[width = 7.6cm]{../Figure-files/Richardson/explicitErrMap.pdf}
% \hspace{5mm}
% \includegraphics[width = 7.6cm]{../Figure-files/Richardson/implicitErrMap.pdf}
% 
% \noindent
% \includegraphics[width = 7.8cm]{../Figure-files/Richardson/af_explicitErrMap.pdf}
% \hspace{5mm}
% \includegraphics[width = 7.8cm]{../Figure-files/Richardson/af_implicitErrMap.pdf}
% 
% \caption{
%   \label{fig_temp_all_4_fig_iso_compare}
%   Temp figures for comparison\\
%   (a) Left top: linear kinematic hardening explicit.  
%   (b) Right top: linear kinematic hardening implicit. 
%   (c) Left bottom: nonlinear kinematic hardening explicit. 
%   (d) Right bottom: nonlinear kinematic hardening implicit. 
% }
% \end{figure}
% %******************************************************************************************
% %*********************** TEMP For Comparison      END            *************************
% %******************************************************************************************




% =================================================================================
% =================================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\subsection{Numerical  Accuracy      Error   and  Uncertainty  Estimation  for Drucker-Prager
\clearpage
\subsection{Numerical  Accuracy  Estimation    for    Associated  Drucker-Prager
Linear Kinematic Hardening Material Model}
\label{subsec_example_associative}


% 
% The previous three Sections~\ref{subsec_example_vonMises_yf},
% ~\ref{subsec_example_DruckerPrager_linearHarding}, and 
% \ref{subsec_example_DruckerPrager_nonlinear_hardening} show the
% error map of the von-Mises and hyperbolic Drucker-Prager yield surfaces.
% %
% However, the numerical uncertainty is still unknown.
% %
% %
% In the following sections, the numerical uncertainty of the elastoplastic algorithms 
% is quantitatively studied and evaluated. 
%


Presented  in  this  section are estimates of numerical accuracy of constitutive
integration for a Drucker-Prager material model with linear kinematic hardening.
%
Three groups of tests are conducted using Drucker-Prager associative material model.
%
For  each  of  numerical  experiments,  the  total strain increment is the same,
however,  the number of increments to reach this strain are different. This means
that the size of strain sub-increments for each test differs.
%
Table~\ref{table_richardson_DP_vary_pressure}  presents  results  of  Richardson
extrapolation   for  three different starting stress point on the Drucker-Prager
yield  surface,  and  corresponds  to  three confining pressures 
($p=10^4{\rm Pa}; 10^5{\rm Pa}; 10^6{\rm Pa} $).
% represents the different groups of tests 
%on the pressure-dependent Drucker-Prager material model. 
%

\begin{table}[!htbp]
\centering
\captionsetup{justification=centering}
\caption{
  \label{table_richardson_DP_vary_pressure}
  Richardson  extrapolation  results  for  shear  stress  $q$ for Drucker-Prager
  linear kinematic hardening material with different confinements.
}
  \begin{tabu}{l|l|l|l}
%    \tabucline[1.2pt]{-}
    Confinement  &    $10^4$Pa  &   $10^5$Pa  &  $10^6$Pa    \\ 
    \hline
    \hline
%    \tabucline[1.2pt]{-}

    Accuracy order $\beta$ & 1.0003 &  1.0017   &  0.9954 \\ \hline 

    GCI  &  0.01\% &  0.01\% &  0.00\% \\ \hline 

    Resulting q  &  32207.705Pa &  41866.802Pa  &  127602.165Pa \\ \hline 

    Richardson exact result for q & 32205.934Pa & 41864.904Pa & 127600.876Pa \\ \hline \hline

    Range of error in q & $\pm$2.213Pa &  $\pm$2.373Pa &  $\pm$1.611Pa \\ 

%    \tabucline[1.2pt]{-}
  \end{tabu}
\end{table}


Figure~\ref{fig_convergen_DP_vary_pressure},   shows   results   for  asymptotic
convergence  of  Drucker-Prager  associative linear kinematic hardening material
model with different starting points (different confinment levels).

\begin{figure}[!htbp]
\centering
\includegraphics[width = 9cm]{../Figure-files/DPAFdiffpressureConverge.pdf}
\captionsetup{justification=centering}
\caption{
  \label{fig_convergen_DP_vary_pressure}
  Asymptotic convergence of Drucker-Prager linear kinematic hardening material with 
  different confinements using implicit, backward Euler algorithm.
}
\end{figure}

%


Based on Figure~\ref{fig_convergen_DP_vary_pressure}, 
A number of observations can be made:
\begin{itemize}
  \item When the strain increment size is greater than $1e-3$,
    the relative error oscillates 
    since the discretization size is too large 
    to represent the actual continua space 
    in the discretized domain.
    %
  \item When the strain increment size is between $1e-13$ and $1e-3$, 
    the relative error decreases asymptotically. 
    This region is called the asymptotic regime of convergence.
    %
  \item When the strain increment size is smaller than $1e-13$,
    the relative error cannot decrease anymore since
    round-off error contributes to the solution error.
    %
  \item When the strain increment is $1e-16$ (less than the machine epsilon $2.22e-16$),
    the stress result oscillates significantly 
    since the round-off error dominates the solution error.
\end{itemize}

Besides, the relative error decreases when the confining pressure increases.
%
This is due to the fact that high pressure states of stress are far from the 
nonlinearities  in the modified hyperbolic function. 
%
This nonlinearity leads to 
increase  in  integration  error,  while  moving  away  from  such  nonlinearity
decreases of the relative error.















% =================================================================================
% =================================================================================
% =================================================================================
% =================================================================================
% =================================================================================
% =================================================================================
% =================================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Numerical  Accuracy  Estimation  for  Non-Associated  Drucker-Prager
Nonlinear Kinematic Hardening Material Model}
\label{subsec_example_non_associative}

%OVDE

In  this section, a non-associated Drucker-Prager material model with nonlinear,
Armstrong-Frederick kinematic hardening is investigated.
%
Both explicit, forward Euler and implicit, backward Euler algorithms are tested with 
various strain incremental sizes.
%
Stress solution, in terms of deviatoric stress 
($q= \sqrt{3/2 s_{ij}s_{ji}}, s_{ij}= \sigma_{ij}- \sigma_m \delta_{ij}$) 
is  used  to illustrate change of accuracy of algorithms for various incremental
step sizes.
%
Figure~\ref{fig_non_associate_plastic_flow_small_increments} shows 
the comparison of deviatoric stress solution for a variation in 
incremental strain size for implicit and explicit algorithms.
%
\begin{figure}[!htbp]
\centering
\captionsetup{justification=centering}
\includegraphics[width = 9cm]{../Figure-files/GCI-Uncertainty/SmallNonAssociate.pdf}
\caption{
\label{fig_non_associate_plastic_flow_small_increments}
Deviatoric  Stress  $q$  results  for a non-associative Drucker-Prager nonlinear
kinemastic  hardening  material  using  implicit,  backward  Euler and explicit,
forward Euler algorithms.}
\end{figure}

%
It is noted that implicit algorithm behaves as expected, that is stress solution
is very accurate even for large strain increments, while explicit algorithm does
produce significant error for larger step sizes.
%
For  explicit  algorithm,  error  is  significantly  reduced  when the step size
becomes small.
%
This argument is used for advocating sub-stepping approach in order to improve 
explicit algorithm accuracy.
%






In  addition  to  stress  result, yield surface values of the stress results are
also used to illustrate numerical accuracy of the algorithms. 
%
Ideally,   constitutive  integration algorithm should return the stress state to
the yield surface, yielding yield surface values to (almost) zero.
%
If      stress    is      not  returned  to  the  yield  surface, equilibrium of
internal stress and external forces might not be achievable.
%
In  addition,    calibration    of    the elastic-plastic material model is done
assuming that stress state is on the yield surface during plastification.
%
If    stress   state is not on the yield surface, then the calibrated parameters
for material model will not be able to replicate real response.
%



Figure~\ref{fig_non_associate_plastic_flow}, shows yield surface values for both
implicit and explicit algorithms for different incremental step sizes.
%
\begin{figure}[!htbp]
\centering
\captionsetup{justification=centering}
\includegraphics[width = 9cm]{../Figure-files/GCI-Uncertainty/NonAssociate.pdf}
\caption{
\label{fig_non_associate_plastic_flow}
Yield surface values of returned stress of non-associative plastic material
using implicit, backward Euler and explicit, forward Euler algorithms.
}
\end{figure}



It  is  noted that, as expected, for explicit algorithm, stress solutions drifts
away  from  the yield surface for larger incremental step sizes, hence producing
erroneous stress solution with calibrated material parameters.
%
On  the other hand, implicit algorithm performs much better, since user controls
stress solution drift from the yield surface by prescribing drift tolerance.

While  explicit  integration  algorithm is simpler to develop and implement than
the  explicit  integration algorithm, it does produce larger error for realistic
incremental strain steps. 
%
Use of sub-incrementation is possible at a cost of increasing computational times.
%
On  the  other  hand,  implicit  integration  algorithm  is more complicated and
requires  iterations,  thus  increasing computational time, however, accuracy is
superior.



%
%  The  reasons  for  the  opposite  trends  in  explicit, forward Euler and 
%  implicit, backward Euler algorithms are discussed as follows:
%  
%  \begin{itemize}
%  
%    \item 
%    In explicit forward Euler algorithm, the stress state at the cross-point
%    between  the  yield surface and the stress prediction is used as the direction
%    of the return-mapping algorithm.
%  %
%    When the elastic predictor is relatively large, the destination of the plastic
%    corrector  lies  outside  the  yield surface, which leads to the larger
%    yield surface values of the final stress.
%  %
%  
%    \item  
%    In implicit backward Euler algorithm, the stress state at the predictor
%    stress  state  is  used to calculate the magnitude of the plastic corrector in
%    the return-mapping algorithm, which is more accurate.
%  %
%    Also, there is a tolerance check to guarantee that the yield surface values
%    satisfy the tolerance requirements. 
%  %
%  
%  
%  \end{itemize}
%  



% =================================================================================
% =================================================================================
% =================================================================================
% =================================================================================
% =================================================================================
% =================================================================================
% =================================================================================
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \section{Pitfalls and Issues to Consider in Verification Process}
% \label{sec_pitfall_verification}
% 
% 
% Finally, we want to mention some points which have to be considered 
% in any code and solution verification studies:
% %
% \begin{enumerate}
% \item 
% Although there is no predefined fixed procedure for verification.
% %
% In general, modular thinking to error/limitation sources is crucial 
% to understand, narrow down, segregation and remove the errors. 
% %
% In-depth knowledge of software quality assurance (SQA), 
% numerical analysis, and physics of the phenomenon 
% could be helpful in those regards. 
% 
% \item
% Given that the concept that mesh convergence study is 
% deviated form Lax equivalence theorem, when it is violated in 
% shock formation situation in nonlinear PDEs, the theorem is not valid. 
% 
% \item  
% Difficulties induced by the scale of the problem: if the test problem 
% for checking a system of PDEs solver is chosen in such a way that one 
% of the including phenomena is becoming significantly larger than the 
% others, the imperfection in the low effect phenomenon may be concealed 
% by the dominant effect.
% 
% \item
% In mesh convergence study, when norms/integrated values are used 
% instead of the original state variables, we have to consider the smoothing effect 
% by using them. 
% %
% For example study of $L_2$ norm (energy norm) always yields a better result 
% as compared to the state variable or $L_\infty$ norm of the error.
% 
% \item
% The problem of accuracy order of other pre/post-processing algorithm(s): 
% any extra pre-processed or post-processed value which is used 
% for verification must be computed via a high-order method to avoid 
% external pollution of the results.
% 
% \item
% The problem of PDE identicalness: mesh convergence test has to be done 
% on an exactly identical differential equation which means the exactly 
% same geometry, initial condition, and boundary conditions, 
% and solution/solver parameters. 
% %
% This makes the problem of mesh convergence study difficult 
% when there are embedded switches in the algorithm 
% (like a material model switch from plasticity to elasticity).
% 
% \item
% Richardson Extrapolation (RE), like all other extrapolations in mathematics 
% might not end up to better results in the presence of high nonlinearities 
% in the system. 
% %
% It is strongly recommended that RE is being used carefully and 
% at least on three consecutive mesh sizes. 
% 
% \item 
% Sometimes a system of PDEs are solved by iterations between different 
% equations (for example Navier-Stokes solvers) instead of a
% fully implicit solution in one shot, it those cases the number of iterations 
% might induce extra numerical error, other than discretization, 
% round-off and algebraic iteration error. 
% 
% \end{enumerate}
% 
% In summary, verification techniques, including PSF, Richardson extrapolation, and 
% grid convergence index are a necessary but not sufficient condition 
% of the soundness of 
% a numerical simulation. 
% 
% %
% Verification studies, as well as design, 
% implementation, and interpretation of verification tests, 
% are complex tasks which require a deep understanding of mathematics, physics 
% and computer science. 
% 
% %
% Sometimes finding a problem could only achieved 
% with a state-of-art innovative numerical technique.    
% 
% % =================================================================================
% % =================================================================================
% % =================================================================================
% % =================================================================================
% =================================================================================
% =================================================================================
% =================================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusion and Final Remarks}


This paper deals with verification of constitutive integration algorithms as 
applied to elastic-plastic materials models used in civil engineering.
%
Both   explicit,  forward  Euler  and  implicit,  backward  Euler  elastoplastic
algorithms are verified on the constitutive level.
%
Three  methods  are  used  in the verification process, including the prescribed
solution  forcing,  Richardson extrapolation, and the extension grid convergence
index (GCI).

In   the   particular   case   of   elastic-plastic  constitutive  integrations,
verification was performed using:
%
\begin{enumerate}
 \item 
 Prescribed  solution  forcing  in  terms of analytical solutions, to verify the
 perfectly plastic material with a hyperbolic Drucker-Prager yield surface.
 
 \item
 Richardson extrapolation, to verify the 
 nonlinear kinematic hardening material models. 

 \item
 GCI method, to estimate the numerical accuracy in 
 the elastoplastic algorithm, especially for Drucker-Prager material models with
 non-associate plastic flow.

\end{enumerate}


It was shown that constitutive algorithms that were subject to verification
perform well and as expected, thus increasing confidence in elastic plastic modeling.

It  is  also  important to note that use of numerical modeling for prediction of
behavior of solids and structures has to be done using verified numerical tools.
%
This   is   rarely   demonstrated   in   project  documentation  and  scientific
publications, which lowers confidence in presented results.




\begin{acknowledgements} ~

Funding from the United States Department of Energy is acknowledged.

\end{acknowledgements}

% BibTeX users please use one of
%\bibliographystyle{spbasic}      % basic style, author-year citations
%\bibliographystyle{spmpsci}      % mathematics and physical sciences
%\bibliographystyle{spphys}       % APS-like style for physics
\bibliographystyle{abbrvnat}
%\bibliographystyle{unsrt}
\bibliography{refmech,refcomp}   % name your BibTeX data base


\end{document}
% end of file template.tex






















































% % For tables use
% \begin{table}
% % table caption is above the table
% \caption{Please write your table caption here}
% \label{tab:1}       % Give a unique label
% % For LaTeX tables use
% \begin{tabular}{lll}
% \hline\noalign{\smallskip}
% first & second & third  \\
% \noalign{\smallskip}\hline\noalign{\smallskip}
% number & number & number \\
% number & number & number \\
% \noalign{\smallskip}\hline
% \end{tabular}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \section{Estimation of Numerically Induced Uncertainty}

% Theory of Roache GCI by Kaveh
% Example of GCI by Yuan + a picture of mesh size versus uncertainty bands
% \end{table}

% =================================================================================
% =================================================================================
% =================================================================================
% =================================================================================
% =================================================================================
% =================================================================================
% =================================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \section{Code Verification -- Linear kinematic Hardening Materials with Prescribed Forcing Solutions}

% \label{sec_linear_hardening_MMS}
% Given that, when the elastic-plastic materials have hardening properties, there is 
% no exact analytical solutions for conducting the code verification. 
% In this section, prescribed solutions forcing (PSF) \citep{Roache2002}
% \citep{Salari2000} are applied to verify the linear kinematic hardening materials.
% Prescribed solutions forcing are applicable to the uniaxial linear kinematic hardening materials, 
% where the stiffness is known before the elastic-plastic algorithms.
% However, it is not suitable to the nonlinear kinematic hardening materials 
% because the stiffness tensor is calculated after the integration algorithms. \textcolor{red}{I do not understand this part!} 
% % ==============================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \subsection{Advantages of Prescribed Solutions Forcing}

% Prescribed forcing solutions \citep{Roache2002} \citep{Leng2013} \citep{Petri2015} 
% has been widely used in computational fluid dynamics (CFD). 
% The idea of prescribed forcing solution is to explicitly assume a solution for the target 
% equations. Then, the assumed solution is substituted 
% into the equation to obtain the initial and boundary conditions.
% The procedures above are outside the program we want to verify. 
% Then, the initial and boundary conditions are transformed as the input of
% the solver. Finally, the output of the solver is compared to the 
% assumed solutions in the beginning, to find the evolution of the error. 

% In the context of elastic-plastic algorithms, the initial and final stress
% are assumed in the beginning. Then, the incremental stress is used to 
% generate the incremental strain through the compliance tensor. 
% The procedures above are outside the program we want to verify. 
% Then, the incremental strains are input to the elastic-plastic algorithms.
% The output contains the stress, which is compared to the assumed stress
% in the beginning. These procedures are illustrated in Figure~\ref{fig_mms_linear}.


% \begin{figure}[!htbp]
%   \centering
%   \includegraphics[width=5cm]{../Figure-files/mms_linear.pdf}
%   \caption{Prescribed Forcing Solutions for Linear kinematic Hardening Materials}
%   \label{fig_mms_linear}
% \end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \subsection{Limitation of Prescribed Forcing Solutions}


% Prescribed forcing solutions are not applicable to the nonlinear kinematic hardening materials. 
% In a typical elastic-plastic algorithms, 
% the stiffness tensor is calculated after the elastic-plastic algorithms. 
% \begin{figure}[!htbp]
%   \centering
%   \includegraphics[width=8cm]{../Figure-files/algorithm_input_output.pdf}
%   \caption{Input and Output of Elastoplastic Algorithms}
%   \label{fig_elastic-plastic_algo_input_output}
% \end{figure}

% Since the compliance tensor is the inverse of the stiffness tensor, 
% compliance tensor is also calculated after the elastic-plastic algorithms. 
% However, according to the data flow of prescribed forcing solutions in Figure~\ref{fig_mms_linear}, 
% the compliance tensor is required as the input to the elastic-plastic algorithms. 
% This becomes a chicken-and-egg problem: 
% If we want to obtain the input strain increment, we have to obtain the 
% compliance tensor first, which is the inverse of the stiffness tensor. 
% However, according to Figure~\ref{fig_elastic-plastic_algo_input_output}, 
% if we want to know the stiffness tensor, we have to obtain the strain increment 
% first. 

% In the special case of uniaxial linear kinematic hardening materials, 
% the stiffness and compliance tensors are known before the elastic-plastic algorithms. 
% A deliberate selection of the linear kinematic hardening parameters guarantees that 
% compliance tensor is known before the numerical experiments 
% such that prescribed solution forcing  are applied to the uniaxial linear kinematic hardening materials.

% \textcolor{red}{we have to rewrite this part in Davis}




% The two-step predictor-corrector algorithms are made of an elastic predictor 
% and a plastic corrector. 
% Based on the returning directions, Euler algorithms are divided into forward Euler 
% and backward Euler algorithms \citep{local-64}. 
% On the side of the Euler algorithms, 
% the generalized trapezoidal return mapping applies the interpolated direction between 
% forward Euler and backward Euler.
% Likewise, the generalized midpoint return mapping applies the interpolated stress between 
% the previous stress states and the initial return stress in the current step.