司芳芳，袁先旭，賀旭照，謝昱飛，葉友達(dá)

(1.北京流體動(dòng)力科學(xué)研究中心， 北京 100120；2.中國空氣動(dòng)力研究與發(fā)展中心， 四川 綿陽 621000)

1 引言

未來戰(zhàn)斗機(jī)需要大幅度提升空戰(zhàn)性能，目前先進(jìn)戰(zhàn)斗機(jī)如F-22和F-35均采用脊形前體，以降低可探測性和提升空戰(zhàn)性能。脊形前體主要是由切拱弧和脊形邊緣組成。攻角變化時(shí)，脊形邊緣使得背風(fēng)流場的分離點(diǎn)固定，與尖前緣高后掠三角翼的流場十分相似。脊形前體高度混合外形良好的隱身性能和超聲速飛行性能使其成為未來高性能戰(zhàn)斗機(jī)和巡航導(dǎo)彈比較理想的前體布局形式。

與具有光滑前體的傳統(tǒng)戰(zhàn)機(jī)相比，脊形前體通常會(huì)產(chǎn)生更強(qiáng)的前體渦流。某些特定飛行條件下，脊形前體產(chǎn)生的前體渦與機(jī)翼的前緣渦以及飛機(jī)其他結(jié)構(gòu)之間存在強(qiáng)干擾，所以國外研究人員對(duì)脊形前體、三角翼等結(jié)構(gòu)之間的渦干擾進(jìn)行了大量的試驗(yàn)研究。目前，脊形前體數(shù)值模擬方面的研究主要有Ravi和Mason等用歐拉方程和雷諾方程數(shù)值模擬研究了大攻角、有側(cè)滑情況下脊形前體的方向穩(wěn)定性。Agosta-Greenman等用薄邊界層雷諾平均納維-斯托克斯方程(navier-stokes equations,N-S方程)數(shù)值模擬了切向吹氣對(duì)脊形前體氣動(dòng)力的影響。在20世紀(jì)90年代中期，美國國家航空航天局蘭利研究中心、洛克希德·沃斯堡公司、洛克希德·馬丁公司和波音公司都致力于用歐拉數(shù)值模擬方法來研究中高攻角時(shí)脊形前體類飛機(jī)的氣動(dòng)力特性，一致的結(jié)論是歐拉解對(duì)飛機(jī)的初期設(shè)計(jì)有一定的指導(dǎo)意義，雖然歐拉方程可以準(zhǔn)確的模擬主渦結(jié)構(gòu)，但是不能研究由粘性產(chǎn)生的二次渦和三次渦以及分離渦的粘性干擾等非線性情況，所以只有采用湍流三維N-S方程來模擬非定常分離渦，以及前體渦與機(jī)翼渦之間的干擾等現(xiàn)象。Jeans等采用DDES(delayed detached-eddy simulation)湍流模型數(shù)值模擬研究了脊形前體類簡化飛機(jī)模型的大攻角非定常以及動(dòng)態(tài)滾轉(zhuǎn)現(xiàn)象，模擬出了試驗(yàn)中的非線性氣動(dòng)力。

盡管已進(jìn)行大量工作，但由于脊形前體大攻角非定常大范圍湍流分離流的復(fù)雜性，目前對(duì)脊形前體飛行器的大攻角流動(dòng)機(jī)理認(rèn)識(shí)尚不充分，深入了解其大攻角氣動(dòng)特性，尤其是脊形前體與翼型相互干擾的氣動(dòng)特性有著重要的工程應(yīng)用價(jià)值以及學(xué)術(shù)意義。

針對(duì)脊形前體類飛行器大攻角非定常湍流大分離流動(dòng)特點(diǎn)，本文采用IDDES(improved delayed detached-eddy simulation)方法來研究脊形角7.5°的前體添加平板機(jī)翼組成的脊形前體飛行器在不同攻角和雷諾數(shù)下的氣動(dòng)特性，尤其是大攻角氣動(dòng)特性和前體渦與機(jī)翼渦耦合后空間渦流場結(jié)構(gòu)演化的影響規(guī)律，為更好地進(jìn)行先進(jìn)高機(jī)動(dòng)飛行器設(shè)計(jì)提供參考。

2 數(shù)值模擬方法

2.1 湍流模型

