高速沖擊下鋼索波動行為SPH模擬研究

趙明基,王發(fā)展,黃克鵬

(西安建筑科技大學 機電工程學院,西安 710055)

鋼索的高速沖擊是一種在極短時間內(nèi)物體間發(fā)生的碰撞響應行為。沖擊速度、碰撞角度和材料性能的變化使其波動特性以及沖擊應力特征更為復雜,以現(xiàn)有實驗和有限元等手段分析鋼索高速碰撞響應機制尚存不足。光滑粒子流體動力學(Smoothed particle hydrodynamics, SPH)在沖擊動力學領(lǐng)域處理大變形特征的動態(tài)響應問題中,具有不可替代的優(yōu)勢[1-3]。因此,利用SPH算法對鋼索高速沖擊問題進行深入探究,研究結(jié)果具有潛在的軍用價值。

近年來,國內(nèi)外學者采用SPH算法對高速沖擊涉及的高應變率、大變形以及流動界面等問題進行了大量研究。Zhang和Liu等[4]將自適應光滑粒子流體動力學方法(ASPH)應用到具有材料強度的動力學問題中,并采用核函數(shù)修正法(KGC)對彈丸侵徹靶板問題進行模擬,但緩沖區(qū)的邊界粒子的核函數(shù)可能會被邊界截斷,產(chǎn)生截斷誤差仍具有不穩(wěn)定性。Sun等[5]將強可壓縮多相SPH體積自適應從笛卡爾坐標系擴展到援助極坐標系來解決軸對稱問題,但其無法計算拉伸力和剪切力。Wang等[6]采用FEM-SPH耦合方法對強夯引起的直流大變形問題進行數(shù)值分析,但該方法在三維計算中較為復雜,精度不高。Markauskas和Kruggel-emden[7]采用DEM-SPH耦合方法對濕法連續(xù)篩分進行了數(shù)值分析,提出了一種新的表征屏蔽表面的簡單模型,但其無法解決因為粒子缺失而帶來的能量損失問題。Riccardo等[8]采用SPH-FE耦合法對陶瓷基彈道盾牌高速沖擊問題進行模擬,但無法解決表面張力不穩(wěn)定問題。Young等[9]提出了一種新的SPH全拉格朗日公式與歐拉公式耦合方法,實現(xiàn)了全拉格朗日公式自適應轉(zhuǎn)換為歐拉公式,并采用此方法對斑片進行二維和三維高速沖擊模擬,但此方法未能解決拉伸不穩(wěn)定問題。Varas等[10-12]基于SPH-ALE耦合法建立有限元模型,研究了不同沖擊速度、不同充液比例對油箱毀傷效應的影響,但此方法不能較好的解決三維高速碰撞問題。以上作者的研究,均為SPH算法在高速沖擊問題中提供了有效的數(shù)值模擬手段,但仍存在計算效率低、拉伸不穩(wěn)定、邊界條件處理難和易產(chǎn)生零能模式等缺陷。但是,國內(nèi)外基于SPH法對鋼索在高速沖擊下的瞬態(tài)響應研究未見報道。

為了利用SPH算法對高速沖擊下的鋼索波動行為進行數(shù)值模擬,選取Monaghan型人工粘性方程有效解決了粒子飛濺問題,采用鏡像粒子法消除了粒子缺失的現(xiàn)象,建立了鋼索高速沖擊的動力學模型,揭示并掌握了其動態(tài)特性和瞬態(tài)響應規(guī)律,通過與FEM算法對比表明SPH算法所建立模型的高效性與優(yōu)越性。理論結(jié)果驗證了仿真所建模型的準確性,為鋼索高速沖擊中應用提供了理論參考。

1 SPH方法

