考虑风-沙双向耦合作用大型风力机体系气动力分布研究

董依帆， 柯世堂， 杨 青

(南京航空航天大学 土木工程系，南京 210016)

近年来随着国内能源结构转型步伐的加快，风力发电设施陆续兴建。其中我国西部风能资源充沛[1]，为更好的捕捉风能，风力机必然朝着高大化方向发展。然而此地区戈壁沙漠面积广大，沙粒受到近地边界层风场的驱动，在一定高度范围内随机运动，两者之间相互耦合并发生能量迁移，使得沙粒以较大速度冲击至风力机表面，形成附加荷载及风蚀效应，对风力机的安全性能和使用寿命产生显著影响[2-3]。因此，开展气-固两相流理论研究，探讨不同风速和沙粒粒径组合对风力机气动载荷的作用机理，具有重要的工程意义和理论价值。

当前国内外对风-沙共同作用的研究多集中在风沙现象和风沙运动本身，包括沙粒运动力学过程[4]、输沙量的垂直分布[5]、输沙率模型[6]、临界起沙风速[7]和风速廓线规律[8]。相较于风荷载作用，风-沙共同作用下建筑物或构筑物荷载效应更加显著[9]。文献[10]通过风洞测力试验，研究了不同沙浓度和风速条件下风沙对低矮建筑整体受力影响，主要考察风沙环境中建筑物表面的风压特性；文献[11]模拟了大气边界层中单体及两相邻对称高层建筑的风沙流场，并指出随着沙粒相体积分数的增大，风压系数不断增大；文献[12]对风沙环境下水平轴风力机进行数值模拟，研究了在某一来流风速和不同沙粒直径组合工况下的风力机风轮转矩。已有研究均未基于风力机整机和沙粒真实环境分布特性展开多参数条件下大型风力机风沙特性的系统研究。

鉴此，本文以南京航空航天大学自主研发的5 MW水平轴三叶片风力机[13]为研究对象，采用改进的k-ε湍流模型和粉尘释放模型(Dust Production Model，DPM)，设置不同风速和不同沙粒粒径形成9种组合计算工况，对处于典型停机状态的塔架-叶片耦合体系展开数值模拟，以获取沙粒在风力机不同部位的分布密度及撞击速度，并对比分析风沙共同作用下风力机体系绕流特性、表面压力分布随风沙流速度和沙粒浓度的变化规律。

1 风-沙双相耦合算法

1.1 沙尘暴等级划分

我国西北地区沙尘暴天气频发，主要集中在春季4～5月，其发生时间集中、频率高、强度大。根据探空资料，沙尘暴发生时形成的沙尘壁高达300 m，沙尘暴影响高度在2 100 m以内[14-15]。沙尘暴天气根据最大风速与最小能见度划分为4个等级，如表1所示。本文选取中、强、特强沙尘暴天气典型风速18 m/s，22 m/s，25 m/s作为风场计算风速。

表1 沙尘暴天气强度划分标准

1.2 风沙流密度

风沙流密度是指空气运动速度在大于起沙风速情况下，单位体积空气气流中所含沙的质量，单位为kg/m3。在近地层中，风受地面摩擦阻力的影响而速度降低，一般而言，摩擦力随着高度的增加而减小，故风速随着高度的增加而增大，风沙流密度与高度、风力等均有密切关系。

据文献[16]统计，当风速在40 m/s以下时，风沙流密度的数量级一般约为10-5，且随高度的增加而下降，基本以3 m高度作为分界点。在3 m高度以下，风沙流密度随高度的增加而大幅下降，在3 m高度以上风沙流密度随高度增加逐步下降，且沙粒运动以飞扬为主。基于工程安全性考虑，采取极端工况进行数值模拟，故风场入口处3 m以下高度风沙流密度取为10-4kg/m3，速度与风速相同；3 m以上高度风沙流密度取为10-5kg/m3，速度为0。

1.3 计算模型与颗粒运动方程

