杨素珍

(漳州职业技术学院 机械工程学院，福建 漳州363000)

台风往往给人类的生产和生活环境带来巨大破坏[1]，造成大量的经济损失和人员伤亡，通过对人们活动环境内重要目标，如建筑、桥梁、塔架和电线杆等，进行台风风载荷作用测试，研究物体风致破坏的机理，从而找出减小台风破坏的措施。常规基于风洞试验和现场实测的风载荷作用检验手段成本高、耗时长；通过模拟台风的方法仿真实施目标风载荷作用[2-3]，可以较好的解决上述问题。因此研究台风模拟手段意义重大，引起国内外学者的广泛关注[4-9]。

风的变化过程可以近似为一种平稳随机过程，一般可将其分解为稳态平均风和动态脉动风。平均风的风速值和风向不随时间变化；脉动风则随时间无规律地改变其风速和风向。脉动风动力成分是强(台)风的重要组成[5]，集中反映其随机特性，该特性可以用紊流强度和功率谱密度函数等来描述。文献[4]借助不同坡度三维山体模型风洞试验，研究了山体各位置脉动风速谱的变化规律。文献[5]利用台风“罗莎”实测数据，研究其高空风参数特性，结合Kolmogrove理论，提出水平方向和垂直方向脉动风功率谱密度函数。上述基于实际风速的方法成本高、耗时长。文献[6] 基于谐波合成法进行海上风机的脉动风速度时程模拟，考虑了过程顺风方向和垂直方向计算点的二维相干性。文献[7]基于标准达文普(Davenport)模拟脉动风的风速和风压时程曲线，并将模拟风谱应用于巨型框架的减振结构设计。文献[8]提出一种结合谐波合成法和本征正交分解的脉动风场数值模拟方法，来解决计算机内存超限问题。文献[9]采用线性滤波法模拟脉动风速，生成作用于上海中心大厦结构上的脉动风载荷时程，为该大厦的风结构设计奠定基础。上述基于脉动风风速时程模拟方法所用的谱密度函数是通过大量实测数据的统计分析所获得的经验式，但在表示特定地区风时，由于各地区风移动环境存在较大差异，经验式不能完整反映该地区的实际风特性。

为此，本文提出一种基于特定地区实际风采集数据的风模拟方法，该方法利用实测台风风速数据和最小二乘法对标准Davenport功率谱密度函数进行修正，再基于一类改进型谐波合成法进行风模拟，在较完整反映特定地区风特性的同时，保证模拟的快速性和精确性。本文首先分析福建某地区实测台风特性，然后详细阐述修正Davenport模型数学建模和基于改进型谐波合成法的多点脉动风模拟算法，最后给出实验情况。

1 特定地区实测风特性分析

为了获取特定地区实际台风的特性信息，通过专用风采集装置采集特定高度实际风速样本，以台风“杜鹃”为例，采集点为福建某城市，取得实际风速样本10000个，采样频率为0.1s，时间总长度为1000s，其反映能量随频率分布的频域峰值谱、功率谱分别如图1和图2所示。由图1可以看出，台风“杜鹃”的风能主要集中在0～3Hz的范围内，超过该范围基本没有风能的分布，实际风幅值的集中范围情况与标准Davenport十分相似。

图1 实测风频域幅值谱(截止频率4Hz)

图1对比结果表明，通过采集装置采集到的风速数据是真实可信的。

图2 实测风频域功率谱

从图2的功率谱上看，二者明显偏离过大，吻合度不够，可见，如果采用标准Davenport进行风模拟，无法完全体现采集所在地区风的特点，需要结合实测风对标准Davenport进行修正，才能更好地反映该地区风的真实特性。

2 实测风修正Davenport模型建模

平稳随机过程可通过自功率谱密度函数来表征其能量在不同频率上的分布，风模拟的关键就是选取能够反映风时程数值特征的风功率谱密度函数。基于Kolmogrov瑞流理论，风工程领域学者通过对自然风的大量观测和实测数据统计，已建立了许多满足瑞流模型的风功率谱密度函数，常用的有Kaimal谱、Simitu谱、ESDU谱和Davenport谱。由于观测地区和数据统计方式的不一致性，不同谱密度函数所模拟的风场频率结构各有差异，一般无法与特定地区的实际风特性完全匹配。

