陳成軍， 陳小偉， 柳明

(1.中國(guó)工程物理研究院 總體工程研究所，四川,綿陽(yáng) 621900；2.北京理工大學(xué) 前沿交叉科學(xué)研究院，北京 100081；3.北京理工大學(xué) 爆炸科學(xué)與技術(shù)國(guó)家重點(diǎn)實(shí)驗(yàn)室，北京 100081)

高效、健壯、精確的侵蝕接觸算法是彈體穿甲/侵徹、高速碰撞等沖擊動(dòng)力學(xué)問(wèn)題數(shù)值模擬的關(guān)鍵，是拉格朗日顯式有限元方法研究與發(fā)展中的重要課題之一. 所謂“侵蝕接觸”是指在接觸-碰撞過(guò)程中，由于材料損傷破壞而刪除相應(yīng)單元導(dǎo)致新物質(zhì)界面產(chǎn)生，新界面繼而與附近物質(zhì)發(fā)生接觸作用. 與一般的動(dòng)態(tài)接觸問(wèn)題相比，侵蝕接觸計(jì)算存在很大挑戰(zhàn)：接觸點(diǎn)與接觸片隨時(shí)間動(dòng)態(tài)變化，原有接觸點(diǎn)/接觸片可能隨單元?jiǎng)h除而消失，伴隨產(chǎn)生新的接觸點(diǎn)/接觸片，接觸對(duì)在空間和時(shí)間上很難預(yù)測(cè)；接觸界面不平滑，接觸片不規(guī)則，點(diǎn)-面、邊-面、邊-邊等多種接觸形式共存，接觸對(duì)和接觸方向確定困難，數(shù)值模擬中易發(fā)生非物理的界面穿透.

迄今已有很多學(xué)者開(kāi)展了動(dòng)態(tài)接觸算法研究[1-2]，但現(xiàn)有研究大多針對(duì)一般動(dòng)態(tài)接觸問(wèn)題開(kāi)展，侵蝕接觸算法的研究還很少. BELYTSCHKO等[3]針對(duì)彈體穿甲問(wèn)題以向量組裝法自動(dòng)識(shí)別/生成接觸界面而構(gòu)建了侵蝕接觸算法，但沒(méi)有討論侵蝕接觸中常遇到的接觸對(duì)遺漏、多重接觸等問(wèn)題. WHIRLEY等[4]基于單元拓?fù)潢P(guān)系自動(dòng)識(shí)別接觸界面，采用“園管”和“小球”處理邊、角接觸，建立了自動(dòng)接觸搜索算法，但沒(méi)有分析算法在侵蝕接觸中應(yīng)用. HEINSTEIN等[5]詳細(xì)討論了接觸對(duì)確定、接觸點(diǎn)計(jì)算及接觸力施加方向不準(zhǔn)確等問(wèn)題，建立了用于接觸對(duì)確定的靜、動(dòng)態(tài)檢測(cè)算法；針對(duì)侵蝕接觸給出了一種接觸面自動(dòng)生成/識(shí)別算法；采用接觸強(qiáng)度方法處理多重接觸問(wèn)題，但該方法在處理尖角附近的接觸時(shí)仍然可能會(huì)失效. 文中作者基于接觸片間的面-面幾何關(guān)系建立了可用于侵蝕接觸計(jì)算的接觸對(duì)確定方法[6]，但在接觸深度計(jì)算中引入了很大的簡(jiǎn)化導(dǎo)致在某些復(fù)雜問(wèn)題中計(jì)算精度降低. ZAVARISE等[7]討論了二維準(zhǔn)靜態(tài)接觸問(wèn)題中接觸對(duì)確定的三種歧義性并建議了多種處理算法. 然而，這些算法用于三維問(wèn)題分析存在很大困難. 綜之，目前為止侵蝕接觸算法的研究還不充分，侵蝕接觸計(jì)算精度仍需提高.

