Simulation of Nonstationary and Nonhomogeneous Wind Velocity Field by Using Frequency–Wavenumber Spectrum

The time history analysis is used to estimate the peak responses of structures subjected to stationary and nonstationary winds. The time histories of the fluctuating wind processes at multiple points can be simulated based on the spectral representation method for given target auto and cross power spectral density (PSD) functions. As the number of the processes of interest increases, the computation time for the simulation increases drastically. For the stationary homogeneous or nonhomogeneous wind fields, this problem can be overcome by using the frequency-wavenumber PSD function to simulate the stochastic propagating waves or fields. In the present study, a technique to simulate the amplitude modulated and frequency modulated nonstationary and nonhomogeneous stochastic propagating wind fields is presented. The technique relies on representing the nonstationary wind velocity by amplitude modulating a process that is time transformed from a stationary process. It is based on the established relations between the PSD functions of the nonstationary and of the stationary wind velocity. Simple to use and implement equations to carry out the simulation for one-dimensional line wind velocity field and two-dimensional nonstationary and nonhomogeneous wind velocity field are presented. The use of the developed technique and its adequacy is illustrated through numerical examples.


INTRODUCTION AND BACKGROUND
Structures such as tall buildings, long bridges, and wind turbines are sensitive to wind actions (Simiu and Scanlan, 1996;Strommen, 2010). The time history of wind velocity at a point on a structure can be represented by the sum of the mean and fluctuating wind components. The fluctuating winds can be characterized by using the power spectral density (PSD) function such as the Davenport spectrum, Kaimal spectrum, and von Karman spectrum. A common characteristic of these spectra is that they can be expressed in terms of reduced frequency or Monin coordinate, which is proportional to the turbulence length scale and inversely proportional to the mean wind velocity. The relation between two fluctuating wind velocity time histories at two spatial points is modelled using the coherence function (Davenport, 1967;Mann, 1994;Krenk, 1996;Simiu and Scanlan, 1996;Solari and Piccardo, 2001;Peng et al., 2018).
The estimation of the structural responses to the wind loading can be carried out in the frequency domain or time domain. It is well recognized that the use of the frequency domain approach to estimate the peak responses for simple structures can be very efficient for linear elastic systems (Simiu and Scanlan, 1996). The use of the time domain approach becomes increasingly popular because of the availability of powerful time-history structural analysis software and computer science and engineering advances. The use of time-domain analysis also facilitates the evaluation of nonlinear inelastic structural responses and reliability by considering wind loading (Vickery, 1970;Hong, 2004;Hong, 2016) or earthquake loading (e.g., Hong and Hong, 2007). For the time history analysis, the synthetic wind velocity fields need to be generated. The simulation of wind velocity can be carried out using one of the several available methods, including the spectral representation method (SRM) (Shinozuka and Jan, 1972;Li and Kareem, 1991;Shinozuka and Deodatis, 1996;Di Paola, 1998;Spanos and Zeldin, 1998;Huang and Chen, 2009), auto-regressive moving-average (ARMA) model (Samara et al., 1985); the empirical mode decomposition (Xu and Chen, 2004); proper orthogonal expansion (Chen and Kareem, 2005;Carassale et al., 2007;Solari et al., 2007), and application of wavelet transform (Gurley and Kareem, 1999).
For the fluctuating wind modelled as a stationary process defined by its PSD function, the most popular method to simulate the time histories of wind velocity is SRM. As the PSD function of fluctuating wind velocity is a function of the mean wind velocity (Simiu and Scanlan, 1996), the PSD function becomes dependent on both the frequency and time. In such a case, the coupling of time and frequency in the PSD function is such that the factorization of the time-frequency-dependent PSD function as a function of time and a function of frequency could be difficult (Hong, 2016) if it is not impossible. To simplify the task of simulating the nonstationary wind velocity, some researchers (e.g., Chay et al., 2006;Chen and Letchford, 2007;Li et al., 2017) assumed that the square root of the whole evolutionary PSD (EPSD) function can be treated as the amplitude modulation function. The assumption is very practical but may not agree with the requirement that the amplitude modulation for an evolutionary process must be a slowly varying function (Priestley, 1965). Hong (2016) proposed the use of a uniformly amplitude modulated and frequency modulated (AM/FM) process to model nonstationary fluctuating wind. In other words, the nonstationary fluctuating wind velocity is represented by amplitude modulating a process that is nonlinearly time transformed from a stationary process. For the fluctuating winds at multiple points in space, this model describes the fluctuating winds with advection according to their corresponding time-varying mean wind velocity.
Given a PSD or EPSD matrix of the fluctuating winds at multiple points in space, the time histories of turbulent winds at these points can be modelled using SRM (Shinozuka and Deodatis, 1996). As the number of the processes defined by the matrix increases, the use of SRM for the simulation of the processes is inefficient since it involves the decomposition of the matrix at many frequencies and the application of the fast Fourier transform (FFT) to improve the computational efficiency may not be possible unless approximations are considered (Li et al., 2017).
Rather than simulating nonstationary wind processes at multiple points, Deodatis and Shinozuka (1989) proposed an efficient technique to simulate the stochastic waves (or fields) for a given frequency-wavenumber PSD (FW-PSD) function, which is an extension of SRM (Di Paolo, 1998). The use of the FW-PSD function based simulation for multi-dimensional random field is much more efficient than the use of multiple points based simulation for vector processes since the former avoids the decomposition of the power spectral density function at many frequencies. For the simulation of the one-dimensional homogeneous wind field, Benowitz and Deodatis (2015) showed that the FW-PSD function could be obtained from the auto PSD and cross PSD (XPSD) Functions of homogeneous fluctuating winds for two points in space. They modified the technique given in Deodatis and Shinozuka (1989) by using FFT to gain additional computational efficiency and showed the adequacy of the simulated wind field for a bridge deck represented as a horizontal line-like exposed structure. Chen et al. (2018) and Song et al. (2018) investigated the simulation of homogeneous and nonhomogeneous wind fields in two spatial dimensions by considering that the inhomogeneity is originated from the height varying mean wind velocity. In their study, the coherence is considered to depend only on the separations in the horizontal and vertical directions or Euclidian distance but independent of the coordinate of the points in space. It is noted that the simulation of nonstationary and nonhomogeneous wind fields by using the FW-PSD function has not been addressed in the literature. The nonstationary fluctuating wind could arise from time-varying mean wind speed due to thunderstorm events, while the inhomogeneity could be caused by the topographic effects or characteristics of the winds that vary with the height or the coherence that depends on the spatial coordinate and separation. These problems are often encountered in evaluating wind-induced responses of bridges (King, 2003), transmission tower and tower-line systems (Momomura et al., 1997;Yang and Hong 2016;Yang et al., 2017), and tall buildings (Irwin, 2009).
The main objective of the present study is to develop a simple technique to simulate the AM/FM and nonhomogeneous stochastic propagating wind velocity fields. It is considered that the nonstationary fluctuating winds at multiple points are represented by AM/FM processes (Hong, 2016). The remaining sections of the present study are organized as follows. Simulation of an Amplitude Modulated Evolutionary Stochastic Field section summarizes the essential theoretical formulation given in Deodatis and Shinozuka (1989) to simulate an amplitude modulated evolutionary stochastic propagating waves or fields. Simulation of Nonstationary and Nonhomogeneous Propagating Stochastic Field section describes the theoretical development leading to the formulation for simulating AM/FM and nonhomogeneous propagating stochastic fields. Numerical Examples section presents the numerical example applications.

