固體炸藥爆轟與惰性介質(zhì)相互作用的一種擴(kuò)散界面模型*

于 明

（1. 北京大學(xué)應(yīng)用物理與技術(shù)中心，北京 100871；2. 北京應(yīng)用物理與計(jì)算數(shù)學(xué)研究所計(jì)算物理重點(diǎn)實(shí)驗(yàn)室，北京 100083）

爆轟是復(fù)雜的流體力學(xué)與劇烈的化學(xué)動(dòng)力學(xué)耦合的一種物理現(xiàn)象。固體炸藥的爆轟能夠提供強(qiáng)大的能量來驅(qū)動(dòng)和壓縮材料，在國(guó)防科技和國(guó)民經(jīng)濟(jì)中獲得廣泛應(yīng)用，因此，固體炸藥爆轟與惰性介質(zhì)相互作用問題是工程應(yīng)用中十分重要的問題。固體炸藥爆轟與惰性介質(zhì)相互作用問題本質(zhì)上是一種可壓縮多物質(zhì)流動(dòng)問題，除了可壓縮單物質(zhì)流動(dòng)中出現(xiàn)的激波、稀疏波、接觸間斷外，一種新的間斷，即物質(zhì)界面出現(xiàn)了，它隔開熱力學(xué)性質(zhì)或狀態(tài)方程不一樣的兩種物質(zhì)，這種物質(zhì)界面實(shí)質(zhì)上也是一種接觸間斷，它滿足物質(zhì)界面壓力相等與速度相等的條件。對(duì)可壓縮多物質(zhì)流動(dòng)，傳統(tǒng)的數(shù)值模擬方法多是拉格朗日方法[1-3]，它不能有效地處理物質(zhì)大變形運(yùn)動(dòng)，為了較好地處理物質(zhì)大變形運(yùn)動(dòng)，歐拉方法成為一個(gè)很好的選擇[3-4]。根據(jù)對(duì)處理物質(zhì)界面的不同方式，歐拉方法大致可以分為4 類：陣面追蹤（front tracking）法[5-6]、流體體積（volume of fluid）法[7-8]、水平集合（level set）法[9-10]以及擴(kuò)散界面（diffuse interface）法[11-17]等，其中，擴(kuò)散界面法由于具有邏輯結(jié)構(gòu)簡(jiǎn)單、物理量守恒性良好、對(duì)界面形狀沒有幾何拓?fù)湎拗?、?duì)多物理問題適應(yīng)性較強(qiáng)的特點(diǎn)，近年來獲得了越來越多的關(guān)注。

擴(kuò)散界面法將離散網(wǎng)格視為包含有多種組分物質(zhì)的混合網(wǎng)格，物質(zhì)界面被認(rèn)為是具有一定厚度的虛擬混合區(qū)，即把無厚度的物質(zhì)界面當(dāng)作有一定厚度的擴(kuò)散界面，擴(kuò)散界面內(nèi)部采用“虛擬狀態(tài)方程”或混合規(guī)則來描述。根據(jù)對(duì)多物質(zhì)混合狀態(tài)的不同處理方式，擴(kuò)散界面模型可以分為兩種[18]：基于單物質(zhì)流動(dòng)歐拉方程組的擴(kuò)展歐拉模型（augmented Euler model）以及基于多相流動(dòng)Baer-Nunziato 方程組的多相流模型（multiphase flow model）。

在擴(kuò)展歐拉模型中，混合規(guī)則通常采用基于力學(xué)平衡和熱學(xué)平衡的假設(shè)。事實(shí)上，熱學(xué)平衡僅適用于均勻混合狀況，而包含有物質(zhì)界面的混合網(wǎng)格內(nèi)各組分物質(zhì)通常不是均勻混合的，由此將熱學(xué)平衡的假設(shè)用于計(jì)算包含界面接觸間斷的混合物質(zhì)的熱力學(xué)性質(zhì)時(shí)，往往使得物質(zhì)界面附近的壓力和速度等物理量出現(xiàn)非物理振蕩現(xiàn)象[19]。為了消除物質(zhì)界面附近的非物理振蕩現(xiàn)象，不得不采取另外的技術(shù)方案進(jìn)行修正處理，如附加狀態(tài)方程參數(shù)演化方程[20]、總能量調(diào)整[21]、守恒與原始變量轉(zhuǎn)換[22]等。

