%%%%%%%%%%%%%%%%%%%%%%% file template.tex %%%%%%%%%%%%%%%%%%%%%%%%%
%
% This is a template file for the global option of the SVJour class
%
% Copy it to a new file with a new name and use it as the basis
% for your article
%
%%%%%%%%%%%%%%%%%%%%%%%% Springer-Verlag %%%%%%%%%%%%%%%%%%%%%%%%%%
%
% First comes an example EPS file -- just ignore it and
% proceed on the \documentclass line
\begin{filecontents*}{example.eps}
%!PS-Adobe-3.0 EPSF-3.0
%%BoundingBox: 19 19 221 221
%%CreationDate: Mon Sep 29 1997
%%Creator: programmed by hand (JK)
%%EndComments
gsave
newpath
  20 20 moveto
  20 220 lineto
  220 220 lineto
  220 20 lineto
closepath
2 setlinewidth
gsave
  .4 setgray fill
grestore
stroke
grestore
\end{filecontents*}
%
% Choose either the first of the next two \documentclass lines for one
% column journals or the second for two column journals.
\documentclass[global,referee]{svjour}
%\documentclass[global,twocolumn,referee]{svjour}
% Remove option referee for final version
%
% Remove any % below to load the required packages
%\usepackage{latexsym}
\usepackage{graphics}
%\usepackage{times}
\usepackage{graphicx}
\usepackage{amsmath}
\usepackage{color}
\usepackage{ulem}
\usepackage[round,colon,authoryear]{natbib}
%\usepackage{subeqn}
%
% NOTE: Don't include the \NameTag{} if you have selected
%       the NoPageNumbers option: this leads to an inconsistency.
%\NameTag{Kuhn}
%

% etc
%
% Insert the name of "your" journal with the command below:
\journalname{Acta Geotechnica}
%
\begin{document}

%
%\title{Probabilistic Elasto-Plasticity:
%Formulation of Evolution Equation of Probability Density Function}
\title{Probabilistic Elasto-Plasticity: Formulation in 1D}
%\subtitle{Do you have a subtitle?\\ If so, write it here}
\author{
Boris Jeremi{\'c}\inst{1}
\and
Kallol Sett\inst{2}
\and
M. Levent Kavvas\inst{3}
}


\institute{
Associate Professor, Department  of  Civil  and  Environmental  Engineering,  University  of
California, Davis, CA 95616. \texttt{jeremic@ucdavis.edu}
\and 
Graduate Research Assistant, Department  of  Civil  and  Environmental  Engineering,  University  of
California, Davis, CA 95616.
\and 
Professor, Department  of  Civil  and  Environmental  Engineering,  University  of
California, Davis, CA 95616.
 }

% \author{First author\inst{1} \and Second author\inst{2}% etc
% % \thanks is optional - remove next line if not needed
% \thanks{\emph{Present address:} Insert the address here if needed}%
% }                     % Do not remove
% %
%\offprints{}          % Insert a name or remove this line
%
% \institute{Insert the first address here \and the second here}
% %
% \date{Received: date / Revised version: date}
% The correct dates will be entered by the editor
%
\maketitle
%
%
% You will need to make the title all-caps
%
% I've found that the \and command doesn't quite work, so just use "and"
%  such as the following (and don't forget the ending curly brace `}').
%\\
%and
%Younyee Kuhn%
%\thanks{Flourishing wife of same.},%
%\ Not a Member, ASCE
 

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

\begin{abstract}
 
A    second-order    exact  expression  for  evolution  of Probability Density  
Function  (PDF)  of  stress  is  derived  for  general,  one  dimensional  (1-D)
elastic-plastic  constitutive rate equations with uncertain material parameters.
The  Eulerian--Lagrangian (EL) form of Fokker--Planck--Kolmogorov (FPK) equation
is  used for this purpose. It is also shown that by using EL form of FPK, the so
called  "closure  problem"  associated with regular perturbation methods used so
far,  is  resolved  too.  The use of EL form of FPK also replaces repetitive and
computationally  expensive deterministic elastic-plastic computations associated
with Monte Carlo technique.

The  derived general expression are  specialized to the particular      cases of
point location scale linear elastic and elastic--plastic constitutive equations,
related to associated Drucker-Prager with linear hardening

In  a  companion  paper,  the  solution  of  FPK  equations for 1D is presented,
discussed and illustrated through number of examples.

\end{abstract}

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


Advanced elasto--plasticity constitutive models, when properly calibrated,
are  very accurate in capturing important aspects of material behavior. However,
all  materials',  and in particular geomaterials' (soil, rock, concrete, powder,
bone  etc.)  behavior is  uncertain due to inherent spatial and
point-wise  uncertainties. These
uncertainties  in  material  properties could  outweigh  the  advantages  gained by using advanced constitutive
models. For example, Fig.~\ref{figure:AnticipatedInfluence} shows a schematic of
anticipated  influence  of  material uncertainties on a bi-linear elastic-plastic
stress-strain behavior. Depending on uncertainties in material properties and
interaction  between  them,  the  behavior  of  the  same material could be very
different. 
%This could be even more complicated for non-linear materials.

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.7\textwidth]{AnticipatedInfluence.eps}
\caption{Anticipated Influence of Material Fluctuations on Stress-Strain Behavior}
\label{figure:AnticipatedInfluence}
\end{center}
\end{figure}


The  uncertainties in material properties are inevitable in real materials  and 
it is best to account for them in
modeling   and   simulation.   In   traditional   deterministic
constitutive  modeling,  material models are calibrated against
set  of  experimental  data.  Although  those  experimental  data sets
generally  exhibit  statistical  distribution,  the  models are
usually  calibrated  against  the  mean of the data and all the
information about  uncertainties  is  neglected.  

The  modeling  and
simulation  of  solids  and  structures with uncertain material
properties   involves   two   steps:   (a)  classification  and
quantification   of   uncertainties   and  (b)  propagation  of
uncertainties through governing differential equations.


The  uncertainties  can  be broadly classified into
aleatory   and  epistemic  types.  Aleatory  uncertainties  are
associated with the inherent variabilities of nature. This type
of   uncertainty   can   not   be   reduced.  Highly  developed
mathematical  theory  is  available  for  dealing with aleatory
uncertainty.  On  the other hand, epistemic uncertainties arise
due to our lack of knowledge.
%
This type of uncertainty can be reduced by collecting more data
but  the  mathematical  tools  to deal with them are not highly
developed  (e.g.  fuzzy logic \citet{Zadeh:1983}, convex models
\citet{book:Elishakoff}, interval arithmetic \citet{book:Moore}
etc.).
%
Hence,  it  proves  useful to trade epistemic uncertainties for
aleatory  uncertainties in order to facilitate their propagation 
through  the  governing  equations  using advanced mathematical
tools.  
%
It  is  important to note that in trading-off epistemic
uncertainties  for  aleatory  uncertainties, one doesn't reduce
the  total  uncertainties  in  the  system, but assumes that the
uncertainties   in   the  system  are  irreducible.  Under  the
framework  of probability theory, uncertain material parameters
are    modeled   as   random   variables   or   random   fields
\citep{book:Vanmarcke}.
%    depending   on   whether   they   are
%specialized  to  a  fixed  location  in  their  continuum  or a
%function  of location in their continuum. 
We note recent works  in  quantifying 
the  uncertainties  in  material  (soil) properties for geotechnical engineering
applications,
%
\citet{Lumb:1966},
\citet{Vanmarcke:1977},
\citet{Mayerhoff:1993},
\citet{Baecher:1993}, 
\citet{Popescu:1995}, 
\citet{Lacasse:1996}, 
\citet{Popescu:1998}, 
\citet{Kulhawy:1999A, Kulhawy:1999B}, 
\citet{Fenton:1999A, Fenton:1999B}, 
\citet{Duncan:2000},
\citet{Rackwitz:2000},
\citet{Marosi:2004}, and
\citet{Stokoe:2004}
%
The   issue  of  uncertain  material  properties  becomes  very
pronounced  when  one  starts  dealing  with the boundary value
problems   with   uncertain  material  properties  (elastic  or
elastic--plastic).




In  mechanics,  the equilibrium equation,$ A \sigma = \phi(t)$,
together   with  the  strain  compatibility  equation, $B  u  =
\epsilon$,  and  the constitutive equation, $\sigma = D \epsilon$,
are    sufficient\footnote{Generalized  stress  is $ \sigma $, $ \phi(t) $
is  generalized  forces  that  can  be time dependent, $ u $ is
generalized  displacements,  $ \epsilon $ is generalized strain
and  $A$,  $B$,  and $D$ are operators which could be linear or
non-linear.}    to    describe   the   behavior   of   the
solid.
%
%
Rigorous  mathematical  theory  has been developed for problems
where   the   only  random  parameter  is  the  external  force
$\phi(t)$.  In this case, the probability distribution function
(PDF)  of  the  response  process will  satisfy  FPK  partial
differential    equation  \citep{Soize1994}   .   With
appropriate  initial and boundary conditions the FPK PDE can be
solved  for  PDF  of  response variable. The numerical solution
method  for  FPK equation by finite element method (FEM) is described
by   number   of   researchers   e.g.  \citet{Langtangen:1991},
\citet{Bergman:2005}.


The  other  extreme case, which is of interest in this work, is
when  the stochasticity of the system is purely due to operator
uncertainty.  Exact  solution  of  the problems with stochastic
operator was attempted by \cite{Hopf:1952} using characteristic
functional   approach.   Later,   \cite{Lee:1974}  applied  the
methodology  to the problem of wave propagation in random media
and  derived  a  FPK  equation  satisfied by the characteristic
functional  of  the  random  wave  field.  This  characteristic
functional approach is very complicated for linear problems and
becomes  even  more  intractable  (and possibly unsolvable) for
nonlinear  problems  and problems with irregular geometries and
boundary conditions.


Monte   Carlo   simulation   technique  is  an  alternative  to
analytical  solution  of  partial  differential  equation  with
stochastic coefficient.
%
Nice  descriptions of different aspects of formulation of Monte
Carlo  technique  for stochastic mechanics problem is described
by \citet{Schueller:1997}.
%
Monte Carlo method is very popular tool with the advantage that
accurate  solution  can  be  obtained  for  any  problem  whose
deterministic  solution  (either  analytical  or  numerical) is
known.
%
Monte  Carlo  technique  has been used by a number of researchers in
obtaining probabilistic solution of geotechnical boundary value
problems,  e.g.
\citet{Griffiths:1996, Griffiths:2002, Griffiths:2003, Griffiths:2005}.
\citet{Popescu:1997},
\citet{Mellah:SFEMgeotechApplication},
\citet{DeLima:SFEMgeotechApplication},
\citet{Deodatis:2002},
\citet{Nobahar:2003}.
%  
The  major  disadvantage  of  Monte  Carlo analysis is the
repetitive  use  of  the deterministic model until the solution
variable become statistically significant. The computational cost
associated with it could be very high especially for non--linear
problems with multiple uncertain material properties.



Various  difficulties   in  finding  analytical   solutions   and   the  high
computational      cost    associated  with  Monte  Carlo  technique  instigated
development  of  numerical method for the solution of
stochastic   differential  equation with random coefficient. For
stochastic  boundary  value  problems Stochastic Finite Element
Method   (SFEM)  is  the  most  popular such method.  There  exist  several
formulations     of     SFEM,    among    which    perturbation
(\citet{book:Kleiber};                \citet{DerKiureghian:1988};
\citet{Mellah:SFEMgeotechApplication};  \citet{DeBorst:1999}) and
Spectral         (\citet{book:Ghanem};        \citet{keese:SPDE};
\citet{Karniadakis:2003};               \citet{Debusschere:2003};
\citet{Hori:2000})  methods  are  very popular. 
%
A nice review on
advantages  and disadvantages of different formulations of SFEM
was  provided  by \citet{Matthies:1997reviewSFEM}. Mathematical
issues  regarding  different formulations of SFEM was addressed
by    \cite{babuska:GalerkinSPDE}  and  \cite{Babuska:SPDE}.  It is important to
note that most  of  the  formulations described in the above mentioned references
are    for  geometrically linear or nonlinear problems with linear
elastic material.
%

A limited number of references is also available related to geometric
non--linear         problems, \citet{DerKiureghian:1991};
\citet{keese:SPDE} and  \citet{keese:2003ReportSFEM}). 
%
Similarly,  there    exist  only        few    published   references related to
material   non--linear   (elastic--plastic)  problems  with  uncertain  material
parameters.
%
The major difficulty in extending the available formulations of
SFEM to general elastic--plastic problem is the high non--linear
coupling in the elastic--plastic constitutive rate equation.
%
First    attempt    to    propagate    uncertainties    through
elastic--plastic   constitutive  equations  considering  random
Young's    modulus    was   published   only   recently,   e.g.
\citet{Hori:1999,  Hori:2000}. 
%
The perturbation expansion
at  the  stochastic mean behavior (considering only the first
term    of  the expansion) was used for the material nonlinear part in the above
mentioned references. In computing
the mean behavior the Authors took the
advantage of bounding media approximation. Although this method
doesn't  suffer  from  computational difficulty associated with
Monte Carlo method for problems having no closed-form solution,
it  inherits  "closure  problem"  and the "small coefficient of
variation" requirements for the material parameters.
%
Closure problem refers to the need for higher order statistical
moments  in  order to calculate lower order statistical moments
\citet{Kavvas:2003}.
%
The  small  COV requirement claims that the perturbation method
can  be  used  (with  reasonable  accuracy)  for  probabilistic
simulations  of solids and structures with uncertain properties
only if their COV $ < 20~\%$ \citep{Kiureghian:2000ReportSFEM}.
For  soils  and  other natural materials, COVs are rarely below
$20~\%$       (\cite{Lacasse:1996};       \cite{Kulhawy:1999A,
Kulhawy:1999B}).
%
Furthermore,  with  bounding  media  approximation,  difficulty
arises  in  computing  the  mean  behavior  when  one considers
uncertainties  in  internal  variable(s) and/or direction(s) of
evolution of internal variable(s).
%
%These difficulties, associated with probabilistic simulation of
%the elastic--plastic constitutive equation, prevent application
%of SFEM to general geomechanics problems.



The focus of present work is on development of methodology for the probabilistic simulation of
constitutive  behavior  of elastic--plastic materials with
uncertain   properties. 
%In addition to that,    and   on  resolving  the  difficulties
%associated with them for general applicability.
%
Recently,     \citet{Kavvas:2003}     obtained     a    generic
Eulerian--Lagrangian  (EL)  form   of   FPK   equation,   exact   to
second-order,   corresponding   to   any  non--linear  ordinary
differential  equation  with  random  coefficients  and  random
forcing.  The approach using EL form of the  FPK  equation doesn't suffer from the
drawbacks of Monte Carlo method and perturbation technique.
%
In  this  paper  the  authors applied developed EL form of the  FPK  equation to
obtain probabilistic formulation  for   a   general,   one-dimensional
incremental  elastic--plastic constitutive equation with random
coefficient.  The solution methodology is designed with several
applications in mind, namely to
%
\begin{itemize}
%
\item  obtain  probabilistic stress--strain behavior from spatial
average  form  (upscaled  form)  of constitutive equation, when
input   uncertain   material  properties  to  the  constitutive
equation are random fields; and
%
\item    obtain   probabilistic   stress-strain   behavior   from
point-location   scale   constitutive   equation,   when  input
uncertain  material properties to the constitutive equation are
random variables.
%
\end{itemize}
%The  objective  of  the second
%application  is to integrate  the  development as
%probabilistic  constitutive  integrator  with  stochastic finite element method,
%where  random  fields  are  discretized  to random variables by standard methods
%(e.g.  by  spectral  method  \citet{book:Ghanem}).  
%
Application  of  the developed methodology is demonstrated on a
particular  point-location  scale  one-dimensional constitutive
equation, namely Drucker--Prager  associative  linear  hardening
elastic--plastic  material  model. 
%
In    this  paper, derivation is made of the EL form of FPK equation that govern
the 1D  probabilistic  elastic--plastic  material  models  with  uncertain
material   parameters.  This  general  formulation  is  then  specialized  to  a
particular 1D Drucker--Prager associative linear
hardening material model.
%
It is noted  that  the  solution  of  FPK equations can be quite computationally
intensive.  However,  in  presented  case,  with  1D constitutive problem, this 
solution poses no significant computational burden. 
%
It  is also noted that the solutions of FPK in stress--strain space
inherently   represents  a  Markov  process. 
In other words, there is no probability of switching to alternate
equilibrium    branch   (say  softening branch with localization of deformation)
as  Markov process requires that paths in probability density space are strictly
followed from initial conditions all through the loading process.
%




In the companion paper the solution methodology of the FPK
equation  corresponding  to  Drucker--Prager associative linear
hardening   material   model  is  described, along with  illustrative
examples.
%
The  methodology  is general enough that it allows extension to
three-dimensions  and  incorporation  into a general stochastic
finite  element  framework.  This  work is underway and will be
reported in future publications.




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


\section{General Formulation}
\label{GeneralFormulation}

