Running FAST based on CFD turbulence fields

Dear all,

I am attempting to do load calculations on wind turbines based on turbulence fields that I generate using large-eddy simulation (LES). Right now I am able to do my LES simulations and to collect two-dimensional turbulent flow fields (basically I am saving the velocity field in the y-z plane with a certain time interval). I would like to use FAST/AeroDyn for simulating the dynamic response of a wind turbine (probably the NREL 5MW standard turbine).
Now my question: How should my data look in order to use it in FAST/AeroDyn? Is there a certain file format that I should save my data in?

Any ideas will be greatly appreciated!

Best regards,
Bastian.

Dear Bastian,

If I read this correctly, you want to produce wind inflow (turbulence) files that can be read in a FAST/AeroDyn simulation. in FAST, these files are actually read by a module called InflowWind, which accepts several different formats (or allows you to create your own format in a UserWind routine, assuming you code the subroutine to read the format you created.)

Most likely, you would create binary full-field turbulece files, such as the ones TurbSim generates (in a y-z plane at constant time intervals; with a .bts extension). The format is described in appendix D of the TurbSim User’s Guide. I also have pulled out the TurbSim code (Fortran) that creates the wind files, and I made some modifications so that it should work outside TurbSim; you can find it here: http://wind.nrel.gov/public/bjonkman/WriteTurbSimBinary.f90 You will need to use the NWTC Library to compile/use it. I’ve tried to make all the variables you’ll need to define as inputs. Note that TurbSim’s V array is only 3 dimensions: dimension 1 is time, dimension 2 is grid point, and dimension 3 is velocity component (U,V,W – see figure 3 in the TurbSim User’s Guide). The top of the routine indicates the numbering convention of the grid points (2nd component).

TurbSim also produces a Bladed-style FF binary file (.wnd, with a corresponding .sum file). This format is described in appendix E in the user’s guide; you can also read this post for help in generating them: http://forums.nrel.gov/t/binary-wnd-specification/180/1

You can use the InflowWind module’s sample driver/test programs (one in Matlab, the other in Fortran) to verify that the files are generated correctly.

Dear Bonnie,

Thank you very much for your extensive answer. You read correctly, I indeed want to produce wind inflow files that can be read by a FAST/Aerodyn simulation. I will try to incorporate the TurbSim subroutine into my CFD code directly, which should be the smoothest way of producing the files, I assume.
I will keep you updated on my progress.

Thanks again.

/Bastian

Hi All,

Just as an FYI, I was struggling to do this a while ago and wrote this Matlab script. It’s very basic and was used to generate wave velocities (water waves) in time and space (hence the ‘Airy_’ functions). If it can help anyone else doing something unorthodox then feel free to use it, although I don’t make any guarantees about it’s correctness!

Pascal

