基于交互模式的柔性体接触碰撞动力学建模方法∗

2017-04-26王检耀刘铸永洪嘉振

物理学报 2017年15期

王检耀刘铸永洪嘉振

(上海交通大学工程力学系,上海 200240)

1 引言

工程中存在大量含接触碰撞的非连续动力学过程,如航天器的交会对接[1]、高速列车受电弓接触网间的碰撞振动[2]、机械系统中铰链间隙引起的反复撞击[3]、轮轨之间的碰撞振动[4]、有偏心转子的打夯机与地面的碰撞[5]等.在这些非连续过程中,系统的拓扑发生突变,撞击力峰值大、作用时间短,会给系统造成剧烈扰动,引起系统动力学性态的突变.因此采用高效高精度的方法对柔性体接触碰撞问题进行准确的动力学仿真是尤为重要的.从物理上来看,接触双方的界面在接触过程中需满足互不嵌入的约束条件.对任意形状柔性体的接触碰撞问题,一般采用有限元方法离散,根据处理约束条件方式的不同,目前柔性体接触碰撞问题的建模方法可以归纳为两类:附加约束方法[6−9]和罚函数方法[10−13].

附加约束方法(ACM)通过引入接触约束方程,与带有Lagrange乘子的动力学方程联立形成封闭的方程,可求得运动学变量和接触力.该方法优点是接触界面互不嵌入的约束条件得到精确的满足.不足之处在于位置约束方程与动力学方程构成指数-3微分代数方程,使得数值求解变得复杂[14−16],并且从未接触到接触时,速度发生突变[15,17],数值方法上也遇到困难.另一个缺点是接触约束是单边约束,方程的维数增加且需随接触状态的变化而改变.

罚函数方法(PFM)将接触作用视为弹簧阻尼力元,当检测到接触发生时,在接触域施加与嵌入量和嵌入速率相关的接触力.与附加约束法相比较,罚函数法的优缺点正好相反.它的优点是不需要求解接触约束方程,且方程维数保持不变;不足之处在于接触约束条件只能近似地被满足.理论上,只有罚因子(碰撞刚度系数)取的足够大时,结果才可靠.但是过大的罚因子会引起系数矩阵的病态,同时降低了数值积分时间步长的临界值.因此使用罚函数法时,罚因子的选取要权衡精度和数值稳定性,往往需要通过多次试算,才能获得合适的罚因子.为了克服罚函数法依赖于罚因子的困难,增广的Lagrange方法被提出并广泛使用.增广的Lagrange方法同样把接触作为力元,通过迭代更新力元大小,直至嵌入量小于嵌入容许量.增广的Lagrange方法减小了对罚因子的依赖,但是Seifried等[18]指出,嵌入容许量的大小同样需要人为设定,过小的嵌入容许量会使得迭代次数大大增加,甚至有可能造成不收敛,而过大的嵌入容许量会使得精度不够.

针对现有解决接触碰撞问题方法的不足,本文提出基于交互模式的建模方法(IMM):将整个模型分为局部静力学模块及主体动力学模块,局部静力学模块用于求解接触力,主体动力学模块用于求解运动学变量,两个模块在计算中进行位移和力的交互.在每一个积分步内,通过主体动力学模块首先计算出局部区域的位移边界,局部静力学模块据此计算出接触力大小,再反馈到主体动力学模块进行下一积分步的计算.局部静力学模块类似附加约束法,施加约束方程使得接触局部满足约束条件;局部静力学模块计算出的接触力以力元的形式施加到主体动力学模块,因此主体动力学模块的数值积分类似于罚函数法.交互模式方法兼具附加约束法和罚函数法的优点,同时避免了二者的不足:在接触力计算上,接触局部满足约束条件保证了接触力的精度,且无需人为选取参数;在数值方面,避免了求解复杂的微分代数方程.本文利用圆柱杆撞击方板的实验,分别使用附加约束法、罚函数法和交互模式法进行数值仿真,并与实验结果进行比较,验证了交互模式法的准确性及高效性.此外,通过对滑块在有间隙滑槽内滑动的多点碰撞问题仿真,进一步验证了所提方法在处理一般性碰撞问题中的有效性.

2 柔性体的动力学方程

