close all
clear all
clc
x=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 ];
x=x';
y=[4.02725 3.3011 2.8372 2.4234 2.06155 1.8402 1.7439 1.645 1.546 1.44 1.3349 1.2369 1.1397 1.0046];
y=y';
z=[3.75 9.75 13.75 17.25 21.5 25.25 27.825 30.475 33.125 35.95 38.775 41.4 44 46.05];
z=z';
h=14;%节点个数
tqd=1;
v10=29.665;%10米处风速
k=1:h;
v(k)=v10*(z(k)/10).^0.16; %任意高度处的风速指数风速廓线
Tstep=6000;%步长的个数
dt=0.1;%步长的大小
f=0.01:0.01:10;%脉动风频率
u=0.005*v10^2;
P=4;
t=(1:6000)*dt;
nt=length(t);
%以下是求RO~R4
R0=zeros(h);
for i=1:h
for j=i:h
dx=x(i)-x(j);
dy=y(i)-y(j);
dz=z(i)-z(j);
h1=z(i)*f/v10/(z(i)/10)^0.16;
h2=z(j)*f/v10/(z(j)/10)^0.16;
w1=f.*(1+50*h1).^(5/3);
w2=f.*(1+50*h2).^(5/3);
s1=200*h1*u./w1; %求kaimal谱
s2=200*h2*u./w2;
xiangguanhanshu=exp(-sqrt(dx.^2./50^2+dy.^2./50^2+dz.^2./60^2)).*cos(2*pi.*f*0.*dt); %求相关函数
y=sqrt(s1.*s2).*xiangguanhanshu; %求互功率谱水平顺风向和竖向互谱
R0(i,j)=trapz(f,y);
R0(j,i)=R0(i,j);
end
end
R1=zeros(h);
for i=1:h
for j=i:h
dx=x(i)-x(j)