考虑叶片停机位置大型风力机塔架风-沙致结构响应分析

柯世堂　董依帆

摘要： 强风停机状态下叶片位置会显著影响风力机塔架的绕流及稳定性能，尤其在沙尘暴极端天气条件下，沙粒的附着也会影响风的湍流特征并对塔架产生附加冲击力，现有工作均缺乏风?沙耦合作用对大型风力机体系气动性能及结构响应的探索。以南京航空航天大学自主研发的5 MW水平轴风力机体系为对象，以风?沙双向耦合算法为核心，基于CFD技术分别采用连续相和离散相模型进行典型风场和沙粒组合的同步迭代模拟。考虑叶片单个旋转周期内8个停机位置，对比分析了风力机塔架表面等效压力系数和气动载荷分布特性;结合有限元方法系统探讨了不同停机位置下风力机体系动力特性、风?沙致结构响应和屈曲稳定性能。研究表明：风?沙耦合环境下沙粒在风力机塔架迎风面底部产生的荷载效应最为显著，部分区域沙致压力系数可达0.55，沙粒冲击荷载与风荷载的比值最高可达23.7%;不同停机位下风荷载均对塔架结构响应起主导作用，但沙粒冲击作用在塔架底部的放大效应不可忽视。

关键词： 大型风力机体系; 风?沙耦合环境; 结构响应; CFD数值模拟; 停机位置

中图分类号： TK83    文献标志码： A    文章编号： 1004-4523（2021）01-0060-12

DOI：10.16385/j.cnki.issn.1004-4523.2021.01.007

引  言

中国西部风能资源充沛，然而此地区戈壁沙漠面积广大，沙粒受到近地边界层风场的驱动，在一定高度范围内随机运动^[1?2];同时，沙粒的附着也会影响风的湍流特征并对风力机体系产生附加冲击力，两者之间相互耦合并发生能量迁移^[3]，使得沙粒以较大速度冲击至风力机表面，形成附加荷载及风蚀效应^[4]。在沙尘暴极端天气条件下，风力机组一般处于停机状态，此时叶片不同停机位置将显著影响整个风力机塔架?叶片耦合体系的气动力分布和结构稳定性能^[5]。因此，对风?沙耦合环境下大型风力机体系气动性能及结构响应开展探索性研究具有重要的工程与理论意义。

大型风力机体系属于典型的风敏感结构，现有研究主要集中在流场作用机理、结构气动力分布^[6]、气弹耦合效应^[7]和体系风致响应特性^[8]等内容。针对风?沙环境下风力机的气动性能研究较少，文献[9]通过低角度加速冲蚀磨损试验研究了风沙环境下不同速度、冲角对风力机叶片的冲刷磨损行为特性，分析发现叶片涂层的冲蚀磨损量随着冲蚀速度的增大而增加;文献[10]采用虚拟风洞的方法研究了风沙环境下模型风轮的转矩特性，获得了风力机风轮转矩随沙尘颗粒粒径和体积分数变化的规律，分析发现，风轮转矩随沙尘体积分数的增大先增大后减小，随沙尘颗粒直径的增大，先增大后減小，最后趋于平稳;综上所述，风?沙环境对于风力机气动性能的影响不可忽略，尤其对于沙尘暴极端天气下大型风力机体系的风?沙耦合作用效应缺乏定量和定性的研究。此外，国内外风电机组设计规范和标准^[11?15]也尚未考虑此类极端环境下风力机的气动和结构性能评估。

鉴于此，本文以南京航空航天大学自主研发的5 MW水平轴三叶片风力机为研究对象，基于风?沙双向耦合算法，分别采用连续相和离散相模型进行典型风场和沙粒组合的同步迭代模拟。在此基础上主要研究：（1）不同停机位置下风力机塔架表面沙致压力系数和气动载荷分布特性;（2）不同停机位置下风力机体系动力特性、风?沙致结构响应和屈曲稳定性能。研究结论可为此类极端天气条件下大型风力机风?沙冲击荷载取值和结构设计提供参考。

1 风-沙双相耦合算法

