摩擦數(shù)值算法簡(jiǎn)介

在過去的三十年間，計(jì)算接觸力學(xué)領(lǐng)域發(fā)展了多種用于求解摩擦接觸問題的算法。在對(duì)摩擦接觸問題進(jìn)行數(shù)值模擬時(shí)，一個(gè)主要困難是摩擦力與切向滑動(dòng)量之間的本構(gòu)關(guān)系是非光滑的，本構(gòu)函數(shù)在某些點(diǎn)上不可微分，從而造成數(shù)值計(jì)算中迭代收斂困難，這個(gè)問題可以通過對(duì)摩擦本構(gòu)關(guān)系的規(guī)則化來解決。

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

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

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

數(shù)學(xué)規(guī)劃方案在摩擦的模擬中也有較為成功的應(yīng)用。此外，其他的一些途徑，例如內(nèi)點(diǎn)法、屏障法、連續(xù)近似法等，也被嘗試用于摩擦接觸分析中。

摩擦接觸分析中，需要區(qū)分滑動(dòng)或粘結(jié)兩種狀態(tài)。對(duì)于粘結(jié)狀態(tài)，通?？梢允褂靡粋€(gè)參數(shù)（罰值）εT將其引入到罰函數(shù)或者拉格朗日格式，εT相當(dāng)于切向滑動(dòng)的剛度。但是這種處理使得靜摩擦力的大小依賴εT的數(shù)值。對(duì)于過小的εT，可能導(dǎo)致物理上有問題的計(jì)算結(jié)果，即接觸界面出現(xiàn)類似粘滑效應(yīng)的行為，與實(shí)際情況不符。雖然可以采用較大的εT來解決這個(gè)問題，但是εT取值過大將導(dǎo)致方程組出現(xiàn)病態(tài)或者迭代流程失去健壯性。

8

摩擦定律與摩擦定律的正則化

 8.1 庫倫摩擦定律

兩接觸表面的切向力一旦超過了某一個(gè)極限值，兩表面不再保持粘結(jié)狀態(tài)，而是發(fā)生相對(duì)的滑動(dòng)。需要使用適當(dāng)摩擦模型來描述這種切向行為。在實(shí)際工程中，庫倫定律因其簡(jiǎn)單和適用性被廣泛采用。庫侖定律的表達(dá)式如下：

其中μ為摩擦系數(shù)，在經(jīng)典庫侖定律中被視為常數(shù)，由兩接觸面的材料種類決定；pN為法向接觸壓力；?T為相對(duì)滑動(dòng)速度；tT為切向摩擦力。

庫侖定律表明，當(dāng)摩擦力小于μpN時(shí)，無相對(duì)滑動(dòng)，處于粘結(jié)狀態(tài)；當(dāng)有相對(duì)滑動(dòng)時(shí)，摩擦力必然達(dá)到極限值μpN。

摩擦系數(shù)通常受多種因素的影響，例如表面粗糙度、相對(duì)滑動(dòng)速度?T、法向接觸壓力pN、環(huán)境溫度θ等。如果考慮這些影響，可以得到經(jīng)典庫侖定律的改進(jìn)形式：

 8.2 庫侖定律的規(guī)則化

經(jīng)典庫倫摩擦定律是不可微的，在相對(duì)滑動(dòng)速度?T=0處發(fā)生由粘接至滑動(dòng)的突變，當(dāng)?T反向時(shí)，摩擦力tT也立即發(fā)生反轉(zhuǎn)。這種突變將導(dǎo)致數(shù)值計(jì)算中迭代的收斂困難。因此提出了一些規(guī)則化（即光滑化）的模型來代替經(jīng)典庫倫摩擦模型，這些規(guī)則化模型可以表達(dá)為如下形式

式中，函數(shù)φ(?T)用來描述由粘結(jié)狀態(tài)至滑動(dòng)狀態(tài)的光滑過渡。不同的規(guī)則化模型中，φ(?T)可采用多種形式，分別對(duì)應(yīng)于不同的規(guī)則化模型。

• 平方根規(guī)則化模型

• 雙曲正切規(guī)則化模型

• 反正切規(guī)則化模型
• 分段直線模型

在以上模型中，標(biāo)量參數(shù)ε為規(guī)則化參數(shù)，當(dāng)ε→0即為經(jīng)典庫倫摩擦定律。ε不宜取值過大，否則規(guī)則化模型將與實(shí)際的接觸界面粘結(jié)-滑動(dòng)規(guī)律有較大出入。

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

當(dāng)采用規(guī)則化的模型后，分析中不存在粘結(jié)狀態(tài)，而統(tǒng)一的應(yīng)用上面的數(shù)學(xué)描述來進(jìn)行摩擦分析。這在物理上可以更真實(shí)地描述實(shí)際摩擦現(xiàn)象，因?yàn)閮蓚€(gè)接觸面上都存在一定的不平度，因此完全粘結(jié)的接觸狀態(tài)是不存在的。只要有切向摩擦力的存在，總要伴隨發(fā)生一定的相對(duì)滑動(dòng)，而規(guī)則化的摩擦模型正好可以描述此類物理現(xiàn)象。

 8.3 摩擦定律的彈塑性類比