瑞流风速主要包括水平、纵向和竖向三个方向风速，其中水平方向产生的风压对结构的作用最直接，破坏性也最强，目前包括美国、德国、加拿大和中国在内的大多数国家的风载荷规范中，普遍采用水平脉动风谱Davenport谱[10]，该谱不随高度变化，频谱中高频谐波分量明显，参数配置简单，方便进行工程应用。Davenport谱的水平顺向风模型为

(1)

式中：x=1200f/U10；f为频率(Hz)；U10为10m高度处的平均风速；k为地面粗糙度系数(VonKaman常数)。

为使Davenport模型能够更真实反映特定地区的实际风特性，本文提出一种修正Davenport模型，利用特定地区实测数据确定该表达式的待定系数，让修正模型中包含特定地区的客观变化和行为信息。定义修正Davenport模型的表达式为

(2)

式中P1、P2、P3、P4为待定系数。

(3)

式中：F(f)为频域幅值谱，E(·)表示均值，T为风速样本持续的时间。

文中利用实测台风“杜鹃”风速数据，求取风速功率谱，得到求解修正Davenport模型系数所需的频率值和功率值样本，然后基于一类非线性最小二乘法拟合求解得到修正Davenport模型中的待定系数。

定义待定系数初值为P0=[P10，P20，P30，P40]，将该初值代入式(2)得到ψ0，则待定系统的求解步骤如下。

(1)对函数(2)在P0点进行泰勒级数展开，并忽略高阶项目，得到近似一阶表达式

(4)

式中： ϑi=Pi-Pi0，i=1,2,3,4，j=1,2,…,n表示用于最小二乘法拟合数据点个数序号；ψj0为序号j数据点初值函数值。

(2)建立带求解系数非线性函数的近似线性关系表达式

(5)

(3)确定求解bi的最小二乘求解准则为：使n个数据点(xj,yj)与非线性函数y距离ηj的平方和最小，即求取bi使得

(6)

为最小值。

(4)对方程(6)求偏导数，并令其等于零，可得

(7)

再将n(n>4)个数据点代入方程(7)，可建立超定方程组

Φ(X,Y,B)=0

(8)

式中：Y=(y1,y2,…,yn)；B=(b1,b2,b3,b4)；

X=(x11,…,xn1,x12,…,xn2,x13,…,xn3,x14,…,xn4)求解超定方程组中的B值即为最小二乘解。

(5)求得待定系数为

P=[P1,P2,P3,P4]=B+P0

(9)

对非线性函数(2)的近似线性形式(5)，可以采用多元线性拟合方法进行求解，降低了计算复杂度。由于采用近似线性形式，初步计算的值存在一定误差，为提高拟合精度，需将求解后的值重新作为初值，多次进行上述求解步骤，直至满足所需精度要求。

3 改进型多点脉动风谐波合成法

离散傅里叶变换可以用来描述具有时间和空间特性的风速时间序列，谐波合成法(谱表示法)的原理是采用一系列具有随机频率的余弦函数序列来模拟脉动风，并考虑结构风速的特点。由于同一时刻作用于结构上的多个点的风载荷不可能完全一样，因此，在模拟风速时程时，应考虑其相关性，即采用多个互相关的随机过程来描述。

风的实测数据表明，风速的脉动部分可看作是具有零均值的高斯(Gauss)平稳随机过程，对具有不同坐标的风点(1，2，3，…，n)，其风速uj(t)(j=1，2，3，…，n)的谱密度函数矩阵为

(10)

式中：S11(ω)，S22(ω)，…，Snn(ω)为不同坐标点风速功率谱；Sij(ω)(i≠j)为i点风和j点风的相关功率谱。风速功率谱采用式(2)的修正Davenport，以让模拟风包含指定地区的风特性，相关功率谱采用式(11)表达式。

(11)

式中：i≠j，Zi、Zj分别为i点风和j点风的高度值，Vi、Vj分别为i点风和j点风的平均风值。

非负定矩阵[S(ω)]可以用Cholesky方法分解成如下形式

[S(ω)]=[H(ω)][H*(ω)]T

(12)

式中

(13)

此时，随机过程uj(t)可以用式[14]进行模拟(当N趋于无穷大时)

θjm(ωml)+φml]

(14)

式中：N为频率域内的数据采集数目；

