二級輕氣炮內(nèi)彈道過程數(shù)學(xué)建模及數(shù)值仿真

莊 宇，陸 欣

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

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

二級輕氣炮內(nèi)彈道過程數(shù)學(xué)建模及數(shù)值仿真

莊 宇，陸 欣

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

根據(jù)二級輕氣炮的發(fā)射特點(diǎn)，采用經(jīng)典內(nèi)彈道模型描述藥室里火藥燃燒狀況和活塞運(yùn)動，同時采用一維非定?？蓧嚎s流動模型描述輕氣室里的氣體流動狀態(tài)和彈丸運(yùn)動，并通過活塞的運(yùn)動狀態(tài)將兩者耦合，從而建立起二級輕氣炮發(fā)射過程的內(nèi)彈道數(shù)學(xué)模型；以某30 mm/120 mm輕氣炮為基本計算模型，運(yùn)用四階Runge-Kutta法求解藥室方程，運(yùn)用二階MacCormack格式求解輕氣室方程，通過兩部分的交替計算，實(shí)現(xiàn)二級輕氣炮內(nèi)彈道過程的數(shù)值仿真，為二級輕氣炮的參數(shù)設(shè)計和發(fā)射性能的提高提供了理論依據(jù)。

二級輕氣炮；數(shù)學(xué)建模；數(shù)值仿真；發(fā)射性能

二戰(zhàn)后，人們對于超高速發(fā)射的研究越來越感興趣。對傳統(tǒng)火炮而言，由于火藥氣體分子量較大，滯止聲速較小，彈丸獲得的最大速度受到限制，很難達(dá)到3 km/s以上。輕氣炮是一種利用熱輕質(zhì)氣體(如氫和氦)膨脹做功的方式推動彈丸，從而增大氣體逃逸速度，減小氣體聲慣性，使之獲得極高速度的發(fā)射裝置。在1946年世界上第一門輕氣炮誕生后的幾十年里，這種高溫低分子量氣體的發(fā)射系統(tǒng)在實(shí)驗室中進(jìn)行的高速發(fā)射研究已經(jīng)取得了輝煌的成果[1]。近年來，隨著輕氣炮發(fā)射技術(shù)不斷發(fā)展，已經(jīng)可以發(fā)射各種形狀的彈丸，而且彈丸的質(zhì)量、尺寸和材料都有較寬的范圍。同時，輕氣炮應(yīng)用領(lǐng)域也不斷拓寬，已經(jīng)從軍工和航天擴(kuò)展到了氣動力、超高速碰撞、材料力學(xué)性能等方面[2-4]。

目前，國內(nèi)外最為常見的是二級輕氣炮。由于二級輕氣炮涉及火藥氣體、活塞、輕質(zhì)氣體以及彈丸這四者之間的相互作用，發(fā)射過程較為復(fù)雜。而結(jié)構(gòu)尺寸與裝填條件的相互影響和共同作用也使得每門二級輕氣炮都存與自身最匹配的發(fā)射參數(shù)[5]。多年來，國內(nèi)外各研究單位對此進(jìn)行了大量的試驗和研究。美國的Douglas公司、GM通用汽車公司、NOL海軍軍械研究室等單位開展的超高速撞擊試驗為美國的軍事及民用方面的發(fā)展做出了巨大貢獻(xiàn)。其中通用汽車公司(GM)甚至將0.053 3 g的彈丸加速到10.80 km/s的驚人速度[6]。但如果僅僅通過反復(fù)試驗來確定每門炮的最佳發(fā)射參數(shù)，顯然事倍功半。因此，從減少實(shí)驗開支等看，輔以必要的二級輕氣炮的數(shù)值仿真研究顯得尤為重要。

為此，針對國內(nèi)外最常見、應(yīng)用最為廣泛的一種二級輕氣炮，運(yùn)用經(jīng)典內(nèi)彈道理論和氣體動力學(xué)模型，建立二級輕氣炮內(nèi)彈道過程的數(shù)學(xué)模型，并提出一種行之有效的求解方法。

1 二級輕氣炮數(shù)學(xué)模型

