有限元理論基礎及Abaqus內部實現方式研究系列35：接觸求解算法

SnowWave02

瀏覽：5885 評論：21 收藏：20

（原創，轉載請注明出處）

1 概述

本系列文章研究成熟的有限元理論基礎及在商用有限元軟件的實現方式，通過

（1）基礎理論

（2）商軟操作

（3）自編程序

三者結合的方式將復雜繁瑣的結構有限元理論通過簡單直觀的方式展現出來，同時深層次的學習有限元理論和商業軟件的內部實現原理。

有限元的理論發展了幾十年已經相當成熟，商用有限元軟件同樣也是采用這些成熟的有限元理論，只是在實際應用過程中，商用CAE軟件在傳統的理論基礎上會做相應的修正以解決工程中遇到的不同問題，且各家軟件的修正方法都不一樣，每個主流商用軟件手冊中都會注明各個單元的理論采用了哪種理論公式，但都只是提一下用什么方法修正，很多沒有具體的實現公式。商用軟件對外就是一個黑盒子，除了開發人員，使用人員只能在黑盒子外猜測內部實現方式。

有限元理論基礎及Abaqus內部實現方式研究系列35：接觸求解算法的圖2

一方面我們查閱各個主流商用軟件的理論手冊并通過進行大量的資料查閱猜測內部修正方法，另一方面我們自己編程實現結構有限元求解器，通過自研求解器和商軟的結果比較來驗證我們的猜測，如同管中窺豹一般來研究的修正方法，從而猜測商用有限元軟件的內部計算方法。我們關注CAE中的結構有限元，所以主要選擇了商用結構有限元軟件中文檔相對較完備的Abaqus來研究內部實現方式，同時對某些問題也會涉及其它的Nastran/Ansys等商軟。為了理解方便有很多問題在數學上其實并不嚴謹，同時由于水平有限可能有許多的理論錯誤，歡迎交流討論，也期待有更多的合作機會。

通用結構有限元軟件iSolver介紹視頻：

有限元理論基礎及Abaqus內部實現方式研究系列35：接觸求解算法的圖4

http://www.yqgqt.org.cn/college/video/c12884

==第35篇：接觸求解算法==

雖然現在Abaqus的功能很多，但在上世紀80年代左右Ansys和Nastran如日中天的時候，Abaqus還能殺出一條血路，主要靠的就是它的非線性功能，如果說線性問題Ansys和Nastran是標準，那么非線性問題Abaqus就是標準。結構非線性主要分為三類：材料非線性、幾何非線性和邊界非線性。在前面的系列文章和視頻中，我們花了大量的篇幅介紹材料和幾何非線性，但一直沒有涉及邊界非線性。其實邊界非線性的求解的基本算法并不難，難的是碰撞前后的整個方程的變化特別大，不但是接觸力發生突變，而且接觸點的節點數目等也發生突變，在數值分析中最怕的就是一些參數的突然增加或者消失了，也就是說你的方程的形式不是唯一的，還需要加入一些額外的邏輯判斷，同時，邊界非線性往往都是和幾何大變形大轉動耦合在一起的，這些困難都造成邊界非線性極易無法收斂的問題。從做一個自主的帶三維交互的CAE軟件的角度上來說，前面說的困難還不是最難的，更難的是前處理對幾何體的拓撲的表達組織和在整個接觸過程中不同時刻點的的兩個離散網格邊界情況的精確設置和判斷，譬如像Abaqus和Ansys一樣在一個整車的幾萬個零部件模型中自動尋找接觸對，這些都需要精度較高的CAD內核才能做到，好的商用CAE都是采用商用的Parasolid、Acis庫來處理的，但對一般的自主CAE開發者來說價格是難以承受的，這些都嚴重制約了自主CAE向邊界非線性方向發展。

本文我們不討論前處理的邊界非線性的處理，僅討論后臺求解器需要做的接觸求解算法，重點還是在如何將接觸后的約束關系加入到有限元基本方程中，為后面實現帶接觸問題的有限元求解打下基礎。有接觸不一定就是邊界非線性，譬如兩個物體用Tie連接在一起，材料依然是線彈性的，應變也沒超過5%，那么我們可以認為依然是線性問題，也就不存在邊界非線性了。我們認為邊界非線性只是接觸的一種，邊界非線性屬于接觸的一種特殊情況，我們討論的接觸求解算法同樣適用于其它接觸問題。

1.1 整體求解算法

Abaqus Standard的接觸求解的整體流程如下，外層按增量法執行，內層按迭代法執行，其實依然是牛頓迭代的范疇，只不過第二步：Form and solve system of equations與只有幾何非線性的方程不同，此時需要加入接觸的方程的形成和求解。

1.2 包括接觸的有限元方程的組成和求解

無論是否存在接觸，有限元方程的建立都是實際問題的等效。

（1）在沒有接觸力時

如下圖情況，物體在體外力和面外力作用下變化。

有限元方程按照虛功原理求解，在物理上可解釋能量守恒原理，即在某一個時刻點，假定在外力作用下有個虛擬的位移，那么外力在虛擬位移下做的虛功=內部應變能的變化相同。

（2）當存在接觸時

當存在接觸時，體積域V和表面積的域包括多個空間獨立的物體。譬如下圖，兩個體的外表面S1，S2發生接觸。

