張煥好，陳志華，姜孝海，韓珺禮

(南京理工大學(xué) 瞬態(tài)物理國家重點實驗室，江蘇 南京210094)

0 引言

彈丸發(fā)射出膛口瞬間，身管內(nèi)高溫、高壓火藥燃?xì)庠谔趴诩眲∨蛎?，形成以欠膨脹射流為主的瞬態(tài)膛口流場［1］。膛口中流出的高溫、高壓、高速氣流對炮身所產(chǎn)生的反作用力，使炮身后坐增大，并對彈丸飛行產(chǎn)生初始擾動，影響射擊精度。

為了改善武器系統(tǒng)后坐性能，可在大威力炮(槍)膛口上加裝膛口制退器，使流出的高溫高壓火藥燃?xì)庠谥仆似鲀?nèi)膨脹，將內(nèi)能轉(zhuǎn)變?yōu)閯幽?，并使部分火藥燃?xì)鈴闹仆似鱾?cè)孔排出，提供與后坐方向相反的作用力，達到減少后坐能量的目的［2-3］。然而，加裝膛口制退器后會增大炮手區(qū)域的沖擊波峰值，一般對炮手危害程度隨制退效率的提高而增大。因此，為了提高膛口制退器效率，減弱或消除其副效應(yīng)，需對制退器內(nèi)外流場進行研究，完善膛口制退器設(shè)計理論。

目前，由于受測試方法限制，通過高速攝影技術(shù)捕捉膛口外流場，卻無法對制退器內(nèi)流場進行顯示。近年來，國內(nèi)外學(xué)者開始應(yīng)用計算流體力學(xué)對膛口流場與制退器效率進行計算與研究。如Sakamoto等［4］利用MUSCL 格式計算了三維膛口流場。Aibarov 等［5］則對三維帶化學(xué)反應(yīng)的膛口流場進行了數(shù)值模擬。Cler 等［6］分別采用Fluent 6.1 中2 階無粘求解器和基于Galerkin 方法的DGM 代碼以及非結(jié)構(gòu)化網(wǎng)格對不考慮運動彈丸的膛口流場(G3 步槍發(fā)射7.62 mm 北約制式槍彈)進行了數(shù)值模擬，并與試驗陰影對照，結(jié)果表明初始流場與試驗陰影相符，但火藥燃?xì)饬鲌霾⒉惶呛?。隨后，Bin 等［7］基于軸對稱Euler 方程的DRP 格式及結(jié)構(gòu)化網(wǎng)格對膛口(G3 步槍發(fā)射7.62 mm 北約制式槍彈)流場進行了模擬，所得初始流場及火藥燃?xì)饬鲌鼍c試驗陰影相吻。樂貴高等［8］利用Osher 格式及非結(jié)構(gòu)網(wǎng)格對某火炮膛口氣流場進行了三維求解。姜孝海等［1，9］基于ALE 方程的結(jié)構(gòu)化動網(wǎng)格和嵌入網(wǎng)格法以及Roe 格式對彈丸由膛內(nèi)高壓氣體推動射出到完全飛離初始流場的整個過程進行了數(shù)值模擬，并討論了初始流場、火藥燃?xì)饬鲌龅男纬膳c發(fā)展以及與彈丸相互作用的過程。

關(guān)于制退器流場的模擬則相對較少。Yang等［10］利用迎風(fēng)格式及結(jié)構(gòu)化動網(wǎng)格技術(shù)對彈丸飛離炮雙擋板制退器過程中二維膛口流場進行了數(shù)值模擬，并與試驗測得超壓峰值及超壓穩(wěn)定值進行對比。張輝等［3］采用三維任意坐標(biāo)系下Euler 方程對某口徑火炮帶3 排側(cè)孔制退器的膛口流場進行了數(shù)值模擬，并提供有關(guān)膛口流場的數(shù)據(jù)。江坤等［11］采用無粘模型及非結(jié)構(gòu)化網(wǎng)格對裝有炮口3 排側(cè)孔制退器車載炮膛口沖擊波場進行了數(shù)值模擬，并對制退器效率進行估算。

