接觸問題求解的效率和穩定性，不僅依賴于對接觸界面的合理離散化，更依賴于算法的質量。目前已經發展了一系列的接觸算法，分別具有不同的適用范圍，對不同的問題各有優缺點。其中以罰函數法、拉氏乘子法和直接約束法最為通用，已經實現于多種商用有限元分析軟件中。

拉氏乘子技術是將數學約束條件引入一個系統的最完美的數學描述。使用拉氏乘子引入施加接觸體必須滿足的非穿透約束條件，用以描述約束極值問題，可以使界面約束條件得到充分的滿足，確保法向不可貫入性。其不足之處是增加了方程的自由度數，并使系統矩陣主對角線元素為零。這就需要在數值方案中處理非正定系統，數學上將發生困難，需實施額外的操作才能保證計算精度，從而使計算費用增加。另外，由于拉格朗日乘子與質量無關，使得包含慣性項的接觸問題的有限元方程和顯式時間積分的求解格式不協調，導致這種由拉格朗日乘子描述的接觸算法難以直接用于顯式動力撞擊問題分析。

拉格朗日乘子技術還經常用于采用特殊的界面單元（GAP單元）描述接觸的問題分析。該方法需要預先知道接觸發生的確切部位，以便施加界面單元。這樣的額外要求對于撞擊、壓力加工等通常事先并不知道準確接觸區域所在的問題是難于滿足的。

罰函數法是另一種施加接觸約束的數值方法。其原理是一旦接觸區域發生穿透，罰函數便夸大這種誤差的影響，從而使系統的求解(滿足力的平衡和位移的協調)無法正常實現。換言之，只有在約束條件滿足之后，才能求解出有實際物理意義的結果。

用罰函數法施加接觸約束的方法可以類比成在物體之間施加非線性彈簧所起的作用。該方法不增加未知量數從而增加系統矩陣帶寬，數值上實施比較容易，在顯示動力分析中被廣泛應用。罰函數的缺點是，約束條件只能被近似的滿足，界面上有少量的相互貫入產生，由于接觸壓力是假定與貫入量成正比的，因此算得的接觸壓力分布通常是振蕩的。理論上講，增大罰參數的數值可令計算精度提高，但在實際分析中罰函數的大小受到嚴格的限制，過大的罰參數將對系統的數值求解造成不良影響。一是罰參數愈大，顯式解法時間步長的臨界值降低得愈多，二是過大的罰參數可能使接觸體的相對運動發生虛假的反向，從而使求解過程不穩定。

罰函數法的數值實現相對簡潔，因此得到了最廣泛的應用。如果罰參數的取值合理，多數情況下罰函數法能夠提供健壯穩定的求解。但是對于由變形驅動的接觸過程，拉氏乘子法通常能給出更精確更穩定的數值結果。

直接約束法也是常用的一種接觸算法。用直接約束法處理接觸問題是追蹤物體的運動軌跡，一旦探測出發生接觸，便將接觸所需的運動約束（即法向無相對運動，切線可滑動）和節點力（法向壓力和切向摩擦力）作為邊界條件直接施加在產生接觸的節點上，兩接觸體的運動約束轉化成了節點自由度的約束和節點力的約束。這種方法對接觸的描述精度高，具有普遍適應性。不需要增加特殊的界面單元，也不涉及復雜的接觸條件變化。該方法不增加系統自由度數（由于接觸關系的變化會增加系統矩陣帶寬）。但是，界面約束的有效引入強烈依賴于界面上的離散情況。對于變形體與剛性體接觸的情況，直接約束法效果明顯，但是對于兩個變形體接觸的情況，如果兩者在界面上的網格不匹配，該方法無法直接使用，需要進行特殊處理。

攝動拉格朗日法，是罰函數法和拉氏乘子法的混合模式。該方法使用的壓力分布插值函數通常比位移分布插值函數低一階。該方法能夠克服純粹的拉氏乘子法所導致的一系列數值困難。

增廣拉格朗日法可視為罰函數法與拉氏乘子法之間的折中方案，該方法綜合了兩者的優點。該方法可用于分析無摩擦接觸問題，也可用于大變形摩擦接觸問題。增廣拉格朗日技術可以與Uszawa算法結合在一起使用，在計算流程中使用嵌套的雙重循環，內循環用以處理接觸約束條件，外循環用以更新拉格朗日乘子。嵌套迭代方式增加了總的迭代次數，但是使算法的數值實現變得簡潔。

