趙光明,宋順成

(1.安徽理工大學(xué)煤礦安全高效開(kāi)采省部共建教育部重點(diǎn)實(shí)驗(yàn)室,安徽 淮南232001;2.西南交通大學(xué)力學(xué)與工程學(xué)院,四川 成都610031)

1 引 言

動(dòng)力侵徹過(guò)程是力學(xué)研究中的一個(gè)重要課題,它與軍事、科技、國(guó)民經(jīng)濟(jì)的發(fā)展密切相關(guān)。因?yàn)榍謴貙?shí)驗(yàn)速度有限、費(fèi)用比較高,數(shù)值模擬成為侵徹過(guò)程研究的有效方法。由于侵徹過(guò)程涉及高溫、高壓和高速度,壓力梯度、溫度梯度、位移速度梯度很大,侵徹過(guò)程中材料破壞涉及沖塞、侵徹、層裂、崩落、濺飛等各種不同復(fù)雜形式,因此侵徹過(guò)程研究及其數(shù)值模擬是動(dòng)力學(xué)中的難點(diǎn)[1-2]。目前常用的方法有拉格朗日(Lagrange)、歐拉(Euler)、拉格朗日-歐拉耦合(A LE)以及光滑質(zhì)點(diǎn)流體動(dòng)力學(xué)(SPH)法等。歐拉法難以跟蹤界面位置,計(jì)算時(shí)間長(zhǎng),強(qiáng)度(失效)狀態(tài)和位移歷程關(guān)系計(jì)算精度差等;拉氏法當(dāng)網(wǎng)格變形嚴(yán)重時(shí)需要網(wǎng)格重分,導(dǎo)致計(jì)算精度下降,甚至計(jì)算失敗。而SPH 法更適合超高速撞擊等大變形、高應(yīng)變率現(xiàn)象的數(shù)值模擬,是一種無(wú)網(wǎng)格的拉格朗日法,比其他方法具有更大的實(shí)用性[3-5]。W.K.Liu等[6-7]的研究表明,SPH 法不能完全滿(mǎn)足一致性條件,而且采用此法偏微分方程時(shí),邊界點(diǎn)的解不穩(wěn)定,他們通入修正函數(shù),對(duì)SPH 方法進(jìn)行改進(jìn),提出了新型的再生核質(zhì)點(diǎn)法(reproducing kernel particle method,RKPM)。

侵徹過(guò)程的數(shù)值模擬,在很大程度上取決于材料的本構(gòu)模型。本文中,通過(guò)引入Johnson-Cook 本構(gòu)模型及損傷模型,采用RKPM 無(wú)網(wǎng)格法實(shí)現(xiàn)高速侵徹過(guò)程數(shù)值模擬,同理論解的結(jié)果進(jìn)行比較,驗(yàn)證數(shù)值模擬方法的有效性。

2 RKPM 形函數(shù)

由于SPH 對(duì)于邊界點(diǎn)以及不規(guī)則離散點(diǎn)不能精確求解,W.K.Liu 等[6-7]對(duì)SPH 方法作了改進(jìn),提出了再生核質(zhì)點(diǎn)法。在SPH 方法[3]中,任意函數(shù)可以近似為

式中:x 是空間點(diǎn),s 是領(lǐng)域內(nèi)的點(diǎn)。在Dirac Delta 條件下,很難解決離散化問(wèn)題,為此,進(jìn)行修改[6]

式中:uα(x)稱(chēng)為u(x)的再生函數(shù),α是膨脹系數(shù),修正核函數(shù)ˉφα可以表示為

式中:C(x-s)是修正函數(shù)。通過(guò)修正,RKPM 方法可以避免SPH 方法在區(qū)域的適用性

式中

根據(jù)Taylor 級(jí)數(shù)展開(kāi)式(2),可以得到

所以

式中

將再生方程離散為

式中:xi 表示質(zhì)點(diǎn),n 是計(jì)算域內(nèi)質(zhì)點(diǎn)數(shù)量,ΔV i 表示質(zhì)點(diǎn)x i 的體積,Ψi 稱(chēng)為再生核質(zhì)點(diǎn)法的形函數(shù)。