侵蝕接觸計(jì)算可分解為4個(gè)步驟：接觸界面自動(dòng)生成/識(shí)別、接觸對(duì)全局搜索、接觸對(duì)局部搜索和接觸約束施加. 其中，接觸對(duì)全局搜索和接觸約束施加算法與一般動(dòng)態(tài)接觸問(wèn)題沒(méi)有區(qū)別，具體算法可參閱相關(guān)文獻(xiàn)[5,8,9]，這里不再贅述. 文中將重點(diǎn)討論侵蝕接觸計(jì)算中的接觸對(duì)局部搜索方法，并給出接觸界面的自動(dòng)生成/識(shí)別算法. 需要說(shuō)明的是，接觸界面的點(diǎn)-面離散模型[9]仍是目前顯式有限元程序最廣泛采用的界面離散方法，而六面體單元是最常用的單元類(lèi)型，因此文中侵蝕接觸算法面向點(diǎn)-面離散模型、六面體單元間的接觸計(jì)算設(shè)計(jì). 文中接觸約束施加采用罰函數(shù)法.

1 接觸對(duì)確定

1.1 接觸對(duì)預(yù)搜索

接觸對(duì)全局搜索給出當(dāng)前及隨后若干時(shí)間步內(nèi)系統(tǒng)中可能發(fā)生接觸的部位(即接觸測(cè)試對(duì))，因此接觸測(cè)試對(duì)集合通常遠(yuǎn)大于當(dāng)前時(shí)間步內(nèi)真實(shí)發(fā)生接觸的(從點(diǎn),主片)集合. 接觸對(duì)局部搜索的目的就是從接觸測(cè)試對(duì)中精確確定出當(dāng)前時(shí)間步的真實(shí)接觸對(duì). 為提高計(jì)算效率，在接觸對(duì)精確確定前引入“接觸預(yù)搜索”以快速排除當(dāng)前時(shí)間步內(nèi)明顯不會(huì)發(fā)生接觸的測(cè)試對(duì). 這里以“范圍檢測(cè)”和“法線(xiàn)檢測(cè)”實(shí)現(xiàn)接觸預(yù)搜索.

① 范圍檢測(cè). 當(dāng)前時(shí)間步，只有位于主片包圍盒內(nèi)的從點(diǎn)才可能與該主片接觸. 主片的包圍盒定義為以主片面心x0為球心，以長(zhǎng)度l為半徑的球形區(qū)域：