因為接觸條件是不等式約束，我們可以導出接觸問題的等效變分不等式形式，位移場的解u必須滿足該不等式。有限變形接觸問題的弱形式如下：

式中，集合K定義為

如果僅考慮超彈性或者線彈性材料本構關系，對于無摩擦接觸問題，上面的不等式與下面的最優化問題等價

為得到變分方程，需要用到法向間隙g_N和切向滑動量g_T的變分，我們將物體A表面采用參數坐標(ξ₁，ξ₂)描述，可以得到三維情況下法向間隙g_N的變分

計算δg_N時，不僅需要考慮坐標x^A和x^B的變分，尚需考慮x^A點處外法線方向?^A對于曲面坐標(ξ1，ξ2)的變分，如下式：

因為?^A是法向單位矢量，所以必然有

因此可以得到以下結果

如果物體B^A是剛性的，法向間隙的變分可以簡化為

切向滑動g_T的變分也可按類似的方式導出

曲面坐標ξ_α的變分可由下式給出

其中

如果接觸界面已知，即在增量形式的迭代求解中已經探測出接觸約束的有效集，(10)式可以寫為等式形式

式中C_c代表了接觸約束活性集的貢獻。

對于超彈性或者線彈性材料，上面的方程等價于對兩個接觸體總勢能的最小化，

式中，附加的泛函∏_c包含了接觸約束的貢獻。

在已經判斷出具體接觸區域的前提下，有數種技術可用于導出C_c和∏_c的格式，包括拉氏乘子法、罰函數法、直接約束法、攝動拉格朗日法以及增廣拉格朗日法等。

拉氏乘子法

應用拉格朗日乘子是將數學約束引入弱形式的經典方案。使用拉氏乘子技術，可以得到∏_c的列式如下：

式中，λ_N和λ_T分別為法向和切向拉氏乘子。g_N和g_T分別為法向間隙和切向的滑動量。

對∏_c變分，可以得到C_c的列式如下

對于兩接觸表面發生切向滑動的情況，切向面力矢量t_T可以根據摩擦滑動本構關系確定，因此應當改寫λ_T?δg_T→t_T?δg_T，從而得下式

罰函數法

在罰函數法中，通過罰函數項將接觸約束引入能量泛函∏

此處ε_N和ε_T分別代表法向和切向罰參數。

對上式取變分，可得

類似于拉氏乘子法，需要區分接觸界面的純粘接和切向滑動兩種狀態。對于滑動狀態，有

直接約束法

該方法將接觸區域的約束作為位移邊界條件直接引入，因此不會導致未知數數量的增加。在判斷出具體接觸區域后，對于進入接觸結點，不等式約束轉換為等式約束

根據上式可得

將上式作為位移約束條件直接引入系統方程組，即可對問題進行求解。

攝動拉格朗日法

該方案是將拉氏乘子法和罰函數法相結合的一種混合格式。攝動拉格朗日法將接觸的貢獻∏^PL表達為

積分號中的第二項和第四項是對拉氏乘子的規則化，可以看作由拉氏乘子引起的補充能量。

對上式取變分可得

右邊的第一個和第三個積分項與拉氏乘子格式相關。第二項和第四項則給出接觸界面的法向和切向的本構關系：

如果將上式代入式（31）積分號內的第一和第三項，即得到標準的罰函數格式；如果令ε_N和ε_T趨向于無窮大，則可得標準的拉氏乘子格式。

增廣拉格朗日法

增廣拉格朗日格式提供了另外一種將不可微的接觸和摩擦條件規則化的方案。其主要思路是將罰函數法或界面的本構定律與C₁可微的鞍點泛函相結合。在增廣拉格朗日法中，用以下泛函引入法向接觸的貢獻：

其中，?_N=λ_N+ε_Ng_N。

上述泛函的結構形式令其不僅適用于?_N≤0的情況，也適用于?_N>0的情況，后者意味著兩接觸體間的間隙是張開的。對上式取變分，得：

