Spline-Interpolation-Based FFT Approach to Fast Simulation of Multivariate Stochastic Processes

Abstract

The spline-interpolation-based fast Fourier transform (FFT) algorithm, designated as the SFFT algorithm, is proposed in the present paper to further enhance the computational speed of simulating the multivariate stochastic processes. The proposed SFFT algorithm first introduces the spline interpolation technique to reduce the number of the Cholesky decomposition of a spectral density matrix and subsequently uses the FFT algorithm to further enhance the computational speed. In order to highlight the superiority of the SFFT algorithm, the simulations of the multivariate stationary longitudinal wind velocity fluctuations have been carried out, respectively, with resorting to the SFFT-based and FFT-based spectral representation SR methods, taking into consideration that the elements of cross-power spectral density matrix are the complex values. The numerical simulation results show that though introducing the spline interpolation approximation in decomposing the cross-power spectral density matrix, the SFFT algorithm can achieve the results without a loss of precision with reference to the FFT algorithm. In comparison with the FFT algorithm, the SFFT algorithm provides much higher computational efficiency. Likewise, the superiority of the SFFT algorithm is becoming more remarkable with the dividing number of frequency, the number of samples, and the time length of samples going up.

1. Introduction

Monte Carlo technique has widely been employed for simulating the stochastic processes which are either one-dimensional or multidimensional, univariate or multivariate, homogeneous or nonhomogeneous, stationary or nonstationary, and Gaussian or non-Gaussian. The Monte Carlo simulation methods are able to generate the sample functions that accurately provide the probabilistic characteristics of the corresponding stochastic processes. For the simulation of stochastic processes, the following approaches are now available: (1) autoregressive (AR) method, such as Mignolet and Spanos [1, 2], Iannuzzi and Spinelli [3], Deodatis and Shinozuka [4], and Novak et al. [5]; (2) autoregressive moving average (ARMA) method, for example, Gersch and Yonemoto [6], Kozin [7], Kareem and Li [8], and Rossi et al. [9]; (3) spectral representation (SR) method, for instance, Shinozuka and Jan [10], Grigoriu [11], Deodatis [12], Grigoriu [13], and Chen and Letchford [14].

It is known that among foregoing simulation methods, the SR method has very high demands on both the computer memory and speed. Notwithstanding this, the SR does not have the problem of model selection, as does the AR and ARMA methods. Likewise, it is easy to implement and has high accuracy. Hence, the SR method has yet been receiving increasing attention in simulating the multivariate stochastic processes with the target cross power spectral density (CPSD) matrices. The primary concept of the SR method dates back to Rice [15], Goto and Toki [16], and Borgman [17]. For the multidimensional, multivariate, and nonstationary cases, Shinozuka [18, 19] first established the SR approach. Subsequently, Deodatis and Shinozuka [20] extended the application of the SR method to the simulation of stochastic waves. Li and Kareem [21] developed a hybrid discrete Fourier transform and digital filtering approach to simulate the multivariate random process. Likewise, Grigoriu [11] compared two different SR models. But, it is worth pointing out that the generated sample functions with resorting to the early algorithm of the SR method by Shinozuka and Jan [10] are not ergodic. In order to generate the ergodic sample functions, several attempts have been made via modifying the existing algorithms. Shinozuka et al. [22] introduced the idea of double-indexing frequencies. But, the obtained sample functions based on the proposed formula are still not ergodic. Deodatis [12] further extended the SR method to simulate the multivariate ergodic stochastic processes. Likewise, the capabilities and efficiency of the proposed algorithm were demonstrated in detail using a one-dimensional trivariate process as an example. On the other hand, the computational efficiency of the early algorithm of the SR method is yet low. In order to cope with this issue, the fast Fourier transform (FFT) algorithm was introduced into the SR method. Yang [23, 24] showed that the FFT algorithm can remarkably enhance the computational efficiency of the SR method and proposed a formula to simulate the random envelop processes. Shinozuka [25] extended the application of the FFT algorithm to the multidimensional cases. Wittig and Sinha [26] showed that the SR method combined with the FFT algorithm seems to be more than an order of magnitude faster than other simulation methods. Further, Li and Kareem [27] used the FFT-based approach to simulate the multivariate nonstationary Gaussian random processes with the prescribed evolutionary spectral description.

Likewise, the multicorrelated stationary random processes, such as the wind velocity or pressure fluctuations, on structures can be transformed into a set of subprocesses through diagonalizing their covariance or CPSD matrices with resorting to either the Cholesky (lower or upper triangular) or eigenvector decomposition. The eigenvector decomposition offers physically meaningful insight into the process as each eigenvector (eigenmode) may be characterized on the basis of its spatial distributions. Theoretically, the eigenvector decomposition is based on the Karhunen-Loeve expansion, consequently referred also to as the proper orthogonal decomposition (POD) [28]. The POD approach has been widely utilized to reduce the dimensions or variables in the large scale systems to enhance the computational speed, such as Holmes et al. [29], Di Paola and Gullo [30], Rathinam and Petzold [31], Chen and Kareem [32], and Solari and Carassale [33]. In [34], Carassale and Solari proposed some strategies such as interpolation of the factorizations of the CPSD matrix, cutoff of the spectral representation based on POD and parallelization of the code, which provide significant computational advantages enabling the simulation of large wind fields. Likewise, Ding et al. [35] also introduced interpolation approximation of the factorizations of the CPSD matrix and cutoff of the spectral representation based on POD to the SR algorithm for the simulation of wind velocity fields on large scale structures. The numerical advantage of the POD technique, akin to the modal analysis of structural dynamics, relies mainly on the reduced-order representation via the truncation of the higher eigenmodes associated with smaller eigenvalues. Of course, this reduced-order representation must warrant that the important characteristics of the stochastic processes and related quantities remain unchanged, or the modification resulting from the approximate representation is acceptable. However, the truncation of higher modes may not necessarily work effectively in the case of local responses, implying that there exists a possibility of underestimating the local wind loads and corresponding effects [36]. In the case of using the CPSD matrix-based POD technique, the similar observations made by Chen and Kareem [37] once again underscore the foregoing phenomenon.

