CN107943748A - 基于Bathe积分策略的多体动力学方程求解方法 - Google Patents
基于Bathe积分策略的多体动力学方程求解方法 Download PDFInfo
- Publication number
- CN107943748A CN107943748A CN201711171910.6A CN201711171910A CN107943748A CN 107943748 A CN107943748 A CN 107943748A CN 201711171910 A CN201711171910 A CN 201711171910A CN 107943748 A CN107943748 A CN 107943748A
- Authority
- CN
- China
- Prior art keywords
- mrow
- mover
- mtd
- msub
- mtr
- Prior art date
- Legal status (The legal status is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the status listed.)
- Granted
Links
Classifications
-
- G—PHYSICS
- G06—COMPUTING; CALCULATING OR COUNTING
- G06F—ELECTRIC DIGITAL DATA PROCESSING
- G06F17/00—Digital computing or data processing equipment or methods, specially adapted for specific functions
- G06F17/10—Complex mathematical operations
- G06F17/11—Complex mathematical operations for solving equations, e.g. nonlinear equations, general mathematical optimization problems
- G06F17/13—Differential equations
Landscapes
- Engineering & Computer Science (AREA)
- Physics & Mathematics (AREA)
- General Physics & Mathematics (AREA)
- Mathematical Physics (AREA)
- Computational Mathematics (AREA)
- Mathematical Analysis (AREA)
- Mathematical Optimization (AREA)
- Pure & Applied Mathematics (AREA)
- Data Mining & Analysis (AREA)
- Theoretical Computer Science (AREA)
- Operations Research (AREA)
- Databases & Information Systems (AREA)
- Software Systems (AREA)
- General Engineering & Computer Science (AREA)
- Algebra (AREA)
- Management, Administration, Business Operations System, And Electronic Commerce (AREA)
- Complex Calculations (AREA)
- Other Investigation Or Analysis Of Materials By Electrical Means (AREA)
Abstract
本发明公开了一种基于Bathe积分策略的多体动力学方程求解方法，包括在动力学多体系统中建立相应的坐标系，并采用相对坐标法，利用拉格朗日乘子法，根据虚功率原理，得到动力学方程；将约束方程对时间求二阶导数，得到指标‑1的微分‑代数方程组，添加约束稳定项，形成含完整约束的动力学方程的一般形式；对含完整约束的动力学方程一般形式在时间步[t,t+h/2]上进行积分迭代求解，得到t+h/2时刻的运动参数；在时间步[t+h/2,t+h]上进行积分迭代求解，得到t+h时刻的运动参数；重复迭代计算，当计算时间达到设置好的仿真总时间，输出系统广义位移、速度和加速度随时间变化的值，即可分析过程闭环多体动力学系统在一定时间内的运动状态；本方法在较大积分步长的情况下也可以得到较好的精度。
Description
技术领域
本发明属于机械动力学领域，特别涉及到了一种基于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、令系统的广义位移初值和广义速度初值分别为P0和假设t时刻的广义位移、广义速度和广义加速度分别为Pt、和首先求解t+h/2时刻的广义位移Pt+h/2。
在[t,t+h/2]上利用梯形公式可推导得到t+h/2时刻的广义速度和广义加速度的表达式为：
将式(6)代入式(5)进行整理后可得到以广义位移Pt+h/2为未知变量的非线性代数方程：
3.2、利用逆Broyden拟牛顿法迭代求解式(7)
迭代过程如下：
式中yk和sk都为迭代的中间变量，Jt+h/2表示t+h/2时刻的雅克比矩阵，上标k和k+1分别表示第k次迭代和第k+1次迭代，Y的表达式如下：
根据泰勒展开原理，式(7)，式(8)中Pt+h/2和Jt+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时刻的广义位移Pt+h。
利用梯形公式，在[t+h/2,t+h]上可推导得到t+h时刻的广义速度和广义加速度的表达式：
将式(12)代入式(5)进行整理后可得到以广义位移Pt+h为未知变量的非线性代数方程：
4.2、利用逆Broyden拟牛顿法迭代求解式(14)
迭代过程如下：
式中Jt+h表示t+h时刻的雅克比矩阵。
根据泰勒展开原理和式(14)，式(15)中Pt+h和Jt+h迭代的初值可设置为：
当计算得到第k+1次迭代所得的和后，利用式(13)可得到和即可进行第k+2次的迭代。设置迭代收敛允许误差tol，并设置收敛条件：
当当前迭代步满足收敛条件式(18)时，迭代停止，计算得到的第k+1次的广义位移广义速度和广义加速度即为t+h时刻的运动参数。
步骤5、重复步骤3与步骤4，直到计算时间达到设置好的仿真总时间时，完成计算，输出系统广义位移、广义速度和广义加速度随时间变化的值，即可分析过程闭环多体动力学系统在一定时间内的运动状态。
实施例1：
通过下面实施例对本方法做进一步详细阐述。以下为本方法的一个仿真算例，算例为一个闭环多体系统问题，具体说明如下：
此处选择了一个含有3个刚性杆、4个转动铰的闭环多体系统，具体结构如图2所示，第一旋转铰J1连接第一地面刚性支架与第一刚性杆B1，第二旋转铰J2连接第一刚性杆B1与第二刚性杆B2，第三旋转铰J3连接第二刚性杆B2与第三刚性杆B3，第四旋转铰J4连接第三刚性杆B3与第二地面刚性支架；初始时刻，第一刚性杆B1与第二刚性杆B2水平布置，第三刚性杆B3垂直于第二刚性杆B2布置，第三刚性杆B3位于第二刚性杆B2下方；惯性坐标系为右手坐标系，选择水平为x轴，竖直向上为y轴，以第一旋转铰J1处为原点o点，建立o-xyz坐标。
3个刚性杆的形状与密度完全一样，杆长都为2m，质量都为m＝1kg，惯量矩阵都为对角矩阵、主对角元素分别为Ixx＝1kg·m2、Iyy＝1kg·m2和Izz＝1kg·m2，在第一刚性杆B1的质心上施加一个初始状态为沿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与不含广义阻尼矩阵情况下的总迭代步数nm进行了对比，结果如表2所示。
表2雅克比矩阵迭代初值不同时的总迭代步数对比
从表2中可以看出，雅克比矩阵中含有广义阻尼矩阵的情况下，迭代步数相对减少，收敛效率较高，并且随着积分步长的增大，收敛效率提高的比率也随之增大。鉴于Bathe积分策略具有在较大步长下也有较好精度的特点，可以在闭环多体动力学计算时选取适当的积分步长，这样就可以即精准又高效的进行求解。
Claims (7)
1.一种基于Bathe积分策略的多体动力学方程求解方法，其特征在于，包括以下步骤：
步骤1、构建含完整约束的动力学方程：在动力学多体系统中建立相应的坐标系，并采用相对坐标法，利用拉格朗日乘子法，根据虚功率原理，得到动力学方程；
步骤2、将动力学方程整理为含完整约束的动力学方程的一般形式：将约束方程对时间t求二阶导数，从而得到指标-1的微分-代数方程组，在动力学方程中添加Baumgarte约束稳定项，整合相应参数并整理方程，形成含完整约束的动力学方程的一般形式；
步骤3、利用Bathe积分策略对含完整约束的动力学方程一般形式在时间步[t，t+h/2]上进行积分迭代求解，得到t+h/2时刻的运动参数；其中h为积分步长；
步骤4、利用Bathe积分策略对含完整约束的动力学方程一般形式在时间步[t+h/2，t+h]上进行积分迭代求解，得到t+h时刻的运动参数；
步骤5、重复步骤3与步骤4，直到计算时间达到设置好的仿真总时间时，完成计算，输出系统广义位移、广义速度和广义加速度随时间变化的值，即可分析过程闭环多体动力学系统在一定时间内的运动状态。
2.如权利要求1所述的一种基于Bathe积分策略的多体动力学方程求解方法，其特征在于，步骤1构建含完整约束的动力学方程为：
<mrow>
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<mi>M</mi>
<mrow>
<mo>(</mo>
<mi>q</mi>
<mo>,</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mover>
<mi>q</mi>
<mo>&CenterDot;&CenterDot;</mo>
</mover>
<mo>+</mo>
<mi>C</mi>
<mrow>
<mo>(</mo>
<mi>q</mi>
<mo>,</mo>
<mover>
<mi>q</mi>
<mo>&CenterDot;</mo>
</mover>
<mo>,</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mover>
<mi>q</mi>
<mo>&CenterDot;</mo>
</mover>
<mo>+</mo>
<msubsup>
<mi>&Phi;</mi>
<mi>q</mi>
<mi>T</mi>
</msubsup>
<mi>&lambda;</mi>
<mo>=</mo>
<mi>Q</mi>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>&Phi;</mi>
<mrow>
<mo>(</mo>
<mi>q</mi>
<mo>,</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mn>0</mn>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</mrow>
式中t为时间，q、和分别为铰的相对位移、相对速度和相对加速度，Φ(q，t)＝0为系统约束方程，Φq为Φ(q，t)对q求导，()T表示对该矩阵的转置，λ为拉格朗日乘子，M(q，t)为质量矩阵，为阻尼矩阵，Q为系统外力。
3.如权利要求2所述的一种基于Bathe积分策略的多体动力学方程求解方法，其特征在于，步骤2建立含完整约束的动力学方程的一般形式，具体包括以下步骤：
步骤2.1、对式(1)中的约束方程Φ(q，t)＝0对时间t求二阶导数，可以将式(1)写成如下指标-1的微分-代数方程组：
<mrow>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mrow>
<mi>M</mi>
<mrow>
<mo>(</mo>
<mi>q</mi>
<mo>,</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
<mtd>
<msubsup>
<mi>&Phi;</mi>
<mi>q</mi>
<mi>T</mi>
</msubsup>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>&Phi;</mi>
<mi>q</mi>
</msub>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
</mtable>
</mfenced>
<mfenced open = "{" close = "}">
<mtable>
<mtr>
<mtd>
<mover>
<mi>q</mi>
<mo>&CenterDot;&CenterDot;</mo>
</mover>
</mtd>
</mtr>
<mtr>
<mtd>
<mi>&lambda;</mi>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>=</mo>
<mfenced open = "{" close = "}">
<mtable>
<mtr>
<mtd>
<mrow>
<mi>Q</mi>
<mo>-</mo>
<mi>C</mi>
<mrow>
<mo>(</mo>
<mi>q</mi>
<mo>,</mo>
<mover>
<mi>q</mi>
<mo>&CenterDot;</mo>
</mover>
<mo>,</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
<mover>
<mi>q</mi>
<mo>&CenterDot;</mo>
</mover>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mi>&Psi;</mi>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>2</mn>
<mo>)</mo>
</mrow>
</mrow>
式中其中Φqq为Φ(q，t)对q求二阶导数，Φqt为Φ(q，t)依次对q和t求一阶导数，Φtt为Φ(q，t)对t求二阶导数；
步骤2.2、在动力学方程中添加Baumgarte约束稳定项：
根据Baumgarte约束稳定法，在式(2)中加入稳定项并移项，将式(2)写成如下形式
<mrow>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mrow>
<mi>M</mi>
<mrow>
<mo>(</mo>
<mi>q</mi>
<mo>,</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
<mtd>
<msubsup>
<mi>&Phi;</mi>
<mi>q</mi>
<mi>T</mi>
</msubsup>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>&Phi;</mi>
<mi>q</mi>
</msub>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
</mtable>
</mfenced>
<mfenced open = "{" close = "}">
<mtable>
<mtr>
<mtd>
<mover>
<mi>q</mi>
<mo>&CenterDot;&CenterDot;</mo>
</mover>
</mtd>
</mtr>
<mtr>
<mtd>
<mi>&lambda;</mi>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>+</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mrow>
<mi>C</mi>
<mrow>
<mo>(</mo>
<mi>q</mi>
<mo>,</mo>
<mover>
<mi>q</mi>
<mo>&CenterDot;</mo>
</mover>
<mo>,</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
</mtable>
</mfenced>
<mfenced open = "{" close = "}">
<mtable>
<mtr>
<mtd>
<mover>
<mi>q</mi>
<mo>&CenterDot;</mo>
</mover>
</mtd>
</mtr>
<mtr>
<mtd>
<msup>
<mi>&lambda;</mi>
<mo>&prime;</mo>
</msup>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>=</mo>
<mfenced open = "{" close = "}">
<mtable>
<mtr>
<mtd>
<mi>Q</mi>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>&Psi;</mi>
<mo>-</mo>
<mn>2</mn>
<mi>&alpha;</mi>
<mover>
<mi>&Phi;</mi>
<mo>&CenterDot;</mo>
</mover>
<mo>-</mo>
<msup>
<mi>&beta;</mi>
<mn>2</mn>
</msup>
<mi>&Phi;</mi>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>3</mn>
<mo>)</mo>
</mrow>
</mrow>
式中，为速度约束，α和β为Baumgarte系数，取h为积分步长，λ′和λ″分别作为λ对时间的一次积分项和二次积分项；
步骤2.3、整合相应参数并整理方程，形成含完整约束的动力学方程的一般形式
令
<mrow>
<mtable>
<mtr>
<mtd>
<mrow>
<mi>P</mi>
<mo>=</mo>
<mfenced open = "{" close = "}">
<mtable>
<mtr>
<mtd>
<mi>q</mi>
</mtd>
</mtr>
<mtr>
<mtd>
<msup>
<mi>&lambda;</mi>
<mrow>
<mo>&prime;</mo>
<mo>&prime;</mo>
</mrow>
</msup>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>,</mo>
<mover>
<mi>P</mi>
<mo>&CenterDot;</mo>
</mover>
<mo>=</mo>
<mfenced open = "{" close = "}">
<mtable>
<mtr>
<mtd>
<mover>
<mi>q</mi>
<mo>&CenterDot;</mo>
</mover>
</mtd>
</mtr>
<mtr>
<mtd>
<msup>
<mi>&lambda;</mi>
<mo>&prime;</mo>
</msup>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>,</mo>
<mover>
<mi>P</mi>
<mo>&CenterDot;&CenterDot;</mo>
</mover>
<mo>=</mo>
<mfenced open = "{" close = "}">
<mtable>
<mtr>
<mtd>
<mover>
<mi>q</mi>
<mo>&CenterDot;&CenterDot;</mo>
</mover>
</mtd>
</mtr>
<mtr>
<mtd>
<mi>&lambda;</mi>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>M</mi>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mrow>
<mi>M</mi>
<mrow>
<mo>(</mo>
<mi>q</mi>
<mo>,</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
<mtd>
<msubsup>
<mi>&Phi;</mi>
<mi>q</mi>
<mi>T</mi>
</msubsup>
</mtd>
</mtr>
<mtr>
<mtd>
<msub>
<mi>&Phi;</mi>
<mi>q</mi>
</msub>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>,</mo>
<mover>
<mi>C</mi>
<mo>&OverBar;</mo>
</mover>
<mo>=</mo>
<mfenced open = "[" close = "]">
<mtable>
<mtr>
<mtd>
<mrow>
<mi>C</mi>
<mrow>
<mo>(</mo>
<mi>q</mi>
<mo>,</mo>
<mover>
<mi>q</mi>
<mo>&CenterDot;</mo>
</mover>
<mo>,</mo>
<mi>t</mi>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
<mtr>
<mtd>
<mn>0</mn>
</mtd>
<mtd>
<mn>0</mn>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>,</mo>
<mi>Q</mi>
<mo>=</mo>
<mfenced open = "{" close = "}">
<mtable>
<mtr>
<mtd>
<mi>Q</mi>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mi>&Psi;</mi>
<mo>-</mo>
<mn>2</mn>
<mi>&alpha;</mi>
<mover>
<mi>&Phi;</mi>
<mo>&CenterDot;</mo>
</mover>
<mo>-</mo>
<msup>
<mi>&beta;</mi>
<mn>2</mn>
</msup>
<mi>&Phi;</mi>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
</mrow>
</mtd>
</mtr>
</mtable>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>4</mn>
<mo>)</mo>
</mrow>
</mrow>
式中P、和分别为动力学方程的广义位移项、广义速度项和广义加速度项，和分别为动力学方程的广义质量矩阵和广义阻尼矩阵，为系统广义力；
则式(3)整理为：
<mrow>
<mover>
<mi>M</mi>
<mo>&OverBar;</mo>
</mover>
<mover>
<mi>P</mi>
<mo>&CenterDot;&CenterDot;</mo>
</mover>
<mo>+</mo>
<mover>
<mi>C</mi>
<mo>&OverBar;</mo>
</mover>
<mover>
<mi>P</mi>
<mo>&CenterDot;</mo>
</mover>
<mo>=</mo>
<mover>
<mi>Q</mi>
<mo>&OverBar;</mo>
</mover>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>5</mn>
<mo>)</mo>
</mrow>
<mo>.</mo>
</mrow>
4.如权利要求3所述的一种基于Bathe积分策略的多体动力学方程求解方法，其特征在于，步骤3得到t+h/2时刻的运动参数具体包括以下步骤：
步骤3.1、令系统的广义位移初值和广义速度初值分别为P0和假设t时刻的广义位移、广义速度和广义加速度分别为Pt、和
在[t，t+h/2]上利用梯形公式推导得到t+h/2时刻的广义速度和广义加速度的表达式为：
<mrow>
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<msub>
<mover>
<mi>P</mi>
<mo>&CenterDot;</mo>
</mover>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
</msub>
<mo>=</mo>
<mfrac>
<mn>4</mn>
<mi>h</mi>
</mfrac>
<mrow>
<mo>(</mo>
<msub>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
</msub>
<mo>-</mo>
<msub>
<mi>P</mi>
<mi>t</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>-</mo>
<msub>
<mover>
<mi>P</mi>
<mo>&CenterDot;</mo>
</mover>
<mi>t</mi>
</msub>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mover>
<mi>P</mi>
<mo>&CenterDot;&CenterDot;</mo>
</mover>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
</msub>
<mo>=</mo>
<mfrac>
<mn>16</mn>
<msup>
<mi>h</mi>
<mn>2</mn>
</msup>
</mfrac>
<mrow>
<mo>(</mo>
<msub>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
</msub>
<mo>-</mo>
<msub>
<mi>P</mi>
<mi>t</mi>
</msub>
<mo>-</mo>
<mfrac>
<mi>h</mi>
<mn>2</mn>
</mfrac>
<msub>
<mover>
<mi>P</mi>
<mo>&CenterDot;</mo>
</mover>
<mi>t</mi>
</msub>
<mo>)</mo>
</mrow>
<mo>-</mo>
<msub>
<mover>
<mi>P</mi>
<mo>&CenterDot;&CenterDot;</mo>
</mover>
<mi>t</mi>
</msub>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>6</mn>
<mo>)</mo>
</mrow>
</mrow>
将式(6)代入式(5)进行整理后可得到以广义位移Pt+h/2为未知变量的非线性代数方程：
<mrow>
<mo>(</mo>
<mfrac>
<mn>16</mn>
<msup>
<mi>h</mi>
<mn>2</mn>
</msup>
</mfrac>
<mover>
<mi>M</mi>
<mo>&OverBar;</mo>
</mover>
<mo>+</mo>
<mfrac>
<mn>4</mn>
<mi>h</mi>
</mfrac>
<mover>
<mi>C</mi>
<mo>&OverBar;</mo>
</mover>
<mo>)</mo>
<msub>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
</msub>
<mo>-</mo>
<mover>
<mi>M</mi>
<mo>&OverBar;</mo>
</mover>
<mo>(</mo>
<mfrac>
<mn>16</mn>
<msup>
<mi>h</mi>
<mn>2</mn>
</msup>
</mfrac>
<msub>
<mi>P</mi>
<mi>t</mi>
</msub>
<mo>+</mo>
<mfrac>
<mn>8</mn>
<mi>h</mi>
</mfrac>
<msub>
<mover>
<mi>P</mi>
<mo>&CenterDot;</mo>
</mover>
<mi>t</mi>
</msub>
<mo>+</mo>
<msub>
<mover>
<mi>P</mi>
<mo>&CenterDot;&CenterDot;</mo>
</mover>
<mi>t</mi>
</msub>
<mo>)</mo>
<mo>-</mo>
<mover>
<mi>C</mi>
<mo>&OverBar;</mo>
</mover>
<mo>(</mo>
<mfrac>
<mn>4</mn>
<mi>h</mi>
</mfrac>
<msub>
<mi>P</mi>
<mi>t</mi>
</msub>
<mo>+</mo>
<msub>
<mover>
<mi>P</mi>
<mo>&CenterDot;</mo>
</mover>
<mi>t</mi>
</msub>
<mo>)</mo>
<mo>=</mo>
<mover>
<mi>Q</mi>
<mo>&OverBar;</mo>
</mover>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>7</mn>
<mo>)</mo>
</mrow>
</mrow>
步骤3.2、利用逆Broyden拟牛顿法迭代求解式(7)
<mrow>
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<msubsup>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</msubsup>
<mo>=</mo>
<msubsup>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mo>-</mo>
<msup>
<mrow>
<mo>&lsqb;</mo>
<msubsup>
<mi>J</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mo>&rsqb;</mo>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mi>Y</mi>
<mrow>
<mo>(</mo>
<msubsup>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>y</mi>
<mi>k</mi>
</msub>
<mo>=</mo>
<mi>Y</mi>
<mrow>
<mo>(</mo>
<msubsup>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</msubsup>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mi>Y</mi>
<mrow>
<mo>(</mo>
<msubsup>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>s</mi>
<mi>k</mi>
</msub>
<mo>=</mo>
<msubsup>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</msubsup>
<mo>-</mo>
<msubsup>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msubsup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msup>
<mrow>
<mo>&lsqb;</mo>
<msubsup>
<mi>J</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</msubsup>
<mo>&rsqb;</mo>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>=</mo>
<msup>
<mrow>
<mo>&lsqb;</mo>
<msubsup>
<mi>J</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mo>&rsqb;</mo>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>+</mo>
<mo>&lsqb;</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>s</mi>
<mi>k</mi>
</msub>
<msup>
<mrow>
<mo>&lsqb;</mo>
<msubsup>
<mi>J</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mo>&rsqb;</mo>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<msub>
<mi>y</mi>
<mi>k</mi>
</msub>
<mo>)</mo>
</mrow>
<msubsup>
<mi>s</mi>
<mi>k</mi>
<mi>T</mi>
</msubsup>
<mo>&rsqb;</mo>
<msup>
<mrow>
<mo>&lsqb;</mo>
<msubsup>
<mi>J</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mo>&rsqb;</mo>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>/</mo>
<mrow>
<mo>(</mo>
<msubsup>
<mi>s</mi>
<mi>k</mi>
<mi>T</mi>
</msubsup>
<msup>
<mrow>
<mo>&lsqb;</mo>
<msubsup>
<mi>J</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mo>&rsqb;</mo>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<msub>
<mi>y</mi>
<mi>k</mi>
</msub>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>8</mn>
<mo>)</mo>
</mrow>
</mrow>
式中yk和sk都为迭代的中间变量，Jt+h/2表示t+h/2时刻的雅克比矩阵，上标k和k+1分别表示第k次迭代和第k+1次迭代，Y的表达式如下：
<mrow>
<mi>Y</mi>
<mo>=</mo>
<mover>
<mi>M</mi>
<mo>&OverBar;</mo>
</mover>
<mover>
<mi>P</mi>
<mo>&CenterDot;&CenterDot;</mo>
</mover>
<mo>+</mo>
<mover>
<mi>C</mi>
<mo>&OverBar;</mo>
</mover>
<mover>
<mi>P</mi>
<mo>&CenterDot;</mo>
</mover>
<mo>-</mo>
<mover>
<mi>Q</mi>
<mo>&OverBar;</mo>
</mover>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>9</mn>
<mo>)</mo>
</mrow>
</mrow>
根据泰勒展开原理，式(7)，式(8)中Pt+h/2和Jt+h/2的迭代初值分别设置为：
<mrow>
<msubsup>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
<mn>0</mn>
</msubsup>
<mo>=</mo>
<msub>
<mi>P</mi>
<mi>t</mi>
</msub>
<mo>+</mo>
<mfrac>
<mi>h</mi>
<mn>2</mn>
</mfrac>
<msub>
<mover>
<mi>P</mi>
<mo>&CenterDot;</mo>
</mover>
<mi>t</mi>
</msub>
<mo>+</mo>
<mfrac>
<msup>
<mi>h</mi>
<mn>2</mn>
</msup>
<mn>16</mn>
</mfrac>
<msub>
<mover>
<mi>P</mi>
<mo>&CenterDot;&CenterDot;</mo>
</mover>
<mi>t</mi>
</msub>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>10</mn>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<msubsup>
<mi>J</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
<mn>0</mn>
</msubsup>
<mo>=</mo>
<mfrac>
<mn>16</mn>
<msup>
<mi>h</mi>
<mn>2</mn>
</msup>
</mfrac>
<mover>
<mi>M</mi>
<mo>&OverBar;</mo>
</mover>
<mo>+</mo>
<mfrac>
<mn>4</mn>
<mi>h</mi>
</mfrac>
<mover>
<mi>C</mi>
<mo>&OverBar;</mo>
</mover>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>11</mn>
<mo>)</mo>
</mrow>
</mrow>
当计算得到第k+1次迭代所得的和后，利用式(6)可得到和即可进行第k+2次的迭代；
设置收敛条件，计算得到的第k+1次的广义位移广义速度和广义加速度即为t+h/2时刻的运动参数。
5.如权利要求4所述的一种基于Bathe积分策略的多体动力学方程求解方法，其特征在于，设置收敛条件为：
<mrow>
<mfrac>
<mrow>
<mo>|</mo>
<mrow>
<msubsup>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</msubsup>
<mo>-</mo>
<msubsup>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msubsup>
</mrow>
<mo>|</mo>
</mrow>
<mrow>
<mo>|</mo>
<msubsup>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</msubsup>
<mo>|</mo>
</mrow>
</mfrac>
<mo><</mo>
<mi>t</mi>
<mi>o</mi>
<mi>l</mi>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>12</mn>
<mo>)</mo>
</mrow>
</mrow>
tol为迭代收敛允许误差。
6.如权利要求4或5所述的一种基于Bathe积分策略的多体动力学方程求解方法，其特征在于，步骤4得到t+h时刻的运动参数，具体包括以下步骤：
步骤4.1、求解t+h时刻的广义位移Pt+h：
利用梯形公式，在[t+h/2，t+h]上可推导得到t+h时刻的广义速度和广义加速度的表达式：
<mrow>
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<msub>
<mover>
<mi>P</mi>
<mo>&CenterDot;</mo>
</mover>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
</mrow>
</msub>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mi>h</mi>
</mfrac>
<msub>
<mi>P</mi>
<mi>t</mi>
</msub>
<mo>-</mo>
<mfrac>
<mn>4</mn>
<mi>h</mi>
</mfrac>
<msub>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
</msub>
<mo>+</mo>
<mfrac>
<mn>3</mn>
<mi>h</mi>
</mfrac>
<msub>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
</mrow>
</msub>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mover>
<mi>P</mi>
<mo>&CenterDot;&CenterDot;</mo>
</mover>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
</mrow>
</msub>
<mo>=</mo>
<mfrac>
<mn>1</mn>
<mi>h</mi>
</mfrac>
<msub>
<mover>
<mi>P</mi>
<mo>&CenterDot;</mo>
</mover>
<mi>t</mi>
</msub>
<mo>-</mo>
<mfrac>
<mn>4</mn>
<mi>h</mi>
</mfrac>
<msub>
<mover>
<mi>P</mi>
<mo>&CenterDot;</mo>
</mover>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
</msub>
<mo>+</mo>
<mfrac>
<mn>3</mn>
<mi>h</mi>
</mfrac>
<msub>
<mover>
<mi>P</mi>
<mo>&CenterDot;</mo>
</mover>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
</mrow>
</msub>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>13</mn>
<mo>)</mo>
</mrow>
</mrow>
将式(12)代入式(5)进行整理后可得到以广义位移Pt+h为未知变量的非线性代数方程：
<mrow>
<mtable>
<mtr>
<mtd>
<mrow>
<mo>(</mo>
<mfrac>
<mn>9</mn>
<msup>
<mi>h</mi>
<mn>2</mn>
</msup>
</mfrac>
<mover>
<mi>M</mi>
<mo>&OverBar;</mo>
</mover>
<mo>+</mo>
<mfrac>
<mn>3</mn>
<mi>h</mi>
</mfrac>
<mover>
<mi>C</mi>
<mo>&OverBar;</mo>
</mover>
<mo>)</mo>
<msub>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
</mrow>
</msub>
<mo>+</mo>
<mover>
<mi>M</mi>
<mo>&OverBar;</mo>
</mover>
<mo>(</mo>
<mfrac>
<mn>1</mn>
<mi>h</mi>
</mfrac>
<msub>
<mover>
<mi>P</mi>
<mo>&CenterDot;</mo>
</mover>
<mi>t</mi>
</msub>
<mo>-</mo>
<mfrac>
<mn>4</mn>
<mi>h</mi>
</mfrac>
<msub>
<mover>
<mi>P</mi>
<mo>&CenterDot;</mo>
</mover>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
</msub>
<mo>+</mo>
<mfrac>
<mn>3</mn>
<msup>
<mi>h</mi>
<mn>2</mn>
</msup>
</mfrac>
<msub>
<mi>P</mi>
<mi>t</mi>
</msub>
<mo>-</mo>
<mfrac>
<mn>12</mn>
<msup>
<mi>h</mi>
<mn>2</mn>
</msup>
</mfrac>
<msub>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
</msub>
<mo>)</mo>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<mo>+</mo>
<mover>
<mi>C</mi>
<mo>&OverBar;</mo>
</mover>
<mrow>
<mo>(</mo>
<mfrac>
<mn>1</mn>
<mi>h</mi>
</mfrac>
<msub>
<mi>P</mi>
<mi>t</mi>
</msub>
<mo>-</mo>
<mfrac>
<mn>4</mn>
<mi>h</mi>
</mfrac>
<msub>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
</msub>
<mo>)</mo>
</mrow>
<mo>=</mo>
<mover>
<mi>Q</mi>
<mo>&OverBar;</mo>
</mover>
</mrow>
</mtd>
</mtr>
</mtable>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>14</mn>
<mo>)</mo>
</mrow>
</mrow>
步骤4.2、利用逆Broyden拟牛顿法迭代求解式(14)
<mrow>
<mfenced open = "{" close = "">
<mtable>
<mtr>
<mtd>
<mrow>
<msubsup>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</msubsup>
<mo>=</mo>
<msubsup>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mo>-</mo>
<msup>
<mrow>
<mo>&lsqb;</mo>
<msubsup>
<mi>J</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mo>&rsqb;</mo>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mi>Y</mi>
<mrow>
<mo>(</mo>
<msubsup>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>y</mi>
<mi>k</mi>
</msub>
<mo>=</mo>
<mi>Y</mi>
<mrow>
<mo>(</mo>
<msubsup>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</msubsup>
<mo>)</mo>
</mrow>
<mo>-</mo>
<mi>Y</mi>
<mrow>
<mo>(</mo>
<msubsup>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msub>
<mi>s</mi>
<mi>k</mi>
</msub>
<mo>=</mo>
<msubsup>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</msubsup>
<mo>-</mo>
<msubsup>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msubsup>
</mrow>
</mtd>
</mtr>
<mtr>
<mtd>
<mrow>
<msup>
<mrow>
<mo>&lsqb;</mo>
<msubsup>
<mi>J</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</msubsup>
<mo>&rsqb;</mo>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>=</mo>
<msup>
<mrow>
<mo>&lsqb;</mo>
<msubsup>
<mi>J</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mo>&rsqb;</mo>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>+</mo>
<mo>&lsqb;</mo>
<mrow>
<mo>(</mo>
<msub>
<mi>s</mi>
<mi>k</mi>
</msub>
<mo>-</mo>
<msup>
<mrow>
<mo>&lsqb;</mo>
<msubsup>
<mi>J</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mo>&rsqb;</mo>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<msub>
<mi>y</mi>
<mi>k</mi>
</msub>
<mo>)</mo>
</mrow>
<msubsup>
<mi>s</mi>
<mi>k</mi>
<mi>T</mi>
</msubsup>
<mo>&rsqb;</mo>
<msup>
<mrow>
<mo>&lsqb;</mo>
<msubsup>
<mi>J</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mo>&rsqb;</mo>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<mo>/</mo>
<mrow>
<mo>(</mo>
<msubsup>
<mi>s</mi>
<mi>k</mi>
<mi>T</mi>
</msubsup>
<msup>
<mrow>
<mo>&lsqb;</mo>
<msubsup>
<mi>J</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msubsup>
<mo>&rsqb;</mo>
</mrow>
<mrow>
<mo>-</mo>
<mn>1</mn>
</mrow>
</msup>
<msub>
<mi>y</mi>
<mi>k</mi>
</msub>
<mo>)</mo>
</mrow>
</mrow>
</mtd>
</mtr>
</mtable>
</mfenced>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>15</mn>
<mo>)</mo>
</mrow>
</mrow>
式中Jt+h表示t+h时刻的雅克比矩阵；
根据泰勒展开原理和式(14)，式(15)中Pt+h和Jt+h迭代的初值设置为：
<mrow>
<msubsup>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
</mrow>
<mn>0</mn>
</msubsup>
<mo>=</mo>
<msub>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
</msub>
<mo>+</mo>
<mfrac>
<mn>2</mn>
<mi>h</mi>
</mfrac>
<msub>
<mover>
<mi>P</mi>
<mo>&CenterDot;</mo>
</mover>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
</msub>
<mo>+</mo>
<mfrac>
<msup>
<mi>h</mi>
<mn>2</mn>
</msup>
<mn>16</mn>
</mfrac>
<msub>
<mover>
<mi>P</mi>
<mo>&CenterDot;&CenterDot;</mo>
</mover>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
<mo>/</mo>
<mn>2</mn>
</mrow>
</msub>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>16</mn>
<mo>)</mo>
</mrow>
</mrow>
<mrow>
<msubsup>
<mi>J</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
</mrow>
<mn>0</mn>
</msubsup>
<mo>=</mo>
<mfrac>
<mn>9</mn>
<msup>
<mi>h</mi>
<mn>2</mn>
</msup>
</mfrac>
<mover>
<mi>M</mi>
<mo>&OverBar;</mo>
</mover>
<mo>+</mo>
<mfrac>
<mn>3</mn>
<mi>h</mi>
</mfrac>
<mover>
<mi>C</mi>
<mo>&OverBar;</mo>
</mover>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>17</mn>
<mo>)</mo>
</mrow>
</mrow>
当计算得到第k+1次迭代所得的和后，利用式(13)得到和即可进行第k+2次的迭代；设置收敛条件，计算得到的第k+1次的广义位移广义速度和广义加速度即为t+h时刻的运动参数。
7.如权利要求6所述的一种基于Bathe积分策略的多体动力学方程求解方法，其特征在于，收敛条件为：
<mrow>
<mfrac>
<mrow>
<mo>|</mo>
<mrow>
<msubsup>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</msubsup>
<mo>-</mo>
<msubsup>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>)</mo>
</mrow>
</msubsup>
</mrow>
<mo>|</mo>
</mrow>
<mrow>
<mo>|</mo>
<msubsup>
<mi>P</mi>
<mrow>
<mi>t</mi>
<mo>+</mo>
<mi>h</mi>
</mrow>
<mrow>
<mo>(</mo>
<mi>k</mi>
<mo>+</mo>
<mn>1</mn>
<mo>)</mo>
</mrow>
</msubsup>
<mo>|</mo>
</mrow>
</mfrac>
<mo><</mo>
<mi>t</mi>
<mi>o</mi>
<mi>l</mi>
<mo>-</mo>
<mo>-</mo>
<mo>-</mo>
<mrow>
<mo>(</mo>
<mn>18</mn>
<mo>)</mo>
</mrow>
<mo>.</mo>
</mrow>
Priority Applications (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711171910.6A CN107943748B (zh) | 2017-11-22 | 2017-11-22 | 基于Bathe积分策略获取闭环多体系统运动状态的方法 |
Applications Claiming Priority (1)
| Application Number | Priority Date | Filing Date | Title |
|---|---|---|---|
| CN201711171910.6A CN107943748B (zh) | 2017-11-22 | 2017-11-22 | 基于Bathe积分策略获取闭环多体系统运动状态的方法 |
Publications (2)
| Publication Number | Publication Date |
|---|---|
| CN107943748A true CN107943748A (zh) | 2018-04-20 |
| CN107943748B CN107943748B (zh) | 2019-04-12 |
Family
ID=61929727
Family Applications (1)
| Application Number | Title | Priority Date | Filing Date |
|---|---|---|---|
| CN201711171910.6A Active CN107943748B (zh) | 2017-11-22 | 2017-11-22 | 基于Bathe积分策略获取闭环多体系统运动状态的方法 |
Country Status (1)
| Country | Link |
|---|---|
| CN (1) | CN107943748B (zh) |
Cited By (6)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108694280A (zh) * | 2018-05-14 | 2018-10-23 | 电子科技大学 | 基于新型杂交应力四面体单元的冲击响应仿真模拟方法 |
| CN109543264A (zh) * | 2018-11-12 | 2019-03-29 | 天津理工大学 | 一种基于多维度重构校正的柔性多体机器人建模与求解方法 |
| CN111159636A (zh) * | 2019-12-04 | 2020-05-15 | 大连理工大学 | 基于绝对节点坐标描述的柔性多体系统动力学半解析灵敏度分析方法 |
| CN111488550A (zh) * | 2020-03-18 | 2020-08-04 | 华东理工大学 | 一种基于时间积分和牛顿法自动联用高效求解稳态微观动力学方程组的方法 |
| CN112560364A (zh) * | 2020-12-22 | 2021-03-26 | 长江航道规划设计研究院 | 一种内河浮标锚链多体系统最小势能求解方法 |
| CN117852322A (zh) * | 2024-03-08 | 2024-04-09 | 西北工业大学深圳研究院 | 一种基于虚功率原理的变体飞行器动力学建模方法及装置 |
Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20150149132A1 (en) * | 2013-11-27 | 2015-05-28 | Battelle Memorial Institute | Time-stacking method for dynamic simulations |
| CN107220421A (zh) * | 2017-05-18 | 2017-09-29 | 北京理工大学 | 一种空间复杂柔性结构多体系统动力学建模与计算方法 |
-
2017
- 2017-11-22 CN CN201711171910.6A patent/CN107943748B/zh active Active
Patent Citations (2)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| US20150149132A1 (en) * | 2013-11-27 | 2015-05-28 | Battelle Memorial Institute | Time-stacking method for dynamic simulations |
| CN107220421A (zh) * | 2017-05-18 | 2017-09-29 | 北京理工大学 | 一种空间复杂柔性结构多体系统动力学建模与计算方法 |
Non-Patent Citations (1)
| Title |
|---|
| 许永生: "闭环多体系统动力学递推建模及其求解方法", 《中国博士学位论文全文数据库-基础科学辑》 * |
Cited By (11)
| Publication number | Priority date | Publication date | Assignee | Title |
|---|---|---|---|---|
| CN108694280A (zh) * | 2018-05-14 | 2018-10-23 | 电子科技大学 | 基于新型杂交应力四面体单元的冲击响应仿真模拟方法 |
| CN109543264A (zh) * | 2018-11-12 | 2019-03-29 | 天津理工大学 | 一种基于多维度重构校正的柔性多体机器人建模与求解方法 |
| CN109543264B (zh) * | 2018-11-12 | 2023-01-20 | 天津理工大学 | 一种基于多维度重构校正的柔性多体机器人建模与求解方法 |
| CN111159636A (zh) * | 2019-12-04 | 2020-05-15 | 大连理工大学 | 基于绝对节点坐标描述的柔性多体系统动力学半解析灵敏度分析方法 |
| CN111159636B (zh) * | 2019-12-04 | 2021-09-24 | 大连理工大学 | 一种柔性多体系统动力学半解析灵敏度分析方法 |
| CN111488550A (zh) * | 2020-03-18 | 2020-08-04 | 华东理工大学 | 一种基于时间积分和牛顿法自动联用高效求解稳态微观动力学方程组的方法 |
| CN111488550B (zh) * | 2020-03-18 | 2023-06-16 | 华东理工大学 | 一种基于时间积分和牛顿法自动联用高效求解稳态微观动力学方程组的方法 |
| CN112560364A (zh) * | 2020-12-22 | 2021-03-26 | 长江航道规划设计研究院 | 一种内河浮标锚链多体系统最小势能求解方法 |
| CN112560364B (zh) * | 2020-12-22 | 2023-10-10 | 长江航道规划设计研究院 | 一种内河浮标锚链多体系统最小势能求解方法 |
| CN117852322A (zh) * | 2024-03-08 | 2024-04-09 | 西北工业大学深圳研究院 | 一种基于虚功率原理的变体飞行器动力学建模方法及装置 |
| CN117852322B (zh) * | 2024-03-08 | 2024-05-10 | 西北工业大学深圳研究院 | 一种基于虚功率原理的变体飞行器动力学建模方法及装置 |
Also Published As
| Publication number | Publication date |
|---|---|
| CN107943748B (zh) | 2019-04-12 |
Similar Documents
| Publication | Publication Date | Title |
|---|---|---|
| CN107943748A (zh) | 基于Bathe积分策略的多体动力学方程求解方法 | |
| Brüls et al. | Lie group generalized-α time integration of constrained flexible multibody systems | |
| Terze et al. | Lie-group integration method for constrained multibody systems in state space | |
| Mian et al. | Numerical investigation of structural geometric nonlinearity effect in high-aspect-ratio wing using CFD/CSD coupled approach | |
| Morton et al. | Hopf-bifurcation analysis of airfoil flutter at transonic speeds | |
| CN105574809B (zh) | 基于矩阵指数的电磁暂态仿真图形处理器并行计算方法 | |
| Beran et al. | Reduced-order modelling of limit-cycle oscillation for aeroelastic systems | |
| Zupan et al. | On conservation of energy and kinematic compatibility in dynamics of nonlinear velocity-based three-dimensional beams | |
| Triantafyllou et al. | Small and large displacement dynamic analysis of frame structures based on hysteretic beam elements | |
| Arnold | DAE aspects of multibody system dynamics | |
| Han et al. | Data‐guided model predictive control based on smoothed contact dynamics | |
| Kotikalpudi | Robust flutter analysis for aeroservoelastic systems | |
| Haug | Multibody dynamics on differentiable manifolds | |
| Howcroft et al. | On the geometrically exact low-order modelling of a flexible beam: formulation and numerical tests | |
| Kim et al. | Unified mechanism synthesis method of a planar four-bar linkage for path generation employing a spring-connected arbitrarily sized rectangular block model | |
| Brugnoli et al. | Port-Hamiltonian flexible multibody dynamics | |
| Zhang et al. | Mathematical modeling and dynamic characteristic analysis of a novel parallel tracking mechanism for inter-satellite link antenna | |
| Chen et al. | Development of an object-oriented framework for the vibration characteristic computation of multibody systems | |
| CN117171864B (zh) | 一种梁结构线性振动预测方法 | |
| Hjelmstad et al. | Variational basis of nonlinear flexibility methods for structural analysis of frames | |
| Valente et al. | OpenFSI interface for strongly coupled steady and unsteady aeroelasticity | |
| Williams et al. | RPI-MATLAB-Simulator: A Tool for Efficient Research and Practical Teaching in Multibody Dynamics. | |
| Zupan et al. | Velocity-based approach in non-linear dynamics of three-dimensional beams with enforced kinematic compatibility | |
| CN112084592A (zh) | 折叠式桁架动力学分析系统、方法、装置和存储介质 | |
| Liang et al. | Symbolic integration of multibody system dynamics with the finite element method |
Legal Events
| Date | Code | Title | Description |
|---|---|---|---|
| PB01 | Publication | ||
| PB01 | Publication | ||
| SE01 | Entry into force of request for substantive examination | ||
| SE01 | Entry into force of request for substantive examination | ||
| GR01 | Patent grant | ||
| GR01 | Patent grant |