基于離散單元法的發射裝藥擠壓破碎模擬實驗*

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

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

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

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

國內外大部分膛炸事故發生于中大口徑火炮，其主裝藥普遍采用花邊十九孔發射藥，因此，探究花邊19孔發射藥的動態力學性能對于研究火炮發射裝藥的發射安全性至關重要。本文中，以硝胺花邊十九孔發射藥為研究對象，運用EDEM 軟件對低溫下硝胺發射裝藥擠壓破碎和破碎發射裝藥密閉爆發器進行模擬，結合硝胺發射裝藥動態擠壓破碎實驗和密閉爆發器實驗，對比模擬和實驗情況下的底部擠壓應力、破碎發射裝藥燃氣生成規律和起始動態活度比，實驗驗證數值模擬系統和參數選取的正確性。

1 發射裝藥擠壓破碎數值模擬系統

1.1 離散單元法

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

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

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

圖1 發射藥顆粒替換黏結過程Fig.1 Replacement and forming bond process of propellant particle

1.1.1 Hertz-Mindlin 無滑動接觸模型

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

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

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

其中：

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

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

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

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

1.1.2 Hertz-Mindlin 黏結接觸模型

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

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

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

其中：

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

當法向和切向剪切應力超過某臨界值時，黏結鍵斷裂。定義法向應力和切向應力的最大值如下：

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

Hertz-Mindlin 黏結接觸模型中參數的計算公式如下[9]：

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

1.2 建立發射裝藥擠壓破碎數值模擬系統

1.2.1建立發射裝藥擠壓破碎力學模型

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

圖4 發射裝藥離散單元力學模型Fig.4 Discrete element mechanical model of propellant charge

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

1.2.2數值模擬參數選取

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

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

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

表1 數值模擬參數與結果Table 1 Model parameters and simulation results

2 發射裝藥擠壓破碎模擬實驗驗證

2.1 發射裝藥動態擠壓破碎實驗

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

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

圖7 實驗發射裝藥Fig.7 Experiment propellant charge

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

表2 實驗條件和結果Table2 Experiment conditionsand results

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

2.2 發射裝藥擠壓破碎模擬結果

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

表3 數值模擬結果Table 3 Simulation results

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

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

2.3 發射裝藥擠壓應力實驗驗證

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

3 破碎發射裝藥燃氣生成規律實驗驗證

運用密閉爆發器實驗系統對圖8中的3發破碎發射裝藥分別進行起始動態活度比實驗，測試獲得相應破碎發射裝藥的燃氣生成規律，依據起始動態活度比理論[1]處理得到起始動態活度比R0。起始動態活度比等于相同壓力下破碎發射裝藥與原未破碎發射裝藥被點燃時刻的面積比，即表征發射裝藥的破碎程度。將圖9所示的3發破碎發射裝藥分別進行密閉爆發器實驗模擬，運用基于離散單元法的破碎發射藥燃燒函數[13]模擬獲得破碎發射裝藥的燃氣生成規律及起始動態活度比。

數值模擬與實驗獲得的破碎發射裝藥密閉爆發器壓力時間對比曲線如圖12所示，起始動態活度比如表4所示。可以看出，對應的壓力時間曲線吻合較好，起始動態活度比誤差最大值為6.05%，誤差較小。這說明建立的發射裝藥擠壓破碎模型和擬合研究獲得的黏結接觸模型參數具有較好的適用性。同時發現，第2發和第3發破碎后的壓力時間曲線前半段實驗的燃氣生成規律大于數值模擬，而后半段則相反，且破碎程度越嚴重，曲線差別越大。這是由于數值模擬獲得的破碎發射裝藥的最小顆粒半徑大于實驗獲得的粉末狀破碎發射藥，這些細小的粉末會導致密閉爆發器實驗初始階段燃氣生成量較大，實驗獲得的發射裝藥破碎程度越嚴重，細小的粉末越多，實驗與模擬獲得的壓力時間曲線誤差越大。當實驗獲得的細小粉末燃燒完后，剩下的大塊實驗破碎發射裝藥的燃燒面積將小于模擬獲得的破碎發射裝藥，因而實驗的燃氣生成規律后半段曲線小于數值模擬。當R0<2.0時[2]，壓力時間曲線吻合較好，說明本文的離散元小球顆粒半徑滿足計算要求。

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

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

4 結 論

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

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

（3）對比實驗和數值模擬情況下的底部擠壓應力、破碎發射裝藥燃氣生成規律和起始動態活度比，實驗驗證了數值模擬系統的正確性和黏結接觸模型參數擬合的合理性，研究結果為再現由發射裝藥破碎導致的膛壓異常奠定了理論和模擬基礎。