基于连续型风速风向联合概率分布的风致结构疲劳分析方法研究

林涛涛，陈 隽1,，李 想，张晓琴

(1.土木工程防灾国家重点实验室，上海 200092； 2.同济大学，上海 200092)

风致结构疲劳是细长高耸钢制结构如天线、电视塔、灯柱、桅杆等重要的设计问题。据统计，既有结构80%以上的破坏由疲劳产生[1]。结构风致疲劳分析方法主要有时域法、频域法和时-频域法。时域法通过结构动力分析首先获得考察点处的应力响应时程，再由循环计数法(如雨流法)来确定损伤值，计算可靠但效率低。频域法基于随机振动理论直接获得关键点应力响应的统计特性，常用的有等效窄带法[2]和上下限法[3]，计算效率较高。时-频域法综合了时域法的精确性和频域法的高计算效率[4]，先通过随机振动理论得到关键点的应力响应功率谱，再用Monte Carlo方法模拟应力过程，最后由雨流计数法获得损伤值。

已有研究表明，风向的分布对结构疲劳寿命的预测有重要影响，忽略风向会低估结构的疲劳寿命。目前一般采用风速风向分离的离散区格方式考虑风向的影响，即先将风向分为若干区间(如常用的16风向区间)，区间内再按风速大小分为若干区格。区格内采用一般算法计算风致疲劳损伤值，并以区格内风速发生概率作为损伤权重系数，再将区格损伤值累加得到总的损伤值。这种处理方式虽然考虑了风向的影响，然而各区格之间的风速函数在交界处不连续，同时区格大小(数量)的划分受风速资料多少的限制。这些都会影响计算结果的准确性。因而，如何在风致结构疲劳寿命计算中更合理准确地考虑风速风向的联合作用是一个值得深入研究问题。

本文利用对风速风向联合概率密度函数(Joint Probability Density Function of wind speed and wind direction，简称JPDF)的研究成果[5]，并结合等效窄带法，将区格离散型的疲劳损伤累加求和公式推广为连续型积分公式，理论上使得风致疲劳寿命计算更为合理和精确。通过将该方法用于一个圆形钢质天线模型的风致疲劳寿命计算,证明了方法的可行性。

1 基于连续型JPDF的风振疲劳频域分析方法

本节通过将JPDF和等效窄带法结合，推导新的风致疲劳寿命计算方法。

1.1 微分区间的损伤值