D={x|(x-x0)2

(1)

式中,l=max(l31,l42)，l31、l42為主片對(duì)角線(xiàn)的長(zhǎng)度，下標(biāo)1～4為主片節(jié)點(diǎn)的局部編號(hào).

② 法線(xiàn)檢測(cè). 當(dāng)從點(diǎn)與主片的外法線(xiàn)方向滿(mǎn)足如下關(guān)系時(shí)，從點(diǎn)與主片才可能發(fā)生接觸：

nm·ns<0

(2)

這里主片外法線(xiàn)方向nm近似為

(3)

從點(diǎn)法線(xiàn)方向ns定義為包含該節(jié)點(diǎn)的所有接觸片外法線(xiàn)方向的加權(quán)平均值，即

(4)

式中：m為與從點(diǎn)相連的單元外表面數(shù)量;Ai為第i個(gè)單元外表面的面積.

1.2 接觸對(duì)的確定

對(duì)于一個(gè)接觸測(cè)試對(duì)，如果從點(diǎn)與主片滿(mǎn)足：

① 從點(diǎn)在主片上的投影點(diǎn)位于該主片內(nèi)；

② 從點(diǎn)與主片間的間隙函數(shù)值≤0.

則該測(cè)試對(duì)在當(dāng)前時(shí)間步內(nèi)為真實(shí)接觸對(duì)；反之，從點(diǎn)與主片不會(huì)發(fā)生接觸.

理想狀態(tài)(單元規(guī)則且變形小)下，從點(diǎn)在主片上的投影點(diǎn)計(jì)算很簡(jiǎn)單，可利用牛頓迭代法求解下列方程組得到[4,10]：

(5)

但在沖擊動(dòng)力學(xué)問(wèn)題中，單元常發(fā)生很大的變形、轉(zhuǎn)動(dòng)甚至扭曲畸變，接觸片常變得很不規(guī)則，此時(shí)該方法存在計(jì)算量大、不穩(wěn)定、收斂困難甚至不收斂等問(wèn)題[8].

為克服牛頓迭代法計(jì)算投影點(diǎn)的困難，文中建立非迭代計(jì)算方法. 首先將接觸從點(diǎn)、接觸主片都投影到某個(gè)平面上，即將三維空間問(wèn)題轉(zhuǎn)化為二維平面問(wèn)題. 投影面的構(gòu)造如圖1所示，投影面的基為

圖1 投影面的構(gòu)造

e′2=e′3×e′1

(6)

假定從點(diǎn)在投影面上的坐標(biāo)為(x′,y′)，主片節(jié)點(diǎn)在投影面上的坐標(biāo)為(x′i,y′i)，(i= 1～4). 以雙線(xiàn)性插值描述主片，并建立局部坐標(biāo)系統(tǒng)(r,s). 如果從點(diǎn)位于主片內(nèi)，則從點(diǎn)坐標(biāo)可表達(dá)為

(7)

(8)

式中:(r,s)為從點(diǎn)的局部坐標(biāo)，(ri,si)為主片第i個(gè)節(jié)點(diǎn)的局部坐標(biāo). 將式(7)中消去s，可得關(guān)于r的二次方程：

Ar2+Br+C=0

(9)

這里A、B、C為與節(jié)點(diǎn)坐標(biāo)相關(guān)的常數(shù)

A=a2b4-a4b2

B=a2b3-a3b2+a4(y′-b0)-b4(x′-a0)

C=a3(y′-b0)-b3(x′-a0)

(10)

以上各式中

a2=(-x′1+x′2+x′3-x′4)/4

a3=(x′3+x′4-x′1-x′2)/4

a4=(x′1-x′2+x′3-x′4)/4

(11)

式中,b0、b2等常數(shù)類(lèi)似地由y′i給出. 利用二次方程的求根公式可計(jì)算出r,然后將r代入式(7)便得到s. 若計(jì)算出的從點(diǎn)局部坐標(biāo)r、s位于區(qū)間[-1,1]，則從點(diǎn)對(duì)主片的投影就位于該主片內(nèi)部.

知道了從點(diǎn)相對(duì)于主片的局部坐標(biāo)(r,s)，便可容易地計(jì)算出從點(diǎn)與主片間的距離，即間隙函數(shù)值：

g=[xs-xc(r,s)]·nm

(12)

式中：xs為從點(diǎn)的坐標(biāo);xc為接觸點(diǎn)(即投影點(diǎn))的坐標(biāo).

1.3 特殊情況的處理

如果接觸片規(guī)則、接觸界面足夠光滑，上述接觸對(duì)確定方法可精確地給出接觸對(duì)及其接觸方向. 然而在侵蝕接觸問(wèn)題中，單元失效刪除前通常會(huì)發(fā)生很大的變形或畸變，導(dǎo)致接觸片很不規(guī)則；單元?jiǎng)h除后會(huì)產(chǎn)生新的不規(guī)則接觸片并存在大量尖角，導(dǎo)致接觸界面極不平滑. 因此，侵蝕接觸問(wèn)題的數(shù)值模擬常存在接觸對(duì)遺漏或接觸方向不準(zhǔn)確等問(wèn)題，致使數(shù)值計(jì)算精度下降甚至錯(cuò)誤.

1.3.1接觸對(duì)遺漏

