趙光明,宋順成

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

1 引 言

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

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

2 RKPM 形函數(shù)

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

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

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

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

式中

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

所以

式中

將再生方程離散為

式中:xi 表示質(zhì)點(diǎn),n 是計算域內(nèi)質(zhì)點(diǎn)數(shù)量,ΔV i 表示質(zhì)點(diǎn)x i 的體積,Ψi 稱為再生核質(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)方程一般可以通過Grüneisen 狀態(tài)方程表示

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

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

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

4 滑移面算法

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

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

從滑移面質(zhì)點(diǎn)A 移動后可能位于主滑移面的3 種位置,即A 可能在Ⅰ、Ⅱ或Ⅲ區(qū)域內(nèi)點(diǎn)′ 或。如果在Ⅰ內(nèi)位置,說明從質(zhì)點(diǎn)移動后在滑移面之外,不發(fā)生接觸,此時質(zhì)點(diǎn)間的各種物理量不需要發(fā)生交換;如果在Ⅱ、Ⅲ內(nèi),說明從質(zhì)點(diǎn)已侵入主滑移面內(nèi),按不侵入條件,則立即將質(zhì)點(diǎn)A垂直拉到滑移面上。對于區(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ì)量和線段長度,v+、v-分別是拉動后、前的質(zhì)點(diǎn)的法向速度?？梢宰C明,式(18)～(20)既保證了滑移面質(zhì)點(diǎn)的動量守恒,也保證了動量矩守恒。如果質(zhì)點(diǎn)A 位移后在區(qū)域Ⅲ內(nèi)A′3位置,此時還須檢查相鄰主滑移面質(zhì)點(diǎ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ā)生改變。

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

5 速度配點(diǎn)法

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

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

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

簡記為

則實現(xiàn)方法為

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

6 彈丸侵徹靶板算例

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

鋼質(zhì)圓柱型彈丸直徑15 mm,長度50 mm,周邊固定的鋼板厚度300 mm,彈丸以1 km/s 速度、與鋼板表面法線方向傾斜角θ=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。

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

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

圖4 ～5 為計算給出的彈丸和靶板不同時刻應(yīng)力分布云圖?？梢钥闯?t=30 μs 時彈丸端頭部出現(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)過侵徹后水平方向的速度基本減至0。圖7 計算給出了彈丸端頭部中心點(diǎn)的速度變化和加速度變化圖。

選擇初始速度分別為0.6、1.0、1.2 和1.5 km/s 進(jìn)行計算,并將彈體入射后的剩余速度與理論解[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)法計算模擬了侵徹過程。在侵徹階段的計算中,將彈丸和靶板劃分核質(zhì)點(diǎn),不需要劃分成單元,利用RKPM 方法對侵徹過程進(jìn)行大應(yīng)變、高應(yīng)變率計算,避免了有限單元法網(wǎng)格重分問題,減少了計算的難度,提高了計算的精度和速度。此外,為了方便解決RKPM 方法中的滑移面接觸和邊界條件問題,提出了新的滑移面算法和配點(diǎn)法。數(shù)值模擬結(jié)果表明,本方法易于編程,計算精度較高。

[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] 宋順成,李國斌,才鴻年,等.戰(zhàn)斗部對混凝土先侵徹后爆轟的數(shù)值模擬[J].兵工學(xué)報,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] 趙光明,宋順成.高速沖擊過程數(shù)值分析的再生核質(zhì)點(diǎn)法[J].力學(xué)學(xué)報,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] 馬曉青.沖擊動力學(xué)[M].北京:北京理工大學(xué)出版社,1991.