吳夢陽,周雲波,趙 豐,王利輝,李進軍,張 明,周 迪
(1.南京理工大學機械工程學院, 南京 210094;2.陸軍某項目管理中心, 北京 100073;3.中國人民解放軍32379部隊, 北京 100073;4.中國人民解放軍32381部隊, 北京 100073)
摘要:以某淺埋衝擊試驗臺架為研究對象,引入新型爆炸粒子算法PBM模擬衝擊波,用離散單元法DEM模擬土壤,構建近距淺埋衝擊試驗臺架仿真模型,分析不同DEM粒子半徑對仿真結果的影響;建立了任意拉格朗日-歐拉算法的完整模型,對比兩種模擬方法的仿真結果;通過開展實爆試驗表明PBM-DEM粒子算法計算結果與試驗結果更接近,計算準確度更高。
關鍵詞:淺埋爆炸;臺架試驗;任意拉格朗日-歐拉算法;PBM-DEM粒子算法;數值計算
現代戰爭中防護型車輛時常面臨著簡易爆炸裝置IED(Improvise Explosive Device)的威脅,其巨大的爆炸衝擊有可能擊穿車輛底部,給車輛及車內乘員帶來毀滅性的損傷。因此,近年來針對軍用防護型車輛底部抗衝擊技術的研究逐步發展,通過分析爆炸衝擊應力波在車輛結構及乘員系統中的傳播路徑,逐層降低應力波強度,以降低乘員損傷指標作為最終目的,系統化提升整車抗爆炸衝擊能力,致力於降低車輛及乘員的損傷[1]。
自北約STANAG 4569標準發布後,各國逐步開展車輛抗爆炸試驗研究,試驗技術也逐漸走向成熟。整車爆炸試驗是評估車輛抗爆炸衝擊能力最直接的手段,為車輛抗爆炸防護設計提供最有效的素材。但是,由於整車爆炸試驗周期長、可重複性不高等原因,利用CAE仿真技術,結合試驗結果,可以有效預測爆炸環境下車輛及乘員的損傷效應,為後期車輛防護能力優化設計提供方向[2]。
近年來,用任意拉格朗日-歐拉算法(ALE)和流固耦合算法模擬爆炸源、衝擊波傳遞以及結構響應是公認的爆炸數值模擬算法之一。但是由於ALE算法需要建立大量的空氣和土壤單元,造成計算效率不高,並且歐拉算法與結構耦合的魯棒性較低,可能會造成計算結果不穩定。
在此基礎上,基於粒子的無網格法逐漸被開發,如SPH、PBM(Particle-Blast Method)、DEM(Discrete-Element Method)等,用來協調計算的魯棒性。李利莎[3]使用Lagrangian算法、ALE算法和SPH算法分別模擬了炸藥衝擊混凝土過程,對比了三種不同算法的優劣性。Trajkovski[4]使用了ALE算法和SPH算法模擬爆炸衝擊靶板,比較兩種算法的差異,發現SPH算法可有效避免衝擊波滲漏問題。Olovsson等[5]首次使用PBM算法進行爆炸模擬,之後,滕海龍[6]對算法進行了改進,使用PBM算法進行了一系列的近場爆炸模擬計算,計算結果與試驗結果具有較好的一致性,同時發現PBM算法在計算經濟性上優於ALE方法。Børvik等[7]用離散單元法DEM模擬了土壤介質,用來研究球形炸藥在土壤中的爆炸特性。Jensen等[8]使用PBM-DEM粒子算法對地雷爆炸進行了模擬研究,旨在了解不同土壤密度、爆炸物形狀等變量下車體的結構響應。
本文主要研究內容是根據某淺埋爆炸試驗臺架試驗要求,分別使用PBM-DEM爆炸算法以及ALE算法進行仿真分析;比較了兩種算法結果中試件背板最大殘餘變形量、臺架的運動量、試件背板的總能量、泡沫鋁內能以及土壤動能;最終進行實爆試驗,驗證兩種算法的精度,對比分析得出PBM-DEM粒子算法計算結果與試驗結果更接近,其計算時間也比ALE算法計算時間短。
1 PBM-DEM數值模擬算法理論1.1 PBM粒子爆炸算法PBM粒子算法是針對基於分子動力學理論KMT(Kinetic Molecular Theory)的微粒法CPM(Corpuscular Method)進行了改進。CPM主要是用於模擬氣囊的爆炸驅動,通過將大量的分子整合為一個粒子,整合比例接近1018∶1。CPM避免了時間尺度與長度尺度的矛盾性,但CPM不適用於熱平衡假設失效的爆炸模擬。
CPM是建立在KMT基礎上的一種粒子法,允許對氣體體積進行數值處理時存在一定偏差。在微觀層面上動力學分子理論表現為研究氣體分子及其相互作用,宏觀層面上表現為理想氣體定律。該理論基於以下假設[9]:
1) 分子之間的平均距離比它們自身尺寸要大。
2) 存在熱力學平衡,即分子處於隨機運動狀態。
3) 分子遵循牛頓運動定律。
4) 分子-分子和分子-結構的相互作用是完全彈性的。
動力學分子理論最早可追溯於1738年,當時Daniel Bernoulli提出了一個理論,即空氣對活塞的壓力是通過離散的分子碰撞建立起來的。以動力學理論為出發點,James Clerk Maxwell導出了熱平衡時分子速度分布的一個非常精確的表達式[10]:
式中: f(v)是分子速度的概率密度函數; M是摩爾質量; R是通用氣體常數; T是溫度。
使用CPM做爆炸模擬時,每克爆轟產物的分子量大約在1023 ~1024之間。所以,在整個爆炸過程中,不可能做到對每個分子進行模擬。在權衡了數值計算精度要求和計算效率之後,針對分子動力學基礎的微粒法進行了以下假設[11]:
1) 粒子均被定義為剛性球體;
2) 每個粒子包含了1015~1020個分子,具體數值取決於實際情況;
3) 對於每個單獨的粒子,平動動能Wt和轉動動能Ws之間是平衡的。
但在實際計算過程中,需要注意真實分子雲的平動能和轉動能的比值是一個統計值,不是對每個分子都有效。在熱平衡假設的基礎上,確定了Wt和Ws比例。對於大多數安全氣囊應用來說,這是一個有效的假設,但在氣流速度極高的爆炸模擬中則不是這樣。故PBM算法針對CPM假設進行如下補充[12]:
1) Ws是每個粒子的集中標量變量。因此,粒子沒有被賦予任何旋轉/振動的自由度。
2) 當粒子撞擊結構時,Wt值將發生變化,而Ws假定不受影響。因此,Ws和Wt的比值發生變化。
3) 當兩個粒子碰撞時,它們的能量被重新分配,以恢復或保持Ws和Wt之間的正確比例。這樣做的方式確保了總能量平衡和動量平衡。
1.2 DEM離散單元法在模型前處理過程中,DEM方法是將模型離散為剛性球體,並通過罰函數定義接觸。DEM是根據牛頓運動定律分別計算每個粒子的運動,採用六自由度系統來描述運動。圖1描繪了DEM粒子間相互作用,其中包括法向和切向的剛度和阻尼力、靜態摩擦力、滾動摩擦力和表面張力。由dint定義的兩粒子相互作用距離。
dint=|x1-x2|-(r1+r2)
其中: xn是球體中心的坐標矢量; rn是球形粒子的半徑。
圖1 粒子相互作用示意圖
粒子碰撞狀態分為3種,如圖2所示,當dint小於0時,表面張力和機械接觸都存在,這可能會導致粒子之間相互滲透;當兩個粒子不直接接觸時,只存在表面張力;當兩個dint足夠大時,粒子間不再存在相互作用力。
圖2 粒子不同接觸狀態示意圖
法向彈簧力Fkn表達式為:
Fkn=Kndint
式中: Kn為法向接觸剛度; KNorm為法向彈簧常數的比例因子;kn為根據彈性模量En和泊松比vn計算得出的壓縮模量。
切向剛度表達式為:
Kt=KnKShear
法向阻尼力Fdn決定了兩個粒子之間的阻尼特性,表達式為:
Fdn=Dndint
式中: mn為質點質量;Dn定義法向臨界阻尼比;dnorm為法向阻尼係數。正像法向阻尼力一樣,切向阻尼力由切向阻尼係數dtan定義。
為了驗證某鋼板複合泡沫鋁夾心結構的抗爆炸性能,設計了一款試驗臺架,三維模型如圖3所示。臺架由外廓支架、試驗臺架以及配重塊組成,均由Q235製造。試驗臺架部分與外廓支架間設立四根滑動導杆,以確保其可以垂向運動。試件安裝在試驗臺架底部,兩端由M24螺栓固定。
試件為鋼板複合泡沫鋁三明治夾芯結構,總體尺寸為600 mm×600 mm×40 mm。面板及背板為6 mm防彈鋼板,芯層結構為28 mm泡沫鋁,鋼板與泡沫鋁使用膠結。
圖3 試驗臺架三維模型示意圖
LS-DYNA是由Livermore Software Technology Corporation(LSTC)開發的高級通用多物理場仿真軟體包。該軟體包被用來計算許多複雜的現實世界問題,其起源和核心競爭力在於使用顯式時間積分的高度非線性瞬態動態有限元分析(FEA)。LS-DYNA多用於汽車,航空航天,建築和土木工程,軍事,製造和生物工程行業。
2.1 PBM-DEM仿真計算根據試驗臺架模型建立1∶1爆炸仿真模型,如圖4所示。臺架模型劃分拉格朗日六面體實體單元網格,網格單元尺寸為10 mm,臺架材料使用選擇關鍵字Mat_Plastic_Motional來定義;試件建立的是5 mm尺寸的實體單元,由於爆炸衝擊的響應時間短,需要考慮材料的應變率效應,因此,防彈鋼材料使用Mat_Johnson_Cook材料模型來描述,防彈鋼材料參數進行霍普金森杆進行動態拉伸試驗獲得,拉伸試樣如圖5所示;泡沫鋁材料模型使用Mat_Crushable_Foam來表示,壓縮曲線由準靜態壓縮試驗獲得,準靜態壓縮試驗如圖6所示。鋼材料模型的主要參數如表1、表2所示(由實驗數據擬合計算得出)。防彈鋼應力應變曲線如圖7,泡沫鋁應力應變曲線如圖8。
根據試驗內容要求,試件迎爆面離地高度為300 mm,炸藥埋深100 mm。根據模型坐標,計算得出炸藥位置,通過Define_Pblast_Geometry關鍵字定義炸藥形狀特徵,並使用Particel_Blast關鍵字定義TNT材料模型,其定義的爆轟波參數如爆速、初始內能以及密度均來源於JWL狀態方程,同時該關鍵字還定義了爆炸粒子與DEM粒子和結構的耦合。
圖4 PBM-DEM模型
圖5 霍普金森拉伸試樣
圖6 準靜態壓縮試驗
圖7 防彈鋼應力應變曲線
圖8 泡沫鋁應力應變曲線
表1 Q235材料參數
表2 防彈鋼材料參數
根據前期大量仿真經驗發現,DEM粒子半徑需要保持在合理的範圍區間內。半徑過大會導致粒子攜帶動能過大,造成擊破模型現象;半徑過小會導致粒子容易穿透拉格朗日網格從而造成滲漏,同時粒子數增加,影響計算效率。一般設置粒子半徑在耦合的拉格朗日網格尺寸的1~1.5倍,由於耦合的試件網格尺寸為5 mm,故將土壤粒子半徑設置為5/6/7/8 mm分別進行計算。通過Disc_Sphere_Generation關鍵字構建DEM粒子。使用Mat_Rigid材料模型定義DEM粒子,粒子間相互作用關係由關鍵字Define_discrete定義,粒子與結構間的相互作用關係由Define_De_To_Surface_Coupling定義,部分參數見表3[8]。
表3 Define_De_To_Surface_Coupling參數
2.2 ALE仿真計算為了與PBM-DEM模型比較,根據試驗內容要求同時建立了基於ALE算法的臺架有限元模型,模型如圖9所示。數值流場長為2 000 mm,空氣場高1 800 mm,土壤場深400 mm,流場的網格尺寸為15 mm。空氣場材料採用Mat_Null模型與線性多項式狀態方程定義[13]:
P=C1+C2μ+C3μ2+(C4+C5μ+C6μ2)E
P表示氣體壓力,Ci(i=0、1、2、3、4、5、6)表示多項式係數,E表示單位體積空氣的初始內能。C0=-0.1,C1=C2=C3=C6=0,C4=C5=0.4;E0=0.253 mJ/mm3;V0=1[14]。土壤的超塑性材料模型由Mat_Soil_And_Foamble_Fail定義。
採用關鍵詞Mat_High_Explosive_Burn來描述TNT的材料模型。爆炸產物的壓力、體積和內能之間的關係由EOS_JWL定義。
式中ω、AJWL、BJWL、R1和R2是與爆炸材料有關的參數; V是相對體積;E0初始內能。TNT的材料性能和JWL參數見表4[14]。
圖9 ALE模型
表4 JWL參數
3 計算結果通過分析粒子半徑在5/6/7/8 mm工況下的計算結果發現,7/8 mm粒子半徑下依然出現擊穿迎爆面的現象;6 mm粒子半徑工況下未出現擊穿情況,但是出現了網格不平滑現象;5 mm粒子半徑工況下,迎爆面呈現了較平滑的耦合狀態,同時未出現粒子穿透現象。具體如圖10~圖13所示。
圖10 半徑8 mm計算結果 圖11 半徑7 mm計算結果
圖12 半徑6 mm計算結果 圖13 半徑5 mm計算結果
當DEM土壤粒子半徑為5 mm時,試件面板未出現被DEM粒子擊穿現象,且粒子與試件有了較好的接觸耦合,所以選用該工況仿真計算結果作為最終PBM-DEM粒子算法計算結果。
3.2 兩種數值模擬方法比較分析對比兩種算法背板殘餘變形量X1、背板的總能量E1、泡沫鋁內能E2、土壤動能Esoil,臺架運動量X2及計算時長t。具體仿真工況結果見表5,其背板變形雲圖和有關曲線見圖14~圖19。
表5 仿真工況結果
圖14 PBM-DEM計算結果
1) 背板最大殘餘變形量曲線對比可以看出,通過PBM-DEM粒子算法計算得出的背板最大殘餘變形量為28.5 mm,ALE算法計算得出的背板最大殘餘變形量為31.3 mm。
2) 背板總能量包括動能以及變形產生的內能,峰值上兩種算法的差異在5%左右,但是最終能量保持階段中,ALE算法中由於使用殼單元模擬背板,模型計算仍在震蕩未平穩,所以殘餘能量較高,但是從曲線中可以看出殘餘能量正在緩慢下降,相比下PBM-DEM粒子算法中背板最終殘餘能量曲線已經基本水平。
3) 泡沫鋁內能曲線對比中可以看出,當爆炸衝擊來臨時,泡沫鋁急劇壓縮吸收能量,內能曲線迅速升高,曲線拐點處能量差別在6%左右。由於ALE算法中忽略了土壤的動能效應,所以當爆炸衝擊完畢後,ALE算法泡沫鋁內能曲線基本水平;PBM-DEM粒子算法中考慮到土壤與結構的耦合效應,土壤衝擊時域上晚於衝擊波,所以在試件經歷了衝擊波之後,土壤繼續對試件產生能量傳遞,造成泡沫鋁吸收能量繼續增長。
圖15 ALE計算結果
圖16 土壤能量曲線
圖17 背板能量曲線
圖18 泡沫鋁內能曲線
圖19 臺架運動量曲線
4) 臺架運動量上對比中可以看出,PBM-DEM粒子算法臺架最大跳動量為68.2 mm,ALE算法臺架最大跳動量為64.1 mm,時域上峰值出現時間PBM-DEM粒子算法較ALE算法滯後。
5) 土壤能量曲線對比中可以看出,峰值上差異性在4%,總體趨勢類似,但是由於ALE算法在土壤飛濺超出計算域時會被刪除,而PBM-DEM粒子算法中土壤粒子不會被刪除,只是由於粒子與結構碰撞引起土壤能量降低,所以其土壤剩餘能量大於ALE算法。
6) 兩種算法計算相同的時長上,PBM-DEM粒子算法需要的時間更短,計算效率更高。
4 臺架試驗4.1 試驗內容本次臺架試驗目的是測試某型鋼板-泡沫鋁複合材料抗爆炸衝擊性能,同時驗證前期爆炸仿真精度。
根據試驗要求製作臺架如圖20所示。外廓支架均採用10 mm厚的Q235鋼材製造;試驗臺架由四根圓柱導杆及Q235實心扁鋼製作;配重塊按照要求配至1 t;試件兩端採用M24螺栓連接在活動支架上。
根據前期仿真邊界,將試件迎爆面調節至離地高300 mm處,地雷埋深100 mm,實現中心引爆。
圖20 臺架試驗裝置
圖21-圖22中顯示了爆炸後臺架及試件的整體狀況。根據試驗前對背板的網格處理,針對各點變形進行測量,最終擬合背板殘餘變形曲面,如圖23所示。試驗前在導杆上塗抹黃油,根據測量試驗後黃油在導杆上的痕跡長度,來描述試驗臺架的垂向運動量,如圖24所示。
圖21 臺架試驗後狀態
圖22 試件試驗後狀態
圖23 變形曲面擬合
圖24 臺架運動量測量
試驗結果表明,背板最大殘餘變形量為29.3 mm,坐標為(200,250),臺架活動部件最大垂直位移為66.7 mm。
4.2 仿真試驗對比分析根據試驗結果測出的各點變形量,擬合變形曲面,發現最大變形區域偏向於試件無約束段,有螺栓約束段變形較小。對比兩種算法計算結果與試驗實際測量結果如表6所示。
表6 仿真結果與試驗結果
PBM-DEM粒子算法計算結果最大變形點與試驗值差別在30 mm左右,最大變形量誤差為3%,臺架運動量差別為2.5%;ALE算法計算結果最大變形點與試驗值差別在133 mm左右,最大變形量誤差為3%,臺架運動量差別在4%;對比分析得出背板變形情況與PBM-DEM粒子算法計算結果更接近。
1) 採用了PBM-DEM粒子算法模擬了淺埋地雷爆炸衝擊試驗臺架試驗,分析了DEM粒子半徑對仿真計算的影響,發現粒子半徑和耦合結構網格尺寸一致時,DEM粒子與結構耦合效果較好,且未產生粒子滲透的現象。
2) 採用ALE算法模擬了淺埋地雷爆炸衝擊試驗臺架試驗,分析了兩種不同算法計算結果中的背板殘餘變形量、背板的總能量、泡沫鋁內能、土壤動能,臺架運動量及計算時長。發現兩種算法在前15 ms內的計算結果吻合較好,15 ms後由於PBM-DEM粒子算法中土壤粒子與結構耦合效應,導致結果差生差異。
3) 和ALE算法相比,PBM-DEM算法仿真計算結果與實爆試驗結果更接近,擬合度更高,且計算效率更高。