董占琢 廖 晖

（东方汽轮机有限公司，四川 德阳，618000）

0 引言

高耸结构横风向振动的机理较为复杂，影响因素很多，在工程结构中较为常见且机理相对清楚的横向风振内容包括：涡激振动[1]、驰振[2]、颤振[3]。本文所研究的塔筒横截面为规则圆形，不存在攻角问题，风绕流塔筒时不会发生驰振和颤振，主要是由卡门涡街的漩涡发放引起的垂直于来流方向的涡激振动。风绕流塔筒产生的卡门涡街以及升阻力方向如图1所示。

图1 卡门涡街与升阻力方向

1 塔筒涡激振动CFD计算

1.1 雷诺相似准则

风绕塔筒的流动主要受粘性力、压力和惯性力的作用。从力学相似的观点看，若两个流场在对应点作用的同种力方向相同、大小成同一比例，则满足动力相似。在几何相似的前提下，两个流动只要在对应点满足代表粘性力与惯性力比值的雷诺数相等，则表示压力与惯性力之比的欧拉准则必然相等，因此风绕流塔筒满足雷诺相似条件。只需给出不同雷诺数下的力系数的大小即可表示不同直径、不同风速下的受力。

1.2 二维圆形CFD计算

圆柱绕流问题是典型的钝体大分离问题，本文采用SST湍流模型进行计算。

1.2.1 计算网格

采用图2所示的计算网格。计算域大小为10D×20D，计算网格总数34302。

图2 计算网格与边界条件

1.2.2 边界条件

入口设置为速度入口条件，按照不同的雷诺数进行计算；出口设置为背压出口，出口为大气压；壁面设置为无滑移壁面条件；计算域上下边界为对称边界，用以模拟无穷大空间。

1.2.3 计算结果

主要计算塔筒表面受力情况，图3为随时间变化圆柱表面的升力系数与阻力系数。

图3 Re=3.34×106时的圆形升、阻力系数变化

由图3可以看出，升力系数和阻力系数都随时间周期按正弦或余弦变化，阻力系数变化幅值很小，可视为不变。升力系数呈关于0线的余弦变化，升力系数以其幅值Cl给出，频率以斯特劳哈尔数的形式给出。

1.2.4 CFD二维计算模型与三维计算模型的比较

我们对圆柱的三维和二维建模都进行了分析，结果对比如表1所示：

表1 二维、三维计算结果的比较

三维和二维的计算结果中，阻力系数和斯特劳哈尔数基本相同，而升力系数中三维的结果明显比二维要小[3]。二维结果更加保守[4]，以二维圆形计算结果来作为塔筒载荷。

1.2.5 CFD计算结果同DIN1055-4标准的比较

GL规范[5]规定横向振动载荷可按DIN1055-4[6]来计算。图4给出了DIN1055-4标准查表数值和本文CFD计算数值的比较。

图4 升、阻力系数DIN标准和CFD计算结果比较图

本文的计算忽略了转捩等复杂因素的影响，因而趋势与实验趋势有所偏差。将在下一节具体分析偏差产生的原因及修正方法。

1.2.6 二维CFD计算结果的修正

旋涡脱落引起的力是复杂的流体力学现象的结果，对描述流体和结构物理特性的许多参数敏感。下面对转捩、表面粗糙度、来流湍流度、实验影响等因素对结果的影响进行分析。

采用关联转捩的模型之后，计算所得的圆柱的阻力系数与Schlichting[7]实验曲线的趋势能够吻合，如图5所示。但由于采用的是经验关联式，与实验值的绝对值有一定的差别。但从图6中可以看到，考虑转捩之后的结果比全湍流的结果偏小。因此，尽管全湍流计算模型有一定局限，为了分析的安全性，我们仍然采用全湍流的计算结果作为以后的分析基础。

图5 考虑转捩后的阻力系数计算结果与实验值比较

图6 转捩对圆柱阻力系数的影响

粗糙表面对计算结果产生了较大的影响，从图7可以看出，总体趋势是使阻力系数和升力系数提高，而使斯特劳哈尔数降低。