接觸對(duì)遺漏包含兩種情況. 第一種情況如圖2(a)所示,即邊界盲區(qū)問(wèn)題，從點(diǎn)與主片邊界接觸，但數(shù)值計(jì)算出的投影點(diǎn)位于主片外. 這種情況主要是由于局部搜索算法引入的某些近似或者數(shù)值截?cái)嗾`差導(dǎo)致的. 此時(shí)可采用處理方法：引入一定的容差ε，適當(dāng)放寬式(12)中投影點(diǎn)局部坐標(biāo)(r,s)的取值范圍：

圖2 特殊情況的二維示意

|r|<1+ε,|s|<1+ε

(13)

數(shù)值計(jì)算經(jīng)驗(yàn)表明，取ε=0.01～0.02可很好地解決絕大多數(shù)的接觸對(duì)遺漏問(wèn)題. 該方法實(shí)質(zhì)上是對(duì)接觸片引入一定的擴(kuò)展域，抵消由于截?cái)嗾`差或投影點(diǎn)計(jì)算中近似處理(如不考慮接觸片的翹曲、接觸方向的近似等)引入的誤差.

接觸對(duì)遺漏的第二種情況如圖2(b)所示，就是所謂的內(nèi)盲區(qū)問(wèn)題[9]. 由于采用罰函數(shù)法施加接觸約束時(shí)，只有當(dāng)從點(diǎn)穿透主片時(shí)才會(huì)施加接觸力；另外，由于對(duì)接觸片引入了式(13)形式的擴(kuò)展，內(nèi)盲區(qū)問(wèn)題轉(zhuǎn)變?yōu)槎嘀亟佑|問(wèn)題.

1.3.2多重接觸

多重接觸即外盲區(qū)問(wèn)題[9]如圖2(c)所示，從點(diǎn)在多個(gè)主片上的投影點(diǎn)都位于相應(yīng)主片內(nèi)部，且其間隙函數(shù)值小于0. 此時(shí)，從點(diǎn)與多個(gè)主片都有可能接觸，傳統(tǒng)接觸搜索算法[9-10]很難確定出真實(shí)的接觸主片.

多重接觸是侵蝕接觸計(jì)算中的常見(jiàn)問(wèn)題. 已有接觸搜索算法[5,9,10]通常僅利用當(dāng)前幾何狀態(tài)量確定從點(diǎn)與主片的接觸關(guān)系，即使引入接觸強(qiáng)度[5]也不能確定出真實(shí)接觸對(duì). 如圖3所示尖角附近的接觸，從點(diǎn)上一時(shí)間步由位置A穿透主片S1運(yùn)動(dòng)到位置B，當(dāng)前時(shí)間步處于位置C. 如僅利用當(dāng)前構(gòu)形易得出從點(diǎn)穿透主片S2的錯(cuò)誤結(jié)論. 因此，文中引入一個(gè)狀態(tài)量跟蹤從點(diǎn)對(duì)主片的穿透歷史以便更好地確定接觸對(duì).

對(duì)每個(gè)從點(diǎn)I定義一個(gè)歷史狀態(tài)量state(I)記錄該從點(diǎn)前一時(shí)間步接觸的主片. 利用該狀態(tài)量可以很容易地處理圖3所示接觸對(duì)確定的困難：如果state(I) = 1，表示從點(diǎn)穿透主片S1，即從點(diǎn)與S1構(gòu)成接觸對(duì). 如果state(I) = -1，則表示前一時(shí)間步從點(diǎn)未與任何主片發(fā)生接觸，而在當(dāng)前時(shí)刻運(yùn)動(dòng)到位置C，如圖4所示. 此時(shí)可按以下邏輯順序和準(zhǔn)則確定接觸主片：

圖3 基于當(dāng)前構(gòu)形確定接觸對(duì)失效

① 計(jì)算從點(diǎn)對(duì)每個(gè)主片的侵徹深度(即間隙函數(shù)的負(fù)值)，選擇侵徹深度小的接觸片為該從點(diǎn)的接觸主片；

② 如果從點(diǎn)對(duì)每個(gè)主片的侵徹深度相等，進(jìn)一步計(jì)算接觸強(qiáng)度nm·ns，并選擇接觸強(qiáng)度小的接觸片為該從點(diǎn)的接觸主片；