目前國際上各種復(fù)雜的分離流動(dòng)模擬中常用的是DES(detached-eddy simulation)類方法。DES類方法的基本思路是用統(tǒng)一的湍流模型，以網(wǎng)格尺度和模型中的特征尺度隱式地劃分大渦模擬(LES)和雷諾平均Navier Stokes(RANS)區(qū)域。DES 方法最初是基于SA(spalart allmaras)模型建立起來的，用LES方法模擬大尺度運(yùn)動(dòng)占主導(dǎo)地位的非定常分離湍流流動(dòng)區(qū)域，用RANS方法模擬小尺度運(yùn)動(dòng)占主導(dǎo)地位的近壁區(qū)域。隨著DES方法研究的深入，當(dāng)邊界層比較厚或者分離區(qū)較窄時(shí)發(fā)現(xiàn)了模型應(yīng)力損耗(modeled stress depletion，MSD)問題。Menter和Kuntz最早在2004年基于SST(Shear Stress Transport)模型的混合RANS/LES方法中應(yīng)用了消除MSD的方法，以避免網(wǎng)格精度不夠高時(shí)，在邊界層內(nèi)DES轉(zhuǎn)換為LES模式。2006年Spalart和Deck參照其思想發(fā)展了更通用的DDES方法。為了避免出現(xiàn)對(duì)數(shù)層不匹配問題，節(jié)省計(jì)算量，Shur和Spalart提出IDDES(improved delayed detached-eddy simulation)方法，其結(jié)合了DDES和壁面模型大渦模擬(wall-modeled large eddy simulation,WMLES)方法。針對(duì)翼身組合體飛行器大攻角流場的復(fù)雜性，本文計(jì)算中采用基于SA模型的IDDES類方法。

WMLES模型主要應(yīng)用于非定常和有湍流的流動(dòng)中，它通過長度尺度耦合RANS和LES，其湍流長度尺度定義為

=(1+)+(1-)

(1)

其中，為RANS模型湍流長度尺度，為LES模型湍流長度尺度。混合函數(shù)從0～1變化時(shí)，模型快速的從RANS模式(=1.0)過渡到LES模式(=0)?；旌虾瘮?shù)是用于修正由于RANS和LES交界面相互作用而損耗過多的雷諾應(yīng)力。

IDDES采用了新的亞格子尺度的定義，它同時(shí)依靠網(wǎng)格大小和壁面距離，其湍流長度尺度定義為

(2)

(3)

其中:為延遲函數(shù)，參數(shù)是模型長度尺度與壁面的距離的比值，為到壁面的距離，是渦粘性，, 是速度梯度，是馮卡門常數(shù)。

2.2 數(shù)值計(jì)算方法

由于迎風(fēng)型NND格式形式簡單，數(shù)值耗散較小，且計(jì)算量小，穩(wěn)定性較好，本文數(shù)值模擬中對(duì)無粘項(xiàng)的空間離散都采用迎風(fēng)型NND格式，非定常時(shí)間推進(jìn)采用雙時(shí)間步隱式迭代法。為了提高非定常的時(shí)間計(jì)算精度，同時(shí)又兼具較高的計(jì)算效率，時(shí)間推進(jìn)采用Jameson提出的雙時(shí)間步(dual-time-step)隱式迭代推進(jìn)方法。

對(duì)于本文采用的IDDES混合湍流模型和非定常算法，已通過大量典型算例進(jìn)行了考核和驗(yàn)證，典型算例包括細(xì)長圓錐有攻角繞流、Hummel單三角翼分離流動(dòng)、低速圓柱繞流、跨聲速方腔流動(dòng)、超聲速圓柱底部流動(dòng)、NACA0012翼型俯仰振蕩、NACA0015翼型深失速分離渦模擬等。圖1、圖2、圖3給出了低速圓柱繞流的數(shù)值模擬結(jié)果。

圖1 用密度著色的瞬時(shí)渦量等值面模擬圖(Q=0.01)Fig.1 Isosurfacescoloured by the density of Q=0.01

圖2 壁面平均壓力系數(shù)曲線Fig.2 Average wall pressure coefficient

圖3 流向速度在尾跡中心線分布曲線Fig.3 Velocity profiles along central streamline of wake region

圖2是壁面壓力系數(shù)與實(shí)驗(yàn)和Kravchenko計(jì)算結(jié)果的比較，實(shí)驗(yàn)數(shù)據(jù)來自Norberg，實(shí)驗(yàn)中雷諾數(shù)為4 020。