3 侵徹材料的本構(gòu)關(guān)系

根據(jù)Johnson-Cook 本構(gòu)方程[8],材料的Mises 等效應(yīng)力

式中:ε為等效塑性應(yīng)變,﹒ε*=﹒ε/﹒ε0,為量綱一等效塑性應(yīng)變,通常取﹒ε0=1 s-1;T*=T/T m,為量綱一溫度,T m 為材料熔化溫度;A、B、C、n 和m 為材料常數(shù)。

彈丸材料的狀態(tài)方程一般可以通過(guò)Grüneisen 狀態(tài)方程表示

式中:K 1、K 2 和K 3 是材料常數(shù),μ=ρ/ρ0-1。

根據(jù)Johnson-Cook 模型,單元內(nèi)材料的損傷演化定義為

式中:D ≤1,當(dāng)D=1 時(shí)材料出現(xiàn)斷裂;Δε為每個(gè)積分步長(zhǎng)中等效塑性應(yīng)變的增量,εf為等效斷裂應(yīng)變,該應(yīng)變與應(yīng)變率、溫度及壓力有關(guān)

4 滑移面算法

通?；泼嫣幚頃r(shí),滑移面條件(不可侵入條件、法向接觸力為壓力的條件及切向摩擦力的條件)都是不等式約束,需要通過(guò)拉格朗日乘子法或罰參數(shù)法將約束條件引入變分原理。具體處理時(shí),事先通過(guò)試探的方法將不等式約束改為等式約束引入方程并進(jìn)行方程求解,計(jì)算后利用不等式條件對(duì)解的結(jié)果進(jìn)行檢查,這樣反復(fù)迭代求解,直至滿(mǎn)足條件。傳統(tǒng)的滑移面處理過(guò)程異常復(fù)雜,計(jì)算量非常大,精度低,對(duì)于RKPM 不適合,也不易實(shí)現(xiàn),為此引入了新的滑移面處理法[9]。

圖1 給出了二維問(wèn)題滑移面處理方法。計(jì)算過(guò)程中,首先檢查從滑移面質(zhì)點(diǎn)A 在下一時(shí)間步時(shí)處于主滑移面的哪兩個(gè)質(zhì)點(diǎn)之間,并檢查從滑移面質(zhì)點(diǎn)與主滑移面的關(guān)系。在圖2 中,ED和F N 是兩條鉛垂線(xiàn),E F 與DF、與F H 相互垂直。

從滑移面質(zhì)點(diǎn)A 移動(dòng)后可能位于主滑移面的3 種位置,即A 可能在Ⅰ、Ⅱ或Ⅲ區(qū)域內(nèi)點(diǎn)′ 或。如果在Ⅰ內(nèi)位置,說(shuō)明從質(zhì)點(diǎn)移動(dòng)后在滑移面之外,不發(fā)生接觸,此時(shí)質(zhì)點(diǎn)間的各種物理量不需要發(fā)生交換;如果在Ⅱ、Ⅲ內(nèi),說(shuō)明從質(zhì)點(diǎn)已侵入主滑移面內(nèi),按不侵入條件,則立即將質(zhì)點(diǎn)A垂直拉到滑移面上。對(duì)于區(qū)域Ⅱ而言,)應(yīng)垂直拉到EF 面的處,并且重新確定從質(zhì)點(diǎn)A(A′)和主質(zhì)點(diǎn)E、F 的法向速度分量,并保持切向速度分量不變

圖1 滑移面Fig.1 Sliding surface

圖2 二維滑移面分析Fig.2 Analysis on 2D sliding surface

