\documentclass[dvips]{nagauth}
%\documentclass[12pt,letterpaper]{article}
%\usepackage{times}
\usepackage{graphicx}
%\usepackage{amsmath}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ovo je predefinisanje datuma za cirilicu ( vidi Kopka93 strana 175 )
\def\today
{\number\day . \space \ifcase\month\or
January\or
February\or
March\or
April\or
May\or
June\or
July\or
August\or
September\or
October\or
November\or
December\fi,\space
\number\year}

%time of day as found in layout.tex
% ----------------------------------------------------------------------
%
% TIME OF DAY
%
\newcount\hh
\newcount\mm
\mm=\time
\hh=\time
\divide\hh by 60
\divide\mm by 60
\multiply\mm by 60
\mm=-\mm
\advance\mm by \time
\def\hhmm{\number\hh:\ifnum\mm<10{}0\fi\number\mm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%za obicne ps fajlove
%\graphicspath{{/home/jeremic/tex/works/psfigures}
%              {/home/jeremic/oofep/fem/ASCE95}
%              {/home/jeremic/oofep/mathematica/BModel}
%              {/home/jeremic/oofep/fem/LStest}
%              {/home/jeremic/oofep/fem/LSTconst}
%              {/home/jeremic/oofep/MGM}
%              {/home/jeremic/oofep/mathematica/BModel/HardFunc}
%              {/home/jeremic/tex/works/Proposals/1999/LDGeomechanics}
%              {/home/jeremic/BG/amblemi}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\newcommand{\nDarray}
%{{\bf \lowercase{n}\kern-0.05em \lower0.3ex\hbox{D}\lowercase{array}}}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% FEMtools LOGO definition
%\newcommand{\FEMtools}
%{{\bf \uppercase{F}\kern-0.25em\uppercase{E}\kern-0.25em\uppercase{M}\lowercase{tools}}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% FEI  LOGO definition
%\newcommand{\FEI}
%{{\sf \uppercase{F}\kern-0.20em\uppercase{E}\kern-0.20em\uppercase{I}}}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% draft special command for postscript       IZBACI ZA SLUZBENU VERZIJU
%
%
%  \special{!userdict begin
%  /bop-hook{
%  gsave
%  105 120 translate     %start position
%  45 rotate            % orientation
%  /Times-Roman findfont %font
%  150 scalefont         % scaling of font
%  setfont
%  0 0 moveto
%  0.95 setgray          % gray level ( 1 -> white ;  0 black )
%  (Draft Paper)    % text you want to see
%  show    % or: true charpath for hollow letters
%  %true charpath
%  stroke grestore}def end}
%  %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % draft special command for postscript       IZBACI ZA SLUZBENU VERZIJU
% %
% %
% \special{!userdict begin
% /bop-hook{
% gsave
% 200 100 translate     %start position
% 60 rotate            % orientation
% /Times-Roman findfont %font
% 250 scalefont         % scaling of font
% setfont
% 0 0 moveto
% 0.80 setgray          % gray level ( 1 -> white ;  0 black )
% (Draft)    % text you want to see
% show    % or: true charpath for hollow letters
% %true charpath
% stroke grestore}def end}
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%JB
%%% from UCthesis.sty
%%\setcounter{topnumber}{4}
%%\def\topfraction{.8}
%%\setcounter{bottomnumber}{1}
%%\def\bottomfraction{.8}
%%\setcounter{totalnumber}{4}
%%\def\textfraction{.2}
%%\def\floatpagefraction{.8}
%%\setcounter{dbltopnumber}{2}
%%\def\dbltopfraction{.7}
%%\def\dblfloatpagefraction{.5}
%%%JB
%%
%%
%%


\def\baselinestretch{2}


\begin{document}



\NAG{1}{6}{01}{28}{01}

\runningheads{Sett, Jeremi{\'c} and Kavvas}
{Nonlinear Hardening/Softening and Probabilistic Elasto-Plasticity}

\title{The Role of Nonlinear Hardening/Softening in Probabilistic Elasto-Plasticity}

\author{Kallol~Sett, Boris~Jeremi{\'c}\corrauth and M.~Levent~Kavvas}

\address{Department  of  Civil  and  Environmental  Engineering,  University  of
California, \\
Davis, CA 95616 \\~\\ {\bf Accepted for publication}}
                                 

\corraddr{ Boris Jeremi{\'c}, Department of Civil and Environmental Engineering,
University of California, 
One Shields Ave., Davis, CA 95616, \texttt{jeremic@ucdavis.edu}}

%\footnotetext[2]{Please ensure that you use the most up to date class file,
%available from the CFM Home Page at\\
%\texttt{http://www.interscience.wiley.com/jpages/1082-5010/}}

%\cgsn{NSF}{EEC-9701568}

\noreceived{}
\norevised{}
\noaccepted{}

%\begin{center}
%{\sc Version: \today{} , \hhmm{} }
%\end{center}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\newpage
%\tableofcontents
%\newpage


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

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

%\begin{abstract}
                        
\section{ABSTRACT}

The  elastic--plastic  modeling and simulations have been studied extensively in
the  last  century.  However, one crucial area of material modeling has received
very  little  attention. The uncertainties in material properties probably have the
largest influence on many aspects of structural and solids behavior. Despite its
importance, effects of uncertainties of material properties on overall response of
structures  and  solids  have  rarely  been studied. Most of the small number of
studies  on  effects  of material variabilities have used repetitive deterministic
models  through  Monte--Carlo type simulations. While this approach might appear
sound,  it cannot be both computationally efficient and statistically accurate
(have statistically appropriate number of data points).


Recently,      we      have      developed      a methodology   to   solve  the 
probabilistic  elastic--plastic differential equations. The methodology is based
on  Eulerian--Lagrangian  form  of  the Fokker--Planck--Kolmogorov equation and
provides  for  full  description  of  the probability density function (PDF) of
stress response for a given strain.


In        this        paper   we   describe our development in some details. In
particular,  we  investigate the effects of nonlinear hardening on predicted PDF
of  stress.  As  it  will  be  shown,  the  nonlinear  hardening  will  create a
discrepancy  between  the  most  likely  stress  solution  and the deterministic
solution.  This  discrepancy,  in fact, means that the deterministic solution is
not  the  most  likely  outcome  of the corresponding probabilistic solution if
material parameters are uncertain (and they always are, we just tend to simplify
that fact and use, for example, mean values for deterministic simulations).


A  number  of  examples  will  be  presented,  illustrating methodology and main
results, some of which are quite surprising as mentioned above.



\keywords{Elasto-Plasticity, Probability Theory, Fokker-Planck-Kolmogorov Equation}

%\end{abstract}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{INTRODUCTION}

Numerical  simulations  of  mechanical  behavior  of solids and
structures  made  of geomaterials are inherently uncertain. The
natural  conditions with which geotechnical engineers deal with
are  unknown  and  must  be inferred from limited and expensive
observations           (Beacher and Christian \cite{book:Christian}).           
For example
Fig.~\ref{figure:SoilVariability} shows measured soil properties in Mexico City.
%
%
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.7\textwidth]{SoilVariability}
%\vspace*{-0.8cm}
\caption{Measured values of mechanical properties of soil from Mexico city, typical soft spot
(after Baecher and Christian \cite{book:Christian}}
\label{figure:SoilVariability}
\end{center}
\end{figure}
%
%
Large variability in material properties can be noted as
a  function  of  depth. 
%This is   especially true for unconfined compressive
%strength,  preconsolidation  load,  and friction angle obtained
%from consolidated quick triaxial tests. 
It  is  important  to  note    that this  soil  formation  is  considered  very 
homogeneous  (Baecher  and Christian \cite{book:Christian}). The usual procedure
in this case is to choose
either  constant properties or fit smooth (usually linear) curves to properties
of  interest  (for example friction angle and/or unconfined compressive strength
and/or  preconsolidation  pressure (refer Fig.~\ref{figure:SoilVariability})) and
use the assumed (much simplified) properties in modeling and simulations related
to  this soil profile. This approach will completely neglect the variations which
might  have a large effect on response of this soil profile. The usual remedy is
to  apply  a  large  factor  of  safety,  which,  it is hoped, will cover up for
neglected  information.  Use of large factors of safety is becoming unacceptable
as  it leads to design solutions that are not only uneconomical but also sometimes not even
safe (Duncan \cite{Duncan:2000}). More appropriate methods for analyzing the random fields
of    material   properties   have   been   described   by   Ghanem and Spanos \cite{book:Ghanem},
Matthies  and  Keese \cite{Matthies2005}, Roberts and Spanos \cite{Roberts1986},
Zhu et al, \cite{Zhu1994}, Soize \cite{Soize1994}.

In      recent  years, civil engineering practice, and particularly geotechnical
engineering  practice  has  seen  an  increasing  emphasis  on quantification of
uncertainties. Quantifications or mathematical descriptions of uncertainties are
usually  done  within  the  framework  of  probability theory, where random soil
parameters  are  modeled as random variables (if they are specialized to a fixed
location  in  the  domain)  or  random  fields  (if  they are specialized to the
function  of  location  in  the domain). Nice descriptions of different types of
soil  parameter uncertainties and their quantification methods were presented by
Phoon   \&   Kulhawy  (\cite{Kulhawy:1999A},  \cite{Kulhawy:1999B})  and  Fenton
(\cite{Fenton:1999A},  \cite{Fenton:1999B}).  The  uncertainty  in $G/G_{max}$
curve was recently quantified by Stokoe et al. \cite{Stokoe:2004}.

 
Even-though  there is a  widely acknowledged necessity to model soil properties
using  theory  of  probability,  there  is  very  little  work on propagation of
uncertainties in soil properties through elastic--plastic
constitutive   equation.   Advanced   elasto--plasticity  based
constitutive  models,  when  properly  calibrated,  although very
accurate,  are, in general, highly sensitive to fluctuations in
model  parameters  (cf. Borja \cite{borja:OpinionPaper}).  In  the field of
geomechanics,  this sensitivity aspect of constitutive model is
of  great  importance as the uncertainties associated with soil
properties could outweigh the advantages gained by advanced and
sophisticated  models.  

First  attempt  to propagate randomness
through    the    elastic--plastic    constitutive   equations,
considering  random Young's modulus was published only recently
by  Anders  and Hori \cite{Hori:1999}. They based their work on perturbation
expansion at  the stochastic mean behavior.
In computing the mean behavior
they  took advantage of bounding media analysis. However, because of
the  use  of Taylor series expansion, this approach is limited to
problems  with small coefficients of variation. According to a recent report by
Sudret  and Der Kiureghian \cite{Kiureghian:2000ReportSFEM}, the coefficient    
of  variation  should  not  exceed  20~\%  if  this approach is to be used. This
limitation  restricts  the  use  of  perturbation  method to general geomechanics
problems, where the coefficients of variation of soil properties are rarely less
than  20\% (eg. Phoon and Kulhawy \cite{Kulhawy:1999A}). 
Another disadvantage of the perturbation method
is   that   it   inherits   the   so-called   "closure   problem"   (cf.  Kavvas
\cite{Kavvas:2003}), where information on higher--order moments is always needed
to  calculate  lower--order  moments.  More  recently,  Fenton   and  Griffiths 
\cite{Griffiths:2005} used a Monte--Carlo type method to propagate uncertainties
through
elastic--plastic   $c-\phi$   soil.   

Conventional    Monte--Carlo simulation,          although      very   accurate,
might   be    computationally   expensive   for  general  elastic--plastic
simulations. The method typically requires a statistically appropriate number of
realizations per random variable (often a very large number) in order to satisfy
statistical   accuracy.   This   repetitive  use  of  computationally  expensive
deterministic   model   renders   the   Monte--Carlo  simulation  impractical.  
Monte--Carlo   simulations,  however,  find  a  great  use  in  verification  of
analytical developments (cf. Oberkampf et al, \cite{Oberkampf2002}).

In   order  to  overcome the drawbacks associated with Monte--Carlo technique and
perturbation  method in dealing with non-linear stochastic differential equations
(SDEs),  Kavvas  \cite{Kavvas:2003} derived a generic Eulerian-Lagrangian form of
Fokker--Planck--Kolmogorov   equation   (FPKE)   for   the  second--order  exact
probabilistic  solution  of  any  non-linear ODE with stochastic coefficient and
stochastic  forcing.  The  development  was  focusing  on  applications in water
resources  engineering  but,  it  is  general  enough to apply to any
non-linear ODE with stochastic coefficient and stochastic forcing.
%
Using    the above mentioned  Eulerian--Lagrangian  form of FPKE, Jeremi{\'c} et
al. \cite{Jeremic:2005a} and Sett et al. \cite{Jeremic:2005b} recently developed
formulation  and  solution  for  the general  1-D  elastic-plastic constitutive rate
equation with random material properties and random strain rate.
%


The    main  advantage      of      FPKE-based    approach  is  that  a complete
probabilistic  description (probability density function (PDF)) of stress can be
obtained  exactly  to  second-order  accuracy (to covariance of time), given random material
property(ies) and/or random strain. In addition to that, even if the FPKE is not
solvable  in  closed-form  solution,  the deterministic linearity of the FPKE in
terms  of the state variable, the probability density of stress, considerably simplifies
the numerical solution process.



In    this    paper   we  present  a  general expression, in the form of a partial
differential  equation,  for  evolution  of  the probability  density of stress with
strain   for   any   general  3-D  point-location  scale/local-average  form  of
elastic-plastic  constitutive  law  with  material  properties  and  strain rate
modeled  as  random  variables/fields. The expression for the PDF of stress is
developed  in a general way, making it usable in any dimension (1--,2-- or 3--D)
and  for any type of incremental elastic--plastic material model. Application of
the  developed  methodology is demonstrated through 1-D point-location scale Cam-Clay
shear  constitutive  law.  In particular,  the effects of non--linear hardening and
softening on the predicted PDFs of stresses are investigated.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{GENERAL EXPRESSION FOR 3--D PROBABILISTIC 
CONSTITUTIVE RATE EQUATION IN LOCAL-AVERAGE FORM}
\label{MasterKey}

The    incremental   form  of  spatial-average  3-D elastic-plastic constitutive
rate equation can be written as:
%
\begin{equation}
\label{eqno_S2}
\frac{d\sigma_{ij}(x_t, t)}{dt}
= 
D^{ep}_{ijkl}(\sigma_{ij}, D^{el}_{ijkl}, f, U, q_*, r_*;x_t,t) 
\frac{d\epsilon_{kl}(x_t,t)}{dt}
\end{equation}
%
%
where,      $D_{ijkl}^{ep}$  is  the  {\it  random,  non-linear  elastic-plastic
coefficient  tensor}  which is a function of random stress tensor ($ \sigma_{ij}
$), random elastic moduli tensor ($D^{el}_{ijkl}$), random yield function ($f$),
random  potential  function  ($U$), random internal variables ($q_*$) and random
direction  of  evolution  of  internal  variables  ($r_*$).  The random internal
variables ($q_*$) could be scalar (for perfectly plastic and isotropic hardening
models),  or  second-order  tensor  (for  translational and rotational kinematic
hardening models), or fourth-order tensor (for distortional hardening models) or
any  combinations  of  the  above. The same classification applies to the random
direction of evolution of internal variables ($r_*$). By denoting all the random
material parameters by a parameter tensor as:
%
\begin{equation}
\label{eqno_S2a}
D_{ijkl} = \left[D^{el}_{ijkl}, f, U, q_*, r_*\right]
\end{equation}
%
and   introducing   a   random  operator  tensor,  $\eta_{ij}$,  one  can  write
Eq.~(\ref{eqno_S2}) as,
%
\begin{equation}
\label{eqno_S13}
\frac{d\sigma_{ij}(x,t)}{dt}
=
\eta_{ij}(\sigma_{ij},D_{ijkl}, \epsilon_{kl};x,t)
\end{equation}
%
with an initial condition,
%                                                                                                                                    
\begin{equation}
\label{eqno_S14}
\sigma_{ij}(x,0)=\sigma_{{ij}_0}
\end{equation}
%
%                                                                                                 
In  Eq.~(\ref{eqno_S13}), the stress tensor $\sigma_{ij}$ can be considered to represent a point
in  a $ 9 $-dimensional stress ($\sigma$)-space and hence Eq.~(\ref{eqno_S13}) determines
the  velocity  for the point in this space. It is possible to visualize this, by
imagining  an  initial  point  in  stress  space, given by its initial condition
$\sigma_{{ij}_0}$,  with  a  trajectory    starting  out  that describes  the 
corresponding solution of the
non-linear    stochastic    ordinary    differential   equation   (ODE)   system
(Eq.~(\ref{eqno_S13})).  Let us now consider  a  cloud  of  initial  points,
described      by    a   density  $\rho(\sigma_{ij},0)$  in  the $\sigma$-space,
and   with   movements      of      these    points    dictated    by  
Eq.~(\ref{eqno_S13}), the phase
density  $\rho$ of $\sigma_{ij}(x,t)$ varies in time according to a continuity
equation   which   expresses  the  conservation  of  all  these  points  in  the
$\sigma$-space.

Expressing  this  continuity  equation in mathematical terms, one obtains Kubo's
stochastic Liouville equation (Kubo \cite{Kubo:1963}):
%                                                                                                                 
\begin{equation}
\label{eqno_S15}
\displaystyle \frac{\partial \rho (\sigma_{ij}(x,t),t)}{\partial t}
=
-\displaystyle 
\frac{\partial}{\partial \sigma_{mn}} \eta_{mn} 
\left(\sigma_{mn}(x,t), D_{mnpq}(x),\epsilon_{pq}(x,t)\right)
\rho(\sigma_{ij}(x,t),t)
\end{equation}
%                                                                                                                                                             
with initial condition,
%                                                                                                                                                             
\begin{equation}
\label{eqno_S16}
\rho(\sigma_{ij},0)=\delta(\sigma_{ij}-\sigma_{{ij}_0})
\end{equation}





%
\noindent where   $\delta(\cdot)$ is the  Dirac  delta  function.  Eq.~(\ref{eqno_S16})  is  the
probabilistic   restatement   in   the  $\sigma$-phase  space  of  the  original
deterministic  initial  condition (Eq.~(\ref{eqno_S14})). One can then apply Van Kampen's
Lemma (Van Kampen \cite{VanKampen:1976}) to obtain:
%                                                                                                                                                             
\begin{equation}
\label{eqno_S17}
<\rho(\sigma_{ij},t)>=P(\sigma_{ij},t)
\end{equation}
%                                                                                                                                                             
where,    the  symbol  $<\cdot>$  denotes  the  expectation  operation,  and
$P(\sigma_{ij},t)$   denotes  evolutionary  probability  density  of  the  state
variable tensor $\sigma_{ij}$ from the constitutive equations.
                                                                                                                                                             
                                                                                                                                                             
Therefore,  in  order  to  obtain  the multivariate probability density function
(PDF), $P(\sigma_{ij},t)$, of  the state variable tensor $\sigma_{ij}$, it is
necessary to obtain the deterministic partial differential equation (PDE) of the
$\sigma$-space  mean  phase  density $ < \rho (\sigma_{ij},t)> $ from the linear
stochastic   PDE  system  (Eqs.~(\ref{eqno_S15})  and  ~(\ref{eqno_S16})).  This
necessitates the derivation of the ensemble average form of Eq.~(\ref{eqno_S15})
for      $      <\rho   (\sigma_{ij},t)>   $.   The  ensemble  average  form  of
Eq.~(\ref{eqno_S15}) was derived by Kavvas \cite{Kavvas:2003} as follows:
% 
\begin{eqnarray}
\nonumber
\lefteqn{\frac{\partial \left<\rho (\sigma_{ij}(x_t,t), t)\right>}{\partial t} = -\displaystyle \frac{\partial}{\partial \sigma_{mn}}
\left[ \left\{\left< \vphantom{\int_{0}^{t}} \eta_{mn}(\sigma_{mn}(x_t,t), D_{mnrs}(x_t), \epsilon_{rs}(x_t,t))\right> \right. \right.} 
\\
\nonumber
&-& \left. \left. \int_{0}^{t} d\tau Cov_0\left[\vphantom{\frac{\partial}{\partial}}\eta_{mn}(\sigma_{mn}(x_t,t),D_{mnrs}(x_t), 
\epsilon_{rs}(x_t,t)); \right. \right. \right. \\
\nonumber
& & \left. \left. \left. \displaystyle \frac{\partial \eta_{ab} (\sigma_{ab}(x_{t-\tau}, t-\tau), D_{abcd}(x_{t-\tau}), 
\epsilon_{cd}(x_{t-\tau}, t-\tau)}{\partial \sigma_{ab}}\right]\right\}\left<\rho (\sigma_{ij}(x_t,t), t)\right>\right] \\
\nonumber
&+& \displaystyle \frac{\partial}{\partial \sigma_{mn}}\left[ \left\{\int_{0}^{t} d\tau Cov_0\left[\eta_{mn}(\sigma_{mn}(x_t,t), 
D_{mnrs}(x_t), \epsilon_{rs}(x_t,t)); \right. \right. \right. \\
\label{eqno_S18}
& & \left. \left. \left. \eta_{ab} (\sigma_{ab}(x_{t-\tau}, t-\tau), D_{abcd}(x_{t-\tau}), \epsilon_{cd}(x_{t-\tau}, t-\tau)) \right]
\vphantom{\int_{0}^{t}} \right\}  \displaystyle \frac{\partial \left<\rho (\sigma_{ij}(x_t,t), t)\right>}{\partial \sigma_{ab}} \right ]
\end{eqnarray}
%                                                                                                                                                             
to  exact  second  order  (to  the order of the covariance time of $ \eta $). In
Eq.~(\ref{eqno_S18}),  $  Cov_0[\cdot] $ is the time ordered covariance function
defined by,
%                                                                                                                                                             
\begin{equation}
\label{eqno_S19}
Cov_0\left[\eta_{mn}(x,t_1), \eta_{ab}(x,t_2)\right]
=
\left<\eta_{mn}(x,t_1) \eta_{ab}(x,t_2)\right>-\left<\eta_{mn}(x,t_1)\right> 
\cdot \left<\eta_{ab}(x,t_2)\right>
\end{equation}
%
By  combining  Eq.~(\ref{eqno_S18})  with Eq.~(\ref{eqno_S17}) and rearranging
the        terms   yields   the   following   Eulerian-Lagrangian   form   of   
Fokker--Planck--Kolmogorov equation (FPKE) (cf. Kavvas \cite{Kavvas:2003}):
%                                                                                                                                                             
\begin{eqnarray}
\nonumber
\lefteqn{\displaystyle \frac{\partial P(\sigma_{ij}(x_t,t), t)}{\partial t} = \displaystyle \frac{\partial}{\partial \sigma_{mn}} 
\left[ \left\{\left< \vphantom{\int_{0}^{t}} \eta_{mn}(\sigma_{mn}(x_t,t), D_{mnrs}(x_t), \epsilon_{rs}(x_t,t))\right> \right. \right.} \\
\nonumber
&+& \left. \left. \int_{0}^{t} d\tau Cov_0 \left[\displaystyle \frac{\partial \eta_{mn}(\sigma_{mn}(x_t,t), D_{mnrs}(x_t), 
\epsilon_{rs}(x_t,t))} {\partial \sigma_{ab}}; \right.  \right. \right. \\
\nonumber
& & \left. \left. \left. \eta_{ab} (\sigma_{ab}(x_{t-\tau}, t-\tau), D_{abcd}(x_{t-\tau}), \epsilon_{cd}(x_{t-\tau}, t-\tau) 
\vphantom{\int_{0}^{t}} \right] \right \} P(\sigma_{ij}(x_t,t),t) \right] \\
\nonumber
&+& \displaystyle \frac{\partial^2}{\partial \sigma_{mn} \partial \sigma_{ab}} \left[ \left\{ \int_{0}^{t} d\tau Cov_0 \left[ 
\vphantom{\int_{0}^{t}} \eta_{mn}(\sigma_{mn}(x_t,t), D_{mnrs}(x_t), \epsilon_{rs}(x_t,t)); \right. \right. \right. \\
\label{eqno_S20}
& & \left. \left. \left. \eta_{ab} (\sigma_{ab}(x_{t-\tau}, t-\tau), D_{abcd}(x_{t-\tau}), \epsilon_{cd}(x_{t-\tau}, t-\tau)) 
\vphantom{\int_{0}^{t}} \right] \vphantom{\int_{0}^{t}} \right\} P(\sigma_{ij}(x_t,t),t) \right]
\end{eqnarray}
%
to  exact  second  order.  The  solution  of  this  deterministic  linear  FPKE
(Eq.~(\ref{eqno_S20}))  in  terms  of  $  P(\sigma_{ij},t)  $  under appropriate
initial  and boundary conditions will yield the PDF of the state variable tensor
$  \sigma_{ij}  $.  It  is  important  to  note that while the original equation
(Eq.~(\ref{eqno_S2}))  is  non-linear, the FPKE (Eq.~(\ref{eqno_S20})) is linear
in  terms  of  its  unknown, the probability density $ P(\sigma_{ij},t) $ of the
state  variable  tensor  $  \sigma_{ij}  $.  This  linearity,  in turn, provides
significant  advantages  in  probabilistic  solution  of  the  constitutive rate
equation.


One  should  also  note  that  Eq.~(\ref{eqno_S20}) is a mixed Eulerian-Lagrangian
equation, since, while the real space location $x_t$ at time $t$ is known, the
location  $x_{t-\tau}$  is  an  unknown. If one assumes, for example, 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_S20a}
d\epsilon=\dot \epsilon \tau = \frac{x_t-x_{t-\tau}}{x_t}
\end{equation}
%
or, by rearranging
%
\begin{equation}
\label{eqno_S20b}
x_{t-\tau}=(1-\dot \epsilon \tau)x_t
\end{equation}

%Once  the probability density function is obtained it can be used to obtain the
%mean of state variable (H) by usual expectation operation as follows,

%\begin{equation}
%\label{eqno_S20c}
%<H(t)>=\int H(t) P(H(t)) dH(t)
%\end{equation}

%\noindent  Another  possible way to obtain the mean of state variable is to use
%the equivalence between FPKE and It\^o stochastic differential equation.
%Kavvas  \cite{Kavvas:2003}  obtained  the  equaivalent  operator  vector  It\^o
%stochastic differential equation as:


Another  interesting  aspect  to  note, possibly of theoretical interest, is the
It\^o  stochastic  differential  equation  (SDE)  form  of the constitutive rate
equation    with   random   material   properties   and   random   strain   rate
(Eq.~(\ref{eqno_S13})).  Utilizing  the  connection  between  FPKE and It\^o SDE
(cf. Gardiner \cite{book:Gardiner}), one can write Eq.~(\ref{eqno_S20}) as:
%
\begin{eqnarray}
\nonumber
\lefteqn{d\sigma_{ij} (x,t) =} \\
\nonumber
&& \left\{ 
\left< \eta_{ij} (\sigma_{ij}(x_t,t), D_{ijkl}(x_t), \epsilon_{kl}(x_t,t)) \vphantom{\int_0^t} \right> 
+ \int_0^t d \tau Cov_0 \left[\displaystyle 
\frac{\partial \eta_{ij} (\sigma_{ij}(x_t,t), D_{ijkl}(x_t), \epsilon_{kl}(x_t,t))}{\partial \sigma_{ab}}; \right. \right. \\
\label{eqno_S20d}
&& \left. \left. \eta_{ab}(\sigma_{ab}(x_{t-\tau},t-\tau), D_{abcd}(x_{t-\tau}), \epsilon_{cd} (x_{t-\tau},t-\tau)) 
\vphantom{\displaystyle \frac{\partial \eta_j (H(x_t,t), A(x_t,t), \Phi(x_t,t))}{\partial \sigma_M}}\right] \right\} dt + 
C_{ijm}(\sigma_K,t) dW_m(t)
\end{eqnarray}
%
where,
%
\begin{eqnarray}
\nonumber
\lefteqn{C_{ijm}(\sigma_{ij},t) C_{abm}(\sigma_{ab},t)  =  2 \int_0^t d \tau Cov_0 \left[\vphantom{\int_0^t} \eta_{ij} (\sigma_{ij}
(x_t,t), D_{ijkl}(x_t), \epsilon_{kl}(x_t,t)); \right.} \\
\label{eqno_S20e}
& & \left. \eta_{ab}(\sigma_{ab}(x_{t-\tau},t-\tau), D_{abcd}(x_{t-\tau}), \epsilon_{cd} (x_{t-\tau},t-\tau)) \vphantom{\int_0^t} \right]
\end{eqnarray}
%
and,  $  dW_i(t) $ is an increment of the vector Wiener process $W_i$ having the
following properties:
%
\begin{equation}
\label{eqno_S20f1}
<dW_i(t)> =0  
\end{equation}
%
and, 
%
\begin{eqnarray}
\nonumber
<dW_i(t)dW_j(t)> &=& \delta_{ij}dt \  \ \mbox{for $\tau =t$} \\
\label{eqno_S20f}
                 &=& 0 \ \ \mbox{for $ \tau \not= t$}
\end{eqnarray}
%
It  is  also  interesting  to  note that all the stochasticity of the original
system  of equations (Eq.~(\ref{eqno_S13})) are lumped together in the last term
(Wiener increment term) of the r.h.s of Eq.~(\ref{eqno_S20d}).



In  theory,  one  could  use  the  It\^o  form of the constitutive rate equation
(Eq.~(\ref{eqno_S20d}))  in  obtaining the mean of stress tensor ($\sigma_{ij}$)
as (cf. Kavvas \cite{Kavvas:2003}):
% 
\begin{eqnarray}
\nonumber
\lefteqn{\displaystyle \frac{<d\sigma_{ij} (x,t)>}{dt} =} \\
\nonumber
&&  \left< \vphantom{\int_0^t} \eta_{ij} (\sigma_{ij}(x_t,t), D_{ijkl}(x_t), \epsilon_{kl}(x_t,t)) \right>  
+  \int_0^t d \tau Cov_0 \left[\displaystyle \frac{\partial \eta_{ij}(\sigma_{ij}(x_t,t), D_{ijkl}(x_t), \epsilon_{kl}(x_t,t))}
{\partial \sigma_{ab}}; \right. \\ 
\label{eqno_S20g}
& & \left. \eta_{ab}(\sigma_{ab}(x_{t-\tau},t-\tau), D_{abcd}(x_{t-\tau}), \epsilon_{cd}(x_{t-\tau},t-\tau)) \vphantom{\int_0^t} \right]
\end{eqnarray} 
%
but the difficulty is that the stress tensor appearing within $\eta_{ij}(\cdot)$
in the covariance term on the r.h.s of Eq.~(\ref{eqno_S20g}) is random and needs
to  be  treated  accordingly.  One  possible way is perturbation with respect to
mean,  however  the  ``closure problem'' (cf. Kavvas \cite{Kavvas:2003}) will make this
type  of  solution  problematic.  Hence  in this study the mean of stress tensor
($\sigma_{ij}$)  will  be  computed  by  standard operation on the PDF of stress
tensor    ($\sigma_{ij}$),    which    will   be   obtained   by   solving   the
Fokker--Planck--Kolmogorov PDE (Eq.~(\ref{eqno_S20})).


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

\section{APPLICATION  TO PROBABILISTIC BEHAVIOR OF 1-D
CAM--CLAY ELASTIC-PLASTIC CONSTITUTIVE EQUATION}
\label{CamClay}

In  this  section  the  applicability  of  the  presented  general expression in
obtaining  PDF  and  subsequently  mean  and variance behaviors of stress will be
demonstrated    for    1-D,  elastic--plastic,  Cam--Clay    constitutive  rate 
equation.  
%Probabilistic solution of stochastic   constitutive   rate  equation 
%could be useful in carrying out
%probabilistic  integrations  at  the  constitutive  level  for stochastic finite
%element  method.  Also  it  could  be used in performing sensitivity analysis of
%constitutive  rate  equation.  
The    general    incremental    form  of  constitutive  equation  with single
scalar internal variable can be written as:
%
\begin{equation}
\label{eqno_CC1}
\frac{d\sigma_{ij}(t)}{dt}= D^{ep}_{ijkl}(t) \frac{d\epsilon_{kl}(t)}{dt}
\end{equation}
%
where  the  stiffness  tensor  $D_{ijkl}^{ep}$  is given as:
%
\begin{eqnarray}
D_{ijkl}^{ep} = \left\{\begin{array}{ll}
%
           D^{el}_{ijkl}
%
%
        \;\;\; & \mbox{when~~~} f <  0 \vee (f = 0 \wedge df <  0) \\
%
           \\
%
           D^{el}_{ijkl}
           -
           \displaystyle \displaystyle \frac{D^{el}_{ijmn}
           \displaystyle \frac{\partial f}{\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 f}{\partial \sigma_{tu}}
           -
           \displaystyle \frac{\partial f}{\partial p_0}\bar p_0}
%
        \;\;\; & \mbox{when~~~} f = 0 \vee df = 0
%
           \end{array} \right.
\label{el_or_elpl}
\end{eqnarray}
%
where $ D^{el}_{ijkl} $ is the elastic stiffness tensor and
$D^{ep}_{ijkl} $ is the elastic--plastic stiffness tensor. The internal variable
$p_0$  in  this  case  depends  on  plastic volumetric strain and has a physical
meaning  (for  Cam--Clay)  as  the  maximum  hydrostatic  stress  the  soil  has
experienced   during previous loading history). It was assumed that the material
obeys associative flow rule (as is it is usually assumed for Cam--Clay material
model).  Hence,  both the  yield  and  the plastic  potential  functions  $f$  are given by
(Cam--Clay material model):
%
\begin{equation}
\label{eqno_CC2}
f=\hat f(p,q) = p^2 - p_0p + \frac{q^2}{M^2} = 0
\end{equation}
%
Material parameter $ M $ represents the slope of
the    critical state line in the $p-q$ space. The stress invariants $p$ and $q$
are given as:
%
\begin{equation}
\label{eqno_CC3}
p 
= 
- \frac{\sigma_{kk}}{3}
=
-\frac{\sigma_{11}+\sigma_{22}+\sigma_{33}}{3} 
\end{equation}
%
\begin{equation}
\label{eqno_CC$}
q
=
\sqrt{\frac{3}{2}{s}_{ij}{s}_{ij}} 
=
\frac{1}{\sqrt{2}}\left[ 
(\sigma_{11}-\sigma_{22})^2 + 
(\sigma_{22}-\sigma_{33})^2 + 
(\sigma_{33}-\sigma_{11})^2 + 
6 \sigma_{12}^2 + 
6 \sigma_{23}^2 + 
6 \sigma_{31}^2 \right]^{1/2} 
\end{equation}
%
%
%\begin{samepage}
%\begin{eqnarray}
%p= -\frac{1}{3} I_{1}\;\;\;\;
%q= \sqrt{3 J_{2D}} \;\;\;\;
%\cos 3\theta= \frac{3\sqrt{3}}{2}\frac{J_{3D}}{\sqrt({J_{2D})^3}}
%\label{invdef1}
%\end{eqnarray}
%\end{samepage}
%
%
%\begin{samepage}
%\begin{eqnarray}
%I_1 = \sigma_{kk}\;\;\;\;
%J_{2D} = \frac{1}{2}{s}_{ij}{s}_{ij} \;\;\;\;
%J_{3D} = \frac{1}{3} {s}_{ij}{s}_{jk}{s}_{ki}\;\;\;\;
%{s}_{ij} = \sigma_{ij} - \frac{1}{3}\sigma_{kk}\delta_{ij}
%\label{invdef2}
%\end{eqnarray}
%\end{
%


%where,    $I_1$  is  the  first  invariant of the stress tensor and $J_2$ is the
%second invariant of the deviatoric stress tensor ${s}_{ij} =
%\sigma_{ij} - {1}/{3}\sigma_{kk}\delta_{ij}$. One
%may  also note that in terms of conventional stress invariants $ I_1 $ and $ J_2
%$, Eq.~(\ref{eqno_CC2}) can be written as:
%
%\begin{equation}
%\label{eqno_CC5}
%f=\hat f(I_1, \sqrt{J_2}) = I_1^2 - I_1 I_{01} + 27 \frac{J_2}{M^2} = 0
%\end{equation}
%                                                                                                                                                             
%\noindent  where  $  I_{01}  $  is analogous to $ p_0 $ of the original cam clay
%formulation and is the value of $ I_1 $ at the intersection of the yield surface
%with  the $ I_1 $ axis in the $I_1-\sqrt{J_2}$ space. 

One can work on components of the Eq.~\ref{el_or_elpl} to obtain:
%                                                                                                                                                             
\begin{align}
\nonumber
A_{kl}
&=
\frac{\partial f}{\partial \sigma_{pq}}D_{pqkl}^{el} 
= \nonumber \\
& \frac{\partial f}{\partial p}\left[2G\left(\frac{\partial p}{\partial
\sigma_{11}}\delta_{1l}\delta_{1k}
+\frac{\partial p}{\partial\sigma_{22}}\delta_{2l}\delta_{2k}
+\frac{\partial p}{\partial
\sigma_{33}}\delta_{3l}\delta_{3k}\right)+\left(K-\frac{2}{3}G
\right)\frac{\partial p}{\partial \sigma_{cd}}\delta_{cd}\delta_{kl}\right]
\nonumber \\
\label{eqno_A6}
&+\frac{\partial f}{\partial q}\left[2G\frac{\partial q}{\partial \sigma_{ij}}\delta_{ik}\delta_{jl}+\left(K-\frac{2}{3}G \right)
\frac{\partial q}{\partial \sigma_{ab}}\delta_{ab}\delta_{kl}\right]
%\label{eqno_A3A}
%&= -3K\alpha \delta_{kl} + \displaystyle \frac{G}{\sqrt{J_2}}S_{kl} 
\end{align}
%                                                                                                                                                             
and
%
\begin{align}
\nonumber
B
&=
\frac{\partial f}{\partial \sigma_{rs}}D_{rstu}^{el}\frac{\partial f}{\partial \sigma_{tu}} 
= \nonumber \\
& 
\left(\frac{\partial f}{\partial p}\right)^2\left[2G\left\{\left(\frac{\partial p}{\partial
\sigma_{11}}\right)^2+\left(\frac{\partial p}{\partial\sigma_{22}}\right)^2+\left(\frac{\partial p}{\partial
\sigma_{33}}\right)^2\right\} + \left(K-\frac{2}{3}G \right)\left\{\frac{\partial p}{\partial \sigma_{ij}}\delta_{ij}\right\}^2\right]
\nonumber \\
\label{eqno_A7}
&+
\left(\frac{\partial f}{\partial q}\right)^2 \left[2G\frac{\partial q}{\partial \sigma_{ij}}\frac{\partial
q}{\partial \sigma_{ij}} +\left(K-\frac{2}{3}G \right)\left\{\frac{\partial q}{\partial
\sigma_{ij}}\delta_{ij}\right\}^2 \right]
%\label{eqno_A4A}
%&= 9K\alpha ^2 + G
\end{align}
%                                                                                                                                                             
where,  $K$  and $G$ are the elastic bulk and shear modulus, respectively.
For the  Cam--Clay formulation the material parameter $ p_0 $ represents an internal
variable  and its evolution depends on the plastic volumetric strain. Hence, the
plastic modulus $K_P$ becomes:
%                                                                                                                                                             
\begin{equation}
\label{eqno_A8}
K_P
=
- \frac{\partial f}{\partial p_0}\bar p_0
= 
(1+e_0)\frac{p_0}{\lambda-\kappa}\frac{\partial f}{\partial p}
\end{equation}
%
One  may  note  that  the  differentiations  appearing  in Eqs.~(\ref{eqno_A6}),
~(\ref{eqno_A7}),  and  ~(\ref{eqno_A8}),  involve  stochastic process and hence
cannot    be    carried    out    in   a   deterministic   sense.   Substituting
Eqs.~(\ref{eqno_A6}),     ~(\ref{eqno_A7}),     and     ~(\ref{eqno_A8})    into
~(\ref{el_or_elpl})  one  may  obtain  the  random operator tensor ($\eta_{ij}$)
specialized to Cam--Clay model as:
%                                                                                                                                                             
\begin{eqnarray}
\eta_{ij} = \left\{\begin{array}{ll}
%
          \left[2G\delta_{ik}\delta_{jl}+\left(K-\displaystyle \frac{2}{3}G\right) \delta_{ij}\delta_{kl}\right] \displaystyle 
          \frac{d\epsilon_{kl}}{dt} 
%
%
         & \mbox{when~~~} f <  0 \vee (f = 0 \wedge df <  0) \\
%
           \\
%
           \left[2G\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}\right] \displaystyle \frac{d\epsilon_{kl}}{dt}
%
        \;\;\; & \mbox{when~~~} f = 0 \vee df = 0
%
           \end{array} \right.
\label{els_or_elpls}
\end{eqnarray}
%
where  $A_{ij}$,  $B$, and  $K_P$ are defined by Eqs.~(\ref{eqno_A6}),
~(\ref{eqno_A7}), and ~(\ref{eqno_A8}) respectively.


Previously  developed  random operator tensor $\eta_{ij}$ for Cam--Clay material
model    will    be   now  specialized  to  1-D shear behavior, resulting in the
following equation for the PDF of shear stress:
%solving for relating PDF of shear
%stress $\sigma_{12}$  and random shear strain $\epsilon_{12}$:
% one can specialize
%the   3-D   general   local-average   form   master-key  relation  presented  in
%Section~\ref{MasterKey}  to cam clay point-location scale 1-D shear constitutive
%rate equation as:
%
% 
\begin{eqnarray}
\nonumber
& &\displaystyle \frac{\partial P(\sigma_{12}(t), t)}{\partial t}=
- \displaystyle \frac{\partial}{\partial \sigma_{12}} \left[ \left\{\vphantom{\int_{0}^{t}}
\left<\eta (\sigma_{12}(t), D_{1212}, \epsilon_{12}(t))\right> \right. \right. \\
\nonumber
&+& \left. \left. \int_{0}^{t} d\tau Cov_0
\left[\displaystyle \frac{\partial \eta (\sigma_{12}(t), D_{1212}, \epsilon_{12}(t))}{\partial \sigma_{12}};
\eta (\sigma_{12}(t-\tau), D_{1212}, \epsilon_{12}(t-\tau) \right)] \right \} P(\sigma_{12}(t),t) \right] \\
\nonumber
&+& \displaystyle \frac{\partial^2}{\partial \sigma_{12}^2} \left[ \left\{ \int_{0}^{t} d\tau Cov_0 \left[ \vphantom{\int_{0}^{t}} 
\eta (\sigma_{12}(t), D_{1212}, \epsilon_{12}(t)); \eta \sigma_{12}(t-\tau), D_{1212}, \epsilon_{12}(t-\tau) \vphantom{\int_{0}^{t}} 
\right] \vphantom{\int_{0}^{t}} \right\} P(\sigma_{12}(t),t) \right] \\
\label{eqno_A9}
\end{eqnarray}
%
%
One  may  note that while solving Eq.~(\ref{eqno_A9}) for probability
density ($P(\sigma_{12})$) of  each fixed value of $\sigma_{12}$ in its domain,
the  $\sigma_{12}$  terms  appearing within $\eta (\cdot)$ are deterministic and
hence  the  differentiations involving $\sigma_{12}$ can be carried out in a usual
sense. In Eq.~(\ref{eqno_A9}) the random operator function $\eta$ is different
for    pre--yield  (elastic)  and  post--yield (elastic--plastic) responses. For
pre--yield,  linear  elastic  response, one can obtain the random operator function
$\eta$ in the form:
%
\begin{equation}
\label{eqno_A10}
\eta = 2 G \frac{d\epsilon_{12}}{dt}
\end{equation}
%
while,  for post-yield, non-linear, elastic--plastic response $\eta$ obtains the
form:

\begin{equation}
\label{eqno_A11}
\eta = \left[2G - \displaystyle \frac{\left(36 \displaystyle \frac{G^2}{M^4} \right) \sigma_{12}^2}
{\displaystyle \frac{(1+e_0)p\left(p-\displaystyle \frac{3\sigma_{12}^2}{pM^2}\right)^2}{\kappa} + \left(18 \displaystyle \frac{G}{M^4}\right) \sigma_{12}^2
+ \displaystyle \frac{1+e_0}{\lambda-\kappa} p p_0 \left(p-\displaystyle \frac{3\sigma_{12}^2}{pM^2}\right)} \right] \displaystyle \frac{d\epsilon_{12}}{dt}
\end{equation}                                                                                                                       


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
\section{RESULTS AND DISCUSSIONS} 
\label{Results}


In    this    section   the  Fokker--Planck--Kolmogorov (FPK) equation is solved
using  appropriate initial and boundary conditions. The FPK equation is  
specialized to 1-D (shear) Cam--Clay constitutive relationship 
with  different  degrees of uncertainties (Eqs.~(\ref{eqno_A9}),
~(\ref{eqno_A10}),  and  ~(\ref{eqno_A11})) and the solution will be in terms 
of the  evolution  of  Probability Density Function (PDF) of shear stress
($\sigma_{12}$)  with  time  (and/or strain), given random material property(ies) and
deterministic  shear  strain.  It  should  be  noted that the strain can be
random    as   well, and the equations can be solved in that case too, but  for 
the  sake  of  simplicity, strain (or time) is preset to deterministic value. In
other words, the solution is driven by deterministic strain.

%Five different cases 
%were considered - (a) Low overconsolidation ratio (OCR) response with random normally distributed shear modulus ($G$); (b) Low OCR 
%response with random normally distributed shear modulus ($G$) and random normally distributed slope of critical state line ($M$); 
%(c) Low OCR response with random normally distributed shear modulus ($G$), random normally distributed slope of critical state 
%line ($M$), and random normally distributed overconsolidation pressure ($p_0$); (d) Low OCR response with random normally distributed 
%shear modulus ($G$) and random normally distributed overconsolidation pressure ($p_0$); (e) High OCR response with random normally 
%distributed shear modulus ($G$) and random normally distributed slope of critical state line ($M$).    

\subsection{Initial and Boundary Conditions}

The   Fokker--Planck--Kolmogorov   PDE   specialized   to  1--D Cam--Clay  shear
constitutive  rate equation (Eq.~(\ref{eqno_A9}) can be written in the following
simplified form:
%
\begin{eqnarray}
\nonumber
\displaystyle \frac{\partial P(\sigma_{12}, t)}{\partial t}=
&=& -\displaystyle \frac{\partial}{\partial \sigma_{12}} \left\{P(\sigma_{12},t) N_{(1)} \right\} + \displaystyle \frac{\partial^2}{\partial 
\sigma_{12}^2} \left\{P(\sigma_{12},t)N_{(2)}\right\} \\
\label{eqno_R1a}
&=& -\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_R1b}
&=& -\displaystyle \frac{\partial \zeta}{\partial \sigma_{12}} 
\end{eqnarray}
%
%
where,      $N_{(1)}$      and      $N_{(2)}$      are   coefficients   of   the
Fokker--Planck--Kolmogorov   PDE   (Eq.~(\ref{eqno_A9}))   and   represent   the
expressions within the curly braces of the first and second term respectively on
the  right hand side  of  Eq.~(\ref{eqno_A9}). They are called the advection and diffusion
coefficients,  respectively as the form of Eq.~(\ref{eqno_A9}) closely resembles
advection-diffusion  equation  (cf. Gardiner \cite{book:Gardiner}). The symbol $\zeta$ in
Eq.~(\ref{eqno_R1b})  can  be  considered  to  be  the  probability current (cf.
Gardiner \cite{book:Gardiner}),  since  Eq.~(\ref{eqno_R1b}) is a continuity equation and
the  state  variable  of  the  equation  (Eq.~(\ref{eqno_R1b}))  is  probability
density.


The  initial  condition  for    pre-yield   (elastic) solution for PDF of shear 
stress,  will  be  deterministic  i.e. at time $t=0$ all the probability mass is
concentrated  at  $\sigma_{12}=0$  or  at  some fixed value of $\sigma_{12}$ (if
there was some initial deterministic stress to begin with). Mathematically, this
is represented using Dirac delta function $\delta(\cdot)$,
%
\begin{equation}
\label{eqno_R1}
P(\sigma_{12},0)=\delta(\sigma_{12})
\end{equation}
%

On    the  other  hand,  for  solving for post-yield  (elastic--plastic)  PDF of
shear  stress,   the  initial  condition  will  be  represented  by  the  PDF of
shear stress corresponding to the pre--yield solution (PDF of shear stress at yield).


The initial probability mass, dictated by Eq.~(\ref{eqno_A9}), will advect and
diffuse  into  the  domain  of  the  system  throughout  the  evolution  of  the
simulation.    As for the boundary condition, a  reflecting barrier boundary  is
chosen since conservation of the probability mass within the system
is  required (i.e. no leaking is allowed at the boundaries).
Mathematically this condition is expressed as (cf. Gardiner \cite{book:Gardiner})
%
\begin{equation}
\label{eqno_R2}
\zeta(\sigma_{12},t)|_{AtBoundaries}=0
\end{equation}
%


By making the assumption that the  material parameters are normally distributed,
theoretically  the  domain  of  stress  could be from -$\infty$ to +$\infty$ and
hence one can write the boundary conditions as:
%
\begin{equation}
\label{eqno_R3}
\zeta(-\infty,t)=\zeta(\infty,t)=0
\end{equation}
%
It   should also be noted that any other probabilistic distribution for material
parameters  can  be  handled.  For example, by assuming that material parameters
obey  exponential distribution, theoretically the domain of stress could be from
0 to +$\infty$ and hence one can write the boundary conditions as:
%
\begin{equation}
\label{eqno_R3b}
\zeta(0,t)=\zeta(\infty,t)=0
\end{equation}

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


\subsection{Numerical Scheme}
\label{NumericalScheme}

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



The    FPKE    specialized   to 1--D  Cam--Clay   constitutive    equation
(Eq.~(\ref{eqno_A9})  with  Eqs.~(\ref{eqno_A10})  or  ~(\ref{eqno_A11}))  with
initial  and  boundary  conditions  (Eqs.~(\ref{eqno_R1})  and (\ref{eqno_R3})
respectively)    were   solved   numerically   by   \textit{method   of   lines}
(using  Mathematica  \cite{book:Mathematica}).
The Fokker--Planck--Kolmogorov PDE is first semi-discretized  into a set of simultaneous
ODE systems and then solved by using an ODE solver.


To        this        end,    differentiating  by parts the  right hand side  of
Eq.~(\ref{eqno_R1b}) yields:
%
\begin{eqnarray}
\nonumber
\frac{\partial P}{\partial t} &=&  -N_{(1)} \frac{\partial P}{\partial \sigma_{12}} - P \frac{\partial N_{(1)}}{\partial \sigma_{12}} + N_{(2)} \frac{\partial^2 P}{\partial \sigma_{12}^2}
+ 2 \frac{\partial P}{\partial \sigma_{12}}\frac{\partial N_{(2)}}{\partial \sigma_{12}} + P \frac{\partial^2 N_{(2)}}{\partial \sigma_{12}^2} \\
\label{eqno_R4}
&=& P\left(\frac{\partial^2 N_{(2)}}{\partial \sigma_{12}} - \frac{\partial N_{(1)}}{\partial \sigma_{12}}\right) + \frac{\partial P}{\partial
\sigma_{12}} \left( 2\frac{\partial N_{(2)}}{\partial \sigma_{12}}-N_{(1)}\right) + \frac{\partial^2 P}{\partial \sigma_{12}^2} N_{(2)} 
\end{eqnarray}
%
then,     using     central    difference    discretization,    as    shown    in
Fig.~\ref{figure:SemiDiscretization}
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.8\textwidth]{StressDiscretization}
\vspace*{-0.8cm}
\caption{Spatial Discretization of the Fokker--Planck--Kolmogorov PDE}
\label{figure:SemiDiscretization}
\end{center}
\end{figure}
%
one  can  write semi-discretized form of     Eq.~(\ref{eqno_R4})    at    any   
intermediate node $i$ to as:
%
\begin{eqnarray}
\nonumber
\frac{\partial P^{(i)}}{\partial t} = && P^{(i-1)} \left[ \frac{N^{(i)}_{(1)}}{2 \Delta \sigma_{12}} + \frac{N^{(i)}_{(2)}}{\Delta \sigma_{12}^2} - \frac{1}
{\Delta \sigma_{12}} \frac{\partial N^{(i)}_{(2)}}{\partial \sigma_{12}} \right] - P^{(i)} \left[ \frac{\partial N^{(i)}_{(1)}}{\partial \sigma_{12}} 
+ 2 \frac{N^{(i)}_{(2)}}{\Delta \sigma_{12}^2} - \frac{\partial^2 N^{(i)}_{(2)}}{\partial \sigma_{12}^2} \right] \\ 
\label{eqno_R5}
&& + P^{(i+1)} \left [ -\frac{N^{(i)}_{(1)}}{2
\Delta \sigma_{12}}  + \frac{N^{(i)}_{(2)}}{\Delta \sigma_{12}^2} +  \frac{1} {\Delta \sigma_{12}} \frac{\partial N^{(i)}_{(2)}}{\partial \sigma_{12}} 
\right] 
\end{eqnarray}
%
which    forms    an    initial value problem in the time dimension. It is worth
noting  that  central finite difference technique might not be the most efficient
in  solving  n-dimensional  Fokker--Planck--Kolmogorov  PDE. Work is underway to
improve  the  efficiency of solution method through adaptivity and reduced order
modeling.

%Using forward difference technique, one can introduce the boundary 
%condition at the left boundary (node 1 in Fig.~\ref{figure:SemiDiscretization}) as:

%\begin{equation}
%\label{eqno_R6}
%N^{(1)}_{(1)}P^{(1)} - N^{(1)}_{(2)} \displaystyle \frac{P^{(2)}-P^{(1)}}{\Delta \sigma_{12}} - P^{(1)}\displaystyle \frac{\partial 
%N^{(1)}_{(2)}}{\partial \sigma_{12}}=0
%\end{equation}

%\noindent Or rearranging,

%\begin{equation}
%\label{eqno_R7}
%P^{(1)}= P^{(2)} \left[\displaystyle \frac{\displaystyle \frac{N^{(1)}_{(2)}}{\Delta \sigma_{12}}}{N^{(1)}_{(1)} + \displaystyle \frac{
%N^{(1)}_{(2)}}{\Delta \sigma_{12}} - \displaystyle \frac{\partial N^{(1)}_{(2)}}{\partial \sigma_{12}}} \right]
%\end{equation}

%\noindent Similarly, using backward difference technique, one can introduce the boundary condition at the right boundary (node $n$ in 
%Fig.~\ref{figure:SemiDiscretization}) as:

%\begin{equation}
%\label{eqno_R8}
%P^{(n)}= P^{(n-1)} \left[\displaystyle \frac{\displaystyle \frac{N^{(n)}_{(2)}}{\Delta \sigma_{12}}}{-N^{(n)}_{(1)} + \displaystyle \frac{
%N^{(n)}_{(2)}}{\Delta \sigma_{12}} - \displaystyle \frac{\partial N^{(n)}_{(2)}}{\partial \sigma_{12}}} \right]
%\end{equation}

It is also noted that, even--though in theory the shear stress ($\sigma_{12}$)
ranges  from  -$\infty$  to +$\infty$, for simulation purpose the domain between
$-0.1~{\rm  MPa}$   and  $+0.1~{\rm  MPa}$  is chosen. This choice is based upon
the material properties of the
example  problem, and by considering the practical range of shear stress ($\sigma_{12}$).
%
In    addition   to  that, the Dirac delta function used in initial condition of
pre--yield (elastic) solution (Eq.~(\ref{eqno_R1}))
is numerically  approximated  using  a Gaussian function with zero mean and
standard deviation of $0.00001~{\rm MPa}$, as shown in Fig.~\ref{figure:DiracDelta}.
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.8\textwidth]{ElasticInitCond-view2}
\vspace*{-0.8cm}
\caption{Approximation  of  Dirac  delta Initial Condition (Eq.~(\ref{eqno_R1}))
with a Gaussian Function of Zero Mean and Standard Deviation
of 0.00001 MPa}
\label{figure:DiracDelta}
\end{center}
\end{figure}


While this approximation introduces some error in solution (PDF of shear stress)
the error quickly disappears as soon as the solution moves away from the initial condition.





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


\subsection{Simulation Plots}

We  begin  by  presenting  the  results of elastic and elastic--plastic probability
density functions for shear stress with only one random material parameter.
%
Fig.~\ref{figure:LowOCR_RandomG}(a)  shows  the evolution of probability density
function   (PDF)   of   shear   stress   with   time/shear   strain  for  a  low
over--consolidation ratio (OCR) Cam--Clay material with normally
distributed  random  shear  modulus  $G$  (mean value of $2.5~{\rm MPa}$ and a
standard  deviation  of  $0.5~{\rm MPa}$). The other material properties are considered
deterministic    and    are  given as follows:  slope of critical state line $M=0.6$;
overconsolidation  pressure  $p_0=0.2~{\rm MPa}$; applied confinement pressure
$p=0.1~{\rm  MPa}$;    slopes of the unloading--reloading  and normal compression
lines in $e-\ln(p^{\prime})$ space $\kappa=0.05$ and $\lambda=0.25$, respectively;
and initial  void  ratio  $e_0=2.18$.  

As    mentioned    earlier,  the  strain is assumed deterministic (strain driven
algorithm)  while  the  strain  rate is assumed constant ($d\epsilon_{12}/dt$ in
Eqs.~(\ref{eqno_A10})  and  ~(\ref{eqno_A11})).  The  numerical  value of assumed
constant  strain  rate  can  be  substituted  in equations ~(\ref{eqno_A10}) and
~(\ref{eqno_A11})  for  the  entire  simulation  of  the  evolution  of  PDF. As
mentioned    earlier    at   the   beginning   of   this   section,   the   FPKE
(Eq.~(\ref{eqno_A9}))   describes   the   evolution   of  PDF  of  shear  stress
($\sigma_{12}$)    with    time    ($t$),    while   the   shear   strain   rate
($d\epsilon_{12}/dt$)  describes the evolution of shear strain ($\epsilon_{12}$)
with  time  ($t$)  also. Combining the two, the evolution of PDF of shear stress
($\sigma_{12}$) with shear strain ($\epsilon_{12}$) is easily obtained. Time has
been brought in this simulation as an intermediate dimension to help in solution
process,    and    hence,   the   numerical   value   of   shear   strain   rate
($d\epsilon_{12}/dt$)  could  be any arbitrary value. This value of shear strain
rate  will  cancel  out  once  one converts the evolution of PDF of shear stress
($\sigma_{12}$) as a function of time, to evolution of PDF of shear stress
($\sigma_{12}$)  as a function of shear strain ($\epsilon_{12}$). For simulation
purpose    we    assumed    the   arbitrary   value   of   shear   strain   rate
($d\epsilon_{12}/dt$) as $0.054/{\rm second}$.


Once  the  evolution  of  PDF  of  shear stress ($\sigma_{12}$) is obtained, the
evolution  of  mean,  mode  and  standard deviations were calculated by standard
operations          on     the     PDF.     The    results    are    shown    in
Fig.~\ref{figure:LowOCR_RandomG}.          In         particular,        
Fig.~\ref{figure:LowOCR_RandomG}(a) shows the evolution (surface) of PDF  of  shear  stress
($\sigma_{12})$ with time ($t$)/shear strain ($\epsilon_{12}$). The evolution of PDF
of shear stress is given in more detail in Fig.~\ref{figure:LowOCR_RandomG}(b)
showing  contours  of the PDF of shear stress ($\sigma_{12}$) and also the mean,
mode, and standard deviations. The deterministic
solution,        obtained   using  the  mean  value  of  shear  modulus ($G$) as
fixed, deterministic input is also shown. As can be seen
from Fig.~\ref{figure:LowOCR_RandomG}(b), the standard deviation of shear stress
increases  in  elastic  regime  (before  approx. $0.056~{\rm s}$) until the yield
point and then decreases as the solution approaches the critical
state  line.  In  other  words, close to critical state the uncertainty in shear
modulus  ($G$)  didn't have much effect on the uncertainty in shear stress as it
had  close  to  yielding.  But on the other hand if one follows the evolution of
deterministic  solution  one  could  see  that even--though the deterministic solution
coincides  with mean and mode for pre-yield (elastic) behavior, it deviates for post-yield
behavior, suggesting that the deterministic solution is neither the mean nor the
most probable solution for post-yield response.
%
%
%\pagebreak
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.7\textwidth]{RandomG_ElasticPlasticPDF-View1-m}
\
\vspace*{-0.1truecm}
\
\begin{center}
(a)
\end{center}
\
\vspace*{-0.1truecm}
\
\includegraphics[width=0.55\textwidth]{RandomG_ContourElasticPlasticPDF-m}
\
\vspace{-0.1truecm}
\
\begin{center}
(b)
\end{center}
\caption{Low  OCR  Cam  Clay  Response  with  Random  Normally Distributed Shear
Modulus  ($G$): (a) Evolution of PDF and (b) Evolution of Contours of PDF, Mean,
Mode,   Standard   Deviations,   and  Deterministic  Solution  of  Shear  Stress
($\sigma_{12}$) with Time (t)/Shear Strain ($\epsilon_{12}$)}
\label{figure:LowOCR_RandomG}
\end{center}
\end{figure}


When           compared        with       Monte--Carlo       simulation, shown in
Fig.~\ref{figure:LowOCR_RandomG_MC}, it can be seen that FPKE approach
predicted  the mean behavior exactly but it slightly over-predicted the standard
deviation  behavior.  This difference is attributed to the approximation used to represent the
initial  deterministic  condition $\sigma_{12}=~0$ (Dirac  delta function). One
may   note   that   at   $\epsilon_{12}=~0$,  the  probability  of  shear  stress
($\sigma_{12}$)  at  $\sigma_{12}=~0$  should  theoretically  be  1  i.e. all the
probability  mass  should  theoretically  be concentrated at $\sigma_{12}=~0$ and
would  be  best  described  by  a  Dirac delta condition. However, for numerical
simulation  of  FPKE, that Dirac delta initial condition was approximated with a
Gaussian  function of mean zero and standard deviation of 0.00001 $MPa$ as shown in
Fig.~\ref{figure:DiracDelta}.  This initial error in standard deviation advected
and  diffused  into  the  domain during the simulation of the evolution process.
This  error  could  be minimized by better approximating the Dirac delta initial
condition (but at the expense of computational cost).
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.7\textwidth]{RandomG_ContourElasticPlasticPDF_MC-m}
\caption{Comparison  of  FPK  Approach  and Monte--Carlo Approach for Low OCR Cam
Clay  Response  with Random Normally Distributed Shear Modulus ($G$) in terms of
Evolution  of  Mean  and Standard Deviation of Shear Stress ($\sigma_{12}$) with
Time ($t$)/Shear Strain ($\epsilon_{12}$)}
\label{figure:LowOCR_RandomG_MC}
\end{center}
\end{figure}



Next few examples examine  results  (PDF  of  shear  stress) by the introduction of
uncertainties in other material parameters (other  than shear modulus $G$).
%
In  examples  where  we  use  more  than  one random material property (see next
examples)   they   are  assumed  uncorrelated  and  independent.  The  presented
development  can  also  deal  with  correlated  random variables, as described by
Fokker--Planck--Kolmogorov  (FPK)  equation  \ref{eqno_S20}.  The  advection and
diffusion  coefficients $N_{(1)}$ and $N_{(2)}$ (in Eq.~(\ref{eqno_R1a})) will
be changed accordingly for any correlation between random soil properties.


Additional uncertainty is introduced for the slope of critical state line $M$
in  terms  of a normal  distribution with mean value  of  0.6  and  a  standard 
deviation  of  0.1. The shear modulus $G$ is again treated as random, as before.
The  resulting  PDF  of  shear  stress  $\sigma_{12}$    is  now  exhibiting   a
non-symmetric distribution in the post-yield (elastic--plastic) region.
This  non--symmetry  is  evident from different  post-yield evolution of mean and
mode of shear stress $\sigma_{12}$ as seen in Fig.~\ref{figure:LowOCR_RandomG_RandomM}(b).  
The  non--symmetry is also evident from the trace of the surface of PDF of shear
stress  at  time ($t) =~0.3{\rm  s}$  (or  shear  strain ($\epsilon_{12})=~1.62\%$) in
Fig.~\ref{figure:LowOCR_RandomG_RandomM}(a).
%
Another       interesting       aspect       to       note       (refer       to
Fig.~\ref{figure:LowOCR_RandomG_RandomM}(b))    is   that   for   this   assumed
combination  of  material  properties (random $G$ and random $M$) the post-yield
deterministic  shear  stress ($\sigma_{12}$) under-predicted the post-yield most
probable  (mode)  shear stress ($\sigma_{12}$) but over-predicted the mean shear
stress ($\sigma_{12}$).
%
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.7\textwidth]{RandomG_RandomM_ElasticPlasticPDF-View1-m}
\
\vspace*{-0.1truecm}
\
\begin{center}
(a)
\end{center}
\
\vspace*{-0.1truecm}
\
\includegraphics[width=0.55\textwidth]{RandomG_RandomM_ContourElasticPlasticPDF-m}
\
\vspace{-0.1truecm}
\
\begin{center}
(b)
\end{center}
\caption{Low  OCR  Cam  Clay  Response  with  Random  Normally Distributed Shear
Modulus  ($G$)  and  Random  Normally  Distributed  Slope of Critical State Line
($M$):  (a)  Evolution  of PDF and (b) Evolution of Contours of PDF, Mean, Mode,
Standard  Deviations, and Deterministic Solution of Shear Stress ($\sigma_{12}$)
with Time (t) /Shear Strain ($\epsilon_{12}$)}
\label{figure:LowOCR_RandomG_RandomM}
\end{center}
\end{figure}


The  effect  of  addition  of  further  uncertainty into the system, in terms of
random normally distributed overconsolidation pressure ($p_0$) with a mean value
of  $0.2~{\rm  MPa}$  and standard deviation of $0.07~{\rm MPa}$, in addition to
previously  introduced  random  shear modulus ($G$) and random slope of critical
state line ($M$), is shown in Fig.~\ref{figure:LowOCR_RandomG_RandomM_RandomP0}.
%
It can be seen  by  comparing  Figs.~\ref{figure:LowOCR_RandomG_RandomM}  and
\ref{figure:LowOCR_RandomG_RandomM_RandomP0},     the     additional
uncertainty   in   overconsolidation   pressure   ($p_0$)   didn't affect much  the
probabilistic  response  of  shear  stress  ($\sigma_{12}$) when compared to its
response  to random normally distributed shear modulus ($G$) and random normally
distributed       slope       of       critical       state      line      ($M$)
(Fig.~\ref{figure:LowOCR_RandomG_RandomM}).    That    is,   the   responses   in
Figs~\ref{figure:LowOCR_RandomG_RandomM}                                     and
\ref{figure:LowOCR_RandomG_RandomM_RandomP0}   are  quite   similar, at  least
qualitatively while quantitatively there are some small differences.
%
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.7\textwidth]{RandomG_RandomM_RandomP0_ElasticPlasticPDF-View1-m}
\
\vspace*{-0.1truecm}
\
\begin{center}
(a)
\end{center}
\
\vspace*{-0.1truecm}
\
\includegraphics[width=0.55\textwidth]{RandomG_RandomM_RandomP0_ContourElasticPlasticPDF-m}
\
\vspace{-0.1truecm}
\
\begin{center}
(b)
\end{center}
\caption{Low  OCR  Cam  Clay  Response  with  Random  Normally Distributed Shear
Modulus  ($G$),  Random Normally Distributed Slope of Critical State Line ($M$),
and   Random   Normally  Distributed  Overconsolidation  pressure  ($p_0$):  (a)
Evolution  of  PDF  and  (b)  Evolution of Contours of PDF, Mean, Mode, Standard
Deviations, and Deterministic Solution of Shear Stress ($\sigma_{12}$) with Time
(t)/Shear Strain ($\epsilon_{12}$)}
\label{figure:LowOCR_RandomG_RandomM_RandomP0}
\end{center}
\end{figure}


The  uncertainty  in  overconsolidation  pressure  ($p_0$) didn't even have much
effect  on  the  probabilistic  response (PDF) of  shear  stress ($\sigma_{12}$) when
considered    separately  with  random normally distributed shear modulus ($G$).
This can  be  seen by  comparing  Fig.~\ref{figure:LowOCR_RandomG_RandomP0},  where  the
probabilistic  behavior  of  shear  stress  ($\sigma_{12}$) for low OCR cam clay
model  with  randomly  distributed  shear  modulus  ($G$)  and  random  normally
distributed    overconsolidation    pressure    ($p_0$)   was   presented,   and
Fig.~\ref{figure:LowOCR_RandomG}.  Similar  to  the case where the shear modulus
was  the only random parameter (Fig.~\ref{figure:LowOCR_RandomG}), here also the
post-yield   deterministic   shear   stress  ($\sigma_{12}$)  overpredicted  the
post-yield   mean   and   post-yield   most   probable   (mode)   shear   stress
($\sigma_{12}$).
%
%\pagebreak
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.7\textwidth]{RandomG_RandomP0_ElasticPlasticPDF-View2-m}
\
\vspace*{-0.1truecm}
\
\begin{center}
(a)
\end{center}
\
\vspace*{-0.1truecm}
\
\includegraphics[width=0.55\textwidth]{RandomG_RandomP0_ContourElasticPlasticPDF-m}
\
\vspace{-0.1truecm}
\
\begin{center}
(b)
\end{center}
\caption{Low  OCR  Cam  Clay  Response  with  Random  Normally Distributed Shear
Modulus   ($G$)  and  Random  Normally  Distributed  Overconsolidation  pressure
($p_0$):  (a) Evolution of PDF and (b) Evolution of Contours of PDF, Mean, Mode,
Standard  Deviations, and Deterministic Solution of Shear Stress ($\sigma_{12}$)
with Time (t)/Shear Strain ($\epsilon_{12}$)}
\label{figure:LowOCR_RandomG_RandomP0}
\end{center}
\end{figure}


The        probabilistic  Cam--Clay  elastic--plastic responses (low OCR)  with 
different  degrees  of  randomnesses  (as  presented  above)  were  compared  in
Fig.~\ref{figure:LowOCR_PDFComparison}.
In that figure the PDFs of shear stress ($\sigma_{12}$) at shear strain
($\epsilon_{12})~=1.62\%$ are  shown and compared with the deterministic value of
shear stress obtained by choosing mean value (which in this case is also mode as
we used normal distribution) of material properties. 
%
It  is very interesting to note that the
deterministic  value  of shear stress ($\sigma_{12}$) is different from the most
probable  (mode)  value  of  shear  stress ($\sigma_{12}$) for all the cases. It
either under-predicted or over-predicted the most probable (mode) value.
%
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.85\textwidth]{LowOCR_ComparisonPDF-m}
\caption{Comparison  of  Shear Stresses at 1.62\% Shear Strain Obtained from Low
OCR Cam Clay Model with different degrees of Randomnesses}
\label{figure:LowOCR_PDFComparison}
\end{center}
\end{figure}


Presented  methodology  is  applicable  to  both hardening (as shown in examples
above) and softening material response.
%
A  high  OCR  Cam--Clay  material sample was analyzed. Following parameters were
used: random normally distributed shear modulus
($G$) with  mean  of  $2.5~{\rm MPa}$  and standard deviation of $0.5~{\rm MPa}$; random
normally  distributed  slope of critical state line ($M$) with mean of $0.6$ and
standard    deviation    of    $0.1$;  deterministic  overconsolidation pressure
$p_0=0.8~{\rm   MPa}$;     deterministic    applied    confinement    pressure  
$p=0.02~{\rm MPa}$; deterministic 
slopes of the unloading--reloading  and normal compression
lines in $e-\ln(p^{\prime})$ space $\kappa=0.05$ and $\lambda=0.25$, respectively;
deterministic initial void ratio $e_0=2.18$. The (arbitrary)  value  of 
strain rate is still the same, $d\epsilon_{12}/dt=0.054/{\rm s}$.          
The          resulting PDF of shear stress is     shown         in
Fig.~\ref{figure:HighOCR_RandomG_RandomM}(a). Fig.~\ref{figure:HighOCR_RandomG_RandomM}(b) 
gives  a more detailed view of results,
in terms of evolution of contours of PDF, mean, mode, standard deviation and
deterministic  solution  of shear stress ($\sigma_{12}$). Similar to the low OCR
simulation  with random $G$ and random $M$, the high OCR simulation also yielded
a  symmetric  distribution  of PDF shear stress ($\sigma_{12}$) in the pre-yield
(elastic)  regime  and  a  non-symmetric  distribution of PDF shear stress
($\sigma_{12}$)   in  the  post-yield  (elastic-plastic)  regime.  Similarly  to
probabilistic response for low  OCR  (with  random  $G$  and  random  $M$),  the
deterministic  shear  stress  ($\sigma_{12}$)  coincided  with the mean and most
probable  (mode)  shear  stress  in the pre-yield (elastic)
regime  but  deviated  in  the post-yield non-linear elastic-plastic regime. The
deterministic shear stress ($\sigma_{12}$) under-predicted the mean shear stress
($\sigma_{12}$)  in  the  entire  post-yield  regime but over-predicted the most
probable  (mode)  shear  stress ($\sigma_{12}$) in the region close to the yield
point,  though  close to critical state line it under-predicted the most probable
(mode) shear stress ($\sigma_{12}$).
%
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.65\textwidth]{HighOCR_RandomG_RandomM_ElasticPlasticPDF-View1-m}
\
\vspace*{-0.1truecm}
\
\begin{center}
(a)
\end{center}
\
\vspace*{-0.1truecm}
\
\includegraphics[width=0.55\textwidth]{HighOCR_RandomG_RandomM_ContourElasticPlasticPDF-m}
\
\vspace{-0.1truecm}
\
\begin{center}
(b)
\end{center}
\caption{High  OCR  Cam  Clay  Response  with  Random Normally Distributed Shear
Modulus  ($G$)  and  Random  Normally  Distributed  Slope of Critical State Line
($M$):  (a)  Evolution  of PDF and (b) Evolution of Contours of PDF, Mean, Mode,
Standard  Deviations, and Deterministic Solution of Shear Stress ($\sigma_{12}$)
with Time (t) /Shear Strain ($\epsilon_{12}$)}
\label{figure:HighOCR_RandomG_RandomM}
\end{center}
\end{figure}


Verification  of  probabilistic  simulations for high OCR sample were again done
using   Monte--Carlo  approach.  Fig~\ref{figure:HighOCR_RandomG_RandomM_MC}
compares the evolution of mean and
standard deviation of shear stress ($\sigma_{12}$) obtained using FPKE approach
and   Monte--Carlo   approach.  It  can  be  seen  that  FPKE  approach  slightly
overpredicted  the  standard deviation behavior because of the reasons discussed
earlier.  
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.7\textwidth]{HighOCR_RandomG_RandomM_ContourElasticPlasticPDF_MC-m}
\caption{Comparison of FPK Approach and Monte--Carlo Approach for High OCR Cam Clay Response with Random Normally Distributed
Shear Modulus ($G$) and Random Normally Distributed Slope of Critical State Line ($M$)
in terms of Evolution of Mean and Standard Deviation of Shear Stress ($\sigma_{12}$) with Time ($t$)/Shear Strain 
($\epsilon_{12}$)}  
\label{figure:HighOCR_RandomG_RandomM_MC}
\end{center}
\end{figure}
%
That  is,  the  FPKE  approach is somewhat ``wider'' since the initial
condition  was  not  a Dirac delta function, but rather an approximation. Again,
this deviation in initial condition (from Dirac delta function) can be controlled
by  choosing  an initial normal distribution for PDF of shear stress with smaller
standard  deviation,  but  this  will increase computational cost in solving the
resulting finite difference system of equations.


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

\section{CONCLUSIONS}

In this paper we  presented    methodology   to   solve  the  probabilistic
elastic--plastic   differential   equations.   The   methodology   is  based  on
Eulerian--Lagrangian   form  of  the  Fokker--Planck--Kolmogorov  equation  and
provides  for  full  description  of the probability density function (PDF) of stress
response for  given strain.
%
Of  particular interest was the modeling of stress--strain constitutive response
for    materials    with    uncertain   material parameters. These uncertainties in
material  properties are present in most  materials, and it is just our
simplification    to    a    single    value    (say mean)  that  renders 
deterministic properties usually used for simulations
%
The     effects  of  nonlinear  hardening/softening  on  predicted  PDF  of  stress  were
investigated  in  some  more  detail. It was shown that the nonlinear hardening/softening
creates  a  discrepancy  between  the  most  likely stress solution and the
deterministic  solution, meaning that the usually obtained deterministic
solution is not the most 
likely   outcome  of  a  constitutive  stress-strain  simulation  if  material
paramaters  are  uncertain.

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

\bibliography{SFEM,refmech}
%\bibliographystyle{acm}
\bibliographystyle{plain}
%\bibliographystyle{alpha}
%\bibliographystyle{agsm}


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

%\include{Appendix}

\begin{appendix}
\section*{NOTATION}

This is a summary of notations used in mathematical derivations:

\begin{itemize}

\item[]  $D_{ijkl}$: Material parameter tensor, function of $D^{el}_{ijkl}$, $f$, $U$,
$q_*$, and $r_*$, 

\item[]  $D^{el}_{ijkl}$:  Elastic stiffness tensor,

\item[]  $D^{ep}_{ijkl}$:  Elastic--plastic  tangent stiffness tensor,

\item[] $\sigma_{ij}$: Cauchy stress tensor, in indicial notation,

\item[] $e$: Void ratio of soil,

\item[] $e_0$: Initial void ratio of soil,

\item[] $f$: Yield surface,

\item[] $G$: Shear modulus of soil,

\item[] $K$: Bulk modulus of soil,

\item[] $K_P$: Plastic modulus,

%\item[] $A_{kl}$: Second-order tensor, equal to 
%        $\displaystyle \frac{\partial f}{\partial \sigma_{pq}} D^{el}_{pqkl}$,
%
%\item[] $B$: Scalar-valued quantity, equal to 
%$\displaystyle\frac{\partial f}{\partial \sigma_{rs}} D^{el}_{rstu}
%\displaystyle\frac{\partial f}{\partial \sigma_{tu}}$,

\item[] $C_{ijm}$: Coefficient tensor in multivariate It\^o equation, in indicial notation,


\item[] $M$: Slope of critical state line in $p-q$ space,

\item[] $N_{(1)}$: Advection coefficient of FPK equation,

\item[] $N_{(2)}$: Diffusion coefficient of FPK equation,

\item[] $P(\cdot, t)$: Probability density at time, t,

\item[] $p$: first invariant of stress tensor,

\item[] $p_0$: Hydrostatic pressure, the soil under analysis has experienced before (internal variable in Cam Clay model),

\item[] $\bar p_0$: Direction of evolution of $p_0$,

\item[] $q$: Invariant of stress tensor, equal to $\sqrt{\frac{3}{2}s_{ij} s_{ij}}$

\item[] $q_*$: Internal variables,

\item[] $r_*$: Directions of evolution of internal variables,

\item[] $s_{ij}$: second invariant of the deviatoric stress tensor,

\item[] $t$, $t_1$, $t_2$: time,

\item[] $U$:  Potential surface,

\item[] $W_i$: Increment of vector Wiener process, 

\item[] $x_t$: Spatial location at time, t,

\item[] $\epsilon_{kl}$: Strain tensor,

\item[] $\dot \epsilon$: Strain rate,

\item[] $\lambda$: Slope of normal consolidation line in $e-ln(p)$ space,

\item[] $\kappa$: Slope of unloading-reloading line in $e-ln(p)$ space,

\item[] $\eta_{ij}$: Spatial-temporal operator tensor, function of $\sigma_{ij}$, $D_{ijkl}$, and $\epsilon_{kl}$,

\item[] $\rho(\cdot, t)$: Phase density at time, t,


\item[] $\tau$: time,

\item[] $\zeta$: Probability current


\end{itemize}


\end{appendix}


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


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

                      
\end{document}
\bye
                                                                                                                                                             

