# Spectral Representation-Based Multidimensional Nonstationary Ground Motion Model for Seismic Reliability Analysis of Frame Structures

Authors: Lawali Moussa Laminou; Xinghua Chen
Journal: Shock and Vibration
 (2021)

---

## Abstract

A framework for a multidimensional nonstationary ground motion model based on spectral representation theory is proposed in this paper. The multidimensional nonstationary ground motion model is built from a local target to fit the multidimensional response spectrum. A four-stage modulation function takes into account the multidimensional intensity correlation and the modified Clough–Penzien (C-P) power spectrum with parameter correlation, which represent the two main aspects, the modulation function and the power spectrum of constructing the multidimensional nonstationary ground motion model. A multidimensional response spectrum constructed according to the standardizing response spectrum is used as the fitting target response spectrum. Samples of random ground motion for random seismic response and dynamic reliability study are finally obtained. The random seismic responses are then combined with the probability density evolution method (PDEM) to carry out the seismic reliability analysis of a randomly base-excited moment-resisting frame structure. In the numerical analysis, the nonlinear seismic responses and reliability of a 10-story reinforced concrete frame structure are carefully investigated in accordance with the Egyptian seismic code. As a result, the effectiveness of the proposed method is fully demonstrated.


---

## Body

## 1. Introduction

After critical earthquake events, structures in the affected zones suffer a considerable amount of damage in most cases. The economic consequences of seismic events can also be catastrophic in the adjacent zones. Although Egypt is affected by only moderate seismic activity compared to other countries, it is exposed to high seismic risk due to the concentration of population as well as important archaeological sites in its seismic zones, together with the lack of appropriate design and construction methods. A damaging earthquake represents a real threat to the safety, social integrity, and economic wellbeing of the region. Thus, seismic studies in Egypt are urgently needed to identify degrees of vulnerability in seismic areas and inform future risk studies, as well as construction codes and land-use planning improvement.Protecting human lives is the leading issue in seismic design [1]. Modern building codes recommend incorporating the seismic response time history under certain conditions in the seismic analysis and design of structures. Reliable ground motion time series are required for accurate nonlinear dynamic analysis of structures [2]. In most cases, there will not be strong motion records for a given site. Even if such recordings are available, there is no basis to expect that a future earthquake will generate the same or similar ground motion. Due to the lack of reliable earthquake records, synthetic earthquake ground motions are often generated instead. Waezi and Rofooei introduce a new time-varying model for generating synthetic nonstationary acceleration records by using a modified Kanai-Tajimi model with time-variant parameters [3]. In other approaches, Li and Wang generate spectrum-compatible multidimensional artificial ground motion via wavelet transformation [4], Zhang et al. present a novel method based on surface fitting to derive sets of rotational components from translational-recorded ground motion [5], and Chen et al. [6] propose a simulation method for fully nonstationary spatially variable ground motion based on the Kameda time-varying power spectrum model. Finally, another recent paper presents a new model of the generalized evolution spectrum of fully nonstationary ground motion. This model is based on the power spectral density function of stationary ground motion [7].The occurrence of earthquakes in time and space is absolutely random in nature, and this leads to a considerable amount of uncertainty in the resulting ground motions [8]. The current basic seismic design codes depend on a uniform hazard spectrum, but this is applied to structures exposed to unpredictable levels of risk [9]. In view of this, probability and statistics concepts and methods are necessary to quantify and evaluate such uncertainty and its effects on engineering systems. In order to examine seismic risk and mitigate potential damage to a structure, it is important to accurately quantify the seismic reliability of a structure [10]. Camacho et al. have developed an alternative reliability-based methodology for structures under seismic loading, with the proposed approach producing adequate results for structural reliability evaluation [11]. Lie and Zio [12] report significant results for their proposed method using a combination of parallel Monte Carlo simulation and the recursive Bayesian method to develop dynamic reliability assessment and failure prognostics with noisy monitored data. In their work, Do et al. [13] present reliability assessment of structures with uncertain-but-bounded parameters under stochastic process excitation. Wang et al. [14] propose a novel reliability-based sensitivity analysis of dynamic random systems with time-dependent random parameters for application to linear structures. In another study, Shittu et al. show that the randomness in design parameters for complex structural systems subject to high uncertainties should be systematically accounted for via stochastic modeling [15].Engineering structures are simultaneously excited by multiple (horizontal and vertical) ground motion components during an earthquake [16]. A growing focus on improving multidimensional seismic analysis is evident in new studies. Among these, Qu et al. develop pushover methods for seismic assessment of buildings under multidimensional earthquakes, with the objective of overcoming the drawbacks of the conventional modal pushover method [17]. In other examples, Farahmand-Tabar and Barghian carry out a seismic assessment of a cable-stayed arch bridge under three-component orthotropic earthquake excitation [18], Lin et al. investigate the seismic behavior of long-span connected structures under multisupported and multidimensional earthquake excitations [19], and Guo et al. study the combined effect of vertical and horizontal ground motions on the failure probability of reinforced concrete chimneys [20]. Wang et al. present a detailed review of the development and evaluation of multidimensional combination rules and their adoption in modern design codes and standards [21]. However, this approach is still not being widely considered and implemented in codes as there remains a lack of necessary knowledge from careful studies to develop a simple and rational approach.Motivated by the necessity for multidimensional analysis and the high level of uncertainty involved in modeling ground motions, the key objective and contribution of this paper is to propose a framework for a multidimensional nonstationary ground motion model based on spectral representation theory and compatible with the Egyptian seismic code. The paper is structured as follows. First, the proposed framework is presented, followed by a detailed discussion about mathematical modeling of the effects of nonstationary multidimensional random ground motion composed of multiple sequences on inelastic structures. Then, the multidimensional representative time history used in the seismic design of building structures is generated according to the Egyptian code, and the efficacy of the proposed method is verified. Next, a brief overview of the PDEM-based dynamic reliability assessment adopted for the modeling is provided. This is followed by a detailed description of the case study, the probability density evolution analysis conducted on a 10-story reinforced concrete frame structure. The results of the analysis are thoroughly discussed, and the paper ends with the conclusions drawn from the analysis based on the key findings.


## 2. Proposed Multidimensional Nonstationary Ground Motion Model

When performing analysis of the stochastic earthquake response of nonlinear structures, it is necessary to have a great number of nonstationary ground motion processes with statistical features identical to the target spectrum. This paper proposes a multidimensional nonstationary ground motion model using a local perspective simplified method to fit the target response spectrum based on the requirements of engineering applications. The proposed model starts from the local simulation of the ground motion process (modulation function, power spectrum, and response spectrum); its control target is to fit the multidimensional response spectrum. The details of the multidimensional nonstationary ground motion model are discussed in the following sections.


### 2.1. Spectral Representation Based on Random Function Method

There are many methods available for generating nonstationary ground motion. In this study, the acceleration time history is expressed as(1)Xt=∑k=1NAt,kXkcoswkt+Yksinwkt,(2)At,k=2SXt,wkΔw,where wk=kΔw, Xk,Yk presents a standard set of orthogonal random variables, and the conditions that it must be met are(3a)EXk=EYk=0,EXjYk=0,(3b)EXjXk=EYjYk=δjk,in which E⋅ denotes the mathematical expectation and δmn is the Kronecker delta.Thus, the accuracy index for controlling the simulation method can be expressed as(4)εN=1−∫0wμ∫0TSxt,wdtdw∫0∞∫0TSxt,wdtdw,where εN represents the mean square relative error; wμ=NΔw indicates the cutoff frequency; and T is the duration of the nonstationary process.The expression of orthogonal random variablesX¯n and Y¯nn=1,2,…,N by a basic random variable Θ is defined as(5a)X¯n=2cosnΘ+α,(5b)Y¯n=2sinnΘ+α,in which Θ is uniformly distributed over the range 0,2π and α is an arbitrary constant, usually taken as α=π/4. Next, after seeking X¯n,Y¯nn=1,2,…,N, the functions rand (“state”) and randperm (N) provided by the Matlab software program are used to map X¯n,Y¯n to Xk,Ykk=1,2,…,N one by one, which is X¯n⟶Xk and Y¯n⟶Yk. And finally, the required Xk,Ykk=1,2,…,N can be uniquely determined.


### 2.2. Evolution Spectrum Model of Nonstationary Ground Motion

According to the construction approach of the nonstationary stochastic process in evolutionary power spectrum theory, the evolution power spectrum of the process of ground motion accelerationXt with nonstationary intensity can be expressed as(6)Sxt,w=ft2⋅Sw,where ft denotes the intensity envelope function and Sw is the stationary process power spectral density. Both will be described in detail in the following sections.


#### 2.2.1. Intensity Modulation Function

This paper sets out to study multidimensional nonstationary seismic response of structures. The commonly used intensity nonstationary model is selected for the single-dimensional component of the model, and the intensity envelope function is used to describe the strength nonstationary characteristics of the model. Furthermore, a four-stage continuous intensity envelope function model, proposed by Li and Lai [22], is considered to more realistically simulate the correlation of three-dimensional ground motion. In contrast to the popular Amin-Ang model [23], the model is obtained by statistical regression of a large number of measured seismic waves with two-level ground motion components using the least-squares method, thus better reflecting the rising phase, stationary phase, and attenuation phase of the ground motion time history. The intensity envelope function ft is expressed as(7)ft=0,0<t≤t0,t−t0t1−t02,t0<t≤t1,1,t1<t≤t2,exp−ct−t2,t>t2,in which c adjusts the decay speed of the falling portion of the model curve, t0 is the start time of the rising segment, t1 and t2 are the first and last moments of the peak plateau, respectively, T1=t1−t0 is the duration of the rising segment, and TS=t2−t1 is the duration of the plateau. The values of parameters t0,T1,TS, and c used in this paper are from Li and Lai [22].Their research on the multidimensional coherence of the components of ground motion shows that the correlation between the seismic components in the frequency characteristics is relatively low and the attenuation is not obvious [22]. Therefore, in order to simplify the complexity of simulating the correlation of each component, it is feasible to only consider the intensity characteristics in the seismic analysis of structures, without considering the effect of three-dimensional frequency correlation. For this reason, the three-dimensional correlation of the frequency characteristics is not considered in this paper.


