非線性侵徹動力過程的再生核質點法*

趙光明,宋順成

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

1 引 言

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

侵徹過程的數值模擬,在很大程度上取決于材料的本構模型。本文中,通過引入Johnson-Cook 本構模型及損傷模型,采用RKPM 無網格法實現高速侵徹過程數值模擬,同理論解的結果進行比較,驗證數值模擬方法的有效性。

2 RKPM 形函數

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

式中:x 是空間點,s 是領域內的點。在Dirac Delta 條件下,很難解決離散化問題,為此,進行修改[6]

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

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

式中

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

所以

式中

將再生方程離散為

式中:xi 表示質點,n 是計算域內質點數量,ΔV i 表示質點x i 的體積,Ψi 稱為再生核質點法的形函數。

3 侵徹材料的本構關系

根據Johnson-Cook 本構方程[8],材料的Mises 等效應力

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

彈丸材料的狀態方程一般可以通過Grüneisen 狀態方程表示

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

根據Johnson-Cook 模型,單元內材料的損傷演化定義為

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

4 滑移面算法

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

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

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

圖1 滑移面Fig.1 Sliding surface

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

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

切向方向速度不發生改變。

對于三維問題計算,首先以質點為頂點,將主滑移面劃分成相連的三角形區域。當主、從滑移面發生接觸時,判斷從滑移面上的某一質點處于主滑移面哪一個三角形區域,并按照二維問題類似的標準判斷,將該從接觸質點與三角形頂點的3 個質點發生速度調整。

5 速度配點法

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

給定質點速度

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

簡記為

則實現方法為

本質邊界條件及加速度邊界條件也可采用與本配點法相同的方法實現。

6 彈丸侵徹靶板算例

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

鋼質圓柱型彈丸直徑15 mm,長度50 mm,周邊固定的鋼板厚度300 mm,彈丸以1 km/s 速度、與鋼板表面法線方向傾斜角θ=10°入射鋼板。狀態方程參數分別為;本構和損傷參數分別為: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。

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

圖3 彈丸與靶板侵徹過程的質點動態圖Fig.3 The dynamic particles' photo of pellet into steel plate

圖4 ～5 為計算給出的彈丸和靶板不同時刻應力分布云圖。可以看出,t=30 μs 時彈丸端頭部出現明顯的破壞斷裂,局部應力達到1.2 GPa 以上。靶板應力成近乎圓形的分布規律向四周擴散。

圖4 彈丸應力分布Fig.4 Nephogram of stress of pellet

圖5 靶板應力分布Fig.5 Nephogram of stress of steel plate

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

圖6 分別給出了彈丸尾部中心點水平方向的速度和加速度變化圖,彈丸經過侵徹后水平方向的速度基本減至0。圖7 計算給出了彈丸端頭部中心點的速度變化和加速度變化圖。

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

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

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

7 結 論

利用再生核質點法計算模擬了侵徹過程。在侵徹階段的計算中,將彈丸和靶板劃分核質點,不需要劃分成單元,利用RKPM 方法對侵徹過程進行大應變、高應變率計算,避免了有限單元法網格重分問題,減少了計算的難度,提高了計算的精度和速度。此外,為了方便解決RKPM 方法中的滑移面接觸和邊界條件問題,提出了新的滑移面算法和配點法。數值模擬結果表明,本方法易于編程,計算精度較高。

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