基于重疊網(wǎng)格的翼型動態(tài)失速數(shù)值模擬

張團元,夏潤澤,劉書巖,員?，|

(1. 海軍航空大學,山東煙臺264001;2.南京航空航天大學機械結(jié)構(gòu)力學及控制國家重點實驗室,江蘇南京210016)

1 引言

動態(tài)失速是指在流場中進行周期振蕩或臨界機動動作等非定常運動的機翼/旋翼周圍發(fā)生氣流分離,從而造成動態(tài)遲滯的現(xiàn)象[1]。該現(xiàn)象在殲擊機的過失速機動、直升機的槳葉旋轉(zhuǎn)、風力機葉片的偏轉(zhuǎn)等流動中非常普遍。其主要特征是翼型周圍發(fā)生非定常分離流動,使得翼型發(fā)生失速直至氣流重新恢復附著流動。動態(tài)失速現(xiàn)象的流動相較于靜態(tài)失速而言,周圍流場特征表現(xiàn)的更加復雜;升力、阻力、力矩特性的變化更為突出。在靜態(tài)失速問題中,翼型迎角達到臨界值時,由于氣流分離氣動力發(fā)生突變。但由于非定常分離渦及大尺度漩渦結(jié)構(gòu)的存在,動態(tài)失速會表現(xiàn)出明顯的氣動力遲滯特性,氣動力系數(shù)也極大的偏離了靜態(tài)值,所以進行動態(tài)失速問題的相關(guān)研究變得尤為重要。

當前主要通過試驗及仿真的方法研究動態(tài)失速問題。試驗方法是在風洞中對縮比模型進行吹風試驗得到氣動力、壓力等,進而直觀認識到動態(tài)過程中的流動延遲現(xiàn)象。但其試驗的前期準備和具體風洞試驗過程非常復雜、成本高昂以及周期時間較長,并只能針對特定工況以及原理性的研究,很難捕捉到模型在振蕩過程中脫體渦的運動形式。數(shù)值模擬計算可以克服以上試驗中存在的不足和難題,針對模型進行多種工況下的動態(tài)失速模擬以及脫體渦整個運動過程的捕捉,該方法在動態(tài)失速問題的研究中已廣泛應用和發(fā)展。

半經(jīng)驗模型方法和計算流體力學(CFD)方法是數(shù)值方法中研究非定常動態(tài)失速問題的兩條重要途徑。國外學者已經(jīng)發(fā)展了一系列的半經(jīng)驗動態(tài)失速理論模型,如Gangwani模型[2]、ONERA模型[3]、Beddoes-Leishman(B-L)模型[4]。其中Gangwani模型是非常早期的半經(jīng)驗模型,而ONERA模型和Leishman-Beddoes(B-L)模型在二維氣動力預測上較為成熟且應用廣泛。以上兩種方法都可以達到工程應用精度,并且大大縮短了計算時間。但半經(jīng)驗模型是基于特定模型試驗數(shù)據(jù)實現(xiàn)的,不具有普適性。此外半經(jīng)驗模型方法僅在氣動力模擬上具有較高的精度和效率,這種方法無法模擬出流場中復雜的分離現(xiàn)象,因此對于復雜流動的機理研究便無法更加深入。

計算流體力學(CFD)方法能夠捕捉到動態(tài)失速過程中渦具體運動形式的細節(jié),是動態(tài)失速問題研究的重要手段。國外學者L.Dubuc[5]等人采用隱式雙時間步長方法求解Euler方程研究了NACA0012翼型跨音速小幅振蕩的非定常運動,準確模擬了翼型淺失速狀態(tài)下的遲滯效應。之后通過TFI方法開發(fā)一種網(wǎng)格變形算法,并通過襟翼的強迫振蕩數(shù)值模擬驗證了該方法的有效性[6]。國內(nèi)學者王軍利[7]采用改進的線性彈簧方法,同樣是求解Euler方程計算了二維翼型以及三維機翼的非定常氣動力問題。張兵[8]等跟據(jù)幾何變形關(guān)系和網(wǎng)格的分塊處理以及通過彈簧類比法和修正的TFI方法分別計算角點位移和子網(wǎng)格內(nèi)部結(jié)點位移,形成了具有彈簧-TFI混合特征的動網(wǎng)格新方法。滿洪海[9]以及雷延生[10]舍棄網(wǎng)格正交性,采用非結(jié)構(gòu)網(wǎng)格通過彈簧近似法和局部網(wǎng)格重構(gòu)方法實現(xiàn)網(wǎng)格變形。目前國內(nèi)外在針對翼型非定常運動的研究中使用的動網(wǎng)格處理方法主要為網(wǎng)格變形法,如代數(shù)法、無限插值方法(TFI)、彈簧法等等[11]。其中TFI方法是一種基于多變量的代數(shù)插值方法,并且為保證網(wǎng)格的正交性只適用于結(jié)構(gòu)網(wǎng)格下的小變形情況;彈簧法即將組成單元的線看作彈簧,通過與邊長相關(guān)的公式計算得到彈簧剛度系數(shù),從而將整個求解域構(gòu)成一個彈簧網(wǎng)格。在邊界節(jié)點運動后,求解網(wǎng)格節(jié)點處的靜力平衡方程來計算新節(jié)點的位置,但其未考慮剪切力的作用,對于大變形來說容易出現(xiàn)網(wǎng)格畸變,從而產(chǎn)生負體積[12]。計算網(wǎng)格的正交性、光滑性以及近翼面貼體網(wǎng)格的布置是評價網(wǎng)格質(zhì)量的重要因素,此外變形網(wǎng)格的方法需要滿足幾何守恒定律,這些都是影響計算精度的要素[13]。

