接觸-碰撞算法研究進(jìn)展

, ,

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

1 引 言

變形體間的接觸-碰撞(或稱(chēng)為動(dòng)態(tài)接觸)廣泛存在于實(shí)際工程中[1-4]，如彈體穿甲/侵徹、安全防護(hù)、汽車(chē)防撞以及沖壓成形等。該類(lèi)問(wèn)題通常涉及幾何、材料與邊界的三重非線(xiàn)性，屬于最困難的非線(xiàn)性問(wèn)題之一[3,4]，數(shù)值方法是解決此類(lèi)問(wèn)題的最有效手段。大量分析表明[5]，接觸-碰撞算法是影響數(shù)值計(jì)算精度的重要因素。另一方面，接觸-碰撞過(guò)程中接觸界面與接觸狀態(tài)的確定非常耗時(shí)，接觸計(jì)算通常占到問(wèn)題總求解時(shí)間的一半以上[6,7,8]。因此，發(fā)展高效、高精度的接觸-碰撞算法對(duì)于實(shí)際工程應(yīng)用非常迫切。

接觸-碰撞數(shù)值算法的研究大致起始于20世紀(jì)70年代，迄今已有很多學(xué)者開(kāi)展了相關(guān)研究，發(fā)展了許多優(yōu)秀算法。Bourago等[9]較全面地列出了拉格朗日、非拉格朗日網(wǎng)格及無(wú)網(wǎng)格方法相關(guān)接觸界面算法的635篇論文和專(zhuān)著，鐘陽(yáng)等[10]總結(jié)了顯式有限元計(jì)算常用的接觸界面搜索算法與接觸約束算法。然而，由于接觸-碰撞問(wèn)題的復(fù)雜性，如接觸界面不光滑、接觸對(duì)動(dòng)態(tài)變化等，現(xiàn)有算法在健壯性、精度及效率等方面仍需發(fā)展和完善。目前尚沒(méi)有一種算法可準(zhǔn)確高效地處理所有接觸-碰撞問(wèn)題，接觸-碰撞模擬仍是一個(gè)具有挑戰(zhàn)性的問(wèn)題。

基于拉格朗日框架的顯式有限元是處理變形體間接觸-碰撞最常用的數(shù)值方法，因此本文主要針對(duì)顯式拉格朗日有限元相關(guān)的接觸-碰撞算法進(jìn)行總結(jié)分析。

2 接觸界面的離散

變形體間接觸-碰撞問(wèn)題的有限元理論與基本數(shù)值方法可參考文獻(xiàn)[1-4]。本文僅概要給出該類(lèi)問(wèn)題的力學(xué)描述、有限元求解模型與基本列式。

接觸-碰撞算法可處理任意多個(gè)物體間的相互作用，多物體間的接觸可歸結(jié)為物體間的兩兩作用。不失一般性地，本文以?xún)蓚€(gè)物體間的接觸為例給出其數(shù)學(xué)描述。如圖1所示，相互接觸的兩個(gè)變形體A和B，t時(shí)刻的構(gòu)形分別標(biāo)識(shí)為ΩA和ΩB，其邊界分別為ΓA和ΓB，邊界的外法向分別為nA和nB，接觸邊界為ΓC。

接觸-碰撞系統(tǒng)除應(yīng)滿(mǎn)足通常的連續(xù)體控制方程外以及接觸界面上還應(yīng)滿(mǎn)足一定的運(yùn)動(dòng)學(xué)條件與動(dòng)力學(xué)條件[4]：不可侵入條件和動(dòng)量守恒?？疾煳挥诮佑|界面上分別隸屬于物體A與物體B的任意兩個(gè)點(diǎn)xA和xB，不可侵入條件可表述為

圖1 兩個(gè)物體間的接觸

Fig.1 A two -body contact system

gN=(xB-xA)·nA≥0

(1)

即兩物體間的間隙要求大于或等于0，0值對(duì)應(yīng)于接觸狀態(tài)，gN稱(chēng)為間隙函數(shù)(其負(fù)值稱(chēng)為侵入深度)。界面的動(dòng)量守恒要求作用于兩物體的面力tA與tB滿(mǎn)足

tA+tB=0

(2)

如果不考慮接觸面的焊接或粘接作用，界面法向作用力pN只能為壓力，即

pN≡tA·nA≤0

(3)

與不含接觸的問(wèn)題相比，接觸-碰撞系統(tǒng)僅需在接觸界面上額外引入約束條件(1～3)。因此，可方便地建立接觸-碰撞系統(tǒng)的控制方程弱形式[11](本文暫不考慮界面間的摩擦)：

(4)

(5)

對(duì)應(yīng)于接觸界面約束的罰函數(shù)法弱形式。若pN考慮為獨(dú)立變量，則

(6)

即為拉格朗日乘子法的弱形式。

對(duì)方程(4)進(jìn)行有限元離散，可得半離散方程：

Man=fextn-fintn+fcn

(7)

式中M為質(zhì)量矩陣，an為加速向量，fextn和fintn分別為節(jié)點(diǎn)外力向量與節(jié)點(diǎn)內(nèi)力向量，fcn為節(jié)點(diǎn)接觸力向量，下標(biāo)n標(biāo)識(shí)時(shí)間迭代步。

中心差分法是顯式有限元最常用的積分方法，此時(shí)方程(7)的時(shí)域離散可概括為

an=M-1(fextn-fintn+fcn)

(8)

vn + 1 / 2=vn - 1 / 2+anΔtn

(9)

un +1=un+vn + 1 / 2Δtn + 1 / 2

(10)

式中 Δtn + 1 / 2=tn + 1-tn，Δtn=tn + 1 / 2-tn - 1 / 2；vn + 1 / 2和vn - 1 / 2分別為tn + 1 / 2和tn - 1 / 2時(shí)刻的速度；un和un + 1 / 2為相應(yīng)時(shí)刻的節(jié)點(diǎn)位移向量。

采用顯式積分求解時(shí)，方程(8)右端項(xiàng)中唯一的未知量為tn時(shí)刻的節(jié)點(diǎn)接觸力向量fcn；接觸-碰撞算法研究的目的就是發(fā)展健壯、精確和高效的方法以求得fcn。

接觸-碰撞問(wèn)題中，接觸面通常不能事前確定，并且隨時(shí)間動(dòng)態(tài)變化，因此接觸力計(jì)算可進(jìn)一步分解為接觸搜索與接觸約束施加兩個(gè)子問(wèn)題。接觸搜索的目的是確定系統(tǒng)中哪些部位發(fā)生了接觸或者哪些原已接觸的部位發(fā)生了分離或滑移，接觸約束施加則根據(jù)系統(tǒng)當(dāng)前接觸狀態(tài)計(jì)算出界面接觸力以滿(mǎn)足接觸約束條件。

有限元方法中，相互接觸的變形體(接觸體)離散為有限單元，接觸界面即為變形體外表面的單元面集合或節(jié)點(diǎn)集合。這樣，接觸體間的作用轉(zhuǎn)化為離散的節(jié)點(diǎn)間、節(jié)點(diǎn)與單元面或單元面間的相互作用。

為便于問(wèn)題描述，本文采用Hallquist等[12,13]提出的主從概念。相互接觸的兩個(gè)物體，一個(gè)為主控體，一個(gè)為從屬體。隸屬于主控體的接觸界面稱(chēng)為主面，構(gòu)成主面的單元外表面稱(chēng)為主片，相應(yīng)的節(jié)點(diǎn)稱(chēng)為主點(diǎn)。類(lèi)似地，定義從屬體上的接觸實(shí)體為從面、從片和從點(diǎn)。

根據(jù)接觸界面離散方式的不同，界面模型可分為三類(lèi)。

(1) 點(diǎn)-點(diǎn)NTN(node -to -node)模型[14]。接觸界面上從點(diǎn)與主點(diǎn)一一對(duì)應(yīng)，接觸過(guò)程中從點(diǎn)與主點(diǎn)幾何位置保持一致。這種離散模型的局限是僅適用于小變形且不存在界面相對(duì)滑動(dòng)的接觸分析，如線(xiàn)性準(zhǔn)靜態(tài)問(wèn)題；主面與從面的網(wǎng)格劃分需要保證主點(diǎn)與從點(diǎn)幾何位置一一對(duì)應(yīng)。由于適應(yīng)性差且誤差偏大，即使對(duì)于小變形問(wèn)題分析，目前也已很少采用這種模型。

(2) 點(diǎn)-面NTS(node -to -segment)模型[11-13]。該類(lèi)模型將兩物體間的接觸離散為從點(diǎn)與主片間的相互作用，通過(guò)強(qiáng)制從點(diǎn)不能侵入主片來(lái)施加接觸約束條件。該模型對(duì)接觸界面網(wǎng)格離散幾乎沒(méi)有限制，并且可應(yīng)用于接觸界面任意滑移等復(fù)雜問(wèn)題，是目前最常用的界面離散方式。

(3) 面-面STS(segment-to -segment)模型[15,16]。面-面接觸是近年發(fā)展起來(lái)的一種計(jì)算模型，它通過(guò)在從片與主片間施加約束來(lái)限制物體間的相互侵入。該類(lèi)模型可應(yīng)用于有限變形和界面滑移分析，計(jì)算精度高，但計(jì)算量大。無(wú)論算法，還是具體實(shí)現(xiàn)上，該類(lèi)模型都存在較大困難，目前仍處于發(fā)展階段。

