%%%%%%%%%%%%%%%%%%%%%%% 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: Solution and Verification in 1D}
%\subtitle{Do you have a subtitle?\\ If so, write it here}
\author{
Kallol Sett\inst{1}
\and
Boris Jeremi{\'c}\inst{2}
\and
M. Levent Kavvas\inst{3}
}


\institute{
Graduate  Research  Assistant,  Department    of    Civil    and  Environmental 
Engineering, University of California, Davis, CA 95616.
\and 
Associate  Professor,  Department   of  Civil  and  Environmental  Engineering, 
University of California, Davis, CA 95616. \texttt{jeremic@ucdavis.edu}
\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}

In   this paper,  a  solution  is presented for evolution of Probability Density
Function  (PDF)  of  elastic--plastic  stress--strain  relationship for material
models  with  uncertain  parameters.  Developments  in  this  paper are based on
already  derived  general  formulation  presented  in  the  companion paper. The
solution   presented   is   then   specialized  to  a  specific  Drucker  Prager
elastic--plastic material model.


Three    numerical   problems  are  used  to  illustrate the developed solution.
The  stress  strain  response  (1D) is given as a PDF of stress as a function of
strain.  The  presentation of the stress strain response through the PDF differs
significantly  from  the  traditional  presentation  of  such results, which are
represented  by  a  single, unique curve in stress--strain space. In addition to
that  the  numerical  solutions are verified against closed form solutions where
available  (elastic).  In  cases  where  the closed form solution does not exist
(elastic--plastic), Monte--Carlo simulations are used for verification.


\end{abstract}

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

\section{Introduction}
\label{Introduction}

An    elastic-plastic    constitutive  law  can be  represented by a set of linear or
non-linear   ordinary   differential   equations   (ODEs),   which  relate  rate
(increments) of stress with the rate of strain
through linear or nonlinear material modulus:
%
\begin{equation}
\label{eqno_I1}
\frac{d \sigma_{ij} (t)}{dt}=D^{ep}_{ijkl} \frac{d \epsilon_{kl} (t)}{dt}
\end{equation}
%
where  $D^{ep}_{ijkl}$  could  be  linear  or a non-linear function of stresses,
strains and internal variables. If either the material modulus or the forcing term
(strain  rate)  becomes  random, this  set of linear or nonlinear ODEs becomes a
set of linear or non-linear stochastic differential equations (SDEs).  
%
%These  types of equations are called SDEs with
%random  coefficient  and  SDEs with random forcing respectively. 
%
The uncertainty
associated   with   the   coefficient  (stiffness)  term  is  generally  attributed  to  the
inherent  variability of the material. Geomaterials are particularly notorious for their
variability,  sampling  and  testing errors and in general, uncertainty in their
properties. 
%
The uncertainty in the
forcing  term  arises when the material is subjected to uncertain loads (usually
dynamic) like wind,
waves  or earthquakes. Due to randomness in the parameter and/or forcing term
the response variable of the elastic-plastic constitutive rate equation (stress)
will  then  be  a  random  process.  There exist several methods to estimate the
probabilistic characteristics of the response variable \citep{book:Gardiner}.

For  the  case  where  the  material modulus is linear and deterministic and the
forcing  is  uncertain (Gaussian), the response is known to be Gaussian and can be estimated
by standard methods \citep{book:Gardiner}. General linear SDE with random forcing
can  be  solved  by cumulant expansion method \citep{VanKampen:1976}. 
%
One possible
way  to  solve  non-linear  SDE  with  random forcing is to write its equivalent
Fokker--Planck--Kolmogorov (FPK) form.  The
advantage  of writing the FPK form is that it is linear and deterministic even--though the
original  equation  is non-linear and stochastic. 
%
The general
solution   method   for   FPK equation   can   be   found   in  any  standard  textbook
(\citet{book:Gardiner},    \citet{book:Risken}).   A  solution  scheme  for  FPK
equation with
application    to   structural   reliability   was   presented   by  \citet{Langtangen:1991}


For    the    particular    case    where the forcing is deterministic and  the 
material  linear  elastic (but still uncertain),  Eq.~\ref{eqno_I1}  simplifies 
to a linear set of algebraic equations of the form,
%
\begin{equation}
\label{eqno_I2}
\sigma_{ij}=D^{el}_{ijkl} \epsilon_{kl}
\end{equation}
%
where  $D^{el}_{ijkl}$  is the elastic stochastic moduli tensor and hence the statistics of
the  response  variable (stress) can be easily obtained by transformation method
of  random  variable.  For general linear SDEs with random coefficients cumulant
expansion  method could be used \citep{VanKampen:1976}. 
%
While   the  solution  for  the stochastic linear elastic stress--strain problem
is    readily    available,  the nonlinear (elastic--plastic) stochastic problem
presents  itself  as  much  harder  to solve. This problems involves finding the
solution for a non-linear SDEs with random coefficients. 
%
Fortunately such solution recently developed by \citet{Kavvas:2003}. 
%
The developed solution is presented as a generic Eulerian--Lagrangian form of the
Fokker--Planck--Kolmogorov  equation.    That  probabilistic solution is  second
order  exact  for  any stochastic  nonlinear ODE or PDE with random coefficients
and random forcing.


The  probabilistic  solution  developed by \citet{Kavvas:2003}  was used in  the
companion paper \citet{Jeremic2005a} to develop the probability density function
(PDF)  of  a  general  local-average  form  of elastic-plastic constitutive rate
equation.  This  Eulerian--Lagrangian  FPK  equation was then specialized to the
particular  cases  of  point-location  scale  linear  elastic and Drucker-Prager
associative linear hardening elastic-plastic constitutive rate equations to show
the applicability of the general formulation. In this paper the solution process
of those particular FPK equations will be presented.



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

\section{Fokker--Planck--Kolmogorov Equation for 
Probabilistic  Elasticity and Elasto--Plasticity in 1-D}
\label{FokkerPlanckEquations} 

In     the companion paper, \cite{Jeremic2005a} applied the Eulerian--Lagrangian
form of 
Fokker--Planck--Kolmogorov (FPK)  equation   to  the  description  of  the  probabilistic 
behavior of elastic and elastic--plastic (Drucker-Prager associative
linear  hardening) 1D constitutive equations with random
material  parameters  and  random  strain  rate. 


By focusing  attention to the
randomness  of  material  properties  only  (i.e.  assuming the forcing function
(strain  rate)  as  deterministic),  partial differential
equation (PDE) describing the evolution of probability density function (PDF) of stress 
can  be  simplified. In particular, for 1D case, and for linear elastic material
(but  still  with  probabilistic material properties, in this case shear modulus
$G$) one can write the following PDE
%(Eq. (23) of \cite{Jeremic2005a}) 
%
\begin{eqnarray}
\nonumber
\displaystyle \frac{\partial P(\sigma_{12}(t))}{\partial t}  
= 
&-& 
\left<G \displaystyle \frac{d\epsilon_{12}}{dt} 
\right> \displaystyle \frac{\partial P(\sigma_{12}(t))}{\partial \sigma_{12}} 
\\ \nonumber
&+& 
\left\{\int_{0}^{t} d \tau Cov_0 
\left[G \displaystyle \frac{d\epsilon_{12}}{dt}; G 
\displaystyle \frac{d \epsilon_{12}}{dt} \right] \right\} 
\frac{\partial^2 P(\sigma_{12}(t))}{\partial \sigma_{12}^2} 
\\
\label{eqno_FPE1}
\end{eqnarray}

                                   
Similarly,  for  elastic--plastic  state, again by  neglecting the randomness in
strain rate, one can write the PDE for evolution of PDF of stress  in 1D as
%zone for materials obeying Drucker-Prager associative
%linear   hardening   type   shear  constitutive  rate  equations  
%%(Eq.  (35)  of \cite{Jeremic2005a})
%  
\begin{eqnarray}
& &
\frac{\partial P(\sigma_{12}(t))}{\partial t}
= 
-\left< \left(G^{ep}(t) \right) \frac{d\epsilon_{12}}{dt} \right> \
\displaystyle \frac{\partial P(\sigma_{12}(t))}
{\partial \sigma_{12}}
\nonumber \\
&+&  
\left\{\int_{0}^{t} d \tau Cov_0
\left[
G^{ep}(t)
\frac{d\epsilon_{12}}{dt};
%\left(
G^{ep}(t-\tau)
 \displaystyle \frac{d\epsilon_{12}}{dt}
 \right] \right \} \frac{\partial^2 P(\sigma_{12}(t))}{\partial \sigma_{12}^2}
\nonumber \\
\label{eqno_FPE2}
\end{eqnarray}
%
where  $G^{ep}(a)$ is the probabilistic elastic--plastic tangent stiffness,
(given 
%(Eq. (34) 
in \cite{Jeremic2005a})
%
\begin{equation}
G^{ep}(a) 
= 
G
-
\displaystyle 
\frac{G^2}{G+9K\alpha^2+\displaystyle \frac{1}{\sqrt{3}}I_1(a)\alpha'}
\label{eqno_FPE2a}
\end{equation}
%
where in the previous equation (\ref{eqno_FPE2a}), $a$ assumes values $t$ or $t-\tau$.
%It  is  important to note that the pre-yield, linear elastic
%probabilistic  behavior  of shear  stress  for this elastic--plastic material
%is governed  by  Eq.~(\ref{eqno_FPE1}).  
With  appropriate initial and boundary
conditions    as    described    in    \cite{Jeremic2005a},    one   can   solve
Eqs.~(\ref{eqno_FPE1})  and  ~(\ref{eqno_FPE2})  for  evolution  of PDF of shear
stress with shear strain.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Example Problem statements}
\label{ExampleProblemStatement}



The applicability of proposed  FPK equations
(Eqs.~(\ref{eqno_FPE1})  and  ~(\ref{eqno_FPE2}))  in  describing  probabilistic
elasto-plastic behavior, is verified using the following three example problems.
%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{\bf Problem I.}
%
Assume the material is linear elastic, probabilistic, with probabilistic shear  modulus  ($G$)
given by a normal distribution at a point--location scale with mean of $2.5$~MPa
and standard deviation of $0.707$~MPa. The aim is to calculated the evolution of
PDF  of  shear  stress ($\sigma_{12}$) with shear strain ($\epsilon_{12}$) for a
displacement-controlled  test  with  deterministic  shear  strain increment. The
other  parameters  are  considered  deterministic  and are as follows: Poisson's
ratio ($\nu = 0.2$, and confining pressure $ I_1 =0.03$~MPa.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
{\bf Problem II.}
%
Assume        elastic--plastic    material  model, composed of linear    elastic
component   and      Drucker--Prager   associative  isotropic  linear  hardening
elastic--plastic    component.   The  probabilistic shear modulus ($G$) is given
through a normal distribution at a point--location scale with mean of $2.5$~MPa 
and standard deviation of $0.707$~MPa. The aim is to  calculate the evolution of the
PDF of shear stress ($\sigma_{12}$) with shear strain
($\epsilon_{12}$) for a displacement-controlled test with deterministic shear
strain  increment.  The other parameters are considered deterministic and are as
follows:  Poisson's  ratio  $ \nu = 0.2$, confining pressure $I_1 = 0.03$~MPa, 
yield  parameter\footnote{The yield parameter $\alpha$ is an internal   variable
and is a   function   of   friction   angle $\phi$ given by ($\alpha  =  2
\sin(\phi)/(\sqrt(3)(3-sin  \phi))$ (e.g. \citep{book:Chen})}  $\alpha  = 0.071$, 
plastic slope\footnote{The plastic slope  $\alpha'$  is a rate of change of friction angle
 governing linear  hardening.} $ \alpha' =5.5$.

{\bf Problem III.}
%
Assume  elastic--plastic  material  model,  with  linear   elastic component and
Drucker--Prager      associative  isotropic  linear  hardening  elastic--plastic
component.  The  probabilistic yield   parameter  ($\alpha$)  is given through a
normal distribution  at  a
point-location  scale with mean of $0.52$ and standard deviation of $0.1$.
The aim is to calculate the  evolution of the PDF of shear stress ($\sigma_{12}$) with
shear  strain  ($\epsilon_{12}$)  for  a  displacement-controlled  test with
deterministic  shear  strain  increment.  The  other  parameters  are considered
deterministic  and  are as follows: shear modulus $G = 2.5$~MPa, Poisson's ratio
$\nu    =    0.2$,  confining  pressure  $I_1  =0.03$~MPa, and the plastic slope
$\alpha' =5.5$.


The  above three problems  will  be solved using the proposed FPK
equation  approach.  In  addition  to  that,  the solution will verified  using 
either    variable    transformation    method,  for  linear    elastic  case or
repetitive Monte Carlo type simulations for elastic-plastic case.


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


\section{Numerical Scheme for Solving Fokker--Planck--Kolmogorov Equation}
\label{NumericalScheme}

For  probabilistic elastic  and  elastic--plastic  constitutive rate equations, 
the PDEs (Eqs.~(\ref{eqno_FPE1}) and (\ref{eqno_FPE2}))  which describe the
evolution of probability densities of $ \sigma_{12} $ have the following general
form:
%
\begin{equation}
\label{eqno_NS1}
\frac{\partial P}{\partial t}
=
-
N_{(1)}\frac{\partial P}{\partial \sigma_{12}}
+
N_{(2)}\frac{\partial^2 P}{\partial \sigma_{12}^2}
\end{equation}
%
with  an  appropriate initial condition, which depends on the type of
problem, and boundary conditions in the form
%
\begin{equation}
\label{eqno_NS2}
\zeta(-\infty,t)=\zeta(\infty,t)=0
\end{equation}
%
where  $N_{(1)}$ and $N_{(2)}$ are called advection and diffusion
coefficients,  respectively.  The  above  PDE  system  (Eqs.~(\ref{eqno_NS1}) and
~(\ref{eqno_NS2}) with appropriate initial condition) were solved numerically by
Method  of Lines \cite{book:Mathematica} using commercially available software
Mathematica  \cite{software:Mathematica}. 

%In order to cone must first define the
%domain of the state variables. 
The  stress  (state) variable $ \sigma_{12} $ theoretically spans space from  $-
\infty$ to $+ \infty$. However, for simulation
(and    practical)   purposes,  this theoretical   domain is reduced to between 
$-0.1$~MPa and $+0.1$~MPa. This reduction is based on the material properties of
the   example   problems   and   span  the  practical  range  of  shear  stress,
$\sigma_{12}$. The Fokker--Planck--Kolmogorov PDE was semi-discretized
(Fig.~\ref{figure:StressDiscretization})  in  stress ($ \sigma_{12} $) domain by
finite difference technique to obtain a set of linear simultaneous ODE systems.
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.9\textwidth]{StressDiscretization.eps}
\caption{Stress Domain Discretization of Fokker--Planck--Kolmogorov PDE}
\label{figure:StressDiscretization}
\end{center}    
\end{figure}
%
This  set  of  linear  simultaneous  ODEs  is  solved  using  central difference
technique.    By    referring    to   Fig.~\ref{figure:StressDiscretization}   a
semi-discretized form of Eq.~ (\ref{eqno_NS1}) can be written at any
intermediate node $ i $ as,
%
\begin{equation}
\frac{\partial P_i}{\partial t} 
= 
P_{i-1} \left( \frac{N_{(1)}}{2 \Delta \sigma_{12}} 
+ 
\frac{N_{(2)}}{\Delta \sigma_{12}^2}\right) 
- 
P_i \left( \frac{2N_{(2)}}{\Delta \sigma_{12}^2}\right) 
+ 
P_{i+1} \left( -\frac{N_{(1)}}{2 \Delta \sigma_{12}} 
+ 
\frac{N_{(2)}}{\Delta \sigma_{12}^2}\right)
\label{eqno_NS3}
\end{equation}
%
Previous  discretized  system of equations forms an  initial  value  problem  in 
the time dimension. By using forward
difference  technique,  one  can  introduce  the  boundary condition at the left
boundary (node 1 in Fig.~\ref{figure:StressDiscretization}) as,
%
\begin{equation}
\label{eqno_NS4}
N_{(1)}P_1-N_{(2)}\frac{P_2-P_1}{\Delta \sigma_{12}}=0
\end{equation}
%
or, after rearranging,
%
\begin{equation}
\label{eqno_NS5}
P_1
=
P_2\left( \displaystyle \frac{\displaystyle \frac{N_{(2)}}{\Delta \sigma_{12}}}{N_{(1)} 
+ 
\displaystyle \frac{N_{(2)}}{\Delta \sigma_{12}}}\right) 
=
0
\end{equation}
%
Similarly,  using  backward difference technique, one can introduce the boundary
condition       at      the      right      boundary      (node      $n$      in
Fig.~\ref{figure:StressDiscretization}) as,
%
\begin{equation}
\label{eqno_NS6}
N_{(1)}P_n-N_{(2)}\frac{P_n-P_{n-1}}{\Delta \sigma_{12}}=0
\end{equation}
% 
which, after rearranging becomes
%
\begin{equation}
\label{eqno_NS7}
P_n
=
P_{n-1}\left( \displaystyle \frac{\displaystyle \frac{N_{(2)}}
{\Delta \sigma_{12}}}{\displaystyle \frac{N_{(2)}}{\Delta \sigma_{12}}
-
N_{(1)}}\right) 
=
0
\end{equation}

The   initial   condition  depends  on  the  type  of the problem and  it  could  be
deterministic  or  random.  For  elastic  constitutive rate equation with random
shear  modulus  (Problem  I)  and for pre-yield elastic-plastic linear hardening
constitutive  rate  equation  with random shear modulus (elastic part of Problem
II)  the  initial  condition  is  deterministic.  It  will,  therefor,   be best
represented as Dirac delta function of the form,
%
\begin{equation}
\label{eqno_NS8}
P(\sigma_{12}) = \delta(\sigma_{12})
\end{equation}

For simulation purpose the Dirac delta initial condition was approximated with a
Gaussian function. That is, for Problem I, the initial condition was approximated with
a  Gaussian  function  with  mean  of 0 and standard deviation of $0.00001$~MPa as
shown in Fig.~\ref{figure:ElasticInitialCondition}
%                 
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.9\textwidth]{ElasticInitCond-view2}
\caption{Approximation  of  Dirac  delta function, used  as an initial condition
for Problem I}
\label{figure:ElasticInitialCondition}
\end{center}
\end{figure}

For     post-yield,  probabilistic  elastic-plastic  behavior (plastic  part  of
Problem  II),  the  initial  condition  is  random  and  it  corresponds  to the
probability  density  function  of $\sigma_{12}$ just prior to yielding, that is
obtained from elastic part of Problem II. For Problem-III, the pre-yield elastic
part  is  deterministic but initial condition for the post-yield elastic-plastic
response  is  random  and  corresponds  to  the  assumed  distribution  in yield
strength.
                                     


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%  
\section{Determination of Coefficients for Fokker--Planck--Kolmogorov Equation}
\label{CoefficientEstimation}

To  solve  Problems  I,  II,  and III, the advection and
diffusion    coefficients $ N_{(1)} $ and $ N_{(2)} $ must be determined for all
three  problems.  For  sake of simplicity, a  constant  strain  rate is assumed 
and hence, terms containing $d \epsilon_{12} / dt$ in coefficients of Eqs.~ (\ref{eqno_FPE1}) and
~(\ref{eqno_FPE2})  can  be  substituted  by  a constant numerical value for the
entire  simulation  of  the  evolution  of  PDF.  It should be noted  that  the 
FPK equation 
(Eqs.~(\ref{eqno_FPE1})  or  ~(\ref{eqno_FPE2})) describes the evolution of PDFs
of  stress with time, while, similarly, strain rate describes the evolution of strain with time.
Combining the two, the  evolution of PDF of stress with strain can be obtained.
Time has been brought in this simulation as an intermediate dimension to help in
solution  process,  and  hence,  the  numerical  value of strain rate could be any
arbitrary  value,  which will cancel out once the time evolution of
PDF  of  stress  is converted to strain evolution of PDF of stress. For simulation of all the
three        example    problems,    an  arbitrary   value  of strain rate of 
$d \epsilon_{12} / dt =  0.054 1/s$ is assumed.


It  should  also  be  noted  that since  the material properties are
assumed as random
variables  at  a  point-location scale, the covariance terms appearing within the
advection  and diffusion coefficients become variances of random variables.
For  estimations  of  means and variances of functions of random variables (e.g.
for  Problems  II  and  III)  from basic random variables, commercially available
statistical software mathStatica \cite{book:Rose} was used.

Substituting  the values of deterministic and random material properties and the
strain  rate, coefficients $N_{(1)}$ and $N_{(2)}$ of the
FPK equations can be obtained for all problems:
%
\begin{description}
%
\item{Problem I}
%
\begin{eqnarray}
\nonumber
N_{(1)} &=& \left<G \frac{d \epsilon_{12}}{dt} \right> \\
\nonumber
&=& 2 \frac{d \epsilon_{12}}{dt} \left< G \right > \\
\nonumber
&=& 0.27 ~~\rm{MPa/s}
\end{eqnarray}
%
\begin{eqnarray}
\nonumber
N_{(2)}  &=& \int_0^t d \tau Var\left[G \frac{d \epsilon_{12}}{dt}\right] \\
\nonumber
&=& 4 t \left(\frac{d \epsilon_{12}}{dt}\right)^2 Var[G] \\
\nonumber
&=& 0.0058 t ~~(\rm{MPa/s})^2
\end{eqnarray}

\item{Problem II}

For  pre-yield linear elastic case, the coefficients $N_{(1)}$ and $N_{(2)}$
will be the same as those for Problem I. For post-yield elastic-plastic case the
coefficients are 

\begin{eqnarray}
\nonumber
N_{(1)}
&=&
\left<\left(G-\displaystyle 
\frac{G^2}{G+9K\alpha^2+\displaystyle 
\frac{1}{\sqrt{3}}I_1\alpha'}\right) \displaystyle 
\frac{d \epsilon_{12}}{dt} \right> 
\\ \nonumber
&=& 
\displaystyle \frac{d \epsilon_{12}}{dt} \left<G-\displaystyle \
\frac{G^2}{G+9K\alpha^2+\displaystyle 
\frac{1}{\sqrt{3}}I_1\alpha'} \right> 
\\ \nonumber
&=& 
0.147 ~~\rm{MPa/s}
\end{eqnarray}
%
\begin{eqnarray}
\nonumber
N_{(2)} 
&=& 
t \left(\displaystyle 
\frac{d \epsilon_{12}}{dt}\right)^2 Var\left[G-\displaystyle 
\frac{G^2}{G+9K \alpha^2
+
\displaystyle \frac{1}{\sqrt{3}}I_1\alpha'}\right] 
\\ \nonumber
&=& 
0.00074t~~(\rm{MPa/s})^2
\end{eqnarray}

\item{Problem III}
%


For   post-yield  elastic-plastic  simulation  the  coefficients  $N_{(1)}$  and
$N_{(2)}$ are
%
\begin{eqnarray}
\nonumber
N_{(1)}
&=&
\left<\left(G-\displaystyle \frac{G^2}{G+9K\alpha^2
+
\displaystyle \frac{1}{\sqrt{3}}I_1\alpha'}\right) 
\displaystyle \frac{d \epsilon_{12}}{dt} \right> 
\\ \nonumber
&=& 
\displaystyle \frac{d \epsilon_{12}}{dt} 
\left<G-\displaystyle \frac{G^2}{G+9K\alpha^2+\displaystyle 
\frac{1}{\sqrt{3}}I_1\alpha'} \right> 
\\ \nonumber
&=& 
0.2365~~\rm{MPa/s}
\end{eqnarray}
%
\begin{eqnarray}
\nonumber
N_{(2)} 
&=& 
t \left(\displaystyle 
\frac{d \epsilon_{12}}{dt}\right)^2 Var\left[G-
\displaystyle \frac{G^2}{G+9K
\alpha^2+\displaystyle \frac{1}{\sqrt{3}}I_1\alpha'}\right] 
\\ \nonumber
&=& 
0.0001t~~(\rm{MPa/s})^2
\end{eqnarray}

\end{description}


It should be noted that for Problem III, since the shear modulus is deterministic, the
pre-yield elastic case is deterministic.


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%\input{Application2-D}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
%\input{NumericalSchemeLinear}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Results and Verifications of Example Problems} 


In      this  section  results  are  presented  for elastic and elastic--plastic
probabilistic  1D  problem.  The  results  are  obtained  by  using FPK equation
approach   described   in   previous   sections   and   in   the  companion  paper
\citep{Jeremic2005a}. In addition to that, the Monte Carlo based verification of
developed  solutions  (results)  is  presented.  The effort to  verify  developed
solutions    (that    are   based    on  FPK  approach)  plays a crucial role in
presented  development  of  probabilistic  elasto--plasticity  as  there  are no
previously  published solutions which could have  been used for verification. 
In    addition  to  that,  verification  and  validation efforts should always 
be included in any modeling and simulations work \citep{Oberkampf2002}.
%
%



For  linear elastic constitutive rate  equations  (Problem-I  and pre-yield case
of  Problem-II)  the  verification  is performed by comparing solutions obtained
through  the  use  of FPK equation approach with high accuracy (exact) solution,
using  a  transformation method of random variables \citep{book:Montgomery}. This
method  is  applicable  as  for rate-independent linear elastic case the 1D shear
constitutive equation simplify to a linear algebraic equation of the form,
%
\begin{equation}
\label{eqno_R1}
\sigma_{12}=G \epsilon_{12}=u(G,\epsilon_{12})
\end{equation}
%
Using the definition of strain rate, the above equation can be written 
in terms of time $t$ as,
%
\begin{equation}
\label{eqno_R2}
\sigma_{12}=G (0.054 t)=v(G,t)
\end{equation}
%
where,  $0.054$~1/s  is  the  arbitrary  strain-rate  assumed  for  this example
problem.  
%
According  to  the  transformation method of random variables
\citep{book:Montgomery}, and,  given  the  continuous  random  variable (shear modulus)
$G$,   with   PDF  $g(G)$  and  Eqs.~(\ref{eqno_R1})  or  ~(\ref{eqno_R2})  as
one-to-one  transformations  between  the  values of random variables of $G$ and
$\sigma_{12}$,  one  can  obtain  the   PDF of shear stress ($\sigma_{12}$),
$P(\sigma_{12})$ as,
%
\begin{equation}
\label{eqno_R3}
P(\sigma_{12})
=
g(u^{-1}(\sigma_{12},\epsilon_{12}))\left|J\right|
\end{equation}
%
which  will  allow  for  predicting  the  evolution of PDF of $\sigma_{12}$ with
$\epsilon_{12}$ or,
%
\begin{equation}
\label{eqno_R4}
P(\sigma_{12})
=
g(v^{-1}(\sigma_{12},t))\left|J\right|
\end{equation}


Eq.~(\ref{eqno_R4}) will predict the evolution of PDF of $\sigma_{12}$ with $t$.
In     Eqs.~(\ref{eqno_R3})    and    ~(\ref{eqno_R4}),   functions  
$G=u^{-1}(\sigma_{12},\epsilon_{12})$   or  $G=u^{-1}(\sigma_{12},t)$  are  the 
inverse of functions
$\sigma_{12}=u(G,\epsilon_{12})$   or   $\sigma_{12}=v(G,t)$   respectively  and
$J=du^{-1}(\sigma_{12},\epsilon_{12})/d\sigma_{12}$ and
$J=dv^{-1}(\sigma_{12},t)/d\sigma_{12}$   are   their  respective  Jacobians  of
transformations.


For  non-linear elastic-plastic constitutive rate equations (post-yield cases of
Problems  II  and  III) the verification is done using Monte-Carlo simulation
technique by generating sample data for material properties from standard normal
distribution  and  by  repeating solution of the deterministic elastic-plastic constitutive
rate  equation  for  each  data generated above. The probabilistic
characteristics  of  resulting  random stress variable for each time (or strain)
step  are  then  easily  computed.  A  relatively  large  number  of  data points
(1,000,000) were generated for each material constant random
variable for this simulation purpose.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Problem I}


The  evolution  of  PDF  of  shear  stress  with  time and shear strain is shown
in Figures~\ref{figure:ElasticPDF} and ~\ref{figure:ActualElasticPDF}.
%
Presented  PDFs are for linear  elastic  material with random shear modulus, and
were  obtained  using  \textit{FPE  approach} (Fig.~~ref{figure:ElasticPDF}) and
\textit{transformation method} (Fig.~\ref{figure:ActualElasticPDF}).
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=\textwidth]{ElasticPDF-m}
\caption{Evolution  of    PDF    of  shear stress versus strain  (or  time)  for
linear  elastic  material  model  with random shear modulus (Problem I) obtained
using FPK equation approach.}
\label{figure:ElasticPDF}
\end{center}
\end{figure}
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=\textwidth]{ActualElasticPDF-m}
\caption{Evolution  of PDF of shear stress versus strain (or  time)  for  linear
elastic  material  model  with  random  shear modulus (Problem I) obtained using
transformation method.}
\label{figure:ActualElasticPDF}
\end{center}
\end{figure}


The contours      of      evolution      of      PDFs   are  compared     in
Fig.~\ref{figure:ComparisonContourElastic_Actual_FPE}. 
%
Similarly, comparison of the evolution of mean
and         standard         deviations        are  shown  in
Fig.~\ref{figure:ComparisonElasticMeanStdDeviation_Actual_FPE}.  
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=\textwidth]{ComparisonContourElastic_Actual_FPE-m}
\caption{Comparison of Contours of Time (or Strain) Evolution of Probability Density Function for Shear Stress for Elastic Constitutive Rate Equation 
with Random Shear Modulus (Problem--I) for FPE Solution and Variable Transformation Method Solution.}
\label{figure:ComparisonContourElastic_Actual_FPE}
\end{center}
\end{figure}
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=\textwidth]{ComparisonElasticMeanStdDeviation_Actual_FPE-m}
\caption{Comparison of Mean and Standard Deviation of Shear Stress for Elastic Constitutive Rate Equation 
with Random Shear Modulus (Problem--I) for FPE Solution and Variable Transformation Method Solution.}
\label{figure:ComparisonElasticMeanStdDeviation_Actual_FPE}
\end{center}
\end{figure}
%
It  can be seen
from  the comparison figure that even-though the FPK approach predicted the mean behavior
exactly, it slightly over-predicted the standard deviation. This is because of
the  approximation used to represent the Dirac delta function, which was used as
the initial condition for the FPK. 
%
One may note that at $\epsilon_{12}=0$, the
probability   of shear stress $\sigma_{12}$ should theoretically
be  $1$  i.e. all the probability mass should theoretically be concentrated
at  $\sigma_{12}=0$. As such, it would be best described by the  Dirac delta function. 
%
However, 
for  numerical  simulation  of  FPK,  Dirac  delta function as initial  condition  was
approximated  with  a  Gaussian  function of mean zero and standard deviation of
$0.00001$~MPa, as shown in Fig.~\ref{figure:ElasticInitialCondition}. 
%
This error in the initial condition advected and diffused into the domain with the
simulation  of the evolution  process.  This  error  could  be  minimized by better
approximating    the    Dirac    delta    initial    condition  (but  at  higher
computational cost). 
%
The effect of approximating the initial condition of the
PDF   of   shear   stress   at   $\epsilon_{12}=0.0426~\%$   is   shown   in
Fig.~\ref{figure:CrudeComparisonAtYield_Actual_FPE}.  
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=\textwidth]{CrudeComparisonAtYield_Actual_FPE-m}
\caption{Effect of Approximating Function of Dirac Delta Initial Condition : PDF of Stress at Yield for Different Approximation of 
Initial Condition with Actual (Variable Transformation Method) Solution).}
\label{figure:CrudeComparisonAtYield_Actual_FPE}
\end{center}
\end{figure}
%
In  this figure the actual
PDF  at  $\epsilon_{12}=0.0426~\%$  obtained  using the transformation method was
compared  with  the  PDFs  at  $\epsilon_{12}=0.0426~\%$  obtained  using the FPK
approach  with  three different approximate initial conditions - all having zero
mean but standard deviations of $0.01$~MPa, $0.005$~MPa and $0.00001$~MPa.



One  may  also  note  that finer approximation of initial condition necessitates
finer  discretization  of  stress  domain  close to (or at) $\sigma_{12}$ =0. The finite
difference  discretization  scheme  adopted  here  uses the same fine discretization
uniformly      all  throughout  the  entire  domain. It is noted that that fine,
uniform discretization is not needed (and is quite expensive) in later stages of
calculation  of  evolution  of PDF, but is kept the same for simplicity sake. 
%
In  presented  examples, to   properly   capture   the   approximate   initial  
condition (as shown in Fig.~\ref{figure:ElasticInitialCondition}),  
the    stress   domain between $-0.1$~MPa  and $+0.1$~MPa was discretized with a
uniform step size of $0.000005$~MPa and hence there is a total  of $40,000$~nodes. 
%
This not only requires large computational effort
but is also  very memory sensitive. 
%
An adaptive discretization technique will be a
much better  approach to solving this  problem. 
Current work is going on in formulating an
adaptive algorithm for the solution of this type of problem.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Problem II}

The solution to this problem involves the solving two FPK equations, one corresponding to
the  pre-yield  elastic  part  and  the  other  corresponding  to the post-yield
elastic-plastic  part.  
%
%
The  elastic  part  of  this  problem  is  identical  to
Problem--I.  The initial condition for the post-yield elastic-plastic part of the
problem  is random and is shown in Fig.~\ref{figure:PlasticInitialCondition}. 
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=\textwidth]{PlasticInitCond-m}
\caption{Initial condition for FPK equation for elastic--plastic zone (Problem--II).}
\label{figure:PlasticInitialCondition}
\end{center}
\end{figure}
%
It  may    be  noted that this initial condition corresponds to the PDF of shear
stress ($P(\sigma_{12}$)) at yield obtained from the solution of FPK equation of
the pre-yield elastic part.
%
A view of the surface of evolution of the PDF of shear stress versus shear strain
(time)  is  shown  in Fig.~\ref{figure:ElasticPlasticPDF-View1}. Another view to
the     PDF     of     stress--strain     surface    is    shown    in
Fig.~\ref{figure:ElasticPlasticPDF-View2}.
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=\textwidth]{ElasticPlasticPDF-m}
\caption{Evolution of PDF of shear stress versus strain (time) for elastic-plastic material with 
random shear modulus (Problem--II). View 1.}
\label{figure:ElasticPlasticPDF-View1}
\end{center}
\end{figure}
%
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=\textwidth]{ElasticPlasticPDF-View2-m}
\caption{Evolution   of   PDF   of   shear  stress  versus  strain  (time)  for
elastic-plastic  material  with random shear modulus (Problem-II). View 2.}
\label{figure:ElasticPlasticPDF-View2}
\end{center}
\end{figure}
%
It  is noted  that the yielding of this   material   occurred   at  $t$=0.00789 
second (which is equivalent to $\epsilon_{12}$= 0.0426 $\%$).
%
The  evolution contours  for PDF of shear stress versus strain (time) along with
the      mean      and      standard      deviations      are      shown      in
Fig.~\ref{figure:ContourElasticPlasticPDF-Detailed}.
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=\textwidth]{ContourElasticPlasticPDF-Detailed-m}
\caption{Contour of evolution of PDF for shear stress versus strain (time) for elastic-plastic material  
with random shear modulus (Problem--II).}
\label{figure:ContourElasticPlasticPDF-Detailed}
\end{center}
\end{figure}
%
It  can  be  seen  from  that figure that, as expected, the evolution of mean of
shear  stress  changes  slope  after  the  material yielded. Another interesting
aspect to note is the relative slope of the evolution of standard deviation with
respect  to  the  evolution of mean. The relative slope in the pre-yield elastic
zone  increases at a higher rate during the evolution process when compared with
that in the post-yield elastic-plastic zone.
%
In  other words, in the evolution process the post-yield elastic--plastic
constitutive  rate  equation  did not amplify the initial uncertainty as much as
the  pre-yield elastic constitutive rate equation did. This can be easily viewed
from  Fig.~\ref{figure:ComparisonPDF-ExtendedElastic-ElasticPlastic}  where  the
post-yield  elastic-plastic  evolution  of PDF of shear stress was compared with
fictitious  extension of  elastic evolution of PDF. 
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=\textwidth]{Plastic_ExtendedElastic_PDFComparison-m}
\caption{Comparison  of  evolution  of  PDF  for  elastic-plastic  material  and
extended elastic material cases for random shear modulus.}
\label{figure:ComparisonPDF-ExtendedElastic-ElasticPlastic}
\end{center}
\end{figure}
%
Comparing  the  PDF  of  shear  stress  at    $\epsilon_{12}=0.0804\%$ (which is
equivalent  to  $t=0.01489  s$), one can conclude that the variance of predicted
elastic-plastic shear stress is much smaller (i.e. prediction is less uncertain)
as compared to the same if the material were modeled as completely elastic.


            