#### 2.2.2. Power Spectrum

For the power spectrum models of ground motion processes, two main categories are considered in this paper: (1) intensity-frequency power spectrum model; (2) power spectra converted from standard design response spectrum.(1)Modified C-P Model. The Clough–Penzien (C-P) model [24] is an improvement on the classic Kanai–Tajimi model [25]. It can correct the energy of the low-frequency part and is widely used. Accordingly, in order to consider the frequency characteristics in the power spectrum model parameters, a modified C-P model is proposed, which can be expressed as(8)St,w=wg4t+4ζg2twg2tw2w2−wg2t2+4ζg2twg2tw2⋅w4w2−wf2t2+4ζf2twf2tw2⋅S0t,the modified C-P model is able to reflect the nonstationary intensity of the ground motion, while the parameter duration and frequency spectrum reflect the nonstationary characteristics of frequency and can be expressed as(9)wgt=w0−atT,ζgt=ζ0+btT,(10)wft=0.1wgt,ζft=ζgt,where w0, ζ0, a, and b are determined according to the site category and T is the ground motion duration. Similar to the C-P model of the steady ground motion acceleration process, wgt and ζgt can be called the dominant circular frequency and damping ratio of the site soil, respectively, and wft and ζft are the corresponding filter parameters. As we can see from equations (9) and (10), the site parameters and filtering parameters are both linear functions of time, which indicates that wgt and ζgt of the site soil are time-varying parameters that change within a certain range.In equation (8), the spectral parameter S0t reflecting the intensity of ground motion can be expressed as(11)S0t=a¯max2r¯2πwgt2ζgt+1/2ζgt,in which a¯max is the mean value of the peak ground acceleration or the direct ground motion acceleration, according to the specification; r¯ is the equivalent peak value factor.(2)Power Spectrum Converted from Standard Design Response Spectrum. Alternatively, when generating the ground motion time history by standard design response spectrum, an approximate conversion relationship between the acceleration response spectrum and power spectrum is required, defined as(12)Sw=−ζπwSa2w,ζ·1ln−π/wTln1−p,where Sw denotes the bilateral power spectral density function; Saw,ζ represents the absolute acceleration response spectrum; T is the ground motion duration; and p is the transcendence probability of the response spectrum, which can be taken as 0.5.In the process of conversion, the discrete control points are selected first on the design response spectrum, and the logarithmic coordinates are then selected uniformly over the period, expressed as(13)lgT0,k=lgT0,1+1−lgT0,1M−1k−1,k=1,2,…,M,where T0,k is the kth control point of the period and M is the total number of discrete control points in a 10-second period. The total number of discrete control points M is determined by equation (13), and the power spectrum value of the corresponding control point is determined by equation (12). The other power spectrum values are determined by linear interpolation, expressed as(14)Swi=Swk+wi−wkwk−1−wkSwk+1−Swk,i=1,2,…,N,k=1,2,…,M,where Swk and Swk+1 represent the power spectrum values of the kth and k+1 th points and Swi denotes the power spectrum value at the ith point, wi=iΔw.In order to make the conversion power spectrum more accurately fit the code’s design response spectrum, the fitting error is defined as the error between the design response spectrum and the average response spectrum [26] and is expressed as(15)εwk=Sawk−SkwkSawk×100%,where Sawk and Skwk represent the design response spectrum value and the average response spectrum value for the kth point. If the fitting error εwg<10%, the accuracy meets requirements; otherwise, the power spectrum requires iteration and the iteration formula [27] is expressed as(16)Si+1wk+1=SiwkSa2wkSk2wk.Due to the random characteristics of ground motion, the design response spectrum is obtained by statistical averaging and smooth processing of the response spectrum of a large number of actual earthquakes. Therefore, the converted power spectrum has practical engineering application value, and the average meaning can be combined with the randomness method.The converted power spectrum obtained by the response spectrum conversion can be used for the simulation of multidimensional ground motion. Taking the seismic zones 5B site and the response spectrum type I as an example, the duration isT = 30 s, and the peak acceleration value is 0.3 g. As can be seen in Figure 1, for the converted power spectrum obtained from the design response spectrum, horizontal x and horizontal y partially overlap in the low-frequency part, and the changing trend of the three-dimensional conversion power spectrum is basically the same on the whole. However, the peak and fall rates of each component are different, the peak is close to the low-frequency part, and the variation range is 10∼50 Hz. Since the response spectra are all response spectra in an average sense, the conversion power spectrum has an average meaning.Figure 1
Converted power spectrum.


### 2.3. Multidimensional Nonstationary Ground Motion Model

The main aspects of establishing a multidimensional ground motion model include considering the correlation of each seismic component from the whole; extending the one-dimensional power spectrum to a multidimensional cross-power spectrum matrix; and studying the correlation between various seismic component parameters (such as response spectrum, amplitude, or duration) locally.The key to establishing the multidimensional ground motion model is to determine the cross-spectral density between the seismic components in the cross-power spectrum matrix. However, the research on the coherence function between the seismic components of a single point is still insufficient, and cross-spectral density is difficult to accurately express. Accordingly, this paper proposes a simplified method to build a multidimensional nonstationary ground motion model from a local perspective. Specifically, based on the local simulation of the ground motion process (modulation function, power spectrum, and response spectrum), the three-dimensional response spectrum of the fitted structure is generally used as the control target, and the correlation between the intensity characteristics of the multidimensional seismic components is considered on the basis of the one-dimensional ground motion model.On one hand, the modified four-segment intensity envelope function is applied to consider the multisegment correlation. On the other hand, from the power spectrum perspective, the power spectrum can use the one converted by the three-dimensional response spectrum, and the same conversion method can be used for both horizontal and vertical components. For the bidirectional horizontal component, the power spectrum is the same value and is completely related. For the vertical component, the same formula is expressed as the bidirectional horizontal component. The parameters are specified to be different. The model uses the modified C-P power spectrum; then(17)Sxt,w=Syt,w,(18)Szt,w=wgV4t+4ζgV2twgV2tw2w2−wgV2t2+4ζgV2twgV2tw2⋅w4w2−wfV2t2+4ζfV2twfV2tw2⋅S0t,(19)wgVt=EGwgHt=1.58wgHt,(20)ζgVt=ζgHt,wfVt=0.1wgHt.At the same time, assuming that the vertical component acceleration peak is 2/3 of the horizontal component [28], then the relationship between the vertical and horizontal component spectral intensity factors is(21)S0Vt=0.281S0Ht.


## 2.1. Spectral Representation Based on Random Function Method

There are many methods available for generating nonstationary ground motion. In this study, the acceleration time history is expressed as(1)Xt=∑k=1NAt,kXkcoswkt+Yksinwkt,(2)At,k=2SXt,wkΔw,where wk=kΔw, Xk,Yk presents a standard set of orthogonal random variables, and the conditions that it must be met are(3a)EXk=EYk=0,EXjYk=0,(3b)EXjXk=EYjYk=δjk,in which E⋅ denotes the mathematical expectation and δmn is the Kronecker delta.Thus, the accuracy index for controlling the simulation method can be expressed as(4)εN=1−∫0wμ∫0TSxt,wdtdw∫0∞∫0TSxt,wdtdw,where εN represents the mean square relative error; wμ=NΔw indicates the cutoff frequency; and T is the duration of the nonstationary process.The expression of orthogonal random variablesX¯n and Y¯nn=1,2,…,N by a basic random variable Θ is defined as(5a)X¯n=2cosnΘ+α,(5b)Y¯n=2sinnΘ+α,in which Θ is uniformly distributed over the range 0,2π and α is an arbitrary constant, usually taken as α=π/4. Next, after seeking X¯n,Y¯nn=1,2,…,N, the functions rand (“state”) and randperm (N) provided by the Matlab software program are used to map X¯n,Y¯n to Xk,Ykk=1,2,…,N one by one, which is X¯n⟶Xk and Y¯n⟶Yk. And finally, the required Xk,Ykk=1,2,…,N can be uniquely determined.


## 2.2. Evolution Spectrum Model of Nonstationary Ground Motion

According to the construction approach of the nonstationary stochastic process in evolutionary power spectrum theory, the evolution power spectrum of the process of ground motion accelerationXt with nonstationary intensity can be expressed as(6)Sxt,w=ft2⋅Sw,where ft denotes the intensity envelope function and Sw is the stationary process power spectral density. Both will be described in detail in the following sections.


### 2.2.1. Intensity Modulation Function

