王方劍，宋玉輝，劉 金，秦 漢，陳 蘭

（中國航天空氣動力技術(shù)研究院，北京 100074）

0 引 言

在機載導彈方面，作為攻防對抗主體的空空導彈，需要比目標飛行器具有更高的機動性和敏捷性，導彈的攻擊方式也由尾后追擊發(fā)展成全向攻擊方式。新型空空導彈面對下一代飛機、下一代空空導彈系統(tǒng)等高機動目標，必須具備全向攻擊能力。采用大攻角和超大攻角飛行，最大限度利用空氣動力實現(xiàn)飛行器飛行軌跡和姿態(tài)變化是直接提高機動性和敏捷性的最優(yōu)手段，由此也帶來超大攻角流動分離問題，包括復雜的非定常旋渦、渦脫落、旋渦非對稱等流動現(xiàn)象，提高對超大攻角復雜非定常流動及其流動機理的認識顯得尤為重要。

早在20世紀50年代，對細長體（錐柱體）分離流動特性的探索就已經(jīng)開始。20世紀80年代以前，由于缺乏強大的計算軟件和先進的風洞試驗方法，學者們關(guān)注點主要集中在附著流和穩(wěn)定分離區(qū)上。隨著技術(shù)的發(fā)展，一些非定常流動特性被揭示出來。Zeiger等對細長旋成體的流動特性進行流動顯示水洞試驗，認為在大攻角時，細長體旋成體上的流動可主要分為3種流動形態(tài)：第1種流動形態(tài)是頭尖部附近的產(chǎn)生的較為集中的渦旋；隨著流動沿軸向向后發(fā)展，流動結(jié)構(gòu)逐漸體現(xiàn)為旋渦脫落流動形態(tài)（第2種流動形態(tài)），渦脫落方向與旋成體呈一定角度；當流動發(fā)展至接近底部時，脫落的旋渦渦軸與旋成體平行，形成第3種流動形態(tài)。Degani等的研究表明，其非定常流動主要包括類卡門渦街的低頻流動、剪切層失穩(wěn)引起的高頻流動，以及介于以上兩種頻率之間的旋渦干擾流動。

本文針對細長體在超大攻角下存在的旋渦非定常性、旋渦脫落、旋渦非對稱等復雜流動現(xiàn)象，采用基于結(jié)構(gòu)網(wǎng)格的雷諾平均/大渦模擬（Reynold averaged Navier-Stokes/large eddy simulation，RANS/LES）混合數(shù)值方法，開展細長體復雜氣動非定常特性研究，分析新型導彈在超大攻角下的流動機理，為超大攻角氣動設(shè)計提供理論支撐。

1 超大攻角數(shù)值模擬方法

細長體在較大攻角時能夠誘導較為明顯的非定常旋渦結(jié)構(gòu)。傳統(tǒng)的RANS 模型可以在不苛刻的網(wǎng)格要求下較好地計算附著流和小分離流動，但是對脫體的大尺度分離流的應用一直存在難以克服的困難。其根源在于空間旋渦流場中RANS 渦粘性模型假設(shè)導致耗散過大，對旋渦非定常流動帶來過沖的阻尼抑制作用。大渦模擬（LES）在模擬湍流的細節(jié)和流動分離等方面都有很好的精準度，但其對高雷諾數(shù)下物面流動的計算量是當前計算機水平所不能企及的，故其與工程應用還有相當大的距離。近年來興起的多種RANS-LES 混合方法綜合了RANS 方法對近壁附著流刻畫和LES 方法對空間大分離模擬的各自優(yōu)點，其中應用最廣泛的是脫體渦模擬（detached eddy simulation，DES）類方法。原始的DES 方法將SST模型的湍動動能輸運方程耗散項改寫為與網(wǎng)格尺度和湍流長度相關(guān)的DES 形式。但該方法存在模型雷諾應力損耗（modeled stress depletion，MSD）問題，進而帶來網(wǎng)格誘導分離等非物理問題。因此，本文發(fā)展了延遲脫體渦模擬（DDES）方法，通過設(shè)置延遲函數(shù)，克服了網(wǎng)格加密時可能形成的附面層內(nèi)LES 啟動過早的問題。

