# Non-Gaussian Stochastic Equivalent Linearization Method for Inelastic Nonlinear Systems with Softening Behaviour, under Seismic Ground Motions

Authors: Francisco L. Silva-González; Sonia E. Ruiz; Alejandro Rodríguez-Castellanos
Journal: Mathematical Problems in Engineering
 (2014)

---

## Abstract

A non-Gaussian stochastic equivalent linearization (NSEL) method for estimating the non-Gaussian response of inelastic non-linear structural systems subjected to seismic ground motions represented as nonstationary random processes is presented. Based on a model that represents the time evolution of the joint probability density function (PDF) of the structural response, mathematical expressions of equivalent linearization coefficients are derived. The displacement and velocity are assumed jointly Gaussian and the marginal PDF of the hysteretic component of the displacement is modeled by a mixed PDF which is Gaussian when the structural behavior is linear and turns into a bimodal PDF when the structural behavior is hysteretic. The proposed NSEL method is applied to calculate the response of hysteretic single-degree-of-freedom systems with different vibration periods and different design displacement ductility values. The results corresponding to the proposed method are compared with those calculated by means of Monte Carlo simulation, as well as by a Gaussian equivalent linearization method. It is verified that the NSEL approach proposed herein leads to maximum structural response standard deviations similar to those obtained with Monte Carlo technique. In addition, a brief discussion about the extension of the method to muti-degree-of-freedom systems is presented.


---

## Body

## 1. Introduction

The method of stochastic equivalent linearization (SEL) is one of the most common methods within the approximate approaches for stochastic dynamic analysis of nonlinear systems. The SEL method has proven to be a versatile and efficient technique from the computation point of view. A great impulse to the method was given by Atalik and Utku [1], who demonstrated that for Gaussian responses the calculation of linearization coefficients can be performed in a very simple way. However, an important deficiency of the method has been found when the response of systems with nonlinear inelastic behavior subjected to intense excitations has been considered as Gaussian [2–6]. The authors [7] have found that when the hypothesis related to the response of single-degree-of-freedom (SDOF) systems with hysteretic behavior has Gaussian distribution, the standard deviation of the displacement response can be underestimated in the order of 45%; consequently, the probability of failure widely deviates from the Monte Carlo simulation results, especially for systems with high displacement ductility demands. This is mainly due to the fact that a Gaussian probability distribution is assumed for all variables; however, the restoring force should lie on a finite region, implying that its probability density function (PDF) is non-Gaussian. Several attempts have been made in the last years to overcome the problem mentioned above [8–10]. For example, Asano and Iwan [11] obtained alternate linearization results for a hysteretic system by using equations that were identical to the usual form on all subsets having a nonzero probability density in the nonlinear system but differed in the zero probability subsets. An ample review related to applications of statistical and equivalent linearization in the analysis of structure and mechanical nonlinear stochastic dynamic systems is found in [12].A non-Gaussian equivalent linearization criterion (NSEL) that uses mathematical expressions for the linearization coefficients is presented herein. In order to derive the non-Gaussian equivalent linearization coefficients, a model appropriately representing the time evolution of the joint probability density function of the structural response of softening systems subjected to seismic ground motions is used. The parameters of the joint PDF depend on a function which measures the non-Gaussianity of the response process.The structural softening models are of particular interest in earthquake engineering because they are commonly used to analyze buildings made with softening materials such as steel, reinforced concrete or masonry, which present saturation of the restoring force when deformation increases [13]. The nonlinearities considered in the present study are only due to material behavior.It is noticed that the NSEL methods mentioned in the first paragraph of this section do not take into account the time evolution character of the joint PDF of the structural response in the way that is formulated here.The mathematical expressions corresponding to the non-Gaussian equivalent linearization coefficients derived here are applied to calculate the response of hysteretic SDOF systems with different vibration periods and design displacement ductility demands, subjected to seismic ground motions that are representative of a nonstationary ergodic random process. The comparison of the results with those obtained by means of Monte Carlo simulation shows that the proposed approach predicts with acceptable accuracy the maximum standard deviation of the structural response of hysteretic SDOF systems.An extension of the proposed method to multi-degree-of-freedom (MDOF) systems is outlined at the end of the paper.


## 2. Equation of Motion

The equation of motion of hysteretic SDOF system studied here is(1)
x
¨
+
2
ξ
ω
x
˙
+
α
2
ω
2
x
+
1
-
α
2
ω
2
z
=
-
a
t
,
where x
¨, x
˙, and x represent the relative acceleration, velocity, and displacement, respectively. ξ
=
c
/
2
m
ω is the fraction of viscous critical damping; ω
=
k
/
m is the circular frequency of vibration of the system; m is the mass; c is the coefficient of viscous damping; k is the initial stiffness. α
2 is the ratio between the postyielding stiffness and the initial stiffness; a
(
t
) represents the ground acceleration; t is the time and z is the hysteretic component of displacement represented here by the Bouc-Wen model [14, 15] which is defined by the following nonlinear differential equation:
(2)
z
˙
=
h
x
˙
,
z
=
A
x
˙
-
β
z
x
˙
z
n
-
1
-
γ
x
˙
z
n
,
where A, β, γ, and n are parameters that define the size and the shape of the hysteretic cycles, n controls the smoothness of the transition from elastic to plastic structural behavior in such a way that the smaller n, the smoother the transition, and n
→
∞ corresponds to a bilinear model. β and γ define the softening or the hardening of the system; the first case corresponds to β
+
γ
>
0 while the second corresponds to β
+
γ
≤
0. Casciati and Faravelli [16], recommend considering β
=
γ for steel and β
=
-
2
γ for reinforced concrete. In this study softening systems with smooth transition represented by n
=
1 are considered.


