在過去的三十年間，計算接觸力學領域發展了多種用于求解摩擦接觸問題的算法。在對摩擦接觸問題進行數值模擬時，一個主要困難是摩擦力與切向滑動量之間的本構關系是非光滑的，本構函數在某些點上不可微分，從而造成數值計算中迭代收斂困難，這個問題可以通過對摩擦本構關系的規則化來解決。

目前已有多種迭代方案用于帶摩擦的接觸分析，可大致分為以下幾類：試探-校核算法、基于塑性理論中的彈塑性類比方法、基于優化理論的數學規劃方法。后兩類方法以嚴密的數學理論為基礎，其可靠性高于試探-校核算法。

試探-校核算法通常是應用于小變形情況。對于摩擦接觸問題，解的唯一性和存在性均不能保證，因此試探-校核算法在很多情況下并不可靠。但試探-校核算法仍然得到了成功的應用，在顯式有限元分析中能夠獲得合理的結果。

近幾十年，Coulomb定律和其他摩擦本構關系被納入塑性理論的研究范疇，基于彈塑性理論的返回映射方案已成功應用于有限元摩擦接觸分析，該方案使算法的收斂行為和可靠性產生了本質的提高。返回映射方案最初用于材料非線性問題，用以積分彈塑性本構關系。將摩擦定律類比為彈塑性本構關系，就可以將返回映射方案應用于帶摩擦的接觸分析。由于摩擦本構關系的非關聯性，返回映射方案所得到的切線剛度矩陣通常不對稱，增加了數值求解的難度。除返回映射方案之外，其他幾種來源于塑性理論的方法，例如屈服極限拉氏乘子法，也已用于摩擦接觸問題的數值分析中。

數學規劃方案在摩擦的模擬中也有較為成功的應用。此外，其他的一些途徑，例如內點法、屏障法、連續近似法等，也被嘗試用于摩擦接觸分析中。

摩擦接觸分析中，需要區分滑動或粘結兩種狀態。對于粘結狀態，通常可以使用一個參數（罰值）ε_T將其引入到罰函數或者拉格朗日格式，ε_T相當于切向滑動的剛度。但是這種處理使得靜摩擦力的大小依賴ε_T的數值。對于過小的ε_T，可能導致物理上有問題的計算結果，即接觸界面出現類似粘滑效應的行為，與實際情況不符。雖然可以采用較大的ε_T來解決這個問題，但是ε_T取值過大將導致方程組出現病態或者迭代流程失去健壯性。

兩接觸表面的切向力一旦超過了某一個極限值，兩表面不再保持粘結狀態，而是發生相對的滑動。需要使用適當摩擦模型來描述這種切向行為。在實際工程中，庫倫定律因其簡單和適用性被廣泛采用。庫侖定律的表達式如下：

其中μ為摩擦系數，在經典庫侖定律中被視為常數，由兩接觸面的材料種類決定；p_N為法向接觸壓力；?_T為相對滑動速度；t_T為切向摩擦力。

庫侖定律表明，當摩擦力小于μp_N時，無相對滑動，處于粘結狀態；當有相對滑動時，摩擦力必然達到極限值μp_N。

摩擦系數通常受多種因素的影響，例如表面粗糙度、相對滑動速度?_T、法向接觸壓力p_N、環境溫度θ等。如果考慮這些影響，可以得到經典庫侖定律的改進形式：

經典庫倫摩擦定律是不可微的，在相對滑動速度?_T=0處發生由粘接至滑動的突變，當?_T反向時，摩擦力t_T也立即發生反轉。這種突變將導致數值計算中迭代的收斂困難。因此提出了一些規則化（即光滑化）的模型來代替經典庫倫摩擦模型，這些規則化模型可以表達為如下形式

式中，函數φ(?_T)用來描述由粘結狀態至滑動狀態的光滑過渡。不同的規則化模型中，φ(?_T)可采用多種形式，分別對應于不同的規則化模型。

• 平方根規則化模型

在以上模型中，標量參數ε為規則化參數，當ε→0即為經典庫倫摩擦定律。ε不宜取值過大，否則規則化模型將與實際的接觸界面粘結-滑動規律有較大出入。

上述規則化模型的前三種能夠提供光滑的摩擦力-滑動量曲線，但第四種模型類似于理想彈塑性應力-應變關系，便于采用彈塑性力學中的相關方案處理。

當采用規則化的模型后，分析中不存在粘結狀態，而統一的應用上面的數學描述來進行摩擦分析。這在物理上可以更真實地描述實際摩擦現象，因為兩個接觸面上都存在一定的不平度，因此完全粘結的接觸狀態是不存在的。只要有切向摩擦力的存在，總要伴隨發生一定的相對滑動，而規則化的摩擦模型正好可以描述此類物理現象。

