\documentclass[sort&compress, review,12pt]{elsarticle}

\usepackage{amsmath,amsfonts,amsthm,amssymb}
\usepackage{graphics,graphicx}
\usepackage{color,xcolor}
\usepackage{rotating}
\usepackage{verbatim}
\usepackage{pdftexcmds}
\usepackage{mathtools}
\usepackage{hyperref}
\usepackage{lineno}
\modulolinenumbers[5]
\usepackage[thinlines]{easytable}
\usepackage{draftwatermark}

\newcommand\pder[2][]{\ensuremath{\frac{\partial#1}{\partial#2}}} 
\newcommand\der[2][]{\ensuremath{\frac{d#1}{d#2}}} 
\SetWatermarkText{Draft rev. 2}
\SetWatermarkAngle{90}
\SetWatermarkHorCenter{1in}

\journal{Computer Methods in Applied Mechanics and Engineering }

%% `Elsevier LaTeX' style
%\bibliographystyle{elsarticle-num}
%\bibliographystyle{elsarticle-harv}
\bibliographystyle{model1-num-names}
\biboptions{square,sort&compress,numbers}

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

\begin{document}

\begin{frontmatter}

\title{Fokker-Planck Linearization for non-Gaussian Stochastic Elastoplastic Finite Elements}

\author[ucd_address]{Konstantinos Karapiperis\corref{mycorrespondingauthor}}
%\ead[url]{www.elsevier.com}

\author[buffalo_address]{Kallol Sett}
\author[ucd_address]{Levent M. Kavvas}
\author[ucd_address,lbnl]{Boris Jeremi\`c}

\cortext[mycorrespondingauthor]{Corresponding author}
%\ead{support@elsevier.com}

\address[ucd_address]{Department of Civil and Environmental Engineering,
 University of California, Davis, CA, USA}
\address[buffalo_address]{Department of Civil, Structural and Environmental Engineering,
 University at Buffalo, The State University at New York, Buffalo, NY, USA}
\address[lbnl]{Earth Science Division, Lawrence Berkeley National Laboratory, Berkeley, CA, USA}

\begin{abstract}
Presented here is   a finite element framework for the solution of stochastic
elastoplastic boundary value problems. The elastic/linearized part is based on a
non-Gaussian  stochastic  Galerkin formulation, where the stiffness random field
is  decomposed  using  Polynomial  Chaos expansion. In the constitutive level, a
Fokker-Planck-Kolmogorov (FPK) plasticity  framework  is  utilized.  A  linearization
procedure  is  developed that serves to update the Polynomial Chaos coefficients
of  the  expanded  random  stiffness  in  the elastoplastic regime, leading to a
nonlinear   least-squares   optimization  problem.  The  proposed  framework  is
illustrated  in a static shear beam example of elastic-perfectly plastic as well
as isotropic hardening material.
\end{abstract}

\begin{keyword}
Fokker-Planck equation, Elastoplasticity, Stochastic Finite Elements, Linearization, Polynomial Chaos
\end{keyword}

\end{frontmatter}

\linenumbers



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


So far, in  constitutive  modeling, material parameters have been defined in a
deterministic   fashion  by  usually  extracting  the  mean  from  a  number  of
experiments.  However,  the  behavior  of  all  engineering materials, let alone
geomaterials,  is  inherently  uncertain,  as  portrayed  by  various  researchers
\citep{Fenton1999b,  Phoon1999,  Jones2002}. The uncertain response follows from
inherent uncertainty of material behavior or spatial non-uniformity of material
distribution.  In  addition,  the nonlinear material behavior present in several
engineering  applications  is usually described using elastoplastic constitutive
relations.  The  physical  or  phenomenological  components  of such a model are
ideally  described  by  random  fields, most of which are non-Gaussian. Modeling
them  as  Gaussian  fields  can  induce  both  inaccuracy and instability to the
solution  of a boundary value problem. For example, a Gaussian representation of
the  material  stiffness  results  in  inaccurate  higher  order  moments  while
physically  allowing  negative realizations of the process to occur (softening).
To  realistically  approximate  such  a  physical  quantity, a strictly positive
definite field is required. 

It  is  generally  accepted  that  intrusive  uncertainty  quantification  (UQ)
frameworks,  in  which  the  uncertainty  is  propagated  through  the governing
differential  equations,  are  more  efficient than non-intrusive ones. However,
most  researchers  have  focused  on  non-intrusive methods, which are easier to
develop  and  utilize  existing  computational  tools,  or  have  limited  their
attention  to  intrusive  UQ  for  simpler  problems.  The simplest example of a
non-intrusive  method  is  Monte  Carlo  Simulation (MCS) \citep{Schueller:1997,
Matthies:1997reviewSFEM},  which  may  be seen as a direct integration method in
which  the  integration  points  are chosen randomly over the probability space.
Depending  on  the application, the latter approach can prove so computationally
demanding   that   any   practical   application   is  hindered,  at  least  for
elasto-plastic models. Lately, more sophisticated sampling-based approaches have
been developed including stochastic collocation \citep{Babuska2007, Xiu2008} and
non-intrusive  Galerkin  techniques  \citep{Giraldi2014}.  The  applicability of
those methods is not affected by the complexity of the problem since they act as
wrappers on a deterministic solver which in turn acts as a "black box".


Several  researchers have dealt so far with intrusive uncertainty quantification
in  computational mechanics with an emphasis in linear problems. A comprehensive
review       of       such       methods       may       be       found       in
\citet{Matthies:1997reviewSFEM,keese:2003ReportSFEM,  Matthies2004},  where  the
authors also provide insight to the well-posedness and structure of a stochastic
boundary  value  problem. So far, the most popular method for the quantification
of   uncertainty   has   been   the  Stochastic  Finite  Element  Method  (SFEM)
\citep{book:Ghanem},  which  relies  on  a  spectral decomposition of parametric
uncertainties  and  a  Polynomial Chaos \citep{Wiener:1938} approximation of the
output  random  field.  It  is  one  of  the  first developments of a stochastic
Galerkin method, where the problem is formulated in a variational form and holds
in   a  weak  sense.  This  class  of  methods  allows  an  explicit  functional
representation  of  the  solution  in  terms of independent random variables. An
overview     of    stochastic    Galerkin    methods    may    be    found    in
\citep{babuska:GalerkinSPDE,    Matthies2004,Matthies2005}.    The    curse   of
dimensionality  associated  with  these  methods  is  one  of  the  topics  that
researchers  have  attempted  to  address lately. \citet{Xiu2002} introduced the
generalized  Polynomial  Chaos  expansion which guarantees optimum (exponential)
convergence  rates  for  different classes of non-Gaussian processes. An optimum
basis from the Askey family of orthogonal polynomials was utilized which reduces
the  dimensionality  of the system. Other researchers have focused on developing
sparse   approximations   by   applying   low-rank   tensor  product  techniques
\citep{Matthies2012},    proper   generalized   decompositions   and   separated
representations \citep{Nouy2010}.


The  first  attempt  to  extend  SFEM  to  nonlinear  material  behavior was by
\citet{Hori:2000},  who  used  a  perturbation  expansion at the stochastic mean
behavior.  In computing the mean behavior they took advantage of bounding media
approximation  by  introducing two fictitious bounding bodies providing an upper
and  a  lower  bound  for  the  mean. This method, however, inherits the "closure
problem"  (essentially the need for higher order statistical moments in order to
calculate   lower  order  statistical  moments)  and  suffers  from  the  "small
coefficient  of  variation"  requirement  for  the  material  parameters. Later,
\citet{Jeremic2005a}  derived  a second-order exact expression for the evolution
of  the  probability  density  function of stress for elastoplastic constitutive
rate    equations    with    uncertain   material   parameters.   Utilizing   an
Eulerian-Lagrangian  form of the Fokker-Planck equation \citep{Kavvas:2003}, the
aforementioned "closure problem" associated with regular perturbation methods is
resolved.  Afterwards,  \citet{Jeremic2009b}  modified their approach to account
for probabilistic rather than expected yielding and incorporated their developed
FPK-based  elastoplastic  model in a Gaussian spectral stochastic finite element
framework  \citep{Sett2010a}. Finally, \citet{Rosic2009} presented in detail the
variational  theory  behind  the  mixed-hardening  stochastic plasticity problem
along  with stochastic versions of relevant established computational plasticity
algorithms.

In this paper, we utilize a Fokker-Planck-Kolmogorov plasticity framework at the
constitutive  level  and  a  stochastic Galerkin framework at the finite element
level.  Non-Gaussian  parametric  uncertainty  is  considered through a combined
Karhunen-Lo{\ `e}ve/Polynomial  Chaos (KL/PC) expansion. The above are coupled through an FPK
linearization  scheme that updates the coefficients of the polynomial chaos (PC)
approximation  of  the  random stiffness. This method may be tailored to provide
varying  order  of  accuracy  counterbalanced  by  computational  efficiency  by
appropriately  choosing  the  KL/PC spaces in which the constitutive integration
procedure   is  performed.  First,  the  stochastic  approximation  schemes  are
discussed  followed  by  the  finite  element  formulation. Next, the underlying
Fokker-Planck-Kolmogorov   framework  is  introduced  along  with  the  proposed
linearization  procedure  and  the  complete  framework  illustrated in a simple
static shear beam example.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Stochastic discretization}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Elastic stiffness}
%
Any  arbitrary non-stationary stiffness random field may be approximated using a
combined  Karhunen-Lo{\ `e}ve/Polynomial  Chaos  methodology. This technique involves
evaluation  of  an  arbitrary  stochastic  process as a polynomial of a suitable
underlying  Gaussian field, whose covariance structure is characterized by means
of  the  Karhunen-Lo{\ `e}ve  expansion  (KLE).  Following  \citet{Sakamoto2002},  we
represent  the  uncertain elastic constitutive tensor field with the help of the
polynomial chaos expansion (PCE):
%
\begin{equation}
D(\mathbf{x}, \theta) = \sum\limits_{i=0}^{M} r_{i}(\mathbf{•}{x}) \Phi_i[\{\xi_r(\theta)\}]
 \label{nonGaussian1}
\end{equation}
%
where  $  \Phi_i[\{\xi_r(\theta)\}]$  is  a  set  of  Hermite  polynomials of an
underlying  Gaussian  set  ${\xi_r(\theta)}$.  It  can  be shown that the PCE is
convergent  in  $L_2(\Omega)$  where  $\Omega$  denotes  the space of the random
variables. A relevant convergence rate study can be found in \citep{Benth1998}.

The spatially dependent coefficients $r_i$ may be computed via simple projection
but this kind of expansion is defined without any reference to the random field
$D(\mathbf{x},\theta)$ and the expected accuracy is low. Therefore,
a correlation structure is endowed to the underlying field by considering 
the following representation:
%
\begin{equation}
D(\mathbf{x}, \theta) = \sum\limits_{i=0}^{M} D_{i}(\mathbf{x}) \Gamma_i[\mathbf{x},\theta)]
 \label{nonGaussian2}
\end{equation}
%
where  $  \Gamma_i[\mathbf{x},\theta]$ is now a set of Hermite polynomials of an
underlying correlated Gaussian field $\gamma(\mathbf{x},\theta)$.
%
The  orthogonality  of the polynomials is employed to calculate the coefficients
$D_i(\mathbf{x})$ as:
%
\begin{equation}
D_i(\mathbf{x}) = \frac{\langle D \Gamma_i\rangle}{\langle \Gamma_i^2 \rangle}
\label{nonGaussian3}
\end{equation}
%
where  the  numerator  can be evaluated using some type of numerical quadrature,
for  example  Monte  Carlo  (MC)  or  Quasi  Monte  Carlo  (QMC) techniques. The
correlation     function     of     $\gamma(\mathbf{x},\theta)$    induced    on
$D(\mathbf{x},\theta)$  is  given  as  the  solution to the following polynomial
equation \citep{Sakamoto2002}:
%
\begin{equation}
\rho_{D}(\mathbf{x}_1,\mathbf{x}_2) = \sum\limits_{i=1}^M D_i(\mathbf{x}_1) D_i(\mathbf{x}_2) i! \left< 
\gamma(\mathbf{x}_1) \gamma(\mathbf{x}_2)\right>^i
 \label{nonGaussian4}
\end{equation}
%
This  equation  is  solved by discretizing the domain into a number of nodes and
solving  the  resulting  system  of equations. Knowing the above, the correlated
random  field  $\gamma(\mathbf{x},\theta)$  may  be  expanded  in  the following
Karhunen-Lo{\ `e}ve form:
%
\begin{equation}
\gamma(\mathbf{x},\theta) = \sum\limits_{i=1}^Q \sqrt{\lambda_i} f_i(\mathbf{x}) \xi_i(\theta)
\label{nonGaussian5}
\end{equation}
%
subject to the following constraint deriving from the unit variance condition imposed on $\gamma(\mathbf{x})$:
%
\begin{equation}
\sum\limits_{i=1}^Q (\sqrt{\lambda_i} f_i(\mathbf{x}))^2 = 1
\end{equation}
%
Thus, it is required that we re-normalize to a unit variance as follows:
%
\begin{equation}
\gamma(\mathbf{x},\theta) = \sum\limits_{i=1}^Q \frac{\sqrt{\lambda_i} f_i(x)}{\sqrt{\sum\limits_{m=1}^Q 
(\sqrt{\lambda_m} f_m(\mathbf{x}))^2}} \xi_i(\theta)
\end{equation}
%
By    equating    the   two   representations   of   $D(\mathbf{x},\theta)$   in
Eq.~(\ref{nonGaussian1}),  (\ref{nonGaussian2})  we  can  find  the coefficients
$r_i(x)$ as
%
\begin{equation}
r_i(\mathbf{x}) = \frac{p!}{\left<\Phi_i^2\right>} D_p(\mathbf{x}) \prod_{j=1}^p 
\frac{\sqrt{\lambda_{k(j)}} f_{k(j)}(\mathbf{x})}{\sqrt{\sum\limits_{m=1}^Q (\sqrt{\lambda_m} 
f_m(\mathbf{x}))^2}}
\end{equation}
%
where $p$ is the order of the polynomial $\Phi_i$ and $k$ is an index on at least one of the 
$\xi_k$ making up $\Phi_i$.
%
Note  that the accuracy of the synthesized marginal probability density function
depends  mainly  on  the  order  $M$  of the PC expansion, while the correlation
accuracy  depends  on  the  dimension  $Q$  of  $\xi(\theta)$. Finally a similar
methodology  may  be  applied  in the case of generalized Polynomial Chaos (gPC)
expansion \citep{Xiu2002}.


This study assumes a strictly positive definite lognormal random stiffness field
in conjunction with classical PCE, which admit the analytical computation of the
respective  coefficients. Assuming an underlying Gaussian field $g(\mathbf{x})$,
the actual stiffness field is given by:
%
\begin{equation}
 D(\mathbf{x}) = e^{g(\mathbf{x})}
 \label{Lognormal1}
\end{equation}
%
with the following mean and variance relations:
%
\begin{eqnarray}
 \bar{D} &=& e^{\bar{g}}\\
 \sigma_D &=& e^{\sigma_g}
 \label{Lognormal2}
\end{eqnarray}
%
The process $g(\mathbf{x})$ is expanded in the Karhunen-Lo{\ `e}ve sense as:
%
\begin{equation}
 g(\mathbf{x}) = \bar{g} + \sum\limits_{i=1}^N g_i \xi_i = \bar{g} + \sum\limits_{i=1}^N \sqrt{\lambda_i} 
f_i(\mathbf{x}) \xi_i
 \label{Lognormal3}
\end{equation}
%
and projection into polynomial chaos yields analytical coefficients $r_i(\mathbf{x})$:
%
\begin{equation}
 r_i(\mathbf{x}) = \frac{\langle e^g \Phi_i \rangle}{\langle \Phi_i^2 \rangle} = \frac{\prod\limits_{j=1}^p 
\sqrt{\lambda_{k(j)}} f_{k(j)}(\mathbf{x})}{\langle \Phi_i^2\rangle} \,e^{\displaystyle{\bar{g} + \frac{1}{2}\sum\limits_{
j=1}^N g_j^2}}
 \label{Lognormal4}
\end{equation}

\vspace{0.1in}
\noindent
Fig.~\ref{Log_Stiffness_PDF}  shows  how  the  synthesized  marginal probability
density  function  using  this  methodology  converges  to  the target lognormal
distribution  for  a  case  of  COV = 30\% for an increasing order of polynomial
chaos  approximation. Finally, Fig.~\ref{Log_Stiffness_Corr} compares the target
and approximated correlation structure for varying KL dimensionality.

\begin{figure}[htbp!]
\centering
\includegraphics[scale=0.7, trim = 0.4in 0in 0in 0in]{figures/Log_Stiffness_PDF.pdf}
\caption{Convergence of the PC approximation (blue) to the target (red) lognormal distribution.}
\label{Log_Stiffness_PDF}
\end{figure}

\begin{figure}[htbp!]
\centering
\includegraphics[scale=0.7, trim = 0.2in 0in 0in 0in]{figures/Log_Stiffness_Corr.pdf}
\caption{Comparing the approximated (blue) to the target (red) correlation structure.}
\label{Log_Stiffness_Corr}
\end{figure}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Shear strength}

In    the    case    of    a    non-Gaussian   shear   strength   random   field
$S_u(\mathbf{x},\theta)$,  the  above methodology may be applied considering the
two  fields  to  be  independent of each other. Alternatively, for computational
efficiency and since the stability of the simulation remains unaffected, one may
choose to approximate the shear strength field simply by KLE as follows:
%
\begin{equation}
S_u(\mathbf{x},\theta) = \bar{S_u}(\mathbf{x}) + \sum\limits_{i=1}^N \sqrt{\lambda_i} f_i(\mathbf{x}) \xi_i(\theta)
\end{equation}
%
by  considering  the  following  Fredholm  integral  equation of the second kind
\citep{Courant1989} with the covariance function $C_{S_u}$ as a kernel:
%
\begin{equation}
\int_D C_{S_u}(\mathbf{x}_1, \mathbf{x}_2) f_k(\mathbf{x}_1)d\mathbf{x}_1 = 
\lambda_k f_k(\mathbf{x}_2)
\label{Fredholm}
\end{equation}
%
This  expansion  is  optimal in the sense that it is the best approximation that
may be achieved in the $L_2(D) \otimes L_2(\Omega)$ norm.


In  some  cases (e.g. triangular, exponential kernel) the above eigenproblem may
be solved analytically, but in the general case a numerical approximation scheme
is  required. In that sense, a number of methods have been applied including FEM
\citep{book:Ghanem},  wavelet-Galerkin \citep{Phoon2002}, $\mathcal{H}$-matrices
\citep{Khoromskij2009} and meshless methods \citep{Rahman2005}.

In a standard finite element setting, each eigenfunction $f_k$ of the kernel is
approximated as:
%
\begin{equation}
f_k(\mathbf{x}) = \sum\limits_{i=1}^N d_{ik} h_i(\mathbf{x})
\end{equation}
%
Utilizing the above representation and requiring the error to be orthogonal to
the approximating space, transforms Eq.~\ref{Fredholm} to the following weak form:
%
\begin{equation}
\hspace{-0.2in}
\sum\limits_{i=1}^N d_{ik}\!\left[ \int_D\!\int_D C_{S_u}(\mathbf{x}_1,\mathbf{x}_2) h_i(\mathbf{x}_2) h_j(\mathbf{x}_1) 
d\mathbf{x}_1 d\mathbf{x}_2 - \lambda_k \int_D h_i(\mathbf{x}) h_j(\mathbf{x}) d\mathbf{x}\right] = 0
 \label{KLGalerkin}
\end{equation}

The  required  discretization  (mesh  size)  depends  on  the correlation length
describing  the  rate  of  fluctuation  of  the  random field. It has been shown
\citep{DerKiureghian:1988, DerKiureghian:1993} that 2-4 elements per correlation
length are usually enough to capture the structure of the random field. In cases
where  the  correlation  structure  is  approximated by long-tailed kernels (e.g
Gaussian),  the  resulting  generalized  "stiffness"  matrix in the eigenproblem
looses  its  sparsity  resulting  in  an  inefficient  numerical solution. It is
therefore common to modify (truncate) the kernels so as to increase the sparsity
of  the  representation.  \citet{Melink2014}  studied  this  problem in terms of
numerical  integration  and  loss  of  positive  definiteness  of the covariance
matrix.
%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Spatial and Stochastic discretization of the solution}

