馮軍紅，魯軍勇，李 開(kāi)，李湘平

(海軍工程大學(xué) 艦船綜合電力技術(shù)國(guó)防科技重點(diǎn)實(shí)驗(yàn)室， 湖北 武漢 430033)

電磁發(fā)射超高速?gòu)椡枋侵咐秒姶拍芗夹g(shù)發(fā)射一體化彈丸至超高速狀態(tài)，其具有初速高(可達(dá)6Ma以上)、射程遠(yuǎn)(可達(dá)200 km以上)、威力大以及發(fā)射能量精確可調(diào)等優(yōu)勢(shì)，是一種先進(jìn)的動(dòng)能殺傷武器[1-3]。作為大空域飛行的超高速?gòu)椡瑁錃鈩?dòng)減阻以及飛行穩(wěn)定特性是關(guān)鍵，值得深入研究。隨著計(jì)算流體力學(xué)的發(fā)展，采用數(shù)值方法計(jì)算超高速?gòu)椡璧臍鈩?dòng)力及力矩逐漸被國(guó)內(nèi)外學(xué)者普遍接受，成為研究彈丸氣動(dòng)特性的重要手段[4-9]。

相對(duì)于傳統(tǒng)的低速?gòu)椡?，超高速?gòu)椡璧娘w行空域和速域范圍較大，幾乎覆蓋了大氣層的不同層段，涉及稠密大氣層減阻、稀薄氣體效應(yīng)以及黏性干擾效應(yīng)等氣動(dòng)現(xiàn)象，使得超高速?gòu)椡铓鈩?dòng)力、力矩以及壓心變化較大，存在大空域氣動(dòng)減阻、氣動(dòng)熱防護(hù)和飛行穩(wěn)定性等諸多問(wèn)題。近年來(lái)，國(guó)內(nèi)南京理工大學(xué)等[9-10]單位開(kāi)展了超高速?gòu)椡璧臍鈩?dòng)特性研究，重點(diǎn)分析了彈丸的氣動(dòng)力和氣動(dòng)熱特性，但缺乏對(duì)彈丸再入飛行段大攻角下的氣動(dòng)特性分析。其他研究單位也較少涉及這一方面的研究。

本文以電磁發(fā)射超高速?gòu)椡铻檠芯繉?duì)象，采用雷諾平均數(shù)值方法計(jì)算了超高速?gòu)椡璧娜S氣動(dòng)流場(chǎng)，并通過(guò)風(fēng)洞試驗(yàn)數(shù)據(jù)驗(yàn)證了數(shù)值計(jì)算精度。重點(diǎn)分析了超高速?gòu)椡柙偃攵未蠊ソ窍碌臍鈩?dòng)特性，得到了彈丸的非線(xiàn)性氣動(dòng)參數(shù)，為超高速?gòu)椡柙偃攵物w行穩(wěn)定性研究提供了氣動(dòng)數(shù)據(jù)，也為彈丸外形設(shè)計(jì)提供了參考。

1 數(shù)值仿真方法

1.1 控制方程

不考慮體積力，采用Favre平均定常張量形式的N-S方程為：

(1)

(2)

(3)

1.2 湍流模型

對(duì)于彈丸的外部流場(chǎng)計(jì)算，最為常用的湍流模型是S-A(Spalart-Allmaras)單方程模型，詳細(xì)的模型說(shuō)明可見(jiàn)文獻(xiàn)[11]。但是S-A湍流模型在預(yù)測(cè)分離流方面誤差較大，石磊等[12]對(duì)比分析了S-A和k-ωSST湍流模型在預(yù)測(cè)高速旋轉(zhuǎn)彈丸氣動(dòng)特性的精度，指出SST湍流模型在可壓縮分離流預(yù)測(cè)方面優(yōu)于S-A模型。SST湍流模型適合用于高速流動(dòng)和由于逆壓梯度造成的分流流動(dòng)[13]，并得到了廣泛的應(yīng)用[9,14]，有關(guān)k-ωSST湍流模型的詳細(xì)說(shuō)明可見(jiàn)文獻(xiàn)[15]。本文擬采用這兩種湍流模型計(jì)算超高速?gòu)椡璧娜S流場(chǎng)，并將計(jì)算結(jié)果與風(fēng)洞試驗(yàn)進(jìn)行對(duì)比，同時(shí)，采用k-ωSST湍流模型計(jì)算超高速?gòu)椡柙偃攵未蠊ソ窍碌臍鈩?dòng)流場(chǎng)。

