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

胡清元， 沈莞薔， 蔣芳芳

(江南大學(xué) 理學(xué)院，無(wú)錫 214122)

1 引 言

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

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

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

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

2 等幾何接觸問(wèn)題

2.1 等幾何分析簡(jiǎn)述

等幾何分析本質(zhì)上是采用NURBS為形函數(shù)的等參有限元方法，采用如下方式建模幾何場(chǎng)，

(1)

(2)

式中uA為控制點(diǎn)位移自由度。

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

a(u,v)=l(v)

(3)

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

2.2 接觸問(wèn)題與接觸條件

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

(4)

定義接觸間隙函數(shù)

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

(5)

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

將接觸面記作Γ，接觸條件可寫(xiě)為

(6)

式中σn為其在n1方向上的分量。接觸條件1表示接觸面不能相互穿透，接觸條件2表示接觸力方向?yàn)閷⒔佑|面向外推，接觸條件3為互補(bǔ)條件。由于考慮的是無(wú)摩擦接觸，故令切向應(yīng)力σt=0。

3 基于Nitsche方法的接觸列式

3.1 Nitsche方法簡(jiǎn)述

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

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

(7)

式中Nitsche部分通?？杀硎緸檫吔绶e分的形式

(8)

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

3.2 Nitsche接觸列式

從Nitsche統(tǒng)一列式[12]出發(fā)，接觸列式可寫(xiě)為

(9)

式中Γ為所有可能產(chǎn)生接觸的界面，[b]R-為標(biāo)量b在R-上的投影

(10)

此外，γ為罰系數(shù)，在文獻(xiàn)[16,17]的基礎(chǔ)上，本文采用經(jīng)驗(yàn)公式(11)，

(11)

式中E為楊氏模量，p為形函數(shù)階次，d為問(wèn)題維度，h1和h2表示接觸面兩側(cè)的網(wǎng)格尺寸。

在可能的接觸面Γ上，Nitsche列式通過(guò)投影算子判斷積分點(diǎn)上的函數(shù)值來(lái)還原接觸條件[9]。

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

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

(3) 當(dāng)g<0時(shí)，違背接觸條件1，通過(guò)大數(shù)γ將接觸反力放大，從而將接觸面彈回。

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

4 擬牛頓迭代求解格式

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

接觸問(wèn)題為強(qiáng)非線(xiàn)性問(wèn)題，接觸問(wèn)題非線(xiàn)性列式的Newton-Raphson迭代求解格式為

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

(12)

式中KT為切線(xiàn)剛度陣，通常需要通過(guò)對(duì)列式進(jìn)行線(xiàn)性化處理得到；余量計(jì)算方式為

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

(13)

Newton-Raphson迭代格式需要進(jìn)行線(xiàn)性化推導(dǎo)，更重要的是，每次迭代都需要求解切線(xiàn)剛度陣的逆

(14)

這嚴(yán)重增加了計(jì)算量和分析時(shí)間。

(15)

并通過(guò)BFGS公式直接更新矩陣的逆，

(16)

4.2 迭代初始化與接觸修正

(17)

(18)

其中

(19)

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

5 算 例

5.1 Taylor接觸

Taylor接觸算例[18]是接觸問(wèn)題的分片實(shí)驗(yàn)，該算例重點(diǎn)關(guān)注非協(xié)調(diào)網(wǎng)格能否均勻地傳遞接觸壓力，用于檢驗(yàn)接觸列式的正確性。如圖1所示，兩物體楊氏模量E=109Pa，泊松比υ=0.3，上方物體受均布?jí)毫=108N/m2，兩模型網(wǎng)格劃分不同。算例解析解為

圖1 Taylor接觸算例

ux=0.03xm,uy=-0.1ym

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

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

圖2 位移云圖

圖3 應(yīng)力相對(duì)誤差

表1 Taylor接觸迭代記錄

5.2 Hertz接觸

Hertz接觸算例是接觸問(wèn)題中的經(jīng)典點(diǎn)接觸算例，原問(wèn)題要求兩個(gè)楊氏模量不同和半徑不同的三維圓柱體接觸。本算例將原問(wèn)題退化，考慮一個(gè)二維彈性圓面和剛性平臺(tái)(即彈性模量和半徑均為無(wú)窮大)的接觸問(wèn)題。如圖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網(wǎng)格計(jì)算該算例，接觸迭代記錄列入表2，其中在第13次迭代時(shí)采用了割線(xiàn)剛度陣修正，修正后的余量范數(shù)大幅度降低，加速了收斂過(guò)程。計(jì)算所得豎向位移uy云圖如圖5所示，結(jié)果表明下方剛體可以有效地支撐上方彈性體，二者之間產(chǎn)生一段平直的接觸面。在不同網(wǎng)格劃分下計(jì)算所得接觸面寬度與接觸應(yīng)力如圖6所示，其中黑色實(shí)線(xiàn)為理論解，結(jié)果表明隨著網(wǎng)格加密，本文列式可以得到較為精確的解答。

圖5 位移云圖

圖6 不同網(wǎng)格下接觸面寬度與接觸應(yīng)力

表2 Hertz接觸迭代記錄

仍采用8×8網(wǎng)格，使用單線(xiàn)程計(jì)算對(duì)比。等幾何方法計(jì)算總耗時(shí)323 s，接觸應(yīng)力最大誤差 9.91%；使用同樣的網(wǎng)格劃分，有限元計(jì)算總耗時(shí)213 s，接觸應(yīng)力最大誤差28.7%。傳統(tǒng)有限元基函數(shù)構(gòu)造簡(jiǎn)單，而等幾何分析中的基函數(shù)需遞推計(jì)算，因此等幾何分析耗時(shí)較長(zhǎng)。但值得指出的是，有限元前處理時(shí)劃分網(wǎng)格的時(shí)間并沒(méi)有統(tǒng)計(jì)在內(nèi)，而等幾何分析直接采用幾何模型計(jì)算。并且，在本算例中等幾何方法僅用較為粗糙的網(wǎng)格即可精確描述接觸邊界，迭代過(guò)程穩(wěn)定，最終得到較精確的結(jié)果，而這在傳統(tǒng)有限元中是難以做到的。

本算例的接觸表面是光滑的，對(duì)粗糙表面的接觸問(wèn)題，其接觸表面為折線(xiàn)狀。因?yàn)闊o(wú)法體現(xiàn)等幾何方法描述曲線(xiàn)邊界時(shí)的精度優(yōu)勢(shì)，此時(shí)的等幾何方法作為一種高階有限元方法，其計(jì)算精度應(yīng)略?xún)?yōu)于傳統(tǒng)有限元。

6 結(jié) 論

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

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