基于SPH的橢圓雙極線性聚能藥包控制爆破數(shù)值模擬研究

吳 波，韋 漢，徐世祥，許 杰，黃宗輝

(1.廣西大學(xué) 土木建筑工程學(xué)院，廣西 南寧 530004；2.教育部工程防災(zāi)與結(jié)構(gòu)安全重點(diǎn)實(shí)驗(yàn)室，廣西 南寧 530004；3.廣西防災(zāi)減災(zāi)與工程安全重點(diǎn)實(shí)驗(yàn)室，廣西 南寧 530004)

0 引 言

聚能裝藥廣泛應(yīng)用于軍事領(lǐng)域，比如爆炸成型彈丸、破甲彈設(shè)計(jì)等，此外聚能裝藥也被廣泛應(yīng)用于工業(yè)生產(chǎn)，尤其是在實(shí)現(xiàn)巖石的定向斷裂控制爆破技術(shù)中。線性聚能裝藥是通過擠壓聚能藥型罩形成聚能射流的，從而達(dá)到定向斷裂控制爆破的效果。近年來，人們對線性聚能裝藥爆破技術(shù)進(jìn)行了廣泛的研究，由于爆破的瞬時性，實(shí)驗(yàn)和理論難以捕捉爆炸的整個過程，而數(shù)值模擬可以再現(xiàn)爆炸的過程，使之成為續(xù)實(shí)驗(yàn)和理論之后的一種研究手段。MOLINARI[1]采用有限元來研究聚能裝藥中銅射流的形成、破碎和侵徹過程。王成等[2]通過ALE研究不同起爆方式對線性聚能射流形成的影響進(jìn)行分析。張先鋒等[3]也通過ALE對3種不同典型線性裝藥的聚能射流形成及侵徹過程進(jìn)行了數(shù)值模擬。李必紅[4]從理論、實(shí)驗(yàn)與數(shù)值模擬研究橢圓雙極線性聚能藥包的爆破機(jī)理，提出了有效裝藥邊界函數(shù)并通過理論推導(dǎo)確定。

以上實(shí)現(xiàn)的數(shù)值分析都是基于網(wǎng)格的拉格朗日方法、歐拉方法或者兩者的結(jié)合，然而實(shí)現(xiàn)的拉格朗日方法由于大變形的存在往往使得計(jì)算終止，而歐拉方法或ALE方法對物質(zhì)邊界的捕捉存在一定的困難。由于無網(wǎng)格法具有處理大變形和破壞問題的特點(diǎn)，使得其成為處理極端變形條件下有限元無法比擬的優(yōu)勢，尤其是無網(wǎng)格SPH方法在線性聚能裝藥中的應(yīng)用。SPH方法最早于1977年由LUCY等[5]提出并將其天體物理學(xué)，隨后廣泛應(yīng)用于各個領(lǐng)域的研究，SWEGLE等[6]首次將其應(yīng)用于水下爆炸數(shù)值研究，實(shí)現(xiàn)SPH在爆炸領(lǐng)域的應(yīng)用，LIU等[7-8]等基于SPH對無藥型罩聚能裝藥爆轟過程及聚能射流形成過程進(jìn)行研究。YANG等[9]基于二維模型對線性聚能射流進(jìn)行SPH模擬研究，并對比分析不同數(shù)值分析方法對聚能射流的影響。李磊等[10-11]通過SPH實(shí)現(xiàn)聚能裝藥的三維數(shù)值模擬，并將SPH方法與有限元耦合實(shí)現(xiàn)聚能射流對鋼板的侵徹效果。楊建輝[12]通過SNSYS/LS-DYNA對切縫聚能管參數(shù)進(jìn)行優(yōu)化，分析不同參數(shù)下對巖石應(yīng)力狀態(tài)的影響。陳帥志等[13]基于LS-DYNA分析切縫藥包與普通藥包爆炸在巖石的傳播特性。岳中文等[14]結(jié)合數(shù)值模擬及試驗(yàn)研究炮孔間距對切縫藥包爆破引起的裂紋擴(kuò)展的影響。張盛等[15]將橢圓雙極線性聚能藥包在巷道中的應(yīng)用，并結(jié)合現(xiàn)場試驗(yàn)對參數(shù)進(jìn)行優(yōu)化研究。

