Orthogonal expansion of Gaussian wind velocity field and PDEM-based vibration analysis of wind-excited structures

J. Wind Eng. Ind. Aerodyn. 99 (2011) 1207–1220

Contents lists available at SciVerse ScienceDirect

Journal of Wind Engineering and Industrial Aerodynamics journal homepage: www.elsevier.com/locate/jweia

Orthogonal expansion of Gaussian wind velocity field and PDEM-based vibration analysis of wind-excited structures Zhang-Jun Liua, Jian-Bing Chenb, Jie Lib,n a b

College of Hydraulic and Environmental Engineering, China Three Gorges University, Yichang 443002, PR China State Key Laboratory of Disaster Reduction in Civil Engineering, Tongji University, 1239 Siping Road, Shanghai 200092, PR China

a r t i c l e i n f o

abstract

Article history: Received 12 September 2009 Received in revised form 13 September 2011 Accepted 16 September 2011 Available online 26 October 2011

An effective procedure for simulation of random wind velocity field by the orthogonal expansion method is proposed in this paper. The procedure starts with decomposing the fluctuating wind velocity field into a product of a stochastic process and a random field, which represent the time property and the spatial correlation property of wind velocity fluctuations, respectively. By an innovative orthogonal expansion technology, the stochastic process for wind velocity fluctuations may be represented as a finite sum of deterministic time functions with corresponding uncorrelated random coefficients. Similarly, the random field can be expressed as a combination form with only a few random variables by the Karhunen–Loeve decomposition. This approach actually simulates the wind velocity field with stochastic functions other than methods such as spectral representation and proper orthogonal decomposition. In the second part of the paper, the probability density evolution method (PDEM) is employed to predict the stochastic dynamic response of structures subjected to wind excitations. In the PDEM, a completely uncoupled one-dimensional partial differential equation, the generalized density evolution equation, plays a central role in governing the stochastic responses of structures. The solution of this equation will give rise to instantaneous probability density function of the responses. Finally, the accuracy and effectiveness of the approach in representing the random wind velocity field and PDEMbased dynamic response of wind-excited building are investigated. & 2011 Elsevier Ltd. All rights reserved.

Keywords: Gaussian wind velocity field Stochastic processes Probability density evolution method Karhunen–Loeve decomposition Orthogonal expansion

1. Introduction In the wind engineering and structural dynamics, it is very common to deal with wind loads as random functions of time and space. In many applications, it is necessary to transform a stochastic time-dependent wind field into a multivariate onedimensional stochastic vector process. For the simulation of multivariate vector processes, the following three approaches are now available: (1) the spectral representation method (SRM) (see, e.g. Shinozuka, 1971; Shinozuka and Deodatis, 1996; Deodatis, 1996; Cao et al., 2000); (2) the autoregressive moving average (ARMA) and autoregressive (AR) algorithms (see, e.g. Spanos and Mignolet, 1990; Li and Kareem, 1990; Di Paola and Gullo, 2001); and (3) the proper orthogonal decomposition (POD) (see, e.g. Di Paola, 1998; Chen and Kareem, 2005; Chen and Letchford, 2005). All the above methods for generating multivariate stochastic processes require the decomposition of covariance or cross power spectral density matrix as a product of two

n

Corresponding author. Tel.: þ86 21 65983526; fax: þ 86 21 65986345. E-mail addresses: [email protected] (Z.-J. Liu), [email protected] (J.-B. Chen), [email protected] (J. Li). 0167-6105/$ - see front matter & 2011 Elsevier Ltd. All rights reserved. doi:10.1016/j.jweia.2011.09.008

or three matrices, where the Cholesky (lower or upper triangular) or eigenvector decomposition procedure is usually adopted. However, in all of these methods a large number of random variables are needed to represent wind field processes, which will lead to quite low computing efficiency in predicting the dynamic response of wind-excited structures. The Karhunen–Loeve (K–L) decomposition is also known as the POD. The K–L decomposition that can handle an extensive class of stationary and non-stationary, Gaussian and non-Gaussian processes is widely investigated, providing approximations for stochastic processes consisting of finite sums of deterministic time functions combined with random variables (see, e.g. Spanos and Ghanem, 1989; Zhang and Ellingwood, 1994; Huang et al., 2001; Phoon et al., 2002a, 2002b; Phoon et al., 2004a, 2004b; Stefanou and Papadrakakis, 2007). However, the efficiency of the K–L decomposition for representing stochastic processes depends crucially on the availability of accurate eigenvalues and eigenfunctions of the integral K–L equation. Generally, the eigenvalues and eigenfunctions of the integral K–L equation need to be found numerically. This situation persists even when dealing with stationary processes if the analysis is based on practical problem. This paper proposed a new procedure for simulation of random wind velocity field with only a few random variables.

1208

Z.-J. Liu et al. / J. Wind Eng. Ind. Aerodyn. 99 (2011) 1207–1220

The procedure starts with decomposing the fluctuating wind velocity field into a product of a stochastic process and a random field, which represent the time property and spatial correlation property of the wind velocity fluctuation respectively. Then the stochastic process for wind velocity fluctuations is represented as a finite sum of deterministic time functions with corresponding uncorrelated random coefficients by an innovative orthogonal expansion technology. Similarly, the random field is expressed as a combination form with only a few random variables by the Karhunen–Loeve decomposition. The proposed technique provides a basis to employ the probability density evolution method (PDEM) which was developed by Li and Chen (2004, 2006, 2008) and Chen and Li (2005) in the past few years to study the stochastic responses of general linear/nonlinear structural systems. In fact, using this method, the instantaneous probability density function and its evolution of arbitrary physical quantities of interest of the linear/nonlinear stochastic systems can be captured successfully and the reliability evaluation of MDOF systems could be obtained (see, e.g. Zhang et al., 2008). A numerical example, which deals with a MDOF frame structure subjected to wind loads, is given for the purpose of illustrating the proposed approach.

correlation matrix R is defined as R ¼ ½rij NN

ð5Þ

Clearly, the correlation matrix R is a real symmetrical matrix. Let us define a random vector C as C ¼ ½c1 ðyÞ,c2 ðyÞ, . . ., cN ðyÞT , the superscript T denoting transpose. Generally C is a correlated random vector because R is usually a non-diagonal matrix. In this context, a correlation matrix decomposition can be employed to transform the correlated random vector C into a vector of standard uncorrelated random variables fxj ðyÞ, j ¼ 1,2, . . ., Ng by C¼

N X

Uj

qffiffiffiffi

lj xj ðyÞ

ð6Þ

j¼1

where lj’s and Uj’s are, respectively, the eigenvalues and standard eigenvectors of the correlation matrix R, i.e. RUj ¼ lj Uj

ð7Þ

From Eq. (6), the zero-mean random variables ck ðyÞ can also be expressed as ck ðyÞ ¼

N qffiffiffiffi X

lj xj ðyÞjjk

k ¼ 1,2, . . ., N

ð8Þ

j¼1

2. The orthogonal expansion of stochastic processes using the Hartley base Consider a real-valued stochastic process uðy,tÞ defined on a probability space ðO,A,PÞ and indexed on a bounded domain ð0,TÞ. Assume the process has a mean uðtÞ and a finite variance, E½uðy,tÞuðtÞ2 , which is bounded for all t A ð0,TÞ. Without loss of generality, the process with zero-mean and second-order moment is considered. It is known that such a stochastic process can be expanded on a basis that involves any complete set of deterministic functions with corresponding random coefficients, similar to the fact that a continuous deterministic function can be expressed on a basis that involves any complete set of deterministic functions with corresponding deterministic coefficients. The series is usually truncated after N terms as follows:

where jjk is the kth element of the standard eigenvector Uj . Substituting Eq. (8) in Eq. (1) yields ^ y,tÞ ¼ uð

N X N qffiffiffiffi X

N qffiffiffiffi X

k¼1j¼1

j¼1

lj xj ðyÞjjk fk ðtÞ ¼

^ y,tÞ ¼ uð

ck ðyÞfk ðtÞ

ð1Þ

k¼1