膛口流場包含高速運動彈丸、射流、沖擊波、湍流與燃燒及爆轟等相互作用，另外，由于高溫、高壓火藥燃?xì)獾募眲∨蛎浺约芭c制退器內(nèi)壁面的激烈碰撞，使制退器內(nèi)流場結(jié)構(gòu)同樣具有復(fù)雜高瞬態(tài)特點。因此，膛口內(nèi)、外流場以及制退器效率的準(zhǔn)確計算，對數(shù)值模擬要求很高，首先數(shù)值格式需準(zhǔn)確捕捉彈丸發(fā)射過程中激波，另外還要對彈丸運動以及二次爆燃等進行計算與描述，雖然非結(jié)構(gòu)化網(wǎng)格在此方面有優(yōu)點，但因耗散過大，會使格式計算精度下降?；诖?，本文采用三維非定常Euler 方程，結(jié)合能準(zhǔn)確捕捉激波的高精度Roe 格式和結(jié)構(gòu)化動網(wǎng)格技術(shù)，對彈丸飛離開腔式膛口制退器過程中，制退器內(nèi)與膛口流場結(jié)構(gòu)的變化過程進行數(shù)值模擬。同時根據(jù)數(shù)值結(jié)果對制退器的計算效率進行計算，并與文獻［12］的試驗結(jié)果進行對比。

1 數(shù)值方法及計算模型

1.1 控制方程

采用三維非定常Euler 方程

式中:U=［ρ，ρu，ρv，ρw，E］T;F=［ρu，ρu2+p，ρuv，ρuw，(E+ p)u］T;G=［ρv，ρuv，ρv2+ p，ρvw，(E+p)v］T;H=［ρw，ρuw，ρvw，ρw2+p，(E+p)w］T.其中:ρ 為氣體密度;u、v、w 分別為x、y 和z 方向的速度分量;E 為總能量，其表達式為

式中r 為理想氣體絕熱指數(shù)。理想氣體的狀態(tài)方程為

式中R 為通用氣體常數(shù)。方程(1)式～(3)式組成封閉方程組。

1.2 數(shù)學(xué)方法及驗證

采用有限體積法對方程進行離散。時間推進采用2 階精度Runge-Kutta 法，對流項則用Roe 格式離散為

式中:ARoe，j+1/2為Roe 矩陣;為Roe 通量。

為了檢驗本文算法對膛口流場復(fù)雜波系的捕捉能力，取7.62 mm JS 狙擊步槍膛口流場為算例。圖1為t=0.43 ms 時，初始膛口流場的計算陰影(上)與壓力(下)等值線和文獻［12］的試驗陰影的對比圖。由圖1可見，計算結(jié)果能成功捕捉膛口流場中瓶狀激波等典型射流結(jié)構(gòu)，并與試驗結(jié)果吻合。

圖1 計算結(jié)果與試驗結(jié)果對比Fig.1 Comparison of calculated and experimental results

1.3 計算模型及動網(wǎng)格方法

帶開腔式膛口制退器的37 mm 火炮計算模型與計算域網(wǎng)格分布如圖2所示。初始時，彈丸底部到炮膛端部距離為2 300 mm，外流場計算域取半徑為85 mm、高為384 mm 的圓柱體。為了減少網(wǎng)格帶來的格式耗散，將整個計算域劃分為96 萬個6 面體結(jié)構(gòu)化網(wǎng)格系統(tǒng)。而計算膛口制退器效率則采用200 萬個6 面體結(jié)構(gòu)化網(wǎng)格。

圖2 帶膛口制退器的火炮結(jié)構(gòu)示意圖Fig.2 Schematic diagram of cannon with muzzle brake

為了模擬彈丸相對炮管的運動過程，把整個計算域劃分為靜網(wǎng)格區(qū)和動網(wǎng)格區(qū)，如圖2(b)中彈丸前后區(qū)域為動網(wǎng)格，而運動分界面將動靜網(wǎng)格區(qū)分開。計算過程中，彈丸以867 m/s 的恒定速度沿x軸飛離膛口制退器，彈丸前后的動網(wǎng)格相對于靜網(wǎng)格滑動，其前端網(wǎng)格被不斷合并，而后端則不斷生成新網(wǎng)格，以彌補彈丸運動造成的尾部空白。

計算所需的初始條件如彈底身管內(nèi)的初始壓力、速度與火藥燃?xì)鉁囟鹊雀鶕?jù)內(nèi)彈道計算得到，其結(jié)果如圖3所示，膛口外流場為常溫常壓大氣條件。

圖3 彈丸離開膛口瞬間的膛內(nèi)燃?xì)馑俣群蛪毫皽囟确植糉ig.3 The velocity，pressure and temperature distribution of the combustion gas inside the barrel at the moment of the projectile ejecting from the muzzle