The unknown displacement random field is semi-discretized in the stochastic
dimension using PCE:
%
\begin{equation}
u(\mathbf{x},\theta) = \sum\limits_{i=0}^P d_i(\mathbf{x}) \Psi_i[\xi_r(\theta)]
\end{equation}
%
with the component $d_i(\mathbf{x})$ being further discretized in the spatial
sense using finite element shape functions:
%
\begin{equation}
d_i(\mathbf{x}) = \sum\limits_{j=1}^N d_{ij} N_j(\mathbf{x})
\end{equation}
%
This results in a final expression for the random displacement field:
%
\begin{equation}
u(\mathbf{x},\theta) = \sum\limits_{i=0}^P \sum\limits_{j=1}^N d_{ij} N_j(\mathbf{x}) \Psi_i[\xi_r(\theta)]
\label{DispStochApprox}
\end{equation}

%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Finite Element Formulation}
%
Employing   the   Galerkin   weak   formulation   of   linearized   static   FEM
\citep{local-86}, we have the following simplified form :
%
\begin{equation}
\sum\limits_e \bigg[ \int\limits_{D_e} \nabla\!N_m (\mathbf{x}) D(\mathbf{x},\theta) 
\nabla\!N_n (\mathbf{x}) d\Omega \, u_m - \int_{D_e} f_m (\mathbf{x},\theta) d\Omega
\bigg] = 0
\label{SpatialGalerkin}
\end{equation}
%
where $\sum\limits_e$ denotes the assembly procedure over all finite elements of
the  discretized  domain $\Omega$ and $f_m(\mathbf{x})$ incorporates the various
elemental contributions to the global force vector.