1.3 計(jì)算模型及網(wǎng)格

為減小氣動(dòng)阻力，采用無(wú)翼式布局，流線(xiàn)型外形，見(jiàn)圖1。以彈丸總長(zhǎng)L為特征長(zhǎng)度，彈丸長(zhǎng)細(xì)比約為7，尾翼呈“×”型結(jié)構(gòu)，翼厚度為2 mm，斜置安裝于彈丸底部。由于尾翼斜置的原因，使得彈丸并不完全對(duì)稱(chēng)，因此，建立了三維全彈模型，同時(shí)考慮超聲速來(lái)流，前場(chǎng)流場(chǎng)區(qū)域?yàn)?.5倍彈長(zhǎng)，后場(chǎng)區(qū)域?yàn)?倍彈長(zhǎng)，徑向?yàn)?.5倍彈長(zhǎng)。采用ICEM CFD軟件進(jìn)行網(wǎng)格劃分，在彈丸和尾翼壁面處進(jìn)行網(wǎng)格加密處理，保證壁面y+值位于30～35之間，網(wǎng)格過(guò)渡比為1.2，如圖1所示。

(a) 計(jì)算模型及對(duì)稱(chēng)面網(wǎng)格(a) Calculation model and grid of symmetry plane

(b) 三維網(wǎng)格結(jié)構(gòu)(b) Three-dimensional grid (c) 尾部剖面網(wǎng)格(c) Tail section grid圖1 計(jì)算模型及網(wǎng)格示意圖Fig.1 Computation model and grids distribution

利用有限體積法、基于隱式方法數(shù)值求解定常的N-S方程，空間上采用二階迎風(fēng)格式進(jìn)行離散，采用AUSM+通量分裂格式，彈丸和彈翼表面采用無(wú)滑移的絕熱壁面條件，在來(lái)流入口、出口以及周向邊界上均采用壓力遠(yuǎn)場(chǎng)邊界條件。整個(gè)計(jì)算工況見(jiàn)表1。基準(zhǔn)工況的大氣壓力為101 325 Pa，溫度為288.15 K，再入工況可取12 km海拔高度的大氣參數(shù)，大氣壓力為19 399 Pa，溫度為216.65 K。

表1 計(jì)算工況

2 計(jì)算結(jié)果分析

2.1 數(shù)值驗(yàn)證

2.1.1 網(wǎng)格無(wú)關(guān)性驗(yàn)證

為驗(yàn)證計(jì)算網(wǎng)格無(wú)關(guān)性，在現(xiàn)有模型基礎(chǔ)上劃分了3套不同數(shù)量的網(wǎng)格，網(wǎng)格加密方式相同，在彈體壁面附近進(jìn)行網(wǎng)格加密，粗網(wǎng)格數(shù)量為160萬(wàn)，最小網(wǎng)格尺寸為0.05 mm，最大網(wǎng)格尺寸為4 mm；適中網(wǎng)格數(shù)量為260萬(wàn)，最小網(wǎng)格尺寸為0.02 mm，最大網(wǎng)格尺寸為1.5 mm；精細(xì)網(wǎng)格數(shù)量為330萬(wàn)，最小網(wǎng)格尺寸為0.01 mm，最大網(wǎng)格尺寸為1 mm。數(shù)值計(jì)算結(jié)果可見(jiàn)表2，升力系數(shù)的相對(duì)誤差較小。俯仰力矩系數(shù)的相對(duì)誤差從210萬(wàn)加密到330萬(wàn)時(shí)計(jì)算結(jié)果相差0.06%，這說(shuō)明網(wǎng)格影響趨于收斂，本文選用適中網(wǎng)格數(shù)量滿(mǎn)足計(jì)算要求。