光滑粒子流體動力學是一種流動質(zhì)點無網(wǎng)格方法,它將計算域離散成一系列承載著各種物理量的粒子,且粒子間通過核函數(shù)相互作用,使數(shù)值模擬時不會受到空間任意分布粒子的影響。相對于傳統(tǒng)的網(wǎng)格數(shù)值方法,在處理大變形、交界面位置捕捉、解決不連續(xù)性問題方面具有先天優(yōu)勢[13]。

1.1 控制方程

高速沖擊時固體材料會表現(xiàn)出流體的性質(zhì),故基于核近似和粒子近似原理對描述沖擊問題的控制方程進行粒子近似,可以得到具有材料強度的流體動力學控制方程的SPH形式[14],即:

(1)

1.2 人工粘性

為有效防止粒子間靠近時的非物理穿透現(xiàn)象,減少數(shù)值耗散,采用Monaghan型人工粘性方程[15],具體方程如下:

(2)

(3)

其中:

(4)

式中:A、B為標準常數(shù),且與A相關(guān)的項得到的是體積粘度,與B相關(guān)的項是用于防止高馬赫數(shù)粒子的相互穿透;c為聲速的速度矢量;v為粒子的速度矢量;φ用于防止粒子間相互靠近時產(chǎn)生的數(shù)值發(fā)散。

1.3 邊界條件

為解決臨近邊界粒子點積分區(qū)域被截斷的問題,本文引入鏡像粒子法來處理邊界問題,使其滿足物面邊界條件。

粒子鏡像處理過程如圖1所示。其中,i粒子為正在進行鏡像操作的粒子,j粒子為鏡像粒子。通常,i粒子與j粒子密度和壓力相同,使求解域內(nèi)實粒子在固壁邊界處法向不穿透和切向可滑移。

圖1 鏡像粒子法

二者的速度滿足以下關(guān)系:

(5)

在鏡像粒子法中關(guān)于應力張量的組成部分f的邊界修正式為

(6)

式中:fbc為邊界粒子在x=xi處預設的邊界張量;Wij為插值核函數(shù)。

2 模型建立

2.1 網(wǎng)格模型

本文采用LS-PrePost軟件對直徑10 mm,長度1 000 mm的鋼索和直徑30 mm,高50 mm的圓柱形剛體進行建模,利用SPH-Generation在鋼索寬、高方向各設10個粒子,長度方向設置10 000個粒子,粒子間距為1.2 mm,模型共80 000個粒子,使用SolidMesher對剛體進行網(wǎng)格自由劃分,如圖2所示。

圖2 網(wǎng)格模型

2.2 材料模型

由于鋼索實際上是由長細比很大的鋼絲繩和其他復合材料擰結(jié)而成,綜合性能等效于鋼材料,故采用Elastic-Plastic-Hydro材料模型結(jié)合Mie-Gruneisen狀態(tài)方程來表征其性能,設置密度為7.80 g/cm3,彈性模量E為210 GPa,剪切模量G為80 GPa。其次采用RIGID表征剛體材料性能,設置密度為7.83 g/cm3。

對于金屬大變形可引用Mie-Gruneisen狀態(tài)方程[16-17],定義壓力為

(7)

式中:pH為沖擊Hugoniot曲線上的點的壓力;Γ為Gruneisen參數(shù)。

(8)

對于不同金屬材料的狀態(tài)方程參數(shù)如表1所示。

表1 不同金屬材料Mie-Gruneisen狀態(tài)方程中的參數(shù)

2.3 邊界設置

將SPH粒子設為從節(jié)點,有限元與SPH粒子連接處的單元表面設為主表面。剛體速度設置為260 m/s、270 m/s 、280 m/s、290 m/s和300 m/s,沖擊角度設置為75°、80°、85°和90°,作用時間設置為200 μs,沖擊點與鋼索的摩擦因數(shù)設置為0.5、0.6、0.7和0.8,剛體與鋼索接觸類型為侵蝕接觸。

2.4 仿真結(jié)果與分析

2.4.1 主應力矢量分布

