基于Nitsche方法與擬牛頓求解的二維接觸問題等幾何分析

胡清元， 沈莞薔， 蔣芳芳

(江南大學 理學院，無錫 214122)

1 引 言

接觸問題是工業工程中常見的一類力學問題[1]。在接觸問題分析方法中，解析方法精度高但適用范圍受限。有限元方法適用范圍廣，但采用的近似網格為接觸邊界帶來嚴重離散誤差，迭代產生的單元畸變可能導致迭代難以收斂。等幾何分析[2]方法采用非均勻有理B樣條(NURBS)為基函數建模，可精確地描述結構邊界，適合分析接觸問題。

然而，當前等幾何接觸分析常見方法存在一定不足。等幾何罰函數法[3]采用類似彈簧剛度的罰系數支撐接觸面，但問題求解受系數取值影響很大。等幾何分析中廣泛采用的Mortar方法[4]通過增加虛擬公共邊界施加接觸條件，該方法將接觸面虛擬自由度增加到系統總自由度中，因此總自由度數隨著接觸面演化而不斷變化。此外，等幾何配點法[5]也用于接觸問題分析。半解析方法[6]將基本解與快速數值求解相結合，具有速度快和計算精度高的特點，在等幾何分析中的應用有待進一步發展。

近年來，基于Nitsche方法的等幾何分析受到持續關注。Nitsche方法最早用于有限元位移邊界條件的施加[7]，之后應用到域分解問題[8]和接觸問題[9]，以及后續的摩擦大變形接觸問題[10]。Apostolatos等[11]綜合對比了多種域分解方法的數值性能，結果表明Nitsche方法的綜合表現最為穩定。Hu等[12]首次采用Nitsche方法還原接觸邊界條件，該列式的優勢是不增加系統自由度，且有完備的數學理論支撐[9]。然而，標準的Nitsche列式需指定罰系數，該系數通常通過求解界面上的廣義特征值問題得到。但對于接觸問題，接觸面隨著迭代不斷變化，若需在每個迭代步前都求解廣義特征值問題，將大幅度增加計算量。

在方程求解方面，由于接觸問題的強非線性，一般采用的Newton-Raphson方法需要對列式進行線性化處理，推導切線剛度陣，并在每一步迭代中對大規模矩陣求逆。這不僅需要繁雜的推導，也加重了迭代計算量。實際上，擬牛頓方法也可用于非線性問題的求解，其核心在于采用了近似的切線剛度陣，并可采取BFGS逆更新格式以避免矩陣求逆。Matthies等[13]結合擬牛頓格式與步長搜索，計算了材料和幾何非線性問題。Laursen等[14]進一步采用擬牛頓方法求解接觸問題，為了保證收斂，當迭代不穩定時需計算切線剛度陣。Gabriel等[15]將第一次迭代求得的解作為第二次迭代的初始近似解。然而，Nitsche方法直接將所有可能產生接觸的界面納入其接觸列式，接觸面的更新會直接反映在Nitsche接觸列式當中。因此，如何在接觸面劇烈演化和擬牛頓迭代不穩定時進行簡單有效地修正，是有待進一步研究的問題。

2 等幾何接觸問題

2.1 等幾何分析簡述

等幾何分析本質上是采用NURBS為形函數的等參有限元方法，采用如下方式建模幾何場，

(1)

(2)

式中uA為控制點位移自由度。

由虛位移原理，得有限元弱形式

a(u,v)=l(v)

(3)

式中a(u,v)為雙線性項，l(v)為線性項。為了便于推導，式(3)忽略了位移邊界條件。

2.2 接觸問題與接觸條件

不失一般性，考慮兩彈性體Ω1與Ω2發生接觸，定義接觸面分別為Γ1和Γ2。假定Γ1上點x1與Γ2上點x2發生接觸，定義接觸面外法向單位矢量為

(4)

定義接觸間隙函數

g(u)=(x2-x1)·n1-[[u]]n

(5)

式中(x2-x1)·n1表示當前接觸間隙，在小變形問題中認為其保持不變，[[u]]n=(u1-u2)·n1表示兩點間的相對位移。

將接觸面記作Γ，接觸條件可寫為

(6)

式中σn為其在n1方向上的分量。接觸條件1表示接觸面不能相互穿透，接觸條件2表示接觸力方向為將接觸面向外推，接觸條件3為互補條件。由于考慮的是無摩擦接觸，故令切向應力σt=0。

3 基于Nitsche方法的接觸列式

3.1 Nitsche方法簡述

Nitsche方法的核心思想是將邊界條件以功的形式引入有限元弱形式，其列式為

a(u,v)-Nitsche =l(v)

(7)

式中Nitsche部分通常可表示為邊界積分的形式

(8)

式中f(σ)和f(v)分別為關于面力和位移的函數，具體取法取決于施加的邊界條件類型。

3.2 Nitsche接觸列式

從Nitsche統一列式[12]出發，接觸列式可寫為

(9)

式中Γ為所有可能產生接觸的界面，[b]R-為標量b在R-上的投影

(10)

此外，γ為罰系數，在文獻[16,17]的基礎上，本文采用經驗公式(11)，

(11)

式中E為楊氏模量，p為形函數階次，d為問題維度，h1和h2表示接觸面兩側的網格尺寸。

在可能的接觸面Γ上，Nitsche列式通過投影算子判斷積分點上的函數值來還原接觸條件[9]。

(1) 當接觸未發生即g>0時，接觸條件3要求σn=0，此時[σn+γg]R-=0。

(2) 當接觸發生即g=0時，接觸條件2要求σn≤0，此時[σn+γg]R-=σn。

(3) 當g<0時，違背接觸條件1，通過大數γ將接觸反力放大，從而將接觸面彈回。