圖1是一個常見的二級輕氣炮的結(jié)構(gòu)示意圖。它的第一級工質(zhì)是火藥氣體，第二級工質(zhì)是輕質(zhì)氣體，一二級之間通過一個自由活塞進(jìn)行能量傳遞。其發(fā)射過程始于固體火藥的點(diǎn)火和燃燒，火藥燃燒產(chǎn)生的氣體驅(qū)動活塞運(yùn)動，從而急劇壓縮泵管中的輕質(zhì)氣體達(dá)到高壓高溫，直到用于隔離彈丸和輕氣的膜片破裂，使得彈丸從靜止開始加速穿過發(fā)射管[7]。需要注意的是，彈丸運(yùn)動后，活塞會繼續(xù)壓縮輕質(zhì)氣體，直至進(jìn)入高壓段開始變形后，由于受到錐形壁摩擦和擠壓以及前端氣體壓力所形成的阻滯力而不斷減速至運(yùn)動停止(如二級輕氣炮參數(shù)匹配得當(dāng)，可使活塞恰好停在錐段附近)。

圖1 二級輕氣炮結(jié)構(gòu)示意圖

考慮上述發(fā)射過程，藥室和泵管的求解方程應(yīng)該單獨(dú)建立并通過活塞的運(yùn)動耦合，從而建立起二級輕氣炮的數(shù)學(xué)模型。本文用一組常微分方程，即經(jīng)典內(nèi)彈道模型描述藥室內(nèi)火藥燃燒和活塞的運(yùn)動規(guī)律，同時用一維非定常可壓縮流體模型描述泵管內(nèi)輕氣的壓縮過程和彈丸的運(yùn)動規(guī)律。鑒于這樣的分析，提出如下假設(shè)：火藥燃燒服從幾何燃燒定律；輕氣視為恒定比熱容的無黏氣體，且為一維非定常變截面流動；活塞運(yùn)動終了時，在緩沖限制器的作用下速度為零；熱損失通過火藥力f和絕熱指數(shù)k修正，不做直接計算；其他假設(shè)等同于經(jīng)典內(nèi)彈道假設(shè)。

1.1 藥室里的基本方程

1) 火藥顆粒形狀函數(shù)：

(1)

式(1)中，z和ψ分別是火藥顆粒的相對燃燒厚度和相對燃燒體積，χ、λ、μ、χs、λs和Zk都是火藥的形狀特征量。

2) 火藥顆粒燃速方程：

(2)

式(2)中，ps是藥室的平均壓力，u1和n分別是火藥的燃速系數(shù)和燃速指數(shù)，e1的值是火藥顆粒初始厚度的一半。

3) 活塞運(yùn)動方程：

(3)

式(3)中，mp、vp和φ分別是活塞的質(zhì)量、速度和次要功系數(shù)，ppf是作用在活塞面上的氦氣壓力，S是活塞的橫截面積。

4) 活塞速度方程：

(4)

式(4)中，lp是活塞的行程。

5) 能量守恒方程：

(5)

式(5)中，f和ω分別是火藥力和火藥的裝藥量，θ是與火藥氣體相關(guān)的常數(shù)，lψ是藥室自由容積縮徑長。

1.2 輕氣室中的控制方程

1) 連續(xù)性方程：

(6)

式(6)中，ρ和u是分別是泵管里氦氣的密度和速度，A是泵管的橫截面積。

2) 動量方程：

(7)

式(7)中，P是泵管中氦氣的壓力，D和L分別是輕氣室的直徑和長度，cf和k分別是Darcy阻力系數(shù)和壓力損失因子[7]。其中，cf可用以下關(guān)系式表示：

(8)

式(8)中，Re為雷諾數(shù)，Λ是對于高速流動十分重要的可壓縮修正因子，可用半經(jīng)驗公式表示

(9)

其中，Ma為當(dāng)?shù)亓鲃拥鸟R赫數(shù)，Ω為恢復(fù)因子。

除非管道面積變化很快，或在流動中隔板阻礙相當(dāng)嚴(yán)重，壓力損失才比較明顯。通常情況下，式(7)中的k都比較小，可忽略不計。

3) 能量方程：

(10)

式(10)中，e是氦氣的內(nèi)能，λh和Nu是氦氣的導(dǎo)熱系數(shù)和努塞爾數(shù)，T和Tw分別是氦氣的溫度和壁溫。

4) 輕氣的狀態(tài)方程：

(11)

式(11)中，γ和α分別是氦氣的比熱容比和余容。

2 數(shù)學(xué)模型的求解

2.1 數(shù)值計算方法

藥室的基本方程是由常微分方程和代數(shù)方程組成，因此可采用4階Runge-Kutta法數(shù)值求解。經(jīng)典Runge-Kutta法與藥室基本方程相容，解無條件收斂和穩(wěn)定[8]，因而時間步長Δt僅由輕氣室所采用的差分格式的穩(wěn)定性來確定。