φ1l、φ2l、…、φml为相互独立的[0，2π]上均匀分布的随机相位角序列；Δω=ωu/N，ωu为截断频率，其大小通常由功率谱密度矩阵中各项与ω的函数关系而定。为避免频率混淆，模拟的时间不长，Δt须满足式(15)。

Δt≤2π/(2ωu)

(15)

采用谐波合成法脉动风模拟尽管计算精度高，但需要在每个频率上进行大量运算，极大地降低了算法的计算速度。本文通过引入快速傅里叶变换FFT[11]改进常规谐波合成法来简化计算量，从而提高算法的计算效率。改进后式(14)的表达式为

(16)

式中：η=1,2,3,…,2N×n-1；j=1,2,3,…,n；λ=1,2,3,…,n-1，为η/2N余数；Qjm(λΔt)的表达式为

(17)

式中

(18)

从式(18)可以知道Qjm(λΔt)是Djm(sΔω)的傅里叶变换，因此，可以用FFT进行相关计算。

4 实验研究

对于修正Davenport模型建模，本文以某时段台风“杜鹃”的风速为样本，计算得到相应的频率值和功率值作为最小二乘法求解的数据点，选取待定系数初值P0=[1，1.5，1.5，2]，利用式(4)～式(9)，经过多次循环求解，得到系数值为

P1=0.4696，P2=1.2721，P3=0.1023，P4=1.8459。

系数代入式(2)可得到修正Davenport模型。

将所建立的修正前后Davenport模型与实测风的功率谱进行对比，结果如图3所示。

图3 台风杜鹃与修正前后Davenport模型功率谱对比

由图3可知，模型修正前，台风杜鹃的功率谱幅值显然比标准Davenport大，二者的吻合度差，且频率越大，标准Davenport偏离台风杜鹃功率谱的趋势加大；而采用修正后的Davenport模型，二者重合度得到了较大的提高。

进一步以某时段台风“苏迪罗”的风速为样本，修正Davenport模型采用与台风“杜鹃”实验相同的系数值，实测风与修正前后Davenport模型的功率谱，对比结果如图4所示。

由图4可以看出，修正后的Davenport模型与修正前相比，同样具有更好的吻合度，因此，本文基于修正Davenport模型进行风速模拟，可以更好地体现实际风速样本采集地的风特性。值得一提的是，通过增加风速样本数量，提高最小二乘法的拟合精度，可以进一步提高吻合度。

图4 台风苏迪罗与修正前后Davenport模型功率谱对比

为了说明本文修正Davenport模型对台风的模拟效果，这里以12级台风为例，进行台风时间序列模拟。地面粗糙度取B类，查风速等级表得到Z1=10m高度的平均风速为V1=30.9m/s，取风速点数为n=2，第二点高度Z2=7m，根据指数风剖面平均风速沿高度变化的规律，可求取相关第二点的平均风速值为

(19)

取风速时程点数N=2048，截断频率ωu=10π，利用式(13)和式(16)，风模拟结果如图5和图6所示。

从图5可以看出，台风随时间变化过程具有明显的脉动风特性，且风速变化跨度大，符合实际测量强台风的变化特点。从图6可以看出，模拟台风的功率谱与修正Davenport目标谱吻合度高，说明风速时间模拟序列较好地体现了指定地区实际台风的风特性。

常规与改进型多点脉动风谐波合成法的上述台风时间序列模拟过程耗时对比如表1所示。

图5 模拟台风风速时程谱

图6 模拟台风功率谱与目标谱比较

表1 常规与改进型多点脉动风谐波合成法耗时对比s

由表1可见，基于快速傅里叶变换FFT的改进型谐波合成法有效减小了算法的计算量，整个模拟过程耗时仅为常规谐波合成法的1/8，算法运行实时性得到了较大的改善。

5 结束语

台风模拟是一种快速有效的分析台风环境内物体风致时程响应的手段，通过实测风数据分析表明，采用常规功率密度函数进行脉动风模拟的方法，不能完全反映出特定地区的风特性，本文通过利用福建某城市实测风数据和最小二乘法对标准Davenport模型进行再建模，获得更能表征该地区风特性的修正Davenport模型，在此基础上，采用多点脉动风谐波合成法进行台风模拟，并通过FFT方法改进其计算效率，提高了与目标谱的吻合度。基于本文所提方法，通过其他地区的实测风数据，就可以实现能反映该地区风特性的台风模拟，为准确研究该地区风致破坏奠定良好基础。