Computer Simulation of Stochastic Wind Velocity Fields for Structural Response Analysis: Comparisons and Applications

Abstract

The digital simulation of wind velocity fields, modeled as multivariate stationary Gaussian processes, is a widely adopted tool to generate the external input for response analysis of wind-sensitive nonlinear structures. The problem does not entail any theoretical difficulty, existing already a large number of well-established techniques, such as the accurate weighted amplitude wave superposition (WAWS) method. However, reducing the computational effort required by the WAWS method is sometimes necessary, especially when dealing with complex structures and high-dimensional simulation domains. In these cases, approximate formulas must be adopted, which however require an appropriate tuning of some fundamental parameters in such a way to achieve an acceptable level of accuracy if compared to that obtained using the WAWS method. Among the different techniques available for this purpose, autoregressive (AR) filters and algorithms exploiting the proper orthogonal decomposition (POD) of the spectral matrix deserve a special attention. In this paper, a properly organized way for implementing stochastic wind simulation algorithms is outlined at first. Then, taking the WAWS method as a reference from the viewpoint of the accuracy of the simulated samples, a comparative study between POD-based and AR techniques is proposed, with a particular attention to computational effort and memory requirements.

1. Introduction

The simulation of wind velocity fields has been one of the main topics of wind engineering for the last decades. In this framework, wind velocity is usually idealized as the sum of a mean part, assumed as constant within a conventional time interval, and a fluctuating part representing the atmospheric turbulence. This last is usually modeled as a stationary zero-mean Gaussian random process [1, 2].

Several techniques were proposed in the literature in order to simulate Gaussian wind velocity fields to be employed in structural analysis [3–5]. Among those, the classic WAWS method, based on the pioneering work by Shinozuka and Jan [6] and then modified by Deodatis [1] in such a way to achieve ergodic realizations and to be efficiently implemented through fast Fourier transform (FFT) algorithms, has proved to guarantee the best quality of the obtained results [7]. Nevertheless, such a procedure requires the Cholesky factorization of the spectral matrix, which unfortunately leads to high computational expenses, especially when dealing with complex structures and high-dimensional simulation domains. These difficulties are mainly related to memory allocation and time consuming operations, thus requiring, on the one hand, the reduction of the problem size. On the other hand, an accurate wind simulation is essential for predicting the wind-induced response of flexible structures, such as transmission power lines, tall buildings, suspension, and cable-stayed bridges. Less demanding, yet approximate, procedures may be obtained by exploiting the properties of the POD decomposition of the spectral matrix, proposed in the papers by Li and Kareem [8], Di Paola [9], and Solari and Carassale [10]. A POD-based technique, in particular, was recently applied to simulate the wind velocity field on a domain representing a long-span suspension bridge [11] with significantly low computational efforts. A third well-established class of simulation formulas is represented by autoregressive (AR) and autoregressive moving average (ARMA) filters, early introduced in the paper by Samaras et al. [12] and applied to wind simulation by Mignolet and Spanos [13], by Li and Kareem [8] and, more recently, by Di Paola and Gullo [14].

This paper addresses, at first, the problem of the efficient implementation of wind simulation algorithms by proposing a two-steps approach which allows to properly handle large numbers of simulations, which is usually the case for performing structural response analysis. Then, the WAWS method is adopted as a reference for calibrating the parameters of POD-based and AR techniques. These last methods are finally compared in terms of computational effort and memory requirements. To this end, two numerical examples are presented. The former, represented by the tower of a suspension bridge, is considered for models calibration. The latter, represented by an entire suspension bridge, is a more demanding case, which is here chosen as a benchmark in order to test the capability of AR and POD-based techniques in handling complex simulations and to compare their computational efficiency. Finally, the convergence rates of the accuracies of the two simplified methods are also investigated, in order to provide indications for tuning the parameters that affect their results.

The considered formulas for wind simulation are only some of those currently available in the literature and this work does not attempt in any way to give a comprehensive literature review on this large topic. In particular, it must be underlined that non-Gaussian techniques were also studied and applied in the literature with the aim of directly simulating wind pressure fields on structures. As examples, this problem was analyzed by Grigoriu [15, 16], by Popescu et al. [17], by Gioffrè et al. [18] and by Borri and Facchini [19], among others. Non-Gaussian responses of wind-excited structures were also studied (see, e.g., Gusella and Materazzi [20, 21]). A quite exhaustive review of all the approaches currently available for simulating wind velocity fields and structural responses was recently proposed by Kareem [22] also considering the nonstationary and non-Gaussian cases.

2. Stochastic Wind Modeling

Let 𝑂𝑥₁𝑥₂𝑥₃ be the global Cartesian reference system with the origin 𝑂 lying on the ground. In the following developments, without loss of generality, it is assumed that the 𝑥₃ axis is parallel to the gravity direction and the mean wind velocity 𝐔 is parallel to the x₁ axis. Following these definitions, the wind velocity field 𝐕(𝐱,𝑡) is idealized as the sum of a mean value 𝐔(𝑥₃), function of the elevation from the ground 𝑥₃, and a stationary zero-mean fluctuation 𝐮(𝐱,𝑡)=(𝑢₁,𝑢₂,𝑢₃), that depends on the position 𝐱=(𝑥₁,𝑥₂,𝑥₃) and varies in time 𝑡. The three components 𝑢₁, 𝑢₂, 𝑢₃ of the vector 𝐮(𝐱,𝑡) represent therefore the longitudinal, lateral and vertical components of turbulence, respectively. As customary in wind engineering, the classic logarithmic law [4] is assumed to represent the variability of the mean wind velocity 𝐔 with 𝑥₃.