where the symbol y indicates the random nature of the corresponding quantity and ffk ðtÞ, k ¼ 1,2, . . ., 1g is a set of realvalued normalized orthogonal bases satisfying ( Z T 1, k¼l fk ðtÞfl ðtÞdt ¼ dkl ¼ /fk , fl S ¼ ð2Þ 0, otherwise 0 Here /US denotes the inner product, dkl is the Kronecker’s delta. In Eq. (1) ck(y)’s are zero-mean real random variables given by the projection of uðy,tÞ on the base function fk ðtÞ Z T ck ðyÞ ¼ /u, fk S ¼ uðy,tÞfk ðtÞdt ð3Þ 0

It should be noted that the right hand side of Eq. (3) is a random integral and should be understood in the sense of mean square. Using Eq. (3), we obtain the correlation coefficient Z T Z T E½ci ðyÞcj ðyÞ ¼ E½ uðy,t 1 Þfi ðt 1 Þdt 1 uðy,t 2 Þfj ðt 2 Þdt 2  0

Z

T

0

Z

T

¼ 0

¼ rij ,

0

Ru ðt 1 ,t 2 Þfi ðt 1 Þfj ðt 2 Þdt 1 dt 2

i,j ¼ 1,2, . . ., N

ð4Þ

where the auto-covariance function Ru ðt 1 ,t 2 Þ ¼ E½uðy,t 1 Þuðy,t 2 Þ is bounded, symmetric and positive definite. Accordingly, the

ð9Þ

P where f j ðtÞ ¼ N k ¼ 1 jjk fk ðtÞ. It is easy to prove that the functions ff j ðtÞ, j ¼ 1,2, . . ., Ng are a set of normalized orthogonal functions within time domain ð0,TÞ satisfying Z T N X f i ðtÞf j ðtÞdt ¼ jjk jik ¼ UTi Uj ¼ dij ð10Þ /f i ,f j S ¼ 0

k¼1

^ y,tÞ in Eq. (9) converges to the Obviously, the process uð original process uðy,tÞ as N-N, i.e. 1 qffiffiffiffi X ^ y,tÞ ¼ uðy,tÞ ¼ lim uð lj xj ðyÞf j ðtÞ ð11Þ N-1

N X

lj xj ðyÞf j ðtÞ

j¼1

The parameters xj(y)’s in Eq. (11) are a set of uncorrelated random variables which can be expressed as Z T 1 1 xj ðyÞ ¼ qffiffiffiffi /u,f j S ¼ qffiffiffiffi uðy,tÞf j ðtÞdt ð12Þ

lj

lj

0

with the mean and covariance given, respectively, by E½xj ðyÞ ¼ 0

ð13Þ

E½xi ðyÞxj ðyÞ ¼ dij

ð14Þ

Then the mean-square absolute error e1 in representing the exact Eq. (11) by the approximate Eq. (9) is given by Z T  1 X ^ y,tÞÞ2 dt ¼ ðuðy,tÞuð lk ð15Þ e1 ¼ E 0

k ¼ Nþ1

where E½U denotes the mathematical expectation. The series expansion in Eq. (11), referred to as the orthogonal expansion of a stochastic process, provides a second-moment characterization in terms of uncorrelated random variables and deterministic orthogonal functions. It is noted that the orthogonal expansion of a stochastic process is similar to the K–L decomposition of the process in expression (Phoon et al., 2002b). In many practical applications where the random quantities vary smoothly with respect to time, only a few terms are necessary to capture the major part of the random fluctuation of the process. Its major advantage is the reduction from a large

Z.-J. Liu et al. / J. Wind Eng. Ind. Aerodyn. 99 (2011) 1207–1220

number of correlated random variables to several most important uncorrelated ones. Thus, the series in Eq. (9) is further approximated by a reduced-order terms, say r (r 5N), such that ~ y,tÞ ¼ uð

In a word, the orthogonal expansion scheme using the Hartley orthogonal bases can be briefly summarized as follows: ^ y,tÞ ¼ uð

r qffiffiffiffi X

lj xj ðyÞf j ðtÞ

1209

N qffiffiffiffi pffiffiffi X 2 lj xj ðyÞf j ðtÞ

ð23aÞ

j¼1

ð16Þ

j¼1

The corresponding auto-covariance function can be evaluated theoretically by R~ u ðt 1 ,t 2 Þ ¼

r X

lj f j ðt1 Þf j ðt2 Þ

ð17Þ

j¼1

The mean-square absolute error e2 in substituting Eq. (16) for Eq. (9) yields e2 ¼ E½

Z

T

N X

^ y,tÞuð ~ y,tÞÞ2 dt ¼ ðuð

0

lk

ð18Þ

f j ðtÞ ¼

N 1 X

jj,k þ 1 fk ðtÞ

ð23bÞ

k¼0

where deterministic functions fk(t)’s are the Hartley orthogonal basic functions, i.e.   1 2kpt , k ¼ 0,1,2, . . ., ðN1Þ ð24Þ fk ðtÞ ¼ pffiffiffi cas T T where the function casðUÞ is defined as casðtÞ ¼ cosðtÞ þ sinðtÞ

ð25Þ

k ¼ rþ1

Considering Eqs. (15) and (18), the mean-square absolute error e in substituting Eq. (16) for Eq. (11) is given by

3. Orthogonal expansion of wind velocity fluctuations

Z e ¼ E½

3.1. Power spectrum model of wind velocity fluctuations

T

~ y,tÞÞ2 dt ðuðy,tÞuð

0

Z T ^ y,tÞÞ þ ðuð ^ y,tÞuð ~ y,tÞÞÞ2 dt ððuðy,tÞuð ¼ E½ 0 1 X

¼ e1 þ e2 ¼

lk

ð19aÞ

k ¼ r þ1

The mean-square relative error e in substituting Eq. (16) for Eq. (11) yields

e ¼ RT 0

e1 þ e2 E½u2 ðy,tÞdt

¼ RT 0

e

ð19bÞ

Ru ðt,tÞdt

2

To implement orthogonal expansion of stochastic processes efficiently, the key issue is to select N base functions, say fi ðtÞ, i ¼ 0,1, . . ., ðN1Þ. Here we assume that N is a large odd number. For normalized trigonometric basis functions, functions up to (N  1)/2 harmonics over the interval ð0,TÞ are defined as pffiffiffi   1 2kpt 2 pffiffiffi cos ð20aÞ fð1Þ fð1Þ 0 ðtÞ ¼ pffiffiffi , 2k1 ðtÞ ¼ T T T pffiffiffi   2kpt 2 p ffiffiffi , fð1Þ ðtÞ ¼ sin 2k T T

k ¼ 1,2, . . ., ðN1Þ=2

ð20bÞ

For the Hartley orthogonal basis function, functions up to (N  1)th order over the interval ð0,TÞ are defined as (Bracewell, 1986) 1

fð2Þ 0 ðtÞ ¼ pffiffiffi ,

T      1 2kpt 2kpt pffiffiffi sin þ cos , fð2Þ k ðtÞ ¼ T T T

In wind engineering, there are usually two approaches to obtain the power spectral density (PSD) function of wind velocity fluctuation processes based on the statistical result of field observed data. One approach is to get the PSD curve as the output of a super low frequency filter exposed to the strong wind data. In the other approach, however, the correlation function curve of the wind velocity is firstly obtained from the field observed data, then the Fourier transform is performed to yield the PSD function. According to the first approach, the PSD of wind velocity fluctuation processes was proposed by Davenport (1961) as follows:

k ¼ 1,2, . . ., ðN1Þ

ð21Þ

Obviously, the relationship between the trigonometric basis function and Hartley basis function are presented by pffiffiffi ð2Þ ð1Þ fð1Þ 2fk ðtÞ, k ¼ 1,2, . . ., ðN1Þ=2 ð22Þ 2k1 ðtÞ þ f2k ðtÞ ¼ It is noted that, although the Hartley basis function is not a complete set of orthogonal series, no information is lost by using the Hartley basis function for a real-valued time function instead of the trigonometric basis function. For the purpose of reference and completion, a summary of the relationship between the Hartley and Fourier transforms is given in Appendix A, which is mainly based on Rodriguez (2003).

Sv ðoÞ ¼ 4K 0 V 10

x2

ð26Þ

oð1 þx2 Þ4=3

where Sv ðoÞ is the one-sided power spectrum of the wind velocity fluctuations (m2/s); K 0 is the non-dimensional surface roughness coefficient; V 10 is the mean wind speed at 10 m height (m/s); o is the frequency in radian per second; x is the integral length scale of turbulence, i.e. x¼

co pV 10

ð27Þ

where the parameter c ¼ 600 m. The variance of the Davenport power spectrum can be approximated by a constant (Li et al., 1993) 2

s2v 6ZK 0 V 10

ð28Þ

where the parameter Z ¼ 6:0. On the other hand, the normalized autocorrelation function of wind pressure fluctuations was given by ]apvte{h as follows (Li et al., 1993): R w ðtÞ ¼ ea9t9 ðcos bt þ m sin b9t9Þ

ð29Þ

The inverse Fourier transform can be written as Z 1 1 R w ðtÞ expðiotÞdt S w ðoÞ ¼ 2p 1 2

¼

1 ðambÞo2 þ ða þ mbÞða2 þ b Þ

p o4 þ2ða2 b Þo2 þ ða2 þ b Þ 2

2 2

ð30Þ

where R w ðtÞ is the normalized autocorrelation function of wind pressure fluctuations, accordingly, S w ðoÞ is the normalized PSD function; a, b, m satisfy m r ða=bÞ. It is noted that the ]apvte{h power spectrum may be unequal to zero at frequency o ¼ 0. In order to remove the influence of the

1210

Z.-J. Liu et al. / J. Wind Eng. Ind. Aerodyn. 99 (2011) 1207–1220

zero frequency part, the following relationship must be satisfied:

m¼

25

a

ð31Þ

b

Considering the relationship of the PSD function between the wind velocity and the wind pressure fluctuation, the two-sided PSD of wind velocity fluctuations can be expressed as

po

4 þ 2ð 2 b2 Þ

a

o2 þða2 þ b2 Þ2

ð32Þ

where the parameter A is a constant.

Equivalent spectrum

PSD (m2/s)

2ao2

A S v ðoÞ ¼

Davenport spectrum

20

15

10

3.2. Equivalent power spectrum for wind velocity fluctuations 5

0

1

2

30

3 4 Frequency (rad/s)

5

6

5

6

Davenport spectrum Equivalent spectrum

25

2

Similarly, the predominant frequency and the peak value of the power spectrum (two-sided) in Eq. (32) are, respectively, qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi op ¼ a 2 þ b 2 ð34aÞ

0

35

PSD (m /s)

An equivalent power spectrum can be defined based on the expression in Eq. (32). The principle of equivalence is that the equivalent power spectrum and the target power spectrum have the same variance, the same predominant frequency and the same peak value, respectively. Thus, the parameters A, a, b in Eq. (32) can be determined according to the equivalence principle. For example, the predominant frequency and the peak value of the Davenport power spectrum (one-sided) (cf. Eq. (26)) are, respectively, pffiffiffiffiffiffi 15 pV 10 ð33aÞ op ¼ 5 c pffiffiffiffiffiffi ffiffiffi 15K 0 cV 10 p 3 Sv ðop Þ ¼ 5 ð33bÞ 4p

20 15 10

S v ðop Þ ¼

A

ð34bÞ

2pa

Therefore, according to the principle of equivalence, the following equations hold by letting the counterparts in Eqs. (33) and (34) be equal: qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi pffiffiffiffiffiffi 15 pV 10 ð35aÞ a2 þ b2 ¼ 5 c pffiffiffiffiffiffi ffiffiffi 15K 0 cV 10 p A 3 2 ¼ 5 ð35bÞ 2pa 4p In addition, employing the variance equivalence principle, the parameterA can be easily defined as 2

2

A ¼ s2v ¼ ZK 0 V 10 ¼ 6:0K 0 V 10

ð36Þ

In turn, the parameters a, b can be solved from Eqs. (35) and (36), i.e.

a ¼ 3:6239

V 10 c

V 10 i ¼ gi c pffiffiffiffiffiffiffi where i ¼ 1 is the imaginary unit, g is given by

b ¼ 2:6853

g ¼ 2:6853

V 10 c

ð37aÞ

ð37bÞ

ð38Þ

In this way, the equivalent power spectrum of wind velocity fluctuations in Eq. (32) is determined via the Davenport power spectrum. The comparison between the Davenport’s spectrum and its equivalent spectrum at the mean wind speed 20 m/s and 30 m/s

5 0

0

1

2

3 4 Frequency (rad/s)

Fig. 1. Comparison between Davenport spectrum and its equivalent spectrum: (a) the mean wind speed is 20 m/s; (b) the mean wind speed is 30 m/s.

are shown in Fig. 1, respectively. Clearly they agree very well. From the spectrum shape, it is easy to see that the PSD function of wind velocity fluctuations decreases rapidly with the increase of frequency and the predominant frequency is approximately close to 0:004V 10 (rad/s). Hence, the energy of wind velocity fluctuations is distributed mainly in the bounded domain (0–2p) rad/s. The energy beyond this part is negligible.

3.3. Pseudo-stochastic process of wind displacement fluctuations The physical meaning of the wind displacement fluctuation is not definitely defined. But the physical relationship between the displacement and the velocity is clear. For convenience, the integral of the wind velocity fluctuation is defined as the pseudostochastic process of wind displacement fluctuation. The PSD of the pseudo-stochastic process of wind displacement fluctuation uðy,tÞ can be deduced from the equivalent PSD of the wind velocity fluctuation vðy,tÞ given by Eq. (32). In fact, the two-sided PSD function of the pseudo-stochastic process uðy,tÞ is

Z.-J. Liu et al. / J. Wind Eng. Ind. Aerodyn. 99 (2011) 1207–1220

given as 2a

1 A Su ðoÞ ¼ 2 S v ðoÞ ¼

p o4 þ 2ða2 b2 Þo2 þ ða2 þ b2 Þ2 A 2a ¼ p o4 þ 2ða2 þ g2 Þo2 þða2 g2 Þ2 o

ð39Þ

where the parameters A, a, g are the same as those in Eqs. (36)–(38). The corresponding autocorrelation function of the pseudostochastic process uðy,tÞ deduced from the Fourier transform of Eq. (39) is given by Z 1 Ru ðtÞ ¼ Su ðoÞ expðiotÞdo 1   A expððagÞ9t9Þ expðða þ gÞ9t9Þ ¼  ð40Þ 2g ag aþg

fluctuation. Generally speaking, the value N ¼ 60021000 can meet the corresponding error eN 51.0. In most situations, only a small number of terms associated with the large eigenvalues are dominant. This attribute helps to realize a reduced-order modeling of the wind displacement fluctuation process in which higher terms associated with small eigenvalues may be further truncated as ~ y,tÞ ¼ uð

r qffiffiffiffi pffiffiffi X 2 lj xj ðyÞf j ðtÞ

In wind engineering, the stochastic process of wind velocity fluctuations is usually assumed as a zero-mean stationary ergodic Gaussian process. In this case, the complete probabilistic characterization is ensured by the PSD function such as the Davenport spectrum. Usually it is difficult to represent the stochastic process of wind velocity fluctuations with only a few series expansion terms (Li and Liu, 2006). In order to capture the main probabilistic characteristics of wind velocity fluctuation process, it is proposed to carry out the orthogonal expansion in terms of the pseudo-wind displacement fluctuation process instead of the wind velocity process. According to the orthogonal expansion scheme of stochastic processes using the Hartley orthogonal base, the pseudo-stochastic process of wind displacement fluctuation uðy,tÞ can be approximately represented as N qffiffiffiffi pffiffiffi X ^ y,tÞ ¼ 2 lj xj ðyÞf j ðtÞ uð

where the non-dimensional random variable set fx1 ðyÞ, x2 ðyÞ, . . ., xr ðyÞg satisfies E½xj ðyÞ ¼ 0,

j ¼ 1,2, . . ., r

N 1 X

jj,k þ 1 fk ðtÞ

ð41aÞ

ð41bÞ

k¼0

where fk(t)’s pffiffiffi are the Hartley orthogonal basic functions, i.e. fk ðtÞ ¼ ð1= T Þcasð2kpt=TÞ, k ¼ 0,1,2, . . ., ðN1Þ, the duration of wind velocity fluctuations T¼600 s; lj’s and Uj ’s are, respectively, the eigenvalues and standard eigenvectors of the wind displacement correlation matrix R, i.e. RUj ¼ lj Uj ; and fj,k þ 1’s are the (kþ1)th element of the standard eigenvector Uj , i.e. Uj ¼ ½jj,1 , jj,2 , . . ., jj,k þ 1 , . . ., jj,N T . The wind displacement correlation matrix R can be defined as R ¼ ½rij NN

ð42Þ

ð43Þ

0

where the wind displacement auto-covariance function Ru ðt 2 t 1 Þ ¼ E½uðy,t 1 Þuðy,t 2 Þ, i.e. Eq. (40). In Eq. (41), the following criterion is usually used to estimate the number of the expanded term N: Z ou Z 1 Sv ðoÞdo ¼ ð1eN Þ Sv ðoÞdo ð44aÞ 0

ou ¼ 2pN=T

i,j ¼ 1,2, . . ., r

ð46bÞ

where dij is the Kronecker’s delta. In Eq. (45), r is the number of the remaining terms. The cumulative percentage of the energy located in the first r terms may be defined as follows: Pr

PðrÞ ¼ PNj ¼ 1

k¼1

lj

lk

 100%

ð47Þ

Generally speaking, we require PðrÞ Z80%, thus the number r can be determined. Considering the relationship between the wind velocity fluctuation and the pseudo-wind displacement fluctuation process, the wind velocity fluctuation process can be approximately expanded as ~ y,tÞ ¼ vð

r qffiffiffiffi pffiffiffi X 2 lj xj ðyÞGj ðtÞ

ð48aÞ

j¼1

Gj ðtÞ ¼

N 1 X

wk þ 1 jj,k þ 1 f_ k ðtÞ

ð48bÞ

where the overdot denotes the time derivative, the undetermined numbers wk þ 1 are defined as harmonic modulated coefficients. Practically, it may be assumed that the wind velocity fluctuation is a zero-mean Gaussian process, and then the random variables xj ðyÞ in Eq. (48a) are independent standard Gaussian variables. From Eq. (48), the autocorrelation function of the wind velocity fluctuation process can be written as ~ y,t 1 ÞUvð ~ y,t 2 Þ ¼ 2 R~ v ðt 1 ,t 2 Þ ¼ E½vð

r X

lj Gj ðt1 ÞGj ðt2 Þ

ð49aÞ

j¼1

or R~ v ðt,t þ tÞ ¼ 2

r X

lj Gj ðtÞGj ðt þ tÞ

ð49bÞ

j¼1

where the correlation coefficient rij is as follows: Z TZ T rij ¼ Ru ðt 2 t 1 Þfi ðt 1 Þfj ðt 2 Þdt 1 dt 2 , i,j ¼ 0,1, . . ., ðN1Þ 0

ð46aÞ

k¼0

j¼1

f j ðtÞ ¼

ð45Þ

j¼1

E½xi ðyÞxj ðyÞ ¼ dij , 3.4. Orthogonal expansion of wind velocity fluctuations

1211

0

ð44bÞ

where eN 51:0 (e.g. eN ¼ 0:10) is the tolerable relative truncated error, ou is the cut-off frequency of the PSD for the wind velocity

Performing Fourier transform in terms of the variable t on both sides of Eq. (49b) yields Z 1 1 R~ v ðt,t þ tÞexpðiotÞdt ð50Þ S~ v ðt, oÞ ¼ 2p 1 The time average of S~ v ðt, oÞ in terms of t in the time period ð0,TÞ is thus Z 1 T ~ Sv ðoÞ ¼ S v ðt, oÞdt T 0 2 3 2   1  r X 2 NX 2k p 2kp 2 2 4 ð51Þ wk þ 1 lj jj,k þ 1 5d o ¼ Tk¼1 T T j¼1 where dðUÞ is the Dirac’s function.

1212

Z.-J. Liu et al. / J. Wind Eng. Ind. Aerodyn. 99 (2011) 1207–1220

From Eqs. (49)–(51), it is easy to see that Z Z 1  1 T  1 T ~ ~ y,tÞUvð ~ y,tÞ dt ¼ E vð Sv ðoÞdo R v ðt,tÞdt ¼ T 0 T 0 0 Z

ð52Þ

Thus, the function Sv ðoÞ may be defined as the mean-power spectral density of the wind velocity fluctuation process during ð0,TÞ. According to the principle of energy equivalence, i.e. Sv ðoÞ ¼ Sv ðoÞ, the values of the coefficients wk þ 1 , k¼0, 1, 2, y, (N 1), can be solved through the following equations:   Z 3p=T r 2 2p 2 2 X w1 ¼ 0, w2 lj j2j,2 ¼ Sv ðoÞdo ð53aÞ T T 0 j¼1  2 Z ð2k þ 1Þp=T r X 2 2kp w2k þ 1 lj j2j,k þ 1 ¼ Sv ðoÞdo, T T ð2k1Þp=T j¼1

4. Orthogonal expansion of wind velocity field Let x, y, z be a point in the space, where z is the height from the ground and x is the lateral wind direction, y is the along-wind direction. The velocity at a given point can be regarded either as a one-variate four-dimensional (1V-4D) random field or a time-dependent 1V-3D stochastic field process. In practice, in order to simplify the concepts, the wind velocity fluctuation can be characterized in the plane y¼0. Then the time-dependent 1V-2D stochastic field process can be examined. Thus, the wind velocity field can be written as

where a is the coefficient of ground roughness degree, V 10 is the mean wind speed at 10 m above the ground. In the case the PSD of wind velocity fluctuation does not vary with height, the fluctuating wind velocity field V~ ðy; x,z; tÞ can be decomposed as (Ou and Wang, 1998) ð56Þ

where vðy,tÞ is the fluctuating wind velocity stochastic process whose complete probabilistic characterization is ensured by the PSD function. The wind velocity fluctuation process vðy,tÞ can be represented using the orthogonal expansion method in Eq. (48). Uðy; x,zÞ is the spatial correlation random field, of which the correlation function can be expressed as (Ou and Wang, 1998) 2 !1=2 3 2 2 ðx x Þ ðz z Þ 2 1 2 1 5 ð57Þ þ RU ðx1 ,x2 ; z1 ,z2 Þ ¼ exp4 L2x L2z where Lx ¼ 50 m, Lz ¼ 60 m. In Eq. (57), the correlation function RU ðx1 ,x2 ; z1 ,z2 Þ can also be approximately expressed as (Simiu and Scanlan, 1996) RU ðx1 ,x2 ; z1 ,z2 Þ6RU x ðx1 ,x2 ÞURUz ðz1 ,z2 Þ

ð58Þ

r x pffiffiffiffiffiffi X

lxi xxi ðyÞf xi ðxÞ

ð61Þ

i¼1

where lxi and f xi ðxÞ are the eigenvalues and corresponding eigenfunctions of the correlation function RU x ðx1 ,x2 Þ, xxi ðyÞ is a set of uncorrelated random variables and r x is the number of K–L terms. The deterministic eigenvalues lxi and eigenfunctions f xi ðxÞ are obtained from the solution of the homogeneous Fredholm integral equation of the second kind Z RUx ðx1 ,x2 Þf xi ðx1 Þdx1 ¼ lxi f xi ðx2 Þ ð62Þ Dx

where Dx is the bounded domain of the random field. The corresponding correlation function can be approximated by the eigen-decomposition RUx ðx1 ,x2 Þ6

rx X

lxi f xi ðx1 Þf xi ðx2 Þ

ð63Þ

i¼1

Similarly, the random field U z ðy,zÞ and its correlation function RUz ðz1 ,z2 Þ in Eq. (60) can be decomposed as

ð54Þ

where the symbol y indicates the random nature of the corresponding quantity and Vðy; x,z; tÞ is the wind speed at the point ðx,zÞ (unit: m/s), VðzÞ is the mean wind speed at the height z (unit: m/s) and V~ ðy; x,z; tÞ is the fluctuating component of the wind speed (unit: m/s), which is assumed to be a normal zero-mean stationary stochastic process. As is known, the mean wind speed VðzÞ can be represented by the power-law (Ou and Wang, 1998) z a VðzÞ ¼ V 10 ð55Þ 10

V~ ðy; x,z; tÞ ¼ Uðy; x,zÞvðy,tÞ

For the one-dimensional homogeneous zero-mean random field U x ðy,xÞ with the exponential correlation function in Eq. (59), a finite Karhunen–Loeve (K–L) series can be employed (Spanos and Ghanem, 1989)

k ¼ 2,3, . . ., ðN1Þ

The sign of the coefficientswk þ 1’s are given empirically by Liu (2007).

ð59Þ

Likewise, RU z ðz1 ,z2 Þ is the vertical correlation function, i.e.   jz2 z1 j RUz ðz1 ,z2 Þ ¼ exp  ð60Þ Lz

U x ðy,xÞ6

ð53bÞ

Vðy; x,z; tÞ ¼ VðzÞ þ V~ ðy; x,z; tÞ

where RU x ðx1 ,x2 Þ is the lateral correlation function, i.e.   jx2 x1 j RUx ðx1 ,x2 Þ ¼ exp  Lx

U z ðy,zÞ6

rz qffiffiffiffiffiffi X

lzj xzj ðyÞf zj ðzÞ

ð64Þ

lzj f zj ðz1 Þf zj ðz2 Þ

ð65Þ

j¼1

RUz ðz1 ,z2 Þ6

rz X j¼1

Thus, the two-dimensional homogeneous random field Uðy; x,zÞ and its correlation function RU ðx1 ,x2 ; z1 ,z2 Þ in Eq. (58) can be approximately expressed as Uðy; x,zÞ6

rx pffiffiffiffiffiffi X

r z qffiffiffiffiffiffi X

i¼1

j¼1

lxi xxi ðyÞf xi ðxÞ

RU ðx1 ,x2 ; z1 ,z2 Þ6

rx X i¼1

lzj xzj ðyÞf zj ðzÞ

lxi f xi ðx1 Þf xi ðx2 Þ

rz X

lzj f zj ðz1 Þf zj ðz2 Þ

ð66Þ

ð67Þ

j¼1

where the zero-mean uncorrelated random variables xxi ðyÞ and

xzj ðyÞ satisfy   E xxi ðyÞxzj ðyÞ ¼ 0

ð68aÞ

  E xxi ðyÞxxk ðyÞ ¼ dik

ð68bÞ

  E xzj ðyÞxzl ðyÞ ¼ djl

ð68cÞ

in which dUU is the Kronecker’s delta; i,k ¼ 1,2, . . ., r x ; j,l ¼ 1,2,. . .,r z . Furthermore, the cross-power spectral density (cross-PSD) function of the along-wind fluctuations ðx1 ,z1 Þ and ðx2 ,z2 Þ can be expressed as SV~ ðx1 ,x2 ,z1 ,z2 ; oÞ ¼ RU ðx1 ,x2 ; z1 ,z2 ÞSv ðoÞ

ð69Þ

where RU ðx1 ,x2 ; z1 ,z2 Þ is the autocorrelation function, Sv ðoÞ is the Davenport’s spectrum.

Z.-J. Liu et al. / J. Wind Eng. Ind. Aerodyn. 99 (2011) 1207–1220

1213

5. PDEM-based dynamic response analysis of wind-excited structures

Eqs. (74)–(76) can be numerically solved in convenience by the following steps:

5.1. Probability density evolution equation of dynamic response

(i) Select representative point set in the domain OH and denote it as hq ¼ ðy1,q , y2,q , . . ., ys,q Þ (q ¼ 1,2, . . ., Nsec ), where Nsec is the total number of the selected points. Simultaneously, determine the corresponding assigned probability Pq. (ii) For the prescribed H ¼ hq , solve the deterministic ordinary differential equation (70) under the initial condition (71) to evaluate the value of X_ ðhq ,t m Þ, where t m ¼ mDt (m ¼0, 1, 2, y), Dt is the time step. (iii) Employ hq in place of h, and X_ ðhq ,t m Þ in place of X_ ðh,tÞ in Eq. (74). Solve the partial differential equation (74) under the initial condition (75) with the finite difference method to obtain the discretized value pX H ðxj , hq ,t k Þ, where xj ¼ jDx (j ¼ 0, 7 1, 7 2, . . .), Dx is the space step, t k ¼ kUDt^ (k ¼ 0,1,2, . . .), Dt^ is the time step in the finite difference

Without loss of generality, the equation of motion of a multidegree-of-freedom (MDOF) structural system under wind loads can be expressed as € þCXðtÞ _ þKXðtÞ ¼ FðH,tÞ MXðtÞ

ð70Þ

subjected to the deterministic initial condition XðtÞ9t ¼ 0 ¼ x0 ,

_ XðtÞ9 ¼ x_ 0 t¼0

ð71Þ

€ X _ and X are the n  1 acceleration, velocity and where X, displacement vector with the overdot denoting the derivative with respect to time t; M, C and K are the n  n mass, damping and stiffness matrix, respectively; F is the n  1 wind excitation vector, which is studied and represented by a finite sum of randomly modulated time functions in the preceding sections (cf. Eqs. (48), (61), (64) and (66)), H ¼ fx1 ðyÞ, x2 ðyÞ, . . ., xnH ðyÞg is the nH  1 random vector with known probability density function (PDF) pH ðhÞ. Obviously, the dynamic response XðtÞ is a random vector process dependent on and determined by H, and can be expressed as XðtÞ ¼ HðH,tÞ

ð72Þ

where H, existent and unique for a well-posed problem, is a deterministic operator. Its any component can be written as XðtÞ ¼ HðH,tÞ

ð73Þ

Here the subscript denoting the components is omitted for the simplicity of writing. Obvious, the response component X(t) is a random function. According to Li and Chen (2004, 2006, 2008) and Chen and Li (2005), the joint PDF of the augmented state vector ðX, HÞ satisfies the governing partial differential equation @pX H ðx, h,tÞ _ @p ðx, h,tÞ þ X ðh,tÞ X H ¼0 @t @x

ð74Þ

where X_ ðh,tÞ is the velocity of the response for a prescribed h. It is easy to specify the initial condition of Eq. (74), i.e. pX H ðx, h,tÞ9t ¼ 0 ¼ dðxx0 ÞpH ðhÞ

ð75Þ

where x0 is the component of x0 as a deterministic value. Once the initial-value problem (Eqs. (74) and (75)) is solved, the PDF of XðtÞ could then be evaluated by Z pX ðx,tÞ ¼ pX H ðx, h,tÞdh ð76Þ OH

where OH is the distribution domain of H. Eq. (74) is referred to as the generalized density evolution equation. Remarkably, it tells that even for a multi-dimensional dynamical system a one-dimensional partial differential equation governing the joint PDF exists. In addition, the physical meaning of Eq. (74) says that as long as the velocity of a physical quantity is known, the PDF of the quantity could be easily captured because varying of the PDF results from varying of the state of the quantity.

method. (iv) Carry out the numerical integration in Eq. (76) to obtain the numerical solution of the PDF of the response. In step (i), special techniques are needed in the strategy of selecting representative points in OH . For instance, the numbertheoretical-method-based algorithm developed by Li and Chen (2007) is employed in the present paper (see Appendix B). In step (iii), the TVD schemes are appropriate in dynamic response analysis while a hybrid difference scheme is preferred in evaluation of the extreme value distribution and the reliability (Chen and Li, 2005). From the above steps, it is seen that a traditional deterministic dynamic response analysis process is embedded in the procedure. Additional finite difference method is then employed. This makes it convenient to use the existent commercial software.

6. Model analysis and validation In order to verify the feasibility and validity of the proposed orthogonal expansion of the wind velocity fluctuation process, the equivalent power spectrum of the Davenport’s spectrum is analyzed. Here K0 ¼0.0033, the mean wind speed V 10 is 20 m/s and 30 m/s, respectively. The parameter values of the orthogonal expansion for wind velocity fluctuation processes are given in Table 1. Fig. 2 shows the eigenvalues of the correlation matrix R of the pseudo-wind displacement fluctuations process. It is seen that the eigenvalues decrease rapidly with the increase of mode number. The former 10 eigenvalues are dominant, making up about 80% total energy of the process. According to the above proposed methodology, the wind velocity time histories can be generated independently using the number theoretical method (NTM) (Li and Chen, 2007) to pick out the representative point sets in the multi-dimensional random parameters space. Fig. 3 shows two typical time histories of wind velocity fluctuation via the orthogonal expansion method when the mean wind speed is 20 m/s and 30 m/s, respectively.

Table 1 Parameter values of the orthogonal expansion of wind velocity fluctuation processes.

5.2. Numerical implementation procedures As is pointed out above, for a general dynamical system, the explicit expression of X_ ðh,tÞ is usually unknown, which makes Eq. (74) analytically intractable for most problems. Nonetheless,

V 10 (m/s)

N

eN (%)

r

PðrÞ (%)

20 30

600 600

6.5 8.5

10 10

89.7 80.5

1214

Z.-J. Liu et al. / J. Wind Eng. Ind. Aerodyn. 99 (2011) 1207–1220

to modify the sample sets " # Nsec M X X k2 u2k ðt m ÞP k Dt ¼ s2v T

2.0x107

Eigenvalues

m¼1

M ¼ T=Dt,

1.5x107

ð77aÞ

k¼1

t m ¼ Dtm

ð77bÞ u2k ðt m Þ

where k is the sample modified coefficient; is the variance value of the kth sample at the time t m (noting the mean value of the samples is zero); s2v is the variance of the Davenport’s spectrum; Nsec is the number of the samples or the number of selected points, in the present example Nsec ¼610; Dt is the time step, Dt¼0.1 s. Thus, the orthogonal expansion of the wind velocity fluctuation process in Eq. (48) can be further expressed as r qffiffiffiffi pffiffiffi X vðy,tÞ ¼ 2 lj xj ðyÞGj ðtÞ ð78aÞ

1.0x107

5.0x106

0.0 0

10

20

30 40 Mode number

50

60

j¼1

Gj ðtÞ ¼

N 1 X

Zk þ 1 jj,k þ 1 f_ k ðtÞ

ð78bÞ

k¼0

Eigenvalues

4x105

Zk þ 1 ¼ kwk þ 1

ð78cÞ

where Zk þ 1 (k ¼ 0,1,2, . . ., ðN1Þ) are defined as the energy equivalence coefficients.

3x105

102

2x105

Simulation Target

1x105

0 0

10

20

30 40 Mode number

50

60

PSD (m2/s)

101

100

Fig. 2. Eigenvalues of the pseudo-wind displacement process: (a) the mean wind value is 20 m/s; (b) the mean wind value is 30 m/s.

wind velocity (m/s)

10-1 20

10-2

10

10-2

0 -10 -20

10-1 Frequency (Hz)

100

102 0

100

200

300

400

500

Simulation

600

time (sec)

Target

20

PSD (m2/s)

wind velocity (m/s)

101

10 0

100

-10 -20

0

100

200

300 time (sec)

400

500

600

Fig. 3. Time history of wind velocity fluctuations: (a) the mean wind speed V 10 is 20 m/s; (b) the mean wind speed V 10 is 30 m/s.

Generally speaking, the total energy of samples in the duration T and the total energy of the original process might be inconsistent. In order to make them equal, a coefficient k is introduced

10-1

10-2

10-2

10-1

100

Frequency (Hz) Fig. 4. Comparison between the samples ensemble PSD and target PSD: (a) the mean wind speed is 20 m/s; (b) the mean wind speed is 30 m/s.

Z.-J. Liu et al. / J. Wind Eng. Ind. Aerodyn. 99 (2011) 1207–1220

Table 2 Eigenvalues of the first 5.

lz1

lnz1

lz2

lnz2

lz3

62.18

18.20

6.86

3.40

1.99

Story node j

Height zj (m)

Mean wind speed (m/s)

The tributary area Bj (m2)

10 9 8 7 6 5 4 3 2 1

100 90 80 70 60 50 40 30 20 10

39.91 38.66 37.32 35.86 34.24 32.41 30.31 27.81 24.62 20.00

160 300 300 300 300 300 300 300 300 300

Table 4 The integer vector ðn,h1 ,h2 , . . ., hs Þ for s ¼ 15 (h1 ¼ 1). 1139 691 1131 480 1106 904 1066 142 1009 487 937347 850242 748799 633750 505923 366239 215705 55406 1026 186 849882

2422 957 2398 094 2323 761 2200 720 2030 234 1814 052 1554 392 1253 920 915717 543256 140357 2134 112 1683 011 1214 641 733806

4395 774 4364 102 4269 316 4112 097 3893 578 3615 335 3279 371 2888 108 2444 365 1951 338 1412 580 831972 213699 3957 988 3277 986

14271 038 14168 215 13860 486 13350 069 12640 642 11737 315 10646 597 9376 347 7935 718 6335 088 4585 990 2701 027 693780 12849 750 10642 098

The along-wind loading generated by the orthogonal expansion method and the response of a wind-excited 100 m tall building with 10 DOFs are performed to illustrate proposed approach. The lumped mass of each story from the bottom to the top are 2.5, 2.5, 2.2, 2.2, 2.2, 2.2, 2.2, 2.2, 1.8, 1.8 (  108 kg). The stiffness of each story from the bottom to the top are 3.48, 2.93, 2.92, 2.94, 2.92, 2.91, 2.92, 2.92, 2.31, 0.43 (  104 kN/mm). The first and second natural periods of the building are 1.764 and 0.601 s, respectively. Rayleigh damping is employed, i.e.C ¼ aM þbK, where a¼0.2657Hz, b ¼0.0071 s, M is the mass matrix and K is the stiffness matrix. The modal damping ratio for each mode is assumed to be 0.05. The mean wind speed at the height zj above the ground is given by the power-law, i.e., Vðzj Þ ¼ V 10 ðzj =10Þ0:3 , where V 10 is the mean wind speed at 10 m above the ground, in the present example takes as the value 20 m/s. The one-sided cross-power spectral density (cross-PSD) of the along-wind fluctuations at the height zi and zj above the ground is given by    zj zi  x2 2 exp  ð79Þ SV~ ðzi ,zj ; oÞ ¼ 4K 0 V 10 Lz oð1þ x2 Þ4=3 where x ¼ 600o=ðpV 10 Þ, K 0 ¼ 0:0033 and Lz ¼ 60 m. The one-dimensional homogeneous zero-mean random field U z ðy,zÞ with correlation function in Eq. (60) can be approximated as (Spanos and Ghanem, 1989) ^ z ðy,z0 Þ ¼ U

3 qffiffiffiffiffiffi X

2 qffiffiffiffiffiffi X n

j¼1

j¼1

lzj xzj ðyÞf zj ðz0 Þ þ

lzj xnzj ðyÞf nzj ðz0 Þ

20 15

0

100

200

300 time (sec)

400

500

600

0

100

200

300 time (sec)

400

500

600

0

100

200

300 time (sec)

400

500

600

40 35 30 25

ð80Þ

where z is the height above the ground; z0 is the local coordinates, i.e., z0 ¼ zH=2, H is the height of building; the eigenvalues lzj and

25

wind velocity (m/s)

139489 138484 135476 130487 123553 114724 104063 91647 77566 61921 44825 26401 6781 125597 104019

Fig. 4 shows the comparison between the ensemble PSD obtained from the 610 samples of the generated wind velocity fluctuation time histories and the target PSD when the mean wind speed is 20 m/s and 30 m/s, respectively. Note that they accord very well, which shows the orthogonal expansion process and the original stochastic process are consistent in the sense of second order statistics.

7. Numerical example

Table 3 The mean wind speed and the tributary area.

n h2 h3 h4 h5 h6 h7 h8 h9 h10 h11 h12 h13 h14 h15

1215

50 40 30

Fig. 5. Samples of simulated time histories for wind velocity.

1216

Z.-J. Liu et al. / J. Wind Eng. Ind. Aerodyn. 99 (2011) 1207–1220

n

n

lzj and the eigenfunctions f zj ðz0 Þ and f zj ðz0 Þ are given in Appendix C, which is mainly based on (Ghanem and Spanos, 1991). Table 2 lists the first 5 eigenvalues which are arranged in a descending order. The cumulative percentage of the energy located in the first 5 eigenvalues is up to 92.63% of the total energy. The external loading vector is FðH,tÞ ¼ fF 1 ðH,tÞ,F 2 ðH,tÞ, . . ., F 10 ðH,tÞgT , where F j ðH,tÞ is the along-wind fluctuating force at the jth story node, which can be represented by F j ðH,tÞ ¼ mD 

1 rV 2 Bj 2

ð81Þ

10

Simulation Target

1

100

10-1

10-2

10-2

10-1 Frequency (Hz)

100

Simulation Target

1

100

10-1

10-2

10-2

10-1 Frequency (Hz)

Cross-PSD between z=30 and z=70 (m2/s)

Aoto-PSD of z=70 (m2/s)

10-2

10-2

10-1 Frequency (Hz)

100

Fig. 6. Comparison of the auto-power spectral density between the ensemble of samples and the original process: (a) the height at z¼ 30 m; (b) the height at z ¼50 m; (c) the height at z ¼70 m.

10-2

10-1 Frequency (Hz)

100

Simulation Target

101

100

10-1

10-2

10-1 Frequency (Hz)

102

Simulation Target

10-1

10-1

10-2

102

100

100

102

100

101

Simulation Target

101

10-2

Cross-PSD between z=50 and z=70 (m2/s)

Aoto-PSD of z=50 (m2/s)

102

10

102 Cross-PSD between z=30 and z=50 (m2/s)

Aoto-PSD of z=30 (m2/s)

102

where mD is the drag coefficient and assumed as 1.2; r is the air density; Bj is the tributary area as listed in Table 3. Employing the PDEM, one will first select representative points in the random vector space H ¼ fx1 ðyÞ, x2 ðyÞ, . . ., x15 ðyÞg using the number theoretical method (NTM) (Li and Chen, 2007). Some tables for values of the integer vector ðn,h1 ,h2 , . . ., h15 Þ are available in (Hua and Wang, 1981). Table 4 gives part of the tables. In the present paper, the parameter values of the above NTM are l ¼ 3:5, s ¼ 15, n¼ 2 422 957, r0 ¼ 1.58, Nsel ¼1038, gs ¼ 0:0428%. The representative time history of the wind velocity field fluctuation can be generated using the orthogonal expansion method in the preceding

100

Simulation Target

101

100

10-1

10-2

10-2

10-1 Frequency (Hz)

100

Fig. 7. Comparison of the cross power spectral density between the ensemble of samples and the original process: (a) the height at z¼ 30 m and 50 m; (b) the height at z¼ 50 m and 70 m; (c) the height at z ¼ 30 m and 70 m.

Z.-J. Liu et al. / J. Wind Eng. Ind. Aerodyn. 99 (2011) 1207–1220

sections. For the prescribed value of H, deterministic structural analysis is carried out to capture the velocity value of the response quantity and then the generalized density evolution Eq. (74) is solved and the numerical integration in Eq. (76) is performed.

Mean (mm)

150 100 50 0

0

100

200

300 400 Time (sec)

500

600

0

100

200

300 400 Time (sec)

500

600

Std.D (mm)

15 10 5 0

Fig. 8. The mean and the standard deviation.

Fig. 5 shows typical samples of the wind velocity time histories generated by the orthogonal expansion method at height z ¼10 m, 50 m and 100 m, respectively. The time step is chosen as 0.1 s. Fig. 6(a), (b) and (c) shows the comparison of the auto-PSDs based on the orthogonal expansion models with the target PSD functions at height z¼30 m, 50 m and 70 m, respectively. Fig. 7(a), (b) and (c) pictures the cross-PSDs based on the orthogonal expansion models with the target PSD functions at height z¼30 m and 50 m, 50 m and 70 m and 30 m and 70 m, respectively. Results shown in Figs. 6 and 7 demonstrate the effectiveness and accuracy of the orthogonal expansion model in describing the wind velocity field. The probabilistic information of the top displacement of the structure is shown in Figs. 8–11. Pictured in Fig. 8 are the mean and the standard deviation. Fig. 9 shows the typical PDF at certain instants of time. Fig. 10 is the surface constructed by the PDF at different instants of time while Fig. 11 is the contour of the surface, respectively. Some short remarks: It is noted that application of the proposed method is limited to analysis of a tall structure where the PSD variation with height can be neglected (considering the PSD of the wind velocity fluctuations does not vary with height only in the present paper). In addition, the proposed methodology, like the K–L decomposition, provides sample functions that are not ergodic (in mean and autocorrelation) in a sample-by-sample

0.14

140 PDF at 480.00 sec PDF at 540.00 sec PDF at 600.00 sec

Displacement (mm)

0.12

0.08 0.06 0.04 0.02 100

110 120 Displacement (mm)

130

130 120 110 100 90 588

140

590

Fig. 9. Typical PDFs at different instants of time.

592

594 Time (sec)

0.08 0.06 0.04 0.02 0 140 Dis

120 pla c

100 em ent (mm )

596

Fig. 11. The contour of the PDF surface.

0.1

PDF

PDF

0.1

0 90

1217

80

588

590

592

594

(s Time

Fig. 10. The PDF evolution against time.

ec)

596

598

600

598

600

1218

Z.-J. Liu et al. / J. Wind Eng. Ind. Aerodyn. 99 (2011) 1207–1220

To exploit this Hermitian symmetry and redundancy in the Fourier transform of real-valued functions, Hartley introduced an integral operation that permits to transform a real-valued temporal function into a real-valued frequency function. The Hartley transform of a real-valued function xðtÞ can be expressed as (Bracewell, 1986) Z 1 Hx ðf Þ ¼ xðtÞcasð2pf tÞdt ðA:4Þ

sense in contrast to the spectral representation method. However, how important is the property of ergodicity is still in argument because the practical stochastic processes are usually not ergodic themselves. Finally, as demonstrated in the preceding sections, the efficiency of orthogonal expansion method for simulating random processes is far superior to the currently existing methodologies, such as the spectral representation in which usually many random variables are involved and thus large computational efforts are needed, the proper orthogonal decomposition or the Karhunen–Loeve decomposition method where usually the Fredholm integral equation need to be solved to find the eigenvalues and eigenfunctions which is usually quite nontrivial.

and the inverse Hartley transform may be defined as Z 1 xðtÞ ¼ Hx ðf Þcasð2pf tÞdf

8. Conclusions

where the kernel transform function, known as the cosine-andsine or cas function, is defined as

1

casðtÞ ¼ cosðtÞ þsinðtÞ This paper proposes an orthogonal expansion model of the Gaussian wind velocity fluctuation. By decomposing the Gaussian wind velocity field into a product of a stochastic process and a random field, the stochastic process of the wind velocity fluctuation may be represented as a finite sum of deterministic time functions with corresponding uncorrelated random coefficients by the orthogonal expansion, and the random field can be expressed as a combination form with only a few random variables by the Karhunen–Loeve decomposition. In the procedure, 15 independent Gaussian random variables are selected, among them 10 are used to represent the time property of the wind velocity fluctuations and the other 5 are employed to represent the spatial correlation property of the random field. It is found that the stochastic response of structures can be predicted by combining the PDEM and the orthogonal expansion model of the Gaussian wind velocity field. An example, which deals with a tall building subjected to wind loads, is investigated to validate the above approach.

ðA:5Þ

1

ðA:6Þ

It is interesting to note that the Hartley transform exhibits a self-inverse property. That is, the direct and inverse transformations use the same integral operation. Note that the sum of the sine and cosine functions is just another sine function shifted by p=4. However, this is the key for the symmetry between the transform and its inverse exhibited by the Hartley transform. Some important relationships existing between the Hartley and Fourier transform are presented as follows (Bracewell, 1986): Z 1 Hx ðf Þ þ Hx ðf Þ Hex ðf Þ ¼ ¼ xðtÞcosð2pf tÞdt ¼ ReðF x ðf ÞÞ ðA:7Þ 2 1 where Hex ðf Þ is the even part of the Hartley transform; ReðUÞ is the real part of the Fourier transform. Similarly, the odd part of the Hartley transform, Hox ðf Þ, is the imaginary part, ImðUÞ, of the Fourier transform Z 1 Hx ðf ÞHx ðf Þ ¼ Hox ðf Þ ¼ xðtÞsinð2pf tÞdt ¼ ImðF x ðf ÞÞ ðA:8Þ 2 1 Then, the Fourier and Hartley transforms of xðtÞ can be readily expressed as

Acknowledgments

F x ðf Þ ¼ Hex ðf ÞiHox ðf Þ

The supports of the National Natural Science Foundation of China for Innovative Research Groups (Grant nos. 50321803 and 50621062), the National Natural Science Foundation of China (Grant nos. 50808113 and 10872148) and the Natural Science Foundation of Hubei Province of China (Grant nos. 2009CDA013 and D20111206) are gratefully appreciated.

Appendix A. Relationship between the Hartley and Fourier transforms (Rodriguez, 2003) The well known integral Fourier transform of a continuous function of time xðtÞ is given by (Bracewell, 1986) Z 1 F x ðf Þ ¼ xðtÞexpði2pf tÞdt ðA:1Þ

Hx ðf Þ ¼ ReðF x ðf ÞÞImðF x ðf ÞÞ

ðA:9Þ ðA:10Þ

It should be noted that the Fourier and Hartley transforms are very similar and that both transforms may be expressed as combinations of the sine and cosine transforms. They are related to one another by Eqs. (A.9) and (A.10). Another important fact is that the Hartley and Fourier transforms are invertible and consequently they carry the whole information about the original function but in a different way. Nevertheless, since the Fourier transform can be recovered from the Harley transform, no information is lost by using the Hartley transform for a realvalued time function instead of the Fourier transform.

Appendix B. The number theoretical method

1

The inverse transform can be written as Z 1 xðtÞ ¼ F x ðf Þexpði2pf tÞdf

ðA:2Þ

1

where the kernel transform function is expð 7i2pf tÞ ¼ cosð2pf tÞ 7isinð2pf tÞ

ðA:3Þ

For a real-valued functionxðtÞ, it is noteworthy that the Fourier transform is Hermitian; that is F x ðf Þ ¼ F nx ðf Þ, where the superscript n denotes complex conjugate. This implies that the Fourier transform is overspecified because dependence exists between transform values for positive and negative frequencies. Therefore, half of the transform is redundant.

According to Hua and Wang’s (1981), the points ðx1,k ,x2,k , . . ., xs,k Þ in a high-dimensional hypercube ½0,1s (1 rs r 18) can be determined by an integer generator vector ðn,h1 ,h2 , . . ., hs Þ, i.e.   hj k hj k int , k ¼ 1,2, . . ., n; j ¼ 1,2, . . ., s ðB:1Þ xj,k ¼ n n where n is the total number of the points, intðUÞ takes the remaining integer after trimming the fraction part. For convenience of application, some tables for values of the integer vector ðn,h1 ,h2 , . . ., hs Þ are available in Hua and Wang (1981). Table B1 is part of the tables fors ¼ 10. Now, suppose H ¼ fx1 ðyÞ, x2 ðyÞ, . . ., xs ðyÞg are mutually independent standardized normally distributed random variables.

Z.-J. Liu et al. / J. Wind Eng. Ind. Aerodyn. 99 (2011) 1207–1220

1219

Table B1 the integer vector ðn,h1 ,h2 , . . ., hs Þ for s ¼ 10 (h1 ¼ 1). n

h2

h3

h4

h5

h6

h7

h8

h9

h10

85633 103661 115069 130703 155093

37677 45681 65470 64709 90485

35345 57831 650 53373 20662

3864 80987 95039 17385 110048

54821 9718 77293 5244 102308

74078 51556 98366 29008 148396

30354 55377 70366 52889 125399

57935 37354 74605 66949 124635

51906 4353 55507 51906 10480

56279 27595 49201 110363 44198

The truncated boundary of random variables xj ðyÞ (j ¼ 1,2, . . ., s) is assumed as l, i.e. 9yj 9 r l. Thus, a translation and scale transformation of the uniformly scattered points in the hypercube ½l, ls can be given as

yj,k ¼ 2lðxj,k 0:5Þ, k ¼ 1,2, . . ., n;

j ¼ 1,2, . . ., s

where n is the total number of the points. Accordingly, the assigned probability ðy1,k , y2,k , . . ., ys,k Þ reads P~ k ¼ pH ðy1,k , y2,k , . . ., ys,k ÞV k ¼

1 ð2pÞs=2

exp 

y21,k þ y22,k þ    þ y2s,k 2

ðB:2Þ of

the

point

V k,

k ¼ 1,2, . . ., n

ðB:3Þ

where V k is the representative volume of the point ðy1,k , y2,k , . . .,

y21,q þ y22,q þ    þ y2s,q rðr0 lÞ2 , q ¼ 1,2, . . ., Nsel

ðB:4Þ pffiffi where r 0 ð0 o r 0 r sÞ is the bounded radius coefficient, and N sel is the number of the selected points. To measure the effect of the sieving condition (B.4), define a sieving ratio

H ¼100 m

H ¼120 m

H ¼150 m

b1 bn1 b2 bn2 b3 bn3

0.01607151 0.03941662

0.0143389 0.03381264

0.01240998 0.028085

0.06766213 0.09762945

0.05709365 0.08188635

0.04647865 0.066124

0.12824834 0.159166275

0.1072883 0.132977764

0.08631895 0.10678415

Spanos, 1991)  qffiffiffiffiffiffi 1 qffiffiffiffiffiffi X U z ðy,z0 Þ ¼ lzj xzj ðyÞf zj ðz0 Þ þ lnzj xnzj ðyÞf nzj ðz0 Þ

