刚柔耦合接触动态粘合算法研究

高 浩，干 锋，戴焕云

(西南交通大学 牵引动力国家重点实验室，成都 610031)

刚柔耦合多体系统中包含了刚体、弹性体以及体与体之间的连接力元和铰，刚柔耦合仿真技术已经基本成熟，几乎所有通用多体动力学软件都支持刚柔耦合仿真，一般采用模态叠加法近似描述弹性体的动态响应。模态叠加法是基于小变形假设的线性处理方法，忽略了结构的非线性特性，同时模态截断造成了高频振动信息的丢失。对于非线性弹性变形或局部模态对仿真结果影响很大的系统，弹性体需要直接采用有限元瞬态动力学分析。可以通过一定的算法将多系统中多刚体和弹性体分为不同的子系统，分别在不同的软件中进行仿真。这是一种松耦合方式，关键技术是如何通过数学工具建立严密的接口关系，通过接口传递界面力和运动信息包括位移、速度和加速度，并保证接口界面处同时满足平衡条件和协调条件。刚柔耦合接触问题是刚柔耦合系统动力学仿真的难题之一，同时也有很多工程应用如齿轮接触动力学和车辆系统弓网接触动力学研究。文献[4]中采用有限元预计算的方法考虑齿轮的弹性变形，研究刚柔耦合特性对齿轮啮合传动特性的影响规律。文献[5]采用直接有限元积分的方法对弓网接触动力学进行研究。Wang[1－2]提出了一种刚柔耦合仿真的松耦合方法称为粘合算法，子系统界面之间传递的信息通过协调器即粘合算法处理并传递给与之相连的其他子系统，这种方法将子系统的运动信息作为输入，输出的界面力作为子系统的下一步积分的边界条件。该方法可以用于多体系统刚柔耦合仿真分析。Run[3]在此基础上改进了粘合算法，针对多体动力学中弹性铰节点提出了弹性粘合算法，减小了连接界面的数值振荡，提高了算法的稳定性。但是粘合算法是针对连接界面固定的子系统提出的，不适合连接界面变化接触问题。为此本文在粘合算法的基础上提出了动态粘合算法，即粘合算法中的装配矩阵和柔度矩阵是根据接触位置动态变化的，以保证接触界面平衡条件和协调条件同时满足。并将其用于处理刚体与弹性体或弹性体与弹性体之间的动态接触问题。当接触发生在没有节点的位置时，如何处理不匹配的网格非常关键，一般采用形函数插值法[6]，Kim[7－8]提出界面单元法，将无网格法中的移动最小二乘近似引入到有限元中。本文在粘合矩阵中考虑单元形函数使得动态粘合算法可以处理非节点位置接触的情况。

1 多体系统粘合算法[2]

多体系统一般有刚体或弹性体多个部件通过力元、铰接或约束组合而成，假设部件i和部件i+1通过界面节点p和p'连接，如图1所示。

图1 多体系统部件粘合Fig.1 Components gluing of multibody system

要保证平衡条件和协调条件同时满足，首先假设合适的界面力向量，使得每个子系统满足平衡条件，在仿真过程中只需要考虑协调条件即界面连接点的运行信息，包括位移、速度和加速度。界面力是运动信息的函数，如式所示:

式(1)可表示为:

接触界面应当满足协调条件，界面运动信息装配向量为U，则:

其中表示一对连接点的位移，假设e为协调条件的误差向量，界面应满足协调条件:

求解式(3)中的界面力向量F，使其满足协调条件式(5)，利用牛顿法－拉夫逊法求解该非线性方程组式，迭代关系式为:

其中:k为迭代次数，Λ称为粘合矩阵，Fk和Fk+1分别为第k步和第k+1步的界面力，具体计算方法如下:

其中:B和C称为装配矩阵，其元素由0、－1和1组成;G称为界面柔度矩阵，式(10)中m和n分别为G的第m行，第n列元素，利用式(1)和式(10)计算获得。对于线性系统可近似的通过在连接界面各自由度施加单位载荷获取柔度矩阵，则其中为在连接界面自由度m施加单位载荷时产生的位移为不加任何载荷时的初始位移。i－th表示i部件对应的元素。