接觸界面離散方式在很大程度上決定了接觸界面搜索方式和接觸力的計(jì)算方法。鑒于NTS模型的廣泛應(yīng)用，本文將針對(duì)基于該類(lèi)模型發(fā)展的相關(guān)算法進(jìn)行討論。

3 接觸搜索算法

如上文所述，接觸界面的確定是接觸-碰撞問(wèn)題數(shù)值計(jì)算中最耗時(shí)的部分。如采用最直接的暴力搜索方法確定接觸對(duì)，其計(jì)算量將為O(N2)，這里N為系統(tǒng)中的節(jié)點(diǎn)數(shù)。因此，暴力搜索對(duì)于實(shí)際工程問(wèn)題(通常需要較多的單元離散，即N很大)幾乎沒(méi)有實(shí)用價(jià)值，必須尋求更高效的接觸對(duì)確定算法。

為加速接觸對(duì)確定，接觸搜索通常分全局搜索和局部搜索兩個(gè)步驟進(jìn)行。全局搜索利用某些簡(jiǎn)單的準(zhǔn)則快速排除掉系統(tǒng)中不可能發(fā)生接觸的部位，或挑選出最可能發(fā)生接觸的從點(diǎn)與主片對(duì)，即確定出接觸測(cè)試對(duì)。局部搜索則從接觸測(cè)試對(duì)中精確確定出真實(shí)發(fā)生接觸的從點(diǎn)與主片對(duì)，即接觸對(duì)。一般而言，全局搜索決定了接觸計(jì)算效率；局部搜索主要影響接觸計(jì)算精度，同時(shí)對(duì)計(jì)算效率也有一定影響。

典型的全局搜索算法有主從面算法[12]、桶排序法[6,8,13,17]、線(xiàn)性位置碼算法[11]、級(jí)域算法[1,18,19]及樹(shù)形算法[20-22]。局部搜索方面，常用算法有點(diǎn)面算法[6,12,13]、內(nèi)外算法[22]、小球算法[23,24]及曲面構(gòu)造法[25,26]。

3.1.1 主從面算法

主從面算法[12]是最早用于求解大規(guī)模接觸-碰撞問(wèn)題的接觸界面搜索方法，現(xiàn)仍然應(yīng)用于通用有限元程序，如LS -DYNA[13]、ABAQUS/Explicit[27]。與暴力搜索方法相比，主從面算法僅在事前指定的有可能發(fā)生接觸的從節(jié)點(diǎn)集與主片集內(nèi)進(jìn)行接觸搜索，即對(duì)每一個(gè)從節(jié)點(diǎn)在接觸主片集內(nèi)確定出與其最近的主節(jié)點(diǎn)，與該主節(jié)點(diǎn)相連的主片視為從點(diǎn)的潛在接觸主片。

主從面算法的計(jì)算復(fù)雜度仍為O(N2)。為減小計(jì)算量，通常引入界面跟蹤技術(shù)[13,28,29]，即假定當(dāng)前時(shí)間步從節(jié)點(diǎn)的最近主節(jié)點(diǎn)為上一時(shí)間步確定的最近主節(jié)點(diǎn)或與該主節(jié)點(diǎn)直接相連的節(jié)點(diǎn)，以減少節(jié)點(diǎn)間距離的計(jì)算次數(shù)。

主從面算法對(duì)于含有少量接觸，且事前可大致確定出接觸部位的問(wèn)題非常有效。但該算法只能處理兩個(gè)表面間的接觸，不能用于單一曲面接觸，也不能用于含材料斷裂破壞的侵蝕接觸分析，甚至不能用于相對(duì)初始構(gòu)形有嚴(yán)重變形的接觸問(wèn)題[6,13]。

3.1.2 桶排序算法

桶排序法是目前應(yīng)用最廣泛的全局搜索算法，該類(lèi)算法的基本思想是[6,13,17],首先，根據(jù)接觸系統(tǒng)單元的某種特征尺寸(如最大單元長(zhǎng)度)將接觸界面空間分割為若干規(guī)則的子區(qū)域(稱(chēng)為桶)；然后，將所有可能參與接觸的節(jié)點(diǎn)根據(jù)其坐標(biāo)排列到這些桶中；接著，根據(jù)節(jié)點(diǎn)所處的桶序號(hào)進(jìn)行排序，形成節(jié)點(diǎn)的有序表；最后，對(duì)每個(gè)從節(jié)點(diǎn)，利用該有序表在其所處桶及相鄰?fù)暗墓?jié)點(diǎn)中確定出距離最近的節(jié)點(diǎn)，該節(jié)點(diǎn)相連接的單元外表面與從節(jié)點(diǎn)構(gòu)成接觸測(cè)試對(duì)。

桶排序法通過(guò)將節(jié)點(diǎn)分組(同一桶中的節(jié)點(diǎn)視為一組)極大縮小了全局接觸搜索范圍，從而大幅提高計(jì)算效率。該算法不需要預(yù)先指定接觸主從面，因此具有良好的適用性，可以用于多體接觸[13]、單曲面接觸[6]和侵蝕接觸[17]等不同接觸類(lèi)型。

對(duì)三維空間問(wèn)題，經(jīng)典桶排序法[6]采用三層嵌套的一維桶排序?qū)崿F(xiàn)。一維桶排序的基本原理如圖2所示，根據(jù)單元特征尺寸將接觸界面空間(xmin,xmax)分割為若干桶，接觸只能發(fā)生在節(jié)點(diǎn)處于同一桶或相鄰?fù)爸械膯卧g。因此，為搜索到所有可能的接觸對(duì)，桶的尺寸應(yīng)不小于l/3，l為系統(tǒng)中最大單元長(zhǎng)度。顯然，這對(duì)于單元網(wǎng)格劃分很不均勻的系統(tǒng)，桶最優(yōu)特征尺寸的確定存在很大難度。桶越小，需要進(jìn)行距離判斷的節(jié)點(diǎn)也就越少，搜索速度越快，但桶過(guò)小可能會(huì)導(dǎo)致接觸搜索遺漏。經(jīng)典桶排序法另外一個(gè)問(wèn)題就是內(nèi)存需求量大[7]。

Heinstein等[7,8]引入包圍盒(the bounding box)改進(jìn)了桶排序算法，利用主片當(dāng)前構(gòu)形與預(yù)測(cè)的下一時(shí)步構(gòu)形定義包圍盒，根據(jù)從節(jié)點(diǎn)運(yùn)動(dòng)速度適當(dāng)擴(kuò)展包圍盒而形成捕捉盒(the capture box)，如圖3 所示；捕捉盒內(nèi)的節(jié)點(diǎn)視為該主片的潛在接觸從點(diǎn)。算法改進(jìn)后，桶的尺寸可取為系統(tǒng)中最小單元的特征長(zhǎng)度，算法效率僅依賴(lài)于從點(diǎn)和主片的數(shù)量，而與桶的尺寸無(wú)關(guān)。另外，該算法通過(guò)定義統(tǒng)一桶編號(hào)而使排序和查找在一維空間進(jìn)行，從而提高了算法效率，降低了內(nèi)存需要。

3.1.3 級(jí)域算法

級(jí)域算法[1,18,19]建立在接觸實(shí)體層級(jí)、層域及接觸域等概念的基礎(chǔ)上。根據(jù)幾何關(guān)系，接觸系統(tǒng)可分為不同層級(jí)的實(shí)體，即接觸體、接觸面、接觸片、接觸邊和接觸點(diǎn)。接觸系統(tǒng)由一個(gè)或多個(gè)接觸體組成，一個(gè)接觸體可包含若干接觸面，一個(gè)接觸面由多個(gè)接觸片構(gòu)成。

層域是某層級(jí)實(shí)體占據(jù)的空間，為避免每個(gè)迭代步都進(jìn)行全局搜索，可將此區(qū)域適當(dāng)擴(kuò)展而形成擴(kuò)展域。當(dāng)從點(diǎn)位于某主點(diǎn)、線(xiàn)或片的接觸域內(nèi)時(shí)，認(rèn)為該從點(diǎn)與相應(yīng)的主點(diǎn)、線(xiàn)或片接觸。點(diǎn)、線(xiàn)與片接觸域的定義如圖4所示。片的接觸域由其邊、外法線(xiàn)方向及某種特征厚度確定，線(xiàn)的接觸域由其外法線(xiàn)方向及其所在的片厚度確定，點(diǎn)的接觸域則由其所在線(xiàn)的外法線(xiàn)方向及相應(yīng)片厚度確定。

接觸搜索時(shí)，先在較高級(jí)的層域間進(jìn)行，若兩個(gè)域存在重疊，則繼續(xù)進(jìn)行下一級(jí)層域的相交檢查，直至確定出接觸測(cè)試對(duì)；若兩域不存在重疊，則兩接觸實(shí)體間不存在接觸關(guān)系，接觸搜索停止。

級(jí)域算法利用了接觸體的層級(jí)結(jié)構(gòu)，在較高級(jí)實(shí)體中排除了不必要的接觸檢查，理論上可大幅提高接觸搜索效率，這對(duì)于同級(jí)域重疊部分較小的接觸問(wèn)題非常有效。但對(duì)于同級(jí)實(shí)體存在較多重疊的接觸問(wèn)題，如單曲面接觸，計(jì)算效率不如位置碼算法[30]。另外，程序?qū)崿F(xiàn)方面，由于該類(lèi)算法涉及不同層級(jí)接觸實(shí)體的定義及相應(yīng)的數(shù)據(jù)管理，程序設(shè)計(jì)與實(shí)現(xiàn)都很繁瑣。