隨著剛體高速沖擊鋼索時,導致鋼索應力分布不均,通過分析主應力的變化趨勢,得到粒子受力規(guī)律。圖3為在沖擊速度為260 m/s沖擊下不同時刻鋼索的主應力矢量分布圖。

圖3 不同時刻下鋼索的主應力矢量分布

由圖3a)可知,當t=14 μs時,由于受到高速沖擊力作用,鋼索中心開始產(chǎn)生形變,鋼索內(nèi)部應力較大,鋼索中心處內(nèi)外兩側(cè)部分粒子產(chǎn)生與其他粒子相反的主應力;由圖3b)可知,當t=34 μs時,鋼索中心外側(cè)應力較大并開始產(chǎn)生斷裂,波峰處粒子產(chǎn)生的主應力比其他部位粒子大,且相切于應力傳遞方向;由圖3c)可知,當t=70 μs時,鋼索內(nèi)所有粒子均受到應力,鋼索中心處內(nèi)側(cè)部分粒子產(chǎn)生反應力,且絕對值大于其他部位粒子;由圖3d)可知,當t=112 μs時,鋼索完全斷裂,粒子產(chǎn)生的主應力減小,中心處粒子主應力方向指向沖擊方向。

2.4.2 結(jié)構(gòu)強度分布

在鋼索高速沖擊過程中,由于沖擊力的瞬間作用,導致鋼索受力不均,粒子分布離散,使得鋼索結(jié)構(gòu)強度改變。圖4為在沖擊速度為260 m/s沖擊下不同時刻鋼索的整體結(jié)構(gòu)強度分布。

圖4 不同時刻下鋼索的整體結(jié)構(gòu)強度分布

由圖4a)可知,當t=14 μs時,鋼索粒子開始受沖擊力作用,大量粒子由內(nèi)向外發(fā)散,導致粒子相互擠壓,部分粒子產(chǎn)生反向作用力,使中心結(jié)構(gòu)強度增加;由圖4b)可知,當t=34 μs時,鋼索外側(cè)開始斷裂,鋼索中心帶粒子呈圓弧狀分布,隨著應力波向兩端傳遞,粒子從中心向兩端分散,且分布對稱,使中心結(jié)構(gòu)強度銳減;由圖4c)可知,當t=70 μs時,鋼索內(nèi)側(cè)未斷裂,鋼索中心處粒子受力增大,且向內(nèi)外兩側(cè)發(fā)散,反作用力由外向內(nèi)逐漸增大,鋼索強度逐漸減小;由圖4d)可知,當t=112 μs時,鋼索內(nèi)外兩側(cè)完全斷裂,鋼索中心處粒子受力銳減,且沿中心向四周發(fā)散,鋼索強度降至最低。

2.4.3 影響因素分析

在摩擦因數(shù)μ=0.7不變條件下,分別選取不同沖擊角度θ,在沖擊角度θ=85°不變條件下,分別選取不同摩擦因數(shù)μ,分析不同摩擦因數(shù)與沖擊角度對鋼索沖擊應力的影響。圖5為不同沖擊角度下沖擊速度與鋼索沖擊應力大小的關(guān)系,圖6為不同摩擦因數(shù)下沖擊速度與鋼索沖擊應力大小的關(guān)系。

圖5 不同沖擊角度下沖擊速度與鋼索沖擊應力大小的關(guān)系

圖6 不同摩擦因數(shù)下沖擊速度與鋼索沖擊應力大小的關(guān)系

如圖5所示,滑動沖擊下鋼索第一次沖擊應力隨著沖擊角度的增大而減小,其差值隨著沖擊速度的增加而增大。如圖6所示,滑動沖擊下鋼索第一次沖擊應力大小隨著摩擦因數(shù)的增大而增大,其差值隨著沖擊速度的增加而增大,且增幅較為明顯。

3 SPH與FEM算法對比分析

3.1 系統(tǒng)總能量