图7 表面粗糙度对圆柱升、阻力系数及斯特劳哈尔数的影响

湍流度对圆柱绕流有一定影响，如图8所示，湍流度增大会使阻力系数降低。

图8 湍流度对圆柱阻力系数的影响

目前，有关圆柱绕流所受阻力的实验值一般取Schlichting[7]的数值，然而，来流湍流度、圆柱表面粗糙度、压力测量方式等都会对实验结果产生影响，实验结果也有很大误差。

考虑表面粗糙度、来流湍流度等的变化和影响，升、阻力系数的修正以偏安全计算作为准则，以防止实际中由于各种参数等的变化产生危险工况。由于本文后续谐响应计算方法，对于斯特劳哈尔数的精度要求不高，因此对其不作修正。对于升、阻力系数，本文推荐修正系数为1.5。经过修正后的升阻力系数如表2所示 （篇幅原因，此表仅为部分计算结果）。

表2 修正后的圆柱升力系数、阻力系数与斯特劳哈尔数

2 塔筒模态分析

为了分析FD82E风力机塔筒在正常运行与吊装过程中的气流激振安全性，需对其正常运行与各吊装环节进行模态分析。分为四种工况，具体定义见表3。

表3 FD82E风力机塔筒模态分析计算工况

下面以工况1为例，介绍模态分析流程及方法。

2.1 工况1塔筒有限元模型及边界条件

塔筒按图纸建模。塔筒顶部用一集中质量点模拟机舱、轮毂和叶片的总质量，坐标 （m）为（-0.5957， 0， 67.867）， 总质量为 106566 kg， Z方向转动惯量Iz=3540000 kg·m2。塔筒底部同样用一集中质量点模拟地基， 坐标为 （0，0，-1.5896），总质量为870040 kg，X、Y方向转动惯量Ix=Iy=9165000 kg·m2。集中质量点分别与塔筒顶部和底部面建立刚性区域。地基质量点周围建立X、Y方向两个扭转弹簧模拟土壤约束，弹簧刚度系数取3×1010N·m/rad，弹簧末端节点坐标为（5， 0， -1.5895）， （0， 5， -1.5895）。整个塔筒网格均为六面体，壁厚方向划分三层网格。模型总单元数为38866个，总节点数为189358个。有限元模型中焊缝壁厚突跳处以及门框和筒体之间采用bonded接触。地基质量点加UX，UY，UZ，ROTZ四个方向约束，扭转弹簧末端采用全约束，如图9所示。

图9 工况1条件下塔筒模态分析有限元模型

2.2 工况1模态计算结果及安全性校核

工况1模态计算结果见表4，具体振型如图10所示。

表4 塔筒前6阶模态计算结果

图10 工况1塔筒X、Y向前两阶弯曲振动振型

已知发电机转速为1100～2000 r/min，风轮转速为 10.58～19.23 r/min，叶片 1P转动频率为0.176～0.321 Hz， 3P 转动频率为 0.528～0.962 Hz。满足GL规范[5]要求的5%避开率。

2.3 所有工况模态分析结果

工况2、3、4的模态结果如表5所示，无风轮激振，风致振的校核会在后续章节分析。

表5 所有工况塔筒前6阶模态计算结果

3 塔筒涡激振动谐响应分析

3.1 谐响应分析载荷简化

实际中由于塔筒外径及风速随高度不断变化，因此雷诺数沿高度方向不为定值，造成实际载荷状况十分复杂，难以分析。本文对风绕流塔筒载荷做如下简化假设：即假设塔筒沿高度方向作用的激振力频率、幅值和相位均相同，激振力幅值取共振频率附近CFD计算得到的最大值。此种简化方法使激振力对塔筒的作用放大。

实际塔筒在风绕流作用下，某一高度塔筒横截面压力分布不均，且随时间变化。CFD中的载荷结果是将横截面的压力换算成作用在中心轴上的集中力。考虑到载荷简化，沿高度方向塔筒载荷相当于作用于中心轴的均布载荷，进一步简化为塔筒顶部的集中载荷和弯矩。

