%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Wind Generator for Aero-elastic simulation purposes %
% Author: Francesco Perrone %
% E-mail: francesco.perrone@forwind.de %
% Date: 18.03.2013 %
% Release: v.1.0 %
% %
% Inputs: %
% 1) Hub Height [m] %
% 2) Grid Width [m] %
% 3) Grid Height[m] %
% 4) Number of points along Y %
% 5) Number of points along Z %
% 6) Number of FFT points (2^12 recommended) %
% 7) Analysis Time (generally 600s) %
% 8) Hub Average Wind Speed [m/s] %
% 9) Coherence Model (TurbSim or Bladed) %
% %
% Outputs: %
% 1) U: Turbulent time series for v-component [m/s] %
% 2) V: Turbulent time series for v-component [m/s] %
% 3) W: Turbulent time series for v-component [m/s] %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clc,clear all,close all
%% INITIALIZATION
tic
HubHt = 90; % Hub Height
GridWidth = 150; % Grid length along Y axis
GridHeight = 150; % Grid length along Z axis
Ny = 35;
Nz = 33;
N = Ny*Nz;
Nfft = 2^12; % It has to be a power of 2 number (Generally 2^12 = 4096 or 2^13 = 8192 are sufficient)
time = 620; % Simulation tim
dt = time/Nfft; % time resolution
Uhub = 8; % Average wind speed at hub height
Coherence = 'Bladed'; % It can be 'TurbSim' or 'Bladed'
%% TURBULENCE PARAMTERS
%set delta1 according to guidelines (chap.6)
if HubHt < 60
lambda1 = 0.7*HubHt;
else
lambda1 = 42;
end
L = 0.8*lambda1;%IEC, (B.12)
%IEC, Table 1, p.22
%IEC, sect. 6.3.1.3
b=5.6;
sigma1 = 1.0;
sigma2 = 0.8;
sigma3 = 0.5;
%derived constants
l=0.8*lambda1; %IEC, (B.12)
xLu = 8.1*l; % Integral length scale along X axis
xLv = 2.7*l; % Integral length scale along Y axis
xLw = 0.66*l; % Integral length scale along Z axis
%% FREQUENCY
% Frequency defition
t = 0:dt:(time-dt);
df = 1/(time);
fn = Nfft/time;
df_5 = df/2;
f = (0:length(t)/2-1)*df;
f_fit = (-length(t)/2:length(t)/2-1)*df;
fit_u = normpdf(f_fit,0,1);
%% GRID DEFINITION
dy = GridWidth/(Ny-1); % Spacing along Y axis
dz = GridHeight/(Nz-1); % Spacing along Z axis
if isequal(mod(Ny,2),0)
iky = [(-Ny/2:-1) (1:Ny/2)];
else
iky = -floor(Ny/2):ceil(Ny/2-1);
end
if isequal(mod(Nz,2),0)
ikz = [(-Nz/2:-1) (1:Nz/2)];
else
ikz = -floor(Nz/2):ceil(Nz/2-1);
end
% define Y and Z coordinates for each grid-point
[Y Z] = ndgrid(iky*dy,(ikz*dz + HubHt));
% define distance matrix
dist = pdist2([Y(:) Z(:)],[Y(:) Z(:)]);
%% COHERENCE MATRIX
f_u = (f)*(xLu/Uhub);
f_v = (f)*(xLv/Uhub);
f_w = (f)*(xLw/Uhub);
% Define normalized power spectra
psd = zeros(length(f),3);
psd(:,1) = (4.*(sigma1.^2).*(xLu/Uhub))./((1 + 6.*f_u).^(5/3));
psd(:,2) = (4.*(sigma2.^2).*(xLv/Uhub))./((1 + 6.*f_v).^(5/3));
psd(:,3) = (4.*(sigma3.^2).*(xLw/Uhub))./((1 + 6.*f_w).^(5/3));
% Coherence constant parameters
a = 12;
Lc = 5.67*min(60,HubHt);
% Generation of Coherence matrix and FFT terms
H_u = zeros(N,length(f));
H_v = zeros(size(H_u));
H_w = zeros(size(H_u));
nn_u = complex((normrnd(0,1.0,size(H_u))),(normrnd(0,1.0,size(H_u))));
nn_v = complex((normrnd(0,1.0,size(H_u))),(normrnd(0,1.0,size(H_u))));
nn_w = complex((normrnd(0,1.0,size(H_u))),(normrnd(0,1.0,size(H_u))));
% TurbSim and Bladed manage the Coherence topic in diametrically opposite
% way, especially when concerning correlation for v- and w- component.
% TurbSim implies correlation only for the u- component of the turbulent
% field, which is in accordance with what stated within the IEC 61400-1:
% indeed, the above standard only mention Coherence for u-component.
% Bladed, instead, introduce a weak correlation for v- and w- components as
% well. This correlation is taken into account by using the standard
% Davenport relationship.
% It must be noted that choosing Bladed coherence is computationally more
% expensive, because of requiring the Cholesky factorization of each of the
% turbulent components.
switch Coherence
case 'TurbSim'
for ii = 1:numel(f)
Coh_u = exp(-a.*dist.*sqrt((f(ii)/Uhub).^2 + (0.12/Lc).^2)).*(df.*psd(ii,1)); % Coherence matrix
Coh_v = eye(N).*sqrt(df.*psd(ii,2));
Coh_w = eye(N).*sqrt(df.*psd(ii,3));
HH_u = chol(Coh_u,'lower'); % Cholesky factorization of the Coherence matrix for u-component
H_u(:,ii) = HH_u*nn_u(:,ii); % Computation of the FFT terms in x direction
H_v(:,ii) = Coh_v*nn_v(:,ii); % Computation of the FFT terms in y direction
H_w(:,ii) = Coh_w*nn_w(:,ii); % Computation of the FFT terms in z direction
end
case 'Bladed'
for ii = 1:numel(f)
if ~isequal(f(ii),0)
Coh_u = exp(-a.*dist.*sqrt((f(ii)/Uhub).^2 + (0.12/Lc).^2)).*(df.*psd(ii,1)); % Coherence matrix
Coh_v = exp(-a.*dist.*(f(ii)/Uhub)).*(df.*psd(ii,2));
Coh_w = exp(-a.*dist.*(f(ii)/Uhub)).*(df.*psd(ii,3));
HH_u = chol(Coh_u,'lower'); % Cholesky factorization of the Coherence matrix for u-component
HH_v = chol(Coh_v,'lower'); % Cholesky factorization of the Coherence matrix for v-component
HH_w = chol(Coh_w,'lower'); % Cholesky factorization of the Coherence matrix for w-component
H_u(:,ii) = HH_u*nn_u(:,ii); % Computation of the FFT terms in x direction
H_v(:,ii) = HH_v*nn_v(:,ii); % Computation of the FFT terms in y direction
H_w(:,ii) = HH_w*nn_w(:,ii); % Computation of the FFT terms in z direction
end
end
end
H_u(:,1) = complex(0,0);
H_v(:,1) = complex(0,0);
H_w(:,1) = complex(0,0);
%% TIME SERIES
U = zeros(N,size(H_u,2)*2);
V = zeros(N,size(H_u,2)*2);
W = zeros(N,size(H_u,2)*2);
for ii = 1:N
tmp_u = (fft((H_u(ii,:)))); % FFT
tmp_v = (fft((H_v(ii,:))));
tmp_w = (fft((H_w(ii,:))));
% Up to this point only the one-sided spectrum has been considered to
% lighten the computational effort, since the double-sided spectrum is
% an even function about the frequency axis.
% Therefore, we can now merge real and imaginary part of the fft to get
% the ultimate time series.
U(ii,:) = [real(tmp_u) imag(tmp_u)];
V(ii,:) = [real(tmp_v) imag(tmp_v)];
W(ii,:) = [real(tmp_w) imag(tmp_w)];
end
sec2timestr(toc)
%% POST PRO
% Scale velocities
% The time series are to be scaled to the right statistics and Gaussian
% shape for a load assumption purpose. Indeed, each time series, must be
% well representive of the given turbulent intensity, so that any present
% smoothing effect has to be carefully avoided.
for ii = 1:N
r_u = tiedrank( U(ii,:) );
p_u = r_u / ( length(r_u) + 1 ); %# +1 to avoid Inf for the max point
r_v = tiedrank( V(ii,:) );
p_v = r_v / ( length(r_v) + 1 ); %# +1 to avoid Inf for the max point
r_w = tiedrank( W(ii,:) );
p_w = r_w / ( length(r_w) + 1 ); %# +1 to avoid Inf for the max point
U(ii,:) = fftshift(norminv( p_u, 0, 1 ));
V(ii,:) = fftshift(norminv( p_v, 0, .8 ));
W(ii,:) = fftshift(norminv( p_w, 0, .5 ));
end