风力机叶轮单向流固耦合工程分析方法研究

关 新，马 骁，陈 旭

（1.沈阳工程学院 新能源学院，辽宁 沈阳 110136；2.沈阳金山能源股份有限公司 新能源分公司，辽宁 沈阳 110000）

目前，国内外对风力机叶片根部的研究较少，对叶片根部的静力学分析并未涉及流固耦合分析和模态分析，所以应通过各种方法加强对风力机叶片根部的分析［1-2］。此外，风力机的叶片根部作为连接部位，所受应力最集中，也应作更深入的研究。

1 计算参数

1.1 风轮直径

风力机的风轮直径可用下式进行估算：

式中，P为风力机输出功率，本文P取5×106W；ρ为空气密度（ρ=1.225 kg/m3）；v为设计风速，v=12 m/s（6级风）；η为风力机效率；Cp=0.4；D为风轮直径，计算后得D≈79 m。

1.2 设计风速

根据中国气象信息中心提供的东北和西北风能资源数据，该地区平均风速为6 m/s。因我国风资源分布遵循威布尔分布函数，故设计风速选取12 m/s。

1.3 风力机翼型选择

翼型的气动性对风力机的运行和寿命起着关键的作用。目前，研发且相对成熟的翼型有瑞典的FFA-W 系列、美国的SERI 系列、NACA44xx 系列、NACA63-4xx系列、NREL-S809和WA等。本文选取国内自主研发的WA 大厚度钝尾缘翼型作为研究对象，该翼型主要是针对兆瓦级风力机进行设计的，厚度分别为45%、50%、55%及60%的翼型轮廓［3-4］，如图1所示。

图1 4种厚度翼型的几何轮廓

1.4 叶片的弦长和安装角

由于WA 大厚度钝尾缘翼型属我国自主研发的新型翼，暂时无叶片截面的详细参数，故参考NACA634 系列翼型对叶片各截面的弦长（C）和安装角（θ）等进行模拟计算，其参数如表1所示。

表1 叶片各个截面参数分布

2 风力机叶片建模

根据已确定的翼型弦长、安装角等参数，确定翼型各叶素离散坐标点数据(x,y,z)，对于两叶素之间的翼型进行二次样条延展。各个叶素空间离散坐标点保存为文本格式（.txt）作为SolidWorks 软件建模的基础数据。用曲线将同一截面的点连接起来，通过放样和扫掠建立风力机叶片的模型，风轮叶片投影轮廓如图2所示。

图2 叶片投影轮廓

最终确定叶片叶素截面及叶片截面放样如图3所示。

图3 叶片叶素截面及放样

将叶片与轮毂以面接触方式进行装配。图4为风轮流固耦合分析三维模型。

图4 风轮流固耦合分析三维模型

3 流固耦合控制

流固耦合分析是一门流体力学和固体力学相互交叉的学科，主要研究流体与固体之间的相互关系。流固耦合须同时定义流体域和固体域，将描述流体变化的变量和固体变化的变量归结为未知变量，其控制方程主要包括流体、固体和流固耦合的3 个控制方程。

3.1 控制方程

3.1.1 流体力学控制方程

流固耦合遵循流体力学的基本假设和相关定律，流体域遵循质量守恒、动量守恒和能量守恒：

式中，t为时间；ff为体积力；ρf为流体密度；v为流体速度矢量；τf为剪切力张量；E为单位质量内能；qf为单位体积热量损失。

3.1.2 固体控制方程

依据牛顿第二定律可以推导出固体的控制方程：

式中，ρs为固体密度；δs为柯西应力张量；fs为体积力矢量；为固体域当地加速度矢量。

3.1.3 流固耦合方程

根据守恒定律，在流固耦合的接触交界面位置上，流体和固体的应力、变形都应该是一一对应守恒的：

式中，τf为流体的应力；τs为固体的应力；nf为流体的单位方向向量；ns为固体的单位方向向量；df为流体的位移；ds为固体的位移［5-7］。

3.2 流固耦合分析流程

流固耦合分析的流体域与固体域共用一个模型，在流体场中求解后，将结果作为输入条件传递给固体场，流程如图5所示。

图5 流固耦合分析流程

将流体域和固体域分别命名为fluid 和solid，搭建流固耦合分析流程，如图6所示。

图6 流固耦合分析流程

4 耦合场建立工程求解

基于工程角度考虑，采用1:1 叶轮实体三维网格进行分析，其单元数和节点体量庞大，不利于工程计算求解，故将风轮实体三维模型缩小为实物的，其流固耦合模型如图7所示。

图7 叶轮流固耦合三维分析模型

4.1 流体域分析

进入Fluent 分析模块对固体域（solid）进行抑制，并分别设置Inlet（进口）、outlet（出口）、sym（对称面）、wall-ground（地面）、wall-fsi（耦合面），同时采用Navier-Stokes方程和SST湍流模型对其求解。因叶轮实体模型缩小为实物的，故边界条件也需进行相应调整（流体入口的风速为12 m/s×40=480 m/s），求解迭代至150步收敛，最终所得流场云图如图8所示。

通过流体分析的结果云图可知，叶片中心的流线呈一定的涡旋状。

4.2 固体域分析

因风力机叶轮为高分子复合材料，故需对材料进行单独设置并赋予叶片，其参数如表2所示。

表2 风力机叶片材料参数

图8 风轮流场分析

为了提高工程求解精度，对流体域进行压缩。网格控制为10 mm，划分后的网格质量为0.834 57，超过了0.7，满足工程要求，单元网格如图9所示。

图9 叶片网格划分

在Static Structural 下Imported Load 后，进行固体域分析。固体域计算结果如图10所示。

求解后可知，风力机叶片的最大应力分布在根部，应力值为111.24 MPa，小于叶片材料的屈服应力，叶片末梢的最小应力值近似于0 MPa。

图10 叶轮压力

5 结论

从计算结果进行分析，叶片根部动力响应能够影响整个叶片的性能和寿命。本文基于工程角度进行叶片实体建模并对叶轮进行流固耦合分析及动力响应研究，结论如下：

1）采用实际工程可使用的方法优化建模，并结合流固耦合分析的基础理论，采取单向流固耦合分析技术，对风力机风轮的流场和结构场进行耦合分析，具有一定的工程应用价值。

2）从风轮结构场分析结果来看，叶轮的叶根位置承载的疲劳载荷较大，故在叶轮实际生产制造过程中须加强叶片根部的结构强度；同时，叶片攻角对叶轮总体应力值影响较大（与叶片厚度等参数相比较），在风力机生产运行过程中以及风况变化情况下，适时调整桨距角对提高风轮运行稳定性有至关重要的作用。

3）由于缺乏试验数据，分析误差率在微观上难以估算；但从宏观结果上分析，符合风力机运行的工程实际情况，可利用3D 打印技术开展实验室风洞试验来修正计算模型，提高计算结果的可靠性。