This paper sets out to study multidimensional nonstationary seismic response of structures. The commonly used intensity nonstationary model is selected for the single-dimensional component of the model, and the intensity envelope function is used to describe the strength nonstationary characteristics of the model. Furthermore, a four-stage continuous intensity envelope function model, proposed by Li and Lai [22], is considered to more realistically simulate the correlation of three-dimensional ground motion. In contrast to the popular Amin-Ang model [23], the model is obtained by statistical regression of a large number of measured seismic waves with two-level ground motion components using the least-squares method, thus better reflecting the rising phase, stationary phase, and attenuation phase of the ground motion time history. The intensity envelope function ft is expressed as(7)ft=0,0<t≤t0,t−t0t1−t02,t0<t≤t1,1,t1<t≤t2,exp−ct−t2,t>t2,in which c adjusts the decay speed of the falling portion of the model curve, t0 is the start time of the rising segment, t1 and t2 are the first and last moments of the peak plateau, respectively, T1=t1−t0 is the duration of the rising segment, and TS=t2−t1 is the duration of the plateau. The values of parameters t0,T1,TS, and c used in this paper are from Li and Lai [22].Their research on the multidimensional coherence of the components of ground motion shows that the correlation between the seismic components in the frequency characteristics is relatively low and the attenuation is not obvious [22]. Therefore, in order to simplify the complexity of simulating the correlation of each component, it is feasible to only consider the intensity characteristics in the seismic analysis of structures, without considering the effect of three-dimensional frequency correlation. For this reason, the three-dimensional correlation of the frequency characteristics is not considered in this paper.


### 2.2.2. Power Spectrum

For the power spectrum models of ground motion processes, two main categories are considered in this paper: (1) intensity-frequency power spectrum model; (2) power spectra converted from standard design response spectrum.(1)Modified C-P Model. The Clough–Penzien (C-P) model [24] is an improvement on the classic Kanai–Tajimi model [25]. It can correct the energy of the low-frequency part and is widely used. Accordingly, in order to consider the frequency characteristics in the power spectrum model parameters, a modified C-P model is proposed, which can be expressed as(8)St,w=wg4t+4ζg2twg2tw2w2−wg2t2+4ζg2twg2tw2⋅w4w2−wf2t2+4ζf2twf2tw2⋅S0t,the modified C-P model is able to reflect the nonstationary intensity of the ground motion, while the parameter duration and frequency spectrum reflect the nonstationary characteristics of frequency and can be expressed as(9)wgt=w0−atT,ζgt=ζ0+btT,(10)wft=0.1wgt,ζft=ζgt,where w0, ζ0, a, and b are determined according to the site category and T is the ground motion duration. Similar to the C-P model of the steady ground motion acceleration process, wgt and ζgt can be called the dominant circular frequency and damping ratio of the site soil, respectively, and wft and ζft are the corresponding filter parameters. As we can see from equations (9) and (10), the site parameters and filtering parameters are both linear functions of time, which indicates that wgt and ζgt of the site soil are time-varying parameters that change within a certain range.In equation (8), the spectral parameter S0t reflecting the intensity of ground motion can be expressed as(11)S0t=a¯max2r¯2πwgt2ζgt+1/2ζgt,in which a¯max is the mean value of the peak ground acceleration or the direct ground motion acceleration, according to the specification; r¯ is the equivalent peak value factor.(2)Power Spectrum Converted from Standard Design Response Spectrum. Alternatively, when generating the ground motion time history by standard design response spectrum, an approximate conversion relationship between the acceleration response spectrum and power spectrum is required, defined as(12)Sw=−ζπwSa2w,ζ·1ln−π/wTln1−p,where Sw denotes the bilateral power spectral density function; Saw,ζ represents the absolute acceleration response spectrum; T is the ground motion duration; and p is the transcendence probability of the response spectrum, which can be taken as 0.5.In the process of conversion, the discrete control points are selected first on the design response spectrum, and the logarithmic coordinates are then selected uniformly over the period, expressed as(13)lgT0,k=lgT0,1+1−lgT0,1M−1k−1,k=1,2,…,M,where T0,k is the kth control point of the period and M is the total number of discrete control points in a 10-second period. The total number of discrete control points M is determined by equation (13), and the power spectrum value of the corresponding control point is determined by equation (12). The other power spectrum values are determined by linear interpolation, expressed as(14)Swi=Swk+wi−wkwk−1−wkSwk+1−Swk,i=1,2,…,N,k=1,2,…,M,where Swk and Swk+1 represent the power spectrum values of the kth and k+1 th points and Swi denotes the power spectrum value at the ith point, wi=iΔw.In order to make the conversion power spectrum more accurately fit the code’s design response spectrum, the fitting error is defined as the error between the design response spectrum and the average response spectrum [26] and is expressed as(15)εwk=Sawk−SkwkSawk×100%,where Sawk and Skwk represent the design response spectrum value and the average response spectrum value for the kth point. If the fitting error εwg<10%, the accuracy meets requirements; otherwise, the power spectrum requires iteration and the iteration formula [27] is expressed as(16)Si+1wk+1=SiwkSa2wkSk2wk.Due to the random characteristics of ground motion, the design response spectrum is obtained by statistical averaging and smooth processing of the response spectrum of a large number of actual earthquakes. Therefore, the converted power spectrum has practical engineering application value, and the average meaning can be combined with the randomness method.The converted power spectrum obtained by the response spectrum conversion can be used for the simulation of multidimensional ground motion. Taking the seismic zones 5B site and the response spectrum type I as an example, the duration isT = 30 s, and the peak acceleration value is 0.3 g. As can be seen in Figure 1, for the converted power spectrum obtained from the design response spectrum, horizontal x and horizontal y partially overlap in the low-frequency part, and the changing trend of the three-dimensional conversion power spectrum is basically the same on the whole. However, the peak and fall rates of each component are different, the peak is close to the low-frequency part, and the variation range is 10∼50 Hz. Since the response spectra are all response spectra in an average sense, the conversion power spectrum has an average meaning.Figure 1
Converted power spectrum.


## 2.2.1. Intensity Modulation Function

This paper sets out to study multidimensional nonstationary seismic response of structures. The commonly used intensity nonstationary model is selected for the single-dimensional component of the model, and the intensity envelope function is used to describe the strength nonstationary characteristics of the model. Furthermore, a four-stage continuous intensity envelope function model, proposed by Li and Lai [22], is considered to more realistically simulate the correlation of three-dimensional ground motion. In contrast to the popular Amin-Ang model [23], the model is obtained by statistical regression of a large number of measured seismic waves with two-level ground motion components using the least-squares method, thus better reflecting the rising phase, stationary phase, and attenuation phase of the ground motion time history. The intensity envelope function ft is expressed as(7)ft=0,0<t≤t0,t−t0t1−t02,t0<t≤t1,1,t1<t≤t2,exp−ct−t2,t>t2,in which c adjusts the decay speed of the falling portion of the model curve, t0 is the start time of the rising segment, t1 and t2 are the first and last moments of the peak plateau, respectively, T1=t1−t0 is the duration of the rising segment, and TS=t2−t1 is the duration of the plateau. The values of parameters t0,T1,TS, and c used in this paper are from Li and Lai [22].Their research on the multidimensional coherence of the components of ground motion shows that the correlation between the seismic components in the frequency characteristics is relatively low and the attenuation is not obvious [22]. Therefore, in order to simplify the complexity of simulating the correlation of each component, it is feasible to only consider the intensity characteristics in the seismic analysis of structures, without considering the effect of three-dimensional frequency correlation. For this reason, the three-dimensional correlation of the frequency characteristics is not considered in this paper.


## 2.2.2. Power Spectrum