針對網(wǎng)格變形及動態(tài)失速特性等問題,為避免網(wǎng)格變形帶來的正交性下降和負體積的產(chǎn)生,本文采用重疊網(wǎng)格的方法及穩(wěn)定的Roe格式進行動態(tài)失速數(shù)值模擬。驗證了網(wǎng)格及計算方法的有效性;最后,根據(jù)數(shù)值模擬結(jié)果,分析了大迎角動態(tài)失速過程中脫體渦產(chǎn)生、分離以及再附著的運動過程以及翼型的氣動特性。

2 模型及計算網(wǎng)格

2.1 物理模型

本文選用NACA0012對稱翼型作為研究對象,其外形及俯仰運動形式如圖1所示。

圖1 NACA0012翼型運動形式示意圖

翼型的運動形式如圖1所示,繞翼型1/4弦線進行俯仰振蕩,其迎角運動變化的規(guī)律為

α=α0+α1sin(kτ)

(1)

式中,α0是平均迎角,α1是幅值迎角。k是減縮頻率,其值定義為

k=πfc/u∞

(2)

τ是無量綱化時間——減縮時間:

τ=2u∞t/c

(3)

式中,u∞——遠前方來流速度;t——時間;f——振蕩頻率;c——翼型參考弦長?；谙议L的雷諾數(shù)由下式定義:

Re=ρu∞c/μ

(4)

式中,ρ——遠前方來流空氣密度;μ——動力黏度。

式(1)中,平均迎角α0是遠前方來流的角度,幅值迎角α1為翼型的俯仰振蕩角度,兩項相加形成總迎角α的變化。

2.2 重疊網(wǎng)格

目前應用比較廣泛的彈簧網(wǎng)格方法以及TFI插值方法雖然可以有效的避免負體積和畸形網(wǎng)格的出現(xiàn),其對于小幅度的振蕩運動是具有一定的準確度的,但當大變形和大幅度振蕩的翼型運動是很難保證網(wǎng)格正交性和光滑性,產(chǎn)生網(wǎng)格畸變引起負體積導致無法計算。重疊網(wǎng)格在不改變網(wǎng)格形狀的前提下,允許不同區(qū)域的網(wǎng)格進行獨立并行計算,并通過插值的方式進行耦合。所以,本文采用重疊網(wǎng)格技術(shù)并通過對翼型動態(tài)失速問題進行研究與分析。

正因如此，GE Digital的“三步走”發(fā)展策略，即GE For GE、GE For Customers、GE For World是很務實的。在邏輯上，前兩步都正確，但第三步踩了“知識壁壘”的紅線，“For World”哪是那么容易？這中間恐怕省略了“For Other Industries”“For USA”等重要步驟，以及為了驗證和優(yōu)化這些重要步驟所必需的長期沉淀與反復打磨。GE Digital多跨出去的那一步，就是導致其業(yè)務停擺的短板——知識壁壘。

如圖2所示,重疊網(wǎng)格由前景網(wǎng)格和背景網(wǎng)格兩部分組成,通過數(shù)據(jù)插值將兩部分網(wǎng)格合并后進行數(shù)值模擬分析。翼型周圍的前景網(wǎng)格采用C-H型結(jié)構(gòu)網(wǎng)格,以保證其良好的貼體性和正交性。背景網(wǎng)格采用正交性非常好的H型結(jié)構(gòu)網(wǎng)格。計算網(wǎng)格形式及邊界條件處理如圖2所示。對翼型物面周圍網(wǎng)格進行了加密處理,同時通過設(shè)置第一層網(wǎng)格高度保證y+≈1,為求解雷諾平均N-S方程提供捕捉粘性和大迎角氣流分離提供良好的條件。

圖2 NACA0012翼型重疊網(wǎng)格及局部放大示意圖