Fig.~\ref{figure:ComparisonElasticPlasticMeanStdDeviation_Actual_FPE}   compares
the  evolution  of  means  and  standard  deviations  of  predicted shear stress
obtained  using  FPK  equation  approach  and  transformation  method (pre-yield
behavior)/Monte-Carlo approach (post-yield behavior).
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=\textwidth]{ComparisonElasticPlasticMeanStdDeviation_Actual_FPE-m}
\caption{Comparison of mean and standard deviation of shear stress for plastic constitutive rate equation 
with random shear modulus (problem-ii) for FPK equation solution and Monte Carlo simulation solution.}
\label{figure:ComparisonElasticPlasticMeanStdDeviation_Actual_FPE}
\end{center}
\end{figure}
%
Although in the pre--yield
response  the  FPK equation  approach  over-predicted the evolution of standard deviations
because  of  reasons  discussed  earlier, in the post-yield response it matched
closely    at    regions  further  from the yielding region. The somewhat larger
difference between FPK equation solution and the verification one (Monte Carlo solution) close to
the yielding  region  is attributed to the fact that the initial condition for solution
of  post-yield  elastic-plastic  FPK equation was obtained from the solution of pre-yield
elastic FPK equation. 
%
One  way  to  better predict the overall probabilistic elastic-plastic behavior,
would  probably  be  to  obtain  the  pre-yield  elastic  behavior  through  the
transformation  method  and  then  use  the  FPK  approach to predict post-yield
elastic-plastic behavior.


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