In order to illuminate solely the effectiveness of interpolation in fast simulation of multivariate stochastic processes without POD technique, the spline function employing the nature of both the smoothness and continuance is utilized to reduce the number of the Cholesky decomposition of power spectral matrix by means of interpolation. The spline-interpolation-based FFT algorithm, designated as the SFFT algorithm, then is proposed in the present paper. Therefore, different from [34, 35], the main purpose of the present study is to evaluate the computational efficiency and accuracy of the SFFT-based SR method with reference to the FFT-based SR method through simulating the multivariate stationary longitudinal wind velocity fluctuations with the phase angles.

2. Simulation of Multivariate Stochastic Processes

According to Deodatis [12] and Shinozuka [19], the multivariate stationary stochastic process 𝑓^0𝑗(𝑡)(𝑗=1,2,…,𝑛) with the mean value equal to zero can be simulated by the following series:𝑓_𝑗(𝑡)=2√Δ𝜔^𝑛_𝑁_𝑚=1_𝑙=1||𝐻_𝑗𝑚𝜔_𝑚𝑙||cos𝜔_𝑚𝑙𝑡−𝜗_𝑗𝑚𝜔_𝑚𝑙+Φ_𝑚𝑙(𝑗=1,2,…,𝑛),(2.1) where 𝐇(𝜔) is the Cholesky decomposition of CPSD matrix 𝐒⁰(𝜔) (see the appendix) and a lower triangular matrix; 𝑁 is the sufficiently large dividing number of frequency; Δ𝜔=𝜔_up/𝑁 is the circular frequency increment; 𝜔_up refers to the upper cutoff circular frequency, with the condition that, when 𝜔>𝜔_up, the value of 𝐒⁰(𝜔) is trivial and negligible; Φ_𝑚𝑙 is the sequences of independent random phase angles, distributed uniformly over the interval [0,2𝜋]; the double-indexing frequency 𝜔_𝑚𝑙=𝑙Δ𝜔−((𝑁−𝑚)/𝑁)Δ𝜔=(𝑙−1)Δ𝜔+(𝑚/𝑁)Δ𝜔.

It is noted that the simulated stochastic process using the SR method is asymptotically Gaussian as 𝑁→∞ (Deodatis [12]). Shinozuka and Deodatis [38] provided rigorous derivations and elaborations about asymptotic Gaussian of the simulated stochastic process according to the central limit theorem. It can be considered approximately Gaussian for most practical applications if 𝑁 is greater than approximately 100 (Deodatis and Micaletti [39]).

In order to avoid aliasing, the time interval Δ𝑡 has to obey the condition as follows:Δ𝑡≤2𝜋2𝜔_up.(2.2) Likewise, the entire period of the sample functions in (2.1) can be calculated as follows:𝑇₀=2𝜋𝑛Δ𝜔=2𝜋𝜔𝑛_up𝑁.(2.3)

It has been proved by Deodatis [12] that the obtained results based on (2.1) possesses the ergodicity. Apparently, once the CPSD matrix is determined and the associated parameters, such as 𝑁,𝜔_up, and Δ𝜔, are properly chosen, the stationary one-dimensional multivariate Gaussian stochastic process can then be simulated quite well with resorting to (2.1).

3. Spline-Interpolation-Based FFT (SFFT) Algorithm

For the simulation of the multivariate stationary stochastic processes, the spline-interpolation-based FFT algorithm, referred to as the SFFT algorithm, is proposed in the present paper.

3.1. Spline Interpolation Technique

Since the Cholesky decomposition has to be conducted separately for each frequency 𝜔_𝑚𝑙, the number of performing the Cholesky decomposition with respect to (2.1) then is 𝑛×𝑁. The computational effort of the SR method is tremendous in the large scale system, though the FFT algorithm can significantly enhance the computational speed.

For the multivariate stationary stochastic processes with the phase angles, the CPSD matrix 𝐒⁰(𝜔) is the complex-valued matrix. Then, the element 𝐻_𝑗𝑚(𝜔)=|𝐻_𝑗𝑚(𝜔)|𝑒^{𝑖𝜗𝑗𝑚} of the obtained lower triangular matrix 𝐇(𝜔) based on the Cholesky decomposition is the complex value. It is noted that |𝐻_𝑗𝑚(𝜔)| varies continuously with regard to the circular frequency. Therefore, as long as |𝐻_𝑗𝑚(𝜔)| at some appropriate circular frequency points are calculated, |𝐻_𝑗𝑚(𝜔)| at other circular frequency points can then be obtained by the cubic spline interpolation. Since the spline function has the nature of both the smoothness and continuance, it has widely been utilized to interpolate and fit data in engineering.

The circular frequency interval [0,𝜔_𝑢] is evenly divided into 𝑟 subintervals by 𝑟−1 frequencies  𝜔₁,𝜔₂,…,𝜔_𝑟−1(0=𝜔₀<𝜔₁<𝜔₂<⋯<𝜔_𝑟−1<𝜔_𝑟=𝜔_𝑢). The corresponding |𝐻_𝑗𝑚(𝜔_𝑖)|(𝑖=0,1,2,…,𝑟) can then be calculated. It is assumed that 𝐻_𝑗𝑚(𝜔_𝑖)=|𝐻_𝑗𝑚(𝜔_𝑖)|(𝑖=0,1,2,…,𝑟) and second continuous derivative 𝐻^{𝑗𝑚}(𝜔_𝑖)=𝑃_𝑖(𝑖=0,1,2,…,𝑟), then 𝐻_𝑗𝑚(𝜔) at other frequency points can be obtained by the cubic spline interpolation. The flow scheme of building spline interpolation 𝐻_𝑗𝑚(𝜔) is shown in Figure 1.

Figure 1 

Flow scheme of building spline interpolation function 𝐻𝑗𝑚(𝜔).