彈塑性類比方案的關鍵思想在于將切向滑動g_T分離為彈性滑動部分g_T^e和塑性滑動部分g_T^p，如圖2。

圖2：切向滑移量的分離

這種分離方可視為對摩擦定律的規則化。但是它也有實際的物理含義：g_T^e代表接觸界面上的彈性微小位移，卸載后g_T^e和由g_T^e引起的變形將恢復為0；g_T^p代表不可恢復的位移，在卸載后仍然存在。

對于切向摩擦的彈性部分，可采用以下各向同性線彈性模型，

式中，c_T為界面彈性常數。

也可采用各向異性線彈性模型

式中，c_T為彈性系數張量。

切向摩擦的塑性部分也有類似的本構方程，可根據彈塑性理論的標準概念導出。與彈塑性理論進行類比，可以建立屈服函數f_s(t_T)和相應的屈服準則，如下：

屈服函數f_s(t_T)可以簡單的表達為：

上式等價于經典庫倫摩擦定律。

利用最大耗散原理，得

然后可以導出塑性滑動的流動法則

此外，還可以得到加載-卸載準則

根據上式可確定參數?.

ABAQUS缺省使用拉氏乘子法施加接觸約束，用戶也可選擇使用罰函數或增廣拉格朗日法。在ABAQUS接觸模擬中，要在各個構件上建立表面，定義可能會發生接觸的一對表面（接觸對），然后按照本構關系求出表面間的法向和切向相互作用。

ABAQUS可以自動確定二維和三維實體單元的自由面，并利用它們建立一個表面。對于殼單元和剛性單元，必須指明是正法向的面還是負法向的面。

ABAQUS應用單純的主控-從屬接觸算法，執行單邊接觸檢查：從屬表面的節點不能穿透主控表面的任何部分，這種算法對于主控表面沒有限制，它可以在從屬表面的節點之間穿透從屬表面。因此必須正確選擇從屬與主控表面，要遵循的簡單原則是：

（1）從屬表面應該是網格劃分的更精細的表面。

（2）如果網格密度相近，從屬表面應當由更柔軟的材料組成。

定義接觸對時，須確定兩表面之間的相對滑動是小量還是有限值。默認設置是有限滑動狀態。當應用小滑動公式時，從模擬開始就建立從屬表面和主控表面之間的關系，確定主控表面哪個部分與從屬表面的結點發生作用。這種相互關系在分析過程中保持不變。對于有限滑動情況，兩個變形表面間的有限滑動模擬僅能用于二維問題，變形面和剛性面的有限滑動計算可用于三維問題。

圖3為ABAQUS接觸分析的流程。ABAQUS在增量步開始前檢查所有接觸對的狀態，判斷從屬結點是否進入接觸，對于接觸的結點，判斷是滑動還是粘接，然后施加相應的約束。

圖3：ABAQUS接觸算法流程

在檢查力或力矩的平衡前，首先檢查從屬結點接觸條件的變化，如果檢測到當前迭代步的任何接觸狀態改變，則標識為嚴重不連續迭代，不進行平衡檢驗，而是改變接觸約束進入第二次迭代。直至接觸狀態不再變化，才進行平衡收斂檢查，如果收斂性不滿足，則進行下一次迭代。

MARC軟件提供了基于直接約束的接觸算法，是解決所有接觸問題的通用方法。特別是對大面積接觸，以及事先無法預知接觸發生區域的接觸問題，程序能根據物體的運動約束和相互作用自動探測接觸區域。當發生接觸時，使用邊界條件直接約束運動體，兩者的運動約束轉化成了節點自由度的約束和節點力的約束。

MARC中的接觸檢查默認采用雙邊檢查，為提高檢測速度，用戶也可選擇單邊檢查。MARC不僅可以處理變形面和剛性面之間的有限滑動，對于變形面與變形面間的有限滑動也能很好的處理。

MARC軟件接觸算法基本流程如圖4所示。從圖4可見，整個過程包括：

• 定義接觸體。

• 探測接觸。

• 施加接觸約束。

• 模擬摩擦。

• 修改接觸約束。

• 檢查約束的變化。

• 判斷分離和穿透。

此外，對于接觸區域事先已知的情況，用戶還可使用間隙單元(gap單元)方案，或者基于罰函數方法，通過用戶子程序USPRNG對接觸界面施加非線性彈簧。

圖4：MARC接觸算法流程

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