因此，本文采用DDES 方法對流動控制方程進行求解，算法設(shè)置如下：

1）空間項離散：離散格式為基于守恒律的上游中心格式（monotone upstream-centered schemes for conservation laws，MUSCL）方法插值方法的FDSRoe格式。

2）時間項推進：時間推進格式采用的是隱式LUSGS 方法，保證了較高的計算效率與精度，非定常計算采用雙時間步方法，并利用當?shù)貢r間步長加速收斂。

3）湍流模擬：湍流模擬中采用基于剪切應力運輸（shear stress transport，SST）湍流模型的DDES方法。

2 細長體外形與網(wǎng)格

2.1 細長體外形

計算采用的外形如圖1，彈身采用的圓柱外形，鴨式舵面為“X”型布局，舵面為梯形翼，彈身頭部采用的是類半球型頭部設(shè)計，導彈質(zhì)心為/=（0.55，0，0）。坐標設(shè)置為軸向后，軸向上。

圖1 細長體外形Fig.1 Slender body configuration

2.2 數(shù)值計算網(wǎng)格

計算網(wǎng)格采用結(jié)構(gòu)網(wǎng)格進行求解，在靠近標模物面處采用O 型網(wǎng)格拓撲，在遠離物面處采用H 型網(wǎng)格拓撲，外邊界距離的設(shè)置以能夠滿足自由來流條件為判斷準則。貼近物面的第一層網(wǎng)格厚度保持～1，以確保邊界層的準確模擬。為了能夠精確地捕捉細長體外形背風側(cè)的分離渦流動，在模型的法向進行了整體加密，在模型背風側(cè)展向網(wǎng)格進行適當加密（如圖3），總網(wǎng)格量為1 500萬。

圖2 物面分布情況Fig.2 Surface grid

圖3 周向網(wǎng)格分布情況Fig.3 Circumferential direction grid

3 超大攻角非定常氣動力特性

非定常數(shù)值模擬的計算工況參照靜態(tài)風洞測力試驗的試驗工況，=0.6、0.8、1.15，=0°～180°，△=15°，共13 個攻角，側(cè)滑角=0°，雷諾數(shù)=2.93×10/m。計算采用非定常DDES方法，每個計算物理時間步為1×10，共計算1 200物理時間步。

圖4～7 為細長體外形在=0.6 氣動力系數(shù)隨攻角變化情況，圖中實線是模型非定常氣動力的平均值，陰影部分代表非定常脈動量的量值范圍。圖4是法向力系數(shù)C隨攻角變化情況，圖中顯示，試驗模型的非定常脈動量較大的攻角范圍是=75°～150°。圖5是俯仰力矩系數(shù)隨攻角變化情況，非定常氣動力的脈動影響較大，在某些攻角下（如=120°），其俯仰力矩瞬時值甚至會發(fā)生變號的現(xiàn)象。圖6～7 是橫側(cè)向力矩系數(shù)、C隨攻角變化情況，在攻角=30°～150°都體現(xiàn)出了較大幅值的非定常特性。從以上圖中可以看出，試驗模型的氣動力在攻角范圍45°～165°內(nèi)具有較為強烈的非定常性，側(cè)向力非定常脈動瞬時值能夠達到法向力平均的1/4～1/2。

圖4 法向力系數(shù)隨攻角變化（Ma=0.6）Fig.4 Normal force coefficient vs α（Ma=0.6）

圖5 俯仰力矩系數(shù)隨攻角變化（Ma=0.6）Fig.5 Pitch moment coefficient vs α（Ma=0.6）

圖6 滾轉(zhuǎn)力矩系數(shù)隨攻角變化（Ma=0.6）Fig.6 Roll moment coefficient vs α（Ma=0.6）

圖7 偏航力矩系數(shù)隨攻角變化（Ma=0.6）Fig.7 Yaw moment coefficient vs α（Ma=0.6）

