For a more robust and time-efficient Matlab implementation, see https://se.mathworks.com/matlabcentral/fileexchange/68632-wind-field-simulation-the-fast-version.
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 co-coherence is a simple exponential decay as done by Davenport . 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 windSim.m has an optional parameter to include two inputs (cf. Example3.mlx): 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 submission 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 getSamplingPara.m which get the time and frequency vectors
- 3 examples files Example1.m, Example2.m, Example3.m
- The function coherence.m that computes the co-coherence.
- Simulating the wind field in a high number of points with a high sampling frequency may take a lot of time.
- This code aims to be highly customizable
A more straightforward version of the present submission has been used to simulate the turbulent wind load on a floating suspension bridge .
 Shinozuka, M., Monte Carlo solution of structural dynamics, Computers and Structures, Vol. 2, 1972, pp. 855 – 874
 Deodatis, G., Simulation of ergodic multivariate stochastic processes, Journal of Engineering Mechanics, ASCE, Vol. 122 No. 8, 1996, pp. 778 – 787.
 Davenport, A. G. (1961), The spectrum of horizontal gustiness near the ground in high winds. Q.J.R. Meteorol. Soc., 87: 194–211
 Wang, J., Cheynet, E., Snæbjörnsson, J. Þ., & Jakobsen, J. B. (2018). Coupled aerodynamic and hydrodynamic response of a long span bridge suspended from floating towers. Journal of Wind Engineering and Industrial Aerodynamics, 177, 19-31.
Cheynet, E. Wind Field Simulation (Text-Based Input). Zenodo, 2020, doi:10.5281/ZENODO.3823864.
Hi Chen Guang,
As you know, the aerodynamic admittance functions (AAFs) account for the deformation of the eddies, as they impinge the deck, which modify the wind load. So AAFs need to be accounted for if the dimensions of the deck are not small compared to the size of the gusts.
I think the AAFs are generally estimated in wind tunnel tests. If you don't have access to them, you can assume that they are equal to one at all frequency. It will probably lead to an overestimated response, but it is OK for a first estimate. If I remember well, AAFs are estimated using pressure sensors. If you work only with velocity data, it can get tricky to implement the AAFs obtained with pressure records as the coherence of pressure fluctuations is generally larger than the coherence of velocity fluctuations. In ref , which relies only on full-scale measurements of velocity fluctuations only, we used a simple second-order filter function to model the AAFs for the lift force. I suppose a similar approach can be used for the overturning moment. We had not used any AAfs for the drag. The reason is that the lateral dynamic motion of the deck is located in the low-frequency range, where turbulence length scales are much larger than the deck dimensions.
 Cheynet, E., Jakobsen, J. B., & Snæbjörnsson, J. (2019). Flow distortion recorded by sonic anemometers on a long-span bridge: Towards a better modelling of the dynamic wind load in full-scale. Journal of Sound and Vibration, 450, 214-230.
Hi E. Cheynet,
Thanks for your kind and quick reply. I would like to ask you another question. For a rigid bluff body, the aerodynamic force is obtained under the action of turbulence wind. If we know the wind spectrum and aerodynamic spectrum, how do we generally get the aerodynamic admittance function (side force, lift force, and overturning moment)?
Hi Chen Guang,
Thank you for the feedback! You raise an interesting question: the answer will depend on the way the integral length scales (ILS) are estimated. I am using here the definition of ILS based on the von Karman spectral model, where the ILS is usually estimated by fitting the normalized von Karman spectral model to the estimated spectra. This method provides different estimates from the one calculated using the cross-correlation function. If I remember well, this lack of consistency is discussed in ESDU 85020, which proposes a modification of the von Karman model. I am aware that the ILS is a crucial parameter in wind tunnel test. However, in the atmospheric surface layer, the estimation of the ILS is associated with large uncertainties. In page 178 of the book of Panofsky and Dutton (1984) , the authors even recommend avoiding using the ILS in "application to atmospheric data". From my (moderately short) experience with in-situ atmospheric measurements, I agree with Panofsky and Dutton.
 Panofsky Hans, A., & Dutton, J. A. (1984). Atmospheric Turbulence, Models and Methods for Engineering Applications. Jhon Willy & sons, New York.