在柔性多体动力学中,浮动坐标系方法是最常用的建模方法.该方法将柔性体的运动分解为刚体运动和相对的弹性变形,其中刚体运动由浮动坐标系描述,相对变形由模态坐标或有限元节点坐标描述.接触碰撞是典型的非线性高瞬态过程,为了精确描述接触的作用,本文采用有限元方法对柔性体进行离散.

如图1所示,柔性体Bi用集中质量有限元离散,er是惯性基,eb是柔性体浮动基.柔性体上任一点P的绝对位置可表示为[19]

式中,ri为浮动基基点的位置坐标阵,A为eb关于er的方向余弦阵,分别为点P变形前后相对浮动基原点的位置坐标阵,为点P相对变形坐标阵,(∗)′表示(∗)在浮动基下的坐标阵.

图1 柔性体的运动学描述Fig.1.Kinematics of a f l exible body.

(1)式对时间t求导一次和两次,可得点P速度、加速度的表达式:

根据速度变分原理,柔性体Bi的动力学方程为

式中,mP是P点质量,fP是P点受到的外力在惯性基上的坐标阵,C和K分别为有限元阻尼阵和刚度阵.

将(2)和(3)式代入(5)式中得到式中mi为质量阵,ςi为与vi二次项相关的惯性力阵,分别为外力阵和内力阵.

将上述矩阵均写为分块矩阵形式,如下:

式中Mi=TTmiT,fi=TT(miT˙q˙i+ςi−fie+fiσ).

3 接触碰撞模型

有限元方法通过将接触域进行细致的空间离散化,在离散后的接触域上形成多个接触对,而接触力元或约束是施加在各个接触对上的.碰撞过程中,接触点对的形成与分离表现了接触区域的时变过程.最广泛使用的接触域离散形式为点-面接触对(node to segment).如图2所示,把撞击物体Bi称为从物体(slave body),被撞击物体Bj称为主物体(master body),从物体上的一点S与主物体单元表面上的最近投影点M构成接触对.接触对间距矢量有如下形式:

接触点间的法向距离可以写为

式中n为主物体表面上点M位置的法向矢量.当0,表示没有接触发生;当=0,表示碰撞刚刚开始;当0,则表示接触点对之间有互相穿透.假设共有m组接触点对,所有接触点对的法向穿透量可组集为

3.1 附加约束法

附加约束法中,接触力虚功项可写为[20]

式中,λ为代表接触力的Lagrange乘子,Gq=∂G/∂q为G的雅可比矩阵.

结合(11)式,系统的总虚功为

式中M=diag,f=

图2 点-面接触对Fig.2.A node-to-segment contact pair.

由于(15)式对任意的δq,δλ都成立,系统的动力学方程可写为

(16)式是指数-3形式的微分代数方程,无论采用何种数值积分方法,求解该类型方程都有以下困难:雅可比矩阵的条件数过大、变步长算法中的局部误差过大及数值稳定性等问题[14].一般在求解之前需要先对其进行处理,如降低指数[14]、通过增广法[15]或缩并法[16]将微分代数方程转变为常微分方程.

附加约束法的另一个困难是,从未接触到接触时(碰撞发生瞬时),接触局部区域的速度是不连续的.如果不考虑速度突变,会引起数值违约及接触力异常振荡[17].Dong等[15]利用连续介质力学的间断面理论建立了速度突变条件.

附加约束法中由于施加了接触约束方程,约束条件是被严格满足的,但是需要求解指数-3形式的微分代数方程,并且要考虑未碰到碰瞬时的速度突变处理,给数值求解带来困难.另一方面,当接触状态改变,即有新的点对进入接触或已有接触点对发生分离时,整个系统的系数矩阵的维数是要随着改变的,在数值计算上要进行多次的大维数矩阵稀疏化处理,会耗费大量的计算时间.

3.2 罚函数法

罚函数法中,接触力的虚功项可写为[20]

其中,εN为罚因子(碰撞刚度参数).

那么系统动力学方程为

罚函数法没有增加方程的维数,并且无需求解微分代数方程,数值求解方便,但是罚因子εN的选取大多依赖经验,往往需要多次试算才能获得可靠结果.

3.3 交互模式法

