考虑摩擦的柔性多体系统斜碰撞理论与实验研究

陈 鹏， 刘锦阳， 洪嘉振

(上海交通大学 船舶海洋与建筑工程学院， 上海 200240)

在工程应用中，接触碰撞是一种普遍存在的动力学行为，大到航天器之间的交会对接，小到机械系统零部件之间的间隙摩擦，接触碰撞对系统的整体动力学响应以及局部构件的可靠性都有着很重要的影响。

多年来，对接触碰撞问题的研究一直是工程与力学领域中的一个重点课题。不同于其他力学行为，接触碰撞是一种持续时间短，响应频率高并且很难用力学模型描述的动力学过程。通过理论与实验的手段，学者们对接触碰撞问题进行了深入的研究并取得了很多成果。

在理论研究方面，对接触碰撞问题的力学描述可以分为刚性碰撞和弹性碰撞两大类，对于刚性碰撞，认为碰撞时间极短，不考虑碰撞过程，常用的方法为基于恢复系数假定[1-3]的动量冲量法；而对于弹性碰撞需要考虑碰撞的持续时间，常用的方法有基于Hertz碰撞理论[4]的连续碰撞力模型，以及基于有限元理论的罚函数和拉格朗日乘子法等。在工程中，对于只关心接触碰撞对系统整体大范围运动产生影响的动力学问题，用刚性接触的处理方法就能高效的解决问题。而在很多情况下，人们需要了解在特定工况下局部碰撞力的大小以及持续时间，因此只能使用弹性碰撞的方法来建模求解。对于几何形状简单的多刚体系统以及部分对局部碰撞仿真精度要求不高的柔性多体系统，连续碰撞力模型具有很高的求解效率，而对于一般的柔性体碰撞，在需要考虑局部接触面的动态变化以及局部接触区域内应力应变响应时，则只能采用有限元方法对柔性体进行空间离散，并使用罚函数方法和拉格朗日乘子法等对接触碰撞进行描述。对于摩擦的描述，最经典且使用最多的是库伦摩擦模型以及学者们针对不同的对象进行改进的修正库伦模型。

多年来国内外对接触碰撞问题的理论研究取得了很多的进展。Schiehlen[5]等在非线性有限元的基础上用罚函数方法对柔性多体系统接触碰撞问题进行研究，深入探讨了碰撞导致的弹塑性效应以及多次碰撞问题。用多尺度研究方法，Seifried等[6]通过对碰撞过程的精细分析得到恢复系数，并将恢复系数用于多体系统大尺度的动力学分析，建立了能高效求解全局碰撞问题的动力学模型。韩石磊等[7]为了缩减碰撞物体的自由度，采用了分区域的多变量方法对柔性多体系统接触碰撞问题进行建模，并用显式算法和隐式算法分别进行求解。Chen等[8]引入固定界面子结构方法，将碰撞物体分为碰撞区域和非碰撞区域，将碰撞区域的节点和铰节点同时作为边界节点，而对内部节点自由度进行缩减，建立了求解柔性多体系统接触碰撞问题高效的动力学模型。

在实验研究方面，由于技术限制，在极短时间内测量碰撞过程中的响应比较困难。Hu等[9]用激光测振仪对钢球和铝杆的轴向碰撞进行了实验测量，在速度响应上与仿真结果非常吻合。董富祥等[10]设计了类似的杆和杆之间的轴向碰撞实验，对应变响应和速度响应都进行实验测量，通过与仿真结果的对比验证了实验结果的可靠性。由于实验条件的限制，国内外对于考虑摩擦的斜碰撞问题的实验研究开展较少。此外，对于斜碰撞问题，碰撞后柔性体上各点的切向速度不为零，不再作直线运动，无法用一维激光测振仪进行速度测量。

本文基于模态综合法建立了柔性多体系统斜碰撞的动力学模型，并采用罚函数方法和修正库伦摩擦模型对接触碰撞进行建模。在实验技术上，分别采用传统电阻应变片和3D-DIC[11]数字图像相关技术对碰撞响应进行测量，3D-DIC技术的使用不仅能够同时测量一个区域内的动力学响应，还可以测量法向和切向速度。最后，对一个钢杆和PVC圆盘的碰撞系统分别采用仿真计算和实验测量的方法进行研究，通过对比验证了本文模型的准确性。斜碰撞问题的实验研究也验证了修正库伦摩擦模型的可靠性，并发现了摩擦因数随相对速度增加而逐步减小的规律。

