65°三角翼亞音速復(fù)雜流場計算和數(shù)據(jù)可視化1)

李 立

(中國航空工業(yè)集團(tuán)公司西安航空計算技術(shù)研究所，西安710065)

應(yīng)用研究

65°三角翼亞音速復(fù)雜流場計算和數(shù)據(jù)可視化1)

李 立2)

(中國航空工業(yè)集團(tuán)公司西安航空計算技術(shù)研究所，西安710065)

提出一種基于非結(jié)構(gòu)混合網(wǎng)格和有限體積法的有效計算策略，對第二期國際渦流試驗項目(second international vortex flow experiment,VFE-2)的尖前緣65°三角翼在馬赫數(shù)0.4，迎角20.3°，雷諾數(shù)2×106條件下的亞音速復(fù)雜流場結(jié)構(gòu)進(jìn)行數(shù)值模擬，重點探討了基于計算數(shù)據(jù)進(jìn)行該類型復(fù)雜渦系干擾表面和空間流場關(guān)鍵特征提取和數(shù)據(jù)可視化問題.通過與相關(guān)試驗類比，建立了與先進(jìn)試驗流動顯示技術(shù)相比擬的定性和定量分析方法，為三角翼這類復(fù)雜流場結(jié)構(gòu)的精細(xì)分析奠定了技術(shù)基礎(chǔ).采用上述方法，細(xì)致分析了亞音速三角翼的大迎角復(fù)雜旋渦流場結(jié)構(gòu)，得到了與試驗一致的結(jié)論.研究證實：在大迎角條件下，三角翼流動物理復(fù)雜，黏性效應(yīng)耦合嚴(yán)重，只有通過N-S方程計算才能準(zhǔn)確地捕捉主渦和二次渦的發(fā)展.

數(shù)據(jù)可視化，混合網(wǎng)格，有限體積法，旋渦運動，VFE-2(second international vortex flow experiment)三角翼

Key wordsflow visualization,hybrid mesh,f i nite volume method,vortex flow,VFE-2(second international vortex flow experiment)delta wing

三角翼是現(xiàn)代戰(zhàn)機(jī)及無人機(jī)的常見布局形式，具有深刻的工程應(yīng)用背景.在亞音速條件下，繞三角翼的流動是典型的幾何簡單，但物理形成機(jī)制復(fù)雜的旋渦流動.在不大的迎角下，三角翼將在上翼面形成前緣分離渦，并隨著迎角增加不斷增強(qiáng)；當(dāng)迎角超過一定限度后，旋渦從穩(wěn)定發(fā)展為不穩(wěn)定，直至從后緣逐步發(fā)生渦破裂.渦破裂對飛機(jī)氣動特性影響非常明顯，嚴(yán)重情況下，會對飛機(jī)的操縱性和穩(wěn)定性帶來致命影響.為此，在過去數(shù)十年間，各國都設(shè)立了相關(guān)項目對三角翼的旋渦結(jié)構(gòu)及形成機(jī)理開展研究，其中尤以美國和歐盟聯(lián)合發(fā)起的第二期國際渦流試驗項目(second international vortex flow experiment,VFE-2)最為知名[1].VFE-2項目的重要價值在于通過發(fā)展光學(xué)壓敏測量技術(shù)(pressure-sensitive paint technique,PSP)、粒子影像測速儀技術(shù)(particle image velocimetry technique,PIV)等先進(jìn)試驗手段，測量得到VFE-2系列三角翼模型的表面流場和空間流場精細(xì)試驗數(shù)據(jù)[2]，為CFD軟件的對比分析和結(jié)果確認(rèn)提供了豐富的資源.

