Elements for numerical simulation of wind time series Elementos para la simulación numérica de series temporales de velocidad de viento

Dynamic characteristics of structures determine the methods to be used for their analysis and design against environmental loads such as wind. In some structures, the effect of the energy contained in gusts sequences must be attended, emphasizing the fluctuating components of the gusts that resonate with them. In these cases, dynamic methods in time domain are frequently used. Numerical simulation of wind loads is a tool that has been widely developed and applied. The aim of this study is to perform a literature review to determine the necessary parameters for simulating wind time series as well as the existing numerical methods for this purpose, emphasizing in tropical cyclone-prone region.


Introduction
The wind is represented by its movement in two different time scales: 1) movements related to the variations of the planetary weather systems and 2) gusts.Van der Hoven (Van der Hoven, 1957) explained this behavior by identifying two phases within the wind power spectrum.The first phase is associated to the first and highest peak within the spectrum, which he calls macro or mesometeorological peak, which shows an energy concentration period of approximately four days corresponding to the typical transit time of a fully developed weather system.The second important peak, called micrometeorological, contains energies with periodic fluctuations around 1 minute up to tenths of seconds, and it is caused by turbulences caused by topographic effects: land roughness or obstacles around the site.The frequency range between the two peaks lacks energy, and it is known as a spectral gap.The greatest relevance of this gap is that enables to represent the wind in terms of a mean speed depending on the height above the ground level, which reflects the synoptic variations only, with the superposition of time-dependent fluctuations in three directions, the so-called gusts..The type of analyzed structure also determines how that load is going to be assigned on the structure, either as a static force or as a time-dependent dynamic force, given the fluctuating nature of the action and, consequently, of the response.Some structures respond much less to the intensity of an individual gust than to the energy contained in gust sequences, particularly regarding the fluctuating components of the gusts that resonate with them.Therefore, dynamic methods in the time or frequency domains are recommended for the analysis and design of these structural types; and for those having also a strong non-linear behavior, the time domain is more reliable.
The numerical simulation of synthetic time series is a tool that has been widely developed and applied in the generation of wind time functions that can subsequently be applied as load in the structures.Two wind components, mean and fluctuation, have to be studied for the simulation, and their elements have to be very carefully defined to accurately reflect the large and small-scale movements of the wind.The objective of this paper is to define a theoretical framework to determine the simulation parameters and methods of wind speed time series, differentiated according to the two components, with emphasis on hurricane-prone regions.

Mean wind speed. Vertical profiles
The wind's mean speed component is mainly conditioned by temperature and frictional interactions with the ground, due to the roughness and orography, which try to slow down the airflow by forming a vertical profile.
Several analytic models can describe the vertical profile of speeds and the three most known and used are: the logarithm law, the power law, and the Deaves and Harris model.According to Tamura and Kareem (Y. Tamura and Kareem, 2013), the log law is the most adequate for strong wind systems; however, it offers some mathematical complications that engineers, for practical reasons, prefer to avoid by using the power law.Studies by Yaojun et al (Yaojun Ge, Xinyang Jin and Cao, 2010) compare the Chinese standard with a total of 15 codes and design standards addressing wind loads in several countries of the region, referred by them as "Asian Pacific" economies.They reaffirm that 87% of the studied codes prefer and use the power law (Canada, China, Hong Kong, India, Indonesia, Japan, Korea, Malaysia, Philippines, Taiwan, Thailand, United States and Vietnam).
The adjustment of the behavior of the mean wind to previous profiles in tropical cyclone-prone regions has been studied.Among the main questionings of the scientific community is the comparison of profiles close to the cyclone eyewall with those of the outer regions of the center vortex of the eye.In 2009, Vickery et al. (Vickery, Wadhera, Powell and Chen, 2009) studied this radial dependence and determined that, in lower areas close to the ground level (above 100 m), the variation of the wind's mean speed is approximately logarithmic, and as the height increases, the applicability of the log law is broken and wind velocities start to reduce; a similar behavior is observed in the inner and outer areas of the cyclonic vortex.The authors propose a new empirical model, thereby recognizing that it properly adjusts to the sea surface, but they recommend extending the studies in order to apply them inland.In 2012 and 2013, Giammanco et al. (Ian M. Giammanco, Schröder and Powell, 2012;Ian M. Giammanco and Schröder, 2013) obtained results that are quite adjustable to those of Vickery et al. (Vickery et al., 2009).Cao et al. (Yukio Tamura, Cao and Giang, 2012) studied the vertical profiles during the path of the Maemi typhoon across Japan, concluding that the height variation of the mean wind speed in the regions close to the surface can be obtained through the classic log law or power law.They also confirmed that no significate differences were observed in the profile's characteristics on the inner or outer area of the eye of the typhoon.Based on researches of Vickery et al. (Vickery et al., 2009) and Giammanco et al. (Ian M. Giammanco et al., 2012;Ian M. Giammanco and Schröder, 2013), Tse et al. (Tse, Li, Chan, Mok and Weerasuriya, 2013) use the power law, the log law and the empirical profile obtained by Vickery et al. (Vickery et al., 2009) to adjust the mean profiles of typhoons observed from 2007 to 2009 in Hong Kong, at averages of 1 h and 10 min.According to their results, they do not recommend using the empirical law of Vickery et al. (Vickery et al., 2009) for 10min averages and heights over 300 m; on the contrary, both the power law and the log law gave adequate results for these same conditions.

