File Exchange

image thumbnail

Wind field simulation (the user-friendly version)

version 7.3 (950 KB) by E. Cheynet
Simulation of spatially correlated wind velocity time-histories


Updated 03 Mar 2019

View License

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 [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 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
- 4 examples files Example1.m, Example2.m, Example3.m and Example4.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.
- 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
- A faster version of the present submission has been used to simulate the turbulent wind load on a floating suspension bridge [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] 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.

Cite As

E. Cheynet (2019). Wind field simulation (the user-friendly version) (, MATLAB Central File Exchange. Retrieved .

Comments and Ratings (32)

E. Cheynet

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. [1]), 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.

[1] 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?

E. Cheynet

Hi Sebastian,
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: 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.

E. Cheynet

Hi Sebastian,
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?

E. Cheynet

Hi Frane,
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.

Hello there,
I am having trouble starting the example 1. Do I copy this code

clearvars;close all;clc;
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?

E. Cheynet

Hi Ifigenia,
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".

Hey everyone,

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..

Med El Ouni


Xia Qi

ni gi

E. Cheynet

@Alan Barile

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.

E. Cheynet

@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.

Excellent code..!!!
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

Alan Barile

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?



Pengyuan Qi


E. Cheynet

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

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

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



Lei Tong

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

E. Cheynet

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).


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

Hi djr,
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?


E. Cheynet

@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.




Added the project website

The example files have been re-written as live-script

Updated the title


Major update on windSim2.m: The integral length scale in the von Karman spectral model is now a function of the height. its value is defined using the asymptotic behaviour of the spectrum at high frequencies. + New Example file.

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.


- new Example, and updated description

updated picture

summary updated

new title

- new example included

- Updated example (need 90 sec to run instead of 2-3 sec ) but it is more complete

-description updated

-description updated

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

- The cross-wind component v has been added !
- Example updated

- update following @Paolo's feedback

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


-correct update now



- title of example file is changed

- New definition of the vector frequency to optimize the computation time
- introduction of the co-coherence function
- The co-coherence is used instead of the magnitude squared coherence to be more rigorous

-dexcription updated



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

-picture updated

- erratum at the end of Run_Me is corrected

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


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

- description updated

-description updated

MATLAB Release Compatibility
Created with R2018a
Compatible with any release
Platform Compatibility
Windows macOS Linux

Inspired: Wind field simulation (the fast version)