For the power spectrum models of ground motion processes, two main categories are considered in this paper: (1) intensity-frequency power spectrum model; (2) power spectra converted from standard design response spectrum.(1)Modified C-P Model. The Clough–Penzien (C-P) model [24] is an improvement on the classic Kanai–Tajimi model [25]. It can correct the energy of the low-frequency part and is widely used. Accordingly, in order to consider the frequency characteristics in the power spectrum model parameters, a modified C-P model is proposed, which can be expressed as(8)St,w=wg4t+4ζg2twg2tw2w2−wg2t2+4ζg2twg2tw2⋅w4w2−wf2t2+4ζf2twf2tw2⋅S0t,the modified C-P model is able to reflect the nonstationary intensity of the ground motion, while the parameter duration and frequency spectrum reflect the nonstationary characteristics of frequency and can be expressed as(9)wgt=w0−atT,ζgt=ζ0+btT,(10)wft=0.1wgt,ζft=ζgt,where w0, ζ0, a, and b are determined according to the site category and T is the ground motion duration. Similar to the C-P model of the steady ground motion acceleration process, wgt and ζgt can be called the dominant circular frequency and damping ratio of the site soil, respectively, and wft and ζft are the corresponding filter parameters. As we can see from equations (9) and (10), the site parameters and filtering parameters are both linear functions of time, which indicates that wgt and ζgt of the site soil are time-varying parameters that change within a certain range.In equation (8), the spectral parameter S0t reflecting the intensity of ground motion can be expressed as(11)S0t=a¯max2r¯2πwgt2ζgt+1/2ζgt,in which a¯max is the mean value of the peak ground acceleration or the direct ground motion acceleration, according to the specification; r¯ is the equivalent peak value factor.(2)Power Spectrum Converted from Standard Design Response Spectrum. Alternatively, when generating the ground motion time history by standard design response spectrum, an approximate conversion relationship between the acceleration response spectrum and power spectrum is required, defined as(12)Sw=−ζπwSa2w,ζ·1ln−π/wTln1−p,where Sw denotes the bilateral power spectral density function; Saw,ζ represents the absolute acceleration response spectrum; T is the ground motion duration; and p is the transcendence probability of the response spectrum, which can be taken as 0.5.In the process of conversion, the discrete control points are selected first on the design response spectrum, and the logarithmic coordinates are then selected uniformly over the period, expressed as(13)lgT0,k=lgT0,1+1−lgT0,1M−1k−1,k=1,2,…,M,where T0,k is the kth control point of the period and M is the total number of discrete control points in a 10-second period. The total number of discrete control points M is determined by equation (13), and the power spectrum value of the corresponding control point is determined by equation (12). The other power spectrum values are determined by linear interpolation, expressed as(14)Swi=Swk+wi−wkwk−1−wkSwk+1−Swk,i=1,2,…,N,k=1,2,…,M,where Swk and Swk+1 represent the power spectrum values of the kth and k+1 th points and Swi denotes the power spectrum value at the ith point, wi=iΔw.In order to make the conversion power spectrum more accurately fit the code’s design response spectrum, the fitting error is defined as the error between the design response spectrum and the average response spectrum [26] and is expressed as(15)εwk=Sawk−SkwkSawk×100%,where Sawk and Skwk represent the design response spectrum value and the average response spectrum value for the kth point. If the fitting error εwg<10%, the accuracy meets requirements; otherwise, the power spectrum requires iteration and the iteration formula [27] is expressed as(16)Si+1wk+1=SiwkSa2wkSk2wk.Due to the random characteristics of ground motion, the design response spectrum is obtained by statistical averaging and smooth processing of the response spectrum of a large number of actual earthquakes. Therefore, the converted power spectrum has practical engineering application value, and the average meaning can be combined with the randomness method.The converted power spectrum obtained by the response spectrum conversion can be used for the simulation of multidimensional ground motion. Taking the seismic zones 5B site and the response spectrum type I as an example, the duration isT = 30 s, and the peak acceleration value is 0.3 g. As can be seen in Figure 1, for the converted power spectrum obtained from the design response spectrum, horizontal x and horizontal y partially overlap in the low-frequency part, and the changing trend of the three-dimensional conversion power spectrum is basically the same on the whole. However, the peak and fall rates of each component are different, the peak is close to the low-frequency part, and the variation range is 10∼50 Hz. Since the response spectra are all response spectra in an average sense, the conversion power spectrum has an average meaning.Figure 1
Converted power spectrum.


## 2.3. Multidimensional Nonstationary Ground Motion Model

The main aspects of establishing a multidimensional ground motion model include considering the correlation of each seismic component from the whole; extending the one-dimensional power spectrum to a multidimensional cross-power spectrum matrix; and studying the correlation between various seismic component parameters (such as response spectrum, amplitude, or duration) locally.The key to establishing the multidimensional ground motion model is to determine the cross-spectral density between the seismic components in the cross-power spectrum matrix. However, the research on the coherence function between the seismic components of a single point is still insufficient, and cross-spectral density is difficult to accurately express. Accordingly, this paper proposes a simplified method to build a multidimensional nonstationary ground motion model from a local perspective. Specifically, based on the local simulation of the ground motion process (modulation function, power spectrum, and response spectrum), the three-dimensional response spectrum of the fitted structure is generally used as the control target, and the correlation between the intensity characteristics of the multidimensional seismic components is considered on the basis of the one-dimensional ground motion model.On one hand, the modified four-segment intensity envelope function is applied to consider the multisegment correlation. On the other hand, from the power spectrum perspective, the power spectrum can use the one converted by the three-dimensional response spectrum, and the same conversion method can be used for both horizontal and vertical components. For the bidirectional horizontal component, the power spectrum is the same value and is completely related. For the vertical component, the same formula is expressed as the bidirectional horizontal component. The parameters are specified to be different. The model uses the modified C-P power spectrum; then(17)Sxt,w=Syt,w,(18)Szt,w=wgV4t+4ζgV2twgV2tw2w2−wgV2t2+4ζgV2twgV2tw2⋅w4w2−wfV2t2+4ζfV2twfV2tw2⋅S0t,(19)wgVt=EGwgHt=1.58wgHt,(20)ζgVt=ζgHt,wfVt=0.1wgHt.At the same time, assuming that the vertical component acceleration peak is 2/3 of the horizontal component [28], then the relationship between the vertical and horizontal component spectral intensity factors is(21)S0Vt=0.281S0Ht.


## 3. Simulation and Verification of Proposed Model

### 3.1. Target Multidimensional Design Response Spectrum

According to the Egyptian code of loads (ECP 201 2012) [29], there are two different spectra, namely, type (1) and type (2), based on expected surface-wave magnitude at site. However, the response spectrum curve type (2) is only considered for coastal zones on the Mediterranean Sea (40 km distance from shore). For all other zones throughout Egypt, the response spectrum curve type (1) is specified. In this paper, only the latter case is adopted. ECP 201-2012 specifies the horizontal response spectrum SeT and the vertical response spectrum SveT, which are expressed as(22)SeT=agγ1S1∗TTB2.5η−1,0≤T≤TB,2.5agγ1Sη,TB≤T≤TC,2.5agγ1SηTCT,TC≤T≤TD,2.5agγ1SηTCTDT2,TD≤T≤4s,(23)SveT=avgγ11∗TTB3ηv−1,0≤T≤TB,3avgγ1ηv,TB≤T≤TC,3avgγ1ηvTCT,TC≤T≤TD,3avgγ1ηvTCTDT2,TD≤T≤4s,where T is the vibration period; ag and avg represent the horizontal and vertical ground acceleration, respectively; γ1 is the importance factor; η is the design damping factor depending on the type of structure; ηv is the vertical design damping factor depending on the type of structure; and S is the soil factor. TB,TC, and TD limits apply depending on the soil type. The earthquake impact coefficient curve is developed as depicted in Figure 2. A study by Hernández and López shows that the response spectrum of the bidirectional horizontal ground motion component is not exactly the same, and the horizontal response spectrum component and vertical ground motion component do not have the same proportional coefficient at each period point [30]. Therefore, the maximum value of the response spectrum curve, the characteristic period, and the descent index of the descent segment are all different.Figure 2
Target response spectrum.


### 3.2. Generation of Multidimensional Nonstationary Ground Motion

In order to validate the proposed method, the multidimensional representative time history of the multidimensional ground motion used in the seismic design of the building structure is generated according to the Egyptian code of loads (ECP 201-2012). The nonstationary ground motion process simulation starts from the local simulation of the multidimensional ground motion process (modulation function, power spectrum, and response spectrum), and the three-dimensional response spectrum of the fitted structure is generally used as the control target.The process is described as follows (Figure3):Step 1. Select the discrete representative point of the basic random variable Θ, which is uniformly distributed in 0,2π.Step 2. Generate the orthogonal random variables X¯n,Y¯n.Step 3. Obtain a modified C-P power spectrum with coherent physical significance and then parameter correlation is considered to filter the wave signal.Step 4. Synthesize a four-stage envelope model that considers the multidimensional intensity correlation and eliminate the spillover effect.Step 5. Develop the representative time histories of the nonstationary ground motion model with 30 s length.Step 6. Calculate the average response spectrum.Step 7. Compare the average response spectrum to the generated site-specific response spectrum under consideration. Generally, interpolation is required in the iteration operation until a satisfactory result is achieved. According to the considered design response spectrum we considered, the interpolation can be carried out by dividing the average response spectrum into four segments, and then polynomial least-squares fitting can be applied for each segment.Step 8. Repeat the same process for different representative time histories input and a new time record is generated.Figure 3
Schematics of artificially generated accelerograms.In the method simulated by the proposed model, the time durationT = 30 s, the seismic zones 5B site, and the response spectrum type (1) are considered, and the design basic seismic acceleration is 0.3 g. The values of parameters t0 = 0.3 s, T1 = 6.9 s, Ts = 5.5 s, and c = 0.1 are used for the four-stage continuous intensity envelope function model. The parameter values of the modified C-P power spectrum model are shown in Table 1.Table 1
The parameter value of the power spectrum model.
ParameterHorizontal componentXHorizontal componentYVertical componentZa888b0.150.150.15γ¯2.743.302.00w0S−1211816ζ00.650.750.85At the same time, the relevant values of the input model are as follows: (1) the time stepΔt = 0.01 s; (2) the cutoff frequency wμ = 240 rad/s; (3) the number of expansion terms N = 1600; (4) the frequency interval is 0.15 rad/s; (5) the total number of discrete points of the power spectral density function in the response spectrum conversion power spectrum is equal to the number of expansion terms N = 1600; (6) the total number of control points M = 200, and (7) the lower cutoff frequency in the iterative fitting is wm = 0.6 rad/s.Based on the spectral representation method, the modified C-P model and the multidimensional four-segment intensity envelope function (Amin-Ang model) are used to simulate the multidimensional nonstationary ground motion time history of the complete probability set. The values of the parameters under rare earthquakes are considered as follows:Tg takes the value 0.40 s, and αmax takes the values 1.0 and 0.8 for the two horizontal components; Tg takes the value 0.29 s, and αmax takes the value 0.6 for the vertical component; and the seismic acceleration values are all 0.4 g.As an example, Figures4 and 5 indicate the 41st generated representative samples and the comparison with the target response spectrum before and after one iteration of correction. In the curved part of the response spectrum, the average response spectrum before and after the iteration is able to well meet the target response spectrum. Furthermore, through further comparison of the mean, standard deviation, and mean response spectrum of each component of the sample with the target value (Figures 6 to 8), it can be seen that the overall characteristics of the representative time history are ideally fitted to the target value. This shows that the proposed construction method can simulate a good representative sample of the multidimensional ground motion acceleration process.Figure 4
(a) Three-dimensional representative time history and (b) comparison with three-dimensional target response spectrum before iteration.
(a)(b)Figure 5
(a) Three-dimensional representative time history and (b) comparison with three-dimensional target response spectrum after iteration.
(a)(b)Figure 6
Comparison of the horizontalX component of the sample with the target value. (a) Representative time history, (b) mean, (c) standard deviation, and (d) average response spectrum.
(a)(b)(c)(d)Figure 7
Comparison of the horizontalY component of the sample with the target value. (a) Representative time history, (b) mean, (c) standard deviation, and (d) average response spectrum.
(a)(b)(c)(d)Figure 8
Comparison of the horizontalZ component of the sample with the target value. (a) Representative time history, (b) mean, (c) standard deviation, and (d) average response spectrum.
(a)(b)(c)(d)In addition, as can be seen from Figures4 to 8, the multidimensional representative time history can reflect the intensity nonstationary characteristics of multidimensional ground motion and the multidimensional correlation between the components, while the representative time history of each component shows that the nonstationary correlation of the strength of the multidimensional ground motion is inconsistent: (1) the correlation between the two-way horizontal components is stronger than the correlation between the vertical and horizontal components, and the nonstationary change trend of the two-way horizontal components is almost consistent, showing a strong correlation; (2) for the vertical and horizontal components, the vertical component rises significantly faster and earlier, but the amplitude is much smaller, which is in line with the three-dimensional characteristics of actual ground motion.