式中:m、l 分別表示質(zhì)點(diǎn)質(zhì)量和線(xiàn)段長(zhǎng)度,v+、v-分別是拉動(dòng)后、前的質(zhì)點(diǎn)的法向速度?？梢宰C明,式(18)～(20)既保證了滑移面質(zhì)點(diǎn)的動(dòng)量守恒,也保證了動(dòng)量矩守恒。如果質(zhì)點(diǎn)A 位移后在區(qū)域Ⅲ內(nèi)A′3位置,此時(shí)還須檢查相鄰主滑移面質(zhì)點(diǎn)線(xiàn)段F G 的位置情況。如果FG 在F H 以下,如FG 在FG2 位置,則將垂直拉到FG2上的,速度可按上述原理發(fā)生交換;如果FG 在F H 以上(如F G在FG1位置)或與F H 在同一位置,質(zhì)點(diǎn)A 只與質(zhì)點(diǎn)F 在垂向發(fā)生速度調(diào)整,即A 和F 調(diào)整后的垂向速度

切向方向速度不發(fā)生改變。

對(duì)于三維問(wèn)題計(jì)算,首先以質(zhì)點(diǎn)為頂點(diǎn),將主滑移面劃分成相連的三角形區(qū)域。當(dāng)主、從滑移面發(fā)生接觸時(shí),判斷從滑移面上的某一質(zhì)點(diǎn)處于主滑移面哪一個(gè)三角形區(qū)域,并按照二維問(wèn)題類(lèi)似的標(biāo)準(zhǔn)判斷,將該從接觸質(zhì)點(diǎn)與三角形頂點(diǎn)的3 個(gè)質(zhì)點(diǎn)發(fā)生速度調(diào)整。

5 速度配點(diǎn)法

對(duì)于給定的速度邊界條件及滑移面質(zhì)點(diǎn)速度的調(diào)整,本文中提出了速度配點(diǎn)法。在計(jì)算域Ω內(nèi)的N P個(gè)質(zhì)點(diǎn)可以分為NB個(gè)給定速度質(zhì)點(diǎn)和NN個(gè)非給定速度質(zhì)點(diǎn),。RKPM 法中,質(zhì)點(diǎn)速度可類(lèi)似給出

給定質(zhì)點(diǎn)速度

將式(23)代入式(22),可以得到

簡(jiǎn)記為

則實(shí)現(xiàn)方法為

本質(zhì)邊界條件及加速度邊界條件也可采用與本配點(diǎn)法相同的方法實(shí)現(xiàn)。

6 彈丸侵徹靶板算例

根據(jù)理論編制了非線(xiàn)性侵徹動(dòng)力過(guò)程數(shù)值分析的RKPM 法,對(duì)具體算例進(jìn)行計(jì)算,驗(yàn)證上述算法。

鋼質(zhì)圓柱型彈丸直徑15 mm,長(zhǎng)度50 mm,周邊固定的鋼板厚度300 mm,彈丸以1 km/s 速度、與鋼板表面法線(xiàn)方向傾斜角θ=10°入射鋼板。狀態(tài)方程參數(shù)分別為;本構(gòu)和損傷參數(shù)分別為:A=790 M Pa,B=510 M Pa,C=0.015,n=0.27,m=1.05,=0.47×10-5,D2=-9.0,D3=3.0,D4=0,D5=0.78,溶化溫度Tm=1.8 kK,初始溫度T0=293 K。

由于具有對(duì)稱(chēng)性,取彈丸和靶板的一半作為分析模型,并劃分成質(zhì)點(diǎn),作為計(jì)算基點(diǎn)。圖3 給出彈丸侵徹前后的質(zhì)點(diǎn)動(dòng)態(tài)變化圖,可以看出靶板在彈丸侵徹過(guò)程中破碎和彎曲變形。

圖3 彈丸與靶板侵徹過(guò)程的質(zhì)點(diǎn)動(dòng)態(tài)圖Fig.3 The dynamic particles' photo of pellet into steel plate

圖4 ～5 為計(jì)算給出的彈丸和靶板不同時(shí)刻應(yīng)力分布云圖?？梢钥闯?t=30 μs 時(shí)彈丸端頭部出現(xiàn)明顯的破壞斷裂,局部應(yīng)力達(dá)到1.2 GPa 以上。靶板應(yīng)力成近乎圓形的分布規(guī)律向四周擴(kuò)散。

圖4 彈丸應(yīng)力分布Fig.4 Nephogram of stress of pellet

圖5 靶板應(yīng)力分布Fig.5 Nephogram of stress of steel plate

