張瑞華，芮筱亭，趙宏立，王瓊林，靳建偉

（1.西安近代化學(xué)研究所，陜西 西安 710065；2.南京理工大學(xué)發(fā)射動力學(xué)研究所，江蘇 南京 210094）

隨著國防現(xiàn)代化的不斷推進(jìn)，對武器性能提出了大口徑、遠(yuǎn)射程、大威力的發(fā)展需求，高能量高安全是彈藥發(fā)展的永恒主題。目前，形成共識的發(fā)射裝藥引起膛炸的機(jī)理為：“低溫發(fā)射裝藥-擠壓-破碎-增面-增燃-增壓-膛炸”，即發(fā)射藥破碎極易導(dǎo)致膛炸事故[1-2]。因此，如何準(zhǔn)確描述相應(yīng)裝藥結(jié)構(gòu)下發(fā)射裝藥的擠壓破碎過程，成為高能量高安全發(fā)射藥研究的技術(shù)瓶頸問題。

學(xué)者們通過試驗或數(shù)值模擬對發(fā)射藥及發(fā)射裝藥擠壓破碎進(jìn)行了研究。Gazanas等[3]利用落錘和高速液壓伺服裝置研究了M30、JA2發(fā)射裝藥的壓縮和撞擊力學(xué)性能，用7.8 mL 小型密閉爆發(fā)器測量破碎藥粒的燃燒規(guī)律，得出破碎藥粒的燃面是未破碎藥粒的6 倍。陳言坤等[4]綜述了用于發(fā)射藥和發(fā)射裝藥力學(xué)性能研究的落錘沖擊試驗、空氣炮試驗、分離式霍普金森壓桿試驗和動態(tài)擠壓物理仿真試驗等，但是對發(fā)射裝藥的破碎模型和仿真還需進(jìn)一步研究。離散單元法（discrete element method, DEM）是Cundall[5]于1971年分析巖石等散體的力學(xué)行為時提出的，已被廣泛應(yīng)用于結(jié)構(gòu)工程、巖土工程、散體力學(xué)以及爆炸力學(xué)等領(lǐng)域。Jiang 等[6]建立了基于離散單元法的發(fā)射裝藥擠壓破碎程序，模擬了發(fā)射裝藥動態(tài)破碎過程并給出了三維動畫顯示；王燕等[7]對不同初始堆積發(fā)射裝藥開展了擠壓破碎實驗，并運(yùn)用離散單元法進(jìn)行了實驗驗證；以上研究采用彈簧-球單元破碎模型，為發(fā)射裝藥擠壓破碎模擬提供了研究手段。彈簧-球模型簡化了接觸力的計算和程序的編制，但是模型中彈性系數(shù)的確定較大程度上依賴于實驗值和經(jīng)驗值。EDEM軟件在計算顆粒破碎方面精確且高效，具有接觸模型參數(shù)由材料物性給定（如密度、彈性模量、泊松比等）、對實驗值和經(jīng)驗值依賴較少、易于確定等優(yōu)點。目前采用專業(yè)的離散單元法模擬軟件如EDEM[8]以及Hertz-Mindlin 黏結(jié)接觸破碎模型[9]等對發(fā)射裝藥擠壓破碎開展的模擬研究較少。

國內(nèi)外大部分膛炸事故發(fā)生于中大口徑火炮，其主裝藥普遍采用花邊十九孔發(fā)射藥，因此，探究花邊19孔發(fā)射藥的動態(tài)力學(xué)性能對于研究火炮發(fā)射裝藥的發(fā)射安全性至關(guān)重要。本文中，以硝胺花邊十九孔發(fā)射藥為研究對象，運(yùn)用EDEM 軟件對低溫下硝胺發(fā)射裝藥擠壓破碎和破碎發(fā)射裝藥密閉爆發(fā)器進(jìn)行模擬，結(jié)合硝胺發(fā)射裝藥動態(tài)擠壓破碎實驗和密閉爆發(fā)器實驗，對比模擬和實驗情況下的底部擠壓應(yīng)力、破碎發(fā)射裝藥燃?xì)馍梢?guī)律和起始動態(tài)活度比，實驗驗證數(shù)值模擬系統(tǒng)和參數(shù)選取的正確性。