圖2 一維桶排序示意圖

Fig.2 Illustration of bucket sorting in 1D

圖3 主片的包圍盒與捕捉盒[8]

Fig.3 Bounding box and capture box for a master segment[8]

3.1.4 位置碼算法

Oldenburg等[11]在桶排序和級(jí)域算法基礎(chǔ)上建立了位置碼算法。該算法對(duì)每一個(gè)主片首先確定位于其片域內(nèi)的接觸從點(diǎn)；對(duì)位于主片域內(nèi)的從點(diǎn)進(jìn)一步檢查是否處于主片接觸域內(nèi)，接觸域內(nèi)的從點(diǎn)與主片構(gòu)成接觸測(cè)試對(duì)。接觸域和片域的二維圖如圖5所示。

定義了主片域后，位于主片域內(nèi)從點(diǎn)的三維空間搜索利用位置碼轉(zhuǎn)換到一維空間中進(jìn)行：接觸界面空間分割成若干桶；根據(jù)節(jié)點(diǎn)所在的桶，對(duì)每一個(gè)節(jié)點(diǎn)賦予一個(gè)位置碼；根據(jù)位置碼與節(jié)點(diǎn)坐標(biāo)將三維空間中的節(jié)點(diǎn)排序到一維數(shù)組(稱(chēng)為位置碼向量)；對(duì)每一個(gè)主片，確定出主片域占據(jù)的桶，利用二分查找算法在位置碼向量中找出這些桶中的節(jié)點(diǎn)(潛在接觸從點(diǎn))。同時(shí)，位置碼算法采用層級(jí)概念定義了三層接觸實(shí)體，即接觸面、接觸片和接觸點(diǎn)，每一個(gè)接觸面定義一個(gè)位置碼向量。

位置碼算法將三維空間的節(jié)點(diǎn)排序和搜索映射到一維數(shù)組中進(jìn)行，并利用二分查找法加速接觸從點(diǎn)的搜索速度，具有很高的計(jì)算效率(計(jì)算復(fù)雜度為O(Nlog2N))。并且，該算法無(wú)需區(qū)分主從面接觸和單曲面接觸，通用性強(qiáng)。

位置碼算法的效率有一定的網(wǎng)格方向依賴(lài)性，這是由桶的編碼方式?jīng)Q定的。Diekmann等[30]以空間填充曲線(xiàn)法SPC(Space Filling Curve algorithm)代替逐行掃描的桶編碼方式，改進(jìn)了位置碼算法。在處理細(xì)長(zhǎng)結(jié)構(gòu)或復(fù)雜結(jié)構(gòu)間的接觸時(shí)，SPC算法具有更高的計(jì)算效率。

3.1.5 樹(shù)形算法

樹(shù)形算法是新近發(fā)展的一類(lèi)全局搜索方法。

圖4 點(diǎn)、線(xiàn)和片的接觸域

Fig. 4 Contact territories of a node,an edge and a segment

圖5 接觸域和片域的定義

Fig.5 Conception of contact territory and segment territory

該方法的基本思想是，采用樹(shù)形數(shù)據(jù)結(jié)構(gòu)，如二叉樹(shù)和四叉樹(shù)，存儲(chǔ)接觸片幾何信息，利用樹(shù)的快速遍歷方法確定出有可能參與接觸的從點(diǎn)與主片或者從片與主片對(duì)。

針對(duì)沖壓成形中板料與模具間的動(dòng)態(tài)接觸，Bruneel等[20]以四叉樹(shù)存儲(chǔ)模具(模具表面離散為多邊形幾何平面)的M個(gè)接觸片，建立了計(jì)算復(fù)雜度為O(Nlog4M)的QuickTrace算法。Yang等[21]針對(duì)mortar接觸中的面-面接觸對(duì)確定，基于層次包圍盒BVH(Bounding Volume Hierarchy)發(fā)展了適用于大變形滑移接觸問(wèn)題的二叉樹(shù)接觸搜索算法；柳陽(yáng)等[31]采用類(lèi)似方法實(shí)現(xiàn)顯式有限元中動(dòng)態(tài)接觸界面的快速搜索。

與經(jīng)典桶排序算法[6,13]相比，樹(shù)形算法的計(jì)算量?jī)H與接觸測(cè)試對(duì)的數(shù)量相關(guān)，而與接觸系統(tǒng)中潛在接觸對(duì)的空間分布無(wú)關(guān)[21]。另外，樹(shù)形算法中子節(jié)點(diǎn)與父節(jié)點(diǎn)有層級(jí)關(guān)系，這類(lèi)似于級(jí)域算法[18]中級(jí)的屬性。因此，對(duì)于大規(guī)模復(fù)雜接觸問(wèn)題，樹(shù)形算法在內(nèi)存需求與計(jì)算效率方面具有較大優(yōu)勢(shì)。在程序?qū)崿F(xiàn)方面，樹(shù)形算法可采用遞歸方式設(shè)計(jì)，程序設(shè)計(jì)與代碼實(shí)現(xiàn)都相對(duì)簡(jiǎn)單。

3.1.6 其他方法

Papadopoulos等[32]基于類(lèi)似于小球算法的思想建立了所謂的球形排序算法(spherical-sorting algorithm)。該算法對(duì)于小規(guī)?；蛘邇H含有少量接觸對(duì)的問(wèn)題計(jì)算效率很高，但不適用于復(fù)雜工程問(wèn)題。因?yàn)樵撍惴ㄔ诿總€(gè)求解步都需要對(duì)所有的潛在接觸片進(jìn)行距離比較和判斷，當(dāng)接觸片較多時(shí)計(jì)算效率會(huì)遠(yuǎn)低于桶排序和位置碼等接觸搜索算法。

Mahadevaiah[33]耦合球形排序算法與桶排序方法建立了球形桶排序算法(sphere-bucket-sort algorithm)。該算法對(duì)每個(gè)可能參與接觸的節(jié)點(diǎn)建立一個(gè)虛擬的包圍球，將包圍球排序到桶中，然后利用經(jīng)典桶排序方法[6]確定出可能的接觸對(duì)。該算法集成了桶排序算法與小球算法的優(yōu)勢(shì)，理論上具有較高的計(jì)算效率。

近年，針對(duì)離散粒子接觸系統(tǒng)，發(fā)展了大量?jī)?yōu)秀的接觸搜索算法，雖然不能直接用于有限元接觸計(jì)算，但為有限元接觸搜索算法的設(shè)計(jì)提供了重要借鑒。王福軍等[25,34]借鑒NBS算法[35]，利用鏈表技術(shù)對(duì)桶排序算法進(jìn)一步改進(jìn)，將節(jié)點(diǎn)位置碼排序到鏈表結(jié)構(gòu)中，接觸搜索轉(zhuǎn)變?yōu)閺膬蓚€(gè)一維數(shù)組中取數(shù)據(jù)的過(guò)程，計(jì)算量降低為O(N)。Chen等[36]借用CGRID算法[37]的不同尺寸粒子遷移策略，引入了層、列和方格概念，建立了計(jì)算效率不受網(wǎng)格尺寸影響以及內(nèi)存需求與計(jì)算成本呈線(xiàn)性復(fù)雜度的全局搜索算法。

與全局搜索相比，局部搜索算法的研究相對(duì)較少，目前主要采用的方法有點(diǎn)面算法和小球類(lèi)算法兩大類(lèi)。另有少量搜索算法主要是針對(duì)兩類(lèi)算法的不足而提出的改進(jìn)方案，如內(nèi)外算法、曲面構(gòu)造法及靜動(dòng)態(tài)檢測(cè)法等。

3.2.1 點(diǎn)面算法

點(diǎn)面算法[6](node -to-segment algorithm)是目前最常用的一類(lèi)局部搜索算法，廣泛應(yīng)用于通用有限元程序[13,27]。該算法的計(jì)算過(guò)程為[6]，對(duì)每個(gè)從點(diǎn)，在與其最近主點(diǎn)相連的主片中確定出可能接觸的主片；計(jì)算從點(diǎn)在該主片上的投影點(diǎn)(或稱(chēng)為接觸點(diǎn))；利用投影點(diǎn)確定從點(diǎn)對(duì)主片的侵入深度。

假定從點(diǎn)ns最近的主點(diǎn)為nm，與nm相連的主片為s1，s2，s3和s4，如圖6所示。如果nm與ns不重合，當(dāng)滿(mǎn)足式(11)時(shí)，ns可能與主片si接觸。

(ci×s)·(ci×ci +1)>0

(11a)

(ci×s)·(s×ci +1)>0

(11b)

式中ci和ci +1為過(guò)點(diǎn)nm的主片兩邊定義的向量,s為向量g在主片si上的投影，而g為連接nm和ns的向量，

s=g-(g·m)m

(12)

(13)

(14)

式中ni為投影方向(或侵徹方向)，

(15)

若l<0，表明從點(diǎn)對(duì)主片有侵入，則根據(jù)接觸約束算法在從點(diǎn)與主片間施加接觸約束力。

當(dāng)接觸界面光滑且網(wǎng)格質(zhì)量較好時(shí)，點(diǎn)面搜索算法非常有效，這也是該方法受到廣泛采用的原因。但算法存在如下不足。

