圓錐前體大攻角繞流的數(shù)值計(jì)算與分析

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

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

0 引言

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

1 數(shù)值方法

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

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

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

2.1 風(fēng)洞實(shí)驗(yàn)壓力分布

細(xì)長(zhǎng)圓錐前體在大攻角時(shí)測(cè)得的40個(gè)滾轉(zhuǎn)角的壓力系數(shù)存在一些內(nèi)在的規(guī)律，以自由來流速度30m/s、攻角35°、第三截面測(cè)得的數(shù)據(jù)為例進(jìn)行分析。圖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沿相反的趨勢(shì)進(jìn)行變化，基本成鏡像狀態(tài)，而它們?cè)诿恳粋€(gè)滾轉(zhuǎn)角下的平均值接近于一個(gè)常數(shù)，近似為 －0.85。通過對(duì)其它風(fēng)速和攻角下的實(shí)驗(yàn)數(shù)據(jù)進(jìn)行分析，也發(fā)現(xiàn)了同樣的規(guī)律:40個(gè)滾轉(zhuǎn)角下的Cp，min，ave接近于一個(gè)常數(shù)或在一個(gè)很小的范圍內(nèi)波動(dòng)，但不同工況下得到的平均值的量值有所差異。

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

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

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

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

圖2 30m/s時(shí)，對(duì)不同周向網(wǎng)格第三截面的壓力預(yù)測(cè)

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

2.3.1 壓力分布比較

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

圖3 30m/s，計(jì)算與實(shí)驗(yàn)壓力分布對(duì)比，第三截面

2.3.2 云圖分析

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

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

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

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

3 結(jié)論

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

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

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

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

4)攻角大于10°時(shí)，在左右舷處出現(xiàn)主渦，每個(gè)主渦在渦核處都有一個(gè)總壓系數(shù)最小值和一個(gè)無量綱的軸向渦量絕對(duì)值的最大值。該最小總壓系數(shù)隨著攻角的增大而減小，隨著自由來流速度的增大而增大，其縱坐標(biāo)隨著攻角的增大而增大。該軸向渦量最大絕對(duì)值一般會(huì)隨著攻角和自由來流速度的增大而減小，其位置與最小總壓系數(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.