SIMULATION OF AN AMPLITUDE MODULATED EVOLUTIONARY STOCHASTIC FIELD
Consider that a stochastic propagating wave or field in time and n-dimensional space, Y( x . ), can be expressed based on May 2021 | Volume 7 | Article 636815 the evolutionary theory (Priestley, 1965;Deodatis and Shinozuka, 1989), is a 1 × (n + l) vector containing the frequency f (Hz) and n wave numbers k i (that has the unit of the inverse of the distance); f corresponds to time, ) is a random process in f . and has orthogonal increments; and the FW- Based on this amplitude modulated stochastic field model, Deodatis and Shinozuka (1989) (see also Shinozuka and Deodatis, 1996) × cos 2π f . j t , j 1 , ..., j n , I t , I 1 , ..., I n · x . + ϕ j t , j 1 , ..., j n , I t , I 1 , ..., I n where S( f . ) in this equation is two-sided PSD, f . j t , j 1 , ..., j n , I t , I 1 , ..., I n j t Δf , j 1 Δk 1 , ..., j n Δk n Θ(I t , I 1 , ..., I n ) Θ is an element by element multiplication operator, Δ f . Δf Δk 1 /Δk n , a symbol following Δ represents the increment of the symbol, and ϕ(j t , j 1 , ..., j n , I t , I 1 , ..., I n ) is the random phase angle between 0 and 2π.