(1) 經(jīng)典點(diǎn)面算法[6,13]采用雙線(xiàn)性參數(shù)化方程表征主片，利用Newton-Raphson迭代求解偏微分方程組(13)的方式計(jì)算接觸點(diǎn)。該方法在大多數(shù)情況下是有效的，但當(dāng)主片發(fā)生較大的翹曲變形時(shí)可能存在收斂困難和不穩(wěn)定等問(wèn)題[22]。

(2) 以主片外法線(xiàn)方向作為從點(diǎn)的投影方向，存在投影點(diǎn)不唯一或無(wú)投影點(diǎn)等奇異性問(wèn)題[38]，即所謂的盲區(qū)(或稱(chēng)為死區(qū))，如圖8所示。

(3) 接觸面網(wǎng)格質(zhì)量較差時(shí)，如細(xì)長(zhǎng)比過(guò)大，可能會(huì)發(fā)生接觸搜索遺漏，如圖9所示。

(4) 由于該方法侵入檢查是針對(duì)從點(diǎn)與主片的關(guān)系進(jìn)行的，因此不能處理如圖10所示的邊-面和面-面接觸情況。

圖6 從點(diǎn)在主片上的投影

Fig.6 Projection of the slave node onto master segment

圖7 接觸點(diǎn)的確定

Fig.7 Determination of contact point

圖8 兩種接觸搜索盲區(qū)

Fig.8 Examples of blind zones

圖9 網(wǎng)格細(xì)長(zhǎng)比過(guò)大造成接觸搜索遺漏[13]

Fig.9 Failure to find contact pairs due to poor aspect ratio in mesh[13]

3.2.2 內(nèi)外算法

面向金屬成形模擬，針對(duì)經(jīng)典點(diǎn)面算法[6]接觸盲區(qū)與計(jì)算效率的不足，Wang等[22]改進(jìn)了從點(diǎn)對(duì)主片的投影方向以及接觸點(diǎn)的計(jì)算方法，從而建立了內(nèi)外算法(inside -outside -algorithm)。

內(nèi)外算法中，以從點(diǎn)的外法線(xiàn)方向作為從點(diǎn)對(duì)主片的投影方向(同時(shí)也是接觸力施加方向)。從點(diǎn)外法向定義為與該節(jié)點(diǎn)相連的接觸片法線(xiàn)方向的平均值，如圖11所示。

(16)

從點(diǎn)相對(duì)于主片一條邊的內(nèi)或外狀態(tài)，通過(guò)從點(diǎn)與該邊構(gòu)成的三角形在主片上的投影來(lái)確定，如圖12所示，

Δa=n·(ra b×ra d)

(17)

若Δa≤0，則從點(diǎn)處于邊內(nèi)部。如果從點(diǎn)處于主片所有邊的內(nèi)部或者外部，則從點(diǎn)位于該主片內(nèi)部，否則位于主片外部。

圖10 邊-面接觸和面-面接觸

Fig.10 Edge -to -surface contact and surface -to -surface contact

圖11 從點(diǎn)S的外法線(xiàn)方向

Fig.11 Definition of normal vector of nodeS

圖12 從點(diǎn)的內(nèi)外狀態(tài)檢查

Fig.12 Inside/outside checking on a 4-node segment

如果從點(diǎn)位于主片內(nèi)部，則利用式(18)計(jì)算從點(diǎn)對(duì)主片的侵入量，

g=n·(x-xs)

(18)

式中x為投影點(diǎn)，利用式(17)可以解析給出；xs為從節(jié)點(diǎn)位置向量。

內(nèi)外算法的優(yōu)勢(shì)在于，接觸點(diǎn)可以唯一確定，避免了從點(diǎn)對(duì)主片的多重接觸，并且接觸點(diǎn)計(jì)算不需要迭代。因此，該算法具有很好的健壯性與較高的計(jì)算效率。由于內(nèi)外算法假定接觸點(diǎn)處主從面的法向共軸(但方向相反)，當(dāng)從點(diǎn)位于兩個(gè)低階單元的交界線(xiàn)上時(shí)，尤其是當(dāng)網(wǎng)格不夠精細(xì)時(shí)，內(nèi)外算法給出的接觸方向可能會(huì)與實(shí)際存在較大偏差[26]。

3.2.3 構(gòu)造曲面法

點(diǎn)面算法存在投影奇異性的根本原因是低階有限元近似使得接觸界面僅C0連續(xù)，接觸片交界處接觸界面外法線(xiàn)方向存在間斷。因此，接觸片光滑化是解決盲區(qū)問(wèn)題的一種有效方式。構(gòu)造曲面法(或稱(chēng)為曲面光滑化法)正是基于這種思想建立的，即基于離散的接觸界面拓?fù)湫畔?，采用高階插值函數(shù)重構(gòu)接觸面，使新接觸面外法線(xiàn)方向空間連續(xù)分布。

構(gòu)造曲面法不僅可以克服接觸盲區(qū)，也有利于提高接觸力計(jì)算精度，消除接觸滑移時(shí)接觸力的非物理震蕩。關(guān)于構(gòu)造曲面法，已有很多學(xué)者開(kāi)展相關(guān)研究，提出了一些光滑化算法，如三次樣條法[39]、NURBS法[40]、Gregory片法[41]和FFS算法[25,26](Free -Formed-Surface algorithm)等。其中，F(xiàn)FS法建立的曲面有C1連續(xù)性，可用于三維問(wèn)題，且無(wú)需引入額外邊界定義或額外自由度。

FFS方法利用接觸片節(jié)點(diǎn)坐標(biāo)構(gòu)造參數(shù)化光滑曲面片。如圖13所示，任意一個(gè)由四條邊構(gòu)成的接觸片，參數(shù)化光滑曲面可表示為

(19)

式中x為曲面的整體坐標(biāo)，u和w為曲面的局部坐標(biāo)，ci j為參數(shù)向量，可根據(jù)相關(guān)節(jié)點(diǎn)坐標(biāo)確定。兩節(jié)點(diǎn)間的邊界線(xiàn)利用相鄰節(jié)點(diǎn)坐標(biāo)以分段三次Hermite曲線(xiàn)構(gòu)造，從而保證了相鄰接觸片間的連續(xù)。

FFS算法構(gòu)造的曲面具有C1連續(xù)性，接觸點(diǎn)與接觸方向計(jì)算精度高，很大程度上克服了接觸盲區(qū)問(wèn)題。算法的不足是計(jì)算量增加較多，尤其是接觸片畸變較大時(shí),需要將接觸片進(jìn)行多次子片剖分才能得到滿(mǎn)足精度的結(jié)果。

3.2.4 靜動(dòng)態(tài)檢測(cè)法

經(jīng)典點(diǎn)面局部搜索算法[6]僅利用主片和從點(diǎn)的當(dāng)前位置信息確定接觸狀態(tài)，對(duì)于發(fā)生穿透的從點(diǎn)，沿主片外法線(xiàn)方向施加接觸力將其外推至主片上。這可能會(huì)導(dǎo)致接觸主片判斷錯(cuò)誤或從點(diǎn)回退方向(也就是接觸力方向)不正確等問(wèn)題，為此Heinstein等[7,8]提出了靜動(dòng)態(tài)檢測(cè)算法。該算法對(duì)于上一時(shí)間步已經(jīng)處于接觸狀態(tài)的從點(diǎn)采用靜態(tài)檢測(cè)法確定當(dāng)前接觸狀態(tài)；對(duì)上一時(shí)間步未發(fā)生接觸，但在當(dāng)前時(shí)間步發(fā)生接觸的從點(diǎn)采用動(dòng)態(tài)接觸檢測(cè)。

(1) 動(dòng)態(tài)檢測(cè)。根據(jù)從點(diǎn)和主片的當(dāng)前位置及下一時(shí)間步的預(yù)測(cè)速度，利用從點(diǎn)和主片接觸時(shí)刻的共面關(guān)系，確定出接觸時(shí)刻與接觸點(diǎn)。三維情況下，這種接觸檢查簡(jiǎn)化為判斷一個(gè)移動(dòng)三角形和一個(gè)移動(dòng)點(diǎn)之間的接觸問(wèn)題(對(duì)于四邊形主片，可以利用主片中心點(diǎn)與四個(gè)邊將主片分解為四個(gè)三角形)，如圖14所示。求得接觸時(shí)刻與接觸點(diǎn)后，進(jìn)一步判斷接觸時(shí)刻是否滿(mǎn)足步長(zhǎng)穩(wěn)定性要求，及接觸點(diǎn)是否位于主片內(nèi)部，如滿(mǎn)足則發(fā)生接觸。

3.2.5 小球算法與分裂小球算法

小球算法(pinball algorithm)是Belytschko等[23]為適應(yīng)高速撞擊和侵徹接觸計(jì)算向量化需求

圖13 自由曲面的構(gòu)造

Fig.13 Construction of FFS

圖14 基于速度的動(dòng)態(tài)接觸檢查

Fig.14 Contact checking based on velocity

提出的一種高效接觸算法。該算法的核心思想是以單元內(nèi)嵌小球(稱(chēng)為等效小球)代替可能參與接觸的表面單元，單元間的侵入檢查與接觸約束轉(zhuǎn)化為小球之間的距離判斷與約束施加。等效小球的體積為相應(yīng)單元的體積，并假定變形過(guò)程中保持不變，球心為單元各節(jié)點(diǎn)坐標(biāo)的平均值。小球概念的二維示意如圖17所示。

引入小球概念后，單元間的侵入檢查變得極為簡(jiǎn)單。假定主面單元與從面單元等效小球的半徑分別為R1和R2，球心位置矢量為C1和C2。當(dāng)兩個(gè)小球球心距離d滿(mǎn)足式(20)時(shí)，認(rèn)為兩單元發(fā)生接觸。

