翼型動態(tài)失速的非定常模擬方法

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

摘要：? 為探究翼型動態(tài)失速的高可信非定常模擬方法，以FFA-W3-241翼型為研究對象，采用開源計(jì)算流體動力學(xué)求解器OpenFOAM開展翼型動態(tài)失速下的流動模擬。研究重疊網(wǎng)格和滑移網(wǎng)格2種不同網(wǎng)格運(yùn)動形式、2種不同時(shí)間步長、2種不同計(jì)算周期和OpenFOAM默認(rèn)湍流模型與修正的k-ω SST湍流模型對動態(tài)失速過程中翼型氣動力的模擬精度，并對流場結(jié)構(gòu)進(jìn)行分析。結(jié)果表明：修正模型預(yù)測的翼型氣動力和流場特征與實(shí)驗(yàn)值更接近;重疊網(wǎng)格在翼型的動態(tài)失速模擬中更具優(yōu)勢。

關(guān)鍵詞：? OpenFOAM; 動態(tài)失速; 湍流模型; 風(fēng)力機(jī); 翼型; 重疊網(wǎng)格

中圖分類號：? V211.41; TK83文獻(xiàn)標(biāo)志碼：? 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

-基金項(xiàng)目：? 國家數(shù)值風(fēng)洞工程項(xiàng)目（NNW2018-ZT7B14）;國家自然科學(xué)基金（51876063）

作者簡介： 于佳鑫（1993—），女，遼寧建昌人，博士研究生，研究方向?yàn)镃FD可信度分析和不確定性方法，（E-mail）820113965@qq.com

通信作者： 王曉東（1979—），男，北京人，教授，博導(dǎo)，研究方向?yàn)楹Ｉ巷L(fēng)電機(jī)組設(shè)計(jì)，（E-mail）wangxd@ncepu.edu.cn0引言

翼型失速分為靜態(tài)失速和動態(tài)失速。靜態(tài)失速是指迎角超過特定迎角值時(shí)，翼型氣動性能下降的現(xiàn)象。動態(tài)失速是指在迎角快速變化時(shí)，迎角超過靜態(tài)失速角導(dǎo)致翼型產(chǎn)生周期性振蕩的現(xiàn)象[1]。在定常流動中，對于給定的翼型幾何形狀，失速角基本為固定值。然而，當(dāng)翼型在包含靜態(tài)失速角的迎角范圍內(nèi)運(yùn)動時(shí)，最大升力的角度大大增加，且強(qiáng)烈依賴于振蕩的速率和振幅。渦沿翼型由前緣運(yùn)動至尾緣直至離開翼型，此時(shí)升力突然下降，升力和俯仰力矩曲線都產(chǎn)生較大的遲滯回環(huán)。動態(tài)失速研究對翼型的氣動設(shè)計(jì)具有重要意義。

動態(tài)失速的研究方法主要有風(fēng)洞實(shí)驗(yàn)[2]和數(shù)值模擬2種。風(fēng)洞試驗(yàn)要求復(fù)雜、成本高，而且難以捕捉動態(tài)過程中的流體細(xì)節(jié)。受實(shí)驗(yàn)條件限制，研究人員大多采用數(shù)值模擬獲取更詳細(xì)的動態(tài)失速信息。動態(tài)失速是強(qiáng)非定常流動過程，準(zhǔn)確的計(jì)算流體動力學(xué)（computational fluid dynamics， CFD）模擬難度較大。張彥軍等[3]計(jì)算不同雷諾數(shù)對動態(tài)失速特性的影響。朱呈勇[4]采用非定常雷諾時(shí)均（unsteady reynolds average navier stocks， URANS）方法研究動態(tài)入流和翼型振蕩對失速特性的影響。

CFD求解器OpenFOAM是一款開源軟件，可以自行修改求解器以應(yīng)對不同物理問題的需求，通過C語言代碼編寫腳本，可使復(fù)雜操作過程更加自動便捷。[5]目前用OpenFOAM求解翼型動態(tài)失速的研究較少，模擬方法有待驗(yàn)證。k-ω SST湍流模型是翼型數(shù)值模擬中廣泛采用的湍流模型。[6]OpenFOAM中默認(rèn)的k-ω SST湍流模型對黏性系數(shù)的求解有一定簡化，針對風(fēng)力機(jī)翼型動態(tài)失速問題，如何適當(dāng)修改湍流模型的系數(shù)值得研究。

