\documentclass[Journal,letterpaper, SectionNumbers]{ascelike-new}
%% Please choose the appropriate document class option:
% "Journal" produces double-spaced manuscripts for ASCE journals.
% "NewProceedings" produces single-spaced manuscripts for ASCE conference proceedings.
% "Proceedings" produces older-style single-spaced manuscripts for ASCE conference proceedings.
%
%% For more details and options, please see the notes in the ascelike-new.cls file.
% Some useful packages...
\usepackage[utf8]{inputenc}
\usepackage[T1]{fontenc}
\usepackage{lmodern}
\usepackage{graphicx}
\usepackage[figurename=Fig.,labelfont=bf,labelsep=period]{caption}
% \usepackage{subcaption}
\usepackage{amsmath}
%\usepackage{amsfonts}
%\usepackage{amsbsy}
\usepackage{mathtools}
\usepackage{amssymb}
\usepackage{mathrsfs}
\usepackage{newtxtext,newtxmath}
\usepackage[colorlinks=true,citecolor=red,linkcolor=black]{hyperref}
\usepackage{subfig}
\usepackage{multirow}
\usepackage{braket}
\mathchardef\mhyphen="2D
%
% Please add the first author's last name here for the footer:
\NameTag{Hexiang Wang, \today}
% Note that this is not displayed if the NoPageNumbers option is used
% in the documentclass declaration.
%
\begin{document}
% You will need to make the title all-caps
\title{Time Domain Intrusive Probabilistic Seismic Risk Analysis of Nonlinear Shear Frame Structure}
\author[1]{Hexiang Wang}
\author[2]{Fangbo Wang}
\author[1]{Han Yang}
\author[1]{Yuan Feng}
\author[1]{Jeff Bayless}
\author[1]{Norman A. Abrahamson}
\author[1,3]{Boris Jeremi{\'c}*}
\affil[1]{Department of Civil and Environmental Engineering, University of California, Davis, CA, USA}
\affil[2]{School of Civil Engineering, Tianjin University, Tianjin, China}
\affil[3]{Earth and Environmental Sciences Area, Lawrence Berkeley National Laboratory, Berkeley, CA, USA, Email:jeremic@ucdavis.edu}
\maketitle
% Please include an abstract:
\begin{abstract}
Presented is a time domain intrusive framework for probabilistic seismic risk
analysis.
%
Seismic source characterization is mathematically formulated.
%
Methodology for simulating non-stationary seismic motions for given source, path
and site is proposed.
%
Both uncertain motions and uncertain structural parameters are characterized as
random process/field and represented with Hermite polynomial chaos.
%
Intrusive modeling of Armstrong-Fredrick kinematic hardening based on Hermite
polynomial chaos is formulated and incorporated into Galerkin stochastic
elastic-plastic FEM.
%
Time-evolving probabilistic structural response is solved through developed
stochastic elastic-plastic FEM.
%
Following that, formulation for seismic risk analysis is derived.
The framework is illustrated by seismic risk analysis of an eight-story shear
frame structure.
%
Uncertainties are propagated from earthquake source into uncertain structural
system.
%
Difficulties of choosing intensity measure in the conventional framework are
avoided since all the uncertainties and important characteristics (e.g.,
spectrum acceleration $Sa$ and peak ground acceleration $PGA$) of seismic
motions are directly carried by the random process excitations in time domain.
%
Stochastic dynamic equations are solved in an intrusive way, circumventing
non-intrusive Monte Carlo simulations.
%
\end{abstract}
\begin{KeyWords}
Probabilistic seismic risk analysis \enspace Stochastic ground motion \enspace Stochastic elastoplastic FEM \enspace Time domain \enspace Fourier amplitude spectra
\end{KeyWords}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Introduction}
Performance-based Earthquake Engineering (PBEE) \cite{Allin2000} has been a
successful framework that allows for objective and quantitative decision-making
through seismic risk analyses.
%
% These decisions may be related to choosing between alternative designs, retrofitting options, or rehabilitation strategies, and the decision may be arrived at following economic, social or environmental considerations, or some combination of these.
%
% Frameworks for PBEE, such as that of the Pacific Earthquake Engineering Research (PEER) center \cite{Allin2000}, allow this decision-making process to be followed in an objective and quantitative manner.
%
Equation \ref{eq:seismic_risk_analysis} demonstrates state of the art
methodology of seismic risk analysis:
\begin{equation}
\lambda(EDP>z) = \int \underbrace{|\frac{d\lambda(IM>x)}{dx}|}_\text{PSHA} \underbrace{G(EDP>z|IM=x)}_{\text{fragility}} dx
\label{eq:seismic_risk_analysis}
\end{equation}
\noindent
where $\lambda(EDP>z)$ is the annual rate of engineering demand parameter (EDP,
i.e., performance target) exceeding specific level $z$.
%
EDP hazard is computed as convolution of probabilistic seismic hazard analysis
(PSHA) results and structural fragility with respect to intensity
measure (IM) of ground shaking.
%
PSHA, usually done by engineering seismologist, estimates exceedance rate of
intensity measure $\lambda(IM>x)$ considering all possible faults and scenarios
near the engineering site.
%
Structural fragility $G(EDP>z|IM=x)$ defines the exceeding probability of EDP
given ground motion with particular IM level $x$.
%
With properly defined damage measure (DM) as a function of EDP(s), seismic risk
of damage state can be calculated.
%
The choice of IM is crucial in seismic risk analysis, as it serves as proxy of
damaging ground motions and all the uncertainties in ground motion are assumed
could be represented by the variability of IM.
%
Spectral acceleration $Sa(T_0)$ is commonly adopted as IM for building
structures.
%
Many ground motion predictions equations (GMPEs) are developed to quantify the
median and aleatory variability of $Sa(T_0)$ \cite{Gregor2014}.
%
However, the problem is that the scalar spectral acceleration cannot fully
describe the influence of ground-motion variability upon engineering objects.
%
\citeN{Stafford2010} investigated different intensity measures and found that
they are generally not strongly correlated, which indicates that knowledge of
just one IM distribution is not sufficient to describe any of the other
ground-motion characteristics.
%
In addition, $Sa(T_0)$ as IM for surface building structures, is based on frequency domain, linear dynamic analysis of single degree of freedom system.
%
% Spectral acceleration, put forward based on dynamic response of single degree of freedom (SDOF) oscillator in frequency domain, is a very useful damaging proxy for linear elastic response of surface buildings, while it does not correlate well with time domain nonlinear response and are not suitable at all for any other types of engineering objects.
%
% It might not be appropriate for widespread nonlinear problems.
%
When nonlinear inelastic and/or higher mode response is expected, use of $Sa(T_0)$ is not
appropriate.
%
Nonlinear response history analysis (RHA) with spectrum-matched ground motion is
found to give un-conservatively biased estimates \cite{Iervolino2010,Huang2009}.
%
\citeN{Grigoriu2016} showed that generally $Sa(T_0)$ is weakly dependent with
engineering demand parameters for realistic structures and fragilities defined
as functions of $Sa(T_0)$ have large uncertainties and of limited practical use.
%
% The bias phenomenon is attributed to asymmetric relationship between conditional spectral acceleration ordinates at periods affecting inelastic response \cite{Seifried2016}.
%
Furthermore, for many other engineering objects (e.g., dams, deeply embedded
structures, etc.), it is very difficult to find a proper IM in engineering
practices.
%
For example, choice of IM among peak ground acceleration (PGA), peak ground
velocity (PGV), Arias intensity (AI) and cumulative absolute velocity (CAV) has
been contentiously argued for deformation analysis of dam embankment
\cite{Davoodi2013}.
%
Though Vector-valued PSHA \cite{Baker2007} was put forward to mitigate this
issue, it is rarely performed in practices.
%
The difficulty lies in fragility computation.
%
The fragility becomes a function of vector IMs (e.g., a fragility surface for
two IMs), which requires a large number of structural analyses to be quantified.
%
Properly choosing multiple IMs is also a problem.
%
Many times, even if proper IMs, such as AI and CAV, are identified, additional
efforts are still needed to develop GMPE for these IMs and their correlation.
%
An effective solution to the aforementioned problems would be to remove
intensity measure (IM) as an intermediate proxy from risk calculation.
%
% To this end, we established a time domain stochastic framework for seismic risk analysis.
%
With this in mind, a time domain intrusive framework for probabilistic seismic
risk analysis is developed and described here.
%
The framework is based on the progress of Fourier amplitude spectrum (FAS)
modeling of seismic motions over last several decades \cite{Brune1970,Boore1983,Boore2003,Boore2015}.
%
Recent advances in inter-frequency correlation of FAS \cite{Stafford2017b,Bayless2017} and Fourier phase derivative modeling \cite{Baglio2017} are also
taken into account.
%
Uncertain motions are simulated from stochastic FAS and Fourier phase spectrum
(FPS), and are modeled as non-stationary random process in time domain.
%
With the proposed framework, engineering seismologists do not need to
interpret/simplify ground motion into IM(s).
%
Correspondingly, structural engineers do not need to compute fragility curve
based on IM.
%
Instead, all the important characteristics and uncertainties in seismic motions
are captured through the random process and propagated into uncertain
engineering system with direct ``communication'' between engineering
seismologists and structural engineers.
%
Another feature of the proposed framework is the circumvention of Monte Carlo
(MC) simulation.
%
MC approach is non-intrusive in the sense that no modifications to the
underlying deterministic solver are required.
%
The state of probabilistic space is characterized by large, statistically
significant number of deterministic samplings of system random parameters.
%
In conventional seismic risk analysis, structural fragility curve is developed
by incremental dynamic analysis (IDA) \cite{Vamvatsikos2002}.
%
IDA, though theoretically straightforward, is numerically demanding because of
the slow convergence rate that is inherent in MC approach.
%
Hundreds of structural analysis need to be performed with deterministic
sampling of uncertain material properties and uncertain ground excitations at
different IM levels.
%
The same issue of MC approach also limits the application of physics-based
seismic waveform modeling techniques \cite{Graves2010,Maechling2007} into
hazard/risk analysis.
%
Millions of MC earthquake scenarios over regional geology have to be simulated
using deterministic wave propagation programs, such as CyberShake
\cite{Graves2011} considering uncertain kinematic sources, crustal geology and
site conditions.
%
\cite{Maechling2007} estimated that ``it would require 300 million CPU-hours
and well over 100 years to complete all the simulations needed to calculate a
PSHA hazard curve''.
%
% Even with reciprocity technique developed by \citeN{Zhao2006},
% \citeN{Maechling2007} states that ``it is still prohibitively computationally
% expensive to produce a waveform-based PSHA hazard map.''
%
To avoid non-intrusive MC simulation, Galerkin stochastic elastic-plastic finite element method
(SEPFEM) has been developed within the authors' research group over the years
\cite{Jeremic2005a,Sett2005a,Sett2010a,Karapiperis2016,Wang2016,Wang2019}.
%
Galerkin SEPFEM is an intrusive approach, requiring new developments based on
variational formulation of the underlying stochastic partial differential
equations (SPDE).
%
Using appropriate choice of orthogonal polynomial chaos basis, intrusive
Galerkin SEPFEM guarantees optimal convergence rates, and is more efficient than
non-intrusive MC approach \cite{Xiu2010,Elman2011}.
%
Both random field structural parameters and random process seismic motions are
represented by Hermite polynomial chaos (PC) \cite{Sakamoto2002} with
correlation structure characterized by Karhunen-Lo{\`e}ve (KL) expansion
\cite{Zheng2017}.
%
Using Galerkin SEPFEM, probabilistic dynamic response of uncertain structural
system driven by uncertain seismic motions is represented by unknown PC
coefficients.
%
Deterministic linear system equations of these unknown temporal-spatial PC
coefficients, equivalent to the original stochastic PDE, are derived from
Galerkin projection technique in weak sense.
%
Seismic risk is then computed from probabilistic dynamic structural response.
%
The organization of this paper is as follows:
%
The proposed time domain intrusive framework for probabilistic seismic risk
analysis is formulated in section \ref{sec:general_description}.
%
Next, the proposed methodology is illustrated by a numerical example in section
\ref{sec:illustrative_example} with conclusions drawn in section
\ref{sec:conclusions}.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Time Domain Intrusive Framework for Seismic Risk Analysis}
\label{sec:general_description}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The proposed framework consists of four components, as shown in Figure
\ref{figure_risk_framework}: seismic source characterization (SSC), stochastic
ground motion modeling, stochastic finite element analysis and seismic risk
computation.
%
\begin{figure}[!htbp]
\centering
\includegraphics[width=0.75\textwidth]{Figures/risk_framework_illustration.pdf}
\caption{\label{figure_risk_framework} Time domain intrusive framework for
seismic risk analysis.}
\end{figure}
In the first step, SSC quantifies the uncertainty in earthquake scenarios so
that the probabilistic scenario space $\lambda(M, R, \Theta)$ for a given
engineering site can be discretized into $N$ mutually exclusive events as
follows:
%
\begin{equation}
\lambda(M, R, \Theta) = \bigcup_{i=1}^{N} \lambda_i(M_i, R_i, \Theta_i)
\label{eq_discretize_scenario}
\end{equation}
\noindent
where $\lambda(\cdot)$ is the annual occurrence rate, $M$ is the magnitude and
$R$ is the distance metric, which could be either rupture distance $R_{rup}$,
hypocenter distance $R_{hyp}$, or Joyner-Boore distance $R_{jb}$.
%
$\Theta$ denotes any other scenario metrics that are required for stochastic
ground motion modeling, for example, style of fault, hanging wall identifier,
etc.
%
$N$ is the total number of seismic scenarios considering all the active faults
in the region.
%
Basic relations for seismic source characterization are formulated in Section
\ref{sec:SSC}.
%
For each scenario event $S_i(M_i, R_i, \Theta_i)$, section
\ref{sec:stochastic_ground_motion_modeling} presents the procedure to simulate
time domain uncertain motions from stochastic FAS and FPS using inverse Fourier
transform.
%
The simulated ground motion population for event $S_i$ is denoted as
$\{\boldsymbol \Gamma_i\}$.
%
At the third step, both uncertain motions and uncertain structural parameters
are represented by Hermite PC-KL expansion as formulated in section
\ref{sec:pc_kl_expansion}.
%
Two choices are provided here: (1) Random process characterization
(i.e., PC-KL expansion) is performed for each individual motion population $\{\boldsymbol
\Gamma_i\}$ and conduct further Galerkin stochastic FEM analysis for each
scenario $S_i$.
%
(2) Seismic motion population from different scenarios is first combined
as an ensemble population $\{\boldsymbol \Gamma\}$ following Equation \ref{eq_combine_population}:
\begin{equation}
\{ \boldsymbol \Gamma \} = \bigcup_{i=1}^{N} \{ w_i \otimes \boldsymbol \Gamma_i \}
\label{eq_combine_population}
\end{equation}
\noindent
with
\begin{equation}
w_i = \frac{\lambda_i}{\sum_{i=1}^{N} \lambda_i}
\label{eq_weight_population}
\end{equation}
\noindent
where $\bigcup_{i=1}^{N} \{w_i \otimes \boldsymbol \Gamma_i \}$ denotes the weighted combination of population
$\{\boldsymbol \Gamma_i \}$ with weight $w_i$ defined as Equation
\ref{eq_weight_population}.
%
The annual occurrence rate of the ensemble population $\{\boldsymbol \Gamma\}$
is $\overline \lambda = \sum_{i=1}^{N} \lambda_i$.
%
The weighted combination can be performed by aggregating individual population $\{\boldsymbol{\Gamma}_i\}$ of different size $n_i$, $i=1, 2, ..., N$ such that $n_i$ is proportional to $w_i$, i.e., $w_i=n_i/\sum_{i=1}^{N}n_i$.
%
Clearly, size $n_i$ for all $i=1, 2, ..., N$ should be large enough to represent the random process motions from individual seismic scenario.
%
As a result, weighted ensemble population $\{\boldsymbol{\Gamma}\}$ with occurrence rate $\overline{\lambda}$ is statistically equivalent to the aggregation of motion population $\{\boldsymbol{\Gamma}_i\}$ from individual scenario with rate $\lambda_i$.
%
Then the ensemble population $\{\boldsymbol \Gamma\}$ can be characterized as a
single random process and single stochastic FEM analysis is performed with
PC-represented random process motions.
%
Compared with PC-KL representation for each individual population $\{\boldsymbol{\Gamma}_i\}$, the consequence of PC-KL expansion for ensemble population $\{\boldsymbol{\Gamma}\}$ is that larger dimension of PC is required since underlying random process of population $\{\boldsymbol{\Gamma}\}$ is more uncertain and less correlated among different times.
%
If both individual population $\{\boldsymbol{\Gamma}_i\}$ and ensemble population $\{\boldsymbol{\Gamma}\}$ are accurately characterized by PC-KL expansion and propagated into uncertain structure through SFEM, EDP hazard can be calculated by either Equation 5 or Equation 6:
%
% Time-evolving probabilistic structural response is obtained from stochastic FEM analysis.
%
% for chosen EDP(s), failure probability $P_i(EDP>z| \boldsymbol \Gamma_i)$
% conditioned on individual scenario $S_i(M_i, R_i, \Theta_i)$ or $P(EDP>z|
% \boldsymbol \Gamma)$ conditioned on combined scenarios can be determined.
% %
% Finally, EDP hazard can be calculated by Equation
% \ref{eq:time_domain_risk_individual} or \ref{eq:time_domain_risk_combine}
% without using any IM(s) and performing Monte Carlo fragility simulations.
\begin{equation}
\lambda(EDP>z) = \sum_{i=1}^{N}\lambda_i(M_i, R_i, \Theta_i)P(EDP>z|\boldsymbol \Gamma_i)
\label{eq:time_domain_risk_individual}
\end{equation}
%
\begin{equation}
\lambda(EDP>z) = \overline \lambda P(EDP>z|\boldsymbol \Gamma)
\label{eq:time_domain_risk_combine}
\end{equation}
%
where $P_i(EDP>z| \boldsymbol \Gamma_i)$ is the failure probability conditioned on individual population $\{\boldsymbol{\Gamma}_i\}$ and $P(EDP>z| \boldsymbol \Gamma)$ is the failure probability conditioned on ensemble population $\{\boldsymbol{\Gamma}\}$.
%
Both Equation 5 and Equation 6 give consistent result for EDP hazard.
%
The difference is that by using Equation 5, many more less expensive SFEM analyses are performed while using Equation 6 requires a single, yet more expensive SFEM analysis.
%
When the number of scenarios $N$ is small, it is practical to perform stochastic
FEM analysis for each scenario and compute EDP hazard by Equation
\ref{eq:time_domain_risk_individual}.
%
The advantage is that controlling scenario can be identified through EDP hazard
de-aggregation.
%
However, when there are many seismic scenarios, quantifying ensemble population as a single random process through PC-KL expansion and performing single stochastic FEM analysis can be computationally more efficient.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
\subsection{Seismic Source Characterization}
\label{sec:SSC}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Seismic source characterization (SSC) and earthquake rupture forecast (ERF) are
complex scientific issues.
%
Earthquake occurrence rate tends to be comprehensively evaluated by multiple
approaches, for example, using historical seismicity, geological information
(e.g., long term slip rates and paleoseismic recurrence intervals) and geodetic
information \cite{Field2017}.
%
Assuming Poisson process of earthquake occurrence, annual occurrence rate
$\lambda^{f}$ of earthquakes on a fault can be estimated based on seismic moment
balance \cite{McGuire2004}:
%
\begin{equation}
\small
\lambda^{f} = \frac{\mu A \mathcal{S}}{\int_{0}^{M_{max}}E(M)f(M) \, dM}
\label{eq:activity_rate_moment_balancing1}
\end{equation}
%
\noindent
where $\mathcal{S}$ is annual slip rate, $\mu$ is shear modulus of crust and $A$
is fault area, $f(M)$ is the probabilistic model of magnitude distribution, which could be
truncated exponential model, Young's and Coppersmith characteristic model
\cite{Youngs1985}, truncated Gaussian model, etc.
%
The seismic moment of earthquake, $E(M)$ with magnitude $M$ is given as:
%
\begin{equation}
\small
E(M) = 10^{1.5M+16.05}
\label{eq:activity_rate_moment_balancing2}
\end{equation}
In engineering practices, only earthquakes greater than certain magnitude
$M_{min}$ are considered, whose annual occurrence rate $\overline{\lambda}^{f} $
is:
%
\begin{equation}
\overline \lambda^{f} = \lambda^{f} \int_{M_{min}}^{M_{max}}f(M) dM
\end{equation}
%
Using probabilistic models of rupture area conditioned on magnitude $f(A|M)$,
rupture width conditioned on rupture area $f(W|A)$ \cite{Leonard2010}, rupture
location along strike (AS) $f(Y)$ and down-dip (DD) $f(Z)$, distance metric $R$
and other scenario metrics $\Theta$, for example, depth to the top of rupture
plane $Z_{tor}$, can be geometrically characterized as $g(R, \Theta| M)$ for a
given engineering site \cite{Hale2018}.
%
The discretized mutually exclusive scenarios $\lambda_i(M_i, R_i, \Theta_i)$ in
Equation \ref{eq_discretize_scenario} is then quantified as:
%
\begin{equation}
\lambda_i(M_i, R_i, \Theta_i) = \sum_{j=1}^{m} \overline{\lambda_j}^{f} \int_{\boldsymbol{\Lambda_i}} f_j(M) g_j(R, \Theta|M) \,dM\,dR\,d\Theta
\end{equation}
\noindent
where $m$ is the total number of active faults, subscript $j$ denotes the
probabilistic models and quantities specific to the $j^{th}$ fault,
$\boldsymbol{\Lambda_i}$ is the integral domain for the $i^{th}$ discretized
scenario with magnitude step $\Delta M$, distance step $\Delta R$ and $\Delta
\Theta$ for any other scenario metrics $\Theta$ if required:
%
\vspace{-0.4cm}
\begin{equation}
\boldsymbol \Lambda_i = [M_i-\frac{\Delta M}{2}, M_i + \frac{\Delta M}{2}] \times [R_i-\frac{\Delta R}{2}, R_i + \frac{\Delta R}{2}] \times [\Theta_i - \frac{\Delta \Theta}{2}, \Theta_i + \frac{\Delta \Theta}{2}]
\label{eq:intergal_domain}
\end{equation}
%
Many PSHA programs could perform SSC, e.g., HAZ45 \cite{Hale2018}.
%
It is noted that presented above are fundamental relations for seismic characterization of fault sources.
%
% The formulations for area sources are similar and omitted.
%
For regions with unknown fault locations or having background seismicity, areal source should also be considered and characterized.
%
See references \cite{Coppersmith2012,Moschetti2015} for more details on seismic source characterization of areal source.
%
Epistemic uncertainties in slip rate, magnitude distribution models and other
parameters, which are typically considered with logic tree approach
\cite{Musson2012}, are not considered here for simplicity.
%
In addition, for some sites, authoritative estimates of magnitude, location and
rate of earthquake ruptures could be determined from established regional
earthquake rupture forecast (ERF) models, for example, UCERF3 \cite{Field2017}
for California region.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%\su
\subsection{Time Domain Stochastic Ground Motion Modeling}
\label{sec:stochastic_ground_motion_modeling}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Stochastic modeling of FAS and FPS are two indispensable components to simulate
% time domain uncertain motions.
%
Time domain uncertain motions can be simulated from stochastic FAS and Fourier phase derivative \cite{Boore2003b,Boore2003}.
%
Specifically, uncertain FAS of seismic motions is modeled as Log-normal
distributed random field \cite{Bora2015,Stafford2017b} in frequency space,
whose marginal median behavior is simulated by the stochastic method of
\citeN{Boore2003}.
%
It is referred to as Boore03 approach hereafter.
%
Boore03 approach simulates FAS using $w^2$ radiated source spectrum
\cite{Brune1970} with modification for path and site effects, as shown in
Equation \ref{eq:FAS_SIMSM}:
%
\begin{equation}
FAS(f) = A_0(M_0, f)Z(R)exp(-\pi f R/Q\beta)S(f)exp(-\pi \kappa_0 f)
\label{eq:FAS_SIMSM}
\end{equation}
%
\noindent
where $M_0$ is the seismic moment; $\beta$ is the source shear wave velocity;
%
$Z(R)$ and $exp(-\pi f R/Q\beta)$ represent the contribution from path effects:
$Z(R)$ is the geometrical spreading term as a function of distance $R$.
%
Term $exp(-\pi f R/Q\beta)$ quantifies the anelastic attenuation as the inverse of
the regional quality factor, $Q$.
%
The site effects including site amplification through crustal velocity gradient
and near surface attenuation are demonstrated by $S(f)$ and
$\kappa_0$ filter $exp(-\pi \kappa_0 f)$, respectively.
%
Term $A_0$ represents the radiated acceleration source spectrum, which
could be characterized by single-corner-frequency model:
% ) or double-corner-frequency model \cite{Boore2014b}.
%
\begin{equation}
A_0(M_0, f) = CM_0\bigg[\frac{(2\pi f)^2}{1+(f/f_0)^2}\bigg]
\label{eq:single_corner_frequency_model}
\end{equation}
\noindent
where $f_0$ is the corner frequency, which in Brune's model \cite{Brune1970} is
related to source stress drop $\Delta \sigma$ as follows:
%
% $C$ is a scaling constant, typically taken as $R_{\Theta \phi}VF/(4\pi\rho\beta^3)$, in which $R_{\Theta \phi}$ is the average radiation coefficient 0.55 for shear waves \cite{Boore1984}, $V=1/\sqrt{2}$ is the partitioning of seismic energy into horizontal components.
%
% $F=2$ is factor for free surface reflection.
\begin{equation}
f_0 = 4.9\times10^6\beta (\Delta \sigma/M_0)^{1/3}
\end{equation}
Boore03 approach is well-recognized for its simplicity and effectiveness to
capture the marginal mean behavior of stochastic FAS.
%
% It is noted that several GMPEs of FAS have been derived recently
% \cite{Bora2015,Bora2018,Bayless2018b} and are also suitable to model the
% marginal behavior of FAS within the proposed framework.
%
\citeN{Bayless2018c} pointed out that the inter-frequency correlation structure
of FAS random field is also important for seismic risk analysis.
%
Misrepresenting-representing the correlation structure, e.g., assuming
inter-frequency independence, would lead to underestimation of seismic risk.
%
% Recently, inter-frequency correlation models for $\epsilon$, which is defined as residual of uncertain $FAS(f)$ with respect to its median estimate normalized by standard derivation, have been developed \cite{Stafford2017b,Bayless2017}.
%
Therefore, inter-frequency correlation model for stochastic FAS developed
recently \cite{Stafford2017b,Bayless2017} is adopted here.
%
% Adopted here is the FAS correlation structure identified by \cite{Bayless2017}.
%
Though the behavior of FAS was well studied, modeling Fourier phase
angles is still challenging.
%
% In the original Boore03 approach \cite{Boore2003}, random phase info is
% contained in the stationary Gaussian white noise modulated by an envelope
% function.
%
Conventionally random phase info is simulated using stationary Gaussian white noise modulated by an envelope function.
%
However, \citeN{Montaldo2003} stated that conventional Gaussian white noise approach could not reliably reproduce the non-stationarity of ground motions.
%
For this reason, the use of phase difference $\Delta \Phi$ was suggested by \citeN{Ohsaki1979}.
%
Using California strong ground motion data, \citeN{Thrainsson2002} modeled phase
differences as Beta distribution.
%
However, the established phase difference models are affected by the signal
length of each record.
%
It is more stable to normalize phase difference by signal length and study the
probabilistic model of phase derivative $\dot{\Phi}$ defined as
\cite{Boore2003b} :
%
\begin{equation}
\dot{\Phi} = \frac{\Delta \Phi}{\Delta f}
\label{eq:phase_derivative_model}
\end{equation}
%
Based on 3551 ground motion records from PEER NGA-West 1 database,
\citeN{Baglio2017} found that the distribution of phase derivative is
leptokurtic and fits well to Logistic model:
%
%%%% new equation
\begin{equation}
f(\dot{\Phi}; \mu, \sigma) = \frac{1}{4\sigma} sech^{2}(\frac{\dot{\Phi}-u}{2\sigma})
\label{eq:phase_derivative_logistic_model}
\end{equation}
\noindent
where $\mu$ and $\sigma$ are the mean and scale parameter of the Logistic
distribution $f(\dot{\Phi}; \mu, \sigma)$, $sech(\cdot)$ is the hyperbolic
secant function.
%
Following \citeN{Baglio2017}, the mean value $\mu$ is a fixed parameter to
position the distribution along the signal length.
%
For example, setting mean parameter $\mu$ equal to $\pi/df$ would align the peak of uncertain seismic motions to the center of simulated signal length.
%
The prediction equation of scale parameter $\sigma$ is correlated to earthquake
magnitude $M$, rupture distance $R_{rup}$, $V_{s30}$ and directivity index
$D_{Dir}=R_{hyp} - R_{rup}$ with coefficients $\alpha_1$, $\alpha_2$, $\beta_1$
$\sim$ $\beta_4$, $\gamma_1$ and $\gamma_2$ determined from maximum likelihood
estimation:
\begin{equation}
log(\sigma/\pi) = \alpha_1 + \alpha_2 log[\beta_1 + 10^{\beta_2M} + \beta_3 R_{rup} + \beta_4 log(V_{s30})+ \gamma_1 + \gamma_2 D_{Dir}]
\label{eq:GMPE_phase_derivative}
\end{equation}
Phase derivatives $\dot{\Phi}(f)$ among frequency coordinates is modeled as
Logistic distributed random field following exponential correlation
with correlation length $l_f=0.05Hz$:
%
\begin{equation}
\displaystyle
{Cov(\dot{\Phi}(f_1), \dot{\Phi}(f_2)) = e^{-\displaystyle{\frac{|f_1-f_2|}{l_f}}}}
\label{eq:correlation_structure_phase_derivative}
\end{equation}
% \begin{equation}
% \displaystyle
% {Cov(\dot{\Phi}(f_1), \dot{\Phi}(f_2)) = e^{-\frac{|f_1-f_2|}{l_f}}}
% \label{eq:correlation_structure_phase_derivative}
% \end{equation}
%
The methodology of time domain stochastic ground motion modeling is summarized
below:
%
\begin{enumerate}
\item Compute marginal median of Log-normal distributed random field $FAS(f)$
following Boore03 approach.
%
\item Generate realizations of Log-normal distributed random field $FAS(f)$
according to the marginal estimation in step 1 and inter-frequency correlation
model by \citeN{Bayless2017}.
%
\item Determine the scale parameter $\sigma$ of marginal Logistic model for
phase derivative random field with Equation \ref{eq:GMPE_phase_derivative}.
%
Set mean value $\mu$ to $\pi/df$ for central peak \cite{Baglio2017}.
%
\item Generate realizations of Logistic distributed random field $\dot{\Phi}(f)$
with marginal distribution from step 3 and exponential correlation structure.
%
\item Multiply realization of phase derivative in step 4 by frequency interval
$df$ to get realizations of phase difference $\Delta \Phi(f)$.
%
Compute realizations of phase angles from $\Delta \Phi(f)$. First phase angle
can be randomly set and would not affect the synthesis.
%
\item Inverse Fourier transformation to time domain with realizations of FAS and
FPS.
\end{enumerate}
Time domain stochastic ground motion for a single scenario with magnitude $M=7$,
distance $R_{rup}=15km$ has been simulated following the above methodology.
%
The detailed modeling parameters for source, path and site effects of this
scenario are given in Table
\ref{table:stochastic_ground_motion_modeling_parameter}.
%
The marginal median FAS is computed using program SMSIM developed by
\citeN{Boore2005}.
%
% Finite faults effects are considered through equivalent point source distance $R_{PS}$ following \cite{Boore2015}.
%
With reference to recent GMPE studies of FAS
\cite{Bora2015,Bora2018}, marginal lognormal standard deviation of
FAS has been adopted as total $\sigma=0.8 \,ln$ units.
%
The maximum modeling frequency $f_{max}$ is 20Hz.
%
It is noted that ergodic assumption was used in developing these GMPEs of FAS.
%
A smaller value of marginal standard deviation can be used for non-ergodic probabilistic seismic risk analysis if additional source, path or site specific information is available.
%
\begin{table}[!htbp]
\footnotesize
\centering
% \caption{\label{table:stochastic_ground_motion_modeling_parameter} Input parameters for median FAS modeling with SMSIM \cite{Boore2005}}
\caption{\label{table:stochastic_ground_motion_modeling_parameter} Source, path
and site parameters for stochastic ground motion modeling of seismic scenario
$M=7$, $R_{rup}=15km$.}
\begin{tabular}{ccc}
\hline
\multicolumn{1}{l}{Parameter type} & Name & Value \\ \hline
\multirow{6}{*}{Source} & Magnitude & M=7 \\
& Source density & $\rho_s=2.8g/cm^3$ \\
& Source velocity & $\beta =3.6km/s$ \\
& $w^2$ source spectrum & Single corner frequency with $\Delta \sigma=8.0 MPa$ \\
& Fault type & Reverse fault $F_{rv}=1$ \\
& Dip angle & $45^{\circ}$ \\ \hline
\multirow{4}{*}{Path} & Distance metrics & \begin{tabular}[c]{@{}c@{}}$R_{rup}=15km$, $R_{hyp}=18km$\\ $R_{jb}=12km$, $R_x=-12km$\end{tabular} \\
& Finite faults effects & \begin{tabular}[c]{@{}c@{}}Equivalent point source model \cite{Boore2015} \\ with $R_{PS}=22.18km$\end{tabular} \\
& Geometrical spreading & \begin{tabular}[c]{@{}c@{}}Hinged line segments model \cite{Atkinson1995} \end{tabular} \\
& Anelastic attenuation Q & Three line segments model by \cite{Boore2003} \\ \hline
\multirow{2}{*}{Site} & Site amplification & \begin{tabular}[c]{@{}c@{}}$V_{s30}=620m/s$\\ Table 4 of \cite{Boore2015} \end{tabular} \\
& $\kappa_0$ attenuation & $\kappa_0=0.03s$ \\ \hline
\end{tabular}
\end{table}
Combining stochastic FAS with uncertain Fourier phase info, 500 realizations of time domain stochastic motions are synthesized.
%
Figure \ref{figure_acc_realizations} shows three different synthesized accelerations.
%
Large variability is observed, for example, peak ground acceleration could vary
from 1.8$m/s^2$ to 5.5$m/s^2$.
\begin{figure}[!htbp]
\centering
\includegraphics[width=0.31\textwidth]{Figures/Acc_time_series_realiztion70.pdf}
\includegraphics[width=0.29\textwidth]{Figures/Acc_time_series_realiztion100.pdf}
% \includegraphics[width=0.31\textwidth]{Figures/Acc_time_series_realiztion350.pdf}
\includegraphics[width=0.325\textwidth]{Figures/Acc_time_series_realiztion367.pdf}
\caption{\label{figure_acc_realizations} Realizations of uncertain acceleration
time series population.}
\end{figure}
Spectral acceleration ($Sa$) of 500 synthesized realizations are calculated and
compared with weighted average prediction of five NGA West-2 GMPEs
\cite{Gregor2014} with weights 0.22 for ASK14, 0.22 for BSSA14, 0.22 for CB14,
0.22 for CY14 and 0.12 for I14.
%
From figure \ref{figure_SA_verification_GMPE}(a), median spectral acceleration
$Sa$ from simulated stochastic ground motion is in very good agreement with GMPE
predictions for all period ordinates.
%
No systematic bias is observed.
%
Figure \ref{figure_SA_verification_GMPE}(b) shows that the standard deviation of simulated response spectra is around 0.65 $ln$ units, which is
consistent with aleatory variability of $Sa$ given by GMPEs.
%
In other words, time domain stochastic ground motions simulated with
aforementioned methodology could not only characterize the median behavior of
$Sa$ very well, but also carry desired amount of uncertainty that is consistent
with empirical GMPEs.
\begin{figure}[!htbp]
\centering
\subfloat[GMPE verification of $Sa$]{
\includegraphics[width=0.51\textwidth]{Figures/SA_GMPE_verification_std_08_no_smooth.pdf}}
\subfloat[Bias check]{
\hspace{-0.3cm}
\includegraphics[width=0.5\textwidth]{Figures/Standard_deviation_std_08_no_smooth.pdf}}
\caption{\label{figure_SA_verification_GMPE} Verification of simulated
stochastic motions: (a) Check median spectral acceleration $Sa$ with NGA West-2 GMPE
(b) Check aleatory variability of simulation spectral acceleration $Sa$.}
\end{figure}
%
The marginal distribution of simulated accelerations at all the time instances is observed to be Gaussian.
%
Similar observation is also made by \citeN{Wang2016} from statistical analysis of seismic records.
%
Therefore, time domain stochastic ground motions is modeled as a Gaussian distributed non-stationary random process.
%
The random process would be represented with Hermite polynomial chaos as
formulated in the next section.
%
It is worthwhile to mention that the random process incorporates much more
information about uncertain ground motions than GMPE used in conventional PBEE.
%
GMPE only quantifies the variability of selected IM, such as $Sa$, while the
random process carries not only consistent variability of $Sa$ but also any
other important characteristics, e.g., PGA, CAV, etc.
%
Realistic inter-frequency correlation of FAS \cite{Bayless2017} is captured.
%
Non-stationarity of seismic motions is quantified through phase derivative modeling without using any modulation function.
%
Compared with existing ground motion modeling techniques commonly adopted by
reliability community, e.g., evolutionary power spectrum and white noise random
phase spectrum, presented methodology is directly compatible with
state-of-the-art seismic source characterization.
%
It could explicitly account for specific source, path and site condition in
both stochastic modeling of FAS and FPS.
%
Many reliability analysis methods, such as probabilistic density evolution
method \cite{Xu2019} can be readily combined with the presented methodology and
incorporated into the proposed risk analysis framework for PBEE.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Hermite Polynomial Chaos Karhunen-Lo{\`e}ve expansion}
\label{sec:pc_kl_expansion}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
This section formulates Hermite polynomial chaos Karhunen-Lo{\`e}ve (PC-KL)
expansion for general heterogeneous random field $D(\boldsymbol{x}, \theta)$ of
arbitrary marginal distribution.
%
Both uncertain motions and uncertain structure parameters can be represented
with PC-KL expansion.
%
To achieve this, we first represent heterogeneous random field
$D(\boldsymbol{x}, \theta)$ of arbitrary marginal distribution through Hermite
polynomial chaos of underlying Gaussian heterogeneous random field
$\gamma(\boldsymbol{x}, \theta)$ up to order $P$ \cite{Sakamoto2002,Wang2016}:
%
% The populations of simulated time series seismic motions are regarded as
% realizations of underlying non-stationary random process.
%
% This random process can be quantified using Polynomial Chaos Karhunen-Lo{\`e}ve (PC-KL) expansion.
%
% Any heterogeneous random field of arbitrary marginal
% distribution (Gaussian or non-Gaussian) can be represented through Hermite polynomial chaos (PC) up to
% order $P$ as \cite{Sakamoto2002,Wang2016}:
%
\begin{equation}
D(\boldsymbol{x}, \theta)
=
\sum_{i=0}^{P} D_i(\boldsymbol{x})\Omega_i(\gamma(\boldsymbol{x}, \theta))
\label{eq:PC_expansion}
\end{equation}
%
\noindent
% where $\boldsymbol{x}$ is the general coordinate in
% space $V$ that can span space and/or time.
%
where $\theta$ denotes the uncertainties.
%
Functions $\{\Omega_i\}= \{1, \gamma, \gamma^2-1, \gamma^3-3\gamma, ...\}$ are
orthogonal, zero mean ($i\geq 1$) Hermite PC bases constructed from zero mean, unit variance kernel Gaussian
random field $\gamma(\boldsymbol{x}, \theta)$.
%
Then at the second step, Gaussian random field $\gamma(\boldsymbol{x}, \theta)$ can be further decomposed by Karhunen-Lo{\`e}ve
(KL) theorem \cite{Zheng2017}.
%
The deterministic PC coefficient field $D_i(\boldsymbol{x})$ can be calculated
through marginal distribution of $D(\boldsymbol{x}, \theta)$, as shown
in Equation
\ref{eq:PC_coefficient}, where $\braket{\cdot}$ is the expectation operator.
%
\begin{equation}
D_i = \frac{\braket{D\Omega_i}}{\braket{\Omega_i^2}}
\label{eq:PC_coefficient}
\end{equation}
%
% The main idea behind Equation \ref{eq:PC_expansion} is to represent arbitrary heterogeneous random field of any type of marginal distribution with a kernel Gaussian random field.
% %
% In this way, deterministic PC coefficients field can be determined and
%
The covariance structure of the original random field $Cov_D(x_1, x_2)$ is
mapped to the Gaussian covariance kernel $Cov_{\gamma}(x_1, x_2)$ as:
%
\begin{equation}
Cov_D(x_1, x_2) = \sum_{i=1}^{P} D_i(x_1) \; D_i(x_2)\; i\,! \; Cov_{\gamma}(x_1, x_2)
\label{eq:correlation_transformation} \;
\end{equation}
%
The Gaussian covariance kernel $Cov_{\gamma}(x_1, x_2)$ can be eigen-decomposed
into probabilistic spaces up to dimension $M$, according to Karhunen-Lo{\`e}ve
(KL) theorem \cite{Zheng2017}:
%
\begin{equation}
\gamma(\boldsymbol{x}, \theta)
=
\sum_{i=1}^{M} \sqrt{\lambda_i} f_i(\boldsymbol{x}) \xi_i(\theta)
\label{eq:KL_expansion}
\end{equation}
%
%
where $\{\xi_i(\theta)\}$ are the multidimensional, orthogonal, zero mean and
unit variance Gaussian random variables, and
$\lambda_i$ and $f_i(\boldsymbol{x})$ are the eigen-values and eigen-vectors of
the covariance kernel $Cov_{\gamma}(x_1, x_2)$ that satisfy Fredholm's integral
equation of the second kind \cite{Sakamoto2002}.
%
% \begin{equation}
% \int_V Cov_{\gamma}(x_1, x_2) f_i(x_1) \, dx_1 = \lambda_{i} f_i(x_2)
% \end{equation}
%
% Due to the unit variance of Gaussian random field $\gamma(\boldsymbol{x},
% \theta)$, eigen-values $\lambda_i$. and eigen-vectors $f_i(\boldsymbol{x})$,
% of
% the covariance kernel $Cov_{\gamma}(x_1, x_2)$,
% are constrained as follows:
% %
% \begin{equation}
% \sum_{i=1} (\sqrt{\lambda_i} f_i(\boldsymbol{x}))^2 = 1
% \label{eq:constrained_KL_condition}
% \end{equation}
% %
Combining Equations \ref{eq:PC_expansion} and \ref{eq:KL_expansion}, the
resultant PC-KL representation of general random field $D(\boldsymbol{x},
\theta)$ is obtained as,
%
%
\begin{equation}
D(\boldsymbol{x}, \theta) = \sum_{i=0}^{K} d_i(\boldsymbol{x}) \Psi_i(\{\xi_j(\theta)\})
\label{eq:PC-KL expansion}
\end{equation}
%
%
% \noindent
% {\huge \bf NOTE:
% where $\{\Psi_i\}$ are multi-dimensional orthogonal Hermite PC bases of order
% $P$, dimension $M$, probabilistic space (i.e., \{$\xi_j(\theta)$\}).}
% %
% %HEXIANG upper sentence needs to be rewritten, last part is not clear
%
%===== Rephrased as =====
where $\{\Psi_i\}$ are multi-dimensional orthogonal Hermite PC bases of order
$P$ constructed from $M$ dimensional probabilistic space (i.e., \{$\xi_j(\theta)$\}, $j=1, 2, ..., M$).
%
%
The total number of multidimensional Hermite PC bases $K$ is related to order $P$
and dimension $M$ as $K = 1+\sum_{s=1}^{P}\frac{1}{s!}\prod_{j=0}^{s-1}(M+j)$.
%
% \begin{equation}
% \displaystyle K = 1+\sum_{s=1}^{P}\frac{1}{s!}\prod_{j=0}^{s-1}(M+j)
% \end{equation}
%
By equating two representations of $D(\boldsymbol{x}, \theta)$ in Equations
\ref{eq:PC_expansion} and \ref{eq:PC-KL expansion}, the coefficients of
multi-dimensional Hermite PC are derived as:
%
% \begin{equation}
% d_i(\boldsymbol{x}) = \frac{P!}{\braket{\Psi_i^2}} D_P(\boldsymbol{x}) \prod_{j=1}^{P} \frac{\sqrt{\lambda_{k(j)}}f_{k(j)}(\boldsymbol{x})}{\sqrt{\sum_{m=1}^{M}(\sqrt{\lambda_m}f_m(\boldsymbol{x}))^2}}
% \end{equation}
%
\begin{equation}
d_i(\boldsymbol{x}) = \frac{p!}{\braket{\Psi_i^2}} D_p(\boldsymbol{x}) \prod_{j=1}^{p} \frac{\sqrt{\lambda_{k(j)}}f_{k(j)}(\boldsymbol{x})}{\sqrt{\sum_{m=1}^{M}(\sqrt{\lambda_m}f_m(\boldsymbol{x}))^2}}
\end{equation}
%
where $p$ is the order of the polynomial $\Psi_i$.
%
From Equation \ref{eq:PC-KL expansion}, PC synthesized marginal mean and variance of the original heterogeneous random field can be calculated as:
%
\begin{equation}
\braket{D(\boldsymbol{x}, \theta)} = d_0(\boldsymbol{x})
\label{eq:PC_marginal_mean}
\end{equation}
%
\begin{equation}
Var(D(\boldsymbol{x}, \theta)) = \sum_{i=1}^{K} d_i^2(\boldsymbol{x}) \braket{\Psi_i^2}
\label{eq:PC_marginal_variance}
\end{equation}
%
PC-synthesized correlation structure can also be computed as:
\begin{equation}
\displaystyle
Cov_D(x_1, x_2) = \frac{\sum_{i=1}^{K} d_i(x_1)d_i(x_2)\braket{\Psi_i^2}}{\sqrt{Var(D(x_1))Var(D(x_2))}}
\label{eq:PC_correlation_structure}
\end{equation}
%
Equations \ref{eq:PC_marginal_mean}, \ref{eq:PC_marginal_variance} and
\ref{eq:PC_correlation_structure} can be used to compare the PC-synthesized
statistics with statistics of original random field $D(\boldsymbol{x}, \theta)$
and check the goodness of PC-KL expansion.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Galerkin Stochastic Finite Element Method}\label{sec:SFEM}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Stochastic Galerkin approach intrusively solves the stochastic PDE with optimal convergence \cite{Ghanem1996,Sett2010a,Wang2016}.
%
% This is different from the non-intrusive methods, such as Monte Carlo and stochastic collocation \cite{Babuska2007b}.
%
Stochastic Galerkin approach intrusively solves the stochastic partial
differential equations (PDE) with optimal convergence
\cite{Sett2010a,Wang2016}.
%
Compared to deterministic finite element method (FEM), Galerkin stochastic
FEM introduces spectral discretization of probabilistic domain in addition to
the spatial and temporal discretization.
%
Using standard spatial FEM discretization, unknown displacement random field
$u(\boldsymbol{x}, t, \theta)$ can be expressed with shape function
$N_i(\boldsymbol{x})$ and uncertain displacement $u_i(t, \theta)$ at nodes:
%
\begin{equation}
u(\boldsymbol{x}, t, \theta) = \sum_{i=1}^{N} N_i(\boldsymbol{x}) u_i(t, \theta)
\label{eq:spatial_discretization}
\end{equation}
%
Uncertain displacement at nodes $u_i(t, \theta)$, can be further represented with
aforementioned multidimensional Hermite PC basis $\phi_j(\{\xi_r(\theta)\})$ of
dimension $M^{u}$, order $P^{u}$:
%
\begin{equation}
u_i(t, \theta) = \sum_{j=0}^{K^{u}} u_{ij}(t) \phi_j(\{\xi_r(\theta)\})
\label{eq:probabilistic_discretization}
\end{equation}
%
Combining Equations \ref{eq:spatial_discretization} and
\ref{eq:probabilistic_discretization}, spatial-probabilistic discretized
expression of $u(\boldsymbol{x}, t, \theta)$ is given:
%
\begin{equation}
u(\boldsymbol{x}, t, \theta) = \sum_{i=1}^{N}\sum_{j=0}^{K^{u}} N_i(\boldsymbol{x}) u_{ij}(t) \phi_j(\{\xi_r(\theta)\})
\end{equation}
%
Galerkin weak formulation of stochastic partial differential equations of motion
can then be written in the following form:
%
\begin{center}
\begin{eqnarray}
\sum_e \bigg[
& &
\int_{D_e} N_m(\boldsymbol x) \rho(\boldsymbol x) N_n(\boldsymbol x) d\,V \, \ddot{u}_n(t, \theta)
+
\nonumber
\\
&+&
\int_{D_e} B_m(\boldsymbol{x}) E(\boldsymbol{x, \theta}) B_n(\boldsymbol x) d\,V\,u_n(t, \theta) - f_m(t, \theta) \bigg]
=
0
% \end{split}
\label{eq:weak_form_FEM}
\end{eqnarray}
\end{center}
%
where $\displaystyle \sum_e$ denotes the assembly procedure over all finite
elements, while $\rho(\boldsymbol{x})$ is deterministic material density field.
%
The shape function gradient function $B_n(\boldsymbol x)$ is given as:
%
\begin{equation}
B_n(\boldsymbol x) = \nabla N_n(\boldsymbol x)
\end{equation}
%
In Equation~\ref{eq:weak_form_FEM},
$E(\boldsymbol{x}, \theta)$ is uncertain tangential stiffness matrix, while
$f_m(t, \theta)$ is uncertain global force vector that incorporates various
elemental contributions.
%
% This uncertain force vector can be used to input uncertain seismic excitations,
% using Domain Reduction Method \cite{Bielak2001,Yoshimura2001}.
% %
Expansion of uncertain stiffness matrix $E(\boldsymbol{x}, \theta)$, and
uncertain force vector $f_m(t, \theta)$ into Hermite PC bases
$\Psi_k(\{\xi_r(\theta)\})$ and $\psi_l(\{\xi_r(\theta)\})$ of dimension
$M^{E}$, order $P^{E}$ and dimension $M^{f}$, order $P^{f}$, respectively, yields:
%
%
\begin{equation}
E(\boldsymbol{x}, \theta)
=
\sum_{k=0}^{K^{E}} E_k(\boldsymbol{x}) \Psi_k(\{\xi_r(\theta)\})
\label{eq:stiffness_discretization}
\end{equation}
%
%
\begin{equation}
f_m(t, \theta)= \sum_{l=0}^{K^{f}} f_{ml}(t) \psi_l(\{\xi_r(\theta)\})
\label{eq:force_discretization}
\end{equation}
%
By combining equations \ref{eq:probabilistic_discretization},
\ref{eq:stiffness_discretization} and \ref{eq:force_discretization} and
equation \ref{eq:weak_form_FEM}, one obtains:
%
%
\begin{equation}
\begin{split}
& \sum_e \bigg[ \int_{D_e} N_m(\boldsymbol x) \rho(\boldsymbol x) N_n(\boldsymbol x) d\,V \, \sum_{j=0}^{K^u} \ddot{u}_{nj} \phi_j(\{\xi_r(\theta)\} - \sum_{l=0}^{K^{f}} f_{ml} \psi_l(\{\xi_r(\theta)\}) \\
& + \int_{D_e} B_m(\boldsymbol{x}) \sum_{k=0}^{K^{E}} E_k(\boldsymbol{x}) \Psi_k(\{\xi_r(\theta)\}) B_n(\boldsymbol x) d\,V\, \sum_{j=0}^{K^u} u_{nj} \phi_j(\{\xi_r(\theta)\}\bigg] = 0
\end{split}
\label{eq:weak_formulatio_expand}
\end{equation}
%
By performing Galerkin projection of Equation \ref{eq:weak_formulatio_expand}
onto PC bases $\phi_i(\{\xi_r(\theta)\})$, to minimize the residual, system of
deterministic ordinary differential equations (ODE) involving temporal
derivative of unknown PC coefficients $u_{nj}$, is developed:
% with Einstein notation:
%
\begin{equation}
M_{minj} \ddot{u}_{nj} + K_{minj} u_{nj} = F_{mi}
\label{eq:deterministic_ODE}
\end{equation}
%
where mass tensor/matrix $M_{minj}$ is given by equation~\ref{eq:Mass_stochastic_matrix}:
%
% {\huge \bf NOTE:
% Hexiang, since material density is determinstic, this mass matrix is then
% deterministic! Is this right?}
% %
%%%%%% === Hexiang Reply: Yes, because of the determinisitc material density, the expanded mass matrix is only determinisitc FEM mass matrix multiplied by the double product of PC bases. And because of the orthogonality of PC bases, it ends up with diagonally duplicate of deterministic mass matrix.
\begin{equation}
M_{minj} = \sum_{e} \int_{D_e} N_m(\boldsymbol x) \rho(\boldsymbol x) N_n(\boldsymbol x) d\,V \, \braket{\phi_i \phi_j}
\label{eq:Mass_stochastic_matrix}
\end{equation}
%
stochastic stiffness tensor/matrix $K_{minj}$ is given by equation~\ref{eq:stiffness_stochastic_matrix}:
%
\begin{equation}
K_{minj} = \sum_{k=0}^{K^{E}} \sum_{e} \int_{D_e} B_m(\boldsymbol{x}) E_k(\boldsymbol{x}) B_n(\boldsymbol x) d\,V\, \braket{\Psi_k \phi_i \phi_j}
\label{eq:stiffness_stochastic_matrix}
\end{equation}
%
and stochastic force tensor/vector $F_{mi}$ is given by
equation~\ref{eq:force_stochastic_matrix}
\begin{equation}
F_{mi} = \sum_{l=0}^{K^{f}} f_{ml} \braket{\psi_l \phi_i}
\label{eq:force_stochastic_matrix}
\end{equation}
In Equations \ref{eq:Mass_stochastic_matrix},
\ref{eq:stiffness_stochastic_matrix} and \ref{eq:force_stochastic_matrix},
terms $\braket{\phi_i \phi_j}$, $\braket{\psi_l \phi_i}$ and
$\braket{\Psi_k\phi_i\phi_j}$ are the ensemble average tensors of double-product
and tri-product of different PC bases.
%
These ensemble average tensors could be pre-computed and used to construct the
stochastic mass matrix $M_{minj}$ and stochastic stiffness matrix $K_{minj}$.
%
It is noted that Einstein's notation for
tensor indices summation is assumed throughout \cite{local-63}.
%
The deterministic system of ordinary differential equations (ODE) from Equation
\ref{eq:deterministic_ODE}, can be integrated in time using
dynamic integrator algorithms, for example Newmark method \cite{Newmark1959},
or Hilber-Hughes-Taylor $\alpha$-method \cite{Hilber1977}.
%
Result of such time marching solution will be time histories of displacement PC
coefficients $u_{nj}$.
%
Those time evolving displacement PC coefficients $u_{nj}$ can then be used
to develop complete probabilistic dynamic finite element response.
%
With resulting complete probabilistic dynamic finite element response,
any damage measure, in fact all damage measures related to EDP(s) can be
applied to trace the failure probability $P_i(EDP>z|\boldsymbol{\Gamma}_i)$ or
$P(EDP>z|\boldsymbol{\Gamma})$.
%
EDP hazard can then be computed according to Equations
\ref{eq:time_domain_risk_individual} and \ref{eq:time_domain_risk_combine}.
The above formulation of Galerkin stochastic FEM is complete for linear elastic
problem with constant uncertain elastic stiffness matrix $E(\boldsymbol{x},
\theta)$.
%
For nonlinear, inelastic problems, additional formulation of stochastic
elastic-plastic FEM (SEPFEM) is required and relies on recent developments
\cite{Jeremic2005a,Sett2005a,Sett2009b,Sett2010a,Arnst2012,Rosic2014,Karapiperis2016}.
%
One of the challenges of the SEPFEM lies in the development of
the probabilistic elastic-plastic stiffness at the constitutive level that is to be
used for finite element level computations.
%
Eulerian-Lagrangian form of the Fokker-Planck-Kolmogorov (FPK) equation has been
successfully used to obtain probabilistic stress solutions
\cite{Jeremic2005a,Sett2005a,Sett2010a}.
%
It is noted in order to produce uncertain stiffness, least square
optimization and linearizion techniques
\cite{Sett2010a,Karapiperis2016} are used.
%
To this end, in one dimension (1D), elastic plastic material
model with vanishing elastic region is used in conjunction with
Armstrong-Fredrick nonlinear kinematic hardening \cite{Armstrong66,Dettmer2004}.
%
This approach simplifies modeling, as elastic plastic response directly follows
Armstrong-Frederick nonlinear equation.
%
% For approach proposed here, direct modeling of nonlinear, probabilistic inelastic
% inter-story drift, the probabilistic nonlinear response between
% inter-story restoring force
% $F^{R}$ and inter-story drift $\eta$ is formulated through direct PC-based
% Galerkin intrusive
% probabilistic modeling of Armstrong-Fredrick hysteretic behavior.
%
For the approach proposed here, probabilistic nonlinear response between inter-story restoring force $F^{R}$ and inter-story drift $\eta$ is formulated through direct PC-based Galerkin intrusive probabilistic modeling of Armstrong-Fredrick hysteretic behavior.
In incremental form, Armstrong-Fredrick kinematic hardening relation
\cite{Armstrong66} between inter-story restoring force $F^{R}$ and
inter-story drift $\eta$ can
be written as:
\begin{equation}
dF^{R} = H_a\, d\eta - C_r F^{R} |d\eta|
\label{eq:AF_constitutive}
\end{equation}
%
where $H_a$ and $C_r$ are model parameters.
%
By setting $dF^{R} =0$, the ultimate inter-story restoring force
becomes $F^{R}_{max} =H_a/C_r$.
%
The tangential stiffness $E(F^{R})$ is a function of restoring force
$F^{R}$:
%
\begin{equation}
E(F^{R}) = \frac{d F^{R}}{d \eta} = H_a - C_r F^{R}\,sgn(d\eta)
\label{eq:tangential_stiffness}
\end{equation}
%
where $sgn(\cdot)$ is the sign function. Equation~\ref{eq:tangential_stiffness}
can be written as:
\begin{equation}
E(F^{R}) = H_a \pm C_r F^{R}
\label{eq_simplified_tangential_stiffness}
\end{equation}
%
where $+$ sign is taken for negative inter-story drift $d\eta$ and $-$ sign is
taken for positive inter-story drift $d\eta$.
%
In the general probabilistic setting, model parameters $H_a$ and $C_r$ can be
uncertain and
modeled as random fields $H_a(\boldsymbol x, \theta)$ and $C_r(\boldsymbol x,
\theta)$.
%
By applying PC expansion with Hermite PC bases $\varphi_i(\{ \xi_r(\theta)\})$ to
those two model parameters, the following equations are obtained:
%
%
\begin{equation}
H_a(\boldsymbol{x}, \theta) = \sum_{i=0}^{P} H_{ai}(\boldsymbol x) \varphi_i(\{ \xi_r(\theta)\})
\label{eq:Ha_expanded}
\end{equation}
%
%
\begin{equation}
C_r(\boldsymbol{x}, \theta) = \sum_{i=0}^{P} C_{ri}(\boldsymbol x) \varphi_i(\{ \xi_r(\theta)\})
\end{equation}
%
%
The inter-story drift increments $d \eta(\boldsymbol{x}, \theta)$, that
represent input to the to constitutive driver (Equation
\ref{eq:AF_constitutive}) are also uncertain due to the probabilistic
structural response $u(\boldsymbol x, t, \theta)$:
%
%
\begin{equation}
d \eta(\boldsymbol{x}, \theta) = \sum_{i=0}^{P} d \eta_i(\boldsymbol x) \varphi_i(\{ \xi_r(\theta)\})
\end{equation}
%
As a result, probabilistic incremental restoring force $dF^{R}(\boldsymbol x,
\theta)$ and probabilistic tangential stiffness $E(\boldsymbol{x}, \theta)$ are
then:
%
%
\begin{equation}
dF^{R}(\boldsymbol x, \theta) = \sum_{i=0}^{P} dF^{R}_i(\boldsymbol x) \varphi_i(\{ \xi_r(\theta)\})
\end{equation}
%
%
\begin{equation}
E(\boldsymbol{x}, \theta) = \sum_{i=0}^{P} E_i(\boldsymbol x) \varphi_i(\{ \xi_r(\theta)\})
\label{eq:E_expanded}
\end{equation}
%
%
Substituting Equations \ref{eq:Ha_expanded} $\sim$ \ref{eq:E_expanded}
into Equations \ref{eq:AF_constitutive} and
\ref{eq_simplified_tangential_stiffness} and applying Galerkin projection on
PC basis $\varphi_i\{ \xi_r(\theta)\}$ yields:
%
%
\begin{equation}
\sum_{m=0}^{P} dF_m^{R} \braket{\varphi_m \varphi_i} = \sum_{j=0}^{P} \sum_{k=0}^{P} H_{aj}d\eta_{k} \braket{\varphi_j \varphi_k \varphi_i} \pm \sum_{l=0}^{P} \sum_{n=0}^{P} \sum_{s=0}^{P} C_{rl} F^{R}_{n} d\eta_{s} \braket{\varphi_l \varphi_n \varphi_s \varphi_i}
\end{equation}
%
%
\begin{equation}
\sum_{i=0}^{P} E_m \braket{\varphi_m \varphi_i} = \sum_{j=0}^{P} H_{aj} \braket{\varphi_j \varphi_i} \pm \sum_{l=0}^{P} \sum_{n=0}^{P} C_{rl} F^{R}_{n} \braket{\varphi_l \varphi_n \varphi_i}
\end{equation}
%
%
% By using the orthogonality of Hermite PC bases $\braket{\varphi_i \varphi_j}=
% \delta_{ij}$, solutions to the unknown PC coefficients of incremental
% inter-story force $dF^{R}(\boldsymbol x,
% \theta)$ and inter-story stiffness $E(\boldsymbol{x}, \theta)$ can be
% written as:
%
By using the orthogonality of Hermite PC bases $\braket{\varphi_i \varphi_j}=
0$ for $i\neq j$, solutions to the unknown PC coefficients of incremental
inter-story force $dF^{R}(\boldsymbol x,
\theta)$ and inter-story stiffness $E(\boldsymbol{x}, \theta)$ can be
written as:
%
% \begin{equation}
% dF_i^{R} = \frac{1}{\braket{\varphi_i^2}} \bigg[ H_{aj}d\eta_{k} \braket{\varphi_j \varphi_k \varphi_i} \pm C_{rl} F^{R}_{n} d\eta_{s} \braket{\varphi_l \varphi_n \varphi_s \varphi_i} \bigg]
% \end{equation}
%
\begin{equation}
dF_i^{R} = \frac{1}{\rm Var[\varphi_i]} \bigg[ H_{aj}d\eta_{k} \braket{\varphi_j \varphi_k \varphi_i} \pm C_{rl} F^{R}_{n} d\eta_{s} \braket{\varphi_l \varphi_n \varphi_s \varphi_i} \bigg]
\end{equation}
%
% \begin{equation}
% E_i = H_{ai} \pm \frac{1}{\braket{\varphi_i^2}} C_{rl} F^{R}_{n} \braket{\varphi_l \varphi_n \varphi_i}
% \end{equation}
%
\begin{equation}
E_i = H_{ai} \pm \frac{1}{\rm Var[\varphi_i]} C_{rl} F^{R}_{n} \braket{\varphi_l \varphi_n \varphi_i}
\end{equation}
%
where $\braket{\cdot}$ is the expectation operator.
%
$\rm Var[\varphi_i]$ is the scalar variance of PC basis $\varphi_i\{
\xi_r(\theta)\}$, which equals to $\braket{\varphi_i^2}$.
%
It is noted that in the above equations, Einstein's tensor summation notation is used with index $i$ as a free index.
%
Terms $\braket{\varphi_j \varphi_k \varphi_i}$, $\braket{\varphi_l \varphi_n
\varphi_i}$ and $\braket{\varphi_l \varphi_n \varphi_s \varphi_i}$ are the
expectation of triple and quadruple product of PC basis $\varphi_i\{
\xi_r(\theta)\}$.
%
The above 1D formulation for SEPFEM is implemented in the context of
explicit, forward Euler algorithm,
%
The expanded stiffness matrix $K_{minj}$ is constructed using stiffness PC
coefficients $^{(n)}E_i$ at step $n$ following Equation
\ref{eq:stiffness_stochastic_matrix}.
%
Displacement PC coefficients $^{(n+1)}u_{nj}$ of step $n+1$ are then solved
by applying force vector $^{(n+1)}F_{mi}$ and using stiffness matrix $K_{minj}$
within Equation \ref{eq:deterministic_ODE}.
%
Following that, incremental inter-story drift PC coefficients $^{(n+1)}d\eta_i$
are calculated from displacement response $^{(n+1)}u_{nj}$ and incremental
uncertain restoring force $^{(n+1)}dF_i^{R}$ can be quantified as:
%
% %
% \begin{equation}
% {}^{(n+1)}dF_i^{R}
% =
% \frac{1}{\braket{\varphi_i^2}} \bigg[ H_{aj}\,^{(n+1)}d\eta_{k}
% \braket{\varphi_j \varphi_k \varphi_i} \pm C_{rl} \, {}^{(n)}F^{R}_{n}
% \,
% {}^{(n+1)}d\eta_{s} \braket{\varphi_l \varphi_n \varphi_s \varphi_i} \bigg]
% \end{equation}
% %
%
\begin{equation}
{}^{(n+1)}dF_i^{R}
=
\frac{1}{\rm Var[\varphi_i]} \bigg[ H_{aj}\,^{(n+1)}d\eta_{k}
\braket{\varphi_j \varphi_k \varphi_i} \pm C_{rl} \, {}^{(n)}F^{R}_{n}
\,
{}^{(n+1)}d\eta_{s} \braket{\varphi_l \varphi_n \varphi_s \varphi_i} \bigg]
\end{equation}
%
Updating the restoring force $^{(n+1)}F_i^{R}$ is then:
\begin{equation}
{}^{(n+1)}F_i^{R} = {}^{(n)}F_i^{R} + {}^{(n+1)}dF_i^{R}
\end{equation}
%
%
while new stiffness PC coefficients $^{(n+1)}E_i$ at step $n+1$ are then:
%
%
% \begin{equation}
% {}^{(n+1)}E_i
% =
% H_{ai} \pm \frac{1}{\braket{\varphi_i^2}} C_{rl}
% \,
% {}^{(n+1)}F^{R}_{n} \braket{\varphi_l \varphi_n \varphi_i}
% \end{equation}
%
\begin{equation}
{}^{(n+1)}E_i
=
H_{ai} \pm \frac{1}{\rm Var[\varphi_i]} C_{rl}
\,
{}^{(n+1)}F^{R}_{n} \braket{\varphi_l \varphi_n \varphi_i}
\end{equation}
%
%
%Other iterative algorithms, such as Backward Euler \cite{Jeremic96b}, could
%also be formulated similarly as in deterministic elastic-plastic FEM and are not
%repeated for conciseness.
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Illustrative Example}
\label{sec:illustrative_example}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
To illustrate the proposed framework, seismic risk of a typical eight story
shear frame structure that has been studied by many researchers \cite{Li2006b,Mitseas2018,Papazafeiropoulos2017,Xu2019} , is developed.
%
% Structure used in this study has been studied earlier by \citeN{Xu2019}.
%
The frame structure is shown in Figure \ref{figure_eight_story_building}.
%
%
\begin{figure}[!htbp]
\centering
\includegraphics[width=10.5cm]{Figures/shear_frame_illustration_update.pdf}
\caption{\label{figure_eight_story_building} Eight-story shear frame structure with uncertain floor stiffness under non-stationary seismic motions.}
\end{figure}
The hysteretic restoring force versus inter-story drift behavior is described by
Armstrong-Fredrick model presented in section \ref{sec:SFEM}.
%
Material parameter $H_a$ of Armstrong-Fredrick model is assumed to be
Gaussian distributed random field with $15\%$ coefficient of variation.
%
Means of material parameter $H_a$ are given for different floors as:
$H_{a1} \sim H_{a2}$ $1.59\times 10^7 N/m$,
$H_{a3} \sim H_{a6}$ $1.66\times 10^7 N/m$ and
$H_{a7} \sim H_{a8}$ $1.76\times 10^7 N/m$.
%
The correlation structure of parameter $H_a$ is assumed to be exponential
between different
floors, with correlation length of $l_c=10$ floors.
%
Material parameter $C_r$ is assumed to be $C_r = 17.6$~1/m.
%
The resultant mean hysteretic behavior of first floor is also shown in Figure \ref{figure_eight_story_building}.
%
Floor masses are assumed to be deterministic.
%
Rayleigh damping $\boldsymbol C=\alpha \boldsymbol M + \beta \boldsymbol K$ is
used with parameters $\alpha$ and $\beta$ chosen to be $\alpha=0.22$Hz and
$\beta=0.008s$.
%
Other structure modeling parameters are given in Table
\ref{tab_structure_parameter}.
%
Those parameters are determined from \citeN{Xu2019}.
%
Parameters $H_a$ and $C_r$ are calibrated to match the hysteretic behavior shown in \citeN{Xu2019}.
\begin{table}[!htbp]
\caption{\label{tab_structure_parameter}Parameters of the eight-story shear
frame structure.}
\centering
\small
\begin{tabular}{cccccc}
%\hline
$h_0[m]$ & $h[m]$ & $m_1$$\sim$$m_2[kg]$ & $m_3$$\sim$$m_4[kg]$ & $m_5$$\sim$$m_6[kg]$ & $m_7$$\sim$$m_8[kg]$ \\
\hline
$4$ & $3$ & 2$\times 10^{5}$ & 2.2$\times 10^{5}$ & 2.4 $\times 10^{5}$ & 2.5 $\times 10^{5}$ \\
%\hline
\end{tabular}
\end{table}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Seismic Source Characterization}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
The structure is
located at coordinate ($10km$, $40km$), $50km$ away from a strike slip fault
with $90^{\circ}$ dip angle, as shown in Figure
\ref{figure_single_fault_problem_setup}.
%
%
\begin{figure}[!htbp]
\centering
\includegraphics[width=0.60\textwidth]{Figures/single_fault_problem_setup.pdf}
\caption{\label{figure_single_fault_problem_setup} Seismic risk analysis of an
eight-story shear frame structure (red triangle) with uncertain stiffness $K$
subjected to earthquakes from a strike slip fault (black line).}
\end{figure}
%
The fault length is $250km$ with annual slip rate of $20mm/yr$.
%
Detailed geometry and model parameters for SSC of the strike slip fault are
given in Table \ref{table:parameter_single_fault}.
%
%
\begin{table}[!htbp]
\caption{\label{table:parameter_single_fault} Parameters for seismic source
characterization (SSC) of the strike slip fault.}
\centering
\footnotesize
\begin{tabular}{cc}
%\hline
Parameter & Value \\
\hline
Fault length & 250km \\
Fault width & 15km \\
Dip angle & $90^{\circ}$ \\
Slip rate $\mathcal{S}$ & 20mm/yr \\
Style of faulting & Strike slip \\
% Coordinate & (60,-100), (60, 150) \\
$f(M)$ & Truncated normal with $\sigma$=0.2 $n\sigma_{max}$=2 \cite{Hale2018}\\
$f(A|M)$ & Delta function at $log(A) = M-4$ \\
$f(W|A)$ & Delta function at $W=\sqrt{1.5A}$, limited to fault width \\
$f(Y)$ & Uniform distribution \\
$f(Z)$ & Uniform distribution \\
%\hline
\end{tabular}
\end{table}
%
Mean characteristic magnitude of the fault $\overline{M}$ is 7.6, and
is related to fault area $A$ \cite{Leonard2010} as:
%
\begin{equation}
\small
\overline{M} = log(A) + 4
\label{eq:mean_characteristic}
\end{equation}
%
%
Only earthquakes with magnitude greater than 5 (i.e. $M_{min}=5$) are
considered.
%
Following the procedure of SSC in section \ref{sec:SSC}, annual rate of
earthquakes occurring on the fault is $\overline{\lambda}=0.0067$/yr.
%
Probabilistic scenario space $\lambda(M, R, \theta)$ is discretized into four
mutually exclusive scenario events $S_i(M_i, R_i, \Theta_i)$ as shown in Table
\ref{table:seismic_scenarios}.
%
The computation is performed with probabilistic seismic hazard analysis program
HAZ45 \cite{Hale2018} using 0.2 for magnitude step $\Delta M$ and 2km for
distance step $\Delta R$.
%
\begin{table}[!htbp]
\caption{\label{table:seismic_scenarios} Seismic scenarios for the strike slip
fault.}
\small
\centering
\begin{tabular}{cccc}
\hline
\ \ Scenario ID & $M$ & $R_{rup}$ [km] & Annual rate $\lambda(M, R_{rup})$ \\ \hline
\ \ \ 1&\ \ \ \ 7.3\ \ \ \ & 56 & 9.54$\times 10^{-4}$ \\
\ \ \ 2&\ \ \ \ 7.5\ \ \ \ & 56 & 2.40$\times 10^{-3}$ \\
\ \ \ 3&\ \ \ \ 7.7\ \ \ \ & 56 & 2.40$\times 10^{-3}$ \\
\ \ \ 4&\ \ \ \ 7.9\ \ \ \ & 56 & 9.54$\times 10^{-4}$ \\ \hline
\end{tabular}
\end{table}
%
%
% Activity rate $N(M>M_{min})$ is computed as by balancing moment rate on the single fault.
%
% Combining the above activity rate $N(M>M_{min})$ with probabilistic models of rupture dimension and hypocenter location, full seismic source characterization can be performed with many available hazard analysis programs.
%
% A list of seismic scenarios with different magnitude $M$, distance $R$ and annual occurrence rate $\lambda(M, R)$ can be obtained.
%
% In this study, program HAZ45 \cite{Hale2018} was used for SSC of the strike slip fault.
%
% Developed are 4 different scenarios as shown in Table \ref{table:seismic_scenarios}.
%
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Time Domain Stochastic Ground Motion Modeling and Representation}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%
For each characterized seismic scenario $S_i(M_i, R_i, \Theta_i)$, $500$
realizations $\{\boldsymbol \Gamma_i\}$ of time domain uncertain motions are
simulated using methodology described in section
\ref{sec:stochastic_ground_motion_modeling}.
%
Figure \ref{figure_time_series_colletction} shows the first 200 realizations of simulated motions for earthquake scenario 1 with $M=7.3$, $R_{rup}=56km$.
\begin{figure}[!htbp]
\centering
\includegraphics[width=0.6\textwidth]{Figures/Acc_realization_200.pdf}
\caption{\label{figure_time_series_colletction} Realizations of uncertain seismic motions for scenario $M=7.3$, $R_{rup}=56km$.}
\end{figure}
In this study, ground motion populations from four different scenarios are
combined into a single population $\boldsymbol{\Gamma}$ using Equation
\ref{eq_combine_population} and modeled as a non-stationary random process.
%
The random process is represented by multi-dimensional Hermite polynomial chaos
(PC) following the technique formulated in section \ref{sec:pc_kl_expansion}.
%
Since marginal distribution of the random process is observed to be Gaussian
(section \ref{sec:stochastic_ground_motion_modeling}), theoretically, PC
representation with order 1 is sufficient.
%
The dimension of PC basis needs to be carefully chosen to reconstruct the
correlation structure of the original random process.
%
To ensure the accuracy of PC-KL representation, following error measurements are
defined and evaluated:
%
\begin{itemize}
\item
The absolute error on marginal mean of the random process:
\begin{equation}
% \varepsilon_m = \frac{1}{N_t} \sum_{k=1}^{N_t}|\frac{\mu(t_k)-\hat{\mu}(t_k)}{\sigma(t_k)}|
\varepsilon_m = \frac{1}{N_t} \sum_{k=1}^{N_t}|\mu(t_k)-\hat{\mu}(t_k)|
\end{equation}
%
\item
The absolute error on marginal standard deviation of the random process:
\begin{equation}
\varepsilon_{std} = \frac{1}{N_t} \sum_{k=1}^{N_t}|\sigma(t_k)-\hat{\sigma}(t_k)|
\end{equation}
%
\item
The absolute error on correlation of the random process:
\begin{equation}
\varepsilon_{corr} = \frac{1}{N_t^2} \sum_{k=1}^{N_t} \sum_{l=1}^{N_t}|Cov(t_k, t_l)-\hat{Cov}(t_k, t_l)|
\end{equation}
%
\end{itemize}
\noindent
where $\mu(t_k)$, $\sigma(t_k)$ and $Cov(t_k, t_l)$ are the marginal mean,
marginal standard deviation and correlation field of simulated ground motion
population $\boldsymbol \Gamma$.
%
Terms $\hat{\mu}(t_k)$, $\hat{\sigma}(t_k)$ and $\hat{Cov}(t_k, t_l)$ are statistics
calculated from PC representation of the random process from Equations
\ref{eq:PC_marginal_mean}, \ref{eq:PC_marginal_variance} and
\ref{eq:PC_correlation_structure}.
%
Term $t_k$ denotes the $k^{th}$ time instance and $N_t$ is the total number of
time instances.
%
Hermite PC bases of order 1, dimension 20, 70, 150 and 300 are examined for
PC-KL expansion of random process motions.
%
The errors for different PC bases are given in Table
\ref{tab_error_measurements}.
%
%
%
\begin{table}[!htbp]
\footnotesize
\caption{\label{table_motion_characterization_error} Error in probabilistic
characterization of non-stationary acceleration and displacement random process
using PC-KL expansion with different dimensions.}
\centering
\label{tab_error_measurements}
\begin{tabular}{ccccc}
%\hline
Dimension of PC & Dim. 20 & Dim. 70 & Dim. 150 & Dim. 300 \\ \hline
Displacement mean error $\varepsilon_m$ & 8.63$\times 10^{-9}$ & 8.63$\times 10^{-9}$ & 8.63$\times 10^{-9}$ & 8.63$\times 10^{-9}$ \\
Displacement S.D. error $\varepsilon_{std}$ & 1.28$\times 10^{-7}$ & 1.28$\times 10^{-7}$ & 1.28$\times 10^{-7}$ & 1.28$\times 10^{-7}$ \\
Displacement correlation error $\varepsilon_{corr}$ & 0.059 & 2.26$\times 10^{-4}$ & 8.27$\times 10^{-6}$ & 3.06$\times 10^{-7}$ \\ \hline
Acceleration mean error $\varepsilon_m$ & 9.84$\times 10^{-9}$ & 9.84$\times 10^{-9}$ & 9.84$\times 10^{-9}$ & 9.84$\times 10^{-9}$ \\
Acceleration S.D. error $\varepsilon_{std}$ & 1.23$\times 10^{-7}$ & 1.23$\times 10^{-7}$ & 1.23$\times 10^{-7}$ & 1.23$\times 10^{-7}$ \\
Acceleration correlation error $\varepsilon_{corr}$ & 0.185 & 0.091 & 0.053 & 0.028 \\
%\hline
\end{tabular}
\end{table}
%
%
It can be observed that in all the four cases marginal behavior of the random
process motions is well captured with very small magnitudes of errors $\varepsilon_m$ and
$\varepsilon_{std}$.
%
As shown in Figure \ref{figure_marginal_verification}, synthesized marginal mean
and marginal standard deviation from PC representation match very well with
statistics of simulated motions.
%
% Due to the transient nature of uncertain seismic motions, the marginal mean of the random process is close to 0.
%
\begin{figure}[!htbp]
\centering
\subfloat[Displacement]{
\hspace{-1.4cm}
\includegraphics[width=0.6\textwidth]{Figures/KL_dis_verification.pdf}}
\subfloat[Acceleration]{
\hspace{-0.8cm}
\includegraphics[width=0.6\textwidth]{Figures/KL_acc_verification.pdf}}
\caption{\label{figure_marginal_verification} Comparison between PC-synthesized
(black line) marginal mean and marginal standard deviation (SD) and statistics
of simulated ground motion realizations (red line).}
\end{figure}
As the dimension of PC increases, the relative error of correlation structure
decreases while the computational cost in stochastic FEM increases.
%
It is noted that PC dimension 70 is already adequate to capture the relatively
smooth random displacement correlation field.
%
However, acceleration correlation field synthesized from PC dimension 70 is
overestimated among many time steps.
%
PC dimension 150 and 300 approximate acceleration correlation structure much
better.
%
Eventually, considering both accuracy and efficiency, Hermite PC of order 1,
dimension 150 is used to spectrally discretize the random process seismic
motions in stochastic FEM analysis.
%
The comparison between the exact correlation structure and the PC synthesized
correlation structure is shown in Figure \ref{figure_motion_correlation}.
%
%
\begin{figure}[!htbp]
\centering
\vspace{-5mm}
\subfloat[Exact displacement correlation field]{
\includegraphics[width=0.50\textwidth]{Figures/KL_corr_exact_dis.pdf}}
%\subfloat[PC synthesized displacement correlation field]{
\subfloat[Synthesized displacement correlation field]{
\includegraphics[width=0.50\textwidth]{Figures/KL_corr_PC_dis.pdf}}
\vspace{-5mm}
\\
%\vspace{-5mm}
\subfloat[Exact acceleration correlation field]{
\includegraphics[width=0.50\textwidth]{Figures/KL_corr_exact_acc.pdf}}
%\subfloat[PC synthesized acceleration correlation field]{
\subfloat[Synthesized acceleration correlation field]{
\includegraphics[width=0.50\textwidth]{Figures/KL_corr_PC_acc.pdf}}
\vspace{-2mm}
\caption{\label{figure_motion_correlation} Verification of PC synthesized
acceleration and displacement correlation field with PC dimension 150.}
\end{figure}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\subsection{Stochastic Galerkin FEM Analysis and Seismic Risk}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
In order to perform stochastic Galerkin FEM analysis, it is also necessary to
characterize the randomness of stiffness of the structural system.
%
In order to do that, Hermite PCs of dimension 2, 4 and 6 is used for
capturing the exponential correlation structure of random field parameter
$H_a(\boldsymbol x, \theta)$.
%
%
%
%
It can be observed from Figure \ref{figure_KL_stiffness} that PC dimension 4 can
reasonably well reconstruct the correlation of $H_a(\boldsymbol x, \theta)$.
%
\begin{figure}[!htbp]
\centering
\includegraphics[width=0.5\textwidth]{Figures/KL_verification_stiffness.pdf}
\caption{\label{figure_KL_stiffness} Characterization of exponential correlation
(correlation length $l_c=10$ floors) of uncertain structural parameter
$H_a(\boldsymbol x, \theta)$ using PCs of different dimensions.}
\end{figure}
%
% As mentioned before, parameter $C_r$ is assumed to be deterministic.
%
% With PC characterized $H_a$, the probabilistic hysteretic behavior of restoring
% force versus inter-story drift can be intrusively modeled following the
% stochastic Galerkin technique formulated in section \ref{sec:SFEM}.
%
With PC characterized structural parameters, the probabilistic hysteretic behavior of restoring
force versus inter-story drift can be intrusively modeled following the
stochastic Galerkin technique formulated in section \ref{sec:SFEM}.
%
Figure \ref{figure_probabilisitc_constitutive_relation} shows the probabilistic
response of restoring force of the first floor under cyclic loading.
%
%
% \begin{figure}[!htbp]
% \centering
% \includegraphics[width=0.6\textwidth]{Figures/constitutive_relation_MC_verification.pdf}
% \caption{\label{figure_probabilisitc_constitutive_relation} Intrusive inelastic
% probabilistic modeling of Armstrong-Frederick hysteretic behavior and verification with Monte Carlo simulation: Gaussian
% distributed $Ha$ with mean 1.76 $\times 10^{7} \ N/m$ and 15\% coefficient of
% variation (COV), $C_r = 17.6$.}
% \end{figure}
%
\begin{figure}[!htbp]
\centering
\subfloat[Uncertain $H_a$]{
\hspace{-0.8cm}
\includegraphics[width=0.53\textwidth]{Figures/constitutive_relation_uncertainHa_certainCr_MC_verification.pdf}}
\subfloat[Uncertain $H_a$ and $C_r$]{
\hspace{-0.2cm}
\includegraphics[width=0.53\textwidth]{Figures/constitutive_relation_uncertainHa_uncertainCr_MC_verification.pdf}}
\vspace{-2mm}
\caption{\label{figure_probabilisitc_constitutive_relation} Intrusive probabilistic modeling of Armstrong-Frederick hysteretic behavior and verification with Monte Carlo simulation: (a) Gaussian distributed $Ha$ with mean 1.76 $\times 10^{7} \ N/m$ and 15\% coefficient of variation (COV), $C_r = 17.6$. (b) Gaussian distributed $Ha$ with mean 1.76 $\times 10^{7} \ N/m$ and 15\% coefficient of variation (COV), Gaussian distributed $C_r$ with mean 17.6 and 15\% COV.}
\end{figure}
%
%
Verification of developed constitutive modeling is performed using 10,000 Monte
Carlo simulations and shown in Figure \ref{figure_probabilisitc_constitutive_relation}
as well.
%
It can be seen that PC-based intrusive probabilistic hysteresis modeling
produces almost the same response as Monte Carlo simulations.
%
It is noted that intrusive probabilistic approach is approximately 2000 times
faster than corresponding Monte Carlo modeling.
%
% From the quantiles plot shown in Figure \ref{figure_probabilisitc_constitutive_relation}(a),
With both uncertain seismic motions (dimension 150) and uncertain structural
parameters (dimension 4) represented by Hermite PCs, probabilistic structural
displacement is described in 154 dimensional probabilistic space of Hermite PCs.
%
The unknown time varying PC coefficients, that contain all the information
about the probabilistic evolution of structural response, are intrusively solved
using developed Galerkin SEPFEM (section \ref{sec:SFEM}).
%
With these solved PC coefficients, a polynomial chaos based surrogate model
is analytically established~\cite{Sudret2008}.
%
After that, any probabilistic structural dynamic response can be
easily reconstructed.
%
Time evolving mean, standard deviation (SD) and correlation field
of any resulting field of interest
can be directly evaluated through Equations \ref{eq:PC_marginal_mean},
\ref{eq:PC_marginal_variance} and \ref{eq:PC_correlation_structure}.
%
By efficiently sampling the PC surrogate model, marginal or joint PDF of any
structural response of interest can also be obtained through kernel density
estimation.
For example, Figure \ref{figure_time_evolving_statistics} shows the time
evolving mean and standard deviation (SD) of the first and top floor deformation
relative to the ground.
%
%
\begin{figure}[!htbp]
\centering
\subfloat[Displacement mean]{
\includegraphics[width=0.55\textwidth]{Figures/mean_disp_longer.pdf}}
\subfloat[Displacement SD]{
\includegraphics[width=0.55\textwidth]{Figures/std_disp_longer.pdf}}
\caption{\label{figure_time_evolving_statistics} Time evolving mean and standard
deviation (SD) of the first and top floor deformation relative to the ground.}
\end{figure}
%
%
Due to inelastic, elastic-plastic response, uncertain permanent deformation is observed
in both mean and standard deviation of floor deformation.
%
It is noted that the deformation of top floor presents much larger variability
than that of the first floor.
Two typical engineering demand parameters (EDPs) are selected for seismic risk
analysis: Maximum inter-story drift ratio (MIDR) and Peak floor acceleration
(PFA) \cite{Miranda2005,Miranda2006}.
%
We define MIDR as a function of probabilistic dynamic floor displacement:
%
\begin{equation}
MIDR_i(\theta) = \max_{t\in [0, T]} \Big\{\frac {|u_i(t, \theta) - u_{i-1}(t, \theta)|}{H_i} \Big\}
\label{eq:MIDR_floor}
\end{equation}
%
\begin{equation}
MIDR(\theta) = \max_{i\in [1, 8]}\max_{t\in [0, T]} \Big\{ \frac {|u_i(t, \theta) - u_{i-1}(t, \theta)|}{H_i} \Big\}
\end{equation}
%
%
%
where $MIDR_i(\theta)$ and $u_i(t, \theta)$ are the probabilistic MIDR and
displacement of the $i^{th}$ floor, respectively, and $H_i$ is the floor height,
while probabilistic MIDR of the whole shear frame structure is given as
$MIDR(\theta)$.
%
Probabilistic floor accelerations are defined as :
%
%
%PFA is defined as a function of probabilistic dynamic floor acceleration:
%
\begin{equation}
PFA_i(\theta) = \max_{t \in [0, T]} \{|\ddot{u}_i(t, \theta)|\}
\end{equation}
%
\begin{equation}
PFA(\theta) = \max_{i\in [1, 8]} \max_{t \in [0, T]} \{|\ddot{u}_i(t, \theta)|\}
\label{eq:PFA_overall}
\end{equation}
%
%
where $PFA_i(\theta)$ and $\ddot{u}_i(t, \theta)$ are the probabilistic PFA and
acceleration of the $i^{th}$ floor, respectively, while $PFA(\theta)$ is the
probabilistic PFA of the whole structure.
%
Since both probabilistic displacements $u_i(t, \theta)$ and probabilistic
accelerations $\ddot{u}_i(t, \theta)$ are well defined through resulting PC
coefficients, probabilistic response of MIDR and PFA are readily
available through Equations \ref{eq:MIDR_floor} to \ref{eq:PFA_overall}.
%
For example problem, the probability density evolution of $MIDR(\theta)$
is shown in Figure \ref{figure_probability_evolution}.
%
\begin{figure}[!htbp]
\subfloat[PDF surface]{
\hspace{-1cm}
\includegraphics[width=0.65\textwidth]{Figures/MIDR_PDF_evolution.pdf}}
\subfloat[PDF contour]{
\includegraphics[width=0.50\textwidth]{Figures/Contour_PDF_evolution_MIDR.pdf}}
\caption{\label{figure_probability_evolution} Time evolving probability density
function (PDF) of MIDR for frame structure.}
\end{figure}
%
%
At $t=0s$, the structure is deterministically at rest, therefore, the PDF of
$MIDR$ tends to infinity, i.e., a delta function centered at zero and as such is
not shown in Figure~\ref{figure_probability_evolution}.
%
Figure \ref{figure_evolving_PDF_combine} shows typical PDFs at three different
times.
%
\begin{figure}[!htbp]
\centering
\includegraphics[width=0.6\textwidth]{Figures/PDF_MIDR_combine.pdf}
\caption{\label{figure_evolving_PDF_combine} PDF of MIDR at different times:
$t=15$s, $21$s and $29$s.}
\end{figure}
%
It can be observed that PDF of MIDR is dispersing during first half of the
seismic loading, while toward the end of the loading, it shows high kurtosis,
due to reduced variation in input excitations.
%
The PDFs of MIDR of several different floors (1st, 3rd, 6th and top
floor) and the whole frame structure are shown in Figure
\ref{figure_MIDR_different_stories}(a).
%
%
It is observed that the mean of MIDR increases along with larger dispersion,
from the top to the bottom floor.
%
This is expected considering the increase of shear force from the top floor to
the base.
%
The MIDR PDF of the first floor almost overlaps with that of the whole
structure, which indicates that the maximum inter-story drift happens at first
floor.
%
From the probabilistic distribution of MIDR, exceeding probability
$P(EDP>z|\boldsymbol \Gamma)$ can be obtained.
%
Combining exceeding probability and scenario rate, EDP hazard of MIDR is
calculated using Equation \ref{eq:time_domain_risk_combine} and is shown in
Figure \ref{figure_MIDR_different_stories}(b).
%
%
\begin{figure}[!htbp]
\centering
\subfloat[PDF MIDR]{
\hspace*{-0.7cm}
\includegraphics[width=0.53\textwidth]{Figures/MIDR_distribution_different_floors.pdf}}
\subfloat[EDP hazard of MIDR]{
\hspace*{-0.7cm}
\includegraphics[width=0.53\textwidth]{Figures/Risk_MIDR.pdf}}
\caption{\label{figure_MIDR_different_stories} PDF and annual exceedance rate of
MIDR between different story over the whole loading history.}
\end{figure}
%
It can be seen that the demand of MIDR is dominantly controlled by lower floors,
e.g., the 1st and 3rd floor.
%
In addition to PDFs of MIDR, PDFs of PFA for different floors and the whole
frame structure are developed and shown in Figure
\ref{figure_PFA_different_stories}(a).
%
The distributions of PFA of the 1st, 3rd and 6th floor are close to each other,
while the PFA of the top floor shows larger mean and variability.
%
The PFA distribution of the top floor is very close to that of the whole
structure, which indicates the top floor tends to experience the maximum
acceleration.
%
EDP hazard of PFA is shown in Figure \ref{figure_PFA_different_stories} (b).
%
The demand of PFA is dominantly controlled by the top floor.
%
\begin{figure}[!htbp]
\centering
\subfloat[PDF PFA]{
\hspace*{-0.7cm}
\includegraphics[width=0.53\textwidth]{Figures/PFA_distribution.pdf}}
\subfloat[EDP hazard of PFA]{
\hspace*{-0.7cm}
\includegraphics[width=0.53\textwidth]{Figures/Risk_PFA.pdf}}
\caption{\label{figure_PFA_different_stories} PDF and annual exceedance rate of
PFA of different stories and the whole frame structure.}
\end{figure}
By assuming that damage measure (DM) is a step function of EDP, seismic risk for damage
states using different levels of MIDR and PFA exceedance can be directly
determined from the EDP hazard curve.
%
As shown in Table \ref{tab_seismic_risk_summary}, seismic risk for MIDR$>1\%$ is
$3.83\times 10^{-3}$ and the risk for PFA$>1m/s^2$ is $1.92\times10^{-3}$.
%
\begin{table}[!htbp]
% \small
\footnotesize
% \scriptsize
\caption{\label{tab_seismic_risk_summary}Seismic risk of damage state for
different levels of MIDR and PFA exceedance.}
\resizebox{0.98\hsize}{!}{
\begin{tabular}{cccccc}
%\hline
MIDR\textgreater{}0.5\% & MIDR\textgreater{}1\% & MIDR\textgreater{}2\% & PFA\textgreater{}0.5${\rm m/s^2}$ & PFA\textgreater{}1${\rm m/s^2}$ & PFA\textgreater{}1.5${\rm m/s^2}$ \\
\hline
6.66$\times 10^{-3}$ & 3.83$\times 10^{-3}$ & 9.97$\times 10^{-5}$ & 6.65$\times 10^{-3}$ & 1.92$\times 10^{-3}$ & 9.45$\times 10^{-5}$ \\
%\hline
\end{tabular}}
\end{table}
As noted earlier, complete probabilistic structural response, including both marginal
distribution and correlation information, is contained in PC coefficients, any
other EDP or other DM defined on multiple EDPs can also be developed with little
additional effort.
%
Figure \ref{figure_joint_PFA_MIDR} shows the 2D joint PDF,
$f(\text{MIDR}, \text{PFA} | \, \boldsymbol \Gamma)$ of two EDPs,
MIDR and PFA, evaluated from
the PC-based surrogate model of probabilistic structure response.
%
%
\begin{figure}[!htbp]
\centering
\subfloat[Joint PDF of MIDR and PFA]{
\hspace*{-0.7cm}
\includegraphics[width=0.57\textwidth]{Figures/2D_EDP_PDF_1e5.pdf}}
\subfloat[Contour of 2D joint PDF]{
\hspace*{-0.7cm}
\includegraphics[width=0.53\textwidth]{Figures/2D_EDP_PDF_downview_1e5.pdf}}
\caption{\label{figure_joint_PFA_MIDR} 2D joint PDF of MIDR and PFA of the whole
shear frame structure.}
\end{figure}
%
It can be observed that in this case MIDR and PFA are positively correlated.
%
The correlation coefficient is $0.64$.
%
For damage measure (DM) defined on multiple EDPs, for example, $DM:
\{\text{MIDR}>z_1\, \vee\, \text{PFA}>z_2\}$, seismic risk can be evaluated as:
\begin{equation}
\lambda(\text{MIDR}>z_1 \vee \text{PFA}>z_2)
=
\overline{\lambda}\int_{\mathscr{D} } f(\text{MIDR}, \text{PFA} |\, \boldsymbol \Gamma) \, d \mathscr{D} % \bigcup
\end{equation}
%
%
where $\overline{\lambda}$ is the annual occurrence rate of seismic scenario
that would induce ground motion population $\boldsymbol \Gamma$, while
$\mathscr{D}$ is the integral domain $(\text{MIDR}, \text{PFA}) \in $ $[z_1,
+\infty] \bigcup [z_2, +\infty]$ according to the definition of damage measure.
%
%
% Using such approach, seismic risk for damage state $DM:\{\text{MIDR}>1\%\,
% \vee\,\text{PFA}>1{\rm m/s^2}\}$ can be calculated as $4.20 \times 10^{-3}$, while the
% risk for damage state $DM:\{\text{MIDR}>1\%\, \wedge \,\text{PFA}>1{\rm m/s^2}\}$ is
% 60\% less, equal to $1.71 \times 10^{-3}$.
%
%
% {\huge \bf NOTE:
% need to say in words what is
% $DM:\{\text{MIDR}>1\%\, \vee\,\text{PFA}>1{\rm m/s^2}\}$ is, and the other one too...
%
% so for example:
%
% . . .
% Damage state DM defined for MIDR that is greater than 1\% and PFA that is greater than $1m/s^2$,
% ($DM:\{\text{MIDR}>1\%\, \vee\,\text{PFA}>1{\rm m/s^2}\}$)
% . . .
%
% In a sense all mathematical equations and symbols have to be described in words if in text...
%
% }
%%%% ===== Hexiang Reply: Has changed accordingly ==========
%
Using such approach, seismic risk for damage state DM defined for either
MIDR greater than 1\% or PFA greater than $1{\rm m/s^2}$ (i.e.,
$DM:\{\text{MIDR}>1\%\, \vee\,\text{PFA}>1{\rm m/s^2}\}$), can be calculated as
$4.20 \times 10^{-3}$, while the risk for damage state defined for both MIDR
greater than 1\% and PFA greater than $1{\rm m/s^2}$ (i.e., $DM:\{\text{MIDR}>1\%\,
\wedge \,\text{PFA}>1{\rm m/s^2}\}$) is 60\% less, equal to $1.71 \times
10^{-3}$.
%
Both of these risk values based on joint EDPs are rather different from the
ones calculated using single EDP.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
\section{Conclusions}
\label{sec:conclusions}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
A time domain intrusive probabilistic seismic risk analysis framework
for performance based earthquake engineering was described in some detail.
%
Methodology to simulate non-stationary stochastic seismic motions was
presented.
%
The presented methodology is directly compatible with state-of-the-art seismic
source characterization.
%
Different source, path and site factors are explicitly accounted for in the
stochastic modeling of Fourier amplitude spectrum and Fourier phase derivative.
%
Both uncertain seismic motions and uncertain structural parameters are
characterized as random process/field and represented with Hermite polynomial
chaos (PC) Karhunen-Lo{\`e}ve (KL) expansion.
%
Direct polynomial chaos based Galerkin intrusive modeling of 1D elastic-plastic response
was formulated and applied to simulate the uncertain hysteretic
behavior of restoring force versus inter-story drift for shear frame structure.
%
Formulations for random stiffness polynomial chaos coefficients were derived and incorporated
into stochastic Galerkin elastic-plastic finite element method.
%
% The issue of PC characterization of random stiffness in current probabilistic
% elastic-plasticity based on Fokker-Planck-Kolmogorov (FPK) equation is resolved.
% %
Using developed stochastic elastic-plastic finite element method, probabilistic dynamic response
of uncertain structural system driven by uncertain motions is intrusively
solved.
%
Following that, seismic risk for damage measure defined on single or multiple
engineering demand parameter(s) was calculated.
%
The proposed framework is illustrated within seismic risk analysis of an eight-story
shear frame structure excited by uncertain strike-slip fault earthquakes.
%
Presented new framework avoids the drawbacks of choosing and using intensity measure(s).
%
All the seismic motion characteristics and their uncertainties, for example,
uncertain peak ground acceleration (PGA), spectrum acceleration (Sa) and others, are captured by random process motions and
directly propagated into uncertain structural system.
%
Development of ground motion prediction equations (GMPEs) for potentially new intensity measures (IMs) (e.g., Arias intensity or cumulative absolute velocity) and repetitive Monte Carlo fragility simulations are
circumvented.
%
Though most of current seismic risk analyses are performed for damage measure
defined on single engineering demand parameter, presented framework can also handle joint engineering demand parameters/failure
criteria without much additional effort.
%
It is found that, for different damage measure defined on joint engineering demand parameters,
corresponding seismic risk significantly varies and is rather different from the
risk value for single engineering demand parameter.
%
Therefore, considering damage measure based on joint engineering demand parameters can be of great
interest in seismic risk analysis.
%
Future work will focus on accuracy and efficiency comparison between the
proposed framework and existing intensity measure based, non-intrusive seismic risk analysis and also applying the proposed framework to more realistic engineering structures.
%
% Here's the list of references:
%
% \label{section:references}
\bibliography{refmech}
%
\end{document}