通过式(6)计算获得的界面力，以外载荷的形式施加到各部件的动力学方程中，实现各部件之间的粘合。

2 接触动态粘合法则

两个部件i和i+1相互接触，其中i为接触体，i+1为目标体，假设i结构上接触点为p点，对应的部件i+1上的p'位置，p与p'为接触对，如图2所示。若p'不是节点，p'坐标由该点所在单元的节点坐标和形函数插值得到，假设接触单元为n个节点单元，p'点坐标如式(11)所示:

其中:Ni接触单元中节点i的形函数，ui为节点i的坐标，Xp'为p'在接触单元中的自然坐标向量。利用粘合算法求解两个部件之间的接触问题，首先要确定连接界面，接触体上的接触点p和目标体中对应点p'所在单元相关的节点组成连接界面，两个部件之间的通过非线性接触力元连接，接触刚度为Kc，随着接触状态的变化，接触界面是动态变化的，因此本文发展动态粘合算法来求解接触问题。

接触点之间不是固接的，令接触点之间的距离qc为:

其中界面装配误差向量为:

图2 两接触部件粘合Fig2.Gluing of two contact components

协调方程可写成:

假设接触点之间为非线性接触弹簧，接触力为:

式(7)中的粘合矩阵改写为:

动态粘合算法需要动态的更新装配矩阵和界面柔度矩阵，装配矩阵B和C的元素不再只有1、－1和0，而是和单元的形函数相关，如式(17)所示。

