超聲速彈丸高空流場數(shù)值模擬與稀薄氣體效應研究

王 勇，王學德，許國龍

(南京理工大學 能源與動力工程學院， 南京 210094)

?

【裝備理論與裝備技術(shù)】

超聲速彈丸高空流場數(shù)值模擬與稀薄氣體效應研究

王 勇，王學德，許國龍

(南京理工大學 能源與動力工程學院， 南京 210094)

利用直角坐標網(wǎng)格技術(shù)研究高超聲速稀薄氣體流動的DSMC模擬方法，編制了高空彈丸外流場數(shù)值模擬通用計算程序，將本方法的計算結(jié)果與有關(guān)文獻中公布的理論模型數(shù)據(jù)進行了對比，驗證本方法的正確性。模擬了稀薄流下不同海拔高度尖拱型彈頭和其他外形彈頭的在超聲速來流中的流場，并推廣至不同外形彈頭氣動受力情況的計算，研究其稀薄效應。分析比較發(fā)現(xiàn)了隨著海拔高度的升高，稀薄效應使彈丸的激波層厚度增加，外部流場壓縮性減弱，表面的熱流密度有所降低。相比于尖拱型彈頭，鈍頭彈丸在駐點處所受的壓力較小，尾部較強，符合高超聲速飛行器鈍頭彈體的受力規(guī)律。

超聲速；數(shù)值模擬；稀薄流；稀薄效應

現(xiàn)代戰(zhàn)爭迫使彈箭的飛行高度越來越高。對于高空飛行的飛行器而言，氣體稀薄效應越發(fā)顯著。對于高空彈箭而言，由于空氣稀薄，經(jīng)典氣體動力學的連續(xù)介質(zhì)假設不成立，無法用傳統(tǒng)的Navier-Stokes方程求解，必須采用一種新的方法。DSMC方法使用大量的模擬分子代表真實氣體并模擬真實分子的運動和碰撞，通過計算機來直接模擬流場本身，摒除了數(shù)學求解的困難，并且將模擬分子的運動信息進行存儲，目前已經(jīng)相當成功地應用到了各個行業(yè)的諸多方面。如美國航天飛機升阻比的計算、長征-4B火箭推進劑在軌排放擾動力矩的計算、美國通訊衰減飛行試驗中的電子密度分布計算等等。有效地分析氣動力參數(shù)對于高空超音速彈箭至關(guān)重要，是其安全飛行的基礎(chǔ)，準確模擬高速稀薄氣流如今只能采用DSMC算法。

國外對于高空稀薄流的研究更早也更加深入,如G.A.Bird[1]關(guān)于DSMC算法在航天飛行器上的應用。國內(nèi)的DSMC研究起步較晚但是發(fā)展速度非常快并且基本成型，李學東、黃海龍等人[2-4]采用非結(jié)構(gòu)化三角網(wǎng)格，利用VHS分子模型、CLL反射模型對高超音速氣體繞圓柱流動問題進行了數(shù)值模擬，蔡國飆[5]、葉友達[6]則分別論述了DSMC的時間步長和高空飛行器的氣動特性，但對于高空超音速彈丸的稀薄流場研究較少。

本文采用直角坐標網(wǎng)格技術(shù)研究高超聲速稀薄氣體流動的DSMC模擬方法，選用由Thacker、Gonzalez和Putland在1980年提出的傳統(tǒng)的四叉樹方法生成直角坐標網(wǎng)格,分子碰撞模型采取VHS變徑硬球模型。模擬計算了二維稀薄氣流下60 km和80 km高空的超音速平板繞流，并與公布的文獻數(shù)據(jù)進行對比驗證，計算了不同海拔高度下尖拱型彈丸以及其他外形彈頭在超聲速稀薄來流的流場，發(fā)現(xiàn)隨著海拔高度的升高，由于高空的稀薄效應使外部流場激波層厚度增加，彈體表面熱流度沿軸線降低等現(xiàn)象。尖拱型彈丸頭部由于連續(xù)外凸的表面，其形成的激波尾部明顯寬于拋物線型彈丸，表現(xiàn)為對彈頭尾部流域的影響更為強烈。

1 數(shù)值模擬方法

