Digital generation of multivariate wind field processes

Probabilistic Engineering Mechanics 16 (2001) 1–10 www.elsevier.com/locate/probengmech

Digital generation of multivariate wind field processes M. Di Paola*, I. Gullo Dipartimento di Ingegneria Strutturale e Geotecnica, Universita` degli Studi di Palermo, Viale delle Scienze, 90128, Palermo, Italy Received 1 May 1996; received in revised form 1 September 1999; accepted 1 November 1999

Abstract A very efficient procedure for the generation of multivariate wind velocity stochastic processes by wave superposition as well as autoregressive time series is proposed in this paper. The procedure starts by decomposing the wind velocity field into a summation of fully coherent independent vector processes using the frequency dependent eigenvectors of the Power Spectral Density matrix. It is shown that the application of the method allows to show some very interesting physical properties that allow to reduce drastically the computational effort. Moreover, using a standard finite element procedure for approximating the frequency dependent eigenvectors, the generation procedure requires the generation of a limited number of univariate fully coherent processes for describing the entire multivariate velocity processes independently of the number of components of the process. 䉷 2000 Elsevier Science Ltd. All rights reserved. Keywords: Multivariate wind field processes; Wave superposition; Power spectral density matrix

1. Introduction An effective approach for predicting the dynamic response to wind actions of elastic structures such as longspan suspended and cable-stayed bridges, towers, antennas and tall buildings, is the Monte Carlo simulation. It mainly consists of generating sample functions having prescribed probabilistic characteristics, then for each sample function the deterministic analysis is performed, and the probabilistic analysis is made on the response samples. The generation of wind velocity histories, as stationary normal processes is currently based on several probabilistic models. Early simulations were based on superposition of harmonic waves with random phases (CAWS and CAWS methods) (see, e.g. Refs. [1–7]). Recently, within the context of mechanical problems, attention has been devoted to simulation algorithms that synthesize stationary random processes as outputs of digital filters exposed to bandlimited white noise inputs. These numerical schemes are commonly referred to as autoregressive moving average (ARMA) and autoregressive (AR) algorithms (see e.g. Refs. [8–13]); they offer the most versatility in seeking an efficient recursive scheme for simulating power spectral compatible realizations of the random process. All methods for generating multivariate (one-dimensional) stochastic processes, require the decomposition of * Corresponding author. ⫹39-91-656-8111; fax: ⫹39-91-656-8407.

the Power Spectral matrix as a product of two matrices, the Choleski decomposition procedure is usually adopted. Nevertheless, in order to simulate a multivariate (onedimensional) vector processes with AR or ARMA or wave superposition models, the Power Spectral Density (PSD) matrix is decomposed by using Cholesky’s method which results in the product of two frequency dependent triangular matrices. As a result, each new component of the multivariate process contains an additional summation with respect to the previous one in order to account for the cross-correlations between time histories in points with different location in space. It follows that the computational effort increases with the law n…n ⫹ 1†=2; n being the number of components of the vector process; then, for very large number of components the generation of a multivariate wind field becomes prohibitive. In order to simulate multivariate stochastic processes, a very interesting way has been proposed by Li and Kareem [14]; it mainly consists of decomposing the multivariate stochastic vector process as a summation of fully coherent vectors, incoherent with each other. This decomposition is not unique and can be performed using different ways for decomposing the PSD matrix. Here the decomposition of the PSD matrix in the (frequency dependent) orthonormal basis of the eigenvectors of the PSD matrix itself is proposed. In this base at each frequency, the multivariate wind process is decomposed in independent fully coherent elementary vector processes.

0266-8920/01/$ - see front matter 䉷 2000 Elsevier Science Ltd. All rights reserved. PII: S0266-892 0(99)00032-6

2

M. Di Paola, I. Gullo / Probabilistic Engineering Mechanics 16 (2001) 1–10

The application to multivariate wind velocity is discussed in detail showing that for this special process only the first few spectral modes exhibit significant power. It follows that one can retain in the analysis only a reduced number of independent wave vectors drastically reducing the computational effort for both computation of eigen-properties and generation of the wind velocity vectors for every kind of chosen generation method. The physical significance of eigenvalues and eigenvectors of the PSD matrix of wind vector velocity is also discussed, showing that the eigenvectors behave like structural mode shapes, slowly varying with the frequency and then are called blowing mode shapes of wind velocity field. Moreover, it will be shown how the use of a standard finite element procedure for approximating the frequency dependent eigenvectors allows to obtain a generation procedure that requires a limited number of univariate coherent processes and this number is independent of the number of components of the process. A numerical example is given for the purpose of demonstrating the efficiency of the proposed method. 2. Stochastic modelling of multivariate wind field In this section the longitudinal wind field velocity is introduced. Let x, y, z be a point in space, z is the height from the ground and x is assumed to be the along-wind direction. Then velocity field V …x; y; z; t† can be written as  ⫹ V…x; y; z; t† V …x; y; z; t† ˆ V…z†

…1†

 is the mean value, and V…x; y; z; t† is the fluctuwhere V…z† ating component of the wind speed, assumed to be a normal zero mean stationary stochastic process. In the following,  is represented by the well-known logathe mean value V…z† rithmic law  ˆ 1 uⴱ ln …z=z0 † V…z† k

…2†

in which k ˆ 0:4 is the Von Karman’s constant, uⴱ is the shear velocity and z0 …m† is the roughness length. The longitudinal fluctuating components of two wind processes in points located in space are fully characterized from a probabilistic point of view by the cross-spectral density function SVr Vs …v†; Vj ˆ V…xj ; yj ; zj ; t† … j ˆ r; s†: The imaginary part was found to be negligible in comparison with the real part [15] while the latter (co-spectrum) can be expressed in the form q …3† SVr Vs …v† ˆ SVr Vr …v†SVs Vs …v† exp …⫺frs …v†† where frs …v† is given as [16] q 兩v兩 Cy2 …yr ⫺ ys †2 ⫹ Cz2 …zr ⫺ zs †2 frs …v† ˆ ⫹  s ††  r † ⫹ V…z 2p…V…z

…4†

the decay coefficients Cy and Cz being determined experi-

mentally. The functions SVr Vr …v†; SVs Vs …v† in Eq. (3) are the PSDF of the processes Vj …t† … j ˆ r; s†: Different expressions have been proposed in the literature, but in the following, the Solari’s spectrum [17] is adopted, that SVj Vj …v; zj † ˆ

6:868s2V f …zj; v†LV =zj …v=2p†‰1 ⫹ 10:302f …zj; v†LV =zj Š5=3

…5†

 j † is the Monin coordinate, LV in which f …zj ; v† ˆ vzj =2pV…z is the integral length scale of turbulence and s2V is the variance of the longitudinal component of the velocity fluctuations defined by

s2V ˆ b…zj †u2ⴱ

…6†

b…zj † being a non-dimensional function. Putting V1 ; V2 ; …; Vn into an n-vector, in the following denoted as V, the complete characterization of the n-variate (onedimensional) vector V, from a probabilistic point of view, can be made by the PSD matrix SV …v† 3 2 SV1 V1 …v† SV1 V2 …v† … SV1 Vn …v† 7 6 6 SV V …v† SV V …v† … SV V …v† 7 7 6 1 2 2 2 2 n 7 6 …7† SV …v† ˆ 6 7 .. .. .. 7 6 … 7 6 . . . 5 4 SV1 Vn …v† SV2 Vn …v† … SVn Vn …v† This matrix is real (because the quadrature spectrum has been neglected), symmetric and positive definite at each frequency v . It will be noted that the multivariate wind velocity vector V…t† is not coherent; in fact, the coherence function gVr Vs …v† between two points in space, is given as 兩SVr Vs …v†兩 ˆ e⫺frs …v† ⬍ 1; gVr Vs …v† ˆ p SVr Vr …v†SVs Vs …v†

᭙r 苷 s …8†

From this expression one can observe that the coherence function decreases by increasing the inter-distance between the location points r and s and decreases by increasing the frequency v . The n-variate stochastic vector process V…t† can be decomposed into a summation of n n-variate fully coherent normal vectors Y j …t† independent of each other, that is: V…t† ˆ

n X

Y j …t†

…9†

jˆ1

This can be done in different ways; as an example wellargued papers [10–14] propose to decompose the PSD matrix as a product of a frequency dependent matrix of order n × m and its transpose. Here the PSD matrix is decomposed into the basis of the eigenvectors of the PSD matrix itself. In order to do this, let C…v† be the eigenmatrix of S V …v†; whose columns are the eigenvectors (real and orthogonal) of SV …v†; normalized with respect to the identity matrix. Then

M. Di Paola, I. Gullo / Probabilistic Engineering Mechanics 16 (2001) 1–10

the following relationships hold

C …v†S V …v†C…v† ˆ L…v†

…10a†

C T …v†C…v† ˆ I

…10b†

T

I being the …n × n† identity matrix and L…v† being a diagonal matrix whose diagonal elements are the eigenvalues L j …v† …j ˆ 1; 2; …; n† associated with c j …v† ( jth column of the matrix C…v†). The real vectors Y j …t† in Eq. (9) can be written in the form q Z∞ c j …v† Lj …v† e ivt dBj …v† …11† Y j …t† ˆ ⫺∞

where Bj …v†is a zero mean normal complex process having orthogonal increments, that is E‰dBj …v†Š ˆ 0;

3

eigenvectors the appeal in using the proposed procedure will be evident as it will be shown in the next section. 3. Eigen-properties of the PSD matrix of the wind velocities In order to understand the physical meaning of eigenvectors and eigenvalues of the matrix SV …v†; two examples are examined. The first one is the case of two points at the same level. In this case the matrix SV …v† reduces to " # 1 e⫺f12 …v† …18† SV …v† ˆ SVV …v† e⫺f12 …v† 1 the eigenvalues of this matrix are

dBj …v† ˆ dBⴱj …v†;

…12†

L1 …v† ˆ SVV …v†…1 ⫹ e⫺f12 …v† †

…19a†

where the star denotes complex conjugate and dpq is the Kronecker delta …dpq ˆ 0 if p 苷 q; dpq ˆ 1 if p ˆ q†: The PSD matrix of Y j …t†; is simply given as

L2 …v† ˆ SVV …v†…1 ⫺ e⫺f12 …v† †

…19b†

E‰dBj …vr †dBⴱk …vs †Š ˆ dvr vs djk dvr

SY j …v† ˆ c j …v†c jT …v†Lj …v†

…13†

and hence n X

SV …v† ˆ

SY j …v† ˆ C…v†L…v†C…v† T

…14†

jˆ1

that is, in virtue of Eq. (10b), it coincides with Eq. (10a). The discretized version of Eq. (11) can be used for simulation purposes, that is Y j …t† ˆ

N X

q c j …vk † Lj …vk †Dv e ivk t P…j† k

…15†

kˆ⫺N

are zero mean where NDv is the cut-off frequency and normal complex random variable obeying the following orthogonality condition: …s† E‰P…j† k Pr Š ˆ djs dkr ;

ⴱ

…s† P…s† r ˆ P⫺r

…16†

Alternatively Eq. (15) can be rewritten in the real form Y j …t† ˆ 2

N X

1 c 2 ˆ p 2

"

#

1 ⫺1

…20b†

that is for a bivariate wind field velocity the spectral modal matrix C is independent of v . It follows that, according to Eqs. (9) and (11), one can write Z∞ p L1 …v† e ivt dB1 …v† V…t† ˆ c 1 ⫺∞

P…j† k

ⴱ

the corresponding eigenvectors, normalized with respect to the identity matrix, are " # 1 1 c 1 ˆ p …20a† 2 1

q c j …vk † Lj …v†Dv …cos vk tRk…j† ⫺ sin vk tIk…j† †

kˆ1

…17† …s† …s† R…s† r and Ir being the real and the imaginary part of Pr ; …s† …s† …s† …s† …s† …Pr ˆ R…r† ⫹ I…r† †; respectively, that is Rr ; Ir are independent, normal zero mean random variables having variance 1/2. At this stage the appeal of the generation of the wind velocity vector by decomposing the PSD matrix into the basis of the eigenvectors is not so evident with respect to the decomposition of the PSD matrix into the product of two triangular matrices using the Choleski method. However, on the basis of some physical meanings of eigenvalues and

⫹ c2

Z∞ p L…v† e ivt dB2 …v† ⫺∞

…21†

Remembering Eq. (9) one can think of the summation in Eq. (21) as two independent processes whose PSD are L1 …v† and L2 …v†; respectively; Then Eq. (21) can be rewritten in the form V…t† ˆ c 1 W1 …t† ⫹ c 2 W2 …t†

…22†

From Eq. (22) one can recognize that for a bivariate wind field velocity one can first generate two independent (scalar) processes W1 …t† and W2 …t† whose PSD are given in Eqs. (19a) and (19b), respectively, and then using Eq. (21), the vector V…t† is easily generated as the contribution of the two fully coherent vectors c 1 W1 …t†and c 2 W2 …t†: In Fig. 1a the PSD functions SVV …v† and the cross PSD  10 † ˆ 30 m=s; SVV …v† e ⫺f12 …v† are plotted versus v , for V…z z1 ˆ z2 ˆ 5 m; Cz ˆ 10; z0 ˆ 0:25 m and y1 ˆ y2 : In Fig. 1b the two functions L1 …v† and L2 …v† given in Eqs. (19a) and (19b) are plotted. From Fig. 1b one can see that the first eigenvector, having

4

M. Di Paola, I. Gullo / Probabilistic Engineering Mechanics 16 (2001) 1–10

Fig. 1. Power spectra of a bivariate wind field velocity. (a) Elements of the frequency dependent PSD matrix. (b) Frequency dependent eigenvalues L1 …v†; L2 …v†:

higher power at low frequency and shape c 1, is the most important, while the second, having shape c 2 at low frequency exhibits comparatively small power. Moreover at high frequency both W 1 …t†and W2 …t† have the same power. From this very simple example one can make some further considerations.

Let two anemometers be in two point locations in the plane y, z (see Fig. 2), as an example at the same level z. Let D be the inter-distance between the two points. By virtue of Eq. (22), one can decompose the wind field V1, V2 into two independent vector processes c 1 W1 …t† and c 2 W2 …t†, the first one (visualized in Fig. 2b) with the higher

Fig. 2. Decomposition of a bivariate wind velocities. (a) Wind action on two anemometers. (b) First blowing mode shape of wind velocity. (c) Second blowing shape of wind velocity.

M. Di Paola, I. Gullo / Probabilistic Engineering Mechanics 16 (2001) 1–10

5

Fig. 3. Elements of PSD matrix for a six-variate wind field velocity: frequency dependent eigenvalues Lj …v† for a six-variate wind field velocity.

power has for a resultant a vector located at the midpoint D=2; while the second having the smaller power (visualized in Fig. 2c), constitutes a moment about the z axis. If the inter-distance D decreases then, for small values of v , L2 …v† ! 0 and L1 …v† ! 2SVV …v†: On the contrary, if D increases, the two independent processes tend to have the same power since e ⫺f12 …v† vanishes. Moreover from this example one clearly realizes the physical significance of eigenvectors and eigenvalues. The latter are the power of the two independent processes W1 …t† and W2 …t†; while the eigenvectors are mode shapes associated with the wind field velocity. As in the case of a structural vibration decomposed as a sum of independent structural mode shapes, the bivariate wind field blows as a sum of two independent totally coherent processes associated with blowing shapes and therefore they will be hereafter called blowing mode shapes of the wind velocity field. Now, in order to show that the above mentioned properties remain almost unchanged increasing the dimension of vector V…t†; when a second example with a six-dimensional wind field is examined. The points 1–6 are located at different levels having inter-distance of 5 m. Also in this case  10 † ˆ 30 m=s: Cz ˆ 10; V…z In Fig. 3 the corresponding eigenvalues Lj …v†; j ˆ 1; 2; …; 6 are also plotted. The components of the matrix L…v† are ordered in such a way that L1 …v† ⬎ L2 …v† ⬎ … ⬎ L6 …v†: These figures show a significant power of the first blowing mode shape c 1 …v† while the other ones are smaller and smaller at low frequencies in comparison with the first eigenvalue. It follows that one can evaluate only the first few eigenvectors, those associated with the higher power Lj …v†; j ˆ 1; 2; …; m …m p n†: In Fig. 4a and b the eigenvectors c 1 and c 6 are plotted versus v . From this figure one can see that the various components of the eigenvector are very regular functions. It follows that using the simple or simultaneous vector iteration methods (see e.g. Refs. [18,19]) for the computations of

the eigenvectors c 1 …v†; if one uses as the initial vector in the iteration c j …vk⫺1 † as a first attempt, the actual eigenvector c j …vk † can be evaluated by means of very few iterations. In Fig. 5a and b the first and the sixth blowing mode shapes of the wind field are plotted for different values of v , showing a surprising similarity with the corresponding structural modes. It follows that for a slender structure the contribution on the response of each natural structural mode will be denominated by the corresponding blowing mode shape. This explains the efficiency of the equivalent wind spectrum technique [20] for the analysis of structures exposed to wind actions. The Solari method mainly consists of adopting a so-called “equivalent wind structure” defined as a multi-dimensional coherent vector V…t†; whose PSD matrix approximates optimally the corresponding PSD matrix related to the actual turbulence configuration.

4. Polynomial approximation of eigenvectors Following what has been stated in the previous section, the multivariate wind velocity process can be decomposed into a summation of independent fully coherent multivariate vectors. It follows that one can generate at each time instant ts and vector Y j …ts † … j ˆ 1; 2; …; n† and then by summing the contribution of the various vectors Y j …ts † one gets V…ts †. The key to the appeal of the decomposition of the PSD matrix in a basis of eigenvectors lies in the fact that, according to Eq. (11) in the vector Y j …t†; the orthogonal increments dBj …v† are scalar ones. In order to take full advantage of this circumstance for generation purposes by an autoregressive (AR) model, another fundamental step has to be made. Let us define the frequency domain ‰v0 ; vc Š; v0 and v c being the appropriate lower and upper cut-off frequencies, subdivided into M parts v0 ⬅ V 0; V 1 ; …; V m ⬅ vc : Then, we can approximate the eigenvector C j …v† in the generic interval ‰V s⫺1 ; V s Š in a polynomial form of a fixed order; in the

6

M. Di Paola, I. Gullo / Probabilistic Engineering Mechanics 16 (2001) 1–10

Fig. 4. Various components of frequency dependent eigenvectors of six-variate wind field velocity; solid line exact; dashed line approximated by third-order polynomial with M ˆ 1: (a) First eigenvector. (b) Sixth eigenvector.

following a third-order polynomial is adopted because it leads to a good representation of the eigenvector as it will be shown later. Then one can adopt for the sth interval the following approximate jth eigenvector.

c j…s† …v† ˆ N…s† j l…v†;

V s⫺1 ⱕ v ⱕ V s

…23†

where l T …v† ˆ ‰1 v v2 v3 Š: The matrix N…s† j (of order n × 4) can be obtained by imposing the continuity at boundaries v ⬅ V s⫺1 ; and v ⬅ V s and hence one can write  j …V s⫺1 ; V s †Lj⫺1 …V s⫺1 ; V s † Nj…s† ˆ C

…24†

 j …V s⫺1 ; V s † ˆ ‰c j …V s⫺1 †c 0j …V s⫺1 †c j …V s †c 0j …V s †Š C

…26†

0

where the prime denotes differentiation with respect to v . By means of Eq. (23), Eq. (11) can be rewritten in the form Y j …t† ˆ

M X

Nj…s†

ZV s Vs ⫺ 1

sˆ1

q l…v† Lj …v† e ivt dB…s† j …v†

…27†

where the orthogonal increments stochastic processes dB…s† j …v† obey the following relationship ⴱ

…u† E‰dB…s† j …vp †dBj …vq †Š ˆ dvp vq dsu dvp

…28†

Hence Eq. (27) can be rewritten in the form

where 2

1

6 6V 6 s⫺1 Lj …V s⫺1 ; V s Š ˆ 6 6 2 6 V s⫺1 4 3 V s⫺1

0

1

1

Vs

2V s⫺1

V s2

2 3V s⫺1

V s3

0

3

7 1 7 7 7 7 2V s 7 5 3V s2

Y j …t† ˆ

M X

…s† N…s† j U j …t†

…29†

sˆ1

…25†

where U…s† j …t† …s ˆ 1; 2; …; M† are fully coherent four-variate processes independent of each other, by virtue of Eq. (28). This is a very remarkable result because in order to generate

M. Di Paola, I. Gullo / Probabilistic Engineering Mechanics 16 (2001) 1–10

7

powers Lj …v† in the frequency interval V s⫺1 ; V s and zeros elsewhere. The other components are related to the time derivatives of a such process; as an example the second …s† …t†; is a process having power v2 Lj …v† in component Uj;2 the frequency interval V s⫺1 ; V s and hence, because …s† …t† ˆ U_ …s† dBj…s† …v† is a scalar process, then Uj;2 j;1 …t†: It will be noted that for M ⬎ 1 we have to generate samples U…s† j …t† having different Nyquist frequencies, then we have to generate the samples of U…s† j …t†, having different time step and since Dt…s† ˆ p=V j ; the following relationship holds Dt …1† ⬎ Dt…2† ⬎ … ⬎ Dt…M†

…31†

and then we can put V j ˆ jV s : In this way, starting from the smallest Dt, that is Dt …M† ; the summation of the various …s† N…s† j U j …t† can be easily performed. Summarizing, by means of the eigenvectors of the PSD matrix the incoherent multivariate process V…t† is decomposed as the summation of n fully coherent independent n-variate vectors Y j …t† according to Eq. (9). By means of piecewise polynomial approximation of the eigenvectors, each fully coherent vector Y j …t† is decomposed as the summation of M independent fully coherent four-variate vectors according to Eq. (29). In the next section we will take full advantage of the aforementioned properties for generation purposes by digital AR models. 5. AR generation In order to evaluate the generic component Y j …t† of the vector V…t† we have to generate the M fully coherent vectors U…s† j …t† defined in Eq. (30). Using the standard generation via AR model, we can write …s† …tk † ˆ Uj;1

p X

…s† …s† …s† a…s† j;u Uj;1 …tk⫺u † ⫹ sj;1 Wj …tk †

…32a†

…s† …s† …s† b…s† j;u Uj;2 …tk⫺u † ⫹ sj;2 Wj …tk †

…32b†

…s† …s† …s† c…s† j;u Uj;3 …tk⫺u † ⫹ sj;3 Wj …tk †

…32c†

…s† …s† …s† dj;u Uj;4 …tk⫺u † ⫹ s…s† j;4 Wj …tk †

…32d†

uˆ1

Fig. 5. Blowing mode shapes of the six-variate wind velocity. (a) First mode. (b) Sixth mode.

…s† …tk † ˆ Uj;2

p X uˆ1

the entire vector Y j …t† we only need to generate …4 × M† independent univariate processes, instead of an n-variate vector process, each of them given in the form 0 2 31 1 B 6 7C Bq6 v 7C ZV s B 6 7C …s† 7C e i vt B Lj …v†=26 U…s† j …t† ˆ B 6 2 7C dBj …v†; Vs ⫺ 1 B 6 v 7C @ 4 5A v3 …s ˆ 1; 2; …; M†

(30)

The four components of the vector U…s† j …t† are characterized by the same Nyquist frequency Dt…s† ˆ p=V s : Moreover, the …s† …t† is a process having first component of this vector Uj;1

…s† …tk † ˆ Uj;3

p X uˆ1

…s† …tk † ˆ Uj;4

p X uˆ1

…s† …s† where a…s† j;u ; …; dj;u are the parameters of the AR model, sj;r …s† …r ˆ 1; …; 4† are the variances of the input and Wj …tk † are the normal random variables with zero mean and unit variance, and p is the number of parameters of the filter. It is worth noting that in Eqs. (32a)–(32d) the input is the same for all components of U…s† j because in Eq. (30) the input is a …s† …s† scalar quantity. The parameters a…s† j;u ; …; dj;u and sj;r ; …r ˆ 1; …; 4† can be evaluated by the usual Yule–Walker scheme

8

M. Di Paola, I. Gullo / Probabilistic Engineering Mechanics 16 (2001) 1–10

[21,22], as an example for the first component of U…s† j …t† one can write RU…s† …tk ⫺ t` † ˆ j;1

