基于PBM算法的沖擊波載荷數值模擬研究

張紹彥，王顯會，彭 兵，孫魁遠，沙康康，謝 淵

(南京理工大學 機械工程學院，南京 210094)

1 引言

近年來，軍用戰術車輛在戰場上面對來自地雷的威脅日益增加，爆炸產生的沖擊波載荷會對裝甲車輛和乘員造成嚴重傷害[1]，開展裝甲車輛抗爆分析、設計與防護日益成為關注的熱點問題。爆炸試驗具有高成本、耗時等特點，因此數值模擬方法被廣泛應用于爆炸載荷的研究。

目前，多材料任意拉格朗日-歐拉(MM-ALE)算法結合流固耦合(FSI)算法模擬爆炸沖擊波傳播是公認的模擬爆炸的方法，國內外學者也對ALE算法可靠性進行了大量研究[2-4]。該方法兼具了Lagrange算法與Euler算法的優點，對于固體材料，允許節點和材料共同運動、變形；而對于流體材料，網格節點固定不變[5]。該方法能夠處理網格畸變問題，可以更加精確地描述結構之間接觸滑移面，但是在地雷爆炸情況下，高能炸藥通常埋藏在土壤中，土壤與結構物相互作用的影響不容忽視[6-8]。為了提高ALE算法模擬結果的準確性，一般需要建立大量的空氣和土壤單元，會造成計算效率低下，同時，拉格朗日-歐拉算法中的接觸問題極易造成計算的不穩定。

因此，國內外學者開發了基于粒子的算法來提高計算效率與穩定性[9-11]，如光滑粒子法(SPH)、粒子爆炸法(PBM)和離散元法(DEM)等。李利莎等[12]運用ALE和SPH兩種算法對炸藥在鋼筋混凝土板表面爆炸進行了模擬分析，對比分析了2種算法的優缺點，發現SPH算法可以更為逼真地模擬爆轟效果并且可以有效地防止能量泄漏。L.Olovsson等[13]最先將PBM算法應用到爆炸仿真中。隨后，Teng等[14]對PBM算法進行了改進，并在LS-DYNA中實現了該算法，模擬了一系列固支鋼板近距離爆炸，同時將仿真結果與實驗結果進行比較，發現該方法具有高精度，較ALE方法具有更好的計算經濟性。

本文以固支鋼板空爆試驗為對標對象，采集爆炸沖擊下鋼板中心最大變形數據，分別開展PBM算法、ALE算法數值仿真模擬，驗證PBM算法計算結果的準確性以及計算的高效性。同時對PBM算法進行進一步研究，探究炸藥粒子數量、空氣粒子數量對計算結果的影響，為尋求模擬爆炸載荷的最優算法提供理論依據。

2 粒子爆炸法

粒子爆炸法(PBM)的本質是對基于分子動力學理論[15]的粒子法的改進，在該算法中，當粒子發生碰撞時，平動動能和轉動動能之間的平衡不再成立，但對整個系統來說，系統達到平衡態時，平動動能和轉動動能之間滿足統計學上的平衡規律。此外，該算法為了更好地描述氣體在極端壓力和溫度下的行為，考慮了協體積效應，使得爆炸產物在高溫高壓下的氣體行為更加精確。

2.1 分子動力學理論

分子動力學理論起源于1738年Daniel Bernoulli提出的一種理論，即氣體壓強是大量氣體分子與器壁碰撞的結果。以動力學理論為出發點，James Clerk Maxwell得出了一個非常簡潔的表達式，表示熱平衡時的分子速度分布：

(1)

式中：f(v)是分子速率的概率密度函數，其中：v為分子速率；M為摩爾質量；R為通用氣體常數；T為溫度。

分子動力學理論是對氣體分子及其相互作用(在微觀水平上)的研究，該理論揭示了氣體宏觀熱學性質和過程的微觀本質，推導出了理想氣體定律(宏觀關系)，給出了宏觀量與微觀量平均值的關系。該理論基于以下假設：① 分子之間的平均距離遠大于分子自身的尺寸大小；② 存在熱力學平衡，即分子處于隨機運動；③ 分子可以看作質點，運動時遵守牛頓運動定律；④ 分子-分子和分子-結構相互作用是完全彈性的碰撞；⑤ 除碰撞的瞬間外，分子間的相互作用力可忽略不計。圖1 為粒子-粒子相互作用示意圖[16]。

圖1 粒子-粒子相互作用示意圖

2.2 粒子爆炸法在LS-DYNA中的應用

每克爆炸物根據分子量的不同可以達到1022～1023分子的數量級，由于無法對爆炸場中的每一個分子進行獨立建模，為了提升仿真過程中的計算效率，對分子動力學理論進行以下假設：