In the circular frequency interval [𝜔_𝑖−1,𝜔_𝑖],𝐻_𝑗𝑚(𝜔) is expressed as cubic polynomial, and then 𝐻^{𝑗𝑚}(𝜔) is a linear function:𝐻^{𝑗𝑚}(𝜔)=𝑃_𝑖−1𝜔_𝑖−𝜔𝜔_𝑖−𝜔_𝑖−1+𝑃_𝑖𝜔−𝜔_𝑖−1𝜔_𝑖−𝜔_𝑖−1=𝑟𝜔_𝑢𝑃_𝑖−1𝜔_𝑖−𝜔+𝑃_𝑖𝜔−𝜔_𝑖−1.(3.1) Then, the continuous double integral of 𝐻^{𝑗𝑚}(𝜔) can be expressed as 𝐻_𝑗𝑚(𝜔)=𝑟6𝜔_𝑢𝜔^𝑖−𝜔³𝑃_𝑖−1+𝜔−𝜔_𝑖−1³𝑃_𝑖+𝐴^𝑖𝜔_𝑖−𝜔+𝐵_𝑖𝜔−𝜔_𝑖−1,(3.2) in which 𝐴_𝑖 and 𝐵_𝑖 are integral constants. According to the boundary conditions𝐻_𝑗𝑚𝜔_𝑖−1=𝜔^2𝑢6𝑟²𝑃_𝑖−1+𝐴_𝑖𝜔_𝑢𝐻_𝑗𝑚𝑟,𝜔_𝑖=𝜔^2𝑢6𝑟²𝑃_𝑖+𝐵_𝑖𝜔_𝑢𝑟,(3.3) the integral constants 𝐴_𝑖 and 𝐵_𝑖 can be expressed as𝐴_𝑖=𝐻_𝑗𝑚𝜔_𝑖−1−𝜔^2𝑢6𝑟²𝑃_𝑖−1𝜔𝑟_𝑢𝐵_𝑖,=𝐻_𝑗𝑚𝜔_𝑖−𝜔^2𝑢6𝑟²𝑃_𝑖𝜔𝑟_𝑢.(3.4) Substituting (3.4) into (3.2), 𝐻_𝑗𝑚(𝜔)=𝑟6𝜔_𝑢𝜔^𝑖+𝐻_𝑗𝑚−𝜔³𝑃_𝑖−1+𝜔−𝜔_𝑖−1³𝑃_𝑖𝜔_𝑖−1−𝜔^2𝑢6𝑟²𝑃_𝑖−1𝜔𝑟_𝑢𝜔_𝑖+𝐻−𝜔_𝑗𝑚𝜔_𝑖−𝜔^2𝑢6𝑟²𝑃_𝑖𝜔𝑟_𝑢𝜔−𝜔_𝑖−1𝜔_𝑖−1≤𝜔≤𝜔_𝑖;𝑖=1,2,3,…,𝑟.(3.5) Then, the derivative of 𝐻_𝑗𝑚(𝜔) can be written as𝐻^𝑗𝑚(𝜔)=𝑟2𝜔_𝑢𝜔−𝜔_𝑖−1²𝑃_𝑖−𝑟2𝜔_𝑢𝜔_𝑖−𝜔²𝑃_𝑖−1+𝜔^𝑢6𝑟𝑃^𝑖−1+𝑟−𝑃_𝑖𝜔_𝑢𝐻^𝑗𝑚𝜔_𝑖−𝐻_𝑗𝑚𝜔_𝑖−1𝜔_𝑖−1≤𝜔≤𝜔_𝑖;𝑖=1,2,3,…,𝑟.(3.6) In the interval [𝜔_𝑖−1,𝜔_𝑖],𝐻^𝑗𝑚𝜔_𝑖=𝜔^𝑢3𝑟𝑃^𝑖+𝜔^𝑢6𝑟𝑃^𝑖−1+𝑟𝜔_𝑢𝐻^𝑗𝑚𝜔_𝑖−𝐻_𝑗𝑚𝜔_𝑖−1,(3.7) and, in next interval [𝜔_𝑖,𝜔_𝑖+1],𝐻^𝑗𝑚𝜔_𝑖=−𝜔^𝑢3𝑟𝑃^𝑖−𝜔^𝑢6𝑟𝑃^𝑖+1+𝑟𝜔_𝑢𝐻^𝑗𝑚𝜔_𝑖+1−𝐻_𝑗𝑚𝜔_𝑖.(3.8) According to (3.7) and (3.8), we can obtain 𝑃_𝑖+1+4𝑃_𝑖+𝑃_𝑖−1=6𝑟²𝜔^2𝑢𝐻^𝑗𝑚𝜔_𝑖+1−2𝐻_𝑗𝑚𝜔_𝑖+𝐻_𝑗𝑚𝜔_𝑖−1(𝑖=1,2,3,…,𝑟−1).(3.9) In order to determine 𝑃_𝑖(𝑖=0,1,2,…,𝑟), the boundary conditions are assumed as 𝐻^𝑗𝑚(𝜔₀)=0 and 𝐻^𝑗𝑚(𝜔_𝑟)=0. Then, we can obtain the following functions:2𝑃₀+𝑃₁=6𝑟²𝜔^2𝑢𝐻^𝑗𝑚𝜔₁−𝐻_𝑗𝑚𝜔₀,(3.10)2𝑃_𝑟+𝑃_𝑟−1=6𝑟²𝜔^2𝑢𝐻^𝑗𝑚𝜔_𝑟−1−𝐻_𝑗𝑚𝜔_𝑟.(3.11) According to (3.9), (3.10) and (3.11), 𝑃_𝑖(𝑖=0,1,2,…,𝑟) can be determined through the following equation system2𝑃₀+𝑃₁=6𝑟²𝜔^2𝑢𝐻^𝑗𝑚𝑃₂𝜔₁−𝐻_𝑗𝑚𝜔₀,+4𝑃₁+𝑃₀=6𝑟²𝜔^2𝑢𝐻^𝑗𝑚𝑃₃𝜔₂−2𝐻_𝑗𝑚𝜔₁+𝐻_𝑗𝑚𝜔₀,+4𝑃₂+𝑃₁=6𝑟²𝜔^2𝑢𝐻^𝑗𝑚𝑃_𝑟𝜔⋮₃−2𝐻_𝑗𝑚𝜔₂+𝐻_𝑗𝑚𝜔₁,+4𝑃_𝑟−1+𝑃_𝑟−2=6𝑟²𝜔^2𝑢𝐻^𝑗𝑚2𝑃_𝑟𝜔_𝑟−2𝐻_𝑗𝑚𝜔_𝑟−1+𝐻_𝑗𝑚𝜔_𝑟−2,+𝑃_𝑟−1=6𝑟²𝜔^2𝑢𝐻^𝑗𝑚𝜔_𝑟−1−𝐻_𝑗𝑚𝜔_𝑟.(3.12) Then, substituting 𝑃_𝑖(𝑖=0,1,2,…,𝑟) into (3.5), we can obtain 𝐻_𝑗𝑚(𝜔)=⎧⎪⎪⎪⎪⎪⎪⎪⎨⎪⎪⎪⎪⎪⎪⎪⎩𝑟6𝜔_𝑢𝜔¹−𝜔³𝑃₀+𝜔−𝜔₀³𝑃₁+𝐻^𝑗𝑚𝜔₀−𝜔^2𝑢6𝑟²𝑃₀𝜔𝑟_𝑢𝜔₁+𝐻−𝜔_𝑗𝑚𝜔₁−𝜔^2𝑢6𝑟²𝑃₁𝜔𝑟_𝑢𝜔−𝜔₀𝜔₀𝑟≤𝜔≤𝜔₁,6𝜔_𝑢𝜔²−𝜔³𝑃₁+𝜔−𝜔₁³𝑃₂+𝐻^𝑗𝑚𝜔₁−𝜔^2𝑢6𝑟²𝑃₁𝜔𝑟_𝑢𝜔₂+𝐻−𝜔_𝑗𝑚𝜔₂−𝜔^2𝑢6𝑟²𝑃₂𝜔𝑟_𝑢𝜔−𝜔₁𝜔₁⋮𝑟≤𝜔≤𝜔₂,6𝜔_𝑢𝜔^𝑖−𝜔³𝑃_𝑖−1+𝜔−𝜔_𝑖−1³𝑃_𝑖+𝐻^𝑗𝑚𝜔_𝑖−1−𝜔^2𝑢6𝑟²𝑃_𝑖−1𝜔𝑟_𝑢𝜔_𝑖+𝐻−𝜔_𝑗𝑚𝜔_𝑖−𝜔^2𝑢6𝑟²𝑃_𝑖𝜔𝑟_𝑢𝜔−𝜔_𝑖−1𝜔_𝑟−1≤𝜔≤𝜔_𝑟.(3.13)

