侵蝕滑移計算方法的改進*

張鳳國，韓 冰，何長江

(北京應(yīng)用物理與計算數(shù)學(xué)研究所，北京100094)

采用拉格朗日有限元程序計算多塊材料相接觸問題時，如高速碰撞和侵徹問題，在接觸面附近往往形成高壓、大變形區(qū)，滑移線的處理對于接觸問題的數(shù)值計算結(jié)果以及相關(guān)物理現(xiàn)象的顯示至關(guān)重要。二維拉格朗日有限元程序一般將密度較大的物質(zhì)定義為主塊，因為它具有較大的慣性，而相對應(yīng)的物質(zhì)定義為從塊，在接觸面上的節(jié)點則分別按順序分成主點和從點，從而構(gòu)成1 組滑移線。在計算過程中，當(dāng)從點侵入到主塊中時，把該節(jié)點沿與接觸面垂直方向拉回到接觸面上，同時按動量、動量矩守恒以及速度協(xié)調(diào)原則調(diào)整相關(guān)節(jié)點的速度。對于與主點相關(guān)單元的侵蝕計算問題，一般的處理方法是以累積塑性應(yīng)變?yōu)榕袚?jù)，即當(dāng)單元內(nèi)的累積塑性應(yīng)變大于給定值時，單元內(nèi)的應(yīng)力、應(yīng)變重零，單元失效。采用這種簡單的處理方法，一方面常常造成計算數(shù)值的不穩(wěn)定性，即可能有從點侵入到主塊中;另一方面在顯示物理現(xiàn)象時，主塊和從塊存在脫離現(xiàn)象以及接觸面不光滑。

二維彈塑性流體力學(xué)拉格朗日有限元程序LTZ-2D 在計算、分析侵徹問題中得到了較好的應(yīng)用［1-3］，本文中將LTZ-2D 程序中侵蝕滑移的計算方法進行改進，使程序在模擬侵徹問題的物理現(xiàn)象方面有所改進。

1 主、從節(jié)點速度的調(diào)整

目前滑移線主、從節(jié)點間動量轉(zhuǎn)換的方法有很多種，針對LTZ-2D 程序的特點，本文中采用G.R.Johnson 等［4］的方法。如圖1 所示，當(dāng)從點s 由A 點侵入到由主點i 和j 組成的主滑移線內(nèi)部B 點時，沿i、j 線的法線方向?qū) 點拉回到i、j 線上的C 點，并對i、j 和s 節(jié)點的法線方向的速度進行調(diào)整

式中:ri=0，rj和rs分別為節(jié)點j 和s 相對于節(jié)點i在r 方向上的距離，tl和ta分別為原來的動量和角動量，mi、mj和ms為節(jié)點質(zhì)量，vi、vj和vs為速度調(diào)整后3 個節(jié)點上的法向速度。

圖1 滑移線的處理Fig.1 The method for conserving normal momentum

求解方程組(1)，可以得到

式中:a=rs/rj，b=1－a，d=mj+ams，e=mjrj+amsrs，g=mi+bms。

2 單元的侵蝕計算

侵蝕滑移的計算，包括侵蝕單元的處理，是拉格朗日有限元程序在處理大變形問題時所采取的主要方法之一。G.R.Johnson 等［5］給出了基于累積塑性應(yīng)變?yōu)榕袛鄻?biāo)準(zhǔn)的單元侵蝕計算方法，即滑移界面附近的單元的等效塑性應(yīng)變超過某一給定值時，認為該單元是失效元，其正應(yīng)力、剪應(yīng)力置零，等效應(yīng)變率也置零，并從所屬物質(zhì)結(jié)構(gòu)中刪除該單元。目前此方法在有限元程序中應(yīng)用較廣泛，J.J.Pyun等［6］對此方法作了一定的改進。

本文中將滑移線上單元分為3 類，如圖2 所示，1 ～9 為滑移線的主點，A ～M 為與滑移線主點相關(guān)的單元，對不同種類的單元采用不同的處理方法。

首先，單元只包含滑移線上的1 個節(jié)點，如圖2 中的B、D、G、H、L 單元;或者雖然單元包含滑移線上的2 個節(jié)點，但單元的邊與滑移線不重合，如單元I，則與之相對應(yīng)的滑移線節(jié)點包括2、3、5、5、8、5 或7。這類單元不允許被侵蝕掉，即使單元的累積塑性應(yīng)變超過給定值，否則將造成滑移線數(shù)據(jù)的不規(guī)則，如單元G 被侵蝕掉后，滑移線的節(jié)點順序?qū)⒆兂?-2-3-4-5-11-12-5-6-7-8-9，物質(zhì)界面內(nèi)部形成了1 個空腔，這種順序的滑移線在拉格朗日程序計算過程中很可能造成數(shù)值計算的不穩(wěn)定性［7］。

其次，單元包含滑移線上的2 個節(jié)點，并且單元有1 條邊與滑移線重合，如圖2 中的A、C、E、F、K、M單元，與之相對應(yīng)的滑移線節(jié)點為(1，2)、(2，3)、(3，4)、(4，5)、(7，8)和(8，9)，當(dāng)這類單元的累積塑性應(yīng)變超過給定值，并且單元相對于滑移線的角度大于某一值或單元的量綱一硬化因子大于某一值時，單元允許被侵蝕掉。

角度的判據(jù)一般取100°～160°之間，這主要根據(jù)彈體的初始速度所確定，如圖2 中的單元A，相對于滑移線的角1-10-2 大于100°，若單元的累積塑性應(yīng)變超過給定值，則單元A 被侵蝕掉，而單元E 則不能被侵蝕，即使單元的累積塑性應(yīng)變超過給定值，因為其相對于滑移線的角3-11-4 為銳角。