圖3中的實(shí)驗(yàn)1數(shù)據(jù)來自O(shè)ng和Wallace，實(shí)驗(yàn)2的數(shù)據(jù)來自Lourenco 和Shih，可以看出，在1～4(是圓柱直徑)區(qū)域，本文計(jì)算結(jié)果雖然與Ong和Wallace的實(shí)驗(yàn)結(jié)果差異較小，但在2～4區(qū)域與Lourenco 和Shih的實(shí)驗(yàn)結(jié)果基本一致。對(duì)比上述計(jì)算結(jié)果，可以看出本文采用的模型方法模擬精度較高。其他算例驗(yàn)證結(jié)果參見文獻(xiàn)[25-26]。

2.3 網(wǎng)格和研究模型

脊形前體飛行器外形(翼身組合體WB1)取自Hall實(shí)驗(yàn)。其機(jī)頭采用脊形角7.5°的脊形前體外形，全長=133.35 mm，后接脊形截面柱狀機(jī)身；添加厚度1.27 mm的平板為機(jī)翼，嵌在脊形線上，機(jī)翼后掠55°，具體尺寸及幾何外形如圖4所示，單位為mm。圖5給出了實(shí)驗(yàn)?zāi)Ｐ屯庑螖?shù)據(jù)和本文根據(jù)數(shù)據(jù)點(diǎn)擬合的7.5°前體橫截面模型數(shù)據(jù)曲線，其中是脊形前體截面高度，是寬度，單位均為mm。

圖4 飛行器示意圖Fig.4 Planform dimensions(mm)

圖5 7.5°脊形角計(jì)算模型橫截面與水洞模型數(shù)據(jù)曲線Fig.5 Comparison of the cross-section of the 7.5°calculation model and the water tunnel model

文獻(xiàn)[28-29]考察了網(wǎng)格疏密度對(duì)脊形前體計(jì)算結(jié)果的影響，為了后續(xù)計(jì)算中能給出較為準(zhǔn)確的空間渦結(jié)構(gòu)，滿足IDDES模擬需求，且綜合考慮計(jì)算效率，本文網(wǎng)格參照脊形前體模型G2的網(wǎng)格規(guī)模，采用C-H型多塊對(duì)接網(wǎng)格，網(wǎng)格規(guī)模為2 455萬，沿流向、法向和展向分布為373×201×181，其中機(jī)身上沿流線分布193個(gè)點(diǎn)，網(wǎng)格在壁面、機(jī)頭、機(jī)翼前緣、尾緣及翼尖處加密，參見圖6，物面法向第一層網(wǎng)格間距為1×10mm，確保<10。

圖6 計(jì)算網(wǎng)格示意圖Fig.6 Computation mesh

3 隨雷諾數(shù)變化規(guī)律

飛機(jī)做戰(zhàn)術(shù)機(jī)動(dòng)和巡航時(shí)往往會(huì)在不同的海拔高度，因此本文研究了不同海拔高度時(shí)翼身組合體的氣動(dòng)特性，主要是雷諾數(shù)的影響，其中參考長度取為機(jī)翼根弦長，=190.5 mm時(shí)=1.77×10，1.1×10，0.65×10，分別對(duì)應(yīng)海拔高度是0 km，5 km，10 km，取=04，=30°，= 0°，無量綱時(shí)間步長Δ=001，力矩積分關(guān)于重心軸向位置=222.9 mm。

圖7為脊形前體飛行器模型不同雷諾數(shù)時(shí)氣動(dòng)力系數(shù)收斂曲線。從圖7中可以看出，在所考察雷諾數(shù)下，氣動(dòng)力系數(shù)隨雷諾數(shù)變化很小。

圖7 氣動(dòng)力系數(shù)收斂曲線Fig.7 Convergence curves for aerodynamic coefficients

圖8為用壓力著色的等值面云圖，取=0001，=(-)(+)，和分別為旋轉(zhuǎn)張量和應(yīng)變率張量。本文選取Q渦準(zhǔn)則是由于Q渦準(zhǔn)則消除了邊界層中平均切應(yīng)力和靠近分離的早期剪切層的影響，比渦度準(zhǔn)則更能反映主要的渦結(jié)構(gòu)。