Wind speed fluctuations
The atmospheric turbulence  !,  , ,  , is a stochastic process involving three variables (3)(, , ) and four dimensions (4) ,  , ;  .It is usually assumed as a stationary random field, with a media of zero, with Gaussian behavior and time-dependent; this review accepts this behavior.It can be assumed that stochastic components (, , ) are statistically independent from each other, changing from a 3 − 4 process towards a 3(1 − 4) .Through the discretization of the space domain in N points, sites where the wind has an impact on the structure, each component of the turbulence (, , ) can be considered a stochastic vector.Each component vector is a stochastic field (1 − 1) (representing one of the components of the wind speed fluctuation in a specific location, dependent on the deterministic parameter t, which is correlated to other components of the same vector).Thereby, it is possible to change from a 1 − 4 process to three  − 1 process, each one corresponding to one of the space components of the turbulence.This method conduces to a simplified model, which disregards any link with another point close to the one considered; rather than taking into account the space distribution of the wind field, from a full probabilistic characterization of the stochastic field reached by the complex form of the cross power spectral density (CPSD).In general, to determine the dynamic longitudinal response of the action of the wind on the structures, the components that are orthogonal to the direction of the fluid (x) may be disregarded; thus, the CPSD is expressed for two points located on the plane perpendicular to the action of the wind ( − ) through Equation 1, assuming that the imaginary part of the cross spectrum is insignificant (Ubertini and Giuliano, 2010).
Where  (− !"  ) , is the coherence function, which quantifies the cross correlation of the analog components of the turbulence at different points of the stochastic field.
is the mean speed and  !y  ! is the mean speed and k y and k z are experimentally determined decay coefficients (Di Paola and Gullo, 2001).The latter differ for each turbulence component and depend on various factors such as height, surface roughness, atmospheric stability, among others.According to Chen et al. (Chen, Li and Xiang, 2014), factors normally range from 7 to 10, and according to Xu et al, they are between 6 and 26 (Xu et al., 2014).Solari and Piccardo (G. Solari and Piccardo, 2001) analyzed 39 reference sources to study the variability of decay coefficients, which allowed creating the summary Table 1 for different directions of the wind action and the turbulence speed components.
With regard to these coefficients in hurricane-prone regions, there are some works based on field measurements.The study of Tamura et al. (Yukio Tamura et al., 2012) is one of them, which shows values that do not differ greatly from those proposed by Solari and Piccardo (G. Solari and Piccardo, 2001) In Equation 1,  !! and  !! are spectral density functions determined at each point.Meteorology and wind energy use different mathematical formulas for them.
In 1948, as a result of a large experimentation developed in wind tunnels with isotropic and homogenous turbulent airflow, von Kármán (von Kármán, 1948) proposed the first recognized formula for  !() , (Equation 3), specifically formulated for the longitudinal component of fluctuations.
Where  = Based on the expression of von Kármán, other ones have derived.The Davenport spectrum (Davenport, 1961), the Kaimal-type spectrum (Kaimal, Wyngaard, Izumi and Coté, 1972) and the spectrum proposed by Solari and Piccardo (G. Solari and Piccardo, 2001) are examples thereof.Among them, the spectra of Solari and Piccardo (Ciampoli, Petrini and Augusti, 2011;Cottone and Di Paola, 2011;Ubertini and Giuliano, 2010), Kaimal (Hernandez, Bernal and Caracoglia, 2013;J. Li, Li and Chen, 2011;Pinheiro, 2010) and von Kármán (An, Quan and Gu, 2012;Fu, Wu, Xu, Li and Xiao, 2012;Q. S. Li, Lunhai Zhi and Fei Hu, 2010) are the most commonly used by the international scientific community.The standard of Australia and New Zealand (AS/NZS1170.2011) and ISO (International Organization for Standardization, 2009) use the von Kármán spectrum.And according to Tamura and Kareem (Y. Tamura and Kareem, 2013), and endorsed by a group of works which monitored structures subjected to the action of hurricanes (Fu et al., 2012;Q. S. Li et al., 2010;Yukio Tamura et al., 2012), this is the most adequate spectrum for describing the wind turbulence during this type of storms.

Turbulence parameters associated to the spectrum
Certain parameters derived from the spectrum in equation 3 deserve special attention, due to a bibliographic divergence regarding their determination.These parameters are mainly: the integral scale and the turbulence intensity, defining the terms variance, friction velocity and turbulence intensity factor for the latter.
The turbulence integral length scale,  !, defines the position of the spectral content of the turbulence.It is a measure of the average size of vortices in the fluid inside the atmospheric boundary layer.It depends essentially on the roughness length,  !, the height from the ground z, the mean speed,  and the parameters of Coriolis, which depend on the latitude of the observation site (WMO, 2010). ! is the height in relation to the ground, where the wind speed is theoretically zero when the profile is assumed as logarithmic.There are two types of approaches to determine this parameter inside the atmospheric boundary layer (L.Li, Kareem, Xiao, Song and Zhou, 2015;Q. S. Li et al., 2010): 1) the morphometric methods, which use the characteristics of the surface to determine the roughness and 2) the micrometeorological methods.The work of Li et al. (L. Li et al., 2015) reviews some of the main works dealing with each of these methods.
There are several ways of obtaining  ! that can be consulted in the references (L.Li, Xiao, Kareem, Song and Qin, 2012;Wang, Hu and Cheng, 2011).Solari and Piccardo (G. Solari and Piccardo, 2001), based on the study of approximately 20 bibliographical sources, propose equation 4 to determine the turbulence integral scale.Equation 4 has been used by several authors who recognize its validity (Carassale and Solari, 2006;Cottone and Di Paola, 2011;Torrielli, 2011).
Some works obtain  !from field measurements in the presence of typhoons (Shiau, 2000;Tamura et al., 2012;Wang et al., 2011).They report that  !increase as the mean speed increases, being higher than the integral scales for the rest of the turbulence components.
According to a study by Li et al. (L. Li et al., 2005),  ! are higher in typhoons than in hurricanes, and they associate this to the fact that basins and latitudes originating the typhoons and hurricanes are a potential source of variations.The mentioned paper makes a comparison between the experimentally measured values and the results derived from applying the formulation proposed by the US standard ASCE7-10, 2010, given by Equation 5.
Where , is equal to 97.54, 152.4 and 198.12 and  equal to 1 3 , 1 5 and 1 8 for exposure categories B, C and D respectively.They conclude that the formulation of ASCE  Solari and Piccardo (G. Solari and Piccardo, 2001).x, y, z: longitudinal, lateral and vertical directions of the wind action 3.0 3.0 0.5 10.0 6.5 6.5 10.0 6.5 3.0 presents good approximations for speeds lower than 25 m/s and underestimates the value according to the real measurements when the speed is higher than 30 m/s.The Japanese code AIJ-RLB of 2004 recommends the use of Equation 6to determine the integral scale (Fu et al., 2012;Q. S. Li et al., 2010).
According to the work of Li et al.(Q.S. Li, Lunhai Zhi and Fei Hu, 2010), which compares the results of finding  ! in field measurements of an instrumented structure in the city of Beijing and the empirical formulations proposed by AIJ-RLB (6) and ASCE ( 5), it is valid to apply the formulation by AIJ-RLB, due to its good adjustments.An et al. (An et al., 2012) and Fu et al. (Fu et al., 2012)