p X

a…s† j;u RU …s† …tk⫺u ⫺ t` † j;1

uˆ1

t` ˆ tk⫺1 ; tk⫺2 ; …tk⫺p

(33)

…s† …t† that where RU …s† …tk ⫺ t` † is the correlation function of Uj;1 j;1 is clearly evaluated as the Fourier transform of Lj …v† in the interval V s⫺1 ; V s ; that is

RU …s† …t† ˆ

ZV s Vs ⫺ 1

j;1

Lj …v† e ivt dv

…34†

Once the filter parameters a…s† j;u are evaluated by solving the linear system (33), the standard parameter s…s† j;1 of the input can be evaluated by the relationship RU…s† …0† ˆ j;1

p X uˆ1

…s† a…s† j;u RU …s† …tk ⫺ tk⫺s † ⫹ sj;1 j;1

2

…35†

For the other components of the vector U…s† j …t† one can proceed in a similar way. It will be noted that the property already outlined in the previous section, that is that the components of U…s† j …t† are related to each other by differentiation is not taken into account because of the difficulty in the definition of derivation of the random variable Wj…s† …t†: By performing the AR generation by the outlined procedure we have to solve 2…M ⫹ 1† eigenproblems at the end of the interval V s⫺1; V s because we need the derivatives of the eigenvectors at the end of the intervals, and this can be done by evaluating the eigenproblems at the frequencies V 0 ⫹ dV; V 1 ⫺ dV; …; V M ⫺ dV; dV being a very small frequency. Moreover 4M Fourier transforms are needed for the generation of each vector Y j …t† for the evaluation of the correlation function by means of Eq. (34). At least for each vector Y j …t† we also need the solution of 4M linear systems of p equations of the p unknowns in Eq. (33) for evaluating the components of the elementary vectors Uj…s† …t†: Once the vector U…s† j …t† …s ˆ 1; …M† is generated, the entire vector Y j …t† can be evaluated by using Eq. (29) and by using Eq. (9) the entire V…t† vector can be easily computed. 6. Example In this section, a numerical example is given for the purpose of demonstrating the efficacy of the method presented in the previous sections. Let V…t† be a six-dimensional wind field obtained by locating six points at different levels having inter-distance of 5 m. Wind profile and turbulence were simulated using Eqs. (2) and (5), respectively. Wind histories were generated using the following data: T ˆ 150 s; time step Dt ˆ 0:157 s; integration limits of the power spectra v0 ˆ 0:08 rad=s and vc ˆ 20 rad=s;

roughness length z0 ˆ 0:25 m; and reference mean wind speed at the 10 m level V…10† ˆ 30 m=s: The order of autoregression assumed is p ˆ 4 while the number of intervals for approximating the eigenvectors of the matrix SV …v† was M ˆ 1: In Fig. 3 the corresponding eigenvalues Lj …v† … j ˆ 1; 2; …; 6† are plotted; from this figure, it can be seen that the first eigenvalue L1 …v† shows a significant power while the others become are smaller and smaller in comparison to the first eigenvalue. That means that one can truncate the analysis including only very few eigenvalues and eigenvectors; then, Eq. (9) can be rewritten in the form V…t† ˆ

m X

Y j …t†;

mpn

…36†

jˆ1

where Y j …t† can be evaluated by Eq. (29). Using the truncation procedure outlined above, the computational effort can be reduced in a drastic way. In Fig. 4a and b the exact and approximated eigenvectors c 1 …v† and c 6 …v† are plotted. From these figures, one can see that using the polynomial approximation and selecting one interval …M ˆ 1†; the approximated eigenvectors cannot be distinguished from the exact one. The esavariate wind field was first simulated taking into account all of the six eigenvectors and only the first one c 1 …t† and M ˆ 1: The results in terms of correlation are reported in Fig. 6a and b. From these figures one can see that the truncation procedure with m ˆ 1 in Eq. (36) and M ˆ 1 in Eq. (29) is enough for describing the multivariate processes. 7. Conclusions In this paper, the generation procedure of multivariate wind field velocity stochastic processes has been discussed in detail by decomposing the Power Spectral Density (PSD) matrix in the frequency dependent orthonormal base of the eigenvectors of the PSD matrix itself. In this way the vector process (usually incoherent) is reconducted to a sum of elementary vectors fully coherent and independent of each other in virtue of the orthogonality of the eigenvectors. The drawback of the proposed generation procedure is the evaluation of the eigen-properties of the PSD matrix at each frequency. However the application of such procedure for generating multivariate wind velocity shows some interesting properties: (i) only very few spectral modes exhibit significant power; (ii) each component of the eigenvectors are very regular functions, and this fact is independent of the analytical model assumed for the PSD functions. In view of the first property, one can truncate the spectral modal matrix retaining only the first few modes shapes, that is those having significant power, say m …m p n†; in this way, the computational effort can be reduced in a drastic way. Because of the second property, by using a standard

M. Di Paola, I. Gullo / Probabilistic Engineering Mechanics 16 (2001) 1–10

9

Fig. 6. Comparison between correlation function obtained by using only the first eigenvector c 1 …v† and all the eigenvectors c j …v† … j ˆ 1; …; 6† with the target correlation function. (a) Auto-correlation function. (b) Cross-correlation function

finite element procedure for approximating the eigenvectors, the generation procedure, via standard AR model, requires a limited number of samples of univariate coherent processes for describing the multivariate wind field and this number is independent of the number of the components of the process. In order to generate the totally coherent vectors we need to approximate the eigenvectors in a polynomial form. This can be done by using the standard finite element procedure. In virtue of the second property (regularity of eigenvectors), we only need a small number of subdivision frequencies for approximating eigenvectors. Very accurate results have been obtained by taking M ˆ 1 into the entire frequency range …v0 –vc † as has been demonstrated in the example. The physical significance of eigenvectors and eigenvalues has been also discussed in detail showing that the wind field behaves exactly like a structure, so that combining the concepts of structural modal analysis and blowing

mode shapes of the wind field implies that the analysis of MDOF systems becomes simple and immediate. In this frame, the proposed method becomes very competitive with the classical AR procedure based on the decomposition of the PSD matrix into triangular matrices by Cholesky’s method. Acknowledgements Comments and suggestions of Prof. Giovanni Solari are fully acknowledged. References [1] Shinozuka M. Simulation of multivariate and multidimensional random processes. J Acoust Soc Am 1970;49(1):357–67.

10

M. Di Paola, I. Gullo / Probabilistic Engineering Mechanics 16 (2001) 1–10

[2] Shinozuka M, Jan CM. Digital simulation of random processes and its applications. J Sound Vibr 1972;25(1):111–28. [3] Shinozuka M, Yun C-B, Seya H. Stochastic methods in wind engineering. J Wind Engng Ind Aerodynam 1990;36:9. [4] Shinozuka M, Deodatis G. Simulation of stochastic processes by spectral representation. Appl Mech Rev, ASME 1991;44(4):191–204. [5] Li Y, Kareem A. Simulation of multivariate random processes: hybrid DFT and digital filtering approach. J Engng Mech, ASCE 1993;119(5):9. [6] Grigoriu M. Simulation of non-stationary Gaussian processes by random trigonometric polynomials. J Engng Mech, ASCE 1993;119(2):9. [7] Yamazaki F, Shinozuka M. Simulation of stochastic fields by statistical preconclutioning. J Engng Mech, ASCE 1990;116(2):268–87. [8] Deodatis G, Shinozuka M. Autoregressive model for non-stationary stochastic processes. J Engng Mech, ASCE 1988;114(11):1995– 2012. [9] Deodatis G. Simulation of multivariate stochastic processes. In: Spanos PD, editor. Computational stochastic mechanics, Amsterdam: Balkema, 1995. pp. 297–305. [10] Li Y, Kareem A. Stochastic decomposition application to probabilistic dynamics. J Engng Mech, ASCE 1993;121(10):9. [11] Spanos PD, Mignolet M. Z-transform modelling of P-M wave spectrum. J Engng Mech, ASCE 1986;112(8):745–59.

[12] Spanos PD, Mignolet MP. Recursive simulation of stationary multivariate random processes: Part II. J Appl Mech 1987;54:681–7. [13] Nagaguma T, Deodatis G, Shinozuka M. ARMA model for two dimensional processes. J Engng Mech, ASCE 1987;113(2):234–51. [14] Li Y, Kareem A. Simulation of multivariate nonstationary random processes by FFT. J Engng Mech, ASCE 1991;117(5):9. [15] Simiu E, Scanlan RH. Wind effects on structures. New York: Wiley, 1978. [16] Davenport AG. The spectrum of horizontal gustiness near the ground in high winds. Quart J R Mater Soc 1961;87:194–211. [17] Solari G. Gust buffeting. I. Peak wind velocity and equivalent pressure. J Struct Engng, ASCE 1993;119(2):365–82. [18] Argyris J, Mlejnek HP. Dynamics of structures. Amsterdam: Elsevier, 1991. [19] Sehemi NS. Large order structural eigenanalysis for finite element systems. Sussex, UK: Ellis Horwood, 1989. [20] Solari G. Equivalent wind spectrum technique: theory and applications. J Struct Engng, ASCE 1988;114(6):1303–23. [21] Yule GU. On a method of investigating periodicities in disturbed time series, with special reference to Walfer’s sunspot numbers. Phil Trans R Soc, Series A 1927;226:267–98. [22] Walker G. On periodicity in series of related terms. Proc R Soc London, Series A 1931;131:518–32.