多体系统动力学基本理论

多体系统动力学基本理论 第2章 多体系统动力学基本理论 本章主要介绍多体系统动力学的基本理论，包括多刚体系统动力学建模、多柔体系统 动力学建模、多体系统动力学方程求解及多体系统动力学中的刚性（Stiff）问。通过本章的学习可以对多体系统动力学的基本理论有较深入的了解，为具体软件的学习打下良好 的理论基础。 2.1 多体系统动力学研究状况 多体系统动力学的核心问题是建模和求解问题，其系统研究开始于20世纪60年代。从60年代到80年代，侧重于多刚体系统的研究，主要是研究多刚体系统的自动建模和数 值求解；到了80年代中期，多刚体系统动力学的研究已经取得一系列成果，尤其是建模理 论趋于成熟，但更稳定、更有效的数值求解方法仍然是研究的热点；80年代之后，多体系统动力学的研究更偏重于多柔体系统动力学，这个领域也正式被称为计算多体系统动力学， 它至今仍然是力学研究中最有活力的分支之一，但已经远远地超过一般力学的涵义。 本节将叙述多体系统动力学发展的历史和目前国内外研究的现状。 2.1.1 多体系统动力学研究的发展 机械系统动力学与仿真是随着计算机技术的发展而不断成熟的，多体系统动力学 是其理论基础。计算机技术自其诞生以来，渗透到了科学计算和工程应用的几乎每一个领 域。数值分析技术与传统力学的结合曾在结构力学领域取得了辉煌的成就，出现了以 ANSYS、NASTRAN等为代表的应用极为广泛的结构有限元分析软件。计算机技术在机构 的静力学分析、运动学分析、动力学分析以及控制系统分析上的应用，则在二十世纪八十 年代形成了计算多体系统动力学，并产生了以ADAMS和DADS为代表的动力学分析软件。两者共同构成计算机辅助工程（CAE）技术的重要内容。 多体系统是指由多个物体通过运动副连接的复杂机械系统。多体系统动力学的根本目 的是应用计算机技术进行复杂机械系统的动力学分析与仿真。它是在经典力学基础上产生 的新学科分支，在经典刚体系统动力学上的基础上，经历了多刚体系统动力学和计算多体 系统动力学两个发展阶段，目前已趋于成熟。 多刚体系统动力学是基于经典力学理论的，多体系统中最简单的情况——自由质点和 一般简单的情况——少数多个刚体，是经典力学的研究内容。多刚体系统动力学就是为多 个刚体组成的复杂系统的运动学和动力学分析建立适宜于计算机程序求解的数学模型，并 寻求高效、稳定的数值求解方法。由经典力学逐步发展形成了多刚体系统动力学，在发展 过程中形成了各具特色的多个流派。 早在1687年，牛顿就建立起牛顿方程解决了质点的运动学和动力学问题；刚体的概念 最早由欧拉于1775年提出，他采用反作用力的概念隔离刚体以描述铰链等约束，并建立了 机械系统动力学分析及ADAMS应用 经典力学中的牛顿-欧拉方程。1743年，达朗贝尔研究了约束刚体系统，区分了作用力和反 作用力，达朗贝尔将约束反力称为“丢失力”，并形成了虚功原理的初步概念。1788年，拉格朗日发表了《分析力学》，系统地研究了约束机械系统，他系统地考虑了约束，并提出 了广义坐标的概念，利用变分原理考虑系统的动能和势能，得出第二类拉格朗日方程—— 最少数量坐标的二阶常微分方程（ODE）；并利用约束方程与牛顿定律得出带拉格朗日乘子 的第一类拉格朗日方程——最大数量坐标的微分代数方程（DAE）。虚功形式的动力学普遍 方程尚不能解决具有非完整约束的机械系统问题，1908年若丹给出了若丹原理——虚功率 形式的动力学普遍方程，利用若丹原理可以方便地讨论碰撞问题和非完整系统的动力学问 题。 对于由多个刚体组成的复杂系统，理论上可以采用经典力学的方法，即以牛顿-欧拉方法为代表的矢量力学方法和以拉格朗日方程为代表的分析力学方法。这种方法对于单刚体 或者少数几个刚体组成的系统是可行的，但随着刚体数目的增加，方程复杂度成倍增长， 寻求其解析解往往是不可能的。后来由于计算机数值计算方法的出现，使得面向具体问题 的程序数值方法成为求解复杂问题的一条可行道路，即针对具体的多刚体问题列出其数学 方程，再编制数值计算程序进行求解。对于每一个具体的问题都要编制相应的程序进行求 解，虽然可以得到合理的结果，但是这个过程长期的重复是让人不可忍受的，于是寻求一 种适合计算机操作的程式化的建模和求解方法变得迫切需要了。在这个时候，也就是20世纪60年代初期，在航天领域和机械领域，分别展开了对于多刚体系统动力学的研究，并且 形成了不同派别的研究方法。 最具代表性的几种方法是罗伯森-维滕堡（Roberson-Wittenburg）方法、凯恩（Kane）方法、旋量方法和变分方法。 罗伯森与维滕堡于1966年提出一种分析多刚体系统的普遍性方法，简称为R/W方法，这种方法的主要特点是利用图论的概念及数学工具描述多刚体系统的结构，以邻接刚体之 间的相对位移作为广义坐标，导出适合于任意多刚体系统的普遍形式动力学方程，并利用 增广体概念对方程的系数矩阵作出物理解释。R/W方法以十分优美的风格处理了树结构多 刚体系统，对于非树系统，通过铰切割或刚体分割方法将非树系统转变成树系统进行处理。 凯恩方法是在1965年左右形成的分析复杂系统的一种方法，其利用广义速率代替广义 坐标描述系统的运动，直接利用达朗伯原理建立动力学方程，并将矢量形式的力与达朗伯 惯性力直接向特定的基矢量方向投影以消除理想约束力，兼有矢量力学和分析力学的特点， 既适用完整系统，也适用于非完整系统。 旋量方法是一种特殊的矢量力学方法（或牛顿-欧拉方法，简称为N/E方法），其特点是将矢量与矢量矩合为一体，采用旋量的概念，利用对偶数作为数学工具，使N/E方程具有极其简明的表达形式，在开链和闭链空间机构的运动学和动力学分析得到广泛运用。 变分方法是不同于矢量力学或分析力学的另一类分析方法，高斯最小拘束原理是变分 方法的基本原理，保保夫和里洛夫从这一原理出发发展了两种不同风格的计算方法。该方 法有利于结合控制系统的优化进行综合分析，而且由于其不受铰的约束数目的影响，适用 于带多个闭环的复杂系统。 这几种方法构成了早期多刚体系统动力学的主要内容，借助计算机数值分析技术，可 以解决由多个物体组成的复杂机械系统动力学分析问题。但是多体系统动力学在建模与求 第2章 多体系统动力学基本理论 解方面的自动化程度，相对于结构有限元分析的成熟来说相差甚远。正是为了解决多体系 统动力学建模与求解的自动化问题，美国Chace和Haug于80年代提出了适宜于计算机自动建模与求解的多刚体系统笛卡尔建模方法，这种方法不同于以罗伯森-维滕堡方法为代表的拉格朗日方法，它是为以系统中每个物体为单元，建立固结在刚体上的坐标系，刚体的 位置相对于一个公共参考基进行定义，其位置坐标统一为刚体坐标系基点的笛卡尔坐标与 坐标系的方位坐标，再根据铰约束和动力学原理建立系统的数学模型进行求解。 20世纪80年代，Haug等人确立了“计算多体系统动力学”这门新的学科，多体系统 动力学的研究重点由多刚体系统走向侧重多柔体系统，柔性多体系统动力学成为计算多体 系统动力学的重要内容。 柔性多体系统动力学在20世纪70年代逐渐引起人们的注意，一些系统如高速车辆、 机器人、航天器、高速机构、精密机械等其中柔性体的变形对系统的动力学行为产生很大 影响。二十多年来柔性多体系统动力学一直是研究热点，这期间产生了许多新的概念和方 法，有浮动标架法、运动-弹性动力学方法、有限段方法以及最新提出的绝对节点坐标法等， 其中浮动标架法最早是在航天领域研究中提出来的。 计算多体系统动力学是指用计算机数值手段来研究复杂机械系统的静力学分析、运动 学分析、动力学分析以及控制系统分析的理论和方法。相比于多刚体系统，对于柔性体和 多体与控制混合问题的考虑是其重要特征。其具体任务为： 1．建立复杂机械系统运动学和动力学程式化的数学模型，开发实现这个数学模型的软 件系统，用户只需输入描述系统的最基本数据，借助计算机就能自动进行程式化处理。 2．开发和实现有效的处理数学模型的计算方法与数值积分方法，自动得到运动学规律 和动力学响应。 3．实现有效的数据后处理，采用动画显示、图表或其他方式提供数据处理结果。 计算多体系统动力学的产生极大地改变了传统机构动力学分析的面貌，使工程师从传 统的手工计算中解放了出来，只需根据实际情况建立合适的模型，就可由计算机自动求解， 并可提供丰富的结果分析和利用手段；对于原来不可能求解或求解极为困难的大型复杂问 题，现可利用计算机的强大计算功能顺利求解；而且现在的动力学分析软件提供了与其它 工程辅助或分析软件的强大接口功能，它与其它工程辅助设计和分析软件一起提供了 完整的计算机辅助工程（CAE）技术。 2.1.2 多体系统动力学研究活动 自20世纪60年代以来，国内外在多体系统动力学方面多次召开了深具意义的会议。 国际范围内的会议主要有： , 1977年国际理论和应用力学学会（International Union of Theoretical and Applied Mechanics - IUTAM）发起在德国慕尼黑由Magnus主持召开第一次多刚体系统动力 学讨论会。 , 1983年北大西洋公约组织与美国国家科学基金委等（NATO-NSF-ARD）联合组织在 美国爱阿华由Haug主持召开“机械系统动力学计算机辅助分析与优化高级研讨会”。 , 1985年第八届国际车辆动力学协会（International Association of Vehicle System 3 机械系统动力学分析及ADAMS应用 Dynamics - IAVSD）会议，Kortum 和Schiehlen发表了用于车辆动力学仿真的多体软 件。 , 1985年IUTAM与国际机器及机构理论联合会（IFTOMM）联合在意大利Udine由Bianchi和Schiehlen主持举行了第二届国际多体系统动力学讨论会，这次会议了 该领域的进展，标志多刚体系统动力学已趋于成熟。 , 1989年由德国斯图加特大学主持对当时比较先进的大型软件进行测试，编辑出版了 “多体系统手册”；以后几乎每年都有国际的多体系统动力学的会议，并出现了多体 系统动力学的专门的刊物。 , 在国内召开的关于多体系统动力学方面的重要会议主要有： , 1986年由中国力学学会一般力学专业委员会在北京主持召开“多刚体系统动力学” 研讨会。 , 1988年在长春召开“柔性多体系统动力学研讨会”。 , 1992年在上海召开“全国多体系统动力学—理论、计算方法与应用学术会议”。 , 1996年由中国力学学会一般力学专业委员会与中国空间学会空间机械委员会联合在 山东长岛召开“全国多体系统动力学与控制学术会议”。 除了国内外重要的会议之外，国内外还出版了多种关于多体系统动力学的优秀教材或 专著。国外的优秀教材和专著主要有： , 1977年，Witenburg出版了“Dynamics of Systems of Rigid Bodies”，这是第一本关于多体系统的教材，其奠定了多刚体系统动力学拉朗日方法的基础，并将图论方法引入 到多体系统动力学中。 , 1986年，Schiehlen出版了“Technische Dynamik”（德语），其将多体系统作为与有限 元系统及连续系统相当的系统来统一考虑。 , 1988年，Nikravesh出版了“Computer-Aided Analysis of Mechanical Systems”。 , 1988年，Roberson 和Schwertassek出版了“Dynamics of Multibody System”。 , 1989年，Haug出版了“Computer-Aided Kinematics and Dynamics of Mechanical Systems”，提出了多刚体系统建模的笛卡尔方法。 , 1989年，Shabana出版了“Dynamics of Multibody Systems”，这是一本关于多柔体动 力学的专著。 , 1990年，Huston出版了“Multibody Dynamics”，其讨论了多体系统动力学、数值方 法以及柔性多体系统动力学。 , 1992年，Bremer和Pfeiffer出版了“Elastische Mehrkorpersystem”（德语），其主要讨论了柔性体，并给出了诸多工程实例。 , 1992年，Amirouche出版了“Computational Methods for Multibody Dynamics”，重点讨论了矩阵方法。 , 1994年，Garcia 和Bayo出版了“Kinematic and Dynamic Simulation of Multibody Systems”，提出了完全笛卡尔坐标方法，并给出了一种求解效率高的计算方法用于实 时仿真的需要。 , 1994年，Shabana出版了“Computational Dynamics”，主要是有关多刚体系统。 , 1996年，Stejskal和Valasek出版了“Kinematics and Dynamics of Machinery”，从空间 第2章 多体系统动力学基本理论 机构的CAD设计入手，讨论了自由刚体问题，给出了高副和低副运动学约束的描述， 并讨论了动力学分析以及数值计算方面的问题。 国内有影响的教材和专著主要有： , 刘延柱，洪嘉振，杨海兴. 多刚体系统动力学. 北京：高等教育出版社，1989 , 黄文虎 邵成勋 等. 多柔体系统动力学. 北京：科学出版社，1996 , 陆佑方. 柔性多体系统动力学. 北京：高等教育出版社，1996 , 洪嘉振. 计算多体系统动力学. 北京：高等教育出版社，1999 2.1.3 多体系统动力学研究现状 计算多体系统动力学中所研究的多体系统，根据系统中物体的力学特性可分为多刚体 系统、柔性多体系统和刚柔混合多体系统。多刚体系统是指可以忽略系统中物体的弹性变 形而将其当作刚体来处理的系统，该类系统常处于低速运动状态；柔性多体系统是指系统 在运动过程中会出现物体的大范围运动与物体的弹性变形的耦合，从而必须把物体当作柔 性体处理的系统，大型、轻质而高速运动的机械系统常属此类；如果柔性多体系统中有部 分物体可以当作刚体来处理，那么该系统就是刚柔混合多体系统，这是多体系统中最一般 的模型。 1．多体系统建模理论 对于多刚体系统，从20世纪60年代到80年代，在航天和机械两个领域形成了两类不 同的数学建模方法，分别称为拉格朗日方法和笛卡尔方法；20世纪90年代，在笛卡尔方法的基础上又形成了完全笛卡尔方法。这几种建模方法的主要区别在于对刚体位形描述的 不同。 航天领域形成的拉格朗日方法，是一种相对坐标方法，以Roberson-Wittenburg方法为 代表，是以系统每个铰的一对邻接刚体为单元，以一个刚体为参考物，另一个刚体相对该 刚体的位置由铰的广义坐标（又称拉格朗日坐标）来描述，广义坐标通常为邻接刚体之间 的相对转角或位移。这样开环系统的位置完全可由所有铰的拉格朗日坐标阵所确定。其q动力学方程的形式为拉格朗日坐标阵的二阶微分方程组，即 ，，，A(q,t)q,B(q,q,t) (2-1) 这种形式首先在解决拓扑为树的航天器问题时推出。其优点是方程个数最少，树系统的坐 标数等于系统自由度，而且动力学方程易转化为常微分方程组（ODEs - Ordinary Differential Equations）。但方程呈严重非线性，为使方程具有程式化与通用性，在矩阵与中常常AB 包含描述系统拓扑的信息，其形式相当复杂，而且在选择广义坐标时需人为干预，不利于 计算机自动建模。不过目前对于多体系统动力学的研究比较深入，现在有几种应用软件采 用拉格朗日的方法也取得了较好的效果。 对于非树系统，拉格朗日方法要采用切割铰的方法以消除闭环，这引入了额外的约束， 使得产生的动力学方程为微分代数方程，不能直接采用常微分方程算法去求解，需要专门 的求解技术。 机械领域形成的笛卡尔方法是一种绝对坐标方法，即Chace和Haug提出的方法，以 系统中每一个物体为单元，建立固结在刚体上的坐标系，刚体的位置相对于一个公共参考 5 机械系统动力学分析及ADAMS应用 基进行定义，其位置坐标（也可称为广义坐标）统一为刚体坐标系基点的笛卡尔坐标与坐 标系的方位坐标，方位坐标可以选用欧拉角或欧拉参数。单个物体位置坐标在二维系统中 为3个，三维系统中为6个（如果采用欧拉参数为7个）。对于由N个刚体组成的系统，位置坐标阵中的坐标个数为3N（二维）或6N（或7N）（三维），由于铰约束的存在，这q 些位置坐标不独立。系统动力学模型的一般形式可表示为 T，，,,Aq，,,Bq (2-2) ,,(q,t),0, 式中,为位置坐标阵的约束方程，为约束方程的雅可比矩阵，为拉格朗日乘子。,q,q 这类数学模型就是微分-代数方程组（DAEs - Differential Algebraic Equations），也称为欧拉- 拉格朗日方程组（Euler-Lagrange Equations），其方程个数较多，但系数矩阵呈稀疏状，适 宜于计算机自动建立统一的模型进行处理。笛卡尔方法对于多刚体系统的处理不区分开环 与闭环（即树系统与非树系统），统一处理。目前国际上最著名的两个动力学分析商业软件 ADAMS和DADS都是采用这种建模方法。 完全笛卡尔坐标方法，由Garcia和Bayo于1994年提出，是另一种形式的绝对坐标方 法。这种方法的特点是避免使用一般笛卡尔方法中的欧拉角或欧拉参数，而是利用与刚体 固结的若干参考点和参考矢量的笛卡尔坐标描述刚体的空间位置与姿态。参考点选择在铰 的中心，参考矢量沿铰的转轴或滑移轴，通常可由多个刚体共享而使未知变量减少。完全 笛卡尔坐标所形成的动力学方程与一般笛卡尔方法本质相同，只是其雅可比矩阵为坐标线 性函数，便于计算。 至于柔性多体系统，从计算多体系统动力学角度看，柔性多体系统动力学的数学模型 首先应该和多刚体系统与结构动力学有—定的兼容性。当系统中的柔性体变形可以不计时， 即退化为多刚体系统。当部件间的大范围运动不存在时，即退化为结构动力学问题。 柔性多体系统不存在连体基，通常选定一浮动坐标系描述物体的大范围运动，物体的 弹性变形将相对该坐标系定义。弹性体相对于浮动坐标系的离散将采用有限单元法与现代 模态综合分析方法。在用集中质量有限单元法或一致质量有限单元法处理弹性体时，用结 点坐标来描述弹性变形。在用正则模态或动态子结构等模态分析方法处理弹性体时用模态 坐标描述弹性变形。这就是莱肯斯首先提出的描述柔性多体系统的混合坐标方法。即用坐 标阵TTTp,(qa)描述系统的位形，其中q为浮动坐标系的位形坐标，a为变形坐标。考虑到多刚体系统的两种流派，在柔性多体系统动力学中也相应提出两种混合坐标，即浮动 坐标系的拉格朗日坐标加弹性坐标与浮动坐标系的笛卡尔坐标加弹性坐标。 根据动力学基本原理推导的柔性多体系统动力学方程，形式同式（2-1）和（2-2），只是将q用p代替。即，柔性多体系统具有与多刚体系统类同的动力学数学模型。 2．多体系统动力学数值求解 多刚体系统拉格朗日方法产生的形如式（2-1）的动力学数学模型，是形式复杂的二阶 常微分方程组（ODEs），系数矩阵包含描述系统拓扑的信息。对于该类问题的求解，通常 采用符号-数值相结合的方法或者全数值的方法。符号-数值方法是先采用基于计算代数的 符号计算方法，进行符号推导，得到多刚体系统拉格朗日模型系数矩阵简化的数学模型， 第2章 多体系统动力学基本理论 再用数值方法求解ODE问题。鉴于计算机技术的发展，目前全数值方法也较为流行，就是 将多刚体系统拉格朗日数学模型当作一般ODE问题进行求解，这方面的技术已经较为成 熟。 多刚体系统笛卡尔方法产生的形如式（2-2）的动力学数学模型，是著名的微分-代数 方程组（DAEs）。DAE问题是计算多体系统动力学领域的热点问题。 柔性多体系统的动力学数学模型，其形式与多刚体系统相同，可以借鉴多刚体系统数 学模型的求解方法。只是混合坐标中描述浮动坐标系运动的刚体坐标q通常是慢变大幅值的变量，而描述相对于浮动坐标系弹性变形的坐标a却为快变微幅的变量，两类变量出现 在严重非线性与时变的耦合动力学方程中，其数值计算呈病态，将出现多刚体系统中见不 到的数值计算困难。 综上所述，多体系统动力学问题的求解集中于微分-代数方程组的求解，下面将简要地介绍一下DAE问题的求解方法，更具体的介绍，将在（2.5）节进行。 （1）微分-代数方程组的特性 多刚体系统采用笛卡尔方法建模生成的微分-代数方程组为： T，，， (2-3) M(q,t)q，,(q,t,),Q(q,q,t),0q (2-4) ,(q,t),0mn,,R其中，，，，q,R、、分别是系统位置、速度、加速度向量，是拉格朗日乘qq m，nnn，n子，Q,Rt,R是时间，为机械系统惯性矩阵，,,R为约束雅可比矩阵，M,Rqm为外力向量，为位置约束方程。 ,,R 将式(1.4)对时间求一阶和二阶导数，得到速度和加速度约束方程： ，，， (2-5) ,(q,q,t),,(q,t)q,,(q,t),0q ，，，，，，，， (2-6) ,(q,q,q,t),,(q,t)q,,(q,q,t),0q ，，，其中，,,,(,q)q,2,q,,,,,,(q,t)称为速度右项，称为加速度右项。 qqqtttt 给定方程组初始条件： q(0),q,0 (2-7) ,，，q(0)q,0, 微分-代数方程组的特性和需要注意的问题有： , 微分-代数方程问题不是常微分方程（ODE）问题； , 由式(2-3)和(2-4)组成的微分-代数方程组是指标3问题，通过对约束方程求导化为由 式(2-3)-( 2-6)组成的微分-代数方程组后，其指标降为1； , 微分-代数方程数值求解的关键在于避免积分过程中代数方程的违约现象； , 初值式(2-7)与位置约束式(2-4)及速度约束式(25)的相容性； , 微分-代数方程组的刚性问题。 2- 自二十世纪七十年代以来，国际上对微分-代数方程问题作了大量的研究，时至如今， 新的算法仍不断涌现。根据对位置坐标阵和拉格朗日乘子处理技术的不同，可以将微分- 代数方程组问题的处理方法分为增广法和缩并法。 传统的增广法是把广义坐标加速度，，,q和拉格朗日乘子作为未知量同时求解，再对加 7 机械系统动力学分析及ADAMS应用 ，，，速度进行积分求出广义坐标速度及广义坐标位置，包括直接积分法和约束稳定法。qqq 近十年来，在传统增广法的基础上又发展形成了超定微分-代数方程组（ODAEs）方法等新的一类算法。 直接积分法：将式(2-3)和(2-6)联立在一起，同时求出，，，，，与，然后对积分得和。,qqqq该方法未考虑式(2-4)和(2-5)的坐标和速度违约问题，积分过程中误差积累严重，很易发散。 在实际的数值计算过程中，并不直接采用直接积分法，但在直接积分法的基础上发展了一 系列控制违约现象的数值方法。 约束稳定法：将控制反馈理论引入微分-代数方程组的数值积分过程以控制违约现象。 通过把式(2-6)右边量替换为含位置约束和速度约束的参数式，保证位置约束和速度约束在 式(2-3)和(2-6)联立求解时恒满足。该方法稳定性好，响应快，但如何选择参数式中速度项 和位置项适当的系数是一个问题。 超定微分-代数方程组（ODAEs）法：将系统速度作为变量引入微分-代数方程组，从而将原来的二阶DAE化为超定的一阶DAE，再为所得方程组引入未知参数，根据模型的 相容性消除系统的超定性，如此可使数值计算的稳定性明显改变。或者将系统位置、速度、 加速度向量和拉格朗日乘子向量联立作为系统广义坐标，再将由式(2-3)、(2-4)、(2-5)和(2-6)组成的微分-代数方程组及速度与位置、加速度与速度的微分关系式作为约束，化二阶DAE为超定的一阶DAE，再根据系统相容性引入二个未知参数，消除超定性，这样所得的最终 约化模型更为简单，但方程组要多n个。在ODAE方法的基础上产生了一系列新的更为有 效的算法。 解耦ODAE法：在ODAE方法的基础上，发展形成了一类解耦思想，就是在ODAEs基础上，对常用的隐式ODE方法采用预估式，再按加速度、速度和位置的顺序进行求解。 后来进一步发展形成了无需对隐式ODE方法利用预估式的解耦思想，更一步地提高了效 率。 缩并法就是通过各种矩阵分解方法将描述系统的n个广义坐标用p个独立坐标表达，从而将微分-代数方程组从数值上化为与式(2-1)类似的数学模型，如此易于用ODE方法进行求解。传统的缩并法包括LU分解法、QR分解法、SVD分解法以及零空间方法等，后 来在传统缩并法的基础上产生了局部参数化缩并方法等新的算法。缩并法中的这些具体方 法，分别对应着约束雅可比矩阵的不同分解。 LU分解法：又称为广义坐标分块法。把广义位置坐标 uv用相关坐标和独立坐标分q 块表示，再将约束雅可比矩阵,，，，用LU分解法分块，得到广义坐标速度q、加速度q用独q 立坐标速度，，，vv、加速度表达的式子。将这两个表达式代入式(2-3)，就可得到形如式(2-1)的关于独立坐标加速度，，v的二阶微分方程。该算法可靠、精确，并可控制误差，但效率稍 低。 QR分解法：通过对约束雅可比矩阵,正交分解的结果作微分流型分析，得到可选作q 受约束系统独立速度的，，并将微分-代数方程组化作形如式(2.1-1)的关于，，的二阶微分方zz程，如此可保证在小时间间隔内由，积分引起的广义坐标的变化不会导致大的约束违约。 z SVD分解法：把约束雅可比矩阵,作奇异值分解所得结果分别用于式(2-3)和(2-6)，q 得到缩并后的系统动力学方程。在该方法推导过程中没有用到式(2-4)和(2-5)，所以也存在 第2章 多体系统动力学基本理论 位置和速度违约问题，可用约束稳定法改善其数值性态。 可微零空间法：通过Gram-Schmidt正交化过程自动产生约束雅可比矩阵,的可微、qm，n(n,m)，n唯一的零空间基，来对系统方程降阶。具体做法是对由和任意矩阵,,RB,Rqn，n构造的矩阵采用Gram-Schmidt正交化过程，将化为正交非奇异矩阵。再引VP,RP Tn入新的速度矢量，，z,Vq，，使满足，将新速度矢量和加速度矢量按正交化结果，，，z,Rzz分块，得到新的独立速度矢量，，，和加速度矢量。如此可将微分-代数方程组化为关于新zzII 的独立加速度矢量，，的动力学方程。 zI 局部参数化缩并方法：先将式(2-3)～( 2-6)改写为等价的一阶形式，再用微分流形理论 的切空间局部参数化方法将等价的欧拉-拉格朗日方程降为参数空间上的常微分方程。 总的说来，微分-代数方程组数值求解的方法都可归为增广法或缩并法，除了上面所介 绍的这些增广法和缩并法所运用的增广和缩并技术外，近几年来还出现了不少独具特色的 处理算法，它们或者是在数值求解算法中独具匠心，或者针对某些具体情况作了专门研究。 （3）相容性问题和刚性问题 初值相容性问题：在微分-代数方程组的数值求解过程中，给定的位置和速度初始条件 与微分-代数方程组中的位置和速度约束的相容性是值得注意的一个问题。相容性是微分-代数方程组有解的必要条件。 刚性问题：由于现代机械系统的复杂性，会由于系统的耦合而使所得到的微分-代数方程组呈现刚性特性。对于刚性问题的求解，目前最常用的方法是隐式方法，隐式方法不仅 用于求解刚性问题，而且相比于显式方法具有更好的稳定性和计算精度。近几年来，无论 是在LU分解法基础上发展起来的新缩并法，还是基于ODAE方法的增广法，或是基于多体系统正则方程的解法，应用的无不是隐式方法。关于刚性问题的进一步介绍见2.5节。 2.2 多刚体系统动力学建模 计算多体系统动力学分析，首先在于提供一个友好方便的界面以利于建立多体系统的 力学模型，并在系统内部由多体系统力学模型得到动力学数学模型；再者需要有一个优良 的求解器对数学模型进行求解，求解器要求效率高、稳定性好，并具有广泛的适应性；最 后还需要对求解结果提供丰富的显示查询手段。这其中的关键技术就是自动建模技术和求 解器设计，所谓自动建模就是由多体系统力学模型自动生成其动力学数学模型，求解器的 设计则必须结合系统的建模，以特定的动力学算法对模型进行求解。 本节首先给出多体系统动力学的一些基本概念，再介绍计算多体系统动力学建模与求 解的一般过程，然后作为重点，按豪格（Haug）的笛卡尔方法给出多体系统的运动学和动 力学数学模型，这是考虑到目前最负盛名的两个多体系统软件都是采用笛卡尔方法，最后 简要叙述自动建模技术。 9 机械系统动力学分析及ADAMS应用 2.2.1 多体系统动力学基本概念 ：这里也称力学模型，由物体、铰、力元和外力等要素组成并具有一定拓扑 构型的系统。 ：多体系统中各物体的联系方式称为系统的拓扑构型，简称拓扑。根据系统 拓扑中是否存在回路，可将多体系统分为树系统与非树系统。系统中任意两个物体之间的 通路唯一，不存在回路的，称为树系统；系统中存在回路的称为非树系统。 ：多体系统中的构件定义为物体。在计算多体系统动力学中，物体区分为刚性体 （刚体）和柔性体（柔体）。刚体和柔体是对机构零件的模型化，定义为质点间距离保持不变的质点系，定义为考虑质点间距离变化的质点系。 ：对系统中某构件的运动或构件之间的相对运动所施加的限制称为约束。约束分 为运动学约束和驱动约束，运动学约束一般是系统中运动副约束的代数形式，而驱动约束 则是施加于构件上或构件之间的附加驱动运动条件。 ：也称为运动副，在多体系统中将物体间的运动学约束定义为铰。铰约束是运动学 约束的一种物理形式。 ：在多体系统中物体间的相互作用定义为力元，也称为内力。力元是对系统中弹 簧、阻尼器、致动器的抽象，理想的力元可抽象为统一形式的移动弹簧-阻尼器-致动器 （TSDA），或扭转弹簧-阻尼器-致动器（RSDA）。 ：多体系统外的物体对系统中物体的作用定义为外力（偶）。 ：分为静力学数学模型、运动学数学模型和动力学数学模型，是指在相应条 件下对系统物理模型（力学模型）的数学描述。 ：装配在一起并允许作相对运动的若干个刚体的组合。 ：研究组成机构的相互联接的构件系统的位置、速度和加速度，其与产生运动 的力无关。运动学数学模型是非线性和线性的代数方程。 ：研究外力（偶）作用下机构的动力学响应，包括构件系统的加速度、速度和 位置，以及运动过程中的约束反力。动力学问题是已知系统构型、外力和初始条件求运动， 也称为动力学正问题。动力学数学模型是微分方程或者微分方程和代数方程的混合。 ：在与时间无关的力作用下系统的平衡，称为静平衡。静平衡分析一种特殊的 动力学分析，在于确定系统的静平衡位置。 ：逆向动力学分析是运动学分析与动力学分析的混合，是寻求运动学上确 定系统的反力问题，与动力学正问题相对应，逆向动力学问题是已知系统构型和运动求反 力，也称为动力学逆问题。 ：固定在刚体上并随其运动的坐标系，用以确定刚体的运动。刚体上每一 个质点的位置都可由其在连体坐标系中的不变矢量来确定。 ：唯一地确定机构所有构件位置和方位即机构构形的任意一组变量。广义坐 标可以是独立的（即自由任意地变化）或不独立的（即需要满足约束方程）。对于运动系统 来说，广义坐标是时变量。 ：确定一个物体或系统的位置所需要的最少的广义坐标数，称为该物体或系统 的自由度。 第2章 多体系统动力学基本理论 ：对系统中某构件的运动或构件之间的相对运动所施加的约束用广义坐标表 示的代数方程形式，称为约束方程。约束方程是约束的代数等价形式，是约束的数学模型。 2.2.2 计算多体系统动力学建模与求解一般过程 一个机械系统，从初始的几何模型，到动力学模型的建立，经过对模型的数值求解， 最后得到分析结果，其如图2.1所示。 计算多体系统动力学分析的整个流程，主要包括建模和求解两个阶段。建模分为物理 建模和数学建模，物理建模是指由几何模型建立物理模型，数学建模是指从物理模型生成 数学模型。几何模型可以由动力学分析系统几何造型模块所建造，或者从通用几何造型软 件导入。对几何模型施加运动学约束、驱动约束、力元和外力或外力矩等物理模型要素， 形成表达系统力学特性的物理模型。物理建模过程中，有时候需要根据运动学约束和初始 位置条件对几何模型进行装配。由物理模型，采用笛卡尔坐标或拉格朗日坐标建模方法， 应用自动建模技术，组装系统运动方程中的各系数矩阵，得到系统数学模型。对系统数学 模型，根据情况应用求解器中的运动学、动力学、静平衡或逆向动力学分析算法，迭代求 解，得到所需的分析结果。联系设计目标，对求解结果再进行分析，从而反馈到物理建模 过程，或者几何模型的选择，如此反复，直到得到最优的设计结果。 在建模和求解过程中，涉及到几种类型的运算和求解。首先是物理建模过程中的几何 模型装配，图2.1中称为“初始条件计算”，这是根据运动学约束和初始位置条件进行的， 是非线性方程的求解问题；再就是数学建模，是系统运动方程中的各系数矩阵自动组装过 程，涉及大型矩阵的填充和组装问题；最后是数值求解，包括多种类型的分析计算，如运 动学分析、动力学分析、静平衡分析、逆向动力学分析等。运动学分析是非线性的位置方 程和线性的速度、加速度方程的求解，动力学分析是二阶微分方程或二阶微分方程和代数 方程混合问题的求解，静平衡分析从理论上讲是一个线性方程组的求解问题，但实际上往 往采用能量的方法，逆向动力学分析是一个线性代数方程组求解问题，这里面，最复杂的 是动力学微分代数方程的求解问题，它是多体系统动力学的核心问题。 在多体系统动力学建模与求解过程中，还有一个问题是值得注意的——初值相容性问 题，这是在任何正式求解之前必须首先解决的问题，直接影响到问题的可解性。初值相容 性是要求系统中所有的位置、速度初始条件必须与系统运动学约束方程相容。对于简单问 题，初值相容性是易于保证的，但对于大型复杂系统，必须有专门的初值相容性处理算法 以判断系统的相容性或由一部分初值计算相容的其它初值。 在多体系统建模与求解过程，求解器是核心，这其中涉及的所有运算和求解，如初始 条件计算、方程自动组装、各种类型的数值求解等等都由求解器所支持，它提供了所需的 全部算法。 实际上，结果分析是需要有专门的数值后处理器来支持的，以提供曲线和动画显示以 及其它各种辅助分析手段。但相比于多体系统建模与求解，数值后处理器相对简单，不存 在什么理论上的重要问题。 11 机械系统动力学分析及ADAMS应用 建 模 求 解 物理建模 数学建模 数值求解 几何模型 物理模型 数学模型 分析结果 初始条件计算 方程自动组集 分析迭代 求 解 器 图2.1 计算多体系统动力学建模与求解一般过程 2.2.3 多刚体系统运动学 对于多体系统的运动学分析，传统的理论力学是以刚体位置、速度和加速度的微分关 系以及矢量合成原理为基础进行分析的，而计算多体系统动力学中的运动学分析则是以系 统中连接物体与物体的运动副为出发点，所进行的位置、速度和加速度分析都是基于与运 动副对应的约束方程来进行的。 基于约束的多体系统运动学，首先寻求与系统中运动副等价的位置约束代数方程，再 由位置约束方程的导数得到速度、加速度的约束代数方程，对这些约束方程进行数值求解， 可得到广义位置坐标及相应的速度和加速度坐标，最后根据坐标变换就可以由系统广义坐 标及相应导数得到系统中任何一点的位置、速度和加速度。 由于机械系统在二维空间运动时，广义坐标、约束方程、问题规模以及问题求解都相 对简单，故本节先讨论二维多体系统运动学以解释多体系统运动学基本理论，在此基础上 再给出三维多体系统的运动学方程。 1 设一个平面机构由nb个刚性构件组成。在机构所在平面上建立一个全局坐标系，xoy ,,,机构在该坐标系中运动；再为机构上每个构件xoyi建立各自的连体坐标系，可由连体坐iiiT,标系的运动确定构件的运动。选定构件r,[x,y]oi连体坐标系原点的全局坐标和连体iiiiT坐标系相对于全局坐标系的转角q,[x,y,,],组成构件i的笛卡尔广义坐标矢量，如iiiii 图2.2所示。由nbnc,3，nb个刚性构件组成的系统的广义坐标数，则系统广义坐标矢量 TTTT可表示为q,[q,q,...,q]。 nb12 第2章 多体系统动力学基本理论 图2.2 平面笛卡尔广义坐标 一个实际的机械系统，系统中构件与支架或构件与构件之间存在运动副的联接，这些 运动副可以用系统广义坐标表示为代数方程。设表示运动副的约束方程数为，则用系统nh 广义坐标矢量表示的运动学约束方程组为： KKKKT,(q),[,(q),,(q),...,,(q)],0 (2-8) 12nh 这里给出的是定常完整约束情况。如果约束方程与时间相关，则自变量中显含时间项， 这种约束被称为非定常约束；更一般的约束方程含有不可积速度项的不等式或关系式，这 种约束称为非完整约束。一般的运动学约束是定常完整约束。 对于一个有nc个广义坐标和nh个约束方程的机械系统，若nc,nh，且这nh个约束 方程是独立、相容的，则系统自由度DOF,nc,nh。为使系统具有确定运动，可以有二种方法： （1）为系统添加与系统自由度DOF相等的附加驱动约束； （2）对系统施加力的作用。 在（1）情况下，系统实际自由度为零，被称为是在运动学上确定的，在此情况下求解 系统运动过程中的位置、速度和加速度的分析是运动学分析，运动学分析本身不涉及作用 力或反作用力问题。但是对于运动学上确定的系统，可以求解系统中约束反力，即已知运 动求作用力，这是动力学逆问题。 在（2）情况下，系统有着大于零的自由度，但是在外力作用力，对于具有确定构型和 特定初始条件的系统，其动力学响应是确定的，这种情况下求解系统运动过程中的位置、 速度和加速度的分析，称为动力学分析。在这种情况下，特殊地，如果外力与时间无关， 可以求解系统的静平衡位置，这就是静平衡分析问题。 考虑运动学分析，为使系统具有确定运动，也就是要使系统实际自由度为零，为系统 施加等于自由度（nc,nh）的驱动约束： D,(q,t),0 (2-9) 在一般情况下，驱动约束是系统广义坐标和时间的函数。驱动约束在其集合内部及其 与运动学约束合集中必须是独立和相容的，在这种条件下，驱动系统运动学上是确定的， 将作确定运动。 由式（2-8）表示的系统运动学约束和式（2-9）表示的驱动约束组合成系统所受的全 13 机械系统动力学分析及ADAMS应用 部约束： K，，,(q,t),(q,t),,0 (2-10) ,,D,(q,t)，， 式（2-10）为nc个广义坐标的nc个非线性方程组，其构成了系统位置方程。求解式 （2-10），就可得到系统在任意时刻的广义坐标位置。 q(t) 2 对式(2-10)运用链式微分法则求导，得到速度方程： ，，， (2-11) ,(q,q,t),,(q,t)q，,(q,t),0qt 若令,,,,(q,t)，则速度方程为 t ，，， (2-12) ,(q,q,t),,(q,t)q,,,0q 如果,，是非奇异的，可以求解式(2-12)得到各离散时刻的广义坐标速度。 qq 对式(2-11)运用链式微分法则求导，可得加速度方程 ，，，，，，，，，， (2-13) ,(q,q,q,t),,(q,t)q，(,(q,t)q)q，2,(q,t)q，,(q,t),0qqqqttt ，，，若令,,,(,q)q,2,q,,，则加速度方程为 qqqttt ，，，，，，，， (2-14) ,(q,q,q,t),,(q,t)q,,(q,q,t),0q 如果,，，是非奇异的，可以求解式(2-14)得到各离散时刻的广义坐标加速度。 qq 在速度方程(2-12)和加速度方程(2-14)中出现的矩阵,，称为雅可比矩阵，雅可比矩阵q是约束多体系统运动学和动力学分析中最重要的矩阵。如果的维数为m，q维数为n，,那么,(,),,,,qm，n,nc，nc维数为矩阵，其定义为。这里为的方阵。 qq(,)ijijq 对式(2-12)中的,和式(2-14)中的进行计算时，会涉及到二阶导数，在实际的数值求, 解中，并不是实时地调用求导算法来进行计算，而是先根据具体的约束类型，导出二阶导 数以及雅可比矩阵的表示式，在计算中只需代入基本的数据即可。 3 在确定系统中构件上任意点的运动时，常要求将构件上点从连体坐标系变换到全局坐 标系中，现讨论连体坐标系与全局坐标系的坐标变换及构件上任意点运动。 ,设矢量,,,xoys在全局坐标系和某连体坐标系中分别表示为： xoyiii T,s[s,s],xy, (2-15) ,T,s[s,s],,,,xy, 若任意点,,,xoy在全局坐标系和连体坐标系中坐标如图2-3所示，则存在如下xoyPiii坐标变换关系： PPP,r,r，s,r，As (2-16) P,其中，rO为点在全局坐标系中的坐标，为连体坐标系原点在全局坐标系中的rP,,PP,坐标，ssss为矢量在全局坐标系中坐标，为矢量在连体坐标系中的坐标，为旋转A 变换矩阵，其形式为： ,cos,sin,，， (2-17) ,,AA(,),,sin,cos,，， 第2章 多体系统动力学基本理论 对时间的导数为： A ,sin,,cos,，，d，，，， (2-18) ,,,A,A,,B,,cos,sin,,d,，， 根据式(2-16)，我们可以得到以连体坐标系表示的构件上的任一点的全局坐标。 图2.3 二维空间坐标变换 图2.4 三维空间坐标变换 式(2.2-9)对时间求导数，可得任意点的速度变换公式： PPP，，,,，，，r,r，As,r，,Bs (2-19) 式(2.2-12)对时间求导数，可得任意点的加速度变换公式： PPPP2P，，，，，，，,,,,，，，，，，r,r，,Bs，,Bs,r，,Bs,,As (2-20) 对于一个平面机构来说，进行运动学分析时，先是选定最大集的广义坐标，再分别根 据式(2-10)、(2-12)和(2-14）求解机构在各离散时刻的广义坐标位置、广义坐标速度和广义 坐标加速度。对于任意一个由连体坐标系确定的构件上的点，可以根据式(2-16)、(2-19)和 (2-20)求解其位置、速度和加速度。 4 三维多体系统的运动分析与二维多体系统较为相似，只是广义坐标选取复杂一些，约 束方程形式复杂一些，问题规模要大一些。三维多体系统广义坐标与二维相似，也是由位 置坐标和方位（或称为姿态）坐标组成，位置坐标表示较为固定，都是由连体坐标系基点 坐标确定，方位坐标则具有多种形式，如方向余弦矩阵、欧拉角、卡尔丹角、有限转动四 元数、欧拉参数等等，最常用的是欧拉角和欧拉参数。这里先给出三维机械系统广义坐标 的方向余弦与欧拉参数和欧拉角几种形式及其之间的变换，再据此给出系统的约束方程、 速度方程和加速度方程的形式。 （1）坐标变换、欧拉参数与欧拉角 对于三维空间机构，采用固联在构件上的连体坐标系来确定系统运动。构件的广义坐 标，由两个部分组成，一是连体坐标系的原点坐标，二是确定连体坐标系相对于全局坐标 系的方位参数。如图2. 4所示，连体坐标系T,,,,,,,,,r,[x,y,z]ooxyzoxyz原点坐标为， 相对于全局坐标系的方位可用方向余弦矩阵表示，也可用欧拉参数或者欧拉角，这几oxyz 种具有相同几何意义，但数值特性不同。 15 机械系统动力学分析及ADAMS应用 方向余弦矩阵定义为： aaa，，111213,,A,[f,g,h],aaa (2-21) 212223,, ,,aaa313233，， 其中，,,,,,,,,,,、和分别为连体坐标系坐标轴、和的单位矢量。方向oxozfhoxyzoyg 余弦矩阵为正交矩阵，因此，中9个变量受6个独立方程的约束，方向余弦矩阵中只AA 存在说明3个转动自由度的独立变量。 ,如果连体坐标系,,,,和全局坐标系的原点重合，即，则矢量在连体坐oxyzr,0soxyz标系中的表示形式,s和在全局坐标系中的表示形式存在如下变换关系： s , (2-22) s,As 更一般的坐标变换式为： PP,r,r，As (2-23) P其中，,,,,,为点在坐标系中的坐标，r为坐标系原点在坐标系中ooxyzroxyzoxyzP P,的坐标，,,,,,,,,s为点在坐标系中的坐标，为相对于的方向余弦矩阵。 oxyzoxyzoxyzPA 对式(2-23)求时间导数，得速度变换式： PPP~，,，，，rrAsr,s,，,， (2-24) ~其中,,,,,,是的斜对称矩阵（斜对称矩阵定义见式2-38），称为连体坐标系相,oxyz 对于全局坐标系的角速度矢量，表示为： oxyz T~，,,AA (2-25) ~~若将角速度矢量,,,,,,运用式(2-22)的相关导出式变换到坐标系并表示sA,Asoxyz ,为,，则存在： T~，,,AA, (2-26) ~~，,,,A,A,A (2-27) 对式(2-24)求时间导数，得加速度变换式： PP，，,，，，，r,r，As (2-28) 其中： ~~~，，，A,A,,A,， (2-29) 如果定义与位移,r,r,,和角速度对应的虚位移和虚转动，则式(2-24)、(2-25)、(2-26) 和(2-27)存在相应的变分形式： PPP~,,r,r,As,r,,s,，,， (2-30) T~,,,,AA (2-31) T~,,,,A,A (2-32) ~~,A,A,A,,,,, (2-33) ,角速度,,,和虚转动是不可积的，因此角速度也被称为拟坐标。 根据刚体转动的欧拉定理，确定刚体的方位还可以采用欧拉定理中的转动轴和转动角。,如果坐标系,,,,oxyzu与坐标系原点重合，由欧拉定理知可设绕单位轴矢量转动oxyzoxyz,,,,,oxyzu角与重合，现可由和定义一个4欧拉参数组： ,, 第2章 多体系统动力学基本理论 ,,e,cos0,2,e，，1 (2-34) ,,,,,eeu,,sin2,,2,,,e3，，, 欧拉参数用列向量表示为： 4，1 TTTp,[e,e],[e,e,e,e] (2-35) 00123 欧拉参数要满足欧拉参数归一化约束： T2222pp,e，e，e，e,1 (2-36) 0123 故欧拉参数4个分量中存在3个独立分量，描述物体转动。 欧拉参数和方向余弦矩阵都是描述物体方位的参数，它们是等价的，其间存在着pA 变换关系，从欧拉参数到方向余弦矩阵的变换为： 2T~A,(2e,1)I，2(ee，ee) (2-37) 00~其中e为单位矩阵，为的斜对称矩阵，其表示为： 3，3eI 0,ee，，32,,~ (2-38) ,0,eee31,, ,,,0ee，，21 从方向余弦矩阵到欧拉参数的变换为： 2,e,(trA，1)/40,,(,)/(4)eaae,132230 (2-39) ,,(,)/(4)eaae213310, ,,(,)/(4)eaae321120, 上式中，e,0e,0，若由式(2-39)中第一式计算得，则由下列式子确定欧拉参数。 002111,e,(1，2a,trA)/42,222,，,e(12atrA)/4,2,333,，,e(12atrA)/4, (2-40) ,a21，a12,4e1e2, ,a，a,4ee311313,,a，a,4ee,322323 式中trA为矩阵的迹。 A 为研究欧拉参数与角速度之间的关系，定义两个辅助矩阵： 17 机械系统动力学分析及ADAMS应用 ,,eeee，，1032,,~,,E,,e,e，eI,,eee,e (2-41) 02301,, ,,,e,eee3210，， ,eee,e，，1032,,~,,Ge,eeIeeee,,,，,,, (2-42) 02301,, ,,,ee,ee3210，，可得到如下关系式： TA,EG (2-43) T，，A,2EG (2-44) ， (2-45) ,,2Ep ,， (2-46) ,,2Gp 11TT,，p,E,,G, (2-47) 22 上述式(2-44)～(2-47)的变分形式为： TT,A,,EG，E,G (2-48) (2-49 ,,,2E,p , (2-50) ,,,2G,p 11TT,,p,E,,,G,, (2-51) 22 对于上述各式中所涉及的角速度和虚转动，并不象欧拉参数的导数或变分一样是可积 的，所以在积分用角速度或虚转动表示的运动方程时，不可直接积分，需要利用上述公式 将角速度或虚转动变换为欧拉参数的导数或变分，再作积分运算。 根据欧拉定理，可以将刚体的方位分解为连体坐标系从与全局坐标系重合的始点起， 依次绕连体坐标系自身的,,,,,,ozoxoz、和转过有限角度、,和,来确定，即相对体313,转动序列，三个角度坐标、,和,即为欧拉角，其中为进动角，,为章动角，,为自,,转角。 用欧拉角表示的方向余弦矩阵为： ,c,,c,,s,,c,,s,,c,,s,,s,,c,,c,s,,s，， ,,A,,s,,c，,c,,c,,s,,s,,s，,c,,c,,c,,c,,s (2-52) ,, ,,s,,s,s,,c,c,，， 从方向余弦矩阵到欧拉角的变换为： 第2章 多体系统动力学基本理论 2,,,arccos(a),arcsin(1,a)3333,, (2-53) ,arccos(,/),arcsin(/),a,sa,s,2313 ,,as,as,,arccos(/),arcsin(/)3231,, 当章动角时，式（2-53）失效，进动角和自转角不能确定，称为欧,,n,(n,0,1,...) 拉角奇异点。 欧拉角表示的欧拉参数为： ,,,，,e,cos(),cos(),022,,,,,e,,cos(),sin(),1 (2-54) 22,,,,,,e,sin(),sin()222,,,,，,e,sin(),cos()4,22, 从欧拉参数到欧拉角的变换为： ,,arctg(e/e)，arctg(e/e),3021,22,,arccos[2(e，e),1] (2-55) ,03 ,,,arctg(e/e),arctg(e/e)3021, （2）位置、速度和加速度分析 在一个三维多体系统中，构件的广义坐标矢量由其连体坐标系原点坐标和欧位参数i 组成，表示为： r，，i,q (2-56) i,,pi，， 对于由nb个构件组成的系统，其广义坐标矢量组为： TTTTq,[q,q,...,q] (2-57) nb12 系统广义坐标维数为7nb。 采用欧拉参数广义坐标，每个构件的欧拉参数广义坐标必须满足式(2-36)的归一化约 束，即： pT,,pp,1,0i,1,2,...,nb， (2-58) iii 系统的欧拉参数归一化约束方程的矢量形式为： ppppt ,,[,,,,...,,],0 (2-59) nb12 方程组维数为nb。 与二维系统类似，设与运动副等价的约束方程数为nh，则系统运动学约束方程的矢量 19 机械系统动力学分析及ADAMS应用 形式为： KKKKT ,(q),[,(q),,(q),...,,(q)],0 (2-60) 12nh 为使系统具有确定运动，对系统施加个独立的驱动约束，系统驱动约束方(6nb,nh)程的矢量形式为： D,(q,t),0 (2-61) 据此，由系统的欧拉参数归一化约束方程、运动学约束方程及驱动约束方程组成的系 统约束方程，或称位置方程为： p，，,(q),,K,(q,t),,(q),0 (2-62) ,,D,,,(q,t)，， 式(2.2-55)包含个广义坐标的个方程。 7nb7nb 为进行速度和加速度分析，对式(2.2--58)求微分，并运用式(2-51)，得到欧拉参数归一化约束的变分为： pTTT,,,,2p,p,pG,,,0 (2-63) iiiiii 从而得到： pTT (2-64) ,,pG,0iii,,i式（2-63）和式（2-64）利用了Gp=0这一事实。式(2-64)表明，与欧拉参数归一化约束有关并以，,,为变量的速度方程完全得到满足，且以为变量的加速度方程也将完全得,,到满足，故对于速度和加速度分析，当以角速度和角加速度为变量时，不需要考虑欧拉参 数归一化约束，只考虑运动学约束及驱动约束即可。 运动学约束方程和驱动约束方程对时间求导即得系统速度方程： KKKK,,nb，，，，,,,，，，，,,,,,,rtii，,r,，,, (2-65) ,,,,,,,,,,,iiDDDD,,,,,,,,,,,i,1,r,，，，，t，，，，ii,, 由于运动学约束方程不涉及时间，故： KK,,,0,, (2-66) t 对式(2.2-58)微分可得到系统加速度方程： KKKKKKnbnb，，,,,,，，，，，，，，,,,,，，，，,,,,,,,,,,rrttiiii，,,，，，r，,,,,r，,, (2-67) ,,,,,,,,,,,,,,,,,,iiiiDDDDDD，，,,,,,,,,,,i,1i,1,,,,,,,,,,r,r,tt，，，，iiii，，，，，，，，,,,, 同样地，运动学约束中不涉及时间，时间仅可能出现在驱动约束中，驱动约束方程是 仅依赖于广义坐标的函数之和或仅依赖于时间的函数之和，故,,0。 qt 在计算速度方程和加速度方程中的雅可比矩阵时，并不是进行实时的数值计算，而是 基于具体的约束类型进行计算。不管是运动学约束还是驱动约束，都可分为有限的几种类 型，针对每一种类型的运动副计算其雅可比矩阵的代数形式，如此，在速度分析和加速度 分析时只要先进行雅可比矩阵的组装，然后在迭代的每一时刻代入具体的构件特性值即可。 第2章 多体系统动力学基本理论 2.2.4 多刚体系统动力学 对于受约束的多体系统，其动力学方程是先根据牛顿定理，给出自由物体的变分运动 方程，再运用拉格朗日乘子定理，导出基于约束的多体系统动力学方程。与运动学分析类 似，先考虑二维多体系统，再讨论三维多体系统，并对动力学三种类型的分析：正向动力 学、逆向动力学和静平衡分析分别予以讨论。 1 先给出自由刚体的运动方程，再根据拉格朗日乘子定理给出约束多体系统带乘子的运 动方程，并讨论系统动力学分析的三种情况和约束反力问题。 （1）自由物体的变分运动方程 任意一个刚体构件,m，质量为，对质心的极转动惯量为J，设作用于刚体的所有外iii力向质心简化后得到外力矢量,,,,F和力矩n，若定义刚体连体坐标系的原点位于刚oxoyii 体质心，则可根据牛顿定理导出该刚体带质心坐标的变分运动方程： T，，,，，,r[mr,F]，,,[J,,n],0 (2-68) iiiiiiii 其中，,r为固定于刚体质心的连体坐标系原点的代数矢量，,为连体坐标系相对于全局oii坐标系的转角，,r,,与分别为r与,的变分。 iiii 取2.2.3节为构件定义的广义坐标： i TTq,[r,,] (2-69) iii 定义广义力： TTQ,[F,n] (2-70) iii 及质量矩阵： ,M,diag(m,m,J) (2-71) iiii 则可将式(2.2-61)写作虚功原理的形式： T，，,q(Mq,Q),0 (2-72) iiii 式(2-72)是连体坐标系原点固定于刚体质心时用广义力表示的刚体变分运动方程。其中 广义坐标选取、广义力及质量矩阵计算分别按式(2-69)、(2-70)及(2-71)进行。 （2）束多体系统的运动方程 考虑由nb个构件组成的机械系统，对每个构件运用式(2-72)，组合后可得到系统的变分运动方程为： nbT，，,q[Mq,Q],0 (2-73) ,iiiii,1 若组合所有构件的广义坐标矢量、质量矩阵及广义力矢量，构造系统的广义坐标矢量、 质量矩阵及广义力矢量为： TTTTq,[q,q,...,q] (2-74) nb12 M,diag(M,M,...,M) (2-75) 12nbTTTTQ,[Q,Q,...,Q] (2-76) nb12 系统的变分运动方程则可紧凑地写为： T，，,q[Mq,Q],0 (2-77) 21 机械系统动力学分析及ADAMS应用 对于单个构件，运动方程中的广义力同时包含作用力和约束力，但在一个系统中，若 只考虑理想运动副约束，根据牛顿第三定律，可知作用在系统所有构件上的约束力总虚功 为零，若将作用于系统的广义外力表示为： TTTAAAAT (2-78) Q,[Q,Q,...,Q]nb12 其中： TAAAT， (2-79) i,1,2,...,nbQ,[F,n]ii 则理想约束情况下的系统变分运动方程为： TA，，,q[Mq,Q],0 (2-80) 式中虚位移与作用在系统上的约束是一致的。 ,q 系统运动学约束和驱动约束的组合如式(2-10)，为： (2-81) ,(q,t),0 注意，在动力学分析中系统约束方程的维数不需要与系统广义坐标维数相等。如果令 nmm,nq,R，则，，且。 n,3，nb,,R 对式(2-81)微分得到其变分形式为： ,,q,0 (2-82) q 式(2-80)和(2-82)组成受约束的机械系统的变分运动方程，式(2-80)对所有满足式(2.2-75) 的虚位移均成立。 ,q 为导出约束机械系统变分运动方程易于应用的形式，运用拉格朗日乘子定理对式(2-80) 和(2-82)进行处理。 nnm，n拉格朗日乘子定理：设矢量b,Rx,R，矢量，矩阵为常数矩阵，如果A,R有： Tbx,0 (2-83) 对于所有满足式（2-84）的x条件都成立。 Ax,0 (2-84) m则存在满足式（2-85）的拉格朗日乘子矢量,,R。 TTbx，,Ax,0 (2-85) 其中x为任意的。 m，nnAnn，n在式(2-80)和(2-82)中，q,RQ,R，，，,,R，运用拉格朗日M,Rqm乘子定理于式(2-80)和(2-82)，则存在拉格朗日乘子矢量,,R,q，对于任意的应满足： ATTTAT，，，，[Mq,Q,]q,，,,q,[Mq，,,,Q],q,0 (2-86) qq 由此得到运动方程的拉格朗日乘子形式： TA，，Mq，,,,Q (2-87) q 式(2-87)还必须满足式(2-10)、(2-12)和(2-14)表示的位置约束方程、速度约束方程及加速度约束方程，如下： ,(q,t),0 (2-88) ，，，,,,,(q,t),(q,q,t),,(q,t)q,,,0， (2-89) tq，，，，，，，，，，，,,,(,q)q,2,q,,,(q,q,q,t),,(q,t)q,,(q,q,t),0， (2-90) qqqtttq 以上三式其维数同式(2-81)。 第2章 多体系统动力学基本理论 式(2-87)、(2-88)、(2-89)和(2-90)组成约束机械系统的完整的运动方程。 将式(2-87)与(2-90)联立表示为矩阵形式： TA，，，，,Mq，，Q，，q, (2-91) ,,,,,,,0,,,,，，q，，，， 式(2-91)即为多体系统动力学中最重要的动力学运动方程，被称为欧拉-拉格朗日方程（Euler-Lagrange Equation），式(2-91)还必须满足式(2-88)和(2-89)。它是一个微分-代数方程组（Differential Algebraic Equations - DAEs），不同于单纯的常微分方程组（Ordinary Differential Equations - ODEs）问题，其求解关键在于避免积分过程中的违约现象，此外， 还要注意DAE问题的刚性问题。 显然，式(2.2-84)有且仅有唯一解的充要条件是其系数矩阵非奇异，但这一条件不利于 实际中的判断，可以给出更为实用的判断。 如果式(2.2-84)满足如下条件： a. Rank,(q,t),mm,n，； qTb. 对任意a,Ker,(q,t)aM(q,t)a,0且，。 a,0q 则式(2.2-84)中的系数矩阵是非奇异的，且，，和是唯一确定的。这是,q 。 可以据此判断，如果系统质量矩阵是正定的，并且约束独立，那么运动方程就有唯一 解。实际中的系统质量矩阵通常是正定的，只要保证约束是独立的，运动方程就会有解。 在实际数值迭代求解过程中，需要给定初始条件，包括位置初始条件q(t)和速度初始0条件，q(t)。此时，如果要使运动方程有解，还需要满足初值相容条件，也就是要使位置初0 始条件满足位置约束方程，速度初始条件满足速度约束方程。对于由式(2-91)及(2-88)、(2-89)确定的系统动力学方程，初值相容条件为： ,(q(t),t),0 (2-92) 00 ，，，,(q(t),q(t),t),,(q(t),t)q(t),,(q(t),t),0 (2-93) 000q00000 （3）正向动力学分析、逆向动力学分析与静平衡分析 对于一个确定的约束多体系统，其动力学分析不同于运动学分析，并不需要系统约束 方程的维数mnm,n等于系统广义坐标的维数，。在给定外力的作用下，从初始的位置 和速度，求解满足位置约束式(2-88)及速度约束式(2-89)的运动方程式(2-91)，就可得到系统的加速度和相应的速度、位置响应，以及代表约束反力的拉格朗日乘子，这种已知外力求 运动及约束反力的动力学分析，称为正向动力学分析。 如果约束多体系统约束方程的维数mnm,n与系统广义坐标的维数相等，，也就是对系统施加与系统自由度相等的驱动约束，那么该系统在运动学上就被完全确定，由2.2.3节的约束方程、速度方程和加速度方程可求解系统运动。在此情况下，式(2-88)的雅可比矩阵是非奇异方阵，即： ,(q,t),0 (2-94) q 展开式(2-91)的运动方程，为： TA，，Mq，,,,Q (2-95) q 23 机械系统动力学分析及ADAMS应用 ，，,q,, (2-96) q 由式(2-96)可解得，，，再由式(2-95)可求得，拉格朗日乘子就唯一地确定了作用在,,q 系统上的约束力和力矩（主要存在于运动副中）。这种由确定的运动求系统约束反力的动力 学分析就是逆向动力学分析。 如果一个系统在外力作用下保持静止状态，也就是说，如果： ，，， (2-97) q,q,0 那么，就说该系统处于平衡状态。将式(2-97)代入运动方程式(2-87)，得到平衡方程： TA (2-98) ,,,Qq 由平衡方程式(2-98)及约束方程式(2.2-80)可求出状态和拉格朗日乘子。这种求系,q 统的平衡状态及在平衡状态下的约束反力的动力学分析称为（静）平衡分析。 （4）约束反力 对于约束机械系统中的构件，设其与系统中某构件存在运动学约束或驱动约束，约ij 束编号为,,,。除连体坐标系外，再在构件上以某点为原点建立一个新的固定于构kxoyiP 件上的坐标系,,,,,,,,,,,，称为运动副坐标系，设从坐标系到坐标系的变换矩阵xPyxPyxoy为,,,C，从坐标系到坐标系的变换矩阵为A，则可导出由约束产生的反作用力xoykxoyii 和力矩分别为： TkTTkk,,F,,CA,, (2-99) iiiriTTTkPTkkk,,,T,(sB,,,), (2-100) iiir,iikkk,,,,以上两式中，FT,为约束对应的拉格朗日乘子，反作用力和力矩均为运动副坐标kii系,,,,中的量。 xPy 2 三维系统的广义坐标比二维系统复杂得多，使得问题规模更大。这里讨论的是与二维 多体系统动力学分析相应的内容，包括微分-代数混合方程组的建立，三种类型的动力学分 析等。只是根据三维系统情况的不同，给出了角速度表示和欧拉参数表示的两种不同形式 的运动方程。 （1）空间自由刚体的变分运动方程 对于空间任意刚体构件,,,,,oxyzo，令其连体坐标系原点固定于刚体质心，此时连iiiiii ,,,,体坐标系也称为质心坐标系，设刚体质量为moxyz，其相对于质心坐标系的惯性张量iiiii ,,,,,,为oxyzJFn，再设作用在刚体上的总外力，外力相对于质心坐标系原点的力矩为，iiiiiii则相对于刚体质心坐标系的刚体牛顿-欧拉变分运动方程为： TT~，,,,,,,,，，,r[mr,F]，,,[J,，,J,,n],0 (2-101) iiiiiiiiiii333,其中，r,R,r,R,,,R为刚体质心的虚位移，为刚体的虚转动，为刚体质心iii3,,,,,位移，,,Roxyz为刚体的在坐标系中表示的角速度。 iiiii （2）空间约束机械系统的运动方程—角加速度形式 考虑由nb个刚体组成的空间约束机械系统，如2.2.3节所述，系统的广义坐标选为： TTTTr,[r,r,...,r] (2-102) nb12TTTTp,[p,p,...,p] (2-103) 12nb 定义： 第2章 多体系统动力学基本理论 TTTT,r,,[r,,r,...,,r] (2-104) nb12 M,diag(mI,mI,...,mI) (2-105) 1323nb3TTTTF,[F,F,...,F] (2-106) nb12TTTT,,,,,,,,[,,,,,...,,,] (2-107) nb12 ,,,,J,diag(J,J,...,J) (2-108) 12nbTTTT,,,,,,,[,,,...,,] (2-109) nb12TTTT,,,,,n[n,n,...,n] (2-110） nb12~~~~,,,,,diag(,,,,...,,) (2-111) ,12nb 运用上述符号将系统中每个构件的牛顿-欧拉方程式(2-101)总和为： TT~，,,,,,,,，，,r[Mr,F]，,,[J,，,J,,n],0 (2-112) 同样地，在理想约束情况下，对于一个系统只要考虑外力或外力矩，如此，由式(2-112)可得到约束机械系统的变分运动方程： TATA~，,,,,,,,，，,r[Mr,F]，,,[J,，,J,,n],0 (2-113) 此式必须适用于所有运动学上所容许的虚位移和虚转动。 系统中的运动学约束和驱动约束如式(2-60)及(2.-61)，联立表示为： K，，,(r,p),(r,p,t),,0 (2-114) ,,D,(r,p,t)，， 与式(2-60)和式(2-61)不同的是，动力学约束中允许总体上大于零的自由度出现，设其 维数为。 nh(nh,6nb) 此外，采用欧拉参数还须满足欧拉参数归一化约束（其维数为nb）： pTTTT,,[pp,1,pp,1,...,pp,1],0 (2-115) 112233 对式(2-114)微分，得到用虚位移和虚转动表示的变分形式： ,,,r，,,,,0 (2-116) ,r, ,式(2-116)即为系统虚位移,r和虚转动,,在运动学上所容许的条件。式中不包括式 ,(2-60)所表示的欧拉参数归一化约束，因为如式(2-64)所表明的，当以,为变量时，与欧拉 ,参数归一化约束相关的速度方程是自动得到满足的，同样，以,,为变量时，与欧拉参数归一化约束相关的变分方程也是自动得满足的。 对式(2-113)和(2-116)应用拉格朗日乘子定理，则存在一个拉格朗日乘子矢量,，满足： TATTAT~，,,,,,,,，，,r[Mr,F，,,]，,,[J,，,J,,n，,,],0 (2-117) ,,r ,式中,r,,和是任意的，由此导出空间约束机械系统的牛顿-欧拉运动方程为： TA，，Mr，,,,F (2-118) rTA~，,,,,,,，,,,,J,n,J, (2-119) ,, 对式(2-114)分别求一次及两次导数，得到2.2.3节中所述的系统速度方程及加速度方 程： ,，,r，,,,, (2-120) ,r, ，,，，,，,,,,r .(2-121) ,r, 式中速度右项,和加速度右项与式(2-65)和(2-67)中定义相同。式(2-120)和(2-121)维, 25 机械系统动力学分析及ADAMS应用 数与式(2-114)相同，为。 nh 由式(2-118)、(2-119)及(2-121)可以得到方阵形式的系统运动方程为： TA，，0MrF，，，，,，，r,T,,A,~,,，0,,,JnJ (2-122) ,,,,,,,,,,,,,,,,,,,,,0,,r,,,,,，，,，，，， 式(2-122)与约束方程式(2-114)及速度方程式(2-120)一起组成描述nb个刚体系统运动的微分-代数方程组。式(2-122)的维数关系满足： 。 [(6nb，nh)，(6nb，nh)],[(6nb，nh)，1],[(6nb，nh)，1] 由于角速度是不可积的，因此应把式(2-122)看作是速度变量及代数变量的一阶微分- 代数混合方程。 与二维约束机械系统类似，如果满足： a. 矩阵[,,,]行满秩； ,r,T,b. 对任意a[M,J]a,0a,Ker[,,,]且，存在。 a,0,r, 则式(2.2-115)中系数矩阵非奇异，系统运动方程有唯一解。 在实际迭代计算中，需要给定初始位置条件，r(t)p(t)r(t)和，以及初始速度条件和000 ,,(t)，初始条件须分别满足位置约束方程和速度约束方程： 0 ,(r,p,t),0 (2-123) 0 ,，,r(t)，,,(t),,(r,p,t) (2-124) ,r0,00 （3）空间约束机械系统的运动方程—欧拉参数形式 运用2.2.3节所介绍的欧拉参数导数与角速度之间的关系，可以得出用欧拉参数导数表 示的空间机械系统运动方程。 在2.2.3节中欧拉参数导数与角速度关系基础上，可进一步导出欧拉参数二阶导数与角 加速度关系，对于单个刚体的这些运动量间关系，在这里一并给出： ,，,,2Gp (2-125) 1T,，p,G, (2-126) 2 ，,，，,,2Gp (2-127) 11TT，,,,，，p,G,,,,p (2-128) 24 ,,2,G (2-129) ,p, 则用欧拉参数表示式(2-120)所示的系统速度方程为： ，，,r，,p,, (2-130） rp 用欧拉参数表示式(2-121)所示的系统加速度方程为： ，，，，,r，,p,, (2-131) rp 式（2-130）和（2-131）维数同样为nh。 对欧拉参数式(2-115)微分，得到欧拉参数速度方程： p，,p,0 (2-132) p 第2章 多体系统动力学基本理论 其维数关系为，且其中： [nb，4nb],[4nb，1],[nb，1] pTTT (2-133) ,,2diag(p,p,...,p)p12nb 对式(2-132)微分，得到欧拉参数加速度方程： pTTTTp，，，，，，，， (2-134) ,p,,2[pp,pp,...pp],,p123 其维数关系同式(2-132)。 约束方程变分式(2-116)用欧拉参数表示为： ,,r，,,p,0 (2-135) rp 欧拉参数约束方程(2-115)的变分式为： p (2-136) ,,p,0p 对于用角速度表示的约束机械系统的变分运动方程式(2-113)，根据式(2-150) 替换 ，,,,，根据式(2-125)替换，根据式(2-127)替换，并利用2.2.3节系列公式，可得到用,,,, 欧拉参数表示的系统变分运动方程： TATTTTA，，,,,，，，，,r[Mr,F]，,p[4GJGp,8GJGp,2Gn],0 (2-137) 其中： G,diag(G,G,...,G) (2-138) 12nb 式(2-137)必须对满足式(2-135)和(2-136)的所有和成立。 ,r,p p应用拉格朗日乘子定理于式(2-137)及(2-135)和(2-136)，则有拉格朗日乘子,和使下,式对任意和成立。 ,r,pTTATTTTTATpp，，,,,，，，，,rM[r,F，,,]，,p[4GJGp,8GJGp,2Gn，,,，,,],0rpp (2-139) 从而得到欧拉参数形式的约束机械系统变分运动方程： TA，，Mr，,,,F (2-140) rTTTppTAT，，,,,，，4GJGp，,,，,,,2Gn，8GJGp (2-141) pp 由式(2-140)和(2-141)，再加上运动学约束和驱动约束的加速度方程式(2-131)及欧拉参数约束的加速度方程式(2-134)，可得到方阵形式的空间约束机械系统运动方程为： TA，，，，M00,r，，F，，rT,,,,,,TTpTAT，，,，，04GJG,,p,,28GnGJGp，pp,,,,,, (2-142) ,,,,,,,00,,,,rp,,,,,,ppp000,,,,,,,p，，，，，， 式(2-142)与式(2-114)和式(2-115)的运动学约束和欧拉参数归一化约束以及式(2-130)和式(2-132)的相应速度方程一起，构成欧拉参数形式表达的系统微分-代数运动方程的混合方程组。其维数关系为： [(8nb，nh)，(8nb，nh)],[(8nb，nh)，1],[(8nb，nh)，1] 式(2-142)的系数矩阵非奇异的条件同本节二维多体系统的运动方程解存在条件。 实际求解中需要给出初始位置，，r(t)p(t)r(t)p(t)和及初始速度和，且初始位置条0000件必须满足位置约束方程，初始速度条件必须满足速度约束方程，如下： 27 机械系统动力学分析及ADAMS应用 ,(r,p,t),0,0 .(2-143) ,ppt,(,),00, ，，,rtptrpt,()，,(),(,,),rp000 (2-144) ,p，pt,(),0p0, 与角速度形式的运动方程不同，式(2-142)、(2-114)、(2-115)、(2-130)和(2-132)组成是二阶微分-代数混合方程组。 在角速度形式的运动方程式(2-122)中，系数矩阵左上方是常数矩阵,和，而在欧JM T,拉参数形式的运动方程式中，系数矩阵中出现了取决于欧拉坐标的项GJG，这会使两种形式的运动方程的数值性态不尽相同。 （4）逆向动力学分析、平衡分析和运动副约束反力 空间约束机械系统的逆向动力学分析、平衡分析与与平面约束机械系统的相应方法完 全相同。 逆向动力学分析，对于运动学上确定的系统，式(2-122)或(2-142)的系数矩阵非奇异，可以直接求出加速度和拉格朗日乘子，再进一步得到系统运动状态和约束反力。 平衡分析，根据平衡状态的定义由式(2-118)和(2-119)或由式(2-140)和(2-141)可以导出系统平衡方程，再由平衡方程及约束方程可以求出系统的平衡位置及平衡时的拉格朗日乘 子，再进一步得到约束反力。 运动副约束反力的计算也与平面约束机械系统类似。考虑一典型运动副，运动副定k kk义点为,,0,，约束方程为，相应的拉格朗日乘子为，运动副反作用力和反作用力P kkk,,,,矩在运动副定义坐标系,,,,,,FT,Pxyz中表示为和，则由拉格朗日乘子计算运动副反ii 作用力和反作用力矩的公式为： TkTTkk,,F,,CA,, (2-145) iiiriTTkTkPTkk~,,,T,,C(,,sA,), (2-146) ,ii,iirii 式中,,,,,,,,,,C为从运动副定义坐标系Pxyz到连体坐标系oxyz的方向余弦变换矩阵。 i 2.2.5 计算多刚体系统动力学自动建模 系统的力学模型是对实际问题的力学抽象，要进行动力学的求解，需要由系统的力学 模型得到系统的数学模型，也就是要得到系统运动方程(2.2-80)～(2.2-83)或(2.2-115)、(2.2-107)和(2.2-113)的具体形式，这其中的关键在于组装系统运动方程中所有的系数矩阵。 计算多体系统动力学是基于约束的运动学和动力学，不仅指运动的速度方程和加速度 方程是在约束方程的基础上建立，动力学的运动方程在约束方程的约束下形成微分-代数方程，也指在多体系统动力学分析过程中，系统运动方程的各种导数不是实时地采用求导算 法进行计算，而是采用基于约束的计算方法。所谓基于约束的计算方法，是指对于有限的 约束类型，包括运动学约束和驱动约束，针对每一种约束计算出在系统运动方程中所需的 各种导数的相应的代数形式，然后在建立数学模型时组装成系统运动方程中各种导数的组 合式，如此，在计算导数时只需代入广义坐标、时间及其它相关参数即可，避免了导数实 第2章 多体系统动力学基本理论 时计算所花的大量费用。 根据系统运动方程(2-87)～(2-90)的要求，对于所有的各种运动学约束和驱动约束，需 要给出如下几个量的代数形式：约束方程,(q,t)、约束方程雅可比矩阵、速度右,(q,t)q 项及加速度右项，所有基本约束相应参量的表示可以直接给出，这里不再,(q,t),(q,v,t) 详细列出。 现在讨论从系统力学模型到系统数学模型建立过程中的自动建模技术，主要是各种系 数矩阵的组装技术。 如式(2-87)～(2-90)所述，系统动力学运动方程为： T，，，，，， (2-147) F(q,q,,q,,t),M(q,t)q，,(q,t,),f(q,q,t),0q (2-148) ,(q,t),0 ，，，,,,,(q,t)， (2-149) ,(q,q,t),,(q,t)q,,(q,t),0tq，，，，，，，，，，，,,,(,q)q,2,q,,， (2-150) ,(q,q,q,t),,(q,t)q,,(q,q,t),0qqqtttq 则需要组装的矩阵有：,(q,t)、、、、和。f(q,v,t),(q,t),(q,t)M(q,t),(q,v,t)q 如果考虑到求解过程中所涉及到的一些运算，还需给出所有约束的其它更复杂参量的代数 T形式，如：[,(q,t)a][,(q,t)v],(q,v,t),(q,v,t)、、、及等。这些形[,(q,t),]qqqvq 式都可以采用本节的统一方法一致获得。 以二维系统为例，三维系统同理。设二维系统由个构件和个约束（包括运动学nbnj约束和驱动约束）组成，并对系统中构件施加力（矩）和力元的作用，再设系统广义坐标 维数为nmr，约化后的约束方程维数为，定义构件的物理参量和约束（设其联接构件jis与）的运动参量如下： 广义坐标 q i 质量 m(q,t) ii ,极转动惯量 J(q,t) ii 外力矢 F(q,t) ii 外力矩 n(q,t) iij约束方程 ,(q,t) r,sj速度右项 ,(q,t) r,sj加速度右项 ,(q,v,t) r,sr,s ,其中，构件质量和极转动惯量考虑的是一般情况，极转动惯量J是对构件质心的极转动惯i量；外力矢和外力矩是由作用于构件的任意时变外力和力矩向质心简化后所得的合力矢和 合力矩，为时间和构件广义坐标函数。组合构件i的质量矩阵和外力及外力矩矢量为 ,质量矩阵 M(q,t),diag(m,m,J) iiiiiiTT外力（矩）矢量 f(q,t),[F,n] iii 在力学模型运动学约束、驱动约束、力元、外力（矩）等要素中，运动学约束和驱动 约束通过代数等价转化为系统约束方程，外力（矩）通过规范化过程转化为系统外力矢量， 至于力元，其对系统的作用可以等价于作用于其所连接的两个构件上的互反的广义力矢量， 在系统数学模型中的体现是并入到系统外力矢量。故此，系统力矢量是由外力矢量和力元 广义力矢量组成。若设力元srk是构件与构件之间的作用，则其等价的广义力为 29 机械系统动力学分析及ADAMS应用 kr力元广义力r fe(q,v,t)r,sr,skksr力元广义力s fe(q,v,t),,fersrs,, 构件上的力元广义力应等于等价于其上的所有力元广义力的合力，设系统中力元个i 数为。 nf nfkil构件上力元广义力 ife,fe,kli,1,, 再定义并计算与构件、约束和力元相关的雅可比阵为 jki j雅可比阵 ， ,(q,t)i,r,sqr,si 由系统中每个构件的参量矩阵，可组装得到系统的参量矩阵。 n1. 系统广义坐标矢量q,R组装（，） n,3nbq TTTTq,[q,q,...,q] (2-151) nb12n，n2. 系统质量矩阵组装（） M,RM M,diag(M,M,...,M) (2-152) 12nbn3. 系统力矢量f,R组装（） f(q,v,t) 系统力矢量由每个构件的外力矢量和力元广义力组成。 TTTTTT12nbT12nbT (2-153) f,[f,f,...,f]，[fe,fe,...,fe] m4. 系统约束方程矢量组装（，） ,(q,t)m,n,,RTTT12njT (2-154) ,,[,,,,...,,] m5. 系统速度右项矢量,,R,(q,t)组装（） TTT12njT (2-155) ,,[,,,,...,,] m6. 系统加速度右项矢量,,R组装（） ,(q,v,t)TTT12njT (2-156) ,,[,,,,...,,] m，n7. 系统雅可比矩阵,(q,t)组装（,,R） qqj因为约束的约束方程,(q,t)只是所联接的两构件的广义坐标和时间的函数，因jr,s 而系统约束雅可比阵为稀疏矩阵，但不是分块对角的。采用分块贡献法进行组装，就是以 约束为考虑对象，分别计算每个约束对所连接的两个构件广义坐标的导数子矩阵，然后加 入到系统约束对系统广义坐标的导数全矩阵中。设系统约束个数为nj，定义矩阵扩展运算 m,nu,v符号u，vm，next(sub)sub，表示将子矩阵扩展为全矩阵，置子矩阵于指定位置i,j (i,j)，其余补充元素全置零，据此可以写出系统雅可比矩阵组装公式。 njmm,3,3jjmnjmnj,, (2-157) ,,[ext(,)，ext(,)],qrcqscq,,,,rsj,1 j,1 其中，,,m，1mkrc,3(r,1)，1，表示第个约束的约束方程数，，k,kk,1 msc,3(s,1)，1j，表示第个约束的约束方程数。 j 第2章 多体系统动力学基本理论 各种参量矩阵的组装是根据力学模型来进行的，得到需要组装的所有参量矩阵之后， 也就得到了系统数学模型，下一步的工作是选用专门的动力学数值分析算法进行迭代求解， 并以适当方式显示分析结果。 2.3 多柔体系统动力学建模 2.3.1 柔性体上点的位置向量、速度和加速度 1 PBP iui Ps iPs i0Prib e Ci rir e 图2.5 柔性体上节点的位置 P rbee柔性体系统中的坐标系如图2.5所示，包括惯性坐标系（）和动坐标系（）。前 者不随时间而变化，后者是建立在柔性体上，用于描述柔性体的运动。动坐标系可以相对 惯性坐标系进行有限的移动和转动。动坐标系在惯性坐标系中的坐标（移动、转动）称为 参考坐标。 与刚体不同，柔性体是变形体，体内各点的相对位置时时刻刻都在变化，只靠动坐标 系不能准确描述该柔性体在惯性坐标系中的位置，因此，引入弹性坐标来描述柔性体上各 点相对动坐标系统的变形。这样柔性体上任一点的运动就是动坐标系的“刚性”运动与弹 性变形的合成运动。由于柔体上各点之间有相对运动，所以动坐标系的选择不是采用连体 坐标系，而需要采用随着柔性体形变而变化的坐标系，即“浮动坐标系”。 在研究多柔体系统时，合适的坐标系是非常重要的。在确定浮动坐标系时有两点准则： 1、便于方程建立求解；2、柔性体刚体运动与变形运动的耦合尽量小。目前常见的浮动坐 标系大致有如下5种，局部附着框架、中心惯性主轴框架、蒂斯拉德框架、巴克凯恩斯框 架以及刚体模态框架。采用何种需因实际情况而定。 31 机械系统动力学分析及ADAMS应用 2 在分析刚体平面运动的时候，把复杂的刚体平面运动分解为几种简单的运动。在对柔 性体的运动，尤其是在小变形的情况下，也可以采用类似的方法。如某柔性体从位置L1 运动到位置L2，其间运动可以分解为：刚性移动－>刚性转动－>变形运动。对于柔性体上任意一点P，其位置向量为： （2-158） rrAsu,，，()0pp 为P点在惯性坐标系中的向量；r为浮动坐标系原点在惯性坐标系中的向量；根据r0 式2.2- 14，A为方向余弦矩阵；为柔性体未变形时P点在浮动坐标系中的向量;为相supp 对变形向量，可以用不同的方法离散化，与讨论平面问题相同，对于点，该单元的变uPp 形采用模态坐标来描述，有： （2-159） uq,,ppf 式（2-159）中，Φ为点满足里兹基向量要求的假设变形模态矩阵，为变形的广qPp 义坐标。 柔性体上任一点的速度向量及加速度向量可以对式求对时间一阶导数和二阶导数得 到： Prr,，，，,AsuAq() （2-160） pppf0 prrAsuAqAq,，，，,，,()2 （2-161） pppfpf0 2.3.2 多柔体系统动力学方程的建立 u PPs Pb()e,P B rB r()e G 图2.6 柔性体变形模型 第2章 多体系统动力学基本理论 1 在ADAMS软件中，外加载荷包括单点力与扭矩、分布式载荷以及残余载荷三部分。 （1）单点力与扭矩 施加于柔性体上某一标记点的单点力和扭矩必须投影到系统的广义坐标上才能起作 用，力和扭矩以矩阵形式写出，在标记点的局部坐标系下表示为： K TT， （2-162） Ffff,[]Tttt,[]KxyzKxyz 广义力由广义平动力、广义扭矩（以欧拉角表示的广义力）和广义模态力组成，可Q 表示为： TQQQQ,[] （2-163） TRM r平动坐标下的广义力可以通过转换单点力e到全局坐标基下来获得，即： FK QAF, （2-164） TK 其中，为标记点上的坐标系相对于全局坐标的欧拉角变换矩阵，可以表示为： AK AAAA,A，如图2.6所示，为局部坐标系点相对于全局坐标系的转换矩阵，Bgkgbbppkgb 即方向余弦阵；A为因节点的小变形引起的标记坐标的方位变化而引入的转换矩阵；Pbp A为定义在柔性体上标记点处的坐标系相对于点坐标系的常值变换矩阵。 Ppk 作用在柔性体标记点处的合力矩，可用相对于全局坐标的矢量矩阵表达为： TATpAF,，， （2-165） totgkKgkK b e为从局部坐标基到受力点处的位移矢量，用斜方阵表示，上式可改写为： p TATpAF,， （2-166） totgkKgkK 将物理坐标系下的力矩向基于欧拉角的坐标系转换，并利用式,,,B，可得广义合B 力矩为： TQABATpAF,，[][] （2-167） RgkKgkK 通过投影单点力和单点力矩到模态坐标上，可得到点处的广义模态力。 P 如在标记点FAF,TF处施加力和力矩，通过转换到全局坐标上，即有：，KIgkKKK 33 机械系统动力学分析及ADAMS应用 TAT,，并将F投影到平动模态坐标、将投影到角模态坐标上，可得广义模态 TIgkKII 力为： TT*QFT,,，, （2-168） FPIPI *其中，,,和对应于节点处的平动和转动自由度的模态斜方阵。由于模态矩阵P,PP 只定义在节点，故单点力和单点力矩只能施加于节点处。 （2）分布式载荷 在ADAMS软件上，分布式载荷可以通过MFORCE的方式来创建。通常在FEM软件中有运动方程： （2-169） MxKxF，, 其中，x和为柔性体上有限单元的质量和刚度矩阵，而和为物理节点的自由MKF度矢量和载荷矢量。 利用模态矩阵将式（2-169）转换到模态坐标下，有： qΦ TTT,,，,,,,MqKqF （2-170） 上式可简化为： ˆˆMqkqf，, （2-171） ˆˆ其中，MK和分别为广义质量和刚度矩阵，为模态载荷矢量。 f 节点力矢量在模态坐标上的投影为： TfF,, （2-172） 上式中，如果是时间的函数则求解时计算开销太大，一种替代的方法是假设空间依F 赖性和时间依赖性可以分开，把载荷看成是一系列依赖于时间的静态载荷的线性组合，F 即： FtstFstF()()(),，， （2-173） 11nn 载荷在模态坐标上的投影计算可在有限元文件MNF过程中完成，而不必在ADAMS仿真中重复进行。如果定义一系列静载荷为载荷矢量，并使其与系统响应显性相关，即表 示成fqt(,)的形式，则模态力又可表示为： fqtsqtfsqtf(,)(,)(,),，， （2-174） nnn1 （3）残余载荷 先前的假设是施加的载荷向模态坐标的投影已实现，即式（2-172）成立。然而由于模 第2章 多体系统动力学基本理论 态截断，很多情况下施加的力并不能进行投影，这些力被称为残余载荷，可以将其看为已 投影到了被截断的高阶模态坐标上。残余载荷可表示为： T,1,,,,FFf[] （2-175） 与残余载荷相关的是残余矢量，可看成是把残余载荷施加于柔性体时产生的变形。残 余矢量可被当作模态振型加入到Craig-Bampton模态基中，增强的模态基完全能够捕捉残余载荷，否则，残余力的丢失不可避免。 2 （1）动能和质量矩阵 考虑节点变形前后的位置、方向和模态，柔性体的广义坐标可以表示为： P TTξrψq,,,[(1,,)][]xyzqiM,,, （2-176） i 速度表达式（2-160）在系统广义坐标式（2-176）的时间导数ξ中表示为： （2-177） vIAsuBA,,，,[()],PppP 柔性体的动能为： 11TTGBTGBTvvdVmvvI,,，,,, （2-178） PPPPPP,,VP22 GB其中，,,,BmI和分别为节点的节点质量和节点惯性张量；，为点相对PBPPPP 于全局坐标基的角速度在局部坐标基中的斜方阵表示。将式（2.3-20）和关系式,,,BPP 代入式（2-178），得到动能的广义表达式： 1TTM,,,,() （2-179） 2 上式中的质量矩阵33，M()ξ为维的方阵，表示为： MMM，，tttrtmT,,MMMM(),, （2-180） trrrrm,,TTtmrmmm,,MMM，， 其中下标trm,,分别表示平动、旋转和模态自由度。质量矩阵的六个独立分量分别表 示为： 35 机械系统动力学分析及ADAMS应用 1,MI,,tt23,MAqB,,，[],,trjj,3,tmMA,,,7889TT （2-181） ,rrjjjijijMBqqqB,,，,[[]],,,,45,Trmjj,MBq,，[],,6,mmM,,,, 其中9个惯性时不变矩阵列表如下： 表2-1 惯性时不变矩阵列表 惯性时不变矩阵 模态坐标数 维数 N 1,,m P标量 ,P1, N 2τs,m PP, (31)，P1, N 3τΦ,m jPP, jM,1,,(3)，MP1, N 4,,,，,,msI PpPPP, (3)，MP1, N 5,,,,m jPpjP,jM,1,, (3)，MP1, N 6TT,,,,,，,,,mI PPPPPP,()MM， P1, N 7T,,，mssI PppP,(33)， P1, N 8,,,ms jPppj,jM,1,, (33)， P1, NjkM,1,,,(33)， 9τ,,,m jkPpjpk,P1, 式(2.3-24)中可以明显看出质量矩阵与模态坐标显性相关，而且由于引入转换矩阵和A 19,~τ，质量矩阵也与系统的方向坐标显性相关。质量矩阵中的9个惯性时不变矩阵可B 通过计算有限元模型的N个节点信息在预处理过程中一次性得到，从而简化运动微分方程 第2章 多体系统动力学基本理论 的求解。节点信息包括：每个节点的质量，未变形时的位置矢量以及模态矩阵。 msΦPPP （2）势能和刚度矩阵 势能一般分为重力势能和弹性势能两部分，可用下列二次项表示： 1T （2-182） WWK,，(),,,g2 在弹性势能中，是对应于模态坐标的结构部件的广义刚度矩阵，通常为常量。重qK 力势能W表示为： g T （2-183） WrgdWrAsqgdW,,,，，,,,[()]gPBPP,,WW 其中fW表示重力加速度矢量，重力可对求导得； ggg ，，[],dWg,W,, ,W,,,AgTfs,,，[()],ΦqgdW （2-184） ,,gPP,W,,ξξ,, T,,A[]Φg,dWP,,,，，W （3）能量损失和阻尼矩阵 阻尼力依赖于广义模态速度并可以从下列二次项中推导得出： 1T,,qDq （2-185） 2 上式称为Rayleigh能量损耗函数。矩阵d包含阻尼系数，它是常值对称阵。当引Dij入正交模态振型时，阻尼矩阵可用对角线为模态阻尼率c的对角阵来表示。对于每一个正i cr交模态的阻尼率都可以取不同值，而且还能以该模态的临界阻尼c的比值形式给出。 i 3 柔性体的运动方程从下列拉格朗日方程导出： dLL,,,,,Ψ,T()[],，，,,,Q0,dt,,ξξ （2-186） ,,ξξ, ,Ψ=0, 其中，为约束方程； Ψ ,为对应于约束方程的拉氏乘子； 37 机械系统动力学分析及ADAMS应用 为如式（2.3-19）定义的广义坐标； ξ 为投影到上的广义力； Qξ 为拉格朗日项，定义为，和分别表示动能和势能； LTW,,WTL 表示能量损耗函数。 , 将求得的代入式（2.3-29），得到最终的运动微分方程为： TW,,, 1,,,MTTMMKfDQ,,,,,,,，,，，，，,[][] （2-187） g,,2,, 其中，ξ,ξ,ξ为柔性体的广义坐标及其时间导数； MM,为柔性体的质量矩阵及其对时间的导数； ,M 为质量矩阵对柔性体广义坐标的偏导数，它是一个,ξ 维张量，为模态数。 (6)(6)(6)MMM，，，，，M 2.4 多体系统动力学方程的求解 如2.2.2节所述，在建模和求解过程中，涉及到几种类型的运算和求解：初始条件计算、 数学模型自动组装、运动学分析、动力学分析、逆动力学分析和静平衡分析。初始条件计 算是非线性位置方程的求解；数学建模是系数矩阵操作；运动学分析是非线性的位置方程 和线性的速度、加速度方程的求解；动力学分析是二阶微分方程或二阶微分方程和代数方 程混合问题的求解；逆向动力学分析是线性代数方程组求解；静平衡分析从理论上是线性 方程组的求解。总的来说，计算多体系统动力学涉及的基本运算包括线性方程组求解、非 线性方程组求解、常微分方程组（ODE）求解和微分代数方程组（DAE）求解。 线性方程组和常微分方程组的求解是数值分析中的基本内容。线性方程组根据问题规 模可以采用消元法或迭代法，在计算多体系统动力学中，一般采用全主元高斯消元法；常 微分方程组的求解可以采用线性多步法或单步的龙格-库塔法。 在这里先简要地给出求解非线性方程组的牛顿-拉夫森（Newton-Raphson）方法，再详细地介绍计算多体系统动力学中至关重要的微分代数方程组求解技术。 2.4.1 非线性代数方程组求解 这里介绍应用广泛的求解非线性代数方程组的牛顿-拉夫森方法。考虑含有n个变量的 第2章 多体系统动力学基本理论 n个非线性方程组： (2-188) ,(q),0TT其中，；；0为n维列向量。 ,,,,q,q,q,...,q,,,,,,...,,nn1212*i*i令q,q为式（2.3-1）的解，记q为q的近似，将在q,q进行一阶泰勒展开,(q)为： iiii (2-189) ,(q),,(q)，,(q)(q,q)，o(q,q)q 其中,为雅可比矩阵，见2.2.3之第2部分说明。 qi，1ii，1令q,q(q,q)为改进的近似，如果足够小，则有： i，1iii，1i (2-190) ,(q),,(q)，,(q)(q,q),0qii，1i令,q,q,q，式（2.4-3）化为： iii (2-191) ,(q),q,,,(q)qi如果雅可比阵非奇异，则式（2-191）有唯一解，且得到： ,(q)qi，1iiq,q，,q (2-192) 0于是可根据式（2.4-4）与（2.4-5），由某一初始值q开始，逐步迭代，直到收敛，得到有效结果，或者不能收敛，初始值不合适或牛顿-拉夫森方法失效。 2.4.2 微分代数方程组求解 在2.1.3节第2小节中，综述了目前微分代数方程求解的常用技术，根据对广义坐标处 理方式的不同，可以分为增广法和缩并法。缩并法一般是在处理过程中将微分代数方程转 化为二阶常微分方程再进行求解，增广法更多的是将微分代数方程直接转化为代数方程组 的求解序列。 在80年代后期至90年代，微分代数方程求解中缩并法占了主要地位，但90年代增广法展现出新的面貌，体现出其独特的优势。这里介绍F. A. Potra提出的基于超定微分代数方程组（ODAEs - Overdetermined Differential Algebraic Equations）方法的解耦ODAE方法，虽然其数学处理比较复杂，但理解起来是比较简单的，它主要采用微分流形的技术，将微 分代数方程组转化为代数方程组的求解序列。该方法充分体现了求解DAE问题中所采用预估-校正策略和显式-隐式方法的区别。 1 式(2-87)~(2-90)给出了约束多体系统动力学方程的完整形式，在这里重写为： T，，，，，，F(q,q,,q,,t),M(q,t)q，,(q,t,),f(q,q,t),0 (2-193) q ,(q,t),0 (2-194) ，，，,,,,(q,t),(q,q,t),,(q,t)q,,(q,t),0， (2-195) tq，，，，，，，，，，，,,,(,q)q,2,q,,,(q,q,q,t),,(q,t)q,,(q,q,t),0， qqqtttq (2-196) nm其中，，，，q,R,,Rq、、分别为系统位置、速度、加速度向量，是拉格朗日乘q m，nn，n子向量，M(q,t),Rt,R是时间，为机械系统质量矩阵，为,(q,t),,,/,q,RqAnm，，约束雅可比矩阵，f(q,q,t),Q(q,q,t),R,(q,t),R为外力向量，为位置约束方程，,和分别称为速度右项和加速度右项。 , 39 机械系统动力学分析及ADAMS应用 系统运动方程(2-193)～(2-196)被称为欧拉-拉格朗日（Euler-Langrange）方程，它是微分-代数方程（DAE），DAE不是单纯的常微分方程（ODE）问题，其求解的关键在于避免积分过程中的违约现象。 上述的系统运动方程给出的是最一般的情况。虽然是以二维系统的形式给出，但三维 系统如果对各系数矩阵进行适当的组装，可以得到与上述几式一致的运动方程，只是其中 各量的含义与二维不同。对于质量矩阵，考虑的是与位置和时间相关的质量矩阵，这是最 一般的情况；对于外力向量，考虑的是与位置、速度和时间相关的力，即考虑到了弹簧力、 阻尼力等情况；对于约束方程，同时考虑到运动学约束和驱动约束。上述各式给出了明确 的函数相关量。 对于上述各式，为了解的存在，一般作如2.2.4节讨论的解的存在性定理中所述的假设： a. Rank,(q,t),mm,n，； qTb. 对任意a,Ker,(q,t)aM(q,t)a,0且，。 a,0q 2-ODAEs 引入速度矢量va和加速度矢量，记： ，，， (2-197） a,vv,q 则式(2-193)～(2-196)可以另写为： ， (2-198) q,v,0 ， (2-199) v,a,0 T (2-200) F(q,v,,a,,t),M(q,t)a，,(q,t),,f(q,v,t),0q (2-201) ,(q,t),0 ， (2-202) ,(q,v,t),,(q,t)v,,(q,t),0q ，， (2-203) ,(q,v,a,t),,(q,t)a,,(q,v,t),0q 定义增广广义坐标： TTTTTx,[q,v,a,,] (2-204) 在方程组(2-198)～(2-203)中，增广广义坐标x维数为3n+m，方程组维数为3n+3m，这 是一个超定问题，即所谓超定微分-代数方程组（ODAEs），但在前述假设条件下它是相容有解的。 由于系统内在的相容一致性，上述方程组中虽然存在冗余，但可以采取某种方式消除 这种冗余性。根据解的存在性定理中的假设，可以认为冗余不会存在于式(2.4-13)～(2.4-16) 中，只会存在于式(2-198)和(2-199)中，对式(2-198)和(2-199)采用微分流形“投影”技术以消除超定性。 T，A(v,a),0 (2-205) 1T，A(q,v),0 (2-206) 2n，(n,m)式中，A,RAAA和。如果取得有效的和，式式(2-198)和(2-199)从n维降到（n-m）1212维，消除了系统冗余性。可证明，如果取AA和满足： 12 T，，Ai(i,1,2)， (2-207) ,,qt,(,),,qcc，， 第2章 多体系统动力学基本理论 非奇异，则方程组(2-205)、(2-206)、(2-200)～(2-203)与方程组(2-208)～(2-213)在某点(q,t)cc 的某个邻域内同解。的一种较为有效的取法是： AiTT A,K(q,t)M(q,t) (2-208) iccccn，(n,m)其中满足： K,R ,(q,t)K(q,t),0,qcccc (2-209) ,T(,)(,)KqtKqt,Icccc, 作如上处理后，方程组(2-198)～(2-203)可写作： TF(x,t)，，M(q,t)a，,(q,t),f(q,v,t),，，q,,,,，，,(x,t),(q,t)a,(q,v,t),q,,,,，,,,(x,t),,,(q,t)v,,(q,t)qH(x,t),,,0 (2-210) ,,,,1,(x,t),(q,t),,,,TT,,,,AR(x)，A(v,a)111,,,,TT，AR(x)A(q,v),,22，，2，， 3n，m3n，m此时，x,RH,R，，消除了系统超定性。 1 采用线性多步法： rr ，y,,y，h,y (2-211) ,,nin,iin,i10i,i, 令 rr ˆq,,q，h,v (2-212) ,,nin,iin,i11i,i, rr ˆv,,v，h,a (2-213) ,,nin,iin,i11i,i, 在avt,tqva,H(x,t),时，将、、和分别写作变量、、和，则可化为： qnnnnn1 FxtMqta,qtfqvtT(,)，，(,)，,(,),(,,)，，nnqnn,,,,xtqt,aqvt，，,(,),(,),(,,)n,,qnn,,xtqt,vqt，,,,,,(,),(,),(,)nqnnHxt(,),,,0 (2-214) ,,,,nxtqt,(,),(,)nn,,,,TTARxAv,hav,,,,()ˆ11(,,)10n,,T,,TARxAqh,vq2220n()ˆ(,,),,,,，，，，3ODAEs 考虑显式方法。采用显式线性多步法时，式(2-211)中,,0。 0 41 机械系统动力学分析及ADAMS应用 在时，取值为： t,tAn2TTˆˆA,K(q,t)M(q,t) (2-215) 2nnnn 将式(2-214)中第6式和第4式组成的位置方程组： T,ˆAqq(,),02n (2-216) ,qt,(,),0n, m应用式（2-216）第一式和式（2-209）第一式，引入新的未知因子,,R，可有： 1TTTˆˆˆ (2-217) A(q,q)，K(q,t),(q,t),,0nnnqnn21 对式(2-217)应用式（2-215）和式（2-209）第二式，可得： Tˆˆˆ (2-218) M(q,t)(q,q)，,(q,t),,0nnnqnn1 式（2-218）称为式（2-216）第一式的扩展等价式。将式（2.4-31）替换到式（2.4-29） 中，可得： T,,ˆˆˆ,(,)(,)，,(,)Mqtqqqt1nnnqnnˆ,,,(,),,0Pq (2-219) 1n,,,(,)qtn,,由式(2-219)计算出q后，再取值为： An1TTA,K(q,t)M(q,t) (2-220) 1nnnn 将式(2-214)中第5式和第3式组成的速度方程组： T,ˆA(v,v),0,1n (2-221) ,qtv,qt,(,),(,),0,qnnnn, 同式（2-219）的处理，应用式(2-209)将式（2-221）改写为扩展的等价式： T,,ˆMqtvvqt,(,)(,)，,(,)nnnqnn2,,V(v,,),,0 (2-222) n2,,qtv,qt,(,),(,)qnnnn,, m式中,,R为新引入的未知因子。 2 由式(2-222)计算出v后，根据由式(2-214)中第1式和第2式组成的方程组： n TMqtaqt,fqvt,,(,)，,(,),(,,)nnqnnnnn,,Aa,(,),,0 (2-223) n,,qta,qvt,(,),(,,)qnnnnn,,可计算出加速度a,和拉格朗日乘子。 nn 由式(2-219)求解q时，由于式(2-219)是非线性方程，运用Newton-Raphson方法。将n 式(2-219)应用Newton-Raphson方法： k，，qknkˆˆˆJ(q,t),,,P(q,,)k,0,1,...， (2-224) ,,nnnnn1kn,1n，， 其中： 第2章 多体系统动力学基本理论 T，，ˆˆM(q,t),(q,t)nnqnnˆˆJ(q,t), (2-225) ,,nnnˆqt,(,)0,,qnn，， 1kk，k，，，，，，qqqnnn (2-226) ,,,,,,,,,kk，1k,,,111nnn，，，，，， 迭代初始值 00ˆq,q， (2-227) ,,0nn1n 应用式(2-224)进行迭代，如果收敛则可计算出。 qn由式(2-222)计算时， 同样应用Newton-Raphson方法： vn k，，vknkˆˆJ(q,t),,,V(v,,)，k,0,1,... (2-228) ,,nnnnn2kn,2n，， ˆ其中，ˆ表示同式(2-225)。 J(q,t)nnn 1kk，k，，，，，，vvvnnn,,, (2-229) ,,,,,,kk，1k,,,222nnn，，，，，， 迭代初始值： 00ˆv,v， (2-230) ,,0nn2n 由式(2-223)计算a,和时，应用Newton-Raphson方法： nn k，，akknˆˆJ(q,t),,,A(a,,)k,0,1,...， (2-231) ,,nnnnnnk,n，， ˆ其中，ˆ表示同式(2-225）。 J(q,t)nnn kk，1k，，，，，，aaannn,,, (2-232) ,,,,,,kk，1k,,,nnn，，，，，， 迭代初始值： 00a,a,,,， (2-233) nn,1nn,1 由上所述，可以写出基于ODAEs的显式算法。 ODAEs 1. 时间t,t(n,0)qv，给定相容初始位置和初始速度，由式(2-223)计算初始加000 速度a,和初始拉格朗日乘子； 00 2. 进入时间步t(n,1)； n 3. 由式(2-212)和(2-213)计算ˆˆqv和； nn4. 由式(2-219)计算q。利用式(2-224)～(2-227)所述的Newton-Raphson迭代方法； n 5. 由式(2-222)计算v。利用式(2-228)～(2-230)所述的Newton-Raphson迭代方法； n 6. 由式(2-223)计算a,和。利用式(2-231)～(2-233)所述的Newton-Raphson迭代方nn 43 机械系统动力学分析及ADAMS应用 法； 7. 如果到终止时间步，结束；否则进行下一时间步，转3。 t,tnn，14ODAEs- 考虑隐式方法。采用隐式线性多步法时，式(2-211)中。 ,,00TT取AA和为： 12TTTA,A,K(q,t)M(q,t) (2-234) 12nnnn 则系统运动方程式(2-214)可改写为扩展的等价式为： Mqtaq,tfqvtT，，(,)，,(,),(,,)nqnn,,qt,aqvt,(,),(,,),,qnn qt,vqt,,,(,),(,)qnnH,x,t (,,,),,0 (2-235) ,,12nqt,(,)n,,TMqtvh,avqt,,,ˆ(,)(,,)，,(,)02nnnqnn,,TMqtqh,vqqt,01nnnqnnˆ(,)(,,)，,(,),,，， m式中,,R,和为新引入的未知因子。 21m用式(2-209)易证明，当且仅当有唯一确定的,,,,R(q,v,a,,,,,,)使是nn12nnnn12nn 式(2.4-48)的解时，(q,v,a,,)为式(2-214)的解。 nnnn 式(2-235)为隐式方程组，现采用预估-校正方法求解。 引入预估式： rr q,,q，h,v (2-236) ,,nin,iin,i11i,i, rr v,,v，h,a (2-237) ,,nin,iin,i11i,i, 上述两式为预估式，不同于式(2-212)和(2-213)所表示的线性多步法的显式部分和。 TTAA12利用预估式估计和，式(2-234)则为： TTTA,A,K(q,t)M(q,t) (2-238) 12nnnn 如此，形如式(2-219)、(2-222)及(2-223)的系统位置、速度和加速度方程则相应地为： T,,ˆ,,(,)(,,)，,(,)Mqtqhvqqt01nnnqnn,,,(,,),,0Pqv (2-239) 1n,,,(,)qtn,, T,,ˆMqtv,havqt,(,)(,,)，,(,)nn0nqnn2,,V(v,a,,),,0 (2-240) n2,,qtv,qt,(,),(,)qnnnn,, TMqtaqt,fqvt,,(,)，,(,),(,,)nqnn,,Aqva,(,,,),,0 (2-241) n,,qta,qvt,(,),(,,)qnn,, 第2章 多体系统动力学基本理论 据此写出基于ODAEs的隐式算法。 ODAEs- 1. 时间，给定相容初始位置和初始速度，由式(2-241)计算初始加t,t(n,0)qv000速度和初始拉格朗日乘子； a,00 2. 进入时间步； t(n,1)n 3. 由式(2-212) 、(2-213) 、(2-236)和(2-237)分别计算ˆˆ、、和，并设 qvqvnnnn000000q,q，v,v，，，a,a，,,, ； (2-242) ,,0,,0nnnnnn,1nn,112nn 4. 对于，计算 k,0,1,... k，1k，，，，qqknnkkJ(,),,P(q,v,,) (2-243) ,,,,nnnn1k，1kn,,11nn，，，， k1k，，，，，vvknnk，1kkJ(,),,V(q,v,a,,) (2-244) ,,,,nnnnn2k1k，n,,22nn，，，， k，1k，，，，aak，1k，1kknnJ(,),,A(q,v,a,,) (2-245) ,,,,nnnnnnk，1k,,nn，，，， 式中： T，，M(q,t),(q,t)nnqnnJJqt,(,), ； (2-246) ,,nnnqt,(,)0qnn，， 5. 如果到终止时间步，结束；否则进行下一时间步t,t，转3。 nn，1 2.5 多体系统动力学中的刚性（Stiff）问题 刚性问题存在于多刚体系统动力学某些情形，更普遍地存在于多柔体系统动力学中， 是多体系统动力学的一个重要问题。刚性首先是在常微分方程求解理论中提出，并形成了 完整的定义和求解理论。常微分方程刚性理论是多体系统动力学中刚性问题的理论基础， 这里先介绍常微分方程刚性问题，再讨论多体系统动力学刚性问题。 2.5.1 微分方程刚性（Stiff）问题 微分方程的刚性（Stiff）问题是微分方程的一个重要问题，微分-代数方程（DAE）中 同样存在刚性问题。微分-代数方程早期的数值求解中，并没有考虑到这个问题，采用的大 多是显式方法，到了二十世纪八十年代，发现一些隐式方法不仅具有更好的适应性，而且 可用于求解刚性问题。 1 为描述刚性方程性质，先考虑线性系统： 45 机械系统动力学分析及ADAMS应用 dy， (2-247) ,Ay(t)，,(t)t,[a,b]dt TmTm其中y(t),[y(t),...,y(t)],R,(t),,[(t),...,,(t)],R为解向量函数，为已m11mm，m知向量函数，,,,，i,为常系数矩阵，设其特征值为，，相应j,1,...,mA,Rjjj的特征向量为,。 j 定义（2-247） 如果在线性系统(2-248)中，,(j,1,...,m)的特征值满足： Aj a. Re(,),0，； j,1,...,mj b. s,max|Re(,)|/min|Re(,)|,,1jj1,j,m1,j,m p则称(2-247)为刚性方程，比值ssO(10)(p,1)称为刚性比，通常刚性比达到就认为 是刚性的。 对于非线性系统： dy,f(t,y) (2-248) dt 令y(a),y为(2-248)满足初始条件的精确解，在解的邻域内对方程(2.5-2)y(t)y(t)0 作线性逼近： dy,J(t)[y,y(t)]，f(t,y(t)) (2-249) dt 或 dy,J(t)y，F(t) (2-250) dt 其中是在点处的雅可比（Jacobi）矩阵的值，如果矩J(t)(t,y(t))f(t,y),f(t,y)/,y阵,,,(t)(j,1,...,m)J(t)的特征值满足： jj a. Re(,),0，j,1,...,m； j b. st(),max|Re(,)|/min|Re(,)|,,1jj1,j,m1,j,m 则称方程(2-248)为刚性方程，s(t)称为在处的局部刚性比。 t 在刚性方程中，刚性比s,,1J(t)，矩阵或是病态的，故刚性方程也称为病态方程A 或坏条件方程。刚性方程数值积分过程中存在快变分量max|Re(,)|和慢变分量jmin|Re(,)|的差别，快变分量要求积分步长很小，而慢变分量则使得在该步长条件下计j 算步数很多，舍入误差较大，这就是刚性方程数值解的困难所在。 考虑到实际问题中可能会出现单个方程情形，或者矩阵的特征值有实部为零或实部A为很小正数的情形，可以给出与实际问题更为接近的刚性方程的定义。 定义2-248 若线性系统(2-249)满足条件： 1. 矩阵的所有特征值实部小于不大的正数； A 2. 至少有一个特征值的实部是很大的负数； A 3. 对应于具有最大负实部的特征值的解分量变化是缓慢的。 则称(2.5-1)是刚性方程。 第2章 多体系统动力学基本理论 对于刚性方程数值稳定性的讨论，一般是针对试验方程： dy， (2-251) ,,yy(0),y0dt 进行的。 定义2-249 一个数值方法以定步长h解试验方程(2-251)，得到线性差分方程的解yn， 当n,,时，若y,0，则称该方法对步长h是绝对稳定的。 n 定义2-250 一个数值方法称为A稳定的，如果它的绝对稳定域包含整个左半平面 。 Re(h,),0 对于A稳定，存在如下结论： 1. 任何显式线性多步法（包括显式RK方法）不可能是A稳定的。 2. A稳定的隐式线性多步法的阶不超过2。 3. 具有最小误差常数的2阶A稳定隐式线性多步法是梯形法。 A稳定的条件过于苛刻，满足条件的数值方法太少，为了突破这个限制，可以放宽稳 定性条件，给出稳定的定义。再从刚性问题特点出发，给出刚性稳定性定义。 A(,) 定义2-251 一个数值方法称为稳定的，如果它的绝对稳定域包含无限的楔形区A(,) 域W,{h,|,,,,,arg(h,),,}，。 ,,(0,,/2), 定义2-252 一个数值方法称为刚性稳定的，如果它是收敛的，并存在正常数使D,,,, 在区域R,{h,|Re(h,),,D}是绝对稳定的，而在区域1 R,{h,|,D,Re(h,),,,|Im(h,)|,,}上具有高精度且是绝对稳定或相对稳定的。 2 刚性稳定性的定义，是考虑到微分方程(2.5-1)或(2.5-2)的解中具有较大时间常数的项对解有重大影响。 2 对常微分方程组初值问题： ,y,f(t,y),t,[a,b], (2-252) ,y(a)y,0, TmTm其中y,[y,...,y],Rf,[f,...,f],R为解向量，为已知向量函数，mm11my,R为已知初始向量。 0 对常微分方程组初值问题，常用的方法有三种。 （1）线性多步法（LMM） kk y,,y，h,f(t,y) (2-253) ,,,,,niniinini,1,0ii 当,,0,,0k,1k,1时为显式公式，当时为隐式公式。当时称为单步法，当时00 称为多步法。 （2）预估校正法（PECE） kk y,,y，h,f(t,y) (2-254) ,,,,,niniinini,1,1ii 47 机械系统动力学分析及ADAMS应用 kk (2-255) y,,y，,hf(t,y)，h,f(t,y),,,0,,nininninini,1,1ii 式(2-254)为显式预估公式，式(2-255)为隐式校正公式。 （3）龙格-库塔法(RKM) s (2-256) y,y，hbknnii,,1i,1 s k,f(t，ch,y，hak),i,1,...,s (2-257) ininijj,,1,1j,1 如果当a,0a,0时，式(2-256)和(2-257)是显式RK方法。如果当时，而对j,ij,iijij角元素a不全为0，式(2-256)和(2-257)称为半隐式RK方法，若此时a均相等，则称为对iiii角隐式RK方法，简称为DIRK方法。 在这些求解常微分方程初值问题数值方法的基础上，产生了求解刚性常微分方程的几 类方法，分别是以BDF方法为代表的线性多步法和隐式龙格-库塔方法。BDF方法是一类特殊的隐式线性多步法。一般的隐式龙格-库塔方法计算量较大，这里只给出一类特殊的隐 式龙格-库塔方法：对角隐式RK方法。 （1）向后差分公式（BDF） k y,,y，h,f(t,y) (2-258) ,ni,ni0nn,i1 BDF方法是隐式线性多步法，为k步k阶方法。当k=1,2时，BDF方法是A稳定的，当k=3～6时，BDF方法是稳定和刚性稳定的。实用的BDF方法只能取k=1～6。 A(,) （2）对角隐式龙格-库塔方法（DIRK） s y,y，hbk (2-259) nnii,,1i,1 i,1 k,f(t，ch,y，hak，hak),i,1,...,s (2-260) in,1in,1ijjiii,j,1 a,a,...,a (2-261) 1122ss 对角隐式龙格-库塔方法常用的有2级3阶和3级4阶两个公式，都是A稳定的，此外还有A稳定的4级4阶公式，但给不出A稳定的更高阶DIRK公式。DIRK方法解高频振荡的问题比Gear方法（即BDF方法）好，但对高精度问题比不上Gear方法好，且计算量比Gear方法大。 第2章 多体系统动力学基本理论 2.5.2 多体系统动力学中Stiff问题 在多刚体系统运动过程中，可能会由于系统中构件之间过大的差异，如不同物体特性 参数的差异，力元（如弹簧、阻尼器）参数的差异等，致使系统中构件运动速度差别很大， 从而使描述系统运动的微分代数方程呈现出刚性特性。在多柔体（或刚-柔混合）系统动力学中，由于柔性体空间大范围运动和其本身小幅值弹性变形的耦合，更容易出现刚性问题。 图2.7 图2.8 如图2.7所示k1＝10 N/mm，k2＝1000000 N/mm，m1=1kg，m2=1kg，计算可知系统的第一阶频率为11.25Hz，第二阶频率为7117.6Hz，为典型的物理刚性系统。如图2.8所示加入阻尼c1＝0.1 N/mm/s，c2＝100000N/mm/s，系统变为典型的数值刚性系统。物理刚性 系统一定存在大刚度弹簧，系统具有高频；数值刚性系统除在刚度方面存在较大差异外， 还有一个很重要的特征是对应高频的阻尼较大，使较高频率被基本阻尼掉，而低频则处于 未阻尼状态。数值刚性系统用非刚性数值算法会出现数值困难，这时候需采用专用于求解 刚性问题的数值方法。 微分-代数方程的求解，无论是缩并法还是增广法，问题还都是归结为常微分方程初值 问题的数值求解，只是求解常微分方程的公式或是用于微分-代数方程转化为常微分方程之 49 机械系统动力学分析及ADAMS应用 后，或是用于转化过程之中。为了使求解的数值方法具有普遍性，既可用于求解良性问题， 也可用于求解刚性问题，微分-代数方程所用的常微分方程数值方法一般采用的都是上节所 述的求解刚性微分方程的方法。除了最常用的BDF方法，隐式RK方法也被考虑用于求解微分-代数方程问题，此外，预估-校正方法也广泛地用于求解微分-代数方程问题之中，它们都可求解存在刚性问题的微分-代数方程。