1.1 沙尘暴等级划分与风沙流密度

中国西北地区沙尘暴天气频发，主要集中在春季4至5月，其发生时间集中、频率高、强度大。根据探空资料，沙尘暴发生时形成的沙尘壁高达300 m，沙尘暴影响高度在2100 m以内^[16?17]。沙尘暴天气根据最大风速和最小能见度划分为四个等级，如表1所示。本文选取中级沙尘暴天气典型风速18 m/s风场计算风速，双相耦合迭代算法中沙粒粒径取为1 mm^[18]。

2 数值模拟

2.1 工况设置

表2给出了南京航空航天大学自主研发的5 MW上风型水平轴三叶片风力机主要结构设计参数及模型示意图，其中塔架为通长变厚度结构，顶壁厚40 mm，底壁厚90 mm;叶片倾角为5°，切出风速25 m/s，各叶片之间成120°夹角，沿周向均匀分布，叶片长度为60 m，沿翼展各叶素截面的详细参数如表3所示;机舱尺寸为18 m×6 m×6 m，分别对应长、宽、高方向。根据以上设计参数依次建立塔架、叶片和机舱等部件，通过布尔运算形成大型风力机三维实体模型。

该风力机位于B类地貌，根据叶片与塔架的相对位置，考虑到三叶片体系旋转过程中存在的周期性，定义叶片与塔架完全重合时为工况1，顺时针旋转15°为一个工况共计8个。此外定义工况1中与塔架完全重合的叶片编号为Ⅰ，顺时针旋转分别为叶片Ⅱ和叶片Ⅲ。不同工况下风力机叶片单周期旋转方向如图1所示。

2.2 计算域设置

如图2所示，本文研究的风力机均处于0°来流风向角下。为保证风力机尾流能够充分发展，计算域尺寸设置为12D×5D×5D（流向X×展向Y×竖向Z，D为风轮直径），风力机置于距离计算域入口3D处。为兼顾计算效率与精度，同时考虑到叶片表面扭曲复杂，网格划分采用混合网格离散形式，将整个计算域划分为局部加密区域和外围区域。局部加密区域内含风力机模型，采用非结构化网格进行划分，外围区域形状规整，采用高质量的结构化网格进行划分。

为确保数值模拟结果的可信度，本文增加了网格无关性验证，表4给出了不同网格方案下网格质量和迎风面压力系数。由表4可知，随着网格总数的增加，网格质量逐渐增加，网格歪斜度和迎风面风压系数呈现逐渐减小的趋势，而840万网格数和1100万网格数的网格质量和计算结果无明显差异，综合计算精度和效率，本文选取840万网格总数的方案。计算域及具体的网格划分如图2所示。

计算域左侧和顶部边界条件为速度入口，右侧为压力出口，数值计算采用3D单精度、分离式求解器，流场流速为绝对速度，空气模型等效为理想不可压缩流体，计算域入口采用幂指数为0.15的风廓线模型，离地10 m高度处的风速设置18 m/s，采用SIMPLEC算法实现速度与压力之间的耦合，对流项求解格式为二阶，计算过程中设置了网格倾斜校正以提高混合网格计算效果，控制方程的能量计算残差设置为1×10^-6，最后初始化风场进行迭代计算。图3给出了平均风速、湍流度剖面模拟结果与理论值的对比曲线，结果表明平均风速和湍流度剖面均与理论值吻合良好。图4给出了8种工况下塔架30 m高度环向压力系数分布图。可知各工况平均风压系数分布曲线的负压极值点和分离点对应角度与规范曲线^[11]一致，迎风和侧风区域风压系数数值吻合较好，仅在背风区负压系数略大于规范值，对比结果验证了风场数值模拟的有效性。

2.3 沙粒冲击荷载

采用式（6）分别进行8种工况下风力机塔架、叶片沙荷载计算，图5?6分别给出了各叶片、塔架沙荷载。由图可知：风力机不同叶片沙荷载分布呈现不同的规律，叶片Ⅰ沙荷载随叶片高度增加先减小后增大;叶片Ⅱ沙荷载随高度增加先增加后减小，且在工况5（60°）达到极大值;叶片Ⅲ沙荷载随高度减小先减小后增大。

