翼型动态失速的非定常模拟方法

于佳鑫　陈江涛　王晓东　吴晓军　康顺

摘要：  为探究翼型动态失速的高可信非定常模拟方法，以FFA-W3-241翼型为研究对象，采用开源计算流体动力学求解器OpenFOAM开展翼型动态失速下的流动模拟。研究重叠网格和滑移网格2种不同网格运动形式、2种不同时间步长、2种不同计算周期和OpenFOAM默认湍流模型与修正的k-ω SST湍流模型对动态失速过程中翼型气动力的模拟精度，并对流场结构进行分析。结果表明：修正模型预测的翼型气动力和流场特征与实验值更接近;重叠网格在翼型的动态失速模拟中更具优势。

关键词：  OpenFOAM; 动态失速; 湍流模型; 风力机; 翼型; 重叠网格

中图分类号：  V211.41; TK83文献标志码：  B

Unsteady simulation method for airfoil dynamic stall

YU Jiaxin CHEN JiangtaoWANG Xiaodong WU Xiaojun KANG Shun

（1. Key Laboratory of Power Station Energy Transfer Conversion and System（Ministry of Education）， North China

Electric Power University， Beijing 102206， China;

2. China Aerodynamics Research and Development Center， Mianyang 621000， Sichuan， China）

Abstract： To explore the highly reliable unsteady simulation method of airfoil dynamic stall， the flow of airfoil under dynamic stall is simulated using the open source computational fluid dynamics（CFD） solver OpenFOAM taking the FFA-W3-241 airfoil as the research object. The accuracy of the aerodynamic simulation of airfoil during dynamic stall is studied under different conditions， that includes two different mesh motion forms（overlapping mesh and sliding mesh）， two different time steps， two different calculation cycles， and OpenFOAM defaulted turbulence model and modified k-ω SST turbulence model. The flow field structure is analyzed. The results shows that the aerodynamic and flow field characteristics predicted by the modified turbulence model are closer to the experimental value. The overset mesh is more advantageous in the dynamic stall simulation of airfoil.

Key words： OpenFOAM; dynamic stall; turbulence model; wind turbine; airfoil; overset mesh

-基金項目：  国家数值风洞工程项目（NNW2018-ZT7B14）;国家自然科学基金（51876063）

作者简介： 于佳鑫（1993—），女，辽宁建昌人，博士研究生，研究方向为CFD可信度分析和不确定性方法，（E-mail）820113965@qq.com

通信作者： 王晓东（1979—），男，北京人，教授，博导，研究方向为海上风电机组设计，（E-mail）wangxd@ncepu.edu.cn0引言

翼型失速分为静态失速和动态失速。静态失速是指迎角超过特定迎角值时，翼型气动性能下降的现象。动态失速是指在迎角快速变化时，迎角超过静态失速角导致翼型产生周期性振荡的现象[1]。在定常流动中，对于给定的翼型几何形状，失速角基本为固定值。然而，当翼型在包含静态失速角的迎角范围内运动时，最大升力的角度大大增加，且强烈依赖于振荡的速率和振幅。涡沿翼型由前缘运动至尾缘直至离开翼型，此时升力突然下降，升力和俯仰力矩曲线都产生较大的迟滞回环。动态失速研究对翼型的气动设计具有重要意义。

动态失速的研究方法主要有风洞实验[2]和数值模拟2种。风洞试验要求复杂、成本高，而且难以捕捉动态过程中的流体细节。受实验条件限制，研究人员大多采用数值模拟获取更详细的动态失速信息。动态失速是强非定常流动过程，准确的计算流体动力学（computational fluid dynamics， CFD）模拟难度较大。张彦军等[3]计算不同雷诺数对动态失速特性的影响。朱呈勇[4]采用非定常雷诺时均（unsteady reynolds average navier stocks， URANS）方法研究动态入流和翼型振荡对失速特性的影响。

