基于Bathe积分策略的多体动力学方程求解方法

摘要

本发明公开了一种基于Bathe积分策略的多体动力学方程求解方法，包括在动力学多体系统中建立相应的坐标系，并采用相对坐标法，利用拉格朗日乘子法，根据虚功率原理，得到动力学方程；将约束方程对时间求二阶导数，得到指标‑1的微分‑代数方程组，添加约束稳定项，形成含完整约束的动力学方程的一般形式；对含完整约束的动力学方程一般形式在时间步[t,t+h/2]上进行积分迭代求解，得到t+h/2时刻的运动参数；在时间步[t+h/2,t+h]上进行积分迭代求解，得到t+h时刻的运动参数；重复迭代计算，当计算时间达到设置好的仿真总时间，输出系统广义位移、速度和加速度随时间变化的值，即可分析过程闭环多体动力学系统在一定时间内的运动状态；本方法在较大积分步长的情况下也可以得到较好的精度。

说明书

技术领域

本发明属于机械动力学领域，特别涉及到了一种基于Bathe积分策略的多体动力学方程求解方法。

背景技术

多体系统是以一定方式相连接的多个物体组成的系统，通常根据拓扑结构的特点，将没有回路的多体系统称为开环系统，带回路的系统称为闭环系统，而系统往往是存在大量闭环结构的复杂多体系统。对于存在闭环结构的多体系统，由于此时系统广义坐标不再是独立的，建模和求解的复杂性大幅度增加。

由于机械机构的复杂性和系统运行环境的多变性，含闭环的多体系统动力学方程一般具有强耦合性和非线性特征，难以获得解析解，因此需要具有良好数值性态的求解算法，从而获得精确、稳定和高效的数值解，实现机械系统的运动学和动力学仿真分析。对于闭环多体系统动力学方程的求解，常用算法有中心差分法，向后差分法，龙格-库塔(Runge-Kutta)方法，亚当斯(Adams,J.C.)方法，Newmark法，HHT法，广义-α法等，目前这几类方法已经被应用于闭环多体系统动力学分析，然而当求解长时间历程的非线性动力学问题时，传统方法可能会变得不稳定，产生数值耗散。

发明内容

本发明所解决的技术问题在于提供一种基于Bathe积分策略的多体动力学方程求解方法，以解决含闭环的多体系统动力学方程求解算法不稳定且效率不高的问题。

实现本发明目的的技术解决方案为：

一种基于Bathe积分策略的多体动力学方程求解方法，包括以下步骤：

步骤1、构建含完整约束的动力学方程：在动力学多体系统中建立相应的坐标系，并采用相对坐标法，利用拉格朗日乘子法，根据虚功率原理，得到动力学方程；

步骤2、将动力学方程整理为含完整约束的动力学方程的一般形式：将约束方程对时间t求二阶导数，从而得到指标-1的微分-代数方程组，在动力学方程中添加Baumgarte约束稳定项，整合相应参数并整理方程，形成含完整约束的动力学方程的一般形式；

步骤3、利用Bathe积分策略对含完整约束的动力学方程一般形式在时间步[t,t+h2]上进行积分迭代求解，得到t+h/2时刻的运动参数；其中h为积分步长；

步骤4、利用Bathe积分策略对含完整约束的动力学方程一般形式在时间步[t+h2,t+h]上进行积分迭代求解，得到t+h时刻的运动参数；

步骤5、重复步骤3与步骤4，直到计算时间达到设置好的仿真总时间时，完成计算，输出系统广义位移、广义速度和广义加速度随时间变化的值，即可分析过程闭环多体动力学系统在一定时间内的运动状态。

本发明与现有技术相比，其显著优点：

(1)本发明的动力学方程建模时使用了指标-1的微分-代数方程，通过添加Baumgarte违约稳定项，使得位移约束和速度约束不违约。

(2)求解方法利用了Bathe积分策略，Bathe算法属于隐式算法，相比于显式积分算法，Bathe算法在较大积分步长的情况下也可以得到较好的精度。

(3)迭代过程利用了逆Broyden拟牛顿法求解方程，由于将动力学方程推导成了以广义位移项作为未知变量的形式，因此雅克比矩阵的迭代初值无需进行求导计算，可直接由广义质量矩阵和广义阻尼矩阵获得，且由于雅克比矩阵中包含广义质量矩阵和广义阻尼矩阵，收敛效率提高，从而能够对动力学方程进行稳定、高效的求解。

下面结合附图对本发明作进一步详细描述。

附图说明

图1为本发明方法的总体流程图。

图2为实施例1中闭环多体系统结构示意图。

具体实施方式

为了说明本发明的技术方案及技术目的，下面结合附图及具体实施例对本发明做进一步的介绍。

结合图1，本发明的一种基于Bathe积分策略的多体动力学方程求解方法，包括以下步骤：