The atmospheric turbulence 𝐮 is usually modeled as a zero-mean, Gaussian, stationary random field that depends on time. If 𝑃 and 𝑃^′ denote two points located at 𝐱 and 𝐱^′, then, from a probabilistic point of view, the complete characterization of this field is ensured by the knowledge of the correlation function 𝑅_{𝑢𝑖𝑢𝑗}(𝐱,𝐱^′,𝜏) for every pair of turbulence components 𝑢_𝑖,𝑢_𝑗,𝜏 being the time lag. Assuming that 𝐮 is ergodic, 𝑅_{𝑢𝑖𝑢𝑗}(𝐱,𝐱^′,𝜏) is given by the Fourier transform of the cross power spectral density (CPSD) function 𝑆_{𝑢𝑖𝑢𝑗}(𝐱,𝐱^′,𝜔) between 𝑢_𝑖 and 𝑢_𝑗, 𝜔 being the circular frequency. Most of the theoretical models adopted in wind engineering express the CPSD in terms of auto-spectra 𝑆_{𝑢𝑖𝑢𝑖}(𝐱,𝑛), 𝑆_{𝑢𝑗𝑢𝑗}(𝐱^′,𝑛) and coherence function Coh_{𝑣𝑖𝑣𝑗}(𝐱,𝐱^′,𝑛), for 𝑖,𝑗=1,2,3, where 𝑛=𝜔/2𝜋 denotes the frequency [9]. Several models were developed in the literature to give an analytical representation to 𝑆_{𝑢𝑗𝑢𝑗}(𝐱,𝑛), for 𝑗=1,2,3. Here, the model by Solari and Piccardo [23] is adopted for this purpose and Coh_{𝑣𝑖𝑣𝑗}(𝐱,𝐱^′,𝑛) is expressed by means of the classic exponential law [11]. The imaginary part of the CPSD (called quadrature spectrum) introduces a time lag between the simulated velocities, which may become significant for points placed in the along-wind direction 𝑥₁. Without loss of generality, the quadrature spectrum is here neglected, since the considered examples focus on simulation domains located on the 𝑥₂𝑥₃ plane. This, however, does not limit the generality of the presented results, which apply also to the case in which the CPSD has a nonnegative imaginary part.

For simulation purposes, the spatial domain is discretized into 𝑁 points which usually represent significant nodes of the case study structure. The position of the 𝑘th point is identified by the vector 𝐱^(𝑘)=(𝑥₁^(𝑘),𝑥₂^(𝑘),𝑥₃^(𝑘)), with 𝑘=1,2,…,𝑁. This spatial discretization allows to represent the wind field 𝐕^(𝑘)(𝑡)=𝐕(𝐱^(𝑘),𝑡) as the sum of the mean velocity vectors 𝐔^(𝑘)=𝐔(𝐱^(𝑘)) and the turbulence vectors 𝐮^(𝑘)(𝑡)=𝐮(𝐱^(𝑘),𝑡)=(𝑢₁^(𝑘),𝑢₂^(𝑘),𝑢₃^(𝑘)). Following this approach, the time-dependent random field 𝐮(𝐱,𝑡) is transformed into a 3N-variate stationary random process 𝐮(𝑡), where 𝐮(𝑡) is a 3N-order vector containing the components 𝑢₁^(𝑘), 𝑢₂^(𝑘), 𝑢₃^(𝑘) of the turbulence vectors 𝐮^(𝑘)(𝑡) for 𝑘=1,2,…,𝑁. The complete characterization of the process 𝐮(𝑡) is given by its power spectral density (PSD) matrix 𝐒_𝐮𝐮(𝑛) containing PSD and CPSD functions between turbulent velocities. Equivalently, the two-side PSD matrix 𝐆_𝐮𝐮(𝜔) can also be utilized [11]. Having neglected the quadrature spectrum, both matrices 𝐒_𝐮𝐮(𝑛) and 𝐆_𝐮𝐮(𝜔) are real and positive definite.

3. Digital Simulation Formulas

Realizations of the process 𝐮(𝑡) are generated along a sequence 𝑡_𝑗 for 𝑗=1,2,…,𝑁_𝑡, of time instants with constant time step Δ𝑡. The circular frequency domain is limited within the interval [−𝜔_𝑐,𝜔_𝑐] where 𝜔_𝑐 denotes the cut-off circular frequency. Such an interval is discretized by means of an equally spaced sequence 𝜔_𝑘 for 𝑘=1,2,…,𝑁_𝜔 of circular frequencies, with step amplitude Δ𝜔, being thereby 𝜔_𝑐=𝑁_𝜔/2Δ𝜔.

The most accurate and well-established algorithm for simulating samples of Gaussian processes is the spectral representation method, early proposed by Shinozuka and Jan [6]. This last was significantly improved by Deodatis [1] who proposed a simulation formula that ensures the generation of ergodic sample functions and which is suitable to be implemented using efficient FFT algorithms. The Shinozuka-Deodatis formula is based on the following frequency-dependent decomposition of 𝐆_𝐮𝐮(𝜔_𝑘):𝐆_𝐮𝐮(𝜔_𝑘)=𝐓(𝜔_𝑘)𝐓∗(𝜔_𝑘)^𝑇,𝑘=1,2,…,𝑁_𝜔, (1) where * denotes the complex conjugate operator. It is worth noting that the decomposition (1) is not unique. To this regard, in the Shinozuka-Deodatis method the Cholesky factorization algorithm is adopted which provides 𝐓(𝜔_𝑘) in lower triangular form. Based on this approach, the following digital simulation formula was obtained [1]:𝑢_𝑖(𝑡_𝑗)≅√Δ𝜔^𝑖𝑁_𝑟=1^𝜔_𝑘=1𝑇_𝑖𝑟(𝜔_𝑟𝑘)cos(𝜔_𝑟𝑘𝑡_𝑗+𝜑_𝑟𝑘),𝑖=1,2,…,3𝑁,𝑗=1,2,…,𝑁_𝑡, (2) where 𝜔_𝑟𝑘=𝜔_𝑘+𝑟Δ𝜔/3𝑁, 𝑇_𝑖𝑟(𝜔_𝑟𝑘) are the elements of matrix 𝐓(𝜔_𝑟𝑘) and 𝜑_{𝑟𝑘𝑟𝑘} are independentrandom phases uniformly distributed in the interval [0,2𝜋]. Without loss of generality, (2) is written under the assumption of neglecting the quadrature spectrum, therefore 𝐓(𝜔_𝑟𝑘) being a real matrix. This method is sometimes called “weighted amplitude waves superposition” (WAWS) method. As observed by Di Paola [9], the central limit theorem ensures that the process simulated by means of (2) is asymptotically Gaussian as 𝑁_𝜔 becomes large. Equation (2) can be efficiently implemented through a fast Fourier transform (FFT) algorithm, as discussed in [1].

An alternative approach to the WAWS formula is represented by the POD-based technique proposed by Carassale and Solari [11]. With such an approach the matrices 𝐓(𝜔_𝑘), in (1), are calculated as𝐓(𝜔_𝑘)=^3𝑁_𝑟=1𝜽_𝑟^(𝑘)𝛾_𝑟^(𝑘),𝑘=1,2,…,𝑁_𝜔, (3) where 𝛾_𝑟^(𝑘), for 𝑟=1,2,…,3𝑁, are the eigenvalues of 𝐆_𝐮𝐮(𝜔_𝑘) and 𝜽_𝑟^(𝑘) are the corresponding eigenvectors. These last assume the well-known physical meaning of wind blowing mode shapes [9]. The use of (3) leads to the following simulation formula:𝐮(𝑡_𝑗)≅^3𝑁_𝑟=1^𝑁𝜔_𝑘=1𝑒^{𝜄𝜔𝑘𝑡𝑗}Δ𝜔𝛾_𝑟^(𝑘)𝜽_𝑟^(𝑘)𝑝_𝑟^(𝑘),𝑗=1,2,…,𝑁_𝑡, (4) where 𝜄 is the imaginary unit and 𝑝_𝑟^(𝑘) are complex-valued Gaussian random numbers with unit variance. Likewise (2), (4) can be efficiently implemented through a fast Fourier transform (FFT) algorithm as discussed in [11], assuming again 𝑁_𝜔=𝑁_𝑡 and Δ𝜔=2𝜋/𝑁_𝑡Δ𝑡.

