风力发电结构在风沙荷载激励下的动力响应分析

李万润, 黄 杰, 丁明轩, 李德强, 王雪平,3, 杜永峰,2,3

(1.兰州理工大学 防震减灾研究所, 甘肃 兰州 730050; 2. 兰州理工大学 甘肃省土木工程减震隔震国际科技合作基地, 甘肃 兰州 730050; 3.兰州理工大学西部土木工程防灾减灾教育部工程研究中心, 甘肃 兰州 730050; 4. 中国华电集团有限公司甘肃公司)

当前我国正处在工业化的重要阶段，高耗能产业占比大，能源短缺，环境保护形势严峻，节能减排任务艰巨.大力发展风电产业能有效缓解能源短缺和环境污染等问题.风能是一种洁净的可再生能源，同时也是当前最廉价、技术最成熟的可再生资源.我国西北地区因其独特的地理优势，拥有得天独厚的风能资源，风力发电产业得到大力的发展，兆瓦级风力发电机得到广泛应用.由于我国西北地区地表干旱，植被稀少，土地沙化严重，大风频繁，风力强劲，极易形成沙尘暴天气.沙尘暴对建(构)筑物侵蚀严重，其荷载效应明显大于净风荷载作用下[1].李秋胜等[2]利用北京气象塔对沙尘暴天气下城市中心的边界层风剖面进行现场实测，基于15个高度的风速和风向实测数据分析了沙尘暴天气下的平均风和湍流强度随高度的变化规律、大气稳定度对风剖面的影响、零平面位移和地面粗糙度的确定等城市边界层风特性.蒋国富等[3]通过对兰新铁路百里风区戈壁风沙流进行现场观测，提出了风沙流密度的概念，认为建(构)筑物在风沙流作用下的压力由风压和砂颗粒对建(构)筑物的冲击两部分组成.目前，国内外对沙尘暴的研究主要集中于风沙气候现象及风沙运动本身，而对沙尘暴作用于建筑物的风沙荷载、风沙流对建筑物环境影响以及将风沙的运动及其对建筑物的作用效应相结合的研究还很少.

风力发电结构的动力响应分析大多只针对风荷载作用，而关于风沙荷载对风力发电结构的动力响应分析及对其造成的危害缺乏深入的研究，针对风沙荷载的破坏程度和等级尚无定论.柯世堂等[4]以5 MW风力发电机为例采用谐波叠加法和叶素-动量理论生成风力发电机的气动荷载，并对风力机塔架-叶片耦合结构进行了非线性时域分析.赵荣珍等[5]根据多自由度模态分析理论，对水平轴风力发电机塔架的振动模态进行数学建模和仿真.章子华等[6]采用不随高度变化的台风脉动风功率谱，基于线性滤波法与竖向相关性简化表达式模拟了某沿海风电场台风风场，通过建立风机结构-桩基础耦合有限元模型，计算了风力发电结构在上、下限台风风速条件下动力时程，分析了极端风况下风力发电机主体结构可能的破坏模式.目前风荷载的数值计算已经从基础研究进入应用阶段，但关于风力发电结构在风沙荷载作用下的数值计算还很少有人研究.

风力发电结构对风荷载具有极强的敏感性，若在风中增加沙颗粒(例如沙尘暴)，则二者的耦合作用对风力发电结构产生的影响也将逐渐增大.因此本文基于动量守恒定律和风沙流密度理论，建立了风沙荷载力学计算模型，通过MATLAB软件分别对风和风沙荷载进行模拟，并根据西北某风场2 MW三桨叶水平轴风力发电机建立三种不同有限元模型，采用ANSYS软件对不同荷载作用下的风力发电结构动力响应差异进行分析，最终求解出不同等级沙尘暴对风力发电结构的影响.

1 风沙荷载

1.1 砂粒荷载

1.1.1砂粒荷载分类

砂粒荷载的分类比较复杂，本文以沙尘浓度作为砂粒荷载分类的主要依据.沙尘浓度是评估沙尘暴强度的一个重要因素，也是衡量风沙流所能携带砂粒量的重要参数，根据文献和各地观测值的综合分析数据可知，不同强度等级的沙尘天气下，沙尘颗粒浓度的综合值是不同[7]，具体参数如表1所示.

表1 沙尘暴颗粒浓度计算参数

1.1.2砂粒荷载计算

根据动量守恒定律和已有试验研究结论，沙尘暴灾害发生时，单个沙尘颗粒撞击挡风结构的撞击力计算公式为

(1)

式中，Fs(τ)为τ时间内沙尘颗粒对挡风结构的撞击力；m为τ时间内撞击挡风结构的沙尘颗粒的质量；v1、v2分别为强风流中单位体积沙尘颗粒撞击结构前、后的速度；t为时间.

