%风电塔项目—风荷载时程模拟-谐波叠加法-Davenport谱
clear;
clc;
%定义基本信息
nd=21; %模拟点数
Um=37.5; %基本风速,通常来讲是10m,对于风电塔结构基本风速是其风机处的风速
H0=65; %基本高度,对于风电塔结构是风机处的高度
dt=0.1; %模拟风速时程的时间步长
wmax=2*pi*10; %上限角频率
T=6000; %模拟总点数
N=ceil(T/(2*nd)); %频率采样点数
M=2*N; %逆傅里叶变换长度,取采样点数的2倍
nT=M*nd; %一个模拟点的时间点数
t=dt*(1:1:nT); %时间序列
dw=wmax/N; %角频率步长
%定义各类常数
r=0.16; %指数型风剖面,考虑地表粗糙度影响的无量纲幂指数,按中国规范取0.16(B)类,荷载规范中可找到
k=0.00215; %地面粗糙度系数k=0.00129(A)0.00215(B)0.00464(C)0.01291(D),此数据来源于文献
Cz=10; %竖直方向指数衰减系数
%定义模拟点竖向高度(不考虑横向相干性)
z(1)=5.880;z(2)=8.820;z(3)=11.760;z(4)=16.554;z(5)=19.494;
z(6)=22.434;z(7)=25.374;z(8)=28.314;z(9)=31.254;z(10)=34.194;
z(11)=37.109;z(12)=39.809;z(13)=42.514;z(14)=45.224;z(15)=47.939;
z(16)=50.659;z(17)=53.384;z(18)=56.114;z(19)=58.849;z(20)=61.569;
z(21)=63.209;
%求平均风速
U0=Um*(10/H0)^r; %标准高度为10m的平均风速
for i=1:nd
U(i)=U0*(z(i)/10)^r;
end %求指数型分布风剖面各模拟点风速
%定义风速谱(Davenport谱)
w=wmax/N:wmax/N:wmax;
n=w./(2*pi);
x=1200*n./U0;
Sv=4*k*U0^2*x.^2./n./(1+x.^2).^(4/3); %Davenport谱
%谐波叠加法运算,利用FFT优化运算速度
V=zeros(nd,nT);
D=zeros(nd,nd,N);
for j=1:nd
rand('state',0); %模拟随机数序列
thet=2*pi*rand(j,N); %生成随机相位
for i=1:N %生成频率序列
w1(i)=(i-1)*dw+j/nd*dw;
end
n1=w1./(2*pi);
x1=1200*n1./U0; %谱的参数
Sv1=4*k*U0^2*x1.^2./n1./(1+x1.^2).^(4/3); %Davenport谱
for j1=1:nd
for j2=1:nd
for i=1:N
Coh(j1,j2,i)=(exp((-n1(i)*sqrt(Cz^2*(abs(z(j1)-z(j2)))^2)*2)/(U(j1)+U(j2))));
S(j1,j2,i)=Sv1(i)*Coh(j1,j2,i);
end
end
end %生成谱矩阵
for i=1:1:N
H(:,:,i)=chol(S(:,:,i));
H(:,:,i)=H(:,:,i)';
end %Cholsky分解
D(:,j,:)=H(:,j,:);
ic=sqrt(-1);
B1=sqrt(2*dw).*D(j,:,:);
for ii=1:j
for jj=1:N
B2(ii,jj)=B1(1,ii,jj);
end
end
B2=B2.*exp(ic.*thet);
for jj=1:j
G(jj,1:M)=fft(B2(jj,:),M);
for jjj=2:nd
G(jj,((jjj-1)*M+1):(jjj*M))=G(jj,1:M);
end
end
for p=1:nT
for i=1:j
V(j,p)=V(j,p)+real(G(i,p)*exp(ic.*i./nd.*dw.*(p-1).*dt));
end
end
end
%绘制时程图
V1=V(1,:);
V2=V(2,:);
figure
subplot(3,1,1); %三行一列图
plot(t,V1,'k-');
xlabel('t/s');
ylabel('v/(m/s)');
axis([0 600 -30 30]);
title('谐波叠加法');
%功率谱
ft=1/dt; %一个模拟点的频率
[power,freq]=psd(V1,nT,ft,boxcar(512),0,'mean'); %boxcar中最好是2的指数倍,为了结果的协调所以调成800
power=power*2*dt; %psd算法归一化修正
%绘制频谱图
subplot(3,1,2);
loglog(freq,power,'k-',n,Sv,'r-'); %第1点的目标谱与模拟谱比较
xlabel('freq');
ylabel('power');
%绘制相干函数
V5=V(5,:);
[power15,freq15]=cpsd(V1,V5,[],512,512,2); %计算第1,2点的目标互功率谱
Sv15=[]; %第1,2点的目标互功率谱
for i=1:N
s15=S(1,5,i);
Sv15=[Sv15,s15];
end
subplot(3,1,3);
plot(freq15,abs(power15),'k',n,Sv15,'r');
axis([0 1 0 30]);
%转换为荷载
load Area; %各模拟点面积
for i=1:nd-1
v(i,:)=V(i,:)+U(i); %计算总风速
Force(i,:)=(v(i,:).*v(i,:))/1630*1.2.*Area(i,:)*10^3; %计算荷载
end
for i=nd:nd
v(i,:)=V(i,:)+U(i);
Force(i,:)=(v(i,:).*v(i,:))*3.986328;
end
%输出结果到excel
fid2=fopen('Windturbines_WAWSFFT_Davenport_data.csv','w');
TableU=[V(:,1:T);v(:,1:T);Force(:,1:T)];
for i=1:size(TableU,1)
fprintf(fid2,'\r\n'); j=1:size(TableU,2); fprintf(fid2,'%e,',TableU(i,j));
end
fclose(fid2);
评论8