翼型的俯仰振蕩是翼型繞特定軸進(jìn)行的旋轉(zhuǎn)振蕩，因此數(shù)值計(jì)算時(shí)需要網(wǎng)格按照相同的方式運(yùn)動。OpenFOAM實(shí)現(xiàn)網(wǎng)格運(yùn)動的方式主要有3種，分別為動網(wǎng)格（Dynamic Mesh）、滑移網(wǎng)格（Sliding Mesh）和重疊網(wǎng)格（Overset Mesh）。動網(wǎng)格可根據(jù)設(shè)定的運(yùn)動，使運(yùn)動物體附近的網(wǎng)格發(fā)生實(shí)時(shí)變形與重新生成。翼型尾緣幾何尺寸較小，為得到較準(zhǔn)確的計(jì)算結(jié)果，網(wǎng)格尺寸往往也較小。采用動網(wǎng)格時(shí)，尾緣附近的網(wǎng)格距離旋轉(zhuǎn)軸較遠(yuǎn)，變形較大，在翼型運(yùn)動過程中容易出現(xiàn)負(fù)網(wǎng)格，導(dǎo)致計(jì)算發(fā)散?；凭W(wǎng)格需要建立靜止區(qū)域和運(yùn)動區(qū)域，在二者間設(shè)置交界面，運(yùn)動區(qū)域內(nèi)的網(wǎng)格整體運(yùn)動，不會涉及到網(wǎng)格的變形與重新生成，因此不會產(chǎn)生負(fù)網(wǎng)格?；凭W(wǎng)格是動網(wǎng)格的簡化形式，通常用于往復(fù)運(yùn)動和旋轉(zhuǎn)運(yùn)動。重疊網(wǎng)格將復(fù)雜的流場區(qū)域進(jìn)行分解，每個(gè)區(qū)域內(nèi)獨(dú)立生成高質(zhì)量網(wǎng)格，區(qū)域網(wǎng)格之間有重疊和共享部分，可用于任意運(yùn)動方式。重疊網(wǎng)格綜合動網(wǎng)格和滑移網(wǎng)格的優(yōu)點(diǎn)，在保證物體運(yùn)動準(zhǔn)確的同時(shí)又能保證運(yùn)動過程中的網(wǎng)格質(zhì)量。因此，網(wǎng)格運(yùn)動形式采用滑移網(wǎng)格和重疊網(wǎng)格2種。

為探究翼型動態(tài)失速的高可信非定常模擬方法，本文基于OpenFOAM，采用URANS方法，研究網(wǎng)格運(yùn)動形式和湍流模型對翼型氣動力和流場特征的影響。

1數(shù)值模型與研究方法

1.1風(fēng)洞實(shí)驗(yàn)回顧

文獻(xiàn)[7]開展FFA-W3-241翼型的相關(guān)實(shí)驗(yàn)研究，采用的雷諾數(shù)為1.6×10⁶，湍流強(qiáng)度為1%，翼型弦長為0.60 m，整個(gè)翼型圍繞x/c=0.4以正弦形式運(yùn)動，運(yùn)動形式示意見圖1，其中α_mea為平均攻角，α_amp為攻角振蕩幅值。

式中：f為振蕩頻率;c為弦長;U_∞為自由來流速度;ω為角速度。

根據(jù)實(shí)驗(yàn)結(jié)果，下文數(shù)值模擬選擇2個(gè)平均攻角α_mea，分別為1.5°和15.6°，2個(gè)平均攻角對應(yīng)的攻角振蕩幅值α_amp分別為1.5°和1.8°。

1.2數(shù)值方法

1.2.1計(jì)算域和網(wǎng)格

計(jì)算域和邊界條件見圖2。

滑移網(wǎng)格與重疊網(wǎng)格的靜止域和旋轉(zhuǎn)域的幾何尺寸保持一致。靜止域?yàn)檎叫?，其邊長為40倍弦長（40c），并將翼型旋轉(zhuǎn)軸置于計(jì)算域中心。旋轉(zhuǎn)域是以旋轉(zhuǎn)軸為圓心、半徑為3c的圓。