Equation (4) generates 𝐮(𝑡) as the superposition of 3N independent fully coherent stochastic processes, which represent the contributions of the different wind modes. A very convenient way for reducing both computer time and allocated memory using (4) is to calculate eigenvectors and eigenvalues of 𝐆_𝐮𝐮(𝜔_𝑘) only for a small number 𝑁_𝜔≪𝑁_𝜔 of circular frequencies distributed along a reduced sequence 𝜔_𝑘, for 𝑘=1,2,…,𝑁_𝜔, and then use interpolation formulas elsewhere. This approach, proposed by Carassale and Solari [11], is here adopted and briefly recalled in Appendix A.

Autoregressive methods generate the wind velocity field by filtering a 3N-order vector 𝐚(𝑡) of uncorrelated band-limited Gaussian white noises 𝑎_𝑗(𝑡), for 𝑗=1,2,…,3𝑁, with unit variance. The turbulent wind process is thus expressed in the following form:𝐮(𝑡_𝑗)≅^𝑝_𝑖=1Φ_𝑖𝐮(𝑡_𝑗+−𝑖Δ𝑡)_𝑞𝑖=0𝐁_𝑖𝐚(𝑡_𝑗−𝑖Δ𝑡),𝑗=1,2,…,𝑁_𝑡, (5) where Φ_𝑖 and 𝐁_𝑖 are convenient 3𝑁×3𝑁 matrices that can be easily calculated by imposing that the simulated process satisfies the covariance structure of the target one [12]. Equation (5) represents a general ARMA method, in which 𝑝 denotes the order of autoregression and 𝑞 is the order of the moving average component. As it is well known [7] an ARMA(𝑝,𝑞) can be approximated by an AR(𝑝₁), with 𝑝₁≫𝑝, where the AR(𝑝₁) filter is readily obtained by assuming 𝑞=0 in (5). Regarding this point, it must be mentioned that AR(20) and ARMA(5,5) provided good results for turbulent wind simulations in many different environmental conditions [7]. The necessary derivations for calculating the coefficient matrices of an AR(p₁) filter are reported in Appendix A.

4. Efficient Implementation of Wind Simulation Methods

The three simulation algorithms presented in Section 3 are implemented in the MATLAB R2009a [24] environment to generate samples of along-wind and across-wind turbulence velocities 𝑢₁ and 𝑢₃, respectively. Without loss of generality, the turbulent wind field is thus reduced to a 2N-variate stationary Gaussian process.

Equation (2) is implemented as the simulation formula of the WAWS method, while (4) is adopted in the POD-based technique. Both in the WAWS and the POD methods, FFT algorithms are employed, as discussed in [1, 11], to improve their computational efficiency. An AR (𝑝₁) filter, given by (5) with 𝑞=0, is considered to represent autoregressive methods assuming a sufficiently large 𝑝₁. A two-step approach is adopted in the implemented codes, following the diagram reported in Figure 1. Accordingly, an initial phase is established (indicated as “phase 0”), which must be run only once, after which a general number n of simulations (Sim) can be performed (phase 1). Typically, in the “phase 0”, all the data which are needed for the successive simulations (e.g., spectral matrix, factorizations, coefficients, etc.) are calculated and collected. This allows minimizing the number of operations when simulating several wind realizations, which is usually the case for performing Monte Carlo simulations of the structural response. The results of the “phase 0” are deterministic and must be stored in the computer memory.

Figure 1 

Adopted computational scheme for wind simulation.

The two main aspects that should be considered when comparing different wind simulation techniques, from the computational point of view, are the allocated memory in the “phase 0” and the time needed for each single simulation 𝑇_sim in the “phase 1”. The former aspect is crucial since collecting a lot of data in the computer memory may cause overflows and thus it may preclude the use of a certain procedure when the dimension of the simulation domain becomes large. The latter aspect is even more important since it gives the measure of the computational effort required by the considered algorithm. On the other hand, when the accuracy of the different methods is concerned, the most relevant points that should be considered are the agreement between target and simulated characteristics of the stochastic process, mainly spectra and correlation functions, and peak values.

The three methods described above are implemented as schematically outlined in the pseudocodes provided in Appendix B. Accordingly, the main computational steps are summarized in Table 1. The comparison between the operations required by the considered algorithms immediately outlines that the WAWS technique is the heaviest one from the computational viewpoint. On the contrary, both POD and AR techniques significantly reduce the computational effort with respect to the WAWS method. Table 1 also summarizes the memory that is allocated in the “phase 0” in each of the three implemented codes. It is worth underlying that the ratio between the memory allocated in the POD-based method and the one allocated in the WAWS method is almost equal to 𝑁_𝜔/𝑁_𝜔. In the next section it will be shown that the POD-based algorithm gives sufficiently accurate results with 𝑁_𝜔=50, if compared to the corresponding WAWS method with 𝑁_𝜔=2¹⁴. The POD-based technique allows therefore drastically reducing the memory required by the WAWS method. This practically eliminates the risk of memory overflows in many technical conditions. Nonetheless, in the most demanding cases, an AR(p₁) filter can be adopted as it further reduces memory allocation even with respect to the POD-based method. Indeed, the ratio between the memory allocated by the AR method and the one allocated by the POD method is equal to(𝑝₁+1)/𝑁_𝜔 and, usually, 𝑝₁+1≪𝑁_𝜔.

Table 1 

Two-steps implementation of the considered simulation methods and corresponding memory allocation for a 2N-variate process (the symbol ∝ denotes proportionality).

5. Numerical Tests and Discussion

5.1. The Case Studies

The above described wind simulation techniques are commonly employed for the definition of the external inputs for time domain response analyses of slender structures, as described in [25–28], among others. Time domain response analysis are essential for evaluating the structural safety of wind sensitive structures, whose dynamic behavior is strongly nonlinear, for instance, in the case of structural cables and long-span bridges [29, 30].