Combining      Equations~\ref{nonGaussian1},      \ref{DispStochApprox}      and
\ref{SpatialGalerkin} and denoting the shape function gradients as:
%
\begin{equation}
\nabla\!N_n (\mathbf{x}) := B_n(\mathbf{x})
\end{equation}
%
yields:
%
\begin{eqnarray}
&& \hspace{-0.6in}\sum\limits_e \bigg[ \int\limits_{D_e} B_m (\mathbf{x}) \sum\limits_{i=0}^{M} r_{i}(\mathbf{x}) \Phi_i[\{\xi_r(\theta)\}] B_n (\mathbf{x})
\sum\limits_{j=0}^P d_{nj} \Psi_j[\xi_r(\theta)]\, d\Omega \nonumber \\ 
&-& \int_{D_e} f_m (\mathbf{x},\theta) d\Omega \bigg] = 0
\label{SpatialGalerkinExpanded}
\end{eqnarray}



Taking now the Galerkin projection of the discretized equation onto each 
arbitrary polynomial basis of the displacement approximation $\Psi_k[\xi_r(\theta)]$:
%
\begin{eqnarray}
&& \hspace{-0.6in} \sum\limits_e \bigg[ \int\limits_{D_e} B_m (\mathbf{x}) \sum\limits_{i=0}^{M} r_{i}(\mathbf{x}) \Phi_i[\{\xi_r(\theta)\}] B_n (\mathbf{x})
\sum\limits_{j=0}^P \sum\limits_{k=0}^P d_{nj} \Psi_j[\xi_r(\theta)] 
 \Psi_k[\xi_r(\theta)]\, d\Omega \nonumber \\ 
&-& \int_{D_e} \sum\limits_{k=0}^P f_m (\mathbf{x},\theta) 
\Psi_k[\xi_r(\theta)]\, d\Omega \bigg] = 0
\label{StochasticGalerkin}
\end{eqnarray}
%
results in the following system of equations:
%
\begin{equation}
\sum_{n=1}^N \sum_{j=0}^P d_{nj} \sum_{k=1}^M b_{ijk} K_{mni} = F_{m}\langle\Psi_k[\{\xi_r\}]\rangle
\end{equation}
%
where
%
\begin{eqnarray}
K_{mni} = \int_D B_m(\mathbf{x}) r_i(\mathbf{x}) B_n(\mathbf{x}) d\Omega
\end{eqnarray}
%
and
%
\begin{eqnarray}
F_m = \int_D f_m (\mathbf{x},\theta) d\Omega
\end{eqnarray}
%
Symbolic     manipulations     are     carried     out     using     Mathematica
\citep{software:Mathematica}  in  order  to  tabulate  the  coefficients  of the
tensor:
%
\begin{eqnarray}
b_{ijk} = \langle \Phi_i[\{\xi_r\}] \Psi_j[\{\xi_r\}] \Psi_k[\{\xi_r\}] \rangle
\end{eqnarray}
%
The  form  of  the  latter  induces  a  special  block sparsity in the resulting
stiffness matrix that may be exploited to develop an efficient solution scheme.
Several  researchers  have  dealt  with such systems of equations arising in the
context  of  the  the spectral stochastic finite element formulation. One of the
first  attempts  was  made  by  \citet{Ghanem1996}  who  proposed  two  solution
procedures,  a  preconditioned  CG method as well as a hierarchical formulation.
Another  iterative scheme of the family of Krylov-subspace methods that has been
applied   is   the   preconditioned  MINRES  \citep{Ullmann2008}.  In  addition,
researchers  have  developed multi-grid approaches \citep{Maitre2010} as well as
incomplete  block-diagonal  preconditioning  schemes based on the FETI-PD solver
\citep{Ghosh2009}.  A  more  complete  review  of  the  methods  may be found in
\cite{Rosic2009}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Elastoplasticity}

