Non-stationary stochastic vector processes: seismic ground motion applications

f

•

r'

.'All PII:S0266-8920(96)00007-0

ELSEVIER

Probabilistic Engineering Mechanics 11 (1996) 149--168 Copyright © 1996 Elsevier Science Limited Printed in Great Britain. All rights reserved 0266-8920/96/$15.00

Non-stationary stochastic vector processes: seismic ground motion applications George Deodatis Department of Civil Engineering and Operations Research, Princeton University, Princeton, NJ 08544, USA A spectral-representation-based simulation algorithm is used in this paper to generate sample functions of a non-stationary, multi-variate stochastic process with evolutionary power, according to its prescribed non-stationary cross-spectral density matrix. If the components of the vector process correspond to different locations in space, then the process is also non-homogeneous in space (in addition to being non-stationary in time). The ensemble cross-correlation matrix of the generated sample functions is identical to the corresponding target. For the important application of earthquake ground motion simulation, an iterative scheme is introduced to generate seismic ground motion time histories at several locations on the ground surface that are compatible with prescribed response spectra, correlated according to a given coherence function, include the wave propagation effect, and have a specified duration of strong ground motion. Three examples involving simulation of earthquake ground motion are presented in order to demonstrate the capabilities of the proposed methodologies. In the first two examples, acceleration time histories at three points on the ground surface are generated according to a prescribed cross-spectral density matrix, while in the third example, the acceleration time histories are generated to be compatible with prescribed response spectra. Copyright © 1996 Elsevier Science Ltd.

reader is referred to Elishakoff, 1 Kozin, 2 Soong and Grigoriu, 3 Grigoriu, 4 and DeodatisS), the spectral representation method is one of the most widely used today. Although the concept of the method existed for some time, 6 it was Shinozuka 7'8 who first applied it for simulation purposes including multi-dimensional, multivariate and non-stationary cases. Yang 9'1° showed that the Fast Fourier Transform (FFT) technique can be used to dramatically improve the computational efficiency of the spectral representation algorithm, and proposed a formula to simulate random envelope processes. Shinozuka n extended the application of the F F T technique to multi-dimensional cases. Recently, Deodatis and Shinozuka 12 further extended the spectral representation method to simulate stochastic waves, Yamazaki and Shinozuka 13 developed an iterative procedure to simulate non-Gaussian stochastic fields, and Grigoriu 14 compared two different spectral representation models. Four recent review papers on the subject of simulation using the spectral representation method were written by Shinozuka 15 and Shinozuka and Deodatis. 16-18 A spectral-representation-based algorithm 735'16'22 is used in the present paper to simulate non-stationary stochastic vector processes with evolutionary power. If the components of the vector process correspond to

INTRODUCTION Several methods are currently available to solve a large number of problems in mechanics involving uncertain quantities described by stochastic processes, fields or waves. At this time, however, Monte Carlo simulation appears to be the only universal method that can provide accurate solutions for certain problems in stochastic mechanics involving non-linearity, system stochasticity, stochastic stability, parametric excitations, large variations of uncertain parameters, etc., and that can assess the accuracy o f other approximate methods such as perturbation, statistical linearization, closure techniques, stochastic averaging, etc. One of the most important parts of the Monte Carlo simulation methodology is the generation of sample functions of the stochastic processes, fields or waves involved in the problem. The generated sample functions must accurately describe the probabilistic characteristics of the corresponding stochastic processes, fields or waves that may be either stationary or nonstationary, homogeneous or non-homogeneous, onedimensional or multi-dimensional, uni-variate or multivariate, Gaussian or non-Gaussian. Among the various methods that have been developed to generate such sample functions (for a review o f these methods, the 149

150

G. Deodatis

different locations in space (as in the seismic ground motion examples provided in this paper), then the process can also be non-homogeneous in space (in addition to being non-stationary in time). The simulation algorithm is simple and straightforward and generates sample functions of the vector process according to a prescribed nonstationary cross-spectral density matrix. For the important application of earthquake ground motion simulation, acceleration, velocity, or displacement time histories can be generated at several locations on the ground surface according to a target crossspectral density matrix. The only drawback in this approach is that it is often preferable to work with time histories that are compatible with prescribed response spectra, rather than with prescribed power spectral density functions (cross-spectral density matrix). In order to address this issue, the present paper introduces an iterative scheme, based on the proposed simulation algorithm, to generate seismic ground motion time histories at several points on the ground surface that are compatible with prescribed response spectra, correlated according to a given coherence function, include the wave propagation effect, and have a specified duration of strong ground motion. Finally, from the rich bibliography related to the various methods currently available to generate sample functions of stochastic processes, fields and waves, the following representative publications dealing with simulation of multi-variate stochastic processes are mentioned here: Shinozuka and Jan 7 (simulation of multivariate processes by spectral representation), Gersch and Yonemoto 19 (multi-variate ARMA model), Mignolet and Spanos 2°'21 (stability and invertibility aspects of AR to A R M A procedures for multi-variate random processes), Li and Kareem 22 (simulation of non-stationary vector processes using an FFT-based approach), Li and Kareem 23 (simulation of multi-variate processes using a hybrid discrete Fourier transform and digital filtering approach), Ramadan and Novak 24 (asymptotic and approximate spectral techniques to simulate multidimensional and multi-variate processes and fields), and Deodatis 25 (simulation of ergodic multi-variate stochastic processes). From this list, the Li and Kareem 22 and Ramadan and Novak 24 papers considered non-stationary vector processes. At this juncture it should be mentioned that a very promising development with potential applications to simulation of non-stationary multivariate stochastic processes is the use of wavelets (e.g. Gurley and Kareem 26 and Spanos and Zeldin 27). In the remaining part of this paper, the theory, the simulation algorithm, and the iterative scheme are presented for tri-variate non-stationary stochastic vector processes (dimension of vector process is three). This is done for the sake of simplicity in the notation. The simulation algorithm and the iterative scheme can be extended in a straightforward fashion to any other dimension of the vector process.

TRI-VARIATE NON-STATIONARY STOCHASTIC PROCESSES Consider a 1D-3V (one-dimensional, tri-variate) nonstationary stochastic vector process with components fl ° (t), f2° (t) and f o (t), having mean value equal to zero: e[J)°(/)] = 0,

j = 1,2,3

(1)

cross-correlation matrix given by:

R°(t, t + 7-) [R°l(t,t+7 -) R°2(t,t+T) = ]R°l(t,t+r) R°2(t,t+r) [R°l(t,t+r) R°2(t,t+r)

R°3(t,t+~-)] R°a(t,t+r) R°3(t,t+r)

(2)

and cross-spectral density matrix given by: -S°l(W,t) S°(w,t) =

S°2(w,t)

S°3(w,t)] S°3(w,/)

S°l(W,t) S°2(w,t)

(3)

s°,(w,t) Note that because of the non-stationarity of the vector process, the cross-correlation matrix is a function of two time instants: t and t + ~- (t = time and ~-= time lag), while the cross-spectral density matrix is a function of both frequency w and time r Adopting the theory of evolutionary power spectra for non-stationary stochastic processes, 28,29 the elements of the cross-spectral density matrix are defined as:

S°(w,t) = IAj(w,t)12Sj(w),

j = 1,2,3

(4)

l)

Aj(w, t)Ak(w , t) ~/Sj(w) Sk(W) rjk(W), j,k= 1,2,3; j -i k (5) where Aj(w, t), j = 1,2, 3 are the modulating functions off1°(/), f°(t) and f°(t), respectively, Sj(w), j = 1,2, 3 are the (stationary) power spectral density functions of fl°(t), f°(t) and Jr(l), respectively, and £jk(W), j,k = 1,2,3, j - i k are the complex coherence functions between fj°(t) and f°(t). It should be pointed out that eqns (4) and (5) imply that the modulating function Aj(w,t) represents the change in the evolutionary power spectrum, relative to the (stationary) power spectral density function Sj(w). Consequently, for any time instant t, the diagonal elements of the cross-spectral density matrix are real and non-negative functions of w satisfying:

S°(w,t)=S°(-w,t),

j = 1,2,3

and for every t

(6) while the off-diagonal elements are generally complex functions of w satisfying:

S°k(w,t) = Sfk*(-w,t),

j,k=

1,2,3,

j -i k

and for every t

Sfk(w,t) = S°*i(w,t), and for every t

(7)

j,k=

1,2,3;

j-ilk (8)

Non-stationary stochastic vector processes where the asterisk denotes the complex conjugate. Equation (8) indicates that the cross-spectral density matrix S°(w, t) is Hermitian for any value of t. The elements of the cross-correlation matrix are related to the corresponding elements of the cross-spectral density matrix through the following transformations:

R (t, t + 7-) = J°°oo aj(w, t)A;(w, t + 7-) e' rs;(w) dw, j=

1,2,3

(9)

SIMULATION FORMULA

In the following, distinction will be made between the non-stationary stochastic vector processfj°(t), j = 1,2, 3 and its simulation fj(t), j = 1,2, 3. In order to simulate the 1D-3V non-stationary stochastic process)5° (t), j = 1,2, 3, its cross-spectral density matrix S°(w, t) must be decomposed at every time instant t under consideration, into the following product: S°(w, t) = n(w, t)nT*(w, t) for every t under consideration

ROk(t,t+7-)=lo°_ooAj(w,t)Ak(w,t+r ) x eicorV/Sj(oo) Sk(w)rjk(w) dw, j,k=1,2,3,

151

j ¢ k

(10)

(17)

where superscript T denotes the transpose of a matrix. This decomposition can be performed using Cholesky's method, in which case H(w, t) is a lower triangular matrix: o

n(w,t) = [a2,(w,t)

Special case: uniformly modulated non-stationary stochastic vector process

[ a31(w, t)

For the special case of a uniformly modulated nonstationary stochastic vector process, the modulating functions Aj(w,t), j = 1,2,3 are independent of the frequency w:

aj(w,t) =Aj(t),

j=

1,2,3

(11)

e

1,2,3

(12)

(13)

In such a case, the three components of the non-

stationary stochastic vector process are expressed as: j = 1,2,3

(14)

where g°(t), j = 1,2, 3 are the three components of a stationary stochastic vector process, having mean value equal to zero: ~[g°(t)] = 0 ,

j=

1,2,3

(15)

and cross-spectral density matrix given by: SI(CO)

S°(c°) = ~ P 2 1

~

(co)

1~,2(¢,.0) SV/~O)) S3 (~.o)F 13(co) ]

S: (~o)

x/S3(W)Sl(W)I~31(co) ~ P 3 2

Hjj(w,t) = Hjj(-w,t),

j=

1,2,3

~F23(w)[ (¢,,0)

$3(co)

(19)

njk(W , t) = I/-/jk(W, t) ]e i°jk(co'tl, k = 1,2;

j > k

(20)

where

X J?~ 4Sj(W)Sk(w)Pjk(w)eicordW,

fj°(t) = Aj(t)g°(t),

whose diagonal elements are real and non-negative functions of w and whose off-diagonal elements are generally complex functions of w. The following relation is satisfied by the diagonal elements of H(w, t):

j = 2, 3;

R°k(t, t+ r) = Aj(t)Ak(t + 7-)

j ¢ k

(18)

If the off-diagonal elements Hjk(w, t) are written in polar form as:

Ro(t,t+r)=Aj(t)Aj(t+7-)l~_ooSj(w) icordw,

j,k=1,2,3,

n32(w, t)

0]

0 /-/33 (w, t)

and for every t

and eqns (9) and (10) reduce to:

j=

n22(w,t )

/ J

(16) Note that the elements of matrix S°(w) shown in eqn (16) consist of terms that have been defined in eqns (4) and (5).

(Im[Hjk (w, t)] "~ 0jk(w, t) = tan q \ Re[/-/a-k(~, t)] J

(21)

with Im and Re denoting the imaginary and the real part of a complex number, respectively, then the following relations are satisfied: ]//jk(W, t) l = I Hjk(--w, t)[ , j = 2, 3; k = 1,2, j > k and for every t

Ojk(W, t) = --Ojk(--W, t), k = 1,2,

j > k

(22)

j = 2, 3; and for every t

(23)

Once matrix S°(w, t) is decomposed according to eqns (17)-(18), the non-stationary stochastic vector process fj°(t),j = 1,2, 3 can be simulated by the following series a s N ~ oo 3 N

~'(t) = 2 Z ~_~lI-Ijm(Wl,t) lv F ~ m=l l=1

X COS[~t t -- 0jm(Wt, t) + ~'ml],

j=

1,2,3 (24)

G. Deodatis

152

l= 1,2,...,N:

or explicitly N

3

fl(t) = 2 Z I H l l ( W l ,

t) IX/Aw

N

fj(i)(t) = 2 E

1=1

E

m=l

x cos[wit - 011 (wh t) + Oll]

(25a)

II-I]m(wl't) lv/-~-~

1=1

× COS[WIt -- Ojm(W h t) nt- ,'h(i)l WmlJ,

j = 1 2, 3

N

f2( t) =

(29)

2 ~', l H21(wl, t) l ~ I=1

The generated values offj (i) (t) according to eqn (29) are bounded as follows:

x cos[w/t - 021 (Wh t) + ~1l] N

3

+ 2 y ~ I H22(Wl, t) l ~ I=1

j = 1,2,3

(30)

(25b)

N f 3 ( t ) = 2 Z I n31(Wl, t) l X / ~ W l=1

For a given form of the cross-spectral density matrix, the above bound is large enough for all practical applications, even for relatively small values of N. 2s It is obvious that this bound can be easily calculated for any form of the cross-spectral density matrix to be used. It will be shown now that the ensemble expected value e[J)(t)], j = 1,2, 3 and the ensemble auto-/crosscorrelation function Rjk(t, t + 7"), j, k = 1,2, 3 of the simulated non-stationary stochastic vector process ~(t) are identical to the corresponding targets, e[j~°(t)] = 0, j = 1,2, 3 and R°k(t, t + 7"), j, k = 1,2, 3, respectively.

X COS[WIt -- 031(Wl, t) + (I)ll ]

N

I H32 (Wh t) l

1=1

x cos[wlt - 032(Wl, t) + ~21] N

+ 2Z

t) lx/A~w,

rn=l l=1

x cos[wit - 022(Wl, t) + ~21]

+ 2~

N

fj(i)(t) < 2 Z E l I - I j m ( W l ,

I H33 (wt, t) I

1=1

x cos[wit - 033 (Wl, t) + O311

(25c)

e[j~(t)] = e[fj°(t)] = 0,

j = 1,2, 3

Proof: This proof is trivial and similar to the corre-

where

w l = lAw,

(i) Show that:

sponding one for 1D-1V processes. 17

l = 1,2,...,N

(26)

Aw = w___u_u N

(27)

Ojm(Wt, t) = tan -1 ( Im[Hjm(wl'/)])

(28)

\ Re[/-/jm (wl, t)]

In eqn (27), wu represents an upper cut-off frequency beyond which the elements of the cross-spectral density matrix (eqn (3)) may be assumed to be zero for any time instant t. As such, wu is a fixed value and hence Aw ~ 0 as N ~ c~ so that N A w = Wu. A criterion to estimate the value of Wu can be found in Shinozuka and Deodatis. 17 The e~lt, ff2t, /b31, I = 1 , 2 , . . . , N appearing in eqns (24)-(25) are three sequences of independent random phase angles distributed uniformly over the interval [0, 2~]. It should be noted that the simulated non-stationary stochastic vector process fj(t), j = 1,2,3 is asymptotically Gaussian as N ~ c~ because of the central limit theorem. 17 A sample function fj(i)(t), j = 1, 2, 3 of the simulated non-stationary stochastic vector process fj(t), j = 1,2, 3 can be obtained by replacing the three sequences of random phase angles tI~ll , (T~21, ~31, 1._ 1,2 . . . . , N (,) with their respective i-th realizations 4hi, "P'21 .~(i)~ .~(i) W31 ,

Rjk(t , t + 7") = R°k(t, t + 7"),

(ii) Show that:

j,k=

1,2,3

Proof: Utilizing the property that the operations of mathematical expectation and summation are commutative, the ensemble auto-/cross-correlation function of the simulated non-stationary stochastic vector process fj(t), j = 1,2, 3 can be written as: Rjk(t , t + 7") = e[fj(t)fk(t + 7")] 3

=4 E

3

Z

N

N

EEIHjmt(WI'

't)]

ml=l m2=1 I1=1 /2=1

x IHk,,2(wt2,t+w)lAw X /3{COS[Wll / -- Ojmt(Wll , t) "}- t~mall ] X COS[WI2 (t -'1-7") -- Okm2(Wl2, t + 7") 3 3 N N

"q--(I)m212]}= 2 E E E Z [/'/j'ml(W/I, t) [ ml=l m2=l I1=1 12=1

x I Hkm2(Wl2, t + 7") I AW × e{COS[(Wh + Wl=)t + Wl27" -- Ojm I (Wl, , t) -- Okrna(Wl2, t Jr" 7") "q- ~rntl l

+ ~rn,l,] + COS[(Wl, -- Wl, )t + W127"+ Ojm,(Wl,, t) -

Okm2(Wl2, t + 7") -- ~b,n,i~ + ~m212]}

(31)

Non-stationary stochastic vector processes Since the 4~'s are independent random variables distributed uniformly over the interval [0,27r], the expected value appearing in eqn (31) is different from zero only when: m I = m2 = m

and

I l = l2 = l

(32)

Under the conditions of eqn (32), eqn (31) can be written as: 3

Ryk(t , t + 7.) = 2 E

N

Z IHi m(O)t' t) l

m=l l=1

x [Hkm(Wl, t+W) lAo) × COS[O)/7. -1- Ojrn(O)l, l) -- Okra(O) h t -']- 7-)]

(33) Taking into account eqns (19), (22) and (23), the expression for Rjk(t, t + 7) in eqn (33) can be written in the limit as AO) ---. 0 and N ~ o0 (while keeping in mind that Wu = NAO) is constant and that the elements of the cross-spectral density matrix S°(O), t) are zero for [O)] > O)u at any time instant t) in the following way:

Rjk(t't + 7-) = J ~

3my~/lHJ = l m(O)'t)l'Hk'(O)'t +

X e i[wr+Oy"(w' t)--Okm(W' t+7")] d w

(34)

Using now the polar form representation of the elements of the H(o), t) matrix (eqn (20)), eqn (34) is written as:

Rjk(t,t + 7-) =

Z Hjm(O),t)Hi*m(O),t + r)ei~" dO) -oo m=l

(35) Then, using the decomposition shown in eqn (17) in conjunction with eqns (4) and (5), it is straightforward to show that:

Rjj(t,t + 7-) =

Aj(O),t)Aj(O),t + 7-)Sy(O))e'~r dO),

