%%%%%%%%%%%%%%%%main procedure%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%单点模拟%%%%%%%%%%%%%%%%%%%%%%%%%%
N=6000;
domega=0.001;
omegaup=2*pi;
n= domega: domega:omegaup;%%频率区间(0.001~6)
v10=16;
k=0.005;
x=1200*n/v10;
delt=0.1;
s1=4*k*v10^2*x.^2./n./(1+x.^2).^(4/3);%�venport谱
subplot(2,2,1)
loglog(n,s1)%%画谱图
xlabel('freq');
ylabel('S');
for i=1:1:omegaup/domega
H(i)=chol(s1(i));%%Cholesky分解
end
thta=2*pi*rand(N,1);%%介于0和2pi之间均匀分布的随机数
t=0.1:0.1:600;%%时间区间(0.1~600s)
ii=sqrt(-1);
for j=1:1:N
v(j)=sqrt(2*domega) *H(j)*exp(ii*thta(j));%%风荷载模拟
end
Y=fft(v,N);%%对数值解作傅立叶变换
for j=1:1:N
vh(j)=real((Y(j)*exp(ii*j*domega)));
end
%[power1,freq1]= psd(vh,N,2,boxcar(512),0,'mean');
[power,freq]=pwelch(vh,boxcar(3024),10,N,1/delt);
subplot(2,2,2)
plot(vh)%%显示风荷载
xlabel('t(s)');
ylabel('v(t)');
axis([0 1800 -10 10]);
subplot(2,2,3)
loglog(freq,power,'r',n,s1,'b')%%比较
xlabel('freq');
ylabel('S');
%%%%main procedure%%%%%%%%%%%%%%%%%%%%
clc
clear
close all
N=6000;
d=0.001;
omegaup=6;
f=d:d:omegaup;
v10=20;
k=0.005;
delt=0.1;
x=1200*f/v10;
s1=4*k*v10^2*x.^2./f./(1+x.^2).^(4/3); %�venport谱表达式
subplot(2,2,1)
loglog(f,s1)
xlabel('freq');
ylabel('S');
%%%%进行Cholesky分解%%%%%%%%%%%%%%%
for i=1:1:omegaup/d
H(i)=chol(s1(i));%%Cholesky分解
end
%%%%风荷载模拟%%%%%%%%%%%%%%%%%%%
thta=2*pi*rand(N,1);
t=1:1:6000;
ii=sqrt(-1);
for j=1:1:N
v(j)=sqrt(2*d)*H(j)*exp(ii*thta(j));%%风荷载模拟
end
%%%%%对风速时程进行FFT变换%%%%%%%%%%%%
Y=fft(v,N);
for i=1:1:N
vh(i)=real(Y(i)*exp(ii*i*d*0.1));
end
[power,freq]=pwelch(vh,boxcar(3024),10,N,1/delt);
subplot(2,2,2)
plot(t/10,vh)
xlabel('t(s)');
ylabel('Y(t)');
%%%%拟合谱与目标谱比较%%%%%%%%%%%%%%
subplot(2,2,3)
loglog(freq,power,'r',f,s1,'b')
xlabel('freq');
ylabel('S');