圖8～11 為細長體外形在=0.8 氣動力系數(shù)隨攻角變化情況。圖8是法向力系數(shù)脈動量隨攻角的變化情況，在攻角60°～120°之間體現(xiàn)出了較為明顯的非定常脈動特性。圖9是俯仰力矩系數(shù)脈動量隨攻角變化情況，在較小攻角下（0°～30°），非定常脈動量非常小，當攻角達到45°時，非定常脈動量突然增大，在攻角范圍45°～135°范圍內(nèi)，其非定常脈動量都體現(xiàn)得較為明顯。圖10～11 是橫航向氣動力矩非定常脈動量隨攻角變化情況，其非定常脈動量隨攻角變化趨勢是一致的，即在攻角60°及120°附近脈動量達到最大，在攻角90°附近脈動量有一定程度的減小。

圖8 法向力系數(shù)隨攻角變化（Ma=0.8）Fig.8 Normal force coefficient vs α（Ma=0.8）

圖9 俯仰力矩系數(shù)隨攻角變化（Ma=0.8）Fig.9 Pitch moment coefficient vs α（Ma=0.8）

圖10 滾轉(zhuǎn)力矩系數(shù)隨攻角變化（Ma=0.8）Fig.10 Roll moment coefficient vs α（Ma=0.8）

圖11 偏航力矩系數(shù)隨攻角變化（Ma=0.8）Fig.11 Yaw moment coefficient vs α（Ma=0.8）

圖12～15 為細長體外形在=1.15 氣動力系數(shù)隨攻角變化情況，圖中顯示其非定常脈動量與=0.6、0.8 相比有很大程度的減小，反映出了=1.15 下，其非定常流動脈動相對較弱的流動特性。圖12為法向力系數(shù)脈動量隨攻角變化情況，圖中顯示，其非定常脈動量已較弱，在45°～60°攻角范圍內(nèi)，非定常脈動量相對較大。圖13是俯仰力矩系數(shù)非定常脈動量隨攻角變化情況，圖中顯示當攻角在45°～60°區(qū)間內(nèi)以及120°附近有相對明顯的脈動量。圖14～15 是橫航向氣動力矩脈動量隨攻角變化情況，圖中顯示攻角在45°附近以及165°附近體現(xiàn)出了明顯的非定常脈動特性。

圖12 法向力系數(shù)隨攻角變化（Ma=1.15）Fig.12 Normal force coefficient vs α（Ma=1.15）

圖13 俯仰力矩系數(shù)隨攻角變化（Ma=1.15）Fig.13 Pitch moment coefficient vs α（Ma=1.15）

圖14 滾轉(zhuǎn)力矩系數(shù)隨攻角變化（Ma=1.15）Fig.14 Roll moment coefficient vs α（Ma=1.15）

圖15 偏航力矩系數(shù)隨攻角變化（Ma=1.15）Fig.15 Yaw moment coefficient vs α（Ma=1.15）

4 超大攻角非定常頻率特性

細長體非定常流動有非常豐富的頻率特性。圖16中，左側(cè)兩張圖是=60°時法向力與側(cè)向力系數(shù)隨時間變化曲線，右側(cè)是其對應的幅值譜。圖中顯示，法向力與側(cè)向力系數(shù)都存在較為強烈的非定常脈動。從幅值譜來看，法向力沒有發(fā)現(xiàn)明顯的主頻，但是側(cè)向力系數(shù)體現(xiàn)出了明顯的主頻為820 Hz。由于在不同攻角下縱向氣動力系數(shù)脈動都體現(xiàn)出了相似的特征，而不同攻角下的橫側(cè)向氣動力系數(shù)則體現(xiàn)了不同的特征，所以本章著重分析橫側(cè)向氣動力系數(shù)的頻率特性。

圖16 法向力與側(cè)向力系數(shù)隨時間變化（Ma=0.6，α=60°）Fig.16 Normal force and side force coefficient vs time（Ma=0.6，α=60°）

圖17是=0.6，=60°工況下，側(cè)向力系數(shù)時間歷程計算值與試驗數(shù)據(jù)對比，從圖中可以看到計算值與試驗值對比較好，計算值稍大。圖18是側(cè)向力系數(shù)的功率譜密度（power spectral density，PSD）計算值與試驗數(shù)據(jù)對比，圖中顯示，計算值與試驗值獲得的側(cè)向力系數(shù)主頻一致，均為820 Hz。

