[转帖]多体动力学建模及求解过程中的几个关键问题

AaronSpark

　　内容简介：多体动力学应用系列文章题纲 邱向军 xiangjunqiu@yahoo.com 1。多体动力学建模及求解过程中的几个关键问题 2。多体数值求解方法中的“patch test“ 3。矿山机械中的多体问题 4。多体随机碰撞数据处理 5。摩擦磨损与多体方法 6。多体系统中的陀螺系统处理 7。哈密尔顿体系与多体系统 8。将流体离散为”多体“的数值方法 9，流固耦合及水甸塘模拟 10。泥沙淤积问题中的多体方法 11。多体动力学的应用随想 多体动力学建模及求解过程中的几个关键问题

　　邱向军

　　xiangjunqiu@yahoo.com

　　我于1990年起在一家位于美国的矿山设备设计咨询公司从事研发工作，从1997年起从事多体动力学在矿山机械中的应用研究。由于应用领域的特殊性，一般的商业通用软件不能满足要求。于是老板同意让我与一位俄裔科学家合作开发了一套供矿山机械应用的多体动力学软件。历时四年，我们的软件终被市场认可。最后，于2002年一家矿山机械制造公司买断了该软件，同时雇我与俄裔科学家继续从事多体动力学的研究与应用。在多年的开发应用与研究中，我有机会接触了许多多体应用课题。这里，谨谈一下对多体动力学建模及求解过程中的几个关键问题的体会。供同行参考，并想借此机会听取同行意见。

　　1. 多体随机碰撞及基本方程体系

　　当一系统的运动物体超过一定数量(>1000)，且物体之间有随机碰撞，运动控制方程的形式必须具有一般性以容纳随机变化的约束条件。习惯有限元的人都喜欢控制方程只含位矢变量的形式。但该形式不便处理随机变化的约束条件。本人偏爱一种称为Discrete Element Method (DEM)所用的控制方程。在该控制方程中，同时呈现了物体的位矢变量及相互约束力。DEM 用类似有限元罚函数方法来处理碰撞约束力并可以方便地处理摩擦力。早期的DEM只能处理多个刚性圆球物体的随机碰撞运动，用途狭窄。近10年来由DEM派生出许多新方法可以处理柔性任意几和形状的多体运动。由于碰撞约束的随机性，每隔几步时间积分必须预测每个物体下几步可能相撞的对象。

　　2. 柔性多体系统的有效算法

　　若多体中某一柔性体变形属大应变，则必须用建立在格林应变与科希霍夫应力基础上的本构关系来处理该柔性体运动。若柔性体变形属小应变大变形或小应变小变形，用处理大应变的方法就非常不经济了。此时，比较有效的方法是将该柔性体分成几个子结构，每个子结构配备一个运动坐标系。在运动坐标系内，相应子结构的变性属小变形。运动坐标系随相应子结构作刚体运动，而子结构的变形力必须参考其运动坐标系来计算。确定运动坐标系的运动有多种方法，这里不作讨论。熟悉有限元的人都知道“patch test”。 其实对 柔性多体系统，也有一套类似的“patch test”来检验柔性多体系统算法的”保刚性“。限于篇幅这里不做讨论。

　　3. 有效多体并行算法

　　随着视窗操作系统的发展及互联网的出现，人们对SERVER计算机要求越来越高。其中多CPU共享内存的SERVER计算机应运而生。人们已意识到这种类型的计算机将对计算力学带来深远影响：新型并行算法。人们对并行算法的研究已至少有半个世纪，但都是针对多CPU分布式(非共享)内存的计算机而进行的。相对于旧并行算法，新并行算法更有效，编程更容易。只要学一下C++的Multi-thread技巧,余下的将无师自通。对多体计算来说， 新并行算法考虑如何将多体分配到每个CPU时，只需按物体顺序号码来进行。而旧并行算法必须根据物体的空间位置来进行。本人在多体新并行算法上已有4年多经验，但苦于无人交流。(偶尔与熟悉多体旧并行算法的人交流一下。) 在此呼吁更多的同行加入新并行算法研究的行列。其实多CPU共享内存的SERVER计算机并不贵。在美国2000美元就可以买到两CPU(AMD)共享内存的低档SERVER计算机。