輕氣室的控制方程是由一組偏微分方程和一個代數(shù)方程組成，一般采用有限差分法進(jìn)行數(shù)值求解。通過對幾種常見差分格式的比較，查閱相關(guān)資料，本文決定選用2階MacCormack 預(yù)估-校正格式來求解二級輕氣炮輕氣室內(nèi)的一維非定?？蓧嚎s流體模型。此差分格式具有時間上的一階精度和空間上的二階精度，整體計算精度較高。

2.2 邊界條件處理

通過活塞的運(yùn)動，將藥室里的基本方程和輕氣室中的控制方程的求解耦合到一起。一方面藥室方程組的計算結(jié)果(活塞速度vp和行程lp)可為輕氣室提供左邊界條件，另一方面輕氣室方程組的計算結(jié)果又為活塞的運(yùn)動提供阻力，二者交替進(jìn)行，直至整個射擊過程完成為止，從而實(shí)現(xiàn)二級輕氣炮發(fā)射過程的數(shù)值仿真。

計算過程中，輕氣室右邊界速度和位移分別由彈丸的速度和行程決定：

(12)

式(12)中，S2、mj、φ2分別為發(fā)射管的橫截面積、彈丸質(zhì)量和彈丸阻力系數(shù)，Pd為彈底壓力，取輕氣室右邊界壓力。

輕氣室密度、壓力、溫度等其他運(yùn)動邊界條件可直接由計算域內(nèi)鄰近節(jié)點(diǎn)的相應(yīng)參數(shù)經(jīng)過插值得到。而輕氣室的固壁邊界條件則需采用第二網(wǎng)格系中的反射法獲得，具體法則如下[9-10]：

(13)

2.3 網(wǎng)格處理

2.3.1 活塞頂端網(wǎng)格處理(輕氣室左端)

在某個時刻活塞運(yùn)動到某一位置，以活塞所處的的第k個網(wǎng)格的步長作為變步長進(jìn)行計算。在下一個時刻活塞運(yùn)動到下一個位置，此時與前面一個網(wǎng)格點(diǎn)(即第(k+1)個網(wǎng)格點(diǎn))的距離為Δx1。若Δx1>0.5Δx,則活塞頂端仍與第k個網(wǎng)格聯(lián)系，并以此新的網(wǎng)格距離Δx1作為新步長進(jìn)行計算；若Δx1<0.5Δx，則將第k個網(wǎng)格與第(k+1)個網(wǎng)格合并，同時把活塞頂端各參數(shù)賦予給新合并的(k+1)網(wǎng)格。新網(wǎng)格步長為(Δx1+Δx)，并且下一時刻的計算就從第(k+1)個網(wǎng)格開始。隨著活塞的不斷向前運(yùn)動，左端的網(wǎng)格不斷減少，直到活塞運(yùn)動到泵管高壓段附近，停止運(yùn)動。

2.3.2 彈丸底部網(wǎng)格處理(輕氣室右端)

當(dāng)彈丸開始運(yùn)動后，輕氣室右邊界由固壁邊界變?yōu)檫\(yùn)動邊界，此時計算網(wǎng)格的數(shù)量要增加。

假設(shè)在某一時刻彈丸運(yùn)動到某一位置處，以彈底與后面第一個格點(diǎn)k的距離作為變步長進(jìn)行計算。在下一個時刻彈丸運(yùn)動到下一個位置，此時與前面第k個網(wǎng)格點(diǎn)的距離為Δx2。若Δx2<1.5Δx，則彈底仍與第k個網(wǎng)格聯(lián)系，并以此新的網(wǎng)格距離Δx2作為新步長進(jìn)行計算；若Δx2>1.5Δx，則分裂出一個新網(wǎng)格(k+1)，并讓此網(wǎng)格與彈底發(fā)生聯(lián)系，此新網(wǎng)格以(Δx2-Δx)作為新步長進(jìn)行計算。

相對于輕氣室左端的網(wǎng)格處理，增加了對于新格點(diǎn)各流動參量的求解。對于新格點(diǎn)的處理，可以采用前后相鄰格點(diǎn)線性插值的方法，來求得新網(wǎng)格點(diǎn)氣體的各個狀態(tài)參量。此時，所有網(wǎng)格點(diǎn)的各狀態(tài)參量就可以求得，從而為下一個時間層的計算做準(zhǔn)備。

2.4 截面濾波

