基于Davenport风速谱单点时程模拟

转自：http://blog.sina.com.cn/s/blog_a1d5b9ba0102wfrw.html

%%%%%%%%%%%%%%%%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;  %%频率区间（0.001～6Hz）

v10=20;        %米高度风速

k=0.005;       %%地面粗糙度系数

delt=0.1;       %%时间增量

x=1200*f/v10;  %�venport谱中间参数

s1=4*k*v10^2*x.^2./f./(1+x.^2).^(4/3); %�venport谱表达式

subplot(2,2,1)

loglog(f,s1)  %%画Davenport经验谱图

xlabel('freq');

ylabel('S');

%%%%进行Cholesky分解%%%%%%%%%%%%%%%

for i=1:1:omegaup/d

H(i)=chol(s1(i));%%Cholesky分解

end

%%%%风荷载模拟%%%%%%%%%%%%%%%%%%%

thta=2*pi*rand(N,1);       %%介于0和6之间均匀分布的随机数

t=1:1:6000;               %%时间区间（0.1～600s）

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');