风阵模拟

引言

《建筑结构荷载规范》在计算风荷载时提到了风振系数，结合实际工程使用，结构上的风荷载可分为两种成分：平均风和脉动风。对应地，风对结构的作用也有静力的平均风作用和动力的脉动风作用。平均风的作用可用静力方法计算，而脉动风是随机荷载，它引起结构的振动，一般采用随机振动理论对其振动进行分析。

风振系数是指风对建筑物的作用是不规则的，风压随风速、风向的紊乱变化而不停地改变，它综合考虑了结构在风荷载作用下的动力响应，其中包括风速随时间、空间的相关性，结构的阻尼特性等因素。通常把风作用的平均值看成稳定风压或平均风压，实际风压是在平均风压上下波动的。平均风压使建筑物产生一定的侧移，而波动风压使建筑物在该侧移附近左右振动。对于高度较大，刚度较小的高层建筑，波动风压会产生不可忽略的动力效应，在设计中必须考虑。目前采用加大风荷载的办法来考虑这个动力效应，在风压值上乘以风振系数。当房屋高度大于30m、高宽比大于1.5时，以及对于构架、塔架、烟囱等高耸结构，均应考虑风振。（ PS:对于30m以下且高宽比小于1.5的房屋建筑，可以不考虑脉动风压影响，此时风振系数取βz=1.0）

规范公式理解

g---峰值因子，取值为2.5

I10---10m高度名义湍流强度，对应A、B、C和D类地面粗糙度，可分别取0.12、0.14、0.23、0.39

R---脉动风荷载的共振分量因子；

其中:

f1---结构第一阶自振频率（Hz）

kw---地面粗糙度修正系数，对A类、B类、C类和D类地面的粗糙度分别取为1.28、1.0、0.54、0.26

结构阻尼比：钢结构取0.01，有填充墙的钢结构取0.02，钢筋混凝土结构取0.05，对其他的一些结构根据相应的一些工程经验进行取值

Bz---脉动风荷载的背景分量因子

1、对于体型和质量沿高度均匀分布的高层建筑个高耸结构可以按照下面的公式来进行相应的计算（通过后面的附录知道，原公式是一个多重积分，经过大量的试算和回归分析，采用非线性最小二乘法拟合得到下面的简化公式，考虑了迎风面和背风面的风压相关性，结合相应的工程经验取0.7的折减系数）：

2、对迎风面和侧风面的宽度沿高度按直线或接近直线变化，而质量沿高度按连续规律变化的高耸结构，8.4.5.1计算的背景分量因子应乘以修正系数。其中一个为构筑物在z高度处的迎风面宽度B(z)与底部宽度B(0)的比值；另外一个按照表8.4.5-2确定。

脉动风荷载的空间相关系数

竖直方向的相关系数按照下面的式子计算：

H---结构总高度，对于A 、B、C、D类地面粗糙度，H的取值分别不应大于300m，350m，450m和550m

2.水平方向的相关系数计算值：

在上面的式子中，B为结构迎风面宽度

3、对于迎风面宽度较小的高耸结构，水平方向的相关系数可以取1

阵型系数应该根据结构的动力计算确定，对外形、质量、刚度沿高度按连续规律变化的竖向悬臂型高耸结构以及沿高度比较均匀的高层建筑，阵型系数也可以根据相对高度z/H按照荷载规范的附录G确定。

附录详解风阵系数

当房屋高度大于30m、高宽比大于1.5时，频谱比较的稀疏，第一阵型起到绝对的作用，这个时候仅考虑结构的第一阵型并通过下面的式子来计算

这两个公式结合起来就可以得到规范中共振因子的公式（见上面）

结构阵型系数的确定

结构阵型系数应该按照结构动力分析确定，用简化的计算公式，主要是根据相应结构变形特点。

理论知识储备

Davenport（译名：达文波特）根据世界上不同地点、不同高度测得到90多次的强风记录，并假定水平阵风谱中的湍流积分尺度L沿高度不变，取常数值1200m，并取脉动风速谱为不同离地高度实测值的平均值，建立了经验数学表达式。