綜上所述，目前對聚能裝藥形成數(shù)值分析主要基于有限元分析，而有限元模型無法動態(tài)捕捉爆生氣體跟聚能罩的動態(tài)演化規(guī)律，或?qū)Ь勰懿鄣木匦嗡幇蚯锌p藥包進(jìn)行聚能分析，而基于無網(wǎng)格實(shí)現(xiàn)橢圓雙極線性聚能裝藥結(jié)構(gòu)的研究較少。鑒于SPH方法在處理爆炸這種極端變形等復(fù)雜條件下具有有限元無法比擬的優(yōu)勢。筆者基于SPH方法研究聚能裝藥的形成，并與已有試驗(yàn)及模擬做對比，驗(yàn)證本文實(shí)現(xiàn)的SPH方法的有效性；然后基于SPH研究橢圓雙極線性聚能藥包爆破機(jī)理并分析外殼及其厚度、藥型罩及其厚度對橢圓雙極線性聚能藥包實(shí)現(xiàn)聚能射流及其對巖石定向斷裂控制爆破技術(shù)的影響。

1 SPH控制方程理論

SPH方法的核心思想就是通過構(gòu)造一個核函數(shù)對粒子的核估計(jì)來計(jì)算梯度的相關(guān)量，不需要求解偏微分方程，而是將偏微分控制方程轉(zhuǎn)化為積分形式，偏微分方程的轉(zhuǎn)化包括函數(shù)核近似和粒子近似兩步，函數(shù)核〈f(x)〉[11]近似如下：

(1)

式中，W(x-x′，h)為光滑函數(shù)或核函數(shù);Ω為積分域；x為空間坐標(biāo)；x′為x緊支域內(nèi)某一點(diǎn)的坐標(biāo)；h為光滑長度。

對于連續(xù)形式的核近似式，函數(shù)f(x)的粒子近似可對粒子i緊支域內(nèi)所有粒子j進(jìn)行加權(quán)求和[16]，即

(2)

式中：mj，ρj分別為粒子i的臨近粒子j的質(zhì)量和密度。

光滑函數(shù)需要滿足歸一化、δ函數(shù)性質(zhì)和緊支性3個基本條件，常見的光滑函數(shù)有3次B樣條函數(shù)、5次B樣條函數(shù)、高斯核函數(shù)、指數(shù)核函數(shù)等，本文采用的是最常用的3次B樣條核函數(shù)[11]，形式如下：

(3)

在粒子近似過程中，如果采用固定光滑長度，粒子集中時將影響計(jì)算效率，同時對于具有強(qiáng)度的材料可能造成張力不穩(wěn)定的現(xiàn)象，粒子散射時會出現(xiàn)粒子緊支域內(nèi)粒子數(shù)減少造成數(shù)值不穩(wěn)定，因此本文采用變光滑長度法，其光滑長度的時間積分[17]格式如下：

(4)

式中:d為空間維數(shù)，本文取d=2；v為粒子速度；h(t)為粒子光滑長度；div為輔助函數(shù)。

2 模型驗(yàn)證

2.1 試驗(yàn)描述問題

通過選取與文獻(xiàn)[18]相同的幾何參數(shù)模型，模型參數(shù)分別為藥型罩頂角W=80，厚度D=1.35 mm，裝藥長度L=25.4 mm，裝藥高度H=28.45 mm，模型如圖1所示。通過LS-DYNA刪除Lagrange單元信息再添加SPH相關(guān)屬性，實(shí)現(xiàn)的SPH二維分析模型如圖2所示。試驗(yàn)?zāi)Ｐ蛯?shí)現(xiàn)的SPH粒子間距為0.008～0.010 cm，炸藥和藥型罩的最小最大光滑長度分別為初始光滑長度的0.2倍和2倍，由于模型對稱，計(jì)算時取1/2模型進(jìn)行計(jì)算。

圖1 數(shù)值模型Fig.1 Numerical model