Bird發(fā)展的DSMC(Direct Simulation Monte-Carlo)方法[7]直接從流動的物理模擬出發(fā)，根據(jù)分子動力學方法，采取概率論判斷分子間是否發(fā)生碰撞。DSMC方法用大量的模擬分子來模擬真實的氣體，將模擬分子的空間坐標、速度、內(nèi)能等數(shù)據(jù)進行儲存，并隨著分子的運動、碰撞和與邊界的作用而變化，模擬所用的時間應該與實際情況的時間一致。整個過程是非定常的，經(jīng)過足夠長時間的模擬后，可以認為流動狀態(tài)達到了定常流。DSMC對于求解過渡領(lǐng)域流場的總體效應方面和精細結(jié)構(gòu)方面均可以得到與實驗準確相符的結(jié)果，該方法對于復雜外形的航天飛機在過渡區(qū)所得結(jié)果與飛行測量數(shù)據(jù)的一致性很好[1]。

DSMC計算方法主要是模擬分子的微觀運動和分子間的碰撞，而分子的運動和碰撞一般是耦合在一起的，常用的數(shù)值模擬方法很難將二者解耦，因此，對于使用DSMC模擬稀薄流場時分子的運動和碰撞需要特別處理。在一個時間步長Δt內(nèi)，先讓各個分子按照各自的速度運動一段距離，再按照碰撞模型以真實的碰撞頻率選出碰撞對，分配碰撞后的分子速度和內(nèi)能，并確定是否發(fā)生化學反應。Δt的取值應該小于平均碰撞時間，通過空間網(wǎng)格計算碰撞并對氣體分子量進行統(tǒng)計得到宏觀量，網(wǎng)格的尺寸與分子平均自由程相比應該足夠小。碰撞只會在同一網(wǎng)格間發(fā)生，從而實現(xiàn)分子運動和碰撞的解耦。模擬計算初始時需給定計算域內(nèi)模擬分子的位置和速度，而后在每一個時間步長內(nèi)，處理分子運動和碰撞，分子的運動超出所設定的邊界，考慮其與邊界的相互作用，引入相應的反射模型，分子間的碰撞處理方式主要是分子碰撞的取樣，即恰當?shù)倪x擇碰撞對并實現(xiàn)一定數(shù)目的碰撞，這是使Δt時間內(nèi)碰撞與運動匹配，使模擬與真實流動過程一致的關(guān)鍵。

2 直角坐標網(wǎng)格生成

DSMC方法中的網(wǎng)格選取十分關(guān)鍵，關(guān)系著分子碰撞對的選取和流場的正確取樣。對于同樣的流場，如果用不同類型的網(wǎng)格進行模擬，結(jié)果也可能不同。所要求網(wǎng)格的尺度要小于當?shù)仄骄杂沙蹋⒈ＷC網(wǎng)格劃分足夠精細，使網(wǎng)格內(nèi)的參數(shù)均勻的假設充分成立。在本研究中，考慮到物面形狀、計算耗時以及計算精度等原因，選用傳統(tǒng)的四叉樹方法生成網(wǎng)格。

通常在流場的模擬中物面邊界的形狀千差萬別，網(wǎng)格的屬性也不相同，除去流場網(wǎng)格和固體網(wǎng)格還有物面網(wǎng)格，本研究采用限制盒技術(shù)，去掉了外部大部分流場網(wǎng)格，然后使用Painting algorithms[8]方法篩選剩余網(wǎng)格。首先在物面外部選擇一個網(wǎng)格，然后在這個網(wǎng)格所有的相鄰網(wǎng)格單元中進行搜索。如果被搜索到的這個網(wǎng)格與物面不相交，那么將這個網(wǎng)格進行標識，放入“存儲庫”中。如果這個網(wǎng)格和物面相交，那么不做標識。搜索完這個網(wǎng)格相鄰的所有網(wǎng)格后，在“儲存庫”中取出一個網(wǎng)格再進行上述運算，直到庫中沒有網(wǎng)格為止。

對于物面邊界而言，其與網(wǎng)格相交被切割成很多個離散的線段。當分子和物面發(fā)生碰撞時，可以認為一條完整的物面邊界所受到的分子碰撞等于各個離散線段受到的分子碰撞之和。而本研究中所討論的離散線段應起始于網(wǎng)格的一邊，終止與該網(wǎng)格的另一條邊。因此，首先將一整條物面邊界分成若干離散線段，然后統(tǒng)計離散線段的起始位置和終止位置，并對每條線段編號。在與網(wǎng)格相交時，會出現(xiàn)線段長度較小的情況，如圖1所示。這種長度較小的情況不能忽略。在程序的設計中，規(guī)定了可識別的最小長度k。凡是低于長度k的都屬于長度較小的情況。這較短的線段應該添加進它之前的線段，然后再對新線段進行k值判斷，直到所有的線段長度都符合要求。接下來模擬每條離散線段的分子碰撞情況，最后將這些數(shù)據(jù)進行匯總，得到這條完整的物面邊界的模擬結(jié)果。