Introducing the spline interpolation into |𝐻_𝑗𝑚(𝜔)|, (2.1) for the simulation of the multivariate stationary stochastic processes can then be rewritten as follows:𝑓_𝑗(𝑡)=2√Δ𝜔^𝑛_𝑁_𝑚=1_𝑙=1𝐻_𝑗𝑚𝜔_𝑚𝑙cos𝜔_𝑚𝑙𝑡−𝜗_𝑗𝑚𝜔_𝑚𝑙+Φ_𝑚𝑙(𝑗=1,2,…,𝑛).(3.14)

3.2. Fast Fourier Transform (FFT) Algorithm

Taking into account that introducing the FFT algorithm renders very high computation efficiency in simulating the stationary processes (Li and Kareem [27]), with resorting to the FFT algorithm, (3.14) can then be derived as follows:𝑓_𝑗(𝑝Δ𝑡)=Re^𝑗_𝑚₌₁ℎ_𝑗𝑚(𝑝Δ𝑡)exp𝑖𝑚Δ𝜔𝑛(𝑝Δ𝑡),(3.15) where𝑝=0,1𝑀,…=,𝑛𝑇×₀𝑀−1;𝑇𝑛Δ𝑡₀;=2𝜋𝑛Δ𝜔=2𝜋𝜔𝑛_up𝑁ℎ_𝑗𝑚;(𝑝Δ𝑡)=⎧⎪⎨⎪⎩𝑔̃^𝑗𝑚𝑔̃_𝑗𝑚(𝑝Δ𝑡),𝑝=0,1,…,𝑀−1,⋮𝑔̃_𝑗𝑚[(𝑝−𝑀)Δ𝑡],𝑝⋮=𝑀,𝑀+1,…,2𝑀−1,𝑔̃_𝑗𝑚[(𝑝−(𝑛−1)𝑀)Δ𝑡],𝑝=(𝑛−1)𝑀,(𝑛−1)𝑀+1,…,𝑛𝑀−1,(𝑝Δ𝑡)=^𝑀−1_𝑙=1𝐵_𝑗𝑚𝑙exp𝑖𝑙𝑝2𝜋𝐵𝑀(𝑝=0,1,…,𝑀−1),_𝑗𝑚𝑙=√2Δ𝜔𝐻_𝑗𝑚𝑙Δ𝜔+𝑚Δ𝜔0,𝑛exp−𝑖𝜗^𝑗𝑚𝑙Δ𝜔+𝑚Δ𝑛𝜔exp𝑖Φ^𝑚𝑙,0𝑁≤<𝑙𝑙≤≤𝑁𝑀,−1.(3.16)

It is noteworthy that 𝑔̃_𝑗𝑚(𝑝Δ𝑡) can be obtained in terms of the inverse FFT of 𝐵_𝑗𝑚𝑙.

4. Numerical Results

In the following, two cases of performing the simulation of multivariate stationary fluctuating wind velocity are taken into consideration in order to demonstrate the capabilities and computational efficiency of the proposed SFFT algorithm.

4.1. Description of Case 1

