%\documentclass[Journal,letterpaper,InsideFigs]{ascelike-new}
\documentclass[Journal,letterpaper]{ascelike-new}

\usepackage{lipsum}
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{lmodern}
\usepackage{graphicx}
\usepackage[figurename=Fig.,labelfont=bf,labelsep=period]{caption}
\usepackage{subfig}
\usepackage{multirow}
\usepackage{amsmath}
\usepackage{newtxtext,newtxmath}
\usepackage[pdfauthor={Han Yang}, 
            colorlinks=true,
            linkcolor=black,
            citecolor=black,
            urlcolor=black]
            {hyperref}

\usepackage{array}
\usepackage{tikz}

\newcommand*\circled[1]{\tikz[baseline=(char.base)]{
            \node[shape=circle,draw,inner sep=1pt] (char) {#1};}}

\newcolumntype{x}[1]{>{\centering}p{#1}}

\NameTag{Yang et al., \today}


\begin{document}

\title{Plastic Energy Dissipation in Pressure-Dependent Materials}

\author[1]{Han Yang}
\author[2]{Hexiang Wang}
\author[3]{Yuan Feng}
% \author[1]{Fangbo Wang, Ph.D.}
% \author[2]{David B McCallen, Ph.D.}
%\author[4,5]{Boris Jeremi{\'c}, Ph.D.}
\author[4]{Boris Jeremi{\'c}, Ph.D.}

\affil[1]{Department of Civil and Environmental Engineering, University of California, Davis, CA, USA.}
\affil[2]{Department of Civil and Environmental Engineering, University of California, Davis, CA, USA.}
\affil[3]{Department of Civil and Environmental Engineering, University of California, Davis, CA, USA.}

\affil[4]{Department of Civil and Environmental Engineering, University of California, Davis, CA, USA.}
%\affil[5]{Earth Science Division, Lawrence Berkeley National Laboratory, Berkeley, CA, USA.}
\affil[ ]{Email: jeremic@ucdavis.edu}


\maketitle


\begin{abstract}
Presented is a thermodynamics-based energy analysis approach for pressure-dependent materials.
%
Formulation of plastic free energy and plastic dissipation for non-associated Drucker-Prager plasticity model is derived based on thermodynamics.
%
It is proven that the proposed energy computation formulation always gives non-negative incremental plastic dissipation, as required by the second law of thermodynamics.
%
Presented methodology is illustrated using numerical simulations of Toyoura sand and Sacramento river sand under different loading conditions.
%
Multi-directional loading and pressure-dependency effects on plastic dissipation are investigated.
%
The continuous, non-negative dissipation of mechanical energy in pressure-dependent frictional materials under complex 3D cyclic loading is properly modeled.
\end{abstract}


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



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


Energy dissipation analysis has gained popularity in recent studies of response of dynamic/static inelastic systems.
%
Plastic energy dissipation, if correctly modeled, can be used as an effective indicator of material damage.
%
It is important to distinguish and properly model different energy dissipation mechanisms.
%
These include plastic energy dissipation, viscous damping, and algorithmic (or numerical) damping for numerical modeling of inelastic material \cite{Yang2018b}.
%
Plastic energy dissipation is defined as the thermal energy irreversibly transformed from mechanical energy during a dissipative process.
%
A thermodynamics-based framework is developed and illustrated to model the plastic energy dissipation for pressure-dependent inelastic materials.



Although energy dissipation has been used to explain behavior of soil and structural systems, there exists a common misconception about plastic work and plastic energy dissipation, especially in the field of structural and geotechnical engineering.
%
Plastic work is defined as the material work done due to plastic deformation.
%
In the classic elastoplasticity theory, the increment of plastic work over a time step is calculated by multiplying the effective stress tensor with the increment of plastic strain tensor, $d W^{pl} = \sigma_{ij} d \epsilon_{ij}^{pl}$.
%
As pointed out by \citeN{Collins2003}, this misconception is originated from the decades-old view that, in granular materials, all permanent deformation contributes to the frictional slipping between particles, and thus all the plastic work is dissipated \cite{Luong1986,Okada1994}.
%
% This view is apparently supported by a number of experimental and theoretical studies \cite{Thurairajah1961, Luong1986, Okada1994}.
%
However, a closer examination of these publications reveals that plastic dissipation was not quantitatively measured in any of these studies.
%
In other words, the assumption that plastic work equals to plastic dissipation for granular materials has never been validated.


The difference between plastic work and plastic energy dissipation is defined as the plastic free energy, also known as stored plastic work or cold work.
%
Plastic free energy developed in inelastic material during plastic loading has been observed in physical experiments \cite{Farren1925,Taylor1934,Rittel2000a} and discussed in a number of modeling studies \cite{Collins97,Dafalias1975,Rosakis2000,Collins2002b,Veveakis2007,Feigenbaum2007,Yang2017a}.
%
The physical interpretation of plastic free energy was explained in detail through a conceptual example by \citeN{Yang2017a}.
%
In short, plastic free energy is the part of plastic work that results in the rearrangement of particles, or change in material fabric, rather than the frictional slipping between particles.
%
It is a thermodynamically essential form of energy when a physically discontinuous material is modeled as a continuum.
%
This point will be further pursued in this paper.


In order to model the plastic energy dissipation in a thermodynamically appropriate fashion, general theories of thermomechanics of inelastic materials were established by \citeN{Ziegler1987} and \citeN{local-63}.
%
\citeN{Collins97} reiterated the theory and applied it to the modeling of geotechnical materials.
%
Following studies by Collins et al. \cite{Collins2002a,Collins2002b,Collins2003,Collins2003b,Collins2003c} presented detailed procedures for constructing critical state models that are popular for soils using the thermomechanical framework.
%
One important discovery is that Drucker's postulate, which was considered to be an equivalent condition to the second law of thermodynamics, cannot be used to check the thermodynamic validity of constitutive models.
%
Particularly, the original Cam clay model \cite{Schofield1968} was found to violate the second law of thermodynamics \cite{Collins2002b}.



The constitutive models developed in the previously mentioned studies are significant in the sense of their sound thermodynamics basis.
%
The material response, for example the stress and strain parameters, is controlled by a predefined free energy function and a dissipation function.
%
Compared to the classic elastoplasticity models, these thermomechanical models are more complicated, and thus difficult to be implemented and used in engineering designs.
%
An effort to incorporate the thermomechanical formulation into the classic elastoplasticity theory was made by Feigenbaum and Dafalias \cite{Feigenbaum2007} for von Mises type material models, which is pressure-independent.
%
Yang et al. \cite{Yang2017a} extended the formulation to finite element method (FEM) for the energy analysis of inelastic solids.
%
As a continued work, this study focuses on the plastic energy dissipation analysis of pressure-dependent inelastic materials.



In the following section, the theoretical formulations of thermomechanics and classic elastoplasticity are summarized and discussed.
%
Equations of the plastic free energy and plastic energy dissipation for pressure-dependent material, modeled using Drucker-Prager plasticity, are derived.
%
It is proven that the presented formulation upholds the first and second laws of thermodynamics, for both von Mises and Drucker-Prager plasticity with various hardening rules.
%
Numerical examples on constitutive level and finite element level are used to illustrate the presented energy computation methodology.







%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Plastic Work and Plastic Dissipation}


Thermodynamics-based energy computation formulation for pressure-dependent material is presented and discussed in this section.
%
Material parameters and internal variables from thermomechanical theory are very different than those from the classic elastoplasticity theory.
%
The main challenge is to develop energy equations, in terms of the parameters from classic elastoplasticity, that are within the thermomechanical framework, so that the energy results follow the principles of thermodynamics.



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Thermomechanical Framework}

% Fundamentals of the thermomechanical theory.

% First law of thermodynamics, or energy balance.

% Second law of thermodynamics leads to non-negative energy dissipation.

% Strain decoupling assumption: separate strain (increment) into an elastic part
% and a plastic part.

% Plastic free energy decoupling assumption: separate plastic free energy
% (increment) into an isotropic hardening part and a kinematic hardening part.
% This paper does not consider anisotropic elastoplasticity models, but an
% anisotropic hardening component can be directly added, according to
% \cite{Feigenbaum2007}.


The local incremental form of the first law of thermodynamics for an isothermal process is given by \cite{Collins97,Yang2017a}
%
\begin{equation}
\sigma_{ij} d \epsilon_{ij} = d \psi + \phi
\end{equation}
%
where $\sigma_{ij}$ is the effective stress tensor, $\epsilon_{ij}$ is the total small strain tensor, $\psi$ is the (Helmholtz) free energy density function, and $\phi$ is the incremental plastic energy dissipation density function.
%
The sign convention of stress and strain components follows the traditional mechanics of materials, i.e. positive in tension.



The free energy density function is assumed to be decomposed into an elastic part, known as the elastic strain energy density $\psi^{el}$, and a plastic part, defined as the plastic free energy density $\psi^{pl}$
%
\begin{equation}
d \psi = d \psi^{el} + d \psi^{pl}
\end{equation}
%
This decomposition naturally rises when the material is assumed to be of the decoupled type \cite{Collins97}, which means that the strain tensor can also be additively decomposed into elastic and plastic components.
%
Note that this assumption is also used in the classic small deformation elastoplasticity theory.
%
The incremental elastic strain energy density can then be written as
%
\begin{equation}
d \psi^{el} = \sigma_{ij} d \epsilon_{ij}^{el}
\end{equation}



The plastic free energy density $\psi^{pl}$ is related to the evolution of material model internal variables, and thus further decomposed into parts that correspond to different hardening types \cite{Feigenbaum2007}.
%
In this study, isotropic and kinematic hardening rules of various types are considered.
%
The material model internal variables used in this study are scalar parameter $k$, for isotropic hardening, defined as the size of yield surface in stress space, and back stress tensor $\alpha_{ij}$, for kinematic hardening, defined as the center of yield surface in stress space.
%
The incremental plastic free energy density then becomes
%
\begin{equation}
d \psi^{pl} = d \psi^{iso} (\sigma_{ij}, d \epsilon_{ij}^{pl}, k) + d \psi^{kin} (\sigma_{ij}, d \epsilon_{ij}^{pl}, \alpha_{ij})
\end{equation}
%
% The specific form of the plastic free energy function is discussed below.



Next, the plastic energy dissipation density function is expressed in terms of the plastic free energy density:
%
\begin{equation}
\phi = \sigma_{ij} d \epsilon_{ij} - d \psi = \sigma_{ij} d \epsilon_{ij}^{pl} - (d \psi^{iso} + d \psi^{kin})
\label{equation_incr_plastic_dissipation}
\end{equation}
%
According to the second law of thermodynamics, the incremental plastic dissipation must always be non-negative during any loading increment.
%
When the plastic free energy function is defined, as will be shown in the next sections, the plastic energy dissipation can be calculated from Equation~\ref{equation_incr_plastic_dissipation}.









%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Review of Energy Computation for von Mises Plasticity}


In this section, the equations of plastic free energy and plastic dissipation for pressure-independent von Mises material model \cite{Yang2017a} are revisited.
%
Note that von Mises plasticity always uses associated plastic flow rule.
%
In other words, the plastic flow direction in the stress space is normal to the yield surface.
%
% This means that the plastic strain is purely deviatoric.


The yield function of von Mises is expressed in the following form
%
\begin{equation}
f = \sqrt{ (s_{ij} - \alpha_{ij})(s_{ij} - \alpha_{ij})}- \sqrt{\frac{2}{3}} k
\label{equation_vm_yield_function}
\end{equation}
%
where $s_{ij} = \sigma_{ij} - ({1}/{3}) \sigma_{kk}$ is the deviatoric part of the stress tensor.


Once the material yields, plastic strain starts to develop.
%
The incremental plastic strain tensor is calculated as:
%
\begin{equation}
d \epsilon_{ij}^{pl} = m_{ij} d \lambda
\label{equation_plastic_strain_incr}
\end{equation}
%
where $d \lambda$ is the scalar loading index that equals to the magnitude of incremental plastic strain.
%
Since associated plasticity is used, the normalized plastic flow direction tensor $m_{ij}$ is calculated by taking the gradient of the yield function in the stress space:
%
\begin{equation}
m_{ij} = \frac{\partial f}{\partial \sigma_{ij}} = \frac{(s_{ij} - \alpha_{ij})}{\sqrt{ (s_{mn} - \alpha_{mn})(s_{mn} - \alpha_{mn}) }}
\label{equation_vm_plastic_flow}
\end{equation}


Armstrong-Frederick kinematic hardening \cite{Armstrong66} is a nonlinear strain hardening rule commonly used to model the cyclic inelastic behavior of various types of materials, including metals, alloys, soils, and other structural/geotechnical materials.
%
The evolution of the incremental back stress tensor $d \alpha_{ij}$ is defined
\begin{equation}
d \alpha_{ij} = \left[ \frac{2}{3} h_a m_{ij} - c_r \alpha_{ij} \sqrt{ \frac{2}{3} m_{rs} m_{rs} } \right] d \lambda
\label{equation_vm_af_hardening}
\end{equation}
%
where $h_a$ and $c_r$ are the non-negative hardening constants.
%
When $h_a > 0$ and $c_r = 0$, the nonlinear Armstrong-Frederick hardening becomes linear hardening rule.
%
If $h_a = 0$ and $c_r = 0$, the material model becomes perfectly plastic with no internal variable hardening.
%
Note that isotropic hardening can be defined in a similar form.
%
To avoid repetition, the remaining part of this paper will focus on kinematic hardening.



The plastic free energy density function that was given by Feigenbaum and Dafalias \cite{Feigenbaum2007} and modified by Yang et al. \cite{Yang2017a} is
%
\begin{equation}
d \psi^{kin} = \frac{3}{2 h_a} \alpha_{ij} d \alpha_{ij}
\label{equation_vm_plastic_free_energy}
\end{equation}



The plastic dissipation density function can be expressed in terms of the material parameters and internal variables.
%
Substituting Equations~(\ref{equation_vm_plastic_free_energy}) and (\ref{equation_vm_af_hardening}) into Equation~(\ref{equation_incr_plastic_dissipation}):
%
\begin{equation}
\begin{aligned}
\phi &= \sigma_{ij} d \epsilon_{ij}^{pl} - \frac{3}{2 h_a} \alpha_{ij} d \alpha_{ij}\\
&= s_{ij} m_{ij} d \lambda - \alpha_{ij} m_{ij} d \lambda + \frac{3 c_r}{2 h_a} \sqrt{ \frac{2}{3} m_{rs} m_{rs} } \alpha_{ij} \alpha_{ij} d \lambda
\label{equation_vm_incr_plastic_dissipation_1}
\end{aligned}
\end{equation}
%
Note that the last term on the right hand side of Equation~(\ref{equation_vm_incr_plastic_dissipation_1}) is non-negative.
%
Substituting the expression of plastic flow $m_{ij}$, Equation~(\ref{equation_vm_plastic_flow}), into Equation~(\ref{equation_vm_incr_plastic_dissipation_1}):
%
\begin{equation}
\begin{aligned}
\phi &\ge (s_{ij} - \alpha_{ij}) m_{ij} d \lambda \\
&= \frac{(s_{ij} - \alpha_{ij}) (s_{ij} - \alpha_{ij})}{\sqrt{(s_{mn} - \alpha_{mn})(s_{mn} - \alpha_{mn})}} d \lambda \\
&= \sqrt{(s_{ij} - \alpha_{ij})(s_{ij} - \alpha_{ij})} d \lambda \\
&= \sqrt{\frac{2}{3}} k d \lambda \ge 0
\label{equation_vm_incr_plastic_dissipation_2}
\end{aligned}
\end{equation}
%
According to Equation~(\ref{equation_vm_incr_plastic_dissipation_2}), the plastic energy dissipation density in the case of (associated) von Mises plasticity is always non-negative, which means that the energy dissipation computation does follow the second law of thermodynamics.






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Energy Computation for Associated and Non-Associated Drucker-Prager Plasticity}


Drucker-Prager type plasticity is commonly used to model pressure-dependent material behavior.
%
In this section, the plastic free energy and plastic dissipation are derived for both associated and non-associated Drucker-Prager plasticity models.
%
It will be shown that the plastic free energy function needs an additional pressure-related term, so that the plastic dissipation calculated for pressure dependent materials is thermodynamically correct.


% Need an additional term in plastic free energy to ensure non-negative plastic
% dissipation for pressure-dependent material, modeled using Drucker-Prager
% plasticity.

% Show that the plastic work for associated Drucker-Prager is always zero. Use
% both equations and a plot of yield surface and plastic flow direction in true
% stress space. Conclude that associated Drucker-Prager model is
% thermodynamically inappropriate, which concurs with the papers by Collins et
% al. \citep{Collins97, Collins2002b, Collins2003b}

% Then derive the plastic dissipation for non-associated Drucker-Prager model,
% and prove that the incremental plastic dissipation is always non-negative.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Associated Drucker-Prager Plasticity}

The Drucker-Prager yield function is
\begin{equation}
f = \sqrt{(s_{ij} - p \alpha_{ij})(s_{ij} - p \alpha_{ij})} - \sqrt{\frac{2}{3}} k p
\label{equation_df_yield_function}
\end{equation}
%
where $p = - ({1}/{3}) \sigma_{kk}$ is the mean stress, or hydrostatic pressure, applied on the material.
%
The negative sign ensures that the pressure $p$ is positive when the material is under compression.
%
Note that in Drucker-Prager plasticity, the internal variables $k$ and $\alpha_{ij}$ are dimensionless, while they have the dimension of stress in von Mises plasticity.



The associated plastic flow is:
%
\begin{equation}
m_{ij} = \frac{\partial f}{\partial \sigma_{ij}}
=
\frac{(s_{ij} - p \alpha_{ij}) + \frac{1}{3} \delta_{ij} \alpha_{pq} (s_{pq} - p \alpha_{pq})}{\sqrt{(s_{rs} - p \alpha_{rs})(s_{rs} - p \alpha_{rs}) }} + \sqrt{\frac{2}{27}} k \delta_{ij}
\label{equation_dp_associated_plastic_flow}
\end{equation}
%
where $\delta_{ij}$ is the Kronecker delta.
%
Due to the pressure term in yield function, the plastic flow has both deviatoric and volumetric components:
%
\begin{equation}
\begin{aligned}
m_{ij}^{dev} &= \frac{s_{ij} - p \alpha_{ij}}{\sqrt{(s_{rs} - p \alpha_{rs})(s_{rs} - p \alpha_{rs})}}\\
m_{ij}^{vol} &= \frac{\delta_{ij} \alpha_{pq} (s_{pq} - p \alpha_{pq})}{3 \sqrt{(s_{rs} - p \alpha_{rs})(s_{rs} - p \alpha_{rs})} } + \sqrt{\frac{2}{27}} k \delta_{ij}
\label{equation_dp_associated_plastic_flow_decomposition}
\end{aligned}
\end{equation}



Equations~(\ref{equation_df_yield_function}) and (\ref{equation_dp_associated_plastic_flow}) show that the kinematic hardening for associated/non-associated Drucker-Prager plasticity is of the rotational type, which means that the cone representing the yield function in stress space rotates around the origin, as the back stress $p \alpha_{ij}$ evolves.
%
Figure~\ref{figure_associated_drucker_prager_model} illustrates the Drucker-Prager yield surface with associated plastic flow and rotational kinematic hardening in stress space.
%
\begin{figure}[!htb]
\centering
%\includegraphics[width=0.7\columnwidth]{\Figures/Associated_Drucker_Prager.pdf}
\includegraphics[width=0.7\columnwidth]{Figure01.pdf}
\caption{\label{figure_associated_drucker_prager_model}
Drucker-Prager yield surface with associated plastic flow and rotational kinematic hardening in stress space.}
\end{figure}
%
Note that the stress tensor $\sigma_{ij}$ and the plastic flow $m_{ij}$, or incremental plastic strain tensor $\epsilon_{ij}^{pl}$, are always orthogonal.




Using Equations~(\ref{equation_plastic_strain_incr}) and (\ref{equation_dp_associated_plastic_flow_decomposition}), the incremental plastic work for associated Drucker-Prager plasticity becomes
%
\begin{equation}
\begin{aligned}
\sigma_{ij} d \epsilon_{ij}^{pl} &= (s_{ij} - p \delta_{ij})(m_{ij}^{dev} + m_{ij}^{vol}) d \lambda \\
&= (s_{ij} m_{ij}^{dev} - p m_{ii}^{vol}) d \lambda \\
&= \left[ \frac{s_{ij} (s_{ij} - p \alpha_{ij})}{\sqrt{(s_{rs} - p \alpha_{rs})(s_{rs} - p \alpha_{rs})}} - \frac{p \alpha_{pq} (s_{pq} - p \alpha_{pq})}{\sqrt{(s_{rs} - p \alpha_{rs})(s_{rs} - p \alpha_{rs})}} - \sqrt{\frac{2}{3}} k p \right] d \lambda \\
&= \left[ \sqrt{(s_{ij} - p \alpha_{ij})(s_{ij} - p \alpha_{ij})} - \sqrt{\frac{2}{3}} k p \right] d \lambda = 0
\label{equation_dp_plastic_work}
\end{aligned}
\end{equation}
%
Equation~\ref{equation_dp_plastic_work} proves that the incremental plastic work for associated Drucker-Prager plasticity is always zero.
%
This conclusion is consistent with the orthogonality between the stress tensor $\sigma_{ij}$ and the incremental plastic strain tensor $\epsilon_{ij}^{pl}$, as shown in Figure~\ref{figure_associated_drucker_prager_model}.
%
However, plastic work is expected to evolve as plastic strain develops in an inelastic material.
%
Thus, associated Drucker-Prager plasticity is thermodynamically inappropriate to be used to model pressure-dependent inelastic materials.







%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsubsection{Non-Associated Drucker-Prager Plasticity}

%
As stated by \cite{Collins97}, non-associated plastic flow rule comes naturally for a pressure-dependent frictional material.
%
In this section, non-associated Drucker-Prager plasticity, expressed in the classic elastoplasticity form, is discussed from the perspective of energy dissipation.


One form of the non-associated plastic flow for pressure-dependent sand material was given by \cite{Manzari97}:
%
\begin{equation}
m_{ij}
=
\frac{s_{ij} - p \alpha_{ij}}{\sqrt{(s_{rs} - p \alpha_{rs})(s_{rs} - p \alpha_{rs})}} - \frac{1}{3} D \delta_{ij}
\label{equation_dp_nonassociated_plastic_flow}
\end{equation}
%
with
%
\begin{equation}
D = \xi \left( \sqrt{\frac{2}{3}} k_{d} - \frac{\sqrt{s_{mn} s_{mn}}}{p} \right)
\label{equation_dp_nonassociated_plastic_flow_D}
\end{equation}
%
where $\xi$ and $k_{d}$ are material model constants that controls the volumetric part of the plastic flow.
%
Note that the plastic flow becomes purely deviatoric when the constant $\xi = 0$.



During loading, both the deviatoric and volumetric components of the plastic flow will contribute to the evolution of plastic work.
%
It has been discussed by \cite{Palmer1967,Jefferies1997,Collins2003c} that the plastic strain due to isotropic compression leads to the change of fabric in granular materials.
%
This means that the volumetric part of the plastic strain should be related to the rise of plastic free energy.
%
Thus, the incremental plastic free energy for Drucker-Prager plasticity can be calculated from
%
\begin{equation}
d \psi^{pl} = \left( \frac{3}{2 h_a} \alpha_{ij} d \alpha_{ij} - m_{ii}^{vol} d \lambda \right) p
\label{equation_dp_plastic_free_energy}
\end{equation}
%
Note that the main differences between Equation~\ref{equation_dp_plastic_free_energy}, for Drucker-Prager plasticity, and Equation~\ref{equation_vm_plastic_free_energy}, for von Mises plasticity, are the pressure dependency and an additional term for volumetric plastic flow.



Armstrong-Frederick nonlinear kinematic hardening is considered again:
%
\begin{equation}
d \alpha_{ij} = \left( \frac{2}{3} h_a m_{ij}^{dev} - c_r \alpha_{ij} \sqrt{\frac{2}{3} m_{rs}^{dev} m_{rs}^{dev}} \right) d \lambda
\label{equation_dp_af_hardening}
\end{equation}
%
Note that the evolution of the internal variable $\alpha_{ij}$ is only related to the deviatoric part of the plastic flow.
%
As a result, the internal variable $\alpha_{ij}$ is a deviatoric tensor.



Combining Equations~(\ref{equation_incr_plastic_dissipation}), (\ref{equation_dp_plastic_free_energy}), and (\ref{equation_dp_af_hardening}), the incremental plastic dissipation density for non-associated Drucker-Prager plasticity with Armstrong-Frederick kinematic hardening can be calculated as the following:
%
\begin{equation}
\begin{aligned}
\phi &= (s_{ij} m_{ij}^{dev} - p m_{ii}^{vol}) d \lambda - \left( \frac{3}{2 h_a} \alpha_{ij} d \alpha_{ij} - m_{ii}^{vol} d \lambda \right) p \\
&= (s_{ij} - p \alpha_{ij}) m_{ij}^{dev} d \lambda + \frac{3 c_r}{2 h_a} \sqrt{\frac{2}{3} m_{rs}^{dev} m_{rs}^{dev}} \alpha_{ij} \alpha_{ij} p d \lambda \\
&\ge \sqrt{(s_{ij} - p \alpha_{ij})(s_{ij} - p \alpha_{ij})} d \lambda \\
&= \sqrt{\frac{2}{3}} k p d \lambda \ge 0
\label{equation_dp_incr_plastic_dissipation}
\end{aligned}
\end{equation}
%
According to Equation~(\ref{equation_dp_incr_plastic_dissipation}), the plastic energy dissipation density in the case of non-associated Drucker-Prager plasticity is always non-negative.
%
This is important because Equation~(\ref{equation_dp_incr_plastic_dissipation}) shows that the presented energy calculation methodology for pressure-dependent inelastic material follows the second law of thermodynamics.
















%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Numerical Examples}

Constitutive level numerical examples are used to illustrate the presented energy analysis methodology.
%
Simulations presented in this paper were conducted using the Real-ESSI Simulator \cite{Real_ESSI_Simulator}, a software, hardware and documentation system for high performance, time domain, linear or nonlinear/inelastic, deterministic or probabilistic, finite element modeling and simulation of soil, structure, and their interaction (\url{http://real-essi.info/}).
%
Constitutive level integrations were performed using the backward Euler algorithm, ensuring the convergence of stress and yield function.
%
Strain-controlled loading was used in all examples.


Pressure-dependent materials modeled using non-associated Drucker-Prager plasticity under different loading conditions are investigated.
%
The size of the yield cone is chosen to be small so that the material begins to yield, and dissipate energy, at a low shear strain level.
%
Similar assumptions have been made in a number of constitutive models for pressure-dependent materials \cite{Manzari97,Taiebat2007a,Pisano2012}.
%
Armstrong-Frederick kinematic hardening, defined by Equation~\ref{equation_dp_af_hardening}, is used to model strain-hardening/softening behavior.
%






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Undrained and Drained Triaxial Tests}

The first set of numerical tests are conducted for the purpose of model parameter calibration.
%
Two widely-accepted triaxial experiments on Toyoura sand by \citeN{Verdugo1996} and Sacramento river sand by \citeN{Lee1967}, respectively, are used.
%
The experimental data is obtained from a paper by \citeN{Taiebat2007a}.
%
In this section, the strain-controlled triaxial loading scheme proposed by \citeN{Bardet91} is used for all simulation cases.



Toyoura sand is a cohesionless soil consisting of sub-round to sub-angular quartzite particles.
%
\citeN{Verdugo1996} reported a series of undrained and drained triaxial tests on isotropically consolidated Toyoura sand samples.
%
In this study, the undrained test results of samples under 100 kPa, 1000 kPa, and 2000 kPa confining pressure are used for model parameter calibration.
%
Figure~\ref{figure_undrained_triaxial} shows the comparison between experimental and numerical undrained triaxial test results on Toyoura sand.
%
\begin{figure}[!htbp]
\centering
% \subfloat[]{
%\includegraphics[width=0.49\columnwidth]{Figures/Undrained_Triaxial/p-q.pdf}
\includegraphics[width=0.49\columnwidth]{Figure02a.pdf}
%\includegraphics[width=0.48\columnwidth]{Figures/Undrained_Triaxial/q.pdf}
\includegraphics[width=0.49\columnwidth]{Figure02b.pdf}
% }
\caption{\label{figure_undrained_triaxial}
Comparison between experimental and numerical undrained triaxial test results on Toyoura sand. Experimental data after \protect\citeN{Verdugo1996}.}
\end{figure}
%
It can be seen that the calibrated numerical model can represent the material behavior quite well in $p$-$q$ space, while not that good in $q$-$\epsilon_{axial}$ space, due to simplicity of the used models.
%
The transition from compressive behavior to dilative behavior, which is a key feature of granular materials, is properly modeled.
%
It should be mentioned that the material model used here, non-associated Drucker-Prager plasticity with Armstrong-Frederick kinematic hardening, is simplistic and used to illustrate energy dissipation calculations and not to perfectly match material response.




The experimental results of Sacramento river sand conducted by \citeN{Lee1967} have been used to calibrate and validate constitutive models in a number of studies \cite{Bardet91,Taiebat2007a,Ching2016b}.
%
In this paper, the drained triaxial test results of Sacramento river sand under 290 kPa, 590 kPa, and 1030 kPa confining pressure, respectively, are used to calibrate numerical model parameters.
%
Figure~\ref{figure_drained_triaxial} shows the experimental and numerical drained triaxial test results on Sacramento river sand.
%
\begin{figure}[!htbp]
\centering
% \subfloat[]{
%\includegraphics[width=0.47\columnwidth]{Figures/Drained_Triaxial/Volumetric_Strain.pdf}
\includegraphics[width=0.47\columnwidth]{Figure03a.pdf}
%\includegraphics[width=0.49\columnwidth]{Figures/Drained_Triaxial/Deviatoric_Stress.pdf}
\includegraphics[width=0.49\columnwidth]{Figure03b.pdf}
% }
\caption{\label{figure_drained_triaxial}
Comparison between experimental and numerical drained triaxial test results on Sacramento river sand. Experimental data after \protect\citeN{Lee1967}.}
\end{figure}
%
It is observed that the volumetric strain behavior and deviatoric stress response of the numerical tests correspond well with those from the physical experiments.
%
For the calibrated parameters, the numerical model shows particularly good performance for the samples under low confining pressure.
%
When the confining pressure is relatively high, the numerical results are still acceptable, especially for small strains.




Table~\ref{table_model_parameters} shows the calibrated parameters for the material models, which are implemented in Real-ESSI, used in this study.
%
\begin{table}[!htb]
\centering
\caption{Material model parameters of the pressure-dependent materials used in this study.}
\label{table_model_parameters}
\small
\begin{tabular}{lcx{1.8cm}x{1.8cm}}
\hline\hline
\multirow{2}{*}{Parameter} & \multirow{2}{*}{Unit} & \multicolumn{2}{c}{Material} \tabularnewline \cline{3-4} 
                     &                   & Toyoura & Sacramento \tabularnewline \hline
\texttt{mass\_density} ($\rho$) & $kg/m^3$ & 2000  & 2000 \tabularnewline
\texttt{elastic\_modulus} ($E$) & $MPa$ & 25.0  & 150.0  \tabularnewline
\texttt{poisson\_ratio} ($\nu$) & & 0.3  & 0.3  \tabularnewline
\texttt{druckerprager\_k} &  & 0.107  & 0.107  \tabularnewline
\texttt{armstrong\_frederick\_ha} ($h_a$) & $MPa$  & 17.5 & 45.0  \tabularnewline
\texttt{armstrong\_frederick\_cr} ($c_r$) &  & 150  & 300  \tabularnewline 
\texttt{plastic\_flow\_xi} ($\xi$) &  & 1.9 & 0.7  \tabularnewline
\texttt{plastic\_flow\_kd} ($k_d$) &  & 0.92 & 0.90  \tabularnewline 
\hline\hline
\end{tabular}
\end{table}
%
Notice that the main differences between the two material models are the non-associated plastic flow parameters $\xi$ and $k_d$, as well as the hardening parameters $h_a$ and $c_r$.
%
This is because a small elastic region was chosen so that the post-yield behavior of the numerical model is dominated by plastic flow and hardening.
%
These two sets of parameters are used in cases presented in the following sections with various loading conditions, in order to investigate the plastic energy dissipation for pressure-dependent materials.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Uniaxial Monotonic and Cyclic Shear Loading}

For Drucker-Prager plasticity, deviatoric loading will lead to yielding and plastic energy dissipation.
%
Energy dissipation caused by uniaxial monotonic shear and cyclic shear loading are investigated in this section.
%
Since the material model is pressure-dependent, the relationship between initial confining pressure $p_0$, shear stress evolution, and plastic energy dissipation is of particular interest.
%
Note that the shear strain in the loading direction follows a linear monotonic or cyclic path, while all other strain components are fixed to simulate undrained shearing condition, as shown in Figure~\ref{figure_uniaxial_shear}~(a).
%
This means that the volumetric strain of the material remains constant, and dilative material behavior will lead to increase in confining pressure.
%
\begin{figure}[!htbp]
\centering
%\subfloat[]{\includegraphics[width=0.3\columnwidth]{Figures/Uniaxial_Shear.pdf}}
\subfloat[]{\includegraphics[width=0.3\columnwidth]{Figure04a.pdf}}
\hspace{2cm}
%\subfloat[]{\includegraphics[width=0.3\columnwidth]{Figures/Biaxial_Shear.pdf}}
\subfloat[]{\includegraphics[width=0.3\columnwidth]{Figure04b.pdf}}
\caption{\label{figure_uniaxial_shear}
Constitutive level strain-controlled, undrained shear loading setup: (a) Uniaxial shear loading; (b) Biaxial shear loading.}
\end{figure}


Figure~\ref{figure_monotonic_shear} shows the stress-strain responses and energy dissipation results for Toyoura sand and Sacramento river sand material models when uniaxial monotonic shearing is applied.
%
\begin{figure}[!htbp]
\centering
\subfloat[]{
%\includegraphics[width=0.35\columnwidth]{Figures/Monotonic_Shear/Toyoura/Stress-Strain.pdf}
%\includegraphics[width=0.35\columnwidth]{Figures/Monotonic_Shear/Toyoura/Confining_Pressure.pdf}
%\includegraphics[width=0.35\columnwidth]{Figures/Monotonic_Shear/Toyoura/Energy_Results.pdf}
\includegraphics[width=0.35\columnwidth]{Figure05a.pdf}
\includegraphics[width=0.35\columnwidth]{Figure05b.pdf}
\includegraphics[width=0.35\columnwidth]{Figure05c.pdf}
}\\
\subfloat[]{
%\includegraphics[width=0.35\columnwidth]{Figures/Monotonic_Shear/Sacramento/Stress-Strain.pdf}
%\includegraphics[width=0.35\columnwidth]{Figures/Monotonic_Shear/Sacramento/Confining_Pressure.pdf}
%\includegraphics[width=0.35\columnwidth]{Figures/Monotonic_Shear/Sacramento/Energy_Results.pdf}
\includegraphics[width=0.35\columnwidth]{Figure05d.pdf}
\includegraphics[width=0.35\columnwidth]{Figure05e.pdf}
\includegraphics[width=0.35\columnwidth]{Figure05f.pdf}
}
\caption{\label{figure_monotonic_shear}
Stress-strain responses and plastic energy dissipation results of pressure-dependent materials under uniaxial monotonic shear loading: (a) Toyoura sand; (b) Sacramento river sand.}
\end{figure}
%
The two materials share very similar shear stress evolution and energy dissipation patterns.
%
Both materials are dilative at large shear strains, which means that the confining pressure keeps increasing as the shearing progresses.
%
This leads to the observed continuous increase of shear stress even after the kinematic hardening internal variable $\alpha_{ij}$ reaches saturation.
%
Also, as expected, the sample under a higher level of confinement develops a larger shear stress.


The plastic dissipation density plots in Figure~\ref{figure_monotonic_shear} present an interesting relationship between plastic dissipation and confining pressure.
%
When the initial confining pressure is increased from 200 kPa to 1000 kPa, more energy dissipation is observed.
%
% This is not a surprise since a higher shear stress is causing energy dissipation in the 1000 kPa case.
%
However, when the confining pressure is raised to 2000 kPa, the plastic dissipation density at low strain level is observed to be smaller than those in the other two cases.
%
This is because when a pressure-dependent material is under a larger confinement, it can resist a higher level of shear stress before significant yielding.
%
Then, after the shear strain becomes large enough, the more-confined material yields and dissipates energy at a higher rate due to its larger shear stress.
%
Such energy dissipation feature could be important when modeling pressure-dependent materials, including soils, mine tailings, and other granular materials.




Figure~\ref{figure_cyclic_shear} shows the stress-strain responses and energy dissipation results of the Toyoura sand and Sacramento river sand under uniaxial cyclic shear loading.
%
\begin{figure}[!htbp]
\centering
\subfloat[]{
%\includegraphics[width=0.35\columnwidth]{Figures/Cyclic_Shear/Toyoura/Stress-Strain.pdf}
%\includegraphics[width=0.34\columnwidth]{Figures/Cyclic_Shear/Toyoura/Confining_Pressure.pdf}
%\includegraphics[width=0.34\columnwidth]{Figures/Cyclic_Shear/Toyoura/Energy_Results.pdf}
\includegraphics[width=0.35\columnwidth]{Figure06a.pdf}
\includegraphics[width=0.34\columnwidth]{Figure06b.pdf}
\includegraphics[width=0.34\columnwidth]{Figure06c.pdf}
}\\
\subfloat[]{
%\includegraphics[width=0.35\columnwidth]{Figures/Cyclic_Shear/Sacramento/Stress-Strain.pdf}
%\includegraphics[width=0.34\columnwidth]{Figures/Cyclic_Shear/Sacramento/Confining_Pressure.pdf}
%\includegraphics[width=0.34\columnwidth]{Figures/Cyclic_Shear/Sacramento/Energy_Results.pdf}
\includegraphics[width=0.35\columnwidth]{Figure06d.pdf}
\includegraphics[width=0.34\columnwidth]{Figure06e.pdf}
\includegraphics[width=0.34\columnwidth]{Figure06f.pdf}
}
\caption{\label{figure_cyclic_shear}
Stress-strain responses and energy computation results of pressure-dependent materials under uniaxial cyclic shear loading: (a) Toyoura sand; (b) Sacramento river sand.}
\end{figure}
%
The responses of the two materials again have very similar patterns under cyclic shear loading.
%
The magnitude of shear strain is increased for each loading cycle, and as a result, the shear stress keeps growing.


The plastic dissipation density rate is always positive, which means that the presented energy analysis methodology follows the second law of thermodynamics, as proven in Equation~\ref{equation_dp_incr_plastic_dissipation}.
%
Notice that the strain energy density is not zero at the beginning due to initial confinement.
%
For the first few cycles with small shear strain magnitudes, the accumulated plastic dissipation density is smaller than strain energy density.
%
Then as the applied shear strain increases, plastic dissipation density rapidly increases and surpasses strain energy density.
%
This means that the majority of input work will be dissipated due to material inelasticity when an inelastic material is loaded with a number of deviatoric loading cycles.






%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Biaxial Shear Loading}

The next example focuses on the material response and plastic energy dissipation when biaxial shear loading is applied.
%
The strain-controlled loading setup of biaxial shear is shown in Figure~\ref{figure_uniaxial_shear}~(b).
%
Compared to uniaxial loading, biaxial shearing condition is one step closer to the realistic, fully three-dimensional loading condition.
%
In addition, biaxial shearing test is a common type of laboratory experiments on granular materials.



The previous examples have shown that the two materials, Toyoura sand and Sacramento river sand, share similar mechanical and energy responses when loaded in shear.
%
Therefore, only Toyoura sand is investigated in this section.
%
Figure~\ref{figure_biaxial_shear} shows the stress-strain responses and energy computation results for the material under biaxial shear loading.
%
\begin{figure}[!htbp]
\centering
%\includegraphics[width=0.49\columnwidth]{Figures/Biaxial_Shear/Energy_Results_Modified.pdf}
%\hfill
%\includegraphics[width=0.46\columnwidth]{Figures/Biaxial_Shear/Principal_Stress_Modified.pdf}
%\\
%\includegraphics[width=0.49\columnwidth]{Figures/Biaxial_Shear/Stress-Strain_xy_Modified.pdf}
%\includegraphics[width=0.49\columnwidth]{Figures/Biaxial_Shear/Stress-Strain_xz_Modified.pdf}
%
\includegraphics[width=0.49\columnwidth]{Figure07a.pdf}
\hfill
\includegraphics[width=0.46\columnwidth]{Figure07b.pdf}
\\
\includegraphics[width=0.49\columnwidth]{Figure07c.pdf}
\includegraphics[width=0.49\columnwidth]{Figure07d.pdf}
%
\caption{\label{figure_biaxial_shear}
Stress-strain responses and energy computation results of the Toyoura sand material under biaxial shear loading.}
\end{figure}
%
As shown in the shear strain path in Figure~\ref{figure_biaxial_shear}, the full loading cycle consists of four loading branches:
\begin{enumerate}
\item[\circled{1}] Increase shear strain $\epsilon_{xy}$ from 0 to 5\%;
\item[\circled{2}] Increase shear strain $\epsilon_{xz}$ from 0 to 5\%;
\item[\circled{3}] Decrease shear strain $\epsilon_{xy}$ from 5\% to 0;
\item[\circled{4}] Decrease shear strain $\epsilon_{xz}$ from 5\% to 0.
\end{enumerate}


The initial confining stress is 1000 kPa.
%
Due to material dilatancy and the constant volume loading condition, the confining pressure evolves as the shearing progresses.
%
This is the reason why the stress state of the material at the end of a full loading cycle is different than its original hydrostatic state before loading, as can be observed in the stress path in the principal stress space. 
%
From the shear stress-shear strain plots in the $xy$ and $xz$ directions shown in Figure~\ref{figure_biaxial_shear}, it is seen that shear loading in one direction does influence the stress response in the other direction.
%
For example, during loading branch \circled{2}, the increase of shear strain in the $xz$ direction not only leads to an increase of shear stress in the same direction, but also caused the shear stress in the $xy$ direction to drop.



For pressure-dependent frictional material under deviatoric loading, evolution of fabric and dissipative slipping between particles are always occurring, even during unloading.
%
According to the plastic dissipation density plot, positive dissipation rates are observed in all four loading branches.
%
% This means that the proposed energy computation methodology properly captures the continuous dissipation of mechanical energy in pressure-dependent frictional materials.





% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% \subsection{Earthquake Loading}

% The final example is intended to illustrate how the proposed energy analysis methodology is used in realistic large-scale simulations of nonlinear inelastic systems.
% %
% As shown in Figure~\ref{figure_npp_model}, the finite element model used in this section is a nuclear power plant model that consists of a superstructure, a foundation, several soil layers, and a contact layer between soil and structure.
% %
% \begin{figure}[!htb]
% \centering
% \subfloat[]{
% \includegraphics[width=0.7\columnwidth]{Figures/NPP_Model.pdf}
% \label{figure_npp_model}
% }\\
% \subfloat[]{
% \includegraphics[width=0.32\columnwidth]{Figures/Motion_Acceleration_X.pdf}
% \includegraphics[width=0.32\columnwidth]{Figures/Motion_Acceleration_Y.pdf}
% \includegraphics[width=0.32\columnwidth]{Figures/Motion_Acceleration_Z.pdf}
% \label{figure_npp_motion_accel}
% }\\
% \subfloat[]{
% \includegraphics[width=0.32\columnwidth]{Figures/Motion_Displacment_X.pdf}
% \includegraphics[width=0.32\columnwidth]{Figures/Motion_Displacment_Y.pdf}
% \includegraphics[width=0.32\columnwidth]{Figures/Motion_Displacment_Z.pdf}
% \label{figure_npp_motion_disp}
% }
% \caption{
% The nuclear power plant model and 3D seismic motion used in this study: (a) Finite element model from \protect\citeN{Sinha2017b}; (b) Acceleration time series; (c) Displacement time series.}
% \end{figure}
% %
% The 3D seismic motions developed by \citeN{Rodgers2017} in SW4 \cite{Petersson2017} is applied using domain reduction method \cite{Bielak2001}.
% %
% Figure~\ref{figure_npp_motion_accel} and \ref{figure_npp_motion_disp} shows the 3D acceleration and displacement time series of the motion at the bottom center of the model.
% %
% % \begin{figure}[!htbp]
% % \centering
% % \includegraphics[width=0.32\columnwidth]{Figures/Motion_Acceleration_X.pdf}
% % \includegraphics[width=0.32\columnwidth]{Figures/Motion_Acceleration_Y.pdf}
% % \includegraphics[width=0.32\columnwidth]{Figures/Motion_Acceleration_Z.pdf}
% % \\
% % \includegraphics[width=0.32\columnwidth]{Figures/Motion_Displacment_X.pdf}
% % \includegraphics[width=0.32\columnwidth]{Figures/Motion_Displacment_Y.pdf}
% % \includegraphics[width=0.32\columnwidth]{Figures/Motion_Displacment_Z.pdf}
% % \caption{\label{figure_npp_motion}
% % Seismic motion applied to the nuclear power plant model: (a) Acceleration time series; (b) Displacement time series.}
% % \end{figure}
% % %
% More details of the model and simulation results were presented in \citeN{Sinha2017b}.





% This paper focuses on the material response at the constitutive level, so two locations at different depth are selected from the large-scale finite element model, as indicated in Figure~\ref{figure_npp_model}.
% %
% Note that the simulations presented in \citeN{Sinha2017b} used von Mises plasticity model, which is pressure-independent.
% %
% In this study, the Toyoura sand material shown in Table~\ref{table_model_parameters} is used in the constitutive simulations at the two selected locations, in order to illustrate the effects of pressure-dependency on the material responses.



% Figure~\ref{figure_npp_surface}
% %
% \begin{figure}[!htbp]
% \centering
% \includegraphics[width=0.49\columnwidth]{Figures/NPP_Surface/Stress-Strain.pdf}
% \includegraphics[width=0.49\columnwidth]{Figures/NPP_Surface/Energy_Results.pdf}
% \\
% \includegraphics[width=0.46\columnwidth]{Figures/NPP_Surface/Principal_Strain.pdf}
% \includegraphics[width=0.46\columnwidth]{Figures/NPP_Surface/Principal_Stress.pdf}
% \caption{\label{figure_npp_surface}
% Stress-strain responses and energy computation results of the Toyoura sand material near the soil surface of a nuclear power plant model under earthquake loading.}
% \end{figure}
% %
% shows the stress-strain responses and energy computation results of the Toyoura sand material near the surface of a nuclear power plant model under earthquake loading.
% %
% Since this location is near the soil surface, the confinement caused by the self-weight of the model is relatively small, which means that the material has a low yield stress.
% %
% According to the plot of the three shear stresses $\sigma_{12}$, $\sigma_{13}$, and $\sigma_{23}$, the material yields at a small shear strain.
% %
% Numerous fluctuations in the stress responses caused by the earthquake motion are observed.
% %
% This leads to the consecutive small increases in the plastic dissipation density curve.
% %
% As a result, plastic dissipation continuously accumulates throughout the entire seismic event.



% On the other hand, figure~\ref{figure_npp_deep}
% %
% \begin{figure}[!htbp]
% \centering
% \includegraphics[width=0.49\columnwidth]{Figures/NPP_Deep/Stress-Strain.pdf}
% \includegraphics[width=0.49\columnwidth]{Figures/NPP_Deep/Energy_Results.pdf}
% \\
% \includegraphics[width=0.46\columnwidth]{Figures/NPP_Deep/Principal_Strain.pdf}
% \includegraphics[width=0.46\columnwidth]{Figures/NPP_Deep/Principal_Stress.pdf}
% \caption{\label{figure_npp_deep}
% Stress-strain responses and energy computation results of the Toyoura sand material at a depth of 120 m in a nuclear power plant model under earthquake loading.}
% \end{figure}
% %
%  shows the stress-strain responses and energy computation results of the Toyoura sand material at a depth of 120 m below the soil surface.
% %
% In this case, the material is under a much higher level of confining pressure, so that it has a larger shear strength.
% %
% This explains the sharp shape and small area observed in the plot of the shear stresses.
% %
% Meanwhile, plastic energy dissipation only occurs when the seismic motion is large enough, which is when the peak acceleration and displacement are observed at about 10 s.
% %
% Compared with the material near surface, the deep material has a much smaller amount of accumulated plastic dissipation density.


% This difference in plastic dissipation at different depths shows that the presented energy analysis methodology is ideal for pressure-dependent materials, such as soils and other frictional granular materials.
% %
% Unlike the traditional energy calculation approach, which is based on complete stress-strain loops, the method presented in this paper allows continuous computation of energy parameters throughout the entire simulation.
% %
% This is particularly useful in engineering design, since perfect stress-strain loops never happen in realistic engineering problems.
% %
% It is also worth mentioning that the same energy analysis issue for pressure-independent von Mises plasticity models has been discussed and illustrated in an earlier publication \cite{Yang2017a}.


% By calculating energy parameters at all integration points, the presented energy dissipation analysis can be easily applied to the entire finite element model.
% %
% According to the distribution of plastic dissipation density in a large-scale model, locations with high risks of large deformation and material damage can be identified.
% %
% Such information is useful in improving the safety and economy of designs of infrastructure objects.
% %
% This idea is currently being studied and will be presented in future publications.



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


This paper presented a thermodynamics-based energy analysis approach for pressure-dependent materials.
%
Theoretical formulation of plastic energy dissipation in pressure-dependent, non-associated Drucker-Prager plasticity model was derived and discussed.
%
The proposed energy computation method was implemented in the Real-ESSI simulator system, and illustrated on a series of numerical examples.
%
It was also shown that the presented energy analysis approach can be used in large-scale finite element simulation of 3D dynamic inelastic system.


The energy analysis equations were derived based on thermomechanics with proper assumptions.
%
The difference between plastic work and plastic dissipation, as well as the importance of plastic free energy, was highlighted.
%
Pressure-independent von Mises plasticity was examined from the perspective of energy storage and dissipation.
%
It was mathematically proven that the energy formulation from \citeN{Yang2017a} does follow the first and second law of thermodynamics.


Next, Drucker-Prager plasticity model with rotational kinematic hardening was discussed with focus on energy dissipation behavior.
%
A close examination of associated Drucker-Prager plasticity showed that zero plastic work is always obtained due to the orthogonality between stress tensor and incremental plastic strain tensor.
%
Thus, it was concluded that associated Drucker-Prager plasticity is a thermodynamically inappropriate constitutive model.


The plastic free energy function for von Mises plasticity was modified to incorporate the influence of confinement on the energy behavior of non-associated Drucker-Prager material model.
%
Based on experimental and theoretical works by a number of researchers, it was assumed that the volumetric part of the plastic strain is related to the rise of plastic free energy.
%
The energy formulation derived based on this assumption was then proven to always give non-negative incremental plastic dissipation, as required by the second law of thermodynamics.



Presented energy computation approach was illustrated using numerical examples of pressure-dependent material models under different loading conditions.
%
Two sets of model parameters were calibrated using the triaxial test data on Toyoura sand \cite{Verdugo1996} and Sacramento river sand \cite{Lee1967}.
%
Uniaxial shearing examples showed that the proposed analysis method can properly account for the influence of pressure-dependency on the energy dissipation behavior of Drucker-Prager model.
%
In general, it was observed that the majority of input work will be dissipated if significant cyclic deviatoric loading is applied.



In the case of biaxial shear loading, it was observed that the evolution of shear stress in one direction influenced the shear stress response in the other direction, due to the pressure dependency of the material model.
%
This also lead to different energy dissipation behaviors when loads were increased or decreased in different directions.
%
More importantly, it was pointed out that the continuous dissipation of mechanical energy in pressure-dependent frictional materials was properly modeled by the proposed energy analysis methodology.


% The last numerical example presented the application of energy dissipation analysis in realistic earthquake soil structure interaction problems.
% %
% It was shown that the stress-strain behavior and energy dissipation results were very different for the materials at different locations.
% %
% Due to the pressure-dependency of the Drucker-Prager model, the material at larger depth below the soil surface had a higher level of confinement, which lead to a larger yield stress and smaller plastic dissipation density.
% %
% By identifying locations with larger energy dissipation density, which means higher risk of large deformation and material damage, the safety and economy of designs of infrastructure objects can be evaluated and possibly improved.



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

This work was supported in part by the US-DOE.

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


% \section*{References}
\bibliography{refmech}

% \bibliographystyle{abbrvnat}


\end{document}