## 3. Stochastic Model of the Seismic Excitation

The seismic motion is modeled as a nonstationary Gaussian random process generated from a filtered white noise modulated in amplitude. The transfer function proposed by Clough and Penzien [17] is used here; therefore, the power spectral density of the seismic excitation is given by
(3)
S
ω
=
ω
g
4
+
4
ξ
g
2
ω
g
2
ω
2
ω
g
2
-
ω
2
2
+
4
ξ
g
2
ω
g
2
ω
2
×
ω
4
ω
f
2
-
ω
2
2
+
4
ξ
f
2
ω
f
2
ω
2
s
0
.
The square of the amplitude modulating function is [18]
(4)
c
2
t
=
a
t
b
d
+
t
e
Exp
-
k
t
.
In practical applications, the amplitude of the power spectral density of the white noise (
s
0
) and the filter parameters (ω
g, ξ
g, ω
f, and ξ
f) are determined with a nonlinear fitting of the power spectral density estimated by means of the Fourier spectrum of the basic seismic ground motion representative of the excitation process. The parameters a, b, d, e, and k of the modulating function are determined with a nonlinear fitting of the Arias intensity curve of the basic accelerogram [19]. In this study, the excitation process is modulated only in amplitude; to study the effect of the frequency modulation on the nonlinear systems response the reader is referred to [20–23].In the analysis presented below, the seismic record obtained in Mexico City at the Ministry of Communications and Transportation (SCT) during the September 19, 1985 earthquake is used to obtain the parameters of the stochastic model of the seismic excitation. The amplitude of the spectral density of the white noise iss
0
=
7.2776
×
10
-
4 cm2/s3, the filter parameters are ω
g
=
3.1017, ξ
g
=
0.0220, ω
f
=
2.2988, and ξ
f
=
0.0492, and those corresponding to the modulating function are a
=
5.8161
×
10
48 cm2/s4, b
=
-
0.3388, d
=
2.166
×
10
47, e
=
26.461, and k
=
-
0.1258. The bimodal power spectral densities of the seismic motion and its fitted model are shown in Figure 1. The cumulative Arias intensity curve and its fitted model (
I
0
) are in Figure 2.Figure 1
Bimodal power spectral density of the accelerogram and its corresponding fitted model.Figure 2
Cumulative Arias intensity curve(
I
0
) and its corresponding fitted model.


## 4. Modeling the Joint Probability Density Function of the Structural Response

The NSEL approach presented here is based on an appropriate representation of the time evolution of the joint PDF of the structural seismic response of softening systems with smooth transition. Results from Monte Carlo simulation analyses have shown that the PDF of the displacement,f
X
(
x
), and of the velocity, f
X
˙
(
x
˙
), of the structural system is Gaussian, and their evolution in time is as shown in Figures 3 and 4, respectively. However, the probability density of the restoring force and of the displacement hysteretic component, f
Z
(
z
), changes with time because it depends on the nonlinearity of the structural behavior. At the beginning and at the end of the excitation, when the structure is under low seismic intensities, f
Z
(
z
) tends to be Gaussian because the structure presents linear behavior; however, as the seismic intensity increases and the structural behavior becomes nonlinear (hysteretic), then f
Z
(
z
) becomes non-Gaussian, as shown in Figure 5. It can be seen that when the system is vibrating within the nonlinear range z equals its maximum value z
u and its PDF presents concentrations at +z
u and −z
u. Figure 5 also shows that as the excitation intensity grows (up to t
=
59 s) and consequently the displacement ductility demand increases, f
Z
(
z
) turns into a bimodal PDF. It is observed that the support of f
Z
(
z
) is −z
u
≤
z
≤
z
u, where the maximum value of z is given by [13, 19]
(5)
z
u
=
A
β
+
γ
1
/
n
.
Considering the concentration of values of z in the vicinity of its maximum value z
u and the shape that such PDF adopts (see Figure 5 for t
=
59 s), the following mathematical model is used:
(6)
f
Z
z
=
φ
a
Z
z
1
-
2
p
+
φ
b
Z
z
-
z
z
p
+
φ
b
Z
z
+
z
z
p
;
0
≤
p
≤
0.5
,
where φ
a
Z is a Gaussian density function with zero mean and standard deviation σ
z
a, φ
b
Z is a Gaussian density function with mean z
z (or -
z
z) and standard deviation σ
z
b, and p is a weighting factor. The idea is that when the nonlinear behavior becomes important, the weighting of φ
b
Z, functions should increase while that corresponding to φ
a
Z should decrease. It is noticed that (6) uses the Gaussian density functions φ
a
Z and φ
b
Z and not Dirac’s pulses as has been proposed by other authors [9]. The advantage of using Gaussian density functions, instead of Dirac’s pulses, is that it is possible to model more appropriately the time evolution of the joint PDF and, as a consequence, the NSEL approach leads to results that better approximate the Monte Carlo results. Furthermore, it is possible to get mathematical expressions of the equivalent linearization coefficients.Figure 3
Time evolution of the probability density function of the displacement response.Figure 4
Time evolution of the probability density function of the velocity response.Figure 5
Time evolution of the probability density function of the hysteretic component of displacement response.Using the marginal PDFsf
X
(
x
), f
X
˙
(
x
˙
), and the information about the correlation structure of the random variables, it is possible to compute the joint PDF [24, 25]. Here, the joint PDF is obtained using the conditional PDF of x and x
˙ given z, h
(
x
,
x
˙
∣
z
), as follows:
(7)
f
X
X
˙
Z
x
,
x
˙
,
z
=
h
x
,
x
˙
∣
z
f
Z
z
.
Based on the simulation results, it is reasonable to assume that h
(
x
,
x
˙

∣

z
) is jointly Gaussian. The parameters of h
(
x
,
x
˙

∣

z
), that is, conditional means and variances, are estimated assuming that the former are linear functions of z, and the latter do not depend on z, as it happens in a mean square estimation of normal variables [26]. Additionally, it is considered that the response process has zero mean.Based on (6) and (7), the following expression for the joint PDF is obtained:
(8)
f
X
X
˙
Z
x
,
x
˙
,
z
=
1
-
2
p
2
π
3
/
2
σ
x
σ
x
˙
σ
z
a
κ
1
/
2
×
Exp
-
1
2
κ
z
2
σ
z
a
2
1
-
ρ
x
x
˙
2
mmmmmmmmmm
+
z
2
σ
z
2
-
z
2
σ
z
a
2
mmmmmmmmmm
×
ρ
x
z
2
+
ρ
x
˙
z
2
-
2
ρ
x
x
˙
ρ
x
z
ρ
x
˙
z
mmmmmmmmmm
+
x
2
σ
x
2
1
-
ρ
x
˙
z
2
+
x
˙
2
σ
x
˙
2
1
-
ρ
x
z
2
mmmmmmmmmm
-
2
x
z
σ
x
σ
z
ρ
x
z
-
ρ
x
x
˙
ρ
x
˙
z
mmmmmmmmmm
-
2
x
˙
z
σ
x
˙
σ
z
ρ
x
˙
z
-
ρ
x
x
˙
ρ
x
z
mmmmmmmnmm
-
2
x
x
˙
σ
x
˙
σ
x
ρ
x
x
˙
-
ρ
x
˙
z
ρ
x
z