表2 網(wǎng)格無(wú)關(guān)性數(shù)值驗(yàn)證

2.1.2 數(shù)值結(jié)果試驗(yàn)驗(yàn)證

采用風(fēng)洞試驗(yàn)數(shù)據(jù)驗(yàn)證數(shù)值仿真結(jié)果，圖2～4顯示了數(shù)值計(jì)算和風(fēng)洞試驗(yàn)的對(duì)比結(jié)果，分別對(duì)比了高馬赫數(shù)下S-A和k-ωSST湍流模型的預(yù)測(cè)精度。可看出，法向力系數(shù)和俯仰力矩系數(shù)吻合較好，最大誤差為1.8%。同時(shí)兩種湍流模型的預(yù)測(cè)精度相當(dāng)。相對(duì)于S-A湍流模型，k-ωSST湍流模型的最大誤差為1.7%，略高于S-A模型的預(yù)測(cè)精度，這是因?yàn)榉ㄏ蛄ο禂?shù)與俯仰力矩系數(shù)主要由氣體壓差產(chǎn)生，黏性項(xiàng)貢獻(xiàn)較小，而兩種湍流模型對(duì)壓力項(xiàng)預(yù)測(cè)精度較高，所以誤差較小。

圖2 軸向力系數(shù)隨攻角變化曲線(xiàn)(Ma=4)Fig.2 Variation curve of axial force coefficient with angle of attack (Ma=4)

針對(duì)軸向力系數(shù)，兩種湍流模型的預(yù)測(cè)精度較低，尤其是k-ωSST湍流模型，最大誤差達(dá)到14%，而S-A湍流模型的最大誤差為4.7%。這是由于相對(duì)于法向力系數(shù)，軸向力系數(shù)的黏性項(xiàng)貢獻(xiàn)較大，而湍流渦黏模型在預(yù)測(cè)湍流黏性項(xiàng)方面的誤差較大，使得計(jì)算結(jié)果偏差較大。而S-A模型在計(jì)算外流場(chǎng)方面具有優(yōu)勢(shì)，因此得到了廣泛的應(yīng)用。但k-ωSST模型在計(jì)算分離流方面更具有優(yōu)勢(shì)[12-15]。因此，本文在基準(zhǔn)工況下采用S-A模型；在再入工況下，由于彈丸攻角較大，會(huì)產(chǎn)生分離流，采用k-ωSST湍流模型進(jìn)行計(jì)算。

圖3 法向力系數(shù)隨攻角變化曲線(xiàn)(Ma=4)Fig.3 Variation curve of normal force coefficient with angle of attack (Ma=4)

圖4 俯仰力矩系數(shù)隨攻角變化曲線(xiàn)(Ma=4)Fig.4 Variation curve of pitching moment coefficient with angle of attack (Ma=4)

2.2 基準(zhǔn)工況氣動(dòng)分析

采用S-A湍流模型的計(jì)算結(jié)果分析了彈丸氣動(dòng)流場(chǎng)特性。圖5顯示了攻角3°工況下，不同截面的流場(chǎng)壓力分布云圖。圖5(a)表示0°截面，即彈體的鉛垂面，圖5(b)表示45°截面，即繞彈中心軸旋轉(zhuǎn)45°的截面。來(lái)流在彈體頭部以及彈翼前緣形成脫體激波，但由于彈丸鈍角較小，脫體激波距離較小。激波向周?chē)臻g延伸，且強(qiáng)度逐漸降低。針對(duì)高馬赫數(shù)來(lái)流工況，來(lái)流激波角較小，使得激波后的氣流貼近物體壁面。氣流經(jīng)過(guò)彈丸底部先膨脹后壓縮，即底部流場(chǎng)除了由于高速氣流引射作用形成的低壓回流區(qū)外，還會(huì)形成再壓縮波或激波，流場(chǎng)總壓進(jìn)一步損失。

(a) 0° 截面(a) 0° section

(b) 45° 截面(b) 45° section圖5 不同截面的流場(chǎng)壓力分布云圖(基準(zhǔn)工況)Fig.5 Pressure contour in different cross sections (base case)