为简化计算，假设沙尘颗粒撞击挡风结构为弹性碰撞，反弹后的速度与撞击前的速度一致(提高了建(构)筑物的安全储备，偏安全)，即：

v1=-v2=vp

(2)

式中:vp为强风流中单位体积砂砾颗粒的平均飞跃速度.

则式(1)可转化为

(3)

(4)

式中，mi(i=1,2,…,N)为砂粒i的质量；ρi为大于起沙风速下某一风速对应的风沙流密度.

假设单位体积内所有沙尘颗粒都作用在风力发电机上，则某一直径的沙尘颗粒荷载作用力之和为

Fi=F(τ)ni

(5)

式中：Fi为τ时间内砂粒i的荷载作用力之和；ni为砂粒i在τ时间内作用在风力发电机上的数量.

时间τ内经过风力发电机的流体体积为

V=2Al

(6)

式中:l为流体在时间τ内经过的距离；A为风力发电机迎风面的面积.假定风力发电机前后迎风面均遭受砂粒荷载的冲击作用，将结构的受力面积简化为2A，因此沙尘颗粒荷载为

(7)

假定风沙观测系统上游沙源在短时间内是不发生变化的，且在大于起沙风速的同一风速情况下，风沙流密度是相同的.基于以上假定，风沙观测仪中某一高度处集沙总质量为

Qh=∑ρiAvhT

(8)

式中:A为风沙观测系统的集沙器进沙口面积；vh为某一高度风沙流密度所对应的风速；T为在取样时间内起沙风所持续的时间.因此，由式(8)可知：

(9)

则沙尘颗粒荷载的计算公式为

FP(t)=4ρivP(t)2A

(10)

1.2 风荷载计算

在某一高度处，某一点的风速可表示为平均风速与脉动风速之和.平均风速具有随高度变化的规律，脉动风对结构的作用是随机的，应按随机振动理论求解其动力响应，本文采用对数律风剖模拟平均风速，采用随高度变化的Davenport风速谱模拟风力发电机的脉动风速.

1.3 风沙荷载计算

风沙流对挡风结构所产生的压强为净风产生的压强与风沙流中砂颗粒撞击挡风结构产生的压强之和.研究风沙颗粒荷载对风力发电结构的影响，实际上是研究风荷载与砂粒荷载叠加后的荷载对风力发电结构的影响，因此本文将沙尘颗粒荷载计算模型与风荷载模型相叠加，故风沙荷载计算模型为

F=FW+FS

(11)

式中:F为单位体积内风沙荷载；FW为单位体积内风荷载；FS为单位体积内沙尘颗粒荷载.

塔身风荷载计算公式为

(12)

式中:μS为风力发电机的体形系数；VS(t)为风速.

由式(11)和式(12)可得

(13)

2 数值模拟

2.1 风力发电机基本参数及模型

本文选用在西北地区风力发电厂中最常见的2 MW三桨叶风力发电机作为原始模型.其中，风轮半径为44 m，各叶片之间呈120°夹角，额定转速为17 r/min.轮毂质量为12 t，叶片质量14 t，机舱质量为60 t.在分析过程中，忽略机舱内部构件将其视为一个矩形，长度为8.44 m，宽度为3.56 m，高度为3.4 m.钢塔高80 m，塔底直径为4.1 m，塔底壁厚为25 mm.塔顶直径为2.5 m，塔顶壁厚为15 mm.塔壁厚度随高度线性变化.风力发电塔筒钢材的弹性模量为2.1×1011Pa，密度为7850 kg/m3.

采用大型有限元软件ANSYS建立三种不同精细程度的有限元模型，如图1所示.

图1 风力发电结构模型Fig.1 Wind turbine model

模型1：风力发电塔筒采用壳单元(shell181)模拟，风轮与机舱简化为集中质量点，采用mass21单元模拟，采用CERIG命令与塔顶形成刚性区域.塔底节点与地面固接.

模型2：在模型1的基础增加机舱、轮毂及叶片，将机舱及轮毂视为一个整体，采用beam189单元模拟，三个叶片采用有限应变梁单元(beam189)模拟.不同构件单元之间采用多点约束单元实现耦合连接.塔底与地面仍采用固接.

模型3：在模型2的基础上，考虑土-结构相互作用[8](简称SSI)对风力发电结构动力响应的影响.在近域地基面上设置粘弹性人工边界[9]，并将其等效为在土体界面上连续分布的并联弹簧-阻尼系统来模拟SSI作用.

2.2 风荷载模拟

为方便计算，对本文所采用的2 MW三桨叶风力发电机风荷载模拟点简化为17个点，具体模拟点的选取高度与数值见图2和表2.

图2 模型风模拟代表点Fig.2 Simulation representative point of model wind