CFD求解器OpenFOAM是一款开源软件，可以自行修改求解器以应对不同物理问题的需求，通过C语言代码编写脚本，可使复杂操作过程更加自动便捷。[5]目前用OpenFOAM求解翼型动态失速的研究较少，模拟方法有待验证。k-ω SST湍流模型是翼型数值模拟中广泛采用的湍流模型。[6]OpenFOAM中默认的k-ω SST湍流模型对黏性系数的求解有一定简化，针对风力机翼型动态失速问题，如何适当修改湍流模型的系数值得研究。

翼型的俯仰振荡是翼型绕特定轴进行的旋转振荡，因此数值计算时需要网格按照相同的方式运动。OpenFOAM实现网格运动的方式主要有3种，分别为动网格（Dynamic Mesh）、滑移网格（Sliding Mesh）和重叠网格（Overset Mesh）。动网格可根据设定的运动，使运动物体附近的网格发生实时变形与重新生成。翼型尾缘几何尺寸较小，为得到较准确的计算结果，网格尺寸往往也较小。采用动网格时，尾缘附近的网格距离旋转轴较远，变形较大，在翼型运动过程中容易出现负网格，导致计算发散。滑移网格需要建立静止区域和运动区域，在二者间设置交界面，运动区域内的网格整体运动，不会涉及到网格的变形与重新生成，因此不会产生负网格。滑移网格是动网格的简化形式，通常用于往复运动和旋转运动。重叠网格将复杂的流场区域进行分解，每个区域内独立生成高质量网格，区域网格之间有重叠和共享部分，可用于任意运动方式。重叠网格综合动网格和滑移网格的优点，在保证物体运动准确的同时又能保证运动过程中的网格质量。因此，网格运动形式采用滑移网格和重叠网格2种。

为探究翼型动态失速的高可信非定常模拟方法，本文基于OpenFOAM，采用URANS方法，研究网格运动形式和湍流模型对翼型气动力和流场特征的影响。

1数值模型与研究方法

1.1风洞实验回顾

文献[7]开展FFA-W3-241翼型的相关实验研究，采用的雷诺数为1.6×10⁶，湍流强度为1%，翼型弦长为0.60 m，整个翼型围绕x/c=0.4以正弦形式运动，运动形式示意见图1，其中α_mea为平均攻角，α_amp为攻角振荡幅值。

式中：f为振荡频率;c为弦长;U_∞为自由来流速度;ω为角速度。

根據实验结果，下文数值模拟选择2个平均攻角α_mea，分别为1.5°和15.6°，2个平均攻角对应的攻角振荡幅值α_amp分别为1.5°和1.8°。

1.2数值方法

1.2.1计算域和网格

计算域和边界条件见图2。

滑移网格与重叠网格的静止域和旋转域的几何尺寸保持一致。静止域为正方形，其边长为40倍弦长（40c），并将翼型旋转轴置于计算域中心。旋转域是以旋转轴为圆心、半径为3c的圆。

静止域的左侧边界为速度入口，湍流强度与实验设置相同，即为1%。翼型表面为光滑无滑移壁面，运动形式在dynamicMeshDict中设置，滑移网格和重叠网格的dynamicFvMesh分别设置为dynamicMotionSolverFvMesh类型和dynamicOverSet FvMesh类型，运动函数为旋转振荡函数oscillatingRotatingMotion。为实现平行流动，上、下边界速度和压力的边界条件设置为零梯度。

网格划分时保证y⁺<1，滑移网格与重叠网格的静止域外侧、旋转域和翼型周向节点保持一致。静止域外侧各边节点数为113。滑移网格交界面处节点数与翼型周向节点总数保持一致，并且均匀分布。2种运动形式的网格见图3。翼型上、下表面各分布200个节点，法向设置115个节点，壁面网格的法向膨胀率为1.1。FFA-W3-241翼型尾缘为钝尾缘，设置50个节点。

此外，绘制3种拓扑结构相同的网格进行网格独立性验证。各边的节点数与上述网格成1.5倍关系，网格参数见表1。

1.2.2湍流模型

采用k-ω SST湍流模型计算URANS方程。OpenFOAM的标准k-ω SST模型是一种全湍流涡黏模型[8]，湍流黏性μ_t的计算公式为