彈塑性類比方案的關(guān)鍵思想在于將切向滑動(dòng)gT分離為彈性滑動(dòng)部分gTe和塑性滑動(dòng)部分gTp，如圖2。

圖2：切向滑移量的分離

這種分離方可視為對(duì)摩擦定律的規(guī)則化。但是它也有實(shí)際的物理含義：gTe代表接觸界面上的彈性微小位移，卸載后gTe和由gTe引起的變形將恢復(fù)為0；gTp代表不可恢復(fù)的位移，在卸載后仍然存在。

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

式中，cT為界面彈性常數(shù)。

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

式中，cT為彈性系數(shù)張量。

切向摩擦的塑性部分也有類似的本構(gòu)方程，可根據(jù)彈塑性理論的標(biāo)準(zhǔn)概念導(dǎo)出。與彈塑性理論進(jìn)行類比，可以建立屈服函數(shù)fs(tT)和相應(yīng)的屈服準(zhǔn)則，如下：

屈服函數(shù)fs(tT)可以簡(jiǎn)單的表達(dá)為：

上式等價(jià)于經(jīng)典庫倫摩擦定律。

利用最大耗散原理，得

然后可以導(dǎo)出塑性滑動(dòng)的流動(dòng)法則

此外，還可以得到加載-卸載準(zhǔn)則

根據(jù)上式可確定參數(shù)?.

9

商用有限元軟件中的接觸和摩擦算法

 9.1 ABAQUS中的接觸算法

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

ABAQUS可以自動(dòng)確定二維和三維實(shí)體單元的自由面，并利用它們建立一個(gè)表面。對(duì)于殼單元和剛性單元，必須指明是正法向的面還是負(fù)法向的面。

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

（1）從屬表面應(yīng)該是網(wǎng)格劃分的更精細(xì)的表面。

（2）如果網(wǎng)格密度相近，從屬表面應(yīng)當(dāng)由更柔軟的材料組成。

定義接觸對(duì)時(shí)，須確定兩表面之間的相對(duì)滑動(dòng)是小量還是有限值。默認(rèn)設(shè)置是有限滑動(dòng)狀態(tài)。當(dāng)應(yīng)用小滑動(dòng)公式時(shí)，從模擬開始就建立從屬表面和主控表面之間的關(guān)系，確定主控表面哪個(gè)部分與從屬表面的結(jié)點(diǎn)發(fā)生作用。這種相互關(guān)系在分析過程中保持不變。對(duì)于有限滑動(dòng)情況，兩個(gè)變形表面間的有限滑動(dòng)模擬僅能用于二維問題，變形面和剛性面的有限滑動(dòng)計(jì)算可用于三維問題。

圖3為ABAQUS接觸分析的流程。ABAQUS在增量步開始前檢查所有接觸對(duì)的狀態(tài)，判斷從屬結(jié)點(diǎn)是否進(jìn)入接觸，對(duì)于接觸的結(jié)點(diǎn)，判斷是滑動(dòng)還是粘接，然后施加相應(yīng)的約束。

圖3：ABAQUS接觸算法流程

在檢查力或力矩的平衡前，首先檢查從屬結(jié)點(diǎn)接觸條件的變化，如果檢測(cè)到當(dāng)前迭代步的任何接觸狀態(tài)改變，則標(biāo)識(shí)為嚴(yán)重不連續(xù)迭代，不進(jìn)行平衡檢驗(yàn)，而是改變接觸約束進(jìn)入第二次迭代。直至接觸狀態(tài)不再變化，才進(jìn)行平衡收斂檢查，如果收斂性不滿足，則進(jìn)行下一次迭代。

 9.2 MSC.MARC中的接觸算法

MARC軟件提供了基于直接約束的接觸算法，是解決所有接觸問題的通用方法。特別是對(duì)大面積接觸，以及事先無法預(yù)知接觸發(fā)生區(qū)域的接觸問題，程序能根據(jù)物體的運(yùn)動(dòng)約束和相互作用自動(dòng)探測(cè)接觸區(qū)域。當(dāng)發(fā)生接觸時(shí)，使用邊界條件直接約束運(yùn)動(dòng)體，兩者的運(yùn)動(dòng)約束轉(zhuǎn)化成了節(jié)點(diǎn)自由度的約束和節(jié)點(diǎn)力的約束。

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

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

• 定義接觸體。
• 探測(cè)接觸。
• 施加接觸約束。
• 模擬摩擦。
• 修改接觸約束。
• 檢查約束的變化。
• 判斷分離和穿透。

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

圖4：MARC接觸算法流程

作者簡(jiǎn)介

王朋波，清華大學(xué)力學(xué)博士，汽車結(jié)構(gòu)CAE分析專家。重慶市科協(xié)成員、《計(jì)算機(jī)輔助工程》期刊審稿人、交通運(yùn)輸部項(xiàng)目評(píng)審專家。專業(yè)領(lǐng)域?yàn)檎嚻谀途?NVH/碰撞安全性能開發(fā)與仿真計(jì)算，車體結(jié)構(gòu)優(yōu)化與輕量化，CAE分析流程自動(dòng)化等。