王光學(xué)，鄧小剛，劉化勇，王運(yùn)濤

（中國空氣動力研究與發(fā)展中心 空氣動力學(xué)國家重點(diǎn)實(shí)驗(yàn)室，四川 綿陽621000）

0 引 言

高階精度格式作為數(shù)值計算的重要工具一直在諸多學(xué)科領(lǐng)域占有顯要地位，譬如，湍流的直接數(shù)值模擬（DNS）和大渦模擬 （LES）、計算氣動聲學(xué)（CAA）、計算電磁學(xué)（CEM）和磁流體力學(xué)（MHD）等。通過高階精度數(shù)值計算，一方面可以觀察流動細(xì)節(jié)、為機(jī)理分析提供精細(xì)數(shù)據(jù)；另一方面可以彌補(bǔ)實(shí)驗(yàn)條件限制和理論分析的不足，擴(kuò)大研究范圍。特別是，隨著人們對計算精度要求的不斷提高，在工程應(yīng)用領(lǐng)域和研究領(lǐng)域都越來越需要更為細(xì)致準(zhǔn)確的數(shù)值模擬，這為高階精度格式的研究和應(yīng)用創(chuàng)造了新的契機(jī)。

WCNS是鄧小剛自20世紀(jì)90年代提出的高階精度格式。首先，他們通過在中心型緊致格式中加入耗散項，于1996年構(gòu)造了單參數(shù)線性耗散緊致格式（Dissipative Compact Schemes，DCS）［1］，通過調(diào)整單一參數(shù)，精度可分別達(dá)到3至9階。同時，他們又提出了自適應(yīng)插值的概念，構(gòu)造了一類高階緊致非線性格式（Compact Nonlinear Schemes，CNS）［2］，這樣可以解決DCS因?yàn)榫€性格式很難用于含有激波流場的困難。然后通過引入加權(quán)思想，構(gòu)造了一系列加權(quán)緊致非線性格式 WCNS（Weighted Compact Nonlinear Scheme）［3－5］，包括隱式和顯式兩種類型。

50年代以來，對三角翼大攻角繞流的研究一直是流體力學(xué)工作者們所特別關(guān)注的一個重要方向。在這個典型的流場環(huán)境中，包含了非常豐富和復(fù)雜的流體力學(xué)現(xiàn)象。這些現(xiàn)象的研究與探討對空氣動力學(xué)理論的發(fā)展與完善具有重大意義，對提高現(xiàn)代飛行器性能具有重要價值。隨著計算機(jī)技術(shù)和高精度、非線性計算方法的迅速發(fā)展，開展了廣泛的數(shù)值模擬和理論研究，但由于該現(xiàn)象的高度復(fù)雜性，人們對其流動機(jī)理尚未充分掌握，還不能形成一個完整合理的理論以解釋渦破裂現(xiàn)象的產(chǎn)生和發(fā)展過程。

本文通過求解任意坐標(biāo)系下的雷諾平均的N－S方程，采用5階精度的加權(quán)緊致非線性格式（WCNSE－5）和多塊對接結(jié)構(gòu)網(wǎng)格技術(shù)，在與相應(yīng)試驗(yàn)結(jié)果對比的基礎(chǔ)上，詳細(xì)研究了 WCNS－E－5格式以及兩種湍流模型對主渦二次渦相互作用、渦破裂位置和表面壓力分布的影響。本文的研究結(jié)果表明，高階精度格式WCNS能成功應(yīng)用于三角翼的跨聲速大攻角流動，網(wǎng)格規(guī)模的增加進(jìn)一步提高流場分辨率，SST湍流模型相對SA湍流模型在三角翼大攻角流動中具有更好的適用性。

1 計算構(gòu)型和計算方法

本文采用的計算模型是65°后掠三角翼，它是1996年NASA Langley為研究該外形的雷諾數(shù)、馬赫數(shù)影響而完成的試驗(yàn)?zāi)Ｐ?，包括可替換的四個前緣外形，本文選用的是尖前緣，其前緣半徑為0。風(fēng)洞試驗(yàn)主要是在NASA Langley的NTF跨聲速風(fēng)洞中完成的，包括測力和測壓試驗(yàn)。后來該外形成為國際旋渦試驗(yàn)研究的標(biāo)準(zhǔn)模型，在德國的DLR等研究機(jī)構(gòu)完成了PSP、PIV試驗(yàn)。計算構(gòu)型如圖1所示，模型詳細(xì)介紹見文獻(xiàn)［6］。

圖1 65°后掠三角翼構(gòu)型Fig.1 Configuration of the 65°delta wing