1 發(fā)射裝藥擠壓破碎數(shù)值模擬系統(tǒng)

1.1 離散單元法

20 世紀(jì)90年代以來，各種離散單元法商業(yè)軟件相繼出現(xiàn)，如美國Itasca 公司開發(fā)的基于球形的PFC3D軟件、英國DEM-Solutions公司開發(fā)的商業(yè)化三維離散元EDEM 軟件等。相比之下，PFC軟件主要通過程序語言實現(xiàn)命令，使用較復(fù)雜。EDEM 軟件中的API 顆粒替換技術(shù)[8]對計算沖擊破碎更有優(yōu)勢。

API顆粒替換技術(shù)中的大顆粒是為了給替換的小顆粒群占據(jù)位置和提供中心坐標(biāo)，替換成功后該顆粒快速離開計算區(qū)域；小顆粒群以大顆粒模型的中心坐標(biāo)為相對位置進(jìn)行顆粒替換，由若干個同種材料的小球黏結(jié)形成需破碎顆粒的形狀，并在指定時刻形成黏結(jié)鍵，該顆粒具備2種接觸模型：Hertz-Mindlin無滑動接觸模型以及Hertz-Mindlin 黏結(jié)接觸模型，通過小球之間的黏接鍵斷裂來表征物料的破碎。

在發(fā)射裝藥擠壓破碎模擬過程中，如果離散顆粒太大，將難以表征真實發(fā)射藥破碎的小顆粒，若離散顆粒太小，則計算耗時過長。因此本文中將花邊十九孔發(fā)射藥離散為由95顆直徑相同的小球顆粒組成的離散單元模型，小球顆粒半徑為1.1 mm，輸出其坐標(biāo)作為API 顆粒替換技術(shù)中的小顆粒模型，替換前后的顆粒以及黏結(jié)模型如圖1所示，共形成286個黏結(jié)鍵。

圖1 發(fā)射藥顆粒替換黏結(jié)過程Fig.1 Replacement and forming bond process of propellant particle

1.1.1 Hertz-Mindlin 無滑動接觸模型

E D E M 軟件中默認(rèn)使用的接觸模型為Hertz-Mindlin 無滑動接觸模型，如圖2所示，該模型用于顆粒與顆粒之間、顆粒與裝置之間的接觸力計算，在顆粒計算方面精確且高效。

圖2 Hertz-Mindlin 無滑動接觸模型Fig.2 Hertz-Mindlin non-sliding contact model

法向力Fn為法向彈性力Fne和法向阻尼力Fnd之和；切向力Fτ為切向彈性力Fτe和切向阻尼力Fτd之和；庫倫摩擦fsFn限制切向力，其中fs為靜摩擦因數(shù)。具體表達(dá)式分別為：

其中：

式中：δn、δτ分別為法向和切向重疊量，vn、vτ分別為相對速度的法向和切向分量，kn、kτ分別為法向和切向剛度，β 為恢復(fù)系數(shù)e的表達(dá)式，Eb、Gb、Rb、mb分別為等效楊氏模量、等效剪切模量、等效半徑、等效質(zhì)量。具體表達(dá)式為：

式中：Ei、Ej分別為顆粒i和顆粒j的楊氏模量，μi、μj分別為顆粒i和顆粒j的泊松比，Ri、Rj分別為顆粒i和顆粒j的接觸球體半徑。

通過在接觸表面施加一個力矩來計算非常重要的滾動摩擦，表達(dá)式為：

