谐波合成法matlab,基于Kaimal谱采用谐波合成法生成脉动风场.pdfVIP

谐波合成法matlab,基于Kaimal谱采⽤谐波合成法⽣成脉动风 场 %********************* 谐波叠加法模拟风速时程 修改1.0(采⽤ aimal谱)************** clc clear %************************* 风速时程参数设定 ************************* m=32; %模拟风速点数 N=2^ 11; %频率采样点数 dt=0.25; %时间间隔 omegaup=3*pi; %上限频率 %************************ 设定风速谱参数 ************************* L=5 14.5; %斜拉桥跨度 z=45; %风速测点离地⾯⾼度 z0=0.03; %地⾯粗糙度 Uz=20; %平均风速 lambda= 10; %空间相关函数中的衰减系数 =0.4; % aimal常数=0.4 M=2*N; %nfft傅⾥叶变换长度，取采样个数的2倍 %形成风速时程矩阵 v=zeros(m,M*m); %创建m⾏，M*m列的时程矩阵 u=zeros(m,M*m); v 1=zeros(M*m,m); %创建m⾏，M*m列的时程矩阵 u2=zeros(M*m,m); t=dt* (0:1:(M*m- 1));%创建时程横坐标时间点 domega=(omegaup-0.00 1)/N;%频率间距 D=zeros(m,m,N); U1= *Uz/log(z/z0); %测点位置的摩阻速度 %形成 ⽬标谱 omega1=0.0 1:domega:omegaup;%形成频率列表，初始化频率列表，初始频率为omegaup/N Sw 1=200*U1^2.*z/Uz./ (1+50.*omega1.*z./ (2*pi*Uz)).^ (5/3);% aimal谱密度表达式 ⽔平向 Su1=3.36*U1^2.*z/Uz./ (1+ 10.*omega1.*z./ (2*pi*Uz)).^ (5/3);%L-P谱密度表达式 竖直向 delta= 12.9; %模拟风速测点间的距离(第⼀个间隔不取为零，否则会出现S不正定的情况) for j = 1:m %对模拟点风速的循环 rand(state,0); thet=2*pi*rand(j ,N);%⽣成随机相位 for l= 1:N omega(l)=(l- 1)*domega+j /m*domega; end Sw=200*U1^2.*z/Uz./ (1+50.*omega.*z./ (2*pi*Uz)).^ (5/3); Su=3.36*U1^2.*z/Uz./ (1+ 10.*omega.*z./ (2*pi*Uz)).^ (5/3); %计算谱数据库矩阵，功率谱计算， aimal谱 for j 1= 1:m for l= 1:m for k= 1:N Coh(j 1,l,k)=(exp(-lambda*omega(k)*delta/ (2*pi*Uz)))^ (abs(j 1-l));%相关系数计算，采⽤Davenport形式。 S(j 1,l,k)=Sw(k)*Coh(j 1,l,k); U(j 1,l,k)=Su(k)*Coh(j 1,l,k); end end end %进⾏Cholesky分解 for i= 1:1:N H(:,:,i)=chol(S(:,:,i)); H(:,:,i)=H(:,:,i); Hu(:,:,i)=chol(U(:,:,i)); Hu(:,:,i)=Hu(:,:,i); end %填充谱数据矩阵D D(:,j ,:)=H(:,j ,:); E(:,j ,:)=Hu(:,j ,:); i=sqrt(- 1); B1=sqrt(2*domega).*D(j ,:,:); Bu1=sqrt(2*domega).*E(j ,:,:); for ii= 1:j for jj = 1:N B2(ii,jj )=B1(1,ii,jj ); Bu2(ii,jj )=Bu1(1,ii,jj ); end end B2=B2.*exp(i.*thet); Bu2=Bu2.*exp(i.*thet); for jj = 1:j G(jj ,1:M