圖2 SPH離散模型Fig.2 SPH discrete model

炸藥采用B炸藥進(jìn)行分析模擬，并采用MAT_HIGH_EXPLOSIVE_BURN模型描述材料本構(gòu)，狀態(tài)方程為JWL:

(5)

式中：p為爆轟壓力；E為比內(nèi)能；V為相對比容；A、B、R1、R2、ω為狀態(tài)方程系數(shù)，炸藥材料參數(shù)及狀態(tài)方程參數(shù)見表1[18]。

表1 B炸藥材料參數(shù)Table 1 B explosive material parameters

金屬藥型罩同樣采用文獻(xiàn)[18]的材料1006號鋼，本構(gòu)模型為MAT_JOHNSON_COOK模型，狀態(tài)方程為EOS_GRUNEISEN：

(6)

表2 1006號鋼材料參數(shù)Table 2 Material parameters of No.1006 steel

2.2 模擬結(jié)果分析

圖3為聚能射流形成過程中不同時刻射流形態(tài)及速度分布。由圖3可知，基于LS-DYNA實(shí)現(xiàn)的SPH線性聚能射流與文獻(xiàn)[18]基于高斯光滑核函數(shù)的自編程序?qū)崿F(xiàn)的聚能射流在不同時刻基本一致，且都在t=4 μs 時爆轟基本完成，射流初步形成，如圖3a所示。

圖3 聚能射流形成不同時刻速度云圖Fig.3 Cloud chart of velocity at different time of jet formation

圖4為t=6.5μs時聚能射流的速度云圖，此刻射流頭部平均最大速度達(dá)到3 507 m/s，與文獻(xiàn)[19]電磁線圈測得射流頭部最大平均速度值3 520 m/s基本一致。文獻(xiàn)[18]實(shí)現(xiàn)的SPH方法獲取的射流頭部峰值速度為3 294 m/s，文獻(xiàn)[19]通過X光技術(shù)測得的射流頭部峰值速度為3 300～3 500 m/s，本文實(shí)現(xiàn)的SPH聚能射流形成穩(wěn)定之后射流頭部平均速度為3 180 m/s，說明本文實(shí)現(xiàn)的SPH方法和文獻(xiàn)[18]及文獻(xiàn)[19]基本一致，表明本文基于DYNA實(shí)現(xiàn)的SPH方法的有效性。

圖4 t=6.5 μs時刻聚能射流速度云圖Fig.4 Cloud chart of jet velocity at t=6.5 μs

3 橢圓雙極線性聚能藥包爆破數(shù)值分析

基于前文實(shí)現(xiàn)的SPH方法，可知SPH可以處理極端變形復(fù)雜條件下的破壞問題，為此對橢圓雙極線性聚能藥包爆破機(jī)理及其影響因素(外殼及其厚度、藥型罩及其厚度)進(jìn)行分析。

3.1 橢圓雙極線性聚能藥包形成

3.1.1 模型的建立

圖5為橢圓雙極線性聚能藥包四分之一模型，幾何參數(shù)如圖5a所示,長短軸比為11∶15，離散模型如圖5b所示，為考慮薄外殼與藥型罩具有足夠的臨近粒子數(shù)，減小外殼與藥型罩的粒子間距，考慮模型尺寸較小，將炸藥粒子間距設(shè)置為外殼及藥型罩相同的粒子間距0.008～0.010 cm。模型材料參數(shù)取值與上節(jié)材料參數(shù)一致，外殼選取與藥型罩相同的材料及本構(gòu)方程和狀態(tài)方程。整個模型采用變光滑長度的計(jì)算方法，炸藥、外殼和藥型罩的最小最大光滑長度分別為初始光滑長度的0.2倍和2倍，計(jì)算時采用對稱邊界處理。

圖5 四分之一聚能藥包模型Fig.5 Quarter model of shaped-charge

3.1.2 聚能藥包作用機(jī)理分析

