%%%%%%%%%%%%%%%%%%%%%%%% 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}
%!PSAdobe3.0 EPSF3.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{fixcm}
%
\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 elasticplastic
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,
elasticplastic integration algorithms, as used in computational nonlinear solid
mechanics. Both explicit and implicit integration algorithms for elasticplastic
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 elasticplastic integration algorithms is applied to a number of
material models: from simple vonMises perfectly plastic to sophisticated
hyperbolic DruckerPrager with nonlinear ArmstrongFrederick rotational
kinematic hardening. Besides, algorithmic uncertainty is estimated with both
associative and nonassociative 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
elasticplastic constitutive equations. Validation of elasticplastic
models is not within current scope of this paper and is left for future publications.
A general elasticplastic 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
elasticplastic 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 (MohrCoulomb or
DruckerPrager) \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 MohrCoulomb.
\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 nonassociated, 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{local91}), second order tensors (kinematic hardening \citep{Armstrong66})
and fourth order tensor (distortional hardening \citep{local99}).
\end{enumerate}
The four components, using different elastic laws, yield functions, plastic flow
directions and hardening laws, can combined to develop a number of
elasticplastic material models. In this paper two models are used for
verification, namely, von Mises and hyperbolic DruckerPrager elasticplastic
material model, both with linear and nonlinear hardening laws. Hyperbolic
DruckerPrager used here is developed based on the idea from \citet{Abbo1995}.
General numerical integration algorithms for elasticplastic constitutes equations are
stateupdating procedure.
Commonly used are explicit (forward Euler) and implicit (backward Euler) elasticplastic integrations
algorithms.
%
%%
%Since the elasticplastic 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_elasticplastic} shortly introduces the elasticplastic algorithms.
Section~\ref{sec_Verification_theory} presents the theoretical basis of the applied
verification techniques.
Section~\ref{sec_elasticplastic_material_model} briefly reviewed the verified
elasticplastic material model.
Section~\ref{sec_verification_example} shows the verification examples of different
elasticplastic 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 wellestablished finite element method
% computer package \citep{selby1997}.
% In another catastrophe in academia, a software bug destroyed the shining
% career of an outstanding and welldecorated 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 onetime 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
NuclearTestBan 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 elasticplastic constitutive integration algorithms.
%
Additionally, implemented and described is Roache Grid Convergence Index \citep{Roache1997} test
to uncover numerical uncertainty in the elasticplastic simulations.
% ====================================================================
% ====================================================================
% ====================================================================
% ====================================================================
% ====================================================================
% ====================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Elastoplastic Constitutive Integration Algorithm}
\label{sec_intro_elasticplastic}
Elastoplastic integration algorithm is used to solve for the stress increment
from a strain increment.
The elasticplastic integration algorithm is a twostep algorithm including the elastic predictor
and plastic corrector \citep{local64}.
The elastic predictor step is similar for different algorithms.
Depending on the direction of the plastic corrector, the elasticplastic
algorithms are generally divided into explicit, forward Euler and implicit,
backward Euler algorithms.
Both explicit and implicit algorithms are special cases of a general midpoint
algorithm \citep{local18,local64,Crisfield1997}. These two constitutive
algorithms represent the most common constitutive integration algorithms.
Explicit algorithm uses the starting elasticplastic 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 elastoplastic 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 elasticplastic stiffness tensor
($E_{ijkl}^{elpl}$) 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) secondorder 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 elasticplastic 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 nonpositive 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 midpoint 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 elasticplastic state, that is on
the yield surface, to calculate the increment in stress ($d \sigma_{ij}$) and
continuum, tangent elasticplastic stiffness tensor ${}^{cont} E_{ijkl}^{ep}$
\citep{local64,Crisfield1997,Chen94}.
%
The plasticmultiplier
$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{local19,local4}.
%
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 plasticmultiplier
$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 elasticplastic 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 roundoff 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]{../Figurefiles/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 elasticplastic 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 bigO 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 elasticplastic 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 abovementioned 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, roundoff error and
iterative error.
%
First, truncation error represents the error due to solving discretized
equations instead of continuum counterpart.
%
Second, roundoff 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 roundoff 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_elasticplastic}
Before the implementation of Richardson Extrapolation, the accuracy order
($\beta$) of the elasticplastic algorithm needs to be determined.
Both the analytical determination and the numerical estimation of accuracy order for
elasticplastic 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{local18}.
In their study, \citet{local18} proved that all elastic plastic
algorithms within generalized trapezoidal and generalized midpoint return
mapping are at least firstorder accurate.
In the case of CrankNicolson 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 elasticplastic algorithm meets the secondorder accuracy.
\noindent
{\bf Remark:} {Conclusions on algorithm accuracy order
made by \citet{local18} 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_elasticplastic_material_model}
Elastoplastic material models used for numerical experiments are based on von
Mises and hyperbolic DruckerPrager yield surfaces.
%
Material models of vonMises type are used to simulate the pressureindependent
material behavior.
%
On the other hand, hyperbolic DruckerPrager family of models are used to simulate
pressuredependent 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 vonMises 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 nonassociated 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{local99} where yield surface changes shape,
is also possible , but it is rarely used.
%
%
Plastic flow directions are usually developed either based on associated or
nonassociated, or, for more sophisticated model, directly as functions of
stress and internal variable space \citep{Dafalias82,local63,Manzari97}.
Presented below are common examples of yield surface functions for vonMises and
DruckerPrager models, and examples of nonassociated plastic flow as well as a
nonlinear kinematic hardening rule.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{The vonMises Yield Surface}
The von Mises yield surface \citep{local63}, 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 DruckerPrager Yield Surface}
The DruckerPrager yield criterion is a pressuredependent material model, where
yielding of materials is controlled by the deviatoric stress, and also by the
confining pressure.
%
The original DruckerPrager yield surface \citep{local63}
has a singularity point at the location of zero deviatoric stress. Hyperbolic
yield surface \citep{Abbo1995} modifies
the DruckerPrager 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{xd}{a}\right)^2  \left(\frac{y}{b}\right)^2 = 1
% \end{equation}
% %
% Extend the yield surface to include isotropic and kinematic hardening.
The hyperbolic DruckerPrager 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 cutoff
on $paxis$ and the original apex points, as shown in
Figure~\ref{fig_hyperbolic_function_diagram1}.
%
%
%
%Figure~\ref{fig_hyperbolic_function_diagram1} shows trace of DruckerPrager yield surface in stress space.
%
\begin{figure}[!htbp]
\centering
\includegraphics[width = 7cm]{../Figurefiles/hyperbolic_dp.pdf}
\caption{The hyperbolic DruckerPrager yield function.}
\label{fig_hyperbolic_function_diagram1}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Nonassociative Plastic Flow}
In pressuredependent 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 nonassociative 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 stressratio.
The material constants $k_d$ governs the critical stressratio,
which controls the direction of the plastic flow.
The dilative or contractive behavior is controlled by the stressratio
(${\sqrt{s_{mn}s_{mn}}}/{p}$)
\citep{local91}.
When dilatancy parameter $D<0$ plastic compaction takes place, whereas when
dilatancy parameter $D>0$ plastic
dilatancy takes place.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{ArmstrongFrederick 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 backstress,
$h_a$ is the initial rate of change of backstress,
$h_a/c_r$ controls the limit, assymptote of the backstress,
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 vonMises 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 nonassociated
hyperbolic DruckerPrager 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 vonMises Perfectly Plastic Material Model}
\label{subsec_example_vonMises_yf}
In this section, the vonMises 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{local78,Lubarda2002a}.
%
%which can be prescribed for error quantification.
% Both forward and backward Euler integration algorithms are verified
% with vonMises 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]{../Figurefiles/cherrormap/error_relative_norm_1vm_perfect_plastic.pdf}
\caption{
\label{fig_error_map1}
Relative stress norm of vonMises 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
subincrementation technique, using 100 subincrementation, 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, subincrementation 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]{../Figurefiles/cherrormap/error_relative_norm_1vm_perfect_plastic_multiStep.pdf}
\caption{
\label{fig_error_map2}
Relative stress norm of vonMises perfectly plastic material model
with \textit{subincrement 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]{../Figurefiles/cherrormap/error_relative_norm_1vm_perfect_plastic_backward.pdf}
\caption{
\label{fig_error_map3}
Relative stress norm of vonMises 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 subincrementation 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 DruckerPrager Linear Kinematic Hardening
Material Model}
\label{subsec_example_DruckerPrager_linearHarding}
For linear hardening, modified hyperbolic DruckerPrager material model, exact
solution does not exist.
%
Linear kinematic hardening simplifies solution, however the nonlinear part of the modified
DruckerPrager 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 subincrementation,
where the total magnitude of 10 subincrementation is
the same as the 1 elastic predictor in Test 1.
\item {Test 3.} The elastic predictor is applied in 100 subincrementation,
where the total magnitude of 100 subincrementation 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 DruckerPrager material model.
% %
% %******************************************************************************************
% \begin{figure}[!htbp]
% \centering
% \captionsetup{justification=centering}
%
% \noindent
% \includegraphics[width = 9.6cm]{../Figurefiles/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]{../Figurefiles/Richardson/explicitErrMap.pdf}
\caption{
\label{fig_nonlinear_ys_explicit}
Error map of explicit algorithms for
hyperbolic DruckerPrager linear kinematic hardening material model.
% %
% The horizontal axis represents the $paxis$ (mean pressure),
% while the vertical axis represents the $qaxis$ (deviatoric stress).
% The yellow line is the hyperbolic DruckerPrager 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 $pq$ 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 DruckerPrager 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 DruckerPrager material model.
%
\begin{figure}[!htbp]
\centering
\captionsetup{justification=centering}
\includegraphics[width = 9cm]{../Figurefiles/Richardson/implicitErrMap.pdf}
\caption{
\label{fig_nonlinear_ys_implicit}
Error map of implicit algorithms for
hyperbolic DruckerPrager 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]{../Figurefiles/Richardson/explicitErrMap.pdf}
% % \hspace{1mm}
% \includegraphics[width = 6.0cm]{../Figurefiles/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 DruckerPrager ArmstrongFrederick 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 elasticplastic 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 ArmstrongFrederick
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 ArmstrongFrederick 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]{../Figurefiles/Richardson/af_explicitErrMap.pdf}
% \caption{ Error map for explicit constitutive integration algorithm using a
% hyperbolic DruckerPrager ArmstrongFrederick material model.}
% \label{fig_armstrong_frederick_explicit}
% \end{figure}
% %
% \begin{figure}[!htbp]
% \centering
% \captionsetup{justification=centering}
% \includegraphics[width = 9cm]{../Figurefiles/Richardson/af_implicitErrMap.pdf}
% \caption{ Error map for implicit constitutive integration algorithm using a
% hyperbolic DruckerPrager ArmstrongFrederick material model.}
% \label{fig_armstrong_frederick_implicit}
% \end{figure}
%
%
\begin{figure}[!htbp]
\begin{center}
\includegraphics[width = 9.0cm]{../Figurefiles/Richardson/af_explicitErrMap.pdf}
\caption{
\label{fig_DP_nonlinearKinematicHardening_Explicit}
Error map of explicit algorithms for
hyperbolic DruckerPrager nonlinear kinematic hardening material model.
}
\end{center}
\end{figure}
%
\begin{figure}[!htbp]
\begin{center}
\includegraphics[width = 9.0cm]{../Figurefiles/Richardson/af_implicitErrMap.pdf}
\caption{
\label{fig_DP_nonlinearKinematicHardening_Implicit}
Error map of implicit algorithms for
hyperbolic DruckerPrager 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]{../Figurefiles/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]{../Figurefiles/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]{../Figurefiles/Richardson/explicitErrMap.pdf}
% \hspace{5mm}
% \includegraphics[width = 7.6cm]{../Figurefiles/Richardson/implicitErrMap.pdf}
%
% \noindent
% \includegraphics[width = 7.8cm]{../Figurefiles/Richardson/af_explicitErrMap.pdf}
% \hspace{5mm}
% \includegraphics[width = 7.8cm]{../Figurefiles/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 DruckerPrager
\clearpage
\subsection{Numerical Accuracy Estimation for Associated DruckerPrager
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 vonMises and hyperbolic DruckerPrager 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 DruckerPrager material model with linear kinematic hardening.
%
Three groups of tests are conducted using DruckerPrager 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 subincrements for each test differs.
%
Table~\ref{table_richardson_DP_vary_pressure} presents results of Richardson
extrapolation for three different starting stress point on the DruckerPrager
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 pressuredependent DruckerPrager material model.
%
\begin{table}[!htbp]
\centering
\captionsetup{justification=centering}
\caption{
\label{table_richardson_DP_vary_pressure}
Richardson extrapolation results for shear stress $q$ for DruckerPrager
linear kinematic hardening material with different confinements.
}
\begin{tabu}{llll}
% \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 DruckerPrager associative linear kinematic hardening material
model with different starting points (different confinment levels).
\begin{figure}[!htbp]
\centering
\includegraphics[width = 9cm]{../Figurefiles/DPAFdiffpressureConverge.pdf}
\captionsetup{justification=centering}
\caption{
\label{fig_convergen_DP_vary_pressure}
Asymptotic convergence of DruckerPrager 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 $1e3$,
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 $1e13$ and $1e3$,
the relative error decreases asymptotically.
This region is called the asymptotic regime of convergence.
%
\item When the strain increment size is smaller than $1e13$,
the relative error cannot decrease anymore since
roundoff error contributes to the solution error.
%
\item When the strain increment is $1e16$ (less than the machine epsilon $2.22e16$),
the stress result oscillates significantly
since the roundoff 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 NonAssociated DruckerPrager
Nonlinear Kinematic Hardening Material Model}
\label{subsec_example_non_associative}
%OVDE
In this section, a nonassociated DruckerPrager material model with nonlinear,
ArmstrongFrederick 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]{../Figurefiles/GCIUncertainty/SmallNonAssociate.pdf}
\caption{
\label{fig_non_associate_plastic_flow_small_increments}
Deviatoric Stress $q$ results for a nonassociative DruckerPrager 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 substepping 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 elasticplastic 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]{../Figurefiles/GCIUncertainty/NonAssociate.pdf}
\caption{
\label{fig_non_associate_plastic_flow}
Yield surface values of returned stress of nonassociative 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 subincrementation 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 crosspoint
% between the yield surface and the stress prediction is used as the direction
% of the returnmapping 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 returnmapping 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.
% %
% Indepth 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/postprocessing algorithm(s):
% any extra preprocessed or postprocessed value which is used
% for verification must be computed via a highorder 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 NavierStokes 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,
% roundoff 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 stateofart innovative numerical technique.
%
% % =================================================================================
% % =================================================================================
% % =================================================================================
% % =================================================================================
% =================================================================================
% =================================================================================
% =================================================================================
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusion and Final Remarks}
This paper deals with verification of constitutive integration algorithms as
applied to elasticplastic 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 elasticplastic 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 DruckerPrager 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 DruckerPrager material models with
nonassociate 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, authoryear citations
%\bibliographystyle{spmpsci} % mathematics and physical sciences
%\bibliographystyle{spphys} % APSlike 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 elasticplastic 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 elasticplastic 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 elasticplastic 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 elasticplastic 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]{../Figurefiles/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 elasticplastic algorithms,
% the stiffness tensor is calculated after the elasticplastic algorithms.
% \begin{figure}[!htbp]
% \centering
% \includegraphics[width=8cm]{../Figurefiles/algorithm_input_output.pdf}
% \caption{Input and Output of Elastoplastic Algorithms}
% \label{fig_elasticplastic_algo_input_output}
% \end{figure}
% Since the compliance tensor is the inverse of the stiffness tensor,
% compliance tensor is also calculated after the elasticplastic 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 elasticplastic algorithms.
% This becomes a chickenandegg 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_elasticplastic_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 elasticplastic 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 twostep predictorcorrector 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{local64}.
% 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.