2 結(jié)果與討論

2.1 流場分析

圖4為火炮發(fā)射時，膛口及開腔式膛口制退器的初始流場計算紋影與試驗陰影對比圖［12］。此時，彈丸仍在膛內(nèi)，受其底部高壓氣推擠在管內(nèi)加速運動，同時壓縮彈頭空氣，使其在彈丸離開膛口前以欠膨脹射流噴出制退器與炮管膛口，形成初始射流。計算結(jié)果清晰地捕捉到初始流場的發(fā)展過程，其中各出口射流邊界、馬赫盤、入射激波、反射激波等所組成的復(fù)雜瓶狀激波系清晰可見(見圖4(a))，并跟試驗結(jié)果(見圖4(b))相符。

圖4 流場結(jié)構(gòu)計算結(jié)果與試驗結(jié)果對比Fig.4 Comparison of calculated and experimental results

當(dāng)火藥燃?xì)鈬姵雠诳诤?，由于制退器分流作用，在制退器?nèi)形成十分復(fù)雜的流場結(jié)構(gòu)。圖5為本文計算條件下，開腔式膛口制退器內(nèi)壓力和密度等值線隨時間的變化圖，其中每幅圖上、下各半分別為壓力和密度?？芍?，彈丸飛離膛口后，其底部高溫、高壓火藥燃?xì)饩o隨其流出膛口，在制退器內(nèi)迅速膨脹，并與制退器內(nèi)壁面及彈丸底部相互作用與反射(見圖5(a)).隨著制退器內(nèi)壁面上、下兩道反射激波在軸線上的相互碰撞以及它們與膛口高壓燃?xì)獾南嗷プ饔?，最終在制退器內(nèi)部耦合成一道強激波(見圖5(b))，提高當(dāng)?shù)貕毫懊芏?，造成膛?nèi)氣流短暫擁塞。由于膛口高速膨脹火藥燃?xì)馀c制退器內(nèi)部反射流場的相互作用，使激波形狀發(fā)生變化，隨著制退器反射激波的減弱，膛口高壓燃?xì)饬髦饾u呈現(xiàn)強勢，在制退器內(nèi)高壓膨脹使激波變形，并在制退器壁面與強激波間形成二次激波(見圖5(c)).隨后，膛口高壓燃?xì)饬髡紦?jù)統(tǒng)治地位，推動激波在制退器軸向出口內(nèi)形成三維斜激波(見圖5(c)～圖5(e))，同時，二次激波(見圖5(d))因膨脹占優(yōu)而消散(見圖5(e)).最終在制退器內(nèi)部及膛口形成復(fù)雜波系結(jié)構(gòu)(見圖5(f))，并在持續(xù)短暫時間后，最終消失。

圖6則為t=0.2 ms 時，膛口制退器內(nèi)的截面速度矢量圖。在xOy 平面上，由于制退器側(cè)孔分流作用，部分氣流加速向側(cè)孔排出，然而因初期受制退器內(nèi)壁面反射氣流的影響，在制退器內(nèi)形成明顯的速度分界面(見圖6(a)).而從y=43 mm 的xOz 平面速度矢量圖可見，前壁面上氣流除了向軸向反射外，還有部分氣流沿著側(cè)孔內(nèi)壁面發(fā)生旋轉(zhuǎn)(見圖6(b))，因此在xOz 面上形成復(fù)雜渦結(jié)構(gòu)。此外，與膛口高速膨脹氣流的碰撞則形成明顯的速度分界層。由此可知，高溫、高壓火藥燃?xì)庠谥仆似魅S空間內(nèi)的膨脹加速、碰撞、反射以及旋轉(zhuǎn)使膛內(nèi)流場結(jié)構(gòu)具有三維復(fù)雜特征。

圖5 壓力(上部)和密度(下部)等值線分布圖Fig.5 Distributions of the isobars (the upper half)and the isopycnics (the lower half)