本文不涉及翼型結(jié)構(gòu)的彈性變形,所以動網(wǎng)格方面通過Fluent的用戶自定義函數(shù)(User-defined functions)來控制重疊網(wǎng)格中前景區(qū)域的運動,從而實現(xiàn)控制翼型的簡諧振蕩運動。

3 離散格式

3.1 控制方程

為較好的模擬動態(tài)失速中的流動分離現(xiàn)象,控制方程采用對附著邊界層湍流和分離湍流模擬具有良好精度的SST k-ω模型[14]。該模型借鑒了Johnson-King模型的思想,在計算中考慮到流動過程中的不平衡作用,得到的翼型動態(tài)失速的氣動性能更加準確,適合于邊界層流動、分離和轉(zhuǎn)捩,因而采用該模型進行翼型動態(tài)失速的模擬分析。

在該模型中,渦粘性定義為

(5)

式中,Ω是渦量絕對值,Ω=|?U/?y|。F2是混合函數(shù),表示為

(6)

其中,湍動能k運輸方程和湍流比耗散率ω方程的具體形式是

τtijSij-β*ρωk

(7)

(8)

式(8)中,右端最后一項為交錯擴散項。式中的生成項Pω≈γρΩ2。同時,識別函數(shù)F1是能夠?qū)崿F(xiàn)k-ε模型與k-ω模型之間切換,其定義為

F1=tanh{min[max(Γ1,Γ3),Γ2]}

(9)

Γ1,Γ2,Γ3可表示為

(10)

該模型下的參數(shù)β,γ,σk,σω可以表述為

φ=F1φ1+(1-F1)φ2(φ={σk,σω,β,γ})

(11)

其中,φ1和φ2分別適用于處理邊界層內(nèi)外流動,給定的系數(shù)值為

σk1=0.85,σω1=0.5,β1=0.075

σk2=1.0,σω2=0.856,β2=0.0828

(12)

γ可表示為

(13)

其它常數(shù)為

a1=0.31,β*=0.09,κ=0.41

(14)

SST k-ω湍流方程因其在預測動態(tài)失速迎角及氣動力計算方面有較好的精度和效率[15],所以該方程目前應用較為廣泛。本文在Fluent軟件中,采用密度基求解器利用該湍流方程求解基于雷諾平均的N-S方程,從而得到俯仰振蕩過程中的氣動力及流場特性。

3.2 空間離散格式

目前,迎風(FDS、FVS)格式、中心格式、TVD格式和ENO格式是CFD計算中常用的格式。本文所采用的格式就是Godunov類求解器的FDS-ROE格式。

該方法是通過構(gòu)造一個合理的常矩陣將非線性的Riemann間斷分解問題轉(zhuǎn)化為線性問題,其轉(zhuǎn)化后的線性方程為

(15)

(16)

(17)

其中,k為自由參數(shù),代表迎風格式的精度。Δ+、Δ-分別為向前、向后差分算子。

Roe格式一般被稱為通量差分分裂(Flux Difference Splitting)格式,相對于其它計算格式,這種方法對激波和接觸間斷都具有較高的分辨率,同時也具有較好的穩(wěn)定性。該格式已成為目前應用最廣、評價最高的CFD空間離散格式之一。

4 非定常動態(tài)失速數(shù)值模擬

采用上述計算方法,通過小幅振蕩及大迎角失速開展對動態(tài)失速特性的模擬研究,并將數(shù)值模擬結(jié)果與試驗值[5],[17]進行對比分析,驗證了計算方法和動網(wǎng)格的有效性。

計算所涉及的初始條件如下表所示。

4.1 小幅振蕩非定常氣動力計算

翼型做小幅振蕩運動的計算條件為上表中的Case1所示。該狀態(tài)的試驗值在文獻[5]給出。其它計算條件為:基于密度基求解器,采用壓力遠場邊界,重疊區(qū)域邊界設(shè)置為overset,壁面條件為無滑移邊界條件,翼型繞其1/4弦線做正弦運動。計算過程為首先計算翼型在平均迎角下的定常流場,將定常流場的計算值作為非定常計算的初始值,以保證瞬態(tài)計算結(jié)果的準確性。計算結(jié)果如圖3所示。

圖3 NACA0012翼型非定常氣動力系數(shù)曲線

圖3分別給出了翼型運動過程中升力系數(shù)和力矩系數(shù)的變化規(guī)律。從圖中可以看出,動態(tài)過程中升力與力矩系數(shù)與靜態(tài)時明顯不一致,隨著迎角的變化均表現(xiàn)出非定常遲滯環(huán)現(xiàn)象。數(shù)值模擬的計算結(jié)果中除力矩系數(shù)中個別值以外,其它同試驗值基本一致,驗證了該方法對翼型小幅振蕩下非定常氣動力計算結(jié)果的精度。

4.2 大迎角動態(tài)失速特性分析

