由于国内气象台站的风向资料记录不全且记录不够准确,许多学者在研究风机塔筒结构风致响应及疲劳寿命分析时仅考虑风速的影响,将风向对结构风振响应和疲劳损伤的影响加以忽略 ^[1,2,3,4]。根据高耸结构的研究结论,忽略风向的影响是不够准确且偏于保守的。

　　 因此,为了准确预测风机塔筒结构在风荷载作用下的疲劳寿命,必须考虑风向的影响。目前常用的做法是确定风机位置处的风速风向联合分布函数。而确定该函数最直接的方法是利用在风机位置处安装的风速风向仪长期观测以获得准确的气象资料。但由于风速风向仪安装在机舱上,对于上风向风机来说,此处有叶片旋转带来的漩涡系和塔影效应的影响,来流风速的大小和方向均会发生改变,最终导致实测数据局部失真。同时风速风向资料属于风机制造企业核心机密,一般较难得到。目前学术界常利用风机位置附近的气象站观测数据统计得到气象站位置的风速风向联合分布函数,进而得到风机位置处的风速风向联合分布函数,以用于塔筒结构的疲劳寿命分析。由于此种方法可应用于所有高层建筑和高耸结构 ^[5,6],是目前解决风向对疲劳寿命影响问题最有效的方法。

　　 本文基于气象站的实际观测数据,推导得到风机位置处的风速风向联合分布函数,采用时域分析方法对不同风速风向工况的风机塔筒结构的风振响应进行分析,并与实测响应进行对比验证,系统研究了风速风向对风机塔筒结构动力响应及疲劳寿命的影响,为风机结构的经济性设计提供技术参考。

　　 由于国内风机建造地区缺乏长期的风速数据观测资料,因此风速统计存在数据不足及由短期记录推测长期记录的问题。如果以年最大风速抽样得到极值风速样本,由于样本数太少,在有些风向区间低风速出现的概率就为0,这就会导致遗漏很多风向区间小风速工况。研究发现,较低风速也会对风机塔筒结构产生疲劳损伤,不能加以舍弃,因此本文采用阶段极值法进行极值样本抽样。由于抽样时间过短会导致风速极值样本存在相关性,因此建议采用每间隔8d抽取风速极值样本数据。

　　 本文研究所采用的1.25MW风机结构位于河北省张家口市张北县。张家口地区地形总体西北高,东南低,阴山山脉横贯中部,将全区分为坝上和坝下两部分。坝上地区是内蒙古高原的一部分,海拔高度一般在1 400～1 600m,张北县属于坝上地区。坝下地区河谷盆地相间,地形复杂,海拔一般在500～900m。受地形和海拔影响,坝上与坝下地区的平均风速、主导风向差异较大。因此在无法获得张家口地区地面观测站风速风向观测资料的情况下,综合考虑距离的远近、海拔高度、地形特征和主导风向的影响,最终选择山西大同气象站和内蒙古化德气象站的地面气象数据,风速样本数据主要来自于国家气象科学数据共享服务平台。

　　 山西大同气象站记录的是从1955年1月到2017年4月共63年的10m高度处日最大风速(10min平均风时距)和相应的风向数据。内蒙古化德气象站记录的是从1952年12月到2017年4月共66年的10m高度处日最大风速(10min平均风时距)和相应的风向数据。图1为山西大同和内蒙古化德气象站各个风速风向区间上风速极值样本的柱状分布图。两个气象站的风玫瑰图如图2所示。从图2可以看出,山西大同气象站气象资料中N向风(对应于图1的0°风向角)出现的频率最大,为23.07%; 在所有风向区间中出现频率最大的风速区间为8～10m/s,频率为23.47%。内蒙古化德气象站WNW向风出现的频率最大,为37%。在所有的风向区间中出现频率最大的风速区间为10～12m/s,其值为25.50%。从图1和图2的对比分析可以看出,大同气象站和化德气象站在各个风向区间出现的概率趋势是一致的,这也为后续利用气象站的风速风向联合分布函数推导得到风机位置处的风速风向联合分布函数中风向相同的假设提供了合理的基础。