SIMULATION OF NONSTATIONARY AND NONHOMOGENEOUS PROPAGATING STOCHASTIC FIELD Adopted Time-Frequency-Dependent Power Spectral Density Function
The application of Eq. 3 to simulate the wind field is simple, provided that S( x ) is ensured. For the nonstationary and nonhomogeneous propagating wind velocity fields, the problem is to find realistic and justifiable models to represent their spatiotemporally varying characteristics. Many of the commonly used PSD function of the fluctuating wind velocity at height z (m) above the ground surface could be expressed in the following mathematical form (Solari and Piccardo, 2001), where S 0f (f) denotes the PSD function, f (Hz) is the frequency, A 0 , B 0 , a 0 , b 0 , and c 0 are model parameters, ζ 0 fz/U(z) denotes the reduced frequency, and U(z) is the mean wind speed. The integration of S 0f (f) over the frequency domain equals the variance. An example of Eq. 5 is the Kaimal spectrum where a 0 0, b 0 1, and c 0 5/3 (Simiu and Scanlan, 1996). Based on the type of PSD function shown in Eq. 5 such as the Kaimal spectrum and the concept of time transformation (i.e., frequency modulation) (Yeh and Wen, 1990), the AM/ FM model for the nonstationary fluctuating winds proposed in Hong (2016) is adopted in the following. The adopted model considers that the wind velocity is represented by the sum of the mean wind velocity U(p j , t) and fluctuating wind velocity u(p j , t) at a point in space p j , where p j (x 1j , x 2j ) in which x 1j and x 2j represent the coordinates of the horizontal and vertical axes, respectively. u(p j , t) is given by, and, where the time-varying standard deviation of the nonstationary process σ(p j , t) represents the amplitude modulation function, τ j is introduced to represent the nondimensional time that is used to define the frequency modulation, andũ j (τ j ) is a stationary process in the domain of τ j , denoted as τ-domain, (but nonstationary in the t-domain). Based on Eq. 5, it is considered that the PSD function Sũ ,jj (ζ) equal to Sũ(ζ) (one-sided spectrum), where A 0 in this equation is a normalization constant to ensure that the integration of Sũ(ζ) over nonnegative value of ζ equal to one. By considering a 0 0, b 0 1, c 0 5/3, A 0 2K/3, and B 0 K 50, Eq. 8 becomes, which represents the normalized Kaimal spectrum with the variance equal to one. The nondimensional quantity τ j viewed as the "time" variable could also be viewed as a "space" variable, depending on the preference of the user. Viewing as the "space" variable, the wind in this model is treated as a frozen wind with its time variation described by advecting the wind with a constant mean wind velocity in the τ-domain. It is homogeneous except for amplitude modulation. The PSD function of u(p j , t), S u,jj (f , t), is given by, where τ ′ j (t) U(p j , t)/x 2j is the time derivative of τ j (t) and f (Hz) is the frequency. This indicates that u(p j , t) represents a nonstationary process with the amplitude modulation σ(p j , t) and frequency modulation f /τ ′ j (t) ζ, which is based on the time transformation defined in Eq. 7. In other words, u(p j , t) is obtained through modulating the process that is mapped from the stationary processũ j (τ j ). Since a key assumption of the evolutionary process (see Eq. 1 for definition) (Priestley, 1965) is that the amplitude modulation function is a slowly varying function, the treatment of the square root of S u,jj (f , t) shown in Eq. 10 as the amplitude modulation function may be questionable unless such a slowly varying characteristic of S u,jj (f , t) is validated. To overcome this, one could simulate the stationary processũ j (τ j ) and map it to the t-domain through the use of Eqs. 6, 7.
In order to deal with multiple processes, each representing the wind velocity at a spatial point, define two stationary processes, u j (τ) at p j andũ k (τ) at p k , with their XPSD function, Sũ ,jk (ζ), given by, in which Sũ ,jj (ζ) Sũ(ζ) Sũ ,kk (ζ), and C x1 and C x2 and are model coefficients. As will be seen, Eq. 13 represents the use of the weighted mean wind speed rather than the simple average mean wind speed. By mapping [ũ j (τ),ũ k (τ)] to [u(p j , t), u(p k , t)] using Eqs. 6, 7, it was shown that the XPSD function u(p j , t) and u(p k , t), S u,jk (f , t), is given by, where τ jk ′ (t) [τ j ′ (t) + τ k ′ (t)]/2. For j k, Eq. 14 reduces to Eq. 10. It can be shown that the lagged coherence for u(p j , t) and u(p k , t), c u,jk (f , p j , p k ), is given by, This function differs from the lagged coherence function given by Davenport (1967) (see also Simiu and Scanlan, 1996;Hong, 2016) on how the mean wind speed is averaged. For fixed points, the lagged coherence decreases as the mean wind velocity increases. The lagged coherence also decreases with increasing frequency or separation. The vertical coordinate plays a minor role in controlling the lagged coherence.