采用以上所述的數(shù)值方法模擬NACA0012翼型大迎角動態(tài)失速過程,試驗條件和結(jié)果在文獻[17]給出。算例的來流條件在表1的Case2中給出。雷諾數(shù)Re=3.45×106,靜溫T=288.15K,對應的黏度μ=1.7894×10-5Pa·s。其它條件設(shè)置和4.1節(jié)計算參數(shù)保持一致。

表1 NACA0012翼型數(shù)值模擬的計算條件

圖4為計算結(jié)果中升力系數(shù)、阻力系數(shù)以及俯仰力矩系數(shù)隨迎角變化的曲線,為保證計算結(jié)果的準確性同樣與試驗值進行了對比。在圖中可以看出,在峰值處計算結(jié)果的絕對值略高于試驗結(jié)果的絕對值,這是由于在動態(tài)失速問題中求解雷諾平均N-S方程時,對脫落渦所誘導的升力有較高的估算[18]。除峰值外,氣動力等系數(shù)計算值與試驗值吻合較好,能夠較好的模擬低馬赫數(shù)下的翼型大迎角動態(tài)失速狀態(tài)下的氣動特性。

圖4 NACA0012翼型氣動力系數(shù)遲滯回線

圖5為NACA0012翼型在深失速俯仰振蕩中翼型表面處不同迎角的流線圖,較好的反映了翼型動態(tài)失速過程中前緣渦的產(chǎn)生、發(fā)展和脫落。對于NACA0012翼型而言,α在15°左右即會發(fā)生靜態(tài)失速[19]。而動態(tài)失速的顯著特點為失速迎角明顯后移。

圖5 NACA0012翼型的動態(tài)流場流線

通過圖5可以看到,在α=15.47°↑時(↑和↓分別代表翼型上仰和下俯狀態(tài)),翼型周圍的流動為完全附著流動。翼型在上仰過程中,迎角逐漸增大,法向力系數(shù)也隨之增大,同時翼型周圍的流動開始出現(xiàn)變化,當迎角上仰角度大于16°時,如圖5(b)所示開始在翼型前緣開始出現(xiàn)脫體渦,并且伴隨著后緣渦的產(chǎn)生。隨著迎角繼續(xù)增大,流動在翼型前緣處發(fā)生分離,即前緣產(chǎn)生脫體渦,并向下游發(fā)展逐步演化成向下脫落的渦。當翼型運動到25°并開始作下俯運動時,翼型周圍相對來流速度變大,前緣脫體渦和氣流分離現(xiàn)象逐漸加劇,前緣渦沿翼型逐漸向后發(fā)展直至脫落,翼型的上下表面氣流完全分離,導致升力系數(shù)驟然下降。與此同時,由于前緣渦往后緣發(fā)展并發(fā)生分離,導致壓力中心后移,所以低頭力矩系數(shù)會出現(xiàn)急劇上升現(xiàn)象。隨著翼型下俯運動的繼續(xù),前緣再次產(chǎn)生脫體渦,所以圖4(a)中升力系數(shù)會呈現(xiàn)一段上升現(xiàn)象。當運動到小迎角時渦消失,翼型表面恢復為附著流動。俯仰運動過程中正是由于前緣渦的產(chǎn)生、運動、脫落及再附著現(xiàn)象導致翼型在俯仰運動過程中同一迎角時而氣動參數(shù)不同,表現(xiàn)出明顯的遲滯現(xiàn)象,這是動態(tài)失速區(qū)別于靜態(tài)失速的一個顯著特征。

5 結(jié)論

本文基于重疊網(wǎng)格,采用雷諾平均N-S方程對NACA0012翼型進行了動態(tài)失速數(shù)值模擬,并通過動態(tài)流場分析了動態(tài)失速過程中脫體渦的運動特點。通過以上計算分析,主要的出以下結(jié)論:

1)本文所采用的動網(wǎng)格數(shù)值計算方法,由于不涉及網(wǎng)格變形,所以其具有良好的網(wǎng)格正交性,具有一定的精確度,可以有效的模擬翼型動態(tài)失速。該方法除了可以模擬小幅振蕩,同時適合大迎角動態(tài)失速的計算分析。

2)動態(tài)失速過程中氣動特性與靜態(tài)失速中的不同,翼型的氣動力會表現(xiàn)出明顯的遲滯環(huán)現(xiàn)象。通過Case1和Case2可以看出迎角和減縮頻率影響著遲滯環(huán)的形狀和大小。

3)大迎角動態(tài)失速過程中,隨著迎角的增大在前緣產(chǎn)生脫體渦,并隨著上仰運動逐漸發(fā)展直至脫落,氣流完全分離,導致升力驟降,低頭力矩急劇上升。在迎角下俯到一定值時,氣流會在前緣再附。氣流分離與再附著的非定常效應導致翼型振蕩過程中升力不對稱,產(chǎn)生遲滯效應。