程誠，張小兵

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

0 引言

火炮內彈道過程是一種典型的伴有化學反應的高溫、高壓、可壓縮瞬態(tài)兩相流動，其循環(huán)過程中存在多種物理強間斷問題，如激波、火焰波以及爆轟波等，因此火炮膛內兩相氣動力過程可以簡化成帶有非齊次源項的黎曼間斷守恒模型［1］。這種黎曼間斷的存在成為兩相流數(shù)值方法應用的難點。

對于內彈道過程的黎曼間斷問題，一般無法得到解析解，只能通過數(shù)值計算方法求解。早期的Lax-Wendroff 兩步差分格式在強間斷附近產(chǎn)生較強的震蕩，在間斷處常常會產(chǎn)生非物理解。在內彈道兩相流領域得到了廣泛應用的MacCormck 格式在進行兩相流數(shù)值模擬時也容易產(chǎn)生非物理震蕩導致計算中斷，所以在方程中必須要加入顯式人工粘性以及濾波等方式才保證計算的穩(wěn)定性，這樣必然影響計算的精度。而后來引進的TVD 格式，其將氣體—固體兩相分開處理，氣相使用TVD 格式，固相仍然使用MacCormark 格式［2］，這種處理方式不僅給程序的編制帶來了很大困難，而且在固相以及氣固方程耦合方面存在一定的精度損失。而最近被廣泛關注的CE/SE 方法其計算單元以及計算格式相對復雜，計算機時較長，且該“菱形網(wǎng)格”推廣到多維(尤其是三維)及高階問題上時會遇到很多較難解決的問題，特別是在進行高雷諾數(shù)流動問題時，其計算結果也不太理想［3］。

Roe 提出了基于近似黎曼解的守恒Roe 格式，該方法不僅延續(xù)了Godunov 型格式接觸間斷處的高分辨率和激波捕捉性能強的特點，而且具有低耗散、高計算穩(wěn)定性以及計算時間相對較短的優(yōu)點，是當前計算流體力學領域最受歡迎的格式之一［4］。傳統(tǒng)Roe 格式是一種只有1 階精度的迎風格式，為了提高計算精度以及消除耗散帶來的誤差，通過引入帶有TVD 性質的2 階修正項，對格式進行高階重構，可以使其達到2 階精度［5］。

點火管作為火炮裝藥結構中極其重要的部分，其管內兩相流動物理過程具有典型的高溫高壓強瞬態(tài)間斷特點，可用黎曼間斷守恒模型描述該過程的物理現(xiàn)象。本文以中心點火管為例，通過高階近似黎曼重構，獲得了兩相流模型的高階近似黎曼離散格式，并通過數(shù)值計算描述了點火管內復雜兩相流動的物理過程，分析了不同因素對點火管點火性能的影響。

1 兩相流數(shù)學模型

火炮膛內點火管的燃燒與流動是火藥顆粒燃燒產(chǎn)生的氣相與未燃顆粒的兩相流動過程，根據(jù)文獻［2］中的假設及物理模型描述，由于固相被考慮成不可壓縮模型，故可忽略固相能量方程。點火管內一維兩相流質量、動量和能量守恒方程可以描述為

式中U、F、S 分別為守恒矢量、對流項和源項，其具體表達式及相關輔助方程可參見文獻［2］。

2 黎曼近似模型

2.1 Roe 格式通量分解

近似黎曼求解與其他方法的區(qū)別在于通量的求解上，其通量求解時將方程線性化，然后通過近似Roe 平均，求得各個局部黎曼問題上的通量，從而擴大到整個流場。通過Jacobian 矩陣可以實現(xiàn)對雙曲方程的線性化。對(1)式進行Jacobian 轉換后得

式中A(U)為Jacobian 矩陣系數(shù)

對(1)式進行1 階Roe 通量離散

式中Fn為控制體單元界面處的數(shù)值通量，可定義為

2.2 Roe 平均參數(shù)

通常在純氣相中往往使用密度變量ρg作為Roe 平均基本參數(shù)，而在對于兩相流問題來說，建議采用單位氣相質量ρgφg作為基本參數(shù)。這樣就可以通過Roe 平均化得到平均氣相密度、平均氣相速度、平均固相速度、平均空隙率等各變量的Roe 平均值。

2.3 熵修正

Roe 格式并不滿足Lax 熵穩(wěn)定條件，容易產(chǎn)生膨脹波問題以及在高馬赫數(shù)時出現(xiàn)“Carbuncle”現(xiàn)象等非物理現(xiàn)象，從而失去計算穩(wěn)定性。為了滿足熵穩(wěn)定，一般需要引入額外修正來保證音速膨脹區(qū)域的Jacobian 矩陣特征值不為0。熵修正方法的形式較多，對于不同問題的表現(xiàn)性能各異。通過數(shù)值比較，本文使用通過改進的Harten-Hyman 型熵修正可以保證足夠的計算穩(wěn)定性［4］。