q ¼ 1,2, . . ., N sel

ðC:1Þ

where z is the height above the ground, z0 is the local coordinates, i.e. z0 ¼ zH=2, H is the total height of the building; the eigenn functions f zj ðz0 Þ and f zj ðz0 Þ are given by cosðbj z0 Þ f zj ðz0 Þ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi hþ ðsinð2bj hÞ=2bj Þ

ðC:2Þ

n

sinðbj z0 Þ n f zj ðz0 Þ ¼ qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi n n hðsinð2bj hÞ=2bj Þ

ðC:3Þ

where h ¼ H=2, namely, the random field indexed on a bounded n domain Dz0 ¼ ½h,h. The corresponding eigenvalues lzj and lzj are

lzj ¼

lnzj ¼

2c

ðC:4Þ

b2j þ c2 2c

ðC:5Þ

bjn2 þc2 n

ðB:5Þ

In the most general occasions, the points selected in Eq. (B.2) with the assigned probability in Eq. (B.3) can be further sieved through the assigned probability. Let the normalized assigned probability P~ q P q ¼ PN , sel P~ k¼1 k

Parameters

j¼1

Generally speaking, however, in the case that s is relatively large, to achieve an acceptable accuracy the needed n is also relatively large. For instance, if s ¼ 10, n might be in the order of 104 2106 . This number is still too large for general structural analysis. Fortunately, the number of selected points can further be reduced greatly such that Nsel 5 n by limiting the selected points in a hyper-ball or a hyper-ellipsoid rather than in a hypercube (Li and Chen, 2007). The joint PDF distributed in the r-dimensional random parameters space OH is usually spherically symmetric. Accordingly, the selected points in the random parameters space OH should naturally also be spherically symmetric. This will yield the condition

