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

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

《风机叶片三维绕流场数值模拟.docx》由会员分享，可在线阅读，更多相关《风机叶片三维绕流场数值模拟.docx（13页珍藏版）》请在冰豆网上搜索。

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

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

周胡，王强，万德成*

（海洋工程国家重点实验室，船舶海洋与建筑工程学院，上海交通大学，上海,200240,

*通信作者：

dcwan@）

摘要：

本研究利用基于开源代码OpenFOAM所开发的两个求解器，对美国国家可再生能源实验室开发的PhaseVI风机叶片的三维黏性绕流场进行数值模拟。

这两个求解器分别是基于任意网络界面元法（ArbitraryMeshInterface,AMI）的瞬态求解器pimpleDyMFoam和基于多参考系（MultiReferenceFrames,MRF）稳态求解器MRFSimpleFoam。

利用这两个求解器分别对相同桨距角、不同风速下三维风机叶片的复杂流场进行了数值模拟，计算得到叶片表面压力分布，叶片的推力、转矩、尾涡等气动力数据。

这两种求解器的计算结果与实验结果进行对比分析，证明采用这两种求解器数值模拟三维风机叶片复杂粘性流场是有效和可靠的。

关键词：

风机叶片；三维黏性流场；AMI；MRF；OpenFOAM

1引言

风能是清洁、无污染的可再生能源之一，得益于机翼设计理论，材料技术、电力供应、叶片加工制造等技术发展，风力发电技术日益成熟，在可再生能源中成本相对较低，应用前景广阔。

随着海上风能的迅速发展，风能的利用再一次吸引了全球学术界和工业界的目光。

国外关于风机叶片三维数值模拟的起步较早，也取得了一些重要的成果。

例如Hansen等1998年在附着在叶片上的旋转坐标系上建立计算模型，基于压力修正方法求解了不可压RANS（ReynoldsAveragedNavier-Stokes）方程[1]，这是首次对完整的转子叶片的全NS数值模拟；YuweiLi等利用动态重叠网格技术（Oversetgridtechnology），使用CFDShip-lowav4.5通用程序分别求解RANS方程的方法和大涡模拟的方法（DES,DetachedEddySimulation）方法模拟了多种工况下的风机的气动力行为[2]；Zahle等使用基于结构化网格不可压有限体积法的EllipSys3D求解器模拟了风机转子和塔架相互作用，成功捕捉了叶片和塔架伴流的非定常相互作用，研究了塔影效应对风机气动力特性的影响[3]。

相对而言，国内相关的研究工作较少。

如李宇红等采用雷诺平均的三维可压NS方程，对风机叶片在特定风速下的绕流流场进行了三维数值模拟[4]；范忠瑶等采用CFD软件包FINE/Turbo对NRELPhaseIV风机进行了研究，对气动力相关参数的数值解和实验室进行了详细对比，分别分析了粘性流场的流动细节[5]。

目前，国内主要使用FLUENT等商用软件对风机进行数值模拟，还没有使用开源代码OpenFOAM对风机叶片进行数值模拟的。

上海交通大学船舶与海洋工程系已经在OpenFOAM的开发和使用上做了不少工作。

曹洪建、查晶晶、万德成[6-7]基于OpenFOAM求解器interDyMFoam，开发实现数值水池造波（包括推板和摇板造波）和阻尼消波；沈志荣和万德成[8]将六自由度运动模块植入OpenFOAM，开发了naoe-FOAM-SJTU求解器，实现了对船舶在波浪上运动的数值模拟。

本研究利用基于OpenFOAM开发的两个求解器，即基于多参考系（MultiReferenceFrames，MRF）稳态求解器MRFSimpleFoam和基于任意网络界面元法（ArbitraryMeshInterface,AMI）的瞬态求解器pimpleDyMFoam，分别对固定桨距角、不同风速下风机叶片的三维粘性绕流场进行了数值模拟。

探讨分析了叶片表面压力分布，叶片的推力、转矩、尾涡等气动力数据。

采用两个求解器计算出的结果与实验结果进行对比分析，证明OpenFOAM这两种不同求解器数值模拟风机叶片三维黏性流场是有效和可靠的。

本研究为模拟海上风机叶片复杂流场问题奠定了基础，为风力机叶片的设计、改进和研发提供了有效的数值分析手段。

2物理模型和计算工况

2000年，美国国家可再生能源实验室（NREL）在美国国家航空航天局Ames120×80风洞对PhaseVI型风力机进行了大量的实验[9]。

PhaseVI为一台额定功率为19.8kW的两叶片实验风机，转速为72RPM。

叶片采用S809翼型（图1），风轮半径为5.029m。

根据实验数据建立的风轮曲面如图2所示。

本研究主要对SequenceS系列实验下四种工况的风机叶片的气动力特性进行计算。

计算时不考虑塔架、轮毂等对叶片绕流流场的影响。

所选用的实验工况下的初始参数如下表所示：

表1四种不同工况下的风速、空气密度、运动黏度和桨距角

实验编号

风速/（m/s）

空气密度/（kg/m3）

运动黏度/（m2/s）

桨距角/（°）

S0500000

5

1.243

1.424E-5

3

S1000000

10

1.246

1.423E-5

3

S1500000

15

1.224

1.450E-5

3

S2500001

20

1.220

1.454E-5

3

图1S809翼型图及坐标定义图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.3MRF方法和AMI方法说明