Formulation for Simulating AM/FM and Nonhomogeneous Propagating Wind Field
As u(p j , t), j 1,. . ., n, shown in Eq. 6 that represent the nonstationary processes in nonhomogeneous field are obtained by mapping stationary processes in the nonhomogeneous field, they can be simulated based on SRM as shown in Hong (2016). It requires to carry out Cholesky decompose of the n × n spectral matrix formed by Sũ ,jk (ζ) in a sequence of frequencies, identify the sub-processes based on the lower diagonal matrix, simulate the stochastic processesũ j (τ) as the superposition of the subprocesses and mapũ j (τ) into u(p j , t). This can be efficient and effective if n is small and a large number of samples of u(p j , t) is to be simulated because the decomposition of the spectral matrix only needs to be carried out once. For n more than a few hundred or thousand, the computing cost of decomposing the spectral matrix at multiple frequencies increases drastically. In other words, its use to obtain wind histories is inefficient for a very large n such as the case where the high spatial fidelity alongwind velocity field acted normal to a long span bridge deck or a horizontal line-like exposed structure. Note that by assuming that the value of x 2j for all points of interest equals z, Eqs. 11-13 becomes, where the use of the one-sided PSD function ofũ j (τ) at any point in the space equal to Sũ(ζ) is made, and positive ζ is considered. The FW-PSD function, S FWũ (ζ, k 1 ), is obtained by applying the Fourier transform to Eq. 16, (17) By carrying out the integration, this equation simplifies to (Gradshteyn and Ryzhik, 1994;Benowitz and Deodatis, 2015;Chen et al., 2018), If instead of considering the lagged coherence shown in Eqs. 12, 13 for x 2j and x 2k equal to z [i.e., Davenport coherence function cũ ,jk (ζ, p j , p k ) exp(−ζ × |C x1 s/z|) as shown within Eq. 16], one could also use other available coherence functions. In particular, by using the coherence function given by Krenk (1996), where ζ z (2πζ) 2 + (z/1.34L u ) 2 /z and L u is the integral length scale and |s| denotes the distance between P j and P k , Eq. 17 becomes, which can be solved analytically as shown in Benowitz and Deodatis (2015).
The simulation of the stationary homogenous wind field u(p, τ) defined by the FW-PSD function shown in this equation for x 2j z can be carried out using Eq. 3, resulting in, The simulated nonstationary wind field u(p, t) x 2j z is then obtained for x 2j z through the mapping defined by Eqs 6, 7 which is re-written as, where I(p, t) is the time-varying turbulence intensity (i.e., coefficient of variation) defined as σ(p, t)/U(p, t).
As shown in Benowitz and Deodatis (2015), Eq. 21 can be rewritten by using FFT and its inverse for computational efficiency. In such a case, the sampled values ofũ(τ, x 1 ) in the τ-domain are given for evenly spaced τ values. This leads to the mapped values of u(p, t) in the t-domain are given for evenly spaced t values only if the time transformation between t and τ is linear (i.e., the mean wind velocity is time-invariant). Otherwise, if the samples of u(p, t) for evenly spaced t values are of interest and the transformation between t and τ is nonlinear, one could only take advantage of FFT and its inverse for a fixed τ value, which is a function of t. This leads to, where the notation τ(t) instead of τ on the right-hand side of the equation is used to emphasize that it is a function of t, and the operators FFT( ) and IFFT( ) denotes the FFT and its inverse, the subscript associated with FFT and IFFT indicates the domain that the operation is carried out. The random field characterized based on Eqs. 11-13 is nonhomogeneous in space because its lagged coherence depends not only on the separation between p j and p k along the horizontal and vertical directions but also on the vertical coordinates of p j and p k . This is the case for the random field characterized based on Eqs. 14, 15 as well. Also, Eq. 15 indicates that (1/x 2j + 1/x 2k ) represents the sum of the weights used to average the mean wind velocity. Therefore, to avoid the difficulty in dealing with this position-dependent inhomogeneity, as a simple pragmatic approach to simulate two-dimensional wind velocity field, one could replace (1/x 2j + 1/x 2k ) by 2/z av , where z av represents a weighted average height for a considered exposure area of a structure for a wind engineering application where the vertical dimension is much smaller than the horizontal dimension. For example, one could consider z av equal to the height of the geometric center of the exposure area or the height of the point where the resulting static wind load for a structure is applied. By using this weighted average, h(p j , p k ) in Eq. 13 is replaced by, and the velocity field defined based on Eqs 11-13 becomes orthotropic in space. In such a case, the FW-PSD function, S FWũ (ζ, k 1 , k 2 ), is given by, By carrying out the integral, this equation is simplified to (Mantoglou and Wilson, 1982;Song et al., 2018), The simulation of the stationary wind fieldũ(p, τ) in twodimensions is then carried out using, S FWũ j t Δζ, I 1 j 1 Δk 1 , I 2 j 2 Δk 2 ΔζΔk 1 Δk 2 × cos 2π j t Δζτ + I 1 j 1 Δk 1 x 1 + I 2 j 2 Δk 2 x 2 + ϕ j t , j 1 , j 2 , 1, I 1 , I 2 The samples of u(p, t) is then obtained through the mapping defined by Eq. 21. Again, as the time transformation between t and τ is nonlinear for the nonstationary processes, for a given value of t, the evaluation of Eq. 28 can be carried out by taking advantage of FFT resulting in, where B I1,I2 j t , τ(t) S FWũ j t Δζ, I 1 j 1 Δk 1 , I 2 j 2 Δk 2 ΔζΔk 1 Δk 2 e i [ 2π ( jt Δζτ ) +ϕ ( jt ,j1,j2,1,I1,I2 )]

