接触与大变形问题的光滑有限元分析

范亚杰, 李 燕, 李中潘, 陈荟键, 冯志强

(西南交通大学 力学与航空航天学院, 成都 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方向上的位移和速度随时间的变化曲线．可以看出:水平方向上的位移变化很小;竖直方向上的位移随着摩擦因数的增大,小球下降的高度越来越低,竖直方向上的速度也逐渐减小．水平方向上的速度随着摩擦因数的增大,其波动也更加剧烈,竖直方向上的速度随着摩擦因数的增大而不断减小．

(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)为碰撞过程中的能量-时间历程曲线,因为只给小球施加了法向上的速度,在切向上没有发生明显的位移,所以能量损耗很小．图示结果与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) 用橡胶作为缓冲材料,适当增加碰撞面的摩擦因数,能得到较好的抗冲击性能．