式中：fr為滾動摩擦因數(shù)，Rt為接觸點到質(zhì)心的距離，ωt為物體在接觸點處的單位角速度矢量。

1.1.2 Hertz-Mindlin 黏結(jié)接觸模型

Hertz-Mindlin 黏結(jié)接觸模型如圖3所示。Hertz-Mindlin 黏結(jié)接觸模型運(yùn)用一個有限大的“膠粘劑”黏結(jié)在接觸范圍內(nèi)的顆粒，該黏結(jié)鍵可承受切向和法向位移，傳遞黏結(jié)顆粒之間的力和力矩，直到黏結(jié)鍵承受的應(yīng)力達(dá)到最大的臨界切向和法向應(yīng)力，即黏結(jié)鍵斷裂點，此后顆粒作為硬球互相接觸作用[9]。這個模型基于Potyondy等[9]的工作，特別適用于模擬混凝土和巖石結(jié)構(gòu)。該模型無法描述發(fā)射藥破碎后局部能量聚集、熱點形成過程。

圖3 Hertz-Mindlin 黏結(jié)接觸模型Fig.3 Hertz-Mindlin bonding contact model

顆粒在黏結(jié)生成時間tbond黏結(jié)在一起。黏結(jié)以后，顆粒上的黏結(jié)力Fˉn、Fˉt和力矩Mˉn、Mˉt設(shè)置為零，并在每個時間步通過下式逐步調(diào)整：

其中：

式中：Rbond為黏結(jié)半徑，Sˉn、Sˉτ分別為法向和切向黏結(jié)剛度（單位面積上的法向和切向剛度，單位為N/m3），Δt為時間步長，vn、vτ分別為顆粒法向和切向速度，ωn、ωτ分別為顆粒法向和切向角速度。

當(dāng)法向和切向剪切應(yīng)力超過某臨界值時，黏結(jié)鍵斷裂。定義法向應(yīng)力和切向應(yīng)力的最大值如下：

以上黏結(jié)力、黏結(jié)力矩是另外加到Hertz-Mindlin 無滑動接觸力中的。該模型只能用于顆粒與顆粒之間，且可以在顆粒沒有實際接觸時起作用，因此接觸半徑的設(shè)置應(yīng)大于實際半徑。

Hertz-Mindlin 黏結(jié)接觸模型中參數(shù)的計算公式如下[9]：

式中：Eˉb為黏結(jié)鍵彈性模量；Sn/Sτ為法向剛度與切向剛度的比值，與泊松比相關(guān)。

1.2 建立發(fā)射裝藥擠壓破碎數(shù)值模擬系統(tǒng)

1.2.1建立發(fā)射裝藥擠壓破碎力學(xué)模型

應(yīng)用EDEM 軟件，首先生成69個大球顆粒，在重力作用下進(jìn)行自然堆積，堆積穩(wěn)定后進(jìn)行API顆粒替換，顆粒替換的同時形成黏結(jié)鍵，替換后的發(fā)射裝藥模型進(jìn)一步進(jìn)行自然堆積，當(dāng)發(fā)射藥顆粒不再運(yùn)動時，則認(rèn)為發(fā)射裝藥堆積密實。堆積密實的發(fā)射裝藥離散單元力學(xué)模型如圖4所示，將圖5的燃燒室壓力曲線作為沖擊載荷加載到圖4中的上蓋板，進(jìn)行發(fā)射裝藥擠壓破碎數(shù)值模擬。

圖4 發(fā)射裝藥離散單元力學(xué)模型Fig.4 Discrete element mechanical model of propellant charge

圖5 上蓋板壓力時間曲線Fig.5 Pressure-timecurve of upper cover plate

1.2.2數(shù)值模擬參數(shù)選取

對保溫到低溫（–40℃）下的單顆硝胺發(fā)射藥進(jìn)行落錘沖擊試驗，試驗前將發(fā)射藥粒加工成軸向均勻、兩端面垂直軸向的試樣，試驗獲得的彈性模量為899.1 MPa、壓縮強(qiáng)度為98.5 MPa。