靜止域的左側(cè)邊界為速度入口，湍流強(qiáng)度與實(shí)驗(yàn)設(shè)置相同，即為1%。翼型表面為光滑無滑移壁面，運(yùn)動形式在dynamicMeshDict中設(shè)置，滑移網(wǎng)格和重疊網(wǎng)格的dynamicFvMesh分別設(shè)置為dynamicMotionSolverFvMesh類型和dynamicOverSet FvMesh類型，運(yùn)動函數(shù)為旋轉(zhuǎn)振蕩函數(shù)oscillatingRotatingMotion。為實(shí)現(xiàn)平行流動，上、下邊界速度和壓力的邊界條件設(shè)置為零梯度。

網(wǎng)格劃分時(shí)保證y⁺<1，滑移網(wǎng)格與重疊網(wǎng)格的靜止域外側(cè)、旋轉(zhuǎn)域和翼型周向節(jié)點(diǎn)保持一致。靜止域外側(cè)各邊節(jié)點(diǎn)數(shù)為113?；凭W(wǎng)格交界面處節(jié)點(diǎn)數(shù)與翼型周向節(jié)點(diǎn)總數(shù)保持一致，并且均勻分布。2種運(yùn)動形式的網(wǎng)格見圖3。翼型上、下表面各分布200個(gè)節(jié)點(diǎn)，法向設(shè)置115個(gè)節(jié)點(diǎn)，壁面網(wǎng)格的法向膨脹率為1.1。FFA-W3-241翼型尾緣為鈍尾緣，設(shè)置50個(gè)節(jié)點(diǎn)。

此外，繪制3種拓?fù)浣Y(jié)構(gòu)相同的網(wǎng)格進(jìn)行網(wǎng)格獨(dú)立性驗(yàn)證。各邊的節(jié)點(diǎn)數(shù)與上述網(wǎng)格成1.5倍關(guān)系，網(wǎng)格參數(shù)見表1。

1.2.2湍流模型

采用k-ω SST湍流模型計(jì)算URANS方程。OpenFOAM的標(biāo)準(zhǔn)k-ω SST模型是一種全湍流渦黏模型[8]，湍流黏性μ_t的計(jì)算公式為

式中：a₁和b₁為封閉系數(shù)，a₁=0.31，b₁=1.0;k為湍動能;ω為湍流耗散率;S為應(yīng)變率幅值;F₂₃參數(shù)參考OpenFOAM文檔[5]取值。

1994年，MENTER[9]對該模型進(jìn)行改進(jìn)，改進(jìn)模型的湍流黏性項(xiàng)為

本文在OpenFOAM中加入改進(jìn)模型，并與默認(rèn)模型進(jìn)行動態(tài)失速仿真的數(shù)值對比。

2結(jié)果分析

2.1網(wǎng)格無關(guān)性驗(yàn)證

在3種網(wǎng)格精度下，滑移網(wǎng)格和重疊網(wǎng)格計(jì)算得到的升力系數(shù)C_L見圖4。在相同網(wǎng)格精度下，滑移網(wǎng)格與重疊網(wǎng)格的結(jié)果一致。當(dāng)攻角α<12°時(shí)，3種網(wǎng)格精度對升力系數(shù)C_L的預(yù)測一致;當(dāng)攻角α>12°時(shí)，粗網(wǎng)格計(jì)算結(jié)果與其他2種網(wǎng)格預(yù)測結(jié)果有些差別，中等網(wǎng)格與密網(wǎng)格的預(yù)測表現(xiàn)一致?？紤]到計(jì)算成本，采用中等網(wǎng)格進(jìn)行動態(tài)失速計(jì)算。

2.2動態(tài)時(shí)間步長和計(jì)算周期驗(yàn)證

為得到穩(wěn)定的計(jì)算結(jié)果并降低計(jì)算成本，進(jìn)行時(shí)間步長驗(yàn)證。以平均攻角α_mea=15.6°、攻角振蕩幅值α_amp=1.8°為例，設(shè)計(jì)4種時(shí)間步長，一個(gè)振蕩周期內(nèi)的步數(shù)分別為1×360、2×360、3×360和4×360步。第2～5個(gè)運(yùn)動周期內(nèi)各時(shí)間步長的升力系數(shù)CL時(shí)間歷程見圖5。