1 基于模态综合法的柔性多体系统动力学建模

(1)

图1 单柔性体示意图

在有限元离散后，可将柔性体的节点坐标分为边界节点和内部节点，设pI和pJ分别为内部节点坐标阵和边界节点坐标阵，NI和NJ为对应的形函数矩阵，任意一点的变形位移可插值为

(2)

在分析接触碰撞问题时，可将碰撞区域的节点与其他铰约束节点都划为边界节点，而对内部节点可通过模态综合法进行缩减，根据固定界面模态综合法，内部节点和边界节点的变形位移可表示为

(3)

pi=Hiyi

(4)

对于任意柔性体，可定义广义速度阵如下

(5)

(6)

(7)

其中

(8)

(9)

(10)

采用拉格朗日方法，无约束的单柔性体虚功率形式的动力学方程如式(11)所示

(11)

惯性力的虚功率可写为

(12)

其中

(13)

设ki为弹性变形刚度阵，采用模态综合法缩减后的柔性体变形坐标yi表示的内力虚功率为

(14)

设fSi为广义内力阵，内力的虚功率为

(15)

式(11)可写为

(16)

动力学变分方程(16)可用广义坐标表示为

(17)

其中

(18)

(19)

设柔性多体系统由Bi(i=1,…,N)组成，系统的动力学变分方程可写为

(20)

(21)

其中

(21)

(22)

引入物体间的约束方程，则柔性多体系统动力学方程可写为

(23)

式中：Φq为约束方程Φ=0对应的雅可比阵；λ为约束方程对应的拉格朗日乘子。

2 接触碰撞模型

对于接触碰撞模型，分为法向接触和切向接触分别进行建模。在有限元离散的条件下，用罚函数方法建立法向碰撞力模型，用库伦摩擦定律建立切向碰撞力模型。

2.1 法向接触模型

在有限元方法中，接触面被离散为一系列单元和节点，将可能发生接触的两个柔性体分为主物体和从物体，当从物体的节点穿透主物体的单元表面时认为发生接触，并记为一个点面接触对，如图2所示。通过接触检测，在每一时刻用有限个接触对来描述接触界面。

图2 发生接触的单元

对从物体表面的从节点P，在主单元面R上的投影点记为Q，如图3所示，n为主单元面的外法线向量，则P点相对于主单元面的距离函数可写为

gn=nT(rP-rQ)

(24)

图3 接触对的描述

当gn小于0时认为接触发生，记为接触对P-R，并在接触对上施加法向碰撞力

(25)

碰撞力的大小用罚函数方法确定

Fn=k|gn|

(26)

其中，k为罚刚度，与碰撞物体之间的弹性模量以及有限元单元尺寸相关，对于三维实体单元，可以表示为

(27)

式中：f为接触刚度的比例因子，与主单元的尺寸相关，在一般情况下可取0.1，若单元尺寸过小，则需要增加f的值；K为接触主单元的体积模量；A为主单元面的面积；V为接触对中主单元的体积。对于其他的单元类型，碰撞刚度有类似的表达形式，具体的求法可以参考有限元文献[12]。

需要指出的是，在有限元离散条件下，接触面是由很多个接触对组成，因此可以用罚函数这种线性弹簧模型来描述，而对于刚体模型以及只需要考虑加一个弹簧力来描述的弹性体碰撞问题，则需要用连续碰撞力模型中的非线性弹簧模型来描述。

2.2 切向接触模型

对于接触碰撞中的摩擦问题，最常用的模型是库伦模型及各种改进形式。在库伦摩擦定律中，对干摩擦问题根据是否有切向相对速度将摩擦分为静摩擦和滑动摩擦，静摩擦力一般通过施加约束方程来求解，而滑动摩擦力通过法向摩擦力乘以摩擦因数来确定。在有限元中，接触面由很多个接触对来描述，对于每一个接触对不必严格区分静摩擦和滑动摩擦，而用统一的摩擦力来表示。

以接触对P-R为例，在当前时刻n+1t, 已知上一时刻摩擦大小为nFt，当前时刻法向碰撞力为n+1Fn，相对于上一时刻从节点在主面上投影点Q的滑移量为Δe，则计算当前时刻摩擦力的流程为