圖6為聚能射流形成不同時刻速度分布云圖。t=0.8 μs時爆轟波先傳到藥型罩錐角頂端，開始擠壓藥型罩形成射流，隨后到達(dá)90°非聚能方向的外殼，開始沖擊作用于90°非聚能方向的外殼，使之開始發(fā)生壓縮變形與向外運(yùn)動。t=4.0 μs時射流初步形成，射流頭部平均最大速度達(dá)到3 224 m/s，外殼不斷發(fā)生塑性拉伸及向外運(yùn)動，同時外殼與藥型罩連接部位受拉力作用不斷拉伸縮頸。在t=14.0 μs時連接部位首先發(fā)生斷裂，爆生氣體從該位置溢出，此時射流不斷被拉伸，并與存在速度梯度的杵體部分逐步分離，由于射流不再受杵體、爆生氣體和外殼的影響，射流頭部平均速度基本保持不變。但最大速度并未發(fā)生在射流頭部，而是在藥型罩與外殼連接處的炸藥粒子，由于射流形成穩(wěn)定后不再受其他因素的影響，而接觸部位炸藥粒子在脫離藥型罩與外殼的約束作用之后受到內(nèi)部相鄰炸藥粒子的直接加速作用。t=20.0 μs時，脫離外殼與藥型罩的炸藥粒子不斷加速追趕射流頭部，此時炸藥粒子最大速度達(dá)到3 772 m/s，而射流頭部平均最大速度為2 985 m/s，同時在90°非聚能方向的外殼不斷發(fā)生縮頸甚至斷裂。

圖6 聚能射流形成不同時刻速度分布云圖Fig.6 Cloud chart of velocity distribution at different time of jet formation

圖7為不同位置粒子的時間歷程曲線，A′為射流頭部粒子編號為302614，B′為藥型罩與外殼接觸部位內(nèi)側(cè)的炸藥粒子編號為292551，C′為90°非聚能方向外殼粒子編號為302831，D′為90°非聚能方向外殼內(nèi)側(cè)炸藥粒子編號為290904。由圖可知，聚能射流在t=4 μs時射流基本形成，與圖6b表述相吻合。t=14.0 μs時爆生氣體溢出藥型罩與外殼連接位置，最終在t=15.0 μs后粒子以3 440 m/s的速度穩(wěn)定向前運(yùn)動。90°非聚能方向外殼粒子與外殼內(nèi)部炸藥粒子的速度增加的總體趨勢一致，但由于炸藥粒子受到爆轟波的往復(fù)作用，速度波動幅度相對外殼粒子要大。此外，射流頭部速度發(fā)展較快，其他粒子速度存在明顯的時間滯后現(xiàn)象。

圖7 不同位置處粒子速度隨時間變化曲線Fig.7 Curve of particle velocity with time at different positions

3.2 聚能藥包參數(shù)優(yōu)化

3.2.1 問題描述

不同幾何參數(shù)的橢圓雙極線性聚能藥包四分之一模型如圖8所示，在保證長短軸比、聚能槽角度和裝藥量一致的情況下，各個模型采用相同的材料參數(shù)和粒子間距(同上節(jié))。粒子的最大、最小光滑長度分別為0.2倍和2倍初始光滑長度，為保證初始時刻每個粒子的支持域內(nèi)具有相同的質(zhì)量，每個粒子初始臨近粒子數(shù)設(shè)置為600。

圖8 不同幾何參數(shù)的聚能藥包四分之一模型Fig.8 Quarter models of shaped-charge cartridge with different geometric parameters

3.2.2 數(shù)值模擬結(jié)果及分析