1.2.1 正向接觸

接觸分法向壓力和切向摩擦力，在Abaqus中，如果是法向接觸力，在設置Constraint enforce method時就是選擇接觸算法

可以選擇Direct、Penalty、Augmented Lagrange三種。

（1）Direct

Abaqus的Direct直接法的含義將接觸中的約束關系直接加入到能量表達式中，如果接觸較硬（譬如硬接觸或者接觸力隨深度急速增加的問題）采用約束乘以Lagrange因子方法插入的方式，如下所示：

此方程是關于u和p的迭代量的一個線性方程組，看起來是兩個方程求解兩個未知數，但在這邊K22=0，此時無法直接縮聚掉接觸力項，給線性方程組求解增加難度，常用算法中迭代法需要主元非0，高斯消元法對角占優的矩陣更容易求，這也是同等維度下有限元相比邊界元更容易求解方程的原因，因為有限元得到的線性矩陣都是對角線附近的非零元素，而邊界元是滿陣，如果主元=0，必須經過多輪行列互換將主元=0的放最下方。

而罰函數法則無以上問題。罰函數法就和我們平常生活中的懲罰的含義是一樣的，小孩考了100分，獎勵，考了80分，什么都不干，考了70分，罰打屁股10下，考了50分不及格，罰打屁股100下，總之，考的越差，懲罰的越多。

（3）Augmented Lagrange

當然Penalty也有缺點，就是上述的K的取值是人為的，K不同得到的結果不同，而對每個模型K的取值還與接觸面的剛度相關，Augmented Lagrange其實是上面兩者的結合，顯然，這種方法計算效率相對更差，在此不詳述了。

1.2.2 切向接觸

如果是切向摩擦力，Abaqus可選擇6種，其中也包括Lagrange乘子和Penalty方法。至于具體的算法介紹和適用范圍可看USim大神的總結：

http://www.yqgqt.org.cn/content/post/1786826

這兩個方向的接觸力中，共同點是都可以選擇Lagrange法和Penalty法（罰剛度法），這也是接觸分析中最常用的兩種方法。下面我們將以一個簡單的接觸分析的算例來說明理論和Abaqus實際中的應用。

1.3 接觸分析的算例

1.3.1 Abaqus模型介紹

Abaqus中創建一個桿和一個剛性板接觸的例子。

桿沿x軸方向，左端簡支，剛性板與x軸成45度，固定在空間中，開始時桿的右端正好和剛性板接觸，在桿上加一個往右的拉力，由于桿的右端無法穿過剛性板，只能往上沿著剛性面移動，且移動的位移u1和u2必然相等。也就是約束關系是：

參數如下：

>尺寸：Truss 長度L0=1，截面積A=1。

>材料：Young’s Modulus Em=10， Poisson Ratio 0。

>邊界：左側節點固支。

>載荷：右側節點加集中力F=1。

>接觸：設置桿靠近剛性板的一點和剛性板發生接觸

>網格：劃分為一個Truss單元。

>分析步：靜力分析，打開NLGeom幾何非線性。

1.3.2 Lagrange法

1.3.2.1 理論解

關于桿的幾何非線性的應變和非線性方程的求解的詳細推導可以參考下方的視頻深入淺出有限元2-非線性：基礎理論->Abaqus操作->編程實現

1.3.2.2 仿真結果

Abaqus中設置接觸，切向無摩擦，法向采用Direct算法：

分析后得到：

可發現U1=U2=1.024e-1，和理論沒有任何差異。

1.3.3 Penalty法

1.3.3.1 理論解

這個方程組的解和K相關。

當K=10時，可得：

u1=0.1036

u2=0.0956

有興趣的可以把這兩個值帶入上面的方程可發現方程成立，可以看出u1和u2不是絕對相等，明顯差異比較大，也就是說Penalty方法得到的解只是近似解，當然有限元中只要K足夠大，那么這個誤差依然可以接受。

1.3.3.2 仿真結果

Abaqus中設置接觸，切向無摩擦，法向采用Penalty算法，且設置stiffness=10：

分析后得到的位移如下：

U1=1.047e-1，U2=8.968e-2，明顯和Penalty方法理論值有較大差異，要么Abaqus的Stiffness的含義和理論的罰剛度不同，要么就是接觸分析的迭代導致了和理論的差異，哪位大神要了解，歡迎討論。

如果Stiffness=100，分析后如下，顯然U1、U2更加接近0.1024這個實際問題的理論解。

1.4 視頻講解和操作驗證演示

如果覺得上面的文字太復雜，也可以看一下視頻的簡要講解和軟件演示，地址如下：

http://www.yqgqt.org.cn/college/video/c12884 20理論系列文章35-接觸求解算法理論.mp4

1.5 總結

本文從接觸分析時如何將接觸的約束加入到有限元基本方程出發介紹了商軟中最常用也是最基本的兩個算法Lagrange和Penalty方法，可以看到Penalty方法中stiffness剛度的選取不同對最終的求解具有較大的影響，我們想要對標Abaqus實現Abaqus完全一致的接觸分析結果，但苦于找不到Abaqus個罰剛度的選取規則，如果有哪位大神對這個問題比較了解，還請不吝賜教。