近年來，數(shù)值計算已在飛行器設(shè)計中得到廣泛應(yīng)用.采用數(shù)值方法進(jìn)行流場分析的一個巨大優(yōu)勢是能夠獲得比試驗多得多的流場信息.但選取哪些數(shù)據(jù)來進(jìn)行流場分析和進(jìn)行直觀的流動顯示，具有一定的挑戰(zhàn)性.由于不同的流動往往具有不同的流動特征，因而針對不同問題不大可能采用完全相同的流動分析和顯示方法來統(tǒng)一處理.最恰當(dāng)?shù)霓k法應(yīng)是針對不同流動的物體特征，研究建立適當(dāng)?shù)年P(guān)鍵流動特征提取方法.基于這一考慮，本文采用混合網(wǎng)格策略和有限體積法，數(shù)值求解了VFE-2尖前緣65°三角翼在馬赫數(shù)0.4，迎角20.3°，雷諾數(shù)2×106條件下的亞音速流場.以此為基礎(chǔ)，重點探討基于計算數(shù)據(jù)如何進(jìn)行該類型復(fù)雜渦系干擾表面和空間流場關(guān)鍵特征提取和數(shù)據(jù)可視化，獲得與先進(jìn)試驗手段(測壓、油流、PSP、PIV等)相比擬的數(shù)據(jù)可視化方法和結(jié)果，從而為三角翼這類復(fù)雜流場結(jié)構(gòu)的精細(xì)分析奠定技術(shù)基礎(chǔ).

1 數(shù)值方法

1.1 控制方程及求解

控制方程為雷諾平均 N-S(Reynolds-averaged Navier-Stokes,RANS)方程，守恒形式如下

其中

這里，ρ,U=(u,v,w)T,E,p,τij,qi分別表示密度、速度矢量、總內(nèi)能、壓強(qiáng)、黏性剪切應(yīng)力張量和熱流;δij表示Kronecker函數(shù),δii=1,δij=0(i/=j).剪切應(yīng)力張量可根據(jù)Boussinesq假設(shè)，表示為

式中，μ為黏性系數(shù)，為層流部分μl和湍流部分μT之和,μ，κ為湍動能,I為單位矩陣.壓強(qiáng)根據(jù)狀態(tài)方程計算，對理想氣體有

式中，γ為比熱比，對空氣，γ=1.4.

湍流模型方程采用Spalart-Allmaras(SA)一方程模型[3].相對其他復(fù)雜的湍流模型，SA一方程模型具有計算效率高、魯棒性好、對分離流動的模擬精度高等優(yōu)點，是航空工程領(lǐng)域應(yīng)用最廣泛的湍流模型之一.為了改善計算的實際收斂性能，本文還引入了渦量修正技術(shù)，即將模型方程的渦量值修正為

式中，fv1=χ3/[χ3+(0.71)3],χ=ρ?v/μl,?v為SA方程的求解變量.

對方程(1)，本文采用非結(jié)構(gòu)混合網(wǎng)格策略和有限體積法進(jìn)行數(shù)值求解.目前，利用非結(jié)構(gòu)網(wǎng)格進(jìn)行有限體積離散主要有兩種方式.一種直接采用網(wǎng)格單元作為控制體，稱為格心(cell-centered)格式；一種以網(wǎng)格單元頂點為中心，采用對偶網(wǎng)格方法建立控制體，稱為格點(cell-vertex)格式.兩種方式各有優(yōu)劣，一般說來，格點格式效率比格心格式要高[4].本文選用格點格式.如圖1所示，對每個網(wǎng)格單元的節(jié)點，虛框包圍部分是實際控制體.由式(1)，在每個控制體上求體積積分，并運用Gauss公式，可得

具體計算采用Jameson等[5]提出的標(biāo)準(zhǔn)中心差分格式進(jìn)行空間離散，并采用3階總變差Runge-Kutta方法進(jìn)行時間推進(jìn)求解.計算過程中，使用了當(dāng)?shù)貢r間步長、隱式殘值光順及聚合多重網(wǎng)格技術(shù)等進(jìn)行計算加速[6].

圖1 格點格式有限體積法的控制體構(gòu)造

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