　　 在得到极值风速样本在各风速风向区间内的频率分布后,就可以根据统计结果得到气象站位置处的风速风向联合分布函数。

　　 杨咏昕等 ^[7]和王钦华 ^[6]认为风速风向概率分布函数P(x,θ)可以由风向频度函数f(θ)(即各个风向区间出现的频率)和各风向上的极值风速概率分布函数F(x)共同确定,两个函数之间是相互独立的,如式(1)所示:

　　 ${\rm P}(x,\theta )=f(\theta )F(x)(1)$

　　 式中:x为极值风速变量; θ为风向角。

　　 一般极值风速概率分布函数F(x)分为三种:极值Ⅰ型(Gumbel)分布、极值Ⅱ型(Frechet)分布和极值Ⅲ型(Weibull)分布,对应的风速风向联合分布函数P(x,θ)具体表达式如式(2)～(4)所示。

　　 极值Ⅰ型分布:

　　 ${\rm P}(x,\theta )=f(\theta )\exp \left[-\exp \left(-\frac{x-b(\theta )}{a(\theta )}\right)\right](2)$

　　 极值Ⅱ型分布:

　　 ${\rm P}(x,\theta )=f(\theta )\exp \left[-\left(\frac{x-b(\theta )}{a(\theta )}\right){}^{-\gamma (\theta )}\right](3)$

　　 极值Ⅲ型分布:

　　 ${\rm P}(x,\theta )=f(\theta )\left\{1-\exp \left[-\left(\frac{x-b(\theta )}{a(\theta )}\right){}^{\gamma (\theta )}\right]\right\}(4)$

　　 式中:a(θ)为尺度参数; b(θ)为位置参数; γ(θ)为形状参数。

　　 由于风速风向联合分布函数P(x,θ)中,f(θ)中没有参数,故本文仅对极值风速概率分布函数F(x)采用矩法、最小二乘法、Gringorten法、极大似然函数法、可靠性矩法和概率权值法进行参数估计。为了衡量不同参数估计方法的优劣性,需要进行拟合优度检验,本文采用四种参数估计优良性指标,即:拟合标准差、拟合相对偏差、相关系数和柯尔莫哥洛夫检验法(K-S检验)。

　　 根据2个气象站极值风速样本数据,最终得到山西大同气象站极值风速样本更接近极值Ⅰ型分布,采用最小二乘法可以对16个风向区间的极值Ⅰ型分布参数进行估计。内蒙古化德气象站极值风速更接近极值Ⅰ型分布,采用极大似然法可以用于对16个风向区间的极值Ⅰ型分布所含参数进行估计。

　　 本文限于篇幅,此处仅给出后续风致疲劳分析所采用的内蒙古化德气象站各个风向区间下极值Ⅰ型分布中所含参数的估计结果,如表1所示。

　　 Coles等 ^[8]提出采用谐波函数来拟合得到风向频度函数f(θ)和极值风速概率分布函数F(x)所含参数沿整个圆周上的变化,如式(5a)～(5c)所示。

　　 对于极值Ⅰ型分布:

　　 $\begin{array}{l}f(\theta )=L{}^{\text{θ}}+{\displaystyle \sum _{i=1}^{n{}_{\text{θ}}}{\rm M}}{}^{\text{θ}}\cos (i\theta -{\rm N}{}^{\text{θ}})(5\text{a})\\ a(\theta )=L{}^{\text{a}}+{\displaystyle \sum _{i=1}^{n{}_{\text{a}}}{\rm M}}{}^{\text{a}}\cos (i\theta -{\rm N}{}^{\text{a}})(5\text{b})\\ b(\theta )=L{}^{\text{b}}+{\displaystyle \sum _{i=1}^{n{}_{\text{b}}}{\rm M}}{}^{\text{b}}\cos (i\theta -{\rm N}{}^{\text{b}})(5\text{c})\end{array}$

　　 式中:L,M和N分别为谐波函数中的待定参数,上标θ,a和b分别用于区分风向频度函数、尺度参数和位置参数; n_θ,n_a和n_b分别为风向频度函数、尺度参数和位置参数的阶数。