根据圣维南原理,接触作用引起的应力集中等现象只局限在接触局部区域,因此可以将接触局部的单元取出计算接触力的大小,而局部单元所占质量很小,惯性力相比于碰撞力为小量,从而可忽略局部的惯性效应,局部单元的接触可做准静态处理.根据这个思路,将整个模型分为局部静力学模块及主体动力学模块,局部静力学模块用于求解接触力,主体动力学模块用于求解运动学变量,两个模块在计算中进行位移和力的交互.在每一个积分步内,通过主体动力学模块首先计算出局部区域的位移边界,局部静力学模块据此计算出接触力大小,再反馈到主体动力学模块进行下一积分步的计算.交互模式法的计算流程如图3所示.

图3 (网刊彩色)交互模式法的计算流程Fig.3.(color online)Schematic diagram of calculation steps using interactive mode method.

接触力项的虚功写为

式中fc为接触力坐标阵.

那么主体动力学模块的方程为

求解(20)式得到整体的广义坐标q,从而可以得到局部静力学模块区域的广义坐标p=p(q).将接触视为约束,局部静力学模块的方程为

式中f(p)为局部静力学区域的内力阵.

局部静力学模块的线性迭代求解格式为

式中,KT=∂H/∂p为切线刚度阵,p∗,λ∗为迭代的初始值.

求解(22)式得到代表接触力的Lagrange乘子λ,转换为fc=fc(λ)再反馈到(20)式进行下一步的求解,由此形成了两个模块间位移和力的交互过程.

交互模式法中,局部静力学模块类似于附加约束法,通过约束方程和Lagrange乘子求解接触力;而主体动力学模块类似于罚函数法,接触力以力元形式加入方程,积分求解运动学变量.相比于附加约束法,无需求解微分代数方程;相比于罚函数法,无需人为选取罚因子.

4 数值仿真与实验结果比较

4.1 实验描述

图4为杆-板正碰撞实验的示意图.一个圆柱形钢杆从一定高度落下,到最低位置时以一定速度撞击保持静止的铝质方形板.两个物体均用细线悬挂,杆在平衡位置时,恰好接触到方板的中心,且此时杆的轴线垂直于板面.为了测量碰撞持续的时间,将两个物体通过导线连接到直流电源上,碰撞过程中电路连通将会产生电流信号.本次实验的碰撞时间约为1.7 ms,在此时间内可认为碰撞的两个物体是自由的水平运动.

碰撞实验中,碰撞力是难以直接测量的,只能通过测试速度和位移等间接物理量来反映碰撞过程.为了测量物体的速度响应,使用了Polytec GmbH公司生产的PSV-300F型激光测振仪.本次实验中,两台激光测振仪分别测试了板背面的中心点及杆后端的中心点的速度响应.杆和板的几何与材料参数列于表1中.

图4 实验原理示意图Fig.4.Schematic diagram of the impact experiment.

表1 实验几何与材料参数Table 1.Geometrical and material parameters.

4.2 结果比较

本文求解接触碰撞问题的方法均是基于有限元离散的,因此首先须考虑网格划分的影响.接触碰撞问题对网格尺寸的要求非常高,既要保证碰撞引起的高频率的弹性波在柔性体中的传播,又要准确表现接触局部范围内的高应力分布及接触域的时变历程.碰撞问题网格的划分基于以下的原则[18]:对于非接触区域,有限元网格应覆盖撞击产生的最高频率的弹性波的传播,网格尺寸需满足l≤c/20fmax,c为弹性波速,fmax为考虑的最高频率,一般取50—100 kHz.根据该原则本例中的柔性体非接触区域网格约取为5 mm;对于接触局部区域,需要更小的网格,可先通过解析方法如Hertz接触公式来粗略估计最大的接触半径a,再确定接触区域内的网格尺寸,以保证接触过程中有多个接触点对发生接触.通过计算本例中a=0.49 mm,接触域内的网格尺寸取为0.2 mm.

对于交互模式法,局部静力学处理部分的区域需要单独拿出来,该区域的尺寸同样按最大接触半径确定,即把最大接触半径范围内的节点所在单元作为局部静力学处理的区域,如图5所示.

图5 (网刊彩色)杆-板碰撞的局部静力学区域Fig.5.(color online)Local static region of the rodplate impact.