2.4 高精度格式構造

2.1 節(jié)描述的Roe 方法僅有1 階精度，對于間斷問題為了得到更高精度的計算結果，必須引入高階格式，但高階項的引入又帶來了間斷處的數(shù)值震蕩，所以需要對其使用限制因子，以保證計算精度與計算穩(wěn)定性。采用通量修正法，1 階通量形式(5)式可寫成如下的2 階形式:

式中通量限制器采用混合限制器，即對氣相密度項使用Superbee 限制器，其余量使用Van Leer 平均限制器［4］。

2.5 源項處理

針對膛內兩相流這種帶有源項的非線性雙曲型守恒系統(tǒng)，一般有兩種方法處理:一種是將源項與通量項耦合在一個離散過程中進行求解;另外一種方法是通過時間步分裂法將帶有源項的非線性雙曲守恒方程分裂成對流項與源項分別求解，對流項使用近似黎曼求解，源項使用常微分方法求解［5］。本文使用后者，該方法可以在保證源項處理精度的基礎上，對通量項的求解使用不同類型的格式進行求解，保證了計算的穩(wěn)定性。

3 數(shù)值驗證及實驗比對

如何保證數(shù)值模擬的準確性，不僅需要與實驗結果的比對，更需要使用大量已經(jīng)被廣泛認可并帶有精確解的經(jīng)典算例進行數(shù)值驗證。針對本文這種帶有源項的非線性雙曲型守恒系統(tǒng)，下文首先應用激波管和源項的經(jīng)典算例對算法進行數(shù)值驗證，然后通過點火管算例與實驗結果進行了比對，進一步驗證數(shù)值模擬的正確性。

3.1 激波管算例

激波管問題是進行數(shù)值格式驗證的經(jīng)典算例，本文采用的激波管算例總長為2 m，網(wǎng)格數(shù)為100，并在x=0 處發(fā)生初始間斷，左右初值條件分別為:pL=100 000 Pa，ρL=1 kg/m3，uL=0;pR=1 000 Pa，ρR=0.01 kg/m3，uR=0.如圖1所示，1 階、2 階格式在光滑處都可以與解析解有很好的符合，但是在間斷處，1 階格式出現(xiàn)了明顯的過度光滑，抹去了間斷處的真實物理信息。而2 階格式與精確解的符合基本一致，說明了基于Roe 格式的2 階精度格式是擁有足夠計算精度的。

圖1 激波管1 階、2 階精度與精確解的對比Fig.1 Comparison between the numerical solution and exact solution of shock tube

3.2 源項算例

下面選用Lowe 等［6］首次提出的帶有源項的歐拉方程，檢驗分裂源項方法的計算精度。針對(1)式，將式中各復合變量用下式代替，就構成了一組無粘可壓縮純氣相歐拉方程，具體形式為

其源項中的各變量選用參見文獻［6］.

具體的物理模型選取與激波管一樣，為了更好捕捉源項的變化過程，將計算區(qū)域劃分成500 個網(wǎng)格，其初始特征量G0=1 294 301 kg/(m3·s)，N=0.1，初始壓力p0=10 140 Pa，初始當?shù)芈曀賑0=330 m/s，初始氣相速度ug0=0 m/s.

圖2所示為氣相速度的計算值與精確值的對比，通過與解析解［6］比較發(fā)現(xiàn)二者有很好的相似性，其相對誤差均小于1%，說明了分裂格式的源項處理在基于2 階Roe 格式的近似黎曼求解中是值得相信的，并具有相當高的精度。

圖2 源項驗證的氣相速度數(shù)值解與精確解對比Fig.2 Comparison between the numerical and analytical solutions of gas velocity for verification of sources

3.3 點火管模型實驗驗證

本節(jié)采用文獻［2］中描述的金屬點火管物理模型，不考慮點火管內的燃燒流動與火炮膛內過程的相互作用，僅對點火管在一個大氣壓條件下的兩相流動及燃燒過程進行研究?；境跏紬l件為:點火管長200 mm，在100～200 mm 的管壁上均勻對稱分布20 個小孔，氣體—固體相速度均為0，初始溫度為298 K，壓力為一個大氣壓，破孔壓力為20 MPa，其他條件為初始裝填條件。點火管兩端均為固壁，采用鏡面反射法來處理邊界［2］。圖3給出了點火管1/2 處壓力隨時間變化的計算結果與實驗結果對比，二者的變化趨勢一致，計算結果符合較好，說明了近似黎曼數(shù)值模擬的準確性。

