旋翼干擾流場數(shù)值模擬中的新型Blend低速預(yù)處理方法

陳龍，夏健，田書玲

南京航空航天大學(xué) 航空宇航學(xué)院，南京 210016

旋翼/機身干擾和直升機/艦船干擾[1-2]等旋翼干擾非定常流場非常復(fù)雜。直升機機身和艦船流場會對旋翼氣動力和操穩(wěn)特性產(chǎn)生影響，旋翼尾跡對機身氣動特性和艦面流場也有較大影響。流場中艦船及機身附近速度很低，是不可壓流動。旋翼槳尖附近卻處于跨聲速狀態(tài)并可能存在激波，是可壓縮流動。因此在研究旋翼干擾流場時，需要考慮可壓/不可壓非定常流動的統(tǒng)一處理，應(yīng)采用基于可壓縮Navier-Stokes方程的低速預(yù)處理方法。旋翼干擾非定常流場的模擬若不采用低速預(yù)處理，過大的剛性和耗散會帶來嚴重的收斂性和準確性的問題。低速預(yù)處理方法的核心問題是預(yù)處理矩陣的設(shè)計[3-6]。預(yù)處理矩陣對預(yù)處理方法的收斂性影響非常大，理想的預(yù)處理矩陣應(yīng)具有特征值相對集中、消除剛性等特點。一般定常雷諾平均Navier-Stokes方程的預(yù)處理采用定常預(yù)處理矩陣[7-8]，非定常雷諾平均Navier-Stokes(URANS)方程的預(yù)處理采用非定常預(yù)處理矩陣[9-10]。

如果采用定常低速預(yù)處理模擬非定常問題，在雙時間步長內(nèi)迭代的收斂性問題會很嚴重[11]，甚至需要比無低速預(yù)處理更多的迭代。非定常低速預(yù)處理是在定常低速預(yù)處理的基礎(chǔ)上引入矩陣中的偽時間和物理時間項，消除了剛性，在雙時間步長內(nèi)迭代中具有較好的收斂性。但是，非定常預(yù)處理在進行旋翼/機身干擾和直升機/艦船干擾為代表的旋翼干擾流場模擬時，物理時間步長需要按照旋翼轉(zhuǎn)速確定。該時間步長較機身或艦船特征時間小很多，導(dǎo)致非定常低速預(yù)處理矩陣中的非定常截斷參數(shù)大于1，這時非定常低速預(yù)處理自動關(guān)閉[12]。也就是說，在處理此類問題時，非定常低速預(yù)處理會自動退化到無預(yù)處理形式，收斂性和預(yù)測精度均無法得到保證。采用定常和非定常低速預(yù)處理在旋翼干擾流場數(shù)值模擬時遇到了兩難困境。鑒于此，Potsdam等[12]將非定常低速預(yù)處理和定常預(yù)處理特征值矩陣進行組合，建立了Blend低速預(yù)處理方法。該方法解決了非定常低速預(yù)處理會自動退化到無預(yù)處理的問題，減小了流場中不可壓區(qū)域的耗散。Sachdev和Hosangadi[13]針對特征值矩陣中定常預(yù)處理部分采用對應(yīng)的定常預(yù)處理矩陣，改進Potsdam的方法為精確的代數(shù)形式，并推廣到AUSM (Advection Upstream Splitting Method)類的通量分裂格式中。

本文在Potsdam和Sachdev的Blend低速預(yù)處理方法基礎(chǔ)上結(jié)合對流項和黏性項量級分析，將定常低速預(yù)處理和非定常低速預(yù)處理特征值矩陣進行組合，建立出一種新型Blend低速預(yù)處理方法。本方法仍采用非定常預(yù)處理矩陣，對定常預(yù)處理和非定常預(yù)處理的通量格式耗散項進行組合得到新的通量格式耗散部分。在南京航空航天大學(xué)自主開發(fā)的OVERU軟件[14]中應(yīng)用該方法，經(jīng)驗證本文發(fā)展的方法相比傳統(tǒng)方法計算效率大幅提高，也具有較好的預(yù)測精度，適用于旋翼干擾非定常流場的數(shù)值模擬，工程應(yīng)用前景廣泛。

