范亞杰, 李 燕, 李中潘, 陳薈鍵, 馮志強(qiáng)

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

0 引 言

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

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

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

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

1 理 論 基 礎(chǔ)

1.1 時間積分計算

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

(1)

矢量．

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

(2)

(3)

(4)

(5)

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

(6)

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

(7)

ui+1=ui+Δu．

(8)

(9)

(10)

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

(11)

(12)

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

(13)

1.2 動態(tài)接觸模型

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

(14)

(15)

(16)

intKμ和bdKμ表示Coulomb摩擦錐的內(nèi)部和邊界．

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

圖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是從局部坐標(biāo)系到全局坐標(biāo)系的轉(zhuǎn)換矩陣．若考慮初始間隙g,式(20)可以改寫為

x=HΔu+g．

(21)

同理,在局部坐標(biāo)系中,接觸力r被定義為

r=rtt+rnn．

(22)

根據(jù)虛功原理,局部接觸力和全局接觸力有以下關(guān)系:

(23)

可得到

Rc=HTr．

(24)

根據(jù)以上推導(dǎo),得到整個接觸系統(tǒng)的方程為

(25)

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

(26)

定義

(27)

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

(28)

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

(29)

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

(30)

(31)

(32)

由此可得

(33)

1.3 雙勢接觸理論

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

(34)

(35)

(36)

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

(37)

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

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

(38)

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

(39)

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

rα=projKμ(rα*)．

(40)

(41)

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

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

(42)

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

1.4 光滑變形梯度

假設(shè)任意一個質(zhì)點在初始構(gòu)型中的幾何位置為X,在當(dāng)前構(gòu)型中的幾何位置為x,該質(zhì)點的位移為u=x-X,則該質(zhì)點從初始構(gòu)型到當(dāng)前構(gòu)型的變形梯度為

(43)

圖3 四邊形單元劃分光滑域個數(shù)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來標(biāo)記．圖4展示了9個四邊形單元被劃分為36個光滑域的情況．

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

(44)

(45)

(46)

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

(47)

(48)

(49)

其中,I是單位張量．

(50)

對應(yīng)的四階材料張量可以寫為

(51)

其分量形式為

(52)

Blatz-Ko模型適用于描述大變形材料[29],其應(yīng)變能密度函數(shù)的表達(dá)式為

(53)

(54)

(55)

2 數(shù) 值 算 例

如圖5所示,超彈性體的碰撞模型由三部分組成:超彈性體小球和兩個對稱的楔形剛體．其中, 小球的圓心O點坐標(biāo)為(0 m,0 m), 半徑R=0.01 m,楔形體各頂點的坐標(biāo)為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,接觸區(qū)域摩擦因數(shù)μ為0.2．

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

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

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

圖6 3種工況下的von Mises應(yīng)力圖Fig. 6 The von Mises stress contours for 3 conditions

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

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

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

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

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

表2 計算效率對比結(jié)果

如圖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,超彈性小球與超彈性體之間的接觸區(qū)域;接觸帶2,超彈性塊與剛性板之間的接觸區(qū)域,兩個接觸區(qū)域摩擦因數(shù)均為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應(yīng)力云圖．可以看到超彈性小球砸向超彈性塊,小球和塊發(fā)生變形并且最后分離的整個過程．

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

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

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

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

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

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

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

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

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

圖15 Von Mises應(yīng)力圖(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達(dá)到最大變形位移,對應(yīng)的值為-0.004 6 m,兩超彈性球同時與超彈性板發(fā)生碰撞．圖16(c)為碰撞過程的能量-時間歷程曲線,在相同時間內(nèi),系統(tǒng)出現(xiàn)了兩次的動能和勢能值相等的現(xiàn)象．發(fā)現(xiàn)在0.002 s之前CSFEM所計算出的能量與MSC.Marc比較吻合,在0.002 s之后,因為兩個超彈性小球、超彈性板和剛性板均已分離,不再產(chǎn)生摩擦效應(yīng),所以總能量應(yīng)該是保持不變的,MSC.Marc計算出的總能量仍有較小的下降,因此本文使用的算法更加精確．

(a) 節(jié)點A的y方向位移 (b) 節(jié)點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的速度,使其產(chǎn)生非對稱碰撞．

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

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

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

圖18 Von Mises應(yīng)力圖(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最先與小球發(fā)生碰撞,碰撞點A和B分別在0.001 1 s和0.001 7 s時發(fā)生最大位移,對應(yīng)的值為0.004 6 m和0.002 4 m．兩次碰撞的速度不同,而且發(fā)生碰撞的時刻也不同,在左邊超彈性小球與超彈性塊接觸時,右邊小球就與超彈性塊發(fā)生了碰撞．因此,兩次碰撞會相互影響,節(jié)點位移曲線也因第二次碰撞的影響而不平滑,圖示結(jié)果和MSC.Marc所得的數(shù)值結(jié)果吻合良好．

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

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

3 結(jié) 論

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

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

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

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