\subsection{Problem III} 

In  this problem, the pre--yield linear elastic part is deterministic, however, at yield
there  is  a  distribution (with very small standard deviation)  in shear stress due to assumed distribution in yield
parameter  $\alpha$.  The distribution in shear stress corresponds to the PDF of
the  random  variable  $\alpha  I_1$ (first invariant of the
stress  tensor or mean confining stress) and is assumed to be
deterministic.  This  PDF of shear stress at yield was assumed to be the initial
condition  for  the  solution  of post-yield elastic-plastic FPK equation and is shown in
Fig.~\ref{figure:RandomYieldInitialCondition}.
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=\textwidth]{RandomYieldInitialCondition-m}
\caption{Initial condition for FPK equation for elastic--plastic material with random yield strength (Problem--III).}
\label{figure:RandomYieldInitialCondition}
\end{center}
\end{figure}



The  evolution   of   PDF for  shear   stress versus strain (time)  is   shown   in
Fig.~\ref{figure:RandomYieldPlasticPDF-View1}. In addition to that the contours (including mean and standard deviation) of the evolution of
PDF  for  shear  stress  versus strain (time) are shown in
Fig.~\ref{figure:RandomYieldContourPlasticPDF-Detailed}. 
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=\textwidth]{RandomYieldPlasticPDF-m}
\caption{Evolution of PDF for shear stress for elastic--plastic material with 
random yield strength (Problem--III) (only plastic zone is shown).}
\label{figure:RandomYieldPlasticPDF-View1}
\end{center}
\end{figure}
%
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=\textwidth]{RandomYieldContourPlasticPDF-Detailed-m}
\caption{Contour   PDF   for      shear      stress  versus  strain  (time)  for
elastic--plastic material with random yield strength (Problem--III).}
\label{figure:RandomYieldContourPlasticPDF-Detailed}
\end{center}
\end{figure}



