小型垂直轴风力机叶片流固耦合分析

韩兆辉，郭志平，张艳锋

（1.内蒙古工业大学机械工程学院，内蒙古 呼和浩特 010051；2.内蒙古工业大学理学院，内蒙古 呼和浩特 010051）

1 引言

可再生能源在现今能源短缺的时代已经充分显示了本身所具有的潜力及重要性。风能作为可再生能源的重要组成部分，使用风力发电已经成为世界上最具前景的能源开发方式。叶片作为风力发电机的重要组成部件，直接影响风力发电机的工作效率。在风力机运行过程中会引起叶片发生变形和振动，随着叶片变形和振动的产生又会对流场中的气流产生影响，为典型流固耦合问题，所以对叶片与周围流场进行耦合作用分析十分重要。

国外研究人员对风力发电机组的研究早于国内。文献［1］对翼型的流固耦合理论做了系统总结；文献［2］对偏航条件下的大型风力机进行了仿真计算，对该风机的气动性能进行了一系列的计算分析；文献［3］基于流固耦合方法对风力机叶片进行了仿真分析，对单个叶片的应变分布进行了计算分析；文献［4］对大型风力机叶片进行了计算分析，考虑气动弹性下叶片的影响，并对大型风力机叶片性能进行了理论分析。

80年代中期开始，国内研究人员开始对风力发电机流固耦合等问题进行研究［5］。文献［6］针对离网小型垂直轴风力发电机，对风力发电机叶片进行了流固耦合仿真分析，研究了翼型的气动特性；文献［7］对风机建筑一体化中垂直轴风力机叶片和主轴的受力情况，使用ANSYS Workbench运用单项流固耦合的方法，在不同来流条件下对叶片和主轴的静应力进行了分析和比较；文献［8］针对在ANSYS Workbench软件中不能在Fluent模块下直接观察到叶片在流场中所受的压力，通过软件ANSYS-CFX 组合，找到分析流固耦合的方法，使风轮数值计算分析更合理、更准确；文献［9］引入类似结构方程组处理流体区域连续变形，对流体网格节点坐标进行更新，采用强耦合方法对垂直轴风力机叶片进行了流固耦合分析，研究了叶片的气动性能。

流固耦合法是观察分析固体场在流体场作用下的一系列行为及固体场反作用于流体场时流体场的变化情况。应用ANSYS Workbench仿真软件对垂直轴风力机叶片进行双向流固耦合分析，是对仿真叶片进行模拟分析的重要手段。所得结果可以为叶片优化设计、研发以及安全运转提供一定的参考。

2 流固耦合理论基础及方法

2.1 流体运动守恒方程

流体运动必须遵循基本的方程，包括连续性方程、动量方程-牛顿运动定律、能量方程-热力学第一定律。这三个守恒方程是CFD的理论基石。针对这里所涉及的问题，连续性方程和动量方程可通过如下方程进行描述。

连续性方程为：

式中：t—时间；ρ—流体密度；vx、vy、vz—沿X轴、Y轴、Z轴的速度。

任何流体问题都必须满足质量守恒定律，即单位时间内流体微元体中质量的增加，等于同一时间间隔内流入该微元的净质量［10］。动量方程-牛顿运动定律：

式中：F—合力；m—质量；v—速度；ρ—密度；

2.2 流固耦合方程

在风力机运转过程中叶片与流场之间相互作用。叶片受到来自流场压力的同时会发生变形，变形后的叶片反作用于流场，进一步扰动叶片所处的流场。因此，在流固耦合壁面处应满足以下方程。

式中：τ—应力；d—位移；q—热流量；T—温度。

3 三维模型的建立