d

(20)

侵入深度g由式(21)計(jì)算得到，

dTd=(R1+R2)2

(21)

式中d=C1-C2+gn，n=(n2-n1)/‖n2-n1‖，n1和n2分別為主片與從片的外法線(xiàn)方向。

小球算法可用于實(shí)體單元或薄殼單元的任意類(lèi)型接觸分析。但用于薄殼單元時(shí)，若殼元厚度相對(duì)于其特征尺寸很小，則小球算法精度會(huì)非常差；若殼單元間存在初始接觸，小球算法可能會(huì)失效。為此，Belytschko等[24]在小球算法基礎(chǔ)上發(fā)展了分裂小球算法(splitting pinball method)。

與小球算法類(lèi)似，分裂小球算法以小球代替每一個(gè)參與接觸的表面單元，但小球的半徑要足夠大以完全覆蓋相應(yīng)單元，該小球稱(chēng)為父球。父球間若有重疊區(qū)則表示兩單元間可能存在接觸，此時(shí)父球分裂為下一級(jí)的子球進(jìn)一步作重疊區(qū)檢查，直至最后一級(jí)子球的半徑與殼厚度相當(dāng)。四邊形殼元的小球?qū)蛹?jí)如圖18所示，接觸檢查與小球逐級(jí)分裂的一維示意如圖19所示。

圖15 從點(diǎn)對(duì)凹面侵徹

Fig.15 Contact checking for a concave surface

圖16 從點(diǎn)對(duì)凸面的侵徹

Fig.16 Contact checking for a convex surface

小球算法和分裂小球算法以統(tǒng)一的方式自動(dòng)處理邊-邊、邊-面及面-體等接觸類(lèi)型，算法便于向量化和并行計(jì)算；單元間的接觸檢查與侵入量計(jì)算簡(jiǎn)單高效。小球算法和分裂小球算法適用于碰撞和侵徹等問(wèn)題的模擬。由于算法未考慮接觸界面間的摩擦效應(yīng)，不能應(yīng)用于滑移和摩擦為關(guān)鍵因素的接觸分析中。另外，當(dāng)接觸界面單元特征尺寸差別較大時(shí)，接觸力的計(jì)算精度較差，不準(zhǔn)確的幾何近似和小球體積不變假設(shè)限制了該算法的應(yīng)用[7]。

4 接觸約束算法

接觸約束算法，也稱(chēng)為接觸力算法，基本方法有拉氏乘子法和罰函數(shù)法[1,4]。基于這兩類(lèi)方法，還發(fā)展了增廣拉氏乘子法和攝動(dòng)拉氏乘子法。拉氏乘子類(lèi)方法計(jì)算格式與顯式有限元計(jì)算不相容，不能直接用于顯式有限元的接觸分析，罰函數(shù)法是目前顯式有限元程序的主要方法。

罰函數(shù)法的基本原理是，如果從點(diǎn)對(duì)主片沒(méi)有侵入則不做處理；如果從點(diǎn)侵入主片，則在從點(diǎn)與主片間引入法向接觸力，將從點(diǎn)推回到變形后的接觸主片上，以滿(mǎn)足不可侵入條件。

圖17 平面單元的內(nèi)嵌小球

Fig.17 Embedded pinball in two dimensions

圖18 四邊形薄殼單元的小球?qū)蛹?jí)

Fig.18 Pinball hierarchy of 4-node shell element

圖19 一維單元的接觸檢查

Fig.19 Contact detection using pinballs

罰函數(shù)法不增加系統(tǒng)方程的自由度數(shù)，接觸力計(jì)算無(wú)需聯(lián)立求解方程組，程序?qū)崿F(xiàn)簡(jiǎn)單且計(jì)算效率高，廣泛應(yīng)用于顯式有限元程序中，如 LS -Dyna[13]，Abaqus/Explicit[27]和Pronto3D[42]等。罰函數(shù)法最大的不足在于不可侵入條件只能近似滿(mǎn)足，計(jì)算精度嚴(yán)重依賴(lài)于罰因子的選取，而罰因子取值往往依賴(lài)于使用者的經(jīng)驗(yàn)。

通用程序一般會(huì)對(duì)罰因子進(jìn)行規(guī)則化處理，以便用戶(hù)選取相對(duì)合理的罰參數(shù)。以L(fǎng)S -Dyna[13]為例，作用于從點(diǎn)的法向接觸力為

fs=-lkini

(22)

式中l(wèi)為從點(diǎn)對(duì)主片的侵入深度

(23)

式中ni,t和r等變量的定義見(jiàn)3.2.1節(jié)。其中ki即為罰因子。LS -Dyna利用主片面積Ai，主片所在單元體積Vi及材料體積模量Ki對(duì)其規(guī)則化

(24)

規(guī)則化后的接觸罰參數(shù)fSI取值范圍通常為 (0.0,1.0]。得到從點(diǎn)的接觸力后，根據(jù)牛頓第三定律，作用在主片接觸點(diǎn)上的力為-fs，按照式(25)可將該合力分配到主片各節(jié)點(diǎn)上，

(25)

拉氏乘子法中，以接觸對(duì)間的接觸力作為乘子來(lái)限制接觸體間的相互侵入，不可侵入條件可精確滿(mǎn)足。但該方法引入了新的未知變量，增加了方程組的自由度數(shù)；接觸力計(jì)算需要聯(lián)立方程求解，這與顯式有限元計(jì)算不相容。為了在顯式有限元中應(yīng)用乘子法進(jìn)行接觸力計(jì)算，必須對(duì)經(jīng)典的拉氏乘子法做必要的格式修正。

Carpenter等[43]利用上一時(shí)間步的乘子與本步的位移建立接觸約束方程，采用修正高斯-塞得爾迭代法求解拉氏乘子，建立了二維接觸問(wèn)題的前增量拉氏乘子法(forward increment Lagrange multiplier)。Sha等[44]在前增量位移中心差分方法(forward incremental displacement-central difference method)基礎(chǔ)上，建立了接觸約束的線(xiàn)性補(bǔ)充方程，并給出了改進(jìn)的共軛梯度求解算法。Zhong[1]提出了防御節(jié)點(diǎn)法(defense node algorithm)；王福軍等[45]考慮相鄰接觸點(diǎn)間的相互影響，在防御節(jié)點(diǎn)法基礎(chǔ)上建立了局部拉氏乘子方法；Chen等[46]對(duì)初始碰撞和后續(xù)接觸采用不同的接觸約束方程，并利用高斯-塞得爾迭代法求解拉氏乘子。這些方法中，防御節(jié)點(diǎn)法不需要任何迭代，計(jì)算效率高，是乘子類(lèi)方法在顯式有限元中成功應(yīng)用的典范。

防御節(jié)點(diǎn)法的基本思想是，將從點(diǎn)與主片之間的接觸作用轉(zhuǎn)化為從點(diǎn)與防御節(jié)點(diǎn)之間的接觸，防御節(jié)點(diǎn)攜帶了主片的運(yùn)動(dòng)學(xué)與動(dòng)力學(xué)信息，如位移、速度及節(jié)點(diǎn)力等。防御節(jié)點(diǎn)位于接觸點(diǎn)，表征主片的運(yùn)動(dòng)，因此其速度和加速度可以利用主片的節(jié)點(diǎn)速度和加速度描述，F(xiàn)為除接觸力以外的力，f為從點(diǎn)與防御節(jié)點(diǎn)間的接觸力。計(jì)算f時(shí)，施加不可侵入條件。得到f后，接觸力就可以分布到主點(diǎn)上。

(26，27)

Ma=F+f

(28)

5 接觸-碰撞問(wèn)題的并行計(jì)算

接觸-碰撞問(wèn)題通常采用顯式有限元進(jìn)行計(jì)算。顯式雖然避免了內(nèi)存的龐大需求，但積分穩(wěn)定性要求時(shí)間步長(zhǎng)必須足夠小，這導(dǎo)致接觸-碰撞問(wèn)題求解通常需要大量迭代，對(duì)于大規(guī)?；驈?fù)雜幾何構(gòu)形問(wèn)題非常耗時(shí)，因此并行計(jì)算是接觸-碰撞數(shù)值分析的必然選擇。

區(qū)域分解(spatial domain decomposition method)是實(shí)現(xiàn)接觸-碰撞問(wèn)題大規(guī)模并行計(jì)算的最有效方式。區(qū)域分解的基本思想[47]是，將一個(gè)復(fù)雜系統(tǒng)按照一定的準(zhǔn)則(如物理特性或幾何形狀)劃分為一系列子系統(tǒng)(稱(chēng)為子區(qū)域)；將子區(qū)域分配到各個(gè)處理器上分別處理；子區(qū)域間通過(guò)消息傳遞協(xié)同計(jì)算，最終求得原問(wèn)題的解。區(qū)域分解可以為大規(guī)模問(wèn)題提供高度并行且可擴(kuò)展的健壯算法。

