接觸與大變形問題的光滑有限元分析

范亞杰, 李 燕, 李中潘, 陳薈鍵, 馮志強

(西南交通大學 力學與航空航天學院, 成都 610031)

0 引 言

接觸和大變形問題廣泛存在于工程應用當中,例如航天器的對接[1-2]、車輛制動[3]等．碰撞作為瞬時接觸[4],屬于非光滑的動態現象,摩擦是接觸表面結構非線性的最重要來源[5]．近年來,有著良好抗震和緩沖性能的超彈性材料應用廣泛,例如在防彈衣中用來吸收槍支發射的子彈或爆炸碎片的沖擊能量,應用到隔震支座中提高支座的抗震性能[6],用在動車組的車鉤緩沖裝置中以提高裝置的緩沖性能[7-8]．還有在動力學仿真研究中,捕捉物體的運動路徑對相關領域有著重要的指導作用,比如預測物體的運動路徑[9],從而可以有效防止二次碰撞對于物體造成的損傷等情況,再結合材料的性能,探究物體接觸面的應力、變形等,可以對車輛碰撞系統、緩沖材料的性能等領域的研究提供一定的指導,非光滑接觸動力學仿真研究,在工程領域有著廣闊的應用前景．

近年來,求解接觸問題的數值算法發展迅速,例如罰函數法[10-11]和Lagrange乘子法等,雖然這些方法運用廣泛,但卻不能精確滿足接觸邊界條件和摩擦定律,而且很難選擇合適的罰因子．Lagrange乘子法雖然滿足法向方向上的接觸邊界條件,但因為引入了Lagrange乘子變量,Lagrange乘子的迭代速度取決于罰因子,罰因子會根據碰撞間隙值的變化而更新．De Saxcé和Feng等[12-13]提出了雙勢方法來求解接觸問題,該方法在不增加自由度的前提下,既能準確地滿足法向的Signorini接觸條件,又能滿足切向的Coulomb摩擦定律,使編程簡單,在處理接觸問題上具有較高的精度和效率．周洋靖等[14]將雙勢方法應用于非關聯材料的彈塑性本構．Peng等[15]在原始Uzawa算法的基礎上,提出了一種具有更高效率的算法,并結合雙勢方法應用于軟體材料的大變形仿真．

隨著有限元法的深入應用,人們發現經典的有限元法存在一定的局限性和缺陷[16],由于完全兼容的Galerkin弱形式使得線性三角形單元和四面體單元的精度較差;并且在計算過程中,對生成網格的質量要求很高,后來出現的混合有限元法、無網格法[17]等都在致力于解決這些問題．比如混合有限元法,其操作仍然局限在單元內部,無網格法雖然超越了單元的界限,但其通用性不是很強,同時無網格法的編程工作和對計算成本的消耗比有限元法高出很多．針對有限元法的不足,Liu[18]在有限元法的基礎上,結合無網格法思想,提出了光滑有限元法(S-FEM)[19],成功克服了以上難題．避免了映射過程的出現,使得域積分轉化為沿著光滑域邊界的線積分[20],計算過程中不涉及形函數的梯度或導數[21],不涉及Jacobi矩陣的求解．在一定程度上,光滑有限元法允許使用質量較差的單元,允許單元在計算過程中發生較大的變形．Li等[22-23]利用光滑有限元法求解彈性和超彈性的接觸問題．Yue等[24]將該方法應用于接觸和散射問題．

本文采用基于面元的光滑有限元法(CSFEM),在雙勢框架下求解超彈性體的接觸大變形問題．首先介紹了求解動力學問題的隱式時間積分方法．其次利用雙勢方法求解得到局部接觸力,再向全局轉化,將其作為外力代入到整體方程中求解位移,在完全Lagrange框架下,將變形梯度光滑化,建立了超彈性材料的應變能密度函數表達式．最后通過與有限元軟件MSC.Marc的數值算例對比,探究了不同摩擦因數對數值結果的影響,證明了該算法的準確性．

1 理 論 基 礎

1.1 時間積分計算

在考慮接觸力的情況下,動力學控制方程可寫為

(1)

矢量．

