王中一，鄭博睿，高 超，方 洪，熊俊濤，劉 鋒，羅時鈞

(1西北工業(yè)大學翼型/葉柵空氣動力學國防科技重點實驗室，西安 710072;2北京臨近空間飛行器系統(tǒng)工程研究所，北京 100076;3美國加州大學尓灣分校，美國加州爾灣 92697-3975)

0 引言

現(xiàn)代飛行器前體和導彈頭部，通常采用細長尖頭旋成體的設計外形。在大攻角下，細長旋成體背風面會產(chǎn)生一對分離渦，當攻角增大到一定程度時，原本對稱的分離渦突然變得非對稱，同時伴隨有方向和大小均無法預估的側(cè)向力和力矩，這對飛行器的操縱性和穩(wěn)定性有很大影響[1－2]。Keener等發(fā)現(xiàn)側(cè)力大小和方向?qū)旤c處的幾何外形很敏感，受雷諾數(shù)或馬赫數(shù)的影響并不大[3]。研究發(fā)現(xiàn)，細長前體分離渦對尖頭處小的擾動非常敏感[4]，而尖鼻點在形狀上近似于圓錐，當?shù)亓鲌鱿喈斢谝粋€正切圓錐處的流場，所以可以通過研究圓錐體繞流來了解細長體上非對稱流場的特征。楊宇等[5]對圓錐前體流動進行了數(shù)值模擬，小于25°攻角得到的計算結(jié)果與實驗對比較好，但大攻角時差異較大。文中在文獻[5]的基礎上進一步優(yōu)化網(wǎng)格，并引入了Gamma theta轉(zhuǎn)捩模型進行數(shù)值模擬。

1 數(shù)值方法

根據(jù)研究內(nèi)容，在CFX中采用3個計算模型進行分析，分別是層流模型，Gamma theta轉(zhuǎn)捩模型和SST湍流模型。其中層流模型以非定常N-S方程作為控制方程，只適合于層流和低雷諾數(shù)的情況;Gamma theta轉(zhuǎn)捩模型多用于轉(zhuǎn)捩預測，它集合了轉(zhuǎn)捩經(jīng)驗關(guān)系式和低雷諾數(shù)湍流模型的優(yōu)勢，通過求解兩個變量(間歇因子和動量厚度雷諾數(shù))的標準輸運方程來實現(xiàn)該目的;SST湍流模型的控制方程為RANS方程，該模型集合了k-ε模型可較好的模擬充分發(fā)展湍流流動的優(yōu)點和k-ω模型可廣泛應用于各種壓力梯度下的邊界層問題的優(yōu)點，在求解邊界層流動時有很高精度，但其對轉(zhuǎn)捩的預測受到質(zhì)疑，因為該模型中阻尼函數(shù)的標定依據(jù)是再現(xiàn)粘性底層的行為，而未考慮轉(zhuǎn)捩的物理機理。

計算初始流場為均勻來流，來流速度為10～70m/s，攻角為15°～35°，基于圓錐段前體底面直徑的雷諾數(shù)為 0.1 ×106～0.7 ×106，壁面光滑且滿足無滑移條件，計算分析類型采用定常態(tài)。計算結(jié)果與Meng XS[6]的實驗數(shù)據(jù)進行了對比。計算模型與文獻[6]尺寸一致，采用圓錐圓柱組合體模型，依次由圓錐段、過渡段和圓柱段三部分組成。圓錐半頂角為10°，為避免數(shù)值計算時尾跡流動對圓錐前體流場造成影響，圓柱段延伸至遠場邊界。圓錐尖端第一個軸向網(wǎng)格長度是圓錐段長度的0.1%，遠離物面的第一個徑向網(wǎng)格高度是0.001mm，為模型圓柱段半徑的10－5倍，徑向遠場邊界距軸線的距離是圓柱段半徑的40倍。實驗模型上的測壓孔分布在沿軸線方向的9個截面上。由于圓錐繞流在相鄰截面的壓力分布差異很小，而第一和第二截面處的測壓孔較少，不能精細地捕捉到壓力分布趨勢，所以文中選用第三截面的實驗數(shù)據(jù)進行對比。

2 計算結(jié)果及討論

2.1 風洞實驗壓力分布