本文進行的數(shù)值仿真中系統(tǒng)總能量是指內(nèi)能、動能、接觸面滑移能、沙漏能、系統(tǒng)阻尼能和剛性墻能的總和,圖7為在不同算法下不同沖擊速度的系統(tǒng)總能量-時間歷程曲線。

圖7 在不同算法下不同沖擊速度的系統(tǒng)總能量-時間歷程曲線

由圖7a)可知,當t=0～17 μs時,鋼索與圓柱體尚未接觸,故系統(tǒng)總能量不變;當t=17～60 μs時,系統(tǒng)總能量迅速減少;當t=60～180 μs時,系統(tǒng)總能量仍減少,但其速率明顯下降;當t=180～200 μs時,系統(tǒng)總能量趨于穩(wěn)定。

由圖7b)可知,當t=0～17 μs時,系統(tǒng)總能量不變;當t=17～20 μs時,系統(tǒng)總能量驟減;當t=20～80 μs時,系統(tǒng)總能量微小波動,但總體呈下降趨勢;當t=80～110 μs時,系統(tǒng)總能量趨于穩(wěn)定;當t=110～200 μs時,系統(tǒng)總能量緩慢增加。

由圖7a)和圖7b)仿真結(jié)果對比可知,系統(tǒng)在SPH算法與FEM算法下均有能量損失,且沖擊速度越大能量損失越大。由于粒子的相互作用以及粒子與單元間接觸算法,造成系統(tǒng)的能量損失,但SPH算法的能量損失要略少于FEM算法。一方面是由于本文剛體采用了網(wǎng)格劃分的方法,因此模型總粒子數(shù)要比純SPH算法的粒子數(shù)少;另一方面是由于粒子與單元接觸算法造成的能量損失要遠小于粒子間相互作用的能量損失。

3.2 應力波分布

在高速沖擊過程中,鋼索由直線預緊狀態(tài)轉(zhuǎn)為三角形布局的應力傳播狀態(tài)為彎折波。彎折波的傳播速度很大,但鋼索應力的傳播速度要遠大于彎折波,因此鋼索中任一點的運動都是從應力傳遞到該點的時刻開始。由于剛體與鋼索的相互作用等效于兩物體間的剛性碰撞,所以二者在嚙合處存在高頻率和大幅度的接觸力,使各處應力和應變不是均勻分布的,圖8和圖9分別為采用SPH算法和FEM算法在沖擊速度260 m/s沖擊下不同時刻下鋼索的范式應力云圖。

圖8 不同時刻下鋼索的Mises應力分布(SPH算法)

初期的攔阻動力學特性主要體現(xiàn)于鋼索本身。在剛體剛與鋼索接觸的時候,鋼索兩端無法在短時間內(nèi)進行動作,因此在碰撞初期鋼索的應力是由剛體與鋼索系統(tǒng)相互作用產(chǎn)生的,故鋼索邊界狀態(tài)對整個碰撞過程所產(chǎn)生的影響可忽略。

由圖8a)～圖8c)可知,在高速沖擊下,在t=14 μs時,鋼索中心處由內(nèi)至外產(chǎn)生形變,隨著時間推移,應力以波的形式由中心向兩端傳遞;在t=34 μs時,鋼索在外側(cè)開始產(chǎn)生斷裂,并向內(nèi)側(cè)延展,但應力波依舊向兩端傳遞;在t=94 μs時,鋼索內(nèi)外兩側(cè)完全斷裂,但應力波依舊向兩端傳遞,并不斷進行著反射疊加。

由圖9a)～圖9c)可知,在高速沖擊下,在t=14 μs時,鋼索中心處由內(nèi)至外產(chǎn)生形變,隨著時間推移,應力以波的形式由中心向兩端傳遞;在t=34 μs時,鋼索外側(cè)達到斷裂極限,但由于網(wǎng)格的限制,鋼索外次沒有產(chǎn)生斷裂,應力波一直在鋼索內(nèi)傳遞并不斷進行著反射疊加。