圖17 Cy時間歷程對比（Ma=0.6，α=60°）Fig.17 Cy vs time（Ma=0.6，α=60°）

圖18 側(cè)向力PSD對比（Ma=0.6，α=60°）Fig.18 Comparison of PSD of Cy（Ma=0.6，α=60°）

圖19是=0.6 不同攻角下的側(cè)向力系數(shù)隨時間變化曲線及其頻譜，圖中左側(cè)為側(cè)向力系數(shù)隨時間變化圖，右側(cè)為其對應的幅值譜。圖中顯示，當攻角=30°時，側(cè)向力系數(shù)沒有體現(xiàn)出明顯的主頻，主要呈現(xiàn)出一種寬頻的特征。當攻角在=60°～150°范圍內(nèi)，側(cè)向力系數(shù)均體現(xiàn)出了較為明顯的主頻。值得注意的是，當攻角=90°時，體現(xiàn)出了兩個比較集中的頻率，分別為689 Hz 與831 Hz，下文中將對這個問題進行分析。

圖19 不同攻角下側(cè)向力隨時間變化曲線及其頻譜（Ma=0.6）Fig.19 Side force coefficient vs time and frequency spectrum at different α（Ma=0.6）

圖20是=0.8 時不同攻角下的側(cè)向力系數(shù)隨時間變化曲線及其頻譜。圖中顯示，在攻角30°工況下，沒有體現(xiàn)出明顯的主頻，其非定常能量主要集中在500～700 Hz 之間。在攻角60°～150°范圍內(nèi)，均體現(xiàn)出了較為集中的主頻現(xiàn)象。

圖20 不同攻角下側(cè)向力隨時間變化曲線及其頻譜（Ma=0.8）Fig.20 Side force coefficient vs time and frequency spectrum at different α（Ma=0.8）

圖21是=1.15 時不同攻角下的側(cè)向力系數(shù)隨時間變化曲線及其頻譜，從圖中可以看到，與=0.6、0.8 不同的是，其在攻角30°～90°范圍內(nèi)都沒有體現(xiàn)出集中的主頻特性。攻角為30°時，其非定常能量主要集中在500～800 Hz 范圍內(nèi)；攻角為60°時其非定常脈動能量范圍主要集中在1 500 Hz范圍內(nèi)；攻角為90°時，體現(xiàn)出的一種寬頻的頻率特征，攻角為120°、150°時，體現(xiàn)出了較為集中的主頻特性。

圖21 不同攻角下側(cè)向力隨時間變化曲線及其頻譜（Ma=1.15）Fig.21 Side force coefficient vs time and frequency spectrum at different α（Ma=1.15）

圖22是=0.6、0.8、1.15工況下，側(cè)向力系數(shù)主頻隨攻角變化情況。從圖中可以看到，隨著攻角的增加，其脈動頻率主要呈現(xiàn)為先增大后減小的現(xiàn)象。當=0.6時，其主頻范圍在400～1 000 Hz之間。當=0.8時，其主頻范圍在500～1 200 Hz之間，=0.8時在各個攻角下的流動主頻均比=0.6時大。當=1.15時，只有在45°、120°、150°攻角下流動有明顯的主頻，其他攻角下并沒有發(fā)現(xiàn)，其頻率范圍在300～1 400 Hz。

圖22 側(cè)向力主頻隨攻角變化Fig.22 Dominant frequency of Cy vs α

圖23是斯特勞哈爾數(shù)（）隨攻角的變化情況。圖中斯特勞哈爾數(shù)的計算中，速度選取的是橫截面的來流速度sin，圖中顯示，不同攻角下的斯特勞哈爾數(shù)量值范圍為0.08～0.35 左右。=0.6 與=0.8工況下，其斯特勞哈爾數(shù)隨攻角變化趨勢接近。當攻角在75°～120°區(qū)間內(nèi)，=0.6 所對應的斯特勞哈爾數(shù)更大，為0.2 左右；當攻角在135°～165°區(qū)間內(nèi)，兩者斯特勞哈爾較為接近，量值在0.25～0.32。=1.15時，攻角150°下斯特勞哈爾較低，只有0.08左右。