圖1 線段的合并

3 高空稀薄流場數(shù)值模擬的驗證

3.1 計算區(qū)域和網(wǎng)格

選取的計算網(wǎng)格模型為平板，如圖2所示，計算區(qū)域為(0.0 m，1.0 m)×(0.0 m，0.6 m)，網(wǎng)格類型為結(jié)構(gòu)網(wǎng)格。網(wǎng)格單元數(shù)為6 000，網(wǎng)格節(jié)點數(shù)為6 197。流場的下邊界為對稱邊界，左側(cè)和上邊界為自由來流，流場右側(cè)設為真空。

圖2 計算網(wǎng)格

3.2 計算條件

來流為空氣，其中氮氣和氧氣的組分分別為0.79和0.21，其他組分的氣體因其所占比例較少而忽略不計。來流溫度為300 K，來流速度為1 412.5 m/s。氮氣分子的質(zhì)量為4.65×10-26kg黏性系數(shù)的溫度指數(shù)為0.74，氧氣分子的質(zhì)量為5.312×10-26kg，黏性系數(shù)的溫度指數(shù)為0.77。隨著高度的變化，分子數(shù)密度也不同，在60 km高度上，分子數(shù)密度5.995×1021m3。在80 km高度，分子數(shù)密度變?yōu)?.264 2×1020m3。采用不同的分子數(shù)密度模擬了不同高度環(huán)境下的平板繞流流場情況。

3.3 計算結(jié)果分析

圖3是兩種高度下平板繞流的密度分布云圖。超音速來流橫掠平板時，在前緣點形成了斜激波，并且在前緣點附近激波與邊界層相互作用相互影響并且逐漸合并。在距前緣點十幾個分子自由程的位置，二者開始發(fā)生分離。平板上方的流場受到了平板前緣擾動的影響，而超聲速來流在流經(jīng)平板后，在流場的右側(cè)邊界也就是真空處，密度迅速下降。對比60 km和80 km兩種工況，隨著高度的增加，流場密度也逐漸減小，可以看出氣體稀薄程度對流場的影響。圖4是不同高度下流場溫度分布云圖，因為在初始邊界條件設定時，設壁面為等溫邊界，這就使得靠近壁面區(qū)域的溫度最高，不但高于激波后的溫度，也高于壁面溫度。前緣產(chǎn)生的擾動影響到了上游溫度分布，可以看出其溫度逐漸升高。圖5顯示了不同高度下的流場壓力分布，經(jīng)過斜激波后壓力增加，在前緣點附近壓力增加到最大。沿平板表面方向氣體流動速度不斷加快，壓力也隨之減小，因為流場右側(cè)設為真空邊界，壓力在這里迅速降低。高度的增加也會導致流場各處壓力減小。這是因為在微觀上，壓力可以解釋為單位時間內(nèi)無數(shù)個氣體分子熱運動對物面的撞擊作用力。在相同的來流速度下，空氣越稀薄，與壁面發(fā)生碰撞的分子數(shù)越少，壓力越小。

圖3 不同高度平板繞流的密度分布云圖

圖4 不同高度平板繞流的溫度分布云圖

圖5 不同高度平板繞流的壓力分布云圖

圖6是60 km高度下平板表面系數(shù)的計算數(shù)值圖，其中壓力系數(shù)、摩擦因數(shù)和熱流通量都和文獻[7]進行了對比，可以看出在平板中部，熱流量較實驗值偏小，這是由于壁面設定的等溫條件，而實際情況該部位溫度由于分子摩擦溫度會升高。尾部由于平板右側(cè)的真空條件，使得模擬數(shù)據(jù)有小幅波動，但總體的變化趨勢基本一致，說明了此方法可行。

圖6 平板表面各計算量與文獻數(shù)據(jù)的比較

4 超聲速彈丸高空流場數(shù)值模擬

4.1 尖拱型母線彈頭

彈丸外形如圖7所示，彈頭部分為相切尖拱型。計算網(wǎng)格區(qū)域為(0.0 m,1.0 m)×(0.0 m,0.6 m)，計算區(qū)域的左邊界和上邊界為自由來流邊界，右面邊界為真空邊界，下邊界是對稱邊界。流場共2 430個網(wǎng)格，2 576個節(jié)點。計算條件來流為空氣，氮氣氧氣分別占79%和21%，其他氣體所占組分較小不予考慮。來流溫度300 K，氣流速度為1 412.5 m/s，高度不同，來流的分子數(shù)密度也不同。

