查看: 478|回复: 0|关注: 0

[未答复] 基于Kaimal谱采用谐波合成法生成脉动风场

[复制链接]

新手

5 麦片

财富积分


050


1

主题

1

帖子

0

最佳答案
发表于 2018-1-14 17:07:05 | 显示全部楼层 |阅读模式
本帖最后由 锴胖纸 于 2018-1-14 18:47 编辑

各位师兄师姐:       你们好!
       最近在研究Matlab生成脉动风场,查阅了文献以及网上分享的几篇关于Matlab生成脉动风场的代码,都跑过,但是我觉得结果不对,先将代码分享如下(采用pwelch函数以及hamming窗口)

  1. %********************* 谐波叠加法模拟风速时程 修改1.0(采用Kaimal谱)**************
  2. clc
  3. clear
  4. %************************* 风速时程参数设定 *************************
  5. m=32;         %模拟风速点数
  6. N=2^11;       %频率采样点数
  7. dt=0.25;      %时间间隔
  8. omegaup=3*pi; %上限频率
  9. %************************  设定风速谱参数 *************************
  10. L=514.5;      %斜拉桥跨度
  11. z=45;         %风速测点离地面高度
  12. z0=0.03;      %地面粗糙度
  13. Uz=20;        %平均风速
  14. lambda=10;    %空间相关函数中的衰减系数
  15. K=0.4;        %Kaimal常数=0.4
  16. M=2*N;        %nfft傅里叶变换长度,取采样个数的2倍

  17. %形成风速时程矩阵
  18. v=zeros(m,M*m);    %创建m行,M*m列的时程矩阵
  19. u=zeros(m,M*m);
  20. v1=zeros(M*m,m);   %创建m行,M*m列的时程矩阵
  21. u2=zeros(M*m,m);
  22. t=dt*(0:1:(M*m-1));%创建时程横坐标时间点
  23. domega=(omegaup-0.001)/N;%频率间距
  24. D=zeros(m,m,N);   
  25. U1=K*Uz/log(z/z0); %测点位置的摩阻速度

  26. %形成目标谱
  27. omega1=0.01:domega:omegaup;%形成频率列表,初始化频率列表,初始频率为omegaup/N
  28. Sw1=200*U1^2.*z/Uz./(1+50.*omega1.*z./(2*pi*Uz)).^(5/3);%Kaimal谱密度表达式 水平向
  29. Su1=3.36*U1^2.*z/Uz./(1+10.*omega1.*z./(2*pi*Uz)).^(5/3);%L-P谱密度表达式 竖直向
  30. delta=12.9;
  31. %模拟风速测点间的距离(第一个间隔不取为零,否则会出现S不正定的情况)
  32. for j=1:m   %对模拟点风速的循环
  33.     rand('state',0);
  34.     thet=2*pi*rand(j,N);%生成随机相位
  35.     for l=1:N
  36.         omega(l)=(l-1)*domega+j/m*domega;
  37.     end
  38.     Sw=200*U1^2.*z/Uz./(1+50.*omega.*z./(2*pi*Uz)).^(5/3);
  39.     Su=3.36*U1^2.*z/Uz./(1+10.*omega.*z./(2*pi*Uz)).^(5/3);
  40.    

  41. %计算谱数据库矩阵,功率谱计算,Kaimal谱
  42.     for j1=1:m
  43.         for l=1:m
  44.             for k=1:N
  45.                 Coh(j1,l,k)=(exp(-lambda*omega(k)*delta/(2*pi*Uz)))^(abs(j1-l));%相关系数计算,采用Davenport形式。
  46.                 S(j1,l,k)=Sw(k)*Coh(j1,l,k);
  47.                 U(j1,l,k)=Su(k)*Coh(j1,l,k);
  48.             end
  49.         end
  50.     end

  51. %进行Cholesky分解
  52.     for i=1:1:N
  53.         H(:,:,i)=chol(S(:,:,i));
  54.         H(:,:,i)=H(:,:,i)';
  55.         Hu(:,:,i)=chol(U(:,:,i));
  56.         Hu(:,:,i)=Hu(:,:,i)';
  57.     end

  58. %填充谱数据矩阵D
  59.     D(:,j,:)=H(:,j,:);
  60.     E(:,j,:)=Hu(:,j,:);
  61.     i=sqrt(-1);
  62.     B1=sqrt(2*domega).*D(j,:,:);
  63.     Bu1=sqrt(2*domega).*E(j,:,:);
  64.     for ii=1:j
  65.         for jj=1:N
  66.             B2(ii,jj)=B1(1,ii,jj);
  67.             Bu2(ii,jj)=Bu1(1,ii,jj);
  68.         end
  69.     end
  70.     B2=B2.*exp(i.*thet);
  71.     Bu2=Bu2.*exp(i.*thet);
  72.     for jj=1:j
  73.         G(jj,1:M)=fft(B2(jj,:),M);
  74.         Gu(jj,1:M)=fft(Bu2(jj,:),M);
  75.         for jjj=2:m
  76.             G(jj,((jjj-1)*M+1):(jjj*M))=G(jj,1:M);
  77.             Gu(jj,((jjj-1)*M+1):(jjj*M))=Gu(jj,1:M);
  78.         end
  79.     end
  80. %谐波叠加生成模拟点的风速时程
  81.     for p=1:M*m
  82.         for k=1:j
  83.             v(j,p)=v(j,p)+real(G(k,p)*exp(i.*k./m.*domega.*(p-1).*dt));
  84.             u(j,p)=u(j,p)+real(Gu(k,p)*exp(i.*k./m.*domega.*(p-1).*dt));
  85.         end
  86.     end
  87. end %第一个for对应
  88. v1=v';u1=u';
  89. %风速时程绘图
  90. figure(1)
  91. subplot(3,1,1);
  92. plot(t(1:1:4096),v(1,1:1:4096));
  93. xlabel('t(s)');
  94. ylabel('v(t)');
  95. title('Point 1');
  96. subplot(3,1,2);
  97. plot(t(1:1:4096),v(2,1:1:4096));
  98. xlabel('t(s)');
  99. ylabel('v(t)');
  100. title('Point 2');
  101. subplot(3,1,3);
  102. plot(t(1:1:4096),v(m,1:1:4096));
  103. xlabel('t(s)');
  104. ylabel('v(t)');
  105. title('Point 32');

  106. figure(2)
  107. subplot(3,1,1);
  108. plot(t(1:1:4096),u(1,1:1:4096));
  109. xlabel('t(s)');
  110. ylabel('u(t)');
  111. title('Point 1');
  112. subplot(3,1,2);
  113. plot(t(1:1:4096),u(2,1:1:4096));
  114. xlabel('t(s)');
  115. ylabel('u(t)');
  116. title('Point 2');
  117. subplot(3,1,3);
  118. plot(t(1:1:4096),u(m,1:1:4096));
  119. xlabel('t(s)');
  120. ylabel('u(t)');
  121. title('Point 32');

  122. [power1,freq1]=pwelch(v(1,:),hamming(3024),10,M*m,20);
  123. [power2,freq2]=pwelch(v(2,:),hamming(3024),10,M*m,20);
  124. [power5,freq5]=pwelch(v(5,:),hamming(3024),10,M*m,20);

  125. %功率谱检验
  126. figure(3)                                                             %第1点模拟自功率谱和计算自功率谱比较(对数坐标)
  127. semilogx(freq1,freq1.*power1,'r',omega1,omega1.*Sw1,'b')
  128. xlabel('Frequency w(rad)')
  129. ylabel('Power wS(w)(rad.m2/s)')
  130. title('Point 1')

  131. figure(7)
  132. loglog(freq5,freq5.*power5,'r',omega1,omega1.*Sw1,'b')
  133. xlabel('Frequency w(rad)')
  134. ylabel('Power wS(w)(rad.m2/s)')
  135. title('Point 5')


  136. figure(8)
  137. plot(freq5,freq5.*power5,'r',omega1,omega1.*Sw1,'b')
  138. xlabel('Frequency w(rad)')
  139. ylabel('Power wS(w)(rad.m2/s)')
  140. title('Point 5')

复制代码
    跑出来的结果中 figure(1)窗口描述了第一个点、第二个点、最后一个点的风速时程,我看了几篇文献,表明因为第一个点以及第二个点间距为12.9米,所以得到的风速时程图像应当近似,然而这里图像不近似,所以想知道是什么原因造成的? (下图左1 程序得到的时程图像;右一 文献中相近几点的图像)

figure(1)

figure(1)
ewds.png






您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

站长推荐上一条 /3 下一条

快速回复 返回顶部 返回列表