风机叶片三维绕流场数值模拟.doc

第十一届全国水动力学学术会议暨第二十四届全国水动力学研讨会并周培源教授诞辰110周年纪念大会文集风机叶片三维绕流场数值模拟周胡，王强，万德成*(海洋工程国家重点实验室，船舶海洋与建筑工程学院，上海交通大学，上海,200240, *通信作者： ) 摘要：本研究利用基于开源代码OpenFOAM所开发的两个求解器，对美国国家可再生能源实验室开发的Phase VI风机叶片的三维黏性绕流场进行数值模拟。这两个求解器分别是基于任意网络界面元法(Arbitrary Mesh Interface, AMI)的瞬态求解器pimpleDyMFoam和基于多参考系(Multi Reference Frames,MRF)稳态求解器MRFSimpleFoam。利用这两个求解器分别对相同桨距角、不同风速下三维风机叶片的复杂流场进行了数值模拟，计算得到叶片表面压力分布，叶片的推力、转矩、尾涡等气动力数据。这两种求解器的计算结果与实验结果进行对比分析，证明采用这两种求解器数值模拟三维风机叶片复杂粘性流场是有效和可靠的。关键词：风机叶片；三维黏性流场； AMI；MRF ；OpenFOAM1 引言风能是清洁、无污染的可再生能源之一，得益于机翼设计理论，材料技术、电力供应、叶片加工制造等技术发展，风力发电技术日益成熟，在可再生能源中成本相对较低，应用前景广阔。随着海上风能的迅速发展，风能的利用再一次吸引了全球学术界和工业界的目光。国外关于风机叶片三维数值模拟的起步较早，也取得了一些重要的成果。例如Hansen等1998年在附着在叶片上的旋转坐标系上建立计算模型，基于压力修正方法求解了不可压RANS（Reynolds Averaged Navier-Stokes）方程1，这是首次对完整的转子叶片的全NS数值模拟；Yuwei Li等利用动态重叠网格技术（Overset grid technology），使用CFDShip-lowa v4.5通用程序分别求解RANS方程的方法和大涡模拟的方法(DES, Detached Eddy Simulation)方法模拟了多种工况下的风机的气动力行为2；Zahle等使用基于结构化网格不可压有限体积法的EllipSys3D求解器模拟了风机转子和塔架相互作用，成功捕捉了叶片和塔架伴流的非定常相互作用，研究了塔影效应对风机气动力特性的影响3。相对而言，国内相关的研究工作较少。如李宇红等采用雷诺平均的三维可压NS方程，对风机叶片在特定风速下的绕流流场进行了三维数值模拟4；范忠瑶等采用CFD软件包FINE/Turbo对NREL Phase IV风机进行了研究，对气动力相关参数的数值解和实验室进行了详细对比，分别分析了粘性流场的流动细节5。目前，国内主要使用FLUENT等商用软件对风机进行数值模拟，还没有使用开源代码OpenFOAM对风机叶片进行数值模拟的。上海交通大学船舶与海洋工程系已经在OpenFOAM的开发和使用上做了不少工作。曹洪建、查晶晶、万德成6-7基于OpenFOAM求解器interDyMFoam，开发实现数值水池造波（包括推板和摇板造波）和阻尼消波；沈志荣和万德成8将六自由度运动模块植入OpenFOAM，开发了naoe-FOAM-SJTU求解器，实现了对船舶在波浪上运动的数值模拟。本研究利用基于OpenFOAM开发的两个求解器，即基于多参考系(Multi Reference Frames，MRF)稳态求解器MRFSimpleFoam和基于任意网络界面元法(Arbitrary Mesh Interface, AMI)的瞬态求解器pimpleDyMFoam，分别对固定桨距角、不同风速下风机叶片的三维粘性绕流场进行了数值模拟。探讨分析了叶片表面压力分布，叶片的推力、转矩、尾涡等气动力数据。采用两个求解器计算出的结果与实验结果进行对比分析，证明OpenFOAM这两种不同求解器数值模拟风机叶片三维黏性流场是有效和可靠的。本研究为模拟海上风机叶片复杂流场问题奠定了基础，为风力机叶片的设计、改进和研发提供了有效的数值分析手段。2 物理模型和计算工况2000年，美国国家可再生能源实验室（NREL）在美国国家航空航天局Ames 12080风洞对Phase VI型风力机进行了大量的实验9。Phase VI为一台额定功率为19.8kW的两叶片实验风机，转速为72RPM。叶片采用S809翼型（图1），风轮半径为5.029m。根据实验数据建立的风轮曲面如图2所示。本研究主要对Sequence S系列实验下四种工况的风机叶片的气动力特性进行计算。计算时不考虑塔架、轮毂等对叶片绕流流场的影响。所选用的实验工况下的初始参数如下表所示：表1 四种不同工况下的风速、空气密度、运动黏度和桨距角实验编号风速/(m/s)空气密度/(kg/m3)运动黏度/(m2/s)桨距角/()S050000051.2431.424E-53S1000000101.2461.423E-53S1500000151.2241.450E-53S2500001201.2201.454E-53 图1 S809翼型图及坐标定义 图2风轮曲面3 数值计算方法3.1 控制方程及离散方法本研究两种方法都基于雷诺平均的NS不可压黏性方程，在惯性参考系下表达为：质量方程： （1）动量方程： （2）其中：为密度，为脉动平均速度， 为压力，为动力黏性系数，为湍流黏性系数。其中使用多重参考系（MRF）时，坐标系随叶片一起旋转，这时惯性参考系下的NS方程就不再适用。必须建立旋转坐标系与惯性参考系的之间关系。上式的NS方程变为： （3）其中：代表旋转坐标系的位置矢量，代表非惯性参考系绕惯性参考系的旋转角速度。 代表了由于旋转引起的科氏加速度项，代表了离心力项。OpenFOAM使用有限体积法对控制方程进行离散，MRFSimpleFOAM稳态求解器采用SIMPLE算法处理速度和压力之间的耦合。差分形式上时间项采用定常格式，对流项采用一阶迎风格式，扩散项采用TVD格式。而PimpleDyMFOAM采用PIMPLE算法（PISO和SIMPLE耦合算法）处理速度和压力的耦合。差分形式上时间项用一阶欧拉隐式格式，对流项采用线性迎风差分格式，扩散项采用面限制的二阶格式。两种方法的湍流模式选取k- SST模型10。3.2 计算区域的网格划分和边界条件为了比较两种求解器数值模拟风机叶片三维绕流场的能力，采用了相同的网格。本文综合参考了其他作者的设置，也结合自身计算机能力，确定为入口至风机距离为5m，风机至出口距离为20m，计算域半径15m（图3）。结合使用ICEM-CFD和snappyHexMesh工具对背景网格进行了划分，最终的网格量约46万。边界条件主要包括入口处具有相同的风速；出口处假定流动充分发展；两种方法在叶片表面边界条件设置方面不一样，AMI方法需要设置成滑移边界条件（movingwallVelocity），而MRF则设置成壁面（fixedValue）。3.3 MRF方法和AMI方法说明任意网格界面元方法(Arbitary Mesh Interface, AMI)是由最新版本OpenFOAM-2.1.0新采用的处理网格运动的一种方法，它可以完全并行计算。这种方法的原理与在OpenFOAM的开发版本中的通用网格界面元(Generalized Grid Interface, GGI)11基本一致，但也存在一些细微区别。它通过带权重的差值实现网格相对运动时交界面内外数据的交互。多重参考系方法（Multi Reference Framework，MRF）是处理网格相对旋转或移动的一种常用方法。根据网格的相对运动，将求解区域划分为不同的子域，每个子域可以相对不同的参考系建立控制方程，在子域的交界面上则通过速度换算成绝对速度形式进行各子域流场信息的交换。MRF方法和AMI方法实处理旋转问题的区别主要表现在处理动态部分和静态部分的相互作用的方式上。 图3 过转轴和叶片展向的剖面网格 图4叶片表面的网格分布4 风机叶片气流流动分析本研究四种工况下的风速分别是5m/s，10m/s，15m/s和25m/s。其中15m/s和25m/s的风速对于风机叶片大部分展向截面来说已经达到失速条件，故而对这两个风速下的数值计算是十分具有挑战性的。4.1涡的分布四种风速下使用AMI数值模拟的风机尾涡的情况可以参见图5。该图由速度梯度张量第二不变量Q(the second invariant of the velocity gradient tensor)12等值面表达，并由速度值染色。从图中可以看到，在5m/s和10m/s的风速下，叶片附近流体大多还是附着流动。从5m/s到10m/s叶片随边附近的涡有加强。这两个风速下，我们可以观察到稳定的叶尖涡，特别是在5m/s的时候更为明显。此外由于叶片的根部截面从S809翼型迅速转化到圆形，变化剧烈，故这部分区域也观察到了十分明显的涡。此外我们也可以注意到尾涡大约在风轮后约一个叶片长度距离后迅速地消散了，这和之前本研究对风机后面约一个叶片长度区域内的网格加密有关。故需要观察到更长距离的尾涡变化，需要对风机后面的网格做更大范围的加密，这当然也会大大增加计算量。而在15m/s时整个叶片范围内都出现了非定常的流动分离，这在25m/s的时候表现地更加厉害。此外，在这两个高风速下，叶尖涡和叶片吸力面的流动分离的相互作用也非常明显，这种相互作用会导致叶尖涡破碎，从图中我们也可以明显看到在高风速下叶尖涡的长度有明显缩短。 图5 四个风速下的尾涡情况（AMI方法）4.2推力和转矩实验中的推力值是通过对压力在整个叶片表面积分得到的，故而只考虑了压力的成分而忽略了比重较小的粘性成分。故本文在求解推力的时候也将只考虑压力成分。实验的转矩通过安装在实验风机低速刹车盘上的传感器测得值加上一些修正得到。 图中，实验的推力和转矩由几个周期的值平均而得，图中的“I”代表了实验值在一个周期内的变化，即代表一个标准差。本研究的AMI方法推力和转矩值由第二圈和第三圈的数据平均而得，由于MRF方法是定常方法，MRF方法的推力和转矩则是稳定后的平均值。图6展示了AMI方法和MRF方法计算得到的推力和力矩的值与实验值比较的情况。从图中可以看出，在低风速下（5m/s和10m/s）AMI方法和MRF方法对于风力机推力的预测非常准确，基本都在实验值的一个标准差内。但是随着风速的增加，MRF方法对推力的预测结果明显偏高，作者分析由于MRF求解器是定常的，而高风速下风机叶片出现了明显的非定常性失速和分离，故而导致了预测的偏差。相比下使用非定常求解器的AMI方法的计算结果与计算值有更好的吻合。在力矩预测方面，10m/s下两种方法的预测结果都出现了较大偏差。图6 两种方法和实验值的推力和力矩的比较4.3 Cn和Ct由于推力和转矩都是通过将压力在整个叶片表面积分得到的，故只能反映风机的整体性能。进一步研究叶片不同位置处截面的空气动力系数，可以让我们更加深入地探讨风机的空气动力性能，也可进一步验证求解器在计算风机空气动力性能细节上的能力。Cn和Ct依旧是对压力积分得到的，但是和之前不同的是，这里仅对压力在特定的截面上进行积分。图7是5m/s和15m/s下两种方法得到的Cn、Ct与实验值的比较。可以看出，在5m/s的风速下，两种方法得到的Cn都与实验值非常吻合，这和该状态下推力和实验值吻合是一致的，相比于AMI方法，MRF计算得到的Ct值略有偏低，这也导致了5m/s风速下MRF方法计算得到的转矩偏小。15m/s风速下，两种方法计算的Cn除了叶跟处偏低外，其余结果几乎与实验值一样，但是Ct预测结果不是很理想，特别是r/R=0.95的位置处都没有预测到试验中Ct回升现象，这可能是叶片表面强的非定常性引起的。还有一点也需要引起注意，Ct实验值在r/R=0.47的位置处有比较明显的下降，这是由于叶片在处附近流体分离点迅速移动到导边处所致13，但是MRF方法没有精确预测，仍然处于高位，这弥补了r/R=0.95处没有上升的问题，从而15m/s下MRF方法计算的转矩与实验值非常接近，这并不能说明MRF方法的优越性。图7 5m/s和15m/s下两种方法得到的Cn(左)、Ct(右)与实验值比较4.4 压力系数Cp前面讨论的Cn和Ct是Cp在特定截面上积分得到的。为了进一步探讨风机空气动力性能的细节，以及探讨之前所述问题的最终原因，我们需要将叶片各个截面压力分布情况与实验进行对比。相较之前，这将更加具有挑战性，因为压力值是针对某一个特定的点而言的，这 也将更加明显地展示OpenFOAM求解器计算风机空气动力性能的能力。 比较了四种风速三个截面处（r/R=0.3,0.63,0.95） 数值模拟的结果与实验值的比较（图8）。从图中可以看出两个很重要的现象，一是两种方法压力面（曲线下部）的数值模拟结果都非常好，与实验值几乎重合；第二个现象是吸力面靠近导边处误差比较大，10m/s和15m/s下这种现象尤为明显，没有出现其他作者13-14计算结果中的峰值，这在一定程度上造成了前面计算得到的转矩偏低。作者推测可能是导边处的网格较疏所致。图8 四种风速下不同截面处的压力系数分布曲线5 结论本研究利用基于OpenFOAM开发的两个求解器，即瞬态求解器pimpleDyMFoam和稳态求解器MRFSimpleFoam，分别采用任意网格界面元法（AMI）和多参考系方法（MRF）对NREL Phase VI风机在5m/s、10m/s、15m/s和25m/s四个风速下的气动力性能进行了数值模拟。通过计算得到了风轮的尾涡情况、推力曲线及扭矩曲线，不同截面位置的Cn、Ct曲线以及Cp的分布曲线丰富的结果，通过与实验的对比，得出以下结论：（1）基于任意网格界面元法（AMI）的pimpleDyFoam求解器和基于多重参考系方法（MRF）的MRFSimpleFoam求解器模拟风机叶片三维复杂粘性流场是可行和有效的。两种方法得到的结果与实验值基本吻合。（2）相比较而言，由于AMI方法能够更准确模拟风机旋转的真实情形，特别在高风速时，模拟能力更好。在计算成本上，由于MRF方法定常求解器的快速收敛性，需要对风机气动力参数进行粗略预估时，MRF方法不失为一个很好的选择。（3）在失速情况下，叶尖涡与叶片吸力面的分离之间相互作用的严重的非线性性，叶片旋转的三维效应等给数值模拟造成了挑战，因此数值模拟的结果具有一定偏差，需要继续提高。（4）使用的snappyHexMesh网格加密方法在处理大曲率边界上具有一定局限性，这种局限性影响了气动力性能预测的精确性，为了得到更精确的结果，应着重关注叶片导边处网格质量。6 致谢本研究得到国家自然科学基金项目(11072154, 50739004)，海洋工程国家重点实验室自主研究课题基金 (GKZD010053-11)，上海东方学者人才计划基金(2008007)和英国劳氏船级社教育基金(LRET)的资助。参 考 文 献1 Sorenson N, Hansen M O L. Rotor performance predictions using a Navier-Stokes method. In 1998.2 Li Y, Paik K J, Xing T, etc. Dynamic overset CFD simulations of wind turbine aerodynamics. Renewable Energy, 2011.3 Zahle F, Srensen N N, Johansen J. Wind turbine rotortower interaction using an incompressible overset grid method. Wind Energy, 2009, 12(6): 594-619.4 李宇红, 张庆麟. 风力机叶片三维流动特性与气动性能的数值分析. 太阳能学报, 2008, 29(9): 1172-1176.5 范忠瑶, 康顺, 王建录. NREL Phase VI 叶片气动性能数值分析. 工程热物理学报, 2009, (010): 1665-1668.6 查晶晶, 万德成. 用OpenFOAM实现数值水池造波和消波. 海洋工程, 2011, 29(3): 1-12.7 曹洪建, 查晶晶, 万德成. 基于 OpenFOAM 编程数值模拟多圆柱绕流流动A. In 中国力学学会学术大会2009 论文摘要集, 2009.8 Shen Z R, Wan D C. RANS computations of added resistance and motions of ship inhead waves. In Proceedings of the 22nd International Offshore and Polar EngineeringConference(In press), 2012.9 Hand M M, Simms D A, Fingersh L J, etc. Unsteady Aerodynamics Experiment Phase VI : Wind Tunnel Test Configurations and Available Data Campaigns Unsteady Aerodynamics Experiment Phase VI : Wind Tunnel Test Configurations and Available Data Campaigns. National Renewable Energy Laboratory, 2001.10 Strelets M, Aeronautics A I o, Astronautics. Detached eddy simulation of massively separated flows. American Institute of Aeronautics & Astronautics, 2001.11 Beaudoin M, Jasak H. Development of a generalized grid interface for turbomachinery simulations with OpenFOAM. In 2008: 4-5.12 Hunt J, Wray A, Moin P. Eddies, streams, and convergence zones in turbulent flows. In 1988: 193-208.13 Srensen N N, Michelsen J, Schreck S. NavierStokes Predictions of the NREL Phase VI Rotor in the NASA Ames 80 ft*120 ft Wind Tunnel. Wind Energy, 2002, 169151-169.14 Li Y W, Paik K J, Xing T, etc. Dynamic overset CFD simulations of wind turbine aerodynamics. Renewable Energy, 2012, 37(1): 285-298.Numerical Simulations of the 3D Viscous Flows around Wind Turbine BladeZHOU Hu, WANG Qiang, WAN De-cheng *(State Key Laboratory of Ocean Engineering, School of Naval Architecture, Ocean and Civil Engineering, Shanghai Jiao To