　　 内蒙古化德气象站风向频度函数f(θ)和极值Ⅰ型分布尺度参数和位置参数的曲线拟合表达式如式(6a)～(6c)所示 ^[9]:

　　 f(θ)=0.061+0.096cos(θ+1.355)-0.063cos(2θ-0.925)-0.045cos(3θ+0.036)-0.034cos(4θ+1.339)+0.031cos(5θ-0.7) (6a)

　　 a(θ)=2.709+0.999cos(θ-0.075)+0.477cos(2θ-0.004)+0.483cos(3θ-33.062)-0.346cos(4θ)+0.096cos(5θ-0.541) (6b)

　　 b(θ)=8.973+1.880cos(θ+1.248)-0.919cos(2θ-39.246)+0.577cos(3θ+13.371)+0.309cos(4θ-0.465)+0.262cos(5θ-1.110) (6c)

　　 根据第1.1节的研究结论,假设两个气象站的风向是一致的,由于气象站和风机位置距离较近,因此可以认为两个位置处的梯度风速相等。在此基础上根据气象站的风速风向联合分布函数推导得到风机位置处的风速风向联合分布函数。经推导,风机位置处的极值Ⅰ型分布风速风向联合分布函数如式(7)所示 ^[9]:

　　 $\begin{array}{l}{\rm P}(\stackrel{—}{{\displaystyle U}}{}_{\text{W}}>x,\theta )\\ =f(\theta )\{1-\exp \left[-\exp \left(-\frac{\left(\frac{h{}_{\text{Μ}}}{{\rm H}{}_{\text{Μ}}}\right){}^{\alpha {}_{\text{Μ}}}\left(\frac{{\rm H}{}_{\text{W}}}{h{}_{\text{W}}}\right){}^{\alpha {}_{\text{W}}}x-b(\theta )}{a(\theta )}\right)\right]\}\\ =f(\theta )\left\{1-\exp \left[-\exp \left(-\frac{\xi x-b(\theta )}{a(\theta )}\right)\right]\right\}(7)\end{array}$

　　 其中:

　　 $\xi =\left(\frac{h{}_{\text{Μ}}}{{\rm H}{}_{\text{Μ}}}\right){}^{\alpha {}_{\text{Μ}}}\left(\frac{{\rm H}{}_{\text{W}}}{h{}_{\text{W}}}\right){}^{\alpha {}_{\text{W}}}$

　　 式中:$\stackrel{—}{{\displaystyle U}}{}_{\text{W}}$为风机位置任意高度处的极值风速; 变量h_M和h_W分别为气象站和风机位置处的风速高度; H_M和H_W分别为气象站和风机位置处的梯度风高度; α_M和α_W分别为气象站和风机位置处的地面粗糙度指数。

　　 将第1.3节拟合得到的风向频度函数、极值Ⅰ型分布尺度参数和位置参数的函数表达式带入到式(7),即可得到基于山西大同气象站和内蒙古化德气象站气象数据的风机位置处的风速风向联合分布函数,如图3所示。

　　 由于本文所研究的某1.25MW风机结构位于河北张家口市张北县,与内蒙古化德县同属坝上高原,此处的海拔高度、地形条件、平均风速和主导风向和位于坝下地区的山西大同市差别较大,这必将对风机塔筒结构的风致疲劳分析结果产生较大的影响,因此,本文采用以内蒙古化德气象站的气象资料为基础计算得到的风速风向联合分布函数来进行风机塔筒结构风致疲劳分析。

　　 本文研究的风机轮毂高度82.35m,塔筒总高度为80m,塔筒总质量为202t。塔筒为变截面钢锥筒,共分4段,从塔底向上各段高度分别为8,18,24,30m。每个塔段的壁厚不变,从塔筒顶部向下各塔段的壁厚分别为18,30,42,52mm。塔顶高度处塔筒外径为2.58m,塔底高度处塔筒外径为4.2m,其余各截面外径沿高度成线性变化。底段塔筒开设椭圆形门洞并安置门框,门洞高度为3.0m,宽度为0.88m,距离塔底1.2m。另外,风轮直径108.76m,风机结构为三叶片形式,机舱的长为13.6m,宽度为4.7m,高为4.7m,质量为85t。基础采用圆截面钢筋混凝土筏板基础,直径为10m,板厚1.5m。通过耦合和约束方程将不同部件进行连接,最终形成风机结构整体有限元模型,见图4。