從圖8可以看出，所考察的3個(gè)雷諾數(shù)下前體渦和前緣渦破裂點(diǎn)和強(qiáng)度差別不大，可見在所考察雷諾數(shù)范圍內(nèi)雷諾數(shù)變化時(shí)流場結(jié)構(gòu)趨同，對(duì)稱性較好。文獻(xiàn)[9]中Hall研究了低雷諾數(shù)(=14×10)的水洞數(shù)據(jù)和高雷諾數(shù)(=06×10)風(fēng)洞試驗(yàn)數(shù)據(jù)之間的比較，研究對(duì)象為圓形機(jī)身，盡管由于機(jī)身改變了主翼的有效攻角、水洞和風(fēng)洞模型縮比尺寸不同以及雷諾數(shù)的差異，使水洞和風(fēng)洞試驗(yàn)得出的渦破裂攻角不同，但渦結(jié)構(gòu)類似，變化趨勢相同，所以，水洞實(shí)驗(yàn)數(shù)據(jù)有較高的可參考性。同時(shí)， Hall也指出，雷諾數(shù)在一定范圍內(nèi)變化時(shí)，脊形前體由于脊形邊的存在，前體渦更強(qiáng)，風(fēng)洞與水洞渦結(jié)構(gòu)類似。

圖8 不同雷諾數(shù)時(shí)用壓力著色的Q等值面云圖(Q=0.001)Fig.8 Isosurfacescoloured by the pressure of Q=0.001 at different Reynolds number

4 隨攻角變化規(guī)律

攻角是很重要的來流參數(shù)，為了分析攻角的影響規(guī)律，主要研究了攻角從5°～60°變化時(shí)對(duì)翼身組合體升阻特性、橫航向特性的影響。取來流馬赫數(shù)= 04，攻角分別為= 5°～60°，=0°，無量綱時(shí)間步長Δ=0.01，海拔0 km的大氣參數(shù)，=288.16?；谏鲜隼字Z數(shù)研究，本節(jié)計(jì)算數(shù)據(jù)與Hall實(shí)驗(yàn)數(shù)據(jù)作簡單的對(duì)比，重點(diǎn)分析前體渦和機(jī)翼前緣渦隨攻角演化規(guī)律。

圖9(a)、圖9(c)中給出了攻角15°和30°時(shí)用馬赫數(shù)著色的空間流線，圖9(b)和圖9(d)為攻角15°和30°時(shí)Hall的水洞實(shí)驗(yàn)結(jié)果，與計(jì)算結(jié)果對(duì)比可以看出，前體渦和前緣主渦位置和破裂點(diǎn)相似，表明本文計(jì)算方法具有較高的精度。

圖9 空間渦結(jié)構(gòu)示意圖Fig.9 Space vortex structure

圖10為用壓力著色的等值面云圖，取=0001。從圖10中可以看出,在5°攻角已經(jīng)出現(xiàn)前體渦和前緣渦，15°攻角時(shí)前緣渦從機(jī)翼尾部開始破裂，破裂點(diǎn)隨攻角增加逐漸向機(jī)翼前緣靠近。圖11為不同攻角下用馬赫數(shù)著色的空間流線云圖?？梢钥闯觯?0°攻角時(shí)，前體主渦和前緣主渦沒有耦合，在機(jī)翼上產(chǎn)生二次渦，前緣主渦在機(jī)翼尾部破裂。15°攻角時(shí)，前體渦和前緣渦進(jìn)一步增強(qiáng)，前緣渦向機(jī)身內(nèi)側(cè)靠攏，與前體渦耦合。20°攻角時(shí)，單獨(dú)前體主渦下面有明顯二次渦(圖12)，單獨(dú)前體數(shù)值模擬數(shù)據(jù)參見文獻(xiàn)[27-28]，而翼身組合體模型前體沒有出現(xiàn)明顯二次渦，前體主渦在機(jī)翼中部也與前緣渦耦合，前緣渦破裂點(diǎn)比15°攻角時(shí)大幅靠近機(jī)身頂點(diǎn)。隨著攻角繼續(xù)增大，30°攻角時(shí)，前體渦和前體二次渦向機(jī)身內(nèi)側(cè)移動(dòng)，前體渦和前緣渦進(jìn)一步抬升，前體主渦、前緣二次渦和前緣主渦耦合，使得前緣主渦進(jìn)一步增強(qiáng)；單獨(dú)脊形前體在二次渦外側(cè)形成三次渦，而翼身組合體沒有前體三次渦。40°攻角時(shí)，前體渦破裂點(diǎn)向機(jī)頭方向移動(dòng)，前體二次渦外側(cè)形成三次渦，并與前緣渦耦合。50°攻角時(shí)，前體主渦破裂點(diǎn)繼續(xù)前移，與前緣渦解耦，該攻角時(shí)通常認(rèn)為主翼已經(jīng)嚴(yán)重失速，但由于脊形前體二次渦和三次渦與前緣渦耦合，延遲了主翼前緣渦的完全破裂；與單獨(dú)脊形前體相比，前體渦破裂點(diǎn)大幅后移。60°攻角時(shí)，前緣渦完全破裂。