對(duì)于超高速?gòu)椡?，氣?dòng)阻力主要由三部分組成，即摩擦阻力、波系阻力以及底部阻力。從圖5中可清楚看出波系阻力和底部阻力的產(chǎn)生機(jī)理。波系阻力主要是由于氣流經(jīng)過(guò)激波后總壓損失，而底部阻力是由于彈丸底部低壓造成的。圖6顯示了不同攻角下升力系數(shù)和阻力系數(shù)的變化曲線(xiàn)。升力系數(shù)隨攻角接近于線(xiàn)性增長(zhǎng)，而阻力系數(shù)隨攻角的增長(zhǎng)規(guī)律呈非線(xiàn)性。

圖6 升力系數(shù)和阻力系數(shù)隨攻角變化曲線(xiàn)Fig.6 Variation curve of lift and drag coefficient with angle of attack

圖7 壓心系數(shù)和升阻比隨攻角變化曲線(xiàn)Fig.7 Variation curve of pressure center coefficient and lift drag ratio with angle of attack

圖7顯示了彈丸壓心系數(shù)和升阻比隨攻角的變化曲線(xiàn)。隨攻角的逐漸增大，彈丸的壓心系數(shù)逐漸變小，說(shuō)明彈丸的靜穩(wěn)定裕度降低，在大攻角下，彈丸的穩(wěn)定性變差，不利于彈丸姿態(tài)收斂。此外，升阻比隨攻角先增大后減小，從圖6也可看出，升力系數(shù)和阻力系數(shù)增長(zhǎng)趨勢(shì)不同，在大攻角下，阻力系數(shù)增長(zhǎng)更快，而升力系數(shù)相對(duì)較慢，導(dǎo)致升阻比增長(zhǎng)逐漸變小，存在極值，針對(duì)本文研究的彈丸，在Ma=4工況下，攻角在12°左右，升阻比達(dá)到最大值。這是因?yàn)樵诖蠊ソ窍?，氣體流動(dòng)分離導(dǎo)致誘導(dǎo)阻力急劇增大，使得總阻力增大，最終導(dǎo)致升阻比達(dá)到極值，這也符合升阻比變化規(guī)律。

2.3 再入工況氣動(dòng)分析

重點(diǎn)分析了再入工況下彈丸在大攻角下的氣動(dòng)特性。圖8顯示了攻角25°下，不同截面的流場(chǎng)壓力分布云圖(圖8(a)和圖8(b)的截面定義與圖5相同)。大攻角下，氣體橫向效應(yīng)比較明顯，彈丸迎風(fēng)面這一側(cè)的激波強(qiáng)度較大，激波后的靜壓急劇增大，而彈丸的背風(fēng)面由于氣流的分離作用使得壓力逐漸減小。相對(duì)于低攻角工況，流場(chǎng)的波系結(jié)構(gòu)對(duì)彈丸的影響比重較大，使得彈丸的氣動(dòng)阻力和升力均增大。

(a) 0°截面(a) 0° section

(b) 45°截面(b) 45° section圖8 不同截面的流場(chǎng)壓力分布云圖(再入工況)Fig.8 Pressure contour in different cross sections (reentry case)

圖9對(duì)比了攻角3°和攻角25°兩種工況下的底部流場(chǎng)結(jié)構(gòu)。在小攻角下，彈丸來(lái)流邊界層未發(fā)生分離現(xiàn)象，底部回流區(qū)相對(duì)對(duì)稱(chēng)，而在大攻角下，由于彈丸背風(fēng)面的邊界層分離與再附作用，使得來(lái)流較為復(fù)雜，回流區(qū)呈現(xiàn)更為復(fù)雜的形態(tài)，存在滑移線(xiàn)，形成剪切層。這些流場(chǎng)形態(tài)均會(huì)影響彈丸的氣動(dòng)特性和飛行穩(wěn)定性。

(a) 3° 攻角(a) 3° angle of attack