細長圓錐前體在大攻角時測得的40個滾轉(zhuǎn)角的壓力系數(shù)存在一些內(nèi)在的規(guī)律，以自由來流速度30m/s、攻角35°、第三截面測得的數(shù)據(jù)為例進行分析。圖1(a)為壓力系數(shù)隨周向角的變化，圖1(b)為Cp，min，p(左舷最小壓力系數(shù))、Cp，min，s(右舷最小壓力系數(shù))和Cp，min，ave(左右舷最小壓力系數(shù)平均值)隨滾轉(zhuǎn)角的變化，其關(guān)系如下:

可以看出:隨著滾轉(zhuǎn)角從0°變到360°，Cp，min，p和 Cp，min，s沿相反的趨勢進行變化，基本成鏡像狀態(tài)，而它們在每一個滾轉(zhuǎn)角下的平均值接近于一個常數(shù)，近似為 －0.85。通過對其它風速和攻角下的實驗數(shù)據(jù)進行分析，也發(fā)現(xiàn)了同樣的規(guī)律:40個滾轉(zhuǎn)角下的Cp，min，ave接近于一個常數(shù)或在一個很小的范圍內(nèi)波動，但不同工況下得到的平均值的量值有所差異。

圖1 30m/s，α =35°，第三截面

2.2 網(wǎng)格無關(guān)性與相關(guān)性

前期數(shù)值計算發(fā)現(xiàn)，軸向和徑向網(wǎng)格數(shù)的密度不同對計算結(jié)果的影響很小，但圓錐前體的壓力分布復雜，周向網(wǎng)格須分布的足夠密。文中對四種不同的周向網(wǎng)格進行了計算，軸向和徑向網(wǎng)格數(shù)均為81和121，周向網(wǎng)格數(shù)分別為120、180、240和360，稱之為網(wǎng)格1到網(wǎng)格4。當計算步數(shù)為500步時，計算殘差可降低2至3個量級。采用 Gamma theta轉(zhuǎn)捩模型對網(wǎng)格無關(guān)性進行了分析。使用四套網(wǎng)格模擬了30m/s，α=25°和α=35°的流場。四套網(wǎng)格在攻角25°時的收斂解幾乎重合，也就是說25°時的數(shù)值計算存在網(wǎng)格無關(guān)性，如圖2(a)，不同的網(wǎng)格計算結(jié)果與實驗數(shù)據(jù)吻合得都很好。而四套網(wǎng)格在35°時的收斂解相差較大，也就是說35°時不同的計算網(wǎng)格會得到不同的解，其原因可能是使用不同網(wǎng)格計算時產(chǎn)生了不同的截斷誤差，如圖2(b)，計算結(jié)果與實驗數(shù)據(jù)有一定偏差，但4個計算解的左右舷吸力峰處的壓力系數(shù)平均值與實驗數(shù)據(jù)一致，在 －0.85左右(參考2.1)。因此，采用四套網(wǎng)格得到的計算結(jié)果均視為合理的。

綜合考慮計算殘差和網(wǎng)格無關(guān)性等因素，文中采用網(wǎng)格3(81×121×240)進行后續(xù)計算。

圖2 30m/s時，對不同周向網(wǎng)格第三截面的壓力預測

2.3 不同模型計算結(jié)果比較

2.3.1 壓力分布比較

對每一工況下40組不同滾轉(zhuǎn)角的實驗數(shù)據(jù)，從中分別挑選出與相同條件下第三截面計算結(jié)果最為吻合的一組來對比。此處只給出來流速度為30m/s的壓力分布對比圖作為典型進行分析，如圖3。分析得到，層流模型對自由來流速度大于50m/s和攻角大于25°的流場的計算結(jié)果與實驗差異較大;SST湍流模型對自由來流速度小于50m/s和攻角小于25°的流場的計算結(jié)果與實驗差異較大;但對文中所計算的各風速和攻角下的流場，Gamma theta轉(zhuǎn)捩模型的計算結(jié)果與實驗數(shù)據(jù)均吻合較好。這說明，文中各工況下的流動主要處于轉(zhuǎn)捩流動狀態(tài)，適合采用Gamma theta轉(zhuǎn)捩模型來模擬。

圖3 30m/s，計算與實驗壓力分布對比，第三截面

2.3.2 云圖分析