Two numerical examples are considered in order to evaluate the performances of the simplified models, when their parameters are chosen in such a way to obtain a similar accuracy with respect to the WAWS method. In the former example the simulation domain is composed by 9 significant points disposed along the height of the tower of a suspension bridge (see Figure 2(a)). This case is here adopted in order to calibrate the parameters of POD-based and AR techniques. In the latter example the simulation domain is represented by an entire suspension bridge. In this case, the spatial domain is composed by 83 nodes located on the mean plane of the bridge (see Figure 2(b)). In both cases, the direction of the mean wind velocity is orthogonal to the plane of the simulation domain and its modulus is assigned by means of the classic logarithmic profile. In particular, a mean velocity of 48.0 m/s is assumed at the top of the tower in the former example (node number 9 in Figure 2(a)) and a mean velocity of 40.1 m/s is assumed at the mid-span of the bridge in the latter example (node 70 in Figure 2(b)). The exponential decay coefficients that appear in the expression of the coherence function are assumed to be equal to those reported in [11], while the variances 𝜎²_{𝑢𝑗(𝑘)} of 𝑢_𝑗^(𝑘), for 𝑗=1,2,3, are calculated following [23]. A nil coherence is assumed between orthogonal turbulence components 𝑢₁ and 𝑢₃. The following parameters are adopted in the simulations: Δ𝑡=0.05s, 𝑁_𝑡=𝑁_𝜔=2¹⁴, Δ𝜔=0.00767 rad/s, 𝜔_𝑐=62.8 rad/s. The simulations have been performed with the aid of a Core2 Duo Intel processor.

(a)

(b)

Figure 2 

Grids for wind simulation: example 1 (a); example 2 (b).

5.2. Calibration of POD-Based and AR Methods

Before discussing the performances of the simplified formulas, the parameters of POD-based and AR techniques must be properly calibrated in order to achieve similar accuracies with respect to the WAWS method. This task is here accomplished in the case of the first example. Literature studies [7, 11] suggest to choose, in similar cases, a reduced number of factorizations (for POD-based technique only) 𝑁_𝜔=50 and an order of autoregression (for AR filter only) 𝑝₁=20. After performing several tests, it was found that similar values could be also adopted here, as it will be clearer in the rest of the paper by comparing the quality of the results obtained by means of the three methods. Therefore, the values 𝑁_𝜔=50 and 𝑝₁=20 are chosen here, thus also allowing a direct comparability with the results presented in other literature works. As an example, a realization of the turbulent wind field generated by means of the WAWS method is shown in Figure 3.

(a)

(b)

Figure 3 

Example 1: wind velocity field sample generated using the WAWS method: along-wind velocities (a); across-wind velocities (b).

The quality of the generated samples is analyzed in Figures 4–10. To this end, the turbulence velocities simulated in two points (number 5 and 9 in Figure 2(a)) of the spatial domain are considered. First of all, by virtue of the ergodicity of the process simulated by means of (2), it is expected that the WAWS method provides a perfect accuracy when one entire period 𝑇₀=2𝑁⋅2𝜋/Δ𝜔 (for a 2N-variate process) of the process is simulated. This result has been confirmed here as shown, for instance, in the results presented in Figure 4, which have been obtained by assuming 𝑁_𝑡=𝑇₀/Δ𝑡 and 𝑁_𝜔=2¹⁴. In this case, as expected, target and computed spectra and correlation functions are overlapped (some meaningless residual differences are related to numeric FFT calculations). Figures 5 and 6 refer to along-wind and across-wind (vertical) velocities, respectively, simulated by means of the WAWS method assuming 𝑁_𝑡=2¹⁴, that is, 𝑁_𝑡≪𝑇₀/Δ𝑡. The corresponding results obtained by means of POD-based and AR techniques, for the chosen parameters values, are shown in Figures 7, 8, 9, and 10. From the presented results it points out that the WAWS method provides a very good match between target and computed spectra and correlation functions, even for 𝑁_𝑡≪𝑇₀/Δ𝑡. The chosen POD-based and AR methods are obviously less accurate than the WAWS formula, which is especially apparent by looking at the errors on auto- and cross-correlation functions. In particular, the analysis of the results, also extended to other points of the spatial domain, confirms that POD-based and AR techniques produce errors in the correlation functions, although, in general, the agreement with the target functions is improved as the coherence is decreased (as, e.g., in the across-wind components). Clearly, the desired levels of accuracy could be obtained in both cases of POD-based and AR algorithms by choosing larger values of either the reduced number of factorizations 𝑁_𝜔, adopted in the former case, or the finite order 𝑝₁ chosen in the latter one. For the purposes of this study, the parameters 𝑁_𝜔=50 and 𝑝₁=20 are considered as good compromises between accuracy and computational efficiency. Indeed, with such choices, the accuracies of the samples generated using POD-based and AR techniques are essentially comparable to each other and the errors with respect to the WAWS method can be regarded as acceptable in view of a structural response analysis [28]. In any case, the comments reported in the following developments of the work could be readily extended to other possible choices of the models’ parameters.

(a)

(b)

(c)

(d)

Figure 4 

Example 1 using WAWS method for 𝑁𝜔=214 and 𝑁𝑡=𝑇0/Δ𝑡 (i.e., a process duration equal to one period 𝑇0; black lines denote target functions, red lines denote numeric approximations): autospectrum of along-wind velocity at point 5 (a); corresponding autocorrelation function (b); cross-spectrum between along-wind velocities at points 5 and 9 (c); corresponding cross-correlation function (d).

(a)

(b)

(c)

(d)

Figure 5 

Example 1 using WAWS method (black lines denote target functions, gray lines denote numeric approximations): autospectrum of along-wind velocity at point 5 (a); corresponding autocorrelation function (b); cross-spectrum between along-wind velocities at points 5 and 9 (c); corresponding cross-correlation function (d).

(a)

(b)

(c)

(d)

Figure 6 

Example 1 using WAWS method (black lines denote target functions, gray lines denote numeric approximations): autospectrum of across-wind velocity at point 5 (a); corresponding autocorrelation function (b); cross-spectrum between across-wind velocities at points 5 and 9 (c); corresponding cross-correlation function (d).

(a)

(b)

(c)

(d)

Figure 7 

Example 1 using POD method (black lines denote target functions, gray lines denote numeric approximations): autospectrum of along-wind velocity at point 5 (a); corresponding autocorrelation function (b); cross-spectrum between along-wind velocities at points 5 and 9 (c); corresponding cross-correlation function (d).

(a)

(b)

(c)

(d)

Figure 8 

Example 1 using POD method (black lines denote target functions, gray lines denote numeric approximations): autospectrum of across-wind velocity at point 5 (a); corresponding autocorrelation function (b); cross-spectrum between across-wind velocities at points 5 and 9 (c); corresponding cross-correlation function (d).

(a)

(b)

(c)

(d)

Figure 9 