本文模擬的狀態(tài)為：與測力測壓試驗(yàn)［6］對應(yīng)的狀態(tài)為馬赫數(shù) M＝0.85，攻角α＝－1°～27°，雷諾數(shù)Re＝6×106；與PSP［6］、PIV［8］試驗(yàn)對應(yīng)的狀態(tài)為馬赫數(shù)M＝0.8，攻角α＝25.5°，雷諾數(shù)Re＝2×106。力矩參考點(diǎn)距頭部頂點(diǎn)三分之二根弦長。Cp分布為沿流向5個截面，其中X／Cr＝0.2、0.4、0.6、0.8、0.95。

本文采用的計算方法為：對流項離散格式采用5階精度加權(quán)緊致非線性格式（WCNS－E－5），粘性項離散格式為4階精度中心格式，湍流模型采用了SA一方程和SST兩方程模型，計算中采用了并行技術(shù)加速。

其中5階精度加權(quán)緊致非線性格式（WCNS－E－5）對流項離散格式如下：（a）內(nèi)點(diǎn)格式

對流項離散采用原始變量型的 WCNS－E－5格式，設(shè)網(wǎng)格間距為h，以ξ方向?yàn)槔?/p>

（b）邊界格式

邊界和靠近邊界格式如下

2 網(wǎng)格影響分析

本文采用的計算網(wǎng)格均自主生成，網(wǎng)格結(jié)構(gòu)為多塊對接網(wǎng)格（1－to－1），網(wǎng)格分為粗網(wǎng)格、中等網(wǎng)格和細(xì)網(wǎng)格三種，網(wǎng)格的詳細(xì)信息如表1所示，圖2給出了粗網(wǎng)格的表面網(wǎng)格和對稱面網(wǎng)格。由表中可以看出，在該網(wǎng)格序列中，各套網(wǎng)格之間并沒有2倍數(shù)關(guān)系，這也許會給網(wǎng)格收斂性研究帶來一定的影響。

圖2 三角翼的計算網(wǎng)格Fig.2 Computational grid of the 65°delta wing

表1 網(wǎng)格統(tǒng)計列表Table 1 Computational grid size for the case

三套網(wǎng)格的網(wǎng)格拓?fù)潢P(guān)系相同，以“C”型網(wǎng)格為主導(dǎo)，計算區(qū)域的遠(yuǎn)場邊界取為25倍弦長。壁面的第一排網(wǎng)格達(dá)到了1.0×10－6弦長，網(wǎng)格在背風(fēng)區(qū)和各個剪切層附近均進(jìn)行了適當(dāng)?shù)募用埽员ＷC分離區(qū)、附面層內(nèi)和剪切層的數(shù)值模擬精度。

首先回顧尖前緣后掠翼的前緣分離的基本特性，由于是尖前緣，在較小的攻角主渦分離就開始于翼尖，從下至上流動卷起形成主渦分離，攻角增大，主渦直徑和強(qiáng)度均增加，位置向機(jī)翼內(nèi)側(cè)移動。主渦旋轉(zhuǎn)產(chǎn)生展向流動，在旋渦下部、三角翼的上翼面加速，出現(xiàn)負(fù)壓力峰值，該點(diǎn)沿展向到翼邊緣是逆壓梯度，將誘導(dǎo)邊界層分離，產(chǎn)生二次分離和二次旋渦，在主渦負(fù)壓力峰值外側(cè)又出現(xiàn)二次負(fù)壓力峰值。當(dāng)攻角大到一定程度，主渦開始破裂，當(dāng)破裂位置發(fā)展到機(jī)翼表面，將引起突然失速。

圖3是馬赫數(shù)M＝0.85，攻角α＝20.6°，雷諾數(shù)Re＝6×106狀態(tài)下，不同計算網(wǎng)格的不同流向位置橫截面壓力系數(shù)的比較。

在X／Cr＝0.2截面，中、細(xì)網(wǎng)格的上翼面出現(xiàn)了兩次負(fù)壓力系數(shù)峰值，表明在此截面，三角翼上翼面流動出現(xiàn)了主渦和二次渦，而粗網(wǎng)格的二次渦不明顯出現(xiàn)。隨網(wǎng)格加密，計算的壓力系數(shù)曲線與試驗(yàn)結(jié)果更為接近，并且更明顯分辨出主渦和二次渦。從X／Cr＝0.4、0.6、0.8截面看出，中、細(xì)網(wǎng)格的計算結(jié)果比較接近。

圖4是粗、細(xì)網(wǎng)格的X／Cr＝0.2截面的壓力云圖和速度矢量圖。粗網(wǎng)格分辨率較低，其二次渦不明顯，而細(xì)密度網(wǎng)格對流場的描述更細(xì)致，計算結(jié)果與試驗(yàn)結(jié)果更為接近，本文下面的計算都采用細(xì)網(wǎng)格模擬。

