\documentclass[dvips]{nagauth}
%\documentclass[12pt,letterpaper]{article}
%\usepackage{times}
\usepackage{graphicx}
%\usepackage{amsmath}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% ovo je predefinisanje datuma za cirilicu ( vidi Kopka93 strana 175 )
\def\today
{\number\day . \space \ifcase\month\or
January\or
February\or
March\or
April\or
May\or
June\or
July\or
August\or
September\or
October\or
November\or
December\fi,\space
\number\year}
%time of day as found in layout.tex
% 
%
% TIME OF DAY
%
\newcount\hh
\newcount\mm
\mm=\time
\hh=\time
\divide\hh by 60
\divide\mm by 60
\multiply\mm by 60
\mm=\mm
\advance\mm by \time
\def\hhmm{\number\hh:\ifnum\mm<10{}0\fi\number\mm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%za obicne ps fajlove
%\graphicspath{{/home/jeremic/tex/works/psfigures}
% {/home/jeremic/oofep/fem/ASCE95}
% {/home/jeremic/oofep/mathematica/BModel}
% {/home/jeremic/oofep/fem/LStest}
% {/home/jeremic/oofep/fem/LSTconst}
% {/home/jeremic/oofep/MGM}
% {/home/jeremic/oofep/mathematica/BModel/HardFunc}
% {/home/jeremic/tex/works/Proposals/1999/LDGeomechanics}
% {/home/jeremic/BG/amblemi}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\newcommand{\nDarray}
%{{\bf \lowercase{n}\kern0.05em \lower0.3ex\hbox{D}\lowercase{array}}}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% FEMtools LOGO definition
%\newcommand{\FEMtools}
%{{\bf \uppercase{F}\kern0.25em\uppercase{E}\kern0.25em\uppercase{M}\lowercase{tools}}}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% FEI LOGO definition
%\newcommand{\FEI}
%{{\sf \uppercase{F}\kern0.20em\uppercase{E}\kern0.20em\uppercase{I}}}
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% draft special command for postscript IZBACI ZA SLUZBENU VERZIJU
%
%
% \special{!userdict begin
% /bophook{
% gsave
% 105 120 translate %start position
% 45 rotate % orientation
% /TimesRoman findfont %font
% 150 scalefont % scaling of font
% setfont
% 0 0 moveto
% 0.95 setgray % gray level ( 1 > white ; 0 black )
% (Draft Paper) % text you want to see
% show % or: true charpath for hollow letters
% %true charpath
% stroke grestore}def end}
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% % draft special command for postscript IZBACI ZA SLUZBENU VERZIJU
% %
% %
% \special{!userdict begin
% /bophook{
% gsave
% 200 100 translate %start position
% 60 rotate % orientation
% /TimesRoman findfont %font
% 250 scalefont % scaling of font
% setfont
% 0 0 moveto
% 0.80 setgray % gray level ( 1 > white ; 0 black )
% (Draft) % text you want to see
% show % or: true charpath for hollow letters
% %true charpath
% stroke grestore}def end}
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%JB
%%% from UCthesis.sty
%%\setcounter{topnumber}{4}
%%\def\topfraction{.8}
%%\setcounter{bottomnumber}{1}
%%\def\bottomfraction{.8}
%%\setcounter{totalnumber}{4}
%%\def\textfraction{.2}
%%\def\floatpagefraction{.8}
%%\setcounter{dbltopnumber}{2}
%%\def\dbltopfraction{.7}
%%\def\dblfloatpagefraction{.5}
%%%JB
%%
%%
%%
\def\baselinestretch{2}
\begin{document}
\NAG{1}{6}{01}{28}{01}
\runningheads{Sett, Jeremi{\'c} and Kavvas}
{Nonlinear Hardening/Softening and Probabilistic ElastoPlasticity}
\title{The Role of Nonlinear Hardening/Softening in Probabilistic ElastoPlasticity}
\author{Kallol~Sett, Boris~Jeremi{\'c}\corrauth and M.~Levent~Kavvas}
\address{Department of Civil and Environmental Engineering, University of
California, \\
Davis, CA 95616 \\~\\ {\bf Accepted for publication}}
\corraddr{ Boris Jeremi{\'c}, Department of Civil and Environmental Engineering,
University of California,
One Shields Ave., Davis, CA 95616, \texttt{jeremic@ucdavis.edu}}
%\footnotetext[2]{Please ensure that you use the most up to date class file,
%available from the CFM Home Page at\\
%\texttt{http://www.interscience.wiley.com/jpages/10825010/}}
%\cgsn{NSF}{EEC9701568}
\noreceived{}
\norevised{}
\noaccepted{}
%\begin{center}
%{\sc Version: \today{} , \hhmm{} }
%\end{center}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\newpage
%\tableofcontents
%\newpage
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\newpage
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\begin{abstract}
\section{ABSTRACT}
The elasticplastic modeling and simulations have been studied extensively in
the last century. However, one crucial area of material modeling has received
very little attention. The uncertainties in material properties probably have the
largest influence on many aspects of structural and solids behavior. Despite its
importance, effects of uncertainties of material properties on overall response of
structures and solids have rarely been studied. Most of the small number of
studies on effects of material variabilities have used repetitive deterministic
models through MonteCarlo type simulations. While this approach might appear
sound, it cannot be both computationally efficient and statistically accurate
(have statistically appropriate number of data points).
Recently, we have developed a methodology to solve the
probabilistic elasticplastic differential equations. The methodology is based
on EulerianLagrangian form of the FokkerPlanckKolmogorov equation and
provides for full description of the probability density function (PDF) of
stress response for a given strain.
In this paper we describe our development in some details. In
particular, we investigate the effects of nonlinear hardening on predicted PDF
of stress. As it will be shown, the nonlinear hardening will create a
discrepancy between the most likely stress solution and the deterministic
solution. This discrepancy, in fact, means that the deterministic solution is
not the most likely outcome of the corresponding probabilistic solution if
material parameters are uncertain (and they always are, we just tend to simplify
that fact and use, for example, mean values for deterministic simulations).
A number of examples will be presented, illustrating methodology and main
results, some of which are quite surprising as mentioned above.
\keywords{ElastoPlasticity, Probability Theory, FokkerPlanckKolmogorov Equation}
%\end{abstract}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{INTRODUCTION}
Numerical simulations of mechanical behavior of solids and
structures made of geomaterials are inherently uncertain. The
natural conditions with which geotechnical engineers deal with
are unknown and must be inferred from limited and expensive
observations (Beacher and Christian \cite{book:Christian}).
For example
Fig.~\ref{figure:SoilVariability} shows measured soil properties in Mexico City.
%
%
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.7\textwidth]{SoilVariability}
%\vspace*{0.8cm}
\caption{Measured values of mechanical properties of soil from Mexico city, typical soft spot
(after Baecher and Christian \cite{book:Christian}}
\label{figure:SoilVariability}
\end{center}
\end{figure}
%
%
Large variability in material properties can be noted as
a function of depth.
%This is especially true for unconfined compressive
%strength, preconsolidation load, and friction angle obtained
%from consolidated quick triaxial tests.
It is important to note that this soil formation is considered very
homogeneous (Baecher and Christian \cite{book:Christian}). The usual procedure
in this case is to choose
either constant properties or fit smooth (usually linear) curves to properties
of interest (for example friction angle and/or unconfined compressive strength
and/or preconsolidation pressure (refer Fig.~\ref{figure:SoilVariability})) and
use the assumed (much simplified) properties in modeling and simulations related
to this soil profile. This approach will completely neglect the variations which
might have a large effect on response of this soil profile. The usual remedy is
to apply a large factor of safety, which, it is hoped, will cover up for
neglected information. Use of large factors of safety is becoming unacceptable
as it leads to design solutions that are not only uneconomical but also sometimes not even
safe (Duncan \cite{Duncan:2000}). More appropriate methods for analyzing the random fields
of material properties have been described by Ghanem and Spanos \cite{book:Ghanem},
Matthies and Keese \cite{Matthies2005}, Roberts and Spanos \cite{Roberts1986},
Zhu et al, \cite{Zhu1994}, Soize \cite{Soize1994}.
In recent years, civil engineering practice, and particularly geotechnical
engineering practice has seen an increasing emphasis on quantification of
uncertainties. Quantifications or mathematical descriptions of uncertainties are
usually done within the framework of probability theory, where random soil
parameters are modeled as random variables (if they are specialized to a fixed
location in the domain) or random fields (if they are specialized to the
function of location in the domain). Nice descriptions of different types of
soil parameter uncertainties and their quantification methods were presented by
Phoon \& Kulhawy (\cite{Kulhawy:1999A}, \cite{Kulhawy:1999B}) and Fenton
(\cite{Fenton:1999A}, \cite{Fenton:1999B}). The uncertainty in $G/G_{max}$
curve was recently quantified by Stokoe et al. \cite{Stokoe:2004}.
Eventhough there is a widely acknowledged necessity to model soil properties
using theory of probability, there is very little work on propagation of
uncertainties in soil properties through elasticplastic
constitutive equation. Advanced elastoplasticity based
constitutive models, when properly calibrated, although very
accurate, are, in general, highly sensitive to fluctuations in
model parameters (cf. Borja \cite{borja:OpinionPaper}). In the field of
geomechanics, this sensitivity aspect of constitutive model is
of great importance as the uncertainties associated with soil
properties could outweigh the advantages gained by advanced and
sophisticated models.
First attempt to propagate randomness
through the elasticplastic constitutive equations,
considering random Young's modulus was published only recently
by Anders and Hori \cite{Hori:1999}. They based their work on perturbation
expansion at the stochastic mean behavior.
In computing the mean behavior
they took advantage of bounding media analysis. However, because of
the use of Taylor series expansion, this approach is limited to
problems with small coefficients of variation. According to a recent report by
Sudret and Der Kiureghian \cite{Kiureghian:2000ReportSFEM}, the coefficient
of variation should not exceed 20~\% if this approach is to be used. This
limitation restricts the use of perturbation method to general geomechanics
problems, where the coefficients of variation of soil properties are rarely less
than 20\% (eg. Phoon and Kulhawy \cite{Kulhawy:1999A}).
Another disadvantage of the perturbation method
is that it inherits the socalled "closure problem" (cf. Kavvas
\cite{Kavvas:2003}), where information on higherorder moments is always needed
to calculate lowerorder moments. More recently, Fenton and Griffiths
\cite{Griffiths:2005} used a MonteCarlo type method to propagate uncertainties
through
elasticplastic $c\phi$ soil.
Conventional MonteCarlo simulation, although very accurate,
might be computationally expensive for general elasticplastic
simulations. The method typically requires a statistically appropriate number of
realizations per random variable (often a very large number) in order to satisfy
statistical accuracy. This repetitive use of computationally expensive
deterministic model renders the MonteCarlo simulation impractical.
MonteCarlo simulations, however, find a great use in verification of
analytical developments (cf. Oberkampf et al, \cite{Oberkampf2002}).
In order to overcome the drawbacks associated with MonteCarlo technique and
perturbation method in dealing with nonlinear stochastic differential equations
(SDEs), Kavvas \cite{Kavvas:2003} derived a generic EulerianLagrangian form of
FokkerPlanckKolmogorov equation (FPKE) for the secondorder exact
probabilistic solution of any nonlinear ODE with stochastic coefficient and
stochastic forcing. The development was focusing on applications in water
resources engineering but, it is general enough to apply to any
nonlinear ODE with stochastic coefficient and stochastic forcing.
%
Using the above mentioned EulerianLagrangian form of FPKE, Jeremi{\'c} et
al. \cite{Jeremic:2005a} and Sett et al. \cite{Jeremic:2005b} recently developed
formulation and solution for the general 1D elasticplastic constitutive rate
equation with random material properties and random strain rate.
%
The main advantage of FPKEbased approach is that a complete
probabilistic description (probability density function (PDF)) of stress can be
obtained exactly to secondorder accuracy (to covariance of time), given random material
property(ies) and/or random strain. In addition to that, even if the FPKE is not
solvable in closedform solution, the deterministic linearity of the FPKE in
terms of the state variable, the probability density of stress, considerably simplifies
the numerical solution process.
In this paper we present a general expression, in the form of a partial
differential equation, for evolution of the probability density of stress with
strain for any general 3D pointlocation scale/localaverage form of
elasticplastic constitutive law with material properties and strain rate
modeled as random variables/fields. The expression for the PDF of stress is
developed in a general way, making it usable in any dimension (1,2 or 3D)
and for any type of incremental elasticplastic material model. Application of
the developed methodology is demonstrated through 1D pointlocation scale CamClay
shear constitutive law. In particular, the effects of nonlinear hardening and
softening on the predicted PDFs of stresses are investigated.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{GENERAL EXPRESSION FOR 3D PROBABILISTIC
CONSTITUTIVE RATE EQUATION IN LOCALAVERAGE FORM}
\label{MasterKey}
The incremental form of spatialaverage 3D elasticplastic constitutive
rate equation can be written as:
%
\begin{equation}
\label{eqno_S2}
\frac{d\sigma_{ij}(x_t, t)}{dt}
=
D^{ep}_{ijkl}(\sigma_{ij}, D^{el}_{ijkl}, f, U, q_*, r_*;x_t,t)
\frac{d\epsilon_{kl}(x_t,t)}{dt}
\end{equation}
%
%
where, $D_{ijkl}^{ep}$ is the {\it random, nonlinear elasticplastic
coefficient tensor} which is a function of random stress tensor ($ \sigma_{ij}
$), random elastic moduli tensor ($D^{el}_{ijkl}$), random yield function ($f$),
random potential function ($U$), random internal variables ($q_*$) and random
direction of evolution of internal variables ($r_*$). The random internal
variables ($q_*$) could be scalar (for perfectly plastic and isotropic hardening
models), or secondorder tensor (for translational and rotational kinematic
hardening models), or fourthorder tensor (for distortional hardening models) or
any combinations of the above. The same classification applies to the random
direction of evolution of internal variables ($r_*$). By denoting all the random
material parameters by a parameter tensor as:
%
\begin{equation}
\label{eqno_S2a}
D_{ijkl} = \left[D^{el}_{ijkl}, f, U, q_*, r_*\right]
\end{equation}
%
and introducing a random operator tensor, $\eta_{ij}$, one can write
Eq.~(\ref{eqno_S2}) as,
%
\begin{equation}
\label{eqno_S13}
\frac{d\sigma_{ij}(x,t)}{dt}
=
\eta_{ij}(\sigma_{ij},D_{ijkl}, \epsilon_{kl};x,t)
\end{equation}
%
with an initial condition,
%
\begin{equation}
\label{eqno_S14}
\sigma_{ij}(x,0)=\sigma_{{ij}_0}
\end{equation}
%
%
In Eq.~(\ref{eqno_S13}), the stress tensor $\sigma_{ij}$ can be considered to represent a point
in a $ 9 $dimensional stress ($\sigma$)space and hence Eq.~(\ref{eqno_S13}) determines
the velocity for the point in this space. It is possible to visualize this, by
imagining an initial point in stress space, given by its initial condition
$\sigma_{{ij}_0}$, with a trajectory starting out that describes the
corresponding solution of the
nonlinear stochastic ordinary differential equation (ODE) system
(Eq.~(\ref{eqno_S13})). Let us now consider a cloud of initial points,
described by a density $\rho(\sigma_{ij},0)$ in the $\sigma$space,
and with movements of these points dictated by
Eq.~(\ref{eqno_S13}), the phase
density $\rho$ of $\sigma_{ij}(x,t)$ varies in time according to a continuity
equation which expresses the conservation of all these points in the
$\sigma$space.
Expressing this continuity equation in mathematical terms, one obtains Kubo's
stochastic Liouville equation (Kubo \cite{Kubo:1963}):
%
\begin{equation}
\label{eqno_S15}
\displaystyle \frac{\partial \rho (\sigma_{ij}(x,t),t)}{\partial t}
=
\displaystyle
\frac{\partial}{\partial \sigma_{mn}} \eta_{mn}
\left(\sigma_{mn}(x,t), D_{mnpq}(x),\epsilon_{pq}(x,t)\right)
\rho(\sigma_{ij}(x,t),t)
\end{equation}
%
with initial condition,
%
\begin{equation}
\label{eqno_S16}
\rho(\sigma_{ij},0)=\delta(\sigma_{ij}\sigma_{{ij}_0})
\end{equation}
%
\noindent where $\delta(\cdot)$ is the Dirac delta function. Eq.~(\ref{eqno_S16}) is the
probabilistic restatement in the $\sigma$phase space of the original
deterministic initial condition (Eq.~(\ref{eqno_S14})). One can then apply Van Kampen's
Lemma (Van Kampen \cite{VanKampen:1976}) to obtain:
%
\begin{equation}
\label{eqno_S17}
<\rho(\sigma_{ij},t)>=P(\sigma_{ij},t)
\end{equation}
%
where, the symbol $<\cdot>$ denotes the expectation operation, and
$P(\sigma_{ij},t)$ denotes evolutionary probability density of the state
variable tensor $\sigma_{ij}$ from the constitutive equations.
Therefore, in order to obtain the multivariate probability density function
(PDF), $P(\sigma_{ij},t)$, of the state variable tensor $\sigma_{ij}$, it is
necessary to obtain the deterministic partial differential equation (PDE) of the
$\sigma$space mean phase density $ < \rho (\sigma_{ij},t)> $ from the linear
stochastic PDE system (Eqs.~(\ref{eqno_S15}) and ~(\ref{eqno_S16})). This
necessitates the derivation of the ensemble average form of Eq.~(\ref{eqno_S15})
for $ <\rho (\sigma_{ij},t)> $. The ensemble average form of
Eq.~(\ref{eqno_S15}) was derived by Kavvas \cite{Kavvas:2003} as follows:
%
\begin{eqnarray}
\nonumber
\lefteqn{\frac{\partial \left<\rho (\sigma_{ij}(x_t,t), t)\right>}{\partial t} = \displaystyle \frac{\partial}{\partial \sigma_{mn}}
\left[ \left\{\left< \vphantom{\int_{0}^{t}} \eta_{mn}(\sigma_{mn}(x_t,t), D_{mnrs}(x_t), \epsilon_{rs}(x_t,t))\right> \right. \right.}
\\
\nonumber
&& \left. \left. \int_{0}^{t} d\tau Cov_0\left[\vphantom{\frac{\partial}{\partial}}\eta_{mn}(\sigma_{mn}(x_t,t),D_{mnrs}(x_t),
\epsilon_{rs}(x_t,t)); \right. \right. \right. \\
\nonumber
& & \left. \left. \left. \displaystyle \frac{\partial \eta_{ab} (\sigma_{ab}(x_{t\tau}, t\tau), D_{abcd}(x_{t\tau}),
\epsilon_{cd}(x_{t\tau}, t\tau)}{\partial \sigma_{ab}}\right]\right\}\left<\rho (\sigma_{ij}(x_t,t), t)\right>\right] \\
\nonumber
&+& \displaystyle \frac{\partial}{\partial \sigma_{mn}}\left[ \left\{\int_{0}^{t} d\tau Cov_0\left[\eta_{mn}(\sigma_{mn}(x_t,t),
D_{mnrs}(x_t), \epsilon_{rs}(x_t,t)); \right. \right. \right. \\
\label{eqno_S18}
& & \left. \left. \left. \eta_{ab} (\sigma_{ab}(x_{t\tau}, t\tau), D_{abcd}(x_{t\tau}), \epsilon_{cd}(x_{t\tau}, t\tau)) \right]
\vphantom{\int_{0}^{t}} \right\} \displaystyle \frac{\partial \left<\rho (\sigma_{ij}(x_t,t), t)\right>}{\partial \sigma_{ab}} \right ]
\end{eqnarray}
%
to exact second order (to the order of the covariance time of $ \eta $). In
Eq.~(\ref{eqno_S18}), $ Cov_0[\cdot] $ is the time ordered covariance function
defined by,
%
\begin{equation}
\label{eqno_S19}
Cov_0\left[\eta_{mn}(x,t_1), \eta_{ab}(x,t_2)\right]
=
\left<\eta_{mn}(x,t_1) \eta_{ab}(x,t_2)\right>\left<\eta_{mn}(x,t_1)\right>
\cdot \left<\eta_{ab}(x,t_2)\right>
\end{equation}
%
By combining Eq.~(\ref{eqno_S18}) with Eq.~(\ref{eqno_S17}) and rearranging
the terms yields the following EulerianLagrangian form of
FokkerPlanckKolmogorov equation (FPKE) (cf. Kavvas \cite{Kavvas:2003}):
%
\begin{eqnarray}
\nonumber
\lefteqn{\displaystyle \frac{\partial P(\sigma_{ij}(x_t,t), t)}{\partial t} = \displaystyle \frac{\partial}{\partial \sigma_{mn}}
\left[ \left\{\left< \vphantom{\int_{0}^{t}} \eta_{mn}(\sigma_{mn}(x_t,t), D_{mnrs}(x_t), \epsilon_{rs}(x_t,t))\right> \right. \right.} \\
\nonumber
&+& \left. \left. \int_{0}^{t} d\tau Cov_0 \left[\displaystyle \frac{\partial \eta_{mn}(\sigma_{mn}(x_t,t), D_{mnrs}(x_t),
\epsilon_{rs}(x_t,t))} {\partial \sigma_{ab}}; \right. \right. \right. \\
\nonumber
& & \left. \left. \left. \eta_{ab} (\sigma_{ab}(x_{t\tau}, t\tau), D_{abcd}(x_{t\tau}), \epsilon_{cd}(x_{t\tau}, t\tau)
\vphantom{\int_{0}^{t}} \right] \right \} P(\sigma_{ij}(x_t,t),t) \right] \\
\nonumber
&+& \displaystyle \frac{\partial^2}{\partial \sigma_{mn} \partial \sigma_{ab}} \left[ \left\{ \int_{0}^{t} d\tau Cov_0 \left[
\vphantom{\int_{0}^{t}} \eta_{mn}(\sigma_{mn}(x_t,t), D_{mnrs}(x_t), \epsilon_{rs}(x_t,t)); \right. \right. \right. \\
\label{eqno_S20}
& & \left. \left. \left. \eta_{ab} (\sigma_{ab}(x_{t\tau}, t\tau), D_{abcd}(x_{t\tau}), \epsilon_{cd}(x_{t\tau}, t\tau))
\vphantom{\int_{0}^{t}} \right] \vphantom{\int_{0}^{t}} \right\} P(\sigma_{ij}(x_t,t),t) \right]
\end{eqnarray}
%
to exact second order. The solution of this deterministic linear FPKE
(Eq.~(\ref{eqno_S20})) in terms of $ P(\sigma_{ij},t) $ under appropriate
initial and boundary conditions will yield the PDF of the state variable tensor
$ \sigma_{ij} $. It is important to note that while the original equation
(Eq.~(\ref{eqno_S2})) is nonlinear, the FPKE (Eq.~(\ref{eqno_S20})) is linear
in terms of its unknown, the probability density $ P(\sigma_{ij},t) $ of the
state variable tensor $ \sigma_{ij} $. This linearity, in turn, provides
significant advantages in probabilistic solution of the constitutive rate
equation.
One should also note that Eq.~(\ref{eqno_S20}) is a mixed EulerianLagrangian
equation, since, while the real space location $x_t$ at time $t$ is known, the
location $x_{t\tau}$ is an unknown. If one assumes, for example, small strain theory,
one can relate the unknown location $ x_{t\tau} $ from the known location
$ x_t $ by using the strain rate, $\dot \epsilon $ (=$d \epsilon/dt$) as,
\begin{equation}
\label{eqno_S20a}
d\epsilon=\dot \epsilon \tau = \frac{x_tx_{t\tau}}{x_t}
\end{equation}
%
or, by rearranging
%
\begin{equation}
\label{eqno_S20b}
x_{t\tau}=(1\dot \epsilon \tau)x_t
\end{equation}
%Once the probability density function is obtained it can be used to obtain the
%mean of state variable (H) by usual expectation operation as follows,
%\begin{equation}
%\label{eqno_S20c}
%=\int H(t) P(H(t)) dH(t)
%\end{equation}
%\noindent Another possible way to obtain the mean of state variable is to use
%the equivalence between FPKE and It\^o stochastic differential equation.
%Kavvas \cite{Kavvas:2003} obtained the equaivalent operator vector It\^o
%stochastic differential equation as:
Another interesting aspect to note, possibly of theoretical interest, is the
It\^o stochastic differential equation (SDE) form of the constitutive rate
equation with random material properties and random strain rate
(Eq.~(\ref{eqno_S13})). Utilizing the connection between FPKE and It\^o SDE
(cf. Gardiner \cite{book:Gardiner}), one can write Eq.~(\ref{eqno_S20}) as:
%
\begin{eqnarray}
\nonumber
\lefteqn{d\sigma_{ij} (x,t) =} \\
\nonumber
&& \left\{
\left< \eta_{ij} (\sigma_{ij}(x_t,t), D_{ijkl}(x_t), \epsilon_{kl}(x_t,t)) \vphantom{\int_0^t} \right>
+ \int_0^t d \tau Cov_0 \left[\displaystyle
\frac{\partial \eta_{ij} (\sigma_{ij}(x_t,t), D_{ijkl}(x_t), \epsilon_{kl}(x_t,t))}{\partial \sigma_{ab}}; \right. \right. \\
\label{eqno_S20d}
&& \left. \left. \eta_{ab}(\sigma_{ab}(x_{t\tau},t\tau), D_{abcd}(x_{t\tau}), \epsilon_{cd} (x_{t\tau},t\tau))
\vphantom{\displaystyle \frac{\partial \eta_j (H(x_t,t), A(x_t,t), \Phi(x_t,t))}{\partial \sigma_M}}\right] \right\} dt +
C_{ijm}(\sigma_K,t) dW_m(t)
\end{eqnarray}
%
where,
%
\begin{eqnarray}
\nonumber
\lefteqn{C_{ijm}(\sigma_{ij},t) C_{abm}(\sigma_{ab},t) = 2 \int_0^t d \tau Cov_0 \left[\vphantom{\int_0^t} \eta_{ij} (\sigma_{ij}
(x_t,t), D_{ijkl}(x_t), \epsilon_{kl}(x_t,t)); \right.} \\
\label{eqno_S20e}
& & \left. \eta_{ab}(\sigma_{ab}(x_{t\tau},t\tau), D_{abcd}(x_{t\tau}), \epsilon_{cd} (x_{t\tau},t\tau)) \vphantom{\int_0^t} \right]
\end{eqnarray}
%
and, $ dW_i(t) $ is an increment of the vector Wiener process $W_i$ having the
following properties:
%
\begin{equation}
\label{eqno_S20f1}
=0
\end{equation}
%
and,
%
\begin{eqnarray}
\nonumber
&=& \delta_{ij}dt \ \ \mbox{for $\tau =t$} \\
\label{eqno_S20f}
&=& 0 \ \ \mbox{for $ \tau \not= t$}
\end{eqnarray}
%
It is also interesting to note that all the stochasticity of the original
system of equations (Eq.~(\ref{eqno_S13})) are lumped together in the last term
(Wiener increment term) of the r.h.s of Eq.~(\ref{eqno_S20d}).
In theory, one could use the It\^o form of the constitutive rate equation
(Eq.~(\ref{eqno_S20d})) in obtaining the mean of stress tensor ($\sigma_{ij}$)
as (cf. Kavvas \cite{Kavvas:2003}):
%
\begin{eqnarray}
\nonumber
\lefteqn{\displaystyle \frac{}{dt} =} \\
\nonumber
&& \left< \vphantom{\int_0^t} \eta_{ij} (\sigma_{ij}(x_t,t), D_{ijkl}(x_t), \epsilon_{kl}(x_t,t)) \right>
+ \int_0^t d \tau Cov_0 \left[\displaystyle \frac{\partial \eta_{ij}(\sigma_{ij}(x_t,t), D_{ijkl}(x_t), \epsilon_{kl}(x_t,t))}
{\partial \sigma_{ab}}; \right. \\
\label{eqno_S20g}
& & \left. \eta_{ab}(\sigma_{ab}(x_{t\tau},t\tau), D_{abcd}(x_{t\tau}), \epsilon_{cd}(x_{t\tau},t\tau)) \vphantom{\int_0^t} \right]
\end{eqnarray}
%
but the difficulty is that the stress tensor appearing within $\eta_{ij}(\cdot)$
in the covariance term on the r.h.s of Eq.~(\ref{eqno_S20g}) is random and needs
to be treated accordingly. One possible way is perturbation with respect to
mean, however the ``closure problem'' (cf. Kavvas \cite{Kavvas:2003}) will make this
type of solution problematic. Hence in this study the mean of stress tensor
($\sigma_{ij}$) will be computed by standard operation on the PDF of stress
tensor ($\sigma_{ij}$), which will be obtained by solving the
FokkerPlanckKolmogorov PDE (Eq.~(\ref{eqno_S20})).
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{APPLICATION TO PROBABILISTIC BEHAVIOR OF 1D
CAMCLAY ELASTICPLASTIC CONSTITUTIVE EQUATION}
\label{CamClay}
In this section the applicability of the presented general expression in
obtaining PDF and subsequently mean and variance behaviors of stress will be
demonstrated for 1D, elasticplastic, CamClay constitutive rate
equation.
%Probabilistic solution of stochastic constitutive rate equation
%could be useful in carrying out
%probabilistic integrations at the constitutive level for stochastic finite
%element method. Also it could be used in performing sensitivity analysis of
%constitutive rate equation.
The general incremental form of constitutive equation with single
scalar internal variable can be written as:
%
\begin{equation}
\label{eqno_CC1}
\frac{d\sigma_{ij}(t)}{dt}= D^{ep}_{ijkl}(t) \frac{d\epsilon_{kl}(t)}{dt}
\end{equation}
%
where the stiffness tensor $D_{ijkl}^{ep}$ is given as:
%
\begin{eqnarray}
D_{ijkl}^{ep} = \left\{\begin{array}{ll}
%
D^{el}_{ijkl}
%
%
\;\;\; & \mbox{when~~~} f < 0 \vee (f = 0 \wedge df < 0) \\
%
\\
%
D^{el}_{ijkl}

\displaystyle \displaystyle \frac{D^{el}_{ijmn}
\displaystyle \frac{\partial f}{\partial \sigma_{mn}}
\displaystyle \frac{\partial f}{\partial \sigma_{pq}}
D^{el}_{pqkl}}
{\displaystyle \frac{\partial f}{\partial \sigma_{rs}}
D^{el}_{rstu}
\displaystyle \frac{\partial f}{\partial \sigma_{tu}}

\displaystyle \frac{\partial f}{\partial p_0}\bar p_0}
%
\;\;\; & \mbox{when~~~} f = 0 \vee df = 0
%
\end{array} \right.
\label{el_or_elpl}
\end{eqnarray}
%
where $ D^{el}_{ijkl} $ is the elastic stiffness tensor and
$D^{ep}_{ijkl} $ is the elasticplastic stiffness tensor. The internal variable
$p_0$ in this case depends on plastic volumetric strain and has a physical
meaning (for CamClay) as the maximum hydrostatic stress the soil has
experienced during previous loading history). It was assumed that the material
obeys associative flow rule (as is it is usually assumed for CamClay material
model). Hence, both the yield and the plastic potential functions $f$ are given by
(CamClay material model):
%
\begin{equation}
\label{eqno_CC2}
f=\hat f(p,q) = p^2  p_0p + \frac{q^2}{M^2} = 0
\end{equation}
%
Material parameter $ M $ represents the slope of
the critical state line in the $pq$ space. The stress invariants $p$ and $q$
are given as:
%
\begin{equation}
\label{eqno_CC3}
p
=
 \frac{\sigma_{kk}}{3}
=
\frac{\sigma_{11}+\sigma_{22}+\sigma_{33}}{3}
\end{equation}
%
\begin{equation}
\label{eqno_CC$}
q
=
\sqrt{\frac{3}{2}{s}_{ij}{s}_{ij}}
=
\frac{1}{\sqrt{2}}\left[
(\sigma_{11}\sigma_{22})^2 +
(\sigma_{22}\sigma_{33})^2 +
(\sigma_{33}\sigma_{11})^2 +
6 \sigma_{12}^2 +
6 \sigma_{23}^2 +
6 \sigma_{31}^2 \right]^{1/2}
\end{equation}
%
%
%\begin{samepage}
%\begin{eqnarray}
%p= \frac{1}{3} I_{1}\;\;\;\;
%q= \sqrt{3 J_{2D}} \;\;\;\;
%\cos 3\theta= \frac{3\sqrt{3}}{2}\frac{J_{3D}}{\sqrt({J_{2D})^3}}
%\label{invdef1}
%\end{eqnarray}
%\end{samepage}
%
%
%\begin{samepage}
%\begin{eqnarray}
%I_1 = \sigma_{kk}\;\;\;\;
%J_{2D} = \frac{1}{2}{s}_{ij}{s}_{ij} \;\;\;\;
%J_{3D} = \frac{1}{3} {s}_{ij}{s}_{jk}{s}_{ki}\;\;\;\;
%{s}_{ij} = \sigma_{ij}  \frac{1}{3}\sigma_{kk}\delta_{ij}
%\label{invdef2}
%\end{eqnarray}
%\end{
%
%where, $I_1$ is the first invariant of the stress tensor and $J_2$ is the
%second invariant of the deviatoric stress tensor ${s}_{ij} =
%\sigma_{ij}  {1}/{3}\sigma_{kk}\delta_{ij}$. One
%may also note that in terms of conventional stress invariants $ I_1 $ and $ J_2
%$, Eq.~(\ref{eqno_CC2}) can be written as:
%
%\begin{equation}
%\label{eqno_CC5}
%f=\hat f(I_1, \sqrt{J_2}) = I_1^2  I_1 I_{01} + 27 \frac{J_2}{M^2} = 0
%\end{equation}
%
%\noindent where $ I_{01} $ is analogous to $ p_0 $ of the original cam clay
%formulation and is the value of $ I_1 $ at the intersection of the yield surface
%with the $ I_1 $ axis in the $I_1\sqrt{J_2}$ space.
One can work on components of the Eq.~\ref{el_or_elpl} to obtain:
%
\begin{align}
\nonumber
A_{kl}
&=
\frac{\partial f}{\partial \sigma_{pq}}D_{pqkl}^{el}
= \nonumber \\
& \frac{\partial f}{\partial p}\left[2G\left(\frac{\partial p}{\partial
\sigma_{11}}\delta_{1l}\delta_{1k}
+\frac{\partial p}{\partial\sigma_{22}}\delta_{2l}\delta_{2k}
+\frac{\partial p}{\partial
\sigma_{33}}\delta_{3l}\delta_{3k}\right)+\left(K\frac{2}{3}G
\right)\frac{\partial p}{\partial \sigma_{cd}}\delta_{cd}\delta_{kl}\right]
\nonumber \\
\label{eqno_A6}
&+\frac{\partial f}{\partial q}\left[2G\frac{\partial q}{\partial \sigma_{ij}}\delta_{ik}\delta_{jl}+\left(K\frac{2}{3}G \right)
\frac{\partial q}{\partial \sigma_{ab}}\delta_{ab}\delta_{kl}\right]
%\label{eqno_A3A}
%&= 3K\alpha \delta_{kl} + \displaystyle \frac{G}{\sqrt{J_2}}S_{kl}
\end{align}
%
and
%
\begin{align}
\nonumber
B
&=
\frac{\partial f}{\partial \sigma_{rs}}D_{rstu}^{el}\frac{\partial f}{\partial \sigma_{tu}}
= \nonumber \\
&
\left(\frac{\partial f}{\partial p}\right)^2\left[2G\left\{\left(\frac{\partial p}{\partial
\sigma_{11}}\right)^2+\left(\frac{\partial p}{\partial\sigma_{22}}\right)^2+\left(\frac{\partial p}{\partial
\sigma_{33}}\right)^2\right\} + \left(K\frac{2}{3}G \right)\left\{\frac{\partial p}{\partial \sigma_{ij}}\delta_{ij}\right\}^2\right]
\nonumber \\
\label{eqno_A7}
&+
\left(\frac{\partial f}{\partial q}\right)^2 \left[2G\frac{\partial q}{\partial \sigma_{ij}}\frac{\partial
q}{\partial \sigma_{ij}} +\left(K\frac{2}{3}G \right)\left\{\frac{\partial q}{\partial
\sigma_{ij}}\delta_{ij}\right\}^2 \right]
%\label{eqno_A4A}
%&= 9K\alpha ^2 + G
\end{align}
%
where, $K$ and $G$ are the elastic bulk and shear modulus, respectively.
For the CamClay formulation the material parameter $ p_0 $ represents an internal
variable and its evolution depends on the plastic volumetric strain. Hence, the
plastic modulus $K_P$ becomes:
%
\begin{equation}
\label{eqno_A8}
K_P
=
 \frac{\partial f}{\partial p_0}\bar p_0
=
(1+e_0)\frac{p_0}{\lambda\kappa}\frac{\partial f}{\partial p}
\end{equation}
%
One may note that the differentiations appearing in Eqs.~(\ref{eqno_A6}),
~(\ref{eqno_A7}), and ~(\ref{eqno_A8}), involve stochastic process and hence
cannot be carried out in a deterministic sense. Substituting
Eqs.~(\ref{eqno_A6}), ~(\ref{eqno_A7}), and ~(\ref{eqno_A8}) into
~(\ref{el_or_elpl}) one may obtain the random operator tensor ($\eta_{ij}$)
specialized to CamClay model as:
%
\begin{eqnarray}
\eta_{ij} = \left\{\begin{array}{ll}
%
\left[2G\delta_{ik}\delta_{jl}+\left(K\displaystyle \frac{2}{3}G\right) \delta_{ij}\delta_{kl}\right] \displaystyle
\frac{d\epsilon_{kl}}{dt}
%
%
& \mbox{when~~~} f < 0 \vee (f = 0 \wedge df < 0) \\
%
\\
%
\left[2G\delta_{ik}\delta_{jl}+\left(K\displaystyle \frac{2}{3}G\right) \delta_{ij}\delta_{kl}\displaystyle \frac{A_{ij}A_{kl}}
{B+K_P}\right] \displaystyle \frac{d\epsilon_{kl}}{dt}
%
\;\;\; & \mbox{when~~~} f = 0 \vee df = 0
%
\end{array} \right.
\label{els_or_elpls}
\end{eqnarray}
%
where $A_{ij}$, $B$, and $K_P$ are defined by Eqs.~(\ref{eqno_A6}),
~(\ref{eqno_A7}), and ~(\ref{eqno_A8}) respectively.
Previously developed random operator tensor $\eta_{ij}$ for CamClay material
model will be now specialized to 1D shear behavior, resulting in the
following equation for the PDF of shear stress:
%solving for relating PDF of shear
%stress $\sigma_{12}$ and random shear strain $\epsilon_{12}$:
% one can specialize
%the 3D general localaverage form masterkey relation presented in
%Section~\ref{MasterKey} to cam clay pointlocation scale 1D shear constitutive
%rate equation as:
%
%
\begin{eqnarray}
\nonumber
& &\displaystyle \frac{\partial P(\sigma_{12}(t), t)}{\partial t}=
 \displaystyle \frac{\partial}{\partial \sigma_{12}} \left[ \left\{\vphantom{\int_{0}^{t}}
\left<\eta (\sigma_{12}(t), D_{1212}, \epsilon_{12}(t))\right> \right. \right. \\
\nonumber
&+& \left. \left. \int_{0}^{t} d\tau Cov_0
\left[\displaystyle \frac{\partial \eta (\sigma_{12}(t), D_{1212}, \epsilon_{12}(t))}{\partial \sigma_{12}};
\eta (\sigma_{12}(t\tau), D_{1212}, \epsilon_{12}(t\tau) \right)] \right \} P(\sigma_{12}(t),t) \right] \\
\nonumber
&+& \displaystyle \frac{\partial^2}{\partial \sigma_{12}^2} \left[ \left\{ \int_{0}^{t} d\tau Cov_0 \left[ \vphantom{\int_{0}^{t}}
\eta (\sigma_{12}(t), D_{1212}, \epsilon_{12}(t)); \eta \sigma_{12}(t\tau), D_{1212}, \epsilon_{12}(t\tau) \vphantom{\int_{0}^{t}}
\right] \vphantom{\int_{0}^{t}} \right\} P(\sigma_{12}(t),t) \right] \\
\label{eqno_A9}
\end{eqnarray}
%
%
One may note that while solving Eq.~(\ref{eqno_A9}) for probability
density ($P(\sigma_{12})$) of each fixed value of $\sigma_{12}$ in its domain,
the $\sigma_{12}$ terms appearing within $\eta (\cdot)$ are deterministic and
hence the differentiations involving $\sigma_{12}$ can be carried out in a usual
sense. In Eq.~(\ref{eqno_A9}) the random operator function $\eta$ is different
for preyield (elastic) and postyield (elasticplastic) responses. For
preyield, linear elastic response, one can obtain the random operator function
$\eta$ in the form:
%
\begin{equation}
\label{eqno_A10}
\eta = 2 G \frac{d\epsilon_{12}}{dt}
\end{equation}
%
while, for postyield, nonlinear, elasticplastic response $\eta$ obtains the
form:
\begin{equation}
\label{eqno_A11}
\eta = \left[2G  \displaystyle \frac{\left(36 \displaystyle \frac{G^2}{M^4} \right) \sigma_{12}^2}
{\displaystyle \frac{(1+e_0)p\left(p\displaystyle \frac{3\sigma_{12}^2}{pM^2}\right)^2}{\kappa} + \left(18 \displaystyle \frac{G}{M^4}\right) \sigma_{12}^2
+ \displaystyle \frac{1+e_0}{\lambda\kappa} p p_0 \left(p\displaystyle \frac{3\sigma_{12}^2}{pM^2}\right)} \right] \displaystyle \frac{d\epsilon_{12}}{dt}
\end{equation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{RESULTS AND DISCUSSIONS}
\label{Results}
In this section the FokkerPlanckKolmogorov (FPK) equation is solved
using appropriate initial and boundary conditions. The FPK equation is
specialized to 1D (shear) CamClay constitutive relationship
with different degrees of uncertainties (Eqs.~(\ref{eqno_A9}),
~(\ref{eqno_A10}), and ~(\ref{eqno_A11})) and the solution will be in terms
of the evolution of Probability Density Function (PDF) of shear stress
($\sigma_{12}$) with time (and/or strain), given random material property(ies) and
deterministic shear strain. It should be noted that the strain can be
random as well, and the equations can be solved in that case too, but for
the sake of simplicity, strain (or time) is preset to deterministic value. In
other words, the solution is driven by deterministic strain.
%Five different cases
%were considered  (a) Low overconsolidation ratio (OCR) response with random normally distributed shear modulus ($G$); (b) Low OCR
%response with random normally distributed shear modulus ($G$) and random normally distributed slope of critical state line ($M$);
%(c) Low OCR response with random normally distributed shear modulus ($G$), random normally distributed slope of critical state
%line ($M$), and random normally distributed overconsolidation pressure ($p_0$); (d) Low OCR response with random normally distributed
%shear modulus ($G$) and random normally distributed overconsolidation pressure ($p_0$); (e) High OCR response with random normally
%distributed shear modulus ($G$) and random normally distributed slope of critical state line ($M$).
\subsection{Initial and Boundary Conditions}
The FokkerPlanckKolmogorov PDE specialized to 1D CamClay shear
constitutive rate equation (Eq.~(\ref{eqno_A9}) can be written in the following
simplified form:
%
\begin{eqnarray}
\nonumber
\displaystyle \frac{\partial P(\sigma_{12}, t)}{\partial t}=
&=& \displaystyle \frac{\partial}{\partial \sigma_{12}} \left\{P(\sigma_{12},t) N_{(1)} \right\} + \displaystyle \frac{\partial^2}{\partial
\sigma_{12}^2} \left\{P(\sigma_{12},t)N_{(2)}\right\} \\
\label{eqno_R1a}
&=& \displaystyle \frac{\partial}{\partial \sigma_{12}} \left[P(\sigma_{12},t)N_{(1)}  \displaystyle \frac{\partial}{\partial \sigma_{12}}
\left\{P(\sigma_{12},t)N_{(2)}\right\} \right] \\
\label{eqno_R1b}
&=& \displaystyle \frac{\partial \zeta}{\partial \sigma_{12}}
\end{eqnarray}
%
%
where, $N_{(1)}$ and $N_{(2)}$ are coefficients of the
FokkerPlanckKolmogorov PDE (Eq.~(\ref{eqno_A9})) and represent the
expressions within the curly braces of the first and second term respectively on
the right hand side of Eq.~(\ref{eqno_A9}). They are called the advection and diffusion
coefficients, respectively as the form of Eq.~(\ref{eqno_A9}) closely resembles
advectiondiffusion equation (cf. Gardiner \cite{book:Gardiner}). The symbol $\zeta$ in
Eq.~(\ref{eqno_R1b}) can be considered to be the probability current (cf.
Gardiner \cite{book:Gardiner}), since Eq.~(\ref{eqno_R1b}) is a continuity equation and
the state variable of the equation (Eq.~(\ref{eqno_R1b})) is probability
density.
The initial condition for preyield (elastic) solution for PDF of shear
stress, will be deterministic i.e. at time $t=0$ all the probability mass is
concentrated at $\sigma_{12}=0$ or at some fixed value of $\sigma_{12}$ (if
there was some initial deterministic stress to begin with). Mathematically, this
is represented using Dirac delta function $\delta(\cdot)$,
%
\begin{equation}
\label{eqno_R1}
P(\sigma_{12},0)=\delta(\sigma_{12})
\end{equation}
%
On the other hand, for solving for postyield (elasticplastic) PDF of
shear stress, the initial condition will be represented by the PDF of
shear stress corresponding to the preyield solution (PDF of shear stress at yield).
The initial probability mass, dictated by Eq.~(\ref{eqno_A9}), will advect and
diffuse into the domain of the system throughout the evolution of the
simulation. As for the boundary condition, a reflecting barrier boundary is
chosen since conservation of the probability mass within the system
is required (i.e. no leaking is allowed at the boundaries).
Mathematically this condition is expressed as (cf. Gardiner \cite{book:Gardiner})
%
\begin{equation}
\label{eqno_R2}
\zeta(\sigma_{12},t)_{AtBoundaries}=0
\end{equation}
%
By making the assumption that the material parameters are normally distributed,
theoretically the domain of stress could be from $\infty$ to +$\infty$ and
hence one can write the boundary conditions as:
%
\begin{equation}
\label{eqno_R3}
\zeta(\infty,t)=\zeta(\infty,t)=0
\end{equation}
%
It should also be noted that any other probabilistic distribution for material
parameters can be handled. For example, by assuming that material parameters
obey exponential distribution, theoretically the domain of stress could be from
0 to +$\infty$ and hence one can write the boundary conditions as:
%
\begin{equation}
\label{eqno_R3b}
\zeta(0,t)=\zeta(\infty,t)=0
\end{equation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Numerical Scheme}
\label{NumericalScheme}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The FPKE specialized to 1D CamClay constitutive equation
(Eq.~(\ref{eqno_A9}) with Eqs.~(\ref{eqno_A10}) or ~(\ref{eqno_A11})) with
initial and boundary conditions (Eqs.~(\ref{eqno_R1}) and (\ref{eqno_R3})
respectively) were solved numerically by \textit{method of lines}
(using Mathematica \cite{book:Mathematica}).
The FokkerPlanckKolmogorov PDE is first semidiscretized into a set of simultaneous
ODE systems and then solved by using an ODE solver.
To this end, differentiating by parts the right hand side of
Eq.~(\ref{eqno_R1b}) yields:
%
\begin{eqnarray}
\nonumber
\frac{\partial P}{\partial t} &=& N_{(1)} \frac{\partial P}{\partial \sigma_{12}}  P \frac{\partial N_{(1)}}{\partial \sigma_{12}} + N_{(2)} \frac{\partial^2 P}{\partial \sigma_{12}^2}
+ 2 \frac{\partial P}{\partial \sigma_{12}}\frac{\partial N_{(2)}}{\partial \sigma_{12}} + P \frac{\partial^2 N_{(2)}}{\partial \sigma_{12}^2} \\
\label{eqno_R4}
&=& P\left(\frac{\partial^2 N_{(2)}}{\partial \sigma_{12}}  \frac{\partial N_{(1)}}{\partial \sigma_{12}}\right) + \frac{\partial P}{\partial
\sigma_{12}} \left( 2\frac{\partial N_{(2)}}{\partial \sigma_{12}}N_{(1)}\right) + \frac{\partial^2 P}{\partial \sigma_{12}^2} N_{(2)}
\end{eqnarray}
%
then, using central difference discretization, as shown in
Fig.~\ref{figure:SemiDiscretization}
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.8\textwidth]{StressDiscretization}
\vspace*{0.8cm}
\caption{Spatial Discretization of the FokkerPlanckKolmogorov PDE}
\label{figure:SemiDiscretization}
\end{center}
\end{figure}
%
one can write semidiscretized form of Eq.~(\ref{eqno_R4}) at any
intermediate node $i$ to as:
%
\begin{eqnarray}
\nonumber
\frac{\partial P^{(i)}}{\partial t} = && P^{(i1)} \left[ \frac{N^{(i)}_{(1)}}{2 \Delta \sigma_{12}} + \frac{N^{(i)}_{(2)}}{\Delta \sigma_{12}^2}  \frac{1}
{\Delta \sigma_{12}} \frac{\partial N^{(i)}_{(2)}}{\partial \sigma_{12}} \right]  P^{(i)} \left[ \frac{\partial N^{(i)}_{(1)}}{\partial \sigma_{12}}
+ 2 \frac{N^{(i)}_{(2)}}{\Delta \sigma_{12}^2}  \frac{\partial^2 N^{(i)}_{(2)}}{\partial \sigma_{12}^2} \right] \\
\label{eqno_R5}
&& + P^{(i+1)} \left [ \frac{N^{(i)}_{(1)}}{2
\Delta \sigma_{12}} + \frac{N^{(i)}_{(2)}}{\Delta \sigma_{12}^2} + \frac{1} {\Delta \sigma_{12}} \frac{\partial N^{(i)}_{(2)}}{\partial \sigma_{12}}
\right]
\end{eqnarray}
%
which forms an initial value problem in the time dimension. It is worth
noting that central finite difference technique might not be the most efficient
in solving ndimensional FokkerPlanckKolmogorov PDE. Work is underway to
improve the efficiency of solution method through adaptivity and reduced order
modeling.
%Using forward difference technique, one can introduce the boundary
%condition at the left boundary (node 1 in Fig.~\ref{figure:SemiDiscretization}) as:
%\begin{equation}
%\label{eqno_R6}
%N^{(1)}_{(1)}P^{(1)}  N^{(1)}_{(2)} \displaystyle \frac{P^{(2)}P^{(1)}}{\Delta \sigma_{12}}  P^{(1)}\displaystyle \frac{\partial
%N^{(1)}_{(2)}}{\partial \sigma_{12}}=0
%\end{equation}
%\noindent Or rearranging,
%\begin{equation}
%\label{eqno_R7}
%P^{(1)}= P^{(2)} \left[\displaystyle \frac{\displaystyle \frac{N^{(1)}_{(2)}}{\Delta \sigma_{12}}}{N^{(1)}_{(1)} + \displaystyle \frac{
%N^{(1)}_{(2)}}{\Delta \sigma_{12}}  \displaystyle \frac{\partial N^{(1)}_{(2)}}{\partial \sigma_{12}}} \right]
%\end{equation}
%\noindent Similarly, using backward difference technique, one can introduce the boundary condition at the right boundary (node $n$ in
%Fig.~\ref{figure:SemiDiscretization}) as:
%\begin{equation}
%\label{eqno_R8}
%P^{(n)}= P^{(n1)} \left[\displaystyle \frac{\displaystyle \frac{N^{(n)}_{(2)}}{\Delta \sigma_{12}}}{N^{(n)}_{(1)} + \displaystyle \frac{
%N^{(n)}_{(2)}}{\Delta \sigma_{12}}  \displaystyle \frac{\partial N^{(n)}_{(2)}}{\partial \sigma_{12}}} \right]
%\end{equation}
It is also noted that, eventhough in theory the shear stress ($\sigma_{12}$)
ranges from $\infty$ to +$\infty$, for simulation purpose the domain between
$0.1~{\rm MPa}$ and $+0.1~{\rm MPa}$ is chosen. This choice is based upon
the material properties of the
example problem, and by considering the practical range of shear stress ($\sigma_{12}$).
%
In addition to that, the Dirac delta function used in initial condition of
preyield (elastic) solution (Eq.~(\ref{eqno_R1}))
is numerically approximated using a Gaussian function with zero mean and
standard deviation of $0.00001~{\rm MPa}$, as shown in Fig.~\ref{figure:DiracDelta}.
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.8\textwidth]{ElasticInitCondview2}
\vspace*{0.8cm}
\caption{Approximation of Dirac delta Initial Condition (Eq.~(\ref{eqno_R1}))
with a Gaussian Function of Zero Mean and Standard Deviation
of 0.00001 MPa}
\label{figure:DiracDelta}
\end{center}
\end{figure}
While this approximation introduces some error in solution (PDF of shear stress)
the error quickly disappears as soon as the solution moves away from the initial condition.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Simulation Plots}
We begin by presenting the results of elastic and elasticplastic probability
density functions for shear stress with only one random material parameter.
%
Fig.~\ref{figure:LowOCR_RandomG}(a) shows the evolution of probability density
function (PDF) of shear stress with time/shear strain for a low
overconsolidation ratio (OCR) CamClay material with normally
distributed random shear modulus $G$ (mean value of $2.5~{\rm MPa}$ and a
standard deviation of $0.5~{\rm MPa}$). The other material properties are considered
deterministic and are given as follows: slope of critical state line $M=0.6$;
overconsolidation pressure $p_0=0.2~{\rm MPa}$; applied confinement pressure
$p=0.1~{\rm MPa}$; slopes of the unloadingreloading and normal compression
lines in $e\ln(p^{\prime})$ space $\kappa=0.05$ and $\lambda=0.25$, respectively;
and initial void ratio $e_0=2.18$.
As mentioned earlier, the strain is assumed deterministic (strain driven
algorithm) while the strain rate is assumed constant ($d\epsilon_{12}/dt$ in
Eqs.~(\ref{eqno_A10}) and ~(\ref{eqno_A11})). The numerical value of assumed
constant strain rate can be substituted in equations ~(\ref{eqno_A10}) and
~(\ref{eqno_A11}) for the entire simulation of the evolution of PDF. As
mentioned earlier at the beginning of this section, the FPKE
(Eq.~(\ref{eqno_A9})) describes the evolution of PDF of shear stress
($\sigma_{12}$) with time ($t$), while the shear strain rate
($d\epsilon_{12}/dt$) describes the evolution of shear strain ($\epsilon_{12}$)
with time ($t$) also. Combining the two, the evolution of PDF of shear stress
($\sigma_{12}$) with shear strain ($\epsilon_{12}$) is easily obtained. Time has
been brought in this simulation as an intermediate dimension to help in solution
process, and hence, the numerical value of shear strain rate
($d\epsilon_{12}/dt$) could be any arbitrary value. This value of shear strain
rate will cancel out once one converts the evolution of PDF of shear stress
($\sigma_{12}$) as a function of time, to evolution of PDF of shear stress
($\sigma_{12}$) as a function of shear strain ($\epsilon_{12}$). For simulation
purpose we assumed the arbitrary value of shear strain rate
($d\epsilon_{12}/dt$) as $0.054/{\rm second}$.
Once the evolution of PDF of shear stress ($\sigma_{12}$) is obtained, the
evolution of mean, mode and standard deviations were calculated by standard
operations on the PDF. The results are shown in
Fig.~\ref{figure:LowOCR_RandomG}. In particular,
Fig.~\ref{figure:LowOCR_RandomG}(a) shows the evolution (surface) of PDF of shear stress
($\sigma_{12})$ with time ($t$)/shear strain ($\epsilon_{12}$). The evolution of PDF
of shear stress is given in more detail in Fig.~\ref{figure:LowOCR_RandomG}(b)
showing contours of the PDF of shear stress ($\sigma_{12}$) and also the mean,
mode, and standard deviations. The deterministic
solution, obtained using the mean value of shear modulus ($G$) as
fixed, deterministic input is also shown. As can be seen
from Fig.~\ref{figure:LowOCR_RandomG}(b), the standard deviation of shear stress
increases in elastic regime (before approx. $0.056~{\rm s}$) until the yield
point and then decreases as the solution approaches the critical
state line. In other words, close to critical state the uncertainty in shear
modulus ($G$) didn't have much effect on the uncertainty in shear stress as it
had close to yielding. But on the other hand if one follows the evolution of
deterministic solution one could see that eventhough the deterministic solution
coincides with mean and mode for preyield (elastic) behavior, it deviates for postyield
behavior, suggesting that the deterministic solution is neither the mean nor the
most probable solution for postyield response.
%
%
%\pagebreak
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.7\textwidth]{RandomG_ElasticPlasticPDFView1m}
\
\vspace*{0.1truecm}
\
\begin{center}
(a)
\end{center}
\
\vspace*{0.1truecm}
\
\includegraphics[width=0.55\textwidth]{RandomG_ContourElasticPlasticPDFm}
\
\vspace{0.1truecm}
\
\begin{center}
(b)
\end{center}
\caption{Low OCR Cam Clay Response with Random Normally Distributed Shear
Modulus ($G$): (a) Evolution of PDF and (b) Evolution of Contours of PDF, Mean,
Mode, Standard Deviations, and Deterministic Solution of Shear Stress
($\sigma_{12}$) with Time (t)/Shear Strain ($\epsilon_{12}$)}
\label{figure:LowOCR_RandomG}
\end{center}
\end{figure}
When compared with MonteCarlo simulation, shown in
Fig.~\ref{figure:LowOCR_RandomG_MC}, it can be seen that FPKE approach
predicted the mean behavior exactly but it slightly overpredicted the standard
deviation behavior. This difference is attributed to the approximation used to represent the
initial deterministic condition $\sigma_{12}=~0$ (Dirac delta function). One
may note that at $\epsilon_{12}=~0$, the probability of shear stress
($\sigma_{12}$) at $\sigma_{12}=~0$ should theoretically be 1 i.e. all the
probability mass should theoretically be concentrated at $\sigma_{12}=~0$ and
would be best described by a Dirac delta condition. However, for numerical
simulation of FPKE, that Dirac delta initial condition was approximated with a
Gaussian function of mean zero and standard deviation of 0.00001 $MPa$ as shown in
Fig.~\ref{figure:DiracDelta}. This initial error in standard deviation advected
and diffused into the domain during the simulation of the evolution process.
This error could be minimized by better approximating the Dirac delta initial
condition (but at the expense of computational cost).
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.7\textwidth]{RandomG_ContourElasticPlasticPDF_MCm}
\caption{Comparison of FPK Approach and MonteCarlo Approach for Low OCR Cam
Clay Response with Random Normally Distributed Shear Modulus ($G$) in terms of
Evolution of Mean and Standard Deviation of Shear Stress ($\sigma_{12}$) with
Time ($t$)/Shear Strain ($\epsilon_{12}$)}
\label{figure:LowOCR_RandomG_MC}
\end{center}
\end{figure}
Next few examples examine results (PDF of shear stress) by the introduction of
uncertainties in other material parameters (other than shear modulus $G$).
%
In examples where we use more than one random material property (see next
examples) they are assumed uncorrelated and independent. The presented
development can also deal with correlated random variables, as described by
FokkerPlanckKolmogorov (FPK) equation \ref{eqno_S20}. The advection and
diffusion coefficients $N_{(1)}$ and $N_{(2)}$ (in Eq.~(\ref{eqno_R1a})) will
be changed accordingly for any correlation between random soil properties.
Additional uncertainty is introduced for the slope of critical state line $M$
in terms of a normal distribution with mean value of 0.6 and a standard
deviation of 0.1. The shear modulus $G$ is again treated as random, as before.
The resulting PDF of shear stress $\sigma_{12}$ is now exhibiting a
nonsymmetric distribution in the postyield (elasticplastic) region.
This nonsymmetry is evident from different postyield evolution of mean and
mode of shear stress $\sigma_{12}$ as seen in Fig.~\ref{figure:LowOCR_RandomG_RandomM}(b).
The nonsymmetry is also evident from the trace of the surface of PDF of shear
stress at time ($t) =~0.3{\rm s}$ (or shear strain ($\epsilon_{12})=~1.62\%$) in
Fig.~\ref{figure:LowOCR_RandomG_RandomM}(a).
%
Another interesting aspect to note (refer to
Fig.~\ref{figure:LowOCR_RandomG_RandomM}(b)) is that for this assumed
combination of material properties (random $G$ and random $M$) the postyield
deterministic shear stress ($\sigma_{12}$) underpredicted the postyield most
probable (mode) shear stress ($\sigma_{12}$) but overpredicted the mean shear
stress ($\sigma_{12}$).
%
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.7\textwidth]{RandomG_RandomM_ElasticPlasticPDFView1m}
\
\vspace*{0.1truecm}
\
\begin{center}
(a)
\end{center}
\
\vspace*{0.1truecm}
\
\includegraphics[width=0.55\textwidth]{RandomG_RandomM_ContourElasticPlasticPDFm}
\
\vspace{0.1truecm}
\
\begin{center}
(b)
\end{center}
\caption{Low OCR Cam Clay Response with Random Normally Distributed Shear
Modulus ($G$) and Random Normally Distributed Slope of Critical State Line
($M$): (a) Evolution of PDF and (b) Evolution of Contours of PDF, Mean, Mode,
Standard Deviations, and Deterministic Solution of Shear Stress ($\sigma_{12}$)
with Time (t) /Shear Strain ($\epsilon_{12}$)}
\label{figure:LowOCR_RandomG_RandomM}
\end{center}
\end{figure}
The effect of addition of further uncertainty into the system, in terms of
random normally distributed overconsolidation pressure ($p_0$) with a mean value
of $0.2~{\rm MPa}$ and standard deviation of $0.07~{\rm MPa}$, in addition to
previously introduced random shear modulus ($G$) and random slope of critical
state line ($M$), is shown in Fig.~\ref{figure:LowOCR_RandomG_RandomM_RandomP0}.
%
It can be seen by comparing Figs.~\ref{figure:LowOCR_RandomG_RandomM} and
\ref{figure:LowOCR_RandomG_RandomM_RandomP0}, the additional
uncertainty in overconsolidation pressure ($p_0$) didn't affect much the
probabilistic response of shear stress ($\sigma_{12}$) when compared to its
response to random normally distributed shear modulus ($G$) and random normally
distributed slope of critical state line ($M$)
(Fig.~\ref{figure:LowOCR_RandomG_RandomM}). That is, the responses in
Figs~\ref{figure:LowOCR_RandomG_RandomM} and
\ref{figure:LowOCR_RandomG_RandomM_RandomP0} are quite similar, at least
qualitatively while quantitatively there are some small differences.
%
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.7\textwidth]{RandomG_RandomM_RandomP0_ElasticPlasticPDFView1m}
\
\vspace*{0.1truecm}
\
\begin{center}
(a)
\end{center}
\
\vspace*{0.1truecm}
\
\includegraphics[width=0.55\textwidth]{RandomG_RandomM_RandomP0_ContourElasticPlasticPDFm}
\
\vspace{0.1truecm}
\
\begin{center}
(b)
\end{center}
\caption{Low OCR Cam Clay Response with Random Normally Distributed Shear
Modulus ($G$), Random Normally Distributed Slope of Critical State Line ($M$),
and Random Normally Distributed Overconsolidation pressure ($p_0$): (a)
Evolution of PDF and (b) Evolution of Contours of PDF, Mean, Mode, Standard
Deviations, and Deterministic Solution of Shear Stress ($\sigma_{12}$) with Time
(t)/Shear Strain ($\epsilon_{12}$)}
\label{figure:LowOCR_RandomG_RandomM_RandomP0}
\end{center}
\end{figure}
The uncertainty in overconsolidation pressure ($p_0$) didn't even have much
effect on the probabilistic response (PDF) of shear stress ($\sigma_{12}$) when
considered separately with random normally distributed shear modulus ($G$).
This can be seen by comparing Fig.~\ref{figure:LowOCR_RandomG_RandomP0}, where the
probabilistic behavior of shear stress ($\sigma_{12}$) for low OCR cam clay
model with randomly distributed shear modulus ($G$) and random normally
distributed overconsolidation pressure ($p_0$) was presented, and
Fig.~\ref{figure:LowOCR_RandomG}. Similar to the case where the shear modulus
was the only random parameter (Fig.~\ref{figure:LowOCR_RandomG}), here also the
postyield deterministic shear stress ($\sigma_{12}$) overpredicted the
postyield mean and postyield most probable (mode) shear stress
($\sigma_{12}$).
%
%\pagebreak
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.7\textwidth]{RandomG_RandomP0_ElasticPlasticPDFView2m}
\
\vspace*{0.1truecm}
\
\begin{center}
(a)
\end{center}
\
\vspace*{0.1truecm}
\
\includegraphics[width=0.55\textwidth]{RandomG_RandomP0_ContourElasticPlasticPDFm}
\
\vspace{0.1truecm}
\
\begin{center}
(b)
\end{center}
\caption{Low OCR Cam Clay Response with Random Normally Distributed Shear
Modulus ($G$) and Random Normally Distributed Overconsolidation pressure
($p_0$): (a) Evolution of PDF and (b) Evolution of Contours of PDF, Mean, Mode,
Standard Deviations, and Deterministic Solution of Shear Stress ($\sigma_{12}$)
with Time (t)/Shear Strain ($\epsilon_{12}$)}
\label{figure:LowOCR_RandomG_RandomP0}
\end{center}
\end{figure}
The probabilistic CamClay elasticplastic responses (low OCR) with
different degrees of randomnesses (as presented above) were compared in
Fig.~\ref{figure:LowOCR_PDFComparison}.
In that figure the PDFs of shear stress ($\sigma_{12}$) at shear strain
($\epsilon_{12})~=1.62\%$ are shown and compared with the deterministic value of
shear stress obtained by choosing mean value (which in this case is also mode as
we used normal distribution) of material properties.
%
It is very interesting to note that the
deterministic value of shear stress ($\sigma_{12}$) is different from the most
probable (mode) value of shear stress ($\sigma_{12}$) for all the cases. It
either underpredicted or overpredicted the most probable (mode) value.
%
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.85\textwidth]{LowOCR_ComparisonPDFm}
\caption{Comparison of Shear Stresses at 1.62\% Shear Strain Obtained from Low
OCR Cam Clay Model with different degrees of Randomnesses}
\label{figure:LowOCR_PDFComparison}
\end{center}
\end{figure}
Presented methodology is applicable to both hardening (as shown in examples
above) and softening material response.
%
A high OCR CamClay material sample was analyzed. Following parameters were
used: random normally distributed shear modulus
($G$) with mean of $2.5~{\rm MPa}$ and standard deviation of $0.5~{\rm MPa}$; random
normally distributed slope of critical state line ($M$) with mean of $0.6$ and
standard deviation of $0.1$; deterministic overconsolidation pressure
$p_0=0.8~{\rm MPa}$; deterministic applied confinement pressure
$p=0.02~{\rm MPa}$; deterministic
slopes of the unloadingreloading and normal compression
lines in $e\ln(p^{\prime})$ space $\kappa=0.05$ and $\lambda=0.25$, respectively;
deterministic initial void ratio $e_0=2.18$. The (arbitrary) value of
strain rate is still the same, $d\epsilon_{12}/dt=0.054/{\rm s}$.
The resulting PDF of shear stress is shown in
Fig.~\ref{figure:HighOCR_RandomG_RandomM}(a). Fig.~\ref{figure:HighOCR_RandomG_RandomM}(b)
gives a more detailed view of results,
in terms of evolution of contours of PDF, mean, mode, standard deviation and
deterministic solution of shear stress ($\sigma_{12}$). Similar to the low OCR
simulation with random $G$ and random $M$, the high OCR simulation also yielded
a symmetric distribution of PDF shear stress ($\sigma_{12}$) in the preyield
(elastic) regime and a nonsymmetric distribution of PDF shear stress
($\sigma_{12}$) in the postyield (elasticplastic) regime. Similarly to
probabilistic response for low OCR (with random $G$ and random $M$), the
deterministic shear stress ($\sigma_{12}$) coincided with the mean and most
probable (mode) shear stress in the preyield (elastic)
regime but deviated in the postyield nonlinear elasticplastic regime. The
deterministic shear stress ($\sigma_{12}$) underpredicted the mean shear stress
($\sigma_{12}$) in the entire postyield regime but overpredicted the most
probable (mode) shear stress ($\sigma_{12}$) in the region close to the yield
point, though close to critical state line it underpredicted the most probable
(mode) shear stress ($\sigma_{12}$).
%
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.65\textwidth]{HighOCR_RandomG_RandomM_ElasticPlasticPDFView1m}
\
\vspace*{0.1truecm}
\
\begin{center}
(a)
\end{center}
\
\vspace*{0.1truecm}
\
\includegraphics[width=0.55\textwidth]{HighOCR_RandomG_RandomM_ContourElasticPlasticPDFm}
\
\vspace{0.1truecm}
\
\begin{center}
(b)
\end{center}
\caption{High OCR Cam Clay Response with Random Normally Distributed Shear
Modulus ($G$) and Random Normally Distributed Slope of Critical State Line
($M$): (a) Evolution of PDF and (b) Evolution of Contours of PDF, Mean, Mode,
Standard Deviations, and Deterministic Solution of Shear Stress ($\sigma_{12}$)
with Time (t) /Shear Strain ($\epsilon_{12}$)}
\label{figure:HighOCR_RandomG_RandomM}
\end{center}
\end{figure}
Verification of probabilistic simulations for high OCR sample were again done
using MonteCarlo approach. Fig~\ref{figure:HighOCR_RandomG_RandomM_MC}
compares the evolution of mean and
standard deviation of shear stress ($\sigma_{12}$) obtained using FPKE approach
and MonteCarlo approach. It can be seen that FPKE approach slightly
overpredicted the standard deviation behavior because of the reasons discussed
earlier.
%
\begin{figure}[htbp]
\begin{center}
\includegraphics[width=0.7\textwidth]{HighOCR_RandomG_RandomM_ContourElasticPlasticPDF_MCm}
\caption{Comparison of FPK Approach and MonteCarlo Approach for High OCR Cam Clay Response with Random Normally Distributed
Shear Modulus ($G$) and Random Normally Distributed Slope of Critical State Line ($M$)
in terms of Evolution of Mean and Standard Deviation of Shear Stress ($\sigma_{12}$) with Time ($t$)/Shear Strain
($\epsilon_{12}$)}
\label{figure:HighOCR_RandomG_RandomM_MC}
\end{center}
\end{figure}
%
That is, the FPKE approach is somewhat ``wider'' since the initial
condition was not a Dirac delta function, but rather an approximation. Again,
this deviation in initial condition (from Dirac delta function) can be controlled
by choosing an initial normal distribution for PDF of shear stress with smaller
standard deviation, but this will increase computational cost in solving the
resulting finite difference system of equations.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{CONCLUSIONS}
In this paper we presented methodology to solve the probabilistic
elasticplastic differential equations. The methodology is based on
EulerianLagrangian form of the FokkerPlanckKolmogorov equation and
provides for full description of the probability density function (PDF) of stress
response for given strain.
%
Of particular interest was the modeling of stressstrain constitutive response
for materials with uncertain material parameters. These uncertainties in
material properties are present in most materials, and it is just our
simplification to a single value (say mean) that renders
deterministic properties usually used for simulations
%
The effects of nonlinear hardening/softening on predicted PDF of stress were
investigated in some more detail. It was shown that the nonlinear hardening/softening
creates a discrepancy between the most likely stress solution and the
deterministic solution, meaning that the usually obtained deterministic
solution is not the most
likely outcome of a constitutive stressstrain simulation if material
paramaters are uncertain.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\bibliography{SFEM,refmech}
%\bibliographystyle{acm}
\bibliographystyle{plain}
%\bibliographystyle{alpha}
%\bibliographystyle{agsm}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%\include{Appendix}
\begin{appendix}
\section*{NOTATION}
This is a summary of notations used in mathematical derivations:
\begin{itemize}
\item[] $D_{ijkl}$: Material parameter tensor, function of $D^{el}_{ijkl}$, $f$, $U$,
$q_*$, and $r_*$,
\item[] $D^{el}_{ijkl}$: Elastic stiffness tensor,
\item[] $D^{ep}_{ijkl}$: Elasticplastic tangent stiffness tensor,
\item[] $\sigma_{ij}$: Cauchy stress tensor, in indicial notation,
\item[] $e$: Void ratio of soil,
\item[] $e_0$: Initial void ratio of soil,
\item[] $f$: Yield surface,
\item[] $G$: Shear modulus of soil,
\item[] $K$: Bulk modulus of soil,
\item[] $K_P$: Plastic modulus,
%\item[] $A_{kl}$: Secondorder tensor, equal to
% $\displaystyle \frac{\partial f}{\partial \sigma_{pq}} D^{el}_{pqkl}$,
%
%\item[] $B$: Scalarvalued quantity, equal to
%$\displaystyle\frac{\partial f}{\partial \sigma_{rs}} D^{el}_{rstu}
%\displaystyle\frac{\partial f}{\partial \sigma_{tu}}$,
\item[] $C_{ijm}$: Coefficient tensor in multivariate It\^o equation, in indicial notation,
\item[] $M$: Slope of critical state line in $pq$ space,
\item[] $N_{(1)}$: Advection coefficient of FPK equation,
\item[] $N_{(2)}$: Diffusion coefficient of FPK equation,
\item[] $P(\cdot, t)$: Probability density at time, t,
\item[] $p$: first invariant of stress tensor,
\item[] $p_0$: Hydrostatic pressure, the soil under analysis has experienced before (internal variable in Cam Clay model),
\item[] $\bar p_0$: Direction of evolution of $p_0$,
\item[] $q$: Invariant of stress tensor, equal to $\sqrt{\frac{3}{2}s_{ij} s_{ij}}$
\item[] $q_*$: Internal variables,
\item[] $r_*$: Directions of evolution of internal variables,
\item[] $s_{ij}$: second invariant of the deviatoric stress tensor,
\item[] $t$, $t_1$, $t_2$: time,
\item[] $U$: Potential surface,
\item[] $W_i$: Increment of vector Wiener process,
\item[] $x_t$: Spatial location at time, t,
\item[] $\epsilon_{kl}$: Strain tensor,
\item[] $\dot \epsilon$: Strain rate,
\item[] $\lambda$: Slope of normal consolidation line in $eln(p)$ space,
\item[] $\kappa$: Slope of unloadingreloading line in $eln(p)$ space,
\item[] $\eta_{ij}$: Spatialtemporal operator tensor, function of $\sigma_{ij}$, $D_{ijkl}$, and $\epsilon_{kl}$,
\item[] $\rho(\cdot, t)$: Phase density at time, t,
\item[] $\tau$: time,
\item[] $\zeta$: Probability current
\end{itemize}
\end{appendix}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%$
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%$
\end{document}
\bye