我国规范及在风工程应用中一般采用Davenport脉动风速谱。Davenport谱比其它谱偏大，而谱值偏大的范围正好是风频率与结构物自振频率接近的地方，影响较大，故Davenport风速谱可能会高估结构的动力响应，其结果可能会偏于保守，但是在结构抗风的设计角度而言，却提高了结构的安全度。

高斯平稳过程：是随机过程的一种，是一系列服从正态分布的随机变量在一指数集内的组合

引入：一个随机振动过程的特征可以用数学期望、方差和相关函数来描述。在工程技术问题中，广泛采用从频率域来描述一个随机振动过程特征的功率谱函数(理解并加以运用)，功率谱密度函数能够反映随机振动的功率关于频率的分布密度（实质内涵）。

一、频谱密度

频谱密度：设一个能量信号为s(t),则它的频谱密度S（w）可以由傅氏变换求得：

能量信号的频谱密度S（f）和功率信号C（jnw）（比如一个周期信号）的频谱主要区别有：

S（f）是连续谱，而C（jnw）是离散谱；
S（f）单位是幅度/频率，而C（jnw）单位是幅度；（这里都是指其频谱幅度）
能量信号的能量有限，并连续的分布在频率轴上，每个频率点上的信号幅度是无穷小的，只有df上才有确定的非0振幅；功率信号的功率有限，但能量无限，它在无限多的离散频率点上有确定的非0振幅

二、功率谱密度

功率谱：也称功率谱密度（PSD），单位是功率/Hz。针对功率有限信号的(能量有限信号用能量谱密度)，所表现的是单位频带内信号功率随频率的变换情况，即信号功率在频域的分布状况（深度理解）功率谱密度函数反映了某一频域上脉动风的能量分布，是脉动风最重要的统计特性之一（这一部分是需要重点理解的）。

三、计算方法

1、周期图法：它是把随机序列x(n)的N个观测数据视为一能量有限的序列，直接计算x(n)的离散傅立叶变换，得X(k)，然后再取其幅值的平方，并除以N，作为序列x(n)真实功率谱的估计。

2、自相关法：根据维纳-辛钦定理，先估计相关函数，再经傅立叶变换得功率谱估计。功率谱与自相关函数是一个傅氏变换对。功率谱具有单位频率的平均功率量纲，所以标准叫法是功率谱密度。通过功率谱密度函数，可以看出随机信号的能量随着频率的分布情况。像白噪声就是平行于w轴，在w轴上方的一条直线。

四、MATLAB实现功率谱函数计算

1、旧版本中较为简单，直接采用psd函数即可

Fs = 1000; t = 0:1/Fs:.296;
x = cos(2*pi*t*200)+randn(size(t));
h = spectrum.welch; % 创建一个韦尔奇谱估计。
Hpsd = psd(h,x,'Fs',Fs); % 计算PSD 值
plot(Hpsd) % 绘制 PSD.图

2、新版本中主要采用pwelch(来自matlab官方说明文档)

%%每s周期性的计数为1000
fs = 1000;
%% 5000个数据列，总共是5s，每个频率间隔是0.001
t = 0:1/fs:5-1/fs;
x = cos(2*pi*100*t) + randn(size(t));
[pxx,f] = pwelch(x,500,300,500,fs);
plot(f,10*log10(pxx))
xlabel('Frequency (Hz)')
ylabel('PSD (dB/Hz)')
%% x为时程信号，nfft为傅里叶变换(fft)点数
%% windows为窗函数，通常设置为hanning(nfft);
%% nooverlap-重叠的点数，通常为nfft/2；
%% px---------求得的信号x的功率谱密度；
%% fx---------对应px的频率序列。
[px,fx]=pwelch(x,Windows,noverlap,nfft,fs);

[pxx,f] = pwelch(x,window,noverlap,NFFT,fs)

x 是一维的信号数据；

window 是计算功率谱每个窗口的信号长度，关于窗函数的长度选择可以参考公式，选择的窗口越长，越能分辨低频的信号，x_length/fren；谱分析中窗的选取

noverlap 是每个窗口之间重叠的长度，通常取33%~50%。窗口之间重叠得越多，图像越平滑（blurred）；反之则更粗糙（blocky）；