1 控制方程與數(shù)值方法

控制方程為任意拉格朗日歐拉(ALE)形式的URANS：

(1)

(2)

式中：ρ為密度；u、v、w為速度；E為內(nèi)能；W為守恒變量；xj為坐標；t為時間；Fj和Gj為通量。

采用非結(jié)構(gòu)混合網(wǎng)格有限體積方法求解控制方程?？臻g離散采用格點格式，可處理多種單元類型。湍流模型可采用Negative-SA一方程模型和Menter的剪切應(yīng)力輸運(SST)兩方程模型。時間離散采用隱式LU-SGS方法。非定常計算采用雙時間步長方法進行二階精度的時間離散。采用隱式重疊網(wǎng)格方法[14]處理網(wǎng)格間的相對運動。

2 低速預(yù)處理方法

預(yù)處理矩陣的建立首先需要將守恒變量W形式的方程轉(zhuǎn)換成原始變量Q的形式，即

(3)

(4)

(5)

將原方程中的雅可比矩陣替換為預(yù)處理矩陣Γ，則預(yù)處理Navier-Stokes方程可寫為

(6)

本文采用Weiss和Smith的預(yù)處理矩陣[15]：

(7)

(8)

式中：Map為定常截斷參數(shù)；c為聲速。

預(yù)處理參數(shù)采用當?shù)伛R赫數(shù)配合全局截斷的方法：

(9)

(10)

式中：Mai為當?shù)伛R赫數(shù)；Ma∞為來流馬赫數(shù)。

為了和預(yù)處理之前的方程保持相同形式，以便于程序編寫，采用守恒形式的預(yù)處理方程為

(11)

稱ΓM-1為守恒形式的預(yù)處理矩陣。定常預(yù)處理矩陣可以直接應(yīng)用到雙時間步長的計算中[10]，但由于其矩陣特征值的病態(tài)，造成了雙時間步長內(nèi)迭代收斂性很差。因此，需要建立非定常預(yù)處理矩陣，引入偽時間，對于雙時間步長方法有

(12)

對該方程物理時間項采用二階向后歐拉格式離散，偽時間項采用一階離散得到離散方程，t是物理時間，τ是偽時間，有

(13)

式中：n表示第n個物理時間步；m表示第m個偽時間步。方程右端項為

(14)

令

(15)

則ΓuM-1為守恒形式的非定常預(yù)處理矩陣[11]。

令

(16)

(17)

Θ′=bΘ-(b-1)ρp

(18)

非定常預(yù)處理參數(shù)的定義為

(19)

(20)

式中：Mau為非定常截斷參數(shù)；L為參考長度；Δt為時間步長。當高Strouhal數(shù)非定常流模擬時，時間步長Δt非常??；當模型的參考長度非常大時，將造成Mau>1，從而導(dǎo)致Map=1，此時非定常低速預(yù)處理為關(guān)閉狀態(tài)。也就是說在處理此類問題時，非定常低速預(yù)處理會自動退化到無預(yù)處理形式，收斂性和預(yù)測精度均無法得到保證。旋翼機身干擾流場就屬于此種情況。同時，對于復(fù)雜算例，特征長度和特征時間很難準確定義，這種不確定性影響著Mau的取值，也影響著非定常低速預(yù)處理的預(yù)測精度和收斂性。

3 Blend低速預(yù)處理方法

為解決非定常截斷參數(shù)Mau>1造成非定常預(yù)處理自動關(guān)閉的問題，Potsdam組合定常低速預(yù)處理和非定常低速預(yù)處理特征值矩陣，從JST(Jameson-Schmidt-Turkel)格式的預(yù)處理形式出發(fā)建立了Blend低速預(yù)處理方法。

