\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}