Hi E. Cheynet,
It's a really nice job. If we have generated fluctuation wind speed using windSim, How to know whether the turbulence integral scale of the fluctuating wind speed is consistent with the target turbulence integral scale?
Hi Xinlong ,
Thank you very much for the feedback. You are right, there should not be a 20x20 matrix but a 20x1 vector. I have modified the function KaimalSpectrum adequately. I will upload a corrected version as soon as possible.
Hi E. Cheynet,
Thanks for sharing these very helpful codes. I found that the windSim.m file may have a bug in Lines 120-123, where the function KaimalSpectrum.m is called. I think we should use the transpose of nodes.Z instead of nodes.Z itself. Specifically, I think Line 120 should be
If we look at Example2, Su in Line 120 should be a 20x1 vector not a 20x20 matrix. I found that Su in the VonKarmanSpectrum is a 20x1 vector. I also compared the theoretical Kaimal spectrum and the simulated spectrum, which match each other if we modify the code as I said.
It would be great if you can check the code again. Thank you very much!
I apologize for the late reply. You can find the submission for one single point here: https://se.mathworks.com/matlabcentral/fileexchange/76854-one-point-random-process-generation
in a few words: These are empirical coefficients that quantify how strongly correlated is a turbulent wind field. If you want standards values for offshore site, you may want to look at the IEC 61400-3:2009 and IEC 61400-1:2005 standards. However, the values used for the decay coefficients were not found from offshore wind measurements (as far as I know). Regarding Cuz, Cvz and Cwz, you can find some values in  estimated from the FINO1 platform and 2 years of data. In ref  there is a shorter and more pedagogical description of the co-coherence for the same data set. However, the coherence model in ref  is different from ref. , so be careful when interpreting the results. Regarding Cuy, Cvy and Cwy I do not think there is any estimate based on measurements for offshore sites, unfortunately.
 Cheynet E, Jakobsen JB, Reuder J. Velocity spectra and coherence estimates in the marine atmospheric boundary layer. Boundary-layer meteorology. 2018 Dec 1;169(3):429-60.
 Cheynet, E., 2018, September. Influence of the measurement height on the vertical coherence of natural wind. In Conference of the Italian Association for Wind Engineering (pp. 207-221). Springer, Cham.