計算模型是VFE-2尖前緣65°三角翼，模型外形和幾何定義參數(shù)定義見參考文獻(xiàn)[2].VFE-2項目2004年由美國和歐盟聯(lián)合發(fā)起，采用不同前緣鈍度的三角翼模型作為研究對象，旨在為CFD軟件對復(fù)雜渦流場的預(yù)測能力評估提供可靠的試驗數(shù)據(jù).對該系列三角翼模型，NASA在1996年已先期開展了大量試驗研究，提供了在較大范圍雷諾數(shù)和可變馬赫數(shù)狀態(tài)下的表面壓力分布、法向力、俯仰力矩系數(shù)等試驗數(shù)據(jù)[7].與NASA早期的試驗不同，VFE-2項目偏重于進(jìn)行精細(xì)風(fēng)洞試驗.試驗在德國的德-荷風(fēng)洞機(jī)構(gòu)的跨音速風(fēng)洞中進(jìn)行，采用PSP和PIV等更先進(jìn)試驗手段，提供了在不同馬赫數(shù)、不同雷諾數(shù)狀態(tài)、不同迎角范圍內(nèi)的詳細(xì)、可靠的流場數(shù)據(jù)[2].本文選取M∞=0.4,α=20.3°,Re=2.×106的狀態(tài)開展研究.該狀態(tài)接近于渦破裂的臨界狀態(tài)，對軟件能力的要求較高[8].

為了節(jié)省計算量，全部計算基于半翼展模型(半模)展開.采用的計算網(wǎng)格自主生成，物面和空間網(wǎng)格均采用網(wǎng)格自適應(yīng)技術(shù)進(jìn)行了加密以提高邊界層、分離區(qū)及剪切層的模擬精度[9-11].如圖2所示，計算區(qū)域沿流向的遠(yuǎn)場邊界取為50倍弦長，沿展向的遠(yuǎn)場邊界取為20倍弦長.邊界層內(nèi)第1層網(wǎng)格達(dá)到5.0×10-6弦長.半模網(wǎng)格總的網(wǎng)格點數(shù)約為65萬，總的網(wǎng)格單元數(shù)約為170萬.計算中，遠(yuǎn)場均采用無反射特征遠(yuǎn)場邊界條件，物面均采用無滑移絕熱壁面邊界條件.

圖2 三角翼模型的計算網(wǎng)格

2 計算結(jié)果及數(shù)據(jù)可視化分析

2.1 表面流場的流動顯示

對于試驗而言，用于表面流場的試驗流動顯示技術(shù)主要包括油流試驗技術(shù)和PSP等.通過計算數(shù)據(jù)對油流試驗、PSP等試驗技術(shù)進(jìn)行復(fù)現(xiàn)，相對比較容易.

對油流試驗結(jié)果，計算可以通過繪制極限流線的方式來復(fù)現(xiàn).按照定義，極限流線是流線無限接近物面時的一種極限狀態(tài).為了通過計算數(shù)據(jù)進(jìn)行極限流線繪制，本文采取的具體步驟是：首先抽取物面網(wǎng)格、第1層邊界層網(wǎng)格及相應(yīng)的解數(shù)據(jù)(速度場)形成獨立的數(shù)據(jù)集；然后針對該數(shù)據(jù)集進(jìn)行流線繪制即得到極限流線.根據(jù)極限流線結(jié)果可以清晰判定流場的分離和再附.圖3給出采用本文RANS方法計算得到的極限流線結(jié)果.由圖可知，該三角翼構(gòu)型在選定狀態(tài)下發(fā)生明顯二次分離，主渦和二次渦附著線清晰可見.

圖3 計算對油流結(jié)果的復(fù)現(xiàn)

PSP是試驗中實時觀測三維流動現(xiàn)象的有用工具，試驗觀測結(jié)果為表面壓力系數(shù)分布.因此，計算可以采用繪制壓力系數(shù)分布云圖的方式來復(fù)現(xiàn).此外，由于計算不像試驗受物理條件限制，根據(jù)需要還可以顯示其他變量的分布云圖，其效果等同于PSP.對傳統(tǒng)PSP顯示進(jìn)行復(fù)現(xiàn)的方法過程為：首先，根據(jù)式(3)，可以由守恒量直接計算壓強(qiáng)；然后，根據(jù)壓力系數(shù)定義可得

