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

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

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

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

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

1 SPH方法

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

1.1 控制方程

高速?zèng)_擊時(shí)固體材料會(huì)表現(xiàn)出流體的性質(zhì),故基于核近似和粒子近似原理對(duì)描述沖擊問(wèn)題的控制方程進(jìn)行粒子近似,可以得到具有材料強(qiáng)度的流體動(dòng)力學(xué)控制方程的SPH形式[14],即:

(1)

1.2 人工粘性

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

(2)

(3)

其中:

(4)

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

1.3 邊界條件

為解決臨近邊界粒子點(diǎn)積分區(qū)域被截?cái)嗟膯?wèn)題,本文引入鏡像粒子法來(lái)處理邊界問(wèn)題,使其滿足物面邊界條件。

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

圖1 鏡像粒子法

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

(5)

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

(6)

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

2 模型建立

2.1 網(wǎng)格模型

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

圖2 網(wǎng)格模型

2.2 材料模型

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

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

(7)

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

(8)

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

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

2.3 邊界設(shè)置

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

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

2.4.1 主應(yīng)力矢量分布

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

圖3 不同時(shí)刻下鋼索的主應(yīng)力矢量分布

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

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

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

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

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

2.4.3 影響因素分析

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

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

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

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

3 SPH與FEM算法對(duì)比分析

3.1 系統(tǒng)總能量

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

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

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

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

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

3.2 應(yīng)力波分布

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

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

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

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

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

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

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

4.1 鋼索波動(dòng)理論

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

(9)

(10)

(11)

(12)

(13)

4.2 對(duì)比驗(yàn)證

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

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

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

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

5 結(jié)論

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

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

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

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