圖7 彈頭外形

4.2 計算結(jié)果分析

圖8給出了高空環(huán)境下的流場密度云圖，可以看出，在高空超音速的工況下，在彈頭的前部形成了十分顯著的斜激波。在彈頭前部駐點附近密度較大，氣流經(jīng)過激波后密度逐漸減小。在靠近右側(cè)邊界，距離彈體附面層較近的地方，密度迅速減小。并且對比60 km和80 km下的兩種工況可以看出二者的激波形狀略有差異。60 km工況下產(chǎn)生的激波更加貼體，同時受高空氣體稀薄效應影響，80 km工況下流場的各點密度較60 km時小。圖9給出了不同高度下流場溫度分布云圖，駐點處的溫度依然最大，過駐點處，彈體表面承受著主要的氣動熱，這一部位的氣動燒蝕也最為顯著。沿彈體的軸線方向，溫度逐漸降低。圖10給出了不同高度下的流場速度分布云圖，可以清楚看到，經(jīng)過斜激波后，速度逐漸減小，在彈體表面處，氣體有一定的滑移速度。沿彈體表面速度逐漸減小直至亞音速流動，但是不會減小到零。這是因為在高空稀薄氣體環(huán)境下，氣體分子之間的平均自由程要大得多，在氣體分子同彈體表面發(fā)生碰撞以后，相當一部分的分子反射出去后沒有發(fā)生二次碰撞，這就表現(xiàn)為圖中所示的氣流滑移現(xiàn)象。圖11給出了不同高度表面分布量計算結(jié)果，壓力系數(shù)在頭部駐點處最大，沿彈丸軸線逐漸下降，最后趨于平緩，摩擦因數(shù)和熱流通量變化也大致相同。

圖8 不同高度下相切尖拱型彈頭流場密度分布云圖

圖9 不同高度下相切尖拱型彈頭流場溫度分布云圖

圖10 不同高度下相切尖拱型彈頭流場速度分布云圖

圖11 相切尖拱型彈頭不同高度表面分布量計算結(jié)果

4.3 其他外形彈頭

采用DSMC算法對80 km高度下的相割尖拱型彈頭和拋物線型彈頭進行了模擬分析。圖12、圖13、圖14分別給出了不同外形彈頭的密度、溫度和速度分布云圖，可以看出在相同的海拔高度、來流馬赫數(shù)下，彈頭外形的不同直接影響了激波形狀。尖拱型彈頭頭部的激波較為顯著，且對彈體表面的影響遠強于鈍頭彈形，其駐點處附近密度較大，氣流經(jīng)過激波后逐漸減小，拋物線型彈頭由于其鈍頭形狀，只在彈頭前部駐點處承受了很大的壓力。

圖15給出了不同外形彈頭表面分布量計算結(jié)果，可以看出尖拱型彈頭前部所承受的壓力要高于鈍頭的拋物線型，這是由于其尖端的面積較小所造成的，在駐點后逐漸下降并開始低于鈍頭體的表面壓力，與實際情況相符，二種彈頭的壓力系數(shù)在經(jīng)過大約十幾個自由程的長度后均逐漸趨于穩(wěn)定。摩擦因數(shù)和熱流通量也符合這一趨勢。其中摩擦因數(shù)在彈頭尾部出現(xiàn)小幅度的增加，這是因為流場右側(cè)邊界采用了真空邊界條件，這和真實的環(huán)境存在差異，因此彈頭尾部的模擬結(jié)果有偏差。用DSMC程序模擬得到的曲線在彈頭的中部會有微小的上下起伏情況，這是由于DSMC程序固有的統(tǒng)計擾動所致，如果取樣統(tǒng)計很好的話，這種擾動很小。但是這并不妨礙研究，總體的曲線趨勢很明朗。現(xiàn)實中的高超聲速飛行器一般都是鈍頭體，這也與實際情況下彈箭受力受熱主要集中在鈍頭部分的現(xiàn)象吻合。

圖12 不同外形彈頭流場密度力分布云圖

圖13 不同外形彈頭流場溫度分布云圖

圖14 不同外形彈頭流場速度分布云圖

圖15 不同外形彈頭表面分布量計算結(jié)果

5 結(jié)論

