脉动风时程matlab程序.doc
- 文档编号:7393134
- 上传时间:2023-05-11
- 格式:DOC
- 页数:4
- 大小:153KB
脉动风时程matlab程序.doc
《脉动风时程matlab程序.doc》由会员分享,可在线阅读,更多相关《脉动风时程matlab程序.doc(4页珍藏版)》请在冰点文库上搜索。
根据风的记录,脉动风可作为高斯平稳过程来考虑。
观察个具有零均值的平稳高斯过程,其谱密度函数矩阵为:
(9)
将进行Cholesky分解,得有效方法。
(10)
其中,
(11)
为的共轭转置。
根据文献[8],对于功率谱密度函数矩阵为的多维随机过程向量,模拟风速具有如下形式:
(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:
1:
1000;%%时间区间(0.1~100s)
forj=1:
1:
1000
a=abs(H);
b=cos((n*j/10)+thta(:
j)');
c=sum(a.*b);
v(j)=(2*d).^(1/2)*c;%%风荷载模拟
end
subplot(2,2,2)
plot(t/10,v)%%显示风荷载
xlabel('t(s)');
ylabel('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')%%比较
axis([-10015-1001000])
xlabel('freq');
ylabel('S');
对源程序的修改:
z=xcorr(v);
Y=fft(z);%%对数值解作傅立叶变换
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')%%比较
axis([-10015-1001000])
xlabel('freq');
ylabel('S');
楼主的修改使模拟得到的功率谱与源谱的数量级对上了,但是吻合不是太好。
但是好像这样做是不对的。
求信号x(t)的功率谱有两种方法,一是对X(t)做傅立叶变换,再平方
S=abs(fft(x))^2
一是先对X(t)求相关系数,再进行傅立叶变换:
S=fft(xcorr(X))
楼主的方法好像是这两个方法的混合。
欢迎大家拍砖^_^
- 配套讲稿:
如PPT文件的首页显示word图标,表示该PPT已包含配套word讲稿。双击word图标可打开word文档。
- 特殊限制:
部分文档作品中含有的国旗、国徽等图片,仅作为作品整体效果示例展示,禁止商用。设计者仅对作品中独创性部分享有著作权。
- 关 键 词:
- 脉动 风时程 matlab 程序