特別要指出，在計算輕氣室各參數(shù)過程中，由于差分格式的色散性，使得波頭可能產(chǎn)生振蕩，且隨著計算的進(jìn)行，振蕩積累到一定程度可能使得數(shù)值計算無法繼續(xù)進(jìn)行，從而造成停機(jī)現(xiàn)象。另外，截面的變化也會使面積的導(dǎo)數(shù)產(chǎn)生間斷。諸如此類問題，通常需要采用濾波處理。常用的濾波方法是Shumann濾波(主要用于截面變化處)，采用的是加權(quán)平均方法，即：

(14)

3 數(shù)值仿真結(jié)果及其分析

在此基礎(chǔ)上，對一個30 mm/120 mm二級輕氣炮進(jìn)行數(shù)值仿真，從而檢驗上述數(shù)學(xué)模型及其計算方法的可行性。30 mm/120 mm二級輕氣炮的主要參數(shù)如表1所示。

表1 30 mm/120 mm二級輕氣炮主要參數(shù)

數(shù)值仿真結(jié)果如圖2～圖6所示。結(jié)合圖2、圖3和圖4可知，當(dāng)藥室壓力達(dá)到活塞起動壓力后，活塞開始運(yùn)動，并不斷壓縮輕質(zhì)氣體。隨著火藥的不斷燃燒，藥室內(nèi)的壓力開始逐漸增大。在大約4.6 ms時，藥室壓力增加到最大值310 MPa。隨后，壓力開始下降。大約在6 ms時，火藥燃燒結(jié)束。此后，隨著活塞的運(yùn)動，火藥氣體空間增加，壓力持續(xù)下降。

圖2 藥室壓力隨時間的變化

圖3 火藥相對已燃百分比隨時間的變化

結(jié)合圖4和圖5可知，在藥室壓力作用下，活塞速度不斷增大至最大速度728 m/s。此后，隨著輕氣室壓力的增大，活塞速度開始減小。當(dāng)活塞運(yùn)動到高壓錐段附近時，由于錐壁阻力的影響，速度迅速降為零，最終停止運(yùn)動。

圖4 活塞和彈丸的速度隨時間的變化

圖5 活塞速度隨行程的變化

結(jié)合圖4、圖6和圖7可以看出，約在22.5 ms時，彈丸開始啟動，彈底壓力開始迅速上升。在維持一段較為平穩(wěn)的壓力過載后，約在24.8 ms，第一個壓力波到達(dá)彈底，彈底壓力達(dá)到第一個最大值347 MPa。之后，由于輕氣室流場中各種壓力波的運(yùn)動和疊加，使得彈底壓力變化劇烈。隨著壓力波在活塞和彈丸之間不斷傳遞反射，彈丸不斷加速運(yùn)動，直至在出炮口處彈丸速度達(dá)到5 906 m/s?？梢灶A(yù)見，如果將燃燒過程、活塞質(zhì)量、輕氣室結(jié)構(gòu)以及彈丸質(zhì)量之間關(guān)系匹配得當(dāng)，使得彈丸在出炮口有一個大幅值的壓力波位于彈底，這對增加彈丸的速度非常有利。

圖6 彈丸速度隨行程的變化

圖7 彈底壓力隨時間的變化

從圖8可以看出，在初始時刻，活塞不斷壓縮輕質(zhì)氣體，導(dǎo)致輕氣室左端壓強(qiáng)比右端壓強(qiáng)稍大，因為此時壓縮波尚未傳到右端。當(dāng)接近高壓段時，輕質(zhì)氣體分布在一個相對較小的空間里，壓強(qiáng)急劇增大，且在空間上的分布出現(xiàn)極大跳動，產(chǎn)生了激波。隨后輕質(zhì)氣體絕熱膨脹，推動彈丸向前運(yùn)動，活塞停止，氣體壓強(qiáng)逐漸減小。

圖8 不同時刻輕氣室壓強(qiáng)沿軸向分布

結(jié)合圖2～圖8，將其與文獻(xiàn)[11]中提到的中國工程物理研究院某30 mm口徑的二級輕氣炮的實(shí)驗數(shù)據(jù)進(jìn)行比較后發(fā)現(xiàn)，該仿真結(jié)果較為準(zhǔn)確地反映了二級輕氣炮發(fā)射過程的實(shí)際狀況。數(shù)值仿真的彈丸出炮口速度為5.91 km/s，實(shí)驗值為5.56 km/s，誤差為6%，在允許范圍內(nèi)。這也在一定程度上驗證了所建立的二級輕氣炮的數(shù)學(xué)模型和所采用的數(shù)值求解方法。

4 結(jié)論