In this study, the elastoplastic behavior is treated in a spectral fashion by
updating  the  coefficients  of  the  stochastic  approximation of the stiffness
according   to   an   underlying  Fokker-Planck-Kolmogorov  framework.  At  each
integration  point  and  orthogonal  KL/PC  space, the nonlinear FPK equation is
solved  incrementally,  and  an  optimization  procedure  yields  the equivalent
linearized  advection  and diffusion terms. The updated PC coefficients are then
computed  based on these terms. We investigate varying approximation accuracy by
restricting  the  number of spaces in which the integration procedure is carried
out.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Formulation of FPK-based probabilistic elastoplasticity}

%
The  incremental form of spatial-average elastoplastic constitutive equation may
be written as
%
\begin{eqnarray}
\der[\sigma_{ij}(x_t,t)]{t} = D_{ijkl}(x_t,t) \der[\epsilon_{kl}(x_t,t)]{t}
\end{eqnarray}
%
where $D_{ijkl}$ is the continuum stiffness tensor which can be either elastic or 
elasto-plastic:
%
\begin{eqnarray}
D_{ijkl} = \left\{ \begin{array}{ c c }
    D_{ijkl}^{el}  \qquad &\mbox{;}  f < 0 \vee (f = 0 \wedge df < 0) \\
    \\
    D_{ijkl}^{el} - \displaystyle\Large{\frac{D_{ijmn}^{el} \pder[U]{\sigma_{mn}} 
\pder[f]{\sigma_{pq}} D_{pqkl}^{el}}{\pder[f]{\sigma_{rs}} D_{rstu}^{el} \pder[U]{\sigma_{tu}} - 
\pder[f]{q_*} r_*}}  \qquad & \mbox{;} f = 0 \wedge df < 0
  \end{array} \right.
\label{ParameterTensor}
\end{eqnarray}
%
according to the established Karush-Kuhn-Tucker conditions.


In  the  above equation, $D_{ijkl}^{el}$ is the elastic stiffness tensor, $f$ is
the  yield  function,  which  is a function of stress $\sigma_{ij}$ and internal
variables $q_*$, while $U$ is the plastic potential function.
%
In its most general form, the incremental constitutive equation takes the form
%
\begin{eqnarray}
\der[\sigma_{ij}(x_t,t)]{t} = \beta_{ijkl}(\sigma_{ij}, D_{ijkl}, q_*, r_*; x_t,t) 
\der[\epsilon_{kl}(x_t,t)]{t}
\end{eqnarray}
%
or
%
\begin{eqnarray}
\der[\sigma_{ij}(x_t,t)]{t} = \eta_{ijkl}(\sigma_{ij}, D_{ijkl},\epsilon_{kl}(x_t,t) q_*, r_*; 
x_t,t)
\end{eqnarray}
%
where  the stochasticity of the operator $\beta$ is induced by the stochasticity
of  $D_{ijkl},  q_*,  r_*$.  This renders the above equation a linear/non-linear
ordinary   differential   equation   with   stochastic  coefficients.  Similarly
randomness  in the forcing term ($\epsilon_{kl}$) results in a linear/non-linear
ordinary  differential equation with stochastic forcing. Combining the two cases
yields  a  linear/non-linear  ordinary  differential  equation  with  stochastic
coefficients and stochastic forcing.
%
Using  the  Eulerian-Lagrangian form of the FPE equation \citep{Kavvas:2003} the
above equation takes the following form in the probability density space
%
\begin{eqnarray}
\pder[P(\sigma_{ij},t)]{t} &=& -\pder[]{\sigma_{mn}} \left[ \left\{ \left\langle 
\eta_{mn}(\sigma_{mn}(t), D_{mnrs}, \epsilon_{rs}(t))\right\rangle \right. \right.\nonumber \\ \nonumber
&& + \left. \left. \int_0^t d\tau Cov_0 \left[\pder[\eta_{mn}(\sigma_{mn}(t), D_{mnrs}, 
\epsilon_{rs}(t))]{\sigma_{ab}}; \right. \right. \right.\\ \nonumber
&& \left. \left. \left. \eta_{ab} (\sigma_{ab}(t-\tau), D_{abcd}, 
\epsilon_{cd}(t-\tau))\right]\right\} P(\sigma_{ij}(t),t)\right]\\ \nonumber
&& + \frac{\partial^2}{\partial \sigma_{mn} \partial \sigma_{ab}} \left[\int_0^t d\tau 
Cov_0\left[\eta_{mn} (\sigma_{mn}(t), D_{mnrs}, \epsilon_{rs}(t));\right.\right.\\ 
&& \left.\left. \eta_{ab}(\sigma_{ab}(t-\tau), 
D_{abcd}, \epsilon_{cd}(t-\tau))\right] P(\sigma_{ij}(t),t)\right]
\label{FPK}
\end{eqnarray}
%
where    $P(\sigma_{ij},t)$    is    the    probability   density   of   stress,
$\langle\cdot\rangle$   is  the  expectation  operator,  $Cov_0[\cdot]$  is  the
time-ordered
covariance  operator  and  $\eta_{ij}$  is a generalized random tensor operator.
Details of this derivation can be found in \citep{Jeremic:2007}.
%
The above equation is equivalent to the following generalized form:
%
\begin{eqnarray}
\hspace{-0.1in}
\pder[P(\sigma_{ij},t)]{t} 
&=& 
-\pder[]{\sigma_{mn}} \left[N_{(1)_{mn}}^{\sigma^{eq}} P(\sigma_{ij},t) - \pder[]{\sigma_{ab}} \left\{N_{(2)_{mnab}}^{\sigma^{eq}} P(\sigma_{ij},t)\right\}\right]
\label{FPKGeneral}
\end{eqnarray}
%
where   $N_{(1)}$   and  $N_{(2)}$  are  advection  and  diffusion  coefficients
respectively  that  are  particular to the constitutive model. Given the initial
and  boundary  conditions  as  well  as  the second-order statistics of material
properties,the equation may be solved with second-order accuracy.

\vspace{0.1in}
To    account    for    the   uncertainty   in   the   probabilistic   yielding,
\citet{Jeremic2008c}   introduced   the   following  equivalent  advection  and
diffusion coefficients:
%
\begin{eqnarray}
 N_{(1)_{mn}}^{\sigma^{eq}}(\sigma_{ij}) &=& (1 - P[f > 0])\, N_{(1)_{mn}}^{el} + P[f > 0]\, N_{(1)_{mn}}^{ep} 
 \label{ProbYieldingN1}\\
 N_{(2)_{mnab}}^{\sigma^{eq}}(\sigma_{ij}) &=& (1 - P[f > 0])\, N_{(2)_{mnab}}^{el} + P[f > 0]\, N_{(2)_{mnab}}^{ep}
 \label{ProbYieldingN2}
\end{eqnarray}
%
where $(1 - P[f > 0])$ represents the probability of the material being elastic,
while $P[f > 0]$ represents the probability of the material being elastoplastic.
$P[f  >  0]$  is  obtained from the cumulative density function, rendering it an
explicit  function of the stress $\sigma_{ij}$ as well as the internal variables
$q_*$.


Utilizing   Eq.~\ref{FPK},   one  may  compute  the  elastic  and  elastoplastic
coefficients addressed in Eq.~\ref{FPKGeneral} as:
%
\begin{eqnarray}
 N_{(1)_{mn}}^{el} &=& \langle D_{mnrs}^{el} \dot{\epsilon}_{rs} \rangle \\
 N_{(2)_{mnab}}^{el} &=& t~ \text{Cov}_0 \left[ D_{mnrs}^{el} \dot{\epsilon}_{rs};
 D_{abcd}^{el} \dot{\epsilon}_{cd} \right]
\end{eqnarray}
%
and
%
\begin{eqnarray}
 N_{(1)_{mn}}^{ep} &=& \langle D_{mnrs}^{ep} \dot{\epsilon}_{rs} \rangle +
 \int_0^t d\tau Cov_0\left[ \pder[]{\sigma_{ab}} \{D_{mnrs}^{ep} \dot{\epsilon}_{rs}\};
 D_{abcd}^{ep} \dot{\epsilon}_{cd} \right] \\
 \label{N1ep}
 N_{(2)_{mnab}}^{ep} &=& \int_0^t d\tau Cov_0\left[ D_{mnrs}^{ep}(t) \dot{\epsilon}_{rs};
 D_{abcd}^{ep}(t-\tau) \dot{\epsilon}_{cd} \right]
 \label{N2ep}
\end{eqnarray}


The  evolution  of  an  internal  variable  $\mathbf{q}$ of the model is handled
through a coupled FPK equation of the form:
%
\begin{equation}
\hspace{-0.35in}
 \pder[\tilde{P}(q_i,t)]{t} = - \pder[]{q_m}\bigg[ N_{{(1)}_m}^{q^{eq}}(\sigma_{mn}, q_m) \tilde{P}(q_i,t) - \pder[]{q_n} \bigg\{ N_{{(2)}_{mn}}^{q^{eq}}(\sigma_{mn}, q_m) \tilde{P}(q_i,t) \bigg\} \bigg] 
 \label{IV_FPK}
\end{equation}
%
The advection and diffusion coefficients in the above equation are given similarly
to Eq.~\ref{ProbYieldingN1} and \ref{ProbYieldingN1} but with no contributions of any 'elastic' state: 
%
\begin{eqnarray}
 N_{(1)_{m}}^{q^{eq}}(\sigma_{ij},q_i) &=& P[f > 0]\, N_{(1)_{m}}^{q^{ep}} 
 \label{N1eq_IV}\\
 N_{(2)_{mn}}^{q^{eq}}(\sigma_{ij}, q_i) &=& P[f > 0]\, N_{(2)_{mn}}^{q^{ep}}
 \label{N2eq_IV}
\end{eqnarray}
%
The elastoplastic components of the equivalent advection and diffusion terms are
functions of the so-called loading index or plastic multiplier $L$ and the rates
of evolution of the internal variables ${r_i}$:
%
\begin{eqnarray}
 N_{(1)_{m}}^{q^{ep}} &=& \langle L r_i \rangle + \int\limits_0^t
 d\tau \text{Cov}_0 \bigg[\pder[]{q_j} L r_i(t); L r_j(t-\tau) \bigg]
 \\
 N_{(2)_{mn}}^{q^{ep}} &=& \int\limits_0^t
 d\tau \text{Cov}_0 \bigg[L r_i(t); L r_j(t-\tau) \bigg]
\end{eqnarray}
%
where $L$ may be expressed as:
%
\begin{equation}
 L = \frac{\displaystyle\pder[f]{\sigma_{ot}} D_{ijot}^{el} \dot{\epsilon}_{ij}}
 {\displaystyle\pder[U]{\sigma_{ab}} D_{abcd}^{el}\pder[f]{\sigma_{cd}} - 
 \pder[f]{q_m} r_m}
\end{equation}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Linearization for stiffness update}