式中n为单位时间内撞击某区域的沙粒数量，F（τ）为单个沙粒对结构的冲击力，F_W为该区域风荷载。表5给出了沿塔架子午向高度风、沙荷载比值，可以发现：不同叶片停机位置对风力机塔架表面沙荷载分布影响较大，但均在（0?0.05）H高度范围内产生该工况下极大沙荷载，工况1，2沙荷载沿高度分布呈现较为明显的变化趋势，其随高度的变化先减小后增大再减小;其余工况沙荷载在塔架0.1H高度以上分布较为一致。

2.4 沙致压力系数

为定量比较不同停机位置下风力机塔筒表面压力分布，定义断面风压系数、沙致压力系数、等效压力系数，如下式所示：

式中C_pi为i点的平均风压系数;P_i为测点平均压力（Pa）;P_H為参考高度处远前方的静压（Pa）;ρ为空气密度，本文数值模拟过程中均取1.225 kg/m³;V_H为参考高度处远前方的平均风速;C_psi为i点的沙致压力系数;n为单位时间内撞击某区域的沙粒数量;F（τ）为单个沙粒对结构的冲击力;S为计算区域面积;C_pei为i点的等效压力系数。

其等效思路为：①将表面各监控点沙粒撞击荷载转化成沙致压强;②计算监控点沙致压强与对应参考高度处风压比值，即沙压系数;③将风压系数与转化后的沙致压力系数相加，得到等效压力系数。等效压力系数可作为等效目标，定量比较不同停机位下风力机塔架所受风?沙耦合作用大小。

考虑沙粒着点分布特性及叶片与塔架之间的气动干扰效应，选取塔架未干扰区段沙粒分布密度差异较大，干扰区段塔架与叶片重合面积不同的2个断面为典型断面。图7?9给出各工况下风力机塔架4个典型断面风压系数、沙致压力系数、等效压力系数对比曲线，分析可知：1）未干扰区（0.025H，0.30H）塔架断面呈现良好的对称性;显著干扰区（0.70H，0.90H）塔架断面工况1，2，3，6等效压力系数、风压系数分布曲线不再保持对称，塔架迎风面中心点处呈现负压。2）不同工况下塔架0.025H高度塔架断面迎风面两侧0°?60°范围内等效压力系数与风压系数差异明显，沙致压力系数最高可达0.549;其余典型断面（0.30H，0.70H，0.90H）等效压力系数与风压系数的数值和分布规律基本一致。

3 响应分析

3.1 有限元建模及动力特性

以图1显示的8个风力机不同停机位的计算工况为例，本文基于大型通用有限元分析软件ANSYS建立大型风力机塔架?叶片一体化有限元模型如图10所示，塔架和叶片采用Shell63单元，机舱采用Beam189单元，环基采用Solid65单元。建立模型后采用Block Lanczos方法求解风力机自振频率和振型^[22?23];同时，为确保有限元分析的可信度，本文增加了网格无关性验证，如表6给出了工况1下、不同网格方案下风力机基频。由表可知，随着网格总数的增加，风力机基频呈现逐渐减小的趋势，而4680网格数和6240网格数的计算结果无明显差异，综合计算精度和效率，本文选取网格总数为4680的方案。

图11给出了不同停机位置下风力机模型前100阶模态固有频率分布曲线及工况2?8相较于工况1风力机有限元模型基频的增/减量，可以看出大型风力机叶片?塔架耦合模型基频均较小，仅为0.138 Hz左右，且各模态之间频率间隔很小。叶片不同停机位置主要影响风力机体系的低阶固有频率，各工况低阶频率出现较小的差别，高阶频率基本一致。可见停机状态下叶片不同位置对风力机体系的频率和振型影响微弱，低阶振型主要以叶片带动机舱及塔架进行前后挥舞和左右摆动，高阶模态出现塔架及叶片本身的结构变形和失稳形态，如图12所示。

3.2 塔架响应