圖3 不同計算網(wǎng)格的截面壓力系數(shù)的比較Fig.3 Comparison of surface pressures coefficient forα＝20.6°

圖4 X／Cr＝0.2截面的壓力云圖和速度矢量圖Fig.4 Comparison of pressures and velocity distribution at X／Cr＝0.2

雖然隨網(wǎng)格規(guī)模的增加，提高了流場分辨率，中、細(xì)網(wǎng)格的計算結(jié)果也比較接近了，但依然不能稱其為達(dá)到網(wǎng)格收斂性。在高階精度模擬條件下，網(wǎng)格收斂性研究對網(wǎng)格序列的要求以及判斷依據(jù)有待進(jìn)一步研究。

3 湍流模型影響分析

在流場求解的RANS方法中，NS方程求解采用了高階精度，約束計算精準(zhǔn)度的瓶頸之一便是湍流模型了，而不同湍流模型在不同流場特性中又表現(xiàn)不同，為了考核高階精度求解條件下，湍流模型在三角翼大攻角旋渦流動中的表現(xiàn)，本節(jié)選取了一方程SA湍流模型和兩方程SST湍流模型進(jìn)行對比計算。

圖5是本文細(xì)網(wǎng)格計算結(jié)果與文獻(xiàn)［6］試驗(yàn)結(jié)果的比較。在攻角20°之前，兩種湍流模型的法向力系數(shù)計算結(jié)果與試驗(yàn)結(jié)果都相近；在大攻角范圍，SST湍流模型的法向力系數(shù)隨攻角變化的拐折點(diǎn)與試驗(yàn)結(jié)果一致，都在攻角23.6°和24.6°，其中攻角24.6°時主渦破裂，SA湍流模型的計算結(jié)果在本文計算的攻角范圍內(nèi)沒出現(xiàn)拐折，表明主渦沒有破裂。俯仰力矩系數(shù)的計算結(jié)果與實(shí)驗(yàn)結(jié)果還有一定差距，特別是俯仰力矩系數(shù)曲線斜率和在大攻角時的曲線拐折點(diǎn)的模擬差距較大，這是數(shù)值模擬的難點(diǎn)也是進(jìn)一步研究的重點(diǎn)。

圖6是兩種湍流模型在攻角20.6°時三角翼上表面壓力系數(shù)等值線云圖的比較。右圖是采用SST湍流模型計算的結(jié)果，在三角翼的一側(cè)可以仔細(xì)分辨出兩道綠色條紋，從頂部向尾部延伸，說明出現(xiàn)了主渦、二次渦分離；左圖是采用SA湍流模型計算的結(jié)果，其SA綠色條紋僅一道，表明只有主渦分離，二次渦不明顯。這從圖7兩種湍流模型的截面壓力系數(shù)的比較也可看出，SST湍流模型的計算結(jié)果與試驗(yàn)更為接近，沿展向出現(xiàn)了兩次壓力峰值，表明出現(xiàn)了主渦和二次渦，其位置也和試驗(yàn)一致，而SA湍流模型的計算結(jié)果沒有出現(xiàn)兩次壓力峰值。

圖5 氣動力系數(shù)的比較Fig.5 Comparison of aerodynamics coefficient

圖6 表面壓力系數(shù)的比較（SA、SST湍流模型）Fig.6 Comparison of surface pressures distributions

圖7 截面壓力系數(shù)的比較（SA、SST湍流模型）Fig.7 Comparison of pressures coefficient

通過以上分析說明，在本文采用的計算方法中，SST湍流模型在模擬大攻角大分離流動時比SA湍流模型更為適用。

4 渦破裂位置的計算結(jié)果與分析

本節(jié)采用細(xì)網(wǎng)格和SST湍流模型精細(xì)模擬了馬赫數(shù)0.8、攻角25°，雷諾數(shù)200萬狀態(tài)下的三角翼繞流。對于這個狀態(tài)，渦破裂位置是數(shù)值模擬的難點(diǎn)和重點(diǎn)，并且在三角翼上表面出現(xiàn)了超音速區(qū)域，以及橫向流動產(chǎn)生較強(qiáng)激波，這都對高階精度格式數(shù)值模擬提出了挑戰(zhàn)。