1
2
κ
z
2
σ
z
a
2
1
-
ρ
x
x
˙
2
+
p
2
π
3
/
2
σ
x
σ
x
˙
σ
z
b
κ
1
/
2
×
Exp
-
1
2
κ
z
-
z
z
2
σ
z
b
2
1
-
ρ
x
x
˙
2
mmmmmmmmmm
+
z
2
σ
z
2
-
z
-
z
z
2
σ
z
b
2
mmmmmmmmmm
×
ρ
x
z
2
+
ρ
x
˙
z
2
-
2
ρ
x
x
˙
ρ
x
z
ρ
x
˙
z
mmmmmmmmmm
+
x
2
σ
x
2
1
-
ρ
x
˙
z
2
+
x
˙
2
σ
x
˙
2
1
-
ρ
x
z
2
mmmmmmmmmm
-
2
x
z
σ
x
σ
z
ρ
x
z
-
ρ
x
x
˙
ρ
x
˙
z
mmmmmmmmmm
-
2
x
˙
z
σ
x
˙
σ
z
ρ
x
˙
z
-
ρ
x
x
˙
ρ
x
z
mmmmmmmnmm
-
2
x
x
˙
σ
x
˙
σ
x
ρ
x
x
˙
-
ρ
x
˙
z
ρ
x
z

z
-
z
z
2
σ
z
b
2

1
2
κ
z
-
z
z
2
σ
z
b
2
1
-
ρ
x
x
˙
2
+
p
2
π
3
/
2
σ
x
σ
x
˙
σ
z
b
κ
1
/
2
×
Exp
-
1
2
κ
z
+
z
z
2
σ
z
b
2
1
-
ρ
x
x
˙
2
mmmmmmmmmm
+
z
2
σ
z
2
-
z
+
z
z
2
σ
z
b
2
mmmmmmmmmm
×
ρ
x
z
2
+
ρ
x
˙
z
2
-
2
ρ
x
x
˙
ρ
x
z
ρ
x
˙
z
mmmmmmmmmm
+
x
2
σ
x
2
1
-
ρ
x
˙
z
2
+
x
˙
2
σ
x
˙
2
1
-
ρ
x
z
2
mmmmmmmmmm
-
2
x
z
σ
x
σ
z
ρ
x
z
-
ρ
x
x
˙
ρ
x
˙
z
mmmmmmmmmm
-
2
x
˙
z
σ
x
˙
σ
z
ρ
x
˙
z
-
ρ
x
x
˙
ρ
x
z
mmmmmmmnmm
-
2
x
x
˙
σ
x
˙
σ
x
ρ
x
x
˙
-
ρ
x
˙
z
ρ
x
z

z
+
z
z
2
σ
z
b
2

1
2
κ
z
+
z
z
2
σ
z
b
2
1
-
ρ
x
x
˙
2
,
where σ
x, σ
x
˙, and σ
z are the standard deviations of x, x
˙, and z, respectively; ρ
x
z, ρ
x
˙
z, and ρ
x
x
˙ are the corresponding correlation coefficients; and κ
=
1
-
ρ
x
z
2
-
ρ
x
˙
z
2
-
ρ
x
x
˙
2
+
2
ρ
x
x
˙
ρ
x
z
ρ
x
˙
z.


### 4.1. Evaluating the Parametersσ
z
a, σ
z
b, p, and z
z