表2 风荷载模拟点参数

图3 模拟点5的脉动风速时程曲线Fig.3 The wind speed time-history curve of simulation point 5

图4 模拟点5的脉动模拟功率谱与目标功率谱比较Fig.4 Comparison between the pulsating analog turbine spectrum and the target turbine spectrum of simulation point 5

2.3 风沙荷载模拟

由式(13)可知，砂粒荷载计算方法与风荷载一致，均与风速相关，且随高度变化而变化.由于风力发电结构为高耸结构，考虑到工程计算的安全性与误差要求，在实际应用中仅考虑砂粒的顺风向荷载.已知砂颗粒运动速度与净风速、气流动力粘度、砂粒质量及风洞边界尺寸特征值等因素有关，因风洞试验结果为风沙流的理想状态，实际野外风沙流中单位体积砂粒运动的平均速度根据相似准则与量纲和谐原理为[10]

(14)

式中：k为与砂颗粒物理性质(如级配、含水率等)有关的系数，通过采样分析获得；h为不受边界影响的充分发展的紊流区域特征尺寸，取值范围为0.2～5 m.

图5为仅考虑风荷载作用时风力发电机模拟点6的风荷载时程曲线，图6为考虑风沙荷载作用时风力发电机模拟点7的风沙荷载时程曲线.

图5 模拟点6的风荷载时程曲线Fig.5 The wind load time-history curve of simulation point 6

图6 模拟点7的风沙荷载时程曲线Fig.6 The wind-blown-load time-history curve of simulation point 7

3 动力响应分析

基于风速脉动时程与风沙荷载时程理论，利用ANSYS软件分别对三种不同模型下的风力发电结构在风荷载与风沙荷载的动力响应进行了分析，按照6种不同等级的风沙荷载分别模拟了不同风速下的风荷载和风沙荷载，通过对不同等级风荷载单独激励与风沙荷载共同激励下结构的动力响应分析.求解出不同模型在不同沙尘暴等级下风荷载与风沙荷载对风力发电结构的响应.

模型1：根据叶片实际气动外型，利用叶素-动量理论(BEM)计算叶片各截面的气动荷载，最后通过叠加气动荷载，得到风轮的轴向气动推力，如图7所示.将轴向推力施加于模型1的塔顶集中质量点来模拟叶片的气动荷载，塔筒风荷载采用模拟得到的风荷载时程施加.计算不同工况下的风力发电结构动力响应.得到不同荷载工况下的塔顶最大位移，如表3所示.图8为特强沙尘暴作用下结构的顶部位移时程曲线.

图7 风轮额定转速下轴向气动推力Fig.7 Axial aerodynamic thrust at rated speed of the wind wheel

表3 模型1不同荷载作用下塔顶最大位移

图8 模型1特强沙尘暴作用下结构顶部位移时程曲线Fig.8 Displacement time-history curve of the top of structure under the action of super strong sandstorm

由表3可知，考虑风沙荷载作用时，在6种不同等级的沙尘暴作用下，风力发电塔顶的顶点位移均比单一风荷载作用下的位移大，其中当风沙等级为特强沙尘暴时，其风力发电机的塔顶位移最大，为159.85 mm，远大于风荷载单独作用下的塔顶位移.而当风沙等级为弱沙尘暴时，风力发电机的塔顶位移增幅最大，为29.23%.因此对于多风沙灾害地区的风力发电机，考虑风沙荷载作用时，不仅要关注最大等级沙尘暴对风力发电结构的作用，也要注意一般等级沙尘暴对风力发电结构位移增幅.

模型2：采用Davenport相关系数经验公式来考虑叶片和塔架之间的相关性，采用谐波叠加法计算风机叶片的脉动风速时程，基于修正的叶素-动量理论(MBEM)[11]，引入普朗特叶尖损失因子和轮毂损失因子模拟叶片两面存在的气压差导致的气流在叶尖和叶根部沿叶片展向产生的二次流.考虑叶片三维旋转效应，加入动态失速延迟模型对叶片升力系数及阻力系数进行修正.考虑湍流-尾流状态，在轴向诱导因子较大时进行Glauert修正.计算风力发电机叶片在不同风速、不同叶片转速及偏航角情况下的动态荷载，进而得到考虑风轮和塔架相干效应及叶片旋转效应的叶片脉动风速时程[12].

根据BEM理论，叶片相对风速Vrel计算公式为

(15)

式中：Vox、Voy分别为谐波叠加法计算沿顺风向和横风向的来流脉动风速；vrot为叶片旋转导致的线速度；W为诱导速度；Vbx、Vby分别为叶片振动速度，各速度关系参见图9.

图9 叶素气流速度矢量图Fig.9 Vector diagram of leaf airflow velocity