圖7為不同時刻膛口制退器x 軸上壓力分布曲線。圖8為相應(yīng)時刻膛口制退器y 軸上壓力分布?？芍?，高溫、高壓火藥燃?xì)饬鞒鎏趴诤?，在制退器?nèi)急劇膨脹，迅速提高出口附近的氣流速度(見圖7膛口)。隨后，前導(dǎo)激波與膛口制退器內(nèi)壁面相碰撞，提高當(dāng)?shù)貕毫?見圖7，t 為0.2～0.4 ms)形成出口通道阻塞，使彈孔出口(見圖7彈孔)附近氣流速度急劇下降，因而在0.2～0.6 ms 期間，壓力在x 軸上的分布呈現(xiàn)出較為明顯的突變特性。同時，由于受激波在制退器內(nèi)的來回反射影響，火藥燃?xì)饬髟趛 軸上的壓力分布出現(xiàn)間斷，且曲線形狀隨時間的變化劇烈(見圖8).可見，在火藥燃?xì)饬鞒龀跗?，制退器y軸上具有高度時空瞬態(tài)特性。然而，隨著氣流持續(xù)流出，彈孔進口(見圖7彈孔，t=0.7 ms)處的擁塞狀況逐漸得到緩解，1 ms 以后，氣流壓力x 軸及y 軸上的分布趨勢變化呈相對較緩降低趨勢。

圖6 t=0.2 ms 時開腔式膛口制退器內(nèi)速度矢量圖Fig.6 Velocity vector in muzzle brake at t=0.2 ms

圖7 不同時刻膛口制退器在x 軸上的壓力分布曲線Fig.7 Pressure distributions along x axis of muzzle brake at different time

表1為不同時刻，開腔式膛口制退器中3 孔內(nèi)的流場參數(shù)計算值，其中，各參數(shù)均為各孔的面積平均值。初始時刻，火藥燃?xì)饬骶猿羲?馬赫數(shù)＞1)流經(jīng)3 孔。對3 孔上質(zhì)量流量m·進行統(tǒng)計發(fā)現(xiàn)流進膛口的質(zhì)量流量等于彈孔和側(cè)孔流出質(zhì)量流量之和，滿足質(zhì)量守恒定律。由于強激波以及二次激波在膛口附近的影響，使得膛口出口壓力及密度在0.4 ms 時突然上升并隨后下降。在0.2～0.5 ms間，氣流擁塞在彈孔出口處，降低了彈孔的出口速度同時提高了該處壓力及密度，隨后均持續(xù)下降。同時，由表可知，t=1 ms 以前，3 孔出口氣流速度方向均與出口平面呈一定角度，并且隨著時間發(fā)生劇烈變化。t=1 ms 后，氣流逐漸趨向于90°垂直于孔平面噴出?；谝陨戏治稣f明火藥燃?xì)鈬姵雠诳诔跗?，在膛口制退器?nèi)形成具有時空瞬態(tài)特性的三維燃?xì)饬鲌觥?/p>

圖8 不同時刻膛口制退器在y 軸上的壓力分布曲線Fig.8 Pressure distributions along y axis of muzzle brake at different time

2.2 制退效率計算

取37 mm 火炮彈丸質(zhì)量mb=0.764 7 kg，裝藥量ω=0.221 kg，彈丸出炮口初速v0=867 m/s.圖9為火炮膛底、制退器以及合力隨時間變化的數(shù)值計算分布曲線。圖10為火炮膛底沖量Io以及總沖量In隨時間變化的數(shù)值計算分布曲線?？芍尤腴_腔式膛口制退器后，由于部分火藥燃?xì)馔ㄟ^制退器側(cè)孔流出，提供了與膛底方向相反的作用力，從而降低了炮身的合力。由圖10可知，加裝制退器后降低了后坐總沖量，從而改善了火炮的后坐性能。當(dāng)身管內(nèi)氣體基本排空時(管內(nèi)壓力降到2 atm 左右)，火藥氣體作用在身管上的沖量為I0=-169.3 N·s，而帶開腔式膛口制退器后作用在身管上的總沖量則為In=5.75 N·s.定義炮口裝置火藥氣體作用系數(shù):

表1 開腔式膛口制退器3 孔的流場參數(shù)Tab.1 The flow field parameters of the three-way muzzle brake

無炮口制退器的火藥氣體作用系數(shù):

制退器效率:

圖9 火炮膛底和制退器及總受力曲線Fig.9 The curves of chamber bottom force，muzzle brake force and total thrust

此時，計算所得開腔式制退器效率為ηT=34.12%，而文獻［12］自由后坐臺相關(guān)試驗結(jié)果約為33%左右，計算結(jié)果基本與試驗結(jié)果相符合。

圖10 火炮膛底沖量I0及總沖量In曲線Fig.10 Distribution of chamber bottom impulse I0 and total impulse In

3 結(jié)論