## 3.1. Target Multidimensional Design Response Spectrum

According to the Egyptian code of loads (ECP 201 2012) [29], there are two different spectra, namely, type (1) and type (2), based on expected surface-wave magnitude at site. However, the response spectrum curve type (2) is only considered for coastal zones on the Mediterranean Sea (40 km distance from shore). For all other zones throughout Egypt, the response spectrum curve type (1) is specified. In this paper, only the latter case is adopted. ECP 201-2012 specifies the horizontal response spectrum SeT and the vertical response spectrum SveT, which are expressed as(22)SeT=agγ1S1∗TTB2.5η−1,0≤T≤TB,2.5agγ1Sη,TB≤T≤TC,2.5agγ1SηTCT,TC≤T≤TD,2.5agγ1SηTCTDT2,TD≤T≤4s,(23)SveT=avgγ11∗TTB3ηv−1,0≤T≤TB,3avgγ1ηv,TB≤T≤TC,3avgγ1ηvTCT,TC≤T≤TD,3avgγ1ηvTCTDT2,TD≤T≤4s,where T is the vibration period; ag and avg represent the horizontal and vertical ground acceleration, respectively; γ1 is the importance factor; η is the design damping factor depending on the type of structure; ηv is the vertical design damping factor depending on the type of structure; and S is the soil factor. TB,TC, and TD limits apply depending on the soil type. The earthquake impact coefficient curve is developed as depicted in Figure 2. A study by Hernández and López shows that the response spectrum of the bidirectional horizontal ground motion component is not exactly the same, and the horizontal response spectrum component and vertical ground motion component do not have the same proportional coefficient at each period point [30]. Therefore, the maximum value of the response spectrum curve, the characteristic period, and the descent index of the descent segment are all different.Figure 2
Target response spectrum.


## 3.2. Generation of Multidimensional Nonstationary Ground Motion

In order to validate the proposed method, the multidimensional representative time history of the multidimensional ground motion used in the seismic design of the building structure is generated according to the Egyptian code of loads (ECP 201-2012). The nonstationary ground motion process simulation starts from the local simulation of the multidimensional ground motion process (modulation function, power spectrum, and response spectrum), and the three-dimensional response spectrum of the fitted structure is generally used as the control target.The process is described as follows (Figure3):Step 1. Select the discrete representative point of the basic random variable Θ, which is uniformly distributed in 0,2π.Step 2. Generate the orthogonal random variables X¯n,Y¯n.Step 3. Obtain a modified C-P power spectrum with coherent physical significance and then parameter correlation is considered to filter the wave signal.Step 4. Synthesize a four-stage envelope model that considers the multidimensional intensity correlation and eliminate the spillover effect.Step 5. Develop the representative time histories of the nonstationary ground motion model with 30 s length.Step 6. Calculate the average response spectrum.Step 7. Compare the average response spectrum to the generated site-specific response spectrum under consideration. Generally, interpolation is required in the iteration operation until a satisfactory result is achieved. According to the considered design response spectrum we considered, the interpolation can be carried out by dividing the average response spectrum into four segments, and then polynomial least-squares fitting can be applied for each segment.Step 8. Repeat the same process for different representative time histories input and a new time record is generated.Figure 3
Schematics of artificially generated accelerograms.In the method simulated by the proposed model, the time durationT = 30 s, the seismic zones 5B site, and the response spectrum type (1) are considered, and the design basic seismic acceleration is 0.3 g. The values of parameters t0 = 0.3 s, T1 = 6.9 s, Ts = 5.5 s, and c = 0.1 are used for the four-stage continuous intensity envelope function model. The parameter values of the modified C-P power spectrum model are shown in Table 1.Table 1
The parameter value of the power spectrum model.
ParameterHorizontal componentXHorizontal componentYVertical componentZa888b0.150.150.15γ¯2.743.302.00w0S−1211816ζ00.650.750.85At the same time, the relevant values of the input model are as follows: (1) the time stepΔt = 0.01 s; (2) the cutoff frequency wμ = 240 rad/s; (3) the number of expansion terms N = 1600; (4) the frequency interval is 0.15 rad/s; (5) the total number of discrete points of the power spectral density function in the response spectrum conversion power spectrum is equal to the number of expansion terms N = 1600; (6) the total number of control points M = 200, and (7) the lower cutoff frequency in the iterative fitting is wm = 0.6 rad/s.Based on the spectral representation method, the modified C-P model and the multidimensional four-segment intensity envelope function (Amin-Ang model) are used to simulate the multidimensional nonstationary ground motion time history of the complete probability set. The values of the parameters under rare earthquakes are considered as follows:Tg takes the value 0.40 s, and αmax takes the values 1.0 and 0.8 for the two horizontal components; Tg takes the value 0.29 s, and αmax takes the value 0.6 for the vertical component; and the seismic acceleration values are all 0.4 g.As an example, Figures4 and 5 indicate the 41st generated representative samples and the comparison with the target response spectrum before and after one iteration of correction. In the curved part of the response spectrum, the average response spectrum before and after the iteration is able to well meet the target response spectrum. Furthermore, through further comparison of the mean, standard deviation, and mean response spectrum of each component of the sample with the target value (Figures 6 to 8), it can be seen that the overall characteristics of the representative time history are ideally fitted to the target value. This shows that the proposed construction method can simulate a good representative sample of the multidimensional ground motion acceleration process.Figure 4
(a) Three-dimensional representative time history and (b) comparison with three-dimensional target response spectrum before iteration.
(a)(b)Figure 5
(a) Three-dimensional representative time history and (b) comparison with three-dimensional target response spectrum after iteration.
(a)(b)Figure 6
Comparison of the horizontalX component of the sample with the target value. (a) Representative time history, (b) mean, (c) standard deviation, and (d) average response spectrum.
(a)(b)(c)(d)Figure 7
Comparison of the horizontalY component of the sample with the target value. (a) Representative time history, (b) mean, (c) standard deviation, and (d) average response spectrum.
(a)(b)(c)(d)Figure 8
Comparison of the horizontalZ component of the sample with the target value. (a) Representative time history, (b) mean, (c) standard deviation, and (d) average response spectrum.
(a)(b)(c)(d)In addition, as can be seen from Figures4 to 8, the multidimensional representative time history can reflect the intensity nonstationary characteristics of multidimensional ground motion and the multidimensional correlation between the components, while the representative time history of each component shows that the nonstationary correlation of the strength of the multidimensional ground motion is inconsistent: (1) the correlation between the two-way horizontal components is stronger than the correlation between the vertical and horizontal components, and the nonstationary change trend of the two-way horizontal components is almost consistent, showing a strong correlation; (2) for the vertical and horizontal components, the vertical component rises significantly faster and earlier, but the amplitude is much smaller, which is in line with the three-dimensional characteristics of actual ground motion.


## 4. PDEM-Based Reliability Assessment

### 4.1. Fundamentals of PDEM