任意网格界面元方法（ArbitaryMeshInterface,AMI）是由最新版本OpenFOAM-2.1.0新采用的处理网格运动的一种方法，它可以完全并行计算。

这种方法的原理与在OpenFOAM的开发版本中的通用网格界面元（GeneralizedGridInterface,GGI）[11]基本一致，但也存在一些细微区别。

它通过带权重的差值实现网格相对运动时交界面内外数据的交互。

多重参考系方法（MultiReferenceFramework，MRF）是处理网格相对旋转或移动的一种常用方法。

根据网格的相对运动，将求解区域划分为不同的子域，每个子域可以相对不同的参考系建立控制方程，在子域的交界面上则通过速度换算成绝对速度形式进行各子域流场信息的交换。

MRF方法和AMI方法实处理旋转问题的区别主要表现在处理动态部分和静态部分的相互作用的方式上。

图3过转轴和叶片展向的剖面网格图4叶片表面的网格分布

4风机叶片气流流动分析

本研究四种工况下的风速分别是5m/s，10m/s，15m/s和25m/s。

其中15m/s和25m/s的风速对于风机叶片大部分展向截面来说已经达到失速条件，故而对这两个风速下的数值计算是十分具有挑战性的。

4.1涡的分布

四种风速下使用AMI数值模拟的风机尾涡的情况可以参见图5。

该图由速度梯度张量第二不变量Q（thesecondinvariantofthevelocitygradienttensor）[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.3Cn和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方法的优越性。

图75m/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）对NRELPhaseVI风机在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）的资助。

参考文献

1SorensonN,HansenMOL.RotorperformancepredictionsusingaNavier-Stokesmethod.In1998.

2LiY,PaikKJ,XingT,etc.DynamicoversetCFDsimulationsofwindturbineaerodynamics.RenewableEnergy,2011.

3ZahleF,SørensenNN,JohansenJ.Windturbinerotor‐towerinteractionusinganincompressibleoversetgridmethod.WindEnergy,2009,12（6）:

594-619.

4李宇红,张庆麟.风力机叶片三维流动特性与气动性能的数值分析.太阳能学报,2008,29（9）:

1172-1176.

5范忠瑶,康顺,王建录.NRELPhaseVI叶片气动性能数值分析.工程热物理学报,2009,（010）:

1665-1668.

6查晶晶,万德成.用OpenFOAM实现数值水池造波和消波.海洋工程,2011,29（3）:

1-12.

7曹洪建,查晶晶,万德成.基于OpenFOAM编程数值模拟多圆柱绕流流动[A].In中国力学学会学术大会'2009论文摘要集,2009.

8ShenZR,WanDC.RANScomputationsofaddedresistanceandmotionsofshipinheadwaves.InProceedingsofthe22ndInternationalOffshoreandPolarEngineeringConference（Inpress）,2012.

9HandMM,SimmsDA,FingershLJ,etc.UnsteadyAerodynamicsExperimentPhaseVI:

WindTunnelTestConfigurationsandAvailableDataCampaignsUnsteadyAerodynamicsExperimentPhaseVI:

WindTunnelTestConfigurationsandAvailableDataCampaigns.NationalRenewableEnergyLaboratory,2001.

10StreletsM,AeronauticsAIo,Astronautics.Detachededdysimulationofmassivelyseparatedflows.AmericanInstituteofAeronautics&Astronautics,2001.

11BeaudoinM,JasakH.DevelopmentofageneralizedgridinterfaceforturbomachinerysimulationswithOpenFOAM.In2008:

4-5.

12HuntJ,WrayA,MoinP.Eddies,streams,andconvergencezonesinturbulentflows.In1988:

193-208.

13SørensenNN,MichelsenJ,SchreckS.Navier–StokesPredictionsoftheNRELPhaseVIRotorintheNASAAmes80ft*120ftWindTunnel.WindEnergy,2002,169151-169.

14LiYW,PaikKJ,XingT,etc.DynamicoversetCFDsimulationsofwindturbineaerodynamics.RenewableEnergy,2012,37

（1）:

285-298.

NumericalSimulationsofthe3DViscousFlows

aroundWindTurbineBlade

ZHOUHu,WANGQiang,WANDe-cheng*

（StateKeyLaboratoryofOceanEngineering,SchoolofNavalArchitecture,OceanandCivilEngineering,

ShanghaiJiaoTongUniversity,Shanghai200240,China

*Correspondingauthor:

dcwan@）

Abstract：

ThispaperpresentstwosolversdevelopedfromtheopensourcepackageOpenFOAMtonumericallysimulatethethreedimensionalviscousflowsaroundtheNRELPhaseVIwindturbine.TheunsteadysolverpimpleDyMFoambasedonarbitrarymeshinterfaceandthesteadysolverMRFSimpleFoambasedonmultireferenceframesareemployedtosimulatethecomplexviscousflowsaroundthebladeinthesamepitchangleandindifferentwindvelocity.Theaerodynamiccharacteristicssuchasbladesurfacepressurecoefficient,thrust,shafttorqueandsoonareanalyzed.Throughthecomparisonoftheresultsobtainedbythetwosolverswiththeexperimentresults,theabilityofthetwosolverssimulatingthe3Dviscousflowaroundwindturbinebladeisproved.

Keywords：

WindturbineBlade;3Dviscousflows;AMI;MRF;OpenFOAM.