Example 1 using AR method (black lines denote target functions, gray lines denote numeric approximations): autospectrum of along-wind velocity at point 5 (a); corresponding autocorrelation function (b); cross-spectrum between along-wind velocities at points 5 and 9 (c); corresponding cross-correlation function (d).

(a)

(b)

(c)

(d)

Figure 10 

Example 1 using AR method (black lines denote target functions, gray lines denote numeric approximations): autospectrum of across-wind velocity at point 5 (a); corresponding autocorrelation function (b); cross-spectrum between across-wind velocities at points 5 and 9 (c); corresponding cross-correlation function (d).

5.3. Discussion

The results presented so far have permitted to choose the parameters of POD-based and AR methods that guarantee efficient simulations with an accuracy that is essentially comparable to that of the WAWS method. The second example is now worth considering to better compare POD and AR methods from the computational viewpoint.

The accuracy of the simulated samples is analyzed in Figures 11, 12, 13, and 14. To this end, to the along-wind and across-wind velocities at the points number 1 (top of one tower) and 70 (bridge mid-span), indicated in Figure 2(b), are considered. The presented results are analogous to those obtained in the first example, thus indicating that the quality of the simulated process is, as expected, independent on the dimension of the simulation domain. In particular, the results confirm that POD-based and AR methods, for the chosen values of 𝑁_𝜔 and p₁, essentially guarantee similar accuracies. The samples generated by means of the three algorithms are also in good agreement as it concerns the peak values 𝑢₁^(𝑘)_max and 𝑢₃^(𝑘)_maxof 𝑢₁^(𝑘) and 𝑢₃^(𝑘). As examples, some relevant results are summarized and compared in Table 2 for the two numerical examples.

Table 2 

Comparison between peak values of samples generated using the three simulation methods.

(a)

(b)

(c)

(d)

(e)

(f)

Figure 11 

Example 2 using POD method (black lines denote target functions, gray lines denote numeric approximations): auto-spectra of along-wind velocities at points 1 and 70 (a), (b); corresponding autocorrelation functions (c), (d); cross-spectrum of along-wind velocities in points 1 and 70 (e); corresponding cross-correlation function (f).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 12 

Example 2 using POD method (black lines denote target functions, gray lines denote numeric approximations): auto-spectra of across-wind velocities at points 1 and 70 (a), (b); corresponding autocorrelation functions (c), (d); cross-spectrum between along-wind velocities at points 1 and 70 (e); corresponding cross-correlation function (f).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 13 

Example 2 using AR method (black lines denote target functions, gray lines denote numeric approximations): auto-spectra of along-wind velocities at points 1 and 70 (a), (b); corresponding autocorrelation functions (c), (d); cross-spectrum between along-wind velocities at points 1 and 70 (e); corresponding cross-correlation function (f).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 14 

Example 2 using AR method (black lines denote target functions, gray lines denote numeric approximations): auto-spectra of across-wind velocities at points 1 and 70 (a), (b); corresponding autocorrelation functions (c), (d); cross-spectrum between along-wind velocities at points 1 and 70 (e); corresponding cross-correlation function (f).

The computer times required by the considered methods in the two numerical examples are summarized in Table 3. These results emphasize that, as expected, POD and AR methods are significantly more computationally efficient than the WAWS formula. Although it is implicit that these results strongly depend on the environment of simulation, they still allow to derive some conclusions. Particularly, the POD-based method seems to be more computationally efficient than the AR filter, both in performing preliminary calculations (“phase 0”), which could become relevant when computing only a few wind simulations, and wind simulations (“phase 1”). Moreover, the computational improvement in using the POD-based technique instead of the AR filter seems to become more relevant as long as the dimension of the simulation domain becomes larger. On the other hand, as expected, the AR filter is the method that mostly reduces memory allocation, at least in the case of high-dimensional simulation domains, which might be greatly beneficial in practical cases.

Table 3 

Computational efforts required by the three different methods for digital wind simulation (T₀ and 𝑇_𝑠𝑖𝑚 denote the computer times necessary for performing the “phase 0” and one single simulation in the “phase 1”, resp.).

It is also worth mentioning that the eigenvalue analysis of the spectral matrix, performed in the POD-based method, reveals interesting properties of the wind field in view of a structural analysis. In particular, as it is well-known, eigenvectors and eigenvalues of the spectral matrix represent respectively the shapes and the associated powers of the 3N stochastic processes summed in (4). The eigenvalue analysis usually reveals that only a few of these spectral modes have large associated powers and that the eigenvectors corresponding to the largest eigenvalues are very similar to the fundamental mode shapes of the structure. Thus, the response of a slender structure is expected to be dominated by the contribution of these waves [9].

As an example, Figure 15(a) shows the first ten eigenvalues 𝛾_𝑟^(𝑘) (𝑟=1,…,10) multiplied by the frequency 𝑛 and plotted versus the circular frequency 𝜔, in the case of the second example. As inferable from this figure, the first few eigenvalues are, as expected, much greater than the others. It is also worth mentioning that, as shown in Figure 15(b), the eigenvalues of the spectral matrix are distinct and thus the lines 𝑛𝛾_𝑟^(𝑘) versus 𝜔 do not cross. The turbulence eigenvectors corresponding to the first six eigenvalues, calculated for a fixed value of the circular frequency (𝜔=0.07 rad/s), are shown in Figure 16. The figure evidences that the first turbulence eigenvector involves the whole bridge in the out-of-plane (along-wind) direction and it is very similar to the typical first structural mode shape of a suspension bridge. Similar observations can be made for the remaining eigenvectors shown in Figure 16. It is noteworthy that hybrid horizontal/vertical eigenvectors are not detected since a nil coherence has been assumed between along-wind and across-wind turbulence components.

(a)

(b)

Figure 15 

Example 2: first ten nondimensional eigenvalues (POD) of atmospheric turbulence as a function of the circular frequency ω (a); detailed view (b).

(a)

(b)

(c)

(d)

(e)

(f)

Figure 16 

Example 2: first 6 blowing mode shapes evaluated for the circular frequency value 𝜔=0.07 rad/s.