图13?14给出了不同工况下塔架风致径向位移响应、风?沙致径向位移响应分布图。由图可得：叶片不同停机位置塔架主要影响塔架中上部位移，沙粒冲击作用则主要影响塔架中下部位移。径向位移随着塔架高度的增加而逐渐增大，最大正负位移均出现在塔顶处0°和180°位置;不同停机位置下塔架位移极值在迎风面和背风面的分布范围不同，极值区域随叶片与塔架位置相对干涉作用增大而增大;工况1叶片完全遮挡塔架，塔架顶部位移出现“逆向效应”，沙粒冲击作用放大，塔架顶部风?沙致位移为负。

为对比叶片位置停机位置、沙粒冲击作用引起塔架底部弯矩值的差异，分别以工况1在净风、风?沙共同作用下响应值作为初始状态，将各工况下塔底内力响应值与之作差，得到塔架风致内力响应特征值、风?沙致内力响应特征值;再将风致、风?沙致内力响应值作差，得到沙致响应特征值。

如图15所示：塔底径向弯矩只在0°，10°，70°以及320°位置出现差异，工况2在0°位置数值相差最大达到146.82 N·m;风致塔底径向弯矩与风?沙致环向弯矩在迎风面差异显著，沙致塔底径向弯矩最高达12.78 N·m;背风面环向弯矩呈左右对称，风致响应与风?沙致响应分布基本一致，沙粒主要冲击塔架迎风面左右60°，沙致弯矩响应极值均出现在0°位置，背风面及侧风面可忽略沙粒冲击作用。

3.3 屈曲稳定性分析

图16给出了不同工况下风力机塔架?叶片耦合体系屈曲位移和临界风速对比示意图。对比可知：1）叶片停机位置对整体结构屈曲失稳时最大位移影响显著，随叶片对塔架遮挡面积的减小最大位移逐渐增大，在工况5时达到峰值;2）风与风?沙耦合作用下临界风速、屈曲最大位移分布曲線规律基本一致，风?沙耦合作用使得临界风速略微增大;工况3对应的临界风速最低，工况5对应的临界风速最高。

4 结  论

本文系统探讨了风?沙耦合环境下大型风力机体系气动力和结构响应特性，主要涉及风?沙双向耦合数值模拟、风荷载分布与风致响应特性、沙荷载分布与沙致响应特性、风?沙致等效荷载及其效应等。得到如下主要结论：

1）不同叶片停机位置对风力机塔架表面沙荷载分布影响较大，但均在（0?0.05）H高度范围内产生该工况下极大沙荷载，工况6（75°）下该高度范围内沙荷载与风荷载比值达23.70%;且随高度增加，沙荷载与风荷载比值逐渐减小。

2）不同工况下沙粒撞击位置主要集中在塔架0.6H（未干扰段）高度以下，集中于迎风区域两侧各60°范围，沙致压力系数最大值为0.549，发生在工况4的塔架底部迎风面。考虑叶片停机位置影响时，叶片遮挡效应对塔架中上部位移影响显著，沙粒冲击作用则加大塔架中下部位移。

3）塔架径向位移分布规律差异较大，工况1（叶片完全遮挡）径向位移对沙粒冲击作用敏感，径向位移随着塔架高度的增加而逐渐增大，最大正负位移均出现在塔顶处0°和180°位置。沙致弯矩响应极值均出现在0°位置，总体分布大致呈左右对称，背风面及侧风面可忽略沙粒冲击作用。

参考文献：

[1]        Bagnold R A. The movement of desert sand[J]. Geographical Journal， 1935，85（4）：342-365.

[2]        Bagnold R A. The transport of sand by wind[J]. Geographical Journal， 1937，89（5）：409-438.

[3]        Cheng J， Jiang F， Xue C， et al. Computational method for maximum sediment discharge and sand-carrying wind load in the prevention and treatment of wind drift sand for railway in strong wind area[J]. China Railway Science， 2012，33（1）：1-5.

[4]        蒋富强，李  荧，李凯崇，等.兰新铁路百里风区风沙流结构特性研究[J].铁道学报，2010，32（3）：105-110.