j = 1,2, 3

(36)

Rjk(t,t + 7") = I f

Aj(w,t)Ak(W,t + 7-)

x v/Sy(O)) Sk(o)) l"jk(w) e i~r dw,

j,k=1,2,3,

j ~ k

(37)

Comparing finally eqns (9) and (10) to eqns (36) and (37), it is evident that:

Rjk(t, t + 7-) = ROk(t, t + 7-),

j, k = 1,2, 3

(38)

Special case: uniformly modulated non-stationary stochastic vector process For the special case of a uniformly modulated nonstationary stochastic vector process, simulation can be

153

performed on the basis of eqn (14), instead of using eqn (24). The simulation formula corresponding to eqn (14) is:

fj.(t) = Aj(t)gy(t),

j = 1,2, 3

(39)

where gy(t) is the simulation of the stationary stochastic vector process g°(t) having mean value equal to zero and cross-spectral density matrix shown in eqn (16). It should be mentioned that the simulation of stationary and ergodic stochastic vector processes can be performed with great computational efficiency using the Fast Fourier Transform (FFT) technique, as described by Deodatis. 25 Comments on computational efficiency At this juncture, it should be pointed out that it is not possible to take advantage of the F F T technique when using the (non-stationary) simulation formula shown in eqn (24), in contrast to the corresponding formula for simulation of stationary stochastic vector processes. 25 This is due to the fact that the coefficients [nym(o)l, t) [ V / ~ in the double summation of eqn (24) are now functions of both frequency and time. However, this should not be of any great concern computationaUy, since in most cases of practical interest the nonstationary stochastic vector process J)°(t), j = 1,2, 3 is limited to relatively short durations by the modulating functions Aj(w, t), j = 1,2, 3 (e.g. ground motion acceleration time histories). A case of non-stationary stochastic vector processes where the F F T technique can be used in the simulation formula is that of uniformly modulated processes (eqn (39)). At this juncture, it should be noted that Li and Kareem 22 have proposed a methodology to simulate non-stationary vector processes taking advantage of FFT.

S I M U L A T I O N OF SEISMIC G R O U N D M O T I O N COMPATIBLE WITH PRESCRIBED RESPONSE SPECTRA The simulation algorithm presented earlier in this paper can generate sample functions of a general nonstationary stochastic vector process with evolutionary power, according to a prescribed non-stationary crossspectral density matrix. A very important application is simulation of earthquake ground motion time histories. In such a case, acceleration, velocity, or displacement time histories can be generated at several locations on the ground surface according to a target cross-spectral density matrix. This will be demonstrated in the first two numerical examples in the following section. The only drawback in this approach is that it is often preferable to work with time histories that are compatible with prescribed response spectra, rather than with prescribed power spectral density functions (cross-spectral density

154

G. Deodatis

matrix). An obvious reason for this preference is that it is much easier and more reliable to find a response spectrum specified for a given location, rather than a power spectral density function. For example, design codes usually provide response spectra as a function of local site conditions. In order to address this issue, a methodology is proposed now to simulate seismic ground motion compatible with the following three prescribed quantifies: (i) response spectra, (ii) complex coherence functions, and (iii) modulating functions. Although the methodology is presented in the following for accelerations, it is straightforward to modify it to accommodate velocities or displacements. The choice of accelerations is done solely for demonstration purposes. According to the proposed methodology, the acceleration time histories at three points on the ground surface are considered to be a uniformly modulated, trivariate, non-stationary stochastic vector process. In general, the three points correspond to different local soil conditions. Consequently, a different target acceleration response spectrum RSAj(w), j = 1,2, 3 is assigned to each one of the three points. In addition, complex coherence functions Fjk(w), j, k = 1,2, 3, j # k are prescribed between pairs of points, and modulating functions A j ( t ) , j = 1,2,3 are assigned at each point. The simulation of the acceleration time histories is then performed according to the iterative scheme shown in Table 1. This scheme is not expected to perfectly converge at all frequencies as the number of iterations increases (perfect convergence at frequency w = w* is expressed as: R S A j ( w * ) / R S A ( f J ) ( w *) = 1, j = 1,2, 3, where RSAU))(w) is defined in Table 1). This is why no convergence criterion is included in the flowchart of Table 1. As will be shown in the third numerical example in the following section, only a small number of iterations is usually needed (in most cases less than ten) for a sufficiently accurate convergence at every frequency. At this point it should be mentioned that the idea for upgrading the individual power spectral density functions depicted in Table 1 has been suggested by Gasparini and Vanmarcke 3° for one-dimensional and uni-variate stochastic processes. Two alternative methodologies to generate ground motion time histories compatible with prescribed response spectra, coherence function, velocity of wave propagation, and duration of strong ground motion, have been suggested by Hao et al. 31 and by Abrahamson. 32 The Hao et al. methodology and the methodology proposed in this paper in Table 1 are conceptually similar, while the Abrahamson approach is a conceptually different time-domain-based methodology. According to Hao et al., 31 a spectral-representationbased simulation algorithm 7 is first used to generate stationary time histories that are compatible with a coherence function and a velocity of wave propagation,

Table 1. Iterative scheme to simulate response spectrum

compatible acceleration time histories at three points on the ground surface

Read input data: t> target acceleration response spectra: RSAj(w); j = 1,2,3 complex cohereaace functions: Fj/,(w); j , k = 1,2,3; j ~ k modulating functions: Ai(t)! j = 1, 2, 3 initialize power spectral density functions Sl(t0), S2(to), S3(w ) by setting them equal to n (noa-zero) constant value over the entire frequency range .___~ rI

I

I

