（机械制造及其自动化专业论文）多柔体系统动力学数值方法研究.pdf

摘要 多柔体系统动力学模型建立与数值解法是多体系统动力学研究的主要内容之 一，其中动力学数学模型可分为两大类，第一类是单纯的二阶微分方程：第二类 是微分代数方程组。至于第一类方程的解法己发展成熟，然而第二类方程组的解 法目前仍然处于探索阶段，不够完善，属于疑难问题。 本文基于上述背景，以多体系统动力学理论为基础，首先较为全面的介绍了 多柔体系统动力学的各种数值方法，重点论述了微分代数方程组的几种违约校正 法，分析其优缺点及适用范围，提出了一种新的算法：即迭代稳定违约校正法。 该方法从b a u m g a r t e 违约稳定法出发，在每一时间步，利用二阶方程的直接积分法 计算迭代初值，根据这个初值，应用中心差分法反复迭代约束稳定法的稳定方程， 然后联立动力学控制方程和n e w t o n r a p h s o n 迭代公式，得出方程的修正解，保证 了约束方程和系统的动力学控制方程的同时满足。最后，给出一典型算例，应用 m a t l a b 6 5 软件进行数值计算仿真，结果验证了本文方法有效性。 关键词：多柔体系统动力学数值方法微分一代数方程组稳定方程 a b s t r a c t m o d e l i n ga n dn u m e r i c a ls o l u t i o n so ff l e x i b l em u l t i b o d yd y n a m i c sa reo n eo f m a j o rc o n t e n t s m a t h e m a t i c a lm o d e lo ff l e x i b l em u l t i b o d yd y n a m i c si n c l u d e s t w o k i n d so ft h e m f i r s ti so r d i n a r yd i f f e r e n t i a le q u a t i o n s e c o n di sd i f f e r e n t i a la l g e b r a i c e q u a t i o n s w i t ht h ed e v e l o p m e n to ff i r s tk i n do fe q u a t i o n ，i th a si n c r e a s e df a s t l y b u t n u m e r i c a ls o l u t i o no fs e c o n dk i n do fe q u a t i o nh a ss t i l lb e e nr e s e a r c h e d t h i sq u e s t i o ni s d i f f i c u l tt os o l v ea tp r e s e n t i nt h i sp a p e r , b a s e do nt h ef o l l o w i n gp a s s a g e ，w i t ht h et h e o r yo ff l e x i b l e m u l t i b o d yd y n a m i c s ，v a r i o u sn u m e r i c a lm e t h o d sf o rd y n a m i c so ff l e x i b l em u l t i b o d y s y s t e m sa r ed i s c u s s e da n dt h e nf o c u so nt h ed i s c u s s i o no fs o m ec o r r e c t i o na l g o r i t h mo f c o n s t r a i n tv i o l a t i o no fd i f f e r e n t i a l a l g e b r a i c e q u a t i o n s ，t h ea d v a n t a g e s a n d d i s a d v a n t a g e s ，a p p l i e dr a n g eo ft h o s em e t h o d sa r ea n a l y z e d an e wm e t h o df o r c o n s t r a i n ts t a b i l i z a t i o ni s p r e s e n t e d i ti s c a l l e dc o r r e c t i o na l g o r i t h mo fi n t e r a t i v e c o n s t r a i n tv i o l a t i o n a te a c ht i m es t e p ，t h ei n i t i a li t e r a t i v ev a l u e sa r ee v a l u a t e db yd i r e c t i n t e g r a lm e t h o do fs e c o n d - o r d e re q u a t i o n ，t h e nn e wg o u po fi t e r a t i v ev a l u e sa r e e v a l u a t e db yc e n t e rd i f f e r e n c em e t h o da n ds t a b i l i z a t i o ne q u a t i o no fc o n s t r a i n t s t a b i l i z a t i o nm e t h o dw i t ht h ef o r m e ri n i t i a lv a l u e s s oan e wg r o u po fc o r r e c t i o nv a l u e s a r ev a l u e db ya p p l y i n gf o rt h en e w t o n r a p h s o ni t e r a t i v ef o r m u l aa n di n t e g r a t i n g c o n t r o le q u a t i o no fm u l t i b o d yd y n a m i c s i tc a nm e e tt h ec o n s t r a i n ta n dc o n t r o l e q u a t i o n so fm u l t i b o d ys y s t e m ss i m u l t a n e o u s l y f i n a l l y , w i t ht h er e s u l t so fo n et y p i c a l e x a m p l e , t h ee f f e c t i v e n e s so fn e w m e t h o di sp r o v e db yt h es o f t w a r em a t l a b 6 5 k e y w o r d ：d y n a m i c so ff l e x i b l em u l t i b o d ys y s t e m s n u m e r i c a lm e t h o d d i f f e r e n t i a la l g e b r a i ce q u a t i o n s s t a b i l i z a t i o ne q u a t i o n 西安电子科技大学 学位论文独创性( 或创新性) 声明 秉承学校严谨的学风和优良的科学道德，本人声明所呈交的论文是我个人在 导师指导下进行的研究工作及取得的研究成果。尽我所知，除了文中特别加以标 注和致谢中所罗列的内容以外，论文中不包含其他人已经发表或撰写过的研究成 果；也不包含为获得西安电子科技大学或其它教育机构的学位或证书而使用过的 材料。与我一同工作的同志对本研究所做的任何贡献均已在论文中做了明确的说 明并表示了谢意。 申请学位论文与资料若有不实之处，本人承担一切的法律责任。 本人签名： 西安电子科技大学 关于论文使用授权的说明 本人完全了解西安电子科技大学有关保留和使用学位论文的规定，即：研究 生在校攻读学位期间论文工作的知识产权单位属西安电子科技大学。学校有权保 留送交论文的复印件，允许查阅和借阅论文：学校可以公布论文的全部或部分内 容，可以允许采用影印、缩印或其它复制手段保存论文。同时本人保证，毕业后 结合学位论文研究课题再撰写的文章一律署名单位为西安电子科技大学。 ( 保密的论文在解密后遵守此规定) 本人签名：兰= ! 二塑 导师签名：刻盟适 日期 k o ；乜 第一章绪论 第一章绪论帚一早珀y 匕 多体系统是指由多个刚体或刚体与弹性体通过一定的方式相互连接构成的复 杂机械系统。多刚体系统动力学是以系统中各部件均抽象为刚体，但可以计及各 部件联结点( 关节点) 处的弹性、阻尼等影响为其分析模型的。多刚体系统动力学是 在经典力学的基础上发展起来的，从6 0 年代至今，多体系统动力学建立模型已经 形成了许多各具特色的方法，如n e w t o n e u l e r t 2 】方法、l a g r a n g e t 2 】方法、 r o b e r s o n w i t t e n b u r g t l 】方法、k a n e 1 】方法、h u s t o n 1 】方法等，并取得了相当完善的 成果，已经解决了工程上的许多实际应用课题。 另外一个重要方面即动力学方程求解，首先动力学数学模型分为两类，一种 是低阶的纯微分方程，一种是微分代数方程组。近二十年来国内外进行了大量的 研究工作，前者的研究算法已经基本成熟，而后者微分代数方程组的解法仍然没 有得到很好解决。目前这类方程组的研究方法大体可分为两类【1 1 ：一种是从微分 代数方程组本身出发，利用现代数学的研究成果将约束定义为流形，对微分代数方 程组进行降阶处理，将其转化为由约束方程定义的流形上的常微分方程。这种方 法的优点是可以直接应用求解常微分方程的技术，避免约束方程违约，但在求解 过程中必须计算由约束方程定义的流形零空间的基，计算工作量大，对复杂的多 体系统，零空间基的计算缺乏成熟的方法，且有时并不唯一；另一种方法是在动力学 方程中引入附加校正项，当约束方程产生违约时，对动力学方程进行校正【4 】。目前 的校正方法多为间接校正方法，不能对系统的广义坐标进行直接的校正以满足约 束方程。还有，在动力学方程中加入附加校正项需给定校正系数，校正系数太小 校正效果不明显，校正系数太大容易引起动力学方程的破坏。到现在还没有校正 系数的自动选取方法，大都凭经验选取校正系数。 1 1多柔体系统动力学国内外数值计算方法研究概况 1 1 1刚性微分方程的数值解法研究进展 多体系统的动力学方程是强耦合、强非线性方程，这种方程的求解目前只能 通过计算机用数值方法解决。慢变的、大的刚性运动和快变的、小的弹性运动的 耦合使得柔性多体系统的动力学方程成为时变的非线性的刚性方程，其求解过程 易出现病态或不收敛，故多柔体系统的动力学方程求解引起学者们的极大关注。 一般数值解法借助计算机分析，复杂系统计算机辅助分析( c a a ) 是计算机 2多柔体系统动力学数值方法研究 辅助工程( c a e ) 的重要环节，并且日益收到国内外工程界的重视，在机械设计 领域，设计者已经习惯于在生产样机前用多体软件对产品进行仿真，分析和优化， 大大缩短了生产周期，降低了产品的生产成本，多体系统动力学c a a 包括： 1 ) 计算机辅助建立系统运动学与动力学数学模型； 2 ) 数学模型的计算机自动分析； 3 ) 数值结果的自动解释。 经过2 0 多年的发展，多体系统的数学模型问题已经基本解决，根据系统的拓 扑结构及描述系统构形所用坐标，多体系统动力学方程可分为以下两类【3 】： 第一类数学模型： a ( q ，f ) 牙= b ( q ，口，f ) ( 1 - 1 ) 其中，彳，b ，口分别为系统的广义质量阵，广义力列阵，广义坐标矢量，该式 可直接转化为如下状态方程 夕= f ( y ，t ) ( 1 - 2 ) 其中，y = g r ，口7 】r ，f = 口7 ( 彳叫b ) 7 r ，式( 1 1 ) 数值积分方法的研究基本成 熟，目前对于解这类方程的数值算法主要有：牛顿一拉斐逊( n e w t o n r a p h s o n ) 法、中 心差分法、威尔逊0 法、纽马克( n e w m a r k ) 法、帕克( p a r k ) 刚性稳定法、龙格库 塔法( r u n g e k u t a ) 等。由于该类方程的解中包含快变分量和慢变分量，给其数值 计算带来了很大的困难。常用的显示积分法大都是条件稳定的，不适用于求解此 类方程。目前刚性常微分方程初值问题数值计算方法绝大多数是隐式算法，每积 分一步都要用n e w t o n 法解2 n 维非线性代数方程组，因此要计算方程j a c o b i 矩阵，当 方程组的维数n 较大时候，j a c o b i 矩阵的推导工作量都是相当巨大的。虽然通过离 散n e w t o n 法( n e w t o n s t e f f o n s e n ) ，用差商代替偏导数，可避免雅可比矩阵，提高 该矩阵的计算效率，但是数值稳定性方面不够完善。为提高数值计算的稳定性， 克服刚性方程在数值计算中的困难，隐式数值计算方法逐步被采用。一般常用的 算法有龙格库塔法、线性多部法、g e a r 法等。 因此，刚性方程的数值求解一直是计算数学、一般力学和其他工程领域学者 们广泛关注的研究热点。在国内有大连理工大学的钟万勰教授等提出的精细时程 积分法【l4 1 ，其关键是矩阵指数函数在计算机字长范围内精确计算的特点，精细时 程积分法是一种绝对稳定的显示积分法，可以用它求解刚性方程。若是在采用精 细积分算法时候选用积分精度较高的科茨积分格式，对一些动力学问题将会得到 非常精确的数值结果。 针对钟万勰提出的精细积分法【1 9 】，吕和祥等人提出的逐步积分法，借助于线 性动力学方程奇次解的解析表达式，构造出适用于非线性动力学的积分方程，其 方程的积分核为一半的正弦和余弦函数【2 5 1 。对于包含非线性项的非齐次项采用了 显示预测矫正四阶精度插值的算法，其转换矩阵也是解析的，因而不受步长大小 第一章绪论 的限制，其计算也不需要精细积分算法中的2 技巧，同时也减少了许多精细积分 算法中矩阵与向量相乘，该方法适用于强非线性、非保守动力学系统。 1 1 2 微分代数方程组的数值解法研究进展 上节介绍了刚性纯微分方程的表达式以及研究进展，本节主要讨论微分一代数 方程组的表达式以及研究发展状况。第二类动力学方程数学模型t 3 】： m ( q ，f ) 尊+ 幻+ c ；( g ，f ) a = q + g ( 1 - 3 ) c ( q ，f ) = 0 ( 1 4 ) 其中m 为系统广义质量矩阵是关于q 和t 复合函数，q 为系统广义坐标矢量，是 关于t 的函数，五为l a g r a n g e 乘子矢量，q 为雅可比矩阵，绋广义外力，g 速度 二次项有关的广义力，这些参数的具体表达式我们将在后面的章节中具体说明。 这是一个典型的微分代数混合方程组，即d a e ( d i f f e r e n t i a la l g e b r a i c e q u a t i o n s ) ，它与常微分方程不同，存在着数值求解上的困难，这也是本文的数值 方法的具体研究对象。 对于微分代数方程组的求解，从一开始就是关注的焦点。1 9 5 0 年h o u b o l t 首先提出了基于l a g r a n g e 插值公式的步进方法【2 1 。1 9 5 9 年n e w m a r k 给出了在位 移、速度的近似关系中引入两个参数的多步隐式方法f l o 】。1 9 6 6 年j o h s o n 利用单 位圆域变换到z 平面左半平面的方法证明了h o u b o l t 法是无条件稳定的】。1 9 7 3 年w i l s o n 提出了引进一族积分参数臼的单步隐式法】。1 9 7 7 年z i e n k i e w i c z 提出 了一类基于广义加权残值方法的多步时域方法，w o o d 对n = 2 和n = 3 的z i e n k i e w i c z 算法进行了比较，同年h i b l e r 等对线性隐式算法做了系统研究【1 2 】。1 9 8 0 年g r o w n 提出了一种新的p e c 算澍27 1 。1 9 8 9 年c a r d o n a 和g 6 r a d i n 将运动方程映射到方程 考虑的某些项的正切空间上完成积分【1 3 】。1 9 9 1 年p e r e i r a 提出了改进的递推方法。 e n r i g h t 精选了2 5 个题目对不同算法的费用和可靠性进行了比较，h u s t o n 在1 9 9 1 年认为四阶r u n g e k u t t a 法是最好的方法【1 2 】。1 9 9 1 年b a y o 给出了隐式梯形公式 的计算方法和所用的j a c o b i 矩阵的具体表达式，便于数值计算。1 9 9 2 年s i m o 提 出能量贮存( e n e r g e - c o n s e r v i n g ) 法，在1 9 9 5 年又提出能量贮存修正 ( e n e r g e c o n s e r v i n gm o d i f i c a t i o n ) 法【2 7 1 。1 9 9 4 年，y o o n 提出了广义坐标主动校正 方法，即通过对约束可展开天线中柔性多体系统动力学d a e 的数值方法研究方 程进行隐式迭代来减少约束方程的违约，动力学方程的数值计算采用显式算法 a s k o h 给出了向后插分格式的计算方法。2 0 0 1 年i b r a h i m b e g o v i c 等人对非线性 柔性多体系统动力学的计算方法与应用做了研究【2 7 1 。国内自2 0 世纪8 0 年以来不 断出现一些与算法研究相关的文献，1 9 8 0 年任曾勋提出了一种工程解法，通过引 4 多柔体系统动力学数值方法研究 入补充未知量构成双直交基来满足求解条件【5 1 。1 9 8 1 年孙焕纯提出了一个改进的 0 法。同年蔡承文证明了当0 0 ) ，且记g j = q ( t o + 访) ( i = 0 ，1 ，2 ，) 。设 精确解g ，的近似解为玩，一般c ( 玩) = 占0 ，对c ( 玩+ 。) t y l o r 展开得： c ( 玩+ i ) = c ( q i ) + h c ( q i ) + ( h 2 2 ) e ( 玩) + d ( 3 ) ( 2 5 1 ) 若近似解玩满足： c + ( 2 h ) c + ( 2 h 2 ) c = 0 ( 2 5 2 ) 则恒有： c ( 玩+ 1 ) = o ( h 3 ) ( 2 - 5 3 ) 即( 2 5 1 ) 式有对约束方程误差进行自我修正的能力，使得数值积分过程具 有好的稳定性，从而提高计算精度。由此可得： 多柔体系统动力学数值方法研究 一q q , r + q , ：- 一k q ( 2 h ) 6 ( 2 h ：，q ) c 2 5 4 ， 【一q c 一 2 2 ) s ij 一。 i 口= 1 h t ：压 ( 2 - 5 5 ) 即给出了校正系数的自动选取方法。此方法的优点是避免了选择口和的盲 目性，数值计算稳定性较高，缺点是对于各种变步长的积分程序，很难确定调用 修正方程右端函数时的当前积分步长。 2 4 5 其它的微分代数方程组的数值方法 约束违约自动修正法 1 9 9 9 年，孔向东，钟万勰提出了一种非线性系统动力学微分代数方程组约束 违约的自动修正法【1 9 】。该方法运用最优控制思想，推导出了约束误差可达计算机 有效精度i 拘b a u m g a r t e 最佳违约修正系数。 一般地，系统动力学方程为： m ( q ，f ) 辱+ 五= f ( q ，毒，f ) ( 2 5 6 ) c ( q ，f ) = 0 ( 2 5 7 ) 约束方程对时间求导得： cq白+ct=0(2-58) c q 自= 一( c q 4 ) q 自一2 c q i 自一c “( 2 - 5 9 ) 设约束方程的位移速度违约值分别为e o ，毛，即c ( q ，f ) = 气，c ( q ，香，f ) = g 。 根据b a u m g a r t e 最佳违约修正技术，由式( 2 6 0 ) 、( 2 6 1 ) 、( 2 6 2 ) 组成的闭环系 统可写为： c ( q ，雪，f ) + a l c ( q ，口，f ) + a o c ( q ，f ) = 0 ( 2 6 0 ) 其中： a o , 口1 r 小。埘；a o = d i a g ( a ? ) ，a 1 = d i a g ( a ：) ，口量0 o ，a ：0 ( k = l ，”，m ) ，贝0 ： c ( q ，雪，f ) = 一a 0 6 0 一口l g l ( 2 - 6 1 ) 由此可将方程( 2 5 6 ) 和( 2 5 9 ) 化成微分方程组的形式： 香= 1 ， ( 2 6 2 a ) = g 兄 0 o得 m q 河i较 0 与 第二章多柔体系统动力学数值方法 1 9 痧= 六+ w ( a o s o + 口l 占1 ) ( 2 6 2 b ) 其中： w = - m 一喏( q m 叫) 。 ( 2 6 3 ) 正= m f 一 c q m 一1 厂+ ( c g e l ) ge l + 2 c 口+ c “】 ( 2 - 6 4 ) 显然当= o 及g ，= 0 时，式( 2 6 5 ) 与直接积分法的方程组等价。 令a t 为步长，t 洲= t ，+ a t ( f - o ，1 ，n 1 ) 。下面确定【f j ，t 圳 时段的a o , a l 值。设 未加入控制量的状态方程为： “= 移：矿= l(265q 26 5 ) 2 v ；，2 。 l j 根据式( 2 6 2 b ) 得 帝= 铲+ w ( a o s o + 口l s i ) ( 2 - 6 6 ) 在【f ，t m 】时段对式( 2 6 6 ) 积分，得到 1 ，= 矿+ w ( a o g o + 口1 9 1 ) f l a t ( 2 - 6 7 ) 由积分中值定理可知，0 六 l 。又由式( 2 6 5 ) 及式( 2 6 7 ) 可得： =+w(aqq ( a o e o + 口l 占1 ) 卣f 2+ + 口l 占lj 告1 f 同理可得： q = 虿+ w ( a o s o + 口l q ) 彘鼻f 2 其中：o 4 0 l 。采用隐式积分格式，且氢= 孝。= 0 5 ，在t 川时刻有： ( 2 - 6 8 ) ( 2 - 6 9 ) q f + i = 玩+ l + 【( 口o o + a l e l ) 】f + la t 2 4 ( 2 - 7 0 ) h + i = 巧+ l + 【矽( 口。占o + a l 毛) 】“l f 2 ( 2 - 7 1 ) 又在，t 】时段，采用e u l e r 修正法得： c i + i = c + a t ( d i 4 l + e ) 2 ( 2 - 7 2 ) 在初始时刻f 。= o ，c o = c o = c o = 0 。从最优控制角度出发，如能控制口o ，口l ， 使f o ，t l 一，t 。各个时刻都有q o ，e 0 ，则由式( 2 6 0 ) 可知奠0 。这种在各 个时刻的位置和速度违约误差皆近似为0 的状态正是此方法要达到的最理想目 标，其相应的控制量a 。，a 。即为最优控制量。 现在假设t ，时刻c o = c o = 0 ，则f 时刻有： 2 0多柔体系统动力学数值方法研究 由此得： c j + l = a t c i + l 2 ( 岛) f + i = a t ( q ) f + i 2 q m = 玩+ l + ( 口o a t 2 4 + a l a t 2 ) 6 0 f + l 1 ，川= 玩+ i + ( 口o a t 2 4 + a l a t 2 ) e 1 j + l 令u o l = d i a g ( a ? a t 2 4 + a ：a t 2 ) ( 尼= 1 ，聊) ，u o l r 刷脚，则 q f + l = 氟l + ( w u o i ) + l ，+ i = 玩+ l + ( w u o l 占1 ) n 1 ( 2 7 3 ) ( 2 7 4 ) ( 2 7 5 a ) ( 2 7 5 b ) ( 2 7 6 a ) ( 2 7 6 b ) 上两式即为式( 2 6 2 ) 的动力学方程的解，它使得两组控制量a o , a 。被一组控制量 u 。所代替。为了求解简便，高效，将式( 2 - 7 6 ) 的隐式格式写成显示迭代格式。 q f + i = 玩+ i + ( g u o i ) l + l ( 瓦) m