由圖8與圖9對比可知,采用FEM算法進行的仿真形變沒有SPH算法明顯,無法完全仿真出彈性性能。若網(wǎng)格劃分過大,則會因為網(wǎng)格畸變導致計算終止;若網(wǎng)格劃分過小,則會僅為計算單元過多導致計算時間過長,效率比SPH算法低。故采用SPH算法對鋼索高速沖擊動力學進行數(shù)值模擬比傳統(tǒng)網(wǎng)格分析計算精度高,效率高,表明了SPH算法在處理大變形、解決不連續(xù)性問題方面的優(yōu)越性。

4 理論及結(jié)果對比

4.1 鋼索波動理論

由于設置的沖擊角度為90°,且初始應力σ0遠小于由于沖擊而產(chǎn)生的新應力σ,可以將波動模型的理論公式簡化,即:

(9)

(10)

(11)

(12)

(13)

4.2 對比驗證

為了驗證SPH算法在模擬鋼索的波動特性上的準確性,本文通過對鋼索波動參數(shù)的理論值和仿真值進行對比。表2和表3分別為采用SPH算法和FEM算法對鋼索在不同沖擊速度下的波動參數(shù)的仿真值與理論值的對比結(jié)果。

表2 不同沖擊速度下的鋼索波動參數(shù)(SPH算法)

表3 不同沖擊速度下的鋼索波動參數(shù)(FEM算法)

由表2可知,在LS-DYNA中用SPH算法得到不同沖擊速度下的縱波、橫波和彎折波的波速均比理論值要小,這是由于鋼索波動性理論分析中沒有考慮鋼索內(nèi)部阻尼。對比可知,SPH算法下的理論值與仿真值高度一致,誤差均在6%以內(nèi),而FEM算法下的誤差在10%左右。從而驗證了仿真所建模型的準確性與優(yōu)越性。

5 結(jié)論

為研究高速沖擊下鋼索應力波的傳播規(guī)律,建立了特殊連續(xù)彈性體模型,通過人工黏性與鏡像粒子法相結(jié)合,解決了鋼索數(shù)值模擬中粒子飛濺和粒子缺失的問題,分別利用SPH算法和FEM算法對鋼索高速沖擊過程進行數(shù)值模擬,得到如下結(jié)論:

1) 在高速沖擊下,鋼索中心處由內(nèi)至外產(chǎn)生形變,隨時間推移,應力波由中心向兩端傳遞,并產(chǎn)生反射疊加現(xiàn)象。通過SPH算法與FEM算法對比可知,二者系統(tǒng)總能量損耗相似,同時采用SPH算法對高速沖擊下鋼索波動行為進行數(shù)值模擬計算精度更高,效率更高。

2) 在高速沖擊作用下,鋼索內(nèi)外兩側(cè)部分粒子產(chǎn)生反應力,中心處粒子由內(nèi)向外發(fā)散,結(jié)構(gòu)強度增加;鋼索中心處外側(cè)產(chǎn)生斷裂時,波峰處粒子產(chǎn)生較大主應力,且相切于應力傳遞方向,中心帶粒子呈圓弧狀分布,粒子向兩端分散,結(jié)構(gòu)強度銳減;鋼索中心內(nèi)側(cè)未斷裂時,所有粒子均受到應力作用,且內(nèi)側(cè)部分粒子產(chǎn)生絕對值較大的反應力,粒子整體向內(nèi)外兩側(cè)發(fā)散,結(jié)構(gòu)強度減小;鋼索完全斷裂后,中心處粒子主應力方向指向沖擊方向,粒子沿中心向四周發(fā)散,結(jié)構(gòu)強度降至最低。

3) 滑動沖擊下鋼索沖擊應力隨沖擊角度的增大而減小,隨摩擦因數(shù)的增大而增大,其差值隨著沖擊速度的增加而增大。通過理論驗證了模型的正確性,為鋼索的設計與優(yōu)化提供理論參考。