分别使用三种接触模型,对该杆-板碰撞的算例进行仿真.使用罚函数法时,需要人为选取罚因子大小.为了研究罚因子的取值对碰撞结果的影响,分别对不同的罚因子进行仿真.不同罚因子下的仿真结果如图6所示,图6(a)和图6(b)分别是测量点的速度响应和碰撞力曲线.当罚因子εN=105N/m时,由于过小的接触刚度导致碰撞时间明显过长,并且碰撞力和速度响应的峰值偏小,振动的高频特性也没有反映出来.随着罚因子值的增加,速度和碰撞力均逐渐收敛.εN=5×106和εN=107N/m时,结果已经几乎一致,此时的罚因子是合适的;当继续增大罚因子时,有可能因为过大的接触刚度发生数值不稳定.可以看到罚函数法的仿真结果非常依赖于罚因子的选取,须通过多次试算直至结果收敛.

与罚函数法相比,附加约束法和交互模式法均无需选择碰撞参数,一次仿真即可得到结果.将附加约束法和交互模式法的结果与罚函数法的收敛解结合在一起,并与实验结果比较,如图7所示.图7(a)给出了三种仿真方法与实验测得速度响应结果的比较.由于实验无法直接测量碰撞力,图7(b)给出了三种仿真方法碰撞力结果的互相比较.可以看到,三种接触模型的仿真结果均与实验结果符合.从局部放大曲线上看,交互模式法与附加约束法的结果几乎一致,而罚函数法的结果与它们仍有微小的偏离.这是由于附加约束法和交互模式法都是通过约束方程求解得到的接触力,而罚函数法通过施加弹簧刚度只能近似满足约束.

该算例的仿真程序在Matlab平台上编写,计算机配置参数为:Intel Core i7-4770,3.40 GHz,内存16 GB.三种方法的仿真的CPU时间列于表2中,可以看出三种方法求解该实验算例的计算效率有着较大差异.这是由于罚函数法中方程维数不随接触状态改变,数值求解简单,因此计算效率最高;交互模式法中,主体动力学模块的数值积分和罚函数法相同,但是每一步内都增加了局部小维数的非线性代数方程组的迭代求解,因此计算效率低于罚函数法;附加约束法中,由于矩阵维数增加且随接触状态改变,加上数值求解的复杂性,使得计算效率最低.虽然罚函数法的单次求解效率最高,但是考虑到其不可避免的试算过程,其计算代价仍然是很大的.因此,交互模式法相比于附加约束法和罚函数法在求解效率上都有优势.

图6 (网刊彩色)不同罚因子下的罚函数方法仿真结果 (a)板背面P1点的速度响应;(b)碰撞力时间历程Fig.6.(color online)Simulation results of penalty function method with different penalty parameters:(a)Velocity of the back point P1on the plate;(b)time history of contact force.

图7 (网刊彩色)三种接触模型的数值仿真与实验结果 (a)板背面P1点的速度响应;(b)碰撞力时间历程Fig.7.(color online)Comparison of simulation and experiment results:(a)Velocity of the back point P1 on the plate;(b)time history of contact force.

表2 三种接触模型的计算规模Table 2.Computation scale of the three contact methods.

总体来说,基于交互模式的建模方法综合了罚函数法和附加约束法的特点,相比于罚函数法,其优点是避免了人为选取参数及多次试算的过程;相比于附加约束法,优点是避免了指数-3形式的微分代数方程和速度突变的处理,在精度上几乎一致,并且在效率上占优.

5 多点碰撞算例

为了验证所提方法对于求解一般性碰撞问题的有效性,对滑块在有间隙的滑槽中滑动问题进行仿真.如图8所示,一圆头形滑块置于光滑的长形滑槽中,滑块和滑槽的几何尺寸标注于图中,厚度均为1 mm,有限元离散后共有3747个节点.滑块与滑槽均为PVC材质,弹性模量E=3 GPa,泊松比µ=0.3,密度ρ=1800 kg/m3.滑块的浮动坐标系建于尾部中点,浮动坐标系原点在惯性系下的初始位置为[0 mm,−1.7 mm]T,滑块轴线与滑槽成10°角.滑槽固定约束于地面,滑块以沿轴线初始速度v=20 m/s弹出,滑块滑动过程中将与滑槽发生多次多点的碰撞.