由圖5可以看出，采用1×360和2×360步的結(jié)果與采用3×360步與4×360步的結(jié)果相差較大，3×360步與4×360步的結(jié)果幾乎重合，可達(dá)到收斂。

選取時(shí)間步數(shù)為3×360，進(jìn)行計(jì)算周期數(shù)收斂性驗(yàn)證，計(jì)算穩(wěn)定后相鄰2個(gè)周期氣動力升力系數(shù)

C_L和阻力系數(shù)C_D計(jì)算結(jié)果的標(biāo)準(zhǔn)差見表2。由此可以看出，3個(gè)周期后標(biāo)準(zhǔn)差不再變化，可達(dá)到收斂。

綜上所述，考慮到計(jì)算時(shí)間成本，采用3×360步數(shù)對應(yīng)的時(shí)間步長，并選取運(yùn)動3個(gè)周期后的結(jié)果進(jìn)行分析。

2.3穩(wěn)態(tài)特性分析

利用滑移網(wǎng)格和重疊網(wǎng)格分別采用OpenFOAM默認(rèn)的k-ω SST模型和修正的k-ω SST模型計(jì)算升力系數(shù)C_L，結(jié)果見圖6。攻角小于10°時(shí)，空氣流動處于附著流動狀態(tài)，未出現(xiàn)失速現(xiàn)象，2種形式的網(wǎng)格和2種湍流模型預(yù)測的結(jié)果與實(shí)驗(yàn)一致。失速以后，即攻角大于10°后，數(shù)值計(jì)算方法對分離點(diǎn)的預(yù)測較實(shí)驗(yàn)值晚，且升力系數(shù)C_L均高于實(shí)驗(yàn)值。采用同一種湍流模型時(shí)，2種形式網(wǎng)格的預(yù)測結(jié)果基本一致。由此可以看出，修正的k-ω SST湍流模型對靜態(tài)氣動力計(jì)算的影響大于網(wǎng)格運(yùn)動形式的影響，采用修正的k-ω SST湍流模型預(yù)測的結(jié)果更接近實(shí)驗(yàn)值。

2.4動態(tài)特性分析

采用2種振蕩工況對動態(tài)計(jì)算結(jié)果進(jìn)行對比，升力系數(shù)C_L的遲滯效應(yīng)見圖7。未發(fā)生失速時(shí)，翼型攻角變化過程為α=1.5°+1.5°sin（ωt），重疊網(wǎng)格與滑移網(wǎng)格預(yù)測的結(jié)果一致。湍流模型對結(jié)果稍有影響，在相同網(wǎng)格下，修正模型的升力系數(shù)C_L小于默認(rèn)湍流模型的計(jì)算結(jié)果。在失速區(qū)內(nèi)，即攻角變化過程為α=15.6°+1.8°sin（ωt）時(shí)，不同網(wǎng)格運(yùn)動形式與湍流模型的預(yù)測結(jié)果均有不同。具體來說，采用默認(rèn)湍流模型時(shí)，滑移網(wǎng)格與重疊網(wǎng)格僅在最大和最小攻角時(shí)升力系數(shù)C_L一致，其他攻角下滑移網(wǎng)格的升力系數(shù)C_L高于重疊網(wǎng)格。采用修正的湍流模型時(shí)，在上仰過程中，最小攻角升至14°時(shí)，滑移網(wǎng)格與重疊網(wǎng)格預(yù)測結(jié)果一致;在下俯過程中，α<16.5°時(shí)，滑移網(wǎng)格與重疊網(wǎng)格預(yù)測結(jié)果一致;在其余攻角下，重疊網(wǎng)格計(jì)算的升力系數(shù)CL小于滑移網(wǎng)格的結(jié)果。

