%% This is file `elsarticle-template-1-num.tex',
%%
%% Copyright 2009 Elsevier Ltd
%%
%% This file is part of the 'Elsarticle Bundle'.
%% ---------------------------------------------
%%
%% It may be distributed under the conditions of the LaTeX Project Public
%% License, either version 1.2 of this license or (at your option) any
%% later version.  The latest version of this license is in
%%    http://www.latex-project.org/lppl.txt
%% and version 1.2 or later is part of all distributions of LaTeX
%% version 1999/12/01 or later.
%%
%% The list of all files belonging to the 'Elsarticle Bundle' is
%% given in the file `manifest.txt'.
%%
%% Template article for Elsevier's document class `elsarticle'
%% with numbered style bibliographic references
%%
%% $Id: elsarticle-template-1-num.tex 149 2009-10-08 05:01:15Z rishi $
%% $URL: http://lenova.river-valley.com/svn/elsbst/trunk/elsarticle-template-1-num.tex $
%%
%%\documentclass[preprint,12pt]{elsarticle}

%% Use the option review to obtain double line spacing
%%\documentclass[preprint,review,12pt,twocolumn]{elsarticle}

%% Use the options 1p,twocolumn; 3p; 3p,twocolumn; 5p; or 5p,twocolumn
%% for a journal layout:
%% \documentclass[final,1p,times]{elsarticle}
%% \documentclass[review,1p,times,twocolumn]{elsarticle}
%% \documentclass[final,3p,times]{elsarticle}
%% \documentclass[final,3p,times,twocolumn]{elsarticle}
%% \documentclass[final,5p,times]{elsarticle}
%
%
%\documentclass[preprint,review,5p,times,twocolumn]{elsarticle}
\documentclass[preprint,review,12pt]{elsarticle}
%
\graphicspath{
{/home/jeremic/tex/works/Papers/2010/SSEPFEM_static/figures/}
%{/home/kallol/publication/2010/JournalPaper/CMAME/figures/}
}
%% if you use PostScript figures in your article
%% use the graphics package for simple commands
%% \usepackage{graphics}
%% or use the graphicx package for more complicated commands
%% \usepackage{graphicx}
%% or use the epsfig package if you prefer to use the old commands 
\usepackage{epsfig}

%% The amssymb package provides various useful mathematical symbols
\usepackage{amssymb}
%% The amsthm package provides extended theorem environments
%% \usepackage{amsthm}

%% The lineno packages adds line numbers. Start line numbering with
%% \begin{linenumbers}, end it with \end{linenumbers}. Or switch it on
%% for the whole article with \linenumbers after \end{frontmatter}.
%% \usepackage{lineno}

%% natbib.sty is loaded by default. However, natbib options can be
%% provided with \biboptions{...} command. Following options are
%% valid:

%%   round  -  round parentheses are used (default)
%%   square -  square brackets are used   [option]
%%   curly  -  curly braces are used      {option}
%%   angle  -  angle brackets are used    <option>
%%   semicolon  -  multiple citations separated by semi-colon
%%   colon  - same as semicolon, an earlier confusion
%%   comma  -  separated by comma
%%   numbers-  selects numerical citations
%%   super  -  numerical citations as superscripts
%%   sort   -  sorts multiple citations according to order in ref. list
%%   sort&compress   -  like sort, but also compresses numerical citations
%%   compress - compresses without sorting
%%
%% \biboptions{comma,round}

% \biboptions{}

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

\journal{Comput. Methods Appl. Mech. Engrg.}

\begin{document}

\begin{frontmatter}

%% Title, authors and addresses

%% use the tnoteref command within \title for footnotes;
%% use the tnotetext command for the associated footnote;
%% use the fnref command within \author or \address for footnotes;
%% use the fntext command for the associated footnote;
%% use the corref command within \author for corresponding author footnotes;
%% use the cortext command for the associated footnote;
%% use the ead command for the email address,
%% and the form \ead[url] for the home page:
%%
%%\title{Title\tnoteref{label1}}
\title{Stochastic Elastic--Plastic Finite Elements}
%% \tnotetext[label1]{}
%% \author{Name\corref{cor1}\fnref{label2}}
\author[label2]{Kallol Sett\corref{cor1}}
\ead{sett@uakron.edu}

