Simulation of spatially correlated wind velocity timehistories
A method to simulate spatially correlated turbulent wind histories is implemented following [1,2].
Two possible vertical wind profiles and two possible wind spectra are implemented. The user is free to implement new ones. The wind cocoherence is a simple exponential decay as done by Davenport [3]. If the wind field is simulated in a grid, the function windSim.m should be used (cf. Examples 1 and 2). For a more complex geometry, such as a radial grid, the function windSim2.m has to be called with two inputs: The first one contains the wind properties, and the second one contains the coordinates of the nodes where wind histories are simulated (cf. Example 3).
The folder windSim.zip contains:
 1 input file INPUT.txt for Example1.m
 1 input file INPUT_MAST.txt for Example2.m
 2 input files windData.txt and circle.txt for Example3.m
 The function windSIm.m
 The function windSim2.m
 3 examples files Example1.m, Example2.m, and Example3.m
 The function coherence.m that computes the cocoherence.
 The function Lx.m that computes the integral length scales
Notes:
 Simulating the wind field in a high number of points with a high sampling frequency may takes a lot of time.
 The software Turbsim has inspired this file. It is more efficient, but may not be best suited for bridge engineering.
 This code aims to be highly customizable
 For an overview of real wind measurements, see [4].
[1] Shinozuka, M., Monte Carlo solution of structural dynamics, Computers and Structures, Vol. 2, 1972, pp. 855 – 874
[2] Deodatis, G., Simulation of ergodic multivariate stochastic processes, Journal of Engineering Mechanics, ASCE, Vol. 122 No. 8, 1996, pp. 778 – 787.
[3] Davenport, A. G. (1961), The spectrum of horizontal gustiness near the ground in high winds. Q.J.R. Meteorol. Soc., 87: 194–211
[4] Buffeting response of a bridge at the inlet of a fjord, E. Cheynet, J. B. Jakobsen, J. þ. Snæbjörnsson, 14th International Conference on Wind Engineering, 2015
6.0  Major update: The standard deviation of the wind fluctuations is now asked in the input files instead of the turbulence intensity. This offers a more realistic simulation, where the turbulence intensity decreases for increasing mean wind velocities. 

5.0  typo 

5.0   new Example, and updated description 

4.2  updated picture 

4.2  summary updated 

4.2  new title 

4.2   new example included 

4.1   Updated example (need 90 sec to run instead of 23 sec ) but it is more complete 

4.0  description updated 

4.0  description updated 

4.0   removed "2D" from title because it is not actual any more 

4.0   The crosswind component v has been added !


3.101   update following @Paolo's feedback 

3.1   I have replaced the name of the function 'cohere' with 'coherence' since 'cohere' is a Matlab function that already exist. 

3.0  description 

3.0  correct update now 

3.0   

3.0  typo 

3.0   title of example file is changed 

3.0   New definition of the vector frequency to optimize the computation time


2.0  dexcription updated 

2.0   

2.0   

2.0  HTML example introduced + update: one single function to run the whole simulation 

1.6  picture updated 

1.6   erratum at the end of Run_Me is corrected 

1.5   included verification.m for a preview of the output ( spectra and time series) 

1.4   

1.3   A description of the output of the simulation is included at the end of the .m file "Run_Me" 

1.2   description updated 

1.1  description updated 
Pengyuan Qi (view profile)
GuillaumeSab (view profile)
coolclw (view profile)
E. Cheynet (view profile)
Hi Lei Tong,
The frequency resolution (denoted f0) is indeed defined as the inverse of the total duration. The upper boundary of the wind spectrum is limited by the Nyquist frequency and the lower boundary is limited by f0. So your point is correct.
However, I would do the scaling very carefully. I personally do not apply it because the simulation is done by using the wind spectrum as a target. The scaling usually leads to larger discrepancies between the target and simulated spectrum. Here, the target wind spectrum include the turbulence intensity as a parameter, but since the signal has a limited duration, the variance of the signal calculated by integrating the PSD over frequency will be lower than the "true variance", which corresponds to the variance of a signal of infinite duration.
Lei Tong (view profile)
Thanks for your quick answer Cheynet. That explains it.
I think another reason is that we are calculating the time series by accumulation over limited frequency range with randomized phase. Of course some energy will be canceled in this process. Anyway, simple scaling will get what I want without harming the expected PSD characteristic.
Thanks again!
E. Cheynet (view profile)
Hei Lei Tong,
Thanks for the feedback.
The fact that the simulated wind field has a slightly lower turbulence intensity than expected is normal. The reason is that the signal you generate has a finite duration. Consequently, there exist an (unavoidable) bias for the estimation of the statistical moments. The bias will get close to zero for sufficiently large duration of the simulation. If you look at the power spectral density (PSD) of the turbulent fluctuations, you will see that the target and simulated PSD agree well. Therefore, I would use the PSD as a "validation tool" rather than the turbulence intensity.
I hope my message answer your question
Regards,
Etienne
Lei Tong (view profile)
Hi Cheynet,
Thanks for your great code. It really provide a simple and clear way of turbulent wind generation.
One question do arise when I'm trying to to use your script, however. If I choose one certain reference height, say 90m, and the relevant mean wind speed, say 10m/s and a stdU equals 1.2m/s, I would expect an output wind with turbulence intensity of 0.12 at the reference height with mean longitudinal wind speed 10m/s. But this is not the case with the current output, the out come turbulence intensity is always lower than expected.
Have I misused the code somehow. Can you help me to understand this? Thanks!
steven chan (view profile)
E. Cheynet (view profile)
Hi djr,
It is perfectly normal to produce new time series every time we run the simulation again. Wind field simulation belongs to Monte Carlo methods, which work with random input(s). Here, the random input is the phase between each frequency components of the time series. I am using rng(1) in the two examples for reproducibility purpose only. Unless we want to look at the examples files, there is no reason to use rng(1).
djr (view profile)
Thanks a lot! It really helped.
Just one small clarification please. Due to rng(1), the same inputs provide the same final results. However, say I have 3 heights; 10, 20 and 30 m above ground and as a result I'll get three wind time series at the end. Now if I add few more heights, the velocities at 10, 20 and 30 m will not be the same as before when I only had these 3 heights. Why is that? I hope this explanation wasn't confusing. :) Thanks a lot!
E. Cheynet (view profile)
Hi djr,
My best answer is an additional html example file I have joined to the submission. I hope it will help you.
djr (view profile)
Hello,
I am interested to simulate fluctuating component of wind speed, but I am not interested in bridge engineering application. Just a reconstruction of fluctuating component of measured wind speed at a mast.
Can you please help me to configure Grid Generation for that case?
Thanks,
DjR
E. Cheynet (view profile)
@Paolo: Hi ! Thanks for the feedback. I got the same issue if I define Nzz>1 but Zmin = Zmax. However there is no problem if Nzz>1 and Zmin<Zmax. It is not surprising since if Zmin=Zmax, the user specifies an horizontal line. However, if at the same time Nzz>1, it is in the contradiction with the definition of the horizontal line. I have added some lines of code to overcome the issue. Did it solve your problem ?
Paolo (view profile)
Hello,
really nice script. Only one issue: I tried to run it with Nzz > 1 and it stops because of a not positive definite Cholesky matrix.
Regards,
Paolo