The presented results have indicated that the POD-based technique exhibits a significant computational efficiency, while the AR filter is especially convenient for preserving the allocated memory. Another important aspect that needs to be addressed in the viewpoint of practical applications is the rate of convergence of the accuracies of the AR filter and the POD-based technique towards the target, when increasing the relevant control parameters, that is, the order of autoregression, p₁, in the case of the AR filter, and the reduced number of factorizations, 𝑁_𝜔, in the case of the POD-based technique. This aspect is analyzed here, by considering the error, ε, between correlation functions of generated samples, 𝑅^generated_{𝑢1(𝑖)𝑢1(𝑗)} and 𝑅^generated_{𝑢3(𝑖)𝑢3(𝑗)}, and the corresponding targets, 𝑅^target_{𝑢1(𝑖)𝑢1(𝑗)} and 𝑅^target_{𝑢3(𝑖)𝑢3(𝑗)}. This error is here defined as:𝜀=^𝑁𝑁_{𝑖=1𝑗=𝑖}‖𝑅^generated_{𝑢1(𝑖)𝑢1(𝑗)}−𝑅^target_{𝑢1(𝑖)𝑢1(𝑗)}‖+^𝑁𝑁_{𝑖=1𝑗=𝑖}‖𝑅^generated_{𝑢3(𝑖)𝑢3(𝑗)}−𝑅^target_{𝑢3(𝑖)𝑢3(𝑗)}‖ (6) and calculated in the time lag interval 0–10 sec. The error defined in (6) is realization dependent and, so, an averaging of ε among a conveniently large number of generations 𝑁_gen is sought, where 𝑁_gen is here taken as 200.

Following the previous definitions, Figure 17 shows plots of the average errors versus p₁ and 𝑁_𝜔, for the AR filter and the POD-based technique, in the case of the first numerical example. For clarity of presentation, the results in Figure 17 are normalized with respect to the residual error of the WAWS method with 𝑁_𝑡=𝑁_𝜔=2¹⁴. The results presented in Figure 17 confirm the expected trends that the errors decrease with increasing 𝑝₁ and 𝑁_𝜔. Moreover, they show that, in the case of the POD-based technique, the rate of convergence appears to be very fast for small values of 𝑁_𝜔. Namely, already for 𝑁_𝜔=2⁴ the system reaches an accuracy which is only slightly worse than that of the WAWS method with 𝑁_𝑡=𝑁_𝜔=2¹⁴ (the normalized average error is about equal to 1.3), while for larger values of 𝑁_𝜔 the convergence rate becomes smaller with some erratic fluctuations of less than ±5%. This result indicates that the interpolation expressions of eigenvectors and eigenvalues along the unequally spaced sequence 𝜔_𝑘 of circular frequencies, proposed by Carassale and Solari [11] and recalled in Appendix A, is very accurate when applied in the case of turbulence characteristics similar to those considered here. The results presented in Figure 17 might also suggest to choose values of 𝑁_𝜔 that are smaller than the one adopted by Carassale and Solari [11] in a similar application and chosen in the numerical examples presented above (𝑁_𝜔=50), without substantially decreasing the resulting accuracy. Nonetheless, calculations performed aside have shown that the average error in predicting peak values, calculated with respect to the WAWS method and considering 200 wind generations, rapidly decreases for 𝑁_𝜔 up to about 50, where it is almost equal to 0.5%, while remaining essentially constant for larger values of 𝑁_𝜔. This result shows that, in the present case, choosing 𝑁_𝜔≅50 is convenient in terms of accuracy in predicting peak values.

(a)

(b)

Figure 17 

Normalized average errors of simulated samples using AR filter (a) and POD-based technique (b) by increasing the order of autoregression (AR filter) and the reduced number of factorizations (POD-based method).

In the case of the AR filter, the results presented in Figure 17 show that for large values of p₁ the error is still slightly larger than the error of the WAWS method with 𝑁_𝑡=𝑁_𝜔=2¹⁴ (the ratio between the two is about 2.5). On this respect, a value of p₁ equal to 200 can be regarded as large if compared to the value 𝑝₁=20 recommended in the literature in similar conditions [7] and adopted in the numerical examples presented above. Looking at the results presented in Figure 17, the ratio between the average error of the AR filter and the residual error of the WAWS method, for 𝑝₁=20, is about equal to 3.2. Considering the high accuracy of the WAWS method, this result can be regarded as acceptable in many practical applications. It is also worth mentioning that, for 𝑝₁=20, the average error in predicting peak values is also small and approximately equal to 2.5%.

It is important to note that the results presented in Figure 17 do not allow to derive conclusions on the relative accuracy of the two simplified methods one with respect to the other. Indeed, for a fair and complete comparison between the two methods, different turbulence characteristics should be considered, which, however, goes beyond the purposes of the present investigation.

6. Conclusions

Simulating the wind velocity field in the case of high-dimensional domains may become a very demanding computational task due to memory occupation and time consuming simulations. In the present paper, a properly organized way for implementing wind simulation methods is presented at first. Then, the widely adopted method using waves superposition (WAWS method) is taken as a reference for calibrating the parameters of two well-known simplified methods (POD-based and AR formulas). These last are then compared, by devoting a special care to algorithm structure and computational expense.

Two numerical examples, with increasing complexity, are considered in the comparative study. The results indicate that POD-based and AR techniques significantly reduce both computer time and memory allocation with respect to the WAWS method, at the expense of smaller accuracies. The POD-based method appears to be more computationally efficient than the AR filter when the relevant parameters affecting their results are tuned in such a way to achieve similar accuracies. On the other hand, the AR filter is greatly preferable in terms of reduction of allocated memory.

Concerning the accuracy of the simulated samples, the WAWS method is known to guarantee the ergodicity of the generations. As shown in the paper, this means that if an entire period of the process is simulated, the match between target and simulated samples is perfect. Even when considering shorter simulations, the accuracy of the WAWS method is seen to be very high, namely the differences between target and computed correlation functions are, for the considered parameters values, almost unnoticeable. The results presented in this paper also confirm that, although less accurate, both POD-based and AR techniques can provide good results if the models’ parameters are properly chosen. It must be also mentioned that, regardless the chosen simulation formula, the POD decomposition of the spectral matrix is always a worth task to be tackled, as it provides useful information about the stochastic wind process in view of a structural analysis. Finally, the rates of convergence of the accuracies of the two simplified approaches with increasing tuning parameters have also been investigated in order to give preliminary indications for choosing these parameters in practical cases.

Appendices

A. Additional Details on the Adopted Simulation Formulas

The rules proposed by Carassale and Solari [11] for reducing the number of eigenvectors and eigenvalues of matrix 𝐆_𝐮𝐮(𝜔_𝑘) allocated in the computer memory, are described, at first, in this section. Then, the derivations for calculating the coefficient matrices of an AR (𝑝₁) filter are briefly recalled.

