Visualize full-field Turbsim data and anemometer equivalent wind speed

Hi all,

Thank you very much for developing and maintaining the available software tools. In relation to my thesis I’m doing some simulations using the 5MW reference turbine, where I’m doing some condition monitoring. I have a few questions regarding Turbsim and InflowWind, since the various wind phenomena are important to my work.

When simulating the reference turbine, the file named “NRELOffshrBsline5MW_InflowWind_12mps” contains the outputs Wind1VelX", “Wind1VelY” and “Wind1VelZ” by default. I’m unable to find anything in the documentation regarding what these velocities are equivalent to (Where they a physically measured). Are these computed on the nacelle, where Wind1Velx is equal to what an anemometer would measure? If not, how do obtain something equivalent? Edit: On further thought, the wind speed measured by an anemometer could be intrepreted as sqrt(Wind1Velx^2 + Wind1Vely^2 +Wind1Velz^2 ).

I’m able to successfully generate different full-field wind data sets and run various simulations using these. Due to my inexperience with wind turbines and wind in general, I’m trying to figure out a best way to visualize the generated data. Also, since it seems that Turbsim uses various random events, see the figure below, which drastically changes the behaviour of the generated power.

I’m hoping that the wind data can be visualized as some sort of 3D plot, any suggestions on how to achieve a visualization of the full-field wind data that makes sense? And is there a way to receive/output a summary of which randoms events occur at the different time steps?

Any sort of help is appreciated. I’m sorry if the answers to my questions are very obvious, please forgive me if I have overlooked something simple.

Best regards,

Thomas Enevoldsen

Dear Thomas,

Regarding InflowWind outputs Wind1VelX, Wind1VelY, and Wind1VelZ, these are the ambient wind speeds (i.e. without the influence of the wind turbine induction / wake) at the 1st output point defined in the Inflow primary input file, i.e. the location defined in the WindVxiList, WindVyiList, and WindVziList lists. See the InflowWind documentation for more information:

TurbSim uses random phases of the various frequency components (the random phases are established via the seed specified in the TurbSim input file), which determine (randomly) where gusts and other events occur in the turbulent time series. You can visualize the turbulence field by reading the wind data in the TurbSim-generated file into MATLAB via the MATLAB script named readfile_BTS provided in the MATLAB Toolbox:

Best regards,

Dear Jason,

Thanks for the reply.

I have completely overlooked the locations defined by the WindVxiList, WindVyiList, and WindVziLists in the “NRELOffshrBsline5MW_InflowWind_12mps.dat” file. Sorry about that.
Also thank you for the clarification regarding the random events which occur when generating the wind data using Turbsim.

Regarding the readfile_BTS script, I encountered a similar script which was included with the installation of Turbsim. Though, I’m not quite sure what the best approach is to visualizing resulting data, which is in the form of an 4-D array. Ideally I would like some sort of 3-D plot, where it is possible to visually see at what times the gusts and other events appears.

I really appreciate your response, once again thank you.

Best regards,

Thomas Enevoldsen

Dear Thomas,

I asked Bonnie Jonkman of Envision Energy about this, and she created a simple MATLAB function that creates separate 3D color plots (time, Y, Z) for each velocity component (u,v, w; the u-component is likely most influential to the turbine response). She has shared the MATLAB function here: … d_figure.m.

I hope that helps.

Best regards,

Hi Jason,

Thank you and Bonnie for this, this will definitely be helpful. I really appreciate the help.

Best regards,

Thomas Enevoldsen

Hi Thomas,

Did you manage to get Bonnies linked function to work?

I am having difficulty in getting the function to read from my .bts file.



My script assumes you have installed the OpenFAST matlab-toolbox.

Dear Bonnie,

I have obtained the wind file in (.bts) format from Turbsim. How I can plot it in 3D and/or 4D format? like the ones in the manual (Fig.1). I tried using readfile_BTS MATLAB script but there is an error.
I would be appreciated to help with this matter.

Kind regards,

Dear Babak,

What error are you getting?

Best regards,

Dear Jason,

Thanks for your reply.
The error is ‘‘Not enough input arguments’’.
In the Matlab script, I have changed the FileName to the wind file from turbsim in bts format.


Dear Babak,

Can you show how you are calling the script from MATLAB? Did you change the source code in some way; if so, how?

Best regards,

Dear Jason,

This is the ‘‘readfile_BTS.m’’ script:

[code]function [velocity, twrVelocity, y, z, zTwr, nz, ny, dz, dy, dt, zHub, z1,mffws] = readfile_BTS(FileName,fileFmt)
%[velocity, twrVelocity, y, z, zTwr, nz, ny, dz, dy, dt, zHub, z1,mffws] = readfile_BTS(FileName,fileFmt)
% Author: Bonnie Jonkman, National Renewable Energy Laboratory
% Input:
% FileName - string: contains file name (.bts extension) to open
% fileFmt - string: optional, contains format of grid points.
% Output:
% velocity - 4-D array: time, velocity component (1=U, 2=V, 3=W), iy, iz
% twrVelocity - 3-D array: time, velocity component, iz
% y - 1-D array: horizontal locations y(iy)
% z - 1-D array: vertical locations z(iz)
% zTwr - 1-D array: vertical locations of tower points zTwr(iz)
% nz, ny - scalars: number of points in the vertical and horizontal
% directions of the grid
% dz, dy, dt - scalars: distance between two points in the vertical
% [m], horizontal [m], and time [s] dimensions
% zHub - scalar: hub height [m]
% z1 - scalar: vertical location of bottom of grid [m above ground level]
% mffws - scalar: mean hub-height wind speed

if ( nargin < 2 )
fileFmt = ‘int16’;
nffc = 3;

fid = fopen( FileName );

if fid > 0
% get the header information

tmp   = fread( fid, 1, 'int16');        % TurbSim format identifier (should = 7 or 8 if periodic), INT(2)

nz    = fread( fid, 1, 'int32');        % the number of grid points vertically, INT(4)
ny    = fread( fid, 1, 'int32');        % the number of grid points laterally, INT(4)
ntwr  = fread( fid, 1, 'int32');        % the number of tower points, INT(4)
nt    = fread( fid, 1, 'int32');        % the number of time steps, INT(4)

dz    = fread( fid, 1, 'float32');      % grid spacing in vertical direction, REAL(4), in m
dy    = fread( fid, 1, 'float32');      % grid spacing in lateral direction, REAL(4), in m
dt    = fread( fid, 1, 'float32');      % grid spacing in delta time, REAL(4), in m/s
mffws = fread( fid, 1, 'float32');      % the mean wind speed at hub height, REAL(4), in m/s
zHub  = fread( fid, 1, 'float32');      % height of the hub, REAL(4), in m
z1    = fread( fid, 1, 'float32');      % height of the bottom of the grid, REAL(4), in m

Vslope(1)  = fread( fid, 1, 'float32'); % the U-component slope for scaling, REAL(4)
Voffset(1) = fread( fid, 1, 'float32'); % the U-component offset for scaling, REAL(4)
Vslope(2)  = fread( fid, 1, 'float32'); % the V-component slope for scaling, REAL(4)
Voffset(2) = fread( fid, 1, 'float32'); % the V-component offset for scaling, REAL(4)
Vslope(3)  = fread( fid, 1, 'float32'); % the W-component slope for scaling, REAL(4)
Voffset(3) = fread( fid, 1, 'float32'); % the W-component offset for scaling, REAL(4)

    % Read the description string: "Generated by TurbSim (vx.xx, dd-mmm-yyyy) on dd-mmm-yyyy at hh:mm:ss."