③ 如果接觸強(qiáng)度相等，則從點(diǎn)是由區(qū)域Ⅲ中的某個(gè)位置運(yùn)動(dòng)到C的，接觸點(diǎn)位于兩個(gè)主片的公共邊或公共點(diǎn)上，如圖4所示的點(diǎn)M，接觸方向?yàn)镃M方向.

圖4 存在多個(gè)潛在接觸片時(shí)接觸對(duì)的確定

2 算法設(shè)計(jì)

2.1 接觸界面的自動(dòng)生成與識(shí)別

侵蝕接觸中，處于物體外表面的單元都有可能參與接觸，因此每次全局接觸搜索前需要首先確定出處于物體外表面的單元面作為全局搜索的輸入. 顯然，位于物體外表面的實(shí)體單元面只與一個(gè)單元相連，而位于內(nèi)部的單元面與兩個(gè)單元相連. 因此，可以利用單元面關(guān)聯(lián)的單元數(shù)量確定該單元面是否屬于外表面[4].

對(duì)于侵蝕接觸，由于存在單元?jiǎng)h除，單元拓?fù)浣Y(jié)構(gòu)生變化，位于外表面的單元面是動(dòng)態(tài)變化的，因此必須設(shè)計(jì)合理的數(shù)據(jù)結(jié)構(gòu)以滿(mǎn)足接觸片動(dòng)態(tài)增減的需要. 文中利用為靜態(tài)數(shù)組surface(I,i) 記錄系統(tǒng)中所有單元面所連接的單元個(gè)數(shù)，I∈[1,ne]，i∈[1,6]，這里ne為系統(tǒng)中的單元數(shù). 同時(shí)，對(duì)每個(gè)單元I定義一個(gè)狀態(tài)量cellFail(I) 表征該單元的激活狀態(tài)：如cellFail(I) = 1表示該單元失效，不參與接觸計(jì)算. 如某個(gè)單元J失效，則對(duì)應(yīng)的surface數(shù)組的J行所有元素值減1. 這樣根據(jù)surface元素值可以方便地識(shí)別出處于物體外表面的單元面，實(shí)現(xiàn)接觸面的自動(dòng)生成/識(shí)別.

2.2 計(jì)算流程

文中侵蝕接觸算法基于自主研發(fā)的并行顯式有限元程序PANDA-Impact進(jìn)行了實(shí)現(xiàn)，給出了每一個(gè)時(shí)間步內(nèi)接觸計(jì)算的偽代碼.

3 算例分析

本節(jié)利用典型算例驗(yàn)證分析文中所建立的侵蝕接觸算法的有效性. 算例一為經(jīng)典的彈性桿碰撞問(wèn)題[11]，用來(lái)驗(yàn)證接觸點(diǎn)計(jì)算方法與程序?qū)崿F(xiàn)的正確性；算例二為凹凸面間的接觸問(wèn)題[5]，分析文中算法處理多重接觸問(wèn)題的能力；第三個(gè)算例為彈體穿甲問(wèn)題[12]，該問(wèn)題涵蓋了結(jié)構(gòu)沖擊響應(yīng)數(shù)值模擬的所有非線(xiàn)性要素，如幾何大變形、材料非線(xiàn)性等，用于綜合驗(yàn)證文中侵蝕接觸算法的有效性.

3.1 彈性桿對(duì)撞

如圖5所示，方形截面桿A以1.0 mm/ms的速度對(duì)心撞擊靜止的方形截面桿B. 桿A與桿B截面積分別為1.0×106mm2、2.0×106mm2，長(zhǎng)度均為1.0×104mm. 假定兩桿材料均為彈性：桿A的彈性模量E=1.0×10-3MPa，B的彈性模量E=4.0×10-3MPa；兩桿密度均為ρ=1.0×10-3g/mm3，泊松比υ均取為0.

由于兩桿長(zhǎng)度遠(yuǎn)大于其截面尺寸，又泊松比取為0，因此可以利用一維應(yīng)力波理論計(jì)算出兩桿撞擊界面處的接觸力F和界面速度v：

(14)