圖23 側(cè)向力St隨攻角變化Fig.23 St of Cy vs α

為了更進一步地研究非定常頻率特性，以=0.6 為例，采用數(shù)值模擬方法，在模型的表面布置測壓點，其中測壓點的分布如圖24所示。

圖24 模型測壓點布置情況Fig.24 Monitored points location on the surface

圖25～27是選取較為典型的3個攻角下的局部壓力脈動曲線，每張圖中有3張小圖，每張小圖里上面的曲線為當?shù)貕毫ο禂?shù)隨時間變化的情況，下面的曲線為與其相對應的幅值譜。圖25展示的是當攻角=60°時典型測壓點的壓力脈動情況，圖中顯示：流動在頭部附近體現(xiàn)出了兩個頻段，一個是0～300 Hz 的較低頻段，一個是800 Hz 左右的主頻；當流動發(fā)展到舵面上時，舵面上的流動沒有體現(xiàn)出明顯的主頻；當流動繼續(xù)向后發(fā)展，其壓力脈動慢慢恢復為較為明顯的主頻，為820 Hz。圖26展示的是當攻角=90°時典型測壓點的壓力脈動情況，圖中顯示：在模型的頭部附近沒有看到明顯的主頻；當流動發(fā)展到模型中部時，體現(xiàn)出了明顯的主頻，為831 Hz；當流動繼續(xù)向后發(fā)展，出現(xiàn)主頻降低的現(xiàn)象；在模型接近尾部的地方，流動主頻降低為689 Hz。圖27所示為攻角=150°時典型測壓點的壓力脈動情況，可以看出，在模型頭部、中部、舵面部分均體現(xiàn)出了很明顯的主頻，而且其幅值譜非常干凈，并且在模型的舵面還體現(xiàn)出了除主頻外的多倍頻率情況。

圖25 不同測壓點壓力脈動情況（α=60°）Fig.25 Fluctuating pressure of different points（α=60°）

圖26 不同測壓點壓力脈動情況（α=90°）Fig.26 Fluctuating pressure of different points（α=90°）

圖27 不同測壓點壓力脈動情況（α=150°）Fig.27 Fluctuating pressure of different points（α=150°）

從非定常頻率特性分析中可以發(fā)現(xiàn)，試驗模型的橫側(cè)向氣動力（側(cè)向力、偏航力矩、滾轉(zhuǎn)力矩）在攻角45°～165°范圍內(nèi)有比較明顯的主頻，其范圍為400～1 000 Hz；氣動力的主頻特性主要由旋成體彈身帶來，彈翼貢獻很小。

5 超大攻角非定常流動特性

超大攻角的非定常氣動力以及頻率特性的本質(zhì)是非定常的流動，本章將對細長體的非定常流動進行分析。

圖28所示為攻角=60°時的非定常瞬時流動，圖中左側(cè)有采用Q 準則進行提取得到的等值面，采用壓力來著色，右側(cè)7 張圖片是從頭部開始不同截面的流動渦量圖。從圖中可以看到，該外形整體呈現(xiàn)出了非對稱旋渦流動，旋渦流動是從前向后逐步發(fā)展，當渦位飄得足夠高時，其下方又會生成一個新的左側(cè)渦，同時在舵面下游以及尾部下游均出現(xiàn)較為明顯的旋渦脫落現(xiàn)象。右側(cè)的渦量圖顯示：模型流動在頭部時（即=0.04），旋渦還在起始發(fā)展中，旋渦的非對稱現(xiàn)象還不明顯；當流動發(fā)展到舵面時（=0.15），體現(xiàn)出了較大的分離流動；當流動發(fā)展至=0.21 時，由舵面誘導的大分離流動依然存在；隨著流動繼續(xù)沿軸向向后發(fā)展，流動開始逐漸演變成類卡門渦街的旋渦流動，出現(xiàn)很明顯的非對稱旋渦現(xiàn)象，靠近物面的旋渦卷起后逐漸遠離物面，同時渦量逐漸減弱，直至脫落。