The four parametersσ
z
a, σ
z
b, p, and z
z of the joint PDF given by (8) should satisfy that the area under the PDF given by (6) is close to unity:
(9)
∫
-
z
u
z
u
f
Z
z
d
z
=
1
-
2
p
Erf
z
z
2
σ
z
a
+
p
Erf
z
u
-
z
z
2
σ
z
b
mim
+
Erf
z
u
+
z
z
2
σ
z
b
≈
1
,
where Erf
(
·
) is the error function given by
(10)
Erf
u
=
2
π
∫
0
u
Exp
-
t
2
d
t
.
The above mentioned condition is necessary because (6) is not a bounded function. If (9) is satisfied, the second moment (in this case, equal to the variance σ
z
2) can be approximated as
(11)
E
z
2
=
∭
-
∞
∞
z
2
f
X
X
˙
Z
x
,
x
˙
,
z
d
x
d
x
˙
d
z
=
∫
-
∞
∞
z
2
f
Z
z
d
z
=
1
-
2
p
σ
z
a
2
+
2
p
z
z
2
+
σ
z
b
2
=
σ
z
2
.
The four parameters σ
z
a, σ
z
b, p, and z
z could be computed by means of their first four moments; however, here they are determined by means of nondimensional coefficients that depend on the system response non-Gaussianity level which is represented by [13]
(12)
λ
t
=
∫
-
∞
-
z
u
1
2
π
σ
z
t
Exp
-
z
2
2
σ
z
2
t
d
z
=
1
2
1
-
Erf
z
u
2
σ
z
t
.
From (12) it is observed that λ is close to zero when the response is Gaussian, that is, when the system strength is high (where z
u is high) and the excitation or response is low (and as a consequence σ
z is low). On the other hand, if the response is non-Gaussian, that is, if the system strength is low (where z
u is low) and the excitation or response is high (and consequently σ
z is high), then λ is high.Here, the following mathematical expressions for the normalized parametersz
z
/
z
z
u, σ
z
a
/
σ
z, σ
z
b
/
σ
z, and p are proposed:
(13)
z
z
t
z
u
=
a
1
Ln
a
2
λ
t
+
a
3
+
a
4
,
p
t
=
b
1
Ln
b
2
λ
t
+
b
3
+
b
4
,
σ
z
a
t
σ
z
t
=
c
1
Ln
c
2
λ
t
+
c
3
+
c
4
,
σ
z
b
t
σ
z
t
=
d
1
Ln
d
2
λ
t
+
d
3
+
d
4
,
where a
i, b
i, c
i, and d
i, i
=
1,2
,
3 and 4, are constants to be determined. The general behavior of (13) is presented in Figures 6(a)–6(d), where it can be seen that z
z and p increase with λ (Figures 6(a) and 6(b)), meaning that when the nonlinearity of the response increases, the function φ
Z
b moves far from z
=
0 and its weight increases; hence, f
Z
z (given by (6)) becomes a bimodal PDF. On the other hand, when the structural behavior is linear λ is close to zero and the weight of φ
Z
b vanishes; therefore f
Z
z becomes Gaussian. It can also be seen (Figures 6(c) and 6(d)) that φ
Z
a
/
σ
z and φ
Z
b
/
σ
z decrease with λ, but it should be noticed that as σ
z increases, λ grows as well.Behavior of the functionsz
z
/
z
z
u, p, σ
z
a
/
σ
z, and σ
z
b
/
σ
z versus λ.
(a)
z
z
/
z
u versus λ
(b)
p versus λ
(c)
σ
z
a
/
σ
z versus λ
(d)
σ
z
b
/
σ
z versus λEquations (13) show that the parameters σ
z
a, σ
z
b, p, and z
z are time dependent; thus, the proposed approach considers the time evolutionary character of the joint PDF of the response, while other methods proposed in the literature do not take it into account. Figure 7 shows an example of the time evolution of f
Z
(
z
) given by (6).Figure 7
Example of time evolution off
Z
(
z
).


#### 4.1.1. Case of Narrow-Band Seismic Inputs

Parametersa
i, b
i, c
i, and d
i, i
=
1,2
,
3 and 4 of  (13) corresponding to the response of single-degree-of-freedom systems (SDOF) systems subjected to the action of narrow-band seismic inputs, were calculated by means of Monte Carlo simulation analysis. The critical case, where the natural vibration period of SDOF system is equal to the dominant period of seismic input (in this case T
=
2.1 s, see Figure 1) and the SDOF system has a high design ductility demand η
=
4, was considered. The design ductility factor η is defined here as the ratio between the expected maximum and the yield displacement of the system. In this study, the first one is found by means of a stationary Gaussian equivalent linearization analysis, and the latter is estimated by means of the initial stiffness and the yield force of the system.The SDOF system was subjected to the action of 50,000 artificial accelerograms based on the narrow-band motion recorded in SCT on September 19, 1985 (see Section3). For this case, it was found that λ lies within the interval 0
<
λ
≤
0.1 and that σ
z
a
>
σ
z
b leads to good results. The parameter values obtained from the analysis are shown in Table 1.Table 1
Values of the parameters in (13).
a
1
a
2
a
3
a
4
0.187
31.990
0.042
0.590
b
1
b
2
b
3
b
4
0.043
455.126
0.003
0.240
c
1
c
2
c
3
c
4
−0.269
8.247
0.125
0.442
d
1
d
2
d
3
d
4
−0.239
1.408
0.007
−0.391Pointsc
1
–
c
6 in Figures 6(a)–6(d) represent the results of Monte Carlo simulation analyses.It is noticed that parameters in Table1 are proposed for narrow-band seismic inputs and are independent of any other parameter, such as structural vibration period and design ductility factor.


## 4.1. Evaluating the Parametersσ
z
a, σ
z
b, p, and z
z