The  constitutive  integration of the FPK-based plasticity model cannot directly
provide   the  updated  generalized  stiffness  at  the  finite  element  level.
Therefore, a numerical scheme is required in order to compute the stiffness
in a PC expansion form as per Eq.~\ref{nonGaussian1}. In this study we assume an
equivalent  linear  FPK equation involving the updated PC coefficients which are
deduced through a least-squares optimization procedure.



For  each  orthogonal  KL/PC  space  $s$,  Eq.~\ref{FPKGeneral} applies with the
advection and diffusion coefficients taking the form:
%
\begin{eqnarray}
 N_{(1)_{mn}}^{s^{eq}}(\sigma^s_{ij},\mathbf{x}) &=& (1 - P[f > 0])\, 
 N_{(1)_{mn}}^{el}  + P[f > 0]\, N_{(1)_{mn}}^{ep} 
 \label{N1fem_nl}\\
 N_{(2)_{mnab}}^{s^{eq}}(\sigma^s_{ij},\mathbf{x}) &=& (1 - P[f > 0])\, N_{(2)_{mnab}}^{el} + P[f > 0]\, N_{(2)_{mnab}}^{ep}
 \label{N2fem_nl}
\end{eqnarray}

For  the  purposes  of this study, let us consider an isotropic linear elastic -
Mises  isotropic  hardening  model  and  derive  the  equivalent  advection  and
diffusion  coefficients for this case. The isotropic linear elasticity tensor in
Eq.~\ref{ParameterTensor} reads:
%
\begin{equation}
D_{ijkl}^{el} = G \delta_{ik} \delta_{jl} + \bigg(K - \frac{2}{3}G\bigg) \delta_{ij}
\delta_{kl}
\label{LinearElasticityTensor}
\end{equation}
%
where  $K$  and  $G$  denote  the  bulk  and  shear modulus respectively and are
represented as random fields (see for example Eq.~\ref{nonGaussian1}). 