兩桿都以20個(gè)單點(diǎn)積分六面體單元離散，見(jiàn)圖5；兩桿間的接觸-碰撞作用以文中建立的侵蝕接觸算法模擬. 圖6給出了兩桿碰撞過(guò)程中接觸力F隨時(shí)間t的變化曲線(xiàn)，圖7給出了接觸界面的速度曲線(xiàn). 作為對(duì)比，圖中同時(shí)列出了相應(yīng)的理論解(由式(14)計(jì)算得到). 可以看出，雖然接觸力和接觸界面速度存在一定程度的震蕩(這與OLDENBURG等[11]給出計(jì)算結(jié)果一致)，但其均值與理論解吻合良好. 總體上，數(shù)值方法較好地模擬了兩彈性桿的碰撞過(guò)程,初步驗(yàn)證了接觸點(diǎn)計(jì)算方法與程序?qū)崿F(xiàn)的正確性.

圖5 兩彈性桿對(duì)撞的有限元網(wǎng)格模型

圖6 接觸力時(shí)程曲線(xiàn)

圖7 接觸界面速度曲線(xiàn)

3.2 凹凸面接觸

如圖8所示，半徑為5.0 mm的半圓柱體在均布?jí)毫ψ饔孟孪蛳逻\(yùn)動(dòng)擠壓帶凹槽(凹槽半徑為5.0 mm)的方塊體. 方塊體底面及左右表面固支,半圓柱體上表面均布?jí)毫隨時(shí)間t變化：

圖8 凹凸面接觸計(jì)算模型

(15)

兩物體材料均為彈性：半圓柱體密度ρ=8.0×10-3g/mm3，彈性模量E=2.0×105MPa，泊松比υ=0.3；方塊體ρ=3.0×10-3g/mm3，E=7.0×104MPa，υ=0.3.

以單點(diǎn)積分六面體單元離散系統(tǒng)，網(wǎng)格剖分如圖9所示，半圓柱體接觸面(凸面)上的節(jié)點(diǎn)與方塊體接觸面(凹面)節(jié)點(diǎn)一一對(duì)應(yīng)(相應(yīng)節(jié)點(diǎn)坐標(biāo)相同). 該問(wèn)題模擬的難點(diǎn)在于，初始階段凸面上的接觸節(jié)點(diǎn)(頂點(diǎn)處的節(jié)點(diǎn)除外)處于凹面上接觸片的內(nèi)盲區(qū)；而對(duì)于凹面上的節(jié)點(diǎn)，它們則位于凸面上接觸片的外盲區(qū)，每個(gè)節(jié)點(diǎn)存在多個(gè)潛在接觸片. 該問(wèn)題的正確模擬不僅要求接觸搜索算法準(zhǔn)確確定出接觸對(duì)，更需要精確給出接觸點(diǎn)與接觸方向. 否則，系統(tǒng)中某些單元可能發(fā)生嚴(yán)重的畸變，半圓柱體甚至發(fā)生XY面內(nèi)的大幅度轉(zhuǎn)動(dòng)或沿母線(xiàn)方向的滑移[5].

圖9給出了文中算法的數(shù)值模擬結(jié)果(計(jì)算中應(yīng)力的單位為MPa，力的單位為N). 可以看出，半圓柱和方塊體上的單元變形均勻，系統(tǒng)的整體變形與應(yīng)力分布關(guān)于YOZ面完全對(duì)稱(chēng). 0.2 ms時(shí)刻半圓柱體接觸界面上的接觸力分布如圖10所示(力的單位為N)，各節(jié)點(diǎn)接觸力的大小與方向符合物理真實(shí). 加載過(guò)程中，半圓柱體承受的接觸力合力的如圖11所示，接觸力僅有Y向分量，圓柱體不會(huì)發(fā)生Z向滑動(dòng).

圖9 0.2 ms時(shí)刻系統(tǒng)的變形

圖10 0.2 ms時(shí)刻半圓柱體上的接觸力

圖11 半圓柱接觸力時(shí)程曲線(xiàn)

3.3 彈體侵徹