The four parametersσ
z
a, σ
z
b, p, and z
z of the joint PDF given by (8) should satisfy that the area under the PDF given by (6) is close to unity:
(9)
∫
-
z
u
z
u
f
Z
z
d
z
=
1
-
2
p
Erf
z
z
2
σ
z
a
+
p
Erf
z
u
-
z
z
2
σ
z
b
mim
+
Erf
z
u
+
z
z
2
σ
z
b
≈
1
,
where Erf
(
·
) is the error function given by
(10)
Erf
u
=
2
π
∫
0
u
Exp
-
t
2
d
t
.
The above mentioned condition is necessary because (6) is not a bounded function. If (9) is satisfied, the second moment (in this case, equal to the variance σ
z
2) can be approximated as
(11)
E
z
2
=
∭
-
∞
∞
z
2
f
X
X
˙
Z
x
,
x
˙
,
z
d
x
d
x
˙
d
z
=
∫
-
∞
∞
z
2
f
Z
z
d
z
=
1
-
2
p
σ
z
a
2
+
2
p
z
z
2
+
σ
z
b
2
=
σ
z
2
.
The four parameters σ
z
a, σ
z
b, p, and z
z could be computed by means of their first four moments; however, here they are determined by means of nondimensional coefficients that depend on the system response non-Gaussianity level which is represented by [13]
(12)
λ
t
=
∫
-
∞
-
z
u
1
2
π
σ
z
t
Exp
-
z
2
2
σ
z
2
t
d
z
=
1
2
1
-
Erf
z
u
2
σ
z
t
.
From (12) it is observed that λ is close to zero when the response is Gaussian, that is, when the system strength is high (where z
u is high) and the excitation or response is low (and as a consequence σ
z is low). On the other hand, if the response is non-Gaussian, that is, if the system strength is low (where z
u is low) and the excitation or response is high (and consequently σ
z is high), then λ is high.Here, the following mathematical expressions for the normalized parametersz
z
/
z
z
u, σ
z
a
/
σ
z, σ
z
b
/
σ
z, and p are proposed:
(13)
z
z
t
z
u
=
a
1
Ln
a
2
λ
t
+
a
3
+
a
4
,
p
t
=
b
1
Ln
b
2
λ
t
+
b
3
+
b
4
,
σ
z
a
t
σ
z
t
=
c
1
Ln
c
2
λ
t
+
c
3
+
c
4
,
σ
z
b
t
σ
z
t
=
d
1
Ln
d
2
λ
t
+
d
3
+
d
4
,
where a
i, b
i, c
i, and d
i, i
=
1,2
,
3 and 4, are constants to be determined. The general behavior of (13) is presented in Figures 6(a)–6(d), where it can be seen that z
z and p increase with λ (Figures 6(a) and 6(b)), meaning that when the nonlinearity of the response increases, the function φ
Z
b moves far from z
=
0 and its weight increases; hence, f
Z
z (given by (6)) becomes a bimodal PDF. On the other hand, when the structural behavior is linear λ is close to zero and the weight of φ
Z
b vanishes; therefore f
Z
z becomes Gaussian. It can also be seen (Figures 6(c) and 6(d)) that φ
Z
a
/
σ
z and φ
Z
b
/
σ
z decrease with λ, but it should be noticed that as σ
z increases, λ grows as well.Behavior of the functionsz
z
/
z
z
u, p, σ
z
a
/
σ
z, and σ
z
b
/
σ
z versus λ.
(a)
z
z
/
z
u versus λ
(b)
p versus λ
(c)
σ
z
a
/
σ
z versus λ
(d)
σ
z
b
/
σ
z versus λEquations (13) show that the parameters σ
z
a, σ
z
b, p, and z
z are time dependent; thus, the proposed approach considers the time evolutionary character of the joint PDF of the response, while other methods proposed in the literature do not take it into account. Figure 7 shows an example of the time evolution of f
Z
(
z
) given by (6).Figure 7
Example of time evolution off
Z
(
z
).


### 4.1.1. Case of Narrow-Band Seismic Inputs

Parametersa
i, b
i, c
i, and d
i, i
=
1,2
,
3 and 4 of  (13) corresponding to the response of single-degree-of-freedom systems (SDOF) systems subjected to the action of narrow-band seismic inputs, were calculated by means of Monte Carlo simulation analysis. The critical case, where the natural vibration period of SDOF system is equal to the dominant period of seismic input (in this case T
=
2.1 s, see Figure 1) and the SDOF system has a high design ductility demand η
=
4, was considered. The design ductility factor η is defined here as the ratio between the expected maximum and the yield displacement of the system. In this study, the first one is found by means of a stationary Gaussian equivalent linearization analysis, and the latter is estimated by means of the initial stiffness and the yield force of the system.The SDOF system was subjected to the action of 50,000 artificial accelerograms based on the narrow-band motion recorded in SCT on September 19, 1985 (see Section3). For this case, it was found that λ lies within the interval 0
<
λ
≤
0.1 and that σ
z
a
>
σ
z
b leads to good results. The parameter values obtained from the analysis are shown in Table 1.Table 1
Values of the parameters in (13).
a
1
a
2
a
3
a
4
0.187
31.990
0.042
0.590
b
1
b
2
b
3
b
4
0.043
455.126
0.003
0.240
c
1
c
2
c
3
c
4
−0.269
8.247
0.125
0.442
d
1
d
2
d
3
d
4
−0.239
1.408
0.007
−0.391Pointsc
1
–
c
6 in Figures 6(a)–6(d) represent the results of Monte Carlo simulation analyses.It is noticed that parameters in Table1 are proposed for narrow-band seismic inputs and are independent of any other parameter, such as structural vibration period and design ductility factor.


## 4.1.1. Case of Narrow-Band Seismic Inputs