上述分析表明，采用數(shù)值計(jì)算分析翼型動態(tài)失速時(shí)，湍流模型和網(wǎng)格運(yùn)動形式對結(jié)果均有影響，其中湍流模型的影響大于網(wǎng)格運(yùn)動形式。2種湍流模型的不同點(diǎn)在于對湍流黏性模型的選擇不同。OpenFOAM默認(rèn)采用的k-ω SST湍流模型對湍流黏性進(jìn)行簡化，修正k-ω SST湍流模型對黏性系數(shù)求解更準(zhǔn)確，因此計(jì)算結(jié)果更可靠。

在α=15.6°+1.8°sin（ωt）振蕩周期中，網(wǎng)格運(yùn)動形式和湍流模型對翼型表面壓力系數(shù)C_p的影響見圖8。由此可以看出，表面壓力系數(shù)Cp的差異主要在吸力面。采用k-ω SST和滑移網(wǎng)格計(jì)算的結(jié)果具有更高的吸力峰，導(dǎo)致預(yù)測的升力系數(shù)C_L較高。采用修正的k-ω SST湍流模型時(shí)，攻角分別在下行14°、上行14°和上行16°時(shí)，滑移網(wǎng)格和重疊網(wǎng)格預(yù)測的表面壓力系數(shù)C_p一致。采用k-ω SST湍流模型時(shí)，在前緣至分離點(diǎn)處，重疊網(wǎng)格預(yù)測的表面壓力系數(shù)Cp低于滑移網(wǎng)格的。下行16°時(shí)，采用修正的k-ω SST湍流模型和重疊網(wǎng)格計(jì)算的壓力系數(shù)Cp曲線整體上移。采用修正湍流模型和滑移網(wǎng)格時(shí)，分離點(diǎn)較其他工況前移。

在α=15.6°+1.8°sin（ωt）振蕩周期內(nèi)，距尾緣1倍弦長處垂直線上的速度分布見圖9（縱坐標(biāo)為無量綱垂直距離，橫坐標(biāo)為無量綱速度）。下行過程的速度變化大于上行過程;除下行16°外，采用修正k-ω SST湍流模型時(shí)，2種網(wǎng)格運(yùn)動形式對尾渦中心速度的預(yù)測一致。在速度從尾渦中心至周圍流場恢復(fù)的過程中，重疊網(wǎng)格預(yù)測的速度梯度明顯小于滑移網(wǎng)格。

3結(jié)束語

基于OpenFOAM，對FFA-W3-241翼型展開非定常模擬方法研究。對動態(tài)計(jì)算所采用的時(shí)間步長和計(jì)算周期進(jìn)行收斂性分析，認(rèn)為在一個(gè)振蕩周期內(nèi)計(jì)算3×360步能夠保證計(jì)算精度，在該時(shí)間步長下，3個(gè)運(yùn)動周期后可達(dá)到收斂狀態(tài)。

研究OpenFOAM默認(rèn)的k-ω SST湍流模型和修正的k-ω SST模型在采用不同網(wǎng)格運(yùn)動形式時(shí)對翼型氣動性能的影響，發(fā)現(xiàn)在靜態(tài)計(jì)算過程中，修正湍流模型對靜態(tài)氣動力計(jì)算的影響大于網(wǎng)格運(yùn)動形式的影響。在動態(tài)計(jì)算過程中，附著流狀態(tài)下，湍流模型和網(wǎng)格運(yùn)動形式對結(jié)果影響程度較弱;失速區(qū)域內(nèi)，重疊網(wǎng)格計(jì)算的升力系數(shù)小于滑移網(wǎng)格的結(jié)果，差異主要出現(xiàn)在上行過程的中后段和下行過程的開始階段。湍流模型影響較大，修正的k-ω SST湍流模型預(yù)測的結(jié)果更準(zhǔn)確。參考文獻(xiàn)：

[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]張彥軍， 趙軻， 張同鑫， 等. 雷諾數(shù)變化對翼型邊界層發(fā)展及失速特性的影響[J]. 航空工程進(jìn)展， 2019， 10（3）： 319-329. DOI： 10.16615/j.cnki.1674-8190.2019.03.004.

[4]朱呈勇， 王同光. 振蕩翼型和振蕩來流下的動態(tài)失速數(shù)值研究[J]. 太陽能學(xué)報(bào)， 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.（編輯武曉英）