某制導炮彈二維兩相流內彈道性能分析與數值模擬研究

程誠，張小兵

（南京理工大學能源與動力工程學院，江蘇南京210094）

某制導炮彈二維兩相流內彈道性能分析與數值模擬研究

程誠，張小兵

（南京理工大學能源與動力工程學院，江蘇南京210094）

為了研究某制導炮彈二維兩相流內彈道性能，簡化兩相流多維數值模擬中彈底運動邊界處理的復雜性，提高運動邊界處的計算精度，建立了基于任意拉格朗日歐拉方法的某制導炮彈內彈道二維氣-固兩相流模型，空間上采用具有TVD特性的高階MUSCL類型的有限體積法對方程進行離散，時間方向采用4階Runge-Kutta方法進行時間推進。通過擁有解析解的數值驗證算例，驗證了數值格式以及動網格生成方法的準確性。對某大口徑制導炮彈內彈道膛內循環過程進行二維兩相流數值仿真。模擬結果準確地反映了整個內彈道循環膛內兩相流動特性及其發展過程，并與實驗結果有較好的一致性。同時分析了不同點火因素對該制導炮彈內彈道性能的影響，為后續深入優化該制導炮彈內彈道性能及發射安全性提供了理論基礎。

兵器科學與技術；內彈道；兩相流；數值模擬；制導炮彈

0　引言

制導炮彈的內彈道過程是一個伴隨著高溫、高壓、高過載的多相燃燒流動過程，特別是由于精確制導設備的存在，其對膛內異常壓力波動及彈丸的異常壓力過載都有著嚴格的要求，因此對制導炮彈的內彈道兩相流理論研究和設計工作都提出了更高的要求。目前對于制導炮彈的內彈道兩相流理論研究多集中在一維數值模擬階段［1］，對于制導炮彈二維甚至是三維的數值模擬研究還未見報道。

在火炮內彈道循環過程中，膛內氣-固兩相流動存在各種復雜邊界，特別是在多維計算過程中，對這些邊界的計算處理成為是否能夠成功對膛內復雜兩相流動過程進行數值模擬的關鍵技術之一。目前廣泛使用的傳統解決方法是采用歐拉坐標系下的方程組，通過對運動邊界處進行控制體守恒方程推導，從而獲得運動邊界處的各物理參量［1-2］。該方法需要在運動邊界處重新建立一套完整的計算差分方法，這樣不僅使計算程序復雜化，而且也在一定程度上影響了計算精度。

基于拉格朗日坐標系下的兩相流方程組可以將彈底邊界固定下來，從而使計算量大大降低。但在兩相流動流場中，由于連續流體線保持性定律不成立，很難建立基于拉格朗日坐標系下的兩相流動方程。為了解決這個問題，宋明等［3］提出使用擬合坐標系，把彈底邊界固定下來，計算中網格點始終保持不變，這樣避免了插值，提高了精度。但當使用該方法處理大口徑長身管火炮內彈道循環過程時，由于網格數目固定，在彈道循環后期會造成單位網格步長不斷增大，網格質量不斷下降，從而影響了對于流場中復雜激波的捕捉能力，降低了計算精度。

任意拉格朗日歐拉（ALE）方法，將拉格朗日方法和歐拉方法統一起來，允許網格以任意速度運動，特別適合用于計算包含運動邊界的流動［4-5］。對于含有運動邊界的膛內氣-固兩相流動問題，基于ALE方法的有限體積法不僅可以使包含運動邊界在內的所有網格采用統一數值差分方法，無需對運動邊界處的網格建立新的差分格式或進行單獨守恒推導，而且該方法可以較為容易地與其他網格生成方法進行結合，保證了內彈道循環后期網格重構的精度，提高了對于膛內兩相流動復雜波系的捕捉能力。

本文基于ALE方法建立了內彈道二維氣-固兩相流模型，對某型艦載制導炮彈內彈道過程進行了數值模擬，數值仿真結果詳細反映了整個內彈道循環膛內兩相流動特性及其發展過程，同時基于該方法討論了不同點火條件對該制導炮彈內彈道性能的影響。

1　ALE形式的二維氣-固兩相流方程組

對于內彈道膛內氣-固兩相流過程，由于膛內沿身管軸線方向呈對稱形式，因此在內彈道膛內氣-固兩相流二維計算時，往往采用軸對稱模型。二維軸對稱可壓、無粘ALE形式的歐拉氣-固兩相流方程可表示為

式中：U為t時刻的守恒矢變量；F和G分別為徑向r、軸向z方向上的矢通量；S為源項。U、F、G、S可寫為如下矩陣形式：