The PDEM, which is based on the principle of probability preservation and its random event description, was proposed by Li and Chen [31]. As a new perspective on stochastic systems, the PDEM provides a powerful tool for the dynamic analysis, engineering reliability, and structural optimal control of stochastic systems [32]. By integrating the PDEM with the representative sample generated by the multidimensional nonstationary ground motion model proposed in the previous section of this paper, the reliability assessment of engineering structures can proceed.Without loss of generality, the equation of motion for a multidegree-of-freedom structural system under earthquake excitations is considered as follows:(24)MX¨t+CX˙t+fX,X˙,t=ΓX¨gΘ,t,where M, C, and f are mass matrix, damping matrix, and restoring force vector, respectively; X¨t, X˙t, and X denote the acceleration, velocity, and displacement vectors of the structural system, respectively; the mass vector Γ=−MI, and I is the loading influence vector; X¨gΘ,t denotes the earthquake ground motion applied to the structure, and the vector Θ denotes the randomness inherent in both structural parameters and the earthquake excitations X¨gΘ,t.Assume the initial displacement and velocity vectors arex0 and x˙0; then, the initial conditions of equation (24) are(25)Xt0=x0,X˙t0=x˙0.For a well-posed dynamics problem, the solution to equation (24) exists, is unique, and is dependent on a function of Θ as follows:(26)Xt=HΘ,t,X˙t=hΘ,t,where H and h=∂H/∂t are deterministic operators.Generally, the physical quantities of interest in the systemZ (the displacements, velocities, stress, internal forces, etc.) are considered in reliability assessment and can be given as(27)Z=HZΘ,t,Z˙=hZΘ,t.Given that the source randomness factor in the structural dynamic system is completely characterized byΘ, the augmented system Z,Θ is probability preserved. For convenience, the joint probability density function of Z,Θ is recorded as pZΘz,θ,t. According to the principle of conservation of probability, the following formulation is obtained:(28)DDt∫Ωt×ΩΘpZΘz,θ,tdzdθ=0,where Ωt×ΩΘ is the arbitrary domain in Ω×ΩΘ. After a series of mathematical manipulations, the generalized density evolution equations can be derived:(29)∂pZΘz,θ,t∂t+∑ℓ=1mZ˙ℓθ,t∂pZΘz,θ,t∂zℓ=0.In the case ofm = 1, equation (29) reduces to a one-dimensional partial differential equation:(30)∂pΖΘz,θ,t∂t+Z˙θ,t∂pΖΘz,θ,t∂z=0.Then, the initial condition can be written as(31)pZΘz,θ,t|t=t0=δz−z0pΘθ,where z0 denotes the deterministic initial value of Zt and δ⋅ is Dirac’s delta function.Once the initial-value problem (equations (30) and (31)) is solved, the instantaneous probability density function of Zt is then given by(32)pZz,t=∫ΩpZΘz,θ,tdθ.


### 4.2. Assessment of Structural Reliability

Reliability, or conversely failure probability, serves as the critical indexes quantifying the performance of structural engineering. The dynamic reliability of seismic structures has always been a major concern of researchers, and two main theories are well developed. One is diffusion process theory, but it has the disadvantage of not being easy to apply to solving high-dimensional partial differential equations. The other is over-cross theory, which is based on multiple assumptions, making it difficult to control calculation accuracy. However, using PDEM makes it easy to find the probability density function of the equivalent extreme value event, so that the dynamic reliability of the engineering structure can be further obtained, and the shortcomings of the above two theories can be avoided.For the first-passage problem, the reliabilityR against the response index Zt can be generally described as(33)R=PrZΘ,t∈Ωs,t∈0,T,where Ωs is the safe domain and T denotes the duration of time. This means that, during the time duration 0,T, the response will never exceed the boundary of the safety domain. In other words, once the response exceeds the boundary, the structure will fail. If a double boundary condition is considered for most of the practical problems, equation (33) can be given equivalently as(34)R=PrZΘ,t<zb,t∈0,T,where zb is the threshold. Zt can be regarded as a random process, in which the value at each time step can be viewed as a random variable. Thereby, equation (34) will be equivalent to(35)R=Pr∩t∈0,TZΘ,t<zb.SupposeX1,X2,…,XK are k random variables, with Wmax being the maximum value of Xj (1 ≤j ≤ k), which is the so-called the equivalent extreme value. Then, there exists(36)R=Pr∩j=1kXj<zb=PrWmax<zb.Therefore, the equivalent extreme value in equation (35) can be expressed as(37)WmaxΘ,t=maxt∈0,TZt.Finally, the first-passage problem in equation (33) can be replaced by(38)R=PrWmax<zb=∫−∞zbpWmaxwdw,where pWmaxw is the probability density function of the equivalent extreme value. Therefore, as long as the probability density function of the equivalent extreme value defined by equation (37) can be obtained, the first-passage problem can be transformed into a simple one-dimensional integration.


## 4.1. Fundamentals of PDEM

The PDEM, which is based on the principle of probability preservation and its random event description, was proposed by Li and Chen [31]. As a new perspective on stochastic systems, the PDEM provides a powerful tool for the dynamic analysis, engineering reliability, and structural optimal control of stochastic systems [32]. By integrating the PDEM with the representative sample generated by the multidimensional nonstationary ground motion model proposed in the previous section of this paper, the reliability assessment of engineering structures can proceed.Without loss of generality, the equation of motion for a multidegree-of-freedom structural system under earthquake excitations is considered as follows:(24)MX¨t+CX˙t+fX,X˙,t=ΓX¨gΘ,t,where M, C, and f are mass matrix, damping matrix, and restoring force vector, respectively; X¨t, X˙t, and X denote the acceleration, velocity, and displacement vectors of the structural system, respectively; the mass vector Γ=−MI, and I is the loading influence vector; X¨gΘ,t denotes the earthquake ground motion applied to the structure, and the vector Θ denotes the randomness inherent in both structural parameters and the earthquake excitations X¨gΘ,t.Assume the initial displacement and velocity vectors arex0 and x˙0; then, the initial conditions of equation (24) are(25)Xt0=x0,X˙t0=x˙0.For a well-posed dynamics problem, the solution to equation (24) exists, is unique, and is dependent on a function of Θ as follows:(26)Xt=HΘ,t,X˙t=hΘ,t,where H and h=∂H/∂t are deterministic operators.Generally, the physical quantities of interest in the systemZ (the displacements, velocities, stress, internal forces, etc.) are considered in reliability assessment and can be given as(27)Z=HZΘ,t,Z˙=hZΘ,t.Given that the source randomness factor in the structural dynamic system is completely characterized byΘ, the augmented system Z,Θ is probability preserved. For convenience, the joint probability density function of Z,Θ is recorded as pZΘz,θ,t. According to the principle of conservation of probability, the following formulation is obtained:(28)DDt∫Ωt×ΩΘpZΘz,θ,tdzdθ=0,where Ωt×ΩΘ is the arbitrary domain in Ω×ΩΘ. After a series of mathematical manipulations, the generalized density evolution equations can be derived:(29)∂pZΘz,θ,t∂t+∑ℓ=1mZ˙ℓθ,t∂pZΘz,θ,t∂zℓ=0.In the case ofm = 1, equation (29) reduces to a one-dimensional partial differential equation:(30)∂pΖΘz,θ,t∂t+Z˙θ,t∂pΖΘz,θ,t∂z=0.Then, the initial condition can be written as(31)pZΘz,θ,t|t=t0=δz−z0pΘθ,where z0 denotes the deterministic initial value of Zt and δ⋅ is Dirac’s delta function.Once the initial-value problem (equations (30) and (31)) is solved, the instantaneous probability density function of Zt is then given by(32)pZz,t=∫ΩpZΘz,θ,tdθ.


## 4.2. Assessment of Structural Reliability

Reliability, or conversely failure probability, serves as the critical indexes quantifying the performance of structural engineering. The dynamic reliability of seismic structures has always been a major concern of researchers, and two main theories are well developed. One is diffusion process theory, but it has the disadvantage of not being easy to apply to solving high-dimensional partial differential equations. The other is over-cross theory, which is based on multiple assumptions, making it difficult to control calculation accuracy. However, using PDEM makes it easy to find the probability density function of the equivalent extreme value event, so that the dynamic reliability of the engineering structure can be further obtained, and the shortcomings of the above two theories can be avoided.For the first-passage problem, the reliabilityR against the response index Zt can be generally described as(33)R=PrZΘ,t∈Ωs,t∈0,T,where Ωs is the safe domain and T denotes the duration of time. This means that, during the time duration 0,T, the response will never exceed the boundary of the safety domain. In other words, once the response exceeds the boundary, the structure will fail. If a double boundary condition is considered for most of the practical problems, equation (33) can be given equivalently as(34)R=PrZΘ,t<zb,t∈0,T,where zb is the threshold. Zt can be regarded as a random process, in which the value at each time step can be viewed as a random variable. Thereby, equation (34) will be equivalent to(35)R=Pr∩t∈0,TZΘ,t<zb.SupposeX1,X2,…,XK are k random variables, with Wmax being the maximum value of Xj (1 ≤j ≤ k), which is the so-called the equivalent extreme value. Then, there exists(36)R=Pr∩j=1kXj<zb=PrWmax<zb.Therefore, the equivalent extreme value in equation (35) can be expressed as(37)WmaxΘ,t=maxt∈0,TZt.Finally, the first-passage problem in equation (33) can be replaced by(38)R=PrWmax<zb=∫−∞zbpWmaxwdw,where pWmaxw is the probability density function of the equivalent extreme value. Therefore, as long as the probability density function of the equivalent extreme value defined by equation (37) can be obtained, the first-passage problem can be transformed into a simple one-dimensional integration.


## 5. Reliability Assessment of Frame Structures

### 5.1. Analytical Model