differ from this criterion by concluding that the expression of the Japanese code does not offer good approximations of the turbulence scale, compared with the measurements made during the path of the Muifa and Megi typhoons across the Chinese of Shanghai and Guangzhou, respectively.
Table 2 summarizes some of the formulations proposed by four of the main wind codes, in addition to the Japanese and the US codes.
The variance of the longitudinal turbulence component  !!=  * where  is a non-dimensional coefficient called turbulence intensity factor, depends on the same parameters as the integral scale.In order to obtain it, the terms turbulence intensity factor and friction velocity must be specified.
According to a large number of studies, which Solari and Piccardo (G.Solari and Piccardo, 2001) quote from Bietry, Sacré and Simiu of 1978, the former propose equation 7 to determine the value of the turbulence intensity factor in the longitudinal direction of the wind, which has been adopted as valid in several published researches (Ciampoli et al., 2011;Cottone and Di Paola, 2011;Ubertini and Giuliano, 2010).
During the path of tropical cyclones, the turbulence mix is strong, which entails higher turbulence intensity factors.Li et al. (L. Li et al., 2012) propose Equation 8to calculate the  , factor, obtained from measurements of the Hagupit typhoon, which struck the southern coast of China in 2008.This paper compares the results with those of Masters et al. (Masters, Tieleman and Balderrama, 2010), concluding that both studies are quite consistent.
In a more updated paper using a larger number of hurricanes, Li et al. (L. Li et al., 2015) propose a better adjustment of equation 8 through Equation 9.
As for the determination of  * , friction velocity, two approaches are frequently used: 1) extrapolations from the logarithmic profile and 2) the expression based on the turbulent exchange motion (L.Li et al., 2015;Wang et al., 2011).The second method is used when anemometers are recording measurements in the three directions.
Concerning the turbulence intensity,  =  !, it can be perceived that it tends to decrease as the mean speed increases, and in the hurricane-prone regions there is a tendency of greater intensity than in the extra-tropical regions.
While Fu et al. (Fu et al., 2012) studied the results from the instrumentation of a building in Guangzhou, China during the Megi typhoon, they concluded that the turbulence intensity values were consistent with those resulting from applying the Solari and Piccardo model to determine the variance (G.Solari and Piccardo, 2001), as well as the results of Li et al. (Q. S. Li et al., 2010) regarding observations made in Beijing.
The work of Li et al. (L. Li et al., 2015) compares the experimentally obtained turbulence intensity during the path of typhoons and hurricanes, with the expression proposed by ASCE7-10 (2010), given by Equation 10.Where c is 0.30, 0.20 and 0.15 for exposure categories B, C and D of the same standard.Based on the comparison, they concluded that ASCE7-10 ( 2010) offered good approximations up to 60 m above the ground for categories D and C, while for B it overestimated the value.These results differ from those obtained by the previous work of Li et al. (Q. S. Li et al., 2010), which compares real-scale measurements with empirical values resulting from the formulations proposed by ASCE ( 10) and .The comparison was made for the B land category, according to ASCE, and from a height of 47 m above the ground level.
is 0.27 for a urban land.
The authors of that study concluded that the experimental values did not properly adjust to any of the empirical profiles of the ASCE or AIJ-RLB codes; thus, their application to high structures is not valid.