三角形單元的硬化因子是指單元的最小高與最長邊的比值，它是判斷三角形單元形狀好壞的重要判據(jù)，單元的量綱一硬化因子是指單元當(dāng)前的硬化因子與單元初始硬化因子的比值，一般情況下取值為0.1，即當(dāng)單元的累積塑性應(yīng)變超過給定值，并且其量綱一硬化因子小于0.1 時，該單元被侵蝕掉。

最后，單元的3 個節(jié)點依次為滑移線上的節(jié)點，此時單元有2 條邊與滑移線重合，如圖2 中的單元J，其相對應(yīng)的節(jié)點為(5，6，7)，當(dāng)單元的累積塑性應(yīng)變超過給定值，或單元的無量綱硬化因子小于0.1，或相對于滑移線的角5-6-7 小于90°時，這類單元將被允許侵蝕掉。

至此，本文中將滑移線的侵蝕標(biāo)準(zhǔn)分為3 類，相應(yīng)的單元是否被侵蝕掉將根據(jù)實際情況分別考慮。如圖2 中的單元A 和J 被侵蝕掉，則滑移線由原來的1-2-3-4-5-6-7-8-9，改變?yōu)?-10-3-4-5-7-8-9。

圖2 主滑移線上的單元種類Fig.2 Classification of the elements along the slideline

3 侵徹問題的數(shù)值模擬

將新的滑移線侵蝕判斷標(biāo)準(zhǔn)引入到LTZ-2D 有限元程序中，對彈體的侵徹問題進行數(shù)值模擬，并與以前的計算結(jié)果進行比較。彈體材料為OFHC 銅，彈體半徑為12.7 mm，彈體長度為50.7 mm，彈體初速度為2 540 m/s，靶材料為4340 鋼，靶厚50.8 mm，靶體半徑為76.2 mm，計算采用的是軸對稱計算模型，圖3 顯示的是初始計算模型。

圖4(a)～4(b)給出了10 μs 時刻的計算結(jié)果，圖4(a)顯示了采用以前的滑移侵蝕判斷標(biāo)準(zhǔn)所引起的接觸面處彈體表面出現(xiàn)的凹凸不平的不規(guī)則現(xiàn)象，而采用新的標(biāo)準(zhǔn)很好地解決了這一問題，彈體表面顯示了很好的圓滑曲線，見圖4(b)。

圖4(c)顯示了采用以前的滑移侵蝕標(biāo)準(zhǔn)所引起的另一種數(shù)值模擬的物理現(xiàn)象，即計算過程中顯示的物理現(xiàn)象存在彈靶脫離的情況，這顯然與實際物理現(xiàn)象存在較大差別，而采用新標(biāo)準(zhǔn)的計算結(jié)果顯示彈靶可以保持較好的接觸，見圖4(d)。此外，圖4(b)、4(d)顯示，采用新的滑移侵蝕判斷標(biāo)準(zhǔn)，彈靶的接觸表面具有較好的光滑性。

圖3 初始計算模型Fig.3 The initial computational model

圖4 用不同的方法在不同時刻的沖擊計算結(jié)果Fig.4 Impact calculation by different methods at different times

4 結(jié) 論

對LTZ-2D 有限元程序中的滑移侵蝕算法進行了改進。計算結(jié)果顯示，采用新的滑移侵蝕判斷標(biāo)準(zhǔn)，在計算彈體對靶板的侵徹過程中，彈靶之間不僅可以基本保持接觸，而且接觸面更光滑，與真實的物理現(xiàn)象更接近，對侵徹過程的數(shù)值描述更精細，同時也有助于深入分析彈靶之間相互作用的力學(xué)機理。

［1］ 張鳳國，李維新，洪滔，等.超高速鎢合金長桿彈對混凝土侵徹及損傷破壞的數(shù)值分析［J］.彈道學(xué)報，2008，20(3):64-67.ZHANG Feng-guo，LI Wei-xin，HONG Tao，et al.Numerical simulation for damage and penetration of concrete driven by long-rod projectile of tungsten alloy under super-high speed［J］.Journal of Ballistics，2008，20(3):64-67.

［2］ 張鳳國.動載荷作用下混凝土靶板損傷破壞的數(shù)值分析［J］.兵工學(xué)報，2009，30(9):19-22.ZHANG Feng-guo.Numerical analysis of the damage of concrete under dynamic loading［J］.Acta Armamentrii，2009，30(9):19-22.

［3］ 張鳳國，馮其京，郝鵬程，等.聚能裝藥侵徹混凝土靶板的數(shù)值模擬［J］.計算物理，2009，26(6):887-891.ZHANG Feng-guo，F(xiàn)ENG Qi-jing，HAO Peng-cheng，et al.Numerical study on penetration of shaped charge jets into concrete targets［J］.Chinese Journal of Computational Physics，2009，26(6):887-891.

［4］ Johnson G R，Stryk R A.User instructions for the EPIC-2 code［R］.AFATL-TR-86-51，1986.

［5］ Johnson G R，Stryk R A.Eroding interface and improved tetrahedral element algorithm for high-velocity impact computations in three dimensions［J］.International Journal of Impact Engineering，1987，5(1/2/3/4):411-421.

［6］ Pyun J J，Kennedy C M，Hruska D.A new slideline/eroding algorithm for EPIC2［J］.International Journal of Impact Engineering，1990，10(1/2/3/4):473-482.

［7］ Pyun J J，Hamilton C W.A slip/collapse algorithm in a hydro code［C］∥Nichols J A.Proceedings of the Nuclear Explosive Design Physics Conference.New Mexico:Los Alamos National Laboratory，1983:127-150.