對圖1所示的單顆發(fā)射藥離散單元模型進(jìn)行落錘沖擊模擬，模擬邊界條件均與實驗相同，目的是尋找合適的Hertz-Mindlin黏結(jié)接觸模型參數(shù)，使得單顆發(fā)射藥落錘沖擊模擬結(jié)果與單顆發(fā)射藥落錘沖擊試驗結(jié)果符合較好[10]。依據(jù)Potyondy等[9]提出的黏結(jié)模型參數(shù)選取原則，令黏結(jié)鍵彈性模量Eˉb與實驗獲得的發(fā)射藥粒彈性模量相同，（Sn/Sτ）與泊松比相關(guān)，選取的黏結(jié)半徑為1 mm。則由式(22)、式(23)計算得到低溫下黏結(jié)模型的法向黏結(jié)剛度和切向黏結(jié)剛度，令其保持不變；法向臨界應(yīng)力和切向臨界應(yīng)力的比值也與泊松比有關(guān)，改變法向臨界應(yīng)力σc和切向臨界應(yīng)力τc的值，以符合最合適的黏結(jié)接觸模型參數(shù)，使得發(fā)射藥顆粒輸出的壓縮強(qiáng)度與落錘實驗獲得的一致。本文中選取6組Hertz-Mindlin黏結(jié)接觸模型參數(shù)，分別進(jìn)行單顆發(fā)射藥落錘沖擊模擬。

模擬參數(shù)和獲得的最大應(yīng)力如表1所示，其中序號4模擬獲得的最大應(yīng)力98.24 MPa與落錘沖擊試驗獲得的壓縮強(qiáng)度98.5 MPa最接近，選擇序號4的一組參數(shù)作為低溫下參數(shù)符合結(jié)果，用于進(jìn)一步進(jìn)行發(fā)射裝藥擠壓破碎數(shù)值模擬。發(fā)射藥顆粒的泊松比μ=0.3[11]，密度與真實發(fā)射藥密度1 650 kg/m3相同。

表1 數(shù)值模擬參數(shù)與結(jié)果Table 1 Model parameters and simulation results

2 發(fā)射裝藥擠壓破碎模擬實驗驗證

2.1 發(fā)射裝藥動態(tài)擠壓破碎實驗

發(fā)射裝藥動態(tài)擠壓破碎實驗裝置[2]如圖6所示。實驗原理為：在燃燒室中加入點火藥，燃燒生成高壓氣體推動活塞向下運(yùn)動，高速運(yùn)動的活塞沖擊模擬藥室中的發(fā)射裝藥試樣，用以模擬火炮發(fā)射過程中彈底發(fā)射裝藥的擠壓破碎過程以及發(fā)射藥床受到的擠壓應(yīng)力，PVDF應(yīng)力傳感器[12]測試發(fā)射裝藥底部擠壓應(yīng)力時間歷程。取圖7所示的69顆硝胺花邊十九孔發(fā)射裝藥為研究對象，實驗前將發(fā)射裝藥保溫到–40℃，對相同質(zhì)量的發(fā)射裝藥分別開展不同燃燒室壓力下的發(fā)射裝藥動態(tài)擠壓破碎實驗。

圖6 發(fā)射裝藥動態(tài)擠壓破碎實驗裝置及原理圖Fig.6 Experiment deviceand schematic diagram of the dynamic compression and fracture of propellant charge

圖7 實驗發(fā)射裝藥Fig.7 Experiment propellant charge

發(fā)射裝藥動態(tài)擠壓破碎實驗條件和結(jié)果如表2所示，實驗獲得的3發(fā)破碎發(fā)射裝藥如圖8所示，可以看出：在第2發(fā)和第3發(fā)實驗中發(fā)射裝藥出現(xiàn)了粉末狀破碎，燃燒時極易引起壓力異?，F(xiàn)象；且在第3發(fā)實驗中幾乎每顆發(fā)射藥粒都發(fā)生了破碎，破碎形狀各異、大小不一。