Thanks for this great piece of work!
I am interested in using your code to simulate the average wind profile for a lot of points in a grid. I am unsure of how to work out/where to find standard values for the co-coherence parameters: Cuy_1, Cuz_1, Cvy_1, Cvz_1, Cwy_1 and Cwz_1. Could you please shine some light on what exactly these values are and what values to use for an offshore site.
Thanks in advance!
If you could create that submissions I'd be grateful!
Yes, the former approach was easier to understand but so much slower. I got several emails from people trying to use this version to generate a wind field in large domain, which is not a good idea without the FFT-based method. I tried to redirect them toward the "fast version" of the code but apparently the slow-version was still the one they prefered. So I decided to keep the text-based version and implement the FFT algorithm inside windsim. I took the occasion to clean a little the code.
To learn a little more about the fft, you can learn a lot by simply playing with the fft and ifft functions with a Gaussian random noise in Matlab.
I can also create a submission where a random process is generated using the fft in a single point instead of an entire grid.
Really nice code, learned a lot about synthesizing wind speed data by means of a random process.
I could understand well what you did prior 7.4 version, after this update I dont get the IFFT process you did on the code.
Is there a source where i can learn about it? I mean besides the   in the description.
Hi Johannes. Firstly, I would suggest you look at the "fast version" of the turbulence generator, which is also available on File Exchange. It's faster and more robust. Secondly, a turbulence intensity that decreases with the height should already be included here because the variance of each velocity component is independant of the height. If you have a logarithmic profile, changing the variance of the velocity will not affect the mean wind speed profile.
I would like to take up Amanda's question again.
For my purpose, I need a turbulence intensity which decreases by increasing height overground. Hence, I would like to adapt the code by changing the T.I. to a given profile (as it is done for the mean wind profile).
So my question is now whether this small modification will change the computation of the wind profile?
Thanks in advance
In the text file, you can set the standard deviation of the velocity components. The turbulence intensity TI is defined as the ratio of the standard deviation of the velocity component by the mean wind speed. By changing the standard deviation values, you can modify the turbulence intensity of the generated wind field.
Hi! Is there anyway to edit the turbulence intensity?
I am not sure I can reproduce your result. Nevertheless, simulating a wind field in 2000 points in space using the "friendly version" of the turbulence generator is not necessarily a good idea because it will take a huge amount of time. I would suggest, instead, to use the fast version of the turbulence generator, available here: https://se.mathworks.com/matlabcentral/fileexchange/68632-wind-field-simulation-the-fast-version
Even with the fast version, using 2000 points in space with a duration of e.g. only 600 sec and a sampling frequency of 10 Hz will be very time-consuming. I have noticed that you generate your points in a 3D space. Since the turbulence generator relies on Taylor's frozen turbulence hypothesis (TFTH), you actually do not need to generate points in the along wind direction, i.e. along the x-axis. The reason is that the TFHT allows you to convert time separation into space separation. Therefore, I suggest you to simply generate the wind field into a 2D plan upstream of the turbine and to generate the other plan by simply introducing a time delay DT that is transformed into a "space delay" DX = U*DT.
I'm trying to simulate LIDAR measurements upstream of a turbine, and have therefore defined a grid consisting of several planes of measurements. That is to say, my input grid is a file containing 3D-points that I want measurements of that are both correlated in time and space. This is key since I'm trying to figure out how well my LIDAR-based predictions are when using a more realistic turbulence model. The file containing the grid looks like this
0.0000 29.8500 0.0000
0.0000 29.8500 19.8111
0.0000 29.8500 39.6222
0.0000 29.8500 59.4333
0.0000 29.8500 79.2444
0.0000 29.8500 99.0556
0.0000 29.8500 118.8667
0.0000 29.8500 138.6778
0.0000 29.8500 158.4889
0.0000 29.8500 178.3000
0.0000 39.2342 0.0000
177.7778 114.3079 158.4889
177.7778 114.3079 178.3000
177.7778 123.6921 0.0000
177.7778 123.6921 19.8111
177.7778 123.6921 39.6222
177.7778 123.6921 59.4333
177.7778 123.6921 79.2444
177.7778 123.6921 99.0556
177.7778 123.6921 118.8667
177.7778 123.6921 138.6778
... and so on
I'm inputting 2000 points where the system should simulate wind history, but the output I'm getting is just a Uy of size (20,322). This is weird because I gave it a lot more points to simulate, and the runtime of the function gets larger the more points I put in, but it does not match (20 != 200). It's also weird because my input file contains the following:
windSim Input File.
author: E. Cheynet
Valid for windSim2 v 2.0, 13.06.2018. Delimiter: tab
------------ Time definition Options ------------------------
fs 2 sampling frequency (Hz)
Duration 200 Duration of the time series (s)
------------ Turbulent Wind data Options ----------------------
Cuy_1 5 co-coherence decay coefficient for u-component
Cuz_1 8 co-coherence decay coefficient for u-component
Cvy_1 5 co-coherence decay coefficient for v-component
Cvz_1 8 co-coherence decay coefficient for v-component
Cwy_1 5 co-coherence decay coefficient for w-component
Cwz_1 8 co-coherence decay coefficient for w-component
roughness 0.001 roughness of the area (useful only if "log" profile is selected)
u_star 0.6 friction velocity (m/s)
betaU 2.4 Used in the 'Von Karman' model only: Ratio stdU/u_star, where stdU is the standard deviation of the along-wind velocity component
betaV 1.9 Used in the 'Von Karman' model only: Ratio stdV/u_star, where stdV is the standard deviation of the crosswind velocity component
betaW 1.3 Used in the 'Von Karman' model only: Ratio stdW/u_star, where stdW is the standard deviation of the vertical wind velocity component
------------ Mean Wind speed Profile ------------------------
type: VonKarman wind spectrum: choice between 'VonKarman' or 'Kaimal'
profile: power Mean wind profile: choice between : 'power' or 'log'
meanU 14 mean wind speed at reference height (useful only if "power" profile is selected)
zr 119 reference height (useful only if "power" profile is selected)
a 0.15 coefficient for power law (useful only if "power" profile is selected)
NOTE : Do not add or remove any lines in this file!
This should give 100 time points, but the output suggests it has simulated for 322 times.
Do you have any idea what is causing this?
Thanks in advance!
Hi Sebastian, I do not think that it is possible to assume that turbulence is "frozen" past the first wind turbine. Said differently, if you want to model the wind load on a wind turbine in the wake of another one using a "simplistic model", you can try first to model the mean wind profile downstream of the first wind turbine using some empirical model described in the literature (e.g. ), which gives the profiles as a function of the downstream distance. Then you can do the same for the profile of velocity variance. I am not sure there is any empirical model for the profiles of the normalized velocity spectra and coherence in the wake of a wind turbine, so you could try with a standard modelling through the Kaimal-type spectra and exponential-decay coherence functions. But there is no guarantee it gives realistic results.
 Stevens, R. J., & Meneveau, C. (2017). Flow structure and turbulence in wind farms. Annual review of fluid mechanics, 49.
