windsim.m

风力发电系统中的风速仿真模型

function [ux,uy,tv] = windsim(t,fsamp,fcut,u,ak,fik)
%|||||||||||||||||||||||||||
% Output :
% ux = vector of windspeed sampled at a frequency fsamp over a total time t [m=s]
% uy = the timederivates of the windspeed in ux [m=s2]
% tv = the time vector according to ux and uy [s]
% Input :
% t = total time [s]
% fsamp = sampling frequency [Hz]
% fcut = the (upper) frequency at which to cut the spectrum [Hz]
% u = mean windspeed over the time t [m=s]
% ak = vector (length > t*fcut) of Rayleigh distributed random numbers
% fik = vector (length > t*fcut) of random numbers uniformly dist. in [0; 2pi]
%|||||||||||||||||||||||||||-
len = round(t*fsamp);% the length of the output vectors
lenc = round(t*fcut); % the number of frequencies up to fcut at which to calculate the spectrum
ak = ak(1:lenc);
fik = fik(1:lenc); % cutting the vectors of random numbers to appropriate length
kdel = (1:lenc)/t; % vector of frequencies up to fcut at which to calculate the spectrum
z = 30; % the height of the point where the wind is simulated
z0 = 0.3; % a ground roughness factor
s = 16.8/(log(z/z0)^2)*z*u./((1+33*z/u*kdel).^(5/3));
% the spectrum calculated at the frequencies of kdel
ex = sqrt(s/t).*ak.*exp(-i*fik);
% the expression to put into the FFT to calculate ux
ey = 2*pi*kdel.*sqrt(s/t).*ak.*exp(-i*(fik+2*pi/2));
% the expression to put into the FFT to calculate uy
n = 0; % n contains the indices of the output vector
len - 1;
cor1 = exp(-i*2*pi*n/len);
cor2 = kdel*t*pi*(1/len-1/(t*fsamp));
% correlation terms to make the FFT useble for our purposes
ux = u+real(cor1.*(fft(ex,len)+n.*fft(i*2*cor2.*ex,len)-n.^2.*fft(2*cor2.^2.*ex,len)));
% calculation of the output ux
uy = real(cor1.*(fft(ey,len)+n.*fft(i*2*cor2.*ey,len)-n.^2.*fft(2*cor2.^2.*ey,len)));
% calculation of the output uy
tv = n/fsamp; % calculation of the output tv