脊形前體飛行器靜態(tài)繞流數值模擬研究

司芳芳，袁先旭，賀旭照，謝昱飛，葉友達

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

1 引言

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

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

盡管已進行大量工作，但由于脊形前體大攻角非定常大范圍湍流分離流的復雜性，目前對脊形前體飛行器的大攻角流動機理認識尚不充分，深入了解其大攻角氣動特性，尤其是脊形前體與翼型相互干擾的氣動特性有著重要的工程應用價值以及學術意義。

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

2 數值模擬方法

2.1 湍流模型

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

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

=(1+)+(1-)

(1)

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

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

(2)

(3)

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

2.2 數值計算方法

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

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

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

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

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

圖2是壁面壓力系數與實驗和Kravchenko計算結果的比較，實驗數據來自Norberg，實驗中雷諾數為4 020。

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

2.3 網格和研究模型

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

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

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

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

圖6 計算網格示意圖Fig.6 Computation mesh

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

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

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

圖7 氣動力系數收斂曲線Fig.7 Convergence curves for aerodynamic coefficients

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

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

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

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

攻角是很重要的來流參數，為了分析攻角的影響規(guī)律，主要研究了攻角從5°～60°變化時對翼身組合體升阻特性、橫航向特性的影響。取來流馬赫數= 04，攻角分別為= 5°～60°，=0°，無量綱時間步長Δ=0.01，海拔0 km的大氣參數，=288.16。基于上述雷諾數研究，本節(jié)計算數據與Hall實驗數據作簡單的對比，重點分析前體渦和機翼前緣渦隨攻角演化規(guī)律。

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

圖9 空間渦結構示意圖Fig.9 Space vortex structure

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

圖10 不同攻角時用壓力著色的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為軸向位置背風面的表面壓力分布曲線，可見，在所研究攻角范圍內，機翼和前體背風面壓力具有較好的對稱性；由于前體渦和前緣渦的影響，在飛行器背風面存在明顯的低壓區(qū)，隨著攻角增大，前體渦和前緣渦向下拉，表面壓力極值變大，但隨著攻角繼續(xù)增大，前緣渦和前體渦破裂，表面壓力極值變小。

圖13 背風面截面壓力系數分布曲線Fig.13 Pressure distributions at constant axial locations of leeward side

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

圖14 氣動力系數隨攻角變化曲線Fig.14 Aerodynamic coefficients variation with the angle of attack

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

5 結論

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

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