（2）式～（5）式中：φ為空隙率；ρg為氣相密度；ρp為固相顆粒密度；ug和vg分別為軸向和徑向方向氣相速度；up和vp分別為軸向和徑向方向固相速度；ub和vb分別為運動邊界在軸向和徑向方向的速度分量；p為氣相壓力；氣相總能，eg為氣體內能；為單位體積、單位時間內燃燒產生氣體質量；為從點火管流進的氣相、固相單位體積質量流量；uig和vig分別為點火燃氣在軸向和徑向方向氣相速度；uip和vip分別為點火燃氣在軸向和徑向方向固相速度；Hig為點火燃氣的滯止焓；為單位體積內固相顆粒等效比表面積；fsz和fsr分別為軸向和徑向的氣-固相間阻力；q和Rp分別為氣-固兩相相間傳熱比熱流以及顆粒間應力。式中各輔助方程參見文獻［2］.

2　數值計算方法

2.1離散格式

本文基于有限體積法，時間方向采用4階Runge-Kutta法進行時間推進，空間采用具有TVD特性的高階MUSCL格式對控制體對流項進行空間離散，具體計算步驟參見文獻［6-10］.

2.2初始條件和邊界條件

膛內初始計算條件為初始環境常量和火炮初始裝填條件。彈底與膛底的形狀比較復雜，本文假設彈底與膛底為垂直于中心軸線的平面。在彈丸未運動前，膛內所有壁面均考慮為靜止固壁邊界。在彈丸開始運動后，彈底為固壁運動邊界?？紤]邊界的運動速度，其他邊界條件參見文獻［2］.針對內彈道動邊界問題這一特殊問題，由于只考慮彈丸沿軸向方向的運動，因此僅考慮軸向方向的速度分量，則該運動邊界的速度即為彈丸運動速度，可根據彈丸運動方程進行求解。

式中：un為運動邊界外虛擬網格點n處的軸向速度。

2.3動網格

當彈底壓力大于擠進壓力后，彈丸開始運動，其彈后的計算區域也在不斷增加，因此需要對彈后網格進行重新劃分，以保證計算精度以及計算穩定性。本文采用內彈道計算中常用的動態層網格技術，即當彈后網格層高度大于某一給定網格高度時，在彈底網格層與鄰近網格層之間增加一層新網格，新增加網格處的物理量可用2階精度插值求出。

3　結果與分析

3.1數值驗證

為了驗證所建立模型以及數值方法的準確性，在程序編制過程中針對激波捕捉能力、時間方向推進精度、動邊界可靠性等方面，使用了具有精確解的經典算例進行數值驗證。由于篇幅有限下面僅給出針對ALE方法的動邊界可靠性驗證的算例。該算例描述了理想氣體壓力驅動條件下活塞自由運動過程，通過氣體動力學知識，其可以獲得解析解，其具體結構示意圖、初始計算參數以及相關解析解表達式參見文獻［11］.

圖1為動邊界驗證算例數值解與解析解對比圖。從圖中可以看出二者有很好的相似性，數值解與解析解基本重合，說明了了基于ALE的動邊界處理方法的準確性以及所編制程序的可靠性。

圖1　動邊界驗證算例數值解與解析解對比Fig.1 Numerical and analytical solutions of moving boundary

3.2膛內兩相流動過程數值模擬及實驗驗證

以帶有中心點火管結構的某制導炮彈內彈道過程的二維兩相流動為研究對象，其具體結構可參見文獻［1］中所示。為了計算方便，暫不考慮點火管的二維效應，即把點火管作為線性點火源考慮，點火管傳火孔破裂后，通過所噴射質量、動量以及能量與膛內燃燒過程進行耦合，彈底與膛底的復雜形狀也簡化為垂直于中心軸線的平面。

3.2.1實驗比對

現針對該制導炮彈內彈道方案的二維兩相流動過程進行數值模擬及實驗比對，由于內彈道實驗的局限性，實驗中僅獲得了彈丸初速以及最大膛壓，分別為851 m/s和319 MPa（銅柱壓力動態誤差修正后），而數值模擬的彈丸初速與最大膛壓分別為852.5 m/s與315.62 MPa，相對誤差都小于1%.

3.2.2數值模擬結果與分析

圖2和圖3為膛內不同時刻壓力分布圖，其清晰地反映了整個發射過程中不同時刻的壓力發展過程。開始階段，如圖2所示，由于點火管破孔區域點火燃氣不斷向膛內噴射，在徑向造成很大的壓力梯度，與此同時破孔區域藥床首先被點燃，從而該破孔區域膛內壓力逐漸上升，并持續向藥室兩端及膛壁傳播。隨著發射藥床不斷被點燃，膛內壓力持續穩定上升，其膛內的二維徑向效應也出現衰減。

圖2　彈丸未運動前膛內不同時刻壓力分布圖（單位：MPa）Fig.2 Contours of in-bore pressures at different times before projectile starts to move（unit：MPa）