标准k-ε模型在科学研究及工程实际中得到了广泛的应用，但用于强旋流、弯曲壁面流动或弯曲流线流动时，会产生一定的失真。为弥补其缺陷，本文采用Realizablek-ε模型，Realizablek-ε模型中k与ε的模数化输运方程组为

(1)

(2)

与标准k-ε模型比较发现，Realizablek-ε模型主要变化是：湍流黏度计算公式发生了变化，引入了与旋转和曲率有关的内容；ε方程发生了很大变化，方程中的产生项不再包含有k方程中的产生项Gk；ε的倒数第二项不具有任何奇异性，即使值很小或为零，也不会为零。

相比于视为连续相的空气来流，沙粒在空气中所占比例较小，即使在特强沙尘暴(风速≥25 m/s)时，沙粒占空气体积分数仍远小于10%。基于此类实际情况，DPM模型被广泛应用于多相流联合研究。本文在风场计算稳定后，将离散相模型作为第二相，介入至连续相中展开风-沙双向耦合运算。沙粒在风场中运动平衡方程为

(3)

(4)

式中：μ为流体黏性系数；dp为颗粒直径；Re为相对雷诺数； 可表示为

(5)

考虑沙粒离散相影响后，风连续相基本控制方程可表示为

(6)

(7)

(8)

式中：I为单位张量； 等式右边第二项为体积膨胀作用。

1.4 颗粒碰撞模型

沙粒冲击到风力机表面过程服从动量守恒定律，求解冲击力的关键在于碰撞时间。计算中忽略沙粒在冲击过程中可能发生的破裂现象，认为沙粒与结构间相互作用遵循牛顿第二定律并假设反弹后的速度与撞击前的速度一致。由动量定理

(9)

式中：f(t)为单个沙粒冲击力矢量， N；vs为沙粒速度矢量。

沙粒在单位时间内对结构的冲击力F(τ)为

(10)

将撞击风力机的沙粒近似看作球体，则

(11)

由于沙粒直径在1 mm及以下，撞击前水平末速度相对较大，可将风力机表面简化为无限大平面[17]，将碰撞时间τ取为

(12)

(13)

式中：m为沙粒的质量；φ为速度恢复系数取1；k1，k2取值与沙粒、风力机材料杨氏模量、泊松比有关。

则沙粒对结构的冲击力可简化为

(14)

2 数值模拟

2.1 工程简介

表2给出了该5 MW风力机的主要结构尺寸及示意图，叶片长度为60 m，机舱尺寸18 m×6 m×6 m，各叶片沿周向成120°夹角分布，塔架通长变厚，高为124 m，塔底半径3.5 m，塔顶半径3.0 m，塔顶厚度90 mm，塔底厚度200 mm。

表2 5 MW风力发电结构体系主要参数列表

2.2 工况设置

该风力机位于B类地貌，综合考虑我国沙尘暴天气等级划分，对比研究3种风速和3种沙粒粒径组合下极端工况对风力机整机的作用规律。其中，小风、中风和大风分别以沙尘暴等级中、强和特强的基本风速进行划分；沙粒粒径分别取为最大统计概率区间内的0.10 mm，0.25 mm以及极端工况下的1 mm，3种风速和沙粒粒径组合共计9种工况，如表3所示。

表3 风力机风-沙耦合计算工况划分表

2.3 计算域设置

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

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

表4 不同网格方案下网格质量和迎风面压力系数

图1 计算域及网格划分示意图Fig.1 Schematic diagram of computational domain and encrypted mesh

计算域左侧和顶部边界条件为速度入口，右侧为压力出口，数值计算采用3D单精度、分离式求解器，流场流速为绝对速度，空气模型等效为理想不可压缩流体，计算域入口采用幂指数为0.15的风廓线模型，离地10 m高度处的风速分别设置为“2.2”节中3种基准风速。流场求解采用SIMPLEC算法实现速度与压力之间的耦合，对流项求解格式为二阶，计算过程中设置了网格倾斜校正以提高混合网格计算效果，控制方程的计算残差设置为1×10-6，最后初始化风场进行迭代计算。图2给出了平均风速、湍流度剖面模拟结果与理论值的对比曲线，结果表明平均风速和湍流度剖面均与理论值吻合良好，风场模拟标准满足工程要求。