Generate gl(t), g2(t), 93(4) as a stationary, tri-variate, stochastic vector process with cross-spectral density matrix shown in Eq. (16). Refer to Eq. (39). Compute acceleration time histories as: fl(t) = At(t)gl(t), f2(t) = A2(t)g2(t), and fa(t) = A3(t)g3(t). Refer to Eq. (39). Compute acceleration response spectra RSA(ID(w), RSA(I*)(w), RSA(Ia)(w), corresponding to fl(t), f2(t), fs(t), respectively.

Yes

.-®

Upgrade power spectral density functions as: Sl(,~) ~ Sl(~)

[ RSAI(w) ]z [RSA(t,)(~)J

s" " f RSA,(~) 1~

[ Rsa~(~) ]' &(w) ~ S~(.~) LRSA(I~)(w) J

but not with the prescribed response spectrum. After multiplying these stationary time histories by an appropriate envelope function to introduce non-stationarity, Hao et al. are then adjusting each non-stationary time history independently to make it compatible with the prescribed response spectrum. This adjustment is performed by Fourier transforming each non-stationary time history to the frequency domain, multiplying its frequency domain Fourier transform by the ratio of the prescribed response spectrum over the computed response spectrum of the non-stationary time history, and then inverse transforming the product back to the time domain. This adjustment procedure is usually repeated a few times. On the other hand, the methodology proposed in this paper in Table 1 is starting by using a different spectral-representation-based simulation algorithm 25 to generate ergodic, stationary time histories (the Shinozuka and Jan 7 algorithm is generating time histories that are not ergodic) that are again compatible with a coherence function and a velocity of wave propagation, but not with the prescribed response spectrum. After multiplying these stationary time histories by an appropriate envelope function to introduce non-stationarity, the methodology proposed in this paper upgrades the power spectral density functions of the components of

Non-stationary stochastic vector processes the vector process as indicated in Table 1, generates new stationary time histories according to the upgraded cross-spectral density matrix, and multiplies them again by the envelope function to introduce non-stationarity. This upgrading procedure is repeated a few times to make the time histories response spectrum compatible to the desired degree. NUMERICAL EXAMPLES In order to demonstrate the capabilities of the proposed algorithms to simulate non-stationary stochastic vector processes, three examples involving simulation of earthquake ground motion are selected. In the first example, ground motion time histories are modeled as a uniformly modulated non-stationary stochastic vector process, and sample functions are generated according to a target cross-spectral density matrix. In the second example, ground motion time histories are modeled as a non-stationary stochastic vector process with amplitude and frequency modulation, and sample functions are again generated according to a target cross-spectral density matrix. Finally, in the third example, ground motion time histories are modeled as a uniformly modulated non-stationary stochastic vector process, but unlike the first example, sample functions are generated to be compatible with prescribed response spectra.

Example 1 In this numerical example, ground motion time histories are modeled as a uniformly modulated non-stationary stochastic vector process, and sample functions will be generated according to a target cross-spectral density matrix. Specifically, the acceleration time histories at three points on the ground surface along the line of main wave propagation are considered to be a tri-variate nonstationary stochastic vector process. The configuration of the three points is shown in Fig. 1 where the arrow indicates the direction of wave propagation. The simulated earthquake ground motion at points 1, 2 and 3 will have the following characteristics: (i) Points 1, 2 and 3 correspond to different local soil conditions. Specifically, point 1 is characterized by stiff soil conditions, point 2 by intermediate soil conditions and point 3 by soft soil conditions. This is a very desirable capability of the method, being able to simulate ground motion at neighboring points on the ground surface with different local soil conditions and consequently different frequency contents. Such a case is encountered, for example, at the supports of intermediate to long-span bridges located in areas with abrupt changes in soil conditions along the axis of the bridge.

t~

50 m

155

--~t-- --50 m

•

g,

O

O

-_~

O

Fig. 1. Configuration of points 1, 2 and 3 on the ground surface along the line of main wave propagation, for Example 1. Point 1 corresponds to rock or stiff soil conditions, point 2 corresponds to deep cohesionless soils and point 3 corresponds to soft to medium clays and sands. (ii) In addition to different frequency constants at points 1, 2 and 3, the acceleration time histories at these three points will also be correlated according to a specified coherence function and they will reflect the wave propagation effect according to a specified velocity of wave propagation. The method is therefore capable to simulate ground motion time histories that, at the same time, are spatially correlated, include the wave propagation effect and correspond to different local soil conditions. (iii) Finally, the acceleration time histories at points 1, 2 and 3 will reflect the non-stationary characteristics of ground motion according to specified modulating functions Aj(~o, t), j = 1,2, 3. The three components of the tri-variate non-stationary stochastic vector process describing the acceleration time histories at points 1, 2 and 3 (see Fig. 1) are denoted by J]°(t), f°(t) and f°(t), respectively. The mean value of the process is equal to zero (see eqn (1)), while the elements of its cross-spectral density matrix (see eqn (3)) are defined as follows:

S°(w, t) = I Aj(w, t)[2Sj(~), t)

= Aj( , xexp[-i~],

j = 1,2, 3

(40)

t) v/Sj( ) j,k=l,2,3,

j Ck (41)

where ")'jk(w) are the (stationary) coherence functions between fj°(t) and f°(t), and exp[-i(w~jk/v)] is the wave propagation term, with ~jk being the distance between points j and k, and v being the velocity of wave propagation. At this juncture, it should be pointed out that the expressions shown in eqns (40) and (41) describe earthquake ground motion that is non-homogeneous in space (since $1(~) ¢ $2(~) ¢ $3(~)) and non-stationary in time. The Clough-Penzien acceleration spectrum 33 is selected to model the (stationary) power spectral density functions Sj(w); j = 1,2, 3 of the acceleration time histories ~.°(t);

G. Deodatis

156 j = 1,2, 3, respectively: 1 + 4 ~gj [~-~gj]

Sj(w)=Soj

1 - [~_]2~+4ffg2j I. gj

a) 2

)