\author[label3]{Boris Jeremi\'{c}}
\ead{jeremic@ucdavis.edu}
%% \ead[url]{home page}

\author[label3]{M. Levent Kavvas}
\ead{mlkavvas@ucdavis.edu}

\cortext[cor1]{corresponding author}

\address[label2]{Department of Civil Engineering, The University of Akron, 302 Buchtel Common, Akron, OH 44325}
\address[label3]{Department of Civil and Environmental Engineering, University of California, One Shields Avenue, Davis, CA 95616}
%\address[label3]{Department of Civil and Environmental Engineering, University of California, One Shields Avenue, Davis, CA 95616}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





\begin{abstract}
%% Text of abstract

A computational framework has been developed for simulations of the behavior of 
solids and structures made of stochastic elastic--plastic materials. 
%
Uncertain elastic--plastic material properties are modeled as random fields, which
appear as the coefficient term in the governing partial differential equation of mechanics. 
%
A spectral stochastic elastic--plastic finite element method with Fokker--Planck--Kolmogorov equation
based probabilistic constitutive integrator is proposed
for solution of this non--linear (elastic--plastic) partial differential equation with
stochastic coefficient.
%
To this end, the random field material properties are discretized, 
in both spatial and stochastic dimension, into finite numbers of
independent basic random variables, using
Karhunen--Lo\`{e}ve expansion.
%
Those random variables are then propagated through the elastic--plastic
constitutive rate equation using Fokker--Planck-Kolmogorov equation approach, to
obtain the evolutionary material properties, as the material plastifies.
%
The unknown displacement (solution) random field is then assembled, as a function of
known basic random variables with unknown deterministic coefficients, using polynomial chaos expansion.
%
The unknown deterministic coefficients of polynomial chaos expansion are obtained,
by minimizing the error of finite representation, by Galerkin technique.



\noindent
%
\rule[.5mm]{\textwidth}{.5mm}


\small{Published as:  Kallol Sett, Boris Jeremi{\'c} and M. Levent Kavvas. 
Stochastic elastic-plastic finite elements. 
Computer Methods in Applied Mechanics and Engineering, 
200(9-12):997-1007, February 2011.}

\noindent
%
\rule[.5mm]{\textwidth}{.5mm}



\end{abstract}

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

\begin{keyword}
%% keywords here, in the form: keyword \sep keyword

%% MSC codes here, in the form: \MSC code \sep code
%% or \MSC[2008] code \sep code (2000 is the default)

elasto--plasticity \sep 
Fokker--Planck--Kolmogorov equation \sep 
Karhunen--Lo\`{e}ve expansion \sep 
polynomial chaos expansion \sep 
stochastic finite elements


\end{keyword}

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

\end{frontmatter}

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

In  mechanics,  simulations  of  static or dynamic behavior of
solids  and  structures  involve  solutions  of  boundary value
problems, which are comprised of the equilibrium equation, $A \sigma
= \phi(t)$, together with the strain compatibility equation, $B
u  =  \epsilon$,  and  the  constitutive  equation, $\sigma = D
\epsilon$,  along with a set of additional restraints (boundary
conditions).  In  the above, $ \sigma $ is generalized stress, $
\phi(t)  $ is generalized force that can be time (t) dependent,
$ u $ is generalized displacements, $ \epsilon $ is generalized
strain  and  $A$,  $B$,  and  $D$  are operators which could be
linear or non-linear.
%

Rigorous  mathematical  theory  has been developed for problems
where   the   only  random  parameter  is  the  external  force
$\phi(t)$.  In this case, the probability distribution function
(PDF)  of  the  response  variable satisfies a  Fokker--Planck-Kolmogorov (FPK)  partial
differential    equation  (cf. Soize \cite{Soize:1994}).   With
appropriate  initial and boundary conditions the FPK PDE can be
solved  for  PDF  of  response variable. The numerical solution
method  of  FPK equation corresponding to structural dynamics problems was described
by   number   of   researchers   (e.g.,  Langtangen \cite{Langtangen:1991},
Masud and Bergman \cite{Bergman:2005}).


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


Monte   Carlo   simulation   technique  is  an  alternative  to
analytical  solution  of  partial  differential  equation  with
stochastic coefficient.
%
A thorough review of different aspects of formulation of Monte
Carlo  technique  for stochastic mechanics problem was presented
in a state-of-the-art report edited by Schueller \cite{Schueller:1997}.
%
Monte Carlo technique has been widely used for probabilistic solution
of uncertain boundary value problems (Paice et al. \cite{Griffiths:1996},
Popescu et al. \cite{Popescu:1997},
Mellah et al. \cite{Mellah:SFEMgeotechApplication},
DeLima et al. \cite{DeLima:SFEMgeotechApplication},
Koutsourelakis et al. \cite{Deodatis:2002},
Griffiths et al. \cite{Griffiths:2002},
Nobahar \cite{Nobahar:2003},
Fenton and Griffiths \cite{Griffiths:2003, Griffiths:2005}).
%
It has the advantage that
accurate  solution  can  be  obtained  for  any  problem  whose
deterministic  solution  (either  analytical  or  numerical) is
known.
%
However, the  major  disadvantage  of  Monte  Carlo analysis is the
repetitive  use  of  the deterministic model until the solution
variable becomes statistically significant. The computational cost
associated with it could be very exorbitantly high (and probably intractable) for three-dimensional 
and/or non--linear problems with multiple uncertain material properties.


The   difficulties with  analytical   solution   and   the  high
computational  cost  associated with Monte Carlo technique lead
to  the  development  of  numerical methods for the solution of
stochastic  differential  equation with random coefficient. 
%
For   stochastic  boundary  value  problems,  stochastic  finite
element method (SFEM) is the most popular.
%
There   exist   several   formulations  of  SFEM,  among  which
perturbation
%
(Kleiber and Hien \cite{book:Kleiber},                
Der Kiureghian and Ke \cite{DerKiureghian:1988};
Mellah et al. \cite{Mellah:SFEMgeotechApplication},  
Gutierrez and de Borst \cite{DeBorst:1999}) 
%
and spectral
%
(Ghanem and Spanos \cite{book:Ghanem},
Keese and Matthies \cite{keese:SPDE},
Xiu and Karniadakis \cite{Karniadakis:2003}, 
Debusschere et al. \cite{Debusschere:2003},
Anders and Hori \cite{Hori:2000})  
%
methods  are  the most common. 
%
Nice  reviews  on  the advantages  and the disadvantages of different
formulations   of   SFEM   were  provided  by  Matthies  et  al.
\cite{Matthies:1997reviewSFEM} and recently, by Stefanou \cite{Stefanou:2009}.
%
Mathematical  issues  regarding  different formulations of SFEM
were  addressed by Deb et al. \cite{babuska:GalerkinSPDE} and by
Babuska and Chatzipantelidis \cite{Babuska:SPDE}.
%
Computational issues were discussed by Ghanem \cite{Ghanem:1999} and by Stefanou \cite{Stefanou:2009}.
%
Recently, a hybrid treatment of various forms of epistemic uncertainties,
including modeling error, was proposed by Soize and Ghanem \cite{SoizeGhanem:2009}.
%
However,  most  of  the  existing formulations  are  for  linear elastic
problems.
%
Though   there  has  been  some  published  works  on  geometric
non--linear     problems     (Liu     and     Der    Kiureghian
\cite{DerKiureghian:1991},       Keese       and       Matthies
\cite{keese:SPDE},  Keese  \cite{keese:2003ReportSFEM}),  there
exist   only  few  published  papers  on  material  non--linear
(elastic--plastic) problems with uncertain material parameters.

The major difficulty in extending the available formulations of
SFEM to general elastic--plastic problem is the high non--linear
coupling in the elastic--plastic constitutive rate equation.
%
First    attempt    to    propagate    uncertainties    through the
elastic--plastic   constitutive  equation  considering  random
Young's    modulus    was   published   only   recently, by
Anders and Hori \cite{Hori:1999,  Hori:2000}. They took perturbation expansion
at  the  stochastic mean behavior and considered only the first
term of the expansion. In computing the mean behavior they took
advantage of bounding media approximation. Although this method
doesn't  suffer  from  computational difficulty associated with
Monte Carlo method for problems having no closed-form solution,
it  inherits  ``closure  problem"  and the ``small coefficient of
variation" requirements for the material parameters.
%
Closure problem refers to the need for higher order statistical
moments  in  order to calculate lower order statistical moments
(cf. Kavvas \cite{Kavvas:2003}).
%
The  small  COV requirement claims that the perturbation method
can  be  used  (with  reasonable  accuracy)  for  probabilistic
simulations  of solids and structures with uncertain properties
only if their COVs are less than 20\% (cf. Sudret and Der Kiureghian \cite{DerKiureghian:2000ReportSFEM}).
For  soils  and  other natural materials, COVs are rarely below
$20\%$       (cf. Baecher and Christian \cite{book:Christian}).
%
Furthermore,  with  bounding  media  approximation,  difficulty
arises  in  computing  the  mean  behavior  when  one considers
uncertainties  in  internal  variable(s) and/or direction(s) of
evolution of internal variable(s).
%
These difficulties, associated with probabilistic simulation of
the elastic--plastic constitutive equation, prevent application
of SFEM to general geomechanics problems.

%Rephrased following Kavvas
Recently,  a  special  nonlocal Eulerian-Lagrangian form of the
Fokker-Planck-Kolmogorov  (FPK)  equation was derived by Kavvas
\cite{Kavvas:2003}   in   order   to  model  the  probabilistic  behavior  of
nonlinear/linear   physical   systems   that   have   uncertain
parameters,  uncertain  forcing functions and uncertain initial
conditions, and that are described by conservation equations in
the  form  of  nonlinear/linear stochastic partial differential
equations  or  stochastic ordinary differential equations. This
Eulerian-Lagrangian  form  of  the  FPK  was  later  applied by
Jeremic  et  al. \cite{Jeremic:2007} to  the  probabilistic  description  of
elasto-plasticity,  by  writing the generalized elastic-plastic
constitutive rate equation in the probability density space.
%
%
%Recently,  in  dealing  with uncertainties in soil constitutive
%behavior,  Jeremi\'{c}  et  al.  \cite{Jeremic:2007}  developed
%Eulerian--Lagrangian  form  of Fokker--Planck--Kolmogorov (FPK)
%equation      (cf.     Kavvas     \cite{Kavvas:2003})     based
%\textit{probabilistic  elasto--plasticity},  by writing the the
%generalized  elastic--plastic constitutive rate equation in the
%probability  density  space.  
%
The  non--linear  (in real space)
elastic--plastic  constitutive  rate  equation becomes a linear
partial differential equation in the probability density space,
simplifying  the  solution  process.  The  resulting FPK
partial    differential    equation   is   second-order   exact
(analytical).  Further,  this  FPK  approach can very easily be
specialized   to   any   particular   constitutive  model.  The
probabilistic  monotonic  behavior  of  Drucker--Prager linear
hardening  and  Cam  Clay  models were discussed by Sett et al.
\cite{Sett:2007a,Sett:2007b}.      Jeremi\'{c}     and     Sett
\cite{Jeremic:2009b}   later   extended  the  FPK  equation  based
\textit{probabilistic     elasto--plasticity}     to    include
probabilistic   yielding   and   in   another  paper  Sett  and
Jeremi\'{c} \cite{Sett:2010a} discussed its influence on cyclic
constitutive   behavior   of   geomaterials.   In   the  above
publications,  it was shown that with probabilistic approach to
geomaterial  modeling,  even the simplest elastic--perfectly plastic
model  picks  up  some of the important features of geomaterial
behavior, which deterministically, could only be possible with
advanced constitutive models.

Encouraged by the probabilistic constitutive responses of geomaterials, in  this
paper,  the  authors,  with  the  broad  goal  of understanding the influence of
spatial  uncertainties  in  probabilistic approach to geomaterial modeling, have
extended  the  spectral  formulation of stochastic finite element method (Ghanem
and Spanos \cite{book:Ghanem}) to general stochastic elastic--plastic problems.
% 
Fokker--Planck--Kolmogorov            equation  based  \textit{probabilistic    
elasto--plasticity}  has  been  used as probabilistic constitutive integrator in
updating the material properties as the material plastifies.
%
Both the formulation and numerical solution scheme are discussed.

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

\section{FPK Equation based \textit{Probabilistic Elasto-Plasticity}}
\label{FPKE}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{General Form of Constitutive Rate Equation}

The constitutive behavior of many materials can be modeled by elastic--plastic constitutive
rate equation, which, in most general form, can be written as: 
%
\begin{equation}
\label{eqno_1}
\displaystyle \frac{d \sigma_{ij}}{dt} = D_{ijkl} \displaystyle \frac{d \epsilon_{kl}}{dt}
\end{equation}
%
where $\epsilon$ is the strain, $t$ is the pseudo time, and $D$ is the
modulus, which could be either elastic or elastic--plastic:
%
\begin{eqnarray}
D_{ijkl} = \left\{\begin{array}{ll}
%
           D^{el}_{ijkl}
%
%
         & \mbox{when~elastic} \\
%
           \\
%
           D^{el}_{ijkl}
           -
           \displaystyle \displaystyle \frac{D^{el}_{ijmn}
           \displaystyle \frac{\partial U}{\partial \sigma_{mn}}
           \displaystyle \frac{\partial f}{\partial \sigma_{rs}}
           D^{el}_{rskl}}
           {\displaystyle \frac{\partial f}{\partial \sigma_{ab}}
           D^{el}_{abcd}
           \displaystyle \frac{\partial U}{\partial \sigma_{cd}}
           -
           \displaystyle \frac{\partial f}{\partial q_*}r_*}
%
         & \mbox{when~elastic--plastic}
%
           \end{array} \right.
\nonumber
\end{eqnarray}
%
where,    $D^{el}$,  $f$,  $U$,  $q_*$,  and $r_*$ are elastic  modulus,  yield 
function,  plastic  potential function,  internal  variable(s),  and  rate(s)  of
evolution   of  internal  variable(s)  respectively.  
%
However,  due  to  various uncertainties associated  with  
our measurement of material  properties, $D$  in  
Eq.~(\ref{eqno_1})  becomes uncertain.  This is especially significant for geomaterials, where 
coefficients of variation of measured properties are typically 50\% or more.
%
In traditional deterministic approach to material modeling, one
typically  applies  engineering  judgment  (qualitative)  in  obtaining  the 'most
probable'  material  parameters  and  substitute in Eq.~(\ref{eqno_1}) to obtain the
'most   probable'  material  constitutive  behavior. 
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Constitutive Equation in Probability Density Space}

In  quantifying  the
uncertainties    in    soil    constitutive   behavior,   Jeremi\'{c}   et   al.
\cite{Jeremic:2007}   proposed  a  probabilistic  approach  to  material
modeling  by  writing  the  constitutive  rate  equation (Eq.~(\ref{eqno_1})) in
the probability   density space using an approach based on Eulerian--Lagrangian form
of Fokker--Planck--Kolmogorov
equation (cf. Kavvas \cite{Kavvas:2003}). Eq.~(\ref{eqno_1}), when written
in the probability density space, takes the form (Sett et al. \cite{Sett:2007a}):
%
%

\begin{eqnarray}
\nonumber
\lefteqn{\displaystyle \frac{\partial P(\sigma_{ij}(t), t)}{\partial t} = - \displaystyle \frac{\partial}{\partial \sigma_{mn}} 
\left[ \left\{\left< \vphantom{\int_{0}^{t}} \eta_{mn}(\sigma_{mn}(t), D_{mnrs}, \epsilon_{rs}(t))\right> \right. \right.} \\
\nonumber
&+& \left. \left. \int_{0}^{t} d\tau Cov_0 \left[\displaystyle \frac{\partial \eta_{mn}(\sigma_{mn}(t), D_{mnrs}, 
\epsilon_{rs}(t))} {\partial \sigma_{ab}}; \right.  \right. \right. \\
\nonumber
& & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \left. \left. \left. \eta_{ab} (\sigma_{ab}(t-\tau), D_{abcd}, \epsilon_{cd}(t-\tau) 
\vphantom{\int_{0}^{t}} \right] \right \} P(\sigma_{ij}(t),t) \right] \\
\nonumber
&+& \displaystyle \frac{\partial^2}{\partial \sigma_{mn} \partial \sigma_{ab}} \left[  \int_{0}^{t} d\tau Cov_0 \left[ 
\vphantom{\int_{0}^{t}} \eta_{mn}(\sigma_{mn}(t), D_{mnrs}, \epsilon_{rs}(t)); \right. \right.  \\
\label{eqno_2}
& & \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \left. \left.  \eta_{ab} (\sigma_{ab}(t-\tau), D_{abcd}, \epsilon_{cd}(t-\tau)) 
\vphantom{\int_{0}^{t}} \right] \vphantom{\int_{0}^{t}}  P(\sigma_{ij}(t),t) \right]
\end{eqnarray}

%
\noindent  where,  $P(\sigma(t), t)$ is the probability density
of  stress, $\left< \cdot \right>$ is the expectation operator,
and $Cov_0 \left[ \cdot \right]$ is the time-ordered covariance
operator.  One may also note that in Eq.~(\ref{eqno_2}), $t$ is
the    pseudo-time    of   the   constitutive   rate   equation
(Eq.~(\ref{eqno_1}))  and  $\eta_{ij}$  is  a  random  operator
tensor,  which  is a function of stress tensor ($\sigma_{ij}$),
material  properties  tensor  ($D_{ijkl}$),  and  strain tensor
($\epsilon_{kl}$).  Eq.~(\ref{eqno_2})  can  be  written  in  a
compact form as follows:

\begin{equation}
\label{eqno_3}
\displaystyle \frac{\partial P(\sigma_{ij}, t)}{\partial t} =
- \displaystyle \frac{\partial}{\partial \sigma_{mn}} \left[ N_{{(1)}_{mn}} P(\sigma_{ij}, t)
- \displaystyle \frac{\partial}{\partial \sigma_{ab}} \left\{ N_{{(2)}_{mnab}} P(\sigma_{ij},t) 
\right\} \right]
\end{equation}

\noindent  
where,  $N_{(1)}$  and  $N_{(2)}$  are advection and
diffusion   coefficients   respectively. 
%
With  appropriate  initial  and  boundary
conditions,  and  given the second-order statistics of material
properties,  Eq.~(\ref{eqno_3})  can be solved for evolution of
probability  density  of  stress  response  with second-order
accuracy, following any particular constitutive law. 
%
The   advection  and
diffusion  coefficients  are function of statistical properties
of  material  parameters  and  strain  as  well  as the type of
constitutive  model. For example, following the general derivation in the Appendix of
Jeremi\'{c} et al. \cite{Jeremic:2007}, it can be shown that
the advection and diffusion coefficients for 1--D von Mises elastic-perfectly 
plastic shear stress ($\sigma$) versus shear strain ($\epsilon$) constitutive
relationship take the following form:

\begin{eqnarray}
N_{(1)}^{vM} = \left\{\begin{array}{ll}
%
        \displaystyle \frac{d \epsilon}{dt} \left< G \right >   
%
%
         & \mbox{when~elastic} \\
%
           \\
%
         0
%
         & \hspace*{-0.5truecm}\mbox{when~perfectly plastic}
%
           \end{array} \right.
N_{(2)}^{vM} = \left\{\begin{array}{ll}
%
          t \left(\displaystyle \frac{d \epsilon}{dt}\right)^2 Var[G]
%
%
         & \mbox{when~elastic} \\
%
           \\
%
          0
%
          & \hspace*{-1.2truecm}\mbox{when~perfectly plastic}
%
           \end{array} \right.
\label{eqno_4}
\end{eqnarray}

\noindent  
In  Eq.~(\ref{eqno_4}),  $G$ is the elastic shear modulus,  $\epsilon$  is  the 
shear strain, $\left< \cdot
\right >$ is the expectation operator, and $Var\left[ \cdot
\right]$  is the variance operator. The transition from elastic
to    elastic--plastic    can    be   controlled using mean yield criteria (Sett
\cite{Sett:2007f}) --
shifting from elastic advection and diffusion coefficients to
elastic--plastic advection and diffusion coefficients, when
mean of elastic solution exceeds mean of yield stress. In mathematical
terms, this  translates to:

\begin{eqnarray}
\begin{array}{ll}
\mbox{if}~~~~~~~~~~~ \left<f\right> < 0 \vee \left(\left<f\right> = 0 \wedge d\left<f\right> < 0 \right) ~~~~~& \mbox{then use elastic $N_{(1)}$ and $N_{(2)}$} \\
\mbox{or, if}~~~~~~~~~ \left<f\right> = 0 \vee d\left<f\right> = 0 ~~~~~& \mbox{then use elastic--plastic $N_{(1)}$ and $N_{(2)}$}
\end{array}
\label{eqno_5}
\end{eqnarray}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Consideration for Probabilistic Yielding}
\label{ProbabilisticYielding}

If  the  yield  strength  of  the  material  is  uncertain,  then there  is  always a
possibility  of "triggering" the  elastic--plastic  advection  and diffusion coefficients 
in  the  pre--yield  elastic region and vice versa. Mean yield criteria
(Eq.~(\ref{eqno_5})),  however,  fails  to  account for those possibilities. One
possible  solution  could  be to assign probability weights, based on cumulative
density  function  (CDF)  of  yield  strength  ($\Sigma_y$),  to the elastic and
plastic     advection    and    diffusion    coefficients    and    solve    the
Fokker--Planck--Kolmogorov  equation  corresponding  to  that  weighted  average
advection      and     diffusion  coefficients  to  obtain  full  probabilistic
elastic-plastic response. Mathematically, the equivalent (weighted) advection and
diffusion  coefficients  ($N_{(1)}^{eq}$  and  $N_{(2)}^{eq}$) can be written as
(Jeremi\'{c} and Sett \cite{Jeremic:2009b}):
%
\begin{eqnarray}
\begin{array}{c}
N_{(1)}^{eq} (\sigma) 
= 
(1 - P[\Sigma_y \leq \sigma]) N_{(1)}^{el} + P[\Sigma_y \leq \sigma] N_{(1)}^{pl} \\
\label{eqno_PY02}
N_{(2)}^{eq} (\sigma) 
= 
(1 - P[\Sigma_y \leq \sigma]) N_{(2)}^{el} + P[\Sigma_y \leq \sigma] N_{(2)}^{pl}
\end{array}
\label{eqno_6}
\end{eqnarray}
%
where   $(1  -  P[\Sigma_y  \leq \sigma])$ represents the probability of material
being  elastic,  while  $P[\Sigma_y  \leq \sigma]$ represents the probability of
material being elastic--plastic. The superscripts  $\cdot^{el}$  and  $\cdot^{pl}$  
on  the  advection and diffusion coefficients  refer to pre--yield elastic region 
and  post-yield elastic-plastic  region. As an example, for von Mises elastic--perfectly
plastic material (elastic and elastic--perfectly plastic advection and diffusion
coefficients  are  given  by  Eq.~(\ref{eqno_4}))  with uncertain yield strength
($\Sigma_y$),    the    equivalent    advection   and   diffusion   coefficients
($N_{(1)}^{{eq}^{vM}}$ and $N_{(2)}^{{eq}^{vM}}$) would become:
%
\begin{eqnarray}
\begin{array}{l}
N_{(1)}^{{eq}^{vM}} (\sigma) = (1 - P[\Sigma_y \leq \sigma]) \displaystyle \frac{d \epsilon}{dt} \left< G \right >  \\
N_{(2)}^{{eq}^{vM}} (\sigma) = (1 - P[\Sigma_y \leq \sigma]) t \left( \displaystyle \frac{d \epsilon}{dt}\right)^2 Var[G] 
\end{array}
\label{eqno_7}
\end{eqnarray}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Probabilistic Constitutive Simulation}

Following the above described framework, any plasticity theory 
based constitutive model can be written in probability
density space and the resulting FPK partial differential equation (PDE) can be solved,
using appropriate initial and boundary conditions, to obtain the constitutive 
response probabilistically.
%
In this paper, the FPK PDE is
solved numerically using method 
of lines - by semidiscretizing the PDE
in the stress domain using finite difference method 
and solving the resulting series of ordinary
differential equations incrementally.

A  typical,  1--D  shear  stress versus shear strain
constitutive  behavior  of von Mises elastic--perfectly plastic
material  with  uncertain  shear  modulus  and  uncertain yield
(shear) strength                is               shown               in
Figure~\ref{figure:contourPDF_and_mean_and_SD_monotonic}.  
%
%
\begin{figure}[!htbp]
\begin{center}
\includegraphics[width=0.75\textwidth]{contourPDF_and_mean_and_SD_monotonic-Ed02.eps}
\caption{Contour of evolutionary probability density function (PDF) of shear stress, when von Mises elastic--perfectly
plastic material with uncertain shear modulus and uncertain yield (shear) strength was monotonically strained; The mean, the mode, and 
the mean$\pm$standard deviation of shear stress -- obtained  by  integrating  the evolutionary PDF of shear stress -- are also shown.}
%{\bf Mode response should be here too!}}
\label{figure:contourPDF_and_mean_and_SD_monotonic}
\end{center}
\end{figure} 
%
%
%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%
% BJ
%
%\noindent
%\rule[2mm]{\textwidth}{2mm}
%
%NOTE:  Can  you  add  the  mode (as you mention it later in the text, as well as
%change the dark colors to something a bit lighter...)
%
%\noindent
%\rule[2mm]{\textwidth}{2mm}
% BJ
%
%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%
%
In  particular,  Figure~\ref{figure:contourPDF_and_mean_and_SD_monotonic}  shows
the  contour  of evolutionary probability density function (PDF) of shear stress
with shear strain, as the material is  strained monotonically up to 1.026\%.
%
The evolutionary PDF was obtained as the solution of FPK
equation   (Eq.~(\ref{eqno_3}))   with   equivalent   von  Mises
elastic--perfectly plastic advection and diffusion coefficients
(Eq.~(\ref{eqno_7})).
%
Shear modulus, with normal distribution, having a mean value of 61.4 MPa and
a coefficient of variation of 38.5\% was assumed in this simulation. 
%
Shear
strength was also assumed to follow normal distribution with a mean of 0.2 MPa and
a standard deviation of 0.14 MPa. 
%
Shear modulus and shear strength were assumed independent of each other.
%
The assumed values are typical for clay material under undrained condition.
%
In Figure~\ref{figure:contourPDF_and_mean_and_SD_monotonic}, the mean and
mean$\pm$standard deviation of shear stress, obtained by integrating the
PDF of shear stress by standard techniques, are also shown. It is interesting
to observe, in Figure~\ref{figure:contourPDF_and_mean_and_SD_monotonic}, the non-Gaussian
nature of the evolutionary PDF of shear stress. In other words, the mean shear stress
differs  from  the  most  probable  (mode) shear stress. 
%{\bf Yes, show mode in Figure 1.}



In addition  to   quantifying   the   uncertainties  in  predicted
constitutive    behavior,   the   probabilistic   approach   to
elasto--plasticity could  significantly influence material modeling in general.
%
As  can  be  observed  from  Figure~\ref{figure:mean_cyclic}, with probabilistic
approach to material modeling, even
the  simplest  elastic--perfectly  plastic material model
captures  some  of the features of (geo-)materials'
behavior,  for example,  modulus reduction  and  modulus degradation 
as material undergoes  cyclic loading.
%
In particular, Figures~\ref{figure:mean_cyclic}(a) and (b) show the mean shear stress behaviors that are
obtained   from   evolutionary   PDF   of   shear   stress   when   von  Mises,
elastic-perfectly  plastic  material (same as used for monotonic simulation) is
strained  cyclically.  Figure~\ref{figure:mean_cyclic}(a) shows an
unsymmetric loading--unloading--reloading            cycle  while 
Figure~\ref{figure:mean_cyclic}(b)
shows a symmetric loading--unloading--reloading cycle.
%
Deterministic simulation of modulus reduction and modulus degradation would require
(far    more)    advanced   constitutive   models   than   the   (probabilistic)
elastic--perfectly plastic
model used here.
%
The  influence  of  hardening  rules  on  probabilistic
approach  to  material  modeling  was  discussed  by  Sett  and
Jeremi\'{c} \cite{Sett:2010a}.

\begin{center}
\begin{figure}[!htbp]
\includegraphics[width=0.47\textwidth]{mean_partcyclic-Ed.eps}
\hspace*{0.15truecm}
\includegraphics[width=0.47\textwidth]{mean_fullcyclic-Ed.eps}
\hspace*{2.75truecm}(a)\hspace*{7.0truecm}(b)
\caption{Mean shear stress, obtained by integrating the evolutionary probability density function (solution of FPK equation (Eq.~(\ref{eqno_3}))), when
von Mises elastic--perfectly plastic material with uncertain shear modulus and uncertain yield (shear) strength was cyclically strained: (a)
loading--partial unloading--reloading cycle and (b) loading--complete unloading--reloading cycle; Two complete cycles are shown for the later}
\label{figure:mean_cyclic}
\end{figure} 
\end{center}

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

\section{Formulation of Stochastic Elastic--Plastic Finite Element}
\label{SFEM}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Governing Equation}

The governing partial differential equation in mechanics can be mathematically written,
combining equilibrium, strain compatibility, and constitutive equations, as:
%
%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%
% BJ
%\noindent
%\rule[2mm]{\textwidth}{2mm}
%
%NOTE:   This is all for 1D (since they are all functions of $x$ and not $x_i$ of
%${\bf x}$ ? How do we represent this for a general 3D case?
%
%Can we simply replace all $x$ with  $x_i$ or  ${\bf x}$?
% 
%\noindent
%\rule[2mm]{\textwidth}{2mm}
% BJ
%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%

\begin{equation}
\label{eqno_Dis01}
\Xi({\bf x}) u({\bf x}) = \phi({\bf x})
\end{equation}

\noindent where $\Xi({\bf x})$ is a linear/non-linear differential operator, $\phi({\bf x})$ is the
external force and $u({\bf x})$ is the response.
%
If the material properties are uncertain (and modeled as random field), $\Xi$ in Eq.~(\ref{eqno_Dis01}) becomes a stochastic
linear/non-linear differential operator and as a result, the response, $u$ becomes a random field.
%
Splitting  the stochastic operator ($\Xi$) into a deterministic part ($L$) and a
random  part  ($\Pi$),  whose  coefficients are zero-mean random fields, one can
write Eq.~(\ref{eqno_Dis01}) as:

\begin{equation}
\label{eqno_Dis02}
\left[L ({\bf x}) + \Pi({\bf x},\theta) \right] u({\bf x}, \theta) = \phi ({\bf x})
\end{equation}

%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%
% BJ
%\noindent
%\rule[2mm]{\textwidth}{2mm}
%
%
%NOTE:   This is all for 1D (since they are all functions of $x$ and not $x_i$ of
%${\bf x}$ ? How do we represent this for a general 3D case?
%
%Can we simply replace all $x$ with  $x_i$ or  ${\bf x}$?
% 
%\noindent
%\rule[2mm]{\textwidth}{2mm}
% BJ
%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%

\noindent In the above equation (Eq.~(\ref{eqno_Dis02})), $\theta$ is introduced to 
denote randomness in the variables. For example, the response $u({\bf x}, \theta)$ is
a spatial random function, where the probability distribution of $u$ changes as a function
of the location, ${\bf x}$, in the space continuum.
%
Further,  if one splits the input material properties
random field into a deterministic trend and zero-mean random (uncertain)
residual             about             trend,  $\mathbb{D}  ({\bf x}, \theta)  =
\hat{\mathbb{D}}({\bf x})  +  \mathbb{R}  ({\bf x}, \theta)$, one can re-write
Eq.~(\ref{eqno_Dis02}) as:

\begin{equation}
\label{eqno_Dis03}
\left[L_1 ({\bf x}) \hat{\mathbb{D}}({\bf x}) + L_2 ({\bf x}) \mathbb{R} ({\bf x},\theta) \right] u({\bf x}, \theta) = \phi ({\bf x})
\end{equation}

\noindent  where  $L_1  ({\bf x})$  and  $L_2  ({\bf x})$ are deterministic
differential   operators. 

%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%
% BJ
%\noindent
%\rule[2mm]{\textwidth}{2mm}
%
%NOTE: perhaps  just  a  tiny  intro  on  random  fields,  and  the  meaning of $\theta$
%(stochastic dimension...)
% 
%\noindent
%\rule[2mm]{\textwidth}{2mm}
% BJ
%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Spatial and Stochastic Discretization}

Spectral approach to stochastic
finite element formulation necessities discretization of 
Eq.~(\ref{eqno_Dis03}) in both stochastic and spatial dimensions, as follows:

\begin{enumerate}

\item Karhunen--Lo\`{e}ve (KL) expansion (Karhunen \cite{Karhunen:1947};
Lo\`{e}ve \cite{Loeve:1948};  Ghanem and Spanos  \cite{book:Ghanem}) 
is used to discretize the zero-mean fluctuating part
of the input material properties random field ($\mathbb{R} ({\bf x}, \theta)$) 
into finite number of independent basic random variables, 

\begin{equation}
\label{eqno_c1}
\mathbb{R} ({\bf x}, \theta) = \sum_{n=1}^L \sqrt \lambda_n \xi_n (\theta) f_n({\bf x})
\end{equation}

%
where,   $\lambda_n$   and   $f_n({\bf x})$  are  the  eigenvalues  and  eigenvector,
respectively of the covariance kernel of the zero-mean random field ($\mathbb{R}
({\bf x},  \theta)$).  The  zero-mean  random  variables, $\xi_n(\theta)$ are mutually
independent and have unit variances.

\item    Polynomial    chaos    (PC)    expansion   (Wiener \cite{Wiener:1938};      
Ghanem and Spanos \cite{book:Ghanem}) is used to represent any unknown random variable
in terms of known random variables. For example, any unknown random variable,
$\chi(\theta)$ can be expanded, truncating after $P$ terms, in functional (polynomial chaos, $\psi_i\left[\left\{\xi_r\right\}\right]$) of
known random variables, $\xi_r$ and unknown deterministic coefficients, $\gamma_i$, as:

\begin{equation}
\label{eqno_c2}
\chi_j(\theta)=\sum_{i=0}^P\gamma_i^{(j)}\psi_i\left[\left\{\xi_r\right\}\right]
\end{equation}

Using PC expansion, the unknown partially-discretized -- in spatial dimension, using KL expansion -- displacement 
random field ($u({\bf x}, \theta)$) in Eq.~(\ref{eqno_Dis03}) can be further discretized in the stochastic dimension,
 

\begin{eqnarray}
%\label{eqno_c3}
u({\bf x},\theta)
&=&
\sum_{j=1}^L e_j \chi_j(\theta)b_j({\bf x}) 
\nonumber \\
&=& 
\sum_{j=1}^L \sum_{i=0}^P \gamma_i^{(j)} \psi_i[\{\xi_r\}]c_j({\bf x}) 
\nonumber \\
&=&
\sum_{i=0}^P \psi_i[\{\xi_r\}] \sum_{j=1}^L \gamma_i^{(j)}c_j({\bf x}) 
\nonumber \\
&=& 
\sum_{i=0}^P \psi_i[\{\xi_r\}] d_i({\bf x})
\label{eqno_c6}
\end{eqnarray}

\noindent
where  $c_j({\bf x}) = e_j b_j({\bf x})$ and $d_i({\bf x})=\sum_{j=1}^L\gamma_i^{(j)} c_j({\bf x})$

\item Shape function expansion (Zienkiewicz  and  Taylor  \cite{book:Zienkiewicz},  Bathe
\cite{book:Bathe}, Ghanem and Spanos \cite{book:Ghanem}) is used to discretize the spatial component
-- $d_i  ({\bf x})$,  in  the  above polynomial chaos expansion, Eq.~(\ref{eqno_c6}) -- of the
unknown random field:

\begin{equation}
\label{eqno_c7}
d_i({\bf x})=\sum_{n=1}^N d_{ni} l_n({\bf x})
\end{equation}

\noindent where $l_m$({\bf x}) are the shape functions.

\end{enumerate}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Stochastic Finite Elements}

Discretizing Eq.~(\ref{eqno_Dis03}) in both spatial and stochastic dimensions using KL, PC, and shape function expansions, one can write
Eq.~(\ref{eqno_Dis03}) as, 

\begin{eqnarray}
\nonumber
&& \sum_{i=0}^P \sum_{n=1}^N \left[ \vphantom{\sum_{m=1}^{M}} \psi_i[\{\xi_r\}] L_1 ({\bf x}) \hat {\mathbb{D}} ({\bf x}) l_n({\bf x}) ~ + \right. \\
\label{eqno_c8}
&& \left. \ \ \ \ \ \ \ \ \ \ \ \sum_{m=1}^{M} \xi_m (\theta) \psi_i[\{\xi_r\}] ~L_2 ({\bf x}) 
\sqrt \lambda_m f_m({\bf x}) l_n({\bf x}) \right]  d_{ni} = \phi ({\bf x})
\end{eqnarray}


\noindent 
Galerkin     type     procedure     may     be    applied    to
Eq.~(\ref{eqno_c8})  to solve for the unknown coefficients ($d_{mi}$) of
the PC-expansion of the displacement random field and after some algebra, one may write 
Eq.~(\ref{eqno_c8}) as (cf. Ghanem and Spanos \cite{book:Ghanem}),

\begin{equation}
\label{eqno_c9}
\sum_{n=1}^N K'_{mn}d_{ni} +\sum_{n=1}^N \sum_{j=0}^P d_{nj} \sum_{k=1}^M c_{ijk} K''_{mnk} 
= \Phi_m \left< \psi_i[\{\xi_r\}]\right>
\end{equation}


\noindent 
where,   $\Phi_{k}  =  \int_{D}  \phi  ({\bf x})  l_k({\bf x})  d{\bf x}$  and  $c_{ijk}  =  \left
<\xi_k(\theta)  \psi_i[\{\xi_r\}]  \psi_j[\{\xi_r\}]  \right  >$.  In  the above
equation  (Eq.~(\ref{eqno_c9})),  one  may  note  that  the  expected values of the
polynomial  chaoses  ($\psi_i[\{\xi_r\}]$)  and  the products of orthonormal random
variables     and    polynomial    chaoses    ($\xi_k(\theta)    \psi_i[\{\xi_r\}]
\psi_j[\{\xi_r\}]$)   can easily be pre-calculated in closed form (symbolically,
for example using Mathematica \cite{software:Mathematica}).
%
Eq.~(\ref{eqno_c9}) can be written in more familiar matrix form as,

\begin{equation}
\label{eqno_c10}
\bar K \bar{u} = \bar F
\end{equation}

\noindent 
where,  $\bar{u}$  is  the  generalized  displacement  vector,  $\bar  F$ is the
generalized  force  vector,  and  $\bar  K$ is the generalized stiffness matrix,
which    is    composed  of  two  components, namely  the deterministic stiffness
matrix,  $K'$, defined as:

\begin{equation}
\label{eqno_c11}
K'_{nk} = \int_{D} L_1 ({\bf x}) l_n({\bf x})~\hat {\mathbb{D}}~l_k({\bf x}) d{\bf x}
\end{equation}

\noindent 
and, the stochastic stiffness matrix, $K''$, defined as:


\begin{equation}
\label{eqno_c12}
K''_{mnk} = \int_{D} L_2 ({\bf x}) l_n({\bf x}) \left\{\sqrt \lambda_m f_m({\bf x}) \right\} l_k({\bf x}) d{\bf x}
\end{equation}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Force-Residual Form}



For non-linear problems (as attempted in this paper), Eq.~(\ref{eqno_c10})
can be written in force-residual form and solved incrementally. 
One  can  write Eq.~(\ref{eqno_c10}) in force-residual
form as:

\begin{equation}
\label{eqno_1EP}
\bar r(\bar u, \lambda) = 0
\end{equation}

\noindent 
where,  $\bar  u$  are  the  generalized degrees of freedom and $\lambda$ is the
load control  parameter. Differentiating Eq.~(\ref{eqno_1EP}), one can write the rate
form of the force-residual equation (Eq.~(\ref{eqno_1EP})) as,

\begin{equation}
\label{eqno_2EP}
\dot {\bar r} = \bar K \dot {\bar u} - \bar q \dot \lambda = 0
\end{equation}

\noindent where, $\bar K$ (=$\partial r_i/\partial u_j$) is the tangent stiffness matrix and $\bar q$ (= $-\partial r_i/\partial \lambda$) is
the load vector. At regular points in the $\bar u-\lambda$ space, the tangent stiffness matrix is non--singular and 
hence one can solve the force-residual rate equation for $\dot {\bar u}$ as,


%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%
% BJ
%\noindent
%\rule[2mm]{\textwidth}{2mm}
%
%See  if  changes  that I made below do make sense (changes run until the of this
%subsection, that is until section 3.5
% 
%\noindent
%\rule[2mm]{\textwidth}{2mm}
% BJ
%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%

\begin{equation}
\label{eqno_5EP}
\dot {\bar u} = \left({\bar K}^{-1} \bar q \right) \dot \lambda
\;\;\;\;
\mbox{with}
\;\;\;\;
{\bar u}^{'} = \frac{d \bar u}{d \lambda} =  {\bar K}^{-1} \bar q 
\end{equation} 

%\noindent or,
%
%\begin{equation}
%\label{eqno_6EP}
%{\bar u}^{'} = \frac{d \bar u}{d \lambda} =  {\bar K}^{-1} \bar q 
%\end{equation}

\noindent 
The resulting set of nonlinear ODEs (Eq.~(\ref{eqno_5EP})) can be solved numerically
for $\bar u$ 
with  a  set of initial conditions (at $\lambda$ = 0, $\bar u$ =
$\bar  u_0$)  using  either pure  incremental  method (forward Euler
method) or incremental--iterative method (Newton method). This
paper only considers the forward Euler method by which knowing the 
solution of $\bar u$ at the $n^{th}$ step ($\bar u_n$), the solution 
at $(n+1)^{th}$ step ($\bar u_{n+1}$) can be obtained as (using load control),

\begin{equation}
\label{eqno_7EP}
\bar u_{n+1} = \bar u_n + \Delta \lambda_n {\bar u_n}^{'}
\end{equation}

\noindent Hence, the incremental solution ($\Delta \bar u_n$) can be written as,

\begin{equation}
\label{eqno_8EP}
\Delta \bar u_n = \bar u_{n+1} - \bar u_n = \left( {\bar K_n}^{-1} \bar q_n \right) \Delta \lambda_n 
\end{equation}





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Evaluation of Generalized Tangent Stiffness Matrix}

In Eq.~(\ref{eqno_8EP}), the  generalized tangent stiffness matrix
($\bar K$)  needs  to  be  re-evaluated  at  each step because the
constitutive  properties ($\hat
{\mathbb{D}}$  in  Eq.~(\ref{eqno_c11}) and $\sqrt \lambda_m$ in
Eq.~(\ref{eqno_c12})) of the material change as the material
plastifies.
%
The evolution of these material properties, as the material plastifies, are
governed by the constitutive rate equation (Eq.~(\ref{eqno_1})). 
%
At each  integration   point   (Gauss   point), 
the constitutive equation, written in probability density space (the FPKE, refer Eq.~(\ref{eqno_3})),
can be solved once for each KL-space\footnote{Note that KL-spaces are orthogonal to each other} to obtain 
updated $\hat{\mathbb{D}}$ (refer Eq.~(\ref{eqno_c11})) 
and $\sqrt \lambda_m$ (refer Eq.~(\ref{eqno_c12})), and thereby the generalized tangent stiffness matrix
($\bar K$; refer Eq.~(\ref{eqno_8EP})) at each load step. 
%
For example, for 1--D problem, the (known,  random) strain increment at the $({k-1})^{th}$ global
step can be written as,

\begin{equation}
\label{eqno_13EP}
\Delta \epsilon^{k-1} (x, \theta) 
= \sum_{i=0}^P \psi_i[\{\xi_r\}] (\theta) \sum_{n=1}^N \Delta d^{k-1}_{ni} B_n(x)
\end{equation}

\noindent Hence, the advection and diffusion coefficients for the FPKE, assuming 1--D von Mises elastic--perfectly plastic 
material behavior, at the $k^{th}$ global step can be calculated using Eq.~(\ref{eqno_7}) as follows:

\begin{itemize}

\item[] For the zeroth KL-space (mean space):
%
\vspace*{-0.1truecm}
%
\begin{eqnarray}
\begin{array}{l}
{N_{(1)}^{{eq}^{vM}}}^k (\sigma) 
= 
(1/\Delta t)~\hat{\mathbb{D}}~\left(1 - P[\Sigma_y \leq \sigma]\right) \\
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \sum_{i=0}^P \left <\vphantom{\sum_{i=0}^P} \psi_i[\{\xi_r\}] \right > \sum_{n=1}^N \Delta d^{k-1}_{ni} B_n  
\\
{N_{(2)}^{{eq}^{vM}}}^k (\sigma) 
= 
t~(1/\Delta t)^2~\hat{\mathbb{D}}^2~\left(1 - P[\Sigma_y \leq \sigma]\right) \\
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \sum_{i=0}^P Var \left[\vphantom{\sum_{i=0}^P} \psi_i[\{\xi_r\}] \right] 
\left(\sum_{n=1}^N \Delta d^{k-1}_{ni} B_n \right)^2 
\end{array}
\label{ADcoeffs01}
\end{eqnarray}

\item[] For any other KL-spaces:
%
\vspace*{-0.1truecm}
%
\begin{eqnarray}
\begin{array}{l}
{N_{(1)}^{{eq}^{vM}}}^k (\sigma) 
= 
(1/\Delta t)~\sqrt \lambda~f ~\left(1 - P[\Sigma_y \leq \sigma]\right) \\
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \sum_{i=0}^P \left <\vphantom{\sum_{i=0}^P} \xi ~\psi_i[\{\xi_r\}] \right > 
\sum_{n=1}^N \Delta d^{k-1}_{ni} B_n 
\\
{N_{(2)}^{{eq}^{vM}}}^k (\sigma) 
= 
t~(1/\Delta t)^2~\lambda~f^2~\left(1 - P[\Sigma_y \leq \sigma]\right)  \\
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \sum_{i=0}^P Var \left[\vphantom{\sum_{i=0}^P} \xi ~\psi_i[\{\xi_r\}] \right] 
\left(\sum_{n=1}^N \Delta d^{k-1}_{ni} B_n \right)^2
\end{array}
\label{ADcoeffs02}
\end{eqnarray}

\end{itemize}

\noindent In Eqs.~(\ref{eqno_13EP})-(\ref{ADcoeffs02}), $B$ is the derivative of shape function. 
%
One   may   note   that   the   means   and   the variances   of   the polynomial   chaoses
($\psi_i[\{\xi_r\}]$)  and  the means and the variances of the products of the orthonormal
Gaussian  random  variables  ($\xi_m$,  where  the  subscript  $m$  denotes  the
KL-space)   and   the polynomial   chaoses   ($\psi_i[\{\xi_r\}]$)   can easily be
pre-calculated symbolically using Mathematica \cite{software:Mathematica}.
% 


Assuming  any  value  of pseudo-time increment $\Delta  t$,  the FPKEs corresponding to each KL-space
(including  the  mean  space)  can  be  solved  at  $t = \Delta t$ to obtain the
corresponding  probability  density function of stress for the strain increment
given by Eq.~(\ref{eqno_13EP}).
%
% out as irrelevant at this point
% For  advanced global level solution schemes, such as Newton-type
% incremental-iterative   scheme,  these  probability  densities  of
% stress at  each KL-space (including the mean space) could be used to compute the internal
% forces.  However,  in  this paper, only pure incremental
% solution  scheme  is  considered. 
The  explicit incremental solution
scheme used here   needs  information  on  the  new tangent material
properties  at  each KL-space (for example, $\hat {\mathbb{D}}$
in  Eq.~(\ref{eqno_c11}) and $\sqrt \lambda_m$  in
Eq.~(\ref{eqno_c12}))  to increment forward. 
%
%In this paper, the  new  tangent  material  properties  at each KL-space is calculated from the
%corresponding  PDF of evolutionary stress at each KL-space.
%
To this end, the    constitutive    rate    equation
(Eq.~(\ref{eqno_1})), after taking expectation on both sides, 
can be written as (Kavvas \cite{Kavvas:2003}; Sett et al. \cite{Sett:2007a}):

\begin{eqnarray}
\nonumber
&& \frac{d \left< \sigma_{ij}(t)\right>}{d t}
=
\left<  D_{ijkl}(\sigma_{ij},  t)  \frac{d  \epsilon_{kl}}{d  t}  \right  >  + \\ 
\nonumber
&&  \ \ \ \ \int_{0}^t ds \;
      Cov_0\left[\frac{\partial}{\partial \sigma_{mn}} \{D_{ijkl}({\sigma_{ij}}_{(t)}, t) \displaystyle \frac{d \epsilon_{kl}}{d t}\}; 
      D_{mnpq}({\sigma_{mn}}_{(t-s)}, t-s) \frac{d \epsilon_{pq}}{d t} \right] \\
\label{eqno_14EPa}
\end{eqnarray}

\noindent 
If one assumes $\sigma (t)$ to be a $\delta$-correlated function of (pseudo) time, then the covariance term on the r.h.s of  
Eq.~(\ref{eqno_14EPa}) becomes variance and hence, Eq.~(\ref{eqno_14EPa}) can be re-written in incremental form as:

\begin{equation}
\label{eqno_14EPb}
\frac{\Delta \left< \sigma_{ij}(t)\right>}{\Delta t}
=
\left< D_{ijkl}(\sigma_{ij}, t) \frac{\Delta \epsilon_{kl}}{\Delta t} \right > +
Var\left[D_{ijkl}(\sigma_{ij}, t) \frac{\Delta \epsilon_{kl}}{\Delta t} \right]
\end{equation}

\noindent It may be noted that if relatively large computational step size is used,
then $\delta$-correlation is a fairly good assumption for processes
with finite correlation lengths as most of the information are in the variances
of those processes. 
%
%
Eq.~(\ref{eqno_14EPb}) relates the statistics of the tangent stiffness
with the rate of change of mean stress, which is known, at each global load step
increment, from the solution of FPKE.  
%
%However, that itself doesn't allow for direct solution
%of the unknowns on the r.h.s, unless the l.h.s (rate of change of mean stress) is known at several intermediate
%steps within the global load step increment, in which case the unknowns on the r.h.s of 
%Eq.~(\ref{eqno_14EPb}) can be approximately solved  by using a least square procedure 
%(e.g. Levenberg-Marquardt technique, Levenberg \cite{Levenberg:1944},
%Marquardt  \cite{Marquardt:1963}).
%
However, that itself doesn't allow for direct solution of the
unknowns on the r.h.s., unless the l.h.s. (rate of change of mean stress) is
known at several intermediate steps within the global load step increment.
%
In such cases, the unknowns on the r.h.s. of Eq.~(\ref{eqno_14EPb}) can be approximately
solved by using a least square procedure e.g., Levenberg-Marquardt 
technique (Levenberg \cite{Levenberg:1944}, Marquardt  \cite{Marquardt:1963}).
%
Hence, in this paper, at each Gauss
integration point, the FPKE is
solved in few sub-steps within each  global step increment.
%
%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%
% BJ
%\noindent
%\rule[2mm]{\textwidth}{2mm}
%
%Might be good to explain a sub-step procedure in constitutive part!
%
%and then reference here!
% 
%\noindent
%\rule[2mm]{\textwidth}{2mm}
% BJ
%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%
%
By      knowing   the  l.h.s.  of  Eq.~(\ref{eqno_14EPb})  at those
sub-steps and noting that the incremental
strain   is   given  by  Eq.~(\ref{eqno_13EP}),  the  evolved  tangent  material
properties  of  each  space  are approximately calculated.  For  example,  for  1-D  problem, for zeroth
KL-space (mean space), Eq.~(\ref{eqno_14EPb}) simplifies to:


\begin{equation}
\label{eqno_15EP}
{\left . \frac{\Delta \left< \sigma (t)\right>}{\Delta t}\right |}_{zero^{th} KL-space}
=
\hat{\mathbb{D}} (t) \left< \frac{\Delta \epsilon }{\Delta t} \right > +
{\hat{\mathbb{D}}^2 (t)} Var\left[\frac{\Delta \epsilon }{\Delta t} \right]
\end{equation}

%
%
%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%
% BJ
%\noindent
%\rule[2mm]{\textwidth}{2mm}
%
%Might  be  good to explain  the  choice  of quadratic function instead of linear
%or  exponential  (we  have  discussed  this way back and you might have notes on
%that...
%
% 
%\noindent
%\rule[2mm]{\textwidth}{2mm}
% BJ
%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%

\noindent   
Now, if   one   assumes   that, within  the  global load step, the   mean   tangent  stiffness 
($\hat{\mathbb{D}}  (t)$) evolves quadratically \footnote{Linear evolution
was also tried. However, quadratic evolution was preferred due to global step-size issue, especially
at the transition region between elasticity and plasticity.} as $ \hat{\mathbb{D}}
(t)~=~a_1~+~a_2~t^2   $,   where  $a_1$  and  $a_2$  are  unknown  deterministic
coefficients, then one can write Eq.~(\ref{eqno_15EP}) as,

\begin{eqnarray}
\nonumber
{\left. \frac{\Delta \left< \sigma (t)\right>}{\Delta t} \right |}_{zero^{th} KL-space}
&=&  \left ( a_1 \left< \frac{\Delta \epsilon }{\Delta t} \right > + a_1^2 Var\left[\frac{\Delta \epsilon }{\Delta t} \right]\right) + \\
\nonumber
&&
\left ( a_2 \left< \frac{\Delta \epsilon }{\Delta t} \right > + 2 a_1 a_2 Var\left[\frac{\Delta \epsilon }{\Delta t} \right] \right) t^2 + \\
\label{eqno_16EP}
&& a_2^2 Var\left[\frac{\Delta \epsilon }{\Delta t} \right] t^4
\end{eqnarray}

\noindent 
Using  Eq.~(\ref{eqno_13EP})  and  the same $\Delta t$ that was used for solving
the  FPKE, the mean and variance terms on the r.h.s. of Eq.~(\ref{eqno_16EP}) can
be  easily  evaluated and hence, knowing the l.h.s. of Eq.~(\ref{eqno_16EP}) over
few  sub-steps,  the unknown deterministic coefficients ($a_1$ and $a_2$) can be
determined    using        Levenberg-Marquardt        technique.  Hence,    the  
evolved  mean  tangent  material property ($\hat{\mathbb{D}} (t)$) at the end of
global load step can be evaluated as $ \hat{\mathbb{D}} (t)~=~a_1~+~a_2~t^2 $ by
substituting $t = \Delta t$.

Similarly, for any other KL-spaces, for 1-D problem, Eq.~(\ref{eqno_14EPb}) simplifies to:

\begin{equation}
\label{eqno_17EP}
{\left. \frac{\Delta \left< \sigma (t)\right>}{\Delta t} \right |}_{non-zero KL-space}
=
\sqrt \lambda (t) f \left< \xi \frac{\Delta \epsilon }{\Delta t} \right > +
\lambda (t) f^2 Var\left[\xi \frac{\Delta \epsilon }{\Delta t} \right]
\end{equation}


%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%
% BJ
%\noindent
%\rule[2mm]{\textwidth}{2mm}
%
%Again,might be  good to explain  the  choice  of quadratic function instead of linear or exponential (we
%have discussed this way back and you might have notes on that...
%
% 
%\noindent
%\rule[2mm]{\textwidth}{2mm}
% BJ
%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%

\noindent Again, within the global load step, one may assume that the tangent material property at the $m^{th}$ KL-space ($\sqrt \lambda_m (t)$) evolves
quadratically\footnote{As for the evolution of mean stiffness at the zero$^{th}$ KL-space, here also, linear evolution
was tried. However, quadratic evolution was preferred due to global step-size issue, especially
at the transition region between elasticity and plasticity.} as,
$ \sqrt \lambda (t)~=~a_3~+~a_4~t^2 $, where $a_3$ and $a_4$ are unknown deterministic coefficients, and write  
Eq.~(\ref{eqno_17EP}) as,

\begin{eqnarray}
\nonumber
{\left. \frac{\Delta \left< \sigma (t)\right>}{\Delta t} \right |}_{non-zero KL-space}
&=&  \left ( f a_3 \left< \xi \frac{\Delta \epsilon }{\Delta t} \right > + f^2 a_3^2 Var\left[\xi \frac{\Delta \epsilon }{\Delta t} \right]\right)  + \\
\nonumber
&&
\left ( f a_4 \left< \xi \frac{\Delta \epsilon }{\Delta t} \right > + 2 f^2 a_3 a_4 Var\left[\xi \frac{\Delta \epsilon }{\Delta t} \right] \right) t^2 + \\
\label{eqno_18EP}
&&
f^2 a_4^2 Var\left[\xi \frac{\Delta \epsilon }{\Delta t} \right] t^4
\end{eqnarray}

\noindent  Hence, following  the  same procedure, as described above for the zeroth KL-space, the  evolved  tangent  material  property  at  the
$m^{th}$  KL-space ($\sqrt \lambda_m (t))$ at the end of global load step can be
estimated   from  Eq.~(\ref{eqno_18EP})       using        Levenberg-Marquardt      technique.

\subsection{A Remark on Non-Gaussian Material Properties}

The formulation presented above utilizes the properties of Gaussian distribution.
%
However, the formulation can be generalized to include
non-Gaussian (random field) material properties too. To this end, the usual technique
of representing a non-Gaussian random variable in terms of
Gaussian random variables through non-linear transformation, e.g., Rosenblatt transformation 
(Rosenblatt \cite{Rosenblatt:1952}, Melchers \cite{book:Melchers}, Grigoriu \cite{book:Grigoriu}) can be used.
%
Alternately, as recently used by several researchers, classical polynomial chaos expansion (Ghanem \cite{Ghanem:nonGaussianSFEM}, 
Sakamoto and Ghanem \cite{Sakamoto:2002}) or  generalized polynomial chaos expansion (Xiu and Karniadakis 
\cite{Karniadakis:2002, Karniadakis:2003a}, Foo et al. \cite{Karniadakis:2007}, Lucor et al. \cite{Karniadakis:2004}) 
may also be utilized to represent non-Gaussian random fields.
%
It may be noted that the FPK equation, used to integrate the constitutive
equation at integration points, is general enough to treat both Gaussian
and non-Gaussian random variables. See Sett and Jeremi\'{c} \cite{Sett:2010a}, where
FPK equation was applied to probabilistically simulate cyclic
constitutive behavior of a von Mises material, whose shear strength
was assumed to have a Weibull distribution.
%


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

\section{SIMULATION RESULTS AND DISCUSSIONS}


%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%
% BJ
%\noindent
%\rule[2mm]{\textwidth}{2mm}
%
%Check  numbers  (material  properties) to be consistent with previous mention of
%these same values (seems like there is an error...).
% 
%\noindent
%\rule[2mm]{\textwidth}{2mm}
% BJ
%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%


In  this section, the spectral stochastic elastic-plastic finite element method, proposed
in    Section~\ref{SFEM},   is   applied   to simulate the probabilistic behavior of a  
1--D soil column subjected to monotonic shear loading.
%
The schematic of the soil column is shown in Figure~\ref{figure:schematic}(a).
%
%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%
% BJ
%\noindent
%\rule[2mm]{\textwidth}{2mm}
%
%Would be good if you can  add  initial sigmas for cor length for 0.1 and 1.0 in the same plot (those
%from  Figure  9  and  also spend one paragraph explaining how you developed this
%initial profile
%
%explain how you came up with initial distribution of data properties!!
%
% 
%\noindent
%\rule[2mm]{\textwidth}{2mm}
% BJ
%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%
%
%
\begin{figure}[!htbp]
\begin{center}
%\hspace*{-0.5truecm}
\includegraphics[height=0.37\textwidth]{ProblemSchematic03.eps} 
\hspace*{-0.25truecm}
\includegraphics[height=0.37\textwidth]{mean_Profile_Stiffness.ps.eps}
\hspace*{-0.4truecm}
\includegraphics[height=0.37\textwidth]{SD_Profile_Stiffness.ps.eps}
\hspace*{0.6truecm}
\includegraphics[height=0.37\textwidth]{mean_Profile_Strength.ps.eps}
\hspace*{-0.4truecm}
\includegraphics[height=0.37\textwidth]{SD_Profile_Strength.ps.eps}
\end{center}
~~~(a) \hspace*{7.1truecm} (b)
%\begin{center}
\caption{(a) Schematic of the 1-D stochastic soil shear beam; (b) Soil profile, in terms of means and standard deviations,
of shear modulus and shear strength along the depth. Note that the standard deviation profiles are approximately calculated 
from the first two eigen modes of respective covariance kernels}
\label{figure:schematic}
%\end{center}
\end{figure}
%
The  shear modulus ($G$) of soil is assumed to be a Gaussian random field with a
mean  of  61.4  MPa,  
%({\bf  CHECK}),  
a  standard deviation of 23.7 $MPa$, 
%({\bf CHECK}), 
and a correlation length of 0.1 m (covariance function is assumed to be
exponential).
%
The shear strength ($c_u$) of the soil is also assumed to be a Gaussian random field
with a mean of 0.2 MPa, a standard deviation of 0.14 MPa, and a correlation length
of 0.1 m (again, covariance function is assumed to be exponential).
%
The shear modulus and shear strength are assumed independent of each other.
%
%Please note that the above assumed values for soil are typical of clay.
%
%
Although the above material parameters are assumed, the
statistical properties of the soil strength and stiffness
parameters can be easily estimated by statistically
analyzing commonly used in-situ test data e.g., standard penetration test (SPT) $N$-value
or cone penetration test (CPT) $q_T$ value (cf. Fenton \cite{Fenton:1999A,Fenton:1999B}).
%
The soil profile, in terms of means and standard deviations of the shear modulus and shear strength,
based on the above assumed statistical properties are shown in Figure~\ref{figure:schematic}(b).
%
The standard deviations of shear modulus and shear strength, as shown in Figure~\ref{figure:schematic}(b), are approximately
calculated at Gauss integration points using first two eigen modes of their respective covariance kernels.
%
It is assumed that the soil follows von Mises yield
criteria, which is a common assumption for simulation  of
undrained clay behavior.
%
Further,  it  is  assumed  that  after  yielding the soil behaves as a perfectly
plastic  material.  Typical  1-D  shear  stress versus shear strain constitutive
behavior of von Mises elastic--perfectly plastic material under uncertainty were
shown     in     Figures~\ref{figure:contourPDF_and_mean_and_SD_monotonic}    --
\ref{figure:mean_cyclic}.


The  input  soil properties random fields are first discretized in both spatial
and  stochastic dimension using KL-expansion (Eq.~(\ref{eqno_c1})). To this end,
the  eigenvalues  and  eigenvectors  of their respective covariance kernels are
obtained  numerically, using  finite  element  technique.
%
In this context, it may be noted that for exponential covariance kernel,
as assumed in this example, the corresponding eigenvalue problem 
has a closed form solution and hence, a numerical treatment may be avoided.
% 
However, except for very few standard covariance kernels,
numerical solution is necessary for any arbitrary covariance function.
%
For details of finite
element  formulation of KL eigenvalue problem, the readers may refer to the book
by  Ghanem  and  Spanos  \cite{book:Ghanem}.  In this study, open source FORTRAN
library,  LAPACK  (Anderson et al \cite{book:LAPACK}),  is  used  for KL eigen analysis. For the
stochastic  simulation  of  the above soil column, only the first two eigenmodes
for both shear modulus and shear strength are considered.
%
In  each loading step, the stochastic displacement (solution) of the soil column
is   assembled   using   PC-expansion  (Eq.~(\ref{eqno_c2})).  In  this  study,
Hermite polynomial chaoses up to second order are considered.
%
The unknown deterministic coefficients of the PC-expansion are then obtained by
solving the spectral stochastic finite element system of equations (Eq.~(\ref{eqno_c9})).
UMFPACK (Davis \cite{manual:UMFPACK}), a set of routines for solving a sparse
linear system using multifrontal method, is used for this purpose. 
%
At  the  constitutive  level, FPKE (Eq.~(\ref{eqno_3})) corresponding to the von
Mises  elastic perfectly plastic model (the advection and diffusion coefficients
are  given  by  Eqs.~(\ref{ADcoeffs01})-(\ref{ADcoeffs02}))  is  solved at each
KL space  (including  the  mean  space)  to  obtain  the  statistics of internal
stresses as the soil evolves.
%
To this end, the FPK partial differential equation is first semi-discretized in
the stress domain on a uniform grid by central differences to obtain a series of
ODE.  These  ODEs  are then solved simultaneously, after incorporating boundary
conditions,    using    a    standard    ODE   solver SUNDIALS (Hindmarsh et al.
\cite{manual:SUNDIALS}),  which  utilizes  ADAMS method and
functional  iteration. 
%
Readily available
open source code LMFIT (Mor{\'{e}} et al. \cite{LMFIT}) for Levenberg-Marquardt type least square 
regression analysis is used in estimating the evolution of statistical
properties of material parameters from the evolutionary statistical
properties of internal stress (refer to Eqs.~(\ref{eqno_15EP})-(\ref{eqno_18EP})).

\begin{figure}[!htbp]
\begin{center}
%\includegraphics[width=0.48\textwidth]{mean.eps}
%\hspace*{0.1truecm}
%\includegraphics[width=0.48\textwidth]{sd.eps} \\
%\hspace*{-0.0truecm}(a)~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~(b)\\
%\vspace*{0.5truecm}
\includegraphics[width=0.65\textwidth]{meanPMsd.ps.eps} \\
%(c)
\caption{Evolution of mean and mean$\pm$standard deviation of displacement at the top of the 12 m tall 
stochastic soil column, when a deterministic static shear force is applied (varied from 0 to 0.3 kN) at top}
\label{figure:mean_and_sd_of_displacements}
\end{center}
\end{figure}

Figure~\ref{figure:mean_and_sd_of_displacements} shows the evolutions of mean and mean$\pm$standard deviation of the displacement at the
top of the soil column when a shear force, applied at the top of the soil column, is increased from 0 to 0.3 kN. 
%
It may be noted that, even
though the material model is assumed elastic--perfectly plastic, due to the introduced uncertainties in yielding, the probabilistic
response is non-linear from the beginning. 
%
Depending on the uncertainty in yield strength, there is always a
possibility that the soil becomes elastic--plastic from the very beginning of loading. This
possibility has been quantified, at the Gauss points, from the mean and variance of the yield strength of soil and taken into consideration
implicitly during constitutive simulation using the equivalent advection and diffusion coefficients ($N_{(1)}^{{eq}^{vM}}$ and $N_{(2)}^{{eq}^{vM}}$, refer to Eq.~(\ref{eqno_7})). 
%
These coefficients assign probability weights to the realizations of stress
response based on the probability of material being elastic or elastic-plastic. 
%
As   discussed  in  section  \ref{ProbabilisticYielding},  initially,  at  small
strains,  the  probability  of  material being elastic-plastic is very small and
hence,  the  initial  probabilistic  response  (ensemble of all realizations) is
closer    (but    not    fully)   to  linear,  elastic  response.  However,  as 
deformation
increases,  the  probability  that  the  material   yields  increases and
consequently,  the probabilistic solution gradually becomes more elastic-plastic
(Figure~\ref{figure:mean_and_sd_of_displacements}).
%



Uncertain  tangent shear    modulus  (stiffness)  evolves  as  well at each Gauss point -- 
both in the mean and the probabilistic (KL) spaces -- with the increase in
shear loading. One
such evolution (of mean of tangent shear modulus and of tangent of $\lambda$ at the
first and  the second KL space), at a Gauss point located at a depth of 6.645 m from the top of the
soil  column,  is  shown in Figure~\ref{figure:tangentStiffnessModulus}.

\begin{figure}[!htbp]
\begin{center}
\includegraphics[width=0.48\textwidth]{tangentModulus.ps.eps} 
\hspace*{0.1truecm}
\includegraphics[width=0.48\textwidth]{lambda.ps.eps} \\
(a)\hspace*{6.1truecm}(b)
\caption{Evolution of the statistics of the tangent of the modulus at a depth of
$6.645$m  from  the  top of the soil column, shear force at the top 
increases from 0 to 0.3 kN: (a) mean, (b) $\lambda$ at the 1st and
2nd KL-space}
\label{figure:tangentStiffnessModulus}
\end{center}
\end{figure}

It  is   also interesting to observe the pattern of evolutionary material behavior along  the  
depth  of  the  soil column. 
%
The evolution of mean of tangent
modulus     along    the    depth    of    soil    column    is    plotted    in
Figure~\ref{figure:evolutionTangentModulus}(a).  
The  mean of tangent modulus at
each  Gauss point evolves differently, depending on their respective advection and diffusion
coefficient  (refer to Eq.~(\ref{ADcoeffs01})). As can be
observed   from   Figure~\ref{figure:evolutionTangentModulus}(a),  the  material
around  the  middle of the soil column has evolved more than that around the top and the
bottom  of the soil column. In other words, middle of the soil column has gotten softer
than  the  top and bottom. It may, however, be noted that this stiffness reduction pattern
depends  on the correlation lengths of the random field material properties. For
example,   Figure~\ref{figure:evolutionTangentModulus}(b)  shows  the  stiffness reduction
pattern  of the same soil column, but with a shear modulus correlation length of
1 m.


%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%
% BJ
%\noindent
%\rule[2mm]{\textwidth}{2mm}
%
%SAME SCALE  (plot below)
% 
%\noindent
%\rule[2mm]{\textwidth}{2mm}
% BJ
%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%
\begin{figure}[!htbp]
\begin{center}
\includegraphics[height=0.3\textwidth]{tangentModulusEvolution.ps.eps}
\hspace*{0.1truecm}
\includegraphics[height=0.3\textwidth]{tangentModulusEvolution_MediumCorrelationModulus.ps.eps} \\
(a)\hspace*{5.7truecm}(b)
\caption{Evolution of mean of tangent modulus along the depth of soil column with load: 
(a)  shear  modulus correlation length $0.1$m, and (b) shear modulus correlation
length $1.0$m (all other parameters are kept the same).}
\label{figure:evolutionTangentModulus}
\end{center}
\end{figure}


%OVDE

The  variabilities  in  the above evolved mean behaviors are captured in terms of
$\sqrt{\lambda}f$,  which  is  a measure of standard deviation, at the mutually
perpendicular KL spaces.
%
Figures~\ref{figure:evolutionLambda01}(a)   and   (b)  show  the  evolutions  of
$\sqrt{\lambda}f$  of  the  tangent  modulus  at  the 1st KL-space for the soil
column with correlation length of shear modulus of 0.1m and 1.0 m, respectively.
%
Figures~\ref{figure:evolutionLambda02}(a)  and  (b)  show  the  same for the 2nd
KL-space.
%
As    can    be   observed   from   Figures~\ref{figure:evolutionLambda01}   and
\ref{figure:evolutionLambda02} that at each Gauss point, as in the mean space,
the  material  properties  at  the  KL-spaces also evolved in different fashions, depending upon their
respective     advection     and     diffusion     coefficients     (refer    to
Eq.~(\ref{ADcoeffs02})).
%
The  evolution  of  standard  deviation  of  the  tangent  modulus, which is the
absolute value of the sum of evolutionary $\sqrt{\lambda}f$ over all KL spaces,
is shown in Figure~\ref{figure:evolutionSDTangentModulus}.
%
%
Comparing            Figures~\ref{figure:evolutionTangentModulus}            and
\ref{figure:evolutionSDTangentModulus},  it  might be concluded that the pattern of stiffness reduction 
along the depth of the soil column and its associated uncertainty depend
on the correlation length of the soil modulus.


%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%
% BJ
%\noindent
%\rule[2mm]{\textwidth}{2mm}
%
%SAME SCALE  (plot below)
% 
%\noindent
%\rule[2mm]{\textwidth}{2mm}
% BJ
%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%
\begin{figure}[!htbp]
\begin{center}
\includegraphics[height=0.3\textwidth]{lambda01Evolution.ps.eps}
\hspace*{0.1truecm}
\includegraphics[height=0.3\textwidth]{lambda01Evolution_MediumCorrelationModulus.ps.eps} \\
(a)\hspace*{5.7truecm}(b)
\caption{Evolution  of  tangent  $\sqrt{\lambda_1}~f_1$  along the depth of soil
column   with   load:   
(a)  shear  modulus correlation length $0.1$m, and 
(b) shear modulus correlation length $1.0$m 
(all other parameters are kept the same).}
\label{figure:evolutionLambda01}
\end{center}
\end{figure}


%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%
% BJ
%\noindent
%\rule[2mm]{\textwidth}{2mm}
%
%SAME SCALE  (plot below)
% 
%\noindent
%\rule[2mm]{\textwidth}{2mm}
% BJ
%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%
\begin{figure}[!htbp]
\begin{center}
\includegraphics[width=0.48\textwidth]{lambda02Evolution.ps.eps}
\hspace*{0.1truecm}
\includegraphics[width=0.48\textwidth]{lambda02Evolution_MediumCorrelationModulus.ps.eps} \\
(a)\hspace*{5.7truecm}(b)
\caption{Evolution  of  tangent  $\sqrt{\lambda_2}~f_2$  along the depth of soil
column   with   load:   
(a)  shear  modulus correlation length $0.1$m, and 
(b) shear modulus correlation length $1.0$m 
(all other parameters are kept the same).}
\label{figure:evolutionLambda02}
\end{center}
\end{figure}


%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%
% BJ
%\noindent
%\rule[2mm]{\textwidth}{2mm}
%
%SAME SCALE  (plot below)
% 
%\noindent
%\rule[2mm]{\textwidth}{2mm}
% BJ
%%%%  %%%%%%%%%   %%%%%%%%%%%   %%%%%%%%%%   %%%%%%%%   %%%%%  %%%%%% %%%%%
\begin{figure}[!htbp]
\begin{center}
\includegraphics[height=0.3\textwidth]{SDtangentModulusEvolution.ps.eps}
\hspace*{0.1truecm}
\includegraphics[height=0.3\textwidth]{SDtangentModulusEvolution_MediumCorrelationModulus.ps.eps} \\
(a)\hspace*{5.7truecm}(b)
\caption{Evolution  of  standard deviation of tangent modulus along the depth of
soil  column  with  load:  
(a)  shear  modulus correlation length $0.1$m, and 
(b) shear modulus correlation length $1.0$m 
(all other parameters are kept the same).}
\label{figure:evolutionSDTangentModulus}
\end{center}
\end{figure}

The  influences  of  the  variances  of  shear modulus and shear strength on the
displacement  behavior  at  the  top  of  the  soil  column have also been studied.
Figure~\ref{figure:mean_and_sd_of_displacements_parametric}(a)     shows     the
mean and mean$\pm$standard  deviation  of displacement at the top of the soil column when
COV  of  shear strength is assumed to be 5\%, 
keeping everything else the same as the original model.
%
%{\bf \large perhaps  change this to not refer to figure 3 for data but rather say that
%it  is  the  same  as original model, with changes to COVs... Same for some other
%reference to Figure 3 ...}
%
As expected,
less  non-linear  response  is  obtained  in  the domain of the simulation 
(compare with Fig.~\ref{figure:mean_and_sd_of_displacements}).
%
Low
uncertainty  (COV  =  5\%)  in  shear modulus, on the other hand, yields a more
non-linear                                                              response
(Figure~\ref{figure:mean_and_sd_of_displacements_parametric}(b))  than  the case
where        the        shear       strength       is       less       uncertain
(Figure~\ref{figure:mean_and_sd_of_displacements_parametric}(a)),    but    less
non-linear  response  than  the original model where both the shear strength and the shear modulus
are very uncertain (Figure~\ref{figure:mean_and_sd_of_displacements}).
%
%
\begin{figure}[!htbp]
\begin{center}
\includegraphics[width=0.48\textwidth]{meanPMsd_LowYieldUncertainty.ps.eps}
\hspace*{0.1truecm}
\includegraphics[width=0.48\textwidth]{meanPMsd_LowModulusUncertainty.ps.eps} \\
(a)\hspace*{6.1truecm}(b)
\caption{Evolution  of  mean and mean$\pm$standard deviation of displacement at the top
of  the  12  m  tall stochastic soil column, assuming (a) low uncertainty (COV =
5\%) in shear strength and (b) low uncertainty (COV = 5\%) in shear modulus. 
%
%{\bf \large perhaps  change this to not refer to figure 3 for data but rather say that
%it  is  the  same  as original model, with changes to COVs... Same for some other
%reference to Figure 3 ...}
%
All
other  material  parameters  are assumed to be the same as the original
model, where  both the shear strength and the shear modulus
are very uncertain}
\label{figure:mean_and_sd_of_displacements_parametric}
\end{center}
\end{figure}






%\clearpage

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

\section{CONCLUDING REMARKS}



A      computational formulation    is    developed    for    numerical  solution  of
elastic--plastic      boundary  value  problems  with  stochastic  coefficients. The
framework  is based on spectral formulation of stochastic finite element method.
At   the   Gauss  integration  points,  Fokker-Plank-Kolmogorov  equation  based
probabilistic constitutive integrator is used to update -- as the material plastifies --
the   statistical   properties   of   the  tangent  modulus.  The  advantage  of
Fokker-Plank-Kolmogorov  equation based probabilistic elasto--plasticity is that
it   transforms  the  nonlinear  stochastic  elastic-plastic  constitutive  rate
equation  in  the real  space  into  a  linear  deterministic  partial
differential equation in the probability density space. This deterministic linearity
simplifies  the  numerical  solution  process  of the probabilistic constitutive
equations.   In  addition,  Fokker-Plank-Kolmogorov  approach  to  probabilistic
elasto-plasticity   yields   second-order  accurate  (analytical)  probabilistic
constitutive  solution  and  doesn't  suffer  from  'closure problem' and 'small
variance   requirement',   associated   with  perturbation  technique  and  high
computation  cost, associated with Monte Carlo technique. Further, the presented
formulation   is   general   enough  to  be  used  with  any  elastic-plastic
constitutive  models  with uncertain material parameters. Different constitutive
models  will  only  yield different advection and diffusion coefficients for the
constitutive Fokker-Planck-Kolmogorov equation.



Developed  stochastic framework  is  applied in simulating the behavior of a 1-D  stochastic  soil  shear column
subjected  to a deterministic shear force at the top. The solution
has been  presented  in  terms  of  mean  and standard deviation of displacement. The
evolutionary pattern of stiffness reduction along the depth of the soil column is discussed.
The  influences  of the correlation length and the variance of the soil parameters on
the response have also been addressed.


On a closing note, while, in this paper,  the application  examples  for  the developed framework deal with the uncertain
soil,  the formulation is generic and can be applied to any elastic-plastic or elastic-damage material with uncertain
material properties.






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

\section*{ACKNOWLEDGEMENT}
The   work presented in this paper was supported in part by  a grant from 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).

%\newpage

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

%% The Appendices part is started with the command \appendix;
%% appendix sections are then done as normal sections
%% \appendix

%% \section{}
%% \label{}

%% References
%%
%% Following citation commands can be used in the body text:
%% Usage of \cite is as follows:
%%   \cite{key}          ==>>  [#]
%%   \cite[chap. 2]{key} ==>>  [#, chap. 2]
%%   \citet{key}         ==>>  Author [#]

%% References with bibTeX database:
%\bibliography{SFEM,refmech}
%\bibliography{/home/kallol/publication/Dissertation/SFEM,/home/kallol/publication/Dissertation/refmech}
%\bibliographystyle{acm}
%\bibliographystyle{wileyj}
%\end{thebibliography}
\bibliographystyle{model1-num-names}
\bibliography{SFEM.Dec2010,refmech}
%\bibliography{/home/kallol/publication/Dissertation/SFEM,/home/kallol/publication/Dissertation/refmech}

%% Authors are advised to submit their bibtex database files. They are
%% requested to list a bibtex style file in the manuscript if they do
%% not want to use model1-num-names.bst.

%% References without bibTeX database:

% \begin{thebibliography}{00}

%% \bibitem must have the following form:
%%   \bibitem{key}...
%%

% \bibitem{}

% \end{thebibliography}


\end{document}

%%
%% End of file `elsarticle-template-1-num.tex'.
