張鳳國(guó)，韓 冰，何長(zhǎng)江

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

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

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

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

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

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

圖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 單元的侵蝕計(jì)算

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

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

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

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

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

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

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

至此，本文中將滑移線的侵蝕標(biāo)準(zhǔn)分為3 類，相應(yīng)的單元是否被侵蝕掉將根據(jù)實(shí)際情況分別考慮。如圖2 中的單元A 和J 被侵蝕掉，則滑移線由原來(lái)的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 侵徹問(wèn)題的數(shù)值模擬

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

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

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

圖3 初始計(jì)算模型Fig.3 The initial computational model

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

4 結(jié) 論

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

［1］ 張鳳國(guó)，李維新，洪滔，等.超高速鎢合金長(zhǎng)桿彈對(duì)混凝土侵徹及損傷破壞的數(shù)值分析［J］.彈道學(xué)報(bào)，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］ 張鳳國(guó).動(dòng)載荷作用下混凝土靶板損傷破壞的數(shù)值分析［J］.兵工學(xué)報(bào)，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］ 張鳳國(guó)，馮其京，郝鵬程，等.聚能裝藥侵徹混凝土靶板的數(shù)值模擬［J］.計(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.