不同于擴(kuò)展歐拉模型，多相流模型[23]中的混合物質(zhì)被認(rèn)為由處于熱學(xué)、力學(xué)、化學(xué)非平衡的多種組分物質(zhì)組成，每種組分物質(zhì)具有各自的物理狀態(tài)并按照各自的動(dòng)力學(xué)規(guī)律分別進(jìn)行演化，演化方程含有用于表達(dá)由組分物質(zhì)非平衡引起的質(zhì)量、動(dòng)量和能量相互轉(zhuǎn)化的各種源項(xiàng)，同時(shí)采用組分物質(zhì)體積分?jǐn)?shù)方程來描述物質(zhì)界面的運(yùn)動(dòng)過程。由于考慮了每種組分物質(zhì)各自的變化規(guī)律，能夠保證混合物質(zhì)的熱力學(xué)性質(zhì)自動(dòng)滿足一致性，從物理建模上消除了物質(zhì)界面附近的非物理振蕩，對(duì)諸如爆轟這種帶有化學(xué)非平衡的可壓縮流動(dòng)，消除物質(zhì)界面附近的非物理壓力振蕩非常重要，因?yàn)楸Z化學(xué)反應(yīng)的激發(fā)與壓力等量密切相關(guān)，如果壓力出現(xiàn)非物理振蕩會(huì)激發(fā)虛假的化學(xué)反應(yīng)，進(jìn)一步會(huì)引起錯(cuò)誤的爆轟特性。針對(duì)具體物理過程，熱學(xué)、力學(xué)、化學(xué)非平衡這些狀態(tài)可以全部滿足也可以部分滿足，在全部滿足非平衡條件下即為著名的Baer-Nunziato 模型，在部分滿足非平衡條件下形成各種簡(jiǎn)化多相流模型，如在力學(xué)平衡、熱學(xué)非平衡條件下有Kapila 模型[24]、Allaire 模型[25]、Massoni 模型[26]、Murrone 模型[27]、Grove 模型[28]等多種典型模型，這些不同模型的差異主要在于以不同方式處理組分物質(zhì)混合狀態(tài)，使得組分物質(zhì)體積分?jǐn)?shù)方程有不同的表達(dá)形式。多相流模型由于考慮了更精確的物理機(jī)制，更加符合物理意義，近年來在可壓縮多物質(zhì)流動(dòng)數(shù)值模擬領(lǐng)域獲得了極大的重視，本文的工作正是基于多相流模型。