圖8是該狀態(tài)下的表面壓力系數(shù)的比較，其中左圖為文獻(xiàn)［7］的PSP試驗(yàn)結(jié)果，左側(cè)表面不是太光滑，出現(xiàn)左右不對稱。右圖為本文計算結(jié)果。其計算結(jié)果與實(shí)驗(yàn)結(jié)果右側(cè)比較，主渦、二次渦的位置和主渦破裂位置都吻合較好。圖9、10是X／Cr＝0.6、0.7截面的速度矢量圖，其中左圖為文獻(xiàn)［8］的PIV試驗(yàn)結(jié)果，右圖為本文計算結(jié)果。圖11主渦破裂位置示意圖，以及X／Cr＝0.2、0.4、0.6、0.8截面的壓力等值線圖。從這些圖中看出，在該計算狀態(tài)下，主渦、二次渦從三角翼頂端分離出，主渦渦核位置較高，其旋轉(zhuǎn)產(chǎn)生側(cè)向流動，在渦核下面沿展向是逆壓梯度，與附面層作用，導(dǎo)致附面層分離，產(chǎn)生二次渦；在X／Cr＝0.6截面，主渦還沒出現(xiàn)破裂，從截面速度矢量圖中看出主渦逆時針旋轉(zhuǎn)，渦核清晰；在X／Cr＝0.7截面，主渦破裂，在渦核位置出現(xiàn)反向旋轉(zhuǎn)流動，從渦破裂示意圖看出，在X／Cr＝0.7截面附近沿渦核流線不再緊密聚集，流線從渦核起向外旋轉(zhuǎn)散開，渦結(jié)構(gòu)已經(jīng)變得很松散，而二次渦并沒有因?yàn)橹鳒u破裂而消失或破裂。

同時，從圖12、13的X／Cr＝0.6、0.7截面的側(cè)向速度和馬赫數(shù)云圖中看出，由于本文選取的計算狀態(tài)為跨聲速流動，其三角翼上表面出現(xiàn)了超音速區(qū)域，在主渦下部靠近壁面位置，激波誘導(dǎo)附面層分離，產(chǎn)生二次分離。對于求解含有激波的跨／超音速流動時，對計算格式的要求更加苛刻，不僅要求格式能以極小的數(shù)值耗散來體現(xiàn)小尺度的流場結(jié)構(gòu)，還希望它具有光滑捕捉流場間斷的性能。而本文采用的高階精度格式（WCNS）在這方面具有一定優(yōu)勢，在減小數(shù)值耗散的同時，能光滑捕捉激波，較好地解決了這兩個幾乎矛盾的雙重目標(biāo)。

5 結(jié) 論

本文采用5階精度的加權(quán)緊致非線性格式（WCNS－E－5）和多塊對接結(jié)構(gòu)網(wǎng)格技術(shù)，通過求解任意坐標(biāo)系下的RANS方程，數(shù)值模擬了65°后掠三角翼的復(fù)雜流場，主要研究了 WCNS－E－5格式以及兩種湍流模型對主渦二次渦相互作用、渦破裂位置和表面壓力分布的影響。通過與相應(yīng)的試驗(yàn)結(jié)果相比較，得到以下一些基本結(jié)論：

（1）高階精度格式（WCNS－E－5）能成功應(yīng)用于三角翼的跨聲速大攻角流動。

（2）SST湍流模型相對SA湍流模型在三角翼大攻角流動中具有更好的適用性。

（3）隨網(wǎng)格規(guī)模的增加，提高了數(shù)值模擬精度和流場分辨率，但高階精度數(shù)值模擬條件下的網(wǎng)格收斂性有待進(jìn)一步研究。

（4）本文采用的數(shù)值方法能較準(zhǔn)確地模擬三角翼主渦破裂位置，其結(jié)果與PSP、PIV實(shí)驗(yàn)結(jié)果一致。

［1］DENG X G，MAEKAWA H，SHEN C.A class of high order dissipative compact schemes［R］.AIAA Paper 96－1972.

［2］DENG X G，MAEKAWA H.Compact high－order accurate nonlinear schemes［J］.J.Comput.Phys.，1997，130：77－91.

［3］DENG X G，MAO M L.Weighted compact high－order nonlinear schemes for the Euler equations［R］.AIAA Paper 97－1941.

［4］DENG X G，ZHANG H X.Developing high－order weighted compact nonlinear schemes［J］.J.Comput.Phys.，2000，165：24－44.

［5］DENG X.High－order accurate dissipative weighted compact nonlinear schemes［J］.ScienceinChina（SeriesA），2002，45（3）：356－370.

［6］CHU J，LUCKRING J M.Experimental surface pressure data obtained on 65°delta wing across Reynolds number and mach number ranges［R］.NASA TM 4645，1996.

［7］KONRATH R，KLEIN C，ENGLER R H，et al.Analysis of PSP results obtained for the VFE－2 65°delta wing configuration at subsonic and transonic speed［R］.AIAA Paper 2006－0060，2006.

［8］SCHR?DER A，AGOCS J，F(xiàn)RAHNERT H，et al.Application of stereo－PIV to the VFE－2 65°delta wing configuration at sub－and transonic speeds［R］.AIAA Paper 2006－3486，2006.