这里利用SolidWorks建模软件对风力机叶片进行三维实体模型的建立。风力机叶片三维实体模型与室外实验场地使用的风力机叶片为同一规格型号。这里所使用仿真叶片翼型是NACA0021，该仿真叶片内部由三根铝管加强筋贯穿五片支撑片，叶片两端由挡板固定，外表进行蒙皮，使用铆钉固定组装而成。建立的仿真叶片模型总长h=1200mm，弦长c=265mm，重量G=4.28kg，最后进行叶片表面处理，这样使仿真叶片模型趋近于室外实验所用模型，且对叶片强度以及刚度都有所增强。具体模型，如图1所示。

图1 叶片整体模型Fig.1 Blade Integral Model

4 实验设备、测试方案及结果

4.1 实验设备

该实验使用的风洞隶属于内蒙古自治区风能太阳能利用机理及优化重点实验室［11］。实验对象为小型垂直轴风力机叶片，实验采用位移测试分析系统对叶片进行耦合实验，系统主要包括位移传感器、位移信号采集分析软件、加速度传感器和数据线等配套设备组装而成。

4.2 测试方案

利用位移测试分析系统在垂直轴风力机研究室进对实验叶片进行位移测试，对实验叶片进行监测点布置，在叶片端部尾缘、0.25h尾缘与0.5h尾缘处分别选取一个监测点，定义为监测点一、二、三。在叶片监测点处布置与实验叶片表面方向垂直的位移传感器，对位移传感器编号与其对应的监测点编号一致，选用实验风速为8m/s对叶片进行实验分析。实验所用支架为刚性结构，支架与叶片、地面采用螺栓连接。叶片置于风洞出风口处中间位置，弦长与风速垂直。三维实验测试简图，如图2所示。

图2 三维实验测试简图Fig.2 Three Dimensional Experimental Test Diagram

4.3 实验结果

通过在不同工况下对叶片进行位移实验，得到位移传感器X、Y、Z3个方向上的位移信号，对信号放大后进行分析处理，并对信号进行一系列函数分析验证。考虑实验叶片主要发生Y方向上的位移，所以提取出实验叶片在Y方向上的位移曲线。当风速为8m/s时，Y方向上各监测点实验与仿真曲线对比图，如图3所示。

图3 实验与仿真曲线对比图Fig.3 Comparison of Experimental and Simulation Curves

根据图3对比Y方向上不同监测点曲线图我们可知：当风速为8m/s时，实验曲线绕仿真曲线呈现上下波动，实验稳定值大于仿真稳定值。

出现以上问题为叶片实验环境激励成分复杂，叶片支架对监测点位移产生影响等，使实验极大值和极小值分别大于和小于仿真曲线极值。

当风速为8m/s时，Y方向上各监测点实验与仿真最大值，如表1所示。

表1 实验与仿真各监测点最大值Tab.1 Maximum Value of Experiment and Simulation at Each Monitoring Point

根据表1 可知监测点一、二、三实验与仿真差值分别为0.3098×10-6m、0.4992×10-6m、0.5091×10-5m，所得差值较小，且实验曲线与仿真曲线波动趋势一致，可得当前仿真模型结果有效。

5 静止叶片双向流固耦合分析

这里使用ANSYS Workbench中Fluent、Transient Structure及System Coupling 三个模块进行数值模拟计算，调用系统耦合器使数据在流体场及固体场之间相互迭代计算，完成对叶片的双向流固耦合计算，模块连接方式，如图4所示。该论文研究的叶片是参考某1000W 垂直轴风力机叶片，三维模型与实体叶片比例为1：1。

图4 双向流固耦合模块示意图Fig.4 Schematic Diagram of Two-Way Fluid-Structure Coupling Module

根据所绘制叶片的三维模型，对叶片外流场进行了针对性设计，对流场划分网格后，对耦合面进行优化。选中耦合面插入细节尺寸Sizing，对Element size进行控制以优化耦合面，优化后耦合面，如图5（b）所示。对固体场外表面网格划分精度与流体场内表面网格划分精度一致，保证两侧网格对应关系更好。由于本论文只针对垂直轴风力机叶片进行流固耦合分析，根据叶片本身实际尺寸情况，在制定计算域时将其设置为叶片的数倍以上。设置流体范围为（12×2.4×2.4）m，叶片置于距流场入口3m 处中间位置，叶片弦向与风向垂直。所画计算域网格数为1217788，计算域边界名称及其网格模型，如图5所示。