图2 不同风速下风场速度及湍流度剖面示意图Fig.2 Velocity and turbulence profiles under different wind speeds

2.4 有效性验证

图3给出了3种风速下塔架30 m高度环向压力系数分布图。对比可知各工况平均风压系数分布曲线的负压极值点和分离点对应角度与规范曲线[18]一致，迎风和侧风区域风压系数数值吻合较好，仅在背风区负压系数略大于规范值，对比结果验证了风场数值模拟的有效性。

图3 3种风速下风力机塔架30 m高度压力系数与规范对比曲线Fig.3 Pressure coefficients on the 30 m height tower under three wind speeds and domestic codes

3 结果分析

为研究风-沙荷载特征值随环境参数变化的作用规律，首先分析风-沙双向耦合作用下风场分布特性、沙粒着点及其分布规律；在此基础上，通过计算沙粒冲击荷载、等效外压系数定量分析不同工况下沙致效应差异。

3.1 风场分布特性

图4给出了不同风速下大型风力机体系风速流线图。由图4可知，由于叶片的遮挡效应，塔架迎风面压力明显减小，叶片与塔架重合部位内侧出现附着于塔架迎风面的自由剪切层，并在塔架背风面形成成熟的旋涡脱落。

图4 风力机体系风速流线图Fig.4 Wind speed streamline of wind turbine

图5和图6分别给出了未干扰及干扰区段塔架涡量云图。对比分析可知：未干扰区段气流在风力机塔架迎风面前缘发生流动分离，并在侧边出现加速效应；而在干扰区段由于叶片对塔架的干扰作用，塔架和叶片之间的区域随着风速的增加出现明显的能量积聚，持续发展后在背风面形成尾流涡旋以及回流，且随风速的增加而加强。

图5 风力机塔架未干扰段湍动能图Fig.5 Wind turbine tower undisturbed turbulent kinetic energy

图6 风力机塔架干扰段湍动能图Fig.6 Turbulence kinetic energy of wind turbine tower interference section

3.2 风沙运动特征

风沙流场中沙粒冲击风力机表面，易产生极端荷载效应。为清晰展示沙粒在风力机不同部位的分布密度，基于颗粒水平向速度对沙粒轨迹进行追踪，同时定义单位时间内冲击塔架沙粒的位置分布为沙粒着点，图7给出了9种工况风沙场中沙粒着点图，并对沙粒的密集程度进行了等比例粗化处理。

由图7可知，由于叶片的对塔架的遮挡效应，各工况沙粒撞击位置多集中分布风力机塔架迎风区域0～0.6H高度范围内，H为风力机塔架高度。受气流旋涡驱动作用，塔架遮挡区和背风区壁面仅有少量沙粒附着；塔架、叶片收集到的沙粒均以工况7最多，且随风速的增加而增加，随沙粒粒径的增大逐渐变小。

图7 风沙场中沙粒着点示意图Fig.7 Sand-reflect motion trajectory in sand and wind field

3.3 风沙荷载分布

3.3.1 沙粒空间分布规律

图8和图9分别给出了9种工况下风力机塔架和叶片收集到的沙粒数量、撞击速度及速度占有率(各速度分布区间内沙粒数量与总数量的比值)对比曲线。由图8和图9对比可知：

(1) 当沙粒粒径较小时，风力机塔架收集到的沙粒数量随风速的增大而增加，当粒径达到中等以上时，风力机塔架收集到的沙粒数量随风速的增大而减少；然而，粒径和风速对风力机叶片的沙粒着点影响并不显著，这是由于叶片位置较高，处于风沙流密度影响较小区域。

图8 风力机塔架沙粒数量与水平末速度分布曲线Fig.8 Wind turbine tower sand number and horizontal velocity distribution curve