圖9為不同幾何參數(shù)聚能藥包在t=14.0 μs時速度云圖，由圖可知不同幾何參數(shù)的聚能藥包實(shí)現(xiàn)的定向效果不同。圖9a為2 mm厚的藥型罩1 mm厚外殼的聚能藥包速度分布圖，由圖可知除射流頭部速度比較大之外，外殼向外運(yùn)動的速度也較大，整體速度大于杵體速度，主要是由于外殼比較薄，爆生氣體對外殼做的功大部分轉(zhuǎn)化為外殼動能。對比圖9b無外殼情況可知，圖9a形成的射流速度較大且藥型罩大部分的粒子形成射流，而圖9b藥型罩大部分粒子用于形成杵體，表明外殼對聚能射流的形成產(chǎn)生較大影響。圖9c為1 mm厚藥型罩和2 mm厚外殼的聚能藥包速度分布圖，對比圖6c的2 mm厚藥型罩情況可知，兩者在同一時刻形成的射流形態(tài)相似，但前者形成的射流頭部速度更大，射流頭部平均最大速度達(dá)到3 670 m/s，更有利于巖石的侵徹作用。圖9d由于無藥型罩而無法形成射流，圖中只給出爆生氣體從外殼切縫處溢出形態(tài)，該形態(tài)與申濤等[20]基于ALE算法模擬切縫藥包爆破相似，由圖可見粒子散射比較明顯，實(shí)現(xiàn)的定向效果不顯著。圖9e為1 mm厚藥型罩和1 mm厚外殼的聚能藥包速度分布圖，對比圖6c可知同樣前者形成的射流更明顯，速度更大，表明藥型罩和外殼厚度同時減小時有利于聚能射流的形成。圖9f為無外殼無藥型罩藥包的爆生氣體速度分布圖，由圖可見，并未形成明顯射流，可見外殼與藥型罩對聚能射流的形成起決定作用。

圖10為不同幾何參數(shù)聚能藥包射流頭部速度隨時間變化曲線，通過比較圖5a、圖8a和圖8b所對應(yīng)的曲線可知，隨著外殼厚度不斷減小，最終射流頭部速度不斷減小，但三者的速度變化趨勢基本一致，再次表明外殼對橢圓雙極線性聚能藥包射流頭部的形成有一定的影響，外殼越厚更有利于穩(wěn)定爆轟作用，從而有利于射流的形成，射流速度越大；比較圖5a、圖8c所對應(yīng)的曲線，可知當(dāng)藥型罩為1 mm，外殼2 mm時形成射流頭部速度最大，最大值為3 700 m/s，而圖8d由于無藥型罩而無法形成聚能射流，對應(yīng)圖10的曲線只不過是聚能方向外殼邊緣處粒子的速度時間歷程曲線，表明藥型罩對聚能射流的形成起決定性作用；對比圖5a和圖8e所對應(yīng)的曲線圖可知，隨著藥型罩和外殼厚度一致減小，聚能射流頭部速度增大，但對比圖9f可知當(dāng)兩者減小為0時不形成射流。對比圖8d和圖8f所對應(yīng)的曲線，由于兩者均沒有藥型罩，只能借助炸藥粒子分析，由圖可見，兩者速度曲線重合，即最終穩(wěn)定速度完全一致，表明當(dāng)無藥型罩時并沒有明顯的聚能效應(yīng)。

圖9 不同幾何參數(shù)聚能藥包在t=14.0 μs時速度云圖Fig.9 Velocity nephogram of shaped charge with different geometric parameters at t=14.0 μs

圖10 不同幾何參數(shù)聚能藥包射流頭部速度時間歷程曲線Fig.10 Time history curve of jet head velocity of shaped-charge with different geometric parameters

4 結(jié) 論

1)橢圓雙極線性聚能藥包通過壓垮藥型罩形成聚能射流，聚能射流形成的同時爆生氣體不斷撞擊外殼向非聚能方向運(yùn)動，變形到達(dá)一定程度之后外殼與藥型罩連接部位斷開，爆生氣體從斷開部位溢出，溢出的爆生氣體速度大于射流頭部速度。

2)隨著外殼厚度不斷增大，射流頭部速度也在增大，且不同外殼厚度的聚能藥包射流頭部速度隨時間變化趨勢趨于一致，表明外殼對橢圓雙極線性聚能藥包射流頭部速度有一定影響，外殼越厚越利于穩(wěn)定爆轟，從而利于射流的形成，射流速度越大。

3)當(dāng)外殼厚度一定時，藥型罩厚度不斷減小，爆生氣體用于壓垮藥型罩做的功越少，轉(zhuǎn)換為聚能射流動能越多，聚能射流頭部速度越大，但當(dāng)藥型罩厚度為0時未能形成明顯的聚能效應(yīng)。