與不含接觸的動(dòng)力學(xué)問(wèn)題或準(zhǔn)靜態(tài)接觸問(wèn)題不同，接觸-碰撞問(wèn)題的接觸域通常具有全局性與動(dòng)態(tài)性，計(jì)算節(jié)點(diǎn)之間的通信關(guān)系在系統(tǒng)求解之前不能確定。如圖20所示的一個(gè)簡(jiǎn)單接觸問(wèn)題，開(kāi)始時(shí)刻下方物體上表面的左前部位與上面物體發(fā)生接觸，而在某個(gè)時(shí)刻下方物體上表面的右后部位與上面物體發(fā)生接觸。接觸界面的全局性和動(dòng)態(tài)性對(duì)接觸-碰撞問(wèn)題的區(qū)域分割和計(jì)算節(jié)點(diǎn)間的通信關(guān)系設(shè)計(jì)帶來(lái)極大挑戰(zhàn)。

(1) 對(duì)有限元計(jì)算(主要是內(nèi)力計(jì)算)高效的某種區(qū)域分解在接觸計(jì)算中往往不能取得好的并行性能，通常需要針對(duì)接觸而引入另外的區(qū)域分割。

(2) 動(dòng)態(tài)接觸的全局性導(dǎo)致計(jì)算節(jié)點(diǎn)間存在大量的數(shù)據(jù)交互，最大限度地減小通信量才能保證算法的擴(kuò)展性與并行效率。

(3) 接觸對(duì)動(dòng)態(tài)變化，局部接觸搜索和約束施加在不同時(shí)刻可能由不同的計(jì)算節(jié)點(diǎn)完成，極易造成各節(jié)點(diǎn)機(jī)負(fù)載的不均衡。

基于區(qū)域分解的并行接觸算法研究始于20世紀(jì)90年代。Malone等[48,49]基于體積檢測(cè)，通過(guò)引入接觸域(contact domain)實(shí)現(xiàn)接觸搜索的本地計(jì)算，設(shè)計(jì)了可應(yīng)用于三維殼體結(jié)構(gòu)計(jì)算的并行接觸-碰撞算法。Waston等[50]針對(duì)簡(jiǎn)單接觸問(wèn)題，設(shè)計(jì)了基于單元分解的并行算法；Zhong等[51]針對(duì)Connection Machine機(jī)型，設(shè)計(jì)了基于桶排序的并行接觸算法；Har等[52]針對(duì)薄殼單元結(jié)構(gòu)，設(shè)計(jì)了三維桶排序并行算法，利用IBM SP2的10個(gè)處理器進(jìn)行了接觸并行計(jì)算；陳成軍等[53,54]開(kāi)展了類(lèi)似工作；文獻(xiàn)[55-58]提出了雙層次區(qū)域分解方案(dual-level domain decomposition scheme)，利用8個(gè)處理器實(shí)現(xiàn)了32萬(wàn)自由度的沖擊接觸問(wèn)題并行計(jì)算；金先龍等[59,60]基于坐標(biāo)遞歸二分法RCB(Recursive Coordinate Bisection algorithm)設(shè)計(jì)了接觸均衡分區(qū)算法CBB(Contact Balance Bisection algorithm)，并在曙光4000A上利用16個(gè)CPU進(jìn)行了汽車(chē)碰撞安全性數(shù)值模擬以及輕軌列車(chē)-斜拉橋耦合振動(dòng)響應(yīng)分析，但該算法僅適用于接觸界面事前可以確定的問(wèn)題；Paik等[61]利用不規(guī)則通信重構(gòu)從節(jié)點(diǎn)集，發(fā)展了并行接觸算法，利用256個(gè)處理器進(jìn)行了球撞擊平板的模擬。

上述算法僅適用于少量CPU核的集群環(huán)境，有良好擴(kuò)展性、適用于大型并行系統(tǒng)(數(shù)千CPU核)的并行接觸算法研究還較少。目前該方面的研究主要限于美國(guó)Sandia國(guó)家實(shí)驗(yàn)室[62,63]、Livermore實(shí)驗(yàn)室[64-66]及中國(guó)工程物理研究院[67,68]等少量團(tuán)隊(duì)。Sandia實(shí)驗(yàn)室的Attaway等[62，63]于1998年使用PRONTO3D程序在Sandia/Intel TFLOPS并行機(jī)上實(shí)現(xiàn)了薄壁柱殼動(dòng)態(tài)屈曲至完全壓潰模擬的3600個(gè)計(jì)算節(jié)點(diǎn)并行，2000年利用3504個(gè)計(jì)算節(jié)點(diǎn)模擬了包含六百多萬(wàn)個(gè)單元的包裝箱碰撞和壓潰響應(yīng)。最近，中國(guó)工程物理研究院總體工程研究所陳成軍等[68]利用8192 CPU核實(shí)現(xiàn)了16.1億自由度規(guī)模的并行接觸計(jì)算。需要特別指出的是，近年一些商業(yè)程序，如LS -DYNA等，也開(kāi)始致力于并行版本的開(kāi)發(fā)[69-71]。

圖20 動(dòng)態(tài)接觸示意圖

Fig.20 Illustration of dynamic contacting

基于區(qū)域分解實(shí)現(xiàn)接觸-碰撞問(wèn)題的粗粒度并行有單一區(qū)域分解法和雙重區(qū)域分解法兩種并行方案。雙重區(qū)域分解法中，有限元計(jì)算采用一種子域剖分(稱(chēng)為主分區(qū))，接觸計(jì)算采用另外一套分區(qū)(稱(chēng)為次分區(qū))；而單一區(qū)域分解法，接觸區(qū)域不作單獨(dú)的子區(qū)分割，接觸的并行計(jì)算基于主分區(qū)進(jìn)行。通常，雙重區(qū)域分解法可獲得更好的負(fù)載均衡，但節(jié)點(diǎn)機(jī)間的通信復(fù)雜；單一分區(qū)法通信拓?fù)浜?jiǎn)單，但易導(dǎo)致負(fù)載不均衡。

單一區(qū)域分解方法中，接觸計(jì)算的并行基于主分區(qū)進(jìn)行。接觸算法并行的關(guān)鍵在于，每個(gè)計(jì)算節(jié)點(diǎn)可有效得到本地計(jì)算所需的全部信息。因此，該類(lèi)算法通常引入影像區(qū)(ghost)或接觸域(contact domain)存儲(chǔ)本地接觸計(jì)算需要但駐留在其他節(jié)點(diǎn)機(jī)上的信息。

基于單一區(qū)域分解的并行接觸算法設(shè)計(jì)相對(duì)簡(jiǎn)單，本文以Malone等[49]的工作為例簡(jiǎn)要描述該類(lèi)算法的實(shí)現(xiàn)方式。

(1) 不考慮接觸界面，對(duì)求解域作區(qū)域剖分，分配給各計(jì)算節(jié)點(diǎn)。所有計(jì)算節(jié)點(diǎn)作如下相同的計(jì)算。

(2) 對(duì)任意計(jì)算節(jié)點(diǎn)PA，構(gòu)建覆蓋該計(jì)算節(jié)點(diǎn)單元幾何空間的包圍盒(bounding box)，并傳遞到其他計(jì)算節(jié)點(diǎn)，同時(shí)接收其他計(jì)算節(jié)點(diǎn)的包圍盒信息。

(3) 利用體積檢查(volume checking)判斷PA包圍盒與其他包圍盒是否有重疊。若無(wú)重疊，則兩個(gè)計(jì)算節(jié)點(diǎn)的單元間無(wú)潛在接觸；若有重疊部分，則將與PA包圍盒有重疊的計(jì)算節(jié)點(diǎn)上的單元信息收集到PA的接觸域，同時(shí)將PA單元傳遞給對(duì)方計(jì)算節(jié)點(diǎn)的接觸域。

(4) 對(duì)每個(gè)計(jì)算節(jié)點(diǎn)利用串行接觸算法進(jìn)行接觸計(jì)算(全局搜索、局部接觸和約束施加)。

基于單一靜態(tài)區(qū)域剖分的并行接觸算法概念簡(jiǎn)單，程序?qū)崿F(xiàn)容易。由于引入影像區(qū)存儲(chǔ)本地計(jì)算所需的全部信息，不需要每步計(jì)算都進(jìn)行數(shù)據(jù)通信。但由于區(qū)域剖分時(shí)未考慮接觸計(jì)算的影響，容易導(dǎo)致計(jì)算節(jié)點(diǎn)間嚴(yán)重的負(fù)載不均衡，算法并行擴(kuò)展性差[62]，很難擴(kuò)展到數(shù)千CPU核。

雙重區(qū)域分解方法中，單元計(jì)算和接觸計(jì)算分別基于主分區(qū)和次分區(qū)進(jìn)行。主分區(qū)通常采用圖剖分法進(jìn)行區(qū)域分割，次分區(qū)采用幾何剖分方法?；陔p重分區(qū)的并行接觸算法通信拓?fù)漭^為復(fù)雜，通常需要每步計(jì)算在主分區(qū)與次分區(qū)間進(jìn)行數(shù)據(jù)交互，但該類(lèi)方法可以較好地實(shí)現(xiàn)負(fù)載平衡。

基于雙重區(qū)域分解的并行接觸算法以圣地亞實(shí)驗(yàn)室Attaway等[63,64]的工作最為典型。他們對(duì)有限元計(jì)算采用靜態(tài)圖剖分方法，接觸界面區(qū)域采用動(dòng)態(tài)RCB算法作區(qū)域分割，大致保證了每個(gè)計(jì)算節(jié)點(diǎn)有相同的接觸從點(diǎn)和接觸主片，算法有效擴(kuò)展至數(shù)千CPU核。圖21給出了該算法的基本計(jì)算流程。

1.Send contact data from FE decomposition to old RCB decomposition

2.Perform parallel RCB to rebalance

3.Share RCB cut information with all processors