7.62 mm尖卵形子彈以827.5 mm/ms的初速度侵徹10 mm厚度的靶板[12]，如圖12所示，靶板長(zhǎng)和寬都為200 mm，四周固支. 子彈與靶板材料的動(dòng)態(tài)力學(xué)行為都以Johnson-Cook本構(gòu)模型描述為

圖12 計(jì)算模型

(16)

表1 子彈和靶板的材料參數(shù)[12]

考慮到系統(tǒng)的對(duì)稱(chēng)性，取1/2模型進(jìn)行建模分析. 彈體和靶板都以單點(diǎn)積分六面單元離散，子彈單元的特征尺寸約為1.0 mm，靶板中間區(qū)域單元的特征尺寸約為0.2 mm，系統(tǒng)總單元數(shù)為409 264.

圖13給出了子彈侵徹靶板過(guò)程中典型時(shí)刻的靶板變形及接觸狀態(tài)圖像. 初始時(shí)刻彈尖處節(jié)點(diǎn)與靶板撞擊面中心處節(jié)點(diǎn)接觸(圖13(a))，此時(shí)一個(gè)節(jié)點(diǎn)存在多個(gè)潛在接觸片，需要接觸算法給出合理的接觸方向. 隨著子彈侵入，靶板撞擊面中心附近材料發(fā)生塑性大變形，當(dāng)塑性應(yīng)變達(dá)到失效應(yīng)變時(shí)單元破壞刪除，新生成的單元外表面繼而與子彈接觸. 此后，彈與靶板間的接觸變得更加復(fù)雜：一方面由于部分單元失效破壞導(dǎo)致靶板上的接觸界面存在大量尖角，極不平滑，如圖14所示；另一方面接觸面附近單元變形嚴(yán)重，這些因素導(dǎo)致接觸對(duì)確定、接觸點(diǎn)與接觸方向計(jì)算存在很大難度. 但由圖14可以看出，文中建立的侵蝕接觸算法合理的模擬了該穿甲過(guò)程：侵徹過(guò)程中子彈受力對(duì)稱(chēng)，未發(fā)生偏轉(zhuǎn)；靶板變形與破壞形貌與實(shí)驗(yàn)[12]一致.

圖13 不同時(shí)刻的變形圖像

圖14 靶體接觸界面形狀

圖15給出了侵徹過(guò)程中子彈的速度與加速度時(shí)程曲線(xiàn)，與類(lèi)似子彈侵徹靶板文獻(xiàn)[13]給出的數(shù)據(jù)在規(guī)律上完全一致. 子彈完全穿過(guò)靶板后，剩余速度為709 mm/ms，與實(shí)驗(yàn)數(shù)據(jù)702 mm/ms及Recht-Ipson模型理論理預(yù)測(cè)給出的710 mm/ms很接近.

圖15 彈體的速度與加速度時(shí)程曲線(xiàn)

4 結(jié) 論

針對(duì)沖擊動(dòng)力學(xué)問(wèn)題中存在單元破壞的動(dòng)態(tài)接觸-碰撞，即侵蝕接觸問(wèn)題，文中發(fā)展了一種可處理復(fù)雜接觸界面的局部搜索算法. 新算法通過(guò)引入接觸歷史狀態(tài)量跟蹤從點(diǎn)對(duì)主片的接觸歷史，并結(jié)合接觸強(qiáng)度確定真實(shí)接觸對(duì)，有效克服了侵蝕接觸計(jì)算中常遇到的接觸盲區(qū)問(wèn)題；在接觸點(diǎn)計(jì)算中，發(fā)展了非迭代計(jì)算方法，提高了接觸計(jì)算的穩(wěn)定性.

利用彈性桿碰撞、高曲率凹凸面接觸、彈體穿甲等典型算例驗(yàn)證了文中算法的有效性. 算例分析表明，文中建立的侵蝕接觸算法可有效實(shí)現(xiàn)侵徹接觸問(wèn)題的數(shù)值模擬，可用于復(fù)雜動(dòng)態(tài)接觸問(wèn)題分析并具有較高的計(jì)算精度.