论坛 › K01：MatLab › 高手看看小弟的模拟风荷载时程分析程序的模拟谱与目标谱不吻合

jizhiyangzzu 发表于 2008-11-9 15:24:36

高手看看小弟的模拟风荷载时程分析程序的模拟谱与目标谱不吻合

clear all;
N=512;                                          %采样点数
wu=4*2*pi;                                    %截断频率                                                      
dm=wu/N;                                        %频率步长
dt=2*pi/(2*wu);                                 %时间步长
k=0.005;                                        %地面粗糙度系数
d=0.001;
f=d:d:10;
v10=30;                                       %10米高处风速为30m/s
x=1200*f/v10;
s=4*k*v10*v10.*x.^2./f./(1+x.^2).^(4/3);
z1=10;                                          %取第一点为10米高度
z2=9.5;                                       %取第二点为5米高度
r=0.12;                                       %考虑地表粗糙度影响的无量纲幂指数
v5=v10*(z2/z1)^r;                               %计算n米高处的平均风速
C=10;                                           %指数衰减系数(取经验值)            
v1=zeros(2*N,1);
v2=zeros(2*N,1);
thta1=rand(N,1);
thta2=rand(N,1);                              %随机数
node=1;
for K=1:node
for j=1:2*N
   sum1=0;
   sum2=0;
      for l=1:N
      m1=l*dm-0.5*dm;
      m2=l*dm;
      x1=1200*m1/(2*pi*v10);
      s11=2*pi*4*k*v10*v10*x1*x1./m1./(1+x1*x1).^(4/3);
      x2=1200*m2/(2*pi*v10);
      s22=2*pi*4*k*v5*v5*x2*x2./m2./(1+x2*x2).^(4/3);
      s12=sqrt(s11*s22).*exp(-2*m2*C*abs(z1-z2)./(2*pi*(v10+v5)));
      s21=sqrt(s11*s22).*exp(-2*m1*C*abs(z1-z2)./(2*pi*(v10+v5)));
      S=;
      H=chol(S);
      a1=abs(H(1,1));
      H1=H';
      a21=abs(H1(2,1));
      a22=abs(H1(2,2));
      b1=cos((m1*dt*(j-1))+2*pi*thta1(l,1));
      b2=cos((m2*dt*(j-1))+2*pi*thta2(l,1));
      c1=a1*b1;
      c21=a21*b1;
      c22=a22*b2;
      d1=(dm).^0.5*c1;
      d2=(dm).^0.5*(c21+c22);
      sum1=sum1+d1;
      sum2=sum2+d2;
      end
   sum1=2*sum1;
   sum2=2*sum2;
   v1(j,K)=sum1;
   v2(j,K)=sum2;
end
end
t=(0:2*N-1)*dt;
subplot(2,2,1);
plot(t,v1,'b-');
xlabel('t(s)');
ylabel('v(t)');
hold on
plot(t,v2,'r-');
xlabel('t(s)');
ylabel('v(t)');
axis([-1 120 -60 60]);

Y=fft(v1);                            %对数值解作傅立叶变换
Y(1)=[];                              %去掉零频量
m=length(Y)/2;                        %计算频率个数；
power=abs(Y(1:m)).^2/(length(Y).^2);%计算功率谱
freq=8*(1:m)/length(Y);   %计算频率，因为时间步长为0.125，而不是1，故乘以8
subplot(2,2,2);
loglog(freq,power,'r-',f,s,'g-');   %比较
axis([-50 15 -10 1000]);
set(gca,'xtick',);          %自动定义刻度
set(gca,'ytick',);
xlabel('freq');
ylabel('P');

ismarlia 发表于 2010-7-23 12:27:43

很受启发，谢谢分享，希望能讨论一下

xiaojunjie 发表于 2010-9-14 21:55:54

楼主强大。。。

wuyuan 发表于 2010-9-16 13:56:35

你的功率谱计算是为什么是除以(length(Y).^2)?matlab中的periodogram 指令貌似只是除以length(Y)。D:\Snap1.bmp

zhaosy05 发表于 2011-4-2 20:04:28

把length（Y）的平方号去掉就行了

zhaosy05 发表于 2011-4-3 16:04:35

同时功率谱需要做归一化处理

zhaosy05 发表于 2011-4-3 16:03:38

同时需要对PW作归一化处理

cjling 发表于 2011-4-7 20:13:00

受教了，楼主强人！！

李新11 发表于 2012-3-17 16:29:30

我也学习一下

zhangweina 发表于 2012-4-24 20:52:31

请问一下，您上面的程序中，频率增量为什么选取0.001，有什么依据吗？

zhangweina 发表于 2012-4-24 21:15:11

是呢，这个应该对功率谱做归一化

zhangweina 发表于 2012-4-24 21:27:31

请问一下，程序中的截断频率是如何选取的？

查看完整版本: 高手看看小弟的模拟风荷载时程分析程序的模拟谱与目标谱不吻合