文中分析了不同風速和攻角下，圓錐前體第三截面上的總壓系數(shù)云圖、速度矢量和軸向渦量云圖。這里只給出了來流速度為30m/s，攻角25°和35°時的云圖(圖4、圖5)，分別用Cp0和ωxd/U∞來代表總壓系數(shù)和渦量，計算模型為Gamma theta轉(zhuǎn)捩模型。為了更好的顯示主渦處的渦量，對云圖中渦量的范圍進行了人工選擇，實際渦量范圍要大一些。分析得出:每一個主渦在渦核處均存在一個總壓系數(shù)最小值，但這并不是整個橫截面上的最小總壓系數(shù)，橫截面最小總壓系數(shù)處于緊貼物面的一層很薄的邊界層內(nèi)。同時，每一個主渦在渦核處還存在一個軸向渦量的最大絕對值，該最大值同樣不是整個橫截面上的最大值，橫截面上的最大軸向渦量絕對值也存在于緊貼物面的一層很薄的邊界層內(nèi)。

進一步分析發(fā)現(xiàn)，左右舷主渦渦核處的最小總壓系數(shù)隨著攻角的增大而減小，隨著自由來流速度的增大而增大，其縱坐標則隨著攻角的增大而增大;渦核處的軸向渦量絕對值的最大值一般會隨著攻角和自由來流速度的增大而減小，其位置與最小總壓系數(shù)的位置非常接近。

圖4 30 m/s，α =25°，Gamma theta 轉(zhuǎn)捩模型，第三截面

圖5 30 m/s，α =35°，Gamma theta 轉(zhuǎn)捩模型，第三截面

3 結(jié)論

通過對圓錐前體上的繞流進行數(shù)值模擬，并與實驗數(shù)據(jù)比較，可得出以下結(jié)論:

1)文中涉及的各種流場主要處于轉(zhuǎn)捩區(qū)域，適合于用Gamma theta轉(zhuǎn)捩模型進行模擬。

2)攻角小于25°的流動主要是對稱分離流，計算時存在網(wǎng)格無關(guān)性。攻角大于25°的流動是非對稱的，采用不同的網(wǎng)格進行計算會得到不同的數(shù)值解，其原因可能是采用不同網(wǎng)格計算時產(chǎn)生的截斷誤差不同。在數(shù)值計算中采用不同網(wǎng)格得到的這種效果和實驗中轉(zhuǎn)動滾轉(zhuǎn)角得到的效果類似。

3)對于非對稱流動，不同滾轉(zhuǎn)角下的實驗得到的左右舷最小壓力系數(shù)的平均值接近于一個常數(shù)，而采用不同網(wǎng)格進行計算時得到的左右舷最小壓力系數(shù)的平均值也接近于一個常數(shù)，這兩個值大小相近。因此，可以說文中在大攻角下的數(shù)值模擬結(jié)果與實驗數(shù)據(jù)是相符的。

4)攻角大于10°時，在左右舷處出現(xiàn)主渦，每個主渦在渦核處都有一個總壓系數(shù)最小值和一個無量綱的軸向渦量絕對值的最大值。該最小總壓系數(shù)隨著攻角的增大而減小，隨著自由來流速度的增大而增大，其縱坐標隨著攻角的增大而增大。該軸向渦量最大絕對值一般會隨著攻角和自由來流速度的增大而減小，其位置與最小總壓系數(shù)的位置相近。

[1]Ericsson L E. Sources of high alpha vortex asymmetry at zero sideslip[J]. Journal of Aircraft，1992，29(6):1086－1090.

[2]Lowson M V，Ponton A J C.Symmetry breaking in vortex flows on conical bodies[J]. AIAA Journal，1992，30(6):1576－1583.

[3]Keener E，Chapman G，Cohen L，et al. Side forces on forebodies at high angles of attack and mach numbers from 0.1 to 0.7:Two tangent ogives，paraboloid and cone，NASA TM，X －3438[R].1977.

[4]Zilliac G G，Degani D，Tobak M.Asymmetric vortices on a slender body of revolution[J]. AIAA Journal，1991，29(5):667－675.

[5]Yang Y，Zheng B R，Yang W，et al.Numerical computation of pressure distributions over conical forebody at high angles of attack，AIAA 2012－1237[R].2012.

[6]Meng X S，Jia C，Qiao Z D，et al. Development of flow asymmetry over a 20-Degree circular cone at low speed，AIAA，2007－4118[R].2007.