步骤1 计算试探摩擦力f*=kΔe

步骤2 判断是否超出最大摩擦力限制

若f*>μn+1Fn，则n+1Ft=μn+1Fn

若f*≤μn+1Fn，则n+1Ft=f*

其中，μ为摩擦因数，可用静摩擦因数μs和滑动摩擦因数μd表示为统一形式

(28)

式中：c为衰减系数，为经验参数，在有限元中也可以不区分静摩擦因数和滑动摩擦因数。在通常情况下认为摩擦因数只与物体表面的粗糙程度相关，而随着实验技术的发展，学者们通过摩擦实验研究发现同样的材料在不同工况下会有不同的摩擦因数[13]。

作用在主从物体上的摩擦力可表示为

(29)

式中：τ为单位矢量，沿着P点相对于Q点的切向速度方向(或相对滑动趋势方向)。

作用于P和Q点的法向碰撞力和摩擦力的合力分别为

(30)

法向碰撞力和摩擦力的虚功率为

(31)

(32)

(33)

2.3 碰撞区域的判定

用有限元方法解决碰撞问题时，接触界面以及碰撞局部区域的网格对求解精度的影响很大，本文提出基于模态综合法的柔性多体系统接触碰撞动力学模型中，保留了接触界面和局部碰撞区域的节点自由度，用以提高模型的计算精度。在实际应用中，对离散后的柔性体碰撞区域的判定也很重要，区域取的太小会降低求解精度，取的太大会影响计算效率。本文采取以下两个策略来判定碰撞区域：

区域1预估接触区域。对于几何形状规则的界面接触，可以用经典Hertz理论预估最大接触面积[14]，对于常用的球面接触，最大接触半径有以下公式：

(34)

式中:F为最大碰撞力预估值；R为接触面的曲率半径；E和ν为被撞物体弹性模量和泊松比。由于Hertz理论是针对静态接触问题的解析分析，因此在动力学分析中需要取一个更加保守的值，通常在低速碰撞情况下，取两倍的碰撞半径足以保证接触区域的准确描述。对于不规则的界面接触，只能通过几何形状以及物体大范围运动来判断可能的接触界面。

区域2确定局部碰撞区域。对于弹性体的碰撞，在碰撞方向上的局部区域内通常会有很大的应力应变梯度，如果要保证这部分区域内的应力应变的计算精度，则需要在碰撞方向上取足够多的单元，使得计算结果收敛，并保留这部分的节点变形自由度；若不关注，对于大部分低速碰撞问题，只需在碰撞方向上取一层单元，保留这部分单元的节点自由度即可保证碰撞力以及整体仿真的精度。

3 仿真与实验分析

为了验证本文提出的动力学模型并通过实验方法研究柔性多体系统斜碰撞问题，分别采用数值仿真和实验测量的方法对一个钢杆和PVC圆盘的斜碰撞过程进行研究。钢杆和圆盘的几何与材料参数如表1所示。在考虑摩擦的斜碰撞工况下，对仿真与实验结果进行了对比。

表1 碰撞物体几何与材料参数

考虑摩擦的斜碰撞工况如图4所示，钢杆沿着圆盘半径方向偏转30°，以0.613 m/s的初速度撞击圆盘侧边。其中，钢杆置于一个有弹射机构的滑槽内，圆盘置于三根固定在实验台的尖头螺钉上，如图4所示。在板的两侧分别用三向应变片和3D-DIC数字图像相关技术测量三个点的三向应变以及速度响应。应变片采用可以同时测量0°，45°，90°方向应变的应变花，DIC测量使用两台Photron FASTCAM SA-X2-1000k-M2型号的高速摄像仪。A面三个测点，S2点在初始碰撞点与圆心的连线上且距离初始碰撞点10 mm的位置，S1和S3点分别距S2点15 mm；B面的三个测点，P1，P2，P3分别距离碰撞点10 mm，20 mm，30 mm，如图5所示。

图4 斜碰撞实验示意图

随着高速摄像与数字图像技术的发展，DIC技术也被越来越广泛地应用于工程中。DIC技术的主要原理是通过在物体表面涂上大小相差不大的不规则散斑，用高速摄像机在一定时间内对散斑区域进行拍照，最后通过数字图像相关技术进行处理得到不同照片散斑点之间的位移信息，如图6所示。本文实验所采用的圆盘的散斑区域，如图7所示。