這里，q∞為自由流動壓，p∞為自由流壓強(qiáng),M∞為來流馬赫數(shù).根據(jù)表面壓力系數(shù)分布可以清晰判定三角翼前緣渦在物面留下的印跡，并能定性和定量分析出吸力峰的位置.圖4給出本文計算得到的PSP結(jié)果與試驗結(jié)果[12]的對比.作為比較，同時給出Euler和RANS計算的不同結(jié)果.圖中，低壓區(qū)即是前緣渦在物面留下的印跡.在圖4的結(jié)果中，計算比試驗給出的低壓區(qū)更清晰，主要是由于對計算結(jié)果進(jìn)行數(shù)據(jù)可視化時對壓力系數(shù)分布區(qū)間作了特殊處理.可以看出，Euler計算給出的主渦低壓區(qū)從機(jī)翼頂點一直拓展到約80%弦長位置，之后發(fā)生壓力陡增.物面壓力分布的陡變實際上反映出渦的破裂.與Euler結(jié)果相比，RANS計算與試驗給出的結(jié)果更加接近，其主渦低壓區(qū)一直拓展到尾緣附近，沒有發(fā)生壓力陡增，表明在選定亞音速條件下，主渦沒有發(fā)生渦破裂.推測造成Euler計算與RANS計算結(jié)果差異大的主要原因是，三角翼邊界層黏性效應(yīng)干擾嚴(yán)重，Euler計算由于忽略了物理黏性影響，造成渦遠(yuǎn)離壁面后很快被耗散掉，因而提前發(fā)生渦破裂.從圖4的計算結(jié)果還可看出，Euler和RANS計算存在一個顯著的區(qū)別是后者能夠準(zhǔn)確預(yù)測出前緣和主渦之間的二次渦.圖4中試驗結(jié)果和RANS計算結(jié)果給出的前緣與主渦之間存在模糊的低壓痕跡直接反映的就是二次渦.這再次說明，對三角翼這類復(fù)雜渦流場，RANS計算結(jié)果明顯優(yōu)于Euler計算.與Euler計算結(jié)果相比，RANS計算給出的主渦和壓力吸力峰位置明顯更靠近翼根，與試驗結(jié)果相符.

圖4 計算對PSP結(jié)果的復(fù)現(xiàn)

根據(jù)PSP結(jié)果，可以進(jìn)一步進(jìn)行定量比較.圖5和圖6分別給出沿主渦渦核及沿流向不同橫向截面的壓力系數(shù)分布比較.圖5的結(jié)果進(jìn)一步證實了對圖4的分析，Euler計算會發(fā)生明顯的壓力陡增.圖6的結(jié)果表明本文RANS結(jié)果與試驗結(jié)果、文獻(xiàn)計算結(jié)果[10]符合較好，沿流向60%弦長之后的壓力吸力峰位置、強(qiáng)度結(jié)果均優(yōu)于文獻(xiàn).

圖5 沿主渦渦核的壓力系數(shù)分布比較(c為平均氣動弦長)

圖6 橫向截面壓力系數(shù)分布比較

2.2 空間流場的流動顯示

試驗中，用于空間流場流動顯示的主要技術(shù)是PIV技術(shù).PIV技術(shù)從20世紀(jì)80年代開始發(fā)展，最早用于應(yīng)力測量，但由于它能在不干擾流場的情況下，獲得整個瞬時以及時均的速度場，并可以進(jìn)一步得到渦流場等參數(shù)，很快在流場測量中得到廣泛應(yīng)用.在流場測量中，PIV給出的結(jié)果主要是不同截面位置的速度矢量場.圖7給出本文RANS計算在x/c=0.7站位對PIV測量速度場[13]的復(fù)現(xiàn)結(jié)果.可以看到，計算得到速度場所反映的關(guān)鍵流動特征與試驗基本一致，在機(jī)翼上方存在明顯的高速流動分離區(qū)，而在翼尖呈現(xiàn)明顯的低速回流區(qū).圖8進(jìn)一步給出沿流向不同站位的典型速度流場，清晰反映出三角翼上方主渦分離區(qū)沿流向方向的發(fā)展.與壓力云圖反映出的流場特征一樣，該三角翼在選定狀態(tài)沿流向方向直至80%弦長位置仍沒有發(fā)生渦破裂.

圖8 三角翼沿流向方向不同站位的速度流場

值得指出，與試驗相比，由于計算結(jié)果包含了更完整的流場信息(密度、壓強(qiáng)、速度等)，因此通過這些信息很容易得到用戶所關(guān)心的、更豐富的物理場.對于亞音速三角翼這類復(fù)雜渦流場，主要關(guān)心兩個問題：(1)渦核的計算；(2)設(shè)計何種物理量作為關(guān)鍵特征量以進(jìn)行更清晰、更直觀的空間流場展示.問題的實質(zhì)是，需要進(jìn)行空間流場關(guān)鍵特征的提取.借鑒課題組在復(fù)雜渦流場網(wǎng)格自適應(yīng)探測器設(shè)計方面的經(jīng)驗[9]，對三角翼渦流場關(guān)鍵特征的提取，本文提出以下具體思路.