本文運(yùn)用經(jīng)典內(nèi)彈道理論和氣體動力學(xué)模型建立了二級輕氣炮內(nèi)彈道過程的數(shù)學(xué)模型，給出了一種運(yùn)用四階Runge-Kutta法和二階MacCormack格式分別交替求解藥室方程和輕氣室方程的數(shù)值計算方法，從而實(shí)現(xiàn)了二級輕氣炮的數(shù)值仿真，較為準(zhǔn)確地刻畫了它的發(fā)射過程。由此可見，通過數(shù)值仿真來輔助二級輕氣炮的實(shí)驗研究是可行的，可以為二級輕氣炮的試驗裝填參數(shù)優(yōu)化和發(fā)射性能的提高提供一定的理論依據(jù)。

[1] BERGGREN R E，REYNOLDS R M.The Light-Gas-Gun Model Launcher[R].AD 713915，1970.

[2] 黃潔,梁世昌,李海燕,等.二級輕氣炮發(fā)射過程內(nèi)彈道數(shù)值計算研究[J].空氣動力學(xué)學(xué)報，2013，31(5):657-661.

[3] 柳森，李毅，黃潔，等．彈丸超高速撞擊單層和多層板結(jié)構(gòu)的碎片特征研究[J]．宇航學(xué)報，2010，31(6)：1672-1677.

[4] ZENG X J,YU Z F,BU S Q,et al.Research on the RCS of Hypervelocity Model and Its Plasma Sheath[J].ACTA Aerodynamica Sinica,2010,28(6):645-649.

[5] 管小榮，徐誠.二級輕氣炮發(fā)射過程數(shù)學(xué)模型和計算方法[J].南京理工大學(xué)學(xué)報，2007，31(1)：22-26.

[6] 王東方，肖偉科，龐寶君.NASA二級輕氣炮設(shè)備簡介[J].實(shí)驗流體力學(xué),2014,28(4):99-104.

[7] 金志明，翁春生.高等內(nèi)彈道學(xué)[M].北京：高等教育出版社，2003.

[8] 鄭慧嬈，陳紹林，莫忠息，等.數(shù)值計算方法[M].武漢：武漢大學(xué)出版社，2002.

[9] 翁春生，王浩.計算內(nèi)彈道學(xué)[M].北京：國防工業(yè)出版社，2006.

[10]黃玲智.具有脈沖擴(kuò)散和擴(kuò)散時滯的單種群模型研究[J].重慶工商大學(xué)學(xué)報(自然科學(xué)版),2015,32(9):5-9.

[11]馮福全.二級輕氣炮內(nèi)彈道仿真研究[D].南京：南京理工大學(xué)，2008.

[12]黃滔，張相炎，劉寧.燃燒輕氣炮點(diǎn)火對內(nèi)彈道性能影響的仿真分析[J].兵工自動化，2014(12):25-28.

(責(zé)任編輯周江川)

Mathematical Modeling and Numerical Simulation for Interior Ballistics Process of Two-Stage Light-Gas Gun

ZHUANG Yu, LU Xin

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

According to the launch features of two-stage light gas gun, classical interior ballistic model was used to describe combustion of powder in chamber and piston motion, and one dimensional unsteady compressible fluid model was used to describe flowing of gas in light-gas chamber and projectile motion. The two groups of equations were linked through piston motion to develop a mathematical model for interior ballistics launch process of two-stage light-gas gun. A 30 mm/120 mm light-gas gun was taken for the primary computing model.The equations in powder chamber were solved numerically with the fourth-order Runge-Kutta method, and the control equations in light-gas chamber were solved numerically with the second-order MacCormack scheme.The numerical simulation for interior ballistics process of two-stage light-gas gun was accomplished by alternation between these two computations. The simulation results provide theoretical references for parameter design and improvement of launching performance of two-stage light-gas gun.

two-stage light-gas gun; mathematical modeling; numerical simulation; launching performance

2016-07-27；

江蘇省前瞻性聯(lián)合創(chuàng)新資金(BY2014004-11)

莊宇(1992—)，男，碩士研究生，主要從事新型發(fā)射理論與控制技術(shù)研究。

10.11809/scbgxb2016.12.018

莊宇，陸欣.二級輕氣炮內(nèi)彈道過程數(shù)學(xué)建模及數(shù)值仿真[J].兵器裝備工程學(xué)報，2016(12):75-79.

format:ZHUANG Yu, LU Xin.Mathematical Modeling and Numerical Simulation for Interior Ballistics Process of Two-Stage Light-Gas Gun[J].Journal of Ordnance Equipment Engineering,2016(12):75-79.

TJ012.1

A

2096-2304(2016)12-0075-05

修回日期：2016-08-30