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