In the case study analysis, the PDEM is applied in the reliability assessment of a 10-story reinforced concrete frame structure under stochastic seismic excitation generated by the proposed multidimensional nonstationary ground motion model. The geometric layout of the structure is presented in Figure9. The structure is modeled and analyzed in this paper using SAP2000. A three-dimensional mathematical model of the physical structure is used to represent the spatial distribution of the mass and stiffness of the structure to an extent that is adequate for the calculation of the significant features of the building’s dynamic response. The building is modeled as a three-dimensional frame structure using frame elements for the columns, longitudinal beams and transverse beams, and shell elements for the slabs. The structure is assumed to have a rigid foundation; therefore, the soil foundation interaction and foundation flexibility effects are ignored. The total height of the building is 40 m, and the typical height of all floors is 4 m. The total length of the building is 30 m for both directions (N-S and E-W), each divided into five typical bays, which are 6 m long. The cross sections of the columns and beams are 0.8 m × 0.8 m and 0.35 m × 0.65 m, respectively, with 0.15 m thickness for the roof slabs and 0.2 m for the remaining floors. It is assumed that the diaphragm is rigid. The compressive strengths of concrete, fc′, are 30 MPa for the columns and 25 MPa for the beams and the slab. The yield strength of the reinforcement bars, fy, is 400 MPa. The modulus of elasticity is E = 25 GPa, and the Poisson ratio μ = 0.2. The density of the structures is 24 kN/m3. The structural basic configuration is defined based on a preliminary process using gravity loads only, which are defined based on the building’s function. The natural frequencies of the unyieided building are, respectively, 0.46, 1.45, 3.00, 5.51, and 9.09 Hz.Figure 9
Geometric layout of the reinforced concrete frame structure. (a) Plan view. (b) Elevation.
(a)(b)


### 5.2. Dynamic Characteristic Analysis

After analyzing the structure using a complete set of multidimensional time-series samples, the dimensionless interstory drift of each story is selected for random seismic response analysis of the building. Combined with the probability density evolution method (PDEM), the analysis results of the probability density evolution of the dimensionless interstory drift between structural layers under rare earthquake (X, XY, and XYZ directions) excitations can be obtained.As an example, Figure10 shows the comparison between the PDFs at three typical time instants (at t = 5 s, 7 s, and 9 s) for the third floor of the 10-story building structure (X, XY, and XYZ directions). It can be clearly seen from the figure that the PDFs of the interlayer displacement angle at different times are different, but all the PDFs are relatively discrete, and the maximum interlayer displacement angle response has exceeded 0.004. In addition, the PDFs at different moments have multiple peaks, which do not satisfy the strict Gaussian distribution.Figure 10
The PDFs of the dimensionless interstory drift for the 10-story building at typical time instants. (a)X direction, (b) XY direction, and (c) XYZ direction.
(a)(b)(c)Figure11 presents the evolutionary PDF contour in the time interval [9 s, 10 s], which implies that some strength failure occurs and a probability preserved system is involved. Figure 11(a) is the equal probability density curve, depicting the contour of the probability density function. The shape is exactly like a flowing liquid, from which the probability density changes over time can be seen. Figure 11(b) is a three-dimensional stereogram of probability density evolution and the shape resembles a continuous mountain peak. It can be seen from the figure that the absolute value of the dimensionless interstory drift is basically within 0.004.Figure 11
Probabilistic information of the dimensionless interstory drift for 10-story building underXYZ directions excitation. (a) Contour of PDF surface. (b) PDF evolution surface.
(a)(b)In fact, the application of the PDEM enables rich probability information of any structural dynamic response of interest to be obtained, especially the fine probability density evolution process of the structural dynamic response in any period, which provides an effective way for realizing refined structural dynamic response analysis.


### 5.3. System Reliability Evaluation

Denote the interstory drifts from the floor to the top byX1t,X2t,…,X10t, the heights of the stories by h1,h2,…,h10. The system reliability of the structure can be defined by(39)R=Pr∩j=110Xjthj<φb,t∈0,T,where φb is the threshold of the interstory angle. For clarity, we define the dimensionless interstory drift as(40)X¯jt=Xjthjφb,j=1,2,…,10.Thus, equation (38) becomes(41)R=Pr∩j=110X¯jt<1,t∈0,T=Pr∩j=110X¯j,max<1,where X¯j,max=maxt∈0,TX¯jt. Further, we define an equivalent extreme value by(42)X¯max=max1≤j≤10X¯j,max.The integral on the PDF of this equivalent extreme value random variable will then give the system reliability and the probability of failure; that is,(43)R=PrX¯max<1=∫01pX¯maxxdx,Pf=1−R.The reliability and the probability of failure of each story can be defined by(44)Rj=PrX¯j,max<1=∫01pX¯maxxdx,Pf,j=1−Rj,j=1,2,…,10.The main objective of seismic codes, as indicated in the Egyptian code (ECP-201), is to achieve satisfactory performance of structural systems by avoiding collapse and so injury and loss of life when subjected to earthquake loading. Concrete construction is generally favored in Egypt, and because the design of moment-resisting frames is usually controlled more by drift limitations rather than strength, a reinforced concrete moment-resisting frame building is chosen for the purpose of evaluating code drift requirements. Story drift limitations vary across codes. The ECP 201 specifies three levels of allowable story drift limit, depending on nonstructural elements and their arrangement in the structure. The code specifies limits of 0.5% of the story height for buildings with nonductile nonstructural components, 0.75% of the story height for buildings with ductile nonstructural components, and 1.0% of the story height for buildings with nonstructural components isolated from the lateral deformation of the structure, or buildings with no nonstructural components. As per the code, in this study, the damage limitation requirement is set to satisfy the story drift limit of 1.0% of the story height.Table2 details the dynamic reliability of all floors of the structure and the overall dynamic reliability of the structure. The seismic reliability of the upper part of the structure is about 100%, and the seismic reliability of the bottom of the structure is about 60%. This shows that the damage to the structure will mainly occur in the weak layer at the bottom, while the upper layer of the structure is basically elastic, which can help ensure the safety of the structure. Under rare earthquake excitations in XYZ directions, it also can be seen from the table that the third floor has the lowest reliability, at only 57.65%, followed by the fourth floor, with reliability of 64.88%. Therefore, the third layer is most vulnerable to damage during the earthquake and the overall reliability can be judged based on the reliability of the third weak layer. For an existing engineering structure, the weak point of the structure can be found by the dynamic reliability of each floor, and the weak layer can be reinforced to improve its reliability, thereby improving overall reliability so that the structure can fully meet the seismic performance requirements.Table 2
Reliability results for the 10-story building structure.
Floor no.ReliabilityX-wavesXY-wavesXYZ-waves194.1393.3487.88281.0180.0779.08359.8558.0257.65467.1665.9764.88583.7282.4581.88694.1393.3592.91710010099.978100100100910010010010100100100The PDFs of the equivalent extreme valueX¯max are shown in Figure 12. It can be seen that the PDFs of the displacement angle between floors under the three working conditions do not all follow a regular curve and the overall change trend of the same floor under different working conditions is relatively consistent. The extreme value probability density function of the bottom floor is greater than the extreme value of the upper floor. For example, the probability density function curve of the 3rd layer displacement angle in Figure 12 has a bimodal form. At the same time, the probability density distribution function of the displacement angle between the layers shown in Figure 13 more accurately gives the reliability of each layer. In fact, if different thresholds are selected, the reliability of each layer and the overall reliability of the structure corresponding to the different thresholds can be obtained. According to the Egyptian seismic code, allowable story drift story should not exceed 0.01 times the story height.Figure 12
The probability density function of extreme value of the 10-story building. (a)X direction, (b) XY direction, and (c) XYZ direction.
(a)(b)(c)Figure 13
The CDF of the dimensionless interstory drift story for the 10-story building underXYZ directions excitation.In addition, the overall dynamic reliability of the 10-story frame structure at the third floor for working conditionsXYZ, XY, and X is 57.65%, 58.02%, and 59.85%, respectively. The three-dimensional working condition XYZ is smaller than the two-way working condition XY, and the two-way working condition XY is smaller than the one-way working condition X, meaning that multidimensional seismic action increases the damage to the middle and bottom of the structure. By comparing the maximum drift, the amplification of the second-order moment corresponding to the evolution of plastic deformations is observed and the extent of inelastic response can be easily assessed. The drift between the upper layers is relatively small and the response is only slightly inelastic due to the minor influence of second-order effects. However, the influence of P-Delta effects significantly increases the drift response in a highly inelastic manner at lower layers. The results obtained demonstrate that the dynamic reliability of the moment-resisting frame structures under multidimensional seismic action in those areas is favorable. In general, it is also shown that the seismic design of high-rise moment-resisting frame structures in those high seismic areas requires reasonable consideration of the building height and the effects of multidimensional excitations.


## 5.1. Analytical Model

In the case study analysis, the PDEM is applied in the reliability assessment of a 10-story reinforced concrete frame structure under stochastic seismic excitation generated by the proposed multidimensional nonstationary ground motion model. The geometric layout of the structure is presented in Figure9. The structure is modeled and analyzed in this paper using SAP2000. A three-dimensional mathematical model of the physical structure is used to represent the spatial distribution of the mass and stiffness of the structure to an extent that is adequate for the calculation of the significant features of the building’s dynamic response. The building is modeled as a three-dimensional frame structure using frame elements for the columns, longitudinal beams and transverse beams, and shell elements for the slabs. The structure is assumed to have a rigid foundation; therefore, the soil foundation interaction and foundation flexibility effects are ignored. The total height of the building is 40 m, and the typical height of all floors is 4 m. The total length of the building is 30 m for both directions (N-S and E-W), each divided into five typical bays, which are 6 m long. The cross sections of the columns and beams are 0.8 m × 0.8 m and 0.35 m × 0.65 m, respectively, with 0.15 m thickness for the roof slabs and 0.2 m for the remaining floors. It is assumed that the diaphragm is rigid. The compressive strengths of concrete, fc′, are 30 MPa for the columns and 25 MPa for the beams and the slab. The yield strength of the reinforcement bars, fy, is 400 MPa. The modulus of elasticity is E = 25 GPa, and the Poisson ratio μ = 0.2. The density of the structures is 24 kN/m3. The structural basic configuration is defined based on a preliminary process using gravity loads only, which are defined based on the building’s function. The natural frequencies of the unyieided building are, respectively, 0.46, 1.45, 3.00, 5.51, and 9.09 Hz.Figure 9
Geometric layout of the reinforced concrete frame structure. (a) Plan view. (b) Elevation.
(a)(b)