Nitsche接觸列式的物理意義是利用接觸力在虛位移v上做的虛功修正有限元弱形式。其中，接觸力主要為接觸面Γ上面力σn，以及當接觸面發生穿透時接觸間隙g通過罰系數γ所得接觸反力。

4 擬牛頓迭代求解格式

4.1 擬牛頓方法與BFGS逆更新公式

接觸問題為強非線性問題，接觸問題非線性列式的Newton-Raphson迭代求解格式為

KT(ui)·Δui + 1=R(ui)

(12)

式中KT為切線剛度陣，通常需要通過對列式進行線性化處理得到；余量計算方式為

R(ui)=l(v)-a(ui,v)+

(13)

Newton-Raphson迭代格式需要進行線性化推導，更重要的是，每次迭代都需要求解切線剛度陣的逆

(14)

這嚴重增加了計算量和分析時間。

(15)

并通過BFGS公式直接更新矩陣的逆，

(16)

4.2 迭代初始化與接觸修正

(17)

(18)

其中

(19)

具體修正策略為,在一定迭代次數(本文設置為8次)后，當本次余量范數大于前一次時，認為接觸面演化較為劇烈，BFGS的自修正特性不足以及時反映接觸面變化，此時需將BFGS的迭代內核替換為式(18)所得割線剛度陣，從而保證加速收斂。

5 算 例

5.1 Taylor接觸

Taylor接觸算例[18]是接觸問題的分片實驗，該算例重點關注非協調網格能否均勻地傳遞接觸壓力，用于檢驗接觸列式的正確性。如圖1所示，兩物體楊氏模量E=109Pa，泊松比υ=0.3，上方物體受均布壓力f=108N/m2，兩模型網格劃分不同。算例解析解為

圖1 Taylor接觸算例

ux=0.03xm,uy=-0.1ym

σx x=0 Pa,σy y=-100 Pa,σx y=0 Pa

接觸迭代記錄列入表1，經過6次擬牛頓迭代后系統收斂。采用本文列式計算所得位移uy結果如圖2所示，結果表明,列式可以得到正確的位移響應。應力σy y的相對誤差如圖3所示，最大誤差為0.13%。需要說明的是，Nitsche方法采用積分的形式將邊界條件弱施加在系統中，對接觸問題，接觸力通過高斯積分點傳遞，因此無法嚴格地通過Taylor接觸實驗。但隨著網格加密，誤差減小，本文列式可以弱通過Taylor接觸實驗。

圖2 位移云圖

圖3 應力相對誤差

表1 Taylor接觸迭代記錄

5.2 Hertz接觸

Hertz接觸算例是接觸問題中的經典點接觸算例，原問題要求兩個楊氏模量不同和半徑不同的三維圓柱體接觸。本算例將原問題退化，考慮一個二維彈性圓面和剛性平臺(即彈性模量和半徑均為無窮大)的接觸問題。如圖4所示，上方圓面彈性模量E=0.02×109Pa，泊松比ν=0.1，半徑r=0.2 m，承受體力gz=-1.3×106N/m3。本算例理論解為[19]，接觸力壓力為2.2918×106Pa，半接觸面寬度為0.0454 m。

圖4 Hertz接觸算例

采用8×8網格計算該算例，接觸迭代記錄列入表2，其中在第13次迭代時采用了割線剛度陣修正，修正后的余量范數大幅度降低，加速了收斂過程。計算所得豎向位移uy云圖如圖5所示，結果表明下方剛體可以有效地支撐上方彈性體，二者之間產生一段平直的接觸面。在不同網格劃分下計算所得接觸面寬度與接觸應力如圖6所示，其中黑色實線為理論解，結果表明隨著網格加密，本文列式可以得到較為精確的解答。

圖5 位移云圖

圖6 不同網格下接觸面寬度與接觸應力

表2 Hertz接觸迭代記錄

仍采用8×8網格，使用單線程計算對比。等幾何方法計算總耗時323 s，接觸應力最大誤差 9.91%；使用同樣的網格劃分，有限元計算總耗時213 s，接觸應力最大誤差28.7%。傳統有限元基函數構造簡單，而等幾何分析中的基函數需遞推計算，因此等幾何分析耗時較長。但值得指出的是，有限元前處理時劃分網格的時間并沒有統計在內，而等幾何分析直接采用幾何模型計算。并且，在本算例中等幾何方法僅用較為粗糙的網格即可精確描述接觸邊界，迭代過程穩定，最終得到較精確的結果，而這在傳統有限元中是難以做到的。

本算例的接觸表面是光滑的，對粗糙表面的接觸問題，其接觸表面為折線狀。因為無法體現等幾何方法描述曲線邊界時的精度優勢，此時的等幾何方法作為一種高階有限元方法，其計算精度應略優于傳統有限元。

6 結 論

本文在等幾何分析框架內，基于Nitsche方法推導了二維小變形無摩擦接觸問題列式。非線性列式采用擬牛頓迭代格式求解，由BFGS公式直接更新近似切線剛度陣的逆。本文給出了一個Nitsche列式中罰系數經驗公式，提出了基于Nitsche界面耦合列式的擬牛頓迭代初始化方法，并采用割線剛度陣解決由接觸面變化導致的BFGS更新能力不足問題。

所提出的接觸分析方法具有較強的適用性，在粗糙網格下依然可以精確描述接觸邊界，列式推導簡單，計算量小。算例表明，基于Nitsche方法的接觸列式以積分形式近似地施加接觸條件，基于擬牛頓方法和割線剛度陣修正的迭代過程能夠較好收斂，最終可以得到較精確的接觸面寬度和接觸面壓力解答。