其中:N=[N1(Xp')，N2(Xp')，…，Nn(Xp')]，接触界面柔度矩阵G也需要动态更新，将式(11)代入式(1)中，结合式(10)求得G矩阵。在整个接触过程中装配矩阵和界面柔度矩阵的维数和数值都会发生变化，因此每一积分步都需要对这些矩阵进行更新。

3 实例仿真

使用本文发展的动态粘合算法，对如图3所示的系统进行动力学仿真，该系统包含简支梁和弹簧－质量块两个部件，部件之间通过接触单元连接，mc为质量块的质量，yc为质量块运动自由度，ks为支撑刚度，kc为接触刚度，ys为支撑基础的运动，V0为弹簧－质量块向前运动速度。

选取仿真参数如表1所示。

图3 简支梁和质量－弹簧系统接触模型Fig3.Simple-supported beam and mass-spring systems contact model

表1 动力学仿真参数Tab.1 Parameters of dynamic simulation

支撑基础正弦输入为 ys=0.001+0.002sin(4πt)，在该正弦输入激扰下，简支梁和质量块会出现接触和脱开的状态转换。采用统一方程法和动态粘合算法两种方法求解该系统的响应。

3.1 统一方程法

采用有限元法建立平面梁的动力学方程，梁单元有两个节点，每个节点有水平、垂向和转动三个自由度，选用梁单元的形函数为:

其中:l为单元长度，ξ为梁单元坐标系下的坐标。整个系统的统一动力学方程为:

N下标表示式(18)形函数中对应位置的元素，ξp'表示p'点在梁单元坐标系中坐标。NI的非零元素分别对应节点 i和节点 j的垂向和转动自由度，即 uiy，uiθ，ujy和ujθ。用Newmarker法求解式(19)，每个积分步判断接触状态，当接触脱开时令kc=0。

3.2 动态粘合算法

系统中包含两个部件，平面简支梁和弹簧－质量系统，两个部件只有一个连接界面，即接触对，两个部件的动力学方程分别为:

假设弹簧质量系统为部件1，简支梁为部件2，部件1的界面位移和界面力为，部件 2 的界面位移和界面力为[FiyMiθFjyMjθ]T，其中 uiy，uiθ，ujy和 ujθ分别为接触单元i节点和 j节点的垂向和转动位移，Fiy，Miθ，Fjy和Mjθ分别为作用在接触单元i节点和j节点的垂向节点力和扭矩。

两个部件之间存在一个连接界面，因此界面装配协调条件为一个方程:

式中:N’(ξp')=[N22(ξp') N23(ξp') N25(ξp')N26(ξp')]。接触点之间的距离为q=up－up'，此时协调方程为e－q=0，接触力为FcA=kcq，根据有限元法将集中力转化成等效节点力，即 FcR=－N’(ξp')FcA，取式(4)中的F=FcA，根据式(8)、式(9)和式 (17)计算部件1和部件2的装配矩阵为:

式(24)～式(27)代入式(16)中，粘合矩阵可表示为:

界面力迭代关系式为:

4 仿真结果及讨论

用统一方程法和动态粘合算法对上述系统进行动力学仿真，图4为质量块运动轨迹及简支梁在不同时刻的变形。质量块的垂向位移时程曲线如图5，梁中部节点的垂向位移时程曲线如图6所示，两种方法计算位移结果基本一致。图7为两种算法计算的位移结果的绝对误差，可见误差非常小，最大误差发生在由脱离状态到接触状态变化瞬间，稳定接触过程中绝对误差小于5e－7 m，完全可以满足仿真精度需求。图8为梁和质量块之间的接触力，两种算法在接触力突变时仍保持很好的一致性，说明粘合算法比较稳定。

图4 不同时刻接触系统的状态Fig.4 Status of contact system at different time

图5 质量块垂向位移Fig.5 Vertical displacement of mass block

图6 梁中部垂向位移Fig.6 Vertical displacement of central node of the beam

图7 位移绝对误差Fig.7 Absolute error of displacement

图8 梁与质量块之间的接触力Fig.8 Contact force between beam and mass block

在使用粘合算法时，用上一积分步的界面力作为界面力迭代的初始条件，经过1～2次迭代即可满足仿真精度要求。可见动态粘合算法可以很好的替代统一方程求解，可以将复杂系统分解成多部件，分别求解在不同的软件或计算机上进行求解，实现分布式计算。

对比式(19)和(22)可以发现接触问题为状态非线性问题，使用统一方程求解时，刚度矩阵是变化的需要动态更新，无法进行事先分解。使用动态粘合算法，刚度矩阵不变，对以刚度矩阵规模较大的有限元系统，可进行事先分解，大大提高计算速度。

5 结论

(1)本文发展了一种动态粘合算法，可以用来解决刚性体与弹性体或弹性体与弹性体之间的接触问题，该方法通过动态的构建粘合矩阵，并迭代求解一组界面力，使得两个接触的系统在接触界面同时满足平衡条件和协调条件。

(2)通过简单刚柔耦合接触系统的动力学仿真，验证了动态粘合算法的计算精度和可靠性。

(3)复杂的接触系统如铁路车辆系统中受电弓与接触网滑动接触，轮轨滚动接触等，使用动态粘合算法可以将这些系统分解成若干部件进行分别在不同的软件或平台上求解，通过数据传输和数据交换完成界面力和运动信息的传递，为接触问题的分布式仿真提供了可行方案。

[1]Wang J Z，Ma Z D，Hulbert G M.A distributed mechanical system simulation platform based on a“gluing algorithm”[J].Journal of Computing and Information Science in Engineering，2005，5:71－76.

[2]Wang J Z.Development of advanced methodlogy for networkdistributed simulation[D].The University of Michigan，2005.

[3]Ryu G S.Integration of heterogeneous simulation models for network-distributed simulation [D]. The University of Michigan，2009.

[4]姚廷强，迟毅林，黄亚宇，等.刚柔耦合齿轮三维接触动力学建模与振动分析[J].振动与冲击，2009，28(2):167－171.

[5]周 宁;张卫华.基于直接积分法的弓网耦合系统动态性能仿真分析[J].中国铁道科学，2008，129(16):71－76.

[6]王勖成.有限单元法[M].北京:清华大学出版社，2003.

[7]Kim H G.Interface element method(IEM)for a partitioned system with non-matching interfaces[J].Comput.Methods Appl.Mech.Engrg.，2002，191:3165 －3194.

[8]Kim H G.Interface element method:Treatment of nonmatching nodes at the ends of interfaces between partitioned domains[J].Comput.Methods Appl.Mech.Engrg.，2003，192:1841－1858.