Jiang Fuqiang， Li Ying， Li Chongkai， et al. Study on structural characteristics of Gobi wind sand flow in 100 km wind area along Lan-Xin Railway[J]. Journal of the China Railway Society， 2010，32（3）：105-110.

[5]        楼文娟，余  江，潘小涛.风力机叶片挥舞摆振气弹失稳分析[J].工程力学，2015，32（11）：236-242.

Lou Wenjun， Yu Jiang， Pan Xiaotao. Calculating for aerodynamic stability response of wind turbine blade in flapwise and edgewise direction?[J]. Engineering Mechanics， 2015，32（11）：236-242.

[6]        柯世堂，余  玮，王同光.停机状态叶片位置对风力机体系气动性能影响[J].浙江大学学报（工学版），2016，50（7）：1230-1238.

Ke Shitang， Yu Wei， Wang Tongguang. Impact for blade position on aerodynamic performance of wind turbine system under stopped status[J]. Journal of Zhejiang University （Engineering Science）， 2016，50（7）：1230-1238.

[7]        李德源，莫文威，夏鸿建，等.水平轴风力机柔性叶片气弹耦合分析[J].太阳能学报，2015，36（3）：734-742.

Li Deyuan， Mo Wenwei， Xia Hongjian， et al. Aeroelastic coupling analysis of flexible blades for horizontal axis wind turbines[J]. Chinese Journal of Solar Energy， 2015，36（3）：734-742.

[8]        柯世堂，余  瑋，王同光.基于大涡模拟考虑叶片停机位置大型风力机风振响应分析[J].振动与冲击，2017，36（7）：92-98.

Ke Shitang， Yu Wei， Wang Tongguang. Analysis of wind-induced vibration response of large wind turbine based on large eddy simulation[J]. Journal of Vibration and Shock， 2017，36（7）：92-98.

[9]        张  永，黄  超，刘  召，等.挟沙风作用下风力机叶片涂层冲蚀特性研究[J].材料导报， 2016，30（10）：95-99.

Zhang Yong， Huang Chao， Liu Zhao， et al. Study on erosion characteristics of wind turbine blade coating under the effect of wind and sand[J]. Material Review， 2016，30（10）：95-99.

[10]      李仁年，范向增，李德顺，等.风沙环境下水平轴风力机的转矩特性[J].兰州理工大学学报，2015，41（4）：55-59.

Li Rennian， Fan Xiangzeng， Li Deshun， et al. Torque characteristics of horizontal-axis wind turbine in sand-blowing environment[J]. Journal of Lanzhou University of Technology， 2015， 41（4）： 55-59.

[11]      中华人民共和国住房和城乡建设部.GB50009-2012，建筑结构荷载规范[S]. 北京： 中国建筑工业出版社， 2012.

MOHURD. GB50009-2012， Building structure load specification [S]. Beijing： China Building Industry Press， 2012.

[12]      中华人民共和国国家发展和改革委员会.DL/T 5383-2007， 风力机发电场设计技术规范[S].北京：中国电力出版社， 2007.

National Development and Reform Commision. DL/T 5383-2007， Technical specification of wind power plant design[S]. Beijing： China Electric Power Press， 2007.

[13]      中华人民共和国国家质量监督检验检疫总局，中国国家标准化管理委员会.GB/T 20319-2006， 风力发电机组验收规范[S].北京：中国机械工业出版社， 2006.

General Administration of Quality Supervision， Inspection and Quarantine of the People's Republic of China， Standardization Administration of the People's Republic of China. GB/T 20319-2006. Code for acceptance of wind turbines generator systems[S]. Beijing： Chinese Machinery Industry Press， 2006.

[14]      IEC 61400-21：2001. Wind turbine generator systems-Part21： Measurement and assessment of power quality characteristics of grid connected wind turbines[S]. International Elector-technical Commission， 2001.

[15]      IEC 61400-23：2001. Wind turbine generator systems-Part23：Full-scale structural testing of rotor blades[S]. International Elector-technical Commission， 2001.