圖3 壓力計算曲線與實驗對比Fig.3 Comparison between the calculated and measured pressure-time tracts

4 結果與討論

為了進一步驗證高階黎曼近似模型在火炮兩相流應用中的可靠性及穩(wěn)定性，針對3.3 節(jié)所述點火管算例，通過近似黎曼解數(shù)值模擬，描述了點火管內復雜兩相流動物理過程，并分析了不同因素對點火管點火性能的影響。

4.1 點火管數(shù)值計算結果與分析

圖4為不同時刻壓力分布規(guī)律等值線圖，其清晰地反映了不同時刻的壓力波陣面在點火管內的傳播過程。開始階段點火管左端靠近底火處首先被點燃，而右端還未被點燃，從而形成一右行壓力梯度。隨著點火波陣面的傳播，壓力保持穩(wěn)定上升。當管內壓力達到破孔壓力后，管內氣相與固相顆粒在管內外壓力差的作用下噴射到管外，此時管內破孔處會出現(xiàn)壓力瞬時驟降，如圖4中23.66 MPa 等值線，其出現(xiàn)了很明顯的壓力震蕩。圖4中的兩條23.66 MPa等值線的對比，可以清晰地反映壓力波陣面向右傳播后在固壁處反射，右端壓力急劇增長的過程。隨著流出量的不斷增加，壓力增長也趨于平緩。當管內流出量大于生成量，管內壓力開始逐漸減小，直到全部流完為止。

圖4 不同時刻壓力分布Fig.4 Calculated pressure distribution at x-t

圖5為不同時刻管內空隙率分布規(guī)律的等值線圖。從圖中可以看出，開始階段由于左端首先被點燃產(chǎn)生氣體，左端的空隙率不斷上升。由于右行壓力梯度的作用，結合圖7不同時刻的固相速度分布圖，可以清晰地看出顆粒在不斷地向右端運動，造成了右端顆粒的堆積，空隙率不斷下降。而當出現(xiàn)破孔時，由于氣相和固相的流出，以及破孔區(qū)域的壓力震蕩會造成管內破孔區(qū)域空隙率的紊亂，如圖5中空隙率為0.7163 等值線附近的渦旋。當管內壓力趨于平緩后，空隙率也在趨于穩(wěn)定增長。

圖6、圖7分別為點火管內不同時刻的氣相、固相速度分布圖。固相速度的分布規(guī)律與氣相基本一致，但遠小于氣相速度，這是由于固相顆粒慣性引起的相對運動滯后。未破空前，點火管左端率先被點燃產(chǎn)生氣相速度，并帶動固相產(chǎn)生固相速度。隨著燃燒的進行，右行速度逐漸增加。當壓力梯度沿右端壁反射后，在左行壓力梯度的作用下，產(chǎn)生圖6、圖7中的負向速度區(qū)域。隨著管內火藥顆粒的全面著火，燃燒逐漸趨于穩(wěn)定，管內壓力分布也逐漸平緩，此時管內的氣相、固相速度分布也逐漸穩(wěn)定并趨向于0。當管內壓力達到破孔壓力后，管內物質不斷流出，打破了管內速度的平衡，左端以及右端未破孔區(qū)域的氣相、固相物質分別向破孔區(qū)域流動，產(chǎn)生了較大的流動速度。直到燃燒后期，隨著流量的逐漸減少，管內流動將會逐漸趨于平緩。

圖8為不同時刻破孔區(qū)域的流量分布圖。從圖中可以看出，當破孔出現(xiàn)后，在壓力的作用下管內流量迅速增加，然后隨著壓力的變化，流量開始逐漸減小。結合圖4不同時刻壓力分布等值線圖，當開始破孔時，管內燃氣生成量大于流出量，這時管內壓力還是處于一個緩慢上升的過程，壓力逐漸達到峰值，例如圖4中的2.36 ms 處。但隨著流出量的不斷增加，管內氣體生成量已經(jīng)無法填補由于流出量而造成的壓力空白，從而壓力開始逐漸下降。

圖9為不同時刻顆粒表面溫度及點火波陣面的分布圖。從圖中可以清晰地看出，點火過程中不同時刻不同位置上的火藥顆粒表面溫度變化情況，以及可以清晰地找到一條點火火焰沿軸向的傳播過程曲線，即點火波陣面。從圖中可以看出點火波陣面與圖4中點火階段的壓力梯度有一定的聯(lián)系。一般認為，點火陣面處在壓力梯度最陡峭的地方［1］。

