脉动风时程matlab程序.doc

根据风的记录，脉动风可作为高斯平稳过程来考虑。观察 n 个具有零均值的平稳高斯过 程，其谱密度函数矩阵为 . . . . . . . . . . . . . . . 21 22221 11211 nnnn n n sss sss sss S ( 9) 将 S 进行 Cholesky分解，得有效方法。 THHS * ( 10) 其中， . . . . . . . . . 0. . . 0. . .0 21 2221 11 nnnn HHH HH H H ( 11) TH * 为 H 的共轭转置。 根据文献 8，对于功率谱密度函数矩阵为 S 的多维随机过程向量，模拟风速具有如 下形式 jm Nl mlljmlljmj tHtv 1 1 co s2 nj .,3,2,1 ( 12) 其中，风谱在频率范围内划分成 N 个相同部分， N 为频率增量， ljmH 为 上述下三角矩阵的模 ， ljm 为两个不同作用点之间的相位角， ml 为介于 0 和 2 之间 均匀分布的随机数， ll 是频域的递增变量。 文中模拟开孔处的来流风，因而只作单点模拟。即式( 4)可简化为 Nl lll tHtv 1 c o s2 ( 13) 本文采用 Davenport水平脉动风速谱 3/42 22 10 1 4 xn kxvnS v ( 14) 式中， nSv 脉动风速功率谱； n 脉动风频率 Hz； k 地面粗糙度系数； ;1200 10vnx 10v 标准高度为 10m处的风速 m/s。 Matlab程序 N10; d0.001; nddN;频率区间( 0.01 10) v1016; k0.005; x1200*n/v10; s14*k*v102*x.2./n./1x.2.4/3;Davenport谱 subplot2,2,1 loglogn,s1画谱图 axis-100 15 -100 1000 xlabelfreq; ylabelS; for i11N/d Hichols1i;Cholesky分解 end thta2*pi*randN/d,1000;介于 0和 2pi之间均匀分布的随机数 t111000;时间区间( 0.1 100s) for j111000 aabsH; bcosn*j/10thta,j; csuma.*b; vj2*d.1/2*c;风荷载模拟 end subplot2,2,2 plott/10,v显示风荷载 xlabelts; ylabelvt; Yfftv;对数值解作傅立叶变换 Y1;去掉零频量 mlengthY/2;计算频率个数； powerabsY1m.2/lengthY.2;计算功率谱 freq10*1m/lengthY;计算频率，因为步长为 0.1，而不是 1，故乘以 10 subplot2,2,3 loglogfreq,power,r,n,s1,b比较 axis-100 15 -100 1000 xlabelfreq; ylabelS; 10 -2 10 0 10 0 10 2 f r e q S 0 50 100 - 2 0 - 1 0 0 10 20 t s v t 10 -2 10 0 10 0 f r e q S 对源程序的修改 zxcorrv; Yfftz;对数 值解作傅立叶变换 Y1;去掉零频量 mlengthY/2;计算频率个数； powerabsY1m.2/lengthY.2;计算功率谱 freq10*1m/lengthY;计算频率，因为步长为 0.1，而不是 1，故乘以 10 subplot2,2,3 loglogfreq,power,r,n,s1,b比较 axis-100 15 -100 1000 xlabelfreq; ylabelS; 楼主的 修改使模拟得到的功率谱与源谱的数量级对上了，但是吻合不是太好。但是好像这样 做是不对的。 求信号 xt的功率谱有两种方法，一是对 Xt做傅立叶变换，再平方 Sabsfftx2 一是先对 Xt求相关系数，再进行傅立叶变换 SfftxcorrX 楼主的方法好像是这两个方法的混合。欢迎大家拍砖 _