表2 實驗條件和結(jié)果Table2 Experiment conditionsand results

圖8 低溫下實驗獲得的破碎發(fā)射裝藥Fig.8 Fracture propellant charge obtained by experimentsat low temperature

2.2 發(fā)射裝藥擠壓破碎模擬結(jié)果

在圖4所示的發(fā)射裝藥離散單元力學(xué)模型的上蓋板分別加載圖8中對應(yīng)實驗的燃燒室壓力，分別進(jìn)行發(fā)射裝藥擠壓破碎模擬。模擬結(jié)果如表3所示，模擬獲得的3發(fā)破碎發(fā)射裝藥離散單元模型和黏結(jié)鍵模型分別如圖9、圖10所示。隨著模擬的進(jìn)行，發(fā)射裝藥發(fā)生破碎，發(fā)射藥粒破碎形狀各不相同，但是模擬發(fā)射裝藥的體積總量沒有發(fā)生改變，通過統(tǒng)計模擬結(jié)束時刻的總體黏結(jié)鍵斷裂情況來描述發(fā)射裝藥破碎程度。發(fā)射裝藥初始黏結(jié)鍵個數(shù)為286×69=19 734，3 發(fā)模擬黏結(jié)鍵斷裂比例分別達(dá)到6%、23%、38%，加載的燃燒室壓力越大，模擬結(jié)束時發(fā)射裝藥黏結(jié)鍵連接個數(shù)越少、擠壓應(yīng)力峰值越大，發(fā)射裝藥破碎越嚴(yán)重，與實驗現(xiàn)象一致。

表3 數(shù)值模擬結(jié)果Table 3 Simulation results

圖9 采用離散單元法模型獲得的擠壓破碎模擬結(jié)果Fig.9 Compression and fracture simulation results with thediscrete element model

圖10 采用黏結(jié)鍵模型得到的擠壓破碎模擬結(jié)果Fig.10 Compression and fracture simulation resultswith the bond model

2.3 發(fā)射裝藥擠壓應(yīng)力實驗驗證

不同燃燒室壓力下模擬與實驗獲得的發(fā)射裝藥底部擠壓應(yīng)力對比曲線如圖11所示。由圖11可知，加載的燃燒室壓力越大，擠壓應(yīng)力曲線脈寬越窄，曲線到達(dá)峰值前的上升速度越快；數(shù)值模擬與實驗獲得的擠壓應(yīng)力時間對比曲線存在少許誤差，這是由于發(fā)射藥的離散單元數(shù)值模擬模型不具有小孔，與真實花邊十九孔發(fā)射藥存在差異；且實驗與數(shù)值模擬的發(fā)射裝藥都是隨機(jī)堆積的，存在隨機(jī)誤差。在加載相同燃燒室壓力下，模擬與實驗獲得的底部擠壓應(yīng)力時間曲線總體吻合較好，但模擬與實驗獲得的破碎發(fā)射裝藥的破碎程度和燃?xì)馍梢?guī)律需進(jìn)一步驗證。

3 破碎發(fā)射裝藥燃?xì)馍梢?guī)律實驗驗證

運(yùn)用密閉爆發(fā)器實驗系統(tǒng)對圖8中的3發(fā)破碎發(fā)射裝藥分別進(jìn)行起始動態(tài)活度比實驗，測試獲得相應(yīng)破碎發(fā)射裝藥的燃?xì)馍梢?guī)律，依據(jù)起始動態(tài)活度比理論[1]處理得到起始動態(tài)活度比R0。起始動態(tài)活度比等于相同壓力下破碎發(fā)射裝藥與原未破碎發(fā)射裝藥被點燃時刻的面積比，即表征發(fā)射裝藥的破碎程度。將圖9所示的3發(fā)破碎發(fā)射裝藥分別進(jìn)行密閉爆發(fā)器實驗?zāi)M，運(yùn)用基于離散單元法的破碎發(fā)射藥燃燒函數(shù)[13]模擬獲得破碎發(fā)射裝藥的燃?xì)馍梢?guī)律及起始動態(tài)活度比。