A hundred velocity points to be simulated are evenly distributed from 10 m to 208 m heights along a vertical line (see Case 1 in Figure 2). The process corresponding to these hundred components is denoted by 𝑓^0𝑗(𝑡)[𝑓⁰¹(𝑡),𝑓⁰²(𝑡),…,𝑓⁰¹⁰⁰(𝑡)]. It is assumed that the mean value of the process is equal to zero. Then, the elements of its CPSD matrix may be given as 𝑆𝑆^0𝑗𝑗(𝜔)=𝑆_𝑗(𝜔)(𝑗=1,2,…,100),(4.1)_0𝑗𝑘(𝜔)=Γ_𝑗𝑘(𝜔)𝑆_𝑗(𝜔)𝑆_𝑘(𝜔)𝑒^{−𝑖𝜃𝑗𝑘(𝜔)}(𝑗,𝑘=1,2,…,100,𝑗≠𝑘),(4.2) in which 𝑆_𝑗(𝜔) represents the power spectral density function of 𝑓^0𝑗(𝑡);Γ_𝑗𝑘(𝜔) denotes the coherence function between 𝑓^0𝑗(𝑡) and 𝑓^0𝑘(𝑡); 𝜃_𝑗𝑘(𝜔) refers to the phase angle between 𝑆_𝑗(𝜔) and 𝑆_𝑘(𝜔).

Figure 2 

Configuration of points along a vertical line in Case 1 and Case 2.

It is emphasized here that the stochastic process described by (4.1) and (4.2) is nonhomogeneous in space, since the longitudinal velocity fluctuations 𝑓⁰¹(𝑡),𝑓⁰²(𝑡),…,𝑓⁰¹⁰⁰(𝑡) have different frequency contents, namely, 𝑆₁(𝜔)≠𝑆₂(𝜔)≠⋯≠𝑆₁₀₀(𝜔).

The following expression proposed by Kaimal et al. [40] is selected to model the two-sided power spectral density (PSD) function of the longitudinal wind velocity fluctuations at different height:𝑆(𝑧,𝜔)=12220𝜋0𝑢^2∗𝑧𝑈(𝑧)[1+50(𝜔𝑧1/2𝜋𝑈(𝑧))]^5/3,(4.3) where 𝑧 is the elevation above the ground, in meters; 𝜔 is the circular frequency, in rad/s; 𝑢_∗ is the shear velocity of the flow, in m/s; 𝑈(𝑧) is the mean wind speed at height 𝑧, in m/s.

Taking into account the following coherence function between the velocity fluctuations at two different heights 𝑧_𝑗 and 𝑧_𝑘 suggested by Davenport [41],Γ(Δ𝑧,𝜔)=exp−𝜔2𝜋𝐶^𝑧Δ𝑧(1/2)𝑈𝑧_𝑗+𝑈𝑧_𝑘(4.4) in which 𝑈(𝑧_𝑗) and 𝑈(𝑧_𝑘) are the mean wind speeds at heights 𝑧_𝑗 and 𝑧_𝑘, respectively; Δ𝑧=|𝑧_𝑗−𝑧_𝑘|; 𝐶_𝑧 is a constant that may be set equal to 10 for structural design purposes (Kristensen and Jensen [42] and Simiu and Scanlan [43]).

Furthermore, the expression of the phase angle 𝜃_𝑗𝑘(𝜔) is taken from Di Paola [44] and Simiu and Scanlan [43],𝜃_𝑗𝑘(𝜔)=𝜔𝑧^𝑗−𝑧_𝑘𝜐^(𝑗,𝑘)_app.(4.5) In (4.2), 𝑒^{−𝑖𝜃𝑗𝑘(𝜔)} is a measure of the wave passage delay due to the apparent velocity of waves 𝜐^(𝑗,𝑘)_app.

Likewise, the apparent velocity of waves can be assumed to be the following form suggested by Simiu and Scanlan [43]:𝜐^(𝑗,𝑘)_app=𝜋𝑈𝑧^𝑗+𝑈𝑧_𝑘𝐶_𝜃,(4.6) in which 𝐶_𝜃 is an appropriate coefficient that has to be determined from experimental data.

In the present paper, the expression by Peil and Telljohann [45] is taken into consideration, which has the form 𝜐^(𝑗,𝑘)_app=𝑈[(𝑧_𝑗+𝑧_𝑘)/2]/5.5. The presence of the phase angle turns into a time shift of the peak of the cross-correlation functions from the data of the measured fluctuations at different heights.

Let it be supposed that the mean wind velocity at the first point 1 (𝑧=10m) is 𝑈(10)=26m/s and that the surface rough roughness length is 𝑧₀=0.001266 m (corresponding shear velocity of the flow 𝑢_∗=1.76 m/s). It is worth mentioning that the values 𝑧₀=0.001266mand 𝑢_∗=1.76m/s are taken from Simiu and Scanlan [43]. The upper cutoff frequency is 𝜔_up=2𝜋 rad/s. The dividing number of frequency is 𝑁=8192. The simulation is carried out at 1×𝑀=1×8192×2=16384 time instants with a time step Δ𝑡=0.5s, over a time length equal to 𝑇=8192s. The mean wind speeds at other points are computed in terms of the logarithmic law (i.e., the vertical profile), which has the form as follows:𝑈𝑧_𝑗=ln𝑧^𝑗/𝑧₀ln𝑧_𝑠/𝑧₀𝑈𝑧^𝑠.(4.7)

4.2. Description of Case 2

Fifty velocity points to be simulated are evenly distributed from 10 m to 255 m heights along a vertical line (see Case 2 in Figure 2). The process corresponding to these fifty components is represented by 𝑓^0𝑗(𝑡)[𝑓⁰¹(𝑡),𝑓⁰²(𝑡),…,𝑓⁰⁵⁰(𝑡)]. The dividing number of frequency is 𝑁=2048. The numerical simulation is implemented at 1×𝑀=1×2048×2=4096 time instants with a time step Δ𝑡=0.5s, over a time length equal to 𝑇=2048s. The other conditions are the same as Case  1.

4.3. Numerical Analyses