在固體炸藥爆轟與惰性介質(zhì)的相互作用過程中，爆轟化學(xué)反應(yīng)過程通常簡(jiǎn)化為固相反應(yīng)物轉(zhuǎn)化成氣相生成物，這樣組分混合物質(zhì)通常包括固相反應(yīng)物、氣相生成物、惰性介質(zhì)這3 種成分，由于這3 種組分物質(zhì)的材料物態(tài)性質(zhì)和熱力學(xué)性質(zhì)差異極大，因此它們之間的相互作用過程被認(rèn)為滿足力學(xué)平衡和熱學(xué)非平衡狀態(tài)[13,24,29]，即在流場(chǎng)控制體中每種組分物質(zhì)擁有相同的壓力和速度、以及不同的溫度和內(nèi)能。基于多相流思想，固體炸藥爆轟與惰性介質(zhì)相互作用過程首先由各種組分物質(zhì)的質(zhì)量守恒方程、混合物質(zhì)的動(dòng)量與總能量守恒方程描述，組分物質(zhì)的質(zhì)量守恒方程還需考慮由化學(xué)反應(yīng)引起的質(zhì)量轉(zhuǎn)化的影響，然后需要補(bǔ)充每種組分物質(zhì)的體積分?jǐn)?shù)的控制方程。由混合物質(zhì)能量守恒方程可以分解獲得每種組分物質(zhì)的內(nèi)能變化方程，再利用每種組分物質(zhì)壓力相等的條件并結(jié)合組分物質(zhì)的質(zhì)量守恒方程，可以推導(dǎo)出每種組分物質(zhì)的體積分?jǐn)?shù)控制方程。由于涉及到熱學(xué)非平衡狀態(tài)，組分物質(zhì)之間存在熱量交換，分解獲得的組分物質(zhì)內(nèi)能變化方程還需考慮熱量交換的影響。這樣，推導(dǎo)出的組分物質(zhì)體積分?jǐn)?shù)控制方程同時(shí)包含了化學(xué)反應(yīng)和熱學(xué)非平衡的影響。并且，混合物質(zhì)的壓力演化方程也被納入到模型方程組中，這樣壓力通過直接離散求解演化方程獲得而不是由流動(dòng)守恒變量計(jì)算獲得，這種方案可以增強(qiáng)消除物質(zhì)界面非物理振蕩的效果，也可以避免由狀態(tài)方程非線性形式引起的壓力迭代求解[28]或者由組分體積方程非守恒型引起的壓力松弛求解[13]，還可以避免計(jì)算壓力對(duì)守恒變量的導(dǎo)數(shù)而簡(jiǎn)化高階精度格式的運(yùn)算過程。最終，獲得的擴(kuò)散界面模型方程組包括組分物質(zhì)的質(zhì)量守恒方程、混合物質(zhì)的動(dòng)量及總能量守恒方程、組分物質(zhì)的體積分?jǐn)?shù)演化方程以及混合物質(zhì)的壓力演化方程。獲得的擴(kuò)散界面模型的最主要特點(diǎn)是考慮了化學(xué)反應(yīng)以及熱學(xué)非平衡的影響，因此具有良好的熱力學(xué)一致性，同時(shí)，該擴(kuò)散界面模型能夠適用于任意表達(dá)形式的狀態(tài)方程以及任意數(shù)目的多種惰性介質(zhì)。

對(duì)所獲得的擴(kuò)散界面模型方程組采用一個(gè)具有波傳播性質(zhì)的時(shí)空二階精度的有限體積法進(jìn)行數(shù)值求解，典型算例結(jié)果顯示，數(shù)值模擬圖像與物理規(guī)律符合，物質(zhì)界面附近不會(huì)出現(xiàn)物理量的非物理振蕩現(xiàn)象。

1 擴(kuò)散界面模型的推導(dǎo)

首先給出物理量的定義?？紤]一個(gè)控制體包含有炸藥固相反應(yīng)物、炸藥氣相生成物以及惰性物質(zhì)，固相反應(yīng)物的物理量用下標(biāo)s 表示、氣相生成物的物理量用下標(biāo)g 表示，惰性物質(zhì)設(shè)有K種物質(zhì)，每種物質(zhì)的物理量用下標(biāo)k表示；其中 ρ 表示密度，e表示內(nèi)能，p表示壓力，u表示速度矢量， α 表示體積分?jǐn)?shù)，v表示比容，Q表示單位質(zhì)量的固體炸藥由固相反應(yīng)物轉(zhuǎn)化為氣相生成物時(shí)所釋放的熱量。在力學(xué)平衡條件下，混合物質(zhì)物理量與組分物質(zhì)物理量的關(guān)系式為：

不考慮各種耗散因素及外力做功情況，則在力學(xué)平衡狀態(tài)條件下固體炸藥爆轟與惰性介質(zhì)的可壓縮流動(dòng)方程組可以表達(dá)為：

式(7)為固體炸藥固相反應(yīng)物的質(zhì)量守恒方程，式(8)為固體炸藥氣相生成物的質(zhì)量守恒方程，式(9)為惰性介質(zhì)的質(zhì)量守恒方程，式(10)為混合物質(zhì)的動(dòng)量守恒方程，式(11)為混合物質(zhì)的總能量守恒方程。