The elastoplastic continuum tangent tensor in Eq.~\ref{ParameterTensor} is given
in the following general form:
%
\begin{equation}
D_{ijkl}^{ep} = G \delta_{ik} \delta_{jl} + \bigg(K - \frac{2}{3}G\bigg) \delta_{ij}
\delta_{kl} - \frac{A_{ij} A^*_{kl}}{B + K_P}
\label{ContinuumTangentTensor}
\end{equation}
%
The Mises linear hardening yield function is written as:
%
\begin{equation}
f = \sqrt{J_2} - k = \sqrt{\frac{1}{2}s_{ij} s_{ij}} - k
\end{equation}
%
Under the assumption of associated flow rule we have:
%
\begin{equation}
\pder[f]{\sigma_{ij}} = \pder[U]{\sigma_{ij}}
\end{equation}
%
which results in the following symmetry:
%
\begin{equation}
A_{ij} = A^*_{ij} = D_{ijkl}^{el} \pder[f]{\sigma_{kl}}
\end{equation}
%
After some algebraic manipulations, one can easily derive:
%
\begin{equation}
A_{ij} = \frac{G}{\sqrt{J_2}} s_{ij} \qquad B = G
\end{equation}
%
The  plastic  modulus  $K_P$  is  computed  here on the basis of a deterministic
hardening rule in terms of the equivalent plastic strain:
%
\begin{equation}
k = k(e_{eq}^p)
\end{equation}
%
After imposing the consistency condition, we have:
%
\begin{equation}
K_P = -\pder[f]{k} \bar{k} = -\frac{1}{\sqrt{3}} \pder[f]{k} \der[k]{e_{eq}^p} \pder[f]{\sqrt{J_2}} = \frac{1}{\sqrt{3}} \der[k]{e_{eq}^p}
\end{equation}
%
Combining              Equations~\ref{ProbYieldingN1}-\ref{N2ep}             and
\ref{LinearElasticityTensor}-  \ref{ContinuumTangentTensor},  we  can derive the
final coefficients of the FPK constitutive rate equation as:
%
\begin{eqnarray}
\hspace{-0.2in}
 N_{(1)_{mn}}^{s^{eq}} &=& (1 - P[f > 0])\, \bigg\langle \bigg[G \delta_{mr} \delta_{ns} + \bigg(K - \frac{2}{3}G \bigg) \delta_{mn} \delta_{rs} \bigg] \dot{\epsilon}_{rs}(t) \bigg\rangle \nonumber \\
 && + P[f > 0]\, \bigg\langle \bigg[G \delta_{mr} \delta_{ns} + \bigg(K - \frac{2}{3}G \bigg) \delta_{mn} \delta_{rs} - \nonumber\\
 && \frac{1}{G + \frac{1}{\sqrt{3}} \der[k]{e_{eq}^p}} \bigg(\frac{G}{\sqrt{J_2}}\bigg)^2 s^s_{ij}(t) s^s_{kl}(t) \bigg] \dot{\epsilon}_{rs}(t) \bigg\rangle \nonumber \\
 && + \int_0^t d\tau Cov_0\bigg[ \pder[]{\sigma^s_{ab}} \bigg\{ \bigg[G \delta_{mr} \delta_{ns} + \bigg(K - \frac{2}{3}G \bigg) \delta_{mn} \delta_{rs} \nonumber \\ && - \frac{1}{G + \frac{1}{\sqrt{3}} \der[k]{e_{eq}^p}} \bigg(\frac{G}{\sqrt{J_2}}\bigg)^2 s^s_{ij}(t) s^s_{kl}(t) \bigg] \dot{\epsilon}_{rs}(t)\bigg\}; \nonumber \\ 
 && \bigg[G \delta_{ac} \delta_{bd} + \bigg(K - \frac{2}{3}G \bigg) \delta_{ab} \delta_{cd} - \nonumber \\
 && \frac{1}{G + \frac{1}{\sqrt{3}} \der[k]{e_{eq}^p}} \bigg(\frac{G}{\sqrt{J_2}}\bigg)^2 s^s_{ab}(t-\tau) s^s_{cd}(t-\tau) \bigg] 
 \dot{\epsilon}_{cd}(t-\tau) \bigg]
 \label{ProbYieldingN1J2}
 \end{eqnarray}
 
 \begin{eqnarray}
 \hspace{-0.2in}
 N_{(2)_{mnab}}^{s^{eq}} &=& (1 - P[f > 0])\, t~ \text{Cov}_0 \bigg[ 
 \bigg\{G \delta_{mr} \delta_{ns} + \bigg(K - \frac{2}{3}G \bigg) \delta_{mn} \delta_{rs} \nonumber \\ && - \frac{1}{G + \frac{1}{\sqrt{3}} \der[k]{e_{eq}^p}} \bigg(\frac{G}{\sqrt{J_2}}\bigg)^2 s^s_{ij}(t) s^s_{kl}(t) \bigg\} \dot{\epsilon}_{rs}(t);  
 \nonumber \\ 
 && \bigg[G \delta_{ac} \delta_{bd} + \bigg(K - \frac{2}{3}G \bigg) \delta_{ab} \delta_{cd} - \nonumber \\
 && \frac{1}{G + \frac{1}{\sqrt{3}} \der[k]{e_{eq}^p}} \bigg(\frac{G}{\sqrt{J_2}}\bigg)^2 s^s_{ab}(t) s^s_{cd}(t) \bigg] 
 \dot{\epsilon}_{cd}(t) \bigg] \nonumber \\
 && + P[f > 0]\, \int_0^t d\tau Cov_0\bigg[ \bigg\{ G \delta_{mr} \delta_{ns} + \bigg(K - \frac{2}{3}G \bigg) \delta_{mn} \delta_{rs} \nonumber \\ && - \frac{1}{G + \frac{1}{\sqrt{3}} \der[k]{e_{eq}^p}} \bigg(\frac{G}{\sqrt{J_2}}\bigg)^2 s^s_{ij}(t) s^s_{kl}(t) \bigg\} \dot{\epsilon}_{rs}(t); \nonumber \\ 
 && \bigg\{G \delta_{ac} \delta_{bd} + \bigg(K - \frac{2}{3}G \bigg) \delta_{ab} \delta_{cd} - \nonumber \\
 && \frac{1}{G + \frac{1}{\sqrt{3}} \der[k]{e_{eq}^p}} \bigg(\frac{G}{\sqrt{J_2}}\bigg)^2 s^s_{ab}(t-\tau) s^s_{cd}(t-\tau) \bigg\}
 \dot{\epsilon}_{cd}(t-\tau) \bigg]
 \label{ProbYieldingN2J2}