渦核計算可采用對速度梯度張量進(jìn)行特征值分析的方法.步驟如下：(1)對網(wǎng)格單元的每個節(jié)點進(jìn)行速度梯度張量的計算Aij=(?Ui/?xj),i,j= 1,2,3；(2)進(jìn)行速度梯度張量的特征值分析，得到相應(yīng)的特征值及特征向量(λi,xi),i=1,2,3；(3)判斷該組特征值是否由一個實特征值及一對共軛的復(fù)特征值組成；如是，則把該節(jié)點標(biāo)記為渦核.

關(guān)于三角翼渦流場關(guān)鍵特征變量的定義問題，很顯然，選擇并不唯一.傳統(tǒng)方法中常把渦量作為反映渦強(qiáng)度及發(fā)展的物理量.通過對渦流場特點的分析，結(jié)合反復(fù)實踐和對比研究，發(fā)現(xiàn)，總壓比、熵估計及湍流渦黏性這幾種物理量均是能較好反映渦發(fā)展的空間流場關(guān)鍵特征量.總壓比和熵估計的定義式分別為

式中，M為當(dāng)?shù)伛R赫數(shù).

圖9給出對三角翼流場進(jìn)行特征提取和流動顯示的典型結(jié)果.采用在截面進(jìn)行空間流線投影的方式對渦的發(fā)展進(jìn)行展示，并采用湍流渦黏性、總壓比作為關(guān)鍵特征物理量進(jìn)行空間流動顯示，清晰顯示出空間上翼尖渦的發(fā)展.空間流場結(jié)果顯示，本文計算準(zhǔn)確捕捉到空間上二次渦的發(fā)展.從這一流動顯示結(jié)果可反映出，相對于試驗，計算在數(shù)據(jù)完備性方面具有天然優(yōu)勢.

圖9 三角翼空間流場關(guān)鍵特征提取和空間流場顯示

3 結(jié) 論

本文提出采用非結(jié)構(gòu)混合網(wǎng)格策略和有限體積法對65°VFE-2尖前緣三角翼進(jìn)行數(shù)值計算.結(jié)合與試驗流動顯示技術(shù)對比，開展了該類型復(fù)雜流動的數(shù)值流動顯示方法研究，為三角翼這類復(fù)雜流場結(jié)構(gòu)的精細(xì)分析奠定了技術(shù)基礎(chǔ).通過本文研究，形成主要結(jié)論如下：

(1)與試驗相比，數(shù)值計算在數(shù)據(jù)完備性方面具有天然優(yōu)勢，能夠完美復(fù)現(xiàn)典型試驗流動顯示技術(shù)(如油流、PSP、PIV等)，并能豐富傳統(tǒng)流動顯示技術(shù)的內(nèi)涵.

(2)本文計算準(zhǔn)確預(yù)測了尖前緣三角翼在亞音速大迎角條件下復(fù)雜的旋渦流場結(jié)構(gòu).計算表明，亞音速三角翼渦流場的黏性效應(yīng)嚴(yán)重，采用RANS計算能夠準(zhǔn)確捕捉到主渦、二次渦的發(fā)展，而Euler計算對渦的精細(xì)捕捉能力明顯不足，過度預(yù)測主渦的發(fā)展.

(3)流場關(guān)鍵特征的提取是流動顯示和數(shù)據(jù)可視化分析的關(guān)鍵問題.對三角翼復(fù)雜渦流場而言，除了傳統(tǒng)方法中常用的渦量外，實踐表明，總壓比和熵估計也均能很好地展示渦強(qiáng)度及渦的發(fā)展.

致謝：本文VFE-2項目試驗數(shù)據(jù)由德國宇航院授權(quán)使用，在此表示感謝.

1 Hummel D,Redeker G.A new vortex fl ow experiment for computer code validation.RTO AVT Symposium on Vortex Flow and High Angle of Attack Aerodynamics,Loen, Norway,2003