1)粒子是剛性的，并且呈球形；

2)每個粒子代表許多分子，分子數通常為1015～1020，具體取決于應用；

3)對于每個單獨的粒子，平動動能Wt和轉動動能WS之間的平衡關系由γ=Cp/Cv確定，其中Cp，Cv分別表示在恒定壓力和恒定體積下的熱容。平衡假定為:

(2)

Ws是每個粒子的一個總的標量變量，因此，粒子沒有為轉動/振動分配任何自由度。當粒子撞擊移動的結構時，Wt將發生變化，而Ws則假定不受影響，因此，Ws和Wt之間的比值將發生改變。當2個粒子碰撞時，它們之間的能量將重新分配以維持Ws、Wt之間比值的恒定，以確保總能量和動量的平衡。粒子的半徑可以動態調整以獲得合理的平均自由程，由于粒子質量是恒定的，半徑過小會導致壓力波不能傳遞，半徑過大會使氣體作用偏離理想氣體定律。此外，隨著粒子半徑的增大，粒子之間接觸計算也會使仿真的成本增加。

3 固支鋼板抗爆轟試驗

3.1 試驗布置

為驗證PBM算法仿真精度，開展了固支鋼板在TNT空爆條件下的試驗研究，試驗臺架如圖2所示，固支鋼板長寬尺寸1 500 mm×1 500 mm，厚度4.5 mm，周邊等距離開孔，試驗臺架上、下方框與鋼板對應位置開孔，通過螺栓將鋼板完全固定，臺架上方框厚度為10 mm，下方框厚度為25 mm，在鋼板正下方布置變形梳測試鋼板在沖擊波載荷作用下的中心點最大撓度，變形梳由梳齒和梳背兩部分組成，梳齒長度呈階梯狀分布，梳齒材質為易產生塑性變形的鋁合金，變形梳如圖3所示，變形梳位置如圖4所示。為增強試驗臺架剛度，內部設置多條加強板，以減小試驗過程中臺架發生大的破壞影響試驗結果。試驗中爆炸物為2 kg圓柱形TNT藥柱，高徑比為0.92，藥柱直徑119.2 mm，高109.8 mm，爆炸品位于靶板中心的正上方，藥柱下表面距離鋼板上表面300 mm，采用中心起爆方式，試驗實物布置如圖5所示。

圖2 臺架示意圖

圖3 變形梳示意圖

圖4 變形梳位置示意圖

圖5 試驗實物布置圖

3.2 試驗結果

爆炸試驗后臺架基座無明顯損壞，四周固定螺栓未脫落，鋼板的試驗及仿真結果如圖6所示，回收置于鋼板底部的變形梳如圖7所示，通過對比變形疏變形前與變形后梳齒的變形測量得鋼板中心點最大撓度為60.1 mm。

圖6 鋼板變形示意圖

圖7 試驗后變形梳圖

4 有限元爆炸仿真研究

采用LS-DYNA非線性動力學仿真軟件，分別使用ALE算法和PBM算法、SPH算法對試驗工況進行了數值模擬。

4.1 ALE有限元模型

由于試驗后臺架基座無明顯損壞且固定螺栓未脫落，臺架的約束邊界未遭到破壞，因此在仿真中簡化臺架模型，保留夾具與靶板，建立炸藥、空氣、靶板與夾具模型，如圖8所示。炸藥和空氣單元均為實體并采用Euler算法，空氣域為長方體，尺寸為1.6 m×1.6 m×1.1 m，為研究空氣網格尺寸對仿真計算效率及精度的影響，空氣單元尺寸分別采用10 mm、15 mm、20 mm。夾具、靶板均為實體并采用Lagrange算法，單元尺寸為10 mm。炸藥、臺架包含在空氣域內，模型中炸藥、空氣域臺架的相互作用通過關鍵字*CONSTRAINED_LAGRANGE_IN_SOLID定義，耦合方式采用罰函數法，夾具與鋼板之間的接觸采用*CONTACT_AUTOMATIC_SINGLE_SURFACE關鍵字定義，在空氣域四周設置無反射邊界條件(*BOUNDARY_NON_REFLECTION)模擬空氣的無限空間，防止爆炸沖擊波到達邊界發生反射影響仿真結果。

圖8 ALE有限元模型示意圖

靶板采用高強防彈鋼，其材料的本構關系采用*MAT_SIMPLIFIED_JOHNSON_COOK模型描述，如公式(3)所示，利用拉伸試驗機和霍普金森實驗裝置分析材料的靜、動態力學性能，獲得材料的主要參數見表1所示。

表1 高強防彈鋼主要材料參數

(3)