Simulation of wind speed time for gaussian processes
There are several simulation methods for generating stationary stochastic fields with Gaussian behavior.Neural networks (Wu and Kareem, 2011) and wave transformations (Benowitz and Deodatis, 2015) are currently applied alternatives.However, the most used methods are classified according to the following diagram in Figure 1.
(AR) and (MA) models are particular (ARMA) cases (Rossi, Lazzari and Vitaliani, 2004).The work of Samaras et al. (Samaras, Shinzuka and Tsurui, was one of the first to introduce (ARMA) models in 1985.A special use of these methods is observed in modelling of complex structures and large simulation domains (Di Paola and Gullo, 2001;J. Li and Li, 2012;Rossi et al., 2004).The (ARIMA) models are an extension of the (ARMA) models, and they are used mainly in the nonstationary behavior of wind series.The main strength of the (ARIMA) models lies in their capacity to reveal temporary interdependencies in the series (Zhihua Zhang and Moore, 2015).
The second traditional approach is to use the superposition of trigonometric functions with random phase angles, "spectral representation methods".Basic ideas were originated in the 50's (Chen et al., 2014), but its implementation in structural dynamics begun to be studied in the 70's, with Shinozuka and Jan as precursors (Shinozuka and Jan, 1972).These authors developed their simulation formulas aimed at stochastic processes, including multivariable, multi-dimensional and non-stationary cases; but when the number of variables increased, it turned out to be computationally inefficient.In 1972, with the purpose of improving the computational efficiency, Yang (Yang, 1972) introduced that the sum of the trigonometric functions can be made through the Fast Fourier Transform (FFT) and proposes a formula to simulate envelope processes.In 1996, Deodatis (Deodatis, 1996) developed an algorithm that improved the simulation formulas established by Shinozuka, because it guaranteed the generation of ergodic series by optimizing the implementation of the FFT.However, the use of the FFT increases the computer storage demand, which is a particular problem for simulations of multivariable and multidimensional processes that require long data-processing times (Chen et al., 2014;J. Li, Li, He and Shen, 2015).Therefore, in order to improve the computational efficiency, techniques have arisen that require the decomposition of XPSD matrix.
In the 90's, Li and Kareem (Yosun Li and Kareem, 1991) introduced the concept of XPSD decomposition, with the aim of simulating random stationary processes, with further application possibilities in multivariable and non-stationary random processes (Y.Li and Kareem, 1997).Two decomposition techniques are mainly used: 1) the Cholesky Decomposition and 2) the Proper Orthogonal Decomposition (POD) (Chen et al., 2014;Ubertini and Giuliano, 2010).
The improved procedure, resulting from combining the studies of Shinozuka et al. (Shinozuka and Jan, 1972) and Deodatis (Deodatis, 1996), applies the Cholesky decomposition of the power spectral density matrix.Solari et al. (Giovanni Solari, Carassale and Tubino, 2007) review the use of POD.According to this work, its use in structural dynamics dates back to the 80's.The works of Di Paola (Di Paola and Gullo, 2001) and Carassale and Solari (Carassale and Solari, 2006) are some of the most quoted by the scientific community regarding the application of this technique.
Some studies have used this technique to establish a comparison between the results' precision and the computational efficiency in the same problem.In 2010, Ubertini and Giuliano (Ubertini and Giuliano, 2010) made two numerical examples.In the first, they simulated the wind field to study the behavior of a suspension bridge tower (in 9 nodes); therefore, they analyzed the AR model efficiency, of the Shinozuka and Deodatis method, and the POD, with the improvements proposed by Carassale and Solari (Carassale and Solari, 2006).The second example simulated the wind field for the entire bridge (83 nodes), using only the AR model and the POD.They concluded that the Shinozuka and Deodatis method is the most accurate, but in relation to the first example, it demands 10 times more computing memory and takes approximately 8 times longer.In the first case, the variations of the other two methods were not significant among them.As for the second, they concluded that the POD is more accurate and computationally efficient in relation to time, but the AR consumes 5 times less computing memory.
Concerning the analysis of flexible structures such as towers and tall buildings (Serrano, Mora and Salazar, 2014), transmission towers (Mercanti, Pizzutti, Aguirre, Fank and Möller, 2011;T. Zhang, Lin and Bai, 2013;Zhuoqun Zhang, Li, Li, Wang and Tian, 2013) or guyed towers (Bu, Law and Zhu, 2012;Clobes and Peil, 2011;Gani and Légeron, 2010), and in agreement with these results, a group of authors has been identified who prefer the use of the spectral representation method, through the Cholesky decomposition.However, if they need to simulate large wind fields in connection with structures like suspension or cable-stayed bridges, they select techniques such as POD (Carassale and Solari, 2006;Y. F. Huang and Yang, 2011) or linear filter models (Proppe and Zhang, 2015).
Based on the methods presented, several works have refined them (Chen et al., 2014;Guoqing Huang et al., 2013;J. Li et al., 2011) or combined them.Regarding the latter, the papers of Huang (G. Huang, 2015;Guoqing Huang, 2014) addressing the simulation of non-stationary processes are an example thereof.

Conclusions
From the bibliographic study dealing with the simulation of wind speed time series, it is concluded that, regarding vertical profiles of mean wind speed, both the power law and the log law can be applied to simulate the behavior of the mean wind speed component under storm conditions.As for the formulations to obtain the spectrum, the expression of van Kármán seems to be the most adequate to describe the wind turbulence during the path of strong storms, such as hurricanes.Due to the wide range of criteria associated to the determination the turbulence scale, the model of Solari and Piccardo (G. Solari and Piccardo, 2001) proposed.Moreover, concerning the parameter of turbulence intensity factor, a comparison between the application of formulas 7 and 9 is proposed.
In relation to the synthetic generation methods of wind speed time series, it is concluded that digital filtering methods for generating wind time series does not require large storage memories, unlike those using the FFT, but they are less accurate than the spectral representation method.From the two main XPSD decomposition procedures, in the spectral representation method, the Cholesky decomposition leads to higher computational costs related to memory and time, but it is more accurate in predicting the response induced by the wind.This is an important aspect for the analysis of flexible and strongly non-linear structures such as transmission line towers, guyed towers and cable-stayed bridges.The Proper Orthogonal Decomposition (POD) is adequate when it is necessary to densify the discretization of the speed fields.

;
the height measured from the ground to the point,  ! is the integral scale of the turbulence component,  !! is the variance of the component, and  * ! is the square of the friction velocity.

Table 1 .
Coefficients proposed by