图8 (网刊彩色)多点碰撞问题示意图Fig.8.(color online)Schematic diagram of multi-point impact problem.

图9 (网刊彩色)多点碰撞问题的仿真结果 (a)碰撞力时间历程;(b)滑块P点的位置轨迹Fig.9.(color online)Results of multi-point impact problem:(a)Time history of contact force;(b)position of point P on the slider.

采用交互模式法对该过程进行动力学仿真,当检测到接触位置发生变化时,局部静力学模块处理的区域也随之改变.将交互模式法的结果与附加约束法的结果进行比较,以验证其对多点碰撞问题仿真的准确性.图9(a)为碰撞力的时间历程,图9(b)为滑块端点P的位置轨迹,在整个过程中共发生了四次不同位置的接触碰撞,交互模式法的结果与附加约束法的结果均能较好地符合,说明交互模式法在多点多次碰撞的一般性碰撞问题中有效.在同等条件的仿真效率比较上,交互模式法在接触阶段耗时552 s,附加约束法在接触阶段耗时699 s,交互模式法的效率仍然较高,接触阶段的计算时间减少了21%.

6 结论

针对现有解决柔性体接触碰撞问题方法的不足,本文提出基于交互模式的建模方法.该方法将整个模型分为局部静力学模块及主体动力学模块,局部静力学模块用以求解接触力,主体动力学模块用以求解运动学变量,两个模块在计算中进行位移和力的交互.在每一个积分步内,通过主体动力学模块首先计算出局部区域的位移边界,局部静力学模块据此计算出接触力大小,再反馈到主体动力学模块进行下一积分步的计算.

相比于罚函数法,交互模式法通过局部模块迭代计算接触力,避免人为选取碰撞力参数和多次试算的过程;相比于附加约束法,交互模式法将约束方程(代数方程)和动力学方程(微分方程)分开求解,避免求解指数-3形式微分代数方程的困难,同时也无需处理从未碰到碰瞬时速度的非连续问题.

然后利用一柔性杆撞击方板的实验,通过激光测振仪测试撞击点的响应,验证所提方法的准确性.对该杆-板碰撞问题,分别使用罚函数法、附加约束法和交互模式法进行仿真.结果显示,罚函数法非常依赖罚因子的选取,当罚因子选取合适时才能获得可靠结果.三种方法的仿真结果都能与实验结果符合较好.通过比较数值仿真结果与实验结果,发现本文所提方法是行之有效的,并且所需计算代价要小于另两种方法.此外,通过对滑块-滑槽多点碰撞问题的仿真,进一步验证了所提方法在处理一般性碰撞问题中的有效性.

[1]Klisch T 1998Multibody Syst.Dyn.2 335

[2]Ambrósio J,Pombo J,Rauter F,Pereira M 2009Multibody Dynamics:Computational Methods and Applications(Dordrecht:Springer Netherlands)p231

[3]Flores P,Koshy C,Lankarani H,Ambrósio J,Claro J C P 2011Nonlinear Dynam.65 383

[4]Lundberg O E,Nordborg A,Arteaga I L 2016J.Sound Vib.366 429

[5]Wang X H,Wang Q 2015Chin.J.Theor.Appl.Mech.47 814(in Chinese)[王晓军,王琪2015力学学报47 814]

[6]Tur M,Fuenmayor F J,Wriggers P 2009Comput.Method Appl.M.198 2860

[7]Duan Y C,Zhang D G,Hong J Z 2013Appl.Math.Mech.Engl.34 1393

[8]Chen P,Liu J Y,Hong J Z 2016Acta Mech.Sin.32 1

[9]Weyler R,Oliver J,Sain T,Cante J C 2012Comput.Method Appl.M.205 68

[10]Tian Q,Zhang Y,Chen L,Flore P 2009Comput.Struct.87 913

[11]Qian Z J,Zhang D G 2015J.Vib.Eng.28 879(in Chinese)[钱震杰,章定国2015振动工程学报28 879]

[12]Zhang J,Wang Q 2016Multibody Syst.Dyn.38 367