Visualize full-field Turbsim data and anemometer equivalent wind speed

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.

Regards,

Joseph,

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,
Babak
winds.JPG

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.

regards
Babak

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’;
end
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));
end


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])
       fclose(fid);
       error(['Could not read entire file: at grid record ' num2str( (it-1)*(nv+nvTwr)+cnt ) ' of ' num2str(nt*(nv+nvTwr))]);
   end

    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))]);
        end

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

end %it

fclose(fid);

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

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

return;[/code]
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,
Babak

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

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
(RotorRadius+ShaftLength))
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!
==================================================[/code]

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

Dear Long,

The MATLAB scripts used to read the TurbSim output and plot these figures are included in the forum posts in this topic above.

Best regards,

Dear mr. Jonkman

Hello and sorry to bother you again. I am now trying to run a longer simulation (300sec) and I need to produce my own wind field with TurbSim. I downloaded the library from Github “GitHub - old-NWTC/TurbSim: A stochastic, full-field, turbulence simulator”. And I try to run TurbSim through windows command prompt as I do with FAST. It is bothering me that it doesn’t have the ending “.exe” like the FAST does but the ending “f90” . When I try to run it I get the message:“‘TurbSim’ is not recognized as an internal or external command,operable program or batch file.” Is there any case that the turbsim file is broken? Or am I doing something wrong??

Thank you a lot in advance

Dear Manolis,

The *.f90 file is the source code; you must compile this version of TurbSim before running a simulation. The compiled version should have a *.exe extension on Windows.

You can find an archive of TurbSim v2.00.07a-bjj on my Google drive that contains the TurbSim source code, documentation, sample models, CertTest, and precompiled Windows executable here: drive.google.com/file/d/1eCvgoO … sp=sharing.

Best regards,

Dear mr. Jonkman

Thank you very much for your help. I managed to produce my own wind field.

Have a very good day