高能炸藥的本構關系采用*MAT_HIGH_EXPLOSIVE_BURN描述，狀態方程采用*EOS_JWL，方程如式(4)所示，使用初始體積法(*INITIAL_VOLUME_FRACTION_GEOMETRY)定義炸藥形狀與位置，其材料參數[17]見表2所示。

表2 炸藥主要參數

(4)

式中：P為爆轟壓力，A、B、R1、R2和ω為炸藥特性參數；V為初始炸藥相對體積；E0為炸藥單位體積初始內能。

空氣的本構關系采用*MAT_NULL描述，其狀態方程采用*EOS_LINEAR_POLYNOMIAL，狀態方程多項式如式(5)所示，空氣的主要參數[18]見表3所示。

表3 空氣參數及狀態方程系數

P=C0+C1μ+C2μ2+C3μ3+

(C4+C5μ+C6μ2)E0

(5)

式中：P為空氣壓力，C0、C1、C2、C3、C4、C5、C6為狀態方程多項式參數，E0為空氣單位體積初始內能。

4.2 SPH模型建立

圖9為SPH仿真模型，仿真模型中的鋼板與夾具的單元類型、尺寸、材料與ALE模型一致，區別在于SPH算法中忽略了空氣的影響，同時將TNT高能炸藥離散為光滑粒子，炸藥藥柱直徑119.2 mm，高109.8 mm，一共離散為10 000個SPH粒子，模型中炸藥與臺架之間的相互作用采用關鍵字*CONTACT_AUTOMATIC_NODES_TO_SURFACES定義，為研究SPH粒子數對仿真計算效率及精度的影響，SPH粒子數分別設置為5 000、10 000、50 000。在SPH算法中容易出現滲漏現象(大量的光滑粒子在沒有接觸力的情況下穿透了靶板)，因此采用以下兩點措施以改善這種狀況：① 細化炸藥模型，將高能炸藥離散為更多SPH粒子，并且進一步細化靶板有限元模型；② 通過接觸控制與時間步長控制，增加靶板中單元的接觸罰函數剛度系數和SPH單元屬性中的平滑步長系數以及減小質量縮放系數等方法是SPH粒子與臺架之間的接觸更加嚴格。

圖9 SPH仿真模型示意圖

4.3 PBM模型建立

圖10為PBM仿真模型，高能炸藥和空氣均采用粒子模擬，通過*DEFINE_PBLAST_GEOMETRY創建炸藥的幾何形狀，定義其特征尺寸與局部坐標參數，在此基礎上生成所需要的圓柱形炸藥。利用關鍵字*PARTICLE_BLAST建立TNT的材料模型，同時該關鍵字還定義了爆炸粒子與拉格朗日網格的相互作用方式，定義PBM模型中高能炸藥的爆速、初始內能和密度等參數見表4所示。此外，該關鍵字還定義了符合理想氣體定律的空氣粒子，夾具與鋼板的單元類型、尺寸與ALE模型中相同。

圖10 PBM仿真模型

表4 高能炸藥材料參數

4.4 仿真結果及分析

由于PBM不是基于連續介質力學，在數值模擬的過程中粒子是非可視的，只能看到定義的粒子點云的運動，如圖11所示。

所有仿真模型計算時間均設置為3 ms，使用同一臺計算服務器進行求解計算，ALE模型中研究空氣網格尺寸對計算精度的影響，將空氣單元尺寸分別設置為10 mm、15 mm、20 mm。PBM模型中研究炸藥粒子數量(Np)對計算精度的影響，在空氣粒子數為0的條件下分別設置了5 000、10 000、50 000、100 000個炸藥粒子。同時為了研究空氣粒子數(Na)對計算精度的影響，在炸藥粒子為50 000的條件下分別設置了0、20 000、50 000、100 000、200 000個空氣粒子。各仿真模型計算結果統計見表5所示。

表5 仿真結果與試驗結果

圖12表示了不同空氣網格尺寸的ALE模型仿真結果與PBM模型仿真結果。從圖11和表5中可以得出：ALE模型和PBM模型數值模擬所得到的靶板中心點撓度隨時間變化的趨勢一致，ALE模型的數值模擬所得到的靶板中心點撓度低于試驗測量值，PBM則略高于試驗測量值。值得注意的是，對于網格尺寸較大的ALE_2和ALE_3工況，數值模擬所得的結果明顯低于試驗測量值，誤差分別達到了15.3%、14.5%，而對于網格尺寸較小的ALE_1工況，數值模擬所得到的結果和試驗結果的誤差較小，僅有1.3%，同PBM模型數值模擬得到的結果有較好的一致性。此外，還可以發現2種算法模型中靶板中心點達到最大撓度的時間不同，ALE模型中靶板中心點達到最大撓度的時間在爆炸開始后1 ms左右，而PBM模型則明顯較晚，靶板中心點達到最大撓度的時間在爆炸開始后2 ms左右，2種算法中臺架的初始結構響應一致，相同的結構響應過程出現在了不同的時間，表明在結構響應方面PBM算法的仿真過程具有滯后性。