基于三維Euler 方程，結(jié)合高精度Roe 格式與結(jié)構(gòu)化動網(wǎng)格技術(shù)，對膛口加裝開腔式制退器的彈丸發(fā)射與運動過程進行了數(shù)值模擬。計算結(jié)果表明，雖然膛口加裝了制退器，其膛口流場典型特征不變，但因火藥燃?xì)庠谥仆似鲀?nèi)高速膨脹，并與制退器內(nèi)壁面及彈丸相互作用，在膛內(nèi)形成復(fù)雜波系結(jié)構(gòu)，且具有時空瞬態(tài)特性。另外，由于制退器分流作用，部分氣體由側(cè)孔排出，從而降低后坐總沖量?；跀?shù)值計算結(jié)果，所得開腔式制退器制退效率為34.12%.計算流場結(jié)構(gòu)與制退效率都與文獻［12］試驗結(jié)果相吻，可為未來類似研究提供重要參考。

致謝本文研究與撰寫過程中得到了李鴻志院士的悉心指導(dǎo)與支持，在此表示感謝。

References)

［1］ 姜孝海，李鴻志，范寶春，等.基于ALE 方程及嵌入網(wǎng)格法的膛口流場數(shù)值模擬［J］.兵工學(xué)報，2007，28(12):1512-1515.JIANG Xiao-hai，LI Hong-zhi，F(xiàn)AN Bao-chun，et al.Numerical simulation of muzzle flow field based on ALE equation and chimera grids［J］.Acta Armamentarii，2007，28(12):1512-1515.(in Chinese)

［2］ 張月林.火炮反后坐裝置設(shè)計［M］.北京:國防工業(yè)出版社，1984:2-55.ZHANG Yue-lin.Design of the reverse recoil device of the cannon［M］.Beijing:National Defense Industry Press，1984:2-55.(in Chinese)

［3］ 張輝，譚俊杰，崔東明.帶制退器膛口流場的三維數(shù)值模擬［J］.力學(xué)與實踐，2008，30(4):30-34.ZHANG Hui，TAN Jun-jie，CUI Dong-ming.Numerical simulation of 3D muzzle brake flow field［J］.Mechanics and Practice，2008，30(4):30-34.(in Chinese)

［4］ Sakamoto K，Matsunnaga K，F(xiàn)ukushima J，et al.Numerical analysis of the propagating blast wave in a firing range［C］//19th International Symposium In Ballistics.Interlaken，Switzerland:IBC，2001:289-296.

［5］ Aibarow A V，Babayev D B，Mironov A A.Numerical simulation of 3D muzzle brake and missile launcher flow field in the presence of movable objects［C］//20th International Symposium In Ballistics.Orlando:IBC，2002:225-232.

［6］ Cler D L，Chevaugeon N，Shephard M S，et al.CFD application to gun muzzle blast a validation case study，AIAA-2003-1142［R］.Reno，Nevada:AIAA，2003.

［7］ Bin J H，Kim M W，Lee S.A numerical study on the generation of impulsive noise by complex flows discharging from a muzzle［J］.Int J Numer Meth Engng，2008，75:964-991.

［8］ 樂貴高，馬大為，馮勇，等.某火炮膛口流場數(shù)值仿真［J］.兵工學(xué)報，2004，25(1):19-22.LE Gui-gao，MA Da-wei，F(xiàn)ENG Yong，et al.Numerical simulation of muzzle blast flow fields of large caliber guns［J］.Acta Armamentarii，2004，25(1):19-22.(in Chinese)

［9］ Jiang X H，Chen Z H，F(xiàn)an B C，et al.Numerical simulation of blast flow fields induced by a high-speed projectile［J］.Shock Waves，2008，18:205-212.

［10］ Yang G，Causon D M，Ingram D M.Calculation of compressible flows about complex moving geometries using a three-dimensional Cartesian cut cell method［J］.Int J Numer Meth Fluids，2000，33:1121-1151.

［11］ 江坤，錢林方.某火炮炮口制退器性能的研究［J］.彈道學(xué)報，2006，18(3):55-57.JIANG Kun，QIAN Lin-fang.A study on the muzzle brake performance of guns［J］.Journal of Ballistics，2006，18(3):55-57.(in Chinese)

［12］ 李鴻志，姜孝海.彈道實驗圖集［R］.南京:瞬態(tài)物理國家重點實驗室，2009.LI Hong-zhi，JIANG Xiao-hai.The album of experimental ballistics［R］.Nanjing:National Key Laboratory of Transient Physics，2009.(in Chinese)