圖28 非定常瞬時流動（α=60°）Fig.28 Unsteady transient flow（α=60°）

圖29所示為攻角=90°時的非定常瞬時流動，圖中顯示：模型流動在頭部（=0.04）時，已經(jīng)體現(xiàn)出了大分離、旋渦脫落等現(xiàn)象；當流動發(fā)展至=0.42 時，其流動特性演變?yōu)轭惪ㄩT渦街流動，同時隨著時間的變化進行左右切換。與攻角=60°不同的是，攻角=90°時旋渦渦位較高，并且渦量較低，這主要是由于在此攻角下軸向流動速度很低。值得注意的是，當流動發(fā)展至模型尾部（=0.98）時，旋渦的非對稱性減弱，近似于對稱渦，并且沒有脫落。這是由于流動繞過尾部在背風側(cè)體現(xiàn)出了一定的軸向速度，從而使得旋渦保持一定的穩(wěn)定性。這也在一定程度上解釋了在此攻角下，模型中部的局部壓力脈動主頻（831 Hz）與尾部主頻（689 Hz）不同的原因。

圖29 非定常瞬時流動（α=90°）Fig.29 Unsteady transient flow（α=90°）

圖30所示為攻角=120°時非定常瞬時流動，圖中顯示：當攻角大于90°時，流動是從模型尾部（=0.98）開始向頭部發(fā)展，尾部的流動主要體現(xiàn)為分離流動；隨著流動向頭部發(fā)展，逐漸體現(xiàn)出類卡門渦街的旋渦特征；當流動發(fā)展至舵面（=0.21）時，體現(xiàn)為大分離流動，非對稱旋渦流動逐漸消失；當流動發(fā)展至頭部（=0.04）時，受舵面的影響，其流動是大分離流動，并伴有一定的旋渦脫落現(xiàn)象。

圖30 非定常瞬時流動（α=120°）Fig.30 Unsteady transient flow（α=120°）

圖31所示為攻角=150°時非定常瞬時流動，流動依然從尾部開始發(fā)展，逐漸發(fā)展為非對稱旋渦流動，由于此攻角存在較大的軸向速度，所以非對稱旋渦能夠較好的保持，沒有出現(xiàn)大范圍的旋渦脫落現(xiàn)象。同時，由于旋渦能量較強，當旋渦發(fā)展至舵面處時，其與物面距離較近，仍然有較強的誘導能力，所以在舵面上也能夠誘導出具有主頻的非定常旋渦流動。并且由于攻角效應，旋渦渦位不會太高，在旋渦下方不會誘導出新的旋渦流動，因此在整個彈身的左右渦位是保持一致的，不會出現(xiàn)在不同截面上左右渦位切換的現(xiàn)象。

圖31 非定常瞬時流動（α=150°）Fig.31 Unsteady transient flow（α=150°）

6 結(jié) 論

本文采用非定常數(shù)值模擬DDES 方法，計算了細長體在=0.6、0.8、1.15，攻角=0°～180°，側(cè)滑角=0°狀態(tài)下的非定常流場，并從超大攻角非定常氣動力特性、非定常頻率特性、非定常旋渦流動特性幾個方面進行深入分析，得到以下幾點主要結(jié)論：

1）試驗模型的六分量氣動力在攻角45°～165°范圍內(nèi)均具有較為強烈的非定常性，側(cè)向力的非定常脈動幅值尤為強烈，其瞬時量值為法向力的1/4～1/2。

2）試驗模型的橫側(cè)向氣動力（側(cè)向力、偏航力矩、滾轉(zhuǎn)力矩）在攻角45°～165°范圍內(nèi)有比較明顯的主頻，范圍在400～1 400 Hz，縱向氣動力/力矩無明顯主頻。

3）氣動力的主頻特性主要由旋成體彈身帶來，彈翼貢獻很小。

4）背風側(cè)復雜旋渦流動是非定常氣動力主要原因，其主要體現(xiàn)為旋渦生成、旋渦切換、渦脫落等。