本文采用隱式梯形法則,即Tamma-Namburu[25]方法進行時間積分:

(2)

(3)

(4)

(5)

式中,F=Fext-Fint,u為位移矢量,參數ξ和θ的取值為0≤ξ≤1和0≤θ≤1．采用分離非線性的方法,將材料非線性和接觸非線性進行分離．通過Newton-Raphson迭代方法對動力學非線性方程組進行迭代,定義迭代步為i,當i=1時,F(i)=Ft有

(6)

式中,K為剛度矩陣．根據式(1),得到關于位移的遞歸形式:

(7)

ui+1=ui+Δu．

(8)

(9)

(10)

本文中,Newton-Raphson法求解過程的收斂精度從力和位移兩個方向同時進行控制,在第i次迭代中

(11)

(12)

只有式(11)和(12)同時滿足,才能達到精度要求,其中收斂精度設置為εF=0.001,εu=0.001．跳出迭代步后,進行速度的更新,如果滿足收斂,速度更新為

(13)

1.2 動態接觸模型

假設兩個物體上的點M和點N接觸,M和N組成一個接觸映射關系α,法向矢量為n,Moreau[26]揭示了接觸模型的動態形式,將其擴展到剛體動力學的沖擊問題當中,在動力學問題中,Signorini模型使用的是基于速度的表達式:

(14)

(15)

(16)

intKμ和bdKμ表示Coulomb摩擦錐的內部和邊界．

如圖1所示,兩物體Ω1和Ω2發生碰撞,P1和P2分別為其接觸節點．基于節點P1建立局部直角坐標系,兩點之間的距離可以表示為

圖1 碰撞模型示意圖Fig. 1 Schematic diagram of the collision model

Δu=x(P1)-x(P2)=xtt+xnn,

(17)

式中,xt為接觸位移在t方向上的投影,xn為接觸位移在n方向的投影,

xt=tTΔu,xn=nTΔu,

(18)

t={1 0}T,n={0 1}T．

(19)

則

x={xtxn}={tTΔunTΔu}={tTnT}Δu=HΔu,

(20)

H是從局部坐標系到全局坐標系的轉換矩陣．若考慮初始間隙g,式(20)可以改寫為

x=HΔu+g．

(21)

同理,在局部坐標系中,接觸力r被定義為

r=rtt+rnn．

(22)

根據虛功原理,局部接觸力和全局接觸力有以下關系:

(23)

可得到

Rc=HTr．

(24)

根據以上推導,得到整個接觸系統的方程為

(25)

將式(7)代入到式(25)中,可得到

(26)

定義

(27)

其中,W是基于接觸系統構建的Delassus算子[27]．為了在迭代過程之前降低問題的維數,將式(7)重構為

(28)

其中,Δuc為接觸點的位移矢量,Δur為其他節點的位移矢量,Δur通過以下方式被消除:

(29)

將式(29)代入到式(28),可以得到

(30)

(31)

(32)

由此可得

(33)

1.3 雙勢接觸理論

雙勢理論[12]是De Saxcé和Feng于1991年提出的．他們在Legendre定理的基礎上,通過Fenchel不等式變換,建立了一組能夠處理對偶變量的方程,用來處理與隱式標準材料有關的問題,并在此基礎上,提出了雙勢函數的概念．在接觸摩擦問題中,接觸距離xα和接觸力rα構成一組對偶變量,可以得到

(34)

(35)

(36)

式(34)的隱式形式可寫為

(37)

由式(34)可知,雙勢函數包含了接觸力與切向位移的耦合項,它們不能拆分成兩個獨立的勢函數．用增廣Lagrange方法對式(37)進行處理可得到

ρbc(-xα,rα*)-ρbc(-xα,rα)≥[(rα-ρxα)-rα]·(rα*-rα),

(38)

其中,ρ為雙勢因子,取值可為柔度矩陣對角線元素最大值的倒數．增廣Lagrange接觸力rα*可以表示為

(39)

增廣Lagrange接觸力rα*在接觸迭代中又被稱為預測接觸力,所以在得到rα*后,需要在Coulomb摩擦錐上進行投影修正:

rα=projKμ(rα*)．

(40)

(41)

圖2 Coulomb摩擦錐Fig. 2 Coulomb’s frictional cone

利用Uzawa算法求解局部接觸力,計算過程包含預測-修正兩部分:

(42)

其中,i和i+1為迭代步．Uzawa算法求解出局部接觸力后,通過式(24)將局部接觸力轉換為全局接觸力,將其代入到式(7)中即可求得Δu．

1.4 光滑變形梯度

假設任意一個質點在初始構型中的幾何位置為X,在當前構型中的幾何位置為x,該質點的位移為u=x-X,則該質點從初始構型到當前構型的變形梯度為

(43)

圖3 四邊形單元劃分光滑域個數Fig. 3 Numbers of smooth domains to be divided

在本文中,將每個單元分別劃分為1,2,3,4,8和16個光滑域,用CSFEM-1SD、CSFEM-2SD、CSFEM-3SD、CSFEM-4SD、CSFEM-8SD和CSFEM-16SD來標記．圖4展示了9個四邊形單元被劃分為36個光滑域的情況．

圖4 有限元背景網格被劃分為4個光滑域Fig. 4 The finite element background mesh divided into 4 smooth domains

(44)

(45)

(46)

因此式(43)引入Heaviside型光滑函數,可以得到光滑變形梯度為

(47)

(48)

(49)

其中,I是單位張量．

(50)

對應的四階材料張量可以寫為

(51)

其分量形式為

(52)

Blatz-Ko模型適用于描述大變形材料[29],其應變能密度函數的表達式為

(53)

(54)

(55)

2 數 值 算 例

如圖5所示,超彈性體的碰撞模型由三部分組成:超彈性體小球和兩個對稱的楔形剛體．其中, 小球的圓心O點坐標為(0 m,0 m), 半徑R=0.01 m,楔形體各頂點的坐標為A(0.005 m,0 m),B(0.015 m,0 m),C(0.015 m,0.035 m),D(0.012 m,0.035 m),給小球施加一個垂直的速度Vy=-30 m/s,小球密度ρ=700 kg/m3．用Blatz-Ko模型來描述超彈性小球,其剪切模量G=3 MPa,接觸區域摩擦因數μ為0.2．

圖5 超彈性體的碰撞模型Fig. 5 The collision model for hyperelastic bodies

表1展示了在不同摩擦因數下,超彈性小球到達最低點的不同時刻及其對應的最大應力值．圖6為超彈性小球到達最低點時的von Mises應力云圖．顯然,3種工況下的最大應力值的位置分布有著明顯的差異:工況①小球下降得更低,導致變形更大,應力值更高,然而隨著摩擦因數的增大,下降的高度越來越低,且最大應力的分布也逐漸移動到接觸表面如工況③所示．

表1 不同摩擦因數不同時刻下的3種工況

圖6 3種工況下的von Mises應力圖Fig. 6 The von Mises stress contours for 3 conditions

圖7(a)為點E的碰撞力-時間歷程曲線．CSFEM取不同的光滑域與MSC.Marc進行對比,可以發現,本文所用算法與MSC.Marc結果趨于一致,不過有較小的波動,原因是工業軟件會對數據曲線進行平滑的處理．圖7(b)—(d)展示了3種摩擦因數下的能量-時間歷程曲線．可以看出,無摩擦時,能量是完全守恒的,有摩擦時總能量減少,被摩擦效應耗散．摩擦因數從0.2增加至0.4,總能的損耗幾乎相同．事實上, 隨著摩擦因數的增大, 摩擦力也增大, 切向上的滑移也減小, 所以能量的耗散不僅取決于摩擦力, 還取決于接觸面上的切向滑移．

(a) 節點E的碰撞力 (b) μ=0時能量-時間歷程曲線 (a) Collision forces of node E (b) Energy-time histories for μ=0

圖8為超彈性小球上的O點在摩擦因數μ=0,0.2,0.4這3種情況下,x,y方向上的位移和速度隨時間的變化曲線．可以看出:水平方向上的位移變化很小;豎直方向上的位移隨著摩擦因數的增大,小球下降的高度越來越低,豎直方向上的速度也逐漸減?。椒较蛏系乃俣入S著摩擦因數的增大,其波動也更加劇烈,豎直方向上的速度隨著摩擦因數的增大而不斷減?。?/p>