The comparisons among several different schemes are made and listed in Table 1. The Cholesky decomposition with the SFFT algorithm is implemented on 128 circular frequency points. It can be seen that in Case 1, the elapsed time for the simulation with the FFT algorithm is 215.9 min, while that for the simulation with the SFFT algorithm is only 10.2 min; in the Case 2, the elapsed time for the simulation with the FFT algorithm is 5.9 min, while that for the simulation with the SFFT algorithm is only 0.5 min. Apparently, in the calculation speed the SFFT algorithm has considerable advantage over the FFT algorithm. Likewise, the superiority becomes more remarkable with the increasing of 𝑁,𝑛, and 𝑇. This is because the number of the Cholesky decomposition with the FFT algorithm equals 𝑛×𝑁, while that with the SFFT algorithm keeps constant with the value equal to 128. For a determined SFFT or FFT algorithm, however, the computational efficiency becomes lower with the increasing of the parameters 𝑁,𝑛, and 𝑇.

Table 1 

Comparisons among several different schemes for wind velocity simulation.

Since in Case 1 the number of the stimulated time histories is very large, without any loss of generality, arbitrary six time histories of all are chosen to demonstrate the effort of simulating the longitudinal wind velocity fluctuations by resorting to the proposed SFFT algorithm. The generated samples of longitudinal wind velocity fluctuations at six points 1, 11, 26, 46, 71, and 100 (see Case  1 in Figure 2), denoted by 𝑓₁(𝑡),𝑓₁₁(𝑡),𝑓₂₆(𝑡),𝑓₄₆(𝑡),𝑓₇₁(𝑡), and 𝑓₁₀₀(𝑡), respectively, are displayed in Figures 3(a), 3(b), 3(c), 3(d), 3(e), and 3(f) for a time length 𝑇=8192s. Likewise, the parts in the first 600 s of these six time histories are further extracted out, respectively, as shown Figure 4, in order to better visualize the differences and the similarities among these six time histories.

(a) Height 𝑧=10m

(b) Height 𝑧=30m

(c) Height 𝑧=60m

(d) Height 𝑧=100m

(e) Height 𝑧=150m

(f) Height 𝑧=208m

Figure 3 

Generated sample functions for longitudinal wind velocity fluctuations at six different heights in Case 1.

(a) Height 𝑧=10m

(b) Height 𝑧=30m

(c) Height 𝑧=60m

(d) Height 𝑧=100m

(e) Height 𝑧=150m

(f) Height 𝑧=208m

Figure 4 

Generated sample functions for longitudinal wind velocity fluctuations at six different heights, over first 600 s of generated sample functions displayed in Figure 3.

As far as the correlation degree among these six time histories is concerned, it is clearly seen from Figure 4 that with the decreasing of the distance between arbitrary two points, the loss of coherence between two samples corresponding, respectively, to these two points will become smaller. For example, since the distance between the points 1 and 11 is the smallest in these six points, consequently, there is the smallest loss of coherence between 𝑓₁(𝑡) and 𝑓₁₁(𝑡). On the other hand, there is the biggest loss of coherence between 𝑓₁(𝑡) and 𝑓₁₀₀(𝑡), as the point 100 is located at 198 m height from the point 1 (see Case 1 in Figure 2). This behavior is controlled by the coherence functions (see (4.4)). Furthermore, the phase angle can also be detected clearly at 100 s time instant between arbitrary two samples of these six time histories. Likewise, the phase angle of these six time histories will become larger with the increasing of the distance between arbitrary two time histories. Similarly, this phenomenon can be controlled by the phase angle functions (see (4.5)).

The forgoing elucidations indicate that the proposed SFFT algorithm is able to simulate the longitudinal wind velocity fluctuations that are spatially correlated according to a prescribed coherence function and that possess the phase angles following a prescribed phase angle function.

Plotted in Figure 5 is the temporal autocorrelation function [𝑅_𝑗𝑗(𝜏)(𝑗=1,11,26,46,71,100)] of the generated sample function shown in Figure 3, along with the target autocorrelation function [𝑅^0𝑗𝑗(𝜏)(𝑗=1,11,26,46,71,100)]. It is worth pointing out that the target autocorrelation function [𝑅^0𝑗𝑗(𝜏)] is computed with resorting to the Wiener-Khintchine transformation. As can be seen in Figure 5, the simulated autocorrelation function [𝑅_𝑗𝑗(𝜏)] practically coincides with the target autocorrelation function [𝑅^0𝑗𝑗(𝜏)]. Displayed in Figure 6 is the temporal cross-correlation function [𝑅_𝑗𝑘(𝜏)(𝑗=1,𝑘=11,26,46,71,100)] of the generated sample function presented in Figure 3, along with the target cross-correlation function [𝑅^0𝑗𝑘(𝜏)(𝑗=1,𝑘=11,26,46,71,100)]. The same, the target cross-correlation function [𝑅^0𝑗𝑘(𝜏)] is calculated according to the Wiener-Khintchine transformation. As expected, the simulated cross-correlation function [𝑅_𝑗𝑘(𝜏)(𝑗≠𝑘)] practically coincides with the target cross-correlation function [𝑅^0𝑗𝑘(𝜏)(𝑗≠𝑘)]. Apparently, the temporal auto correlation and cross-correlation functions of any sample function [𝑓_𝑗(𝑡)(𝑗=1,2,…,100)] are, respectively, identical to the target autocorrelation and cross-correlation functions based on the SFFT algorithm. Moreover, Figure 6 demonstrates that the presence of the phase angle turns into a time shift of the peak of the cross-correlation functions. Likewise, the time shift will become larger with the increasing of the distance between arbitrary two points. It is proved once again that the proposed SFFT algorithm can efficiently carry out the simulation of longitudinal wind velocity fluctuations with the phase angles.

(a) Height 𝑧=10m

(b) Height 𝑧=30m

(c) Height 𝑧=60m

(d) Height 𝑧=100m

(e) Height 𝑧=150m

(f) Height 𝑧=208m

Figure 5 

Temporal autocorrelation functions of generated sample function displayed in Figure 3 versus corresponding targets.

(a) 𝑧=10m and 𝑧=30m

(b) 𝑧=10m and 𝑧=60m

(c) 𝑧=10m and 𝑧=100m

(d) 𝑧=10m and 𝑧=150m

(e) 𝑧=10m and 𝑧=208m

Figure 6 

Temporal cross-correlation functions of generated sample functions displayed in Figure 3 versus corresponding targets.