2 Konrath R,Klein C,Schr¨oder A.PSP and PIV investigations on the VFE-2 con fi guration in sub-and transonic fl ow.46th AIAA Aerospace Sciences Meeting and Exhibit, Reno,Nevada,2008

3 Spalart PR,Allmaras SR.A one-equation turbulence model for Aerodynamic fl ows.La Recherche Aerospatiale,1994, (1)：5-21

4 閻超,于劍,徐晶磊等.CFD模擬方法的發(fā)展成就與展望.力學(xué)進(jìn)展,2011,41(5)：562-588

5 Jameson A,Schmidt W,Turkel E.Numerical solutions of the Euler equations by fi nite volume methods using Runge-Kutta time stepping scheme.AIAA Paper 81-1259,1981

6 朱培燁.三維非結(jié)構(gòu)網(wǎng)格的歐拉方程聚合多重網(wǎng)格法.航空計算技術(shù),2004,34(3)：6-9

7 Chu J,Luckring JM.Experimental surface pressure data obtained on 65°delta wing across Reynolds number and Mach number ranges.NASA Technical Memorandum 96-4645,1996

8 Fritz W.What was learned from the numerical simulations for the VFE-2.AIAA Paper 2008-399,2008

9 Bai W,Qiu Z,Li L.Recent e ff orts to establish adaptive hybrid grid computing capability at ACTRI.Computational Fluid Dynamics Journal,2007,15(4)：438-449

10 Pirzadeh SZ.Vortical fl ow predication using an adaptive unstructured grid method.RAT AVT Symposium on Advanced Flow Management：Part A-Vortex Flows and High Angle of Attack for Military Vehicles,Loen,Norway,2001

11 李立,白文,梁益華.基于伴隨方程方法的非結(jié)構(gòu)網(wǎng)格自適應(yīng)技術(shù)及應(yīng)用.空氣動力學(xué)報,2011,29(3)：309-316

12 Konrath R,Klein C,Engler RH,et al.Analysis of PSP results obtained for the VFE-2 65°delta wing con fi guration at sub-and transonic speeds.44th AIAA Aerospace Sciences Meeting and Exhibit,Reno,Nevada,2006

13 Konrath R,Schr¨oder A,Kompenhans J.Analysis of PIV results obtained for the VFE-2 65°delta wing con fi guration at sub-and transonic speeds.24th Applied Aerodynamics Conference,San Francisco,California,2006

(責(zé)任編輯：劉希國)

NUMERICAL SIMULATION AND VISUALIZATION OF COMPLEX SUBSONIC FLOW OVER A 65°DELTA WING1)

LI Li2)
(AVIC Aeronautics Computing Technique Research Institute,Xi’an 710065,China)

The VFE-2 65°delta wing with sharp leading edge in subsonic flows at Mach number of 0.4,angle of attack of 20.3°,and Reynolds number of 2×106is numerically simulated by using an unstructured hybrid mesh based f i nite volume method,with emphasis on how to extract the key flow features for visualization both on surface and in space for such type of complex flows.Approaches for the advanced flow visualization techniques are used for qualitative and quantitative analyses,as a solid foundation for elaborated analysis of complex flow structures of delta wing.With these approaches,the complex vortex flow structure of the subsonic delta wing at high angle of attack is analyzed,and results are consistent with experiments.It is shown that the flow over delta wing at high angle of attack has a complex physical nature with a strong viscous coupling ef f ect,and the evolutions of the primary and secondary vortices can be accurately captured only by a Navier-Stokes equation based simulation.

V211.3

A

10.6052/1000-0879-16-195

2016-06-12收到第1稿，2016-07-07收到修改稿.

1)國家863計劃(2012AA01A304)和航空科學(xué)基金(2015ZA31002)資助項目.

2)李立，高級工程師，研究方向為計算流體力學(xué).E-mail：westlili@163.com

李立.65°三角翼亞音速復(fù)雜流場計算和數(shù)據(jù)可視化.力學(xué)與實踐,2017,39(1)：18-24

Li Li.Numerical simulation and visualization of complex subsonic flow over a 65°delta wing.Mechanics in Engineering,2017,39(1)：18-24