程誠，張小兵

(南京理工大學(xué) 能源與動(dòng)力工程學(xué)院，江蘇 南京210094)

0 引言

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

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

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

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

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

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

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

2 黎曼近似模型

2.1 Roe 格式通量分解

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

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

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

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

2.2 Roe 平均參數(shù)

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

2.3 熵修正

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

2.4 高精度格式構(gòu)造

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

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

2.5 源項(xiàng)處理

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

3 數(shù)值驗(yàn)證及實(shí)驗(yàn)比對

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

3.1 激波管算例

激波管問題是進(jìn)行數(shù)值格式驗(yàn)證的經(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)了明顯的過度光滑，抹去了間斷處的真實(shí)物理信息。而2 階格式與精確解的符合基本一致，說明了基于Roe 格式的2 階精度格式是擁有足夠計(jì)算精度的。

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

3.2 源項(xiàng)算例

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

其源項(xiàng)中的各變量選用參見文獻(xiàn)［6］.

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

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

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

3.3 點(diǎn)火管模型實(shí)驗(yàn)驗(yàn)證

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

圖3 壓力計(jì)算曲線與實(shí)驗(yàn)對比Fig.3 Comparison between the calculated and measured pressure-time tracts

4 結(jié)果與討論

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

4.1 點(diǎn)火管數(shù)值計(jì)算結(jié)果與分析

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

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

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

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

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

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

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

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

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

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

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

4.2 不同因素對點(diǎn)火管性能的影響

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

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

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

5 結(jié)論

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

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

3)研究了合理設(shè)計(jì)傳火孔位置可以改善點(diǎn)火管內(nèi)瞬時(shí)性、均勻性及安全性等點(diǎn)火性能，并充分說明了所建立模型具有較強(qiáng)的適應(yīng)能力。

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

References)

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

［2］ 袁亞雄，張小兵.高溫高壓多相流體動(dòng)力學(xué)基礎(chǔ)［M］.哈爾濱:哈爾濱工業(yè)大學(xué)出版社，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］ 閻超.計(jì)算流體力學(xué)方法及應(yīng)用［M］.北京:北京航空航天大學(xué)出版社，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］ 張會(huì)生.高初速火炮新型點(diǎn)傳火技術(shù)及裝藥優(yōu)化研究［D］.南京:南京理工大學(xué)，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)