面面接觸算法研究

面面接觸算法研究

1接觸算法的計(jì)算方法接觸問題是工程中常見的力學(xué)問題。接觸問題的特點(diǎn)和難點(diǎn)是接觸限制和接觸力的未知。由于接觸界面的區(qū)域大小、位置以及接觸狀態(tài)都是未知的、隨時(shí)間變化的，因此接觸問題表現(xiàn)出顯著的非線性特征。接觸問題的非線性決定了接觸分析過程中需要經(jīng)常插入接觸界面的搜尋判定。有限元法為分析接觸問題提供了有效的工具，在進(jìn)行有限元求解接觸問題時(shí)，接觸體被離散為空間單元集，接觸界面的相互作用被轉(zhuǎn)化為離散的單元表面之間的作用或者離散的節(jié)點(diǎn)與單元表面的作用。接觸問題的動(dòng)力分析主要包括選取時(shí)間積分方案和接觸算法，其中接觸算法包含接觸搜索算法和接觸力算法。接觸搜索的目的是確定整個(gè)系統(tǒng)中有哪些部位發(fā)生了接觸或者有哪些原已接觸的部分發(fā)生了滑移或脫離。在接觸動(dòng)力分析中，接觸搜索占較大的計(jì)算量，因此接觸搜索算法的計(jì)算效率至關(guān)重要。由于接觸搜索是接觸力計(jì)算的基礎(chǔ)，所以其計(jì)算精度也非常關(guān)鍵。在一般情況下接觸搜索需要先進(jìn)行全局搜索，再進(jìn)行局部搜索。通過全局搜索先粗略找到圍繞特殊點(diǎn)所有潛在可能的接觸單元面。全局算法有主從面算法局部搜索目的是計(jì)算節(jié)點(diǎn)到單元面的精確距離，找出接觸投影點(diǎn)所在單元面的相對(duì)位置關(guān)系，從而確定接觸節(jié)點(diǎn)與單元面間的接觸狀態(tài)。局部搜索算法主要有點(diǎn)面算法在完成接觸搜索之后，接觸力求解通常采用拉格朗日乘子法或罰函數(shù)法強(qiáng)震作用下的高混凝土壩系統(tǒng)災(zāi)變過程伴隨著材料和接觸復(fù)合非線性問題，對(duì)其進(jìn)行全過程精細(xì)化數(shù)值模擬和全面深入的抗震安全評(píng)價(jià)，需要求解未知量高達(dá)百萬甚至千萬級(jí)，盡管各種商業(yè)軟件在求解一般常見問題時(shí)可以顯示出所謂的強(qiáng)大功能，但在計(jì)算方法和計(jì)算實(shí)施方案上很難適用于解決高混凝土壩復(fù)雜工況下地震響應(yīng)分析所遇到的復(fù)雜而又特殊的問題，因此，高性能并行計(jì)算是解決該問題的必然選擇。本文擬在研究現(xiàn)有接觸算法的基礎(chǔ)上，提出多體接觸問題快捷高效的面面接觸算法。這種算法既可解決搜索盲區(qū)問題，又能回避求解非線性投影方程。本文算法采用拉格朗日乘子法動(dòng)力接觸方程的增量形式求解接觸力，另外，為了便于并行計(jì)算及其程序編制，本文還根據(jù)接觸面分布特點(diǎn)進(jìn)行預(yù)分區(qū)處理，并基于提出的多體接觸問題面面接觸算法，開發(fā)實(shí)現(xiàn)其算法的FORTRAN源代碼程序，為后續(xù)實(shí)現(xiàn)高混凝土壩和巖體高邊坡靜動(dòng)力分析的并行計(jì)算提供支持。2搜索法的分析2.1全局搜索過程的實(shí)體面求解本文將主從面算法與位碼算法相結(jié)合，實(shí)現(xiàn)接觸面的全局搜索。利用主從面算法本文構(gòu)建了描述“接觸面對(duì)”的主面單元和從面單元的共享實(shí)常數(shù)組：主面單元用正整數(shù)表示，從面單元用其對(duì)應(yīng)的負(fù)整數(shù)表示，當(dāng)絕對(duì)值相同時(shí)主面單元和從面單元構(gòu)成“接觸面對(duì)”；不同的“接觸面對(duì)”通過不同的實(shí)常數(shù)定義；一組實(shí)常數(shù)可以對(duì)應(yīng)多個(gè)邊界面。以下是對(duì)一組“接觸面對(duì)”的全局搜索步驟：（1）搜索預(yù)處理：統(tǒng)計(jì)“接觸面對(duì)”單元面上的主節(jié)點(diǎn)和從節(jié)點(diǎn)個(gè)數(shù)，計(jì)算各單元面外法向量及各節(jié)點(diǎn)外法向量。其中節(jié)點(diǎn)外法向量由包含該節(jié)點(diǎn)的單元法向量取平均得到。（2）根據(jù)“接觸面對(duì)”的位置關(guān)系，構(gòu)造包含接觸單元面的最小長(zhǎng)方體。在接觸問題的有限元計(jì)算過程中，全局搜索是相對(duì)耗時(shí)較大的計(jì)算環(huán)節(jié)，為了避免在每個(gè)時(shí)間步都進(jìn)行全局搜索，適當(dāng)擴(kuò)大該長(zhǎng)方體的區(qū)域，確定接觸搜索范圍。（3）確定立方格的網(wǎng)格尺寸，立方格尺寸與接觸面中平均的單元尺寸相接近或可以稍大一些。（4）建立典型的立方格結(jié)構(gòu)，其三個(gè)坐標(biāo)（x,y,z）方向的網(wǎng)格尺寸一致。按先x方向，再y方向，最后z方向的順序?qū)α⒎礁襁M(jìn)行編號(hào)。（5）對(duì)當(dāng)前“接觸面對(duì)”的所有主從節(jié)點(diǎn)循環(huán)，根據(jù)節(jié)點(diǎn)坐標(biāo)確定其所在的立方格的編號(hào)基于以上全局搜索流程編制的程序框圖如圖1所示。2.2局部搜索流程本文結(jié)合點(diǎn)面算法與內(nèi)外算法完成接觸局部搜索。由全局搜索獲得“接觸面對(duì)”的點(diǎn)面關(guān)系，采用內(nèi)外算法判斷“接觸面對(duì)”主面上的節(jié)點(diǎn)落入哪些從面接觸面單元內(nèi)，為了避免搜索盲區(qū)，本文采用了相關(guān)單元節(jié)點(diǎn)法向“平均向量”。如圖2所示，在接觸單元中，節(jié)點(diǎn)B的法向量n取單元AB的法向量n下面給出一組“點(diǎn)面接觸測(cè)試對(duì)”的局部搜索步驟：（1）對(duì)測(cè)試對(duì)中的面單元數(shù)進(jìn)行循環(huán)，根據(jù)“點(diǎn)”法向量與單元面法向量，進(jìn)一步判斷是否為潛在接觸關(guān)系。（2）針對(duì)潛在接觸單元，計(jì)算“點(diǎn)”到單元面投影點(diǎn)位置及其貫入量。（3）由內(nèi)外算法圖4給出了基于局部搜索流程編制的程序框架。綜上所述，本文采用主從面算法與位碼算法相結(jié)合完成接觸全局搜索，點(diǎn)面算法與內(nèi)外算法相結(jié)合完成接觸局部搜索。同時(shí)在執(zhí)行接觸搜索前預(yù)先考慮了接觸區(qū)域分塊劃分，使得本文算法可以方便地進(jìn)行高混凝土壩和巖體高邊坡靜動(dòng)力穩(wěn)定分析的并行計(jì)算處理。3拉格朗日乘子法根據(jù)可能接觸邊界的分布特點(diǎn)，將計(jì)算域分解成不同的子區(qū)域?？紤]由N個(gè)子區(qū)域Ω式中：ρ為質(zhì)量密度；c為阻尼系數(shù)；u?不同子區(qū)域間的可能接觸邊界條件為分離狀態(tài)、黏著狀態(tài)或滑動(dòng)狀態(tài)。接觸力求解采用拉格朗日乘子法。對(duì)每一個(gè)編號(hào)為i的子區(qū)域，獨(dú)立進(jìn)行網(wǎng)格剖分和有限元離散后，在引進(jìn)拉格朗日乘子λ后，動(dòng)力學(xué)方程式（1）的空間離散形式為：式中：i為子區(qū)域編號(hào)；M為質(zhì)量矩陣；C為阻尼矩陣；K為剛度矩陣；F為荷載項(xiàng)；B采用Newmark直接積分法，將式（2）寫成增量形式：可能接觸的邊界條件的離散形式為：式中：δ為接觸面間隙，會(huì)隨著時(shí)間步變化，應(yīng)根據(jù)可能接觸邊界的幾何位置自動(dòng)判定和修正，法向接觸間隙δ將式（3）和式（4）寫成整體的形式：由式（5）消去未知量dU，整理出關(guān)于乘子力的柔度矩陣D及接觸荷載向量Q:其中：將接觸力的計(jì)算轉(zhuǎn)化為不等式方程組的求解：采用擬高斯迭代法4算法程序?qū)Ρ认旅鎸⒉捎脙蓚€(gè)算例驗(yàn)證本文提出的面面接觸算法。算例選取的一般原則是模型簡(jiǎn)單典型，將計(jì)算結(jié)果與其理論解對(duì)比，再者借助于商業(yè)軟件相應(yīng)的功能對(duì)本文算法程序進(jìn)行補(bǔ)充驗(yàn)證。第一個(gè)為塊體接觸模型，主要驗(yàn)證算法的穩(wěn)定性；第二算例為經(jīng)典的赫茲問題，將數(shù)值計(jì)算的接觸力和接觸半徑結(jié)果與理論解比較，并借助商業(yè)軟件考查接觸面上的變形分布是否與本文算法結(jié)果是否一致，以驗(yàn)證本文算法的正確性和精確性。4.1突加重力作用下的穩(wěn)定性考慮由兩個(gè)塊體組成的系統(tǒng)，如圖5所示，塊體1和塊體2在x方向、y方向和z方向的長(zhǎng)度均為100m，塊體1在3個(gè)尺寸方向均為10等分；塊體2在豎向10等分，在橫向兩個(gè)方向20等分。塊體1和塊體2的密度均取10kg/m圖6給出了兩塊體在突加重力作用下穩(wěn)定后的變形云圖。圖6顯示，最大變形位于塊體2上的上表面，其豎向位移值為-1.5×10圖7和圖8分別給出了位于兩塊體接觸面塊體2底面中點(diǎn)和上表面中點(diǎn)的位移時(shí)程曲線。圖7和圖8顯示，時(shí)間步長(zhǎng)為0.01s時(shí)的變形很快穩(wěn)定下來，其次是步長(zhǎng)為0.04s，然后為0.06s和0.10s。由于本文算法動(dòng)力方程的時(shí)間積分采用了無條件穩(wěn)定的Newmark積分，因此，時(shí)間步長(zhǎng)的大小僅會(huì)影響計(jì)算精度而不會(huì)影響其數(shù)值計(jì)算的穩(wěn)定性。由圖7和圖8可以看出，在計(jì)算穩(wěn)定后，接觸塊體2底面和上表面豎向位移值分別為-1.0×104.2算法的計(jì)算結(jié)果兩平行接觸的圓柱體接觸體在本文中，圓柱彈性模量E=210GPa，圓柱半徑R由本文算法程序得到的沿Y向的應(yīng)力分量云圖如圖10所示，最大接觸應(yīng)力發(fā)生在接觸中心點(diǎn)，其值為4136.6MPa，與赫茲理論值為4256.4MPa相比，相對(duì)誤差僅為0.9%；計(jì)算得到接觸半徑為0.188mm，赫茲理論解值為0.185mm，相對(duì)誤差為0.6%，本文算法的計(jì)算結(jié)果與理論解幾乎完全吻合。該算例結(jié)果表明，本文的接觸算法是正確的和有效的。另外，采用ANSYS軟件對(duì)兩接觸圓柱體變形進(jìn)行補(bǔ)充計(jì)算，得到如圖11所示的兩平行接觸的圓柱體接觸體豎向變形云圖。從圖11可以看出，本文算法程序與ANSYS軟件計(jì)算得到的變形分布一致。對(duì)比表1列出的沿兩圓柱上下接觸面由左到右的豎向位移，二者上接觸面位移的最大誤差為0.71%，計(jì)算結(jié)果幾乎完全一致。5局部搜索算法通常位移邊界條件的引入有3種方法：（1）合成總剛矩陣時(shí)劃去邊界位移約束對(duì)應(yīng)的各行各列元素，緊縮總剛度及荷載列陣，將約束的影響作用轉(zhuǎn)移到荷載列陣，這個(gè)方法降低了位移方程的階數(shù)。但是在面面接觸算法中，若接觸面上含有邊界約束條件，形成接觸力柔度矩陣時(shí)，容易出現(xiàn)矩陣奇異。（2）對(duì)角元素充大數(shù)法，就是把邊界位移為零的那一行對(duì)角元素充一個(gè)大數(shù)，使得非對(duì)角元素相對(duì)地較小，獲得最后求解得到的位移值趨近于零的效果。盡管這種方法較為簡(jiǎn)便，但取值不當(dāng)往往會(huì)影響計(jì)算結(jié)果的精度，有時(shí)還會(huì)使方程變成病態(tài)而得不到解；（3）“充0置1”，即將邊界位移約束的這一行的對(duì)角線元素置1，與這個(gè)對(duì)角線元素相應(yīng)行和列的其它元素都充0，位移方程右端的相應(yīng)行為約束位移值，這樣就保證了邊界行的位移等于約束值。通過本文反復(fù)的數(shù)值試驗(yàn)，認(rèn)為第三種方法的位移約束邊界處理方法最適宜本文面面接觸模型隱式接觸力的求解模式。面面接觸模型的優(yōu)點(diǎn)可用于任意形狀多體接觸面，接觸面可以具有不同的網(wǎng)格。但這樣在接觸分析過程中需要進(jìn)行接觸狀態(tài)判定，以及確定接觸點(diǎn)位置。面面接觸算法通過全局搜索粗略判斷所有可能潛在接觸測(cè)試對(duì)，再由局部搜索找出接觸投影點(diǎn)所在單元面的相對(duì)位置關(guān)系，從而確定接觸節(jié)點(diǎn)與單元面間的接觸狀態(tài)。本文全局搜索算法融合了主從面算法與位碼算法，局部搜索算法繼承了點(diǎn)面算法與內(nèi)外算法的優(yōu)點(diǎn)。本文研究在現(xiàn)有算法作了以下工作：（1）在全局搜索中構(gòu)建了描述“接觸面對(duì)”的主面單元和從面單元的共享實(shí)常數(shù)，可以將主面單元和從面單元相關(guān)聯(lián)實(shí)現(xiàn)快速搜索；（2）在局部搜索中采用了相關(guān)單元節(jié)點(diǎn)法向“平均向量”，從而解決搜索盲區(qū)問題，同時(shí)避開求解了點(diǎn)到面投影的非線性方程；（3）在接觸力求解方面，采用拉格朗日乘子法的動(dòng)力方程的增量法，以便解決材料非線性問題；（4）在接觸搜索預(yù)處理過程中已考慮了接觸