面面接觸算法研究

面面接觸算法研究

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