圖5 不同時刻空隙率分布Fig.5 Calculated porosity distribution at x-t

圖6 不同時刻氣相速度分布Fig.6 Calculated gas velocity distribution at x-t

圖7 不同時刻固相速度分布Fig.7 Calculated solid velocity distribution at x-t

圖8 不同時刻破孔流量分布圖Fig.8 Calculated mass flow rate distribution of vent holes at x-t

圖9 不同時刻顆粒表面溫度及點火波陣面分布Fig.9 Calculated particle temperature and ignition wavefront distribution at x-t

4.2 不同因素對點火管性能的影響

大量的實驗及數(shù)值模擬結果表明［1］，不同的點火性能對壓力波的發(fā)展、火焰波的傳播以及整個內彈道性能有著顯著的影響。點火管點火性能的主要因素包括傳火孔位置、傳火孔面積以及點火位置等，其中第1 排傳火孔的位置是控制點火管內壓力變化以及破孔規(guī)律的重要影響因素。為了進一步驗證基于近似黎曼解模型所編制程序的適應能力，以及研究不同傳火孔位置對點火管性能的影響，下面通過黎曼近似模型對其進行分析。

圖10為點火管內不同破孔位置壓力波對比圖。從圖中可以看出，隨著第1 排傳火孔位置的逐漸后移，點火管內壓力逐漸提高，并且點火管內的最大壓力波幅也會隨著明顯增大。這與理論分析及實驗描述［1］均完全相一致。這種點火壓力的提高，勢必會造成膛內最大壓力的上升，并且引起很大程度上的點火延時，給發(fā)射安全帶來隱患。所以一般在選取第1 排傳火孔位置時盡量靠近左端，但是相關研究表明［7］，當點火管過于靠近左端時往往會造成炮口動能的降低，達不到火炮設計要求。因此需要考慮最大膛壓和炮口速度的綜合變化情況，對點火管第1 排破孔位置進行優(yōu)化選擇。

圖10 不同破孔位置壓力波動對比Fig.10 Comparison of pressure wave at different locations

5 結論

1)基于Roe 的近似黎曼解，通過具有TVD 性質的高階重構以及兩步分裂源項處理，獲得了高階精度的內彈道兩相流近似黎曼解離散格式。

2)通過與激波管、源項驗證等具有解析解的經(jīng)典算例進行對比，驗證了所構造高階方法的準確性與穩(wěn)定性。點火管算例的數(shù)值驗證，準確形象地描述點火管內復雜兩相流動的實際物理過程，實驗驗證了該方法處理火炮兩相流問題的可靠性。

3)研究了合理設計傳火孔位置可以改善點火管內瞬時性、均勻性及安全性等點火性能，并充分說明了所建立模型具有較強的適應能力。

4)以高精度近似黎曼解為基礎，建立的火炮內彈道兩相流模型，無需濾波、粘性修正以及將氣體—固體兩相分開處理，提高了對膛內強間斷以及復雜波系的捕捉能力，改善了數(shù)值模擬的計算精度。

References)

［1］ 金志明，袁亞雄，宋明.現(xiàn)代內彈道學［M］.北京:北京理工大學出版社，1992.JIN Zhi-ming，YUAN Ya-xiong，Song Ming.Modern interior ballistic［M］.Beijing:Beijing Institute of Technology Press，1992.(in Chinese)

［2］ 袁亞雄，張小兵.高溫高壓多相流體動力學基礎［M］.哈爾濱:哈爾濱工業(yè)大學出版社，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］ Chang S C.The method of space-time conservation element and solution element:a new approach for solving the Navier-Stokes and Euler equations［J］.Journal of Computational Physics，1995，119:259-324.

［4］ 閻超.計算流體力學方法及應用［M］.北京:北京航空航天大學出版社，2006.YAN Chao.Method and application of computational fluild dynamics［M］.Beijing:Beijing University of Aeronautics and Astronautics Press，2006.(in Chinese)

［5］ Toro E F.Riemann solvers and numerical methods for fluid dynamics［M］.3rd.Berlin:Springer，2009.

［6］ Lowe C，Clarke J.A class of exact solution for the Euler equations with sources［J］.Mathematical and Computer Modeling.2002，(36):275-291.

［7］ 張會生.高初速火炮新型點傳火技術及裝藥優(yōu)化研究［D］.南京:南京理工大學，1998.ZHANG Hui-sheng.The optimization studies of new pattern ignition and flame spreading system and charge structure technology in high-velocity guns［D］.Nanjing:Nanjing University of Science and Technology，1998.(in Chinese)