在給定每種組分物質(zhì)的狀態(tài)方程條件下，并確定了每種組分物質(zhì)的體積分?jǐn)?shù)的控制方程之后，上述流動(dòng)方程組成為封閉系統(tǒng)進(jìn)而可以獲得定解。組分物質(zhì)體積分?jǐn)?shù)的控制方程可由混合物質(zhì)的能量守恒方程并結(jié)合混合網(wǎng)格內(nèi)每種組分物質(zhì)壓力相等的條件導(dǎo)出，混合物質(zhì)的能量守恒方程可以分解為每種組分物質(zhì)的內(nèi)能變化方程。

上述流動(dòng)方程組式(7)～(11)可以經(jīng)過變換得到混合物質(zhì)的內(nèi)能守恒方程：

式(14)中3 個(gè)括號(hào)中分別表示炸藥固相反應(yīng)物、炸藥氣相生成物、惰性介質(zhì)的內(nèi)能變化，再考慮到這多種物質(zhì)由于處于熱學(xué)非平衡狀態(tài)而存在溫度差異，物質(zhì)之間會(huì)產(chǎn)生熱量交換，則可以將式(14)分裂為各種組分物質(zhì)的內(nèi)能變化方程：

同理可以得到氣相生成物與惰性介質(zhì)的體積分?jǐn)?shù)控制方程，這樣各組分物質(zhì)的體積分?jǐn)?shù)控制方程均可獲得，同時(shí)，混合物質(zhì)的壓力演化方程也成為一個(gè)流動(dòng)控制方程，這樣固體炸藥爆轟與惰性介質(zhì)相互作用的擴(kuò)散界面模型被推導(dǎo)出：

2 數(shù)值方法

為了便于數(shù)值求解，獲得的界面擴(kuò)散模型方程組(29)～(36)在二維直角坐標(biāo)系下可以寫成如下矢量-矩陣形式的偏微分方程組：

這里狀態(tài)變量q=[ρsαs, ρgαg, ρkαk, ρu, ρv, ρE, αs, αk,p]T（不失一般性，惰性介質(zhì)只寫一種物質(zhì)），矩陣A(q) 與矩陣B(q) 是具有類似表達(dá)形式的Jacobi 矩陣，其中

非齊次非守恒型雙曲律方程組(30)可用Strang 分裂方法[30]來求解，求解過程依次為雙曲律方程組求解步與常微分方程組求解步。在雙曲律方程組求解步中，一個(gè)具有波傳播特征的二階精度Godunov型格式[31]被采用，該格式能夠很好地處理雙曲律方程組中的非守恒型對(duì)流項(xiàng)。在常微分方程組求解步中，首先根據(jù)固體炸藥爆轟通常具有快反應(yīng)過程和慢反應(yīng)過程的特點(diǎn)[29]，將爆轟化學(xué)反應(yīng)率分解成快反應(yīng)率和慢反應(yīng)率兩部分，然后一個(gè)具有顯-隱離散特性的二階精度Runge-Kutta 格式[32]被采用，顯式離散用于處理慢反應(yīng)率源項(xiàng)、隱式離散用于處理快反應(yīng)率源項(xiàng)。

3 典型算例

3.1 一維爆轟增長(zhǎng)問題

圖1 一維爆轟的壓力增長(zhǎng)過程Fig. 1 Growth of pressure in one-dimensional detonation

3.2 滑移爆轟約束效應(yīng)問題

考察一平面爆轟在金屬銅的約束下向前傳播情況，計(jì)算域如圖2 所示，固體炸藥固相反應(yīng)物和氣相生成物均采用Jones-Wilkins-Lee 形式狀態(tài)方程[34]；爆轟化學(xué)反應(yīng)率采用Ignition-Growth 模型[34]；銅采用Mie-Grüneisen 形式狀態(tài)方程[18]。起爆區(qū)域設(shè)置有壓力為36.5 GPa 的高壓靜止氣相生成物，其余區(qū)域均設(shè)置為標(biāo)準(zhǔn)狀態(tài)；計(jì)算域入口設(shè)置為定值，其余邊界設(shè)置為無反射邊界條件。