圖12 ALE模型與PBM模型仿真結果曲線

圖13為不同SPH粒子數的SPH模型與PBM模型的仿真結果，由于SPH算法中忽略了空氣的影響，仿真所得到的結果略高于ALE算法與PBM算法的仿真結果。由圖13可得：SPH算法與PBM算法數值模擬所得到的靶板中心點撓度隨時間的趨勢大致相同，SPH算法與PBM算法的數值模擬所得到的靶板中心點撓度都略高于試驗值，對于光滑粒子數較少的SPH_1和SPH_2工況，數值模擬所得到的結果誤差較大，分別達到了15.6%和6.7%，而對于粒子數較多的SPH_1工況，數值模擬的誤差較小，為4.0%，因此提高SPH粒子數有利于提高仿真結果的精確性。

圖13 SPH模型與PBM模型仿真結果曲線

為了研究PBM模型中炸藥粒子數量對仿真結果的影響，在空氣粒子為0的情況下，進行了PBM_1到PBM_4四種工況的數值模擬，4種工況靶板的中心點撓度隨時間的變化趨勢如圖14所示。由圖14可知：4種工況的數值模擬結果中靶板中心點撓度隨時間變化的趨勢大致相同，都是先增大至最大值，然后逐漸減小，但是爆炸粒子數量對靶板中心點撓度達到最大值的時間有一定影響。此外，從圖中還可以得出：炸藥粒子數對數值模擬結果的影響較大，靶板中心點的最大撓度隨著炸藥粒子數的增加而減小，并逐漸趨向于一個穩定值。

圖14 中心點撓度隨炸藥粒子的變化曲線

為了研究PBM模型中空氣粒子數量對數值模擬結果的影響，進行了PBM_3、PBM_5、PBM_6、PBM_7、PBM_8五種工況的數值模擬，五種工況靶板的中心點撓度隨時間的變化趨勢如圖15所示。由圖15可知：從整體上看，五種工況靶板的中心點撓度隨時間變化的趨勢大致相同，靶板的中心點撓度隨著時間先增大后減小。此外，在炸藥粒子數量一定的情況下，靶板中心點撓度的最大值隨著空氣粒子數的增加而增加，并且逐漸趨近于穩定值，在空氣粒子數大于50 000時，PBM模型中空氣粒子數對數值模擬的結果的影響很小。

圖15 中心點撓度隨空氣粒子數的變化曲線

對比PBM算法、ALE算法SPH算法的計算效率，分別統計了各工況計算時長，見表5所示。從整體來看，ALE模型和SPH模型的計算時長遠大于PBM模型的計算時長，網格尺寸對ALE數值模擬結果的影響較大，為提高計算精度減小網格尺寸會造成計算時長大幅增加,SPH粒子數對SPH模型的計算時長影響較大，且SPH算法的計算成本隨著SPH粒子數的增加而增加，PBM模型的計算時長隨粒子數的增加而增加，但遠小于ALE模型與SPH模型，并且粒子數對靶板中心點的最大撓度影響較小。

5 結論

1)ALE模型中網格尺寸對數值模擬結果影響較大，減小網格尺寸可以有效提高計算精度，當空氣網格尺寸與結構網格尺寸接近1:1時，仿真結果與試驗結果誤差較小，與PBM模型計算結果具有良好的一致性，與ALE算法相比，在結構響應方面PBM算法的仿真過程具有滯后性；SPH算法和ALE算法的計算精度相似,并且當SPH粒子數的較多時，仿真誤差較小，與PBM算法的一致性較好。

2)PBM模型中，在空氣粒子數量一定的情況下，炸藥粒子數對數值模擬結果的影響較大，靶板中心點的最大撓度隨著炸藥粒子數的增加而減小，并逐漸趨近于試驗值。

3)PBM模型中，在炸藥粒子數一定的情況下，靶板中心點的最大撓度隨空氣粒子數的增加而增加，并且當空氣粒子數大于50 000時，空氣粒子數量對數值模擬結果的影響較小，靶板中心點撓度的變形狀況趨于穩定。

4)在誤差較小的情況下，與ALE算法、SPH相比，PBM算法可以大大的節省計算成本，此外，PBM模型中，計算成本隨著粒子數的增加而增加并且增加炸藥粒子所增加的計算成本遠大于增加空氣粒子所增加的計算成本。