圓錐前體大攻角繞流的數值計算與分析

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

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

0 引言

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

1 數值方法

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

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

2 計算結果及討論

2.1 風洞實驗壓力分布

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

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

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

2.2 網格無關性與相關性

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

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

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

2.3 不同模型計算結果比較

2.3.1 壓力分布比較

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

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

2.3.2 云圖分析

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

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

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

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

3 結論

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

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

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

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

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

[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.