圖3　彈丸運動后膛內不同時刻壓力分布圖（單位：MPa）Fig.3 Contours of in-bore pressures at different times after projectile starts to move（unit：MPa）

當彈丸開始運動后，如圖3所示，隨著燃燒的不斷進行其壓力也在不斷上升，其沿徑向的二維效應已基本消失，膛內從膛底向彈底沿軸線方向形成了一定的壓力梯度，并推動彈丸逐漸運動。值得注意的是，在坡膛區域的壓力明顯比其周圍區域高，并有一定的徑向效應存在。

為了更為詳細地描述整個發射過程中的壓力發展過程，圖4給出了不同時刻r=0.04 m處沿軸向方向的壓力分布曲線。圖4（a）～圖4（b）為早期點火階段r=0.04 m處沿軸向方向的膛內壓力分布情況，可以看出壓力在破孔區域首先不斷上升，并向膛底和彈底方向逐漸發展。圖4（c）所示，隨著燃燒不斷進行，膛內壓力也在不斷上升，當t= 9 ms左右時，壓力達到最大，并隨著彈丸的運動彈后空間逐漸增大，壓力開始不斷下降，直到彈丸離開炮口。

圖4 r=0.04 m處沿軸向方向的壓力分布Fig.4 Pressure distributions along the axis at r=0.04 m

圖5為不同時刻的氣相速度及矢量分布圖，圖中每個時間子圖中左邊為速度云圖，右邊為該時刻所對應的速度矢量圖。從圖中可以看出在點火初期徑向效應較為明顯，當火焰傳播到膛壁并反射后，徑向效應逐漸減小，軸向流動成為膛內運動的主要形式。如圖5（a）～圖5（d）所示，隨著火藥逐漸被點燃，膛內壓力不斷上升并且向膛底和彈底不斷運動，當到達膛底或彈底后又被反射做反向運動。如圖5（e）所示，當彈丸開始運動后，彈后空間出現了局部低壓區，在壓力梯度作用下，彈底區域的火藥燃氣及固相顆粒立即改變運動方向，逐步向彈底運動，而此時其他向膛底方向運動的火藥燃氣及固相顆粒還未及時改變運動方向，因此在圖中出現了一定程度的渦旋區域。如圖5（f）所示，隨著膛內壓力的逐步上升以及彈丸的不斷運動，膛內氣相和固相都向彈底不斷運動，其速度大小不斷上升，直到彈丸離開炮口。

圖5　不同時刻氣相速度（左）及速度矢量（右）分布Fig.5 Velocity（left）and vector（right）of gas phase at different times

3.2.3不同因素分析

下面基于以上二維兩相流數值模擬方法，討論不同點火條件對該制導炮彈內彈道性能的影響。

圖6為不同傳火孔孔徑D0變化條件下最大膛壓pm和炮口初速v0的變化情況。通過理論分析及已有數值研究結果可以知道，隨著傳火孔徑的增大，其傳火總面積也在增加，因此流向發射藥床的點火燃氣也就相應更多，也就是說點火激勵和點火猛度更強，從而相應的膛壓和初速都會越大。而圖6中所反映的變化規律與已有理論分析結果較為一致。

圖7為不同點火藥量mig對與最大膛壓pm和初速v0的影響情況。從圖中可以看出，隨著點火藥量的增加，點火強度也會增加，對于發射藥床的激勵程度也會越強，所對應的最大膛壓和和炮口初速也會相應增加。數值模擬所展現的物理變化規律與已有理論分析和實驗研究結果較為一致，進一步說明了所建立方法和編制程序為研究該制導炮彈內彈道性能提供了有力的理論工具。

圖7　最大膛壓和初速隨點火藥量的變化分布Fig.7 Change of maximum pressure and muzzle velocity with black powder mass in igniter

4　結論

本文以某制導炮彈內彈道性能研究為背景，建立了該制導炮彈內彈道二維兩相流數學模型并進行數值模擬與實驗比對，通過不同因素分析討論了點火條件對該制導炮彈內彈道性能的影響，具體結論如下：

1）建立了基于ALE方法的內彈道二維氣-固兩相流數學模型，空間采用具有TVD特性的高階MUSCL格式對控制體對流項進行空間離散，時間方向采用4階Runge-Kutta法進行時間推進。并通過具有解析解的經典驗證算例，驗證了所構造數學模型及數值模擬方法的準確性。

2）基于該方法對某制導炮彈內彈道兩相流過程進行了數值模擬與實驗比對，數值模擬準確形象地描述了整個內彈道復雜兩相流動的實際物理過程，實驗對比驗證了該方法處理火炮內彈道兩相流問題的可靠性。