如前所述，目前一般在风速风向区间间隔(v

图1 微分区间和区间间隔处理方式比较示意图

根据等效窄带法的假定，疲劳验算点处的应力过程S(t)为零均值、窄带平稳高斯过程，则微分区间(v

(1)

式中，S是应力范围；σSv,θ是相应于微分区间(v

NSm=K

(2)

式中，m、K是和材料有关的常数；N是相应于常应力范围S时直到发生疲劳破坏的循环次数。

相应于微分区间(v

(3)

(4)

式中，α0v,θ、α2v,θ分别是相应于微分区间(v

由式(3)和式(4)，可得到相应于风速风向的联合微分区间(v

(5)

则当应力范围区间为(0,+∞)时，窄带过程的损伤值为：

(6)

将式(1)代入式(6)，则相应于微分区间(v

(7)

式中，Γ(*)为Gamma函数。

1.2 宽带修正

对于宽带高斯应力过程采用式(7)时需对损伤预测值进行修正。根据文献[2]，可由实际应力功率谱的谱矩确定一个经验上的宽带修正系数λv,θ，以此来修正窄带过程而得到宽带过程的损伤值：

(8)

式中，λv,θ是相应于微分区间(v

λv,θ=a+(1-a)(1-εv,θ)b

(9)

式中，a,b是与材料有关的参数，其分别由式(10)和式(11)确定；εv,θ是相应微分区间(v

a=0.926-0.033m

(10)

b=1.587m-2323

(11)

(12)

1.3 全部风速风向区间内的累积损伤计算

结合风速风向JPDF，在风速风向微分区间(v

Pv,θ=fV,Θ(v,θ)dvdθ

(13)

式中，fV,Θ(v,θ)为风速风向联合概率密度函数。

则在全风速风向区间(0

(14)

令D=1，由式(14)可得连续积分形式的疲劳寿命计算公式为：

(15)

在文献[7]给出了公式(1)-(15)的全部推导过程和细节。

1.4 与离散区格公式的区别及联系

以上推导时若采用风速风向区间间隔(v

(16)

对于式(16)，令D=1，可得求和累加形式的疲劳寿命计算公式为：

(17)

比较式(14)和(16)，式(15)和(17)可知，由于引入了JPDF函数，损伤计算的区格累加形式变成了整个风速风向空间上的连续积分形式，使得疲劳损伤和疲劳寿命的计算在理论上更加合理。双重连续积分形式在数值计算中可以任意设置数值离散积分间隔，可以达到满意的精度；而双重累加形式和风速风向频度统计有关，要想达到满意的精度，又需要重新统计风速风向频度资料。此外，区隔累加形式的风速上限值为实测统计资料中的最大风速，而JPDF可选用任何合理的形式，从而包括更大风速情况的影响。以上这些特点都可提高疲劳损伤计算的精度。

2 JPDF的选取

文献[5]对现有的风速风向的联合分布函数模型以及它们的适用性进行了详细地讨论。主要模型包括：经验模型、区间模型、双高斯模型、角-线性模型、Farilie-Gumbel-Morgensternmoxign (FGM)模型和Simple C-Type Distribution (SC模型)等。本文采用文献[5]推荐的FGM模型。

2.1 FGM模型

直接引入风速随机变量V和风向随机变量Θ，根据文献[9]的研究，FGM模型的通用表达式为：

FV,Θ(v,θ)=FV(v)×

FΘ(θ)[1+δ(1-FV(v))(1-FΘ(θ))]

(18)

式中，v、θ是风速值和角度值；FV,Θ(v,θ)是风速风向联合累积分布函数(JCDF)；FV(v)是风速累积分布函数(CDF)；FΘ(θ)是风向CDF；δ是表示风速、风向之间联合作用大小的相关性系数，其值在-1≤δ≤1之间，|δ|越小说明风速、风向之间的相关性越小，反之则越大。

式(18)两边分别对v和θ求偏导数，经过化简整理，得到风速风向联合概率密度函数(JPDF)的表达式为：

fΘ(θ){1+δ[1-2FV(v))(1-2FΘ(θ)]}

(19)

式中，fV,Θ(v,θ)是风速风向JPDF；fV(v)是风速概率密度函数(PDF)；fΘ(θ)是风向PDF。

2.2 风速PDF模型和风向PDF模型

上述计算中需要风速和风向的概率密度函数模型。其中风速PDF模型取为正态-威布尔分布(Normal-Weibull distribution，NW模型)为：

φ1,φ2)+

-∞≤v≤∞

(20)

(21)

φ1,φ2)dv

(22)

式中，ω0为Normal分布模型所占的权重；α、β分别为Weibull分布模型中的形状参数和尺度参数；φ1、φ2是和随机变量单位一致的参数。

风向PDF取混合von Mises模型[10]：

0≤θ≤2π

(23)

(24)

I0(κj)=

(25)

式中，fΘ(θ)是风向PDF；N是混合von Mises模型的组成个数，本文取N=5；ωj为第j个von Mises模型所占的权重，满足式(24)；κj和μj是两个参数，分别满足κj≥0、0≤μj<2π；I0(κj)是0阶修正的第1类Bessel函数。

3 实例分析

3.1 实例概况

分析实例采用文献[11]D类地区某超高层建筑上的钢质天线结构，该超高层结构的顶部标高(即天线底部标高)为246 m，天线结构总高为87 m，结构分为三段，从下往上三段天线的长度分别为34 m、27 m、26 m，三段天线对应的截面外径分别为2.4 m、1.75 m、1 m，三段天线对应的截面壁厚分别为3 cm、2.5 cm、2 cm(见图2)。天线所用钢材的质量密度为7.85×103kg/m3；结构阻尼采用Rayleigh阻尼，阻尼矩阵由质量矩阵和刚度矩阵决定，即[C]=a[M]+b[K]，a,b为常数，对于钢质天线结构，取前4阶的阻尼比为0.01。结构动力特性分析采用集中质量法。在仅考虑顺风向荷载作用下，圆形截面结构的风致疲劳寿命较低的关键点(即最危险关键点)往往出现在风向发生频度较高的位置(关于圆心对称的关键点的疲劳寿命也较低)，则根据风速风向发生概率，所取关键点位置与X轴夹角分别为90°、120°、130°(见图2)。顺风向脉动风速谱采用Kaimal谱，顺风向空间相关函数采用Shiotani模型。

天线结构的振动有X和Y两个方向，且质量和刚度关于X和Y两个方向都是对称的，这大大简化了整个分析和计算过程。在此规定，奇数阶振型和频率是X方向的，偶数阶振型和频率是Y方向的。前6阶频率见表1。

计算步骤：①整理风速风向实测数据；②使用Mathcad软件拟合风速PDF和风向PDF模型参数，获得模型表达式；③使用Matlab软件拟合JPDF模型中的相关性系数，获得JPDF表达式；④分析实例的动力特性；⑤选择控制截面及控制点；⑥编程计算疲劳寿命。

图2 钢天线模型

表1 前6阶频率表

JPDF模型参数的计算数据为上海龙华气象站1956年1月1日至1990年12月31日全部12784天10米高度处的日最大风速值(10分钟平均时距)和相应的风向记录的统计结果。为了避免不同样本值之间的相关性，按照阶段极值抽样方法，抽取8天时间间隔的极值风速组成极值子样[12]。

采用Mathcad软件进行风速PDF和风向PDF各参数的拟合，采用Matlab软件进行相关性系数δ的拟合。需要注意的是参数拟合属于最优化问题，当拟合初值不同时，参数拟合结果会有所不同，但拟合优度R2接近。

风速PDF和风向PDF各参数的拟合结果分别见表2和表3。FGM模型中的相关性系数δ的拟合结果为0.0993125。拟合优度达到0.970350，可见FGM模型的精度是相当令人满意的。图3和图4分别是实测的JPDF和拟合的JPDF,比较看出，拟合效果相当好。

表2 风速PDF模型参数拟合值

表3 风向PDF模型参数拟合值

图3 实测的JPDF

图4 拟合的JPDF

图5 第二变截面关键点2处不同风速风向角下的应力功率谱

3.3 疲劳寿命估算

在数值计算时，风向角度间隔取为15°，10 m高度处的风速间隔取1 m/s；风速的上限值取10 m高度处对应的20 m/s，此上限值已经足够满足数值积分的精度要求；根据《钢结构设计规范 GB 50017-2003》，假定构件类别为第8类，则参数K取值为4.1×1011，参数m取值为3。

根据风速风向JPDF的分布形式，考虑风向发生概率等综合因素，所选取的疲劳分析关键点位置为天线结构最外边缘，且与X轴方向夹角分别为90°、120°、130°。计算截面分别为底部截面、第一变截面、第二变截面，如图2所示。

图5给出了天线结构底部截面和第二变截面关键点2处在不同风速风向角下的应力功率谱图，风速是10 m高度处的风速，从图中可以看出随着风速的增大，应力功率谱值是增大的；在结构对应的自振频率处，其功率谱值达到峰值。

根据连续积分形式的疲劳寿命计算公式(15)，计算得到底部截面、第一变截面和第二变截面各关键点的疲劳寿命，见表4。最终认为频域计算方法估算的疲劳寿命为：173.89年。

表4 关键点疲劳寿命(单位：年)

4 结 论

本文结合连续型的风速风向联合概率密度函数(JPDF)，推导了基于等效窄带法的连续积分形式的疲劳寿命计算公式，使得疲劳寿命的估算公式在理论上变得更加合理准确。应用于实例工程的疲劳寿命预测表明了连续积分形式计算公式的可行性。为更加精确的疲劳寿命估算提出了一条新思路。

[1] Repetto M P.Wind-induced response and fatigue of slender vertical structures [D].Universita Degli Studi Di Genova，2003.

[2] Wirsching P H，Light M C.Fatigue under wide band random stress[J].Journal of Structural Division.1980，106(7): 1593-1607.

[3] Holmes J D.Fatigue life under along-wind loading-closed-form solutions[J].Engineering Structures.2002，24(1): 109-114.

[4] Gu M，Xu Y L，Chen L Z，et al.Fatigue life estimation of steel girder of Yangpu cable-stayed Bridge due to buffeting[J].Journal of Wind Engineering and Industrial Aerodynamics.1999，80: 383-400.

[5] 张晓琴.风速风向联合分布模型与应用研究[D].上海：同济大学，2011.

[6] 俞载道，曹国敖.随机振动理论及其应用 [M].上海: 同济大学出版社，1988.

[7] 林涛涛.基于风速风向联合概率密度函数的结构风致疲劳计算方法研究 [D].上海：同济大学，2013.

[8] 王钦华.复杂结构风致疲劳的理论和应用研究[D].上海：同济大学,2009.

[9] Erdem E，Shi J.Comparison of bivariate distribution construction approaches for analysing wind speed and direction data[J].Wind Energy，2011，14 (1): 27-41.

[10] Carta J A，Bueno C，Ramirez P.Statistical modelling of directional wind speeds using mixtures of von Mises distributions: Case study[J].Energy Conversion and Management，2008 (49): 897-907.

[11] 李想.基于风速风向联合分布函数的结构风致疲劳寿命计算方法研究[D].上海:同济大学,2011.

[12] 杨咏昕，葛耀君，项海帆.风速风向联合分布的平均风统计分析[J].结构工程师，2002 (3): 29-36.

YANG Yong-xin，GE Yao-jun，XIANG Hai-fan.Statistic analysis of wind speed based on joint distribution of wind speed and wind direction [J].Structural Engineers，2002(3):29-36.