图5 计算域及其网格模型Fig.5 Computational Domain and its Grid Model

5.1 流体场条件分析设置

对三维叶片模型处理完成后，在Fluent中对流场进行设置。该论文对计算域湍流模型的选用为：realizable k-ε模型，设置湍流强度为5%。一般垂直轴风机有效风速介于（3～20）m/s之间，适宜风速为（6～8）m/s，该论文主要考虑中高风速下叶片各位置位移变形情况，所以选取入口风速值依次为8m/s、12m/s、15m/s。

计算前对流体场各个条件进行设置，为使本次计算模拟更加接近风机叶片实际工作情况，我们认为叶片已经安装在风机支架上，即流场四周均不受气流扰动影响，所以流场上下左右四个面均为对称面，前后两端分别为入口及出口，出口设置为自由流出，将流场内表面定义为双向流固耦合壁面。双向流固耦合中固体场会发生变化，进而影响流场，所以在计算仿真中采用动网格进行仿真计算。最后，合理的时间步长对于求解精度以及对计算机资源的利用可以达到预期要求，为求解准确性以及计算机资源的充分使用，设置时间步长为0.005s，每一步保存一次。当残差值达到1.0×10-3收敛。

5.2 瞬态结构计算分析设置

对流体场条件设置完成后，进入Transient Structure中对叶片进行设定。确定叶片材料为6061t6铝。对叶片进行网格划分后，将叶片表面做细化处理，精度与流场内表面一致，使结果更准确。定义约束面1为远端位移，使约束面1沿叶展方向自由移动，其余方向全部固定，约束面2为固定约束。在瞬时分析中的时间步长采取与流体场中一致，在此不单独设置，即为0.005s。

6 计算结果分析

在对垂直轴风力机叶片模型进行非稳态分析前进入System Coupling中对时间步长及截止时间进行合理设定，保证工作站资源的合理有效利用。考虑到叶片结构为上下对称且最大变形出现在靠近尾缘处，我们在叶片端部尾缘、0.25h尾缘与0.5h尾缘处分别选取一个监测点，定义为监测点一、二、三，以便观察叶片在受到流场压力过程中3个位置的网格位移变化情况。观察点位置，如图1所示。

图6（a）为当来流条件为8m/s时流体场叶片位置的风速流线图，图6（b）、图6（c）为叶片在来流方向及背流方向的压力分布。由图6（b）、图6（c）可知，叶片来流方向压力大于叶片背流方向压力。当来流条件设置为12m/s与15m/s时，压力变化趋势与8m/s时相同。不同风速下叶片位移变化曲线及3个监测点的总网格位移曲线，如图7所示。

图6 风场风速流线及叶片来流方向、背流方向压力图Fig.6 Wind Speed Streamline of Wind Field and Pressure Diagram of Flow Direction and Back Flow Direction of Blade

图7 不同风速下叶片位移变化曲线及3个监测点的总网格位移曲线Fig.7 The Blade Displacement Curve Under Different Wind Speeds and the Total Grid Displacement Curve of Three Monitoring Points

由图7 可知：当风速为8m/s 时，网格位移量最小为1.5129×10-23m，最大为7.0861×10-6m，当风速为12m/s 时，网格位移量最小为1.1328×10-22m，最大为1.1328×10-5m。当风速为15m/s 时，网格位移量最小为8.4104×10-19m，最大为2.3977×10-5。

比较3种不同风速条件下网格位移量可得，随着风速的增加位移量随之增加，由监测点一到监测点二到监测点三，最大位移量由大减小再变大，得到叶片产生最大位移的位置。即靠近叶片中部尾缘处。