图9 风力机叶片沙粒数量与水平末速度分布曲线Fig.9 Wind turbine blade grit and horizontal velocity distribution curve

(2) 撞击风力机塔架的沙粒，随着风速增大，小粒径和大粒径沙粒的平均速度均增大，而中等粒径的沙粒平均速度随风速的增大而减小；风力机叶片上追踪的沙粒平均速度特性也呈现相似规律，不同在于风力机叶片上的大粒径沙粒平均速度随风速的增大而减小。

(3) 风力机塔架、叶片收集到的沙粒水平速度占有率分布规律基本一致，速度占有率随着水平末速度(顺风向平均速度)的增大先增加后减小。此外，收集到的沙粒平均水平速度均远小于基准风速，且沙粒水平撞击速度均随着沙粒直径的增加而增大。

3.3.2 叶片、塔架沙粒冲击荷载

定义与塔架重合的叶片为A、迎风面左侧叶片为B、右侧叶片为C。采用式(12)分别进行9种工况下风力机塔架、叶片沙荷载计算，同时给出了塔架、叶片不同高度(叶展长度)范围沙荷载以及风力机不同高度范围沙荷载与该区域风荷载的比值，如图10、图11和表5所示。经对比发现：

(1) 大粒径沙粒与不同风速的组合工况(工况3、工况6、工况9)在风力机表面产生显著荷载，由于位置和形状的不同，风力机塔架、叶片沙荷载分布呈现不同的规律；

(2) 叶片A、叶片B、叶片C沙荷载均沿叶展方向呈减小趋势，由于所处高度不同，风沙流密度影响程度不同，叶片A承受最大沙荷载约为叶片B、叶片C的两倍；

(3) 不同工况对风力机塔架表面沙荷载分布影响较大，但均在0～0.05H高度范围内产生该工况下极大沙荷载，工况3下该高度范围内沙荷载与风荷载比值达22.57%；

(4) 显著工况(工况3、工况6、工况9)下风力机塔架极大沙荷载随风速的变大而增大，沙荷载沿塔架高度先减小后增大再减小，0.2H～0.6H高度范围内沙荷载极大值也随风速的增大而增大。由于叶片A的遮挡效应，塔架0.5H高度以上沙荷载减小。

图10 各工况风力机叶片沙荷载示意图Fig.10 Sand load distribution on different height of wind turbine blades

图11 风力机塔筒沿子午线高度荷载示意图Fig.11 Comparison curves of equivalent pressure in typical meridian of wind turbine tower

表5 工况1～工况9塔架不同高度范围风沙荷载比值列表

3.3.3 沙压系数

沙致压力系数(简称沙压系数)计算方法为：①将表面各监控点沙粒撞击荷载转化成沙致压强；②计算监控点沙致压强与对应参考高度处风压比值，即沙压系数。图12给出了各工况下风力机塔架沙压系数三维分布示意图。由图12可知：

(1) 各工况风力机塔架沙压系数均主要集中于迎风面两侧各60°范围内，其余范围数值基本为0，各工况沙压系数最大值均发生在0～0.05H高度范围内，各工况沙压系数最大值为0.517；

(2) 沙压系数分布呈现明显三维效应，其随塔架高度增加而迅速减小，塔架与叶片重合部分沙压系数趋于0；

(3) 不同组合工况下沙压系数大小差异明显，且沙粒粒径的增大而显著增大，此外，其最大值随风速的增大呈减小趋势，而其总体分布更趋于平均化。

图12 各工况风力机塔架沙致压力系数三维分布示意图Fig.12 Three-dimensional distribution of sand-reflect and sand pressure coefficient on surface of wind turbine tower

3.3.4 风-沙共同作用等效压力系数

为定量比较不同组合工况下的塔架风-沙致压力分布，将沙压系数与风压系数矢量加和，定义等效压力系数作为沙致效应的衡量参数。

考虑沙粒着点分布特性及叶片与塔筒之间的气动干扰效应，选取塔架未干扰区段沙粒分布密度差异较大，干扰区段塔架与叶片重合面积不同的两个断面为典型断面。图13给出各工况下风力机塔架4个典型断面等效压力系数对比曲线，分析可知：