用同樣的方法，對于遵循經典庫侖摩擦定律的切向條件，可以導出類似的增廣拉格朗日格式。采用切向相對滑移增量Δg_T和切向拉氏乘子?_T=λ_T+ε_TΔg_T，構造如下的泛函引入切向約束

其中μ為摩擦系數，?_N(p頂部黑點代表箭頭)為增廣的法向接觸壓力。對于未接觸狀態(?_N>0)，可采用如下的切向泛函：

在增廣拉格朗日法中，拉氏乘子λ是未知的，因此需要在算法流程中采用迭代循環以更新λ的值。

接觸搜索是計算接觸力學中最關鍵的問題之一。一般說來，接觸搜索可分為兩個階段：首先是空間搜索，探測出可能進入接觸的物體或有限單元；然后是接觸探測，確定相互接觸的單元對或者結點對。下面我們對空間搜索和接觸探測的算法作一簡單綜述。

空間搜索實際上是采用一些排序算法對固體表面的單元進行排序操作。計算接觸力學中發展了以下幾種方案用以執行該操作。

網格元法

該方法將一包含全部接觸體的空間劃分為多個形狀規則的網格元。當單元均勻的分布于各網格元時，網格元法效率很高。但如果單元主要聚集于少數網格元中，該算法不具有優勢。因此出現了一種自適應的網格元法，可根據有限單元空間分布的非均勻程度自動調整網格元的劃分。然而，生成自適應網格將導致額外的計算花費。

二叉樹法

二叉樹法同樣是基于一個包含多個規則網格元的空間格柵。但是二叉樹結構中僅保留含有有限單元的網格元。搜索消耗的時間取決于格柵的初始構造。利用一些特殊技術，例如平衡樹的分支以及最小化樹的深度，可以使搜索時間進一步減少。生成二叉樹的時間消耗消耗為O(NlogN)，搜索的時間消耗也為這個量級。

空間堆排序算法

在該算法中，按照沿空間軸的坐標遞增的順序對一個單元鏈表排序。該算法通常與基于接觸體的網格元相結合使用，其搜索時間為O(NlogN)量級。空間堆排序算法的優點在于：不需要特殊的數據結構，對于對象的空間尺度不敏感，而且僅需要一個O(N)量級的數組來存儲必要的數據。

在構建接觸探測算法流程時，需要區分幾種不同的情況：剛性體與剛性體接觸、變形體與剛性體接觸、變形體與變形體接觸。

1. 變形體與剛性體接觸

在變形體與剛性體接觸的情況下，剛性體表面可用一個隱式函數表示

剛性體表面也可用于定義接觸法方向，因此可將其作為主接觸體或參考體。由利用隱函數f給出剛性體表面的解析描述后，判斷一個從結點x_s=(x_s,y_s,z_s)是在剛性體的內部還是外部僅需要考慮f(x_s,y_s,z_s)的數值。如果f<0，x_s在剛性體內部，如果f>0則在外部。剛性體表面法方向也很容易得到：

2. 變形體與變形體接觸

探測變形體與變形體的接觸，遠較探測變形體與剛性體的接觸困難。目前已有數種有效技術用于變形體之間的接觸探測。

Benson和Hallquis發展了一種可用于二維和三維問題的探測方法。該方法把一指定從結點x_s的局部接觸探測分為三步：判斷離x_s距離最近的主結點x_k^A、確定包含x_s在主接觸面的投影?^A(兩點代表橫杠，下同)的單元面和計算?^A點的空間流動坐標。

Kane等則提出了另外一種方法，對物體進行三角剖分，可對有限元離散模型進行貫入探測。其思路是基于以下事實：當有相互貫入發生時必然出現邊界段的交叉。在二維情況下，可通過計算邊界段的面積來確定是否發生交叉。在三維情況下，必須考慮邊界面的交叉，因此需要進行體積檢查。

來源：模態空間

王朋波，清華大學力學博士，汽車結構CAE分析專家。重慶市科協成員、《計算機輔助工程》期刊審稿人、交通運輸部項目評審專家。專業領域為整車疲勞耐久/NVH/碰撞安全性能開發與仿真計算，車體結構優化與輕量化，CAE分析流程自動化等。王朋波私人微信：poplewang；加微請注明：單位+姓名。