圖2 滑移爆轟約束構(gòu)型圖Fig. 2 Configuration of confinement effect

計(jì)算獲得爆轟波達(dá)到定常狀態(tài)時(shí)的密度及壓力的分布，如圖3 所示，可以看出，炸藥與銅界面附近的爆轟波陣面呈彎曲狀態(tài)，爆轟波向銅介質(zhì)折射一激波。數(shù)值計(jì)算獲得定常爆轟狀態(tài)下炸藥與銅界面附近的爆轟波陣面的形態(tài)，如圖4 所示，爆轟定常傳播速度計(jì)算值為7 670 m/s（理論值為7 655 m/s）；爆轟波陣面邊緣角 α 的計(jì)算值為81.3°（理論值為78.8°），偏轉(zhuǎn)角 θ 的計(jì)算值為4.86°（理論值為4.75°），計(jì)算結(jié)果與理論結(jié)果吻合較好。同時(shí)也可以看出，爆轟波向金屬進(jìn)行的折射為正規(guī)折射，界面處沒有出現(xiàn)反射波。

圖3 銅約束爆轟波傳播的密度及壓力分布Fig. 3 Distribution of density and pressure in detonation flowfield under copper confinement

圖4 銅約束下爆轟波陣面形態(tài)Fig. 4 Detonation flowfield nearby explosives under copper confinement

3.3 平面爆轟波繞射問題

定常的平面爆轟波繞一個(gè)由金屬銅構(gòu)成的90°擴(kuò)張通道進(jìn)行傳播，計(jì)算域如圖5 所示。固體炸藥固相反應(yīng)物和氣相生成物均采用Jones-Wilkins-Lee 形式狀態(tài)方程[34]；爆轟化學(xué)反應(yīng)率采用Ignition-Growth 模型[34]；銅采用Mie-Grüneisen 形式狀態(tài)方程[18]。起爆區(qū)域設(shè)置有壓力為36.5 GPa 的高壓靜止氣相生成物，其余均為標(biāo)準(zhǔn)狀態(tài)；計(jì)算域入口設(shè)置為定值，上邊界設(shè)置為固壁條件，其余邊界設(shè)置為無反射邊界條件。

圖5 爆轟波繞射構(gòu)型圖Fig. 5 The configuration for the diffraction of detonation wave

計(jì)算獲得爆轟波在t= 4.25，5.00，5.75，6.50 μs 這4 個(gè)時(shí)刻的密度和壓力分布，如圖6 所示。

圖6 爆轟波繞射流場(chǎng)圖Fig. 6 The flowfield for the diffraction of detonation wave at various simulation times

從密度圖看，爆轟波在銅介質(zhì)中產(chǎn)生正規(guī)折射激波，該折射激波跨過臺(tái)階時(shí)會(huì)預(yù)壓炸藥；從壓力圖看，銅介質(zhì)中折射激波在向炸藥折射過程中產(chǎn)生稀疏波引起銅介質(zhì)出現(xiàn)一個(gè)低壓區(qū)。

4 結(jié) 論

本文中提出一種具有熱力學(xué)一致性的數(shù)值模擬固體炸藥爆轟與惰性介質(zhì)相互作用的擴(kuò)散界面模型，該模型基于混合網(wǎng)格內(nèi)各組分物質(zhì)處于力學(xué)平衡與熱學(xué)非平衡假設(shè)，推導(dǎo)獲得的模型流動(dòng)控制系統(tǒng)包括組分物質(zhì)的質(zhì)量守恒方程、混合物質(zhì)的動(dòng)量及總能量守恒方程，以及組分物質(zhì)的體積分?jǐn)?shù)演化方程和混合物質(zhì)的壓力演化方程。典型算例表明，數(shù)值計(jì)算結(jié)果與物理規(guī)律符合，并且在物質(zhì)界面附近不會(huì)出現(xiàn)物理量的非物理振蕩現(xiàn)象、適用于任意表達(dá)形式的物質(zhì)狀態(tài)方程以及任意數(shù)目的惰性介質(zhì)。