接觸-碰撞有限元計算算法綜述

接觸-碰撞有限元計算算法綜述

1相關(guān)研究現(xiàn)狀實際設(shè)計中存在著變形體之間的接觸沖突（或動態(tài)接觸）。接觸-碰撞數(shù)值算法的研究大致起始于20世紀(jì)70年代,迄今已有很多學(xué)者開展了相關(guān)研究,發(fā)展了許多優(yōu)秀算法。Bourago等基于拉格朗日框架的顯式有限元是處理變形體間接觸-碰撞最常用的數(shù)值方法,因此本文主要針對顯式拉格朗日有限元相關(guān)的接觸-碰撞算法進行總結(jié)分析。2接觸界面的色散2.1處理工具中界面接觸的描述變形體間接觸-碰撞問題的有限元理論與基本數(shù)值方法可參考文獻接觸-碰撞算法可處理任意多個物體間的相互作用,多物體間的接觸可歸結(jié)為物體間的兩兩作用。不失一般性地,本文以兩個物體間的接觸為例給出其數(shù)學(xué)描述。如圖1所示,相互接觸的兩個變形體A和B,t時刻的構(gòu)形分別標(biāo)識為Ω接觸-碰撞系統(tǒng)除應(yīng)滿足通常的連續(xù)體控制方程外以及接觸界面上還應(yīng)滿足一定的運動學(xué)條件與動力學(xué)條件即兩物體間的間隙要求大于或等于0,0值對應(yīng)于接觸狀態(tài),g如果不考慮接觸面的焊接或粘接作用,界面法向作用力p2.2接觸-碰撞問題的拉格朗日乘子法離散方程與不含接觸的問題相比,接觸-碰撞系統(tǒng)僅需在接觸界面上額外引入約束條件(1~3)。因此,可方便地建立接觸-碰撞系統(tǒng)的控制方程弱形式式中σ對應(yīng)于接觸界面約束的罰函數(shù)法弱形式。若p即為拉格朗日乘子法的弱形式。對方程(4)進行有限元離散,可得半離散方程:式中M為質(zhì)量矩陣,a中心差分法是顯式有限元最常用的積分方法,此時方程(7)的時域離散可概括為式中Δt采用顯式積分求解時,方程(8)右端項中唯一的未知量為t接觸-碰撞問題中,接觸面通常不能事前確定,并且隨時間動態(tài)變化,因此接觸力計算可進一步分解為接觸搜索與接觸約束施加兩個子問題。接觸搜索的目的是確定系統(tǒng)中哪些部位發(fā)生了接觸或者哪些原已接觸的部位發(fā)生了分離或滑移,接觸約束施加則根據(jù)系統(tǒng)當(dāng)前接觸狀態(tài)計算出界面接觸力以滿足接觸約束條件。2.3接觸界面離散模型有限元方法中,相互接觸的變形體(接觸體)離散為有限單元,接觸界面即為變形體外表面的單元面集合或節(jié)點集合。這樣,接觸體間的作用轉(zhuǎn)化為離散的節(jié)點間、節(jié)點與單元面或單元面間的相互作用。為便于問題描述,本文采用Hallquist等根據(jù)接觸界面離散方式的不同,界面模型可分為三類。(1)點-點NTN(node-to-node)模型(2)點-面NTS(node-to-segment)模型(3)面-面STS(segment-to-segment)模型接觸界面離散方式在很大程度上決定了接觸界面搜索方式和接觸力的計算方法。鑒于NTS模型的廣泛應(yīng)用,本文將針對基于該類模型發(fā)展的相關(guān)算法進行討論。3接觸搜索的實現(xiàn)如上文所述,接觸界面的確定是接觸-碰撞問題數(shù)值計算中最耗時的部分。如采用最直接的暴力搜索方法確定接觸對,其計算量將為O(N為加速接觸對確定,接觸搜索通常分全局搜索和局部搜索兩個步驟進行。全局搜索利用某些簡單的準(zhǔn)則快速排除掉系統(tǒng)中不可能發(fā)生接觸的部位,或挑選出最可能發(fā)生接觸的從點與主片對,即確定出接觸測試對。局部搜索則從接觸測試對中精確確定出真實發(fā)生接觸的從點與主片對,即接觸對。一般而言,全局搜索決定了接觸計算效率;局部搜索主要影響接觸計算精度,同時對計算效率也有一定影響。典型的全局搜索算法有主從面算法3.1整個搜索方法3.1.1確定接觸部位主從面算法主從面算法的計算復(fù)雜度仍為O(N主從面算法對于含有少量接觸,且事前可大致確定出接觸部位的問題非常有效。但該算法只能處理兩個表面間的接觸,不能用于單一曲面接觸,也不能用于含材料斷裂破壞的侵蝕接觸分析,甚至不能用于相對初始構(gòu)形有嚴(yán)重變形的接觸問題3.1.2節(jié)點分組優(yōu)化桶排序法是目前應(yīng)用最廣泛的全局搜索算法,該類算法的基本思想是桶排序法通過將節(jié)點分組(同一桶中的節(jié)點視為一組)極大縮小了全局接觸搜索范圍,從而大幅提高計算效率。該算法不需要預(yù)先指定接觸主從面,因此具有良好的適用性,可以用于多體接觸對三維空間問題,經(jīng)典桶排序法Heinstein等3.1.3接觸域和接觸搜索級域算法層域是某層級實體占據(jù)的空間,為避免每個迭代步都進行全局搜索,可將此區(qū)域適當(dāng)擴展而形成擴展域。當(dāng)從點位于某主點、線或片的接觸域內(nèi)時,認(rèn)為該從點與相應(yīng)的主點、線或片接觸。點、線與片接觸域的定義如圖4所示。片的接觸域由其邊、外法線方向及某種特征厚度確定,線的接觸域由其外法線方向及其所在的片厚度確定,點的接觸域則由其所在線的外法線方向及相應(yīng)片厚度確定。接觸搜索時,先在較高級的層域間進行,若兩個域存在重疊,則繼續(xù)進行下一級層域的相交檢查,直至確定出接觸測試對;若兩域不存在重疊,則兩接觸實體間不存在接觸關(guān)系,接觸搜索停止。級域算法利用了接觸體的層級結(jié)構(gòu),在較高級實體中排除了不必要的接觸檢查,理論上可大幅提高接觸搜索效率,這對于同級域重疊部分較小的接觸問題非常有效。但對于同級實體存在較多重疊的接觸問題,如單曲面接觸,計算效率不如位置碼算法3.1.4基于位置碼的接觸分析Oldenburg等定義了主片域后,位于主片域內(nèi)從點的三維空間搜索利用位置碼轉(zhuǎn)換到一維空間中進行:接觸界面空間分割成若干桶;根據(jù)節(jié)點所在的桶,對每一個節(jié)點賦予一個位置碼;根據(jù)位置碼與節(jié)點坐標(biāo)將三維空間中的節(jié)點排序到一維數(shù)組(稱為位置碼向量);對每一個主片,確定出主片域占據(jù)的桶,利用二分查找算法在位置碼向量中找出這些桶中的節(jié)點(潛在接觸從點)。同時,位置碼算法采用層級概念定義了三層接觸實體,即接觸面、接觸片和接觸點,每一個接觸面定義一個位置碼向量。位置碼算法將三維空間的節(jié)點排序和搜索映射到一維數(shù)組中進行,并利用二分查找法加速接觸從點的搜索速度,具有很高的計算效率(計算復(fù)雜度為O(Nlog位置碼算法的效率有一定的網(wǎng)格方向依賴性,這是由桶的編碼方式?jīng)Q定的。Diekmann等3.1.5創(chuàng)建接觸片幾何信息樹形算法是新近發(fā)展的一類全局搜索方法。該方法的基本思想是,采用樹形數(shù)據(jù)結(jié)構(gòu),如二叉樹和四叉樹,存儲接觸片幾何信息,利用樹的快速遍歷方法確定出有可能參與接觸的從點與主片或者從片與主片對。針對沖壓成形中板料與模具間的動態(tài)接觸,Bruneel等與經(jīng)典桶排序算法3.1.6接觸搜索算法的發(fā)展Papadopoulos等Mahadevaiah近年,針對離散粒子接觸系統(tǒng),發(fā)展了大量優(yōu)秀的接觸搜索算法,雖然不能直接用于有限元接觸計算,但為有限元接觸搜索算法的設(shè)計提供了重要借鑒。王福軍等3.2主要使用的算法與全局搜索相比,局部搜索算法的研究相對較少,目前主要采用的方法有點面算法和小球類算法兩大類。另有少量搜索算法主要是針對兩類算法的不足而提出的改進方案,如內(nèi)外算法、曲面構(gòu)造法及靜動態(tài)檢測法等。3.2.1點面搜索算法的不足點面算法假定從點n式中c式中m=(c式中t為從點n式中n若l<0,表明從點對主片有侵入,則根據(jù)接觸約束算法在從點與主片間施加接觸約束力。當(dāng)接觸界面光滑且網(wǎng)格質(zhì)量較好時,點面搜索算法非常有效,這也是該方法受到廣泛采用的原因。但算法存在如下不足。(1)經(jīng)典點面算法(2)以主片外法線方向作為從點的投影方向,存在投影點不唯一或無投影點等奇異性問題(3)接觸面網(wǎng)格質(zhì)量較差時,如細(xì)長比過大,可能會發(fā)生接觸搜索遺漏,如圖9所示。(4)由于該方法侵入檢查是針對從點與主片的關(guān)系進行的,因此不能處理如圖10所示的邊-面和面-面接觸情況。3.2.2從點部分植物的侵入量計算面向金屬成形模擬,針對經(jīng)典點面算法內(nèi)外算法中,以從點的外法線方向作為從點對主片的投影方向(同時也是接觸力施加方向)。從點外法向定義為與該節(jié)點相連的接觸片法線方向的平均值,如圖11所示。從點相對于主片一條邊的內(nèi)或外狀態(tài),通過從點與該邊構(gòu)成的三角形在主片上的投影來確定,如圖12所示,若Δ如果從點位于主片內(nèi)部,則利用式(18)計算從點對主片的侵入量,式中x為投影點,利用式(17)可以解析給出;x內(nèi)外算法的優(yōu)勢在于,接觸點可以唯一確定,避免了從點對主片的多重接觸,并且接觸點計算不需要迭代。因此,該算法具有很好的健壯性與較高的計算效率。由于內(nèi)外算法假定接觸點處主從面的法向共軸(但方向相反),當(dāng)從點位于兩個低階單元的交界線上時,尤其是當(dāng)網(wǎng)格不夠精細(xì)時,內(nèi)外算法給出的接觸方向可能會與實際存在較大偏差3.2.3接觸力的計算精度點面算法存在投影奇異性的根本原因是低階有限元近似使得接觸界面僅C構(gòu)造曲面法不僅可以克服接觸盲區(qū),也有利于提高接觸力計算精度,消除接觸滑移時接觸力的非物理震蕩。關(guān)于構(gòu)造曲面法,已有很多學(xué)者開展相關(guān)研究,提出了一些光滑化算法,如三次樣條法FFS方法利用接觸片節(jié)點坐標(biāo)構(gòu)造參數(shù)化光滑曲面片。如圖13所示,任意一個由四條邊構(gòu)成的接觸片,參數(shù)化光滑曲面可表示為式中x為曲面的整體坐標(biāo),u和w為曲面的局部坐標(biāo),cFFS算法構(gòu)造的曲面具有C3.2.4移動點接觸檢查經(jīng)典點面局部搜索算法(1)動態(tài)檢測。根據(jù)從點和主片的當(dāng)前位置及下一時間步的預(yù)測速度,利用從點和主片接觸時刻的共面關(guān)系,確定出接觸時刻與接觸點。三維情況下,這種接觸檢查簡化為判斷一個移動三角形和一個移動點之間的接觸問題(對于四邊形主片,可以利用主片中心點與四個邊將主片分解為四個三角形),如圖14所示。求得接觸時刻與接觸點后,進一步判斷接觸時刻是否滿足步長穩(wěn)定性要求,及接觸點是否位于主片內(nèi)部,如滿足則發(fā)生接觸。(2)靜態(tài)檢測。靜態(tài)檢測基于預(yù)測構(gòu)形進行,從點回退方向根據(jù)接觸面的凹凸性質(zhì)分別處理。對于凹面,接觸點為從點在主片上的投影點,或者是主片頂點,從點對主片的侵徹方向如圖15所示;對于凸面,侵徹方向認(rèn)為沿著上一迭代步主片的法向,即圖16所示的方向m^3.2.5殼單元間的接觸檢查與分裂材料小球算法(pinballalgorithm)是Belytschko等引入小球概念后,單元間的侵入檢查變得極為簡單。假定主面單元與從面單元等效小球的半徑分別為R侵入深度g由式(21)計算得到,式中d=C小球算法可用于實體單元或薄殼單元的任意類型接觸分析。但用于薄殼單元時,若殼元厚度相對于其特征尺寸很小,則小球算法精度會非常差;若殼單元間存在初始接觸,小球算法可能會失效。為此,Belytschko等與小球算法類似,分裂小球算法以小球代替每一個參與接觸的表面單元,但小球的半徑要足夠大以完全覆蓋相應(yīng)單元,該小球稱為父球。父球間若有重疊區(qū)則表示兩單元間可能存在接觸,此時父球分裂為下一級的子球進一步作重疊區(qū)檢查,直至最后一級子球的半徑與殼厚度相當(dāng)。四邊形殼元的小球?qū)蛹壢鐖D18所示,接觸檢查與小球逐級分裂的一維示意如圖19所示。小球算法和分裂小球算法以統(tǒng)一的方式自動處理邊-邊、邊-面及面-體等接觸類型,算法便于向量化和并行計算;單元間的接觸檢查與侵入量計算簡單高效。小球算法和分裂小球算法適用于碰撞和侵徹等問題的模擬。由于算法未考慮接觸界面間的摩擦效應(yīng),不能應(yīng)用于滑移和摩擦為關(guān)鍵因素的接觸分析中。另外,當(dāng)接觸界面單元特征尺寸差別較大時,接觸力的計算精度較差,不準(zhǔn)確的幾何近似和小球體積不變假設(shè)限制了該算法的應(yīng)用4接觸限制算法接觸約束算法,也稱為接觸力算法,基本方法有拉氏乘子法和罰函數(shù)法4.1接觸力的計算罰函數(shù)法的基本原理是,如果從點對主片沒有侵入則不做處理;如果從點侵入主片,則在從點與主片間引入法向接觸力,將從點推回到變形后的接觸主片上,以滿足不可侵入條件。罰函數(shù)法不增加系統(tǒng)方程的自由度數(shù),接觸力計算無需聯(lián)立求解方程組,程序?qū)崿F(xiàn)簡單且計算效率高,廣泛應(yīng)用于顯式有限元程序中,如LS-Dyna通用程序一般會對罰因子進行規(guī)則化處理,以便用戶選取相對合理的罰參數(shù)。以LS-Dyna式中l(wèi)為從點對主片的侵入深度式中n規(guī)則化后的接觸罰參數(shù)f式中ue7884.2顯式有限元接觸力的計算拉氏乘子法中,以接觸對間的接觸力作為乘子來限制接觸體間的相互侵入,不可侵入條件可精確滿足。但該方法引入了新的未知變量,增加了方程組的自由度數(shù);接觸力計算需要聯(lián)立方程求解,這與顯式有限元計算不相容。為了在顯式有限元中應(yīng)用乘子法進行接觸力計算,必須對經(jīng)典的拉氏乘子法做必要的格式修正。Carpenter等防御節(jié)點法的基本思想是,將從點與主片之間的接觸作用轉(zhuǎn)化為從點與防御節(jié)點之間的接觸,防御節(jié)點攜帶了主片的運動學(xué)與動力學(xué)信息,如位移、速度及節(jié)點力等。防御節(jié)點位于接觸點,表征主片的運動,因此其速度和加速度可以利用主片的節(jié)點速度和加速度描述,式中φ為插值形函數(shù)。防御節(jié)點的運動為F為除接觸力以外的力,f為從點與防御節(jié)點間的接觸力。計算f時,施加不可侵入條件。得到f后,接觸力就可以分布到主點上。5接觸界面的全局動態(tài)性接觸-碰撞問題通常采用顯式有限元進行計算。顯式雖然避免了內(nèi)存的龐大需求,但積分穩(wěn)定性要求時間步長必須足夠小,這導(dǎo)致接觸-碰撞問題求解通常需要大量迭代,對于大規(guī)?；驈?fù)雜幾何構(gòu)形問題非常耗時,因此并行計算是接觸-碰撞數(shù)值分析的必然選擇。區(qū)域分解(spatialdomaindecompositionmethod)是實現(xiàn)接觸-碰撞問題大規(guī)模并行計算的最有效方式。區(qū)域分解的基本思想與不含接觸的動力學(xué)問題或準(zhǔn)靜態(tài)接觸問題不同,接觸-碰撞問題的接觸域通常具有全局性與動態(tài)性,計算節(jié)點之間的通信關(guān)系在系統(tǒng)求解之前不能確定。如圖20所示的一個簡單接觸問題,開始時刻下方物體上表面的左前部位與上面物體發(fā)生接觸,而在某個時刻下方物體上表面的右后部位與上面物體發(fā)生接觸。接觸界面的全局性和動態(tài)性對接觸-碰撞問題的區(qū)域分割和計算節(jié)點間的通信關(guān)系設(shè)計帶來極大挑戰(zhàn)。(1)對有限元計算(主要是內(nèi)力計算)高效的某種區(qū)域分解在接觸計算中往往不能取得好的并行性能,通常需要針對接觸而引入另外的區(qū)域分割。(2)動態(tài)接觸的全局性導(dǎo)致計算節(jié)點間存在大量的數(shù)據(jù)交互,最大限度地減小通信量才能保證算法的擴展性與并行效率。(3)接觸對動態(tài)變化,局部接觸搜索和約束施加在不同時刻可能由不同的計算節(jié)點完成,極易造成各節(jié)點機負(fù)載的不均衡?；趨^(qū)域分解的并行接觸算法研究始于20世紀(jì)90年代。Malone等上述算法僅適用于少量CPU核的集群環(huán)境,有良好擴展性、適用于大型并行系統(tǒng)(數(shù)千CPU核)的并行接觸算法研究還較少。目前該方面的研究主要限于美國Sandia國家實驗室基于區(qū)域分解實現(xiàn)接觸-碰撞問題的粗粒度并行有單一區(qū)域分解法和雙重區(qū)域分解法兩種并行方案。雙重區(qū)域分解法中,有限元計算采用一種子域剖分(稱為主分區(qū)),接觸計算采用另外一套分區(qū)(稱為次分區(qū));而單一區(qū)域分解法,接觸區(qū)域不作單獨的子區(qū)分割,接觸的并行計算基于主分區(qū)進行。通常,雙重區(qū)域分解法可獲得更好的負(fù)載均衡,但節(jié)點機間的通信復(fù)雜;單一分區(qū)法通信拓?fù)浜唵?但易導(dǎo)致負(fù)載不均衡。5.1靜態(tài)區(qū)域剖分法設(shè)計單一區(qū)域分解方法中,接觸計算的并行基于主分區(qū)進行。接觸算法并行的關(guān)鍵在于,每個計算節(jié)點可有效得到本地計算所需的全部信息。因此,該類算法通常引入影像區(qū)(ghost)或接觸域(contactdomain)存儲本地接觸計算需要但駐留在其他節(jié)點機上的信息。基于單一區(qū)域分解的并行接觸算法設(shè)計相對簡單,本文以Malone等(1)不考慮接觸界面,對求解域作區(qū)域剖分,分配給各計算節(jié)點。所有計算節(jié)點作如下相同的計算。(2)對任意計算節(jié)點PA,構(gòu)建覆蓋該計算節(jié)點單元幾何空間的包圍盒(boundingbox),并傳遞到其他計算節(jié)點,同時接收其他計算節(jié)點的包圍盒信息。(3)利用體積檢查(volumechecking)判斷PA包圍盒與其他包圍盒是否有重疊。若無重疊,則兩個計算節(jié)點的單元間無潛在接觸;若有重疊部分,則將與PA包圍盒有重疊的計算節(jié)點上的單元信息收集到PA的接觸域,同時將PA單元傳遞給對方計算節(jié)點的接觸域。(4)對每個計算節(jié)點利用串行接觸算法進行接觸計算(全局搜索、局部接觸和約束施加)。基于單一靜態(tài)區(qū)域剖分的并行接觸算法概念簡單,程序?qū)崿F(xiàn)容易。由于引入影像區(qū)存儲本地計算所需的全部信息,不需要每步計算都進行數(shù)據(jù)通信。但由于區(qū)域剖分時未考慮接觸計算的影響,容易導(dǎo)致計算節(jié)點間嚴(yán)重的負(fù)載不均衡,算法并行擴展性差5.2分區(qū)區(qū)域分割方法雙重區(qū)域分解方法中,單元計算和接觸計算分別基于主分區(qū)和次分區(qū)進行。主分區(qū)通常采用圖剖分法