文獻[12]認為不可壓流場定常計算時，無預(yù)處理黏性項的速度場和溫度場耗散過大，壓強場耗散不足。通過定常低速預(yù)處理使得黏性項速度場、溫度場和壓強場均為合理的耗散。非定常計算時，問題就變得復(fù)雜的多，尤其是較小的物理時間步長或高Strouhal數(shù)的問題。若采用定常低速預(yù)處理模擬非定常問題，會造成病態(tài)特征值和雙時間步長內(nèi)很差的收斂性，同時人工黏性項中壓強場耗散也會過大。若采用非定常低速預(yù)處理，時間步長較小造成Mau>1，非定常低速預(yù)處理就會自動退化到?jīng)]有預(yù)處理的形式。為擺脫這個兩難困境，Potsdam通過將定常低速預(yù)處理和非定常低速預(yù)處理特征值矩陣組合建立了Blend低速預(yù)處理方法[12]，方法中特征值矩陣為

(21)

(22)

引入選擇矩陣

(23)

式(22)可表示為

(24)

在此基礎(chǔ)上，Sachdev針對特征值矩陣中定常預(yù)處理部分采用對應(yīng)的定常預(yù)處理矩陣改進了Postsdam的方法(式(25)方框項)，這種形式更適合通量差分格式耗散項的構(gòu)造。文獻[13]經(jīng)過驗證后證明Sachdev的方法和Potsdam的方法計算結(jié)果差別極小。Sachdev的Blend低速預(yù)處理中JST格式的人工黏性項為

D1/2=

(25)

量級分析方法是低速預(yù)處理收斂性研究的一個重要手段[3,16]。文獻[3]從數(shù)值分析的角度，進行了對流項和人工黏性項的量級分析，研究了定常預(yù)處理方法對低馬赫數(shù)時Euler方程收斂性的影響，當Ma→0時，未采用預(yù)處理的對流項的量級為

(26)

黏性項的量級為

(27)

當Ma→0時，預(yù)處理后對流項的量級為

(28)

黏性項的量級為

(29)

可見當Ma→0時，無預(yù)處理的對流項和人工黏性項量級并不匹配，有且僅有速度場的黏性項偏大，而在預(yù)處理后各方程量級保持一致。采用文獻[3，16]的量級準則對Potsdam的Blend方法進行量級分析，在Mau>1時預(yù)處理系統(tǒng)的特征值量級為

(30)

(31)

則本文的Blend低速預(yù)處理中選擇矩陣為

(32)

本文的Blend低速預(yù)處理中JST格式的人工黏性項為

D1/2=

(33)

在非定常預(yù)處理的框架下，當Mau>1時非定常預(yù)處理關(guān)閉，會出現(xiàn)格式耗散偏大的問題。而此時本文方法中對流項和黏性項的量級仍然是匹配的，減小了數(shù)值耗散，從而解決非定常低速預(yù)處理自動退化為無預(yù)處理的問題。同時本文的思想也適用于AUSM類格式的改進[13]。

4 計算結(jié)果及分析

4.1 圓柱非定常擾流

本算例模擬圓柱非定常擾流，取很小物理時間步長模擬非定常低速預(yù)處理自動退化為無預(yù)處理的問題，對本文建立方法的收斂性和正確性進行考核。網(wǎng)格采用三角形和四邊形混合網(wǎng)格，約6萬網(wǎng)格點，如圖1所示。計算條件：來流馬赫數(shù)為0.02，雷諾數(shù)ReD=1 000，采用層流計算，圓柱直徑D=0.01 m，時間步長Δt=5×10-6s，非定常截斷參數(shù)Mau=D/(πΔtc)≈1.87。此時非定常預(yù)處理自動退化為無預(yù)處理的形式。