式中：a₁和b₁为封闭系数，a₁=0.31，b₁=1.0;k为湍动能;ω为湍流耗散率;S为应变率幅值;F₂₃参数参考OpenFOAM文档[5]取值。

1994年，MENTER[9]对该模型进行改进，改进模型的湍流黏性项为

本文在OpenFOAM中加入改进模型，并与默认模型进行动态失速仿真的数值对比。

2结果分析

2.1网格无关性验证

在3种网格精度下，滑移网格和重叠网格计算得到的升力系数C_L见图4。在相同网格精度下，滑移网格与重叠网格的结果一致。当攻角α<12°时，3种网格精度对升力系数C_L的预测一致;当攻角α>12°时，粗网格计算结果与其他2种网格预测结果有些差别，中等网格与密网格的预测表现一致。考虑到计算成本，采用中等网格进行动态失速计算。

2.2动态时间步长和计算周期验证

为得到稳定的计算结果并降低计算成本，进行时间步长验证。以平均攻角α_mea=15.6°、攻角振荡幅值α_amp=1.8°为例，设计4种时间步长，一个振荡周期内的步数分别为1×360、2×360、3×360和4×360步。第2～5个运动周期内各时间步长的升力系数CL时间历程见图5。

由图5可以看出，采用1×360和2×360步的结果与采用3×360步与4×360步的结果相差较大，3×360步与4×360步的结果几乎重合，可达到收敛。

选取时间步数为3×360，进行计算周期数收敛性验证，计算稳定后相邻2个周期气动力升力系数

C_L和阻力系数C_D计算结果的标准差见表2。由此可以看出，3个周期后标准差不再变化，可达到收敛。

综上所述，考虑到计算时间成本，采用3×360步数对应的时间步长，并选取运动3个周期后的结果进行分析。

2.3稳态特性分析

利用滑移网格和重叠网格分别采用OpenFOAM默认的k-ω SST模型和修正的k-ω SST模型计算升力系数C_L，结果见图6。攻角小于10°时，空气流动处于附着流动状态，未出现失速现象，2种形式的网格和2种湍流模型预测的结果与实验一致。失速以后，即攻角大于10°后，数值计算方法对分离点的预测较实验值晚，且升力系数C_L均高于实验值。采用同一种湍流模型时，2种形式网格的预测结果基本一致。由此可以看出，修正的k-ω SST湍流模型对静态气动力计算的影响大于网格运动形式的影响，采用修正的k-ω SST湍流模型预测的结果更接近实验值。

2.4动态特性分析

采用2种振荡工况对动态计算结果进行对比，升力系数C_L的迟滞效应见图7。未发生失速时，翼型攻角变化过程为α=1.5°+1.5°sin（ωt），重叠网格与滑移网格预测的结果一致。湍流模型对结果稍有影响，在相同网格下，修正模型的升力系数C_L小于默认湍流模型的计算结果。在失速区内，即攻角变化过程为α=15.6°+1.8°sin（ωt）时，不同网格运动形式与湍流模型的预测结果均有不同。具体来说，采用默认湍流模型时，滑移网格与重叠网格仅在最大和最小攻角时升力系数C_L一致，其他攻角下滑移网格的升力系数C_L高于重叠网格。采用修正的湍流模型时，在上仰过程中，最小攻角升至14°时，滑移网格与重叠网格预测结果一致;在下俯过程中，α<16.5°时，滑移网格与重叠网格预测结果一致;在其余攻角下，重叠网格计算的升力系数CL小于滑移网格的结果。

上述分析表明，采用数值计算分析翼型动态失速时，湍流模型和网格运动形式对结果均有影响，其中湍流模型的影响大于网格运动形式。2种湍流模型的不同点在于对湍流黏性模型的选择不同。OpenFOAM默认采用的k-ω SST湍流模型对湍流黏性进行简化，修正k-ω SST湍流模型对黏性系数求解更准确，因此计算结果更可靠。