Looking           at Fig.~\ref{figure:RandomYieldContourPlasticPDF-Detailed} and
comparing  the  slopes  of  evolution  of  mean  and standard deviation, one can
conclude  that  the elastic-plastic evolution process didn't amplify the initial
uncertainty  in yield strength significantly. The initial (at yield) probability
density  function  of  shear  stress  just  advected  into the domain during the
elastic--plastic      evolution      process      without     diffusing     much.
Fig.~\ref{figure:RandomYieldPlasticPDF-View1}   clearly   shows  this  advection
process.  The evolution of mean and standard deviations of shear stress obtained
from  the  FPK equation approach was compared with those obtained from the Monte
Carlo           simulation           and           is          shown          in
Fig.~\ref{figure:ComparisonMeanStdDeviation_Actual_FPE}.
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=\textwidth]{RandomYieldComparisonElasticMeanStdDeviation_Actual_FPE-m}
\caption{Comparison of Mean and Standard Deviation of Shear Stress for Elastic-Plastic Constitutive Rate Equation 
with Random Yield Strength (Problem-III) for FPE Solution and Monte Carlo Simulation Solution}
\label{figure:ComparisonMeanStdDeviation_Actual_FPE}
\end{center}
\end{figure}



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



In  this  paper  a  solution  was presented for the evolution of the Probability
Density  Function  (PDF) of elastic--plastic stress--strain relationship, in 1D.
The  solution  was  based  on  expressions  developed  in  a companion paper and
specialized  to  the  Drucker  Prager elasto--plastic material model with linear
isotropic hardening.

