风力机复合材料叶片流固耦合分析*

张 旭， 刘安宇

(天津工业大学 机械工程学院， 天津 300387)

风力发电已经成为世界上最具发展前景的能源开发方式.目前，利用风能最为行之有效的装置是风力发电机，其中最具重要性的部件即为叶片，因为其性能与风力机运行现状及稳定性存在直接关联.在风力机的运行过程中，叶片周围空气的流动对叶片产生的作用力促使叶片出现一定程度的形变，这一变化会导致周边空气的流向发生变化.为了得到精准的叶片动力学响应，首要要考虑叶片及流场间存在的作用，即必须考虑叶片与其周围流场的耦合作用.

国外研究人员对流体和固体风力发电机组进行了大量研究.Tongchitpakdee等对偏航条件下的大型风机进行了数值模拟计算，并对其气动性能进行了系统分析[1]；Rachid等在定常情况及非定常情况下对NACA4415翼型的叶片弯曲及扭转效应进行分析研究，延长了风力机的使用寿命[2]；Kim等考虑了气动弹性对叶片的影响，并对大型风力机叶片性能进行了预测[3].

自80年代中期至今，国内研究人员对实体研究流固耦合问题与风力发电机整机的流固耦合进行了大量的研究.徐建源等利用ADAMS软件建立了风波联合作用海上风力机模型，并对其进行流固耦合分析，得到了复杂工况下的风力机叶片的动态响应参数[4]；李少华等利用Fluent软件并结合SSTk-ω模型对1.2 MW风力机旋转风轮流场进行模拟，但是并没有考虑风力机塔筒与流场间的相互作用对数据造成的影响[5]；胡丹梅等在水平轴风力机模型不同尖速比条件下对风轮下游流场速度进行测量，获得风力机叶片尾迹的流场定量信息，由于该测量方法容易受外界因素的影响，因此，流动情况与实际情况有一定误差[6-7]；田琳琳等结合制动盘理论与CFD方法，并通过Fluent软件对9台风力机的尾流互相干扰情况进行模拟，得到风力机尾流相互影响的结果[8]；陈海萍等研究了叶片周围流场对叶片结构特性的影响[9]；潘萍萍等通过Fluent软件对1.5 MW风力发电机塔筒进行流固耦合分析，得到了风力机塔筒的动态评估方法[10].

本文应用ANSYS软件对改进后的2 MW风力机叶片模型进行流固耦合分析，通过对传统翼型NACA4412的改进，应用铺层设计把叶梢、叶根、展向等厚度改成渐变式.首先得到风力机叶片及风轮的流场计算结果，再将风力机流场的计算结果施加到固体计算中，得到风力机的应力分布，进而得出应力最大值及变形量最大位置，所得结果可以为风力机的设计和安全运行提供指导[11].

1 流固耦合方法

流固耦合问题从数据传输的角度来看，其计算可以分为两种：单向流固耦合和双向流固耦合.单向流固耦合是指通过流体计算获得的压力、速度等数据可以通过连接界面导入固体结构，或者将从固体计算获得的网格位移转移到流体的计算.

计算流体力学简称CFD，主要使用离散点上的一系列变量值代替时域和空间域中的一系列物理量.通过建立可以表示变量之间关系的代数方程，将基本方程离散化，通过CFD的数值模拟，可以显示流动领域的运动.在通过CFD分析的工程应用中，可以深入研究耦合问题以实现实验指南，达到节约成本的目的.流固耦合分析流程图如图1所示.

图1 耦合系统流程示意图Fig.1 Schematic diagram of coupling system

2 三维模型建立

本文利用NASYS APDL建立风力机叶片的三维实体模型.采用传统翼型NACA4412应用复合材料建立叶片有限元模型如图2a所示，该模型总长44.5 m，翼型最大弦长3.1 m，根部为直径2.34 m的圆形，质量为10.6 t.为了加强叶片稳定性，在叶片中部添加两个厚为1.6 cm的加强筋，叶根部也进行了加厚，具体模型如图2b所示.叶片的建模工作主要在传统翼型NACA4412模型下进行改进，应用铺层设计把叶梢、叶根、展向等厚度改成渐变式，这样不仅更加接近实体模型，在强度、刚度上也有所加强，风轮模型如图2c所示.

3 静止单叶片流固耦合分析

3.1 流场计算分析

根据叶片结构在ANSYS软件中绘得实体模型，有针对地对流体实际范围进行了设定，同时进行了网格划分，不过因常规风力机体积过大，因此在制定计算域的过程中需要将其设置为叶片的数倍以上.设置的流体域范围为202 m×103 m×63 m，所划分的网格数为4 601 159，计算域及计算域网格模型如图3所示.