From Figure 7, it is clear that the estimated spectra derived from these sample time histories and the target spectra are in remarkably good agreement with each other. Likewise, it is seen from Figures 8 and 9 that the mean and variance values of these 100 time histories generated with the SFFT algorithm approach those with the FFT algorithm, which implies that the SFFT and FFT algorithms practically achieve the same computation accuracy. More specifically, Figure 8 shows that the variance values of the generated time histories with resorting to both the SFFT and FFT are identical to the target values, respectively. It is observed in Figure 9 that the mean values of the generated time histories with resorting to both the SFFT and FFT are not more than 0.2 m/s, effectively meaning that these mean values are considerably small with respect to the zero target value.

(a) Height 𝑧=10m

(b) Height 𝑧=30m

(c) Height 𝑧=60m

(d) Height 𝑧=100m

(e) Height 𝑧=150m

(f) Height 𝑧=208m

Figure 7 

Power spectrum of generated sample functions displayed in Figure 3 versus corresponding targets.

Figure 8 

Variance values of sample generated by resorting to the SFFT or the FFT with respect to corresponding targets in Case 1.

Figure 9 

Mean values of sample generated with resorting to the SFFT or the FFT with reference to corresponding targets in Case 1.

In order to demonstrate Gaussian of the simulated wind velocity fluctuations, Figure 10 shows the probability density functions of generated samples displayed in Figure 3 versus corresponding target (Gaussian). It is clear that the probability density functions estimated from these samples match closely Gaussian probability density function.

(a) 𝑧=10m

(b) 𝑧=30m

(c) 𝑧=60m

(d) 𝑧=100m

(e) 𝑧=150m

(f) 𝑧=208m

Figure 10 

Probability density functions of generated sample functions displayed in Figure 3 versus corresponding target (Gaussian).

4.4. Accuracy and Time Expense of Various Interpolations

Four kinds of interpolation techniques, such as SFFT, the cubic Lagrangian interpolation based FFT (CL-FFT), square Lagrangian interpolation based FFT (SL-FFT), and neural network interpolation based FFT (NN-FFT), are taken into consideration to simulate the wind velocity field for the evaluation of accuracy and time expense of various interpolations [46]. The root mean square error (RMSE) and error factor (EF) [47] are introduced to evaluate the accuracy of these four interpolation-based FFT techniques:RMSE=⎷1𝑁^𝑁_𝑖=1𝑦_𝑖EF=1−∑_𝑁𝑖=1−𝑌_𝑖²,||𝑦_𝑖−𝑌_𝑖||∑^𝑁𝑖=1||𝑌_𝑖||,(4.8) in which, 𝑁 is the total number of time points in the generated sample; 𝑌_𝑖 and 𝑦_𝑖 signify generated wind velocities at a certain instant employing the FFT and the interpolation-based FFT, respectively. It is noted that, the interpolation based FFT is more precise if the value of RMSE is closer to zero; the interpolation-based FFT is more precise if the value of EF is closer to one.

The data coming from [46] are listed in Tables 2 and 3. Table 2 shows the elapsed time in the simulation using the SL-FFT, CL-FFT, SFFT, NN-FFT, and FFT, respectively, and Table 3 records both the RMSE and EF values. It is from Table 2 seen that the efficiency of the SFFT is the highest due to the minimal time expense. It is from Table 3 observed that the NN-FFT presents the highest accuracy. Worth noting, the accuracy by the SFFT is a little lower than that by the NN-FFT. Forasmuch, the SFFT, generally, is the best.

Table 2 

Comparison of time expense among different simulation schemes (noted that the data come from [46]).

Table 3 

Comparisons of MRSE and EF among different simulation schemes (noted that the data come from [46]).

5. Conclusions

This study presented the SFFT algorithm to further enhance the computational efficiency of simulating the multivariate stochastic processes. More specifically, the main purpose of the proposed SFFT algorithm, simultaneously with resorting to the spline interpolation and FFT algorithms, is to reduce the number of the Cholesky decomposition. The main conclusions of the present research effort on the SFFT algorithm can be drawn as follows.(1)The proposed SFFT algorithm is able to efficiently generate the longitudinal wind velocity fluctuations that are spatially correlated according to a prescribed coherence function and that possess the phase angles following a prescribed phase angle function.(2)The proposed SFFT algorithm has considerable advantage over the FFT algorithm. Likewise, this superiority will become more and more remarkable with the increasing of the parameters 𝑁,𝑛, and 𝑇.(3)Although introducing the spline interpolation approximation in decomposing the cross-power spectral density matrix, the SFFT and FFT algorithms practically achieve the same computational accuracy.(4)Employing the SFFT algorithm, the temporal autocorrelation and cross-correlation functions of any sample function are identical to the target autocorrelation and cross-correlation functions, respectively. Likewise, the estimated spectra derived from these sample time histories and the target spectra are in remarkably good agreement with each other.(5)The SFFT, generally, is the best in many interpolation methods.

Eventually, it is worth emphasizing that since the spline interpolation is able to significantly reduce the number of the Cholesky decomposition of time-varying power spectral matrix, the SFFT algorithm can improve the performance of the SR method for the simulation of the multivariate nonstationary stochastic processes. It is known that the power spectrum of nonstationary stochastic process is time-varying; thus, meaning that the FFT algorithm is very hard to be used in this case. Only when the time-varying power spectrum satisfies some special conditions, the FFT can be utilized. However, it is anticipated that the SFFT algorithm can be used to enhance the computational efficiency of simulating the multivariate nonstationary stochastic processes, since the spline interpolation technique is not affected by the time-varying power spectrum.

Appendix

Cholesky Decomposition