Concerning the POD-based technique, the following sequence of circular frequencies must be introduced:𝜔_ℎ=2𝜋𝑁_𝑡(𝑁Δ𝑡2^𝑡)^ℎ−1/𝑁_𝜔−1,(ℎ=1,…,𝑁_𝜔). (A.1) Eigenvectors and eigenvalues of matrix 𝐆_𝐮𝐮, calculated for the circular frequency values reported in (A.1), are denoted by 𝛾_𝑟^(ℎ) and 𝜃_𝑟^(ℎ), respectively. Then, by defining the sequence 𝑠_𝑘:𝑠_𝑘=ln(𝑘−1)ln(𝑁_𝑡(𝑁/2)^𝜔𝑁−1)+1(𝑘=2,…,^𝑡2+1) (A.2) the eigenvalues 𝛾_𝑟^(𝑘) of matrix 𝐆_𝐮𝐮(𝜔_𝑘) are interpolated as:𝛾_𝑟^(𝑘)=(𝑠_𝑘−round⁺(𝑠_𝑘))𝛾^{(round−(𝑠𝑘}_𝑟⁾⁾+(1−𝑠_𝑘+round⁺(𝑠_𝑘))𝛾^{(round+(𝑠𝑘}_𝑟⁾⁾𝑁(𝑘=2,…,^𝑡2𝛾+1),_𝑟(𝑘)=𝛾^(𝑁𝑡_𝑟^+2−𝑘)𝑁(𝑘=^𝑡2+2,…,𝑁_𝑡), (A.3) where round⁺(⋅) and round₋(⋅) provide the upper and lower integer rounds of the argument respectively. Finally, the eigenvectors are approximated stepwise as:𝜽_𝑟^(𝑘)=𝜽^{(round(𝑠𝑘}_𝑟⁾⁾𝑁(𝑘=2,…,^𝑡2𝜽+1),_{(𝑁𝑡𝑟+1−𝑘)}𝑁(𝑘=^𝑡2+2,…,𝑁_𝑡), (A.4) where round(⋅) returns the closest integer to the argument.

The simulation formula of an AR(p₁) filter is obtained, by assuming 𝑞=0 in (5), as:𝐮(𝑡_𝑗)=^𝑝_𝑖=1Φ_𝑖𝐮(𝑡_𝑗−𝑖Δ𝑡)+𝐁₀𝐚(𝑡_𝑗). (A.5) In order to calculate the coefficient matrices Φ_𝑖 and 𝐁₀, the correlation matrix 𝖱_𝐮𝐮(𝜏) must be defined as:𝖱_𝐮𝐮(𝜏)=𝐸[𝐮(𝑡)𝐮^𝑇(𝑡+𝜏)]. (A.6) Assuming the shorter notation 𝖱_𝐮𝐮(𝑘Δ𝑡)=𝖱_𝐮𝐮(𝑘), the following hyper-system of algebraic equations can be obtained [12]:𝐊Φ₁Φ₂⋯Φ_𝑝1^𝑇𝖱=𝐛,𝐊=^𝐮𝐮(0)𝖱_𝐮𝐮(1)⋯𝖱_𝐮𝐮𝖱(𝑝−1)_𝐮𝐮(−1)𝖱_𝐮𝐮(0)⋯𝖱_𝐮𝐮𝖱(𝑝−2)⋮⋮⋱⋮_𝐮𝐮(1−𝑝)𝖱_𝐮𝐮(2−𝑝)⋯𝖱_𝐮𝐮𝖱(0)𝐛=_𝐮𝐮(−1)𝖱_𝐮𝐮(−2)⋯𝖱_𝐮𝐮(−𝑝)^𝑇 (A.7) Once the correlation matrices are calculated by means of inverse FFT algorithms applied to CPSD functions, the coefficient matrices Φ_𝑖 can be easily determined by solving the algebraic hyper-system (A.7). Then it is possible to calculate matrix 𝐁₀, by post-multiplying (A.1) by 𝐮(𝑡_𝑗)^𝑇 and taking the average, namely:𝐸[𝐮(𝑡_𝑗)𝐮^𝑇(𝑡_𝑗)]=^𝑝_𝑖=1Φ_𝑖𝐸[𝐮(𝑡_𝑗−𝑖Δ𝑡)𝐮^𝑇(𝑡_𝑗)]+𝐁₀𝐸[𝐚(𝑡_𝑗)𝐮^𝑇(𝑡_𝑗)]. (A.8) By definition of correlation matrix, (A.8) may be rewritten in the following form:𝖱_𝐮𝐮(0)=^𝑝_𝑖=1Φ_𝑖𝖱_𝐮𝐮(𝑖)+𝐁₀𝖱_𝐚𝐮(0). (A.9) From (A.9), it follows that the covariance structure of the process is preserved if the following condition is satisfied:𝐁₀𝖱_𝐚𝐮(0)=𝖱_𝐮𝐮(0)−^𝑝_𝑖=1Φ_𝑖𝖱_𝐮𝐮(𝑖). (A.10) Thus, the choice of matrix 𝐁₀ is not unique and a possible strategy is to assume 𝐁₀=𝖱_𝐚𝐮(0)^𝑇 and to obtain it through the Cholesky decomposition.

B. Implemented Algorithms

The wind simulation algorithms implemented in this study are schematically outlined in this section, for a 2N-variate process. The simulated wind field is denoted by the matrix variable field_u. Variables definitions and adopted functions are implicit from the presented algorithms.

B.1. WAWS Algorithm

Step 0. initialize wind simulation parametersfor 𝑟=1∶2𝑁  for 𝑘=1∶𝑁_𝜔  calculate 𝜔_𝑟(𝑘) as a function of r and k    for 𝑖=1∶2𝑁    for 𝑗=𝑖∶2𝑁     𝐺(𝑖,𝑗)=𝐺_{𝑢𝑢𝑖𝑗}(𝜔_𝑟(𝑘))    end   end   matr𝐺=𝐺+upper_triangular_part(𝐺)^𝑇   𝑇=Cholesky_decomposition(matr𝐺)^𝑇   matr𝑇(∶,𝑘,𝑟)=𝑇(∶,𝑟)  endendsave wind simulation parameterssave matr𝑇.

Step 1. load wind simulation parametersload matr𝑇field_𝑢=zero_matrix(2𝑁,𝑁_𝑡)time=[Δ𝑡,2Δ𝑡,…,𝑁_𝑡Δ𝑡]phi=2𝜋⋅uniformly_distributed_random_numbers   (2𝑁,𝑁_𝜔)for 𝑖=1∶2𝑁 for 𝑟=1∶𝑖  for 𝑘=1∶𝑁_𝜔   𝐵_𝑟𝑘(𝑘)=matr𝑇(𝑖,𝑘,𝑟)⋅exp(𝜄⋅phi(𝑟,𝑘))  end  add_field_𝑢=𝑁_𝜔⋅Inverse_FFT(𝐵_𝑟𝑘)  for 𝑗=1∶𝑁_𝑡/𝑁_𝜔  field_𝑢(𝑖,(𝑗−1)⋅𝑁_𝜔+1∶𝑗⋅𝑁_𝜔)  =field_𝑢(𝑖,(𝑗−1)⋅𝑁_𝜔+1∶𝑗⋅𝑁_𝜔)  +√Δ𝜔⋅real(add_field_𝑢.∗exp(𝜄⋅𝑟/2𝑁  ⋅time((𝑗−1)⋅𝑁_𝜔+1∶𝑗⋅𝑁_𝜔)⋅Δ𝜔))  end endendsave field_𝑢.

B.2. POD-Based Algorithm

Step 0. initialize wind simulation parametersfor 𝑘=1∶𝑁_𝜔 for 𝑖=1∶2𝑁  for 𝑗=𝑖∶2𝑁   𝐺(𝑖,𝑗)=𝐺_{𝑢𝑢𝑖𝑗}(𝜔_𝑘)  end end  matr𝐺_hat=𝐺+upper_triangular_part(𝐺)^𝑇 matrPSI_hat(∶,∶,𝑘) =matrix_eigenvectors(matr𝐺_hat) matrGAMMA_hat(∶,𝑘) =eigenvalues(matr𝐺_hat)endsave wind simulation parameterssave matrPSI_hat and matrGAMMA_hat

Step 1. load wind simulation parametersload matrPSI_ℎ𝑎𝑡 and matrGAMMA_hatfield_𝑢=zero_matrix(2𝑁,𝑁_𝑡)𝑅=1/√2⋅normally_distributed_random_numbers (𝑁_𝜔,2𝑁)𝐼=1/√2⋅normally_distributed_random_numbers (𝑁_𝜔,2𝑁)time=[Δ𝑡,2Δ𝑡,…,𝑁_𝑡Δ𝑡]for 𝑘=1∶𝑁_𝜔 if 𝑘≤𝑁_𝑡/2+1  𝑠𝑘=log(𝑘−1)/log(𝑁_𝑡/2)⋅(𝑁_𝜔−1)+1  PSI_𝑘=matrPSI_hat(∶,∶,round(𝑠𝑘))  GAMMA_𝑘=(𝑠𝑘−floor(𝑠𝑘))⋅matrGAMMA_hat    (∶,∶,ceil(𝑠𝑘))  +(1−𝑠𝑘+floor(𝑠𝑘))⋅matrGAMMA_hat  (∶,∶,floor(𝑠𝑘)) end if 𝑘>𝑁_𝑡/2+1  𝑠𝑘2=log(𝑁_𝜔+1−𝑘)/log(𝑁_𝑡/2)⋅(𝑁_𝜔−1)+1  PSI_𝑘=matrPSI_ℎ𝑎𝑡(∶,∶,round(𝑠𝑘2))  GAMMA_𝑘=(𝑠𝑘2−floor(𝑠𝑘2))  ⋅matrGAMMA_hat(∶,∶,ceil(𝑠𝑘2))  +(1−𝑠𝑘2+floor(𝑠𝑘2))⋅matrGAMMA_hat  (∶,∶,floor(𝑠𝑘2)) end add_field_𝑢=zero_matrix(2𝑁,1) for 𝑟=1∶2𝑁  add_field_𝑢=add_field_𝑢  +Δ𝜔⋅GAMMA_𝑘(𝑟,1)  ⋅PSI_𝑘(∶,𝑟)⋅(𝑅(𝑘,𝑟)+√−1⋅𝐼(𝑘,𝑟))  end  field_𝑢(∶,𝑘)=field_𝑢(∶,𝑘)+add_field_𝑢endfor 𝑟=1∶2𝑁field_𝑢(𝑟,∶)=𝑁_𝜔⋅Inverse_FFT(field_𝑢(𝑟,∶))endfield_𝑢=real(field_𝑢)save field_𝑢

B.3. AR Algorithm

Step 0. for 𝑖=1∶2𝑁 for 𝑗=𝑖∶2𝑁  for 𝑘=𝑁_𝜔/2+1∶𝑁_𝜔   matr𝑆(𝑖,𝑗,𝑘−𝑁_𝜔/2)=𝑆_{𝑢𝑢𝑖𝑗}(𝜔_𝑘)  end  matr𝑅(𝑖,𝑗,∶)=Inverse_FFT(matr𝑆(𝑖,𝑗,∶)) endendfor 𝑗=1∶𝑝₁+1matr𝑅(∶,∶,𝑗)=matr𝑅(∶,∶,𝑗)+upper_𝑡𝑟𝑖𝑎𝑛𝑔𝑢𝑙𝑎𝑟_𝑝𝑎𝑟𝑡(𝑚𝑎𝑡𝑟𝑅(∶,∶,𝑗))^𝑇endarrange 𝑚𝑎𝑡𝑟𝐾=𝐊 and 𝑚𝑎𝑡𝑟𝑏=𝐛 from 𝑚𝑎𝑡𝑟𝑅solve hyper-system (A.7) 𝑚𝑎𝑡𝑟𝑝ℎ𝑖=𝑚𝑎𝑡𝑟𝐾⧵𝑚𝑎𝑡𝑟𝑏𝐵0=𝑚𝑎𝑡𝑟𝑅(∶,∶,1)for 𝑖=1∶𝑝₁ 𝑝ℎ𝑖(∶,∶,𝑖)=𝑚𝑎𝑡𝑟𝑝ℎ𝑖(2𝑁(𝑖−1)+1∶2𝑁𝑖,∶)^𝑇 𝐵0=𝐵0−𝑝ℎ𝑖(∶,∶,𝑖)⋅𝑚𝑎𝑡𝑟𝑅(∶,∶,𝑖+1)end𝐵0=𝐶ℎ𝑜𝑙𝑒𝑠𝑘𝑦_𝑑𝑒𝑐𝑜𝑚𝑝𝑜𝑠𝑡𝑖𝑜𝑛(𝐵0)^𝑇save 𝑝ℎ𝑖 and 𝐵0

Step 1. load wind simulation parametersload 𝑝ℎ𝑖 and 𝐵0𝑎=normally_distributed_random_numbers(2𝑁,𝑁_𝑡+𝑝₁)field_𝑢=𝑧𝑒𝑟𝑜_𝑚𝑎𝑡𝑟𝑖𝑥(2𝑁,𝑁_𝑡)for 𝑗=𝑝₁+1∶𝑁_𝑡+𝑝₁ for 𝑖=1∶𝑝₁  field_𝑢(∶,𝑗)=field_𝑢(∶,𝑗)  +𝑝ℎ𝑖(∶,∶,𝑖)⋅field_𝑢(∶,𝑗−𝑖) end field_𝑢(∶,𝑗)=field_𝑢(∶,𝑗)+𝐵0⋅𝑎(∶,𝑗)endfield_𝑢=field_𝑢(∶,𝑝₁+1∶𝑁_𝑡+𝑝₁)save field_𝑢.