在进行流场分析时应先抑制叶片固体域，以实现对流体域的保护，设置流体控制温度为25 ℃.依照入口实际条件对其速度进行控制，并将这一设定速度称为恒速，流体的方向则需要垂直于入口界面，结合实际情况设置流体速度大小为13 m/s，而流体出口处设定为自由流出，风力机叶片边界条件设置为光滑无滑移壁面.风场流速图如图4a所示，叶片在迎风面和背风面的压力分布如图4b、c所示，在背风面叶片受到较小的压力，而迎风面叶片受到压力最大.

图2 叶片及风轮模型Fig.2 Models for blade and wind wheel

图3 计算域及计算域网格模型(静止状态)Fig.3 Models for computational domain and computationaldomain grid (static state)

图4 风场流速及迎风面、背风面压力分布Fig.4 Wind field velocity and pressure distributionof windward and leeward surfaces

3.2 叶片静力分析

叶片翼面材料为玻璃钢，该材料为线性材料，用MP命令定义其材料特性.其各层的线性材料特性为正交异性，具体参数为：展向模量42.6×109Pa，径向模量16.5×109Pa，剪切模量5.5×109Pa，泊松比0.22，密度1 950 kg/m3.

通过Static Structural接收风场的分析结果数据，并实现结构静力的相关研究.首先单独保留叶片，确保流程部分处于制约状态.由于叶片实际网格划分直接影响最后的研究结果，因此在划分过程中要确保其规则性.在Static Structural模块下选择Engineering Data，定义材料参数密度、弹性模量、泊松比和剪切模量，导入固体模型，定义网格划分参数划分网格，设置的网格数为980 916，定义固定约束.划分网格后的模型如图5所示.

图5 叶片网格模型Fig.5 Model for blade grid

4 旋转状态风轮流固耦合分析

4.1 流场计算分析

将Fluent模块的Solution与Static Structural模块的Setup建立连接进行单向流固耦合分析，导入风轮分别建立旋转流体域和外部空气流体域，旋转流体域半径为48 m，高度为6 m.外部空气流体域尺寸为606 m×296 m×296 m，设置网格划分参数及膨胀层，计算域及计算域网格模型如图6所示.

图6 计算域及计算域网格模型(旋转状态)Fig.6 Models for computational domain and computationaldomain grid (rotation state)

设置室温为25 ℃，入口速度为13 m/s，定义旋转轴为y轴，旋转速度为1.894 49 rad/s，风场流速图如图7a所示，通过Contour云图显示的压力云图如图7b所示，可以看出压力从叶根到叶尖逐渐减少，叶根最大，而在叶尖处最小.

图7 风场流速及风轮压力分布Fig.7 Wind field velocity and pressuredistribution of wind wheel

4.2 风轮结构静力分析

在Static Structural模块下选择Engineering Data进行材料参数定义：密度为1 950 kg/m3，弹性模量为42.6×109Pa，泊松比为0.22，剪切模量为5.5×109Pa，导入固体模型，划分网格，设置叶轮网格数为1 393 827，定义固定约束.划分网格后的模型如图8所示.

图8 风轮网格模型Fig.8 Model for wind wheel grid

4.3 结果输出

导入流体分析结果即压力，选择导入压力的面为叶片外表面，得到应变分布如图9所示.由图9可知，在额定速度下风力机运行时，风力机结构受风场的影响出现变形，变形的趋势与在风场中所受到的叶片表面压力分布是基本吻合的，此外在形变量最大的叶尖处，相对于叶根处的约束部分最远，并且此处叶片材料是最薄的，与文献[11]实验结果基本吻合.

图9 风轮应变分布Fig.9 Stain distribution of wind wheel

5 结 论

本文合理应用ANSYS Workbench实现风力系统中风轮及叶片的单向流固耦合分析，得到两者在运行数据上存在的关联.实验中得到的数据为进行叶片疲劳及共振频率测试提供了一定的参考，为之后的研究提供基础.在实际运行工作中，叶片的固有频率为耦合之后整个系统的固有频率，这个固有频率包括了风力机塔架在系统中对叶片的影响，系统中旋转离心力的作用，风场在系统中对叶片的影响，系统中重力作用等一系列的作用，这些影响与单个叶片在静止状态下的模态频率有一定的区别，但一般而言频率的变化很小.在今后的叶片优化设计过程中，同样需要关注该部分问题，以便得到的模型与实际运作更为接近，减小并优化设计周期.