图11 简化载荷雷诺数计算取值说明

3.2 工况1条件下塔筒谐响应分析

此工况为风力机吊装完毕的正常运行条件，包括所有风力机部件。在考虑风绕流塔筒的横向激振力之外还需考虑风轮转动对塔筒的激振作用。

由于CFD计算本身存在一些误差，安全起见取与工况1固有频率0.453Hz最为接近的较大结果 （参见表2）， 此时雷诺数约为 1.95×106， 风速约为7m/s，均布载荷为58.73N/m （依照GL规范，此载荷已乘以安全系数1.35）。假设风来流方向在X向，则风绕流塔筒产生的横向激振力在Y向，由理论力学简化得到塔筒顶部集中载荷幅值为Fy=3884.99N，Mx=128496.04N·m。从计算结果可以看出，风力机正常运行时塔筒共振频率附近的风绕流激振力很小 （仅有不到4kN），对塔筒影响也很小，可以忽略。

3.3 工况2条件下塔筒谐响应分析

3.3.1 工况2塔筒谐响应分析有限元模型

工况2塔筒涡激振动谐响应分析有限元模型除塔顶质量点的设置与工况1模态分析的有限元模型有差别外，其他部分完全相同。工况2谐响应分析塔顶质量点用于模拟机舱质量，坐标 (m)为 (-1.19，-0.02996，67.69)， 质量为69000kg，Z方向转动惯量Iz=359500 kg·m2，X方向转动惯量Ix=91270 kg·m2， Y 方向转动惯量 Iy=393600 kg·m2。塔筒顶部同时建立一无质量点，用于加载集中载荷，坐标为(0，0，66.15)。网格均为六面体，塔筒壁厚方向划分三层网格。模型总单元数为38867个，总节点数为189359个，如图12所示。

图12 工况2条件下塔筒谐响应分析有限元模型

3.3.2 工况2加载方式与边界条件

取与工况2固有频率0.542Hz最为接近的较大结果 （参见表2），此时雷诺数约为2.5×106，风速约为9m/s，均布载荷为95.18N/m(已考虑安全系数1.35)。经简化得到塔筒顶部集中载荷幅值为Fy=6295.83N，Mx=208234.45N·m，阻尼系数值取0.005，相位角均取零，激振频率范围为0.541～0.5425Hz，计算步数为15步。有限元模型中焊缝壁厚突跳处以及门框和筒体之间采用bonded接触。地基质量点加UX，UY，UZ，ROTZ四个方向约束，扭转弹簧末端采用全约束。

3.3.3 工况2谐响应计算结果

由于结构中定义了阻尼，所以结构响应与激振力之间不同步，存在相位差[5]。ANSYS计算的结果会以复数形式存储，而实际的计算结果是由实部和虚部的SRSS值给出。

由振动理论可知，谐响应激励作用下塔筒一阶共振时塔筒顶部位移最大，因此取塔筒顶部节点为观察节点。可观察到的塔筒顶部节点Y向振幅随激振频率变化，如图13所示。

图13 工况2条件下塔筒顶部节点Y向振幅随激振频率变化

由图13可知，在激振力频率为0.5419Hz时，塔筒顶部节点位移最大，最大位移为0.103304m。由相位角-90.7555°可判断此时已基本与共振峰值点重合。进一步观察0.5419Hz下塔筒共振时整体位移和应力分布情况，如图14和图15所示。

图14 工况2条件下塔筒共振时应力分布/Pa

图15 工况2条件下塔筒共振时位移分布/m

由图14和图15可知，共振时实部计算结果为：最大应力3.43MPa，最大位移0.001636m。虚部计算结果为：最大应力25.2MPa，最大位移0.107364m。

综合实部、虚部计算结果可知，塔筒产生的最大应力为25.43MPa，位于壁厚发生突跳的焊缝处。塔筒顶部最大振幅为0.107376m。由此可知，工况2条件下风绕流塔筒产生的横向激振力对塔筒的影响很小，几乎不对塔筒造成危害。