本文利用DSMC技術(shù)研究高超聲速稀薄氣體流動模擬了不同高度下平板擾流流場，通過平板表面各計算值與文獻數(shù)據(jù)的比較，驗證了模擬方法和網(wǎng)格結(jié)構(gòu)的有效性。而后對高空超音速彈丸的氣動特性進行了DSMC模擬計算。建立其他外形的彈頭氣動模型并且使用自主編制的DSMC程序計算了不同外形彈頭的氣動力參數(shù)，模擬了高度變化對于稀薄氣體流場的影響，得出以下結(jié)論：

隨著海拔高度的增加，稀薄效應使得激波的厚度明顯增厚，流場的壓縮性減弱，但駐點處的溫度依然最大，該部位的熱流密度也最強，彈體表面的氣體滑移影響減小，這是由于流動變得稀薄，分子間的碰撞頻率明顯降低，能量交換也不充分。

對于不同外形的彈頭，尖拱型彈丸頭部的激波尾部明顯寬于拋物線型彈丸，直接表現(xiàn)為對彈頭尾部流域的影響更為強烈，同時尖拱型彈丸前部的所受的壓力與阻力均強于鈍型彈頭，但尾部的空氣阻力效應稍弱，因而高超聲速飛行器一般都采用鈍頭體。

[1] BIRD G A.Application of the DSMC method to the full shuttle geometry[J].AIAA,1990-1692.

[2] 李學東.稀薄氣體高超聲速流動的非結(jié)構(gòu)DSMC并行化計算[J].科技回顧,2010,28(4):64-67.

[3] 李學東.非結(jié)構(gòu)化網(wǎng)格下二維鈍頭體繞流DSMC數(shù)值模擬[C]//全國沖擊動力學學術(shù)會議論文集.焦作:中國力學學會爆炸力學專業(yè)委員會沖擊動力學專業(yè)組,2009.

[4] 黃海龍.稀薄氣體中高超速圓柱繞流仿真研究[J].哈爾濱商業(yè)大學學報,2013,29(4):483-487.

[5] 蔡國飆.DSMC當?shù)貢r間步長技術(shù)的理論發(fā)展[J].科學中國,2012,10(2):2750-2756.

[6] 葉友達.高空高速飛行器氣動特性研究[J].力學進展,2009,39(4):387-397.

[7] BIRD G A.Molecular Gas Dynamics and the Direct Simulation of Gas Flow[M].Oxford:Clarendon Press,1994.

[8] AFTOSMIS M J.Solution Adaptive Cartesian Grid methods for Aerodynamic flow with complex geometries[C]//28th Computational fluid dynamics.[S.l.]:[s.n.],1997.

(責任編輯 周江川)

Numerical Simulation and Rarefied Gas Effect of Supersonic Warhead in High Flow Field

WANG Yong, WANG Xue-de, XU Guo-long

(School of Energy and Power Engineering, Nanjing University of Science and Technology, Nanjing 210094, China)

A general computer program for numerical simulation of supersonic warhead in the high flow field by using DSMC simulation method with the Cartesian coordinate’s grid technique researching of hypersonic gas flow was worked out. To prove the validity of this method, the article made a contrast of the results with the data of theoretical model released related literature. Then, rarefied flow was simulated and calculated under different altitude, about the point arch warhead and other kinds of warheads in supersonic flow, the objective extending to the aerodynamic force calculation of different warhead shapes was used to study the effect of rarefied atmosphere. Analysis shows that with the elevation of altitude, the thickness of shock layer of the warhead increases because of rarefaction effect, and the compressibility of external flow field recedes, and the heat flux on the surface of warhead decreases. Compared with the pointed arch warhead, the pressure of blunt projectile is small in the stagnation point and is strong in the rear, which is in accord with the force law of blunt body of hypersonic aircrafts.

supersonic velocity; numerical simulation; rarefied flow; rarefaction effect

2016-05-23；

2016-06-20

江蘇省青藍工程資助項目

王勇(1991—)，男，碩士研究生,主要從事稀薄氣體空氣動力學研究。

10.11809/scbgxb2016.10.003

王勇，王學德，許國龍.超聲速彈丸高空流場數(shù)值模擬與稀薄氣體效應研究[J].兵器裝備工程學報，2016(10):13-18.

format:WANG Yong, WANG Xue-de, XU Guo-long.Numerical Simulation and Rarefied Gas Effect of Supersonic Warhead in High Flow Field[J].Journal of Ordnance Equipment Engineering,2016(10):13-18.

TJ011.+2

A

2096-2304(2016)10-0013-06