步骤1、构建含完整约束的动力学方程：在动力学多体系统中建立相应的坐标系，并采用相对坐标法，利用拉格朗日乘子法，根据虚功率原理，得到动力学方程。

在动力学多体系统中任意位置建立惯性坐标系，并且在动力学多体系统中的每个刚体的质心上建立连体坐标系，每个铰上建立铰坐标系，以相邻铰的坐标系的相对坐标(即铰的相对运动参数)作为系统广义坐标，采用相对坐标法，考虑系统所受的约束，采用拉格朗日乘子法引入约束方程，根据虚功率原理得到动力学方程，其为指标-3的微分-代数方程组(DAEs)：

式中t为时间，q、和分别为铰的相对位移、相对速度和相对加速度，Φ(q,t)＝0为系统约束方程，Φ_q为Φ(q,t)对q求导，()^T表示对该矩阵的转置，λ为拉格朗日乘子，M(q,t)为质量矩阵，为阻尼矩阵，Q为系统外力。

步骤2、将动力学方程整理为含完整约束的动力学方程的一般形式：将约束方程对时间t求二阶导数，从而得到指标-1的微分-代数方程组，在动力学方程中添加Baumgarte约束稳定项，整合相应参数并整理方程，形成含完整约束的动力学方程的一般形式。

2.1、对式(1)中的约束方程Φ(q,t)＝0对时间t求二阶导数，可以将式(1)写成如下指标-1的微分-代数方程组(DAEs)：

式中其中Φ_qq为Φ(q,t)对q求二阶导数，Φ_qt为Φ(q,t)依次对q和t求一阶导数，Φ_tt为Φ(q,t)对t求二阶导数。

2.2、在动力学方程中添加Baumgarte约束稳定项

为稳定求解动力学方程，根据Baumgarte约束稳定法，在式(2)中加入稳定项可将式(2)写成如下形式

式中，为速度约束，α和β为Baumgarte系数，取h为积分步长，λ′和λ″分别作为λ对时间的一次积分项和二次积分项，由于在式(3)中与拉格朗日乘子相关的项里只有λ参与了运算，λ′与λ″并没有进行实际的运算，因此并不改变原方程的形式。

2.3、整合相应参数并整理方程，形成含完整约束的动力学方程的一般形式

令

式中P、和分别为动力学方程的广义位移项、广义速度项和广义加速度项，和分别为动力学方程的广义质量矩阵和广义阻尼矩阵，为系统广义力。

则式(3)可整理为：

上式即为含完整约束的动力学方程的一般形式。

步骤3、利用Bathe积分策略对含完整约束的动力学方程一般形式在时间步[t,t+h/2]上进行积分迭代求解，得到t+h/2时刻的运动参数：

3.1、令系统的广义位移初值和广义速度初值分别为P₀和假设t时刻的广义位移、广义速度和广义加速度分别为P_t、和首先求解t+h/2时刻的广义位移P_t+h/2。

在[t,t+h/2]上利用梯形公式可推导得到t+h/2时刻的广义速度和广义加速度的表达式为：

将式(6)代入式(5)进行整理后可得到以广义位移P_t+h/2为未知变量的非线性代数方程：

3.2、利用逆Broyden拟牛顿法迭代求解式(7)

迭代过程如下：

式中y_k和s_k都为迭代的中间变量，J_t+h/2表示t+h/2时刻的雅克比矩阵，上标k和k+1分别表示第k次迭代和第k+1次迭代，Y的表达式如下：

根据泰勒展开原理，式(7)，式(8)中P_t+h/2和J_t+h/2的迭代初值可分别设置为：

当计算得到第k+1次迭代所得的和后，利用式(6)可得到即可进行第k+2次的迭代。

设置迭代收敛允许误差tol，并设置收敛条件：

当当前迭代步满足收敛条件式(12)时，迭代停止，计算得到的第k+1次的广义位移广义速度和广义加速度即为t+h/2时刻的运动参数。

步骤4、利用Bathe积分策略对含完整约束的动力学方程一般形式在时间步[t+h/2,t+h]上进行积分迭代求解，得到t+h时刻的运动参数：

4.1、求解t+h时刻的广义位移P_t+h。

利用梯形公式，在[t+h/2,t+h]上可推导得到t+h时刻的广义速度和广义加速度的表达式：

将式(12)代入式(5)进行整理后可得到以广义位移P_t+h为未知变量的非线性代数方程：

4.2、利用逆Broyden拟牛顿法迭代求解式(14)

迭代过程如下：

式中J_t+h表示t+h时刻的雅克比矩阵。

根据泰勒展开原理和式(14)，式(15)中P_t+h和J_t+h迭代的初值可设置为：

当计算得到第k+1次迭代所得的和后，利用式(13)可得到和即可进行第k+2次的迭代。设置迭代收敛允许误差tol，并设置收敛条件：