nchar    = fread( fid, 1, 'int32');     % the number of characters in the description string, max 200, INT(4)
asciiINT = fread( fid, nchar, 'int8' ); % the ASCII integer representation of the character string
asciiSTR = char( asciiINT' );

disp( ['Reading from the file ' FileName ' with heading: ' ] );
disp( ['   "' asciiSTR '".' ] ) ;

% get the grid information

nPts        = ny*nz;
nv          = nffc*nPts;               % the size of one time step
nvTwr       = nffc*ntwr;

% velocity = zeros(nt,nffc,nPts);
velocity = zeros(nt,nffc,ny,nz);
twrVelocity = zeros(nt,nffc,ntwr);

if strcmpi(fileFmt,'float32')
    Voffset = 0.0*Voffset;
    Vslope  = ones(size(Vslope));

for it = 1:nt
    %get the grid points

   [v, cnt] = fread( fid, nv, fileFmt ); % read the velocity components for one time step
   if ( cnt < nv ) 
       disp([ it nt ny nz nffc nv cnt])
       error(['Could not read entire file: at grid record ' num2str( (it-1)*(nv+nvTwr)+cnt ) ' of ' num2str(nt*(nv+nvTwr))]);

    ip = 1;
    for iz = 1:nz
        for iy = 1:ny
            for k=1:nffc
                velocity(it,k,iy,iz) = ( v(ip) - Voffset(k))/Vslope(k) ;
                ip = ip + 1;
            end %k
        end %iy
    end % iz

    %get the tower points
    if nvTwr > 0
        [v, cnt] = fread( fid, nvTwr, fileFmt ); % read the velocity components for the tower
        if ( cnt < nvTwr ) 
            error(['Could not read entire file: at tower record ' num2str( (it-1)*(nv+nvTwr)+nv+cnt ) ' of ' num2str(nt*(nv+nvTwr))]);

        for k=1:nffc      % scale the data
            twrVelocity(it,k,:) = (v(k:3:nvTwr) - Voffset(k))/Vslope(k); 

end %it


error(['Could not open the wind file: ’ FileName]) ;

y = [0:ny-1]dy - dy(ny-1)/2;
z = [0:nz-1]*dz + z1;
zTwr = z1 - [0:ntwr-1]*dz;

it is in MATLAB script code so I just call it from MATALB and as as far as I know,I only need to change the FileName to my wind file in .bts format, right?

Dear Babak,

Correct. Here is a sample:

[velocity, twrVelocity, y, z, zTwr, nz, ny, dz, dy, dt, zHub, z1,mffws] = readfile_BTS('90m_12mps_twr.bts');

Is that how you are calling readfile_BTS.m?

Best regards,

Dear Jason,

Thanks for your reply.
The problem was about entering file name. It needs to be entered as 'FileNamenot just FileName.
Now I have a 4D velocity. Could you please tell me how can I plot it to obtain figures like in my previous post?

Kind regards,

Dear Babak,

That is what Bonnie Jonkman’s script named FFwind_figure.m is for. See the link to FFwind_figure.m in my forum post dated Mar 20, 2019 above.

Best regards,

Dear Jason,

Thanks for your help. Now, I have obtained the following graphs for each of the wind speed components, as below ( shown only for u component)

Please met me know how I can have plots as below?
I appreciate for your helps.

Dear Babak,

It looks like you already generated the image on the left (although your TurbSim grid is much coarser than the image from the TurbSim documentation). The image on the right is the same, but without any turbulence.

Best regards,

Dear Jason,

Thanks for your reply.
Do you mean NumGrid_Z and NumGrid_Y ? or other parameters needs to modified?
Here I have attached my Tubsim input file.
Thanks for your time.

[code]TurbSim Input File. Valid for TurbSim v1.06.00, 21-Sep-2012

---------Runtime Options-----------------------------------
2318573 RandSeed1 - First random seed (-2147483648 to 2147483647)
RanLux RandSeed2 - Second random seed (-2147483648 to 2147483647) for intrinsic pRNG, or an alternative pRNG: “RanLux” or “RNSNLW”
False WrBHHTP - Output hub-height turbulence parameters in binary form? (Generates RootName.bin)
False WrFHHTP - Output hub-height turbulence parameters in formatted form? (Generates RootName.dat)
False WrADHH - Output hub-height time-series data in AeroDyn form? (Generates RootName.hh)
True WrADFF - Output full-field time-series data in TurbSim/AeroDyn form? (Generates Rootname.bts)
False WrBLFF - Output full-field time-series data in BLADED/AeroDyn form? (Generates RootName.wnd)
False WrADTWR - Output tower time-series data? (Generates RootName.twr)
False WrFMTFF - Output full-field time-series data in formatted (readable) form? (Generates RootName.u, RootName.v, RootName.w)
False WrACT - Output coherent turbulence time steps in AeroDyn form? (Generates RootName.cts)
True Clockwise - Clockwise rotation looking downwind? (used only for full-field binary files - not necessary for AeroDyn)
0 ScaleIEC - Scale IEC turbulence models to exact target standard deviation? [0=no additional scaling; 1=use hub scale uniformly; 2=use individual scales]

--------Turbine/Model Specifications-----------------------
5 NumGrid_Z - Vertical grid-point matrix dimension
5 NumGrid_Y - Horizontal grid-point matrix dimension
0.05 TimeStep - Time step [seconds]
600 AnalysisTime - Length of analysis time series [seconds] (program will add time if necessary: AnalysisTime = MAX(AnalysisTime, UsableTime+GridWidth/MeanHHWS) )
100 UsableTime - Usable length of output time series [seconds] (program will add GridWidth/MeanHHWS seconds)
34.6189 HubHt - Hub height [m] (should be > 0.5GridHeight)
8 GridHeight - Grid height [m]
8 GridWidth - Grid width [m] (should be >= 2
0 VFlowAng - Vertical mean flow (uptilt) angle [degrees]
0 HFlowAng - Horizontal mean flow (skew) angle [degrees]

--------Meteorological Boundary Conditions-------------------
“IECKAI” TurbModel - Turbulence model (“IECKAI”=Kaimal, “IECVKM”=von Karman, “GP_LLJ”, “NWTCUP”, “SMOOTH”, “WF_UPW”, “WF_07D”, “WF_14D”, “TIDAL”, or “NONE”)
“1-ED3” IECstandard - Number of IEC 61400-x standard (x=1,2, or 3 with optional 61400-1 edition number (i.e. “1-Ed2”) )
“C” IECturbc - IEC turbulence characteristic (“A”, “B”, “C” or the turbulence intensity in percent) (“KHTEST” option with NWTCUP model, not used for other models)
“NTM” IEC_WindType - IEC turbulence type (“NTM”=normal, “xETM”=extreme turbulence, “xEWM1”=extreme 1-year wind, “xEWM50”=extreme 50-year wind, where x=wind turbine class 1, 2, or 3)
default ETMc - IEC Extreme Turbulence Model “c” parameter [m/s]
default WindProfileType - Wind profile type (“JET”;“LOG”=logarithmic;“PL”=power law;“H2L”=Log law for TIDAL spectral model;“IEC”=PL on rotor disk, LOG elsewhere; or “default”)
34.6189 RefHt - Height of the reference wind speed [m]
7 URef - Mean (total) wind speed at the reference height [m/s] (or “default” for JET wind profile)
default ZJetMax - Jet height [m] (used only for JET wind profile, valid 70-490 m)
default PLExp - Power law exponent [-] (or “default”)
default Z0 - Surface roughness length [m] (or “default”)

--------Non-IEC Meteorological Boundary Conditions------------
default Latitude - Site latitude [degrees] (or “default”)
0.05 RICH_NO - Gradient Richardson number
default UStar - Friction or shear velocity [m/s] (or “default”)
default ZI - Mixing layer depth [m] (or “default”)
default PC_UW - Hub mean u’w’ Reynolds stress (or “default”)
default PC_UV - Hub mean u’v’ Reynolds stress (or “default”)
default PC_VW - Hub mean v’w’ Reynolds stress (or “default”)
default IncDec1 - u-component coherence parameters (e.g. “10.0 0.3e-3” in quotes) (or “default”)
default IncDec2 - v-component coherence parameters (e.g. “10.0 0.3e-3” in quotes) (or “default”)
default IncDec3 - w-component coherence parameters (e.g. “10.0 0.3e-3” in quotes) (or “default”)
default CohExp - Coherence exponent (or “default”)

--------Coherent Turbulence Scaling Parameters-------------------
“C:\users\bmehdiza\Desktop\coh-events\EventData” CTEventPath - Name of the path where event data files are located
“Random” CTEventFile - Type of event files (“LES”, “DNS”, or “RANDOM”)
true Randomize - Randomize the disturbance scale and locations? (true/false)
1.0 DistScl - Disturbance scale (ratio of wave height to rotor disk). (Ignored when Randomize = true.)
0.5 CTLy - Fractional location of tower centerline from right (looking downwind) to left side of the dataset. (Ignored when Randomize = true.)
0.5 CTLz - Fractional location of hub height from the bottom of the dataset. (Ignored when Randomize = true.)
30.0 CTStartTime - Minimum start time for coherent structures in RootName.cts [seconds]

NOTE: Do not add or remove any lines in this file!

Dear Babak,

Yes, NumbGrid_Y and NumGrid_Z. Normally you’d want to set these so that spatial increment (dY and dZ):

dY = GridWidth/(NumGrid_Y-1) dZ = GridHeight/(NumGrid_Z-1)
are roughly equal to maximum chordlength of the rotor blade.

Bets regards,

Dear Babak,
Can you tell me how to plot the pictures as the one that you have ploted?

Thank you
Best ragards