圖2給出了雙時間步長內(nèi)迭代的殘差收斂曲線，可見定常預(yù)處理難以收斂，Blend預(yù)處理收斂性較好，本文方法略快于Sachdev的Blend方法，無/非定常預(yù)處理則因其耗散最大所以收斂速度最快。此處“無/非定常預(yù)處理”表示此狀態(tài)Mau>1非定常預(yù)處理會自動退化到無預(yù)處理的形式，無預(yù)處理和采用非定常預(yù)處理計算結(jié)果是一致的。由于同樣采用了非定常預(yù)處理矩陣，本文方法在Mau>1時預(yù)處理同樣會關(guān)閉，但是改進過的黏性項與對流項量級仍然匹配，數(shù)值耗散也更合理。采樣點(距圓心3.4D)的速度歷程如圖3所示，定常預(yù)處理、本文Blend預(yù)處理和Sachdev的方法幾乎重合，而無/非定常預(yù)處理的過大耗散使其結(jié)果與其他計算結(jié)果差別很大。圖2和圖3中也對比了無/非定常預(yù)處理時采用Roe和HLLC (Harten-Lax-van Leer-Contact)迎風格式，可見無預(yù)處理時兩種格式在此算例中無法收斂。圖4和圖5分別給出了無/非定常預(yù)處理和本文Blend預(yù)處理的渦量云圖，可見無/非定常預(yù)處理的過大耗散使流場中的渦強在脫離圓柱后快速降低，而Blend預(yù)處理的渦量云圖則顯示出Blend低速預(yù)處理的耗散合理，得到了較好的渦量結(jié)果。綜上，兩種Blend預(yù)處理計算結(jié)果一致，但收斂效率較定常預(yù)處理有較大提高。而無/非定常預(yù)處理會耗散過大，計算結(jié)果較差。

圖1 圓柱非定常擾流計算網(wǎng)格Fig.1 Computational grid for unsteady flow around cylinder

圖2 圓柱擾流非定常殘差收斂曲線Fig.2 Convergence history of residual of unsteady flow around cylinder

圖3 采樣點速度歷程Fig.3 Velocity history of sampling point

圖4 圓柱擾流無/非定常低速預(yù)處理計算結(jié)果Fig.4 No/unsteady low speed preconditioning result of flow around cylinder

圖5 圓柱擾流本文Blend低速預(yù)處理計算結(jié)果Fig.5 Low speed preconditioning result of flow around cylinder of the proposed Blend method

4.2 旋翼/機身干擾流場

本算例進行直升機前飛時旋翼機身干擾的數(shù)值模擬，對本文Blend預(yù)處理方法收斂性和正確性進行考核。采用直徑2.5 m的單旋翼直升機基準模型，包括4片槳葉的旋翼模型和Robin機身模型。旋翼實度為0.077，翼型采用NACA0012，弦長0.072 m，幾何扭轉(zhuǎn)0°，旋翼的周期變距為θ0=3.1°，θ1c=-0.6°，θ1s=1.5°，槳盤傾角αs=-5°，前進比為0.1。

圖6為本算例所采用的重疊網(wǎng)格系統(tǒng)，共包含7塊子網(wǎng)格，其中2塊為遠場背景網(wǎng)格，它們與繞Robin直升機機身的網(wǎng)格為靜止網(wǎng)格，2塊背景網(wǎng)格單元數(shù)分別為31萬和163萬，Robin直升機機身網(wǎng)格單元數(shù)為677萬，均為六面體網(wǎng)格；繞每片旋翼葉片的網(wǎng)格隨旋翼運動，網(wǎng)格單元數(shù)為107萬。重疊網(wǎng)格系統(tǒng)總網(wǎng)格數(shù)約1 300萬。物面第1層網(wǎng)格y+≤1。

圖7給出了Robin直升機機身表面網(wǎng)格。湍流模型采用Negative-SA一方程模型。時間步長取旋翼轉(zhuǎn)過0.5°方位角，參考長度為機身長度L=2.5 m，則非定常截斷參數(shù)Mau=L/(πΔtc)≈38.8。此時非定常預(yù)處理自動退化為無預(yù)處理的形式。