3.4 所有工况谐响应计算结果

由于工况3和工况4计算过程与工况2类似，不再详述过程。工况3条件下塔筒产生的最大应力为131.535MPa，塔筒顶部最大振幅为0.5027m。工况4塔筒产生的最大应力为265.55MPa，塔筒顶部最大振幅为0.5075m。塔筒仍符合安全性要求。

4 塔筒安装建议

塔筒安装主要关心高于多少风速时不能继续实施吊装，主要考虑风激振频率与已安装塔筒固有频率是否会发生共振。此时研究的工况主要为工况2、3和4，因工况2条件下塔筒与风激振力共振时影响很小，因此忽略此工况。图16为激振力频率随风速的分布。

图16 激振力频率随风速分布

风力机吊装能否实施取决于风激振与这两种工况塔筒固有频率是否会发生共振，计算结果见表6。

表6 吊装工况下塔筒振动避开率

从表6可看出，塔筒在风速为15m/s时不满足避开率大于5%的要求，按照GL规范要求取安全系数1.1，此时吊装风速应为14/1.1=12.73m/s。相关文献[8,9]风速大于12m/s时，风力机应停止吊装，以保证工程安全。结合以上计算结果，建议当风速大于12m/s时FD82E风力机停止吊装。

5 结论

通过对FD82E塔筒涡激振动、固有频率、谐响应等方面分析计算，得出了以下结论：

（1）塔筒附近流动符合雷诺相似准则，升力系数、阻力系数和斯特劳哈尔数等可近似认为仅与雷诺数相关。

（2）塔筒绕流的CFD计算中，采用二维方法比三维方法计算出的升力系数偏高，而阻力系数和斯特劳哈尔数基本相同，因而二维计算方法所得载荷更趋保守。

（3）采用转捩计算模型可以得到与阻力系数实验值趋势符合良好的结果，但转捩模型计算值比全湍流模型计算值偏低。采用全湍流模型计算可得到更加保守的载荷。

（4）表面粗糙度、来流湍流度等对塔筒的升力系数、阻力系数有较大影响，综合考虑这些因素，在实际计算中建议乘以修正系数，本文建议取值为1.5。

（5）在塔筒正常运行的工况1中，塔筒的固有频率和风轮振动频率满足GL规范要求的5%避开率。由工况1条件下的谐响应分析结果可知，风力机正常运行时塔筒共振频率附近的风绕流激振力很小 (仅有不到4kN)，对塔筒影响也很小，可以忽略。

（6）风绕流塔筒谐响应计算中，当激振力与塔筒产生共振时，塔筒位移与激振力的相位角相差约π/2。由工况2、3、4条件下谐响应分析结果可知，风绕流塔筒产生激振力对塔筒影响较小（最大应力265MPa，小于许用应力336MPa），不会产生应力破坏；

（7）通过分析计算不同风速下风绕流塔筒激振力与塔筒的共振情况，得出当风速大于12m/s时，风力机应该停止吊装。

[1]詹昊,李万平,方秦汉,等.不同雷诺数下圆柱绕流仿真计算 [J].武汉理工大学学报,2008,30(12)：129-132

[2]谢宇新.驰振稳定性分析及其工程应用 [D].天津大学：硕士学位论文,2003

[3]陈桂彬,杨超,邹丛青.气动弹性设计基础 [M].北京：北京航空航天大学出版社,2010

[4]王亚玲,刘应中,缪国平.圆柱绕流的三维数值模拟 [J].上海交通大学学报,2001,35(10)：1464-1469

[5]GL规范.GL Wind 2003

[6]DIN1055-4德国风载荷标准，NABau 2005

[7]何鸿涛.圆柱绕流及其控制的数值模拟研究 [D].北京交通大学：硕士学位论文,2009

[8]廖明夫,R.Gasch,J.Twele.风力发电技术 [M].西安：西北工业大学出版社,2009

[9]何显富,卢霞,杨跃进,刘万琨.风力机设计、制造与运行 [M].北京：化学工业出版社,2009