当当前迭代步满足收敛条件式(18)时，迭代停止，计算得到的第k+1次的广义位移广义速度和广义加速度即为t+h时刻的运动参数。

步骤5、重复步骤3与步骤4，直到计算时间达到设置好的仿真总时间时，完成计算，输出系统广义位移、广义速度和广义加速度随时间变化的值，即可分析过程闭环多体动力学系统在一定时间内的运动状态。

实施例1：

通过下面实施例对本方法做进一步详细阐述。以下为本方法的一个仿真算例，算例为一个闭环多体系统问题，具体说明如下：

此处选择了一个含有3个刚性杆、4个转动铰的闭环多体系统，具体结构如图2所示，第一旋转铰J₁连接第一地面刚性支架与第一刚性杆B₁，第二旋转铰J₂连接第一刚性杆B₁与第二刚性杆B₂，第三旋转铰J₃连接第二刚性杆B₂与第三刚性杆B₃，第四旋转铰J₄连接第三刚性杆B₃与第二地面刚性支架；初始时刻，第一刚性杆B₁与第二刚性杆B₂水平布置，第三刚性杆B₃垂直于第二刚性杆B₂布置，第三刚性杆B₃位于第二刚性杆B₂下方；惯性坐标系为右手坐标系，选择水平为x轴，竖直向上为y轴，以第一旋转铰J₁处为原点o点，建立o-xyz坐标。

3个刚性杆的形状与密度完全一样，杆长都为2m，质量都为m＝1kg，惯量矩阵都为对角矩阵、主对角元素分别为I_xx＝1kg·m²、I_yy＝1kg·m²和I_zz＝1kg·m²，在第一刚性杆B₁的质心上施加一个初始状态为沿y轴向上、大小为10N的随体运动的力F，计算0～10s内系统的运动状态，依据上述本发明的具体实施方式，对该闭环多体系统进行分析。

步骤1、建立该闭环多体系统的动力学方程，依据建立的惯性坐标系，并采用相对坐标法，利用拉格朗日乘子法，根据虚功率原理，得到动力学方程：

由于本实施例中有4个z轴方向的旋转铰，因此上式中的运动参数q、和即为这4个旋转铰z轴的相对转角、相对转角速度和相对转角加速度，质量矩阵M(q,t)中包含3个刚性杆的质量与惯量信息，阻尼矩阵中包含相对运动关系，系统外力Q即为施加的力F。

步骤2、将动力学方程整理为含完整约束的动力学方程的一般形式：在动力学方程中添加Baumgarte约束稳定项，整合相应参数并整理方程，形成含完整约束的动力学方程的一般形式：

并进一步整理为

将系统的初始状态作为运动参数的初值输入，设定时间积分步长h与仿真计算总时间10s，并设置迭代收敛允许误差tol，开始进行仿真计算。

步骤3、从t＝0开始计算，依据式(6)～(11)，利用Bathe积分策略对含完整约束的动力学方程一般形式在时间步[t,t+h/2]上进行积分迭代求解，当满足收敛条件(即误差小于收敛允许误差tol时，式(12))得到t+h/2时刻的运动参数。

步骤4、依据式(13)～(17)，利用Bathe积分策略对动力学方程在时间步[t+h/2,t+h]上进行积分迭代求解，当满足收敛条件(即误差小于收敛允许误差tol时，式(18))，得到t+h时刻的运动参数。

步骤5、重复步骤3与步骤4，直到计算时间t达到设置的仿真总时间10s时，完成计算，输出4个旋转铰z轴方向的相对转角、相对转角速度和相对转角加速度随时间变化的值，即可分析该闭环多体系统在10s内的运动状态。

将本发明计算的结果与Adams软件计算的结果进行了对比，表1中列出了本发明与Adams软件分别计算的4个铰在10s时的角速度。

表1本发明与Adams软件分别计算的4个铰在10s时的角速度

从表1的数据中可以看出本发明的计算值与Adams软件的计算结果相吻合，本发明可以在较大积分步长的情况下也可以得到较好的精度。

为了验证雅克比矩阵的迭代初值中含有广义质量矩阵和广义阻尼矩阵后收敛效率提高，将上述模型计算的总迭代步数n与不含广义阻尼矩阵情况下的总迭代步数n_m进行了对比，结果如表2所示。

表2雅克比矩阵迭代初值不同时的总迭代步数对比

从表2中可以看出，雅克比矩阵中含有广义阻尼矩阵的情况下，迭代步数相对减少，收敛效率较高，并且随着积分步长的增大，收敛效率提高的比率也随之增大。鉴于Bathe积分策略具有在较大步长下也有较好精度的特点，可以在闭环多体动力学计算时选取适当的积分步长，这样就可以即精准又高效的进行求解。