翼型動態失速的非定常模擬方法

于佳鑫　陳江濤　王曉東　吳曉軍　康順

摘要：? 為探究翼型動態失速的高可信非定常模擬方法，以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.（編輯武曉英）