在α=15.6°+1.8°sin（ωt）振荡周期中，网格运动形式和湍流模型对翼型表面压力系数C_p的影响见图8。由此可以看出，表面压力系数Cp的差异主要在吸力面。采用k-ω SST和滑移网格計算的结果具有更高的吸力峰，导致预测的升力系数C_L较高。采用修正的k-ω SST湍流模型时，攻角分别在下行14°、上行14°和上行16°时，滑移网格和重叠网格预测的表面压力系数C_p一致。采用k-ω SST湍流模型时，在前缘至分离点处，重叠网格预测的表面压力系数Cp低于滑移网格的。下行16°时，采用修正的k-ω SST湍流模型和重叠网格计算的压力系数Cp曲线整体上移。采用修正湍流模型和滑移网格时，分离点较其他工况前移。

在α=15.6°+1.8°sin（ωt）振荡周期内，距尾缘1倍弦长处垂直线上的速度分布见图9（纵坐标为无量纲垂直距离，横坐标为无量纲速度）。下行过程的速度变化大于上行过程;除下行16°外，采用修正k-ω SST湍流模型时，2种网格运动形式对尾涡中心速度的预测一致。在速度从尾涡中心至周围流场恢复的过程中，重叠网格预测的速度梯度明显小于滑移网格。

3结束语

基于OpenFOAM，对FFA-W3-241翼型展开非定常模拟方法研究。对动态计算所采用的时间步长和计算周期进行收敛性分析，认为在一个振荡周期内计算3×360步能够保证计算精度，在该时间步长下，3个运动周期后可达到收敛状态。

研究OpenFOAM默认的k-ω SST湍流模型和修正的k-ω SST模型在采用不同网格运动形式时对翼型气动性能的影响，发现在静态计算过程中，修正湍流模型对静态气动力计算的影响大于网格运动形式的影响。在动态计算过程中，附着流状态下，湍流模型和网格运动形式对结果影响程度较弱;失速区域内，重叠网格计算的升力系数小于滑移网格的结果，差异主要出现在上行过程的中后段和下行过程的开始阶段。湍流模型影响较大，修正的k-ω SST湍流模型预测的结果更准确。参考文献：

[1]MCCROSKEY W J. Phenomenon of dynamic stall： NASA TM-81264[C/OL]. （2013-08-11）[2021-12-01]. https：//ntrs.nasa.gov/citations/19830012625.

[2]BHAT S S， GOVARDHAN R N. Stall flutter of NACA 0012 airfoil at low Reynolds numbers[J]. Journal of Fluids and Structures， 2013， 41： 166-174. DOI： 10.1016/j.jfluidstructs.2013.04.001.

[3]张彦军， 赵轲， 张同鑫， 等. 雷诺数变化对翼型边界层发展及失速特性的影响[J]. 航空工程进展， 2019， 10（3）： 319-329. DOI： 10.16615/j.cnki.1674-8190.2019.03.004.

[4]朱呈勇， 王同光. 振荡翼型和振荡来流下的动态失速数值研究[J]. 太阳能学报， 2019， 40（9）： 35-42.

[5]OpenFOAM user guide 2.3.1 technical report[EB/OL]. （2021-01-01）[2021-12-01]. http：//www.openfoam.com/documentation/user-guide/.

[6]CHITSOMBOON T， THAMTHAE C. Adjustment of k-ω SST turbulence model for an improved prediction of stalls on wind turbine blades[C]// Proceedings of 2011 Linkping Electronic Conference. Linkping： World Renewable Energy Congress. DOI： 10.3384/ecp110574114.

[7]FUGLSANG P， ANTONIOU I， DAHL K S. Wind tunnel tests of FFA-W3-241， FFA-W3-301 and NACA 63-430 airfoils[EB/OL]. （1998-01-31）[2021-12-01]. https：//www.osti.gov/etdeweb/servlets/purl/314855.

[8]MENTER F. Zonal two equation k-ω turbulence models for aerodynamic flows[C]// Proceedings of 24th Fluid Dynamics， Plasmadynamics， and Lasers Conference， 1993. Orlando： AIAA. DOI： 10.2514/6.1993-2906

[9]MENTER F R. Two-equation eddy-viscosity turbulence models for engineering applications[J]. AIAA Journal， 1994， 32（8）： 1598-1605. DOI： 10.2514/3.12149.（编辑武晓英）