The incremental form of  {spatial-average} elastic-plastic constitutive
equation can be written as
%
\begin{equation}
\label{eqno_S2}
\frac{d\sigma_{ij}(x_t, t)}{dt}
=
D_{ijkl}(x_t,t)
\frac{d\epsilon_{kl}(x_t,t)}{dt}
\end{equation}
%
where    $x_t$  is  current  spatial location at time $t$, and the    continuum 
stiffness tensor $D_{ijkl}(x_t,t)$ can be either elastic or elastic-plastic
%
\begin{eqnarray}
D_{ijkl} = \left\{\begin{array}{ll}
%
           D^{el}_{ijkl}
%
%
        \;\;\; & \mbox{~;~} f <  0 \vee (f = 0 \wedge df <  0) \\
%
           \\
%
           D^{el}_{ijkl}
           -
           \displaystyle \displaystyle \frac{D^{el}_{ijmn}
           \displaystyle \frac{\partial U}{\partial \sigma_{mn}}
           \displaystyle \frac{\partial f}{\partial \sigma_{pq}}
           D^{el}_{pqkl}}
           {\displaystyle \frac{\partial f}{\partial \sigma_{rs}}
           D^{el}_{rstu}
           \displaystyle \frac{\partial U}{\partial \sigma_{tu}}
           -
           \displaystyle \frac{\partial f}{\partial q_*}r_*}
%
        \;\;\; & \mbox{~;~} f = 0 \vee df = 0
%
           \end{array} \right.
\label{el_or_elpl}
\end{eqnarray}
%
and where $ D^{el}_{ijkl} $ is the elastic stiffness tensor, 
$D^{ep}_{ijkl} $ is the elastic--plastic continuum stiffness tensor, $f$  is  the  yield
function, which is a function of stress ($\sigma_{ij}$) and
internal  variables  ($q_*$), $U$ is the plastic potential function (also a function
of  stress  and  internal  variables). The
internal  variables  ($q_*$)  could  be  scalar(s)  (for  perfectly-plastic and
isotropic  hardening  models),  second-order  tensor  (for  translational and
rotational   kinematic  hardening)  or  fourth-order  tensor  (for  distortional
hardening). 
%The same classification applies to the  direction
%of  evolution of internal variables ($r_*$). 
Therefore, the most
general  form  of  incremental constitutive equation in terms of its parameters can be
written as
%
\begin{equation}
\frac{d \sigma_{ij} (x_t,t)}{dt}
=
\beta_{ijkl}(\sigma_{ij}, D_{ijkl}, q_*, r_*; x_t,t)
\frac{d \epsilon_{kl} (x_t,t)}{dt}
\label{eqno_S2a}
\end{equation}
%
Due  to  
{randomnesses}  in  elastic  constants  ($D^{el}_{ijkl}$) and internal
variables ($q_*$) and/or rate of evolution of internal variables ($r_*$) the
material    stiffness   operator $\beta_{ijkl}$ in  Eq.~(\ref{eqno_S2a}) becomes
stochastic.  It  follows that the Equation~(\ref{eqno_S2})  becomes  a 
linear/non-linear ordinary
differential  equations  with  stochastic coefficients. 
%
Similarly,  {randomness}        in    he    forcing  term ($ \epsilon_{kl} $) of
Equation~(\ref{eqno_S2a}) results in  Equation~(\ref{eqno_S2a}) becoming linear/non-linear
ordinary  differential  equations with stochastic forcing. 
%
This  can  be  generalized,  so  that  {randomnesses} in material properties and
forcing function of     Equation~(\ref{eqno_S2a}) results in Equation~(\ref{eqno_S2a})
becoming a  linear/non-linear  ordinary  differential  equation  with stochastic
coefficients and stochastic forcing.



In  order  to  gain  better  understanding {of} the effects of 
{random}
material  parameters and forcing on  response, focus is shifted from a general 3D case to a
1D {case}. 
%
In   what   follows,   the  probabilistic  formulation  for  1-D    constitutive
elastic--plastic incremental equation with stochastic coefficient and
stochastic forcing  is  derived. 
%
In  addition  to  that,  the  probabilistic  formulation for 1-D elastic linear 
constitutive equation is obtained as a
special case of non-linear general derivation.


Focusing on 1-D behavior, the Eq.~(\ref{eqno_S2a}) is written as
%
\begin{equation}
\label{eqno_S3}
\frac{d\sigma(x_t, t)}{dt}
=
\beta(\sigma,D,q,r;x_t,t)
\frac{d\epsilon(x_t,t)}{dt}
\end{equation}
%
which is a non-linear ordinary differential equation with stochastic coefficient
and     stochastic   forcing.   The right hand side of Eq.~(\ref{eqno_S3}) is
replaced with the function $\eta$ as
%
\begin{equation}
\eta(\sigma,D,q,r,\epsilon; x, t)
=
\beta(\sigma,D,q,r;x_t,t) \frac{d\epsilon(x_t,t)}{dt}
\label{eqno_S4}
\end{equation}
%
so  that  now Eq.~(\ref{eqno_S3}) can be written as
%
\begin{equation}
\label{eqno_S5}
\frac{\partial \sigma(x_t, t)}{\partial t}= \eta(\sigma,D,q,r,\epsilon; x, t)
\end{equation}
%
with initial condition,
%
\begin{equation}
\label{eqno_S6}
\sigma(x,0)=\sigma_0
\end{equation}

In    the    above Eq.~(\ref{eqno_S5}) $\sigma$ can be considered to represent a
point  in  the  $\sigma$-space and hence, the Eq.~(\ref{eqno_S5}) determines the
velocity  for  the  point in that $ \sigma $-space. This may be visualized, from
the  initial point, and given initial condition $\sigma_0$, as a trajectory that
describes  the  corresponding  solution  of  the  non-linear stochastic ordinary
differential  equation  (ODE)  (Eq.~(\ref{eqno_S5})). Considering now a cloud of
initial  points  (refer  to  Fig.~\ref{figure:CloudOfPoints}),  described  by  a
density $\rho(\sigma,0)$ in the $\sigma$-space.
%

\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.7\textwidth]{CloudOfPoints01.ps}
\caption{Movements of Cloud of Initial Points, described by density $\rho(\sigma,0)$, in the $\sigma$-space}
\label{figure:CloudOfPoints}
\end{center}
\end{figure}
%
The phase density $\rho$ of $\sigma(x,t)$ (movement  of  any   point   dictated 
by  Eq.~(\ref{eqno_S5})) varies in time according to a continuity equation which
expresses the conservation of all these points in the $ \sigma $-space.
%
This  continuity  equation can be expressed in mathematical terms, using Kubo's
stochastic Liouville equation \citep{Kubo:1963}:
%
\begin{equation}
\label{eqno_S7}
\displaystyle \frac{\partial \rho (\sigma(x,t),t)}{\partial t}
=
-\displaystyle \frac{\partial}{\partial \sigma} \eta 
[\sigma(x,t), D(x), q (x), r (x), \epsilon(x,t)].\rho[\sigma(x,t),t]
\end{equation}
%
with an initial condition,
%
\begin{equation}
\label{eqno_S8}
\rho(\sigma,0)=\delta(\sigma-\sigma_0)
\end{equation}
%
where     $\delta$  is the Dirac  delta  function  and  Eq.~(\ref{eqno_S8})  is 
the  probabilistic  restatement  in  the  $\sigma$-phase  space  of the original
deterministic  initial  condition  (Eq.~(\ref{eqno_S6})).  Here it proves useful
to recall  Van Kampen's
Lemma  \citep{VanKampen:1976}, which states that the ensemble average of a phase
density is the probability density
%
\begin{equation}
\label{eqno_S9}
<\rho(\sigma,t)>=P(\sigma,t)
\end{equation}
%
where,    the    symbol    $<\cdot>$   denotes  the  expectation  operation, and
$P(\sigma,t)$  denotes  evolutionary  probability  density of the state variable
$\sigma$ of the constitutive rate equation (Eq.~(\ref{eqno_S3})).