NUMERICAL EXAMPLES
Three examples are presented in this section. In the first example, the one-dimensional horizontal wind velocity field is simulated by considering the mean wind velocity that varies along the horizontal line because of the topographic effect. The second example is the same as the first one, except that the mean wind velocity is assumed to vary with time as well. Finally, the simulation of two-dimensional nonstationary and nonhomogeneous wind fields is carried out. In some cases, multiple simulation runs are carried out, and the auto and cross PSD functions and lagged coherence function are evaluated by using the simulated wind fields and their expected values are compared to the PSD and target coherence functions.
Example 1: Consider that the simulation of wind along the main span of a bridge across a canyon is of interest. For the simulation, it is considered that U(p, t) is time-invariant and equals U max p(x 1 ), where the maximum mean wind speed U max equals 40 m/s and the wind profile along the horizontal main span, p(x 1 ) [sin(πx 1 /450) + 7]/8, is illustrated in Figure 1. For the numerical analysis z 40 m, K 50 (see Eq. 9), the FW-PSD given by Eq. 18 with C x1 20, and I(p, t) 0.12 (see Eq. 22) are used. For the sampling, it is considered Δζ 1/600 [i.e., z/(TU max ) 1/600, T 600 s], N t 6,000, Δt 0.1 s, Δk 1 1/450 (i.e., 1/L where the main span length L equals 450 m), and N 1 512 (see Eqs. 23,24).
The obtained samples of the one-dimensional horizontal wind field at a few selected points are illustrated in Figure 2. Figure 2A illustrates the time histories for selected points, where the mean wind velocity is included in the plot, showing that the time histories of the wind velocity near the middle of the main span of the bridge are associated with greater fluctuating winds than those near the pylons. The plot of the line wind field for selected times that are depicted in Figure 2B shows that the line wind field for a given time follows the assigned horizontal wind profile. Finally, the propagating line wind fields are shown in Figure 2C.
For assessing the adequacy of the simulated fluctuating wind, the simulation is run 25 times. The average PSD function of the simulated fluctuating wind velocities at three selected points along the bridge deck are evaluated and shown in Figure 3. The PSD functions calculated from the simulated samples compares well to the target PSD function at the considered sites.
In addition to the above, the lagged coherence of the wind velocities for simulated records at three pairs of points is calculated. The evaluation of the lagged coherence is carried out in the τ-domain rather than the t-domain, so the calculated values are independent of time. A comparison of the calculated lagged coherence from the samples to the targets shown in Figure 4 indicates that the match is excellent for the target lagged coherence value greater than 0.4. For cases where the target lagged coherence is less than about 0.3, the lagged coherence from the simulated records is higher than that of the target value. This is expected since it is well-known that the lagged coherence even for the simulated white noise is in the order of about 0.2 (Bendat and Piersol, 1971).
Example 2: This example is the same as the previous example except that the mean wind is considered to be equal to U max p(x 1 )[sin(πt/600) + 5]/6, where U max equals 40 m/s. Typical samples of one-dimensional nonstationary and nonhomogeneous horizontal wind field at a few selected points are illustrated in Figure 5A  Example 1, spectral analysis carried out by using the simulated samples of the propagating wind field indicates that the auto and cross PSD function calculated using the samples in the τ-domain match well their target values.
Example 3. In this example, the two dimensional AM/FM and nonhomogeneous propagating wind fields are simulated for the area shown in Figure 6A. The bottom line of the area is 25 m above the ground surface, and the height of the area equals 30 m, and the width equals 450 m. The profile of the mean wind velocity is defined by, which is shown in Figure 6B.
For the fluctuating winds, the FW-PSD shown in Eq. 27 is applicable with C x1 20 and C x2 16, z av (25 + 30)/2 27.5 m and Sũ(ζ) given by Eq. 9 with K 50. The typical simulated wind time histories and fields for a few selected points and time by applying Eqs. 29, 30 are presented in Figure 7. As expected, the plots of the time histories in each row that are for a fixed horizontal location resemble well since their vertical separations are relatively small. As the separation in the horizontal direction increases the similarity of the time histories decreases, which is expected since the lagged coherence decreases with increased horizontal separation. The plots of wind fields in the figure illustrate that the application of the proposed technique can simulate the nonstationary and nonhomogeneous fluctuating winds with very refined mesh. The total wind velocity of the propagating wind field at a given time is simply equal to the mean wind component plus the fluctuating wind component (as shown in Figures 6, 7).

CONCLUSION
In the present study, a simple technique to simulate AM/FM and nonhomogeneous propagating wind fields is presented. The technique is based on the theoretical development indicating that the nonstationary wind velocity can be represented as amplitude modulated and frequency modulated stationary wind where the frequency modulation is achieved through the nonlinear time transformation.
Simple to use and implement equations to carry out the simulation for AM/FM one-dimensional line wind velocity field and two-dimensional wind velocity field are presented by incorporating the fast Fourier transform for computational efficiency. The use of the developed technique and its adequacy is illustrated through numerical examples.  Frontiers in Built Environment | www.frontiersin.org May 2021 | Volume 7 | Article 636815 9