4.For all my surfaces

IF surface capture box extends beyond my RCB box,THEN

Determine which other processors need it

END IF

5.Send overlapping surfaces to nearby processors

6.Find contacts within corresponding RCB Box

CALL global search and local search routines

7.Send contact searching results to FE owners

圖21 并行接觸計(jì)算的區(qū)域分割

Fig.21 Outline of dynamic decomposition for parallel contacts

6 結(jié) 論

本文總結(jié)分析了顯格式有限元計(jì)算常用的接觸-碰撞算法，包括全局搜索、局部搜索、約束施加算法以及基于區(qū)域分解的并行計(jì)算方法。本文可為接觸-碰撞算法的深入研究提供基礎(chǔ)的理論參考。

接觸-碰撞算法，特別是局部接觸搜索算法，實(shí)現(xiàn)細(xì)節(jié)非常豐富。看似可忽略的細(xì)節(jié)常常決定了算法的健壯性與精度，程序?qū)崿F(xiàn)時(shí)必須予以充分考慮。

經(jīng)30多年的發(fā)展，接觸-碰撞算法已經(jīng)能夠解決大多數(shù)的工程問(wèn)題，但仍然存在一些問(wèn)題需要進(jìn)一步深入研究。

(1) 局部搜索精度。目前，接觸-碰撞算法主要是基于點(diǎn)-面接觸模型發(fā)展起來(lái)的，相關(guān)算法相對(duì)簡(jiǎn)單且效率高，但在處理不連續(xù)、非光滑接觸界面時(shí)易發(fā)生接觸遺漏，死區(qū)仍然是一個(gè)沒(méi)有完全解決的難點(diǎn)。面-面接觸模型有望避免涉及死區(qū)問(wèn)題，基于面-面接觸模型發(fā)展相應(yīng)的接觸搜索算法是接觸-碰撞問(wèn)題數(shù)值算法的一個(gè)研究重點(diǎn)。

(2) 全局搜索效率。全局搜索是影響接觸-碰撞問(wèn)題數(shù)值計(jì)算效率的關(guān)鍵因素。目前已有全局搜索算法仍不能滿(mǎn)足計(jì)算需求，仍然需要發(fā)展高效的搜索算法。

(3) 并行擴(kuò)展性。接觸-碰撞并行算法設(shè)計(jì)的關(guān)鍵在于全局搜索算法。目前，接觸搜索算法大多是針對(duì)串行計(jì)算而設(shè)計(jì)，缺乏內(nèi)在并行性。動(dòng)態(tài)接觸界面的動(dòng)態(tài)性、全局性與事先不確定性特點(diǎn)決定了并行計(jì)算時(shí)常需要大量的全局通信及不規(guī)則的通信拓?fù)洹Ａ硗?，接觸系統(tǒng)中實(shí)際發(fā)生接觸的部位往往是局部的，極易導(dǎo)致節(jié)點(diǎn)機(jī)的嚴(yán)重負(fù)載不平衡。這些因素使得接觸-碰撞并行算法的設(shè)計(jì)極具挑戰(zhàn)性，可用于數(shù)千CPU核的并行接觸算法仍處于發(fā)展中。

:

[1] Zhong Z H.FiniteElementProceduresforContact-ImpactProblems[M].New York:Oxford University Press,1993.

[2] Laursen T A.ComputationalContactandImpactMechanics[M].Berlin:Springer,2002.

[3] Yastrebov V A.NumericalMethodsinContactMechanics[M].John Wiley & Sons,2013.

[4] Belytschko T,Liu W K,Moran B,et al.NonlinearFiniteElementsforContinuaandStructures[M].John Wiley & Sons,2013.

[5] Shukla A，Ravichandran G,Rajapakse Y.DynamicFailureofMaterialsandStructures[M].New York：Springer,2010.

[6] Benson D J,Hallquist J O.A single surface contact algorithm for the post-buckling analysis of shell structures [J].ComputerMethodsinAppliedMechanicsandEngineering,1990,78(2):141-163.

[7] Heinstein M W,Attaway S W,Swegle J W,et al.A General-Purpose Contact Detection Algorithm for Nonlinear Structural Analysis Codes [R].SAND92-2141,Sandia National Laboratories,Albuquerque,1993.

[8] Heinstein M W,Mello F J,Attaway S W,et al.Contact-impact modeling in explicit transient dynamics [J].ComputerMethodsinAppliedMechanicsandEngineering,2000,187(3-4):621-640.

[9] Bourago N G,Kukudzhanov V N.A review of contact algorithms [J].MechanicsofSolids,2005,40(1):35-71.

[10] 鐘 陽(yáng),鐘志華,李光耀,等.機(jī)械系統(tǒng)接觸碰撞界面顯式計(jì)算的算法綜述 [J].機(jī)械工程學(xué)報(bào),2011,47(13):44-58.(ZHONG Yang,ZHONG Zhi-hua,LI Guang-yao,et al.Review on Contact algorithms calculating the contact-impact interface in mechanical system with explicit FEM [J].JournalofMechanicalEngineering,2011,47(13):44-58.(in Chinese))

[11] Oldenburg M,Nilsson L.The position code algorithm for contact searching [J].InternationalJournalforNumericalMethodsinEngineering,1994,37(3):359-386.

[12] Hallquist J O,Goudreau G L,Benson D J.Sliding interfaces with contact-impact in large -scale Lagrangian computations [J].ComputerMethodsinAppliedMechanicsandEngineering,1985,51(1):107-137.

[13] Hallquist J O.LS-DYNATheoryManual[M].Livermore Software Technology Corporation,2006.

[14] Hughes T J R,Taylor R L,Sackman J L,et al.A finite element method for a class of contact-impact problems [J].ComputerMethodsinAppliedMechanicsandEngineering,1976,8(3):249-276.

[15] Puso M A,Laursen T A.A mortar segment-to -segment contact method for large deformation so -lid mechanics [J].ComputerMethodsinAppliedMechanicsandEngineering,2004,193(6):601-629.

[16] Laursen T A,Puso M A,Sanders J.Mortar contact formulations for deformable -deformable contact:past contributions and new extensions for enriched and embedded interface formulations[J].ComputerMe-thodsinAppliedMechanicsandEngineering,2012,205:3-15.

[17] Belytschko T,Lin J I.A three -dimensional impact-penetration algorithm witherosion [J].Computers&Structures,1987,5(1-4):111-127.

[18] Zhong Z H,Nilsson L.A contact searching algorithm for general contact problems [J].Computers&Structures,1989,33(1):197-209.

[19] Zhong Z H,Nilsson L.A unified contact algorithm based on the territory concept [J].ComputerMethodsinAppliedMechanicsandEngineering,1996,130(1-2):1-6.

[20] Bruneel H C J,De Rycke I.QuickTrace:a fast algorithm to detect contact [J].InternationalJournalforNumericalMethodsinEngineering,2002,54(2):299-316.

[21] Yang B,Laursen T A.A contact searching algorithm including bounding volume trees applied to finite sliding mortar formulations [J].ComputationalMechanics,2008,41(2):189-205.

[22] Wang S P,Nakamachi E.The inside -outside contact search algorithm for finite element analysis [J].InternationalJournalforNumericalMethodsinEngineering,1997,40(19):3665-3685.

[23] Belytschko T,Neal M O.Contact-impact by the pinball algorithm with penalty and Lagrangian methods [J].InternationalJournalforNumericalMethodsinEngineering,1991,31(3):547-572.

[24] Belytschko T,Yeh I S.The splitting pinball method for contact-impact problems [J].ComputerMethodsinAppliedMechanicsandEngineering,1993,105(3):375-393.

[25] 王福軍.沖擊接觸問(wèn)題有限元法并行計(jì)算及其工程應(yīng)用 [D].清華大學(xué),2000.(WANG Fu-jun.Parallel Computation of Contact-Impact Problems with FEM and Its Engineering [D].Tshinghua University,2000.(in Chinese))

[26] Wang F J,Cheng J G,Yao Z H.FFS contact sear-ching algorithm for dynamic finite element analysis [J].InternationalJournalforNumericalMethodsinEngineering,2001,52(7):655-672.

[27] Abaqus 6.12TheoryManual[M].Dassault Systemes Simulia Corporation,2012.

[28] Whirley R G,Engelmann B E.Automatic Contact in DYNA3D for Vehicle Crashworthiness [R].Lawrence Livermore National Laboratories,CA,1993.

[29] Whirley R G,Engelmann B E.Automatic contact algorithm in DYNA3D for crashworthiness and impact problems [J].NuclearEngineeringandDesign,1994,150(2):225-233.

[30] Diekmann R,Hungershofer J,Lux M,et al.Efficient contact search for finite element analysis [A].European Congress on Computational Methods in Applied Sciences and Engineering[C].Barchlona,Spain,2000.

[31] 柳 陽(yáng),柳 明,陳成軍.簡(jiǎn)化相鄰依賴(lài)關(guān)系的mortar面面接觸算法研究[A].中國(guó)計(jì)算力學(xué)大會(huì)[C].貴陽(yáng),2014.(LIU Yang,LIU Ming,CHEN Cheng-jun.Investigation of the mortar segment-to-segment contact algorithm with the neighbor-dependent relationship simplified[A].China Computational Mecha-nics Conference[C].Guiyang，2014.(in Chinese))

[32] Papadopoulos P,Taylor R L.A simple algorithm for three -dimensional finite element analysis of contact problems [J].Computers&Structures,1993,46(6):1107-1118.