　　最后，谈谈对多体通用程序的看法。由于多体系统的复杂性多样性远远超过有限元所能处理的结构系统，商业多体程序不可能在有效性与通用性上达到一致。通常是牺牲有效性来满足通用性。对于特殊复杂多体系统，例如矿石球磨机中含百万个物体的多体系统且有随机碰撞，商业多体通用程序根本就无能为力。这时，针对某一类型的多体系统，需要专门研究设计专门的方法自己编程求解。这也是做多体分析的工作者的饭碗所在。以上只是一管之见，与同行交流。

　　多体数值求解方法中的“patch test“

　　邱向军

　　xiangjunqiu@yahoo.com

　　熟悉有限元的人都知道“patch test”。它要求

　　1) 当物体发生刚体位移时单元不产生内力，

　　2) 当物体朝一个方向发生均匀变形时，有限元应得到与理论相符的常应力。

　　由于当有限元划分趋密，每个单元的应力将趋于常数，因此通过patch test是有限元收敛的充分条件。

　　毫无疑问，这种patch test也一定适用于柔性多体方法。只是相应的patch test更复杂。对每个柔性体的patch tes要求如下

　　1) 当物体发生刚体平移时单元不产生内力，

　　2)当物体延一个方向发生均匀变形但处于静平衡时，柔性体方法应得到与理论相符的常应力。

　　3)当初始运动条件(包括初速度及初内应力分布)将使物体处于转动动平衡状态时(例如绕固定轴作匀速转动)，柔性体方法应得到与理论相符的动平衡运动解。

　　比较上述两组patch test条件，可以发现柔性多体方法比有限元方法的patch test多了第三条件。增加这第三条是希望当每个柔性体自由度数趋于无穷大而柔性多体系统趋于动平衡时，柔性多体方法的解应趋于理论解。必须说明的是，这三个条件并不是柔性多体方法对于一般运动求得的数值解趋于理论解的充分条件，而只是对于动平衡运动的数值解趋于理论解的充分条件。

　　由于大多数柔性多体方法采用有限元作为基本元素，因此一般都能满足前两个条件。然而这第三个条件对很对新方法构造者来说是一个挑战。目前绝大多数柔性多体方法都采用显式时间积分方法来求数值解，因此只要采

　　用了前一篇文章提到的运动参考坐标系，就可满足第三个条件。我的一位同事曾经给我讲了他在作博士后时的一段经历。为了课题的需要，他必须自己编制一个柔性多体数值求解程序。他试用了很多常规有限元都无法通过这第三条件试验。原因是他不知道应该同时建立一个运动参考系。(由于该同事是做有限差分法研究的，对有限元了解有限。因此走了弯路。) 最后他用了一种由英国力学家Argyria六十年代提出的“非寻常有限元”才满足了第三条件。尽管“非寻常有限元”不是主流学派，我对此并不陌生。在80年代初，Argyria的弟子陈兆鎏博士来华讲学时，我听过他的报告。因此，我立刻告诉我同事，非寻常有限元的一个特点就是单元参考坐标系总是建立在单元变形后形状上的。所以不需要特地再建立运动参考系也能通过第三条件。不过我同事的这段经历也给了我一个启发，采用非寻常有限元作为柔性多体方法的基本元素也许更便于柔性多体系统的建模与计算。

　　有些柔性多体方法采用隐式时间积分方法来求数值解，这时要满足第三条件困难就更大了。这时，尽管是小应变，你也必须用建立在格林应变与科希霍夫应力基础上的本构关系来建立刚度阵，否则别想通过第三条件这一关。除非该柔性多体系统是病态的，一般不会用隐式时间积分方法来求解。可是，在九十年代中期，我遇见了一群土木工程学者。他们开辟了一个新领域，称为???，专门研究岩石结构的稳定性，包括失稳前和失稳后运动。该方法用隐式时间积分方法来求解，尽管他们所处理的柔性多体系统并不病态。在我看来这是一个标准的柔性多体问题。因此，我在读他们的论文时，自然要验证一下他们的方法是否满足patch test三条件。结论是，除了个别例外，大多数已发表的方法不满足patch test的第三条件。考虑到岩石结构在失稳发生前，大都是小转动，patch test第三条件对岩石结构的失稳前期运动也许不很重要，因此也就得过且过了吧。不过我是不会用他们的方法的。

　　最后想对读者说句话，若你发现我的系列文章中有错，请向我指正。来信请寄xiangjunqiu@yahoo.com。我将表示感谢。

　　来自科研在线