x[/ u/ 03

2 2

2

_

~

2

1.1

'

+ 4 ~ j ~-~fj. 1

j = 1,2,3

(42)

where Soj is a constant determining the intensity of acceleration at point j, OJgj and (gj c a n be thought of as some characteristic frequency and damping ratio of the ground at point j, and wfj and (fj a r e filtering parameters for point j. The Harichandran-Vanmarcke model 35 is chosen to describe the (stationary) coherence functions 7jk(~O), j, k = 1,2, 3, j ~ k between fj°(t) and f°(/):

where al and a 2 are model parameters depending on such factors as earthquake magnitude and epicentral distance. It should be noted that eqn (45) includes the wave propagation effect. At this juncture, it should he pointed out that the models used for Sj(•), 7jk(W) and Aj(w, t) in eqns (42), (43) and (45), respectively, were selected for demonstration purposes only. There are several other models in the literature that can be used for Sj(~), 7jk(W) and Aj(w, t). The next step is to select numerical values for the parameters appearing in eqns (42)-(45). Starting from %j and (gj, j = 1,2, 3 appearing in eqn (42), the values suggested by Ellingwood and Batts 36 for three different soil conditions are used in this study:

Point 1: Rock or stiff soil conditions: Wgl = 87rrad/s,

(gl = 0.60

Point 2: Deep cohesionless soils: Wg2 = 57r rad/s,

(g2 = 0.60

+ (1 - A) exp [ - 2~jk (1 - A + aA)]

L

j , k = 1,2,3,

j ~ k

(43)

where ~jk is the distance between points j and k, O(a)) is the frequency-dependent correlation distance:

j

= k 1 + \w0:

(44)

and A, a, k, ~0 and b are model parameters. Since the acceleration time historiesfj°(t),j = 1,2, 3 are modeled as a uniformly modulated non-stationary stochastic vector process, the modulating functions Aj(~, t), j = 1,2, 3 will be functions of time only. The BogdanoffGoldberg-Bernard model 34 is used for this purpose:

a l ( ~ , t) = .41(t) = altexp(-a2t)

for t > 0 (45a)

A2(w , t) = A2(t ) a

t --

0

for t

>

(2__~1

V

(45b) A3(a;, t) = A3(t) /al

t-

~fl = 0.87rrad/s,

(fl = 0.60

(47a)

Point 2: Deep cohesionless soils: wf2 = 0.57rrad/s, (f2 = 0-60

(47b)

Point 3: Soft to medium clays and sands: ~;f3 = 0.247rrad/s,

0

t -~ 1

fort>~31 _

for 0 < t<

(45c)

(47c)

S02 = 99.7 cm2/s 3,

184.5cm2/s 3

(48)

For the various parameters appearing in eqns (43) and (44), the values estimated by Harichandran and Wang 38 by analyzing data from the SMART-1 seismograph array will be used in this study for demonstration purposes: a = 0"022,

w0 = 12.692 rad/s, ~3--~1 v

(f3 = 0"85

Finally, the last parameter appearing in eqn (42), Soj, j = 1,2, 3, is computed so that the standard deviation of the Kanai-Tajimi part of the (stationary) power spectral density function is equal to 100 cm/s 2 for all three points 1, 2 and 3:

A = 0.626, exp -a2

(46c)

Point 1: Rock or stiff soil conditions:

803 =

/ < ~2..~1

for0<

(g3 = 0.85

The filtering parameter a;fj (eqn (42)) is set equal to 0.10 of the corresponding Wgj value, while the other filtering parameter (fj (eqn. (42)) is set equal to the corresponding (gj value, following the recommendation by Hindy and Novak: 37

S01 = 62.3 cm2/s 3, exp --a2 t --

(46b)

Point 3: Soft to medium clays and sands: O)g3 = 2.47rrad/s,

7 j k ( w ) = A e x p [ - - ~2(jk (1 -- A +etA)]

(46a)

k = 19700m,

b = 3.47

(49)

Parameters a~ and az appearing in the expressions for the modulating functions (eqn (45)) are set equal to: al = 0.906

and

a2 = 1/3

(50)

Non-stationary stochastic vector processes so that the maximum value of A/(t),j = 1,2, 3 occurs at time t = (3 +~/l/V)s, j = 1,2,3 and is equal to unity (from Fig. 1 it is obvious that ~l] = 0, ~2] = 50m, ~31 = 100m). Finally, the velocity of wave propagation v is set equal to: 38 v = 1000m/s

wu = 126rad/s (20Hz)

and

N = 126

(52)

The simulation is performed at 3072 time instants, with a time step At = 6.14 × 10 -3 S, over a length equal to 3072 × 6.14 × 10 -3 = 18.85s. One generated sample function for the acceleration at points 1, 2 and 3 (see Fig. 1), denoted by fl (t), f2 (t) and f3 (t), respectively, is displayed in Fig. 5. The most obvious characteristic of the acceleration time histories in Fig. 5 is their non-stationarity which is described by the modulating functions shown in Fig. 4. The different frequency contents o f f ] (t), f2(t) and f3 (t) (making the acceleration non-homogeneous in space

~

0

¢5-

~

(51)

Based on the numerical values given in eqns (46)-(51), the (stationary) power spectral density functions Sj(w), j = 1,2, 3 are plotted in Fig. 2, the (stationary) coherence functions "y/k(w),j, k = 1,2, 3, j # k are plotted in Fig. 3 and the modulating functions Aj(t), j = 1,2, 3 are plotted in Fig. 4. Since the non-stationary stochastic vector process in this example is a uniformly modulated one, the generation of its sample functions will be performed using eqn (39). It is reminded that the simulation of stationary stochastic vector processes involved in eqn (39) is performed with great computational efficiency using the F F T technique, as described by Deodatis. 25 The upper cut-off frequency Wu and the value of N (see eqn (27)) are set equal to:

157

d"

~(O;)")'3N"~.2 3 1 k

((zJ)

t::;dca

Frequencyw(rad/sec) Fig. 3. Coherence functions ~21(~), ")'31(~) and ~32(~), for Example 1. and described by the (stationary) power spectral density functions shown in Fig. 2) can also be detected in Fig. 5. Specifically, the characteristic frequency of the ground for fl(t) (b3g1 = 87rrad/s) is higher than that for 3~(t) (Wg2 = 57r rad/s), which is then higher than that forj~ (t) (Wg3 = 2.4~r rad/s). This characteristic is not as obvious as the non-stationarity one in Fig. 5, because of the relatively large bandwidths of SI (w), $2 (w) and $3 (w), which are controlled by (el, ~e2 and (g3, respectively. The wave propagation effect and the correlation among fl(t), f2(t) and f3(t) cannot be observed in Fig. 5. For this reason, a segment of the acceleration time histories shown in Fig. 5 (specifically from time t = 3.1 s up to time t = 6.2 s) is magnified and displayed in Fig. 6. In this figure, the wave propagation effect is easily detected by following the movement of peak A in fl (t), f2 (t) and f3 (t), which reflects the v = 1000 m/s velocity of wave propagation. The correlation among Ji(t), f2(t) andf3(t) is also easily detected in Fig. 6, as major peaks in the three time histories (like peak A) retain their

) c:l

o "~-

~,-

i

o-

2'0

,'0

6'o

~,

,~o

Frequencyw(rad/sec)

o

110

Fig. 2. The power spectral density functions at points 1, 2 and 3, for Example 1: Sl(w)~ rock or stiff soil conditions; S2(w)---~deep

cohesionless soils; and S3(w)~soft medium clays and sands.

to

1'5

Time (see) Fig. 4. Modulating functions At(t), A2(t) and Example 1.

.43(t),for

G. Deodatis

158

O

}~ li i~["l l ~ l J

Point I--+ rock or stiff soil conditions

~J

-+,

V,r+'+v,,".",+ I

6

.,,+.+,,+..+,..,-+ ....... ............

I

+

i0

1+

a0

•

E3

•z ,

•

°

o+

~, 0

10

20 Q

0

ta

o

-----'-r

o

"v-----

I+o

(s Time (see)

a0

Fig. 5. Generated sample function for the acceleration at points 1, 2 and 3, over a length equal to 18.85 s (Example 1).

general shape after some loss of coherence described by the (stationary) coherence functions shown in Fig. 3. It is therefore obvious from Figs 5 and 6 that the proposed algorithm is able to simulate non-stationary ground motion time histories that are spatially correlated according to a given coherence function, include the wave propagation effect and are non-homogeneous in space (or equivalently they correspond to different frequency contents). Finally, the ensemble auto-/cross-correlation function Rjk(t , t + "r) is computed from 1000 sample functions at time instant t = 3.14s and plotted as a function of r in Fig. 7 versus the target auto-/cross-correlation function R°k(t, t + "r). As can be seen in Fig. 7, the agreement between R°k(t, t + 7-) and Rjk(t, t + "r) is very good, especially in the vicinity of the dominant peak of the auto-/cross-correlation functions. Some small differences can be observed away from this dominant peak in Fig. 7. However, these differences disappear when the ensemble auto-/cross-correlation function Rjk( t , t + -c) is computed from 100 000 sample functions, as can be seen in Fig. 8 where R?k(t , t + T) and Rjk(t, t + -c) practically coincide. It should be noted that the agreement between R°k(t, t + r) and R/k(t, t + 7-) observed in Figs 7 and 8 for time instant t = 3.14 s is typical for any time instant

from t = 0 up to t = 18.85s (the duration of the time histories shown in Fig. 5). Example 2

Unlike the first example where ground motion time histories were modeled as a uniformly modulated nonstationary stochastic vector process, in this example they are modeled as a non-stationary stochastic vector process with amplitude and frequency modulation. This means that both the amplitude and the frequency content of ground motion change as a function of time (in the first example, it was only the amplitude of ground motion that was varying with time). As in the previous example, sample functions will be generated according to a target cross-spectral density matrix. Figure 9 displays an acceleration record from the 1964 Niigata earthquake involving both amplitude and frequency content variation as a function of time. Specifically, this record shows an abrupt change of its frequency content between approximately 8 and 10 s. During this 2 s period, the frequency content is transformed from one containing a relatively broad band of frequencies to one containing essentially a single low frequency. It is believed that this phenomenon is due to soil liquefaction. The objective of this example will be to reproduce the general

Non-stationary stochastic vector processes

i s

159

A

A --

, I

f~

-

i~.

/V V

Point 1 -* rock or stiff soil conditions

wv

vv v

\

,,q

t:~

3'.5

5:0

415

5:5

6'.0

|! ~

_

°

315 ~-

~i

A

~

I

,~ )

Point 2 --, deep cohesionless soils

415

/| IA

510

515

610

Point 3 --, soft to medium clays and sands

/\

,v

3.5

4.0

4.5

5.0

5.5

6.0

Time(see)

Fig. 6. Generated sample function for the acceleration at points I, 2 and 3 from t = 3.1 s up to t = 6.2 s (Example 1). frequency and amplitude variation characteristics of the acceleration record shown in Fig. 9. The acceleration time histories at three points on the ground surface are again considered to be a tri-variate non-stationary stochastic vector process. The configuration of the three points is shown in Fig. 10, indicating that it is not necessary for the points to be on a straight line, or to be equidistant (compare to previous example). For simplicity, no wave propagation will be considered in this example. The simulated earthquake ground motion at points 1, 2 and 3 will have the following characteristics: (i) Points 1, 2 and 3 correspond to the same local soil conditions. (ii) The frequency content of the acceleration time histories at these three points will change with time as will be described in the following, in order to capture the unique characteristics of the record shown in Fig. 9. (iii) Tbe acceleration time histories at points 1, 2 and 3 will be correlated according to a specified coherence function. (iv) Finally, the amplitude variation as a function of time of the acceleration time histories at these three points will be described by prescribed modulating functions Aj(u3, t), j = 1,2, 3.

The three components of the tri-variate non-stationary stochastic vector process describing the acceleration time histories at points 1, 2 and 3 (see Fig. 10) are denoted by f]°(t),f°(t) andf3°(t), respectively. The mean value of the process is equal to zero (see eqn (1)), while the elements of its cross-spectral density matrix (see eqn (3)) are given by the expressions shown in eqns (40) and (41). It should be pointed out, however, that since no wave propagation is considered in this example, the wave propagation term exp[--i(w~jk/V)] in eqn (41) should be set equal to unity (v ~ c~). The Clough-Penzien acceleration spectrum shown in eqn (42) is selected again to model the power spectral density functions Sj(~), j = 1,2, 3 of the acceleration time historiesfj°(t), j = 1,2, 3, respectively. However, in order to account for the variation of the frequency content with time, the characteristic frequency and damping ratio of the ground are considered now to vary as a function of time as follows:

~gj(t)

=

15.56rad/s f o r 0 _ < t < 4 . 5 s 27'12tp3 - 40"68t2 + 15"56 2.0rad/s

for 4.5 < t < 5.5 s, j = 1 , 2 , 3 for t _> 5.5s

(53)

G. Deodatis

160

R~n ~ continuous line R n ~ dotted line

80

t

R°l --, continuous line

-4

~]°I

I R21 --* dotted line

tO

1

. • ..:...-.'.:.....,....~.'..:..

. ......"

.~..'-...'.. ' . .

01

gq

R~2 --* continuous line R22 --* dotted line

R°I ~ continuous line :. Rsx --* dotted line

O O O tO

O O O tM

g ...

,.

.

.

.

.

.

..

o

.....

0

0

; R~32 -* continuous line --* dotted line

/~ss -+ continuous line

°

I

Rss -* dotted line

O

g. (D

•

'

'

.-"%...,t

,.. ,.._ . . _ . ".."...',

. .. "'....

. ' .. "' . '"-...f",'..

""_-"..,.'.

..

r (see)

r (see)

Fig. 7. Ensemble auto-/cross-correlation functions [R/~(t, t + r)] computed from 1000 sample functions versus the corresponding targets [R°k(t, t + r)]. Both R j k ( t , t + r) and Ryk(t, t + r) are plotted at time instant t = 3.14s as a function of r.

0-64

(gj(t)

for 0 < t < 4.5s

instant: cr 2

=

Soj( t)

1.25tp3 - 1.875tff + 0 . 6 4 for 4.5 < t < 5.5s, j - - 1,2,3 0.015

tp = t -

j=

(fj(t) = (gj(t),

1,2,3 (56)

4.5 s and

~fj(t) = O'lwgj(t),

7rWgj(t)[2~gj(t)+~]'

for t >_ 5.5s (54)

where

=

j = 1,2, 3 (55)

The expressions for Wgj(t) and (g/(t) are taken f r o m Deodatis and Shinozuka 39 and describe a sudden d r o p in the values o f the characteristic frequency and d a m p i n g ratio during the one-second period f r o m t = 4.5 to 5.5 s. The constant determining the intensity o f the acceleration is also going to be a function o f time, so that the standard deviation o f the K a n a i - T a j i m i part o f the spectrum is equal to ~r = 100 cm/s 2 at every time

The H a r i c h a n d r a n - V a n m a r c k e model shown in eqns (43) and (44) is chosen again to describe the coherence functions 7jk(•), j, k = 1,2, 3; j # k between fj°(t) and f°(t). The model p a r a m e t e r s A, a, k, w0 and b are assigned the values f r o m the previous example shown in eqn (49). The B o g d a n o f f - G o l d b e r g - B e r n a r d model is selected again for the m o d u l a t i n g functions:

Aj(~,t)=Aj(t)=altexp(-a2t),

j=1,2,3

(57)

with model p a r a m e t e r s al and a2 set equal to al = 0.680

and

a 2 = 1/4

(58)

so that the m a x i m u m value o f Aj(t), j = 1,2, 3 occurs at time t = 4 s and is equal to unity.

Non-stationary stochastic vector processes

8

R~I1 --* c o n t i n u o u s line

R~ 1 ~ c o n t i n u o u s line

R l l ---, dotted line

R21 ~ dotted line

x

O

O

04

804

O O O"

O O

oo.

R~ 2 --* continuous line

Rs°I --* continuous line

R22 ---* dotted line

Rsl ~ dotted line

i

oo -X

161

F-

"

l

80

~

R~s2 ~ continuous line

R~ss --, continuous line ---, dotted line

O

oo tO O O

z --' dotted line

8

r (sec) Fig. 8. Ensemble auto-/cross-correlation functions [Rjk(!, t + T)] computed from 100 000 sample functions versus the corresponding targets [R°k(t, t + T)]. Both Rjk(t, t + r) and Ryk(t, t + T) are plotted at time instant t = 3.14s as a function of T.

131.8

0.0

-131.8 . 0.0

.

.

6.79

.

.

13.57

.

.

.

20.36

27.15

Time (sec) Fig. 9. Acceleration record from the 1964 Niigata earthquake.

33.93

G. Deodatis

162

,S~j(w,=t7sec)(dividedby170)

Point 2

s e c ) [8Poin

t (/Y" _

50 m - - -

~ ' ~ o i nrlt

l

Fig. 10. Configuration of points 1, 2 and 3 on the ground surface, for Example 2. Points 1, 2 and 3 correspond to the same local soil conditions. At this juncture, it should be pointed out that the definitions for Sj(~), j = 1,2, 3 (see eqns (42), (53), (54), (55), (56)) and Aj(~,, t) (see eqn (57)) do not imply anymore that the modulating function Aj-(~, t) represents the change in the evolutionary power spectrum, relative to the (stationary) power spectral density function Sj(~). (Note that Sj.(~) is now a function of both frequency and time since the characteristic frequency and damping ratio of the ground are functions of time, and Aj(w, t) is a function of time only.) The evolutionary power spectra S°.(w, t), j = 1,2, 3 (see eqn (4)) are plotted in Fig. 11 at three time instants: t = 1, 4 and 7 s. As seen in this figure, the sudden drop in the values of the characteristic frequency and damping ratio of the ground from 4-5 to 5"5 s causes the formation of a very sharp peak in the evolutionary power spectra at frequency 2"0 rad/s (refer to eqns (53) and (54)). Since the non-stationary stochastic vector process in this example is not a uniformly modulated one, the generation of its sample functions can only be performed using eqn (24). The upper cut-off frequency Wu and the value of N (see eqn (27)) are set equal to 0;u = 128rad/s

and

N = 1024

(59)

The reason for using such a high value for N (compare to the corresponding value in eqn (52)) is to describe accurately the very sharp peak in the evolutionary power spectra after time instant 5"5 s (refer to Fig. 11). The simulation is performed at 2000 time instants, with a time step At = 0.01 s, over a length equal to 2000 x 0.01 = 20s. One generated sample function for the acceleration at points 1, 2 and 3 (see Fig. 10), denoted by fl (t), f2 (t) and f3 (t), respectively, is displayed in Fig. 12. It is obvious that the general frequency and amplitude variation characteristics of the acceleration record of the 1964 Niigata earthquake shown in Fig. 9 are reproduced very well in the sample function plotted in Fig. 12. This is accomplished by using the time-varying characteristic frequency and damping ratio of the ground shown in eqns (53) and (54), and the modulating function shown in eqn (57). In addition, it is possible to detect the correlation among fl (t), f2(t) andfa(t) in Fig. 12, as major peaks in the three time histories retain their general shape after

o1'0

2'0

~,

,~

~

Frequency ~(rad/sec)

Fig. 11. The evolutionary power spectra S° (w, t); j = 1,2, 3 at three time instants (Example 2). some loss of coherence described by the coherence functions 7jk(W), j, k = 1,2, 3, j ~ k (refer to eqns (43), (44) and (49)). Finally, it is possible to show that the ensemble auto-/cross-correlation function converges to the target auto-/cross-correlation function as the number of sample functions increases (refer to the first numerical example for a similar demonstration).

Example 3 In the first two numerical examples, sample functions of ground motion time histories were generated according to a target cross-spectral density matrix. The only drawback in this approach is that it is often preferable to work with time histories that are compatible with prescribed response spectra, rather than with prescribed power spectral density functions (cross-spectral density matrix), as indicated earlier in the section 'Simulation of seismic ground motion compatible with pre-scribed response spectra.' In order to address this issue in this numerical example, ground motion time histories are modeled as a uniformly modulated non-stationary stochastic vector process, but unlike the first example, sample functions will be now generated to be compatible with prescribed response spectra. The acceleration time histories at three points on the ground surface along the line of main wave propagation are considered one more time to be a tri-variate non-stationary stochastic vector process. The configuration of the three points is shown in Fig. 13 where the arrow indicates the direction of wave propagation. The simulated earthquake ground motion at points 1, 2 and 3 will have the following characteristics: (i) Points 1, 2 and 3 correspond to different local soil conditions. Specifically, point 1 corresponds to the Uniform Building Code's Soil Type 1 (rock and stiff soils), point 2 to Soil Type 2 (deep cohesionless or stiffclay soils), and point 3 to Soil Type 3 (soft to medium clays and sands).

Non-stationary stochastic vector processes

163 Point 1

v

0

15

20

Point 2

Co,

;

- - - "

'

(5

10

0

" - -

20

Point 3

C o oO4

tl ..%o

b

io

5

15

Time (sec)

2'o

Fig. 12. Generated sample function for the acceleration at points 1, 2 and 3, over a length equal to 20 s (Example 2). (ii) The acceleration time histories at these three points will be correlated according to a prescribed coherence function and they will reflect the wave propagation effect according to a specified velocity of wave propagation. (iii) Finally, the amplitude variation as a function of time of the acceleration time histories at these three points will be described by prescribed modulating functions Aj(co, t), j = 1,2, 3. The three components of the tri-variate non-stationary stochastic vector process describing the acceleration time histories at points 1, 2 and 3 (see Fig. 13) are denoted by fj0 (t), f o (t) and f3° (t), respectively.

The acceleration response spectra specified by the Uniform Building Code 4° are selected for the three points in Fig. 13. For the purposes of this numerical example, the peak ground acceleration is set equal to 200 cm/s2 and the corresponding Uniform Building Code (UBC) acceleration response spectra RSAj(CO), j = 1,2, 3 are plotted in Fig. 14. It is reminded that sample functions of ground motion time histories will be generated to be compatible with these UBC response spectra. The Abrahamson model 32 is chosen to describe the (stationary) coherence functions ~jk(CO), j , k = 1,2,3, j ~ k between fj°(t) a n d f ° ( t ) :

1

"yjk(co) =

CO

6

73

k /

I~

50 m

-,--'-

50 m

,-1

x tanh J

c3 ((jk) CO

CO2

L1 + ~c4(~jk) + 4-~ cT(~jk) Fig. 13. Configuration of points l, 2 and 3 on the ground surface along the line of main wave propagation, for Example 3. Point 1 corresponds to rock and stiff soils (UBC Type 1), point 2 corresponds to deep cohesionless or stiff clay soils (UBC Type 2), and point 3 corresponds to soft to medium clays and sands (UBC Type 3).

+ [4"80 - c3((j~)] exp c6(~jk)~-~ + 0"35 ,

j , k = 1,2,3,

j Ck

(60)

164

G. Deodatis

ff,5% damping

&

ge-

i

.<

o

0'.5

1'.o

5'.0

1~.o

d.o

,oko

F~'equencyw (rad/sec)

Fig. 14. Uniform Building Code acceleration response spectra RSAj(~v),j = 1,2, 3, assigned to points 1, 2 and 3 (see Fig. 13), respectively. Point 1 corresponds to rock and stiff soils (UBC Type 1), point 2 corresponds to deep cohesionless or stiff clay soils (UBC Type 2), and point 3 corresponds to soft to medium clays and sands (UBC Type 3) (Example 3). where ~jk is the distance between points j and k, and 3.95 c3(~jk) -- (1 + 0"0077~yk + 0"000023 ~j2k) + 0"85 exp{--0-000 13~yk}

0.4 1

(61)

3 1+

E1+ \190JJ l+\l--~j] c6(~jk)=3(exp{-~o}-l)-O'O018~/k

(63)

Since the acceleration time histories j~0 (t), j = 1,2, 3 are modeled as a uniformly modulated non-stationary stochastic vector process, the modulating functions A/(~v,t), j = 1,2,3 will be functions of time only. In order to control the duration of strong ground motion, the model suggested by Jennings et ai.4! is selected for the modulating functions A/(t), j = 1,2,3. Figure 16 plots AI (t) and indicates the numerical values chosen for all the parameters (A2(t) and Aa(t) have the same form as Al(t), but they are shifted by ~21/v and ~31/v, respectively, to include the wave propagation effect in the same way as in eqn (45)). The generation of sample functions of the acceleration time histories at the three points shown in Fig. 13 will be performed using the iterative scheme shown in Table 1 (refer to the section entitled 'Simulation of seismic ground motion compatible with prescribed response spectra'). According to Table 1, the power spectral density functions Sy(~v),j = 1,2, 3 must be initialized by setting them equal to a constant value over the entire frequency range. This value can be selected arbitrarily and in this example is set equal to: Sj(o3) : 100 cm2/s 3,

j = 1,2, 3

(68)

while the upper cut-off frequency wu and the value of N (see eqn (27)) are set equal to: ~vu = 128 rad/s

and

N = 128

(69)

The simulation is performed at 6144 time instants, with a time step At = 3.07 x 10 -3 s, over a length equal to 6144 x 3.07 x 1 0 - 3 = 18-85s. One sample function for the acceleration at points 1, 2 and 3 (see Fig. 13), denoted by fl(t), f2(t) and j~(t), respectively, is generated after 10 iterations and displayed in Fig. 17. As in the previous two examples, it is possible to detect the following characteristics in the time histories shown in Fig. 17: (i) their mutual correlation according to the coherence functions defined in eqns (60)-(65); (ii)

C7(~jk) = -0.598 + 0-1061n(~jk + 325)

-- 0"0151 exp{--0"6~yk} C8(~yk) = exp{8.54 -- 1.07 ln(~yk + 200)} + 100 exp{-- (jk}

(64)

O

(65)

Abrahamson's model for the coherence function has the advantage that it can be used for a broad range of soil conditions. 32 The 7jk(~v), j, k = 1,2, 3, j ¢ k defined in eqns (60)-(65) are plotted in Fig. 15. The velocity of wave propagation v is set equal to: v = 2000 m/s

(66)

so that the complex coherence functions can be expressed as: Fjk(0))

=

~/jk(W)exp[--i~L-~,

j,k-- 1,2,3,

j ~k

Frequency w (rad/sec)

(67)

Fig. 15. Coherence functions 721(~v), 731(~v) and 732(~v), for Example 3 (Abrahamson's model32).

N o n - s t a t i o n a r y stochastic vector processes

"go ,~ 0

5

10

15

T i m e t (sec)

Fig. 16. Modulating function Al(t), for Example 3 (Jennings

165

The iterative scheme shown in Table 1 is therefore capable of simulating seismic ground motion time histories that: (i) are compatible with prescribed response spectra (that can be different at different points on the ground surface); (ii) are correlated according to a given coherence function; (iii) include the wave propagation effect; and (iv) have an amplitude variation as a function of time according to a prescribed modulating function.

et al. model4°).

CONCLUSIONS the wave propagation effect according to a velocity of wave propagation v = 2000m/s; (iii) their amplitude variation as a function of time according to the modulating function plotted in Fig. 16. In this example, however, the main objective is to generate acceleration time histories that will be compatible with prescribed response spectra. For this purpose, the acceleration response spectra RSA(g)(w), j = 1,2, 3 computed using the ground motion time histories shown in Fig. 17 should be compared with the target UBC response spectra RSAj(w), j = 1,2,3 plotted in Fig. 14. This comparison is carried out in Fig. 18, where it can be seen that the 10 iterations performed to obtain the sample function shown in Fig. 17 are enough for an excellent match at every frequency. 0

A spectral-representation-based simulation algorithm was used to generate sample functions of a nonstationary, multi-variate stochastic process with evolutionary power, according to its prescribed non-stationary cross-spectral density matrix. If the components of the vector process correspond to different locations in space, then the process is also non-homogeneous in space (in addition to being non-stationary in time). The ensemble cross-correlation matrix of the generated sample functions is identical to the corresponding target and the generated sample functions are Gaussian in the limit as the number of terms in the frequency discretization of the cross-spectral density matrix approaches infinity. Point 1 --~ rock and stiffsoils ( U B C Type I)

~-. o.

" - ~0

"5-

1" 0

15 '

20 '---

Point 2 --* deep cohesionless or stiff clay soils (UBC Type 2) U

- - "

0

5'

o0

;o

;s

20 "--

Point 3 --* soft to medium clays and sands (UBC Type 3)

~

,,e

0

5

•

10

15

20

Time (sec) Fig. 17. Generated sample function for the acceleration at points 1, 2 and 3, after 10 iterations, over a length equal to 18"85s (Example 3).

166

G. Deodatis

if-. ~q 5% damping

5% damping

;-4 t)

;'4

/ 1,4

/

RSA(I~)(w) ~ oscillating line

0:5 s'.o ,~.o Frequency w (rad/sec) '

t)

'

'

O.

RSA2(w)--' smooth line RSA(h)(w) ~ oscillating line

o'.s 5:0 ~.o Frequency w (rad/sec) '

'

'

2~

RS A3 (o~) ---, smooth line RS A (h )(w ) ---, oscillating line

0'.5 '

5:0 '

~.o '

Frequency w (rad/sec)

Fig. 18. Acceleration response spectra RSA (fJ)(w), j = 1,2, 3 computed using the generated acceleration time histories shown in Fig. 17 versus the UBC acceleration response spectra RSAj(w), j = 1,2, 3. For the important application of earthquake ground motion simulation, an iterative scheme was introduced to generate seismic ground motion time histories at several locations on the ground surface that are compatible with prescribed response spectra, correlated according to a given coherence function, include the wave propagation effect, and have a specified duration of strong ground motion. Three examples involving simulation of earthquake ground motion were presented in order to demonstrate the capabilities of the proposed methodologies. In all three examples, the acceleration time histories at three points on the ground surface were considered to be a trivariate, non-stationary stochastic vector process. In the first example, the time histories were modeled as a uniformly modulated non-stationary stochastic vector process, and sample functions were generated according to a target cross-spectral density matrix. In the second example, the time histories were modeled as a nonstationary stochastic vector process with amplitude and frequency modulation, and sample functions were again generated according to a target cross-spectral density matrix. Finally, in the third example, ground motion time histories were modeled as a uniformly modulated non-stationary stochastic vector process, but unlike the first example, sample functions were generated to be compatible with prescribed response spectra. The generated acceleration time histories at the three points were spatially correlated according to prescribed coherence functions (all three examples), included the wave propagation effect (Examples 1 and 3), and were non-homogeneous in space, or equivalently corresponded to different local soil conditions (Examples 1 and 3). The generated ground motion time histories can be directly used as input for the dynamic seismic analysis of elongated structures such as bridges, lifelines, retaining walls, etc.

ACKNOWLEDGEMENTS This work was supported by the National Science Foundation under Grant # BCS-9257900 with Dr Clifford J. Astill as Program Director, and by the N C E E R Highway Project ( F H W A Contract DTFH61-92-C00106).

REFERENCES 1. Elishakoff, I., Probabilistic Methods in the Theory of Structures, John Wiley, New York, 1983. 2. Kozin, F., Auto-regressive moving-average models of earthquake records. J. Probabilistic Engng Mech., 3 (1988) 58-63. 3. Soong, T. T. & Grigoriu, M., Random Vibration of Mechanical and Structural Systems, Prentice-Hall, Hemel Hempstead, 1993. 4. Grigoriu, M., Applied Non-Gaussian Processes, PrenticeHall, Hemel Hempstead, 1995. 5. Deodatis, G., Simulation of stochastic processes and fields to model loading and material uncertainties, Probabilistic Methods for Structural Design, ed. Carlos Guedes Soares. Kluwer Academic, Dordrecht, 1996. 6. Rice, S. O., Mathematical analysis of random noise, Selected Papers on Noise and Stochastic Processes, ed. Nelson Wax, Dover, New York, 1954, pp. 133-294. (paper originally appeared in two parts in the Bell System Technical Journal 23, July 1944 and 24, January 1945). 7. Shinozuka, M. & Jan, C.-M., Digital simulation of random processes and its applications. J. Sound and Vibration, 25 (1972) 111-28. 8. Shinozuka, M., Monte Carlo solution of structural dynamics. Computers and Structures, 2 (1972) 855-74. 9. Yang, J.-N., Simulation of random envelope processes. J. Sound and Vibration, 25 (1972) 73-85. 10. Yang, J.-N., On the normality and accuracy of simulated random processes. J. Sound and Vibration, 26 (1973) 41728. 11. Shinozuka, M., Digital simulation of random processes in engineering mechanics with the aid of FFT technique, Stochastic Problems in Mechanics, eds S. T. Ariaratnam

Non-stationary stochastic vector processes

12. 13. 14. 15.

16. 17. 18. 19.

20. 21. 22. 23. 24. 25. 26.

27.

and H. H. E. Leipholz, University of Waterloo Press, Waterloo, 1974, pp. 277-86. Deodatis, G. & Shinozuka, M., Simulation of seismic ground motion using stochastic waves. J. Engng Mech., ASCE, 115 (1989) 2723-37. Yamazaki, F. & Shinozuka, M., Digital generation of nonGaussian stochastic fields. J. Engng Mech., ASCE, 114 (1988) 1183-97. Grigoriu, M., On the spectral representation method in simulation. J. Probabilistic Engng Mech., 8 (1993) 75-90. Shinozuka, M., Stochastic fields and their digital simulation, Stochastic Methods in Structural Dynamics, eds G. I. Schu~ller and M. Shinozuka, Martinus Nijhoff Publishers, Dordrecht, 1987, pp. 93-133. Shinozuka, M. & Deodatis, G., Stochastic process models for earthquake ground motion. J. Probabilistic Engng Mech., 3 (1988) 114-23. Shinozuka, M. & Deodatis, G., Simulation of stochastic processes by spectral representation. Appl. Mech. Rev., ASME, 44 (1991) 191-204. Shinozuka, M. & Deodatis, G., Simulation of multidimensional Gaussian stochastic fields by spectral representation. Appl. Mech. Rev., ASME, 49 (1996) 29-53. Gersch, W. & Yonemoto, J., Synthesis of multi-variate random vibration systems: a two-stage least squares ARMA model approach. J. Sound and Vibration, 52 (1977) 553-65. Mignolet, M. P. & Spanos, P.-T. D., Recursive simulation of stationary multi-variate random processes. Part I. J. Appl. Mech., 54 (1987) 674-80. Mignolet, M. P. & Spanos, P.-T. D., Recursive simulation of stationary multi-variate random processes. Part II. J. Appl. Mech., 54 (1987) 681-7. Li, Y. & Kareem, A., Simulation of multivariate nonstationary random processes by FFT. J. Engng Mech., ASME, 117 (1991) 1037-58. Li, Y. & Kareem, A., Simulation of multi-variate random processes: hybrid DFT and digital filtering approach. J. Engng Mech., ASCE, 119 (1993) 1078-98. Ramadan, O. & Novak, M., Simulation of spatially incoherent random ground motions. J. Engng Mech., ASCE, 119 (1993) 997-1016. Deodatis, G., Simulation of ergodic multi-variate stochastic processes. J. Engng Mech., ASCE, 122 (1996). Gurley, K. & Kareem, A., On the analysis and simulation of random processes utilizing higher order spectra and wavelet transforms. Proc. Second International Conf. on Computational Stochastic Mechanics, Athens, Greece, 1995, Balkema, pp. 315-24. Spanos, P.-T. D. & Zeldin, B., Random field simulation using wavelet bases. Proc. ICASP-7 Conf., Paris, France, 10-13 July 1995, pp. 1275-83.

167

28. Priestley, M. B., Evolutionary spectra and non-stationary processes. Jr. Royal Statistical Society, Series B, 27 (1965) 204-37. 29. Priestley, M. B., Non-linear and Non-stationary Time Series Analysis, Academic Press, New York, 1988. 30. Gasparini, D. & Vanmarcke, E. H., Simulated earthquake motions compatible with prescribed response spectra. Technical Report, Department of Civil Engineering, Massachusetts Institute of Technology, Publication No. R76-4, 1976. 31. Hao, H., Oliveira, C. S. & Penzien, J., Multiple-station ground motion processing and simulation based on SMART-1 array data. Nuclear Engng and Design, 111 (1989) 293-310. 32. Abrahamson, N. A., Spatial variation of multiple support inputs. Proc. 1st U.S. Seminar on Seismic Evaluation and Retrofit of Steel Bridges, University of California at Berkeley, San Francisco, 18 October, 1993. 33. Clough, R. W. & Penzien, J., Dynamics of Structures, McGraw-Hill, New York, 1975. 34. Bogdanoff, J. L., Goldberg, J. E. & Bernard, M. C., Response of a simple structure to a random earthquaketype disturbance. Bulletin of the Seismological Society of America, 51 (1961) 293-310. 35. Harichandran, R. S. & Vanmarcke, E. H., Stochastic variation of earthquake ground motion in space and time. J. Engng Mech., ASCE, 112 (1986) 154-74. 36. Ellingwood, B. R. & Batts, M. E., Characterization of earthquake forces for probability-based design of nuclear structures. Technical Report BNL-NUREG-51587, NUREG/CR-2945, Department of Nuclear Energy, Brookhaven National Laboratory, Prepared for Office of Nuclear Regulatory Research, U.S. Nuclear Regulatory Commission, 1982. 37. Hindy, A. & Novak, M., Pipeline response to random ground motion. J. Engng Mech., ASCE, 106 (1980) 33960. 38. Harichandran, R. S. & Wang, W., Effect of spatially varying seismic excitation on surface lifelines. Proc. Fourth U.S. National Conf. on Earthquake Engineering, Palm Springs, 20-24 May 1990, pp. 885-94 (Vol. 1). 39. Deodatis, G. & Shinozuka, M., Auto-regressive model for non-stationary stochastic processes. J. Engng Mech., ASCE, 114 (1988) 1995-2012. 40. International Conference of Building Officials, Uniform Building Code. Vol. 2, 1994. 41. Jennings, P. C., Housner, G. W. & Tsai, N. C., Simulated earthquake motions. Technical Report, Earthquake Engineering Research Laboratory, California Institute of Technology, 1968.