Hi E. Cheynet,
Instead of another turbine I meant a Lidar that measures the shear profile lets say 100m before the turbine. I basically want to simulate a gust that first passes by the Lidar and then with a certain delay arrives at the turbine. If I understood correctly, in your model you assume a homogenous wind field in x-direction so I would just add a discrete gust and let it advect with the wind speed given by your model. Is this a somewhat realistic assumption?
A wind turbine located 100m downstream from another one will experience a wind field quite disturbed, both in terms of mean wind speed and in terms of turbuent velocities.
As an example from offshore measurements: https://doi.org/10.1016/j.egypro.2016.09.217. During this short "experiment" scanning lidar measurements were observed to be "disturbed" by the wake of a wind turbine located a couple of kilometers upstream of the scanned area. While empirical relations can be used to model the disturbed mean wind speed profile, it might be more challenging to model the disturbed turbulent field.
Thanks for your quick reply, that makes sense. I have a follow-up question: Let's say I calculate a wind field as defined in Example 2 at a certain location dependent on time and altitude. Is it reasonable to assume that a wind turbine located behind this point (couple of 100m) would experience approximately the same wind field but delayed proportional to the distance and divided by mean wind speed? I essentially want to model advection of a gust.
Don't worry, that is not a stupid question at all. I am not familiar with the Von-Karman turbulence model in the aerospace toolbox in Simulink. From the description, it seems that the limited dimension of the gusts is not accounted for in the aerospace toolbox. Said differently. the modelling of the "wind turbulence" in the aerospace toolbox may be limited to a spatial dimension of a couple of meters only. While this description may be good enough for a small plane, that is not the case for a bigger structure, which has a non-negligible size compared to a typical length scale of turbulence.
I'm new to wind field simulation so that's maybe a "stupid" question: What is the difference between this wind field generator and for instance the Von-Karman turbulence model in the aerospace toolbox in Simulink?
windSim is a function so you will not be able to run it like a script. The scripts files are Example1.m andExample2.m; these are the files you need to run to see how the function windSim works.
I am having trouble starting the example 1. Do I copy this code
rng(1); % to ensure reproducibility of the example.
filename = 'INPUT_Example1.txt';
[u,v,w,t,nodes] = windSim(filename);
or do I open the file named windsim.m?
The nodes represent the points in space where the time series are computed. If u has the dimension [20 x 3601], that means that time series made of 3601 time steps are generated in twenty positions. These positions are provided by the user. You can lso find the coordinates of the points in the variable "nodes".
I am relatively not so acquainted with Matlab,so my question is, what do the nodes represent?
the u velocity is 20x3601. 3601 are the samples according time..But the 20 is not really clear to me..
It is true that the turbulence length scales (Lux, Lvx and Lwx) change with the height. In the new version of windSim2, the length scales are estimated using the known behaviour of the wind spectrum at high frequencies. This leads to a relationship between the height and the length scales. This means also that the range of applicability of the "modified" von Karman model is limited to the atmospheric surface layer (ASL), i.e. the first 10% of the atmospheric boundary layer (ABL).
Regarding the standard deviations (stdU, stdV and stdW), if one assumes that they are proportional to the friction velocity u_star, then they are constant with the height as u_star is constant in the ASL if one uses the Monin-Obukhov similarity theory. This is the assumption I have done in windSIm2, although it is not necessarily verified in full-scale.
To my knowledge, there is no consensus on the height-dependency of the standard deviations (stdU, stdV and stdW) and the turbulence length scales (Lux, Lvx and Lwx) in the ABL.
I have included an example (Example4.m), where the Kaimal model is used. The latter model does not require Lux, Lvx and Lwx or stdU, stdV and stdW as it relies on the friction velocity, the measurement height and the mean wind velocity only.
@Qazi Shahzad Ali: Thank for the feedback. I am not sure I have understood all the parameters you wish to extract, but you can get several of them (delta_t, maxtime, Ntime) without simulating the data, simply based on the knowledge of the sampling frequency and by fixing either the duration or the number of time step.
I wand to take wind simulation data on following (NA,NR, radpol, thetapol, upol, delta_t, maxtime, Ntime) parameters to run my unsteady matlab code in Polar Grid (for 5 meter Dia Turbine). How can I acquire that data by using your code?
Need guidance. Thanks
Hi, Mr. Etienne.
Thanks for the great code.
I have a question. Shouldn't the turbulence length scales (Lux, Lvx and Lwx) and standard desviations (stdU, stdV and stdW) values change with each height used in the geometry file?
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.
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.
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
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!
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).
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!
My best answer is an additional html example file I have joined to the submission. I hope it will help you.
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?
@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 ?
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.
Find the treasures in MATLAB Central and discover how the community can help you!Start Hunting!