Nsel n

n

Values of the parameters bj and bj .

!

ys,k Þ, i.e. V k ¼ ð2lÞs =n.

gs ¼

Table C1

ðB:6Þ

In the presented paper, the parameter values of the above NTM are l ¼4.0, s ¼ 10, n ¼ 155093, r0 ¼1.05, Nsel ¼610, gs ¼0.393%.

Appendix C. The Karhunen–Loeve decomposition of random field with exponential covariance The Karhunen–Loeve decomposition of the homogeneous zero-mean random field with the exponential correlation function defined by Eq. (60) can be represented as (Ghanem and

in which c ¼ 1=Lz ¼ 1=60; bj and bj transcendental equations

can be solved by the

cb tanðbhÞ ¼ 0

ðC:6Þ

bn þ c tanðbn hÞ ¼ 0

ðC:7Þ

In Eqs. (C.6) and (C.7), the transcendental equations can be solved by numeric arithmetic when H is assigned. Thus, the eigenvalue problem of Eqs. (C.4) and (C.5) can be solved with its eigenvalues arranged in descending order to obtain the corresponding eigenfunctions. The cumulative percentage of the energy of the first r z eigenvalues may be defined as follows: Pr z

lj j ¼ 1 lj

Pðrz Þ ¼ Pj1¼ 1

ðC:8Þ

Generally, 90% of the energy of the random field needs to be retained, r z is simply chosen so that Pðr z Þ Z90%. n Table C1 lists the former 6 values of the parameters bj and bj when H assigned as 100 m, 120 m and 150 m, respectively.

1220

Z.-J. Liu et al. / J. Wind Eng. Ind. Aerodyn. 99 (2011) 1207–1220

References Bracewell, R.N., 1986. The Hartley Transform. Oxford University Press, New York. Cao, Y.H., Xiang, H.F., Zhou, Y., 2000. Simulation of stochastic wind velocity field on long-span bridges. Journal of Engineering Mechanics 126 (1), 1–6. Chen, L., Letchford, C.W., 2005. Simulation of multivariate stationary Gaussian stochastic processes: hybrid spectral representation and proper orthogonal decomposition approach. Journal of Engineering Mechanics 131 (8), 801–808. Chen, J.B., Li, J., 2005. Dynamic response and reliability analysis of nonlinear stochastic structures. Probabilistic Engineering Mechanics 20 (1), 33–44. Chen, X., Kareem, A., 2005. Proper orthogonal decomposition-based modeling, analysis, and simulation of dynamic wind load effects on structures. Journal of Engineering Mechanics 131 (4), 325–339. Davenport, A.G., 1961. The spectrum of horizontal gustiness near the ground in high winds. Quarterly Journal of the Royal Meteorolgical Society 87, 194–211. Deodatis, G., 1996. Simulation of ergodic multivariate stochastic processes. Journal of Engineering Mechanics 122 (8), 778–787. Di Paola, M., 1998. Digital simulation of wind field velocity. Journal of Wind Engineering and Industrial Aerodynamics 74-76, 91–109. Di Paola, M., Gullo, I., 2001. Digital generation of multivariate wind field processes. Probabilistic Engineering Mechanics 16 (1), 1–10. Ghanem, R., Spanos, P.D., 1991. Stochastic Finite Element: A Spectral Approach. Springer, Berlin. Hua, L.K., Wang, Y., 1981. Applications of Number Theory to Numerical Analysis. Springer, Science Press, Berlin, Beijing. Huang, S.P., Quek, S.T., Phoon, K.K., 2001. Convergence study of the truncated Karhunen–Loeve expansion for simulation of stochastic processes. International Journal for Numerical Methods in Engineering 52 (9), 1029–1043. Li, J., Chen, J.B., 2004. Probability density of evolution method for dynamic response analysis of structures with uncertain parameters. Computational Mechanics 34, 400–409. Li, J., Chen, J.B., 2006. The probability density evolution method for dynamic response analysis of nonlinear stochastic structures. International Journal for Numerical Methods in Engineering 65 (6), 882–903. Li, J., Chen, J.B., 2007. The number theoretical method in response analysis of nonlinear stochastic structures. Computational Mechanics 39 (6), 693–708. Li, J., Chen, J.B., 2008. The principle of preservation of probability and the generalized density evolution equation. Structural Safety 30 (1), 65–77. Li, J., Liu, Z.J., 2006. Expansion method of stochastic processes based on normalized orthogonal bases. Journal of Tongji University 34 (10), 1279–1283 (in Chinese). Li, Y., Kareem, A., 1990. ARMA systems in wind engineering. Probabilistic Engineering Mechanics 5 (2), 50–59.

Li, G.Q., Cao, H., Li, Q.S., 1993. The Theory of Dynamic Reliability of Structures and its Application. Earthquake Press, Beijing (in Chinese). Liu, Z.J., 2007. Orthogonal Expansion Theory for Stochastic Dynamic Excitations in Engineering and its Applications. Tongji University, Shanghai (in Chinese). Ou, J.P., Wang, G.Y., 1998. Random Vibration of Structures. Higher Education Press, Beijing (in Chinese). Phoon, K.K., Huang, S.P., Quek, S.T., 2002a. Simulation of second-order processes using Karhunen–Loeve expansion. Computers and Structures 80 (12), 1049–1060. Phoon, K.K., Huang, S.P., Quek, S.T., 2002b. Implementation of Karhunen–Loeve expansion for simulation using a wavelet-Galerkin scheme. Probabilistic Engineering Mechanics 17 (3), 293–303. Phoon, K.K., Quek, S.T., Huang, H., 2004a. Simulation of non-Gaussian processes using fractile correlation. Probabilistic Engineering Mechanics 19 (4), 287–292. Phoon, K.K., Huang, H.W., Quek, S.T., 2004b. Comparison between Karhunen–Loeve and wavelet expansions for simulation of Gaussian processes. Computers and Structures 82, 985–991. Rodriguez, G., 2003. Analysis and simulation of wave records through fast Hartley transform. Ocean Engineering 30, 2255–2273. Shinozuka, M., 1971. Simulation of multivariate and multidimensional random processes. Journal of the Acoustical Society of America 49 (1), 357–368. Shinozuka, M., Deodatis, G., 1996. Simulation of multi-dimensional Gaussian stochastic fields by spectral representation. Applied Mechanics Review 49 (1), 29–53. Simiu, E., Scanlan, R.H., 1996. Wind Effects on Structures3rd edition John Wiley & Sons, New York. Spanos, P.D., Ghanem, R., 1989. Stochastic finite element expansion for random media. Journal of Engineering Mechanics 115 (5), 1035–1053. Spanos, P.D., Mignolet, M.P., 1990. Simulation of stationary random processes: two stage MA to ARMA approach. Journal of Engineering Mechanics 116 (3), 620–641. Stefanou, G., Papadrakakis, M., 2007. Assessment of spectral representation and Karhunen–Loeve expansion methods for the simulation of Gaussian stochastic fields. Computer Methods in Applied Mechanics and Engineering 196, 2465–2477. Zhang, J., Ellingwood, B., 1994. Orthogonal series expansions of random fields in reliability analysis. Journal of Engineering Mechanics 120 (12), 2660–2677. Zhang, L.L., Li, J., Peng, Y.B., 2008. Dynamic response and reliability analysis of tall buildings subject to wind loading. Journal of Wind Engineering and Industrial Aerodynamics 96, 25–40.