NFFT，即FFT数据点的个数，可以变化。但是最大长度不能超过每一段的点数。当然，通常设置NFFT为大于每一段的点数的最小2次幂，这样可以得到最高的频域分辨率。NFFT越小，最终会越粗糙；

fs是采样频率，最终的结果，横坐标的最大值为采样频率的一半；

pxx 为计算得到的功率谱数值；

f 为功率谱数值对于频率的位置；

pwelch的方法概括步骤如下

将信号分为多段，每段之间可以有overlapping,也可以没有。
每一段加窗
每一段做谱分析
求平均。

以一个小例子做进一步的分析：准备一个余弦函数，周期为1，以10hz采样

n = 0:0.1:20;
x = cos(2*pi*n);
plot(n,x,'.-');
fs=10;

其中不同的参数会对结果有一定的影响

不同的NFFT值决定频谱的最低频率
不同window大小决定谱的粗糙程度
不同fs影响横坐标的位置

现在用pwelch画一下，功率谱

n = 0:0.1:20;
x = cos(2*pi*n);%+randn(size(n));
%plot(n,x,'.-');
fs=10;  %%采样频率，横坐标的最大值为采样频率的一半
NFFT=50;%%FFT数据点的个数
window=50;%%计算功率谱每个信号的长度
noverlap=30;%%每个窗口之间的重叠长度
[pxx,f] = pwelch(x,window,noverlap,NFFT,fs);
plot(f,10*log10(pxx))
xlabel('Frequency (Hz)')
ylabel('PSD (dB/Hz)')

五、频域分析与时域分析的比较：

频域分析：利用傅里叶变换将风荷载变换为一系列的简谐荷载，最后再将结构在每个简谐荷载下的响应叠加起来得到结构的总响应。这种方法的缺点是不能考虑结构的非线性，计算结果不够精确。

时域分析：直接运用风洞试验的风压时程或数值模拟的风压时程作用于计算模型上，通过在时间域内直接求解运动方程得到结构的响应，能够考虑到结构的非线性影响，对于索网、悬索结构这种非线性柔性结构具有一定的优越性。

六、得到风阵系数的步骤

通过风洞试验或数值模拟来确定风荷载时程：没有风洞试验考虑用数值模拟的方法确定风荷载时程
确定结构的质量M，阻尼C，刚度K等参数，建立有限元模型;
将风荷载时程作为外荷载时程作用到有限元模型上;
采用瞬态分析方法得出结构的动力响应，从而得到风振系数。

备注：详细的一些理论基础参见陈博士的风时程生成程序技术说明

Davenport风速谱的代码详解

利用Davenport风速谱进行时程分析，确定本工程柔性支架风振系数，模拟风速谱、风压谱以及模拟结果与Davenport风速谱对比，在参考资料里面有一篇帖子讲述了比较详细的一份分析过程。

利用matlab来进行风速时程的数值模拟

根据风的记录，脉动风可作为高斯平稳过程来考虑，通过观察n个具有零均值的平稳高斯过程.....(这一块的内容是自己没有接触到的部分，也是理论的一个部分，确实比较的难搞定的一块)

算例一模拟出来的结果如下：

%%%%%%%%%%%%%%%%main procedure%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%进行单点模拟%%%%%%%%%%%%%%%%%%%%%%%%%%
N=6000;%%模拟点数6000个
domega=0.001;%%频率增量0.001
omegaup=2*pi;%%频率区间上限
n= domega:domega:omegaup;%%频率区间（0.001～6.28）Hz
v10=16;%%10米高度的风速为16m/s
k=0.005;%%地面粗糙度系数k=0.005
x=1200*n/v10;%%x为1行200pi列的数组,Dvenport谱中间参数
delt=0.1;%%时间增量0.1s

%%采用Davenport风速度谱，s1为脉动风速功率谱
s1=4*k*v10^2*x.^2./n./(1+x.^2).^(4/3);%%由公式进行注写
subplot(2,2,1)%%在本区域显示两行两列个图像，本图显示在第一个位置
loglog(n,s1)%%画谱图，loglog(n,s1) 在 n 轴和 s1 轴上应用对数刻度来绘制 n 和 s1 坐标