## 5.2. Dynamic Characteristic Analysis

After analyzing the structure using a complete set of multidimensional time-series samples, the dimensionless interstory drift of each story is selected for random seismic response analysis of the building. Combined with the probability density evolution method (PDEM), the analysis results of the probability density evolution of the dimensionless interstory drift between structural layers under rare earthquake (X, XY, and XYZ directions) excitations can be obtained.As an example, Figure10 shows the comparison between the PDFs at three typical time instants (at t = 5 s, 7 s, and 9 s) for the third floor of the 10-story building structure (X, XY, and XYZ directions). It can be clearly seen from the figure that the PDFs of the interlayer displacement angle at different times are different, but all the PDFs are relatively discrete, and the maximum interlayer displacement angle response has exceeded 0.004. In addition, the PDFs at different moments have multiple peaks, which do not satisfy the strict Gaussian distribution.Figure 10
The PDFs of the dimensionless interstory drift for the 10-story building at typical time instants. (a)X direction, (b) XY direction, and (c) XYZ direction.
(a)(b)(c)Figure11 presents the evolutionary PDF contour in the time interval [9 s, 10 s], which implies that some strength failure occurs and a probability preserved system is involved. Figure 11(a) is the equal probability density curve, depicting the contour of the probability density function. The shape is exactly like a flowing liquid, from which the probability density changes over time can be seen. Figure 11(b) is a three-dimensional stereogram of probability density evolution and the shape resembles a continuous mountain peak. It can be seen from the figure that the absolute value of the dimensionless interstory drift is basically within 0.004.Figure 11
Probabilistic information of the dimensionless interstory drift for 10-story building underXYZ directions excitation. (a) Contour of PDF surface. (b) PDF evolution surface.
(a)(b)In fact, the application of the PDEM enables rich probability information of any structural dynamic response of interest to be obtained, especially the fine probability density evolution process of the structural dynamic response in any period, which provides an effective way for realizing refined structural dynamic response analysis.


## 5.3. System Reliability Evaluation

Denote the interstory drifts from the floor to the top byX1t,X2t,…,X10t, the heights of the stories by h1,h2,…,h10. The system reliability of the structure can be defined by(39)R=Pr∩j=110Xjthj<φb,t∈0,T,where φb is the threshold of the interstory angle. For clarity, we define the dimensionless interstory drift as(40)X¯jt=Xjthjφb,j=1,2,…,10.Thus, equation (38) becomes(41)R=Pr∩j=110X¯jt<1,t∈0,T=Pr∩j=110X¯j,max<1,where X¯j,max=maxt∈0,TX¯jt. Further, we define an equivalent extreme value by(42)X¯max=max1≤j≤10X¯j,max.The integral on the PDF of this equivalent extreme value random variable will then give the system reliability and the probability of failure; that is,(43)R=PrX¯max<1=∫01pX¯maxxdx,Pf=1−R.The reliability and the probability of failure of each story can be defined by(44)Rj=PrX¯j,max<1=∫01pX¯maxxdx,Pf,j=1−Rj,j=1,2,…,10.The main objective of seismic codes, as indicated in the Egyptian code (ECP-201), is to achieve satisfactory performance of structural systems by avoiding collapse and so injury and loss of life when subjected to earthquake loading. Concrete construction is generally favored in Egypt, and because the design of moment-resisting frames is usually controlled more by drift limitations rather than strength, a reinforced concrete moment-resisting frame building is chosen for the purpose of evaluating code drift requirements. Story drift limitations vary across codes. The ECP 201 specifies three levels of allowable story drift limit, depending on nonstructural elements and their arrangement in the structure. The code specifies limits of 0.5% of the story height for buildings with nonductile nonstructural components, 0.75% of the story height for buildings with ductile nonstructural components, and 1.0% of the story height for buildings with nonstructural components isolated from the lateral deformation of the structure, or buildings with no nonstructural components. As per the code, in this study, the damage limitation requirement is set to satisfy the story drift limit of 1.0% of the story height.Table2 details the dynamic reliability of all floors of the structure and the overall dynamic reliability of the structure. The seismic reliability of the upper part of the structure is about 100%, and the seismic reliability of the bottom of the structure is about 60%. This shows that the damage to the structure will mainly occur in the weak layer at the bottom, while the upper layer of the structure is basically elastic, which can help ensure the safety of the structure. Under rare earthquake excitations in XYZ directions, it also can be seen from the table that the third floor has the lowest reliability, at only 57.65%, followed by the fourth floor, with reliability of 64.88%. Therefore, the third layer is most vulnerable to damage during the earthquake and the overall reliability can be judged based on the reliability of the third weak layer. For an existing engineering structure, the weak point of the structure can be found by the dynamic reliability of each floor, and the weak layer can be reinforced to improve its reliability, thereby improving overall reliability so that the structure can fully meet the seismic performance requirements.Table 2
Reliability results for the 10-story building structure.
Floor no.ReliabilityX-wavesXY-wavesXYZ-waves194.1393.3487.88281.0180.0779.08359.8558.0257.65467.1665.9764.88583.7282.4581.88694.1393.3592.91710010099.978100100100910010010010100100100The PDFs of the equivalent extreme valueX¯max are shown in Figure 12. It can be seen that the PDFs of the displacement angle between floors under the three working conditions do not all follow a regular curve and the overall change trend of the same floor under different working conditions is relatively consistent. The extreme value probability density function of the bottom floor is greater than the extreme value of the upper floor. For example, the probability density function curve of the 3rd layer displacement angle in Figure 12 has a bimodal form. At the same time, the probability density distribution function of the displacement angle between the layers shown in Figure 13 more accurately gives the reliability of each layer. In fact, if different thresholds are selected, the reliability of each layer and the overall reliability of the structure corresponding to the different thresholds can be obtained. According to the Egyptian seismic code, allowable story drift story should not exceed 0.01 times the story height.Figure 12
The probability density function of extreme value of the 10-story building. (a)X direction, (b) XY direction, and (c) XYZ direction.
(a)(b)(c)Figure 13
The CDF of the dimensionless interstory drift story for the 10-story building underXYZ directions excitation.In addition, the overall dynamic reliability of the 10-story frame structure at the third floor for working conditionsXYZ, XY, and X is 57.65%, 58.02%, and 59.85%, respectively. The three-dimensional working condition XYZ is smaller than the two-way working condition XY, and the two-way working condition XY is smaller than the one-way working condition X, meaning that multidimensional seismic action increases the damage to the middle and bottom of the structure. By comparing the maximum drift, the amplification of the second-order moment corresponding to the evolution of plastic deformations is observed and the extent of inelastic response can be easily assessed. The drift between the upper layers is relatively small and the response is only slightly inelastic due to the minor influence of second-order effects. However, the influence of P-Delta effects significantly increases the drift response in a highly inelastic manner at lower layers. The results obtained demonstrate that the dynamic reliability of the moment-resisting frame structures under multidimensional seismic action in those areas is favorable. In general, it is also shown that the seismic design of high-rise moment-resisting frame structures in those high seismic areas requires reasonable consideration of the building height and the effects of multidimensional excitations.


## 6. Summary and Conclusions

Nonstationary multidimensional random ground motion has an important role to play in the seismic reliability of engineering structures. To address the difficulty of expressing and solving the cross-power spectrum matrix, a simplified approach was taken to establish a multidimensional nonstationary ground motion model. Starting from the local simulation of the ground motion process (modulation function, power spectrum, and response spectrum), this paper proposed determining nonstationary multidimensional ground motion using parameters according to the Egyptian seismic code for buildings. The three-dimensional response spectrum of the combined structure was the control target. Samples of random ground motion for random seismic responses and dynamic reliability were obtained. The random seismic responses were then combined with the PDEM to carry out seismic reliability analysis of randomly base-excited moment-resisting frame structures. In the numerical analysis, the nonlinear seismic responses and reliability of a 10-story reinforced concrete frame structure were carefully investigated.Based on the analysis conducted by this study, the following conclusions can be drawn:(i)
The analyses showed that the simulation based on the proposed methodology can generate nonstationary, multidimensional random ground motion that matches the response spectrum. The resulting analyses of the nonstationary multidimensional random ground motion simulations demonstrated the capability of the proposed method for dynamic seismic response analysis of engineering structure systems.(ii)
PDEM can provide efficient solutions for the reliability problems encountered in engineering structures, and excellent control in stochastic systems as its satisfactory performance here exposes the logical fundamentals of randomness propagation in physical systems.(iii)
The accuracy of reliability assessment can be guaranteed by specifying a reasonable randomness condition in the process, which can be endorsed for the adjustment between accuracy and efficiency.(iv)
The study revealed that the results are significant as the models built proved useful for making a reasonably accurate reliability assessment for unexpected conditions that were not considered in the design development of structures in the specific seismic area. Thus, the proposed model serves to cover a range of conditions in a more effective way and improve the future design in the area. Its suitability and improvement under variable structural systems will be further explored in the future.



---

*Source: 5592249-2021-04-22.xml*