图5 斜碰撞实验装置图

图6 DIC技术的基本原理

图7 实验圆盘散斑区域

首先通过钢材料圆柱球面与PVC平面的切向碰撞实验测量得到在切向相对速度为0.31 m/s时滑动摩擦因数为0.11。在此工况下，用本文提出的柔性多体系统接触碰撞动力学模型进行数值求解，求解参数如表2所示。数值仿真得到的圆盘受到的法向碰撞力和切向摩擦力曲线，如图8所示。

表2 仿真参数

图8 碰撞力时间历程

对于应变响应，S1，S2，S3三点实验直接测量得到的x，y方向应变与仿真结果对比如图9～11所示。可见在考虑摩擦的斜碰撞工况下，仿真结果与实验结果基本吻合。在此工况下S1和S3两点x方向应变的绝对值小于y方向应变，而S2点的x方向应变的绝对值大于y方向的应变。此外，由三组曲线可以看出，在距离圆盘边缘处，S1，S3两点的x方向应变的方向与S2点相反，物理上可以解释为由于碰撞作用，碰撞点附近区域会因为局部挤压而向外变形。

由于圆盘受到摩擦力作用的方向为y轴负方向，且在碰撞过程中碰撞点会沿着y轴负方向滑动，因此S3点的y方向应变的绝对值比对称位置的S1点大。

三向应变花能同时测量得到测点x，y方向以及45°方向的应变，通过换算可以得到测点的x-y方向剪切应变。三点剪切应变的仿真结果与实验结果对比如图12～14所示。

图9 S1点x,y方向正应变仿真结果与实验结果对比

Fig.9 The comparison of normal strain inxandydirection ofS1between simulation result and experiment result

图10 S2点x,y方向正应变仿真结果与实验结果对比

Fig.10 The comparison of normal strain inxandydirection ofS2between simulation result and experiment result

图11 S3点x,y方向正应变仿真结果与实验结果对比

Fig.11 The comparison of normal strain inxandydirection ofS3between simulation result and experiment result

图12 S1点剪切应变仿真结果与实验结果对比

Fig.12 The comparison of shear strain ofS1between simulation result and experiment result

图13 S2点剪切应变仿真结果与实验结果对比

Fig.13 The comparison of shear strain ofS2between simulation result and experiment result

图14 S3点剪切应变仿真结果与实验结果对比

Fig.14 The comparison of shear strain ofS3between simulation result and experiment result

可以看出，另外S2点和S3点的剪切应变与仿真结果基本吻合。由于S1点的三向应变都较小，测量误差比S2点和S3点大，因此折算出的剪切应变与仿真结果也比S2点和S3点大。此外可以看出，S1和S3点的剪切应变方向也是相反，S3点的剪切应变绝对值大于S1点。

P1，P2，P3三点沿着碰撞法向即x方向的速度响应仿真与实验结果对比如图15～17所示。

图15 P1点y方向速度响应对比

Fig.15 The comparison of velocity inydirection ofP1between simulation result and experiment result

图16 P2点y方向速度响应对比

Fig.16 The comparison of velocity inydirection ofP2between simulation result and experiment result

图17 P3点y方向速度响应对比

Fig.17 The comparison of velocity inydirection ofP3between simulation result and experiment result

可以看出，3D-DIC不接触形式测量圆盘表面的速度响应结果与仿真结果也非常吻合，同时也验证了摩擦因数测量实验得到的摩擦因数0.11的可靠性。

为了研究不同相对速度的接触碰撞对摩擦因数的影响，增加另外一种工况即钢杆在0.78 m/s的初始速度下以同样角度撞击圆盘侧面，初始切向速度为0.39 m/s，测得切向速度为0.39 m/s时摩擦因数为0.1。S2点y方向正应变仿真结果与实验结果对比如图18所示，可以看出仿真结果与实验结果吻合较好。

图18 S2点y方向正应变仿真结果与实验结果对比

Fig.18 The comparison of normal strain inydirection ofS2between simulation result and experiment result