(b) 25°攻角(b) 25° angle of attack圖9 彈丸底部流場(chǎng)速度分布云圖Fig.9 Velocity distribution of flow field at the base of projectile

圖10 流向截面的壓力分布云圖(再入工況)Fig.10 Pressure distribution in streamwise section (reentry case)

由于彈丸在大攻角下橫流效應(yīng)明顯，這里詳細(xì)分析了流向截面的壓力分布和流向渦結(jié)構(gòu)分布。圖10顯示了靠近彈丸尾部流向截面的壓力分布云圖，可明顯看出經(jīng)過(guò)激波后的流場(chǎng)壓力突增，類(lèi)似于乘波結(jié)構(gòu)，而在彈丸的背風(fēng)面，由于分離作用彈丸表面壓力較小，使得彈丸上下壓差較大，形成氣動(dòng)升力。圖11對(duì)比了不同攻角下的流向渦結(jié)構(gòu)，重點(diǎn)分析橫流繞過(guò)彈丸表面后，在背風(fēng)面的氣流形態(tài)，可看出，在小攻角下，氣體橫向流速較小，基本貼體流動(dòng)；而在大攻角下，橫向氣流繞過(guò)彈體后形成脫落渦結(jié)構(gòu)及一對(duì)清晰的流向渦，流向渦導(dǎo)致彈丸背風(fēng)面形成低壓回流區(qū)，使得壓力降低。由此可知，大攻角下彈丸的迎風(fēng)面氣動(dòng)流場(chǎng)類(lèi)似于乘波結(jié)構(gòu)，彈丸的背風(fēng)面流場(chǎng)類(lèi)似于低壓脫落渦結(jié)構(gòu)。這兩種流場(chǎng)機(jī)理均對(duì)彈丸的氣動(dòng)特性存在影響。

(a) 3° 攻角(a) 3° angle of attack

(b) 25° 攻角(b) 25° angle of attack圖11 不同攻角流向截面的流向渦分布云圖Fig.11 Streamwise vortex distribution in different streamwise sections

圖12 彈丸表面壓力系數(shù)流向變化曲線(xiàn)Fig.12 Pressure coefficient curve on projectile surface

為進(jìn)一步分析這兩種機(jī)理對(duì)彈丸氣動(dòng)特性的影響效果，分析了彈丸表面的壓力系數(shù)變化規(guī)律，如圖12所示。氣流經(jīng)過(guò)彈丸頭部后，由于激波作用，壓力系數(shù)急劇降低，隨后流經(jīng)迎風(fēng)面的氣流壓力系數(shù)逐漸減小，而背風(fēng)面的氣流壓力系數(shù)變化較小，在靠近尾翼處，由于尾翼的影響，使得壓力系數(shù)產(chǎn)生突變。對(duì)比不同攻角的壓力系數(shù)變化，可看出，隨攻角的增大，迎風(fēng)面的壓力系數(shù)變化較大，而背風(fēng)面的壓力系數(shù)變化相對(duì)較小，這說(shuō)明迎風(fēng)面的波系結(jié)構(gòu)、流場(chǎng)特性對(duì)彈丸的氣動(dòng)特性影響更為顯著。

圖13顯示了再入工況下(Ma=2)彈丸阻力系數(shù)和壓心系數(shù)隨攻角變化曲線(xiàn)。隨著彈丸攻角的增大，阻力系數(shù)迅速增大，尤其是攻角大于15°后，彈丸的阻力系數(shù)增長(zhǎng)較快，這說(shuō)明了大攻角下阻力系數(shù)的非線(xiàn)性特性。同時(shí)，彈丸的壓心系數(shù)也逐漸減小，當(dāng)攻角大于15°后，由于彈丸壓心系數(shù)的前移，使得靜穩(wěn)定裕度降低了4.5%，這說(shuō)明再入工況下，彈丸的再入攻角較大，使得彈丸氣動(dòng)阻力迅速增大，靜穩(wěn)定裕度急劇降低，彈丸的姿態(tài)收斂較慢，從而使彈丸再入段的速度衰減較大。