[33] Mahadevaiah U.Development Andvalidation of New Algorithms to Improve Contact Detection and Robustness in Finite Element Simulations [D].The George Washington University,2009.

[34] Wang F J,Cheng J G,Yao Z H.A contact searching algorithm for contact-impact problems [J].ActaMechanicaSinica,2000,16(4):374-382.

[35] Munjiza A,Andrews K R F.NBS contact detection algorithm for bodies of similar size [J].InternationalJournalforNumericalMethodsinEngineering,1998,43(1):131-149.

[36] Chen H,Lei Z,Zang M Y.LC-Grid:a linear global contact search algorithm for finite element analysis [J].ComputationalMechanics,2014,54(5):1285-1301.

[37] Williams J R,Perkins E,Cook B.A contact algorithm for partitioning N arbitrary sized objects [J].EngineeringComputations,2004,21(21314):235-248.

[38] Zavarise G,De Lorenzis L.The node -to -segment algorithm for 2D frictionless contact:classical formulation and special cases[J].ComputerMethodsinAppliedMechanicsandEngineering,2009,198(41):3428-3451.

[39] EI-Abbasi N,Meguid S A,Czekanski A.On the mo -delling of smooth contact surfaces using cubic splines[J].InternationalJournalforNumericalMethodsinEngineering,2001,50(4):953-967.

[40] Corbett C J,Sauer R A.NURBS -enriched contact finite elements[J].ComputerMethodsinAppliedMechanicsandEngineering,2014,275:55-75.

[41] Puso M A,Laursen T A.A 3D contact smoothing method using Gregory patches[J].InternationalJournalforNumericalMethodsinEngineering,2002,54(8):1161-1194.

[42] Tylor L M,Flanagan D P.PRONTO 3D:A Three -Dimensional Transient Solid Dynamics Program [R].Sandia National Lab.Albuquerque,1989.

[43] Carpenter N J,Taylor R L,Katona M G.Lagrange constraints for transient finite element surface contact [J].InternationalJournalforNumericalMethodsinEngineering,1991,32(1):103-128

[44] Sha D,Tamma K K,Li M.Robust explicit computational developments and solution strategies for impact problems involving friction[J].InternationalJournalforNumericalMethodsinEngineering,1996,39(5):721-739.

[45] Wang F J,Wang L P,Cheng J G,et al.Contact force algorithm in explicit transient analysis using finite -element method [J].FiniteElementsinAnalysisandDesign,2007,43(6):580-587.

[46] Chen H,Zhang Y X,Zang M Y,et al.An explicit Lagrange constraint method for finite element analysis of frictionless 3D contact/impact problems [J].AppliedMechanicsandMaterials,2014,553:751-756.

[47] 金先龍,李淵印.結(jié)構(gòu)動(dòng)力學(xué)并行計(jì)算方法及應(yīng)用[M].北京：國(guó)防工業(yè)出版社,2008.(JIN Xian-long,LI Yuan-yin.ParallelCalculationMethodandApp-licationofStructuralDynamics[M].Beijing:National Defence Industry Press,2008.(in Chinese))

[48] Malone J G,Johnson N L.A parallel finite element contact/impact algorithm for non-linear explicit transient analysis:Part I—The search algorithm and contact mechanics [J].InternationalJournalforNumericalMethodsinEngineering,1994,37(4):559-590.

[49] Malone J G,Johnson N L.A parallel finite element contact/impact algorithm for non-linear explicit transient analysis:Part II—Parallel implementation [J].InternationalJournalforNumericalMethodsinEngineering,1994,37(4):591-603.

[50] Watson B C,Noor A K.Large -scale contact/impact simulation and sensitivity analysis on distributed-memory computers [J].ComputerMethodsinAppliedMechanicsandEngineering,1997,141(3):373-388.

[51] Zhong Z H,Nilsson L.Contact-impact algorithms on parallel computers [J].NuclearEngineeringandDesign,1994,150(2):253-263.

[52] Har J,Fulton R E.Aparallel finite element procedure for contact-impact problems [J].EngineeringwithComputers,2003,19(2):67-84.

[53] 白小勇,何穎波,陳成軍.顯式有限元中的一種并行接觸算法[J].計(jì)算物理,2011,28(3):341-346.(BAI Xiao -yong,HE Ying-bo,CHEN Cheng-jun.A parallel contact detection algorithm for explicit finite element analysis [J].ChineseJournalofComputationalPhysics,2011,28(3):341-346.(in Chinese))

[54] 陳成軍,柳 陽(yáng),張?jiān)?等.基于PANDA的并行顯式有限元程序開(kāi)發(fā)[J].計(jì)算力學(xué)學(xué)報(bào),2011,28(s1):204-207，214.(CHEN Cheng-jun,LIU Yang,ZHANG Yuan-Zhang,et al.Programming of parallel explicit finite element based on PANDA [J].ChineseJournalofComputationalMechanics,2011,28(s1):204-207，214.(in Chinese))

[55] 寇哲君,程建鋼,姚振漢.可擴(kuò)展的沖擊-接觸并行計(jì)算研究[J].計(jì)算力學(xué)學(xué)報(bào),2003,20(3):325-328.(KOU Zhe -jun,CHENG Jian-gang,YAO Zhen-han.Study of scalable parallel computation for contact-impact problems [J].ChineseJournalofComputationalMechanics,2003,20(3):325-328.(in Chinese))

[56] Wang F J,Feng Y T,Owen D R J.Interprocessor communication schemes in parallel finite -discrete element analysis on PC clusters [J].EngineeringComputations,2003,20(8):1065-1084.

[57] Wang F J,Feng Y T,Owen D R J.Parallelisation for finite -discrete element analysis in a distributed-memory environment [J].InternationalJournalofComputationalEngineeringScience,2004,5(1):1-23.

[58] 王福軍,王利萍,程建鋼,等.并行有限元計(jì)算中的接觸算法[J].力學(xué)學(xué)報(bào),2007,39(3):422-427.(WANG Fu-jun,WANG Li-ping,CHENG Jian-gang,et al.A contact algorithm for parallel computation of FEM [J].ChineseJournalofTheoreticalandAppliedMechanics,2007,39(3):422-427.(in Chinese))

[59] 亓文果,金先龍,張曉云.沖擊-接觸問(wèn)題有限元仿真的并行計(jì)算[J].振動(dòng)與沖擊,2006,25(4):68-72.(Qi Wen-guo,JIN Xian-long,ZHANG Xiao -yun.Study on parallel algorithm for finite element simulation of contact-impact problems [J].JournalofVibrationandShock,2006,25(4):68-72.(in Chinese))

[60] 杜新光,金先龍,陳向東.基于并行計(jì)算的大跨度斜拉橋行車(chē)安全分析[J].振動(dòng)與沖擊,2010,29(7):5-8,78.(DU Xin-guang,JIN Xian-long,CHEN Xiang-dong.Simulation analysis for running safety of a light-rail train on a long span cable -stayed bridge based on parallel computation [J].JournalofVibrationandShock,2010,29(7):5-8,78.(in Chinese))

[61] Paik S H,Moon J J,Kim S J,et al.Parallel perfor-mance of large scale impact simulations on Linux cluster super computer [J].Computers&Structures,2006,84(10):732-741.

[62] Attaway S W,Hendrickson B A,Plimpton S J,et al.A parallel contact detection algorithm for transient solid dynamics simulations using PRONTO3D [J].ComputationalMechanics,1998,22(2):143-159.

[63] Brown K,Attaway S,Plimpton S,et al.Parallel stra-tegies for crash and impact simulations [J].ComputerMethodsinAppliedMechanicsandEngineering,2000,184(2-4):375-390.

[64] Hoover C G,Groot A J,Sherwood R J.Parallel Contact Algorithms for Explicit Finite Element Analysis [R].UCRL-JC-130461,Lawrence Livermore National Lab.Livermore,CA,1998.

[65] Pierce T G.A Parallel Algorithm for Contact in a Finite Element Hydrocode [R].Lawrence Livermore National Lab.Livermore,CA,2003.

[66] Pierce T,Rodrigue G.A parallel two -sided contact algorithm in ALE3D [J].ComputerMethodsinAppliedMechanicsandEngineering,2005,194(27-29):3127-3146.

[67] 陳成軍,柳 明.大規(guī)模并行顯式有限元軟件開(kāi)發(fā)[R].GF-A0188466G,中國(guó)工程物理研究院,2015.(CHEN Cheng-jun,LIU Ming.Development of Large Scale Parallel Explicit Finite Software[R].GF-A0188466G,China Academy of Engineering Physics,2015.(in Chinese))

[68] 肖永浩.沖擊-接觸問(wèn)題的并行計(jì)算研究[D].中國(guó)工程物理研究院,2013.(XIAO Yong-hao.Parallel Computing Research on Impact-Contact Problems [D].Chinese Academy of Engineering Physics,2013.(in Chinese))

[69] Lin Y Y,Wang J.Performance of the hybrid LS -DYNA on crash simulation with the multicore architecture[A].7t hEuropean LS -DYNA Conference[C].2009.

[70] Makino M.The performance of 10-million element car model by MPP version of LS -DYNA on fujitsu PRIMEPOWER [A].10t hInternational LS -DYNA Users Conference[C].Michigan，2008.

[71] Kenshiro K,Mitsuhiro M.The performance of car crash simulation by LS -DYNA hybrid parallel version on fujitsu FX1 [A].11t hInternational LS -DYNA Users Conference[C].2010.