Three    numerical   problems were used to illustrate the developed solution and
discuss  the  general  behavior  of  elastic--plastic  materials  which  exhibit
uncertainty  in  material  parameters.  The solutions to numerical problems were
verified  against  closed, analytical forms, where available, while Monte--Carlo
simulations were used for all other verifications.

The approach to solving probabilistic elastic-plastic problems presented here is
quite  unique  and  shows great promise in dealing with general 3D probabilistic
constitutive  problems. Subsequently, the developed methodology is to be used in
solving  general,  probabilistic  elastic--plastic boundary value problems using
the finite element method. Current work is progressing in that direction.




%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
\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)




                    
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\begin{thebibliography}{12}
\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[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[Van~Kampen(1976)]{VanKampen:1976}
N.~G. Van~Kampen.
\newblock Stochastic differential equations.
\newblock \emph{Phys. Rep.}, 24:\penalty0 171--228, 1976.

\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[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[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[Jeremi{\'c} et~al.(2006)Jeremi{\'c}, Sett, and Kavvas]{Jeremic2005a}
Boris Jeremi{\'c}, Kallol Sett, and M.~Levent Kavvas.
\newblock Probabilistic elasto-plasticity: Formulation of evolution of
  probability density function.
\newblock \emph{Acta Geotechnica}, 2006.
\newblock In Review.

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

\bibitem[Wolfram(1991)]{book:Mathematica}
Stephen Wolfram.
\newblock \emph{Mathematica: A System for Doing Mathematics by Computer}.
\newblock Addison-Wesley, Redwood City, second edition, 1991.

\bibitem[Wolfram Research~Inc.(2003)]{software:Mathematica}
.~Wolfram Research~Inc.
\newblock \emph{Mathematica Version 5.0}.
\newblock Wolfram Research, Inc., Champaign, Illonois, 2003.

\bibitem[Rose and Smith(2002)]{book:Rose}
Colin Rose and Murray~D. Smith.
\newblock \emph{Mathematical Statistics with Mathematica}.
\newblock Springer Texts in Statistics. Springer-Verlag, New York, 2002.

\bibitem[Oberkampf et~al.(2002)Oberkampf, Trucano, and Hirsch]{Oberkampf2002}
William~L. Oberkampf, Timothy~G. Trucano, and Charles Hirsch.
\newblock Verification, validation and predictive capability in computational
  engineering and physics.
\newblock In \emph{Proceedings of the Foundations for Verification and
  Validation on the 21st Century Workshop}, pages 1--74, Laurel, Maryland,
  October 22-23 2002. Johns Hopkins University / Applied Physics Laboratory.

\bibitem[Montgomery and Runger(2003)]{book:Montgomery}
Douglas~C. Montgomery and George~C. Runger.
\newblock \emph{Applied Statistics and Probability for Engineers}.
\newblock John Wiley {\&} Sons, 605 Third Avenue, New York, NY 10158, third
  edition, 2003.

\end{thebibliography}

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

%\input{Appendix}

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

                      
\end{document}
\bye
                                                                                                                                                             