圖10 不同攻角時(shí)用壓力著色的Q等值云圖(Q=0.001)Fig.10 Isosurfacescoloured by the pressure of Q=0.001 at different angle of attack

圖11 WB1模型空間流線和截面總壓損失云圖Fig.11 Streamline and total pressure loss at the section of WB1 model

圖12 前體模型空間流線和截面總壓損失云圖Fig.12 streamlines and total pressure loss at the section of chine forebody model

圖13為軸向位置背風(fēng)面的表面壓力分布曲線，可見，在所研究攻角范圍內(nèi)，機(jī)翼和前體背風(fēng)面壓力具有較好的對(duì)稱性；由于前體渦和前緣渦的影響，在飛行器背風(fēng)面存在明顯的低壓區(qū)，隨著攻角增大，前體渦和前緣渦向下拉，表面壓力極值變大，但隨著攻角繼續(xù)增大，前緣渦和前體渦破裂，表面壓力極值變小。

圖13 背風(fēng)面截面壓力系數(shù)分布曲線Fig.13 Pressure distributions at constant axial locations of leeward side

圖14為=04時(shí)氣動(dòng)力系數(shù)隨攻角的變化曲線，氣動(dòng)力系數(shù)值是每50步輸出計(jì)算結(jié)果然后求時(shí)間平均。從圖14可以看出：<15°時(shí)，升力系數(shù)隨攻角線性增加；15°<<25°時(shí)，增加變緩，但仍呈線性增加；25°<<40°之間變化緩慢；>40°時(shí)，開始減小。阻力系數(shù)隨攻角增加而增加。隨攻角增加，側(cè)向力系數(shù)、滾轉(zhuǎn)力矩系數(shù)和偏航力矩系數(shù)都趨于零。

圖14 氣動(dòng)力系數(shù)隨攻角變化曲線Fig.14 Aerodynamic coefficients variation with the angle of attack

可見，機(jī)翼前緣渦對(duì)脊形前體主渦有明顯的卷吸作用，極大地延緩了前體主渦的破裂，這導(dǎo)致機(jī)翼前緣渦破裂后，前體渦誘導(dǎo)的升力仍隨攻角增加而增大，整體失速特性有明顯改善。機(jī)翼前緣渦誘導(dǎo)脊形前體主渦更靠近壁面，因而抑制了前體二次渦、三次渦產(chǎn)生和發(fā)展，相應(yīng)地，也延緩了前體二次渦、三次渦的破裂攻角；前體二次渦、三次渦與機(jī)翼前緣渦先后耦合，延遲了機(jī)翼前緣渦的完全破裂，這導(dǎo)致在25°<α<40°之間整體升力變化緩慢，也改善了失速特性，且較單獨(dú)脊形前體，流場對(duì)稱性保持更好。由于機(jī)翼前緣渦誘導(dǎo)低頭力矩，與脊形前體渦誘導(dǎo)抬頭力矩相互抵消，導(dǎo)致翼身組合體的縱向穩(wěn)定性特性，較單獨(dú)脊形前體有明顯改善。

5 結(jié)論

所研究雷諾數(shù)對(duì)脊形前體飛行器氣動(dòng)力系數(shù)和渦流場結(jié)構(gòu)影響較小。攻角變化時(shí)，機(jī)翼前緣渦對(duì)脊形前體主渦有明顯的卷吸作用，極大地延緩了前體主渦的破裂，并誘導(dǎo)脊形前體主渦更靠近壁面，因而抑制了前體二次渦、三次渦產(chǎn)生和發(fā)展，相應(yīng)地，也延緩了前體二次渦、三次渦的破裂攻角；前體二次渦、三次渦與機(jī)翼前緣渦先后耦合，延遲了機(jī)翼前緣渦的完全破裂，整體失速特性有明顯改善。且較單獨(dú)脊形前體，無側(cè)滑時(shí)流場對(duì)稱性保持更好。由于機(jī)翼前緣渦誘導(dǎo)低頭力矩，與脊形前體渦誘導(dǎo)抬頭力矩相互抵消，導(dǎo)致翼身組合體的縱向穩(wěn)定性特性，較單獨(dú)脊形前體有明顯改善。

高機(jī)動(dòng)飛行器一般不以大攻角平飛，而是機(jī)動(dòng)進(jìn)入或經(jīng)過大攻角狀態(tài)。靜態(tài)大攻角繞流是一種近似模擬，后續(xù)將進(jìn)一步研究高機(jī)動(dòng)飛行器的動(dòng)態(tài)氣動(dòng)特性。