圖6 前飛直升機流場模擬重疊網(wǎng)格系統(tǒng)Fig.6 Overset grids system for flowfield simulation of helicopter in forward flight

圖7 機身表面網(wǎng)格Fig.7 Fuselage surface grid

圖8 定常預(yù)處理參數(shù)對單獨機身氣動力的影響Fig.8 Isolated fuselage forces vs steady preconditioning parameters

圖9 單獨機身壓力系數(shù)分布計算值與試驗值對比Fig.9 Comparison between experimental and computed distribution of isolated fuselage pressure coefficient

圖10 單獨機身壓力云圖對比Fig.10 Comparison of pressure contour of isolated fuselage

圖11 定常預(yù)處理參數(shù)對收斂曲線的影響Fig.11 Convergence history vs steady preconditioning parameter

采用無/非定常低速預(yù)處理、定常預(yù)處理和本文Blend預(yù)處理方法開展直升機前飛時旋翼/機身干擾流場的模擬。隨旋翼方位角的變化機身法向力系數(shù)的變化可以反映旋翼對機身氣動力的干擾，圖12和圖13分別為采用不同雙時間步長內(nèi)迭代數(shù)和通量格式對機身法向力系數(shù)CN的影響。從收斂性方面考慮：無/非定常預(yù)處理時，機身法向力在內(nèi)迭代(Nin為迭代次數(shù))10、20和50次時很接近，僅需要約10步內(nèi)迭代即可收斂；定常預(yù)處理時，從20增加到100次內(nèi)迭代后機身法向力仍然沒有收斂，收斂性非常差；本文Blend預(yù)處理，20次和50次內(nèi)迭代的機身法向力幾乎重合，表明其僅需約20次內(nèi)迭代即可收斂，收斂性較定常預(yù)處理大幅改善。從流場預(yù)測精度考慮：充分收斂后定常低速預(yù)處理和Blend低速預(yù)處理計算結(jié)果表現(xiàn)出較好的一致性，無/非定常低速預(yù)處理由于存在較大耗散，JST格式和HLLC格式均得到較平緩的計算結(jié)果，存在較大的誤差。

圖14給出了不同方法的內(nèi)迭代殘差收斂曲線，可見本文Blend預(yù)處理方法的收斂效率略高于Sachdev的Blend預(yù)處理方法，定常預(yù)處理方法后期收斂速度會變慢，無/非定常預(yù)處理時JST格式由于其耗散最大收斂速度最快，HLLC格式則收斂困難。圖15給出了旋翼方位角0°時各方法計算的機身表面壓力分布，定常預(yù)處理方法和本文Blend預(yù)處理方法結(jié)果符合得很好，而無/非定常預(yù)處理結(jié)果中的壓力分布存在振蕩和較大誤差。圖16給出了流場的Q值等值面，本文Blend方法可以更好地進行旋翼尾跡的捕捉，而無/非定常預(yù)處理方法由于較大的耗散，顯著影響了尾跡捕捉的精度。

圖12 內(nèi)迭代數(shù)對機身法向力的影響Fig.12 Fuselage normal force history vs number of sub-iterations

圖13 通量格式對機身法向力系數(shù)的影響Fig.13 Fuselage normal force coefficient vs flux schemes

綜上，本文建立的Blend低速預(yù)處理方法保持了定常低速預(yù)處理的預(yù)測精度，并大幅改善了收斂性，總體計算效率可提高約5倍，相比定常預(yù)處理方法可提高內(nèi)迭代收斂標準，使計算精度得到保證。

圖14 直升機前飛模擬時內(nèi)迭代殘差收斂曲線Fig.14 Residual value cnvergence history in simulation of helicopter in forward flight

圖15 直升機前飛壓力系數(shù)分布對比Fig.15 Comparison of pressure coefficients distribution of helicopter in forward flight

圖16 直升機前飛流場的Q等值面Fig.16 Q criterion iso-surface for helicopter in forward flight