(a) 節點O的x方向位移 (b) 節點O的y方向位移 (a) The x-direction displacement of node O (b) The y-direction displacement of node O

本文程序運行環境為Visual Studio 2019,編程語言為C++．表2為CSFEM劃分不同光滑域的情況下與傳統有限元和商業軟件MSC.Marc在同一臺計算機上的計算效率對比．能明顯看出,CSFEM方法結合雙勢理論的計算效率高于傳統有限元方法,略低于商業軟件,且隨著光滑域劃分的數量增多,計算時間也越長．

表2 計算效率對比結果

如圖9所示, 超彈性小球半徑R=0.025 m, 小球密度ρ為700 kg/m3, 給它施加一個向下的速度Vy=1 m/s;下面的超彈性塊長L=0.16 m,高度h=0.06 m,最下面是剛性板,用Blatz-Ko模型來描述超彈性小球和超彈性塊,其剪切模量G=3 MPa．接觸帶1,超彈性小球與超彈性體之間的接觸區域;接觸帶2,超彈性塊與剛性板之間的接觸區域,兩個接觸區域摩擦因數均為0.2．

圖9 超彈性體碰撞模型Fig. 9 The collision model for hyperelastic bodies

圖10分別展示了T=0.000 5 s, 0.001 s, 0.001 5 s和0.002 s這4個時刻下的von Mises應力云圖．可以看到超彈性小球砸向超彈性塊,小球和塊發生變形并且最后分離的整個過程．

(a) T=0.000 5 s (b) T=0.001 s

從圖11(a)可知在0.000 95 s時刻,x方向上的變形位移最大值為0.001 38 m,同時還可以觀察到,碰撞導致節點A左側的節點往右拉,右側節點往左拉,因此x方向上的位移時間歷程曲線呈現反對稱分布;從圖11(b)可知,y方向上的變形位移最大值為0.007 2 m,可以清楚地觀察到點A附近的壓痕最為明顯,y方向位移時間歷程曲線呈現對稱分布．

(a) x方向位移 (b) y方向位移 (a) The x-displacement (b) The y-displacement

由圖12(a)可知超彈性塊上的點A和點B的速度時間分布．小球一開始與超彈性塊發生碰撞時, 點A與小球最先接觸, 之后隨著縱波向其底部傳遞, 最先傳遞到節點B, 此時長方體塊還未彈起, 開始產生振動現象．圖12(b)為碰撞過程中的能量-時間歷程曲線,因為只給小球施加了法向上的速度,在切向上沒有發生明顯的位移,所以能量損耗很?。畧D示結果與MSC.Marc吻合得比較好．

(a) 節點A和B的y方向速度 (b) 能量-時間歷程曲線 (a) The y-direction velocities of nodes A and B (b) Energy-time histories

如圖13所示,基于上一個算例,系統由一個超彈性小球碰撞增加為兩個超彈性小球,并以相同初速度Vy=-20 m/s與超彈性板發生碰撞．接觸帶1,左邊超彈性小球與超彈性板接觸的區域;接觸帶2,右邊超彈性小球與超彈性板接觸的區域;接觸帶3,超彈性板與最下層的剛性板接觸的區域．

圖13 兩個小球碰撞模型圖Fig. 13 The model for 2 colliding balls

圖14為接觸帶1-2上的節點位移圖．能明顯看出碰撞節點A和B兩端的節點均被拉向兩側,碰撞點A左側區域向x負方向變形,而碰撞點B右側區域向x正方向變形,因為兩個超彈性小球速度相同,所以碰撞點A和B處區域的壓痕相同,為0.003 43 m．圖15為碰撞產生的von Mises應力圖,與圖14的結果相吻合,碰撞點A和B的應力最大,兩次碰撞產生的應力相互影響,在接觸帶1和2的中間節點發生波的匯聚,波的能量相互影響產生一部分的耗散．

(a) x方向位移 (b) y方向位移 (a) x-displacement (b) y-displacement

