我于1990年起在一家位于美国的矿山设备设计咨询公司从事研发工作，从1997年起从事多体动力学在矿山机械中的应用研究。由于应用领域的特殊性，一般的商业通用软件不能满足要求。于是老板同意让我与一位俄裔科学家合作开发了一套供矿山机械应用的多体动力学软件。历时四年，我们的软件终被市场认可。最后，于2002年一家矿山机械制造公司买断了该软件，同时雇我与俄裔科学家继续从事多体动力学的研究与应用。
在多年的开发应用与研究中，我有机会接触了许多多体应用课题。这里，谨谈一下对多体动力学建模及求解过程中的几个关键问题的体会。供同行参考，并想借此机会听取同行意见。
当一系统的运动物体超过一定数量(>1000)，且物体之间有随机碰撞，运动控制方程的形式必须具有一般性以容纳随机变化的约束条件。习惯有限元的人都喜欢控制方程只含位矢变量的形式。但该形式不便处理随机变化的约束条件。本人偏爱一种称为Discrete Element Method (DEM)所用的控制方程。在该控制方程中，同时呈现了物体的位矢变量及相互约束力。DEM 用类似有限元罚函数方法来处理碰撞约束力并可以方便地处理摩擦力。早期的DEM只能处理多个刚性圆球物体的随机碰撞运动，用途狭窄。近10年来由DEM派生出许多新方法可以处理柔性任意几和形状的多体运动。由于碰撞约束的随机性，每隔几步时间积分必须预测每个物体下几步可能相撞的对象。
若多体中某一柔性体变形属大应变，则必须用建立在格林应变与科希霍夫应力基础上的本构关系来处理该柔性体运动。若柔性体变形属小应变大变形或小应变小变形，用处理大应变的方法就非常不经济了。此时，比较有效的方法是将该柔性体分成几个子结构，每个子结构配备一个运动坐标系。在运动坐标系内，相应子结构的变性属小变形。运动坐标系随相应子结构作刚体运动，而子结构的变形力必须参考其运动坐标系来计算。确定运动坐标系的运动有多种方法，这里不作讨论。熟悉有限元的人都知道“patch test”。其实对 柔性多体系统，也有一套类似的“patch test”来检验柔性多体系统算法的“保刚性”。限于篇幅这里不做讨论。
随着视窗操作系统的发展及互联网的出现，人们对SERVER计算机要求越来越高。其中多CPU共享内存的SERVER计算机应运而生。人们已意识到这种类型的计算机将对计算力学带来深远影响：新型并行算法。人们对并行算法的研究已至少有半个世纪，但都是针对多CPU分布式(非共享)内存的计算机而进行的。相对于旧并行算法，新并行算法更有效，编程更容易。只要学一下C++的Multi-thread技巧,余下的将无师自通。对多体计算来说， 新并行算法考虑如何将多体分配到每个CPU时，只需按物体顺序号码来进行。而旧并行算法必须根据物体的空间位置来进行。本人在多体新并行算法上已有4年多经验，但苦于无人交流。(偶尔与熟悉多体旧并行算法的人交流一下。) 在此呼吁更多的同行加入新并行算法研究的行列。其实多CPU共享内存的SERVER计算机并不贵。在美国2000美元就可以买到两CPU(AMD)共享内存的低档SERVER计算机。
最后，谈谈对多体通用程序的看法。由于多体系统的复杂性多样性远远超过有限元所能处理的结构系统，商业多体程序不可能在有效性与通用性上达到一致。通常是牺牲有效性来满足通用性。对于特殊复杂多体系统，例如矿石球磨机中含百万个物体的多体系统且有随机碰撞，商业多体通用程序根本就无能为力。这时，针对某一类型的多体系统，需要专门研究设计专门的方法自己编程求解。这也是做多体分析的工作者的饭碗所在。以上只是一管之见，与同行交流。
熟悉有限元的人都知道“patch test”，它要求：
-
-
当物体朝一个方向发生均匀变形时，有限元应得到与理论相符的常
应力。
由于当有限元划分趋密，每个单元的应力将趋于常数，因此通过patch test是有限元收敛的充分条件。
毫无疑问，这种patch test也一定适用于柔性多体方法。只是相应的patch test更复杂。对每个柔性体的patch tes要求如下：
-
-
当物体延一个方向发生均匀变形但处于静平衡时，柔性体方法应得到与理论相符的常应力。
-
当初始运动条件(包括初速度及初内应力分布)将使物体处于转动动平衡状态时(例如绕固定轴作匀速转动)，柔性体方法应得到与理论相符的动平衡运动解。
比较上述两组patch test条件，可以发现柔性多体方法比有限元方法的patch test多了第三条件。增加这第三条是希望当每个柔性体自由度数趋于无穷大而柔性多体系统趋于动平衡时，柔性多体方法的解应趋于理论解。也就是说，前两个条件是收敛于静平衡解的充分条件，而第三个条件是收敛于动平衡解的充分条件。必须强调的是，这三个条件并不是柔性多体方法对于一般运动求得的数值解趋于理论解的充分条件，而仅仅是对于动平衡运动的数值解趋于理论解的充分条件。
由于大多数柔性多体方法采用有限元作为基本元素，因此一般都能满足前两个条件。然而这第三个条件对很对新方法构造者来说是一个挑战。目前绝大多数柔性多体方法都采用显式时间积分方法来求数值解，因此只要采用了前一篇文章提到的运动参考坐标系，就可满足第三个条件。
我的一位同事曾经给我讲了他在作博士后时的一段经历。为了课题的需要，他必须自己编制一个柔性多体数值求解程序。他试用了很多常规有限元都无法通过这第三条件试验。原因是他不知道应该同时建立一个运动参考系。(由于该同事是做有限差分法研究的，对有限元了解有限。因此走了弯路。) 最后他用了一种由英国力学家Argyria六十年代提出的“非寻常有限元”才满足了第三条件。
尽管“非寻常有限元”不是主流学派，我对此并不陌生。在80年代初，Argyria的弟子陈兆鎏博士来华讲学时，我听过他的报告。因此，我立刻告诉我同事，非寻常有限元的一个特点就是单元参考坐标系总是建立在单元变形后形状上的。所以不需要特地再建立运动参考系也能通过第三条件。不过我同事的这段经历也给了我一个启发，采用非寻常有限元作为柔性多体方法的基本元素也许更便于柔性多体系统的建模与计算。
有些柔性多体方法采用隐式时间积分方法来求数值解，这时要满足第三条件困难就更大了。这时，尽管是小应变，你也必须用建立在格林应变与科希霍夫应力基础上的本构关系来建立刚度阵，否则别想通过第三条件这一关。除非该柔性多体系统是病态的，一般不会用隐式时间积分方法来求解。
可是，在九十年代中期，我遇见了一群土木工程学者。他们开辟了一个新领域，称为？？？，专门研究岩石结构的稳定性，包括失稳前和失稳后运动。该方法用隐式时间积分方法来求解，尽管他们所处理的柔性多体系统并不病态。在我看来这是一个标准的柔性多体问题。
因此，我在读他们的论文时，自然要验证一下他们的方法是否满足patch test三条件。结论是，除了个别例外，大多数已发表的方法不满足patch test的第三条件。考虑到岩石结构在失稳发生前，大都是小转动，patch test第三条件对岩石结构的失稳前期运动也许不很重要，因此也就得过且过了吧。不过我是不会用他们的方法的。
矿山机械门类繁多，但有一共同点：通过控制某种特定的运动来处理矿石。而这种特定的运动就是矿石与机器相互作用下的多体运动。例如，皮带运输机，SAG与AG球磨机，矿石粉碎机，矿石筛选机，等等。这些机械作矿石处理的目的是将矿石尺寸一步一步变小，最终在30至200微米的范围内，将金属与杂物分离。不言而喻，多体动力学方法在矿山机械的设计与控制方面应起到关键性作用。
说来叫人难以相信，直到八十年代末，多体动力学才开始被矿山机械类大学研究人员重视起来。美国的犹他大学于八十年代末才培养出全美第一个将多体动力学应用于矿石球磨机的博士。但本人发现凡由矿山专业毕业的多体动力学博士力学基础较薄弱。因此，他们的研究进展缓慢，很难深入。
1997年，当我当时的老板决定进军矿山机械多体应用这一领域时，我与我的搭档俄裔科学家轻而易举地超过了当时的对手，并一步走到了世界前沿。这里用一个数据来证明此言不虚。2001年，在“SAG 2001”国际会议上，犹他大学教授作“Key Note"发言时宣布他的研究团队已能模拟含5万个物体的缩小了尺寸的球磨机，并预言两年内达到20万个物体”。他不曾想到，随后我老板发言时宣布，我们已能模拟含百万个物体的全尺寸球磨机。读者若想对矿山机械中的多体运动摸拟有一些感性了解，可以点击阅读原文。
当然，内行要看得是门道。下面作一个简单介绍。
1. 必须有非常有效的随机碰撞处理功能。设想有N个物体，每个物体都可能与其余N-1个物体在下一时刻相撞。这要花多少时间才能完成一次碰撞处理：与N平方成正比！可是我和俄裔科学家所发明的方法仅与N一次方成正比。
2.必须有一套建立在概率论基础上的应力，能量分析方法。同时还要有一套建立在概率论基础上的分析实际矿山机械多体运动随机信号方法以与多体运动数值模拟结果作对比。我们已经发展出了这样一套软件。
3.在模拟矿石粉碎机中的多体运动时，必须非常有效地模拟矿石的破裂。我们已取的成功。
从功能上来说，矿山机械多体运动数值模拟必须为矿山机械设计者提供他门最想要的数据：
对於能耗多大及效率多高问题，多体力学工作者可能觉得不难回答。但寿命多长问题就很难轻松回答了。要知道，大多数矿山机械都是因摩擦磨损或更换部件或寿终正寝的。因此，通过多体运动模拟是否能提供由摩擦磨损决定的寿命数据就摆在了我和我同事的面前了。经过努力，我们最终取得了成功。目前，提供球磨机寿命预测及优化服务已是我们部门最赚钱的业务。关于如何模拟摩擦磨损并预测和优化寿命，可参阅参考文献[4]。
1995年，我参加在美国弗罗里达州召开的土木工程国际会议期间遇见了钟万勰教授，并聆听了他的关于“弹性力学求解新体系”报告。由于报告时间只有20分，这着实难为了听众。我和许多听众一样，如坠雾里。会后，钟教授私下跟我介绍了该新体系的精髓：1)采用了对偶变量 2)上升为哈密尔顿体系。
我当时听不懂他所说的哈氏体系。但对偶变量一词我懂。因为我1991年就用了这一概念解决了重型皮带运输机滚筒在百万磅侧向拉力下的弯曲变形及应力分析。但我不敢将自己的工作与钟教授的相比，因为钟教授的新体系是对铁摩辛科弹性力学旧体系的革命。而我当时解滚筒应力时毫无革人之命的感觉。
半年之后，我终於得到了一本钟教授著的<<弹性力学求解新体系>>，并一口气读完。这才发现，我所用的滚筒应力分析法就是钟教授的新体系方法。具体地说，我用了与钟教授新体系法性质相同的对偶变量，将一个长度坐标作为“时间”变量建立了常微分方程，并用传递矩阵手段求解。当然，差别也是很显著的。钟教授新体系法的抽象性一般性，如分离变量，共轭辛正交，及展开定理等令我看到自己与大师的差别。我的方法只能算作钟教授新体系法在薄壁圆柱壳-变厚度薄板-变厚度平面应力-梁所组成的一个特殊结构系统上的一个具体应用。
这里，我想简单谈一下自己当时怎么会独立想到采用对偶变量这一概念来解滚筒应力问题的。也许我的这段经历对在校学生有一定的启发作用。1983年，在我念硕士研究生期间，听了一门非本专业的“现代控制论”课，对其中的定义在时域上的状态变量概念及传递矩阵方法留下深刻印象。我将其中的状态变量概念及传递矩阵方法与本科生期间上“振动理论”课时学过的一种计算一维结构固有频率方法中的(定义在一维空间域上的)状态变量概念及传递矩阵方法[1]作了对比，发现原来时域上状态变量就是弹性力学一维空间域上的对偶变量。但是，与时域相比，这种一维空间域上的计算结构固有频率方法缺少相应的状态变量的常微分方程。於是，我试着去推导出一维空间域上的状态变量的常微分方程，这实际上就是哈密尔顿方程。
但是，我毕竟不是大师，我即不知道这就是哈密尔顿方程，也看不出它可以推广到三维空间。当时就想，将来可能会有机会用它来解一些一维空间域上的弹性力学问题。这一机会直到1991年才来到。1991年我工作所在的那家美国矿山设备设计咨询公司老板给了我一个项目：找一个能够有效求解皮带运输机滚筒应力的分析方法，并发展成一个软件。当时滚筒制造商不怎么喜欢有限元，因为商业有限元程序在当时的工作站计算机上需要两天的时间才能算得一个滚筒应力解(如今只需在PC机上用1-2小时)。而滚筒制造商在初步设计阶段，要从上百个方案中选出几个方案。因此，对滚筒制造商来说，一两分钟就能完成计算才是理想的有效方法。
由于滚筒结构是薄壁圆柱壳-变厚度薄板-变厚度平面应力-梁所组成轴对称系统，用Fourier级数展开后就成一维空间问题了，我立刻意识到应用对偶变量及传递矩阵方法的绝好机会来了。经过一个多月的方程推导，包括导出特征值及向量，Jordan链等一系列工作，又加上三个月的编程，试算，与已知理论解及有限元数值解反复对比，终于证实这一方法就是滚筒制造商期待的理想方法。1992年，该软件，(命名为PStress)，正式进入市场。1993年为了宣传，老板授权让我发表了该方法的理论[2]。
那么，弹性力学的哈密尔顿体系能否应用于多体动力学系统呢？回答是肯定的，因为我已经在用了。矿山机械中有许多圆柱形滚动部件，工作时与其它部件作接触滚动，有时需要将它们处理成二维柔性体。此时，可以将径向坐标当作“时间”坐标，周向用Fourier级数展开。在得到与每一项Fourier级数相应的传递矩阵后，浓缩到可能的接触边界上形成该柔性体的刚度矩阵。以后的计算与一般柔性体完全一样。
另外，我计划将这种方法用于皮带运输机中的皮带上，期望能获得一种比目前德国科学家提倡的简化方法更有效的新方法。 可惜，这一计划提出八年了，一直没时间实施。 从我自身经验中，我可以预期许多了解弹性力学求解新体系的多体力学工作者一定能找到用武之地的。
前不久，从朋友处借来李连杰主演的<<太极张三丰>>光碟。该故事讲的是两个少林和尚从小一起习武(李连杰饰演师兄)。长大后，两人因违反门规而被同时赶出了少林。来到尘世后，两人因不同的处世态度而发展至反目成仇。但师兄的武功在其弟之下。在倍受挫折后，师兄潜心练功。一次，偶然看见一只球掉入水缸引起水的运动，水在水缸中的运动给了师兄启发而悟出太极真谛。遂将过去练就的硬功改变为柔软如水的太极功夫。最后击败其弟，为民除害。
看了该片后，不知为何使我联想起多体动力学。其实，若将多刚体动力学比喻成“硬功”，则多柔体动力学大概就可比作“太极”了吧。进一步，若有一种方法能将水也处理成”多体“，这种方法岂不是“太极中之太极”？其实，这并不是痴人做梦，确实有一种这样的方法，名叫SPH方法。
1997年在一次商务活动中，我认识了来自澳大利亚的保罗博士。保罗博士向我简单介绍了他所从事的SPH方法原理。该方法屏弃了
流体力学的传统欧拉描述法，而采用拉格朗日描述法，将
流体离散成介于有形与无形之间的SPH颗粒。说它有形是因为每个颗粒在空间都有固定的影响半径。说它无形是因为两颗粒相遇时，不会像固体那样因”边界纠纷“而反弹，而是相互融合。真是你中有我，我中有你。
看了保罗的图片介绍，我立刻被其中的“太极”精神吸引了。 我向老板建议研究该方法。 得到老板同意后，我花了一个多月时间读了十几篇论文并编制了程序(相对于传统
流体力学算法，SPH方法因无需网格而编程非常容易)。 可是当我开始算题时，问题一个接一个冒出来了。 原来，这种方法有非常强的数值不稳定因素。 我仔细研究了两个颗粒间的作用力，发现当改变它们之间的距离时，相互作用会从引力转变成斥力，也会从斥力转变成引力。 个中原因我无法得到物理解释。 这时除了怪自己愚钝而不能领悟其“太极”真谛外，只有听凭老板训斥了。
谁知，老板比我还执着，他派我前往澳大利亚当面请教保罗博士。 在澳大利亚， 我无心欣赏异国风光，每天与保罗博士交谈，期望提高修炼程次。 当我问他为何相互作用会从引/斥力转变成斥/引力时，保罗的回答让我摸不到头脑： "The force is as it is"。 听他的话里好像我根本不该问这问题。 看来我还没入“太极”之门呢。 本来计划拜访住在同一城市的保罗的SPH启蒙老师蒙那汉(Monaghan)教授。 蒙教授是SPH方法的开创人之一，机会难得。 可惜，正好蒙教授度假不在。 此时，也只好自我安慰:避免了一次被“蒙”。 于是，我只得向保罗请教。 白天问保罗问题，晚上编程序。 结果发现只有对特定的问题选择特定的参数才可避免数值不稳定性。 但参数的选择并无可循的规律，完全靠反复试算得到。 这对学术研究来说，可能不是问题。 但对以赢利为目的的商业活动来说，是不能允许的。
或许“太极”本身只许正人君子探索，而不许唯利是图的凡夫俗子糟塌？不过，在澳大利亚的经历也不是一无是处。有一段经历值得一提。我发现，凡数值不稳定发生时，众多颗粒总会因引力而聚成一团而不再分开。于是我凭固体力学的直觉，在SPH颗粒间人为地加入类似固体边界碰撞时应有的约束力以防止SPH颗粒聚团。结果，不稳定现象就消失了。但同时我又觉得这样做毫无道理：流体就是流体，不是固体。于是我也羞于向保罗提及这一尝试，以免授人笑柄。后来回到公司，向老板作了汇报，并告之结论：SPH方法目前还不适合商业应用。老板虽然失望，但也无奈，只得望利兴叹。
谈到此，相信大家和我老板当时的心情一样失望。然而，7年后的今天，问题有了转机。今年一月的一天，我的俄籍同事找到我说，“我读了一篇蒙那汉教授2000年发表的文章。他人为地引入了一种斥力而彻底避免了数值不稳定性......”，我没等他说完，就迫不及待地从他手中拿过那篇文章，一口气读到底。天哪，蒙教授用的方法与我当时在澳大利亚羞于向保罗提及的尝试非常相似。蒙教授还举出许多算例来证明改进方法的可行性。我开始懊悔。当时，自己已进入“太极”之门，却胆怯地离去，而且一去七年。可能这就是凡夫俗子与大师的区别。
在懊恼的同时，我的眼前又浮现出一片利润的风光： 一种基于SPH与多体动力学的新型流固耦合方法，和一种基于SPH与多体动力学的新型流固两相流方法(在我看来流固两相流也是一种流固(多体)耦合问题)，它们可以方便地模拟含水球磨机，泥浆泵设计，水甸塘，甚至有可能模拟河床泥沙淤积，等等。 要知道这些课题背后都有丰厚的利润。 前两个项目直接针对着我们公司的产品。 后两个项目可是目前国内水利界最重要的科研课题。 太诱惑了。 打住，我的凡夫俗子之心又开始蠢蠢欲动了。 不过，动就让它动吧，反正也没指望成“太极”高手了。
钟爱多体动力学的人应当正视一个事实：多体动力学这一名称对某些正从事这一工作的人来说缺乏吸引力。也许，你不同意这一观点，也许你怀疑。不要急，听我慢慢道来。
据我所知，在美国至少有这样几个力学工作者群体从事着名符其实的多体动力学工作，但是却从不称自己是多体力学工作者，而是另立门户。例如，研究颗粒状物体运动(granule flower)的群体称他们所用的多体动力学方法为离散单元法(Discrete Element Method)或DEM。美国三大汽车公司都雇用力学工作者专门研究如何做汽车撞击破坏的数值模拟，他们称自己的研究领域为“Crashworthiness”。研究如何模拟金属在冲床挤压下成形的群体称自己的领域为“金属成形(Metal Forming)“。还有我第二篇文章中提到的研究岩石结构稳定性的群体也自立了门户。
另外值得一提的是橡胶轮胎与路面相互作用的数值模拟也有相当多的多体动力学成分。不过该领域的专家不提多体动力学也无可非厚，毕竟他们关心的更多的是轮胎的其他问题。在以上提到的几个多体动力学分支领域中，数Crashworthiness群体的专业水平最高。他们要处理的是含大塑性变形，后压曲及有断裂现象的多体碰撞。一位该领域专家告诉我，模拟汽车一秒钟的碰撞，需要花个把月的计算时间。仅此即可窥见其难度的一斑。其次，金属成形模拟及橡胶轮胎撞击障碍物破裂前后的运动模拟也有相当的难度。
一个有趣的现象是，以上提到的几个多体动力学分支领域的专家很少参加以多体动力学命名的国际会议。若问他们的专业应属什么分支，他们多数会回答属有限元分支，甚至称属流体力学分支(例如，研究granule flower的力学工作者就是这么说的)。如果要探究为何造成多体动力学“江山破碎”这样的局面，我说不准，但大概可以归因于正宗的多体动力学太难了。像图论树结构等概念让人望而生畏。其实，“太难“是一个褒义词；若一定要我(并非我意愿地)用一个贬义词，那应该是“僵硬”。
正宗的多体动力学，不管是笛卡尔流派还是拉格朗日流派，都不太适合处理随机碰撞问题。Crashworthiness的专家们都不懂图论树结构等概念，但照样做着高水平的多体动力学模拟。我虽然读过罗伯森的多刚体动力学一书，但当我和同事在开发矿山机械多体动力学软件时，选择的却是DEM类的多体动力学方法。其实，旁门左道(非正宗)多体动力学的法宝就是两样：有限元(处理变形)加上罚函数法(处理约束)，就是这么简单。了解罚函数方法的人应该同意这样的通俗说法，旁门左道多体动力学的罚函数法是一种拉格朗日乘子(约束反力)的软处理法，而正宗多体动力学的笛卡尔流派的约束处理法是一种拉格朗日乘子的硬处理法。软处理法比硬处理法简便得不是一点点。这也是我(被逼)选用贬义词“僵硬”的来源。
当然，事情都不是绝对的，若能让正宗的和旁门左道的多体动力学方法相互交流取长补短，我们将会看到另一副景象。旁门左道多体动力学方法的长处是处理随机约束方便，短处是计算精度受到罚函数的限制，另外计算所费时间也与罚函数大小有关。正宗多体动力学(拉格朗日流派)的长处是处理树结构系统的完整约束方法非常漂亮，短处是不便于处理随机约束。如果一个系统可以分成这样一些子系统，每个子系统内部不存在碰撞，但子系统之间可能存在着随机碰撞。那么，不用我说完，大家都明白应该如何将正宗与旁门左道方法结合起来。不过细心的读者也许注意到，在这种正宗与旁门左道方法相结合的方法中，笛卡尔流派下岗了。
但是，也有一些多体动力学的分支，尽管没有另立门户，却与现代多体动力学其他分支缺少交流。去年回国，从同学处要来一本名为<<兵器多体系统发射动力学>>专著。我是为其书名所吸引。但读后才发现，该书即无正宗多体动力学元素，也无旁门左道多体动力学踪影。该书呈现的主要研究手段是将发射系统简化为几个自由度的线性系统。若作者能将现代计算多体动力学处理有限位移，有限转动，随机碰撞，图像演示等手段应用到发射动力学中去，一定能更上一层楼。
近几年，柔性多体动力学发展势头不错，但目前仍缺少受力学界各方敬重的重量级人物将各个分支流派统合起来(这里”统合“的意思是指，至少各个分支都来参加同一个国际会议)。值得一提的是，目前国内的多体动力学前景看好。这体现在两个方面。一是许多重点大学都有可授博士学位甚至博士后工作站的一般力学专业(美国大学就没有该专业)，而多体动力学又是一般力学专业的重点研究方向。二是一般力学这一名称已与航空航天器的重要力学问题挂钩，从而有了资金的保障。这两个条件为国内多体动力学今后的发展提供了良好的环境。
来源：声振之家
当前暂无评论，小编等你评论哦！