4.3 直升機/艦船干擾流場

本算例進行直升機/艦船干擾非定常流場的數(shù)值模擬，對本文Blend預(yù)處理方法收斂性和正確性進行考核。采用美國兩棲攻擊艦LHA的簡化模型作為艦船模型，將4.2節(jié)直升機模型放大到18 m作為模擬的旋翼和機身模型，如圖17所示。直升機位于艦首的第一起降點，槳盤中心距甲板13 m。艦船前進速度為20 m/s。旋翼的周期變距為θ0=10.3°，θ1c=-2.7°，θ1s=2.4°。

圖18為本算例所采用的重疊網(wǎng)格系統(tǒng)，共包含6塊子網(wǎng)格，分別為靜止的艦船網(wǎng)格和直升機機身的網(wǎng)格，網(wǎng)格單元數(shù)分別為1 140萬和677萬。繞每片旋翼葉片的網(wǎng)格隨旋翼運動，網(wǎng)格單元數(shù)為107萬。艦船網(wǎng)格采用非結(jié)構(gòu)混合網(wǎng)格，其余為六面體網(wǎng)格。重疊網(wǎng)格系統(tǒng)總網(wǎng)格數(shù)約2 440萬 。物面第1層網(wǎng)格y+≤1。圖19給出了直升機和艦船的表面網(wǎng)格。湍流模型采用Menter的SST兩方程模型。時間步長取旋翼轉(zhuǎn)過0.5° 方位角，參考長度為艦船甲板長度L=200 m，則非定常截斷參數(shù)Mau=L/(πΔtc)≈489.9。此時非定常預(yù)處理自動退化為無預(yù)處理的形式。

圖17 直升機/艦船干擾流場模擬模型Fig.17 Simulation model for helicopter/ship interaction flowfield

圖18 直升機/艦船干擾流場模擬重疊網(wǎng)格系統(tǒng)Fig.18 Overset grids system for helicopter/ship interaction flowfield simulation

圖19 直升機/艦船干擾流場表面網(wǎng)格Fig.19 Surface grids for helicopter/ship interaction flowfield

圖20 直升機/艦船干擾模擬流程Fig.20 Flowchart for simulation of helicopter/ship interaction

圖21 直升機/艦船干擾流場Q等值面Fig.21 Q criterion iso-surface for helicopter/ship interaction flowfield

圖22 直升機/艦船干擾流場垂向速度Fig.22 Vertical velocity of helicopter/ship interaction flowfield

圖23 直升機/艦船干擾模擬時殘差收斂曲線Fig.23 Residual value convergence history for simulation of helicopter/ship interaction

5 結(jié) 論

本文針對定常和非定常低速預(yù)處理在旋翼干擾流場數(shù)值模擬時遇到的兩難困境，在Potsdam和Sachdev的Blend低速預(yù)處理方法基礎(chǔ)上結(jié)合對流項和人工黏性項量級分析，將非定常低速預(yù)處理和定常低速預(yù)處理特征值矩陣進行組合，建立出一種新型Blend低速預(yù)處理方法。

1) Blend低速預(yù)處理方法通過非定常低速預(yù)處理矩陣和定常低速預(yù)處理矩陣進行組合，解決了在旋翼干擾流場數(shù)值模擬時非定常低速預(yù)處理自動退化為無預(yù)處理的問題。在雙時間步長內(nèi)迭代中相比定常低速預(yù)處理收斂效率大幅提高，可提高內(nèi)迭代收斂標準，使計算精度得到保證。

2) 無/非定常低速預(yù)處理在處理此類流場時，存在耗散過大的問題。充分收斂后定常低速預(yù)處理和Blend低速預(yù)處理計算結(jié)果表現(xiàn)出較好的一致性。

3) 本文建立的Blend低速預(yù)處理方法相比于Potsdam和Sachdev的方法，對流項和人工黏性項量級更匹配，修正了溫度場黏性過小的問題，提高了收斂效率和魯棒性。