为了进一步说明相对切向速度对摩擦因数的影响，当初始切向速度为0.39 m/s时，分别以0.11，0.10和0.09三组摩擦因数进行仿真得到S2点y方向应变响应的三组仿真结果与实验测量结果的对比如图19所示。可以看到，同样两种材料的接触碰撞在不同摩擦因数下的仿真结果也会不同，当摩擦因数为0.11时，理论结果与实验结果差异较大，而当摩擦因数为0.1时，理论结果与实验结果吻合较好。在第一种工况下，初始相对切向速度为0.31 m/s，此时摩擦因数为0.11。由此可见，在一定范围内随着碰撞物体相对速度的增大实际摩擦因数会减小，这也与摩擦因数测量实验中的结论相同。

图19 S2点y方向正应变仿真结果与实验结果对比

Fig.19 The comparison of normal strain inydirection ofS2between simulation result and experiment result

4 结 论

(1) 基于模态综合法和混合坐标法建立了适用于求解柔性多体系统斜碰撞问题的动力学模型，分别罚函数方法和修正库伦摩擦模型描述法向碰撞力和切向摩擦力。该模型在保证局部碰撞区域描述精确的条件下能大幅度缩减非碰撞区域的自由度。

(2) 用传统的电阻应变片和3D-DIC基于高速摄像的数字图像相关技术分别测量碰撞物体段时间内的应变响应和速度响应，通过与仿真结果的对比验证了本文提出的柔性多体系统考虑摩擦的动力学模型的准确性。同时，3D-DIC测量结果也证明了这种前沿的实验技术在测量极短时间内物体局部区域在冲击作用下的动力学响应的可靠性。

(3) 通过对不同工况的仿真与实验对比，发现了摩擦因数的变化规律，随着初始相对速度增大，摩擦因数逐渐减小。

[1] DJERASSI S. Collision with friction, Part A: Newton’s hypothesis[J]. Multibody Syst Dyn, 2009, 21: 37-54.

[2] DJERASSI S. Collision with friction, Part B: Poisson’s and Stornge’s hypothesis[J]. Multibody Syst Dyn, 2009, 21: 55-70.

[3] DJERASSI S. Three-dimensional, one-point collision with friction[J]. Multibody Syst Dyn, 2012, 27: 173-195.

[4] TIMOSHENKO S P, GERE J M. Theory of elastic stability[M]. Westford, MA: Courier Dover Publications, 2009.

[5] SCHIEHLEN W, SEIFRIED R, EBERHARD P. Elastoplastic phenomena in multibody impact dynamics[J]. Computer Methods in Applied Mechanics and Engineering, 2006, 195(50): 6874-6890.

[6] SEIFRIED R, SCHIEHLEN W, EBERHARD P. The role of the coefficient of restitution on impact problems in multi-body dynamics[J]. Proceedings of the Institution of Mechanical Engineers, Part K: Journal of Multi-body Dynamics, 2010, 224(3): 279-306.

[7] 韩石磊, 洪嘉振. 柔性多体碰撞问题的多变量方法[J]. 力学学报, 2011, 43(5):886-893.

HAN Shilei, HONG Jiazhen. Multi-variable method for flexible multi-body systems with contact/impact[J]. Chinese Journal of Theoretical and Applied Mechanics, 2011, 43(5): 886-893.

[8] CHEN P, LIU J Y, HONG J Z. An efficient formulation based on the Lagrangian method for contact-impact analysis of flexible multi-body system[J]. Acta Mechanica Sinica, 2016(2):1-9.

[9] HU B, SCHIEHLEN W, EBERHARD P. Comparison of analytical and experimental results for longitudinal impacts on elastic rods[J]. Journal of Vibration and Control, 2003, 9(1/2): 157-174.

[10] 董富祥. 刚柔耦合多体系统碰撞动力学建模理论与实验研究[D]. 上海：上海交通大学，2010

[11] WANG H, XIE H, WU L, et al. Study on the effect of DIC deformation sensor on mechanical property of substrate[J]. Measurement, 2014, 49: 283-288.

[12] ZIENKIEWICZ O C, TAYLOR R L. The finite element method: solid mechanics[M]. Oxford: Butterworth-Heinemann, 2000.

[13] LIN Y, QIN J, LU F, et al. Dynamic friction coefficient of two plastics against aluminum under impact loading[J]. Tribology International, 2014, 79: 26-31.

[14] POPOV V. Contact mechanics and friction: physical principles and applications[M]. Dordrecht: Springer Science & Business Media, 2010.