脉动风时程matlab程序Word文档格式.doc

脉动风时程matlab程序Word文档格式.doc

《脉动风时程matlab程序Word文档格式.doc》由会员分享，可在线阅读，更多相关《脉动风时程matlab程序Word文档格式.doc（4页珍藏版）》请在冰豆网上搜索。

（12）

其中，风谱在频率范围内划分成个相同部分，为频率增量，为上述下三角矩阵的模，为两个不同作用点之间的相位角，为介于和之间均匀分布的随机数，是频域的递增变量。

文中模拟开孔处的来流风，因而只作单点模拟。

即式（4）可简化为：

（13）

本文采用Davenport水平脉动风速谱：

（14）

式中，脉动风速功率谱；

脉动风频率（Hz）；

地面粗糙度系数；

标准高度为10m处的风速（m/s）。

Matlab程序:

N=10;

d=0.001;

n=d:

d:

N;

%%频率区间（0.01～10）

v10=16;

k=0.005;

x=1200*n/v10;

s1=4*k*v10^2*x.^2./n./（1+x.^2）.^（4/3）;

%%Davenport谱

subplot（2,2,1）

loglog（n,s1）%%画谱图

axis（[-10015-1001000]）

xlabel（'

freq'

）;

ylabel（'

S'

fori=1:

1:

N/d

H（i）=chol（s1（i））;

%%Cholesky分解

end

thta=2*pi*rand（N/d,1000）;

%%介于0和2pi之间均匀分布的随机数

t=1:

1000;

%%时间区间（0.1～100s）

forj=1:

1000

a=abs（H）;

b=cos（（n*j/10）+thta（:

j）'

c=sum（a.*b）;

v（j）=（2*d）.^（1/2）*c;

%%风荷载模拟

subplot（2,2,2）

plot（t/10,v）%%显示风荷载

t（s）'

v（t）'

Y=fft（v）;

%%对数值解作傅立叶变换

Y

（1）=[];

%%去掉零频量

m=length（Y）/2;

%%计算频率个数；

power=abs（Y（1:

m））.^2/（length（Y）.^2）;

%%计算功率谱

freq=10*（1:

m）/length（Y）;

%%计算频率，因为步长为0.1，而不是1，故乘以10

subplot（2,2,3）

loglog（freq,power,'

r'

n,s1,'

b'

）%%比较

对源程序的修改：

z=xcorr（v）;

Y=fft（z）;

楼主的修改使模拟得到的功率谱与源谱的数量级对上了，但是吻合不是太好。

但是好像这样做是不对的。

求信号x（t）的功率谱有两种方法，一是对X（t）做傅立叶变换，再平方

S=abs（fft（x））^2

一是先对X（t）求相关系数，再进行傅立叶变换:

S=fft（xcorr（X））

楼主的方法好像是这两个方法的混合。

欢迎大家拍砖^_^