工程科学学报，第39卷，第9期：1453-1462,2017年9月 Chinese Journal of Engineering,Vol.39,No.9:1453-1462,September 2017 D0l:10.13374/j.issn2095-9389.2017.09.020:htp:/journals..usth.edu.cn 全风向来流非高斯风场风机疲劳寿命可靠性分析 双 妙12)，宋波12) 1)北京科技大学土木与资源工程学院，北京1000832)北京科技大学强震区轨道交通工程抗震研究北京市国际科技合作基地，北京 100083 ☒通信作者，E-mail:b20140020@x.usth.edu.cm 摘要在Hermite矩模型基础上，根据Kaimal谱生成某典型风机结构正常风速条件下，三种不同概率特性风场（高斯、非高 斯硬化和软化)，在考虑来流风向和平均风速联合概率密度条件下，以塔架基础连接处为例，对风机进行疲劳寿命可靠性分 析.由叶片的气动模型和多体动力，计算出风机的动力响应，并对响应的时域和频域特性进行分析.基于线性损伤累积理论 和Pis公式，对来流全风向条件下的裂纹形成寿命和裂纹扩展寿命进行了详细讨论.结果表明，裂纹形成寿命对风荷载的非 高斯性较为敏感，而裂纹扩展寿命对风荷载的非高斯性并不敏感，需要考虑风荷载的非高斯性对风机结构疲劳损伤的影响 此外，在考虑全风向来流条件下，疲劳裂纹形成和扩展阶段的失效位置相同，均在主导风向上 关键词全风向；非高斯性：雨流计数法：裂纹形成；裂纹扩展 分类号TU311.3 Reliability analysis of the fatigue life of wind turbines under a non-Gaussian wind field with a full-direction inflow SHUANG Miao)SONG Bo'2) 1)School of Civil and Resource Engineering,University of Science and Technology Beijing,Beijing 100083,China 2)Beijing International Cooperation Base for Science and Technology-Aseismic Research of the Rail Transit Engineering in the Strong Motion Area,Uni- versity of Science and Technology Beijing,Beijing 100083,China Corresponding author,E-mail:b20140020@xs.ustb.edu.en ABSTRACT Using the Hermite moment model,three types of wind fields with different probability characteristics,namely the Gaussian,hardening non-Gaussian,and softening non-Gaussian processes,were generated via the Kaimal spectrum for a typical wind turbine under operational conditions.Reliability analysis of the fatigue life of the wind turbine was performed by taking tower base con- nections as an example with consideration of the joint probability density distribution of the wind direction and mean wind speed.The dynamic response was calculated by an aerodynamic model of the blade and multi-body dynamics,and the time-and frequency-domain characteristics of the response were analyzed.Using linear damage accumulation theory and the Paris equation,the fatigue crack initia- tion life and crack growth life of the turbine were discussed in detail.Fatigue estimation shows that the crack initiation life of the tur- bine is more sensitive to non-Gaussian winds,whereas its crack propagation life is less sensitive to the non-Gaussian characteristics of wind load.The influence of the non-Gaussian characteristics of wind load on the fatigue damage of the wind turbine should be consid- ered.In terms of the full-direction inflow,the failure positions of the crack initiation and propagation stages are identical and in the dominant wind direction. KEY WORDS full-direction inflow;non-Gaussian characteristic;rain-flow counting method;crack initiation;crack propagation 收稿日期：2017-01-11 基金项目：国家自然科学基金资助项目(51178045)：北京科技大学与台北科技大学专题联合研究计划资助项目(TW201601):高校2016年度 引智项目资助项目(110000201420160126)
工程科学学报,第 39 卷,第 9 期:1453鄄鄄1462,2017 年 9 月 Chinese Journal of Engineering, Vol. 39, No. 9: 1453鄄鄄1462, September 2017 DOI: 10. 13374 / j. issn2095鄄鄄9389. 2017. 09. 020; http: / / journals. ustb. edu. cn 全风向来流非高斯风场风机疲劳寿命可靠性分析 双 妙1,2) 苣 , 宋 波1,2) 1) 北京科技大学土木与资源工程学院, 北京 100083 2) 北京科技大学强震区轨道交通工程抗震研究北京市国际科技合作基地, 北京 100083 苣通信作者, E鄄mail: b20140020@ xs. ustb. edu. cn 收稿日期: 2017鄄鄄01鄄鄄11 基金项目: 国家自然科学基金资助项目(51178045);北京科技大学与台北科技大学专题联合研究计划资助项目(TW201601);高校 2016 年度 引智项目资助项目(110000201420160126) 摘 要 在 Hermite 矩模型基础上,根据 Kaimal 谱生成某典型风机结构正常风速条件下,三种不同概率特性风场(高斯、非高 斯硬化和软化),在考虑来流风向和平均风速联合概率密度条件下,以塔架基础连接处为例,对风机进行疲劳寿命可靠性分 析. 由叶片的气动模型和多体动力,计算出风机的动力响应,并对响应的时域和频域特性进行分析. 基于线性损伤累积理论 和 Paris 公式,对来流全风向条件下的裂纹形成寿命和裂纹扩展寿命进行了详细讨论. 结果表明,裂纹形成寿命对风荷载的非 高斯性较为敏感,而裂纹扩展寿命对风荷载的非高斯性并不敏感,需要考虑风荷载的非高斯性对风机结构疲劳损伤的影响. 此外,在考虑全风向来流条件下,疲劳裂纹形成和扩展阶段的失效位置相同,均在主导风向上. 关键词 全风向; 非高斯性; 雨流计数法; 裂纹形成; 裂纹扩展 分类号 TU311郾 3 Reliability analysis of the fatigue life of wind turbines under a non鄄Gaussian wind field with a full鄄direction inflow SHUANG Miao 1,2) 苣 , SONG Bo 1,2) 1) School of Civil and Resource Engineering, University of Science and Technology Beijing, Beijing 100083, China 2) Beijing International Cooperation Base for Science and Technology鄄Aseismic Research of the Rail Transit Engineering in the Strong Motion Area, Uni鄄 versity of Science and Technology Beijing, Beijing 100083, China 苣Corresponding author, E鄄mail: b20140020@ xs. ustb. edu. cn ABSTRACT Using the Hermite moment model, three types of wind fields with different probability characteristics, namely the Gaussian, hardening non鄄Gaussian, and softening non鄄Gaussian processes, were generated via the Kaimal spectrum for a typical wind turbine under operational conditions. Reliability analysis of the fatigue life of the wind turbine was performed by taking tower base con鄄 nections as an example with consideration of the joint probability density distribution of the wind direction and mean wind speed. The dynamic response was calculated by an aerodynamic model of the blade and multi鄄body dynamics, and the time鄄 and frequency鄄domain characteristics of the response were analyzed. Using linear damage accumulation theory and the Paris equation, the fatigue crack initia鄄 tion life and crack growth life of the turbine were discussed in detail. Fatigue estimation shows that the crack initiation life of the tur鄄 bine is more sensitive to non鄄Gaussian winds, whereas its crack propagation life is less sensitive to the non鄄Gaussian characteristics of wind load. The influence of the non鄄Gaussian characteristics of wind load on the fatigue damage of the wind turbine should be consid鄄 ered. In terms of the full鄄direction inflow, the failure positions of the crack initiation and propagation stages are identical and in the dominant wind direction. KEY WORDS full鄄direction inflow; non鄄Gaussian characteristic; rain鄄flow counting method; crack initiation; crack propagation
·1454· 工程科学学报，第39卷，第9期 风机结构自振频率与高度在60~160m之间的高 轮教高度90m叶片直径126m的5MW三叶片陆上风 层建筑相似，并且风机在设计中不考虑舒适度，因此风 机为例)].该风机采用变速变桨控制系统，额定风速 振响应比一般结构更为显著.同时，由于叶片与来流 为11.4ms1,额定转速12.1rmin.对结构在平稳 之间相互作用且机舱结构接近于刚体质点，动力响应 高斯、硬化和软化非高斯3种不同概率特性风荷载作 按传统有限元计算会产生较大误差.常用方法可分为 用下的动力响应和疲劳可靠性进行分析 多体动力-)和基于叶片与来流的气动模型两大类. 算例采用FAST(fatigue,aerodynamics,structures 目前的风机设计规范中，风场通常简化为平稳高 and turbulence)建模，风机叶片划分为I7个单元，塔筒 斯过程.对于平坦均匀地形的边界层风场，高斯假设 划分为10个单元（如图1），塔底与地面刚接.叶片和 是恰当的.然而，实测数据表明，复杂地形的脉动风往 塔体的材料参数和气动参数可参见文献[15]. 往显示出明显的非高斯特性.同时，非高斯荷载作用 10 20 下的建筑结构和海岸工程具有较大的极值响应和疲劳 损伤.因此需要根据风机结构的动力和控制特性，对 280 310 非高斯荷载作用下结构的广义动力响应和疲劳可靠性 截面2 进行分析.此外，由于来流风向随时间变化，疲劳估计 590 520 中需根据工程所处地区气象统计资料，考虑全风向来 流对疲劳损伤累积的影响. 对于风机等细长结构的风致疲劳问题，Do等[5-] 931 940 95 面1 964 根据Monte Carlo模拟和多体动力分析了塔架基础连 3.87m 09 接处的疲劳寿命，并提出基于性能的设计方法.Da- 270 wood等[)在疲劳损伤累积准则和雨流计数法基础上 5 mm /6=25mm 考虑材料的不确定性，对典型细长高耸结构的风致疲 截面1 截面2 180P 截面1 劳进行研究.Repetto与Solari[s-]根据动力响应的背 图1风机结构和风场示意图 景和共振两部分，采用双峰计数法[]研究了高斯荷载 Fig.1 Sketch of the wind turbine and wind field 作用下高耸结构的风致疲劳. 现有的风荷载非高斯性对风机结构影响的研究多 1.2平均风速和来流风向 集中于极值响应分析-]，对疲劳损伤的研究相对较 风机结构任意高度处的风荷载均可分解为平均和 少.由破坏实例可知，风机失效多为叶片根部和塔架 脉动两部分.对于平均风荷载，可根据指数律由参考 基础连接处.目前，复合材料的广泛使用提高了叶 高度处平均风速换算得到.平均风速()的概率密度 片强度和疲劳寿命，因此本文将控制截面限定于塔架 函数可按Weibull分布简化，即： 基础连接处.结合结构和风荷载特性可知，塔体底部 ee (1) 截面由面内、面外弯曲应力和自重引起的压应力控制， 式中，c和k分别是尺度和形状参数，可根据工程所处 处于单轴应力状态.同时，控制截面不同位置的弯曲 地的风速统计资料按最大似然估计法(MLM)、统计矩 应力与来流风向相关，因此在可靠性分析中需考虑来 法(MM)和能量密度法(PDM)计算得到.考虑到方法 流风向的影响. 的简便性及结果的准确性，本文按能量密度法对两参 基于以上分析，本文在Hermite矩模型基础上，通 数Weibull分布进行估计，即uo]: 过非高斯穿越过程模拟了多维相关平稳高斯、非高斯 3.69 硬化和软化三种不同概率特性的风场，详细研究了三 k=1+ 种风场作用下风机结构广义动力响应的时域和频域特 (Ei (2) 性.结合气象资料，根据来流风向将控制截面划分成 不同控制点，采用线性疲劳损伤累积理论，对全风向来 r小+) 流不同概率特性风场作用下不同控制点的疲劳裂纹形 式中，T(·)是Gamma函数；E是能量系数，定义为： 成寿命进行讨论：并采用Pais公式，计算出疲劳裂纹 扩展寿命 E=7 (3) 1风机与来流特性 如图2所示为NREL135m气象塔88m高度处 超声风速仪测得的2014年全年10min平均风速（即轮 1.1风机结构 毂高度处的平均风速山)的玫瑰图.图3所示为3种 本文以NREL(national renewable energy laboratory) 不同方法Weibull分布参数估计与实测数据的比较
工程科学学报,第 39 卷,第 9 期 风机结构自振频率与高度在 60 ~ 160 m 之间的高 层建筑相似,并且风机在设计中不考虑舒适度,因此风 振响应比一般结构更为显著. 同时,由于叶片与来流 之间相互作用且机舱结构接近于刚体质点,动力响应 按传统有限元计算会产生较大误差. 常用方法可分为 多体动力[1鄄鄄2]和基于叶片与来流的气动模型[3鄄鄄4]两大类. 目前的风机设计规范中,风场通常简化为平稳高 斯过程. 对于平坦均匀地形的边界层风场,高斯假设 是恰当的. 然而,实测数据表明,复杂地形的脉动风往 往显示出明显的非高斯特性. 同时,非高斯荷载作用 下的建筑结构和海岸工程具有较大的极值响应和疲劳 损伤. 因此需要根据风机结构的动力和控制特性,对 非高斯荷载作用下结构的广义动力响应和疲劳可靠性 进行分析. 此外,由于来流风向随时间变化,疲劳估计 中需根据工程所处地区气象统计资料,考虑全风向来 流对疲劳损伤累积的影响. 对于风机等细长结构的风致疲劳问题,Do 等[5鄄鄄6] 根据 Monte Carlo 模拟和多体动力分析了塔架基础连 接处的疲劳寿命,并提出基于性能的设计方法. Da鄄 wood 等[7]在疲劳损伤累积准则和雨流计数法基础上 考虑材料的不确定性,对典型细长高耸结构的风致疲 劳进行研究. Repetto 与 Solari [8鄄鄄9] 根据动力响应的背 景和共振两部分,采用双峰计数法[10]研究了高斯荷载 作用下高耸结构的风致疲劳. 现有的风荷载非高斯性对风机结构影响的研究多 集中于极值响应分析[11鄄鄄13] ,对疲劳损伤的研究相对较 少. 由破坏实例可知,风机失效多为叶片根部和塔架 基础连接处[14] . 目前,复合材料的广泛使用提高了叶 片强度和疲劳寿命,因此本文将控制截面限定于塔架 基础连接处. 结合结构和风荷载特性可知,塔体底部 截面由面内、面外弯曲应力和自重引起的压应力控制, 处于单轴应力状态. 同时,控制截面不同位置的弯曲 应力与来流风向相关,因此在可靠性分析中需考虑来 流风向的影响. 基于以上分析,本文在 Hermite 矩模型基础上,通 过非高斯穿越过程模拟了多维相关平稳高斯、非高斯 硬化和软化三种不同概率特性的风场,详细研究了三 种风场作用下风机结构广义动力响应的时域和频域特 性. 结合气象资料,根据来流风向将控制截面划分成 不同控制点,采用线性疲劳损伤累积理论,对全风向来 流不同概率特性风场作用下不同控制点的疲劳裂纹形 成寿命进行讨论;并采用 Paris 公式,计算出疲劳裂纹 扩展寿命. 1 风机与来流特性 1郾 1 风机结构 本文以 NREL(national renewable energy laboratory) 轮毂高度 90 m 叶片直径 126 m 的 5 MW 三叶片陆上风 机为例[15] . 该风机采用变速变桨控制系统,额定风速 为 11郾 4 m·s - 1 ,额定转速 12郾 1 r·min - 1 . 对结构在平稳 高斯、硬化和软化非高斯 3 种不同概率特性风荷载作 用下的动力响应和疲劳可靠性进行分析. 算例采用 FAST ( fatigue, aerodynamics, structures and turbulence)建模,风机叶片划分为 17 个单元,塔筒 划分为 10 个单元(如图 1),塔底与地面刚接. 叶片和 塔体的材料参数和气动参数可参见文献[15]. 图 1 风机结构和风场示意图 Fig. 1 Sketch of the wind turbine and wind field 1郾 2 平均风速和来流风向 风机结构任意高度处的风荷载均可分解为平均和 脉动两部分. 对于平均风荷载,可根据指数律由参考 高度处平均风速换算得到. 平均风速( u)的概率密度 函数可按 Weibull 分布简化,即: f(u) = k ( c u ) c k - 1 ·e - (u / c) k . (1) 式中,c 和 k 分别是尺度和形状参数,可根据工程所处 地的风速统计资料按最大似然估计法(MLM)、统计矩 法(MM)和能量密度法(PDM)计算得到. 考虑到方法 的简便性及结果的准确性,本文按能量密度法对两参 数 Weibull 分布进行估计,即[16] : k = 1 + 3郾 69 (Epf) 2 ; c = u 祝 ( 1 + 1 ) k ì î í ï ïï ï ïï . (2) 式中,祝(·)是 Gamma 函数;Epf是能量系数,定义为: Epf = u 3 u 3 . (3) 如图 2 所示为 NREL 135 m 气象塔[17] 88 m 高度处 超声风速仪测得的 2014 年全年 10 min 平均风速(即轮 毂高度处的平均风速 uhub )的玫瑰图. 图 3 所示为 3 种 不同方法 Weibull 分布参数估计与实测数据的比较, ·1454·
双妙等：全风向来流非高斯风场风机疲劳寿命可靠性分析 ·1455· ☐25ms<4<30ms1 0.10 NNW ☐20m:sl<<25ms 0.09 ☐15ms<b<20msl 0.08 10ms<15 m.s 0.07 NW ne ☐5ms'<<10ms 0.06 ☐0m·s'<4<5mg 0.05 0.04 WNW ENE 0.03 0.02 0.01 0 W 0.01 0.02 0.03 0.04 WSW ESE 0.05 0.06 0.07 0.08 SW SE 0.09 0.10 SSW SSE 图2风速玫瑰图 Fig.2 Wind rose diagram 图中柱状图是实测数据，曲线分别是按PDM、MM和 0.664.,湍流尺度系数A.=0.7min(30,H),H是 MLM方法计算得到的两参数Weibull分布的概率密度 轮毂高度 函数(PDF). 本文根据EC规范，仅考虑空间中任意i,两点顺 0.16 风向脉动风荷载的相关性，表示为： -PDM --·MM coh()=e. (5) 0.12 --MLM 式中，r表示生成风速网格点间的距离；a,、b,表示衰减 和偏移系数，按EC规范取值，即a.=12和b.=3.527 008 ×10-4 风机正常风速条件下按EC规范中的设计荷载条 0.04 件(DLC)1.1和正常湍流模型(NTM)考虑载荷. 00 10 15 20 25 30 2风场模拟 a/m·g- 2.1高斯风场模拟 图3三种方法Weibull分布参数估计与实测数据的比较 根据谐波叠加法，生成平稳高斯脉动风速第m个 Fig.3 Comparison of parameters estimated by the Weibull distribu- tion using three methods and the measured data 时程在时刻的风速，表示为： N/2-l 由图3可知，3种方法的计算结果相近，但能量密 uc.n(t)=∑F(o,)ey. (6) j.-N2 度法不需通过迭代即可快速估计出Weibull分布的形 式中，F.(ω，)表示第m个待生成脉动风速时程的第j 状和尺度参数，且结果精度高，具有明显的优势 个傅里叶谱，即： 1.3脉动风速 根据脉动风速谱，采用Monte Carlo方法模拟如图 F(o)= √△aHn(o)e (7) 1所示风场，风场中心位于轮毂位置.本文脉动风速谱 式中：ω=ω，△w是圆频率间隔；P,是第r个待生成 IEC 61400-1 International Electronic Technical 时程相位谱9，的第j个随机相位角；H（ω，)是傅里叶 Commission)规范中Kaimal谱模型（式(4）)，但本文方 幅值谱，由脉动风速谱Cholesky分解得到，即： 法并不限于Kaimal谱，对其他风速谱同样成立. S(@)=H(@)H*(@). (8) 40 L/upab 式中，S(w)是目标谱矩阵，H·(w)为H(w)的共轭 S.D=1+6/ua) (4) 矩阵 式中：∫表示角频率；o,是脉动风速的标准差(k取u、v 2.2非高斯风场模拟 和w,分别表示顺风向、横风向和垂直向)，取σ。= 非高斯过程根据峰度系数(y4)的不同分为软化 0.16(0.75u+5.6),0,=0.8o.和0.=0.50.；L是 过程(y:>3)和硬化过程(y,<3).任意标准非高斯过程 湍流积分尺度，分别取L。=8.1A,L.=2.74.和L.= X(t)均可由单调函数变换成标准高斯过程U(t),即：
双 妙等: 全风向来流非高斯风场风机疲劳寿命可靠性分析 图 2 风速玫瑰图 Fig. 2 Wind rose diagram 图中柱状图是实测数据,曲线分别是按 PDM、MM 和 MLM 方法计算得到的两参数 Weibull 分布的概率密度 函数(PDF). 图 3 三种方法 Weibull 分布参数估计与实测数据的比较 Fig. 3 Comparison of parameters estimated by the Weibull distribu鄄 tion using three methods and the measured data 由图 3 可知,3 种方法的计算结果相近,但能量密 度法不需通过迭代即可快速估计出 Weibull 分布的形 状和尺度参数,且结果精度高,具有明显的优势. 1郾 3 脉动风速 根据脉动风速谱,采用 Monte Carlo 方法模拟如图 1 所示风场,风场中心位于轮毂位置. 本文脉动风速谱 采用 IEC 61400鄄鄄 1 ( International Electronic Technical Commission)规范中 Kaimal 谱模型(式(4)),但本文方 法并不限于 Kaimal 谱,对其他风速谱同样成立. Sk(f) = 4滓 2 kLk / uhub (1 + 6fLk / uhub ) 5 / 3 . (4) 式中:f 表示角频率;滓k是脉动风速的标准差( k 取 u、v 和 w,分别表示顺风向、横风向和垂直向),取 滓u = 0郾 16(0郾 75uhub + 5郾 6),滓v = 0郾 8滓u和 滓w = 0郾 5滓u ;Lk 是 湍流积分尺度,分别取 Lu = 8郾 1撰u ,Lv = 2郾 7撰u和 Lw = 0郾 66撰u ,湍流尺度系数 撰u = 0郾 7min(30, Hhub ),Hhub是 轮毂高度. 本文根据 IEC 规范,仅考虑空间中任意 i,j 两点顺 风向脉动风荷载的相关性,表示为: cohij,u (f) = e - au (f·r/ uhub) 2 + (bu r) 2 . (5) 式中,r 表示生成风速网格点间的距离;au 、bu表示衰减 和偏移系数,按 IEC 规范取值,即 au = 12 和 bu = 3郾 527 伊 10 - 4 . 风机正常风速条件下按 IEC 规范中的设计荷载条 件(DLC)1郾 1 和正常湍流模型(NTM)考虑载荷. 2 风场模拟 2郾 1 高斯风场模拟 根据谐波叠加法,生成平稳高斯脉动风速第 m 个 时程在 t 时刻的风速,表示为: uG,m (t) = 移 N/ 2 -1 j = -N/ 2 Fm (棕j)·e i棕j t . (6) 式中,Fm (棕j)表示第 m 个待生成脉动风速时程的第 j 个傅里叶谱,即: Fm (棕j) = 移 m r = 1 驻棕Hmr(棕j)·e i渍rj . (7) 式中:棕j = j驻棕,驻棕 是圆频率间隔;渍rj是第 r 个待生成 时程相位谱 渍r的第 j 个随机相位角;Hmr(棕j)是傅里叶 幅值谱,由脉动风速谱 Cholesky 分解得到,即: S(棕) = H(棕)H *T (棕). (8) 式中,S(棕) 是目标谱矩阵,H * ( 棕) 为 H( 棕) 的共轭 矩阵. 2郾 2 非高斯风场模拟 非高斯过程根据峰度系数( 酌4 ) 的不同分为软化 过程(酌4 >3)和硬化过程(酌4 < 3). 任意标准非高斯过程 X(t)均可由单调函数变换成标准高斯过程 U(t),即: ·1455·
·1456· 工程科学学报，第39卷，第9期 x=g(u)=F:[Φ(u)]. (9) 45m a 式中：u和x分别是U(t)和X()的随机变量；中(·)表 95%C 示标准高斯过程的概率分布函数(CDF);F表示标 30 准非高斯过程CDF的逆函数；g(·)表示非高斯穿越过 程，可根据Hermite矩模型us]或CDF映射等方法 95%CL 得到. 100 200 300 400 500 600 对于非高斯风荷载模拟，首先需对日标谱S(ω) 45 标准化，即： (b) 95%CI o)8 j,k=1,2,…,m. (10) 二 式中，Oc是第j个时程的标准差. 再根据式(6)生成标准高斯过程的风速时程 95%C ((t))后，通过非高斯转化过程（式(9）)生成标准 100 200 300400 500 600 t切s 45 非高斯过程的风速时程(())： (c 95%CI+ u8(t)=g[u2(t)] (11) 为了提高非高斯过程的准确性，可对生成时程的 功率谱密度进行迭代，即20-2]： S(ω)= TS.c(o）1 Is2()s () (12) 100 200 300400500600 式中，B是收敛速度和收敛精度的优化系数，通常取 图4轮毂高度处风速时程.(a)硬化过程：(b)高斯过程：(c)软 B=0.3] 化过程 2.3非高斯风场正常风速模拟 Fig.4 Wind speed time history at the hub height:(a)hardening 在风机运行状态，将切入风速(3ms·)到切出风 process;(b)Gaussian process:(c)softening process 速(25ms1)划分为12个单元，单元宽度为2ms; 在风机停机状态，考虑低于切入风速的1ms和2m· 系统对结构动力响应有较大影响.因此本文采用基于 s风速，共14种不同平均风速的硬化、软化和高斯过 耦合气动模型的FAST程序计算结构的动力响应.为 程风场.本文按Hermite矩模型生成非高斯风场，软化 考虑来流风向影响，按图2所示将来流划分为16个风 过程的偏斜和峰度分别取0和4.5，硬化过程取0和 向角. 2.0.每种情况按相同的随机种子模拟硬化、软化和高 风机结构动力响应分析中，采用叶素动量理论考 斯风场各50条，每条时程的长度为630s,时间间隔为 虑叶片尾流模型，按Beddoes-Leishman模型考虑叶片 0.05$.为了消除启动瞬间的影响，从30s后对结构动 的气动效应.由FAST程序计算出塔体底部控制截面 力响应进行分析 的内力后，根据材料属性和几何尺寸计算出控制截面 图4所示为轮毂高度处平均风速25m·s时，按 不同位置的应力.由于面内弯曲应力远大于面外弯曲 式(6)和式(11)生成的硬化、高斯和软化风荷载.图 应力，因此本文仅对控制截面面内弯曲正应力(S,)的 中虚线表示均值加减1.96倍标准差（对于高斯分布即 动力响应进行分析 95%置信度区间C).图中小图比较了目标谱与生成 图5所示为风机结构在三种不同概率特性风荷载 荷载的模拟谱，右图实线为与风荷载均值和方差相同 作用下，控制截面面内最大弯曲正应力的响应时程 的高斯分布概率密度函数(PDF),柱状图为模拟风荷 可以看出，风机结构应力响应的非高斯性较风荷载的 载PDF.由图可知，三种不同风场功率谱密度与目标 非高斯性大大减弱. 谱均吻合较好，软化过程和硬化过程均表现出明显的 图6所示为不同风速区间内不同概率特性风场各 非高斯性.在相同目标谱、平均风速和湍流度条件下， 50条计算得到的应力响应的均值、标准差、偏态和峰 高斯过程的极值大小和数量介于硬化和软化过程 度的平均值.从图中可知在高斯和非高斯风荷载作用 之间. 下，均值、标准差和偏态系数基本相同.然而，由于高 阶矩的不稳定性，3种荷载作用下动力响应的峰度系 3风机结构的动力响应分析 数有所不同，但整体趋势仍保持一致 在风荷载作用下，叶片与来流、叶片与塔体间均相 图7所示为轮毂高度处平均风速25m·s条件 互耦合，动力响应分析具有高度的非线性.同时，控制 下，3种不同概率特性风荷载作用下，控制截面弯曲正
工程科学学报,第 39 卷,第 9 期 x = g(u) = F - 1 X [椎(u)]. (9) 式中:u 和 x 分别是 U(t)和 X(t)的随机变量;椎(·)表 示标准高斯过程的概率分布函数(CDF);F - 1 X 表示标 准非高斯过程 CDF 的逆函数;g(·)表示非高斯穿越过 程,可根据 Hermite 矩模型[18] 或 CDF 映射[19] 等方法 得到. 对于非高斯风荷载模拟,首先需对目标谱 S T NG (棕) 标准化,即: S 0 NG,jk(棕) = ST,NG,jk(棕) 滓NG,j滓NG,k ,j,k = 1,2,…,m. (10) 式中,滓NG,j是第 j 个时程的标准差. 再根据式 ( 6 ) 生 成 标 准 高 斯 过 程 的 风 速 时 程 (u 0 G (t))后,通过非高斯转化过程(式(9)) 生成标准 非高斯过程的风速时程(u 0 NG (t)): u 0 NG (t) = g[u 0 G (t)]. (11) 为了提高非高斯过程的准确性,可对生成时程的 功率谱密度进行迭代,即[20鄄鄄21] : S (i + 1) G (棕) = [ ST,NG (棕) S (i) NG (棕 ] ) 茁 S (i) G (棕). (12) 式中,茁 是收敛速度和收敛精度的优化系数,通常取 茁 = 0郾 3 [22] . 2郾 3 非高斯风场正常风速模拟 在风机运行状态,将切入风速(3 m·s - 1 )到切出风 速(25 m·s - 1 )划分为 12 个单元,单元宽度为 2 m·s - 1 ; 在风机停机状态,考虑低于切入风速的 1 m·s - 1和 2 m· s - 1风速,共 14 种不同平均风速的硬化、软化和高斯过 程风场. 本文按 Hermite 矩模型生成非高斯风场,软化 过程的偏斜和峰度分别取 0 和 4郾 5,硬化过程取 0 和 2郾 0. 每种情况按相同的随机种子模拟硬化、软化和高 斯风场各 50 条,每条时程的长度为 630 s,时间间隔为 0郾 05 s. 为了消除启动瞬间的影响,从 30 s 后对结构动 力响应进行分析. 图 4 所示为轮毂高度处平均风速 25 m·s - 1时,按 式(6)和式(11)生成的硬化、高斯和软化风荷载. 图 中虚线表示均值加减 1郾 96 倍标准差(对于高斯分布即 95% 置信度区间 CI ± ). 图中小图比较了目标谱与生成 荷载的模拟谱,右图实线为与风荷载均值和方差相同 的高斯分布概率密度函数( PDF),柱状图为模拟风荷 载 PDF. 由图可知,三种不同风场功率谱密度与目标 谱均吻合较好,软化过程和硬化过程均表现出明显的 非高斯性. 在相同目标谱、平均风速和湍流度条件下, 高斯过程的极值大小和数量介于硬化和软化过程 之间. 3 风机结构的动力响应分析 在风荷载作用下,叶片与来流、叶片与塔体间均相 互耦合,动力响应分析具有高度的非线性. 同时,控制 图 4 轮毂高度处风速时程. (a)硬化过程;( b)高斯过程;( c)软 化过程 Fig. 4 Wind speed time history at the hub height: ( a) hardening process; (b) Gaussian process; (c) softening process 系统对结构动力响应有较大影响. 因此本文采用基于 耦合气动模型的 FAST 程序计算结构的动力响应. 为 考虑来流风向影响,按图 2 所示将来流划分为 16 个风 向角. 风机结构动力响应分析中,采用叶素动量理论考 虑叶片尾流模型,按 Beddoes鄄鄄 Leishman 模型考虑叶片 的气动效应. 由 FAST 程序计算出塔体底部控制截面 的内力后,根据材料属性和几何尺寸计算出控制截面 不同位置的应力. 由于面内弯曲应力远大于面外弯曲 应力,因此本文仅对控制截面面内弯曲正应力( Sb )的 动力响应进行分析. 图 5 所示为风机结构在三种不同概率特性风荷载 作用下,控制截面面内最大弯曲正应力的响应时程. 可以看出,风机结构应力响应的非高斯性较风荷载的 非高斯性大大减弱. 图 6 所示为不同风速区间内不同概率特性风场各 50 条计算得到的应力响应的均值、标准差、偏态和峰 度的平均值. 从图中可知在高斯和非高斯风荷载作用 下,均值、标准差和偏态系数基本相同. 然而,由于高 阶矩的不稳定性,3 种荷载作用下动力响应的峰度系 数有所不同,但整体趋势仍保持一致. 图 7 所示为轮毂高度处平均风速 25 m·s - 1 条件 下,3 种不同概率特性风荷载作用下,控制截面弯曲正 ·1456·
双妙等：全风向来流非高斯风场风机疲劳寿命可靠性分析 ·1457· 70 60 60 (a 95%CI+ ( 。硬化过程 口软化过程 一高斯过程 002010 均值 40 0 0 0-0 95%CI -10 100 200 300 400500 600 20 70 标准差 (b) 60 95%CI+ 0000000 50 0 10 15 20 25 2 m·s） 0 。硬化过程 0 95%C 口软化过程 -10 一高斯过程 100 200 300 400 500 600 tis 70 峰度 60 (c) 95%C+ 0888 50 40 30 偏态 20 10 0 95%CI -1 0 100 200 300 400 500 600 10 15 20 25 tis 丛/m·s) 图5风荷载作用下面内弯曲应力时程.(a)硬化过程：(b)高斯 图63种不同概率特性风荷载作用下控制截面面内弯曲正应力 过程；(c)软化过程 的前四阶统计矩.(a)均值与标准差：(b)偏态和峰度 Fig.5 Time history of the in-plane bending stress of the critical sec- Fig.6 First four statistical moments of the in-plane bending stress of tion under a wind load:(a)hardening process;(b)Gaussian the critical section under wind loads with different probability charac- process;(c)softening process teristics:(a)mean and standard deviation;(b)skewness and kurto- sis 应力的功率谱密度(PSD).由于动力响应时程的时间 102 间隔为0.O5s,根据Nyquist抽样定理，PSD的有效频率 硬化过程 1P 高斯过程 范围为[0,10Hz].由图可知，三种不同概率特性风荷 10°F 2P 软化过程 载作用下弯曲应力的PSD基本相同. 4P 图8是3种不同概率特性风场作用下，由控制截 102 1P:0.327Hz 面面内弯曲正应力的功率谱密度计算得到的带宽系数 's 10E 2P:0.578Hz α、a,和q(式(13)).根据随机振动理论，当g≈0或 3P:1.302Hz 4P:2.915Hz α心，≈1时，应力过程是窄带过程.由图8可知，3种应 10 力过程均表现出明显的宽带特性，且带宽系数基本 10-2 10 10 10 相同. f/Hz 入 图73种不同概率特性风荷载作用下控制截面面内弯曲应力的 a=() 功率谱密度 入2 Fig.7 PSDs of the in-plane bending stress of the critical section un- a2= (oA4)迈， (13) der wind loads with different probability characteristics 9=/1- 4疲劳可靠性分析 入。A2 如图9所示为硬化、高斯和软化风荷载作用下，不 结构的使用寿命通常由疲劳裂纹形成寿命和裂纹 同风速区间内由50个应力时程得到的响应最大值. 扩展寿命两部分组成.本文在考虑全风向来流条件 结合图6和图9可知，在均值和方差近似相等的条件 下，采用线性损伤累积理论对裂纹形成寿命进行计算， 下，软化、高斯和硬化过程的最大值依次减小，且软化 并考虑由结构和荷载两方面不确定性引起的疲劳损伤 过程最大值的离散程度最大 不确定性[2a].同时，采用Paris公式计算疲劳裂纹扩展
双 妙等: 全风向来流非高斯风场风机疲劳寿命可靠性分析 图 5 风荷载作用下面内弯曲应力时程. ( a)硬化过程;( b)高斯 过程;(c)软化过程 Fig. 5 Time history of the in鄄plane bending stress of the critical sec鄄 tion under a wind load: ( a ) hardening process; ( b ) Gaussian process; (c) softening process 应力的功率谱密度(PSD). 由于动力响应时程的时间 间隔为 0郾 05 s,根据 Nyquist 抽样定理,PSD 的有效频率 范围为[0, 10 Hz]. 由图可知,三种不同概率特性风荷 载作用下弯曲应力的 PSD 基本相同. 图 8 是 3 种不同概率特性风场作用下,由控制截 面面内弯曲正应力的功率谱密度计算得到的带宽系数 琢1 、琢2和 q(式(13)). 根据随机振动理论,当 q抑0 或 琢2抑1 时,应力过程是窄带过程. 由图 8 可知,3 种应 力过程均表现出明显的宽带特性,且带宽系数基本 相同. 琢1 = 姿1 (姿0姿2 ) 1 / 2 , 琢2 = 姿2 (姿0姿4 ) 1 / 2 , q = 1 - 姿 2 1 姿0姿2 ì î í ï ï ï ï ï ï ï ï . (13) 如图 9 所示为硬化、高斯和软化风荷载作用下,不 同风速区间内由 50 个应力时程得到的响应最大值. 结合图 6 和图 9 可知,在均值和方差近似相等的条件 下,软化、高斯和硬化过程的最大值依次减小,且软化 过程最大值的离散程度最大. 图 6 3 种不同概率特性风荷载作用下控制截面面内弯曲正应力 的前四阶统计矩. (a)均值与标准差; (b)偏态和峰度 Fig. 6 First four statistical moments of the in鄄plane bending stress of the critical section under wind loads with different probability charac鄄 teristics: (a)mean and standard deviation; ( b) skewness and kurto鄄 sis 图 7 3 种不同概率特性风荷载作用下控制截面面内弯曲应力的 功率谱密度 Fig. 7 PSDs of the in鄄plane bending stress of the critical section un鄄 der wind loads with different probability characteristics 4 疲劳可靠性分析 结构的使用寿命通常由疲劳裂纹形成寿命和裂纹 扩展寿命两部分组成. 本文在考虑全风向来流条件 下,采用线性损伤累积理论对裂纹形成寿命进行计算, 并考虑由结构和荷载两方面不确定性引起的疲劳损伤 不确定性[23] . 同时,采用 Paris 公式计算疲劳裂纹扩展 ·1457·