根据风的记录,脉动风可作为高斯平稳过程来考虑。观察n 个具有零均值的平稳高斯过程,其谱密度函数矩阵为:
?
?
???
??
?????=)(...)
()(............)(...)
()()(...)()()(21222
21
11211ωωωωωωωωωωnn n n n n s s s s s s s s s S (9)
将)(ωS 进行Cholesky 分解,得有效方法。
T H H S )()()(*ωωω?= (10)
其中,
?
?
???
?
?
?????=)(...)
()(............0...)
()(0
...
0)
()(2122
21
11ωωωωωωωnn n n H H H H H H H (11)
T H )(*ω为)(ωH 的共轭转置。
根据文献[8],对于功率谱密度函数矩阵为)(ωS 的多维随机过程向量,模拟风速具有如下形式:
[]
∑∑==++???=j
m N
l ml l jm l l jm j t H t v 11)(cos 2)()(θωψωωω
n j ...,3,2,1= (12)
其中,风谱在频率范围内划分成N 个相同部分,N ωω=?为频率增量,)(l jm H ω为上述下三角矩阵的模,)(l jm ωψ为两个不同作用点之间的相位角,ml θ为介于0和π2之间
均匀分布的随机数,ωω??=l l 是频域的递增变量。
文中模拟开孔处的来流风,因而只作单点模拟。即式(4)可简化为:
[]∑=+???=N
l l l l t H t v 1
cos 2)()(θωωω (13)
本文采用Davenport 水平脉动风速谱:
3
/422
210
)1(4)(x n kx v
n S v += (14) 式中,--)(n S v 脉动风速功率谱;
--n 脉动风频率(Hz);
--k 地面粗糙度系数;
;1200
10v n x
--10v 标准高度为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([-100 15 -100 1000]) xlabel('freq'); ylabel('S');
for i=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 )
for j=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([-100 15 -100 1000]) xlabel('freq'); ylabel('S');
10
10
10
10
2
freq
S
50100
-20
-10010
20
t(s)
v (t )
10
-2
100
10
freq
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([-100 15 -100 1000]) xlabel('freq'); ylabel('S');
楼主的修改使模拟得到的功率谱与源谱的数量级对上了,但是吻合不是太好。但是好像这样做是不对的。
求信号x(t)的功率谱有两种方法,一是对X(t)做傅立叶变换,再平方 S=abs(fft(x))^2
一是先对X(t)求相关系数,再进行傅立叶变换:
S=fft(xcorr(X))
楼主的方法好像是这两个方法的混合。欢迎大家拍砖^_^