3D Turbulent Wind Field Simulation by means of Sandia Method for wind energy applications

by

 

27 Mar 2012 (Updated )

3D turbulent wind field generator with Kaimal model

wind_gen.m
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 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

Contact us