Parametersa
i, b
i, c
i, and d
i, i
=
1,2
,
3 and 4 of  (13) corresponding to the response of single-degree-of-freedom systems (SDOF) systems subjected to the action of narrow-band seismic inputs, were calculated by means of Monte Carlo simulation analysis. The critical case, where the natural vibration period of SDOF system is equal to the dominant period of seismic input (in this case T
=
2.1 s, see Figure 1) and the SDOF system has a high design ductility demand η
=
4, was considered. The design ductility factor η is defined here as the ratio between the expected maximum and the yield displacement of the system. In this study, the first one is found by means of a stationary Gaussian equivalent linearization analysis, and the latter is estimated by means of the initial stiffness and the yield force of the system.The SDOF system was subjected to the action of 50,000 artificial accelerograms based on the narrow-band motion recorded in SCT on September 19, 1985 (see Section3). For this case, it was found that λ lies within the interval 0
<
λ
≤
0.1 and that σ
z
a
>
σ
z
b leads to good results. The parameter values obtained from the analysis are shown in Table 1.Table 1
Values of the parameters in (13).
a
1
a
2
a
3
a
4
0.187
31.990
0.042
0.590
b
1
b
2
b
3
b
4
0.043
455.126
0.003
0.240
c
1
c
2
c
3
c
4
−0.269
8.247
0.125
0.442
d
1
d
2
d
3
d
4
−0.239
1.408
0.007
−0.391Pointsc
1
–
c
6 in Figures 6(a)–6(d) represent the results of Monte Carlo simulation analyses.It is noticed that parameters in Table1 are proposed for narrow-band seismic inputs and are independent of any other parameter, such as structural vibration period and design ductility factor.


## 5. Non-Gaussian Equivalent Stochastic Linearization

In the equivalent linearization method, (2) is replaced by a linear differential equation as follows [5, 27]:
(14)
z
˙
=
s
e
x
+
c
e
x
˙
+
k
e
z
,
where s
e, c
e, and k
e are linearization coefficients. In order to minimize the expected value of the square of the difference of (2) and (14), the linearization coefficients must satisfy [5, 28, 29]
(15)
H
e
=
E
Q
Q
T
-
1
E
h
Q
,
where Q is the response vector Q
=
x
x
˙
z
T, H
e is the vector of equivalent linearization coefficients given by H
e
=
s
e
c
e
k
e
T, and h is given by (2).When the joint PDF proposed in this paper (see (8)) is substituted in (15), mathematical expressions for the linearization coefficients are derived (see Appendix A). In particular, for n
=
1, which corresponds to SDOF systems with softening behavior, the equivalent linearization coefficients are as follows:
(16)
s
e
=
0
,
c
e
=
(
1
-
2
p
)
A
-
2
π
β
ρ
x
˙
z
σ
z
a
2
ε
a
-
2
π
γ
σ
z
a
+
2
p
A
-
β
2
π
ρ
x
˙
z
σ
z
b
2
ε
b
Exp
-
z
z
2
ρ
x
˙
z
2
2
ε
b
mmmmmmmmm
+
z
z
Erf
z
z
ρ
x
˙
z
2
ε
b
mmmmm
-
γ
2
π
σ
z
b
Exp
-
z
z
2
2
σ
z
b
2
mmmmmmiimm
+
z
z
Erf
z
z
2
σ
z
b
,
k
e
=
(
1
-
2
p
)
A
ρ
x
˙
z
σ
x
˙
σ
z
a
2
-
σ
z
2
σ
z
3
mmmiimmm
-
2
π
β
σ
x
˙
σ
z
a
2
σ
z
3
ε
a
+
ρ
x
˙
z
2
σ
z
a
2
-
σ
z
2
ε
a
mmmmiimm
-
2
π
γ
ρ
x
˙
z
σ
x
˙
σ
z
a
2
σ
z
a
2
-
σ
z
2
σ
z
3
+
2
p
A
ρ
x
˙
z
σ
x
˙
z
z
2
+
σ
z
b
2
-
σ
z
2
σ
z
3
mmmmm
-
β
σ
x
˙
σ
z
3
mmmmm
×
2
π
Exp
-
z
z
2
ρ
x
˙
z
2
2
ε
b
mmmmmiimmmmm
×
z
z
2
ε
b
+
σ
z
b
2
ε
b
+
ρ
x
˙
z
2
σ
z
b
2
-
σ
z
2
ε
b
mmmmmmmm
+
z
z
ρ
x
˙
z
z
z
2
+
3
σ
z
b
2
-
σ
z
2
Erf
z
z
ρ
x
˙
z
2
ε
b
mmmmm
+
γ
ρ
x
˙
z
σ
x
˙
σ
z
3
miimmm
×
2
π
Exp
-
z
z
2
2
σ
z
b
2
σ
z
b
z
z
2
+
2
σ
z
b
2
-
σ
z
2
mmmmmm
+
z
z
z
z
2
+
3
σ
z
b
2
-
σ
z
2
Erf
z
z
2
σ
z
b
,
where ε
a
=
(
1
-
ρ
x
˙
z
2
)
σ
z
2
+
ρ
x
˙
z
2
σ
z
a
2 and ε
b
=
(
1
-
ρ
x
˙
z
2
)
σ
z
2
+
ρ
x
˙
z
2
σ
z
b
2.


## 6. Covariance of the Response of the Hysteretic System