诱导速度W可由式(16)表示：

(16)

式中:B为叶片数；L为叶片升力；φ为入流角；ρ为空气密度；r为叶素沿叶片轴向的长度；n为推力方向的单位向量；F为普朗特叶尖损失因子.

其中，叶片入流角φ为

(17)

式中：Ω为叶片旋转速度；a为轴向诱导因子；a′为切向诱导因子；v∞为叶片处风速.

(18)

基于上述的风场模拟方法，得到考虑叶片旋转的脉动风速时程，如图10所示.从而得到不同风速下的风荷载及风沙荷载，采用图2的加载方式，得到模型2不同荷载工况下的塔顶最大位移如表4所示.图11为特强沙尘暴作用下结构的顶部位移时程曲线.

图10 叶片模拟点14风速时程Fig.10 The wind speed time-history curve blade simulation point 14

图11 模型2特强沙尘暴作用下结构顶部位移

表4 模型2不同荷载作用下塔顶最大位移

由表4可知，考虑叶片旋转效应时，在6种不同等级的风沙荷载作用下，风力发电结构的顶点位移均比相应等级风荷载作用下的响应大，这与不考虑叶片旋转效应得到的风力发电结构的位移响应变化规律相同.但在相同等级风荷载及风沙荷载作用下，考虑叶片旋转效应的风力发电结构的位移响应均比不考虑叶片旋转效应的风力发电结构大.

模型3：采用模型2中计算得到的叶片及塔筒脉动风速时程以及风力发电机风荷载和风沙荷载时

程进行加载.采用粘弹性人工边界模拟风力发电结构动力响应时的土-结构相互作用.在近域地基面上设置人工边界，将其等效为在人工截断边界上连续分布的并联弹簧-阻尼系统，如图12所示.其中弹簧元件的弹性系数Kb及粘性阻尼器的阻尼系数Cb为

图12 人工边界等效弹簧-阻尼示意图Fig.12 Schematic diagram of artificial boundary equivalent spring-damping

式中：ρ、G分别为介质的质量密度和剪切模量；R为散射波源至人工边界的距离,取地基表面中心点到人工边界的平均距离；土体尺寸为60 m×60 m×40 m；c为介质中的波速；参数α参照表5确定.场地土参数如表6表示.

表6 场地土体参数

表5 三维粘弹性人工边界参数α的取值

Tab.5 The value of three-dimensional viscoelastic

artificial boundary parameterα

类型方向α三维人工边界法向4/3切向2/3

注：各土层的泊松比均为0.2.

基于上述土-结构相互作用的模拟方法及风力发电机所处场地的土参数.采用模型2中得到的风荷载及风沙荷载时程进行有限元动力响应分析，得到模型3在不同沙尘暴等级下风荷载与风沙荷载对风力发电结构塔顶最大位移的影响，如表7所示.

表7 模型3不同荷载作用下塔顶最大位移

由表7可知，考虑土-结构相互作用时，削弱了地基和基础对结构的约束作用，增加了塔体基础与地基连接的柔度，因此风力发电结构在不同等级荷载作用下的动力响应均增大，塔筒顶部的最大位移相较于不考虑土-结构相互作用时大约增加25.6%.

4 结论

本文是以2 MW三桨叶风力发电机为模型，分别建立了塔顶集中质量模型、考虑叶片旋转效应模型以及考虑土-结构相互作用三种不同精细程度的有限元模型.通过数值模拟得到不同风沙等级的风荷载和风沙荷载时程曲线，对上述三种风力发电机有限元模型进行动力响应分析，得到如下结论：

1) 通过对比单独风荷载与风沙荷载作用下三种不同精细程度的风力发电结构的动力响应发现，考虑风沙荷载作用下的风力发电结构比单一风荷载作用下的塔顶位移大.因此虽然风荷载是风力发电结构的主要设计荷载，但沙尘颗粒对风力发电结构的动力响应具有明显的放大效应，因此在西北多风沙地区风沙荷载对风力发电结构的危害应予以关注.

2) 考虑叶片旋转效应计算不同等级风沙荷载以及风荷载动力响应时，风力发电结构的塔顶最大位移增大，因此在进行结构动力响应分析时，采用简化的塔顶集中质量模型，得到的结果偏于不安全，建议在进行动力响应分析时采用考虑叶片旋转效应的整机有限元模型.

3) 考虑土-结构相互作用计算不同等级风荷载及风沙荷载作用下风力发电结构的动力响应时，由于地基对动力响应的放大作用，其风力发电结构塔顶的位移增大，其增幅较为显著.因此，应考虑土-结构相互作用对风力发电结构的影响.

致谢：本文得到兰州理工大学红柳优秀青年基金的资助，在此表示感谢.