%%横纵坐标的标签分别为freq和S
xlabel('freq');
ylabel('S');

for i=1:1:omegaup/domega  %% 这里相当于对列进行相应的循环，有2000pi列，每次按照1来增量  
H(i)=chol(s1(i));%%Cholesky分解，相当于对数组中的每一个数进行开平方处理
end

thta=2*pi*rand(N,1);%%thta为介于0和2pi之间均匀分布的随机数角度，rand(N,1)为介于0和1之间均匀分布的随机数
t=0.1:0.1:600;%%时间区间（0.1～600s），时间增量为0.1
ii=sqrt(-1);%%虚数单位
for j=1:1:N  %%对6000个模拟点进行相应的模拟
v(j)=sqrt(2*domega) *H(j)*exp(ii*thta(j));%%风荷载模拟,这里有三角函数与指数函数的转换的问题
end

Y=fft(v,N);%%对数值解作傅立叶变换
for j=1:1:N
vh(j)=real((Y(j)*exp(ii*j*domega)));%%real函数是求复数的实部
end

%[power1,freq1]= psd(vh,N,2,boxcar(512),0,'mean');

[power,freq]=pwelch(vh,boxcar(3024),10,N,1/delt);

subplot(2,2,2) %%在本区域显示两行两列个图像，本图显示在第二个位置

plot(vh)%%显示风荷载

xlabel('t(s)');

ylabel('v(t)');

axis([0 1800 -10 10]);%%axis([xmin xmax ymin ymax])，设置x和y轴的坐标范围的上下限

subplot(2,2,3) %%在本区域显示两行两列个图像，本图显示在第三个位置

loglog(freq,power,'r',n,s1,'b')%%%%%%拟合谱与目标谱比较%%%%%%%%%%%%%%

xlabel('freq');

ylabel('S');

第二段代码，也是比较类似的

%%%%main procedure%%%%%%%%%%%%%%%%%%%%
clc
clear
close all
N=6000;       %%模拟点数
d=0.001;      %%频率增量
omegaup=6;    %%频率上限
f=d:d:omegaup;  %%频率区间（0.001～6Hz）
v10=20;        %米高度风速
k=0.005;       %%地面粗糙度系数
delt=0.1;       %%时间增量
x=1200*f/v10;  %Dvenport谱中间参数

s1=4*k*v10^2*x.^2./f./(1+x.^2).^(4/3); %Dvenport谱表达式
subplot(2,2,1)
loglog(f,s1)  %%画Davenport经验谱图
xlabel('freq');
ylabel('S');

%%%%进行Cholesky分解%%%%%%%%%%%%%%%

for i=1:1:omegaup/d
H(i)=chol(s1(i));%%Cholesky分解
end

%%%%风荷载模拟%%%%%%%%%%%%%%%%%%%

thta=2*pi*rand(N,1);       %%介于0和6之间均匀分布的随机数
t=1:1:6000;               %%时间区间（0.1～600s）
ii=sqrt(-1);

for j=1:1:N
v(j)=sqrt(2*d)*H(j)*exp(ii*thta(j));%%风荷载模拟
end

%%%%%对风速时程进行FFT变换%%%%%%%%%%%%

Y=fft(v,N);             %%对数值解作傅立叶变换
for i=1:1:N
vh(i)=real(Y(i)*exp(ii*i*d*0.1));
end

[power,freq]=pwelch(vh,boxcar(3024),10,N,1/delt);
subplot(2,2,2)

plot(t/10,vh)                %%显示风荷载
xlabel('t(s)');
ylabel('Y(t)');

%%%%拟合谱与目标谱比较%%%%%%%%%%%%%%
subplot(2,2,3)
loglog(freq,power,'r',f,s1,'b')    %%拟合谱与目标功率谱进行比较
xlabel('freq');
ylabel('S');

未完待续

还有很多的一些东西没有整理好，也会在后面慢慢整理和补充出来......，当然难免会有不足的地方，只能勉励自己后面逐步精进。

参考资料：

《建筑结构荷载规范 GB5009-2012》
https://zhuanlan.zhihu.com/p/59364215
https://blog.csdn.net/deng_sai/article/details/52467619
http://dinochen.com/