The covariance matrix of the response of the hysteretic SDOF system,Σ
t
=
E
u
u
T where u
(
t
)
=
x
x
˙
z
x
f
x
g
x
˙
f
x
˙
g
T, is calculated by solving the following equation [28, 29]:
(17)
d
Σ
t
d
t
=
H
t
Σ
t
+
Σ
t
H
T
t
+
2
π
S
F
,
where H is a matrix depending on the mechanical and geometrical properties of the system, the linearization coefficients, the modulating function coefficients, and the filter parameters. In this study, H is given by(18)
H
t
=
0
1
0
0
0
0
0
-
α
2
ω
2
-
2
ξ
ω
-
1
-
α
2
ω
2
ω
f
2
c
t
-
ω
g
2
c
t
2
ξ
f
ω
f
c
t
-
2
ξ
g
ω
g
c
t
s
e
c
e
k
e
0
0
0
0
0
0
0
0
0
1
0
0
0
0
0
0
0
1
0
0
0
-
ω
f
2
ω
g
2
-
2
ξ
f
ω
f
2
ξ
g
ω
g
0
0
0
0
-
ω
g
2
0
-
2
ξ
g
ω
g
.S
F is a 7
×
7 matrix of zeros, except the element S
F
(
7,7
) which is equal to the amplitude of the white noise spectral density (
s
0
).In order to computeΣ
t, it is essential to start the integration process of (17) when the Clough-Penzien filter response has reached the stationary stage. Thus, the initial condition for Σ
t should correspond to the covariance of the stationary response of the Clough-Penzien filter as described in Appendix B [30].


## 7. Verification of the Accuracy of the Proposed NSEL Approach

In this section the standard deviation of the response of several hysteretic SDOF systems is calculated by means of the proposed NSEL criterion using the mathematical expressions obtained above ((13) and (16) using values of Table 1). The SDOF systems have vibration periods T
=
2.1 s and T
=
4.0 s and design ductility factors η
=
1.5 and η
=
4, and they are subjected to a nonstationary random process with the statistical properties of the seismic motion recorded in Mexico City at the SCT during the September 19, 1985 earthquake (see Section 3).The response of the SDOF hysteretic systems is calculated by means of SEL and alternatively by Monte Carlo simulation using 50,000 simulated accelerograms. Two SEL criteria are used: (a) assuming Gaussian response [1] and (b) assuming non-Gaussian response using the criterion and the mathematical expressions proposed herein.The statistical response of the SDOF system withT
=
2.1 s and η
=
1.5 is shown in Figures 8, 9, and 10, and that corresponding to the system with T
=
2.1 s and η
=
4.0 is in Figures 11, 12, and 13. Figures 8 to 13 show that the proposed approach predicts with reasonable accuracy the maximum response; however, as well as the Gaussian method, it does not predict the permanent drift of the hysteretic system especially for high values of the design ductility (see Figures 8 and 11). The statistical response of the SDOF system with T
=
4.0 s and η
=
1.5 is shown in Figures 14, 15, and 16, and that corresponding to the system with T
=
4.0 s and η
=
4.0 is in Figures 17, 18, and 19. As in the case of SDOF system with T
=
2.1 s, the proposed approach predicts with reasonable accuracy the maximum response of the SDOF system with T
=
4.0 s; however it does not predict with good accuracy the permanent drift.Figure 8
Standard deviation of displacement.T
=
2.1 s and η
=
1.5.Figure 9
Standard deviation of velocity.T
=
2.1 s and η
=
1.5.Figure 10
Standard deviation of hysteretic variable.T
=
2.1 s and η
=
1.5.Figure 11
Standard deviation of displacement.T
=
2.1 s and η
=
4.Figure 12
Standard deviation of velocity.T
=
2.1 s and η
=
4.Figure 13
Standard deviation of hysteretic variable.T
=
2.1 s and η
=
4.Figure 14
Standard deviation of displacement.T
=
4.0 s and η
=
1.5.Figure 15
Standard deviation of the velocity.T
=
4.0 s and η
=
1.5.Figure 16
Standard deviation of the hysteretic variable.T
=
4.0 s and η
=
1.5.Figure 17
Standard deviation of displacement.T
=
4.0 s and η
=
4.0.Figure 18
Standard deviation of the velocity.T
= 4.0 s and η
= 4.0.Figure 19
Standard deviation of the hysteretic variable.T
=
4.0 s and η
=
4.0.The relative error (ε) in the calculation of the maximum standard deviation of the response (σ
x, σ
x
˙ and σ
z) is obtained as follows:(19)
ε
%
=
response
obtained
from
linearization
-
response
obtained
from
simulation
response
obtained
from
simulation
×
100
.The relative errors(
ε
) of the maximum response, corresponding to all the cases analyzed, are shown in Tables 2 and 3 corresponding to structural vibration periods of 2.1 and 4.0 s, respectively. A negative value of ε indicates that the equivalent linearization underestimates the response. Tables 2 and 3 show that the proposed NSEL method leads to results close to those obtained by means of Monte Carlo simulation.Table 2
Relative errors of the maximum response of SDOF system withT
=
2.1 s.
Method
Error in:max
⁡
[
σ
x
]
Error in:max
⁡
[
σ
x
˙
]
Error in:max
⁡
[
σ
z
]
η
=
1.5
η
=
4
η
=
1.5
η
=
4
η
=
1.5
η
=
4
Gaussian
−5.60
−13.59
−5.67
−9.29
−7.08
−19.42
Proposed method
−1.04
2.72
−1.72
6.50
−4.21
−0.31Table 3
Relative errors of the maximum response of SDOF system withT
=
4.0 s.
Method
Error in:max
⁡
[
σ
x
]
Error in:max
⁡
[
σ
x
˙
]
Error in:max
⁡
[
σ
z
]
η
=
1.5
η
=
4
η
=
1.5
η
=
4
η
=
1.5
η
=
4
Gaussian
−4.99
−13.02
−1.99
−3.83
−4.92
−17.95
Proposed method
−3.92
−8.22
−1.21
0.86
−4.07
−1.48