\end{eqnarray}

Next  let  us consider a linearized FPK equation for the stress corresponding to
the  orthogonal  space  $s$  at  the k\textsuperscript{th}-step in the following
form:
%
\begin{eqnarray}
\pder[P^{lin}(\sigma_{ij}^s,t)]{t} &=& -N_{{(1)}_{mn}}^{s^{lin}}\pder[P(\sigma_{ij}^s,t)]{\sigma^s_{mn}} + N_{{(2)}_{mnab}}^{s^{lin}} \pder[^2P(\sigma_{ij}^s,t)]{\sigma^s_{mn}\sigma^s_{ab}}
\label{FPKLinear}
\end{eqnarray}
%
where the linearized advection and diffusion coefficients are given by:
%
\begin{eqnarray}
 N_{{(1)}_{mn}}^{s^{lin}} &=& r_{mnab}^{s^{(k)}} \sum\limits_{i=0}^P \langle \Phi_i \Psi_s \rangle \nonumber \\ 
 && \frac{1}{2\Delta t}\sum\limits_{j=1}^N \bigg[ N_{j,b}(\mathbf{x}) \Delta d_{ija}^{k-1} 
+ N_{j,a}(\mathbf{x}) \Delta d_{ijb}^{k-1} \bigg]
 \label{N1fem_lin}\\
 N_{{(2)}_{mnab}}^{s^{lin}} &=& t\, (r_{mnab}^{s^{(k)}})^2 \sum\limits_{i=0}^P \text{Var} \bigg[\Phi_i \Psi_s \bigg]  \nonumber \\ 
 &&  \frac{1}{4\Delta t^2} \bigg[ N_{j,b}(\mathbf{x}) \Delta d_{ija}^{k-1} + N_{j,a}(\mathbf{x}) \Delta d_{ijb}^{k-1} \bigg]^2
 \label{N2fem_lin}
\end{eqnarray}

In  an  explicit  scheme,  the strain increment at the (k-1)\textsuperscript{th}
step is
utilized,    while    the    fourth-order    tensor    valued   PC   coefficient
$r_{mnab}^{s^{(k)}}$  is  unknown. Depending on the specific constitutive model,
the  above  equations  may  be  simplified to include scalar PC coefficients and
deterministic    bases    in    an    appropriate   tensor   format.   Combining
Equations~\ref{FPKGeneral}   and   \ref{FPKLinear},   one   ends   up   with  an
over-determined residual system of equations in terms of the unknown coefficients
at time step $k$:
%
\begin{equation}
R_i(r_{mnab}^{s^{(k)}}) = \pder[P^{lin}(\boldsymbol{\sigma}_i^s,t)]{t} - \pder[P(\boldsymbol{\sigma}_i^s,t)]{t} = 0, \quad i = 1,\ldots,N
\end{equation}
%
Each  equation corresponds to a single point in the stress domain and the system
of  equations  may  be  solved  in the least squares sense using for example the
Levenberg - Marquardt algorithm \citep{Levenberg:1944}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Varying order of accuracy}

The  outlined  linearization  scheme  is  accurate  to  the  order  of the KL/PC
approximation  of  the  stiffness. However, one can restrict the accuracy of the
method   to  second  order  with  significant  computational  time  savings,  by
considering  the  integration  of  a  single  FPK  equation  at any point in the
discretized   domain.   This   is  achieved  by  considering  the  total  stress
$\boldsymbol{\sigma}$     rather    than    each    KL/PC    stress    component
$\boldsymbol{\sigma}^s$. The linearized equation becomes:

\begin{eqnarray}
\pder[P^{lin}(\sigma_{ij},t)]{t} &=& -N_{(1)}^{lin}\pder[P(\sigma_{ij},t)]{\sigma_{mn}} + N_{(2)}^{lin} \pder[^2P(\sigma_{ij},t)]{\sigma_{mn}\sigma_{ab}}
\label{FPKLinearized}
\end{eqnarray}
%
where:
%
\begin{eqnarray}
\hspace{-0.2in}
 N_{(1)}^{lin} &=& \sum\limits_{s=0}^M r_{mnab}^{s^{(k)}} \sum\limits_{i=0}^P \langle \Phi_i \Psi_s \rangle \nonumber \\
 && \frac{1}{2\Delta t}\sum\limits_{j=1}^N \bigg[ N_{j,b}(\mathbf{x}) \Delta d_{ija}^{k-1} 
+ N_{j,a}(\mathbf{x}) \Delta d_{ijb}^{k-1} \bigg]
 \label{N1fem_lin2}\\
 \hspace{-0.2in}
 N_{(2)}^{lin}(r_s^k, \mathbf{x}) &=& t\, \sum\limits_{s=0}^M (r_{mnab}^{s^{(k)}})^2 \sum\limits_{i=0}^P \text{Var} \bigg[\Phi_i \Psi_s \bigg] \nonumber \\
 && \frac{1}{4\Delta t^2} \bigg[ N_{j,b}(\mathbf{x}) \Delta d_{ija}^{k-1} + N_{j,a}(\mathbf{x}) \Delta d_{ijb}^{k-1} \bigg]^2
 \label{N2fem_lin2}
\end{eqnarray}

Again  the  resulting system of equations is generally over-determined and may be
solved  using  least-squares techniques. Due to the form of the FPK constitutive
integrator,  we  do  not  expect higher order accuracy in the linearized tangent
stiffness.  Indeed,  studying  the  above  system,  it  is  evident  that only 2
coefficients  may  be  determined,  while the polynomial chaos coefficients that
correspond  to  third  and  higher  order  are dependent (negatively correlated)
variables.  It  is  proposed  that  the  higher  order coefficients retain their
elastic  values in order to achieve higher order accuracy during elastic loading
or unloading. 

In a more general sense, one can choose the number of orthogonal KL/PC spaces in
which  the  integration  procedure is carried out based on \textit{a posteriori}
error  estimation  techniques.  However,  a  study  of  the  accuracy  of such a
technique is out of the scope of this study.


%1D:
%\begin{eqnarray}
%\pder[P^{lin}(\sigma_{ij},t)]{t} &=& -N_{(1)}^{lin}\pder[P(\sigma_{ij},t)]{\sigma_{mn}} + N_{(2)}^{lin} \pder[^2P(\sigma_{ij},t)]{\sigma_{mn}\sigma_{ab}}
%\label{FPKLinear}
%\end{eqnarray}
%%
%where:
%%
%\begin{eqnarray}
% N_{(1)}^{lin}(r_s^k, \mathbf{x}) &=& \sum\limits_{s=0}^M r_s^{k} \sum\limits_{i=0}^P \langle \Phi_i \Psi_s \rangle \frac{1}{\Delta t}\sum\limits_{j=1}^N \Delta d_{ij}^{k-1} N_j(\mathbf{x}) 
% \label{N1fem_lin}\\
% N_{(2)}^{lin}(r_s^k, \mathbf{x}) &=& t\, \sum\limits_{s=0}^M (r_s^{k})^2 \sum\limits_{i=0}^P \text{Var} \bigg[\Phi_i \Psi_s \bigg] \frac{1}{\Delta t^2} \bigg(\sum\limits_{j=1}^N \Delta d_{ij}^{k-1} N_j(\mathbf{x})\bigg)^2
% \label{N2fem_lin}
%\end{eqnarray}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Numerical illustrations}

