基於Davenport風速譜兩點時程模擬

風荷載時程曲線生成

從網絡收集而來的程序

基於Davenport風速譜兩點時程

clear all;
clc;
N=500;                                   %采樣點數
omegaup=4*2*pi;                   %截斷頻率                                                      
dm=omegaup/N;                    %頻率步長
dt=3.2*pi/(2*omegaup);          %時間步長0.2
k=0.00464;                             %地面粗糙度系數地面粗糙度等級A B C D：K= 0.00129、0.00215、0.00464、0.01291
d=0.001;
f=d:d:10;                                 %時間從0.001到10s，步進值為0.001
v10=26.8;                               %設計風速 26.8m/s——50年一遇十分鍾平均風速最大值
x=1200*f/v10;
s=4*k*v10*v10.*x.^2./f./(1+x.^2).^(4/3);         %venport風速譜

%x=1200*f/v10;
%Sv_Davenport=4*x.^2./(6*(1+x.^2).^(4/3));   %venport風速譜

%x1=z*f/(v10*(z/10)^2);
%Sv_Simiu=200*x1./(6*(1+50*x1).^(5/3));        %Simiu風速譜

%Sv_Kaimal=200*fn./(6*(1+50*fn).^(5/3));       %Kaimal風速譜

%x2=1800*f/v10;
%Sv_Harris=4*x2./(6.677*(2+x2.^2).^(5/6));    %Harris風速譜

%Lu=100*(z/30)^(1/2);
%x3=Lu*f/v0;
%Sv_Karman=4*x3./((1+70.8*x3.^2).^(5/6));    %Karman風速譜 

z1=10;                                    %取第一點為10米高度
z2=50;                                    %取第二點為50米高度
r=0.22;                                    %考慮地表粗糙度影響的無量綱冪指數，按中國規范取0.22-c類
v5=v10*(z2/z1)^(r);                 %計算n米高處的平均風速——52.8m處平均風速期望值
Cz=10;                                      %Cz表示z方向的指數衰減系數(取經驗值)           
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);            %1點自互功率譜
       x2=1200*m2/(2*pi*v10);
       s22=2*pi*4*k*v5*v5*x2*x2./m2./(1+x2*x2).^(4/3);              %2點自互功率譜
       s12=sqrt(s11*s22).*exp(-2*m2*Cz*abs(z1-z2)./(2*pi*(v10+v5))); %僅考慮豎向相關性互功率譜
       s21=sqrt(s11*s22).*exp(-2*m1*Cz*abs(z1-z2)./(2*pi*(v10+v5))); %僅考慮豎向相關性互功率譜
       S=[s11 s12;s21 s22];   %譜矩陣
       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=0.8*sum1;
    sum2=0.8*sum2;
    v1(j,K)=sum1; 
    v2(j,K)=sum2;
  end
end
u1=v1+ v5;
u2=v2+ v5;
t= (0:2*N-1)*dt;
subplot(2,2,1);
plot(t,u1,'b-');
xlabel('t(s)');
ylabel('v(t)');
axis([-1 180 0 80]);

subplot(2,2,3);
plot(t,u2,'r-');
xlabel('t(s)');
ylabel('v(t)');
axis([-1 180 0 80]);

Y=fft(v2);                           %對數值解作傅立葉變換
Y(1)=[];                             %去掉零頻量
m=length(Y)/2;                       %計算頻率個數；
power=150*abs(Y(1:m)).^2/(length(Y).^2);      %計算功率譜
freq=5*(1:m)/length(Y);       %計算頻率，因為時間步長為0.125，而不是1，故乘以8
subplot(2,2,2);
loglog(freq,power,'r-',f,s,'g-');    %對數顯示，比較plot(t,v2,'r-');
axis([-50 15 -10 1000]);       %畫出坐標軸比例axis([xmin xmax ymin ymax])
set(gca,'xtick',[0.1 1 10]);         %自動定義刻度
set(gca,'ytick',[0.1 1 10]);
grid on
xlabel('頻率');
ylabel('功率');