Script used to create the four-dimensional input ’wind’ velocity files for FAST. The script calculates
wave particle velocity using linear wave theory in 3 dimensions in both space and time.
%Script to create a 4−D array for manipulation into a .wnd binary velocity
%full−field file to be used in FAST. 18/04/12. Pascal Galloway.
clear all
close all
FileName = '05Hz_150mm_09zHub';
% FileName = 'nowaves';
% NOTE − if no waves are desired, but full−field file still required,
% change u, v and w to zero in loops
total_time = 60; %run time (s)
dx = 0.05; %delta x in m
dy = 0.1; %lateral point spacing (m)
dz = 0.1; %vertical point spacing (m)
ny = 20; %number of lateral points
nz = 20; %number of vertical points
zHub = 1.21; %reference hubheight (m) − 0.9 > 0.95, 0.67 > 1.21
MFFWS = 0.9; %mean full field current speed (m/s)
d = 1.880; %water depth (m)
H = 0.150; %wave height (m)
Ti = 2; %intrinsic wave period (s) − 1.33(0.75), 2(0.5) or 2.85(0.35)
L = Airy_Wavelength(d,Ti); %calculating the wavelength
dt = dx/MFFWS; %time difference (s).
y = linspace(−ny/2*dy,ny/2*dy,ny); %equally spaced vector
z = linspace(−nz/2*dz,nz/2*dz,nz)+zHub; %equally spaced vector
z1 = z(1); %vertical location of bottom of grid [m above ground level]
193
time = 10:dt:total_time+10; %time vector
zi = −z; %depth position (m) note that d=−z
nt = (length(time)−1)/2; %half number of time steps
%format according to TurbSim
velocity = zeros(length(time),3,length(ny),length(nz));
for it = 1:length(time)
for iz = 1:nz
for iy = 1:ny
[u,v,w] = Airy_Range_TurbSim(d,H,Ti,MFFWS,time(it),zi(iz),L);
velocity(it,1,iy,iz) = u; %axial velocity (m/s)
velocity(it,2,iy,iz) = v; %lateral velocity (m/s)
velocity(it,3,iy,iz) = w; %vertical velocity (m/s)
end
end
end
%% initialize variables
fileFmt = 'int16'; %file format 16−bit integer
ConvFact = 1.0; %results in meters and seconds (change conversion factor
% for other units)
str = {'HUB HEIGHT','CLOCKWISE','UBAR','TI(U','TI(V','TI(W'};
%MUST be in UPPER case
numVars = length(str);
%% −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
% WRITING THE HEADER TO THE BINARY FILE
% −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
fid_wnd = fopen( [ FileName '.wnd' ],'w');
if ( fid_wnd <= 0 )
error( 'Wind file could not be opened.' );
end
nffc = −99; fwrite( fid_wnd, nffc, 'int16' );
% number of components
fc = 4; fwrite( fid_wnd, fc, 'int16' );
% should be 4 to allow turbulence intensity to be stored in the header
nffc = 3; fwrite( fid_wnd, nffc, 'int32' );
lat = 0.0; fwrite( fid_wnd, lat, 'float32' );
% latitude (deg)
z0 = 0.001; fwrite( fid_wnd, z0, 'float32' );
% Roughness length (m)
fwrite( fid_wnd, zHub, 'float32' );
% Reference height (m) = Z(1) + GridHeight / 2.0
TI_U = 1; fwrite( fid_wnd, TI_U, 'float32' );
% Turbulence Intensity of u component (%)
TI_V = 1; fwrite( fid_wnd, TI_V, 'float32' );
% Turbulence Intensity of v component (%)
TI_W = 1; fwrite( fid_wnd, TI_W, 'float32' );
% Turbulence Intensity of w component (%)
fwrite( fid_wnd, dz, 'float32' ); % delta z in m
fwrite( fid_wnd, dy, 'float32' ); % delta y in m
fwrite( fid_wnd, dx, 'float32' ); % delta x in m
fwrite( fid_wnd, nt, 'int32' ); % half the number of time steps
fwrite( fid_wnd, MFFWS, 'float32'); % mean full−field wind speed
fwrite( fid_wnd, 0, 'float32' );
fwrite( fid_wnd, 0, 'float32' );
194
fwrite( fid_wnd, 0, 'float32' ); %zLu,yLu,xLu:unused variables (for BLADED)
fwrite( fid_wnd, 0, 'int32' );
fwrite( fid_wnd, 0, 'int32' ); %unused variables (for BLADED)
fwrite( fid_wnd, nz, 'int32' ); %number of points in vertical direction
fwrite( fid_wnd, ny, 'int32' ); %number of points in horizontal direction
fwrite( fid_wnd, 0, 'int32' );
fwrite( fid_wnd, 0, 'int32' );
fwrite( fid_wnd, 0, 'int32' );
fwrite( fid_wnd, 0, 'int32' );
fwrite( fid_wnd, 0, 'int32' );
fwrite( fid_wnd, 0, 'int32' );
%other length scales: unused variables (for BLADED)
%% −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
% WRITE THE GRID DATA TO THE BINARY FILE
% −−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
nt = max([nt*2,1]); %number of time steps
nv = nffc*ny*nz; %the size of one time step
SummVars = [zHub; 1; MFFWS; 1; 1; 1];
%data from .sum file − SummVars(2) if '−1' anti−clockwise rotation
Scale = 0.00001*SummVars(3)*SummVars(4:6);
Offset = [SummVars(3) 0 0];
if SummVars(2) > 0 %clockwise rotation
%flip the y direction
y_ix = ny:−1:1;
else
y_ix = 1:ny;
end
%changing the dimension of velocity so that it's 4−D instead of 3−D
velocity_norm = zeros(length(time),3,length(ny),length(nz));
for it = 1:nt
for iz = 1:nz
for iy = y_ix
for k = 1:nffc
velocity_norm(it,k,iy,iz) = velocity(it,k,iy,iz)...
/Scale(k) − Offset(k);
fwrite( fid_wnd, velocity_norm(it,k,iy,iz), fileFmt );
end %for k
end %iy
end % iz
end %it
fclose(fid_wnd);
% Linear Wave Theory accounting for Doppler shift of uniform current
% Outputs feed into 4Darray.m
% 18/04/12. Pascal Galloway.
% This function gives horizontal and vertical velocity in
% time and depth for a uniform current with small amplitude waves
function [u,v,w] = Airy_Range_TurbSim(d,H,Ti,U,ti,zi,L)
% Inputs:
% d = water depth (m)
% H = wave height (m)
% Ti = intrinsic wave period (s)
% U = uniform current speed (m/s)
% zi = depth position (m) note that d=−z
phi = 0; %phase angle i.e. phase shift, use for moving along x−axis
mu = 0; %angle of inclination between wave and current
195
k = (2*pi)/L; %wavenumber
Ta = (2*pi*Ti)/(2*pi + Ti*k*U); %observed wave period at rotor
sigma = ((2*pi)/Ta)*ti + phi; %wave phase
steepness = H/L;
%steepness of waves. when ratio reaches 1:7 waves break i.e. 0.14
if steepness > 0.14
display 'WAVE BREAKING!'
end
u = ((pi*H)/Ti)*cos(mu)*(cosh(k*(d+zi))/sinh(k*d))*cos(sigma);
v = 0; %zero lateral velocity in towing tank
w = ((pi*H)/Ti)*(sinh(k*(d+zi))/sinh(k*d))*sin(sigma);
return

Thanks, Pascal. We always appreciate users helping each other!