（机械制造及其自动化专业论文）求解多柔体系统动力方程的违约修正零空间法.pdf

摘要 建模方法和数值求解方法是多柔体系统动力学两大核心内容，许多的系统控 制问题的数学模型可以用一个微分一代数方程组来描述，微分一代数方程的数值求 解方法已成为近年来的研究热点与难点。 本文主要对数值求解方法进行了论述和研究。首先，以多体系统动力学理论 为基础，介绍了多体系统动力学的各种数值方法，分析了其优缺点及适用范围。 详细介绍了几种基于广义坐标的独立性分解的缩并算法，给出了各种矩阵分解方 法在微分代数方程中的应用，并介绍了零空间法的思想。然后在缩并法的基础上 结合违约修正的思想，提出了一种违约修正零空间法。此方法用种简便零空间 的求解方法，分离出独立广义坐标，把微分代数方程化为纯微分方程，在每一时 刻步，利用纽马克法计算迭代初值，代入违约修正公式，判断违约，进行修正。 最后，通过一个典型算例说明本文算法的稳定性与有效性。 关键词：多柔体系统动力学独立坐标数值方法零空间约束违约零空间法 a b s t r a c t m o d e l i n ga n d n u m e r i c a ls o l u t i o n sa r et h ed y n a m i c so ff l e x i b l e m u l t i b o d ys y s t e m s k e r n e lc o n t e n t s m a n ys y s t e mc o n t r o lp r o b l e mc a r lb ed e s c r i b e dw i t hd i f f e r e n t i a l a l g e b r a i ce q u a t i o n s ，a n dt h en u m e r i c a ls o l u t i o n so ft h ed i f f e r e n t i a l - a l g e b r a i ce q u a t i o n s h a v eb e e nt h es t u d ys p o t l i g h ta n dd i f f i c u l tp r o b l e mi nr e c e n ty e a r s i nt h i sp a p e r , w em a i n l yd e s c r i b e da n ds t u d i e dt h en u m e r i c a ls o l u t i o n sm e t h o d s f i r s t l y , 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 ys y s t e m sa r e d i s c u s s e db a s e do nt h et h ed y n a m i c so ff l e x i b l em u l f i b o d ys y s t e m st h e o r y , t h e a d v a n t a g e sa n dd 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 s e v e r a l c o n d e n s e dm e t h o d sa r ei n t u r d u c e d ，w h i c ha r eb a s e do nt h eg e n e r a l i z e dc o o r d i n a t e i n d e p e n d e n td e c o m p o s i t i o ni nd e t a i l g i v e nt h eu s e so fa l l t h em a t r i xf a c t o r i z a t i o n m e t h o d si nd i f f e r e n t i a l - a l g e b r a i ce q u a t i o n sa n dt h en u l ls p a c em e t h o d s e c o n d l y , c o m b i n et h ea d v a n t a g eo ft h e c o n d e n s e dm e t h o dw i t hc o n s t r a i n tv i o l a t i o n c o r r e c t i o n ，w er a i s e dan e wm e t h o d ，c a l l e d an u l ls p a c em e t h o do fc o n s t r a i n tv i o l a t i o n c o r r e c t i o n t h i sm e t h o du s eas i m p l em e t h o dt og e tt h en u l ls p a c e ，s e p a r a t e di n d e p e n d e n t c o o r d i n a t e ，t r a n s f o r mt h ed i f f e r e n t i a l - a l g e b r a i ce q u a t i o n si n t oo r d i n a r y d i f f e r e n t i a l e q u a t i o n s ，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 sw e r ee v a l u a t e db yn e w m a r k - 1 3 m e t h o d ，w h e ng e tt h ev a l u e s ，t a k et h e mt ot h ec o n s t r a i n tv i o l a t i o nc o r r e c t i o nf o r m u l a ，j u d g e t h ed i s p l a c e m e n t sa r ec o n s t r a i n tv i o l a t i o no rn o t ，m o d i f yt h ed i s p l a c e m e n t s ，s ot h a ti tc a n m e e tt h ed i s p l a c e m e n t sc o n t r o le q u a t i o n so fm u l t i b o d ys y s t e m s f i n a l l y , t h es t a b i l i t y a n de f f i c i e n c yo ft h ep r e s e n t e dm e t h o da r ei l l u s t r a t e db yt h es i m u l a t i o nr e s u l t so fo n e t y p i c a le x a m p l e 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 n u l ls p a c em e t h o dg e n e r a l i z e dc o o r d i n a t e sd e c o m p o s i t i o nn u l ls p a c em e t h o d o fc o n s t r a i n tv i o l a t i o nc o r r e c t i o n 西安电子科技大学 学位论文独创性( 或创新性) 声明 秉承学校严谨的学风和优良的科学道德，本人声明所呈交的论文是我个人在 导师指导下进行的研究工作及取得的研究成果。尽我所知，除了文中特别加以标 注和致谢中所罗列的内容以外，论文中不包含其他人已经发表或撰写过的研究成 果；也不包含为获得西安电子科技大学或其它教育机构的学位或证书而使用过的 材料。与我一同工作的同志对本研究所做的任何贡献均已在论文中做了明确的说 明并表示了谢意。 申请学位论文与资料若有不实之处，本人承担一切的法律责任。 本人签名：瘟空鍪只期翌鲨 西安电子科技大学 关于论文使用授权的说明 本人完全了解西安电子科技大学有关保留和使用学位论文的规定，即：研究 生在校攻读学位期间论文工作的知识产权单位属西安电子科技大学。学校有权保 留送交论文的复印件，允许查阅和借阅论文；学校可以公布论文的全部或部分内 容，可以允许采用影印、缩印或其它复制手段保存论文。同时本人保证，毕业后 结合学位论文研究课题再撰写的文章一律署名单位为西安电子科技大学。 ( 保密的论文在解密后遵守此规定) 本学位论文属于保密，在_ 年解密后适用本授权书。 本人签名： 导师签名： 日期塑! z ：至：! 兰 日期丝丝：呈丝 垒一一 强； 平一仔一 襞 第一章绪论 第一章绪论弟一早 三百t 匕 柔性多体系统动力学研究由刚体和柔性体组成的复杂机械系统在经历大范围 空间运动时的动力学行为，是多刚体系统动力学的自然延伸和发展。它主要研究 柔性体的变形与其大范围空间运动之间的相互作用或相互耦合，以及这种耦合所 导致的动力学效应。柔性体的变形运动与柔性体大范围空间运动的同时出现及其 相互耦合是柔性多体系统动力学的本质特征，这个特征使其动力学模型不仅区别 于多刚体系统动力学，也区别于结构动力学，是两者的结合与推广。柔性多体系 统动力学是与经典动力学、结构动力学、控制理论及计算机技术紧密相联的一门 新兴交义学科，在航空航天、机器人、高速机构及车辆等各个领域有着广泛的应 用，成为目前理论和应用力学最活跃的分支之一。 柔性多体系统的动力学数学模型分为两类，类数学模型可简化为纯微分方 程，另一类则是微分一代数混合方程。近二十年来国内外进行了大量的研究工作， 纯微分方程的数值解法已经基本成熟，而微分一代数方程的数值算法是近年来研究 的热点问题。微分一代数方程的数值算法目前的研究方法大体可分为两类：一种是 从微分一代数方程组本身出发，利用现代数学的研究成果将约束定义为流形，对微 分一代数方程组进行降阶处理，将其转化为由约束方程定义的流形上的常微分方 程。这种方法的优点是可以直接应用求解常微分方程的技术，避免约束方程违约。 但在求解过程中必须计算由约束方程定义的流形零空间的基，计算工作量大，对复 杂的多体系统，零空间基的计算缺乏成熟的方法，且有时并不唯一；另种方法是 在动力学方程中引入附加校正项，当约束方程产生违约时，对动力学方程进行校 正。目前的校正方法多为间接校正方法，不能对系统的广义坐标进行直接的校正 以满足约束方程。另外，在动力学方程中加入附加校正项需给定校正系数，校正 系数太小校正效果不明显，校正系数太大容易引起动力学方程的破坏。目前还没 有校正系数的自动选取方法，大都凭经验选取校正系数。 1 2 多柔体系统动力学研究概况 对多柔体系统的研究有着实际的工程应用背景。随着动力学模拟的深入，人 们发现，系统中某些物体的变形有时会对系统性能产生非常重要的影响。在航天 器、机器人领域和机构设计等方面，部件有向轻质量和高速度发展的趋势，其中 2 求解多柔体系统动力方程的违约修正零空间法 系统中轻质量大尺度部件的高速运动往往引起系统的剧烈振动，达不到高精度要 求，甚至毁坏系统的某些部件。为了解决这个问题，就应该考虑系统中某些部件 的弹性变形，即在抽象物理模型时就要考虑某些部件的柔性效应。同时，必须考 虑柔性体的变形与其大范围空间运动之间的相互作用或相互耦合，以及这种耦合 所导致的动力学效应的研究，这类系统称为柔性多体系统或n - 柔混合多体系统 ( s y s t e mc o n s i s t i n go fr i g i da n df l e x i b l eb o d i e s ) t 1 1 。 对于柔性多体系统来说，它的各构件之间一般都存在着大的相对平动和转动， 而且在运动中要考虑构件的柔性。这就使得系统运动的自由度，各构件互相之间 的运动学关系都大大地复杂化了。同时引起了复杂而变化的离心力和哥氏力力场， 影响了多体之间相互运动的力学条件。这就使得考虑多体运动和柔性效应之间的 耦合显得极其重要。文献 1 认为：柔性多体系统不同于多刚体系统，它含有柔性 的部件，变形不可忽略，其逆运动学是不确定的；它与结构动力学不同，部件在 自身变形运动同时，在空间中经历着大的刚性移动和转动，刚性运动与变形运动 互相影响、强烈耦合，与一般的系统不同，它是一个多输入、多输出的无穷维、 时变、高度耦合、高度非线性的复杂系统。总之，多柔体系统动力学是与经典动 力学、结构动力学、连续介质力学、计算力学、现代控制理论及计算机技术紧密 相联的一门新兴交叉、边缘性学科，在航空航天、机器人、高速机构及车辆等各 个领域有着广泛的应用，成为目前理论和应用力学最活跃的分支之一。它的主要 任务是研究建立系统的适合计算机的动力学模型的建模方法和有效的数值求解方 法。 随着部件尺寸的增大、结构重量的减轻，从而刚度的减弱以及运行速度的提 高，在许多方面都提出了多柔体系统建模的需求。在人造卫星、航天飞机、大型 空间站等的动力学分析中，由于它们的天线和太阳能帆板的伸展尺寸与本体尺寸 相比，可能大到几倍甚至十几倍，此时弹性变形不再可以忽略。1 9 5 8 年美国发射 的第一颗人造卫星“探险者1 号”( e x p l o r e r i ) 。由于在系统的动力学建模时没 有计及4 根鞭状天线的弹性影响，导致卫星入轨后翻滚、失控；1 9 8 2 年美国“陆地 卫星一4 ( l a n d s t a i v ) 的观测仪的旋转部分受到柔性太阳帆板驱动系统的干扰 而产生微小扰动，从而降低了图像质量；“国际通讯卫星v 号”( i n t e r s a t v ) 柔 性太阳帆板扭振频率与驱动系统发生谐振时，导致帆板停转和打滑。各种问题的 提出引发了人们对多柔体系统建模的思索和重视，开始了多柔体建模的探讨和研 究。 七十年代初期开始，p wl i k i n s ，w jb o o k ，jps a d l e r 等人对柔性系统进行了大 量的研究工作。随着有限元方法的出现，发展和成熟，1 9 7 1 年r cw i n f r e y t 2 】和 e r d m a n 3 】先后在不考虑构件弹性变形对其大范围刚体运动影响的情况下，把结构 分析中的有限元法引入弹性系统的分析中，从而为展开弹性多体系统动力学开辟 第一章绪论 了一条新路。w i n f r e y 等人的工作，标志着机构弹性动力学( k i n e t o e l a s t i od y n a m i c a n a l y s i s ，简称k e d ) 近期研究工作的开始。该方法的要点是，不考虑构件的弹性变 形对其大范围运动的影响，而是通过对多刚体系统动力学分析得到构件运动性态， 再加上构件的惯性特性，以惯性力的形式加到构件上，然后根据惯性力和系统的 外力对构件进行弹性变形以及强度分析。这种方法实质上是将柔性多体系统动力 学问题转变成多刚体系统动力学与结构动力学的简单叠加，忽略了二者之间的耦 合。 到了八十年代，柔性多体系统动力学进入了高速发展阶段。k e d 更加精确化， 形成了比较精确的数学模型。在k e d 方法日趋成熟之时，柔性多体动力学，俗称 f m d ( f l e x i b l em u l t i b o d yd y n a m i c s ) 逐渐成为研究的焦点。随着轻质、高速的现代机 械系统的不断出现，k e d 方法的局限性日益暴露出来。为了计及构件弹性变形对 其大范围运动的影响，人们首先对柔性构件建立一个浮动坐标系，将构件的位形 认为是浮动坐标系的大范围运动与相对于该坐标系的变形的叠加，提出了用大范 围浮动坐标系的刚体坐标与柔性体的节点坐标( 或模态坐标) 建立动力学模型的方 法。在具体建模的过程中先将构件的浮动坐标系固化，弹性变形按照某种理想边 界条件下的结构动力学有限元( 或模态) 进行离散，然后仿照多刚体系统动力学的方 法建立离散系统数学模型。这种方法虽然考虑了构件弹性变形对大范围运动的影 响，但在对柔性体离散时并没有考虑大范围运动对其的影响，且在对有限元( 或模 态) 进行离散时具有很大的随意性。实质上这种方法是柔性多体系统的一种零次近 似的耦合动力学。尽管如此，国内外的学者对这种模型的研究持续约十来年，在 建模方法研究的基础上重点解决了数学模型数值解的病态问题，并在工程领域得 到一些应用。近年来的研究表明，采用这种零次近似的耦合方法得到的柔性多体 系统动力学的分析结果，有的和工程实际比较接近，精度较高。b o l a n d ，k a n e ， ws a n a d a ，sd u b o n s k y ，aas h a b a n a ，ejh a n g 等都为柔性多体系统动力学建模 做出了杰出贡献1 4 】。 中国学者对于多体系统动力学的研究始于8 0 年代【5 】，对太阳能电池帆板阵、挠 性机械臂、星载可展开天线等的动力学分析做了许多研究工作。文献 6 对柔性多 体系统动力学研究进行了较为全面的回顾。到目前为止，柔性多体系统动力学的 研究虽然取得了一些成果，在实际工程中得到了应用，但是远没有达到多刚体系 统动力学的研究水平，其主要原因是对物体大范围运动和弹性变形耦合问题的认 识和处理方法上遇到困难。 1 3多体系统动力学数值计算方法研究概况 柔性多体系统的动力学方程是强耦合、强非线性方程，这种方程的求解目前 4 求解多柔体系统动力方程的违约修正零空间法 只能通过计算机用数值方法解决。数值求解方法是多体系统动力学研究的一个重 要方面。慢变的、大的刚性运动和快变的、小的弹性运动的耦合使得柔性多体系 统的动力学方程成为时变的非线性的刚性方程，其求解过程易出现病态或不收敛， 故多柔体系统的动力学方程求解引起学者们的极大关注，形成了各种不同的求解 法。多柔体系统动力学各种方法的数学模型都可以归纳为刚性微分方程组和微分 代数混合方程组两种类型。 刚性微分方程的数值积分方法的研究【7 】已经很成熟，目前对于解这类方程的数 值算法主要有：牛顿一拉斐逊( n e w t o n r a p h s o n ) 法；直接积分法( 中心差分法、威尔 逊口法、纽马克法、帕克( p a r k ) 冈1 性稳定法，龙格库塔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 方法) ，用差商代替偏导数，提高该矩阵的计算效率贝, l j n e w t o n 迭代法将是相当理 想的方法。为提高数值计算的稳定性，克服刚性方程在数值计算中的困难。隐式 数值计算方法逐步被采用【8 】。常用的隐式算法有r u n g e k u t a 法、线性多步法、g e a r 方法等。 根据拉格朗同方程，使用非独立广义坐标的方法建立的多体系统动力学方程 是一组非线性的微分代数混合方程组，又称为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 ) 通常要通过数值方法进行求解。多体系统动力学微分代数方程的解法从 方程的结构形式看可分为两大类【9 】。一类是拉格朗同乘子法，即通过引入拉格朗同 乘子释放系统的约束，因此，控制方程中既有独立的广义坐标又有不独立的广义 坐标，同时又包含特定乘子，系统的运动微分方程与系统的约束方程构成了系统 的动力学控制方程，这是一组微分代数混合方程组。当用欧拉参数描述动参考系 在空间的方位时，在约束方程中，还包含相应的数学约束，这种形式的方程人们 通常称为最大未知量法( 增广法) 。其主要特点是：方程数目多，是微分一代数混 合方程组，存在着求解的数值困难。不过这时方程组的系数矩阵通常呈现稀疏状， 易于应用稀疏矩阵的特点进行快速数值求解，提高数值算法效率。在约束方程以 约束库的形式存入计算机里的情况下，这种形式便于对复杂系统自动建模，并可 根据需要同时求出任何约束的约束反力。但它要求位置和速度的初值严格满足约 束条件，否则容易因误差积累而导致解的发散。另一类是全部用独立的广义坐标 描述的纯微分方程组，通常称为最小未知量法( 缩并法) 。其主要特点是：参与 求解的方程数目最少，是纯微分方程组，但方程呈现强耦合，系数矩阵通常为满 阵且矩阵中各元素都很复杂，其复杂程度还与独立坐标的选取有关，在建模时需 第一章绪论 要人为干预，或需要在程序中作一定的预处理工作，分离并消去不独立的广义坐 标，例如0 r 法，奇异值分解( s v d ) 法等，因而导致复杂的递推公式，使占用o p u 的时间过长。此外，对于那些需要求出约束反力的系统来说，这种形式反而不太 理想。 数值算法方面研究的主要问题包括常微分方程的刚性问题、微分一代数方程的 数值解法问题、多体系统非线性动力学行为的数值分析方法等。微分一代数方程组 的求解方法已成为目前多体系统动力学的难点问题。近二十年国内外进行了大量 的研究工作。目前的研究方法大体上可以总结为以下三种【l o 】： ( 1 ) 通过数值方法，利用计算机自动寻找独立变量个数并选择，将方程缩并 为个数与自由度数相等的方程后再进行数值处理。将微分代数方程转化为常微分 方程o d e ( o r d i n a r yd i f e r e n t i a le q u a t i o n s ) ，( 方程的个数等于系统独立的广义坐标 数) ，该方法称为缩并法( c o n d e n s e dm e t h o do re l i m i n a t e dm e t h o d ) 。 ( 2 ) b a u m g a r t e 1 1 】稳定化方法，将全部广义坐标与拉格朗日乘子作为未知量 统一处理，将方程变为较大变量数的方程，由微分代数方程得到常微分方程( 方程 的个数与系统的广义坐标数相等) ，该方法称为增广法( a u g m e n t e dm e t h o d ) 。 ( 3 ) 利用多体系统动力学方程的特点，直接构造微分代数方程的数值差分格 式。文献【1 2 】作了较全面的综述，在这方面，jwb a u m g a r t e ，raw e h a g e ，tw p a r k ， ejh a u g 等人作出了重要贡献，相应出现的方法有直接积分法( d i r e c ti n t e g r a t i o n ) ， 约束稳定法( c o n s t r a i n ts t a b i l i z a t i o n ) ，广义坐标分离法( g e n e r a l i z e dc o o r d i n a t e p a r t i t i o n i n g ) ，混合数值积分法( h y b r i dn u m e r i c a li n t e g r a t i o n ) 等。所以微分代数方 程的处理是多体系统动力学的重要内容。 1 4 数值算法中的病态问题 柔性多体系统是多自由度的强耦合非线性系统，建立柔性多体系统动力学方 程时，由于系统中存在量级异常差别的慢变大位移与快变弹性变形之间的相互耦 合，使柔性多体系统动力学方程表现出微分方程组中变化快的解分量很快地趋于 它的稳定值，而变化慢的分量则缓慢地趋于它的稳定值，这种性质给这类力学模 型的数值求解带来了s t i f r 特性。在数学上称之为刚性方程，也称为病态方程或坏 条件方程【1 3 】。例如：在柔性多体系统动力学中含有慢变的动参考系运动与快变的 柔性体的振动；在控制系统中，一般来说控制部件反应灵敏，是快变的，具有小 的时间常数，而受控物体一般惯性大，是慢变的，具有大的时间常数，如航天飞 行器的控制问题。 对此种刚性方程求解方法的研究，是近二十年来的热点之一，本文前面提到 的吉尔等人做了许多重要工作，但还有许多问题尚待解决，非线性问题数值稳定 6求解多柔体系统动力方程的违约修止零空间法 性的研究刚开始，在解决各种实际问题时，需要更适用的算法。 造成动力学方程组刚性的原因主要有：( 1 ) 柔性构件大范围运动与弹性变形 间的强耦合；( 2 ) 方程系数的强非线性；( 3 ) 积分变量过多。这些因素表现在 动力学方程的数值积分过程中产生快变量与慢变量，从而造成数值积分受到限制， 积分速度过慢。然而无论哪种形式的柔性多体系统的动力学方程，都存在着同时 含有描述刚性运动的慢变量与描述变形运动的快变量的问题，因此有关数值病态 问题需要特别关注。 1 5 研究数值算法的必要性 柔性多体系统动力学控制方程，由于其强耦合和高度非线性，通常得不到响 应的解析解，必须借助于数值方法。微分一代数方程组的数值求解方法是近年来的 研究热点，也是多柔体系统动力学的研究难点，它的研究进展影响着多柔体系统 动力学的发展。 目前已有许多用于求解微分一代数方程的数值方法，但是这些方法都有其优缺 点及适用范围。如广义坐标分离法【1 4 】，坐标的选取没有一定的准则，需要人为干 预；b a u m g a r t e 约束稳定法，缺点是必须依靠经验选取校正系数，不能将违约有效 地控制在给定的精度范围内；洪嘉振的约束违约直接修正法，缺点是修正解不满 足动力学控制方程。 因此，寻求一个稳定、可靠、精度高、通用性好的数值方法是多体系统动力 学研究的重要任务，数值算法的研究具有很大的必要性，本文的重点工作就是研 究数值算法。 1 6 本文的主要工作及章节安排 1 6 1本文的主要工作 多柔体系统动力学是近年来的研究热点，动力学方程的数值解法是其研究的 主要内容，对多柔体动力学方程准确、迅速的求解是多体系统动力学的主要研究 任务。目前的数值求解的方法很多，但是这些研究方法都有其局限性与适用范围。 本文在研究了各种数值算法的基础之上提出了一种违约修正零空间算法，首先基 于广义坐标的分解，用一种简单方法求得了约束方程雅可比矩阵的零空间，然后 把微分代数方程化为关于独立坐标的纯微分方程，然后用纽马克法进行数值求 解，最后用一种新的违约修正迭代格式对数值解进行修正。本方法的提出既避免 了以前计算零空间的复杂的广义逆的求解问题，又避免因坐标约束方程的违约而 第一章绪论 7 引起的数值解的发散问题，可以得到比较理想的效果。基于本文算法的提出，我 们所作的主要工作如下： ( 1 )对多体系统动力学的数值算法进行了深入研究，着重研究分析了各种 微分一代数方程数值算法的优缺点及适用范围，其中又重点研究了违约稳定法与违 约修正法。 ， ( 2 )深入研究了缩并法的基本思想，然后研究了基于广义坐标的独立性分 解的几种常用微分一代数方程的数值算法，包括l u 分解法、高斯消去法、q r 分解法、 s v d 法，并总结了其优缺点及其适用范围。 ( 3 )深入研究了零空间法【”】在多体系统动力学中的应用情况，从中发现其 优缺点及适用范围，并给出了一种简便的零空间的求解方法。 ( 4 )在缩并法的基础上，提出了一种零空间法与违约修正思想相结合的新 的数值算法一违约修正零空间法。首先确定约束方程的雅可比矩阵的零空间，分 离出独立的广义坐标，用独立的广义坐标把微分一代数方程化为了纯微分方程，然 后用纽马克法对纯微分方程进行迭代求解，用一种新的违约修j 下迭代格式对迭 代初值进行违约修正，从而使所得的数值解既能满足控制方程又能满足约束方程。 ( 5 )通过一个旋转杆一滑块系统，建立了其动力学控制方程和约束方程，用 本文提出的违约修正零空间法对其进行仿真计算，与其它数值算法的仿真结果进 行对比分析，说明了本文提出的算法的正确性与有效性。 1 6 2 本文的章节安排 第一章绪论。介绍了多柔体系统动力学研究概况，数值计算方法的研究概 况以及研究算法的必要性，并对数值算法中的病态问题进行了讨论。 第二章多柔体系统动力学方程的数值方法研究。介绍了各种数值求解方法， 指出其优缺点及适用范围。重点论述了微分代数方程的数值求解方 法，最后详细给出了缩并法的思想。 第三章基于广义坐标分解的数值算法研究。详细介绍了各种基于广义坐标分 解的缩并算法，最后介绍了零空间算法的思想。 第四章违约修正零空间法。在零空间法的基础之上，结合违约修正的思想， 提出了一种新的违约修正零空间法，并给出了具体计算步骤和算法 流程图。 第五章动力学算例仿真与分析。对一个典型动力学模型进行仿真分析，运用 本文提出的新方法对其进行仿真计算，并与其它方法的仿真结果进 行了分析比较。 第六章总结与展望。总结全文，指出下步的研究方向与内容。 第二章多柔体系统动力学方程的数值方法研究 9 第二章多柔体系统动力学方程的数值方法研究 2 1 引言 l a g r a n g e 方法是建立多体系统动力学方程的普遍方法之一，随着多体系统动 力学研究内容的不断深入，研究范围的不断扩展，l a g r a n g e 方法也得到了相应的 发展。传统的l a g r a n g e 方程可分为第一类和第二类，前者适用于广义坐标独立的 质点系，其动力学方程为常微分方程组。后者适用于广义坐标不完全独立的质点 系，其动力学方程为微分代数方程组。 对于具有完整、定常约束的树形结构多体系统，当广义坐标取铰链坐标和变 形体的模态坐标时，系统的广义坐标是独立的，应用第一类l a g r a n g e 方程，该系 统的动力学方程为 m q + 啪+ 娑：q ( 2 1 ) o q 上述表达式是一组二阶常微分方程组，未知量是广义坐标，其个数与方程组的个 数相等，都等于系统的广义自由度数。这是数目最小的封闭微分方程组，其求解 方法可直接沿用一阶常微分方程的数值解法。 根据l a g r a n g e 公式，弓l h l a g r a n g e 乘子，可得第二类l a g r a n g e 方程： 绚+ 幻+ c ：a = 绯+ 9 ( 2 - 2 ) c ( q ，) = 0 ( 2 3 ) 其中m ( q ，t ) 尺删为系统广义质量矩阵，q r ”为系统广义坐标矢量，a r ”为 l a g r a n g e 乘子矢量，c ( q ，t ) r ”为系统约束方程的左部，c 。( g ，t ) = o c ( q ，t ) 幻为 j a c o b i 矩阵，k e 尺是刚度矩阵，q ，广义主动力绋为与速度二次项有关的广义 力。 这是一个典型的微分代数混合方程组，它与常微分方程不同，存在着数值求 解上的困难，这也是本文的数值方法的具体研究对象。本章将主要介绍多柔体系 统动力学方程的各种数值方法，指出其优缺点及适用范围。 2 2 多柔体系统动力学方程的数值方法介绍 通过近二十年国、内外学者的研究和探索，出现了许多求解多体动力学方程 的数值方法，本节将介绍各种数值求解方法，并对其中几个典型的、适用性广的 i o 求解多柔体系统动力方程的违约修止零空间法 数值方法进行详细地论述。 2 2 1二阶微分方程的直接积分法 二阶微分方程直接数值积分就是不对方程作任何变换，将求任意时刻都满足 此式的的广义位移向量g ( f ) ，代以求仅在离散时刻点( 江0 ，1 ，) 满足该方程的 吼。一般来说，吼与g ( ) 总有一些差异，g ，是9 ( ) 的近似解。其误差来源于积分 的每一步的截断误差和舍入误差以及每步所产生的误差在以后各步中的传播。截 断误差与积分步长有关，舍入误差则影响算法本身的稳定性。积分过程是个逐步 数值积分过程。 1 中心差分法 对于图2 1 所示的位移时间历程曲线 - 7 1 a t g _g i 冒“_ 图2 1 位移时间历程曲线 在时间间隔f 中点的速度可写为 嘭+ 等2 警 t 点的加速度为： 自k 一矗k 磊2 将( 2 4 ) 式代入( 2 5 ) 式，有 ( 2 4 ) ( 2 5 ) 玩= 矿1 ( 一2 9 ，+ ) ( 2 6 ) 于是中心差分格式为： 或= 面1 【q ，+ a t - 2 q , + 吼一出】 玩= 古( 一2 吼+ ) 若设二阶微分方程的具体形式为 m 畦+ c 矗+ k q = q t ( 2 7 ) ( 2 8 ) ( 2 9 ) 第二二章多柔体系统动力学方程的数值方法研究ii 把( 2 - 7 ) 、( 2 8 ) 式代入上式后，有 ( 古膨+ 击c ) = q 一( k 一志肘( 古膨一击c ) 乳& ( 2 - 1 0 ) 由( 2 1 0 ) 式可求得纵a f 。而对应t 时刻的速度或和加速度或将由( 2 7 ) 、( 2 8 ) 式求得。 用中心差分格式求吼+ 出时，由( 2 1 0 ) 式可见，涉及到前两步的值，吼及吼一厶f 。 是种显示积分格式。 注意到，当垃 夕乞时，中心差分格式将是不稳定的。式中。为离散系 统的最高频率。 在计算每一步的g i 时，需要知道g 咄的值，但是一般初始条件仅给出q o 及t l 。的 值，即存在一个如何起动的问题。为此，需要由( 2 7 ) 式解出q o + = q 垃，代入( 2 - 8 ) 式，整理得： g ：吼一幻j + 兰譬魏 ( 2 1 1 ) 式中魏在实际计算中可直接给出，或通过方程 m o o + 回o + 脚。= o o ( 2 1 2 ) 求得。 2 威尔逊9 法 威尔逊秒法是假设在时间间隔8 a t ，秒1 o 中，加速度为线性变化，如图2 2 所示。且系统的动力特性在此间隔中保持不变。 图2 2 加速及时间历程曲线 设r 为t 与f + 铋f 间隔内的时间增量( o f o a t ) ，因此，在f 到f + 抛f 间隔内假设 或+ f - 玩+ 志( 或恤一 ，( 2 1 3 ) 对上式积分求得f + f 时刻的速度和位移 ，2 “= 或+ 玩付盏( 玩+ 越一玩) ( 2 - 1 4 ) 1 2 求解多柔体系统动力方程的违约修止零空间法 鸹+ 或f + 三矿+ 盖( 玩+ 甜一玩) ( 2 - 1 5 ) 令f = o a t ，代换后解得用g ，+ 触表示的速度和加速度为： 或+ 触= 志( 纵触一q t ) - 去或_ 2 玩 ( 2 - 1 6 ) 吼+ 触2 两砑吼+ 触一 面吼一孢 6 ：瓦3 ( 吼他一g f ) 2 4 , 一跳 。17)-24,-i-q,(2-1吼+ 触。面( 吼他一g f ) 一 现在考虑t + 0 a t 时刻的运动方程，以此解得t + a t 时刻的位移、速度和加速度。由 于加速度是线性变化的所以力向量也可以作线性假设，即 吮+ 锄+ q ，+ 趾+ 物，+ 弛= q + 甜 ( 2 - 1 8 ) 将( 2 1 6 ) 式、( 2 1 7 ) 式代入( 2 1 8 ) 式，得 ( 志m + 面3c + k ) g ，+ 甜= q 抛+ ( 志m + 丢。吼+ ( 嘉m + 2 c 地+ ( 2 肘+ 警c ) 苗 ( 2 1 9 ) 由( 2 - 1 9 ) 式求得q ，+ 他后，再l h - f 式求得t + a t 时刻的速度和加速度及位移值 q t + t u - - 上0 3 ( a t ) 2 ( g f + 蛐一q t ) - 志州l 一扣 ( 2 - 2 0 ) ( g 惝一瓦g r + 【1 一万坦 2 。2 u 舢= 玩+ 等( 或+ f + 玩)(一21)q 22 1 ，+ 血2 吼+ - = ( g f + f + g ，) 【。 g f + 血：吼+ 增，+ 竽( 玩+ f + 2 或) ( 2 2 2 ) 显然，威尔逊0 法是一种隐式差分格式。对于线性系统，当0 1 3 7 时，方法 是无条件稳定的。而对于非线性系统，常取0 = 1 5 。 注意到一个异常情况，在t + a t 时刻用该法计算，系统从不满足平衡。 3 帕克刚性稳定法 帕克刚性稳定法是综合了吉尔两步法和三步法而得到的，它对于低频和所有 高频成分都是稳定的。这个方法的速度和加速度差分格式为 或+ 2 面1 ( 1 0 q t + - - 1 5 q t + 6 吼一出一吼拙) ( 2 - 2 3 ) 或+ & = 志( 1 0 或+ 州或+ 6 玩矿或拙) ( 2 - 2 4 ) 由纽马宽法的列式求得吼+ 址，再将其代入( 2 2 3 ) 、( 2 2 4 ) 式可得f + f 时刻的 速度和加速度。 帕克法在计算吼讪时，涉及到t 、卜f 、f 一2 a t 时刻的位移和速度，因而需要 第二章多柔体系统动力学方程的数值方法研究 1 3 有专门的起动过程。 本方法的缺点在于需要大量的计算机内存空间用以存储前三个时间步上的位 移和速度。 2 2 2 一阶微分方程的数值积分法 由于高阶微分方程总可以化为一阶微分方程组求解，因此多体系统的二阶微 分方程组也可以变为相应的一阶微分方程组，这不仅是因为一阶微分方程的数值 积分已经有许多成熟的算法可供应用，而且，由于将柔性多体系统的动力学控制 方程转变到状态空间，这样更便于与现代控制方法的仿真兼容。 一阶微分方程的初值问题可表述为给定微分方程： 雪= f ( q ，f ) ( 2 2 5 ) 和初始条件 q ( t o ) = 吼 ( 2 2 6 ) 若f ( q ，t ) 在区域gcr 肿1 上连续并存在连续偏导数，则存在唯一解。 对( 2 2 5 ) 、( 2 2 6 ) 两式的数值积分，就是在已经求得q 的基础上，从推到f 川 的数值解q 川，即 q ( t n + i ) 一g ( 乙) = r ( g ( f ) ，t ) d t ( 2 2 7 ) 当q ( t n ) 为已知时 q ( t n + 1 ) = g ( 乙) + r 1 厂( g ( f ) ，f ) 仍 ( 2 - 2 8 ) 由于任何连续函数可以用足够高次的多项式在任意区间内逼近，因而，在 k ，乙+ 。】中的f ( q ，f ) 可以通过对过去那些点上的z 的插值多项式来逼近，以求得近 似解。因而微分代数方程的解归结为( 2 2 8 ) 式右端的积分问题，用不同的插值多 项式和它们的不同阶次逼近厂( g ( f ) ，t ) ，于是出现了不同的数值积分方法的列式。 1 龙格库塔法 龙格库塔法应用非常广泛，是一种求解常微分方程的方法，最常用的是4 阶4 级定步长法，该方法不求函数的导数。 四阶龙格库塔法的列式为 g 州2 吼+ 姆 ( 2 2 9 ) 式中 1 g = ( 石+ 2 五+ 2 石+ 工) ( 2 - 3 0 ) 1 4 求解多柔体系统动力方程的违约修止零空间法 其中 z = f ( q 。，t 。) 石= 几。+ 宝朋。+ 互h ) 石川参。，+ 争 厶= 厂( g 。+ 坛，f 。+ 冬) 另外一种列式为 g 肘l = q 。+ 鲁( z + 3 l + 3 六+ 丘) ( 2 - 3 1 ) 其中 z = f ( q 。，t 。) 以川”弘”争 工川了hz + h l + 警) 六= f ( q 。+ m h l + 鹿厶，t 。+ h ) 龙格库塔法的截断误差为r = o ( h 5 ) ，其主要缺点是在每个时间步长上函数 f ( q 。，t 。) 都要求解四次，而这些函数值在以后的计算步中不再用到，因而效率较低。 2 刚性微分方程的吉尔解法 在可以用常微分方程来描述的许多实际物理或者化学过程中，往往包含了许 多复杂的子过程以及它们之间的相互作用，其中有的子过程表现为快变化的，另 一些相对来说是慢变化的，并且变化的速度可以相差非常大的量级，描述这些过 程的常微分方程的解中也将包含快变分量和慢变分量。如果在一个过程中的快变 子过程与慢变子过程的变化速度相差非常大，在数学上称这种过程具有刚性，而 描述这种过程的常微分方程则称为刚性方程，也称为病态方程或坏条件方程。 从数值观点看，为获得响应的快变分量，此时应该用小的时间步长积分，因 为一方面在相当长的积分区间内采用如此小的步长的积分过程，往往是不现实的， 另一方面当快变分量已趋稳定或消失时，就应该用较大的时间步长积分。但是理 论与实践说明，很多方法，特别是显式方法的步长仍不能放大，否则便出现数值 不稳定现象和误差急骤增加，以致掩盖了真解。正是由于这种性质，使得传统的 常微分方程的数值积分遇到了极大的困难。 1 9 6 8 年吉尔指出向后差分的数值积分公式在无穷远处具有良好的稳定性质， 并且用下面的公式编制了一个求解刚性微分方程组的通用程序包。 第二章多柔体系统动力学方程的数值方法研究 1 5 g 枞= 口j q 州+ j l l a 以“ ( 2 3 2 ) i = 0 吉尔法是求解刚性方程的许多方法中最好的几种方法之一，至今仍是最有效 的通用算法，它主要具有下面的三个优点：l 容易改变阶次和步长；2 。能够应 用高阶的和高稳定格式；3 。每前进一个步长，解隐式方程组所需要的工作量比 较小。它的缺点是方程阶数增大，初始值难选定，多步时间离散对约束方程影响 较大，导致解的漂移。 2 2 3 其它数值方法 1 精细时程积分法【1 6 】 利用指数函数的加法定理，计算数矩阵丁的精度，从而确定解的精度，是一种 基于指数矩阵的精细计算格式的2 算法。这个算法的优点是显示积分；无条件稳 定，时间步长由精度控制；计算精度高，效率高，适用面广。缺点是非齐次问题 需要进行矩阵求逆和大量满秩矩阵相乘；需大量的存储空间，计算工作量大，精 细积分法适合自由度数较少的简单结构。 2 约束稳定罚函数法 把原动力学方程变为两个独