In      order    to   obtain  the deterministic\footnote{The probability density
function  describes  probability  of  occurrence  of  certain  event,  but  is  a
deterministic quantity.} probability density function (PDF)
$(\sigma,t)$  of  the  state  variable,  $\sigma$, it is necessary to obtain the
deterministic  partial  differential  equation  (PDE) of the $\sigma$-space mean
phase   density   $<\rho(\sigma,t)>$  from  the  linear  stochastic  PDE  system
(Eqs.~(\ref{eqno_S7}) and ~(\ref{eqno_S8})). 
%
This   necessitates   the   derivation  of  the    ensemble    average  form  of
Eq.~(\ref{eqno_S7})  for $<\rho (\sigma,t)>$. 
%
%
%
%
%
%
% NEW REWORKED     
% NEW REWORKED      %
% NEW REWORKED      %
% NEW REWORKED      %
% NEW REWORKED      To this end, for writing convenience, denoting the ensemble average of
% NEW REWORKED      the function $\eta$ at time $t$ by $\Xi (x,t)$ as,
% NEW REWORKED      %
% NEW REWORKED      \begin{equation}
% NEW REWORKED      \label{eqno_S10a}
% NEW REWORKED      \Xi (x,t) = \left< \vphantom{\int_{0}^{t} d\tau} \eta(\sigma(x_t,t), D(x_t),
% NEW REWORKED      q(x_t), r(x_t), \epsilon(x_t,t)) \right>
% NEW REWORKED      \end{equation} 
% NEW REWORKED      %
% NEW REWORKED      the  time-ordered covariance\footnote{Defined for two dummy variables $p(x,t_1)$
% NEW REWORKED      and    $p(x,t_2)$   as  $  Cov_0\left[p(x,t_1);
% NEW REWORKED      p(x,t_2)\right]=\left<p(x,t_1)    p(x,t_2)\right>-\left<p(x,t_1)\right>    \cdot
% NEW REWORKED      \left<p(x,t_2)\right> $ }
% NEW REWORKED      %
% NEW REWORKED      between  the  function  $\eta$  at  time  $t$ and the derivative of the function
% NEW REWORKED      $\eta$ with respect to stress, $\sigma$, at time $t-\tau$ by $\Pi(x,t,\tau)$ as,
% NEW REWORKED      %
% NEW REWORKED      \begin{eqnarray}
% NEW REWORKED      \nonumber
% NEW REWORKED      \lefteqn{\Pi (x,t,\tau) 
% NEW REWORKED      = 
% NEW REWORKED      Cov_0\left[\vphantom{\int_{0}^{t} d\tau} \eta(\sigma(x_t,t), D(x_t), q(x_t), r(x_t),
% NEW REWORKED      \epsilon(x_t,t)); \right.} \\
% NEW REWORKED      \label{eqno_S10b}
% NEW REWORKED      && \left. \displaystyle \frac{\partial \eta (\sigma(x_{t-\tau}, t-\tau), D(x_{t-\tau}), q(x_{t-\tau}), 
% NEW REWORKED      r(x_{t-\tau}), \epsilon(x_{t-\tau}, t-\tau))}{\partial \sigma} \right]
% NEW REWORKED      \end{eqnarray}
% NEW REWORKED      %
% NEW REWORKED      and the time-ordered covariance between  function $\eta$ at time $t$ and function
% NEW REWORKED      $\eta$ at time $t-\tau$ by $\Omega(x,t,\tau)$ as,
% NEW REWORKED      %
% NEW REWORKED      \begin{eqnarray}
% NEW REWORKED      \nonumber
% NEW REWORKED      \lefteqn{\Omega (x, t,\tau) = Cov_0\left[\vphantom{\int_{0}^{t} d\tau} \eta(\sigma(x_t,t), D(x_t), q(x_t), r(x_t),
% NEW REWORKED      \epsilon(x_t,t)); \right.} \\
% NEW REWORKED      \label{eqno_S10c}
% NEW REWORKED      && \left. \eta (\sigma(x_{t-\tau}, t-\tau), D (x_{t-\tau}), q(x_{t-\tau}), r(x_{t-\tau}), \epsilon
% NEW REWORKED      (x_{t-\tau}, t-\tau)) \vphantom{\int_{0}^{t} d\tau} \right]
% NEW REWORKED      \end{eqnarray}
% NEW REWORKED      %
% NEW REWORKED      one can write the  ensemble  average  form of Eq.~(\ref{eqno_S7}) as 
% NEW REWORKED      %
% NEW REWORKED      \begin{eqnarray}
% NEW REWORKED      \nonumber
% NEW REWORKED      \lefteqn{\frac{\partial \left<\rho (\sigma(x_t,t), t)\right>}{\partial t} =} \\
% NEW REWORKED      \nonumber
% NEW REWORKED      & - & \displaystyle \frac{\partial}{\partial \sigma}\left[ \left\{ \Xi (x,t) 
% NEW REWORKED      - \int_{0}^{t} d\tau \Pi(x,t,\tau) \right\} \left<\rho (\sigma(x_t,t), t) \right>\right] \\
% NEW REWORKED      \label{eqno_S10}
% NEW REWORKED      &+& \displaystyle \frac{\partial}{\partial \sigma } \left[\int_{0}^{t} d\tau 
% NEW REWORKED      \Omega(x,t,\tau)  \vphantom{\int_{0}^{t}}  
% NEW REWORKED      \displaystyle \frac{\partial \left<\rho (\sigma (x_t,t), t)\right>}{\partial \sigma} \right ] 
% NEW REWORKED      \end{eqnarray}
% NEW REWORKED      %
% NEW REWORKED      %\begin{eqnarray}
% NEW REWORKED      %\nonumber
% NEW REWORKED      %\lefteqn{\frac{\partial \left<\rho (\sigma(x_t,t), t)\right>}{\partial t}=}  \\
% NEW REWORKED      %\nonumber
% NEW REWORKED      %&-&\displaystyle \frac{\partial}{\partial \sigma}\left\{ \left[\left< \vphantom{\int_{0}^{t} d\tau} \eta(\sigma(x_t,t), D(x_t),
% NEW REWORKED      %q(x_t), r(x_t), \epsilon(x_t,t)) \right> \right. \right. \\
% NEW REWORKED      %\nonumber
% NEW REWORKED      %&-& \left. \left. \int_{0}^{t} d\tau Cov_0\left[\vphantom{\int_{0}^{t} d\tau} \eta(\sigma(x_t,t), D(x_t), q(x_t), r(x_t),
% NEW REWORKED      %\epsilon(x_t,t)); \right. \right. \right. \\
% NEW REWORKED      %\nonumber
% NEW REWORKED      %& & \left. \left. \left. \displaystyle \frac{\partial \eta (\sigma(x_{t-\tau}, t-\tau), D(x_{t-\tau}), q(x_{t-\tau}), r(x_{t-\tau}) \epsilon(x_{t-\tau}, t-\tau)}{\partial \sigma}\right]\right]\left<\rho (\sigma(x_t,t), t)\right>\right\} \\
% NEW REWORKED      %\nonumber
% NEW REWORKED      %&+& \displaystyle \frac{\partial}{\partial \sigma }\left\{ \left[\int_{0}^{t} d\tau Cov_0\left[\eta(\sigma (x_t,t), D(x_t), q(x_t),
% NEW REWORKED      %r (x_t), \epsilon(x_t,t)); \right. \right. \right. \\
% NEW REWORKED      %\nonumber
% NEW REWORKED      %& & \left. \left. \left. \eta (\sigma(x_{t-\tau}, t-\tau), D (x_{t-\tau}), q(x_{t-\tau}), r(x_{t-\tau}), \epsilon
% NEW REWORKED      %(x_{t-\tau},
% NEW REWORKED      %t-\tau)) \right] \vphantom{\int_{0}^{t}} \right]  \displaystyle \frac{\partial \left<\rho (\sigma (x_t,t), t)\right>}{\partial \sigma}
% NEW REWORKED      %\right \} \\
% NEW REWORKED      %\label{eqno_S10}
% NEW REWORKED      %\end{eqnarray}
% NEW REWORKED      %
% NEW REWORKED      to  exact  second  order  (to  the order of the covariance time of $ \eta $). By
% NEW REWORKED      combining  Eqs.~(\ref{eqno_S9})  and ~(\ref{eqno_S10}) and rearranging the terms
% NEW REWORKED      yields  an  Eulerian--Lagrangian  form  of Fokker-Planck equation \footnote{also
% NEW REWORKED      known   as  Forward--Kolmogorov  Equation  or  Fokker--Planck--Kolmogorov  (FPK)
% NEW REWORKED      equation               \citep{book:Risken},               \citep{book:Gardiner},
% NEW REWORKED      \citep{Schueller:1997}}    (FPE)    (for    details,   refer   to   Kavvas
% NEW REWORKED      \citep{Kavvas:2003}),   which,  for  writing  convenience,  after  denoting  the
% NEW REWORKED      covariance between the derivative of the function $\eta$ with respect to stress,
% NEW REWORKED      $\sigma$   at   time   $t$   and   the  function  $\eta$  at  time  $t-\tau$  by
% NEW REWORKED      $\Phi(x,t,\tau)$, as,
% NEW REWORKED      %
% NEW REWORKED      \begin{eqnarray}
% NEW REWORKED      \nonumber
% NEW REWORKED      \lefteqn{\Phi (x, t,\tau) = 
% NEW REWORKED      Cov_0 \left[ \displaystyle \frac{\partial \eta(\sigma(x_t,t), D(x_t), 
% NEW REWORKED      q(x_t), r(x_t) \epsilon(x_t,t))}{\partial \sigma} ; \right.} \\
% NEW REWORKED      \label{eqno_S12a}
% NEW REWORKED      && \left. \eta(\sigma(x_{t-\tau}, t-\tau), D(x_{t-\tau}), q(x_{t-\tau}), r(x_{t-\tau}),
% NEW REWORKED      \epsilon(x_{t-\tau}, t-\tau)) \vphantom{\int_{0}^{t} d\tau} \right],
% NEW REWORKED      \end{eqnarray}
% NEW REWORKED      
% NEW REWORKED      \noindent can be written as,
% NEW REWORKED      
% NEW REWORKED      \begin{eqnarray}
% NEW REWORKED      \nonumber
% NEW REWORKED      \lefteqn{\frac{\partial P(\sigma(x_t,t), t)}{\partial t} =} \\
% NEW REWORKED      \nonumber
% NEW REWORKED      & - & \displaystyle \frac{\partial}{\partial \sigma}\left[ \left\{ \Xi (x,t) 
% NEW REWORKED      + \int_{0}^{t} d\tau \Phi(x,t,\tau) \right\} P(\sigma(x_t,t), t) \right] \\
% NEW REWORKED      \label{eqno_S12}
% NEW REWORKED      &+& \displaystyle \frac{\partial^2}{\partial \sigma^2 } \left[\int_{0}^{t} d\tau 
% NEW REWORKED      \Omega(x,t,\tau)  \vphantom{\int_{0}^{t}}  
% NEW REWORKED      P (\sigma (x_t,t), t) \right ] 
% NEW REWORKED      \end{eqnarray}
% NEW REWORKED      %
% NEW REWORKED      %
% NEW REWORKED      %\begin{eqnarray}
% NEW REWORKED      %\nonumber
% NEW REWORKED      %\lefteqn{\displaystyle \frac{\partial P(\sigma(x_t,t), t)}{\partial t}=} \\
% NEW REWORKED      %\nonumber
% NEW REWORKED      %&-& \displaystyle \frac{\partial}{\partial \sigma} \left[ \left\{\left< \vphantom{\int_{0}^{t} d\tau} \eta(\sigma(x_t,t), D(x_t),
% NEW REWORKED      %q(x_t), r(x_t) \epsilon(x_t,t)) \right> \right. \right. \\
% NEW REWORKED      %\nonumber
% NEW REWORKED      %&+& \left. \left. \int_{0}^{t} d\tau Cov_0 \left[ \displaystyle \frac{\partial \eta(\sigma(x_t,t), D(x_t), q(x_t), r(x_t)
% NEW REWORKED      %\epsilon(x_t,t))}{\partial \sigma}; \right. \right. \right. \\
% NEW REWORKED      %\nonumber
% NEW REWORKED      %& & \left. \left. \left. \eta(\sigma(x_{t-\tau}, t-\tau), D(x_{t-\tau}), q(x_{t-\tau}), r(x_{t-\tau}),
% NEW REWORKED      %\epsilon(x_{t-\tau}, t-\tau) \vphantom{\int_{0}^{t} d\tau} \right] \right \} P(\sigma(x_t,t),t) \right] \\
% NEW REWORKED      %\nonumber
% NEW REWORKED      %&+& \displaystyle \frac{\partial^2}{\partial \sigma^2} \left[ \left\{ \int_{0}^{t} d\tau Cov_0 \left[ \vphantom{\int_{0}^{t}}
% NEW REWORKED      %\eta(\sigma(x_t,t), D(x_t,t), q(x_t,t), r (x_t,t), \epsilon(x_t,t)); \right. \right. \right. \\
% NEW REWORKED      %\nonumber
% NEW REWORKED      %& & \left. \left. \left. \eta_1 (\sigma(x_{t-\tau}, t-\tau), D(x_{t-\tau}), q(x_{t-\tau}), r(x_{t-\tau}),
% NEW REWORKED      %\epsilon(x_{t-\tau}, t-\tau)) \vphantom{\int_{0}^{t}} \right] \vphantom{\int_{0}^{t}} \right\} P(\sigma (x_t,t),t) \right] \\
% NEW REWORKED      %\label{eqno_S12}
% NEW REWORKED      %\end{eqnarray}
% NEW REWORKED      %
% NEW REWORKED      %
% NEW REWORKED      to   exact  second  order.  The  terms  $\Xi(x,t)$,
% NEW REWORKED      $\Phi(x,t,\tau)$,        and        $\Omega(x,t,\tau)$       in
% NEW REWORKED      Eq.~(\ref{eqno_S12})   are   given  by  Eqs.~(\ref{eqno_S10a}),
% NEW REWORKED      (\ref{eqno_S12a}),    and    (\ref{eqno_S10c})    respectively.
% NEW REWORKED      Eq.~(\ref{eqno_S12})   is   the   most   general  relation  for
% NEW REWORKED      probabilistic   behavior   of   inelastic   (non--linear)   1-D
% NEW REWORKED      constitutive  rate equation. 
% NEW REWORKED      %
% NEW REWORKED      %
% NEW REWORKED      %
% NEW REWORKED      %
%
%
%
%       %
%       %
%       %
%       %
%ORIGINAL     
%ORIGINAL      %       %
%ORIGINAL      %       %
%ORIGINAL      %       %
%ORIGINAL      %       %
%ORIGINAL      This       ensemble       average       was       recently       derived   as
%ORIGINAL      (for detailed derivation, refer to \citet{Kavvas:1996})
%ORIGINAL      %
%ORIGINAL      \begin{eqnarray}
%ORIGINAL      \nonumber
%ORIGINAL      \lefteqn{\frac{\partial \left<\rho (\sigma(x_t,t), t)\right>}{\partial t}=}  \\
%ORIGINAL      \nonumber
%ORIGINAL      &-&\displaystyle \frac{\partial}{\partial \sigma}\left\{ \left[\left< \vphantom{\int_{0}^{t} d\tau} \eta(\sigma(x_t,t), D(x_t),
%ORIGINAL      q(x_t), r(x_t), \epsilon(x_t,t)) \right> \right. \right. \\
%ORIGINAL      \nonumber
%ORIGINAL      &-& \left. \left. \int_{0}^{t} d\tau Cov_0\left[\vphantom{\int_{0}^{t} d\tau} \eta(\sigma(x_t,t), D(x_t), q(x_t), r(x_t),
%ORIGINAL      \epsilon(x_t,t)); \right. \right. \right. \\
%ORIGINAL      \nonumber
%ORIGINAL      & & \left. \left. \left. \displaystyle \frac{\partial \eta (\sigma(x_{t-\tau}, t-\tau), D(x_{t-\tau}), q(x_{t-\tau}), r(x_{t-\tau}) \epsilon(x_{t-\tau}, t-\tau)}{\partial \sigma}\right]\right]\left<\rho (\sigma(x_t,t), t)\right>\right\} \\
%ORIGINAL      \nonumber
%ORIGINAL      &+& \displaystyle \frac{\partial}{\partial \sigma }\left\{ \left[\int_{0}^{t} d\tau Cov_0\left[\eta(\sigma (x_t,t), D(x_t), q(x_t),
%ORIGINAL      r (x_t), \epsilon(x_t,t)); \right. \right. \right. \\
%ORIGINAL      \nonumber
%ORIGINAL      & & \left. \left. \left. \eta (\sigma(x_{t-\tau}, t-\tau), D (x_{t-\tau}), q(x_{t-\tau}), r(x_{t-\tau}), \epsilon
%ORIGINAL      (x_{t-\tau},
%ORIGINAL      t-\tau)) \right] \vphantom{\int_{0}^{t}} \right]  \displaystyle \frac{\partial \left<\rho (\sigma (x_t,t), t)\right>}{\partial \sigma}
%ORIGINAL      \right \} \\
%ORIGINAL      \label{eqno_S10}
%ORIGINAL      \end{eqnarray}
%ORIGINAL      %
%ORIGINAL      to  exact  second  order  (to  the order of the covariance time of $ \eta $). 
%ORIGINAL      %
%ORIGINAL      In  Eq.~(\ref{eqno_S10}), $Cov_0[\cdot]$ is the time ordered covariance function
%ORIGINAL      defined by
%ORIGINAL      %
%ORIGINAL      \begin{equation}
%ORIGINAL      \label{eqno_S11}
%ORIGINAL      Cov_0\left[\eta(x,t_1), \eta(x,t_2)\right]
%ORIGINAL      =
%ORIGINAL      \left<\eta(x,t_1) \eta(x,t_2)\right>
%ORIGINAL      -
%ORIGINAL      \left<\eta(x,t_1)\right> \cdot \left<\eta(x,t_2)\right>
%ORIGINAL      \end{equation}
%ORIGINAL      
%ORIGINAL      By  combining  Eqs.~(\ref{eqno_S10})  and ~(\ref{eqno_S9}) and rearranging the terms
%ORIGINAL      yields    the   following   Fokker-Planck   equation   (FPE,   also   known   as
%ORIGINAL      Forward--Kolmogorov Equation or Fokker--Planck--Kolmogorov FPK Equation) 
%ORIGINAL      (\citet{book:Risken}, \citet{book:Gardiner}, \citet{Schueller:1997report}):
%ORIGINAL      %
%ORIGINAL      \begin{eqnarray}
%ORIGINAL      \nonumber
%ORIGINAL      \lefteqn{\displaystyle \frac{\partial P(\sigma(x_t,t), t)}{\partial t}=} 
%ORIGINAL      \\ \nonumber
%ORIGINAL      &-& 
%ORIGINAL      \displaystyle \frac{\partial}{\partial \sigma} 
%ORIGINAL      \left[ \left\{\left< \vphantom{\int_{0}^{t} d\tau} \eta(\sigma(x_t,t), D(x_t),
%ORIGINAL      q(x_t), r(x_t) \epsilon(x_t,t)) \right> \right. \right. 
%ORIGINAL      \\ \nonumber
%ORIGINAL      &+& 
%ORIGINAL      \left. \left. \int_{0}^{t} d\tau Cov_0 
%ORIGINAL      \left[ \displaystyle \frac{\partial \eta(\sigma(x_t,t), D(x_t), q(x_t), r(x_t)
%ORIGINAL      \epsilon(x_t,t))}{\partial \sigma}; \right. \right. \right. 
%ORIGINAL      \\  \nonumber
%ORIGINAL      & & 
%ORIGINAL      \left. \left. \left. \eta(\sigma(x_{t-\tau}, t-\tau), 
%ORIGINAL      D(x_{t-\tau}), q(x_{t-\tau}), r(x_{t-\tau}),
%ORIGINAL      \epsilon(x_{t-\tau}, t-\tau) \vphantom{\int_{0}^{t} d\tau} \right] \right \} 
%ORIGINAL      P(\sigma(x_t,t),t) \right] 
%ORIGINAL      \\ \nonumber
%ORIGINAL      &+& 
%ORIGINAL      \displaystyle \frac{\partial^2}{\partial \sigma^2} 
%ORIGINAL      \left[ \left\{ \int_{0}^{t} d\tau Cov_0 \left[ \vphantom{\int_{0}^{t}}
%ORIGINAL      \eta(\sigma(x_t,t), D(x_t,t), q(x_t,t), r (x_t,t), 
%ORIGINAL      \epsilon(x_t,t)); \right. \right. \right. 
%ORIGINAL      \\ \nonumber
%ORIGINAL      & & 
%ORIGINAL      \left. \left. \left. \eta_1 (\sigma(x_{t-\tau}, t-\tau), 
%ORIGINAL      D(x_{t-\tau}), q(x_{t-\tau}), r(x_{t-\tau}),
%ORIGINAL      \epsilon(x_{t-\tau}, t-\tau)) \vphantom{\int_{0}^{t}} \right] 
%ORIGINAL      \vphantom{\int_{0}^{t}} \right\} P(\sigma (x_t,t),t) \right] 
%ORIGINAL      \\
%ORIGINAL      \label{eqno_S12}
%ORIGINAL      \end{eqnarray}
%ORIGINAL      %
%ORIGINAL      to  exact  second  order.  
%ORIGINAL      %
%ORIGINAL      This  is  the most general relation for probabilistic
%ORIGINAL      behavior  of  inelastic (non--linear, elastic--plastic) 
%ORIGINAL      1-D stochastic incremental constitutive equation. 
%ORIGINAL      %
%ORIGINAL      %
%ORIGINAL      %
%ORIGINAL      %
%ORIGINAL      %
%ORIGINAL      %
%ORIGINAL      %
%ORIGINAL      
%ORIGINAL      
%ORIGINAL      

%       %
%  LATEST REWORKED     
%  LATEST REWORKED 
%  LATEST REWORKED 
%  LATEST REWORKED 
%       %
This       ensemble       average       was       recently       derived   as
(for detailed derivation, refer to \citet{Kavvas:1996})
%

\begin{eqnarray}                                                                              
\nonumber                                                                                     
\lefteqn{\frac{\partial \left<\rho (\sigma(x_t,t), t)\right>}{\partial t} =} \\               
\nonumber                                                                                     
& - & \displaystyle \frac{\partial}{\partial \sigma}\left[ \left\{ \Xi (x,t)                  
- \int_{0}^{t} d\tau \Pi(x,t,\tau) \right\} \left<\rho (\sigma(x_t,t), t) \right>\right] \\   
\label{eqno_S10}                                                                              
&+& \displaystyle \frac{\partial}{\partial \sigma } \left[\int_{0}^{t} d\tau                  
\Omega(x,t,\tau)  \vphantom{\int_{0}^{t}}                                                     
\displaystyle \frac{\partial \left<\rho (\sigma (x_t,t), t)\right>}{\partial \sigma} \right ] 
\end{eqnarray}                                                                                
%
to exact  second  order  (to  the order of the covariance time of $ \eta $). 
In Eq.~(\ref{eqno_S10}), $\Xi (x,t)$ is the ensemble average of
the function $\eta$ at time $t$ and is given as,
%
\begin{equation}
\nonumber
%\label{eqno_S10a}
\Xi (x,t) = \left< \vphantom{\int_{0}^{t} d\tau} \eta(\sigma(x_t,t), D(x_t),
q(x_t), r(x_t), \epsilon(x_t,t)) \right>;
\end{equation}
%
and $\Pi(x,t,\tau)$ is the covariance between  the  function  $\eta$  at  time  $t$ 
and the derivative of the function $\eta$ with respect to stress, $\sigma$, at time $t-\tau$ and is
given as,
%
\begin{eqnarray}
\nonumber
\lefteqn{\Pi (x,t,\tau) 
= 
Cov_0\left[\vphantom{\int_{0}^{t} d\tau} \eta(\sigma(x_t,t), D(x_t), q(x_t), r(x_t),
\epsilon(x_t,t)); \right.} \\
\nonumber
%\label{eqno_S10b}
&& \left. \displaystyle \frac{\partial \eta (\sigma(x_{t-\tau}, t-\tau), D(x_{t-\tau}), q(x_{t-\tau}), 
r(x_{t-\tau}), \epsilon(x_{t-\tau}, t-\tau))}{\partial \sigma} \right];
\end{eqnarray}
%
while $\Omega(x,t,\tau)$ is the covariance between function $\eta$ at time $t$ and function
$\eta$ at time $t-\tau$ and is given as,
%
\begin{eqnarray}
\nonumber
\lefteqn{\Omega (x, t,\tau) = Cov_0\left[\vphantom{\int_{0}^{t} d\tau} \eta(\sigma(x_t,t), D(x_t), q(x_t), r(x_t),
\epsilon(x_t,t)); \right.} \\
\nonumber
%\label{eqno_S10c}
&& \left. \eta (\sigma(x_{t-\tau}, t-\tau), D (x_{t-\tau}), q(x_{t-\tau}), r(x_{t-\tau}), \epsilon
(x_{t-\tau}, t-\tau)) \vphantom{\int_{0}^{t} d\tau} \right]
\end{eqnarray}
%
%
%
%\begin{eqnarray}
%\nonumber
%\lefteqn{\frac{\partial \left<\rho (\sigma(x_t,t), t)\right>}{\partial t}=}  \\
%\nonumber
%&-&\displaystyle \frac{\partial}{\partial \sigma}\left\{ \left[\left< \vphantom{\int_{0}^{t} d\tau} \eta(\sigma(x_t,t), D(x_t),
%q(x_t), r(x_t), \epsilon(x_t,t)) \right> \right. \right. \\
%\nonumber
%&-& \left. \left. \int_{0}^{t} d\tau Cov_0\left[\vphantom{\int_{0}^{t} d\tau} \eta(\sigma(x_t,t), D(x_t), q(x_t), r(x_t),
%\epsilon(x_t,t)); \right. \right. \right. \\
%\nonumber
%& & \left. \left. \left. \displaystyle \frac{\partial \eta (\sigma(x_{t-\tau}, t-\tau), D(x_{t-\tau}), q(x_{t-\tau}), r(x_{t-\tau}) \epsilon(x_{t-\tau}, t-\tau)}{\partial \sigma}\right]\right]\left<\rho (\sigma(x_t,t), t)\right>\right\} \\
%\nonumber
%&+& \displaystyle \frac{\partial}{\partial \sigma }\left\{ \left[\int_{0}^{t} d\tau Cov_0\left[\eta(\sigma (x_t,t), D(x_t), q(x_t),
%r (x_t), \epsilon(x_t,t)); \right. \right. \right. \\
%\nonumber
%& & \left. \left. \left. \eta (\sigma(x_{t-\tau}, t-\tau), D (x_{t-\tau}), q(x_{t-\tau}), r(x_{t-\tau}), \epsilon
%(x_{t-\tau},
%t-\tau)) \right] \vphantom{\int_{0}^{t}} \right]  \displaystyle \frac{\partial \left<\rho (\sigma (x_t,t), t)\right>}{\partial \sigma}
%\right \} \\
%\label{eqno_S10}
%\end{eqnarray}
%
%to  exact  second  order  (to  the order of the covariance time of $ \eta $). 
%%
It is important to note that the covariances appearing in the above equations are time ordered 
covariances, which for two dummy variables $p(x,t_1)$
and    $p(x,t_2)$   can be written as,  
%
\begin{equation}
\label{eqno_S11}  
Cov_0\left[p(x,t_1); p(x,t_2)\right] = \left<p(x,t_1)    p(x,t_2)\right>-\left<p(x,t_1)\right>    \cdot
\left<p(x,t_2)\right> 
\end{equation}
%
%In  Eq.~(\ref{eqno_S10}), $Cov_0[\cdot]$ is the time ordered covariance function
%defined by
%%
%\begin{equation}
%\label{eqno_S11}
%Cov_0\left[\eta(x,t_1), \eta(x,t_2)\right]
%=
%\left<\eta(x,t_1) \eta(x,t_2)\right>
%-
%\left<\eta(x,t_1)\right> \cdot \left<\eta(x,t_2)\right>
%\end{equation}

By  combining  Eqs.~(\ref{eqno_S9})  and ~(\ref{eqno_S10}) and rearranging the terms
yields    the   following   
Eulerian--Lagrangian  form  of second-order accurate Fokker-Planck equation \footnote{also
known   as  Forward--Kolmogorov  Equation  or  Fokker--Planck--Kolmogorov  (FPK)
equation               \citep{book:Risken},               \citep{book:Gardiner},
\citep{Schueller:1997}}    (FPE)    (for    details,   refer   to   Kavvas
\citep{Kavvas:2003}):
%
\begin{eqnarray}
\nonumber
\lefteqn{\frac{\partial P(\sigma(x_t,t), t)}{\partial t} =} \\
\nonumber
& - & \displaystyle \frac{\partial}{\partial \sigma}\left[ \left\{ \Xi (x,t) 
+ \int_{0}^{t} d\tau \Phi(x,t,\tau) \right\} P(\sigma(x_t,t), t) \right] \\
\label{eqno_S12}
&+& \displaystyle \frac{\partial^2}{\partial \sigma^2 } \left[\int_{0}^{t} d\tau 
\Omega(x,t,\tau)  \vphantom{\int_{0}^{t}}  
P (\sigma (x_t,t), t) \right ] 
\end{eqnarray}
%
where, $\Phi(x,t,\tau)$ is the time-ordered covariance between the derivative of the function 
$\eta$ with respect to stress, $\sigma$   at   time   $t$   and   the  function  $\eta$  
at  time  $t-\tau$  and is given as,
%
\begin{eqnarray}
\nonumber
\lefteqn{\Phi (x, t,\tau) = 
Cov_0 \left[ \displaystyle \frac{\partial \eta(\sigma(x_t,t), D(x_t), 
q(x_t), r(x_t) \epsilon(x_t,t))}{\partial \sigma} ; \right.} \\
\nonumber
%\label{eqno_S12a}
&& \left. \eta(\sigma(x_{t-\tau}, t-\tau), D(x_{t-\tau}), q(x_{t-\tau}), r(x_{t-\tau}),
\epsilon(x_{t-\tau}, t-\tau)) \vphantom{\int_{0}^{t} d\tau} \right]
\end{eqnarray}

%\begin{eqnarray}
%\nonumber
%\lefteqn{\displaystyle \frac{\partial P(\sigma(x_t,t), t)}{\partial t}=} 
%\\ \nonumber
%&-& 
%\displaystyle \frac{\partial}{\partial \sigma} 
%\left[ \left\{\left< \vphantom{\int_{0}^{t} d\tau} \eta(\sigma(x_t,t), D(x_t),
%q(x_t), r(x_t) \epsilon(x_t,t)) \right> \right. \right. 
%\\ \nonumber
%&+& 
%\left. \left. \int_{0}^{t} d\tau Cov_0 
%\left[ \displaystyle \frac{\partial \eta(\sigma(x_t,t), D(x_t), q(x_t), r(x_t)
%\epsilon(x_t,t))}{\partial \sigma}; \right. \right. \right. 
%\\  \nonumber
%& & 
%\left. \left. \left. \eta(\sigma(x_{t-\tau}, t-\tau), 
%D(x_{t-\tau}), q(x_{t-\tau}), r(x_{t-\tau}),
%\epsilon(x_{t-\tau}, t-\tau) \vphantom{\int_{0}^{t} d\tau} \right] \right \} 
%P(\sigma(x_t,t),t) \right] 
%\\ \nonumber
%&+& 
%\displaystyle \frac{\partial^2}{\partial \sigma^2} 
%\left[ \left\{ \int_{0}^{t} d\tau Cov_0 \left[ \vphantom{\int_{0}^{t}}
%\eta(\sigma(x_t,t), D(x_t,t), q(x_t,t), r (x_t,t), 
%\epsilon(x_t,t)); \right. \right. \right. 
%\\ \nonumber
%& & 
%\left. \left. \left. \eta_1 (\sigma(x_{t-\tau}, t-\tau), 
%D(x_{t-\tau}), q(x_{t-\tau}), r(x_{t-\tau}),
%\epsilon(x_{t-\tau}, t-\tau)) \vphantom{\int_{0}^{t}} \right] 
%\vphantom{\int_{0}^{t}} \right\} P(\sigma (x_t,t),t) \right] 
%\\
%\label{eqno_S12}
%\end{eqnarray}
%
%to  exact  second  order.  
%

%Eq.~(\ref{eqno_S12}) is  the most general relation for probabilistic
%behavior  of  inelastic (non--linear, elastic--plastic) 
%1-D stochastic incremental constitutive equation. 
%
%
%
%
%
%
%









%   
%   This       ensemble       average       was       recently       derived   as
%   (for detailed derivation, refer to \citet{Kavvas:1996})
%   %
%   \begin{eqnarray}
%   \nonumber
%   \lefteqn{\frac{\partial \left<\rho (\sigma(x_t,t), t)\right>}{\partial t}=}  \\
%   \nonumber
%   &-&\displaystyle \frac{\partial}{\partial \sigma}\left\{ \left[\left< \vphantom{\int_{0}^{t} d\tau} \eta(\sigma(x_t,t), D(x_t),
%   q(x_t), r(x_t), \epsilon(x_t,t)) \right> \right. \right. \\
%   \nonumber
%   &-& \left. \left. \int_{0}^{t} d\tau Cov_0\left[\vphantom{\int_{0}^{t} d\tau} \eta(\sigma(x_t,t), D(x_t), q(x_t), r(x_t),
%   \epsilon(x_t,t)); \right. \right. \right. \\
%   \nonumber
%   & & \left. \left. \left. \displaystyle \frac{\partial \eta (\sigma(x_{t-\tau}, t-\tau), D(x_{t-\tau}), q(x_{t-\tau}), r(x_{t-\tau}) \epsilon(x_{t-\tau}, t-\tau)}{\partial \sigma}\right]\right]\left<\rho (\sigma(x_t,t), t)\right>\right\} \\
%   \nonumber
%   &+& \displaystyle \frac{\partial}{\partial \sigma }\left\{ \left[\int_{0}^{t} d\tau Cov_0\left[\eta(\sigma (x_t,t), D(x_t), q(x_t),
%   r (x_t), \epsilon(x_t,t)); \right. \right. \right. \\
%   \nonumber
%   & & \left. \left. \left. \eta (\sigma(x_{t-\tau}, t-\tau), D (x_{t-\tau}), q(x_{t-\tau}), r(x_{t-\tau}), \epsilon
%   (x_{t-\tau},
%   t-\tau)) \right] \vphantom{\int_{0}^{t}} \right]  \displaystyle \frac{\partial \left<\rho (\sigma (x_t,t), t)\right>}{\partial \sigma}
%   \right \} \\
%   \label{eqno_S10}
%   \end{eqnarray}
%   %
%   to  exact  second  order  (to  the order of the covariance time of $ \eta $). 
%   %
%   In  Eq.~(\ref{eqno_S10}), $Cov_0[\cdot]$ is the time ordered covariance function
%   defined by
%   %
%   \begin{equation}
%   \label{eqno_S11}
%   Cov_0\left[\eta(x,t_1), \eta(x,t_2)\right]
%   =
%   \left<\eta(x,t_1) \eta(x,t_2)\right>
%   -
%   \left<\eta(x,t_1)\right> \cdot \left<\eta(x,t_2)\right>
%   \end{equation}
%   
%   By  combining  Eqs.~(\ref{eqno_S10})  and ~(\ref{eqno_S9}) and rearranging the terms
%   yields    the   following   Fokker-Planck   equation   (FPE,   also   known   as
%   Forward--Kolmogorov Equation or Fokker--Planck--Kolmogorov FPK Equation) 
%   (\citet{book:Risken}, \citet{book:Gardiner}, \citet{Schueller:1997report}):
%   %
%   \begin{eqnarray}
%   \nonumber
%   \lefteqn{\displaystyle \frac{\partial P(\sigma(x_t,t), t)}{\partial t}=} 
%   \\ \nonumber
%   &-& 
%   \displaystyle \frac{\partial}{\partial \sigma} 
%   \left[ \left\{\left< \vphantom{\int_{0}^{t} d\tau} \eta(\sigma(x_t,t), D(x_t),
%   q(x_t), r(x_t) \epsilon(x_t,t)) \right> \right. \right. 
%   \\ \nonumber
%   &+& 
%   \left. \left. \int_{0}^{t} d\tau Cov_0 
%   \left[ \displaystyle \frac{\partial \eta(\sigma(x_t,t), D(x_t), q(x_t), r(x_t)
%   \epsilon(x_t,t))}{\partial \sigma}; \right. \right. \right. 
%   \\  \nonumber
%   & & 
%   \left. \left. \left. \eta(\sigma(x_{t-\tau}, t-\tau), 
%   D(x_{t-\tau}), q(x_{t-\tau}), r(x_{t-\tau}),
%   \epsilon(x_{t-\tau}, t-\tau) \vphantom{\int_{0}^{t} d\tau} \right] \right \} 
%   P(\sigma(x_t,t),t) \right] 
%   \\ \nonumber
%   &+& 
%   \displaystyle \frac{\partial^2}{\partial \sigma^2} 
%   \left[ \left\{ \int_{0}^{t} d\tau Cov_0 \left[ \vphantom{\int_{0}^{t}}
%   \eta(\sigma(x_t,t), D(x_t,t), q(x_t,t), r (x_t,t), 
%   \epsilon(x_t,t)); \right. \right. \right. 
%   \\ \nonumber
%   & & 
%   \left. \left. \left. \eta_1 (\sigma(x_{t-\tau}, t-\tau), 
%   D(x_{t-\tau}), q(x_{t-\tau}), r(x_{t-\tau}),
%   \epsilon(x_{t-\tau}, t-\tau)) \vphantom{\int_{0}^{t}} \right] 
%   \vphantom{\int_{0}^{t}} \right\} P(\sigma (x_t,t),t) \right] 
%   \\
%   \label{eqno_S12}
%   \end{eqnarray}
%   %
%   to  exact  second  order.  
%   %
%   %
%   %
%   %
%   %
%   %
%   %
%   %
%   This  is  the most general relation for probabilistic
%   behavior  of  inelastic (non--linear, elastic--plastic) 
%   1-D stochastic incremental constitutive equation. 
%   %
%   %
%   
%   
%   
%   
%   
%   
%   
%   
%  
%   
%   
%   
%   
%   
%   
%   
%   
%   
%   
%   
%   
%   
%   
%   
Eq.~(\ref{eqno_S12}) is  the most general relation for probabilistic
behavior  of  inelastic (non--linear, elastic--plastic) 
1-D stochastic incremental constitutive equation.
%
The solution
of    this  deterministic  linear  FPE  (Eq.~(\ref{eqno_S12})),  in  terms  of  the
probability density $P(\sigma,t)$,
under appropriate initial and boundary conditions will yield the PDF of the
state variable $\sigma$ of the original 1-D non-linear stochastic constitutive
rate  equation  (Eq.~(\ref{eqno_S3})).  
%
It  is  important  to  note that while the
original     equation     (Eq.~(\ref{eqno_S3}))    is    non-linear,    the    FPE
(Eq.~(\ref{eqno_S12}))  is linear in terms of its unknown, the probability density
$P(\sigma,t)$  of  the  state  variable $ \sigma $. 
%
This linearity, in turn,
provides significant advantages in the solution of the probabilistic behavior of
the incremental constitutive  equation (Eq.~(\ref{eqno_S3})).


One  should  also  note  that  Eq.~(\ref{eqno_S12}) is a mixed Eulerian-Lagrangian
equation.  This  stems from the fact that while the real space location $x_t$ at
time $t$ is known, the
location  $x_{t-\tau}$  is  an  unknown. {If one assumes small strain theory,
one can relate the unknown location $ x_{t-\tau} $ from the known location 
$ x_t $ by using the strain rate, $\dot \epsilon $ (=$d \epsilon/dt$) as}, 
%
%  \begin{equation}
%  \label{eqno_S13}
%  {d\epsilon=\dot \epsilon \tau = \frac{x_t-x_{t-\tau}}{x_t}}
%  \end{equation}
%  %
%  {Or, rearranging}
%
\begin{equation}
\label{eqno_S14}
{x_{t-\tau}=(1-\dot \epsilon \tau)x_t}
\end{equation}
%
Once  the  probability density function $P(\sigma,t)$ is obtained it can be used
to obtain the mean of state variable ($\sigma$) by usual expectation operation
%
\begin{equation}
\label{eqno_S16}
<\sigma(t)>=\int \sigma(t) P(\sigma(t)) d\sigma(t)
\end{equation}


Another  possible  way  to  obtain  the  mean  of  state  variable is to use the
equivalence   between   FPE   and   It{\^o}--operator   stochastic  differential  equation
\citep{book:Gardiner}.    In  this case  It{\^o}--operator stochastic differential equation
equivalent to Eq.~(\ref{eqno_S12}) is
%
\begin{eqnarray}
\nonumber
d\sigma (x,t)
&=&
\left\{ \left< \vphantom{\displaystyle 
\frac{\partial \eta_1 (\sigma(x_t,t), A(x_t,t), 
\epsilon(x_t,t))}{\partial \sigma}} 
\eta (\sigma(x_t,t), D(x_t), q(x_t), r(x_t), \epsilon(x_t,t)) \right> \right. 
\\ \nonumber
&+& 
\left. \int_0^t d \tau Cov_0 
\left[\displaystyle \frac{\partial 
\eta (\sigma(x_t,t), D(x_t), q(x_t), r(x_t), \epsilon(x_t,t))}
{\partial \sigma}; \right. \right. 
\\ \nonumber
& & 
\left. \left. \eta (\sigma (x_{t-\tau},t-\tau), 
D(x_{t-\tau}), q(x_{t-\tau}), r(x_{t-\tau}), 
\epsilon (x_{t-\tau},
t-\tau)) \vphantom{\displaystyle 
\frac{\partial \eta_1 (\sigma(x_t,t), A(x_t,t), 
\epsilon(x_t,t))}{\partial \sigma}}\right] \right\} dt 
\\ 
&+& 
b(\sigma,t)dW(t)
\label{eqno_S17} 
\end{eqnarray}
%
where,
%
\begin{eqnarray}
\nonumber
b^2(\sigma,t) 
&=& 
2 \int_0^t d \tau Cov_0 
\left[\vphantom{\int_0^t d \tau} 
\eta(\sigma(x_t,t), D(x_t), q(x_t), r(x_t),
\epsilon(x_t,t)); \right. 
\\
& & 
\left. \eta(\sigma (x_{t-\tau},t-\tau), 
D (x_{t-\tau}), q(x_{t-\tau}), r(x_{t-\tau}), \epsilon
(x_{t-\tau},t-\tau))
\vphantom{\int_0^t d \tau} \right]
\nonumber \\
\label{eqno_S18}
\end{eqnarray}
%
and,  $dW(t)$  is  an  increment of Wiener process W with $ <dW(t)> = 0 $. It is
also  interesting  to note that all the stochasticity of the original equation
(Eq.~(\ref{eqno_S3})) are lumped together in the last term (Wiener increment term)
of  the  right-hand-side  of  Eq.~(\ref{eqno_S17}).  
%
By taking  advantage of the
independent  increment  property  of the Wiener process ($<dW(t)>=0$), one can
derive  the differential equation which describes the evolution of mean of state
variable  ($\sigma$) of the nonlinear constitutive rate equation in time and
space as, (e.g. \citep{Kavvas:2003})
%
\begin{eqnarray}
\nonumber
\displaystyle \frac{<d\sigma (x,t)>}{dt}
&=&
\left< \vphantom{\int_0^t d \tau} \eta(\sigma(x_t,t), D(x_t), q(x_t), r(x_t),
\epsilon(x_t,t)) \right> 
\\ \nonumber
&+& 
\int_0^t d \tau Cov_0 \left[\displaystyle 
\frac{\partial \eta (\sigma (x_t,t), 
D (x_t), q(x_t), r(x_t), \epsilon (x_t,t))}
{\partial \sigma}; \right. 
\\
& & 
\left. \eta (\sigma (x_{t-\tau},t-\tau), 
D(x_{t-\tau}), q(x_{t-\tau}), r(x_{t-\tau}), \epsilon
(x_{t-\tau},t-\tau)) \vphantom{\int_0^t d 
\tau Cov_0 \left[\displaystyle \frac{\partial \eta (\sigma (x_t,t), A (x_t,t),
\epsilon (x_t,t))}{\partial \sigma}; \right.} \right]
\nonumber \\
\label{eqno_S19}
\end{eqnarray}

Eq.~(\ref{eqno_S19})   is   a   nonlocal   integro-differential  equation  in
Eulerian-Lagrangian  form,  since, although the location $x_t$ at time $t$ is
known,  the  Lagrangian location  $x_{t-\tau}$ is an unknown which is
determined   by   Eq.~(\ref{eqno_S14}).
%\sout{or   it's   first   order   approximation
%(Eq.~\ref{eqno_S15})}.  
It is important to note that the state variable appearing
within  $\eta(\cdot)$ on the right-hand-side of Eq.~(\ref{eqno_S19}) is random
and needs to be treated accordingly.


This  concludes  the  development of  relation  for  probabilistic  behavior  of
1-D  elastic--plastic    constitutive  incremental   equation  with  stochastic 
coefficients  and stochastic  forcing  in  most general form.
%
In the following section the developed general
relation is specialized to two particular types of {point-location scale}
constitutive modeling: a) 1-D (shear) linear elastic constitutive behavior, and b)
1-D (shear) elastic-plastic  Drucker-Prager  associative  linear hardening 
constitutive behavior.


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

\section{Linear Elastic Probabilistic 1-D Constitutive Incremental Equation}
\label{LinearElastic}

The  1-D (shear) {point-location scale} linear elastic  constitutive rate equation
can be written as
%
\begin{equation}
\label{eqno_LE1}
\frac{d \sigma_{12}}{dt}=G \;\; \frac{d \epsilon_{12}}{dt}
\end{equation}
%
so  that   the function $ \eta $,  defined in Eq.~(\ref{eqno_S4}), can be written 
as
%
\begin{equation}
\label{eqno_LE2}
\eta=G \;\; \frac{d \epsilon_{12}}{dt}
\end{equation}
%
and    hence,  considering  both  the  shear modulus, 
%\sout{G (x)} 
{$G$} and the strain rate,
%\sout{$d\epsilon_{12}/dt(x_t,t)$} 
{$d\epsilon_{12}/dt(t)$} as random, one can substitute Eq.~(\ref{eqno_LE2}) in
Eq.~(\ref{eqno_S12})    to  obtain  the  particular  FPK  equation for
the  probabilistic  behavior  of 1-D (shear) {point-location scale} linear elastic
constitutive incremental equation as,
%
\begin{eqnarray}
\nonumber
{
\lefteqn{\displaystyle \frac{\partial P(\sigma_{12}(t),t)}{\partial t}=}
}
\\
&-& {
\displaystyle \frac{\partial}{\partial \sigma_{12}} \left[\left\{\left<G 
\displaystyle \frac{d\epsilon_{12}}{dt} (t) \right> \right. \right.
}
\nonumber                                                                      
\\
&+& {
\left. \left. \int_{0}^{t} d\tau Cov_0 \left[\displaystyle \frac{\partial}{\partial \sigma_{12}}
\left(G  \displaystyle \frac{d \epsilon_{12}}{dt} (t) \right);
G \displaystyle \frac{d\epsilon_{12}}{dt}(t-\tau) \right] \right \} 
P(\sigma_{12}(t),t) \right]
} 
\nonumber
\\
&+& {
\frac{\partial^2}{\partial \sigma_{12}^2}\left [\left\{\int_{0}^{t} d \tau Cov_0 
\left[G \displaystyle \frac{d\epsilon_{12}}{dt} (t); 
G \displaystyle \frac{d\epsilon_{12}}{dt} (t-\tau)\right] \right\}
P(\sigma_{12}(t),t) \right]
}
\nonumber \\
\label{eqno_LE3}
\end{eqnarray} 

The  first random 
%\sout{field} 
{process} in the covariance term of the first coefficient on
the r.h.s of above equation (Eq.~(\ref{eqno_LE3})) is independent of $\sigma_{12}$
and by noting  that  covariance  of  zero  with  any random 
%\sout{field} 
{process} is zero, one can
further simplify the above equation as,
%
\begin{eqnarray}
{
\lefteqn{\displaystyle \frac{\partial P(\sigma_{12}(t), t)}{\partial t}=}
} 
\nonumber
\\
&-& {
\displaystyle \frac{\partial}{\partial \sigma_{12}} \left[\left<G \displaystyle 
\frac{d\epsilon_{12}}{dt} (t) \right> P(\sigma_{12}(t),t) \right]
} 
\nonumber
\\
&+& {
\frac{\partial^2}{\partial \sigma_{12}^2}\left [\left\{\int_{0}^{t} d \tau Cov_0 
\left[G \displaystyle \frac{d\epsilon_{12}}{dt} (t); 
G \displaystyle \frac{d\epsilon_{12}}{dt} (t-\tau)\right] \right\} P(\sigma_{12}(t),t) \right]
}
\nonumber \\
\label{eqno_LE4}
\end{eqnarray}

Given  the  
%\sout{spatial  random  field  of} 
{random} shear modulus and 
%\sout{temporal-spatial random field} 
shear  strain  rate {random process},  with  appropriate  initial and boundary
conditions,  Eq.~(\ref{eqno_LE4}) will predict the evolution of probability density
function   of  shear  stress  with  time  following  1-D  linear  elastic  shear
constitutive  incremental  equation.  
%
Once the PDF of shear stress is obtained, one can
integrate   it   to  obtain  mean  (ensemble  average)  and  covariance  of  the
temporal-spatial  random field of shear stress. 
%
It should be noted that there exist several other exact methods (e.g. cumulant
expansion  method)  to  obtain  probabilistic  behavior  of stochastic processes
driven  by  linear  ordinary  differential  equations. 
% 
The  main  objective  of
FPK   based  approach  presented  here  is  to  deal  with  non-linear
stochastic  processes  (e.g. elastic-plastic constitutive rate equation) and the
linear elastic case is obtained as a special case of that non-linear process.





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
    
\section{Elastic--Plastic Probabilistic 1-D Constitutive Incremental Equation}
%\section{Probabilistic Behavior of Drucker-Prager Associative Linear Hardening Type 1-D Elastic-Plastic Shear Constitutive Rate Equation}
\label{ElasticPlastic}




For  materials  obeying  Drucker-Prager  yield  criteria (without cohesion), the
yield surface can be written as:
%
\begin{equation}
f = \sqrt{J_2} - \alpha I_1
\label{eqno_A1}
\end{equation}
%
where $J_2 = \frac{1}{2}S_{ij}s_{ij}$ is  the second invariant of the deviatoric
stress  tensor  $s_{ij} = \sigma_{ij} - 1/3 \delta_{ij} \sigma_{kk}$, and $I_1 =
\sigma_{ii}$  is  the  first  invariant  of  the stress tensor, and $\alpha$, an
internal   variable,   is   a   function   of   friction   angle  ($\alpha  =  2
\sin(\phi)/(\sqrt(3)(3-sin  \phi))$,  where  $\phi$  is the friction angle (e.g.
\citep{book:Chen}).


By  assuming  associative  flow  rule,  so  that the yield function has the same
derivatives as the
plastic flow function
%
\begin{equation}
\frac{\partial f}{\partial \sigma_{ij}}=\frac{\partial U}{\partial \sigma_{ij}}
\label{eqno_A2}
\end{equation}
%
one  can  expand  parts  of  the  tangent  constitutive tensor given in Eq.
~(\ref{el_or_elpl}) (from  Page  \pageref{el_or_elpl}), to read\footnote{A more 
detailed derivation of this probabilistic differentiation is given in
the Appendix.}
%
\begin{eqnarray}
A_{kl}
=
\frac{\partial f}{\partial \sigma_{pq}} D_{pqkl}
=
& &
A_{kl}
\frac{\partial f}{\partial I_1}\left( G \left( \frac{\partial I_1}
{\partial \sigma_{11}}\delta_{1l}
\delta_{1k}
+
\frac{\partial I_1}{\partial\sigma_{22}}
\delta_{2l} \delta_{2k}
+
\frac{\partial I_1}{\partial \sigma_{33}}
\delta_{3l} \delta_{3k} \right) \right.
\nonumber \\
&+&
\left.
\left(K-\frac{2}{3} G \right)
\frac{\partial I_1}{\partial \sigma_{cd}}
\delta_{cd} \delta_{kl} \right)
\nonumber \\
&+&
\frac{\partial f}{\partial \sqrt{J_2}}
\left( G \frac{\partial \sqrt{J_1}}{\partial \sigma_{ij}}
\delta_{ik}
\delta_{jl}
+
\left(K-\frac{2}{3} G \right)
\frac{\partial \sqrt{J_2}}{\partial \sigma_{ab}}
\delta_{ab}
\delta_{kl}
\right)
%\label{eqno_A3A}
%&= -3K\alpha \delta_{kl} + \displaystyle \frac{G}{\sqrt{J_2}}S_{kl}
\nonumber \\
\label{eqno_A3}
\end{eqnarray}
%
and,
%
\begin{eqnarray}
B
&=&
\frac{\partial f}{\partial \sigma_{rs}}D_{rstu}\frac{\partial f}{\partial \sigma_{tu}}
= 
\nonumber \\ 
& &
\left(
\frac{\partial f}{\partial I_1}
\right)^2
\left(
G \left(
\left(
\frac{
\partial I_1}{\partial \sigma_{11}}
\right)^2
+
\left(
\frac{\partial I_1}{\partial \sigma_{22}}
\right)^2
+
\left(
\frac{\partial I_1}{\partial \sigma_{33}}
\right)^2
\right)
\right.
\nonumber \\
&+&
\left.
\left(
K - \frac{2}{3} G
\right)
\left(
\frac{\partial I_1}{\partial \sigma_{ij}}
\delta_{ij}
\right)^2
\right.
\nonumber \\
&+&
\left.
\left(
\frac{\partial f}{\partial \sqrt{J_2}}
\right)^2
\left(
G \frac{\partial \sqrt{J_2}}{\partial \sigma_{ij}}
\frac{\partial \sqrt{J_2}}{\partial \sigma_{ij}}
+
\left(K - \frac{2}{3} G \right)
\left(
\frac{\partial \sqrt{J_2}}{\partial \sigma_{ij}}
\delta_{ij}
\right)^2
\right)
\right)
\nonumber \\
%\label{eqno_A4A}
%&= 9K\alpha ^2 + G
\label{eqno_A4}
\end{eqnarray}
%
where,    $K$ and $G$ are the elastic bulk modulus and the elastic shear modulus
respectively.



By further assuming that the evolution of internal variable is a function
of        equivalent      plastic      strain\footnote{This  is a fairly common
assumption,  e.g.  \citep{book:Chen}},     $e_{eq}^p ={2}/{3} e_{ij}^p  e_{ij}^p$
then one can write
%
\begin{equation}
K_P
=
- \frac{\partial f}{\partial q_n} \; r_n
=
- \displaystyle
\frac{1}{\sqrt{3}}
\; \displaystyle \frac{\partial f}{\partial \alpha}
\; \displaystyle \frac{d \alpha}{de^p_{eq}}
\; \frac{\partial f}{\partial \sqrt{J_2}}
\label{eqno_A5}
\end{equation}
%
It  should  be  noted  that since material properties are assumed to be  random,  the
resulting stress  tensor  will  also become random and hence the derivatives of the stress
invariants  with  respect  to  stress  tensor  ($\sigma_{ij}$)  will  become
%\sout{stochastic} 
{random}.   Therefore,    differentiations   appearing  in
Eqs.~(\ref{eqno_A3}), ~(\ref{eqno_A4}),  and ~(\ref{eqno_A5}) can not be carried out in
a deterministic sense.
%


The parameter tensor in Eq.~(\ref{eqno_S2}) then becomes
%
\begin{eqnarray}
D_{ijkl}^{ep} =
\left\{\begin{array}{ll}
     G
     \delta_{ik}
     \delta_{jl}
     +
     \left(K-\displaystyle \frac{2}{3}G\right)
     \delta_{ij}
     \delta_{kl}
        \;\;\; & \mbox{~;~} f <  0 \vee (f = 0 \wedge df <  0)
        \\ \\
    G
    \delta_{ik}
    \delta_{jl}
    +
    \left( K - \displaystyle \frac{2}{3} G \right)
    \delta_{ij}
    \delta_{kl}
    -
    \displaystyle \frac{A_{ij} A_{kl}}{B + K_P}
        \;\;\; & \mbox{~;~} f = 0 \vee df = 0
                       \end{array} \right.
\label{StochasticTangetnStiffness}
\end{eqnarray}
%
where  tensor    $A_{ij}$  and  scalars  $B$    and      $K_P$  are  defined  by
Eqs.~(\ref{eqno_A3}),  ~(\ref{eqno_A4}),  and ~(\ref{eqno_A5}) respectively. The
above equation (Eq.~\ref{StochasticTangetnStiffness}) represents a probabilistic
continuum  stiffness  tensor  for  an elastic--plastic material model, in this
case Drucker-Prager isotropic linear hardening material with associated plasticity.
%
By focusing  our      attention   on  one  dimensional  {point-location scale}  shear 
constitutive   relationship   between   $\sigma_{12}$  and  $\epsilon_{12}$  for
Drucker-Prager      material model,      one     can     simplify     the     function
$\eta(\sigma,D,q,r,\epsilon;  x,  t)$ as defined in Eq.~(\ref{eqno_S4}) (on Page
\pageref{eqno_S4}) to read
%
\begin{eqnarray}
\eta
=
\left\{\begin{array}{ll}
    G \displaystyle \frac{d \epsilon_{12}}{dt}
        \;\;\; & \mbox{~;~} f <  0 \vee (f = 0 \wedge df <  0)

        \\
%
    \left( G
           -
           \displaystyle \frac{G^2
            \left(\displaystyle \frac{\partial f}{\partial \sqrt{J_2}}
                  \displaystyle \frac{\sqrt{J_2}}{\partial \sigma_{12}}
            \right)^2}
            {B+K_P}\right)
            \displaystyle \frac{d \epsilon_{12}}{dt}
        \;\;\; & \mbox{~;~} f = 0 \vee df = 0
\end{array} \right.
\end{eqnarray}
%
%
By considering both the material properties (shear modulus $G$, bulk
modulus  $K$,  friction angle $\alpha$, and rate of change of friction angle
(linear  hardening)  $\alpha'    $)   and  the  strain  rate ($d\epsilon_{12}/dt
(t)$) as random,  one  can  substitute $\eta$ as derived in Eq.~(\ref{eqno_S12}) to
obtain  the  particular FPK equation for the probabilistic behavior of
Drucker-Prager  associative  linear  hardening,  1-D  {point-location scale} elastic-plastic shear
constitutive  rate  equation.  In particular, two cases are recognized, one for
elastic (pre--yield) behavior of material ($f <  0 \vee (f = 0 \wedge df <  0$))
%
\begin{eqnarray}
{
\lefteqn{\displaystyle \frac{\partial P(\sigma_{12}(t), t)}{\partial t}=}
}
\nonumber
\\
&-&
{
\displaystyle \frac{\partial}{\partial \sigma_{12}}
\left[\left<G \displaystyle \frac{d\epsilon_{12}}{dt} (t) \right>
P(\sigma_{12}(t),t) \right]
}
\nonumber
\\
&+&
{
\frac{\partial^2}{\partial \sigma_{12}^2}\left [\left\{\int_{0}^{t} d \tau Cov_0
\left[G \displaystyle \frac{d\epsilon_{12}}{dt} (t);
G \displaystyle \frac{d\epsilon_{12}}{dt}
(t-\tau)\right] \right\}P(\sigma_{12}(t),t) \right]
}
%\nonumber
%\\
\label{eqno_A10}
\end{eqnarray}
%
noting  that  this  is  the  same equation as Eq.~(\ref{eqno_LE4}). In addition to
that,  the  case of elastic--plastic behavior ($f = 0 \vee df = 0$) is described
by the following probabilistic equation
%
\begin{eqnarray}
& & {
\frac{\partial P(\sigma_{12}(t), t)}{\partial t}
=
-\displaystyle \frac{\partial}{\partial \sigma_{12}}\left[ \left\{\left<
\left(
G^{ep}(t)
\right)
\frac{d\epsilon_{12}}{dt} (t) \right> \right. \right.
}
\nonumber
\\
&+& {
\left. \left.
\int_{0}^{t} d\tau Cov_0
\left[\frac{\partial}{\partial \sigma_{12}}
\left(
G^{ep}(t)
\frac{d\epsilon_{12}}{dt} (t)\right);
G^{ep}(t-\tau)
%\right. \right. \right.
%\nonumber
%\\
%&&
%\left. \left. \left.
\frac{d\epsilon_{12}}{dt}
(t-\tau) \right] \right \} P(\sigma_{12}(t),t) \right]
}
\nonumber
\\
&+& {
\frac{\partial^2}{\partial \sigma_{12}^2} \left[ \left\{\int_{0}^{t} d \tau Cov_0
\left[
%\left(
G^{ep}(t)
%\right) 
\frac{d\epsilon_{12}}{dt} (t);
%\right. \right. \right. \\
%\nonumber
%&& \left. \left. \left.
%\left(
G^{ep}(t-\tau)
%\right)
%\right. \right. \right.
%\\
%&& \left. \left. \left.
%\left. \left. \left.
%\left(G(x_{t-\tau}) -\frac{G^2 (x_{t-\tau})}{G(x_{t-\tau})+9K(x_{t-\tau})\alpha^2(x_{t-\tau}) + \displaystyle
%\frac{1}{\sqrt{3}}I_1(x_{t-\tau},t-\tau) \alpha'(x_{t-\tau})}\right)
% \right. \right. \right.
 \displaystyle \frac{d\epsilon_{12}}{dt}
(t-\tau) \right] \right \} P(\sigma_{12}(t),t) \right]
}
\nonumber \\
\label{eqno_A10a}
\end{eqnarray}
%
where  
%\sout{$G^{ep}(a,b)$} 
{$G^{ep}(a)$} is defined as probabilistic elastic--plastic kernel and is
introduced to shorten  the writing (but will also have other uses later)
%
\begin{eqnarray}
{
G^{ep}(a) =
\left(G -\frac{G^2
\left( \displaystyle  \frac{\partial  f}{\partial \sqrt{J_2}(a)}
\displaystyle
\frac{\partial \sqrt{J_2}(a)}{\partial  \sigma_{12}(a)}\right)^2}
{B(a)+K_P(a)}\right)
}
\label{ProbElPlKernel01}
\end{eqnarray}
%
and  $a$  
%\sout{and  $b$}  
assumes  values  
%\sout{$x_t$  or  $x_{t-\tau}$  and} 
$t$ or $t-\tau$.
%\sout{respectively}.


It  is  important to note that the differentiations appearing in the coefficient
terms  of the FPK PDE (Eq.~(\ref{eqno_A10a})),
within      the    probabilistic    elastic--plastic  kernel  
%\sout{$G^{ep}(a,b)$} 
{$G^{ep}(a)$} (i.e.
Eq.~(\ref{ProbElPlKernel01})),  are  for  fixed  values of $\sigma_{12}$ and hence
those          differentiations     can  be carried out in a deterministic sense.
After  carrying  out  the  differentiations,  the probabilistic elastic--plastic
kernel becomes
%
\begin{eqnarray}
{
\left. G^{ep}(a)\right|_{\sigma_{12} \rightarrow const.} =
\left(G -\frac{G^2 }{G+9K\alpha^2 
+
\frac{1}{\sqrt{3}}I_1(a)\alpha'}\right)
}
\label{ProbElPlKernel02}
\end{eqnarray}
%
%
%
which,  after  substitution,  result  in  simplification  of  the  FPK
equation~(\ref{eqno_A10a}).
Further  simplification is possible by noting that the first random 
%\sout{field} 
{process} in the
covariance    term    of    the    first    coefficient   on  the  r.h.s  of the
equation~(\ref{eqno_A10a})  is  independent of $\sigma_{12}$. Furthermore, since
the
covariance   of  zero with any random 
%\sout{field} 
{process}  is  zero, the FPK equation~(\ref{eqno_A10a}) is further simplified to
read
%
%
\begin{eqnarray}
{
\lefteqn{\displaystyle \frac{\partial P(\sigma_{12}(t), t)}{\partial t}=}
}
\nonumber
\\
&-& {
\displaystyle \frac{\partial}{\partial \sigma_{12}}\left[\left<
G^{ep}(t)
%\left(G (x_t)-\frac{G^2 (x_t)}{G(x_t)+9K(x_t)\alpha^2 (x_t)
%+\displaystyle \frac{1}{\sqrt{3}}I_1(x_t,t)\alpha'(x_t)}\right)
\displaystyle \frac{d\epsilon_{12}}{dt} (t) \right> P(\sigma_{12}
(t),t) \right]
}
\nonumber
\\
&+& {
\frac{\partial^2}{\partial \sigma_{12}^2} \left[ \left\{\int_{0}^{t} d \tau Cov_0
\left[
G^{ep}(t)
%\left(G (x_t) -\displaystyle \frac{G^2 (x_t)}{G(x_t)+9K(x_t)\alpha^2 (x_t) +
%\displaystyle \frac{1}{\sqrt{3}}I_1(x_t,t)\alpha'(x_t)}\right)
\displaystyle \frac{d\epsilon_{12}}{dt} (t); G^{ep}(t-\tau)
%\left(G(x_{t-\tau}) -\frac{G^2 (x_{t-\tau})}{G(x_{t-\tau})+
%9K(x_{t-\tau})\alpha^2(x_{t-\tau}) + \displaystyle
%\frac{1}{\sqrt{3}}I_1(x_{t-\tau},t-\tau) \alpha'(x_{t-\tau})}\right)
%\right. \right. \right. \\
%&& \left. \left. \left.
%   %
%   \vphantom{\left. \left. \left. \left(G(x_{t-\tau}) -\frac{G^2 (x_{t-\tau})}{G(x_{t-\tau})+9K(x_{t-\tau})\alpha^2(x_{t-\tau}) + \displaystyle
%   \frac{1}{\sqrt{3}}I_1(x_{t-\tau},t-\tau) \alpha'(x_{t-\tau})}\right) \right. \right. \right.}
%   %
\displaystyle \frac{d\epsilon_{12}}{dt}
(t-\tau) \right] \right \} P(\sigma_{12}(t),t) \right]
}
\nonumber \\
\label{eqno_A10c}
\end{eqnarray}
%
where  the  probabilistic  elastic--plastic kernel 
%\sout{$G^{ep}(a,b)$} 
{$G^{ep}(a)$} is given by the
Eq.~(\ref{ProbElPlKernel02}).



The  evolution  of  a    mean value of shear stress $\sigma_{12}$ is obtained by
substituting  $\eta$ (derived for Drucker-Prager material in Eq.~(\ref{eqno_S19}))
as
%
\begin{eqnarray}
& & {
\displaystyle \frac{<d\sigma_{12} (t)>}{dt}=
\left <
G^{ep}(t)
%\left(G (x_t)-\frac{G^2 (x_t) \left( \displaystyle \frac{\partial f(x_t)}{\partial \sqrt{J_2}(x_t,t)} 
%\displaystyle \frac{\partial
%\sqrt{J_2}(x_t,t)}{\partial \sigma_{12}(x_t,t)}\right)^2}{B(x_t,t)+K_P(x_t,t)}\right) 
\displaystyle \frac{d\epsilon_{12}}{dt} (t)
\right>
} 
\nonumber
\\
& & { 
+ 
\int_{0}^{t} d\tau Cov_0 \left[\displaystyle \frac{\partial}{\partial \sigma_{12}}\left(
G^{ep}(t)
%
%\left(G (x_t)- \displaystyle \frac{G^2
%(x_t)\left( \displaystyle \frac{\partial f(x_t)}{\partial \sqrt{J_2}(x_t,t)} \displaystyle \frac{\partial \sqrt{J_2}(x_t,t)}{\partial
%\sigma_{12}(x_t,t)}\right)^2}{B(x_t,t)+K_P(x_t,t)}\right) 
%
\frac{d\epsilon_{12}}{dt} (t)\right); 
%\right. \\
%& & \left. 
G^{ep}(t-\tau)
%
%\left(G(x_{t-\tau}) -\frac{G^2 (x_{t-\tau}) \left( \displaystyle \frac{\partial f(x_{t-\tau})}{\partial \sqrt{J_2}(x_{t-\tau},
%t-\tau)} \displaystyle \frac{\partial \sqrt{J_2}(x_{t-\tau},t-\tau)}{\partial \sigma_{12}(x_{t-\tau},t-\tau)}\right)^2}{B(x_{t-\tau},
%t-\tau)+K_P(x_{t-\tau}, t-\tau)}\right)
%
\frac{d\epsilon_{12}}{dt} (t-\tau) \right]
}
%\nonumber
%\\
\label{eqno_A10d}
\end{eqnarray}


It is important to note that the derivatives appearing in the mean and covariance term of
the  above  Eulerian-Lagrangian  integro-differential  equation 
(Eq.~(\ref{eqno_A10d}) with the probabilistic elastic--plastic kernel defined
through the Eq.~(\ref{ProbElPlKernel02}))
are  
%\sout{stochastic} 
{random}
differentiations    and    need    to   be  treated  accordingly.  One  possible
%\sout{approach} 
{approach} to obtaining these differentiations 
%\sout{is by using} 
{could be} perturbation with respect to mean
\citep{Hori:2000} but the "closure  problem" will appear.  Hence,  in  this  study
the  evolution  of  mean  of $\sigma_{12}$ will be obtained by the expectation
operation on the PDF (Eq.~(\ref{eqno_S16})).


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

\section{Initial and Boundary Conditions for the probabilistic elastic--plastic PDE}
\label{MethodOfSolution}


The   PDE describing  the  probabilistic 
behavior  of constitutive rate equations can be written in the following general
form:
%
\begin{eqnarray}
%\begin{subequations}
\frac{\partial P(\sigma_{12},t)}{\partial t} 
&=& 
- 
\frac{\partial}{\partial
\sigma_{12}}\left\{ P(\sigma_{12},t) N_{(1)} \right\} 
+ 
\frac{\partial^2}{\partial \sigma_{12}^2} \left\{ P(\sigma_{12},t) N_{(2)} \right\} 
%\label{eqno_A11}
\nonumber
\\
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
&=& 
-
\displaystyle   \frac{\partial   }{\partial   \sigma_{12}}\left[P(\sigma_{12},t)
N_{(1)}
- \displaystyle \frac{\partial}
{\partial \sigma_{12}}\left\{ P(\sigma_{12},t) N_{(2)} \right\} \right] 
%\label{eqno_A12}
\nonumber
\\
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
&=& 
-
\displaystyle \frac{\partial \zeta}{\partial \sigma_{12}}
%\label{eqno_A13}
\label{eqno_A11}
%\end{subequations}
\end{eqnarray}
%
where,    $N_{(1)}$  and $N_{(2)}$ are coefficients\footnote{Indices in brackets
are not used in index summation convention.} of the PDE and
represent  the expressions within the curly braces of the first and second terms
respectively       on       the       right--hand--side       of       Eqs.
(~\ref{eqno_LE4}),~(\ref{eqno_A10}), and~(\ref{eqno_A10a})
%
These terms  are called the advection ($N_{(1)}$) and diffusion ($N_{(2)}$) coefficients
as  the form of Eq.~(\ref{eqno_A11}) closely resembles advection--diffusion
equation  \citep{book:Gardiner}.  The symbol $\zeta$ in Eq.~(\ref{eqno_A11}) can be
considered    to    be   the   probability   current.   This   follows   from   
Eq.~(\ref{eqno_A11}),  which    is  a  continuity    equation   and  the  state 
variable of the equation is probability density.


After  introducing     initial   and   boundary   conditions,   one  can  solve 
Eq.~(\ref{eqno_A11}) for probability densities of $ \sigma_{12} $ with evolution
of  time.   The initial condition could be deterministic or stochastic depending
on the   type  of  problem. 
%
For  probabilistic  behavior  of  linear  elastic
constitutive  rate  equation  (Eq.~(\ref{eqno_LE4})),  one can assume that all the
probability  mass at time $t=0$ is concentrated at $\sigma_{12}=0$ or at
some constant value of $\sigma_{12}$ if there were some initial stresses to begin
with (e.g. overburden pressure on a soil mass).
 
In mathematical term, this translates to,
%
\begin{equation}
\label{eqno_A13a}
P(\sigma_{12},0)=\delta(\sigma_{12})
\end{equation}
%
where, $\delta(\cdot)$ is the Dirac delta function. 


For  the  post--yield  behavior  of probabilistic elastic-plastic  constitutive 
rate   Equation\footnote{Specialized   to   Drucker-Prager   associated  linear
hardening model.} (\ref{eqno_A10c}), there
will  be a distribution of $\sigma_{12}$, corresponding to the solution of the
pre-yield probabilistic behavior (Eq.~(\ref{eqno_A10})), to begin with.
%
%
%%THIS SENTENCE NEEDS SOME WORK
%\vspace*{1cm}
%THIS SENTENCE NEEDS SOME WORK
This  probability  mass ($P(\sigma_{12}(t), t)$), dictated
by      Eq.~(\ref{eqno_S12}),      will  advect  and  diffuse  into  the  domain
($\sigma_{12}, t$ space) of the system
throughout the evolution (in time/strain) of the simulation.
%THIS SENTENCE NEEDS SOME WORK
%\vspace*{1cm}
%
%THIS SENTENCE NEEDS SOME WORK
%
Since it is required that the probability mass within the system is conserved i.e.
no  leaking is allowed at the boundaries, a reflecting barrier at the boundaries
will  be  the  preferred  choice.  In  mathematical  term,  one can express this
condition as \citep{book:Gardiner}
%
\begin{equation}
\label{eqno_A13b}
\left. \zeta(\sigma_{12},t) \right|_{{At Boundaries}} =0
\end{equation}
%
In theory, the stress domain could extend from  $-\infty$ to $\infty$
so that boundary conditions are then
%
\begin{equation}
\label{eqno_A13c}
\zeta(-\infty,t)=\zeta(\infty,t)=0
\end{equation}

With these    initial    and    boundary    conditions,    the  probabilistic
differential  equation {(with random material 
properties and random strain)} for  elasto--plasticity,  specialized  in  this  case to
associated  Drucker--Prager  material  model with linear hardening, and by using
FPK  transform described above, can be solved for probability 
%\sout{density function for} 
{densities of} shear  stress  ($\sigma_{12}$) as it evolves 
{with time/shear strain ($\epsilon_{12}$)}.
%\sout{due to random field of material properties and strains}.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%\input{NumericalSchemeLinear}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%\input{Results}

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


\section{Discussions and Conclusions}
\label{Conclusions}

In  this  paper,  a second-order exact expression for evolution of PDF of stress
was   derived   for   any   general,   one-dimensional   {spatial-average   form}
elastic-plastic     constitutive     rate    equation.    The    advantage    of
Fokker--Planck--Kolmogorov based
approach  as  presented  here  is  that it doesn't suffer from "closure problem"
associated  with regular perturbation approach nor it requires repetitive use of
computationally expensive deterministic elastic-plastic model as associated with
Monte   Carlo   simulation  technique.  Furthermore,  the  developed  expression
(Eq.~(\ref{eqno_S12}))  is  linear  and  deterministic  PDE whereas the original
constitutive  rate  equation (Eq.~(\ref{eqno_S3})) was stochastic and non-linear
ODE. 
%
This    deterministic linearity  will  in   turn   provide  great  simplicity 
in  solving  the  PDE  (Eq.~(\ref{eqno_S12}))  for  probability densities and 
subsequently mean and variance   behavior   of  the  state  variable  (stress)  
of  the  original  ODE (Eq.~(\ref{eqno_S3})).  
%
An  associated Drucker--Prager material model with uncertain parameters was used
to  illustrate  the  general  form of probabilistic elastic and elastic--plastic
equations on one dimensions. 


In    a  companion  paper,  solution,  discussion and illustrative examples
are presented.


%In  simulation  of  boundary  value  problems with uncertain material properties
%stochastic  finite  element  method is a popular technique. In stochastic finite
%element  formulation  input material properties random fields are discretized to
%random  variables  by  standard  techniques.  But for non-linear elastic-plastic
%problems  they  need  to  be  re-evaluated  for each load steps. In other words,
%probabilistic  constitutive  integration  needs  to  be carried out at the Gauss
%points.  In  order  to  show  the  applicability of the developed methodology in
%carrying  out  probabilistic  constitutive  integration,  needed  for stochastic
%elastic-plastic   finite   element   method,   the  derived  general  expression
%(Eq.~(\ref{eqno_S12}))    was   specialized   to   the   particular   cases   of
%{point-location  scale} linear elastic (Eq.~(\ref{eqno_LE4})) and Drucker-Prager
%associated  linear  hardening  type  elastic-plastic (Eqs.~(\ref{eqno_A10a}) and
%~(\ref{eqno_A10c})) constitutive equations.


       
%In the companion paper, the solution of the derived PDE (Eq.~(\ref{eqno_S12})) for
%PDF  and  subsequently  mean  and  standard  deviation  will be obtained for two
%particular cases of {point-location scale} linear elastic (Eq.~(\ref{eqno_LE4}))
%and   Drucker-Prager   associative   linear   hardening   type   elastic-plastic
%(Eqs.~(\ref{eqno_A10a}) and ~(\ref{eqno_A10c})) constitutive rate equations with
%input  parameters  modeled  as  random  variable  and deterministic strain rate.
%Future   work   is   going   on   to  extend  the  one-dimensional  relationship
%(Eq.~(\ref{eqno_S12}))  to  a  three-dimensional setting and incorporating it in
%the formulation of three-dimensional
%%\sout{global level (load-deflection behavior) of stochastic
%%continuum mechanics problem} 
%{stochastic elastic-plastic finite element method}.


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



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

%\bibliography{SFEM,refmech}
%\bibliographystyle{unsrtnat}

\begin{thebibliography}{57}
\providecommand{\natexlab}[1]{#1}
\providecommand{\url}[1]{\texttt{#1}}
\expandafter\ifx\csname urlstyle\endcsname\relax
  \providecommand{\doi}[1]{doi: #1}\else
  \providecommand{\doi}{doi: \begingroup \urlstyle{rm}\Url}\fi

\bibitem[Zadeh(1983)]{Zadeh:1983}
J.~Zadeh.
\newblock The role of fuzzy logic in the management of uncertainty in expert
  systems.
\newblock \emph{Fuzzy Sets and Systems}, 11:\penalty0 199--228, 1983.

\bibitem[Ben-Haim and Elishakoff(1990)]{book:Elishakoff}
Y.~Ben-Haim and I.~Elishakoff.
\newblock \emph{Convex Models of Uncertainty in Applied Mechanics}.
\newblock Elsevier, Amsterdam, 1990.

\bibitem[Moore(1979)]{book:Moore}
R.~Moore.
\newblock \emph{Methods and Applications of Interval Analysis}.
\newblock SIAM, Philadelphia, 1979.

\bibitem[Vanmarcke(1983)]{book:Vanmarcke}
Erik Vanmarcke.
\newblock \emph{Random Fields: Analysis and Synthesis}.
\newblock The MIT Press, Cambridge, Massachusetts, 1983.

\bibitem[Lumb(1966)]{Lumb:1966}
P.~Lumb.
\newblock The variability of natural soils.
\newblock \emph{Canadian Geotechnical Journal}, 3:\penalty0 74--97, 1966.

\bibitem[Vanmarcke(1977)]{Vanmarcke:1977}
E.~H. Vanmarcke.
\newblock Probabilistic modeling of soil profiles.
\newblock \emph{Journal of Geotechnical Engineering Division, ASCE},
  109\penalty0 (5):\penalty0 1203--1214, 1977.

\bibitem[Mayerhoff(1993)]{Mayerhoff:1993}
G.~G. Mayerhoff.
\newblock Development of geotechnical limit state design.
\newblock In \emph{Proceedingd of International Symposium on Limit State Design
  in Geotechnical Engineering, Cpenhagen}, pages 1--12. ASCE, ASCE, New York,
  1993.

\bibitem[DeGroot and Baecher(1993)]{Baecher:1993}
Don~J. DeGroot and Gregory~B. Baecher.
\newblock Estimating autocovariance of in-situ soil properties.
\newblock \emph{Journal of Geotechnical Engineering}, 119\penalty0
  (1):\penalty0 147--166, January 1993.

\bibitem[Popescu(1995)]{Popescu:1995}
R.~Popescu.
\newblock \emph{Stochastic variability of soil properties: data analysis,
  digital simulation, effects on system behavior}.
\newblock Doctoral dissertation, Princeton University, Princeton, NJ, 1995.

\bibitem[Lacasse and Nadim(1996)]{Lacasse:1996}
Suzanne Lacasse and Farrokh Nadim.
\newblock Uncertainties in characterizing soil properties.
\newblock In Charles~D. Shackelford and Priscilla~P. Nelson, editors,
  \emph{Uncertainty in Geologic Environment: From Theory to Practice,
  Proceedings of Uncertainty '96, July 31-August 3, 1996, Madison, Wisconsin},
  volume~1 of \emph{Geotechnical Special Publication No. 58}, pages 49--75.
  ASCE, ASCE, New York, 1996.

\bibitem[Popescu et~al.(1998)Popescu, Prevost, and Deodatis]{Popescu:1998}
R.~Popescu, J.~H. Prevost, and G.~Deodatis.
\newblock Spatial variability of soil properties: Two case studies.
\newblock In Panos Dakoulas and Mishac Yegian, editors, \emph{Geotechnical
  Earthquake Engineering and Soil Dynamics III: Proceedings of a Speciality
  Conference held in Seattle, Washington, August 3-6, 1998}, Geotechnical
  Special Publication No. 75, pages 568--579. ASCE and Geo-Institute, ASCE,
  Reston, VA, 1998.

\bibitem[Phoon and Kulhawy(1999{\natexlab{a}})]{Kulhawy:1999A}
Kok-Kwang Phoon and Fred~H. Kulhawy.
\newblock Characterization of geotechnical variability.
\newblock \emph{Canadian Geotechnical Journal}, 36:\penalty0 612--624,
  1999{\natexlab{a}}.

\bibitem[Phoon and Kulhawy(1999{\natexlab{b}})]{Kulhawy:1999B}
Kok-Kwang Phoon and Fred~H. Kulhawy.
\newblock Evaluation of geotechnical property variability.
\newblock \emph{Canadian Geotechnical Journal}, 36:\penalty0 625--639,
  1999{\natexlab{b}}.

\bibitem[Fenton(1999{\natexlab{a}})]{Fenton:1999A}
Gordan~A. Fenton.
\newblock Estimation of stochastic soil models.
\newblock \emph{Journal of Geotechnical and Geoenvironmental Engineering},
  125\penalty0 (6):\penalty0 470--485, June 1999{\natexlab{a}}.

\bibitem[Fenton(1999{\natexlab{b}})]{Fenton:1999B}
Gordan~A. Fenton.
\newblock Random field modeling of cpt data.
\newblock \emph{Journal of Geotechnical and Geoenvironmental Engineering},
  125\penalty0 (6):\penalty0 486--498, June 1999{\natexlab{b}}.

\bibitem[Duncan(2000)]{Duncan:2000}
J.~M. Duncan.
\newblock Factors of safety and reliability in geotechnical engineering.
\newblock \emph{ASCE Journal of Geotechnical and Geoenvironmental Engineering},
  126\penalty0 (4):\penalty0 307--316, April 2000.

\bibitem[Rackwitz(2000)]{Rackwitz:2000}
R.~Rackwitz.
\newblock Reviewing probabilistic soil modelling.
\newblock \emph{Computers and Geotechnics}, 26:\penalty0 199--223, 2000.

\bibitem[Marosi and Hiltunen(2004)]{Marosi:2004}
Karen~T. Marosi and Dennis~R. Hiltunen.
\newblock Characterization of spectral analysis of surface waves shear wave
  velocity measurement uncertainty.
\newblock \emph{Journal of Geotechnical and Geoenvironmental Engineering},
  130\penalty0 (10):\penalty0 1034--1041, October 2004.

\bibitem[Stokoe~II et~al.(2004)Stokoe~II, Darendeli, Gilbert, Menq, and
  Choi]{Stokoe:2004}
K.~H. Stokoe~II, R.~B. Darendeli, R.~B. Gilbert, F.-Y. Menq, and W.~K. Choi.
\newblock Development of a new family of normalized modulus reduction and
  material damping curves.
\newblock In \emph{International Workshop on Uncertainties in Nonlinear Soil
  Properties and their Impact on Modeling Dynamic Soil Response, PEER
  Headquarters, UC Berkeley, March 18-19, 2004}, Plenary Paper,
  http://peer.berkeley.edu/lifelines/
  Workshop304/pdf/Stokoe{\_}PlenaryPaper.pdf, 2004. PEER.

\bibitem[Soize(1994)]{Soize1994}
C.~Soize.
\newblock \emph{The Fokker-Planck Equation for stochastic dynamical systems and
  its explicit steady state solutions}.
\newblock World Scientific, Singapore, 1994.

\bibitem[Langtangen(1991)]{Langtangen:1991}
H.~Langtangen.
\newblock A general numerical solution method for fokker-planck equations with
  application to structural reliability.
\newblock \emph{Probabilistic Engineering Mechanics}, 6\penalty0 (1):\penalty0
  33--48, 1991.

\bibitem[Masud and Bergman(2005)]{Bergman:2005}
Arif Masud and Lawrence~A. Bergman.
\newblock Application of multi-scale finite element methods to the solution of
  the fokker-planck equation.
\newblock \emph{Computer Methods in Applied Mechanics and Engineering},
  194\penalty0 (1):\penalty0 1513--1526, April 2005.

\bibitem[Hopf(1952)]{Hopf:1952}
Eberhard Hopf.
\newblock Statistical hydromechanics and functional calculus.
\newblock \emph{Journal of rational Mechanics and Analysis}, 1:\penalty0
  87--123, June 1952.

\bibitem[Lee(1974)]{Lee:1974}
L.~C. Lee.
\newblock Wave propagation in a random medium: A complete set of the moment
  equations with different wavenumbers.
\newblock \emph{Journal of Mathematical physics}, 15\penalty0 (9):\penalty0
  1431--1435, September 1974.

\bibitem[Schueller(1997{\natexlab{a}})]{Schueller:1997}
G.~I. Schueller.
\newblock A state-of-the-art report on computational stochastic mechanics.
\newblock \emph{Probabilistic Engineering Mechanics}, 12\penalty0 (4):\penalty0
  197--321, 1997{\natexlab{a}}.

\bibitem[Paice et~al.(1996)Paice, Griffiths, and Fenton]{Griffiths:1996}
G.~M. Paice, D.~V. Griffiths, and Gordan~A. Fenton.
\newblock Finite element modeling of settelement on spatially random soil.
\newblock \emph{Journal of Geotechnical Engineering}, 122\penalty0
  (9):\penalty0 777--779, 1996.

\bibitem[Griffiths et~al.(2002)Griffiths, Fenton, and
  Manoharan]{Griffiths:2002}
D.~V. Griffiths, Gordan~A. Fenton, and N.~Manoharan.
\newblock Bearing capacity of rough rigid strip footing on cohesive soil:
  Probabilistic study.
\newblock \emph{Journal of Geotechnical and Geoenvironmental Engineering},
  128\penalty0 (9):\penalty0 743--755, 2002.

\bibitem[Fenton and Griffiths(2003)]{Griffiths:2003}
Gordan~A. Fenton and D.~V. Griffiths.
\newblock Bearing capacity prediction of spatially random $c-\phi$ soil.
\newblock \emph{Canadian Geotechnical Journal}, 40:\penalty0 54--65, 2003.

\bibitem[Fenton and Griffiths(2005)]{Griffiths:2005}
Gordan~A. Fenton and D.~V. Griffiths.
\newblock Three-dimensional probabilistic foundation settlement.
\newblock \emph{Journal of Geotechnical and Geoenvironmental Engineering},
  131\penalty0 (2):\penalty0 232--239, February 2005.

\bibitem[Popescu et~al.(1997)Popescu, Prevost, and Deodatis]{Popescu:1997}
R.~Popescu, J.~H. Prevost, and G.~Deodatis.
\newblock Effects of spatial variability on soil liquefaction: Some design
  recommendations.
\newblock \emph{Geotechnique}, 47\penalty0 (5):\penalty0 1019--1036, 1997.

\bibitem[Mellah et~al.(2000)Mellah, Auvinet, and
  Masrouri]{Mellah:SFEMgeotechApplication}
R.~Mellah, G.~Auvinet, and F.~Masrouri.
\newblock Stochastic finite element method applied to non-linear analysis of
  embankments.
\newblock \emph{Probabilistic Engineering Mechanics}, 15:\penalty0 251--259,
  2000.

\bibitem[De~Lima et~al.(2001)De~Lima, Teixeira, and
  Ebecken]{DeLima:SFEMgeotechApplication}
B.~S. L.~P. De~Lima, E.~C. Teixeira, and N.~F.~F. Ebecken.
\newblock Probabilistic and possibilistic methods for the elastoplastic
  analysis of soils.
\newblock \emph{Advances in Engineering Software}, 132:\penalty0 569--585,
  2001.

\bibitem[Koutsourelakis et~al.(2002)Koutsourelakis, Prevost, and
  Deodatis]{Deodatis:2002}
S.~Koutsourelakis, J.~H. Prevost, and G.~Deodatis.
\newblock Risk assesment of an interacting structure-soil system due to
  liquefaction.
\newblock \emph{Earthquake Engineering and Structural Dynamics}, 31:\penalty0
  851--879, 2002.

\bibitem[Nobahar(2003)]{Nobahar:2003}
A.~Nobahar.
\newblock \emph{Effects of Soil Spatial Variability on Soil-Structure
  Interaction}.
\newblock Doctoral dissertation, Memorial University, St. John's, NL, 2003.

\bibitem[Kleiber and Hien(1992)]{book:Kleiber}
M.~Kleiber and T.~D. Hien.
\newblock \emph{The Stochastic Finite Element Method: Basic Perturbation
  Technique and Computer Implementation}.
\newblock John Wiley {\&} Sons, Baffins Lane, Chichester, West Sussex PO19 1UD
  , England, 1992.

\bibitem[Der~Kiureghian and Ke(1988)]{DerKiureghian:1988}
A.~Der~Kiureghian and B.~J. Ke.
\newblock The stochastic finite element method in structural reliability.
\newblock \emph{Journal of Probabilistic Engineering Mechanics}, 3\penalty0
  (2):\penalty0 83--91, 1988.

\bibitem[Gutierrez and {De Borst}(1999)]{DeBorst:1999}
M.~A. Gutierrez and R.~{De Borst}.
\newblock Numerical analysis of localization using a viscoplastic
  regularization: Influence of stochastic material defects.
\newblock \emph{International Journal for Numerical Methods in Engineering},
  44:\penalty0 1823--1841, 1999.

\bibitem[Ghanem and Spanos(2003)]{book:Ghanem}
Roger~G. Ghanem and Pol~D. Spanos.
\newblock \emph{Stochastic Finite Elements: A Spectral Approach}.
\newblock Dover Publications, Inc., Mineola, New York, 2003.

\bibitem[Keese and Matthies(2002)]{keese:SPDE}
Andreas Keese and Hermann~G. Matthies.
\newblock Efficient solvers for nonlinear stochastic problem.
\newblock In H.~A. Mang, F.~G. Rammmerstorfer, and J.~Eberhardsteiner, editors,
  \emph{Proceedings of the Fifth World Congress on Computational Mechanics,
  July 7-12, 2002, Vienna, Austria},
  http://wccm.tuwien.ac.at/publications/Papers/fp81007.pdf, 2002.

\bibitem[Xiu and Karniadakis(2003)]{Karniadakis:2003}
D.~Xiu and G.~E. Karniadakis.
\newblock A new stochastic approach to transient heat conduction modeling with
  uncertainty.
\newblock \emph{International Journal of Heat and Mass Transfer}, 46:\penalty0
  4681--4693, 2003.

\bibitem[Debusschere et~al.(2003)Debusschere, Najm, Matta, Knio, and
  Ghanem]{Debusschere:2003}
B.~J. Debusschere, H.~N. Najm, A.~Matta, O.~M. Knio, and R.~G. Ghanem.
\newblock Protein labeling reactions in electrochemical microchannel flow:
  Numerical simulation and uncertainty propagation.
\newblock \emph{Physics of Fluids}, 15:\penalty0 2238--2250, 2003.

\bibitem[Anders and Hori(2000)]{Hori:2000}
Maciej Anders and Muneo Hori.
\newblock Three-dimensional stochastic finite element method for elasto-plastic
  bodies.
\newblock \emph{International Journal for Numerical Methods in Engineering},
  51:\penalty0 449--478, 2000.

\bibitem[Matthies et~al.(1997)Matthies, Brenner, Bucher, and
  Soares]{Matthies:1997reviewSFEM}
Hermann~G. Matthies, Christoph~E. Brenner, Christian~G. Bucher, and C.~Guedes
  Soares.
\newblock Uncertainties in probabilistic numerical analysis of structures and
  soilds - stochastic finite elements.
\newblock \emph{Structural Safety}, 19\penalty0 (3):\penalty0 283--336, 1997.

\bibitem[Deb et~al.(2001)Deb, Babuska, and Oden]{babuska:GalerkinSPDE}
Manas~K. Deb, Ivo~M. Babuska, and J.~Tinsley Oden.
\newblock Solution of stochastic partial differential equations using galerkin
  finite element techniques.
\newblock \emph{Computer Methods in Applied Mechanics and Engineering},
  190:\penalty0 6359--6372, 2001.

\bibitem[Babuska and Chatzipantelidis(2002)]{Babuska:SPDE}
Ivo Babuska and Panagiotis Chatzipantelidis.
\newblock On solving elliptic stochastic partial differential equations.
\newblock \emph{Computer Methods in Applied Mechanics and Engineering},
  191:\penalty0 4093--4122, 2002.

\bibitem[Liu and Der~Kiureghian(1991)]{DerKiureghian:1991}
P.~L. Liu and A.~Der~Kiureghian.
\newblock A finite element reliability of geometrically nonlinear uncertain
  structures.
\newblock \emph{Journal of Engineering Mechanics}, 117\penalty0 (8):\penalty0
  1806--1825, 1991.

\bibitem[Keese(2003)]{keese:2003ReportSFEM}
Andreas Keese.
\newblock A review of recent developments in the numerical solution of
  stochastic partial differential equations (stochastic finite elements.
\newblock Scientific Computing 2003-06, Deapartment of Mathematics and Computer
  Science, Technical University of Braunschweig, Brunswick, Germany, 2003.

\bibitem[Anders and Hori(1999)]{Hori:1999}
Maciej Anders and Muneo Hori.
\newblock Stochastic finite element method for elasto-plastic body.
\newblock \emph{International Journal for Numerical Methods in Engineering},
  46:\penalty0 1897--1916, 1999.

\bibitem[Kavvas(2003)]{Kavvas:2003}
M.~Levent Kavvas.
\newblock Nonlinear hydrologic processes: Conservation equations for
  determining their means and probability distributions.
\newblock \emph{Journal of Hydrologic Engineering}, 8\penalty0 (2):\penalty0
  44--53, March 2003.

\bibitem[Sudret and Der~Kiureghian(2000)]{Kiureghian:2000ReportSFEM}
B.~Sudret and A.~Der~Kiureghian.
\newblock Stochastic finite element methods and reliability: A state of the art
  report.
\newblock Technical Report UCB/SEMM-2000/08, University of California,
  Berkeley, 2000.

\bibitem[Kubo(1963)]{Kubo:1963}
R.~Kubo.
\newblock Stochastic liouville equations.
\newblock \emph{J. Math. Phys.}, 4\penalty0 (2):\penalty0 174--183, 1963.

\bibitem[Van~Kampen(1976)]{VanKampen:1976}
N.~G. Van~Kampen.
\newblock Stochastic differential equations.
\newblock \emph{Phys. Rep.}, 24:\penalty0 171--228, 1976.

\bibitem[Kavvas and Karakas(1996)]{Kavvas:1996}
M.~Levent Kavvas and A.~Karakas.
\newblock On the stochastic theory of solute transport by unsteady and steady
  groundwater flow in heterogeneous aquifers.
\newblock \emph{Journal of Hydrology}, 179:\penalty0 321--351, 1996.

\bibitem[Risken(1989)]{book:Risken}
H.~Risken.
\newblock \emph{The Fokker-Planck Equation: Methods of Solution Applications}.
\newblock Springer-Verlag, Berlin Heidelberg, second edition, 1989.

\bibitem[Gardiner(2004)]{book:Gardiner}
C.~W. Gardiner.
\newblock \emph{Handbook of Stochastic Methods for Physics, Chemistry and the
  Natural Science}.
\newblock Springer:Complexity. Springer-Verlag, Berlin Heidelberg, third
  edition, 2004.

\bibitem[Schueller(1997{\natexlab{b}})]{Schueller:1997report}
G.~I. Schueller.
\newblock A state-of-the-art report on computational stochastic mechanics.
\newblock \emph{Probabilistic Engineering Mechanics}, 12\penalty0 (4):\penalty0
  197--321, 1997{\natexlab{b}}.

\bibitem[Chen and Han(1988)]{book:Chen}
W.~F. Chen and D.~J. Han.
\newblock \emph{Plasticity for Structural engineers}.
\newblock Springer-Verlag, 1988.

\end{thebibliography}


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

\section{Appendix}
\label{Appendix}

\subsection{Derivation of $ A_{kl} $ (Eq.~(\ref{eqno_A3}))}


Chain  rule  of  differentiation  for  yield function $f=f(I_1, J_2)$ is used to
develop equation
%
\begin{equation}
\label{eqno_AP1}
\frac{\partial f}{\partial \sigma_{pq}}D_{pqkl} 
= 
\frac{\partial f}{\partial I_1}\frac{\partial I_1}{\partial \sigma_{pq}}D_{pqkl}
+
\frac{\partial f}{\partial \sqrt{J_2}}
\frac{\partial \sqrt{J_2}}{\partial \sigma_{pq}}D_{pqkl}
\end{equation}
                                                                                                                                                             
The first term in the R.H.S of Eq.~(\ref{eqno_AP1}) can be written as:
%
\begin{align}
\nonumber
%\label{eqno_AP2}
\frac{\partial f}{\partial I_1}\frac{\partial I_1}{\partial \sigma_{pq}}D_{pqkl}&=\frac{\partial f}{\partial I_1}\frac{\partial I_1}{\partial 
\sigma_{pq}}\left[G\delta_{pk}\delta{ql}+\left(K-\frac{2}{3}G\right)\delta_{pq}\delta_{kl}\right]\\
&=G\frac{\partial f}{\partial I_1}\frac{\partial I_1}{\partial \sigma_{pq}}\delta_{pk}\delta_{ql}+\left(K-\frac{2}{3}G\right)\frac{\partial 
f}{\partial I_1}\frac{\partial I_1}{\partial \sigma_{pq}}\delta_{pq}\delta_{kl}
\label{eqno_AP2}
\end{align}
%
where now one can write the first term on the R.H.S of Eq.~(\ref{eqno_AP2}) as, 
%
\begin{eqnarray}
\nonumber
\frac{\partial I_1}{\partial \sigma_{pq}}\delta_{pk}\delta_{ql}
&=&
\left(\frac{\partial I_1}{\partial \sigma_{1q}}\delta_{1k}+\frac{\partial 
I_1}{\partial \sigma_{2q}}\delta_{2k}+\frac{\partial I_1}{\partial \sigma_{3q}}\delta_{3k}\right)\delta_{ql}\\
\nonumber
&=&
\left(\frac{\partial I_1}{\partial \sigma_{1q}}\delta_{ql}\delta_{1k}+\frac{\partial I_1}{\partial \sigma_{2q}}\delta_{ql}\delta_{2k}+\frac{\partial 
I_1}{\partial \sigma_{3q}}\delta_{ql}\right)\delta_{3k}\\
\nonumber
&=&
\left(\frac{\partial I_1}{\partial \sigma_{11}}\delta_{1l}+\frac{\partial I_1}{\partial \sigma_{12}}\delta_{2l}+\frac{\partial
I_1}{\partial \sigma_{13}}\delta_{3l}\right)\delta_{1k}\\
\nonumber
& & +\left(\frac{\partial I_1}{\partial \sigma_{21}}\delta_{1l}+\frac{\partial I_1}{\partial \sigma_{22}}\delta_{2l}+\frac{\partial
I_1}{\partial \sigma_{23}}\delta_{3l}\right)\delta_{2k}\\
\nonumber
& & +\left(\frac{\partial I_1}{\partial \sigma_{31}}\delta_{1l}+\frac{\partial I_1}{\partial \sigma_{32}}\delta_{2l}+\frac{\partial
I_1}{\partial \sigma_{33}}\delta_{3l}\right)\delta_{3k}\\
%\nonumber
\label{eqno_AP3}
&=&
\frac{\partial I_1}{\partial \sigma_{11}}\delta_{1l}\delta_{1k}+\frac{\partial I_1}{\partial \sigma_{22}}\delta_{2l}\delta_{2k}+\frac{\partial 
I_1}{\partial \sigma_{33}}\delta_{3l}\delta_{3k}
%\label{eqno_AP3}
\end{eqnarray}
%
It  is  important to note that the differentiation is with respect to random 
stress    so that   terms  in  last  line  of previous equation cannot be simplified
further. 
%
Similarly, from Eq.~(\ref{eqno_AP2}) it follows that
%
\begin{align}
G\frac{\partial f}{\partial I_1}\frac{\partial I_1}{\partial \sigma_{pq}}\delta_{pk}\delta_{ql}&=
G\frac{\partial f}{\partial I_1}\left[\frac{\partial I_1}{\partial \sigma_{11}}\delta_{1l}\delta_{1k}+\frac{\partial I_1}{\partial 
\sigma_{22}}\delta_{2l}\delta_{2k}+\frac{\partial I_1}{\partial \sigma_{33}}\delta_{3l}\delta_{3k}\right]
\label{eqno_AP4}
\end{align}
%
where detailed derivations yield 
%
\begin{eqnarray}
\nonumber
\frac{\partial I_1}{\partial \sigma_{pq}}\delta_{pq}\delta_{kl}
&=&
\frac{\partial I_1}{\partial \sigma_{11}}\delta_{11}\delta_{kl}+\frac{\partial 
I_1}{\partial \sigma_{22}}\delta_{22}\delta_{kl}+\frac{\partial I_1}{\partial \sigma_{33}}\delta_{33}\delta_{kl}\\
\nonumber
& & +
2\left(\frac{\partial I_1}{\partial \sigma_{12}}\delta_{12}\delta_{kl}+\frac{\partial I_1}{\partial \sigma_{23}}\delta_{23}\delta_{kl}+\frac{\partial 
I_1}{\partial \sigma_{31}}\delta_{31}\delta_{kl}\right)\\
\label{eqno_AP5}
&=&
\frac{\partial I_1}{\partial \sigma_{11}}\delta_{kl}+\frac{\partial I_1}{\partial \sigma_{22}}\delta_{kl}+\frac{\partial I_1}{\partial 
\sigma_{33}}\delta_{kl}
\end{eqnarray}
%
so  that  the  second  part  of the  R.H.S. of Eq.~(\ref{eqno_AP2})) can be
written as
%
\begin{equation}
\label{eqno_AP6}
\left(K-\frac{2}{3}G\right)\frac{\partial f}{\partial I_1}\frac{\partial I_1}{\partial \sigma_{pq}}\delta{pq}\delta_{kl} = 
\left(K-\frac{2}{3}G\right)\frac{\partial f}{\partial I_1}\left[\frac{\partial I_1}{\partial \sigma_{11}}\delta_{kl}+\frac{\partial I_1}{\partial 
\sigma_{22}}\delta_{kl}+\frac{\partial I_1}{\partial \sigma_{33}}\right]\delta_{kl}
\end{equation}


First part of the R.H.S. of the Eq.~(\ref{eqno_AP1})) can be written as
%
\begin{align}
\nonumber
\frac{\partial f}{\partial I_1}\frac{\partial I_1}{\partial \sigma_{pq}}D_{pqkl} &= \frac{\partial f}{\partial I_1}G \left(\frac{\partial 
I_1}{\partial \sigma_{11}}\delta_{1l}\delta_{1k}+\frac{\partial I_1}{\partial \sigma_{22}}\delta_{2l}\delta_{2k}+\frac{\partial I_1}{\partial 
\sigma_{33}}\delta_{3l}\delta_{3k}\right)\\
\label{eqno_AP7}
&+\frac{\partial f}{\partial I_1}\left(K-\frac{2}{3}G\right)\left(\frac{\partial I_1}{\partial \sigma_{11}}+\frac{\partial I_1}{\partial 
\sigma_{22}}+\frac{\partial I_1}{\partial \sigma_{33}}\right)\delta_{kl}
\end{align}
%
while the second part of the same Eq.~(\ref{eqno_AP1})) can be written as
%
\begin{align}
\nonumber
\frac{\partial f}{\partial \sqrt{J_2}}\frac{\partial \sqrt{J_2}}{\partial \sigma_{pq}}D_{pqkl}&=\frac{\partial f}{\partial \sqrt{J_2}}\frac{\partial 
\sqrt{J_2}}{\partial 
\sigma_{pq}}\left[G\delta_{pk}\delta{ql}+\left(K-\frac{2}{3}G\right)\delta{pq}\delta_{kl}\right]\\
&=G\frac{\partial f}{\partial \sqrt{J_2}}\frac{\partial \sqrt{J_2}}{\partial 
\sigma_{pq}}\delta_{pk}\delta_{ql}+\left(K-\frac{2}{3}G\right)\frac{\partial 
f}{\partial \sqrt{J_2}}\frac{\partial \sqrt{J_2}}{\partial \sigma_{pq}}\delta_{pq}\delta_{kl}
\label{eqno_AP8}
\end{align}
%
The second term of the previous equation  (Eq.~(\ref{eqno_AP8})) then becomes 
%
\begin{eqnarray}
\frac{\partial \sqrt{J_2}}{\partial \sigma_{pq}}\delta_{pk}\delta_{ql}
= & &
\left(\frac{\partial \sqrt{J_2}}{\partial 
\sigma_{1q}}\delta_{1k}+\frac{\partial 
\sqrt{J_2}}{\partial \sigma_{2q}}\delta_{2k}
+
\frac{\partial \sqrt{J_2}}{\partial \sigma_{3q}}\delta_{3k}\right)\delta_{ql}
\nonumber
\\
= & &
\frac{\partial \sqrt{J_2}}{\partial \sigma_{1q}}\delta_{ql}\delta_{1k}
+
\frac{\partial \sqrt{J_2}}{\partial \sigma_{2q}}\delta_{ql}\delta_{2k}
+\frac{\partial \sqrt{J_2}}{\partial \sigma_{3q}}\delta_{ql}\delta_{3k}
\nonumber
\\
= & &
\left(\frac{\partial \sqrt{J_2}}{\partial \sigma_{11}}\delta_{1l}
+
\frac{\partial \sqrt{J_2}}{\partial \sigma_{12}}\delta_{2l}
+
\frac{\partial \sqrt{J_2}}{\partial \sigma_{13}}\delta_{3l}\right)\delta_{1k}
\nonumber
\\
&+&
\left(\frac{\partial \sqrt{J_2}}{\partial \sigma_{21}}\delta_{1l}
+
\frac{\partial \sqrt{J_2}}{\partial \sigma_{22}}\delta_{2l}
+
\frac{\partial \sqrt{J_2}}{\partial \sigma_{23}}\delta_{3l}\right)\delta_{2k}
\nonumber
\\
&+&
\left(\frac{\partial \sqrt{J_2}}{\partial \sigma_{31}}\delta_{1l}
+
\frac{\partial \sqrt{J_2}}{\partial \sigma_{32}}\delta_{2l}+\frac{\partial
\sqrt{J_2}}{\partial \sigma_{33}}\delta_{3l}\right)\delta_{3k}
\label{eqno_AP9}
\end{eqnarray}
%
which   can   be   used   to   write   the   first   term   on   the  R.H.S.  of
Eq.~(\ref{eqno_AP2})) as
%
\begin{align}
\nonumber
G\frac{\partial f}{\partial \sqrt{J_2}}\frac{\partial \sqrt{J_2}}{\partial \sigma_{pq}}\delta_{pk}\delta_{ql}&=
G\frac{\partial f}{\partial \sqrt{J_2}}\left[\frac{\partial \sqrt{J_2}}{\partial \sigma_{11}}\delta_{1l}\delta_{1k}+\frac{\partial 
\sqrt{J_2}}{\partial 
\sigma_{12}}\delta_{2l}\delta_{1k}+\frac{\partial \sqrt{J_2}}{\partial \sigma_{13}}\delta_{3l}\delta_{1k}\right]\\
\nonumber
&+G\frac{\partial f}{\partial \sqrt{J_2}}\left[\frac{\partial \sqrt{J_2}}{\partial \sigma_{21}}\delta_{1l}\delta_{2k}+\frac{\partial 
\sqrt{J_2}}{\partial 
\sigma_{22}}\delta_{2l}\delta_{2k}+\frac{\partial \sqrt{J_2}}{\partial \sigma_{23}}\delta_{3l}\delta_{2k}\right]\\
\label{eqno_AP10}
&+G\frac{\partial f}{\partial \sqrt{J_2}}\left[\frac{\partial \sqrt{J_2}}{\partial \sigma_{31}}\delta_{1l}\delta_{3k}+\frac{\partial 
\sqrt{J_2}}{\partial 
\sigma_{32}}\delta_{2l}\delta_{3k}+\frac{\partial \sqrt{J_2}}{\partial \sigma_{33}}\delta_{3l}\delta_{3k}\right]
\end{align}
%
and since
%
\begin{align}
\nonumber
\frac{\partial \sqrt{J_2}}{\partial \sigma_{pq}}\delta_{pq}\delta_{kl}&=\frac{\partial \sqrt{J_2}}{\partial 
\sigma_{11}}\delta_{11}\delta_{kl}+\frac{\partial 
\sqrt{J_2}}{\partial \sigma_{22}}\delta_{22}\delta_{kl}+\frac{\partial \sqrt{J_2}}{\partial \sigma_{33}}\delta_{33}\delta_{kl}\\
\nonumber
&+2\left(\frac{\partial \sqrt{J_2}}{\partial \sigma_{12}}\delta_{12}\delta_{kl}+\frac{\partial \sqrt{J_2}}{\partial 
\sigma_{23}}\delta_{23}\delta_{kl}+\frac{\partial 
\sqrt{J_2}}{\partial \sigma_{31}}\delta_{31}\delta_{kl}\right)\\
\label{eqno_AP11}
&=\frac{\partial \sqrt{J_2}}{\partial \sigma_{11}}\delta_{kl}+\frac{\partial \sqrt{J_2}}{\partial \sigma_{22}}\delta_{kl}+\frac{\partial 
\sqrt{J_2}}{\partial 
\sigma_{33}}\delta_{kl}
\end{align}
%
with 
%
\begin{equation}
\label{eqno_AP12}
\left(K-\frac{2}{3}G\right)\frac{\partial f}{\partial \sqrt{J_2}}\frac{\partial \sqrt{J_2}}{\partial \sigma_{pq}}\delta{pq}\delta_{kl} = 
\left(K-\frac{2}{3}G\right)\frac{\partial f}{\partial \sqrt{J_2}}\left[\frac{\partial \sqrt{J_2}}{\partial \sigma_{11}}\delta_{kl}+\frac{\partial 
\sqrt{J_2}}{\partial 
\sigma_{22}}\delta_{kl}+\frac{\partial \sqrt{J_2}}{\partial \sigma_{33}}\right]\delta_{kl}
\end{equation}
%
and with 
%
\begin{align}
\nonumber
\frac{\partial f}{\partial \sqrt{J_2}}\frac{\partial \sqrt{J_2}}{\partial \sigma_{pq}}D_{pqkl} &= 
G\frac{\partial f}{\partial \sqrt{J_2}}\left[\frac{\partial \sqrt{J_2}}{\partial \sigma_{11}}\delta_{1l}\delta_{1k}+\frac{\partial 
\sqrt{J_2}}{\partial 
\sigma_{12}}\delta_{2l}\delta_{1k}+\frac{\partial \sqrt{J_2}}{\partial \sigma_{13}}\delta_{3l}\delta_{1k}\right]\\
\nonumber
&+G\frac{\partial f}{\partial \sqrt{J_2}}\left[\frac{\partial \sqrt{J_2}}{\partial \sigma_{21}}\delta_{1l}\delta_{2k}+\frac{\partial 
\sqrt{J_2}}{\partial 
\sigma_{22}}\delta_{2l}\delta_{2k}+\frac{\partial \sqrt{J_2}}{\partial \sigma_{23}}\delta_{3l}\delta_{2k}\right]\\
\nonumber
&+G\frac{\partial f}{\partial \sqrt{J_2}}\left[\frac{\partial \sqrt{J_2}}{\partial \sigma_{31}}\delta_{1l}\delta_{3k}+\frac{\partial 
\sqrt{J_2}}{\partial 
\sigma_{32}}\delta_{2l}\delta_{3k}+\frac{\partial \sqrt{J_2}}{\partial \sigma_{33}}\delta_{3l}\delta_{3k}\right]\\
\label{eqno_AP13}
&+\left(K-\frac{2}{3}G\right)\frac{\partial f}{\partial \sqrt{J_2}}\left(\frac{\partial \sqrt{J_2}}{\partial \sigma_{11}}\delta_{kl}+\frac{\partial 
\sqrt{J_2}}{\partial 
\sigma_{22}}\delta_{kl}+\frac{\partial \sqrt{J_2}}{\partial \sigma_{33}}\right)\delta_{kl}
\end{align}
%
will finally yield Eq.~(\ref{eqno_AP1})),
%
\begin{align}
\nonumber
\frac{\partial f}{\partial \sigma_{pq}}D_{pqkl} &= 
\frac{\partial f}{\partial I_1}\left[G\left(\frac{\partial I_1}{\partial 
\sigma_{11}}\delta_{1l}\delta_{1k}+\frac{\partial I_1}{\partial\sigma_{22}}\delta_{2l}\delta_{2k}+\frac{\partial I_1}{\partial 
\sigma_{33}}\delta_{3l}\delta_{3k}\right)+\left(K-\frac{2}{3}G
\right)\frac{\partial I_1}{\partial \sigma_{cd}}\delta_{cd}\delta_{kl}\right]\\
\label{eqno_AP14}
&+\frac{\partial f}{\partial \sqrt{J_2}}\left[G\frac{\partial \sqrt{
J_1}}{\partial \sigma_{ij}}\delta_{ik}\delta_{jl}+\left(K-\frac{2}{3}G \right)\frac{\partial \sqrt{J_2}}{\partial 
\sigma_{ab}}\delta_{ab}\delta_{kl}\right]
\end{align}


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

\subsection{Derivation of B (Eq.~(\ref{eqno_A4}))}

%\Large{\bf{Derivation of B (Eq.~\ref{eqno_A4})}}:

Using Eq.~(\ref{eqno_AP14}), one can write:
%
\begin{align}
\nonumber
\frac{\partial f}{\partial \sigma_{rs}}D_{rstu}\frac{\partial f}{\partial \sigma_{tu}} & =  \left[\frac{\partial f}{\partial I_1} G 
\left(\frac{\partial I_1}{\partial 
\sigma_{11}}\delta_{1l}\delta_{1k}+\frac{\partial I_1}{\partial\sigma_{22}}\delta_{2l}\delta_{2k}+\frac{\partial I_1}{\partial 
\sigma_{33}}\delta_{3l}\delta_{3k}\right) \right.\\
\nonumber
&+ \left. \frac{\partial f}{\partial I_1} \left(K-\frac{2}{3}G \right)\left(\frac{\partial I_1}{\partial \sigma_{11}}\delta_{kl}+\frac{\partial 
I_1}{\partial\sigma_{22}}\delta_{kl}
+\frac{\partial I_1}{\partial \sigma_{33}}\delta_{kl}\right) \right. \\
\nonumber
&+ \left. G \frac{\partial f}{\partial \sqrt{J_2}}\left( \frac{\partial \sqrt{J_2}}{\partial \sigma_{11}}\delta_{1l}\delta_{1k}+ \frac{\partial 
\sqrt{J_2}}{\partial 
\sigma_{12}} \delta_{2l}\delta_{1k} + \frac{\partial \sqrt{J_2}}{\partial \sigma_{13}}\delta_{3l}\delta_{1k} \right. \right. \\
\nonumber
&+ \left. \left. \frac{\partial \sqrt{J_2}}{\partial \sigma_{21}} \delta_{1l}\delta_{2k} + \frac{\partial \sqrt{J_2}}{\partial 
\sigma_{22}}\delta_{2l}\delta_{2k} + 
\frac{\partial \sqrt{J_2}}{\partial \sigma_{23}} \delta_{3l}\delta_{2k} \right. \right. \\
\nonumber
&+ \left. \left. \frac{\partial \sqrt{J_2}}{\partial \sigma_{31}} \delta_{1l}\delta_{3k} + \frac{\partial \sqrt{J_2}}{\partial 
\sigma_{32}}\delta_{2l}\delta_{3k} + 
\frac{\partial \sqrt{J_2}}{\partial \sigma_{33}} \delta_{3l}\delta_{3k} \right) \right. \\
\label{eqno_AP15}
&+ \left. \frac{\partial f}{\partial \sqrt{J_2}} \left(K-\frac{2}{3}G \right)\left(\frac{\partial \sqrt{J_2}}{\partial 
\sigma_{11}}\delta_{kl}+\frac{\partial \sqrt{J_2}}
{\partial\sigma_{22}}\delta_{kl} +\frac{\partial \sqrt{J_2}}{\partial \sigma_{33}}\delta_{kl}\right) \right] \frac{\partial f}{\partial \sigma_{kl}}
\end{align}
%
where, by using
%
\begin{eqnarray}
\delta_{1l}\delta_{1k}\frac{\partial f}{\partial \sigma_{kl}} = \frac{\partial f}{\partial \sigma_{11}}
\;\;\;\;\;\; \mbox{;} \;\;\;\;\;\;
\delta_{2l}\delta_{2k}\frac{\partial f}{\partial \sigma_{kl}} = \frac{\partial f}{\partial \sigma_{22}}
\;\;\;\;\;\; \mbox{;} \;\;\;\;\;\;
\delta_{3l}\delta_{3k}\frac{\partial f}{\partial \sigma_{kl}} = \frac{\partial f}{\partial \sigma_{33}}
\nonumber
\\
\delta_{2l}\delta_{1k}\frac{\partial f}{\partial \sigma_{kl}} = \frac{\partial f}{\partial \sigma_{12}}
\;\;\;\;\;\; \mbox{;} \;\;\;\;\;\;
\delta_{3l}\delta_{1k}\frac{\partial f}{\partial \sigma_{kl}} = \frac{\partial f}{\partial \sigma_{13}}
\;\;\;\;\;\; \mbox{;} \;\;\;\;\;\;
\delta_{1l}\delta_{2k}\frac{\partial f}{\partial \sigma_{kl}} = \frac{\partial f}{\partial \sigma_{21}}
\nonumber
\\
\delta_{3l}\delta_{2k}\frac{\partial f}{\partial \sigma_{kl}} = \frac{\partial f}{\partial \sigma_{23}}
\;\;\;\;\;\; \mbox{;} \;\;\;\;\;\;
\delta_{1l}\delta_{3k}\frac{\partial f}{\partial \sigma_{kl}} = \frac{\partial f}{\partial \sigma_{31}}
\;\;\;\;\;\; \mbox{;} \;\;\;\;\;\;
\delta_{2l}\delta_{3k}\frac{\partial f}{\partial \sigma_{kl}} = \frac{\partial f}{\partial \sigma_{32}}
\label{eqno_AP24}
\end{eqnarray}
%
one can write the first part of the R.H.S. of Eq.~(\ref{eqno_AP15}) as
%
\begin{align}
\nonumber
& \frac{\partial f}{\partial I_1} G \left(\frac{\partial I_1}{\partial \sigma_{11}}\delta_{1l}\delta_{1k}+\frac{\partial 
I_1}{\partial\sigma_{22}}\delta_{2l}\delta_{2k}
+\frac{\partial I_1}{\partial \sigma_{33}}\delta_{3l}\delta_{3k}\right)\frac{\partial f}{\partial \sigma_{kl}} 
=
\\
&= G \left[ \left(\frac{\partial f}{\partial \sigma_{11}}\right)^2 + \left(\frac{\partial f}{\partial \sigma_{22}}\right)^2 + \left(\frac{\partial 
f}{\partial 
\sigma_{33}}\right)^2 \right] 
\nonumber
\\
&= G \left(\frac{\partial f}{\partial I_1}\right)^2 \left[\left(\frac{\partial I_1}{\partial \sigma_{11}}\right)^2 + \left(\frac{\partial 
I_1}{\partial \sigma_{22}}\right)^2 
+ \left(\frac{\partial I_1}{\partial \sigma_{33}}\right)^2 \right]
\label{eqno_AP25}
\end{align}
%
while the second part of the same equation (Eq.~(\ref{eqno_AP15})) becomes
%
\begin{eqnarray}
 G \frac{\partial f}{\partial \sqrt{J_2}}
 \left(   
\frac{\partial \sqrt{J_2}}{\partial \sigma_{11}}\delta_{1l}\delta_{1k}
+ 
\frac{\partial \sqrt{J_2}}{\partial \sigma_{12}} \delta_{2l}\delta_{1k} 
+ 
\frac{\partial \sqrt{J_2}}{\partial \sigma_{13}}\delta_{3l}\delta_{1k}
\right.
\nonumber
\\
\left.
+
\frac{\partial \sqrt{J_2}}{\partial \sigma_{21}} \delta_{1l}\delta_{2k} 
+ 
\frac{\partial \sqrt{J_2}}{\partial \sigma_{22}}\delta_{2l}\delta_{2k} 
+ 
\frac{\partial \sqrt{J_2}}{\partial \sigma_{23}} \delta_{3l}\delta_{2k}  
\right. 
\nonumber
\\
\left.
+ 
\frac{\partial \sqrt{J_2}}{\partial \sigma_{31}} \delta_{1l}\delta_{3k} 
+ 
\frac{\partial \sqrt{J_2}}{\partial \sigma_{32}}\delta_{2l}\delta_{3k} 
+ 
\frac{\partial \sqrt{J_2}}{\partial \sigma_{33}} \delta_{3l}\delta_{3k} 
\right) \frac{\partial f}{\partial \sigma_{kl}} 
\nonumber
\\
=  
G \left( \frac{\partial f}{\partial \sqrt{J_2}} \right)^2 
\left[\left(\frac{\partial \sqrt{J_2}}{\partial \sigma_{11}}\right)^2 
+ 
\left(\frac{\partial \sqrt{J_2}}{\partial \sigma_{12}}\right)^2 
+ 
\left(\frac{\partial \sqrt{J_2}}{\partial \sigma_{13}}\right)^2 
\nonumber
\right. 
\\
\left. 
+ 
\left(\frac{\partial \sqrt{J_2}}{\partial \sigma_{21}}\right)^2 
+ 
\left(\frac{\partial \sqrt{J_2}} {\partial \sigma_{22}}\right)^2 
+ 
\left(\frac{\partial \sqrt{J_2}}{\partial \sigma_{23}}\right)^2 
\nonumber
\right. 
\\
\left. 
+ 
\left(\frac{\partial \sqrt{J_2}}{\partial \sigma_{31}}\right)^2
+ 
\left(\frac{\partial \sqrt{J_2}}{\partial \sigma_{32}}\right)^2 
+ 
\left(\frac{\partial \sqrt{J_2}}{\partial \sigma_{33}}\right)^2 \right]
\label{eqno_AP26}
\end{eqnarray}
%
and since 
%
\begin{equation}
\label{eqno_AP27}
\frac{\partial f}{\partial \sigma_{kl}} \delta_{kl} 
= 
\frac{\partial f}{\partial \sigma_{11}} + \frac{\partial f}{\partial \sigma_{22}} + 
\frac{\partial f}{\partial \sigma_{33}}
\end{equation}
%
it follows that
%
\begin{align}
\nonumber
& \frac{\partial f}{\partial I_1} \left(K-\frac{2}{3}G \right)\left(\frac{\partial I_1}{\partial \sigma_{11}}\delta_{kl}+\frac{\partial 
I_1}{\partial\sigma_{22}}\delta_{kl}
+\frac{\partial I_1}{\partial \sigma_{33}}\delta_{kl}\right) \frac{\partial f}{\partial \sigma_{kl}} \\
\label{eqno_AP28}
&=\left(K-\frac{2}{3}G \right) \left(\frac{\partial f}{\partial I_1} \right)^2\left[\frac{\partial I_1}{\partial \sigma_{11}} + \frac{\partial 
I_1}{\partial \sigma_{22}} 
+ \frac{\partial I_1}{\partial \sigma_{33}} \right]^2 
\end{align}
%
and
%
\begin{align}
\nonumber
& \frac{\partial f}{\partial \sqrt{J_2}} \left(K-\frac{2}{3}G \right)\left(\frac{\partial \sqrt{J_2}}{\partial \sigma_{11}}\delta_{kl}+\frac{\partial 
\sqrt{J_2}}
{\partial\sigma_{22}}\delta_{kl} +\frac{\partial \sqrt{J_2}}{\partial \sigma_{33}}\delta_{kl}\right) \frac{\partial f}{\partial \sigma_{kl}} \\
\label{eqno_AP29}
&= \left(K-\frac{2}{3}G \right) \left(\frac{\partial f}{\partial \sqrt{J_2}}\right)^2 \left[\frac{\partial \sqrt{J_2}}{\partial \sigma_{11}} + 
\frac{\partial \sqrt{J_2}}{\partial \sigma_{22}} + \frac{\partial \sqrt{J_2}}{\partial \sigma_{33}} \right]^2
\end{align}

Therefore,  adding  Eqs.~(\ref{eqno_AP24}), ~(\ref{eqno_AP25}), ~(\ref{eqno_AP27}), and
~(\ref{eqno_AP28}) one can write:
%
\begin{align}
\nonumber
B = 
\frac{\partial f}{\partial \sigma_{rs}}D_{rstu}\frac{\partial f}{\partial \sigma_{tu}} 
&=
\left(\frac{\partial f}{\partial I_1}\right)^2
\left[
G\left\{
\left(\frac{\partial I_1}{\partial\sigma_{11}}\right)^2
+
\left(\frac{\partial I_1}{\partial\sigma_{22}}\right)^2
+
\left(\frac{\partial I_1}{\partial \sigma_{33}}\right)^2
\right\} 
\right.
\nonumber 
\\
& + 
\left.
\left(K-\frac{2}{3}G \right)\left\{\frac{\partial I_1}{\partial \sigma_{ij}}\delta_{ij}\right\}^2\right]
\nonumber
\\
&+
\left(\frac{\partial f}{\partial \sqrt{J_2}}\right)^2 \left[G\frac{\partial \sqrt{J_2}}{\partial \sigma_{ij}}\frac{\partial
\sqrt{J_2}}{\partial \sigma_{ij}} +\left(K-\frac{2}{3}G \right)\left\{\frac{\partial \sqrt{J_2}}{\partial
\sigma_{ij}}\delta_{ij}\right\}^2 \right]
\label{eqno_AP30}
\end{align}
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Derivation of $ K_P $ (Eq.~(\ref{eqno_A5}))}

%\Large{\bf{Derivation of $ K_P $ (Eq.~\ref{eqno_A5})}}

>From the consistency condition ($ \dot f = 0 $), one can derive
%
\begin{equation}
\label{eqno_AP32}
K_P = - \frac{\partial f}{\partial q_n}r_n
\end{equation}
%
with
%
\begin{equation}
\label{eqno_AP33}
\dot q_n = <L> r_n
\end{equation}
%
where, $ <\cdot> $ is the McCaulay bracket and $L$ is the loading index given by
%
\begin{equation}
\label{eqno_AP34}
L= \frac{1}{K_P} L_{ij} \dot \sigma_{ij}
\end{equation}
%
where,  $ L_{ij} $ is a vector normal to the loading surface in stress space and
dot represents time derivative.
%
For Drucker-Prager yield criteria (Eq.~(\ref{eqno_A1})) one can simplify the above
equation (Eq.~(\ref{eqno_AP32})) as
%
\begin{equation}
\label{eqno_AP35}
K_P = - \frac{\partial f}{\partial \alpha} \bar \alpha
\end{equation}
%
where,
%
\begin{equation}
\label{eqno_AP36}
\dot \alpha = <L> \bar \alpha
\end{equation}
%
By  assuming  that  the  friction  angle  like  internal  variable $\alpha$ is a
function of equivalent plastic strain ($e_{eq}^p$) i.e.
%
\begin{equation}
\label{eqno_AP37}
\alpha=\alpha(e_{eq}^p)
\end{equation}
%
one    can   write, (using the relationship: $\dot e_{eq}^p = {1}/{\sqrt{3}} <L>
\displaystyle {\partial f}/{\partial \sqrt{J_2}}$)
%
\begin{align}
\nonumber
\dot \alpha = \frac{d \alpha}{d e_{eq}^p}\dot e_{eq}^p 
            = \frac{1}{\sqrt{3}} \frac{d \alpha}{d e_{eq}^p} <L> \frac{\partial f}{\partial \sqrt{J_2}} 
            = <L> \bar \alpha
\label{eqno_AP38}
\end{align}
%
and therefore,
%
\begin{equation}
\label{eqno_AP39}
K_P 
= 
- \frac{\partial f}{\partial q_n}r_n 
= 
- \frac{\partial f}{\partial \alpha} \bar \alpha  
= 
- \displaystyle \frac{1}{\sqrt{3}}\displaystyle \frac{\partial f}{\partial \alpha} 
\displaystyle \frac{d \alpha}{de^p_{eq}}\frac{\partial f}{\partial \sqrt{J_2}}
\end{equation}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%$
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%$

                      
\end{document}
\bye
                                                                                                                                                             