叶片在受力过程中总网格位移曲线并没有呈现明显的上下波动，只出现一个波峰之后便趋于平缓，这是由于仿真叶片在约束选取上并不是使用叶片两端进行约束，而是分别选取了距离中部330mm处与叶片支架连接处进行约束。

监测点一、二、三均是在35步左右趋于平缓，可知叶片在受到流场传递的压力时，中间部位与两端的振动频率是相等的，网格位移量不同，在靠近叶片中部尾缘部位网格量达到最大值。

当设置来流条件为8m/s时监测点一、二、三在0步时的初始值与达到稳定后的值不同，监测点一的位移增量为2.11×10-6m，监测点二的位移增量为1.30×10-6m，监测点三的位移增量为2.21×10-6m。

0.5h尾缘处网格位移量大于叶片端部尾缘网格位移量大于0.25h尾缘处网格位移量，叶片端部尾缘与0.5h尾缘处网格位移量相差较小，说明叶片主要变形区域发生在靠近叶片中部尾缘处。当来流条件为12m/s、15m/s时，位移增量的变化趋势与来流条件8m/s时一致。

不同风速下各个监测点的曲线对比图，如图8所示。由图8可知，同一监测点在不同风速下初始值也不同，随着风速的增加，总位移初始值随之增加，稳定值也越大。

图8 不同风速下各个监测点曲线对比图Fig.8 Curve Comparison of Each Monitoring Point Under Different Wind Speeds

经过计算叶片端部尾缘的位移增量，可知不同来流情况下该点的位移增量值：风速为8m/s时，位移增量为2.11×10-6m；风速为12m/s 时，位移增量为3.36×10-6m；风速为15m/s 时，位移增量为4.88×10-6m。

当来流条件为15m/s 时，叶片端部尾缘处总网格位移增量最大。

由图可知，随着风速的增加，叶片端部网格位移增加量随之增加，监测点曲线波动幅值随之增加。0.25h尾缘与0.5h尾缘处监测点对比曲线有同样的趋势。

沿叶片展向各方向与总位移的曲线对比图，如图9所示。考虑到叶片上下对称，叶片取值范围为（0～0.5）h。

图9 沿叶片展向各方向与总位移的曲线对比图Fig.9 Comparison Diagram of Displacement Curves Along Each Direction of Blade Spanning

由图9可知叶片主要变形出现在（0.25～0.35）h处，沿各个方向上的位移最大值均出现在（0.35～0.4）h处。

沿着叶片展向x，y，z三个方向上的位移先增加后减少，然后再增加减少，另一半叶片曲线变化趋势与此相同。

y方向上的位移大于x方向与z方向的位移，并与总位移的趋势相同。

对比0.5h处的位移值，当风速为8m/s时，y方向的位移值为3.44×10-6m，而总位移的位移值为3.47×10-6m，两者相差不大，即叶片在受力过程中主要发生y方向上的变形。

对比叶片y方向上的位移量，沿叶片展向位移量呈现上下波动，伴随风速的增加位移量也增加，可知风速与叶片在耦合面上的摆振成正比。

7 结论

该论文合理应用ANSYS Workbench实现垂直轴风力机系统中叶片的双向流固耦合作用，使用实验和仿真相互验证的方法，对叶片进行合理分析。根据计算结果对叶片的位移变形进行分析，分析结果表明：（1）随着风速的增加，叶片位移增加明显，0.5h处大于叶片端点处大于0.25h处。（2）风载与叶片位移成正比，最大变形出现在（0.25～0.35）h、（0.65～0.75）h 处，最大位移出现在（0.35～0.4）h、（0.6～0.65）h处，y轴正方向与来流方向相同，风载是影响叶片产生变形的主要载荷，且沿叶片展向位移值上下波动。（3）通过开展风力机叶片流场和结构场双向流固耦合数值模拟计算可知，风载与叶片在耦合面上的摆振成正比，且沿叶片展向上下波动。