## 8. Extension to Multiple-Degree-of-Freedom Systems

In this section an extension of the proposed method to multiple-degree-of-freedom (MDOF) systems is presented.Consider a nonlinear dynamic MDOF structural system subjected to a vector inputP
(
t
) with an m-dimensional response q
(
t
) satisfying the nonlinear differential equation:
(20)
Φ
q
(
t
)
,
q
˙
(
t
)
,
q
¨
(
t
)
=
P
t
.
Here bold letters are related to MDOF systems. The equivalent linear system of (20) is defined as
(21)
M
e
q
¨
t
+
C
e
q
˙
t
+
K
e
q
t
=
P
t
,
where M
e, C
e, and K
e are defined as the equivalent mass, viscous damping, and stiffness matrices of m
×
m. They are determined minimizing the mean square error e
=
E
ε
ε
T, where ε is the m-dimensional vector:
(22)
ε
=
Φ
q
,
q
˙
,
q
¨
-
M
e
q
¨
+
C
e
q
˙
+
K
e
q
.
Denoting by Q and H
e the 3
m-dimensional response vector and the 3
m
×
m equivalent structural matrix given, respectively, by
(23)
Q
=
q
¨
,
q
˙
,
q
T
;
H
e
=
M
e
,
C
e
,
K
e
T
it can be shown that H
e can be obtained from [5, 28, 29]
(24)
H
e
=
E
Q
Q
T
-
1
E
Q
Φ
T
Q
.
Herein, the random excitation P
(
t
) is modelled as a linear combination of the responses of linear filters to white or shot noise. With the aim of calculating the statistical second-order response of the MDOF structure, it is convenient to express (21) in state space form:
(25)
d
Z
t
d
t
=
A
e
Z
t
+
B
Z
F
t
,
where Z
t
=
q
t
,
q
˙
t
T is the 2
m
×
1 state vector and A
e the time-dependent system matrix given by
(26)
A
e
t
=
0
I
-
M
e
-
1
K
e
-
M
e
-
1
C
e
and the matrix B contains the coefficients of the linear combination of the response Z
F of the filters. The dynamics of the latter is governed by
(27)
d
Z
F
t
d
t
=
D
Z
F
t
+
W
′
t
,
where the matrix D contains the coefficients of the filters and the vector W
'
(
t
) the white or shot noise. Equations (25) and (27) can be rewritten as
(28)
d
u
d
t
=
H
t
u
+
W
,
where
(29)
u
t
=
Z
t
Z
F
t
,
(30)
H
t
=
A
B
0
D
.
Matrix W has the same structure as W
' but with higher dimension. Thus, (28) controls the response of an augmented system consisting of the linear filters and the structural system in series [5]. Finally, the covariance matrix of the response of the MDOF system is calculated by solving
(31)
d
Σ
(
t
)
d
t
=
H
t
Σ
t
+
Σ
t
H
T
t
+
2
π
S
F
,
where S
F is a matrix of zeros, except the last element which is equal to the amplitude of the white noise spectral density (
s
0
).The mathematical expressions given in this section correspond to a MDOF structural system withm-dimensional response, while those given in the previous sections correspond to the special case of a SDOF with 2-dimensional response (
q
(
t
)
=
x
,
z
T
). For example (20) took the particular form of (1) and (2); (21) took the particular form of (14); (24) was reduced to (15) while (30) took the particular form of (18).


## 9. Conclusions

A model appropriately representing the time evolution of the joint PDF of the structural response of softening systems subjected to seismic ground motions was proposed. The displacement and velocity was assumed jointly Gaussian and the marginal PDF of the hysteretic component of the displacement was modeled by a mixed PDF which is Gaussian when the structural behavior is linear and turns into a bimodal PDF when the structural behavior is hysteretic.The NSEL method proposed in this study uses mathematical expressions for the equivalent linearization coefficients (16) corresponding to the non-Gaussian response of inelastic nonlinear systems. The mathematical expressions depend on the parameters of the joint PDF of the structural response, and they take into account the time evolutionary character of the PDF. It is noticed that other NSEL methods found in the literature have proposed mathematical expressions for the equivalent linearization coefficients; however, they do not consider the time evolution of the joint PDF.The mathematical expressions derived here (16) can be applied to obtain the non-Gaussian response of any kind of SDOF softening system with smooth transition from elastic to plastic behavior (n
=
1). It is possible to derive expressions of the linearization coefficients for other values of n.The NSEL method proposed was applied to calculate the statistical response of SDOF systems with different vibration periods and different design ductility factors. The results show that, for the SDOF systems analyzed, the proposed criterion estimates with good accuracy the time evolution of the structural response; however, it does not predict the permanent drift exhibited by hysteretic systems with high design ductility values.It is concluded that the proposed NSEL approach is useful to estimate the maximum standard deviation of the structural response of softening systems subjected to seismic ground motions. The relative errors in the calculation of the maximum standard deviation of the displacement, velocity, and hysteretic component of displacement as compared with those obtained by means of Monte Carlo technique are low and in general they are smaller than those obtained by the Gaussian equivalent linearization method.The parametersa
i, b
i, c
i, and d
i, i
=
1,2
,
3,4, in (13) can be found for other types of seismic excitation (e.g., wide-band ground motions).



---

*Source: 539738-2014-11-25.xml*