圖15 Von Mises應力圖(Vy0,L=-20 m/s, Vy0,R=-20 m/s)Fig. 15 The von Mises stress contour (Vy0,L=-20 m/s, Vy0,R=-20 m/s)

如圖16(a)、(b)所示,碰撞點A和B在0.001 11 s達到最大變形位移,對應的值為-0.004 6 m,兩超彈性球同時與超彈性板發生碰撞．圖16(c)為碰撞過程的能量-時間歷程曲線,在相同時間內,系統出現了兩次的動能和勢能值相等的現象．發現在0.002 s之前CSFEM所計算出的能量與MSC.Marc比較吻合,在0.002 s之后,因為兩個超彈性小球、超彈性板和剛性板均已分離,不再產生摩擦效應,所以總能量應該是保持不變的,MSC.Marc計算出的總能量仍有較小的下降,因此本文使用的算法更加精確．

(a) 節點A的y方向位移 (b) 節點B的y方向位移 (a) The y-direction displacement of node A (b) The y-direction displacement of node B

給左邊小球施加Vy=-20 m/s的速度,右邊小球施加Vy=-10 m/s的速度,使其產生非對稱碰撞．

圖17為接觸帶1-2上的節點位移圖, 隨著時間的增加, 接觸區域逐漸增大, 兩處碰撞也相互產生影響, 碰撞點A左側區域向X正方向發生變形, 碰撞點B右側區域也向X正方向變形．點A處壓痕明顯更大, 為0.004 58 m．

(a) x方向位移 (b) y方向位移 (a) x-displacement (b) y-displacement

圖18為超彈性體的von Mises應力云圖,在接觸帶1處最先產生應力,碰撞點A和B處區域應力較大．兩次碰撞產生的應力相互影響,在接觸帶1和2的中間節點發生波的匯聚,波的能量相互影響,產生一部分的耗散．接觸帶1產生的波最快到達接觸帶3,使得超彈性板左側開始脫離剛性板,并向兩側擴展;接觸帶2產生的波傳遞到接觸帶3時,超彈性板右側也逐漸脫離剛性板．

圖18 Von Mises應力圖(Vy0,L=-20 m/s, Vy0,R=-10 m/s)Fig. 18 The von Mises stress contour (Vy0,L=-20 m/s, Vy0,R=-10 m/s)

如圖19(a)所示,明顯看出點A最先與小球發生碰撞,碰撞點A和B分別在0.001 1 s和0.001 7 s時發生最大位移,對應的值為0.004 6 m和0.002 4 m．兩次碰撞的速度不同,而且發生碰撞的時刻也不同,在左邊超彈性小球與超彈性塊接觸時,右邊小球就與超彈性塊發生了碰撞．因此,兩次碰撞會相互影響,節點位移曲線也因第二次碰撞的影響而不平滑,圖示結果和MSC.Marc所得的數值結果吻合良好．

(a) 節點A的y方向位移 (b) 節點B的y方向位移 (a) The y-direction displacement of node A (b) The y-direction displacement of node B

圖19(b)同上一對稱碰撞算例一樣,后續過程中,兩個超彈性小球、超彈性板和剛性板均已分離,不再產生摩擦效應,所以總能量應該保持不變,MSC.Marc計算出的總能量仍有較小的下降,因此本文使用的算法更加準確．

3 結 論

本文發展了一種數值方法來處理超彈性體的接觸碰撞和大變形問題,該方法結合了CSFEM和計算接觸摩擦的雙勢理論．利用光滑有限元法在求解大變形問題上的優點,結合雙勢方法計算接觸力的優勢,從而可以達到較高的精度要求．最后,根據數值結果分析,得到了以下結論:

1) 在分析碰撞體的變形時,本文所用算法數值精度較高,與MSC.Marc計算出的結果基本一致,并且滿足了能量守恒或耗散定律．

2) 碰撞力的大小,能量的耗散,不僅受摩擦因數的影響,也會受到接觸界面結構的影響．

3) 用橡膠作為緩沖材料,適當增加碰撞面的摩擦因數,能得到較好的抗沖擊性能．