[16]      唐国利，巢清尘.近48年中国沙尘暴的时空分布特征及其变化[J].应用气象学报， 2005，16（b03）：128-132.

Tang Guoli， Chao Qingchen. The temporal and spatial distribution characteristics and variation of sandstorm in China in recent 48 years[J]. Journal of Applied Meteorology， 2005，16（b03）：128-132.

[17]      戴雪荣，师育新，薛  滨.兰州现代特大尘暴沉积物粒度特征及其意义[J].兰州大学学报（自科版）， 1995，31（4）：168-174.

Dai Xuerong， Shi Yuxin， Xue Bin. Particle size characteristics and its significance of modern duststorms in lanzhou[J]. Journal of Lanzhou University（Natural Science）， 1995，31（4）：168-174.

[18]      董依帆，柯世堂，杨  青.考虑风-沙双向耦合作用大型风力机体系气动力分布研究[J].振动与冲击， 2019，38（16）：84-92.

Dong Yifan， Ke Shitang， Yang Qing. Wind-sand movement characteristics and aerodynamic distribution of large wind turbine systems considering bidirectional coupling effect?[J]. Journal of Vibration and Shock， 2019，38（16）：84-92.

[19]      李天涯.冲量及动量定理的示例探讨[J].物理教学探讨， 2014，32（12）：36-36.

Li Tianya. An example of the impulse and momentum theorem[J]. Discussion on Physics Teaching， 2014，32（12）：36-36.

[20]      楊庆峰.用牛顿第二定律解决问题的几种方法[J].数理化解题研究，2018，（13）： 78-79.

[21]      罗贝格.弹性及弹-塑性介质中的冲击波[M].北京：科学出版社， 1965.

Robberg. Shock Waves in Elastic and Elastic-Plastic Media [M]. Beijing： Science Press， 1965.

[22]      Ke Shitang， Wang Tongguang， Ge Yaojun， et al. Wind-induced responses and equivalent static wind loads of tower-blade coupled large wind turbine system[J]. Structural Engineering and Mechanics， 2014，52（3）：485-505.

[23]      邓  露，肖志颖，黄民希，等.考虑流固耦合的近海风机动力响应数值计算[J].湖南大学学报（自然科学版）， 2015，42（7）：1-8.

DENG Lu， XIAO Zhiying， HUANG Mingxi， et al. Numerical simulation of dynamic response for offshore wind turbines including fluid-structure interaction[J]. Journal of Hunan University （Natural Sciences）， 2015，42（7）：1-8.

Abstract： The position of the blades in a strong wind-off state will significantly affect the flow around wind turbine towers and stability of them. Especially in the extreme weather conditions of sandstorms， the adhesion of sand will also affect the turbulence characteristics of the wind and generate additional impact force on the tower. The existing researches lack the exploration of wind-sand coupling effect to explore the aerodynamic performance and structural response of large-scale wind turbines. Taking the 5 MW horizontal-axis wind turbine system independently developed by Nanjing University of Aeronautics and Astronautics as the object and wind-sand-wave two-way coupling algorithm as the core， the simultaneous iterative simulation of the typical wind field and sand combination is performed on the continuous phase and discrete phase model respectively by CFD technology. Considering the eight stop positions in a single rotation cycle of the blade， the equivalent pressure coefficient and the aerodynamic load distribution characteristics are compared and analyzed. Then combined with the finite element method， the dynamic characteristics， wind-sand induced structural response and buckling stability of the wind turbine system at different stopping positions are discussed. The results show that the load effect of sand particles on the windward surface of wind turbine tower is the most significant under the wind-sand coupling environment. The pressure coefficient of sand in some regions can reach 0.55， and the ratio of sand impact load to wind load can reach 23.7%. The wind load under different stand positions plays a dominant role in the response of the tower structure， but the effect of sand impact on the bottom of the tower cannot be ignored.

Key words： large wind turbine system; wind-sand coupling environment; structural response; CFD numerical simulation; stop position

作者簡介： 柯世堂（1982-），男，博士，教授。电话：13621581707; E-mail： keshitang@163.com