3）利用該二維兩相流數值方法對不同點火條件下的制導炮彈內彈道性能進行了分析，為后續深入優化該制導炮彈內彈道性能提供了必要的理論基礎。

4）基于ALE方法的內彈道二維氣-固兩相流模型，無需對運動邊界處的網格建立新的差分格式或進行單獨守恒推導，保證了內彈道循環后期網格重構的精度，提高了對于膛內兩相流動復雜波系的捕捉能力。

（

）

［1］程誠，張小兵，Rashad M，等.基于高階黎曼近似解的膛內多相燃燒過程研究［J］.彈道學報，2013，25（3）：79-82. CHENG Cheng，ZHANG Xiao-bing，Rashad M，et al.Study on multi-phase combustion based on high resolution approximate Riemann solver in guns［J］.Journal of Ballistics，2013，25（3）：79-82.（in Chinese）

［2］袁亞雄，張小兵.高溫高壓多相流體動力學基礎［M］.哈爾濱：哈爾濱工業大學出版社，2005. YUAN Ya-xiong，ZHANG Xiao-bing.Multiphase hydrokinetic foundation of high temperature and high pressure［M］.Harbin：Harbin Institute of Technology Press，2005.（in Chinese）

［3］宋明，楊新民.內彈道兩相流動計算中兩類邊界條件的處理［J］.兵工學報，1993，14（3）：6-11. SONG Ming，YANG Xin-min.Treatment of two kinds of boundary conditions in the calculation of two-phase flow in interior ballistics［J］. Acta Armamentarii，1993，14（3）：6-11.（in Chinese）

［4］Trepanier J Y，Reggio M，Zhang H，et al.A finite-volume method for the Euler equations on arbitrary Lagrangian-Eulerian grids［J］. Computers&Fluids.1991，20（4）：399-409.

［5］Hirt C W，Amsden A A，Cook J，et al.An arbitrary Lagrangian-Eulerian computing method for all flow speeds［J］.Journal of Computational Physics，1974，14（3）：227-253.

［6］Van Leer B.Towards the ultimate conservative difference scheme［J］. Journal of Computational Physics，1997，135（2）：229-248.

［7］Toro E F.Riemann solvers and numerical methods for fluid dynamics：a practical introduction［M］.3rd ed.Berlin：Springer，2009.

［8］程誠，張小兵.高階近似黎曼解模型在火炮內彈道兩相流中的應用研究［J］.兵工學報，2011，32（10）：1200-1205. CHENG Cheng，ZHANG Xiao-bing.Research and application of higher-order approximate Riemann solver to two-phase flow in gun［J］.Acta Armamentarii，2011，32（10）：1200-1205.（in Chinese）

［9］Demird?i'c I，Lilek ?，Peri'c M，et al.A collocated finite volume method for predicting flows at all speeds［J］.International Journal for Numerical Methods in Fluids，1993，16（12）：1029-1050.

［10］Rahman M M，Miettinen A，Siikonen T，et al.Modified SIMPLE formulation on a collocated grid with an assessment of the simplified QUICK scheme［J］.Numerical Heat Transfer，1996，30（3）：291-314.

［11］Jacobs P A.Shock tube modeling with L1d［D］.Queensland：The University of Queensland，1998.

Two-dimensional Numerical Simulation on Two-phase Flow Interior Ballistic Performance of a Guided Projectile

CHENG Cheng，ZHANG Xiao-bing
（School of Energy and Power Engineering，Nanjing University of Science and Technology，Nanjing 210094，Jiangsu，China）

To simplify the complexity of moving boundary treatment in the two-phase flow numerical simulation and improve the computational accuracy of moving boundary treatment，a two-dimensional twophase flow model based on ALE method is established.The governing equations are discretized with the TVD-type MUSCL scheme to obtain a second order accuracy in finite volume form.A two-step and fourthorder Runge-Kutta method is used to solve the source terms.The dynamic mesh update method and numerical method are verified by a standard test with the exact solution.The numerical application on a large caliber guided projectile shows the complete development and the two-phase flow characteristics of the interior ballistic process.The calculated results are in good agreement with the experimental results. The effects of different ignition conditions on interior ballistic performance are investigated.The study lays a theoretical foundation for deep optimization of the interior ballistic performance and firing safety.

ordnance science and technology；interior ballistics；two-phase flow；numerical simulation；guided projectile

TJ012.1

A

1000-1093（2015）01-0058-06

10.3969/j.issn.1000-1093.2015.01.009

2013-11-15

江蘇省自然科學基金項目（BK20131348）；瞬態物理國家重點實驗室基金項目（9140C300206120C30110）

程誠（1986—），男，講師。E-mail：chengcheng@njust.edu.cn；張小兵（1968—），男，教授，博士生導師。E-mail：zhangxb680504@163.com