风场模拟的MATLAB程序
clear tic v10=24;z=116;n=0.001:0.001:10; xn=z*n./v10; k=0.003;%地面粗糙长度 ti=0.1;%时间步长s1=200*k*v10^2*xn./n./(1+50*xn).^(5/3);%K谱
%产生空间点坐标 for i=1:20
x(i)=5+i;
z(i)=8+i; end
%求R矩阵 syms f
R0=zeros(20); for i=1:20
for j=i:20
H0=inline('200*k* v10^2*116*f./v10./f./(1+(50*116*f./v10)).^(5/3).*exp(-sqrt(dx^2/50^2+dz^2/60^2))','f','k','dx','dz','v10');
dx=x(i)-x(j);
dz=z(i)-z(j);
R0(i,j)=quadl(H0,0.001,10,0.001,0,k,dx,dz,v10);
R0(j,i)=R0(i,j);
end
end
R1=zeros(20); for i=1:20
for j=i:20
H1=inline('200*k* v10^2*116*f./v10./f./(1+(50*116*f./v10)).^(5/3).*exp(-sqrt(dx^2/50^2+dz^2/60^2)).*cos(2*pi*f*1*ti)','f','k','dx','dz','ti','v10');
dx=x(i)-x(j);
dz=z(i)-z(j);
R1(i,j)=quadl(H1,0.001,10,0.001,0,k,dx,dz,ti,v10);
R1(j,i)=R1(i,j);
end
end
R2=zeros(20); for i=1:20
for j=i:20
H2=inline('200*k* v10^2*116*f./v10./f./(1+(50*116*f./v10)).^(5/3).*exp(-sqrt(dx^2/50^2+dz^2/60^2)).*cos(2*pi*f*1*ti)','f','k','dx','dz','ti','v10');
dx=x(i)-x(j);
dz=z(i)-z(j);
R2(i,j)=quadl(H2,0.001,10,0.001,0,k,dx,dz,ti,v10);
R2(j,i)=R2(i,j);
end
end
R3=zeros(20); for i=1:20
for j=i:20
H3=inline('200*k* v10^2*116*f./v10./f./(1+(50*116*f./v10)).^(5/3).*exp(-sqrt(dx^2/50^2+dz^2/60^2)).*cos(2*pi*f*1*ti)','f','k','dx','dz','ti','v10');
dx=x(i)-x(j);
dz=z(i)-z(j);
R3(i,j)=quadl(H3,0.001,10,0.001,0,k,dx,dz,ti,v10);
R3(j,i)=R3(i,j);
end
end
R4=zeros(20); for i=1:20
for j=i:20
H4=inline('200*k* v10^2*116*f./v10./f./(1+(50*116*f./v10)).^(5/3).*exp(-sqrt(dx^2/50^2+dz^2/60^2)).*cos(2*pi*f*1*ti)','f','k','dx','dz','ti','v10');
dx=x(i)-x(j);
dz=z(i)-z(j);
R4(i,j)=quadl(H4,0.001,10,0.001,0,k,dx,dz,ti,v10);
R4(j,i)=R4(i,j);
end
end
Q1=zeros(20); Q2=zeros(20); Q3=zeros(20); Q4=zeros(20); A=;%80X80矩阵 B=;%80X20矩阵 X=A\B;
%此式相当于A*X=B，X(80×20矩阵)为自回归系数Ψ q1=X(1:20,:);%取X的第一个20×20矩阵 q2=X(20+1:2*20,:);%取X的第二个20×20矩阵 q3=X(2*20+1:3*20,:);%取X的第三个20×20矩阵 q4=X(3*20+1:4*20,:);%取X的第四个20×20矩阵 Q1=q1'; Q2=q2'; Q3=q3'; Q4=q4'; RN=R0+Q1*R1+Q2*R2+Q3*R3+Q4*R4;%计算时出现
L=zeros(20);L=chol(RN);L=L';a=zeros(20,2048);
for i=1:20
for j=1:2048
a(i,j)=normrnd(0,1,1,1);
%产生均值0，方差1的正态随机数矩阵
endend
V(1:20,1)=L*a(:,1); V(1:20,2)=Q1*V(1:20,1)+L*a(:,2); V(1:20,3)=(Q1*V(1:20,2)+Q2*V(1:20,1))+L*a(:,3); V(1:20,4)=(Q1*V(1:20,3)+Q2*V(1:20,2)+Q3*V(1:20,1))+L*a(:,4); for t=5:2048
V(1:20,t)=(Q1*V(1:20,t-1)+Q2*V(1:20,t-2)+Q3*V(1:20,t-3)+Q4*V(1:20,t-4))+L*a(:,t); end
V1=V(2,:);%取第一点的风速
t=(1:2048)*ti; figure subplot(2,1,1); plot(t,V1,'b-'); xlabel('t(s)'); ylabel('v(t)'); axis(); %与目标谱进行比较 %=psd(V1,2048,10,boxcar(1024),0,'mean'); =cpsd(V1,V1,blackman(1024),0.5,2048,10);power=power.*2./10;%规一化修正subplot(2,1,2); loglog(freq,power,'r-',n,s1,'g-'); xlabel('freq'); ylabel('s1/power'); toc
页:
[1]