圖6 彈丸尾部中心點(diǎn)水平速度和加速度Fig.6 Velocity and acceleration at the center point of end part of pellet

圖6 分別給出了彈丸尾部中心點(diǎn)水平方向的速度和加速度變化圖,彈丸經(jīng)過(guò)侵徹后水平方向的速度基本減至0。圖7 計(jì)算給出了彈丸端頭部中心點(diǎn)的速度變化和加速度變化圖。

選擇初始速度分別為0.6、1.0、1.2 和1.5 km/s 進(jìn)行計(jì)算,并將彈體入射后的剩余速度與理論解[10]進(jìn)行比較,結(jié)果如圖8,兩者完全一致。

圖7 彈丸端頭部中心點(diǎn)垂直速度和加速度Fig.7 Vertical velocity and acceleration at the center point of head part of pellet

圖8 彈體剩余速度Fig.8 Residual velocity of pellet

7 結(jié) 論

利用再生核質(zhì)點(diǎn)法計(jì)算模擬了侵徹過(guò)程。在侵徹階段的計(jì)算中,將彈丸和靶板劃分核質(zhì)點(diǎn),不需要?jiǎng)澐殖蓡卧?利用RKPM 方法對(duì)侵徹過(guò)程進(jìn)行大應(yīng)變、高應(yīng)變率計(jì)算,避免了有限單元法網(wǎng)格重分問(wèn)題,減少了計(jì)算的難度,提高了計(jì)算的精度和速度。此外,為了方便解決RKPM 方法中的滑移面接觸和邊界條件問(wèn)題,提出了新的滑移面算法和配點(diǎn)法。數(shù)值模擬結(jié)果表明,本方法易于編程,計(jì)算精度較高。

[1] 宋順成,才鴻年.彈丸侵徹混凝土的SPH 算法[J].爆炸與沖擊,2003,23(1):56-60.SONG Shun-cheng, CAI H ong-nian.SPH algorithm for projectile penetrating into concrete[J].Explosion and Shock Waves,2003,23(1):56-60.

[2] 宋順成,李國(guó)斌,才鴻年,等.戰(zhàn)斗部對(duì)混凝土先侵徹后爆轟的數(shù)值模擬[J].兵工學(xué)報(bào),2006,27(2):230-234.SONG Shun-cheng,LI Guo-bin,CAI Hong-nian,et al.Numerical simulation of penetration-then-detonation of concrete target with projectile[J].Acta Armamentarii,2006,27(2):230-234.

[3] Lucy L B.A numerical approach to the testing of the fission hypothesis[J].The Astronomical Journal, 1977,82:1013-1024.

[4] H ayhurst C J,Clegg R A.Cylindrically symmetric SPH simulations of hypervelocity impacts on thin plates[J].International Journal of Impact Engineering, 1997,20:337-348.

[5] Dancygier A N, Yankelevsky D Z.High strength concrete response to hard projectile impact[J].International Journal of Impact Engineering,1996,18(6):583-599.

[6] Liu W K,Jun S,Zhang Y F.Reproducing kernel particle method[J].International Journal for Numerial Methods in Fluids, 1995,20(8-9):1081-1106.

[7] Liu W K, Jun S,Li S F,et al.Reproducing kernel particle methods for structural dynamics[J].International Journal for Numerial Methods in Engineering,1995,38(10):1655-1679.

[8] Holmquist T J, Johnson G R, Cook W H.A computational constitutive model for concrete subjective to large strain,high strain rates, and high pressure[C]∥14th International Symposium on Ballistics.USA:American Defense Preparedness Association, 1993:591-600.

[9] 趙光明,宋順成.高速?zèng)_擊過(guò)程數(shù)值分析的再生核質(zhì)點(diǎn)法[J].力學(xué)學(xué)報(bào),2007,39(1):63-69.ZH AO Guang-ming, SONG Shun-cheng.The reproducing kernel particle method for numerical analysis of highspeed impact process[J].Chinese Journal of Theoretical and Applied Mechanics, 2007,39(1):63-69.

[10] 馬曉青.沖擊動(dòng)力學(xué)[M].北京:北京理工大學(xué)出版社,1991.