图13 风力机塔架典型断面环向等效压力系数对比曲线Fig.13 Comparison curves of equivalent pressure coefficient in typical section of wind turbine tower

(1) 不同工况下相同高度截面处等效压力系数分布规律及数值基本一致，位于未干扰区(0.025H，0.30H)塔筒断面环向风压系数呈现良好的对称性；而位于显著干扰区(0.70H，0.90H)塔筒断面，其环向表面压力系数分布曲线不再保持对称，塔筒迎风面中心点处呈现负压；

(2) 同一工况下4个典型断面等效压力系数数值均略有差异，其中最大负压随着高度的增加先增大后减小，在叶片遮挡区域附近约为-1.5；

(3) 不同工况下0.025H高度塔筒断面迎风面两侧0～30°内压力系数差异明显，工况3、工况6、工况9大于其余工况下等效压力系数，且最大压力系数随风速的增大而减小，背风面呈负压；

(4) 位于干扰区段塔架断面，随高度增加迎风面、侧风面和背风面负压分别呈增大、减小和增大的趋势。

图14给出了风力机塔架0°，90°，180°及270°四条典型子午线的等效压力系数对比曲线，分析可得：

(1) 风力机塔架90°与270°子午线等效压力系数随高度先增大后减小再增大，均在0.5H，0.8H高度处出现拐点；0°和180°子午线迎风面等效压力系数发生数值分离多集中在沙粒着点数目较多的0.5H高度以下的未干扰区段，0.1H～0.6H高度范围内背风面等效压力系数出现明显差异；

(2) 90°和270°子午线等效压力系数最大值均产生在0.5H高度处，最大值约为为-1.890；受风沙流影响，0°子午线最大等效压力系数发生0.005H高度处，约为1.092；270°子午线方向最大等效压力系数发生在塔架顶部附近，最大值为-0.282；

(3) 0°子午线底部等效压力系数受风沙流影响显著，大粒径沙粒对塔架的正面冲击作用加强，沙粒着点随风速的增大而减少，故等效压力系数随着风速的增大而减小。

图14 风力机塔筒典型子午线等效压力系数对比曲线Fig.14 Comparison curves of equivalent pressure coefficient in typical meridians of wind turbine tower

4 结 论

(1) 仅考虑单向风作用时风力机体系风场和风压系数均呈现明显的三维效应，考虑风-沙双向耦合后沙粒主要冲击塔架迎风面中下部，仅有少量沙粒随气流进入塔筒背风面；沙粒水平方向作用力与分布范围大小受风速和粒径的影响，在塔架和叶片等不同部位呈现不同的规律。

(2) 风力机塔架、叶片沙粒冲击数量随风速的增大而增加，随粒径的增大而减少，沙粒统计数量最多的为工况7(特强沙尘暴发生风速+典型工况下较小沙粒粒径)。塔架上的沙粒速度占有率最大值随风速、粒径的增大而降低，速度分布更加均匀，而冲击叶片沙粒速度占有率分布规律不受风速、粒径等因素干扰，主要表现为峰值右移。

(3) 不同工况下沙粒撞击位置主要集中在塔架0.6H高度以下，受两侧流动分离涡裹挟，集中于迎风区域两侧各60°范围，且大粒径沙粒冲击塔架表面更易形成显著荷载，沙致压力系数最大值为0.517，发生在工况3(中等沙尘暴发生风速+极端工况下沙粒粒径)的塔架底部迎风面。

(4) 不同工况下风力机塔架最大负压值随着高度的增加先增大后减小，背风区域负压值随着高度的增加先减小后增大，在叶片遮挡区域附近约为-1.5。各子午线等效压力系数随着高度的增加变化趋势不同，0°子午线底部等效压力系数受风沙流影响显著，大粒径沙粒对塔架的正面冲击作用加强，最大等效压力系数发生0.005H高度处，数值为1.092。