數(shù)值模擬與實驗獲得的破碎發(fā)射裝藥密閉爆發(fā)器壓力時間對比曲線如圖12所示，起始動態(tài)活度比如表4所示?？梢钥闯?，對應(yīng)的壓力時間曲線吻合較好，起始動態(tài)活度比誤差最大值為6.05%，誤差較小。這說明建立的發(fā)射裝藥擠壓破碎模型和擬合研究獲得的黏結(jié)接觸模型參數(shù)具有較好的適用性。同時發(fā)現(xiàn)，第2發(fā)和第3發(fā)破碎后的壓力時間曲線前半段實驗的燃?xì)馍梢?guī)律大于數(shù)值模擬，而后半段則相反，且破碎程度越嚴(yán)重，曲線差別越大。這是由于數(shù)值模擬獲得的破碎發(fā)射裝藥的最小顆粒半徑大于實驗獲得的粉末狀破碎發(fā)射藥，這些細(xì)小的粉末會導(dǎo)致密閉爆發(fā)器實驗初始階段燃?xì)馍闪枯^大，實驗獲得的發(fā)射裝藥破碎程度越嚴(yán)重，細(xì)小的粉末越多，實驗與模擬獲得的壓力時間曲線誤差越大。當(dāng)實驗獲得的細(xì)小粉末燃燒完后，剩下的大塊實驗破碎發(fā)射裝藥的燃燒面積將小于模擬獲得的破碎發(fā)射裝藥，因而實驗的燃?xì)馍梢?guī)律后半段曲線小于數(shù)值模擬。當(dāng)R0<2.0時[2]，壓力時間曲線吻合較好，說明本文的離散元小球顆粒半徑滿足計算要求。

表4 數(shù)值模擬與實驗起始動態(tài)活度比對比Table 4 Comparision of initial dynamic vivacity ratios in simulation and experiment

圖12 密閉爆發(fā)器數(shù)值模擬與實驗壓力時間對比曲線Fig.12 Comparison of simulated and experiment pressure-time curves in the closed bomb

4 結(jié) 論

（1）加載不同燃燒室壓力下發(fā)射裝藥擠壓破碎模擬獲得的黏結(jié)鍵斷裂比例分別達(dá)到6%、23%、38%，與發(fā)射裝藥動態(tài)擠壓破碎實驗現(xiàn)象定性一致，且模擬與實驗獲得的底部擠壓應(yīng)力時間曲線總體吻合較好。

（2）在不同破碎程度的發(fā)射裝藥密閉爆發(fā)器模擬與實驗中，對應(yīng)破碎發(fā)射裝藥的壓力時間曲線吻合較好，起始動態(tài)活度比誤差最大值為6.05%，誤差較小，驗證了破碎發(fā)射裝藥燃燒規(guī)律的一致性。當(dāng)數(shù)值模擬破碎發(fā)射裝藥的起始動態(tài)活度比R0<2.0時，離散元小球顆粒半徑滿足計算要求。

（3）對比實驗和數(shù)值模擬情況下的底部擠壓應(yīng)力、破碎發(fā)射裝藥燃?xì)馍梢?guī)律和起始動態(tài)活度比，實驗驗證了數(shù)值模擬系統(tǒng)的正確性和黏結(jié)接觸模型參數(shù)擬合的合理性，研究結(jié)果為再現(xiàn)由發(fā)射裝藥破碎導(dǎo)致的膛壓異常奠定了理論和模擬基礎(chǔ)。