In this last section, the proposed framework is applied to the static loading of
a  shear  beam  representing  a  one-dimensional  soil  column  under  undrained
conditions. The simplified numerical model is shown in Fig.~\ref{1d_shear_beam}.
Two  different  cases  are considered to test the methodology against parameters
that  differentiate  the contribution of each orthogonal space, the evolution of
the stress PDF as well as the global response (Table~\ref{Cases1-2}).


\begin{figure}[!htb]
\centering
\includegraphics[scale=0.4, trim = 0.0in 0in 0in 0in]{figures/1d_shear_beam.pdf}
\caption{Realization of the stiffness random field and simplified numerical model.}
\label{1d_shear_beam}
\end{figure}

\begin{table}[htbp!]
\begin{center}
\begin{TAB}(r,0.5cm,0.5cm)[0.5pt]{c|ccccccccc}{c|cc}
\textit{Case} $\quad$ & $\quad l_c^{G,\sigma_y} \quad$ & $\langle G \rangle \quad$ &
$COV_G \quad$ & $\langle \sigma_y \rangle \quad$ & $COV_{\sigma_y} \quad$ & $k_{hard} \quad$ & $n_{KL} \quad$ & $m_{PC_d} \quad$ & $s_{PC_k} \quad$ \\
1 & 0.2 & 50 & 0.4 & 0.8 & 0.1 & 0 & 2 & 2 & 2\\
2 & 1 & 50 & 0.1 & 0.8 & 0.4 & 10 & 2 & 2 & 2\\
\end{TAB}
\end{center}
\caption{Parameters for the examples in this study}
\label{Cases1-2}
\end{table}


\begin{figure}[!htb]
\centering
\includegraphics[scale=0.7, trim = 0.4in 0in 0in 0in]{figures/PDF_evol_Case1.pdf}
\caption{Evolution of the PDF of shear stress at the first 4 orthogonal KL/PC spaces (Case 1).}
\label{PDF_evol_Case1}
\end{figure}

\begin{figure}[!htb]
\centering
\includegraphics[scale=0.43, trim = 0.4in 0in 0in 0in]{figures/PC_coef_evol_Case1.pdf}
\caption{Evolution of PC coefficients of the linearized random shear stiffness (Case 1).}
\label{PC_coef_evol_Case1}
\end{figure}

\begin{figure}[!htb]
\centering
\includegraphics[scale=0.7, trim = 0.4in 0in 0in 0in]{figures/PC_coef_profile_evol_Case1.pdf}
\caption{Evolution of the profile of PC coefficients of the linearized random shear stiffness along the depth of the shear beam (Case 1).}
\label{PC_coef_profile_evol_Case1}
\end{figure}

\begin{figure}[!htb]
\centering
\includegraphics[scale=0.43, trim = 0.4in 0in 0in 0in]{figures/f_d_Case1.pdf}
\caption{Force- mean $\pm$ std. dev. of displacement plot at the top of the shear column (Case 1).}
\label{f_d_Case1}
\end{figure}

Fig.~\ref{PDF_evol_Case1}  shows the evolution of the PDF of stress at the first
4  orthogonal Kl/PC spaces at the top of the shear beam for \textit{case 1}. Due
to  the small correlation length and large coefficient of variation of the shear
modulus,  all  stress  spaces  are  active. The evolution of the PDF in the mean
shear  stress  space  is  initially  diffusive  and  then  sharpens  in  a quick
transition  to the elastoplastic regime. On the other hand, the remaining stress
spaces    exhibit    mostly    diffusion.    At    the   same   spatial   point,
Fig.\ref{PC_coef_evol_Case1}  shows  the evolution of PC coefficients derived by
means   of  the  proposed  linearization  procedure.  After  a  few  steps,  the
optimization  procedure  has converged and the values of the coefficients remain
almost  constant  until  their sharp decline towards zero at the end of loading.
The evolution of the profile of coefficients (along the depth of the shear beam)
is  given  in  Fig.~\ref{PC_coef_profile_evol_Case1},  where  the  shape  of the
initial  profile  (light color) is determined by the underlying KL eigenvectors.
It is evident that the aforementioned profile values ultimately tend to zero due
to the elastic-perfectly plastic nature of the model. Finally, the global force-
(mean   $\pm$   std.dev)   displacement   response   at  the  top  is  shown  in
Fig.~\ref{f_d_Case1},  where  we  can identify a sharp transition to a perfectly
plastic response.


\begin{figure}[!htb]
\centering
\includegraphics[scale=0.7, trim = 0.4in 0in 0in 0in]{figures/PDF_evol_Case2.pdf}
\caption{Evolution of the PDF of shear stress at the first 4 orthogonal KL/PC spaces (Case 2).}
\label{PDF_evol_Case2}
\end{figure}

\begin{figure}[!htb]
\centering
\includegraphics[scale=0.43, trim = 0.4in 0in 0in 0in]{figures/PC_coef_evol_Case2.pdf}
\caption{Evolution of PC coefficients of the linearized random shear stiffness (Case 2).}
\label{PC_coef_evol_Case2}
\end{figure}

\begin{figure}[!htb]
\centering
\includegraphics[scale=0.7, trim = 0.4in 0in 0in 0in]{figures/PC_coef_profile_evol_Case2.pdf}
\caption{Evolution  of  the  profile of PC coefficients of the linearized random
shear stiffness along the depth of the shear beam (Case 2).}
\label{PC_coef_profile_evol_Case2}
\end{figure}

\begin{figure}[!htb]
\centering
\includegraphics[scale=0.43, trim = 0.4in 0in 0in 0in]{figures/f_d_Case2.pdf}
\caption{Force- mean $\pm$ std. dev. of displacement plot at the top of the shear column (Case 2).}
\label{f_d_Case2}
\end{figure}

\textit{Case  2}  involves  a more uncertain initial yield strength along with a
deterministic  hardening  modulus, which results in the characteristic evolution
of  the  PDF of shear strength at the top as shown in Fig.~\ref{PDF_evol_Case2}.
Due  to  the  large correlation length and small coefficient of variation of the
shear  modulus,  the  mean  stress  space  is mostly active as well as the first
stress  space,  which again is mostly diffusive. The associated values of the PC
coefficients  are  shown  in Fig.~\ref{PC_coef_evol_Case2}, which shows a smooth
decline  of the governing coefficient due to the wide range of the elastoplastic
transition.  Fig.~\ref{PC_coef_profile_evol_Case2}  shows  the  evolution of the
profile  of  the  PC  coefficients of the linearized random stiffness similar to
\textit{Case 1}.
Finally  the  global force-displacement response at the top of the shear beam is
shown  in  Fig.~\ref{f_d_Case2},  where we can identify a smooth transition to a
linear hardening response.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusions}
%
We  have proposed a numerical technique to solve inelastic random boundary value
problems    based   on   stochastic   Galerkin   techniques   and   a   nonlocal
Fokker-Planck-Kolmogorov   plasticity   framework.   It   relies  on  a  general
linearization procedure that couples any functional representation of parametric
uncertainty   with   an  underlying  advection-diffusion  model  describing  its
evolution.  Being  an  intrusive  framework  it  has  the  potential  for higher
convergence  rates  than  conventional non-intrusive techniques, especially when
combined  with  sparse PC representations and efficient FPK solution methods. An
additional  advantage  of  the  method  is its potential to balance accuracy and
computation based on error estimates.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
~\newline
\section{Acknowledgements}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
~\newline
\clearpage
\section*{References}

%\bibliography{refmech}
\bibliography{Konstantinos_refmech.bib}

\end{document}