圖13 阻力系數(shù)和壓心系數(shù)隨攻角變化曲線(xiàn)(再入工況)Fig.13 Drag coefficient and pressure center coefficient change with angle of attack (reentry case)

彈丸的阻力系數(shù)是攻角的偶函數(shù)，略去四階以上微分項(xiàng)，可寫(xiě)為：

(4)

(5)

式中，ΔCD=(CDj-CD(j-1))/(αj-αj-1)，Δα=αj-αj-1。計(jì)算結(jié)果見(jiàn)表3。初步評(píng)估在大攻角工況下，彈丸的阻力二階項(xiàng)是小攻角的2倍，這說(shuō)明在大攻角下，阻力系數(shù)二階項(xiàng)或誘導(dǎo)阻力是彈丸氣動(dòng)阻力的主要組成部分。

彈丸升力系數(shù)是攻角的奇函數(shù)，略去高階項(xiàng)，可得：

(6)

可采用一階差分格式評(píng)估彈丸氣動(dòng)升力的非線(xiàn)性特性。計(jì)算公式為：

(7)

計(jì)算結(jié)果可見(jiàn)表3。在小攻角工況下，ΔCL/Δα的值較小，且基本相同，這說(shuō)明升力系數(shù)隨攻角的線(xiàn)性變化特性；在大攻角下，ΔCL/Δα的值是小攻角的2倍，這說(shuō)明一階線(xiàn)性近似誤差較大，而且隨著攻角的逐漸增大，ΔCL/Δα值逐漸增大，這說(shuō)明大攻角下彈丸升力系數(shù)的非線(xiàn)性特性，高階項(xiàng)的比重逐漸增大，需要引入更高階量才能更為準(zhǔn)確地計(jì)算氣動(dòng)參數(shù)。

表3 氣動(dòng)系數(shù)與攻角的差分量

3 結(jié)論

采用數(shù)值計(jì)算方法研究了超高速?gòu)椡璧臍鈩?dòng)流場(chǎng)特性。利用風(fēng)洞試驗(yàn)數(shù)據(jù)對(duì)比了兩種湍流模型對(duì)超高速?gòu)椡铓鈩?dòng)數(shù)據(jù)的預(yù)測(cè)能力，計(jì)算結(jié)果表明兩種湍流模型在計(jì)算法向力和俯仰力矩系數(shù)等方面的預(yù)測(cè)精度較高，可達(dá)到2%以?xún)?nèi)，而在阻力系數(shù)預(yù)測(cè)方面，精度較差，S-A湍流模型的預(yù)測(cè)精度達(dá)到4.6%，而k-ωSST湍流模型的預(yù)測(cè)精度達(dá)到14%，這說(shuō)明S-A湍流模型在預(yù)測(cè)彈丸外流場(chǎng)特性方面具有較高的精度，這與先前的研究結(jié)果一致。

分析了彈丸的底部和橫向流場(chǎng)分布特性。相對(duì)于小攻角的氣動(dòng)流場(chǎng)，大攻角的彈丸底部流場(chǎng)更為復(fù)雜，回流區(qū)存在滑移線(xiàn)等非定常流場(chǎng)結(jié)構(gòu)。大攻角下的彈丸橫流效應(yīng)更為明顯，使得高速流場(chǎng)在彈丸背風(fēng)面形成對(duì)稱(chēng)的脫落流向渦，這些流向渦導(dǎo)致彈丸背風(fēng)面的壓力降低。在迎風(fēng)面，激波效應(yīng)較為明顯，激波后壓力急劇增大，類(lèi)似于乘波結(jié)構(gòu)，而迎風(fēng)面的壓力變化對(duì)彈丸的氣動(dòng)阻力和升力特性貢獻(xiàn)作用更大。

最后，分析了彈丸氣動(dòng)非線(xiàn)性特性，在大攻角下，彈丸的氣動(dòng)升力和阻力系數(shù)呈現(xiàn)明顯的非線(xiàn)性特性，阻力系數(shù)明顯增大，而且彈丸的靜穩(wěn)定裕度也急劇降低，使得彈丸的收斂特性變差，這是引起彈丸再入段速度衰減的主要原因。