　　 由于风机结构是风敏感结构,其非线性较为显著。如果采用频域分析,必须对原结构进行线性化处理,误差较大。时域分析方法是业内唯一认可的风机结构风振响应计算方法。该方法精度较高,可以直接了解结构的动力性能和疲劳信息,非常适合处理响应为非高斯应力过程和非线性结构的疲劳问题。本文采用时域分析方法进行风机塔筒结构的风振响应分析,并与已有现场实测结果进行对比。叶片的旋转效应通过引入旋转傅里叶谱模型来考虑 ^[9],该模型的引入将叶片转动的运动学问题转化成静力学问题,极大地方便了风机结构的动力响应分析。

　　 Thomsen等 ^[10]研究发现,疲劳损伤主要来源于风机正常运行状态下的荷载设计工况,占到总疲劳损伤的99.5%。因此,本文主要对正常运行状态时风机塔筒结构在不同风向和风速下的风振响应进行分析,未考虑风机启动、关机和停机等工况下的风振响应及疲劳损伤。

　　 本文将风机结构的空气动力荷载分为风轮叶片气动力荷载和塔筒结构气动力荷载。作用在塔筒结构上的气动力荷载主要考虑顺风向气动力荷载F_D(z,t)、漩涡脱落导致的横风向升力F_V(z,t)和横向湍流分量引起的横风向气动力荷载F_L(z,t)。作用在风轮叶片上的空气动力荷载主要考虑轴向推力Q_x和切向力Q_y。其中,塔筒结构顺风向气动力荷载F_D(z,t)可以利用准定常方法确定 ^[11]; 漩涡脱落导致的横风向升力F_V(z,t)在风速处于亚临界和跨临界区域时采用简谐升力 ^[11],在风速处在超临界区域时基于Vickery提出的升力功率谱 ^[12]得到; 横向湍流分量引起的横风向气动力荷载F_L(z,t)根据Holmes ^[13]提出的合成理论基于准定常假设得到。作用在风轮叶片上的空气动力荷载主要基于叶素-动量理论(BEM理论)及其相关修正理论 ^[14]得到。

　　 在任意风向角某一特定风速下,作用在风机结构上的空气动力荷载示意如图5所示。图5中X轴和Y轴为塔筒结构的体轴。习惯上,在进行结构风振响应计算时要把顺风向和横风向的气动力荷载投影到结构的体轴上。由于风机上的偏航系统和控制系统相配合,使风力发电机的风轮始终处于迎风状态,因此本文不考虑风向的变化对叶片气动力荷载的影响,仅考虑风向的变化对塔筒结构风振响应及疲劳寿命的影响。

　　 将每个风向角范围[θ_i,θ_i+1]下风速范围$[(\stackrel{—}{{\displaystyle U}}{}_{\text{h}\text{u}\text{b}}){}_j,(\stackrel{—}{{\displaystyle U}}{}_{\text{h}\text{u}\text{b}}){}_{j+1}]$定义为第[i,j]个工况,风向角间隔Δθ_i=22.5°。轮毂高度处的平均风速$\stackrel{—}{{\displaystyle U}}{}_{\text{h}\text{u}\text{b}}$介于切入风速(4m/s)和切出风速(24m/s)之间,本文取区间端点值进行计算,风速间隔为2m/s。

　　 典型风速下,不同风向角时塔筒结构顺风向和横风向位移响应沿高度的变化规律分别如图6和图7所示。从图6,7可以看出,对于任意风向角,塔筒顺风向和横风向位移均值和均方根值均随着塔筒高度的增加而逐渐增大。但是风向角对塔筒顺风向和横风向位移的影响规律不一致。当风向角为22.5°时,顺风向位移均值和均方根值达到最大值,而横风向位移均值和均方根值则在风向角为67.5°达到最大值 ^[15],主要原因为顺风向和横风向气动力荷载的作用机理不同。