For a given CPSD matrix with the element |𝑆^0𝑗𝑘(𝜔)|𝑒^{𝑖𝜃𝑗𝑘}(𝑘=1,2,…,𝑛;𝑗=1,2,…,𝑛), the element of the lower triangular matrix 𝐻(𝜔) can be computed as follows:𝐻₁₁(𝜔)=||𝐻₁₁(𝜔)||𝑒^𝑖𝜗11=||𝑆⁰¹¹(𝜔)||𝑒^𝑖𝜃11^1/2=||𝑆⁰¹¹(𝜔)||^1/2𝑒^{𝑖(1/2)𝜃11}⟹𝜗||𝐻¹¹₁₁(𝜔)||=||𝑆⁰¹¹(𝜔)||^1/2,𝐻=𝜃₁₁=0,_𝑗1(𝜔)=||𝐻_𝑗1(𝜔)||𝑒^{𝑖𝜗𝑗1}=||𝑆^0𝑗1(𝜔)||𝑒^{𝑖𝜃𝑗1}||𝐻₁₁(𝜔)||𝑒^𝑖𝜗11=||𝑆^0𝑗1(𝜔)||||𝐻₁₁(𝜔)||𝑒^{𝑖𝜃𝑗1}⟹⎧⎪⎨⎪⎩||𝐻(𝑗=2,3,…,𝑛)_𝑗1(𝜔)||=||𝑆^0𝑗𝑘(𝜔)||𝜗||𝐻₁₁(𝜔)||(𝑗=2,3,…,𝑛),_𝑗1𝐻=𝜃_𝑗1(𝑗=2,3,…,𝑛),_𝑘𝑘(𝜔)=||𝐻_𝑘𝑘(𝜔)||𝑒^{𝑖𝜗𝑘𝑘}=||𝑆^0𝑘𝑘(𝜔)||𝑒^{𝑖𝜃𝑘𝑘}−^𝑘−1_𝑟=1||𝐻_𝑘𝑟(𝜔)||𝑒^{𝑖𝜗𝑘𝑟}||𝐻_𝑟𝑘(𝜔)||𝑒^{−𝑖𝜗𝑘𝑟}^1/2=||𝑆^0𝑘𝑘(𝜔)||−^𝑘−1_𝑟=1⟹⎧⎪⎨⎪⎩||𝐻||𝐻_𝑘𝑟(𝜔)||⋅||𝐻_𝑟𝑘(𝜔)||^1/2𝑒^{𝑖(1/2)𝜃𝑘𝑘}(𝑘=2,3,…,𝑛)_𝑘𝑘(𝜔)||=||𝑆^0𝑘𝑘(𝜔)||−^𝑘−1_𝑟=1𝜗||𝐻_𝑘𝑟(𝜔)||⋅||𝐻_𝑟𝑘(𝜔)||^1/2(𝑘=2,3,…,𝑛),_𝑘𝑘𝐻=𝜃_𝑘𝑘=0(𝑘=2,3,…,𝑛),_𝑗𝑘(𝜔)=||𝐻_𝑗𝑘(𝜔)||𝑒^{𝑖𝜗𝑗𝑘}=||𝑆^0𝑗𝑘(𝜔)||𝑒^{𝑖𝜃𝑗𝑘}−∑^{𝑘−1𝑟=1}||𝐻_𝑗𝑟(𝜔)||𝑒^{𝑖𝜗𝑗𝑟}||𝐻_𝑟𝑘(𝜔)||𝑒^{−𝑖𝜗𝑘𝑟}||𝐻_𝑘𝑘(𝜔)||𝑒^{𝑖𝜃𝑘𝑘}=||𝑆^0𝑗𝑘(𝜔)||𝑒^{𝑖𝜃𝑗𝑘}−∑^{𝑘−1𝑟=1}||𝐻_𝑗𝑟(𝜔)||⋅||𝐻_𝑟𝑘(𝜔)||𝑒^{𝑖(𝜗𝑗𝑟−𝜗𝑘𝑟)}||𝐻_𝑘𝑘(𝜔)||𝑒^{𝑖𝜃𝑘𝑘}=||𝑆^0𝑗𝑘(𝜔)||𝑒^{𝑖𝜃𝑗𝑘}−∑^{𝑘−1𝑟=1}||𝐻_𝑗𝑟(𝜔)||⋅||𝐻_𝑟𝑘(𝜔)||𝑒^{𝑖𝜗𝑗𝑘}||𝐻_𝑘𝑘(𝜔)||𝑒^{𝑖𝜃𝑘𝑘}=||𝑆^0𝑗𝑘(𝜔)||−∑^{𝑘−1𝑟=1}||𝐻_𝑗𝑟(𝜔)||⋅||𝐻_𝑟𝑘(𝜔)||||𝐻_𝑘𝑘(𝜔)||𝑒^{𝑖𝜃𝑗𝑘}⟹⎧⎪⎪⎨⎪⎪⎩||𝐻(𝑘=2,3,…,𝑛;𝑗=𝑘+1,𝑘+2,…,𝑛)_𝑗𝑘(𝜔)||=||𝑆^0𝑗𝑘(𝜔)||−∑^{𝑘−1𝑟=1}||𝐻_𝑗𝑟(𝜔)||⋅||𝐻_𝑟𝑘(𝜔)||𝜗||𝐻_𝑘𝑘((𝑘𝜔=)||1,2,…,𝑛;𝑗=𝑘+1,𝑘+2,…,𝑛),_𝑗𝑘=𝜃_𝑗𝑘(𝑘=1,2,…,𝑛;𝑗=𝑘+1,𝑘+2,…,𝑛).(A.1)

Subsequently, |𝐻_𝑗𝑘(𝜔)|(𝑘=1,2,…,𝑛;𝑗=𝑘,𝑘+1,𝑘+2,…,𝑛) can be obtained by resorting to the Cholesky decomposition of a real matrix which consists of the element |𝑆^0𝑗𝑘(𝜔)|(𝑘=1,2,…,𝑛;𝑗=1,2,…,𝑛). Likewise, the phase angle contents with the relationship 𝜗_𝑗𝑘=𝜃_𝑗𝑘(𝑘=1,2,…,𝑛;𝑗=1,2,…,𝑛).

Acknowledgments

The work reported in this paper has been supported by the National Science Foundation of China with Grant no. 11162005. This support is gratefully acknowledged. The authors wish to thank the reviewers for their careful, unbiased, and constructive suggestions that significantly improved the quality of this paper.