　　 区别于高耸结构和超高层建筑等风敏感建筑,塔顶横风向位移均值也不为零(图7(a)),其主要原因为在叶片剖面上作用着沿横风向不均匀分布的切向力(图5)。

　　 图8为典型风向角下塔顶顺风向位移均值和均方根值随风速的变化曲线。由于顺风向和横风向位移响应随风速的变化规律相同,即塔顶位移均值和均方根值均随着风速的增大而增大。此处仅给出顺风向位移响应结果。

　　 值得注意的是,由于风机制造商对叶片翼型参数数据的保密性,当风速大于额定风速(18m/s)而小于切出风速(24m/s)时,没有考虑叶片桨距角的变化对塔筒结构动力响应的影响。

　　 本文依据张家口地区该风机结构的现场实测数据,将测得的加速度响应和数值分析得到的风振响应结果进行对比分析,以验证时域分析方法的正确性 ^[9]。基于风环境的高度随机性和风机内部附加振动的影响,目前国内外文献中鲜有学者提出一种合理的评价指标来评价实测响应和数值计算结果的差异性。笔者认为只要实测响应和数值计算结果处在同一量级上,就能较好地反映本文建立的整体有限元模型与风振响应动力学模型的准确性,模型即可用于后续风机塔筒结构疲劳寿命分析。

　　 测点布置在机舱顶部,观测数据包含10min内风速、风向、机舱高度处沿机舱主轴和垂直于主轴两个方向的振动加速度时程。测试采样时间为10min,采样频率为50Hz。本文共选择9种实测工况(对应于9个测试数据样本)进行对比。

　　 由于各测点实测加速度均值均接近于0,与数值模拟得到的相应点加速度响应均值(其均值为0)较为吻合。此处仅给出各工况下典型测点的加速度均方根值与对应数值模拟结果的对比图,如图9所示。从图9可以看出,实测加速度均方根值与数值模拟得到的均方根值总体差距不大,均在同一个数量级上。两者的差距随着风速的增大(对应工况6和工况5)而增大。

　　 典型工况下加速度自相关函数和响应功率谱对比分别如图10和图11所示。实测加速度自相关函数、加速度响应功率谱和数值模拟得到的结果吻合较好,均在同一个量级上。

　　 上述实测响应和数值模拟得到的响应存在差距的原因可能是,结构整体建模与实际结构存在一定的差异,叶片在正常转动过程中不能实时精确的对准风向,导致气动力无法与叶片的质心重合,引起交变荷载,阵风、风速、风向快速变化引起的交变荷载等,这些因素都会对正常运行状态下的风机产生振动影响,最终使实测响应和数值模拟得到的响应值存在不可避免的误差。

　　 采用时域分析方法得到的疲劳寿命T_f的计算公式如下式所示。

　　 $\begin{array}{l}{\rm T}{}_{\text{f}}=\frac{C}{\Omega }(8)\\ \Omega ={\displaystyle \sum _{i=1}^p{\displaystyle \sum _{j=1}^q{\rm P}}}{}_{ij}f{}_{L{}_{ij}}{\displaystyle \sum _{m=1}^{k}f}{}_{\text{Δ}\sigma {}_{\text{m}}}\left(\text{Δ}\sigma {}_{\text{m}}\right){}^{\beta }(9)\end{array}$

　　 式中:Ω为应力参数; C为疲劳强度S-N曲线中与材料性质有关的常数,S-N曲线采用我国《钢结构设计标准》(GB 50017—2017) ^[16]中规定的曲线形式,其中S为疲劳应力幅,N为循环次数; β为S-N曲线的斜率; P_ij为每个风速风向工况出现的概率,可由风速风向联合分布函数确定; f_Δσm为第[i,j]个工况下第m级应力幅Δσ_m出现的频率,可从雨流计数法计算的应力幅频数分布直方图中得到; p为风向角的总数; q为参与风振响应计算的风速的总数;k为疲劳应力幅的总个数;f_Lij为作用在结构上单位时间内应力幅的平均次数,计算公式如下:

　　 $f{}_{L{}_{ij}}=\frac{{\rm N}{}_{L{}_{ij}}}{L{}_{ij}}(10)$

　　 式中N_Lij为第[i,j]工况下应力时程作用时间L_ij内用雨流计数法得到的应力幅的总循环次数。

　　 根据《风力发电机组 塔架》(GBT 19072—2010) ^[17]的相关规定,塔筒门框开口区域和筒体连接区域应进行塔筒的疲劳分析。同时,在实际工程中,塔筒底部与基础的连接区域 ^[18,19]也容易发生疲劳失效。本文选择塔筒与门框连接区域及塔筒底部与基础连接区域这两个易发生疲劳失效的区域来进行风机塔筒结构疲劳寿命分析。塔筒与门框连接区域的构造与德国规范 ^[20]中连接类别526相符合,此时2×10⁶次对应的疲劳强度Δσ_c=90MPa。塔筒底部与基础的连接区域的构造与德国规范 ^[20]中连接类别232相符合,此时2×10⁶次对应的疲劳强度Δσ_c=71MPa。本文以沿塔筒高度方向的正应力σ_z作为风机塔筒结构疲劳分析的名义应力。

　　 根据疲劳寿命时域分析理论,编制MATLAB程序进行风机塔筒结构疲劳寿命分析计算。从而得到风机塔筒结构的疲劳寿命为22.86年。

　　 为了研究疲劳累积损伤沿风向的分布,将各风向上对应不同风速的损伤进行合计,得到疲劳累积损伤随风向角的变化规律,如图12所示。类似地,将各风速上对应不同风向角的损伤进行合计,得到疲劳累积损伤随风速的变化规律,如图13所示。

　　 从图12中可以看出,270°～337.5°风向角下对应的疲劳累积损伤较大,其他风向角下对应的疲劳累积损伤值较小。这与图3中风机位置处各风向工况出现的概率分布特征完全一致。由此说明,风向对结构的疲劳累积损伤影响较大,风向出现概率较大的区间产生的疲劳累积损伤较大。同时,从图13可看出,当风速区间介于16～24m/s时,塔筒结构疲劳累积损伤较大。

　　 为了更全面地描述塔筒结构疲劳累积损伤在各个风速风向工况下的分布规律,图14给出了疲劳寿命T_f期限内各风速风向工况疲劳累积损伤三维柱状分布图。

　　 根据疲劳寿命时域分析理论,得到不考虑风向角影响下的风机塔筒结构疲劳寿命为6.6年。

　　 显然,若不考虑风向角对塔筒结构疲劳累积损伤的影响,计算得到的疲劳寿命不满足标准风机设计时对塔筒最低使用寿命为20年以上的要求(德国规范 ^[20]),且过于保守,从而导致设计人员采取纵向加强筋或者环形加强筋对塔筒结构门框开口区域进行加强,最终影响到风机结构设计的经济性。

　　 图15为不考虑风向角条件下疲劳累积损伤随风速的变化规律。对比图13和图15可以看出,虽然不考虑风向影响的疲劳累积损伤仍然在风速区间16～24m/s时较大,但是其疲劳累积损伤在各风速区间的分布规律有较大改变。

　　 因此,建立精确的风机塔筒结构的疲劳寿命计算模型必须考虑风向角的影响,只有这样才能同时兼顾风机设计的可靠性和经济性要求。

　　 (1)风向角对塔筒位移的影响规律不一致。其中,当风向角为22.5°时顺风向位移均值和均方根值达到最大值,而横风向位移均值和均方根值则在风向角为67.5°达到最大值。

　　 (2)实测加速度响应和数值模拟得到的加速度响应均在一个数量级上,验证了本文所用方法的准确性,可用于风机塔筒结构疲劳寿命分析。

　　 (3)风向对塔筒结构的疲劳累积损伤影响较大,在风向出现概率较大的区间产生的疲劳累积损伤较大。与考虑风向角下得到的疲劳寿命相比,不考虑风向角影响时,计算得到的塔筒结构疲劳寿命偏小,但计算结果过于保守。因此在进行塔筒结构疲劳寿命分析时,需考虑风向对结构疲劳累积损伤的影响。