CN105868425B - 基于精确波效应的多刚体碰撞仿真方法 - Google Patents

Abstract

本发明涉及一种基于精确波效应的多刚体碰撞仿真方法。该方法在多体碰撞系统的基本物理准则约束基础上，增加波效应约束，所述波效应约束是指在碰撞过程中，应力波携带着动量和能量在相互接触的物体间传播；然后根据多体碰撞系统中接触点处能量的瞬时变化过程，采用冲量的积分代替仿真时间的积分，并为接触势能构建关于接触冲量的二次能量方程，通过求解该方程计算出碰撞后刚体的可靠的运动状态，从而实现多体碰撞过程的仿真。本发明能够精确模拟刚体碰撞过程中蕴含的应力波效应，从而仿真出非常精细的刚体多体碰撞现象。

Description

基于精确波效应的多刚体碰撞仿真方法

技术领域

本发明属于计算机图形学、仿真动画技术领域，具体涉及一种基于精确波效应的多刚体碰撞仿真计算及其动画方法。

背景技术

为多刚体碰撞问题进行建模与计算求解是刚体仿真领域中的核心问题。从物理学角度观察，物体碰撞的过程可以依次分为压缩和松弛这两个阶段。然而在计算机仿真应用中，仿真算法所能获得的往往只有离散的物体位姿信息，用离散的刚体信息仿真这两个连续的阶段是一件难以完成的事情。因此，在大多数图形与动画领域的仿真模型中，刚体之间的碰撞被当作瞬间发生的事件来处理。这些模型假设碰撞持续的时间为无穷小，仅关注碰撞前与碰撞后这两个瞬间的状态而忽略其内在的碰撞过程与阶段转换。然而这样的假设会导致一系列不适定的问题，比如碰撞发生次序的问题。同时模型和传播模型是针对这一问题的两种处理方式。同时模型认为所有的碰撞是在同一瞬间发生的，因此应当统一地进行求解。传播模型认为碰撞是次序发生的，因此需要逐个碰撞点按次序求解。

同时和传播模型有不同的适用场合，也都各自存在缺陷。Smith等人[SMITH,B.,KAUFMAN,D.M.,VOUGA,E.,TAMSTORF,R.,AND GRINSPUN,E.2012.Reflections onsimultaneous impact.ACM Trans.Graph.31,4(July),106:1–106:12.]提出了五个物理准则约束，作为衡量刚体仿真算法模型正确性的准则。这五个约束是：分离约束、对称性约束、能量约束、动量约束以及单边冲量约束。

在图形与动画领域，线性互补规划(Linear Complementary Problem，LCP)模型是最常用的刚体仿真模型。然而它无法正确处理牛顿摆的仿真，因此不满足分离约束。LCP错误地采用了碰撞前的相对速度来预测碰撞后的相对速度，这会导致一系列问题。例如当碰撞前的相对速度为零，即两个物体相对静止地接触在一起时，LCP计算得出的碰撞后相对速度依然为零，即使在这两个物体间有碰撞冲击发生。这使得LCP模型无法满足分离约束，并且并不适合用于多体碰撞中。为了使其满足分离约束，近年来有一些针对LCP的改进。Smith等人(参考文献同上)在LCP模型的基础上进行了改进，提出了Generalized Reflection(GR)模型。GR模型可以在保证其他约束不被破坏的前提下，满足分离约束。

然而归根溯源，分离现象产生的根本原因是在碰撞物体中应力波的传递。为了从更本质的角度描述这一现象，本发明提出了一个新的物理准则约束：波效应约束。波效应约束是指：在碰撞过程中，应力波携带着动量和能量在相互接触的物体间传播。牛顿摆中的分离现象就是一种典型的由于波效应造成的现象。波效应也在实际的物理实验中被观察到，比如在图1中的三球碰撞的场景。理想条件下，碰撞后两侧的球交换速度而中间的球依然保持静止。

原则上来说，分离约束是波效应约束的一种特例。由于其对于碰撞后速度的错误预测，LCP无法满足分离约束，更无法满足波效应约束。GR在对于碰撞后速度的预测上采用了与LCP基本相同的策略，因此在处理图1中的情况时，GR会得到与LCP一样的结果。这个结果是与现实不符的。因此，GR也无法满足波效应约束。

综上所述，多刚体碰撞仿真问题是图形仿真与动画领域中一个富有挑战性的问题，其主要难点在于为同时多点接触的情况提出一个精确适用的算法模型和仿真计算方法。理想的用于描述多体碰撞的算法模型应当满足一系列基本的物理准则约束，例如能量约束、对称性约束等。然而现有的算法模型均无法同时满足这些约束。

发明内容

本发明针对多刚体碰撞仿真问题，提出了一种基于精确波效应的多刚体碰撞仿真方法，采用二次接触能量模型(QCE，Quadratic Contact Energy)描述碰撞发生时物体与物体之间表面形成的接触所产生的能量变化，能够精确模拟刚体碰撞过程中蕴含的应力波效应，从而仿真出非常精细的刚体多体碰撞现象。

为实现上述目的，本发明采用如下技术方案：

一种基于精确波效应的多刚体碰撞仿真方法，其步骤包括：

1)在多体碰撞系统的基本物理准则约束基础上，增加波效应约束，所述波效应约束是指在碰撞过程中，应力波携带着动量和能量在相互接触的物体间传播；

2)根据多体碰撞系统中接触点处能量的瞬时变化过程，采用冲量的积分代替仿真时间的积分，并为接触势能构建关于接触冲量的二次能量方程，通过求解该方程计算出碰撞后刚体的可靠的运动状态，从而实现多体碰撞过程的仿真。

进一步地，步骤2)将接触碰撞过程中的势能函数E_i函数的自变量由时间t替换为该点的接触冲量P_i，由于冲量P_i可以表示为接触应力与时间的函数，即dP_i＝F_id，并且相对速度v_i可以表示为广义速度与雅各布矩阵的乘积，即v_i＝J_id，因此接触势能的公式变换为：

P_i是接触点i处产生的接触冲量，它在碰撞的过程中不断增长，J_i是接触点i的雅各布矩阵，

是广义速度；通过这次变换，得以用接触冲量而非接触时间作为自变量来分析接触势能的变化，从而避免由于对远小于仿真时间步长的碰撞接触持续时间做微分而引入的较大误差。

进一步地，步骤2)所述关于接触冲量的二次能量方程为：

其中，E_i是接触点i处的接触势能，P_i是接触点i处产生的接触冲量，J_i是接触点i的雅各布矩阵，J是系统中所有接触点的雅各布矩阵，

是初始的广义速度，M是广义质量矩阵，Rⁱ是表示接触冲量P中每一个元素与P_i关系的比例向量；

通过求解上述一元二次方程，得到使接触势能向量中的所有元素归为零的每个接触点的接触冲量，从而得到整个系统碰撞后的运动状态。

进一步地，对于接触点经历多次压缩-松弛阶段的复杂情况，采用循环迭代的方式将多次压缩-松弛阶段分离开来，在每一次循环中，一个接触点只能经历一次压缩-松弛过程；如果某个点在计算完成后，其相对速度再次变为负数，即需要进入再压缩，该点将在该次循环中被暂时忽略，其压缩阶段将被推迟到下次循环中执行。在为一个接触点解出其接触冲量以后，重新构建所有与其属于同一个物体的接触点的能量方程。

与现有技术相比，本发明的有益效果如下：

本发明精确刻画了多体碰撞系统中接触点处能量的变化过程，并为接触势能构建了关于接触冲量的二次方程，通过求解该方程，计算出可靠的碰撞后刚体的运动状态。该动力学模型和仿真方法可以满足所有的基本物理约束，并且可以仿真出多种复杂碰撞现象，特别是刚体碰撞中产生的波效应。除此之外，本发明的模型还具有高兼容性，可以较为容易地与其它传统模型如LCP(线性互补规划)模型结合起来，结合后的模型可以为带恢复系数和摩擦系数的刚体仿真提供更准确的结果。本发明中二次能量方程可以解析求解，避免了数值求解带来的误差和庞大时间开销，可以快速并且准确地对同时多点接触的情况进行仿真，仿真结果与现实情况高度吻合。

附图说明

图1是刚体碰撞中由于波效应导致的速度交换现象示意图。A，B，C是三个质量相同的球体，碰撞前，B球静止，A球和C球以完全相反的方向运动，A的速率大于C的速率(速度的绝对值)；A球和C球完全同时正碰上B球，碰撞的瞬间状态分成两部分(接触前、接触后)，在接触后，A球和C球的运动状态发生了改变，即A和C的运动方向完全调换，同时A和C的运动速率也发生了调换，也就是速度交换现象。

图2是本发明模拟的经典的牛顿摆碰撞试验效果图。

图3是本发明模拟的大量复杂物体相互碰撞的场景。

图4是本发明仿真得到的震动台沙砾震动的奇特规则模板效果图。

具体实施方式

为使本发明的上述目的、特征和优点能够更加明显易懂，下面通过具体实施例对本发明做进一步说明。

1.波效应理论基础

波效应是刚体碰撞中的一种常见效应，然而对其进行仿真并不容易。在牛顿摆这一例子中，对碰撞后速度的不正确预估使得LCP未能将原本贴在一起的物体分开。考虑一个有n个刚体的系统，其广义坐标可以表示为

广义质量矩阵为

在此系统中有m个接触点，这些接触点组成了一个活跃集A。系统中物体的恢复系数为。对于此系统，其典型的LCP方程为：

公式里

是系统中所有碰撞点的雅各布矩阵，

是包含所有冲量大小的向量。是用于预测碰撞后相对速度的向量，在LCP中它被简单地用

即碰撞前相对速度的反向速度来代替。这一替换与物理实验中观察到的现象并不吻合，比如牛顿摆实验。因此，LCP无法满足分离约束和波效应约束。

为了弥补LCP的缺陷，Smith等人提出了GR模型[SMITH,B.,KAUFMAN,D.M.,VOUGA,E.,TAMSTORF,R.,AND GRINSPUN,E.2012.Reflections on simultaneous impact.ACMTrans.Graph.31,4(July),106:1–106:12.]。在GR使用一个新集合V代替活跃集A，

只有广义法向与广义速度方向相反的接触点才会被添加进V中。直观地来说，V只包含相对速度为负值(即相互接近中)的接触点。在牛顿摆实验中，在应力波的传递过程中，任意时刻只会有一个接触点拥有负相对速度。因此在GR的每次迭代中，V只包含一个接触点，而其他的接触点将在构建LCP方程是暂时被忽略掉。GR的算法核心就是利用集合V进行迭代，不断构建LCP方程并求解。通过这个方法，GR能够满足分离约束。

通过暂时忽略相对速度为零的接触点，GR得以仿真出物体分离的相关现象，然而GR并无法满足波效应约束。举例来说，在图1的场景中，两个接触点都拥有负相对速度，因此这两个接触点都将被添加到集合V中。在这种情况下，V与A相同，因此对于这个例子GR得到的结果与LCP的结果完全一致，两者都无法仿真出由波效应导致的速度交换现象。

2.本发明技术方案的具体实施方法

本发明的贡献在于：提出了一种新的基于精确波效应的二次接触能量模型与仿真计算方法，用以生成准确的碰撞后物体运动状态；该模型满足以上所有的物理约束，包括Smith等人提出的五个约束以及本发明提出的波效应约束。由于多刚体碰撞问题的复杂性，在现有的算法模型中，准确的碰撞后运动状态求解往往被认为是非常耗时并且不必要的。因此一些物理约束，比如波效应约束就被忽略了。这一问题可以通过引入一个二次接触能量模型来解决。本发明通过将接触势能精确建模为关于接触冲量的二次函数，详细分析并计算了碰撞的细节过程，使得仿真结果与现实情况高度吻合。模型中的二次能量方程可以解析求解，避免了数值求解带来的误差和庞大时间开销。下面进行具体说明。

1)刚体接触势能分析与建模

为了分析波的传递，碰撞通常被认为是一个拥有过程的事件而非瞬时的事件，物体之间的碰撞必然伴随着物体表面的接触。在局部形变的作用下，物体的速度随着时间而变化，应力波也藉此在物体间传播。因此波效应可以通过分析物体速度的变化而计算得出。受此启发，本发明从分析碰撞细节着手来进行建模。

根据碰撞点的相对速度，碰撞过程可以被分为两个阶段：压缩阶段和松弛阶段。在压缩阶段中，两个物体相互挤压，它们之间的相对速度从一个负值变为零，动能转化为接触势能。在松弛阶段中，两个物体相互弹开，相对速度从零变为一个正值，弹性势能转化为动能。在这两个阶段中，局部形变是将动能与弹性势能相互转换并提供接触应力的关键。接触点处的接触势能与其局部形变存在正相关关系。在碰撞开始时，接触势能为零。其在压缩阶段逐渐增加，并在松弛阶段逐渐减少。让接触势能重新将为零时，该接触点的碰撞过程结束。因此，本发明首先将每个接触点的接触势能变化作为求解多刚体碰撞问题的关键。

对于任意一个接触点i，它的接触势能E_i可以表示为接触应力所做功的负数：

其中，E_i为接触点i的接触势能，F_i为接触点i的接触应力，v_i为接触点i的相对速度，T为接触的时间。对上式作进一步调整，将E_i函数的自变量替换为该点的接触冲量P_i。由于冲量P_i可以表示为接触应力与时间的函数，即dP_i＝F_idt，并且相对速度v_i可以表示为广义速度与雅各布矩阵的乘积，即

因此接触势能的公式变换为：

P_i是接触点i处产生的接触冲量，它在碰撞的过程中不断增长。通过这次变换，得以用接触冲量而非接触时间(碰撞时间或者碰撞接触时间)作为自变量来分析接触势能的变化。因为碰撞接触持续的时间通常远小于仿真时间步长，对这样微小的量再做微分将不可避免地引入较大误差，而采用接触冲量来描述这一过程可以避免该问题。当整个系统没有能量损耗时，接触势能E_i在碰撞结束时降回至零。因此使E_i降回至零的接触冲量P_i，即是在碰撞过程中该点最终产生的接触冲量。解多体碰撞系统的关键即是求解出冲量向量P，使得接触势能向量E中的所有元素归为零，即

E(P)＝0 (4)

2)接触势能转换成二次接触能量模型及每个接触点的求解

为了求解这一势能方程，本发明提出二次接触能量模型。从公式3开始，进行一系列的推导和变换以进行求解。首先，广义速度

可以分为两部分：初始的广义速度

以及由于接触冲量P引起的速度增量，即

因此接触点i处的接触势能表示为：

在公式6中，除去P和P_i之外的所有量在碰撞过程中都是常量。为了求出这个非常量积分

P中的每一个元素与P_i的关系应当被求解出。

为了求解出这个关系，从接触力学中接触应力-穿透深度的方程出发，推导出法向冲量之间的比例关系：

其中γ_ji＝k_j/k_i，k_j和k_i分别表示第j个点和第i个点处的接触刚度。η表示接触的类型

(η＝3/2代表赫兹接触模型，η＝1代表线性弹簧模型)。E_ji是接触点j和i处接触势能的比例，即E_ji＝E_j(P_j)/E_i(P_i)

通过该比例关系可以完成积分

因此方程6可以被改写成入公式8所示的隐式方程的形式：

E_i(P_i，E)＝0 (8)

然而，直接求解该隐式方程是非常耗时间的。作为一个复杂的一阶的差分方程，它只能通过迭代的方式近似求解。在图形学的仿真中，时间开销是非常关键的衡量标准。

为了能快速求解方程6，本发明提出一种新的用于描述P中每一个元素与P_i关系的方式。与冲量相关比例(impulse correlation ratio，ICR)类似，假设dP_j/dP_i比例是一个常量。ICR是用于描述冲量之间比例的一个常量，

ICR的值仅仅依赖于系统的固有属性以及碰撞前的相对速度。本发明对ICR进行了扩展，假设r_ji＝dP_j/dP_i这一比例是一个常量。与ICR类似，r_ji是碰撞前相对速度的比例。这使得P中每一个元素与P_i的关系变得非常清楚，并可以用一个比例向量Rⁱ来表示，即：

这个比例关系反映了一个事实，即碰撞前相对速度的比例影响了碰撞的结果。现实物理世界中，具有更高碰撞前速度的接触点能够产生更高的接触应力，因此其接触冲量的微分dP_i也更高，这是与本发明提出的这一比例关系相符合的。因此，用碰撞前相对速度的比例来估计P中每一个元素与P_i的关系不仅仅是一种对于ICR的扩充，更是符合直觉、并且有理可循的。需要注意的是，接触点i处的碰撞前相对速度有可能为零。然而作为比例，r_ji并不允许其分母上出现零。为了解决这一问题，本发明实际使用的是碰撞前相对速度，而非它们的比例。这一替换将在下文中详细讨论，这里依然使用比例向量Rⁱ来解释本发明的模型。

有了这一个比例向量，方程6中的E_i就可以被转换为一种更加简单的形式：

可以看出，i点的接触势能被化成了一个关于其接触冲量P_i的二次函数。找出让E_i归为零的P_i不再是一个复杂的隐式方程求解问题，而是一个简单的一元二次方程求解问题。这个一元二次方程可以很容易地进行解析求解，避免了数值方法带来的昂贵的时间开销和数值误差。对方程11进行求解可以得到每个点的接触冲量，因此整个系统的碰撞后运动状态就可以被计算得出。

3)多次压缩-松弛的处理与计算

某些复杂的情况下，接触点将会经历多次压缩-松弛阶段。这将为计算该点相关的比例r带来一些困难。对于一个刚刚进入再压缩阶段的接触点，其相对速度v₀仅仅略大于零，因此，当将这个相对速度代入公式10来计算相应的比例时，算出的比例将会过大(当其位于分母时)或过小(当其位于分子时)。这将为接下来的计算带来误差。

为了解决这一问题，本发明采用循环迭代的方式来将多次压缩-松弛阶段分离开来。在每一次循环中，一个接触点只能经历一次压缩-松弛过程。如果某个点在计算完成后，其相对速度再次变为负数，即需要进入再压缩，该点将在该次循环中被暂时忽略，其压缩阶段将被推迟到下次循环中执行。通过这个方法，避免了由于再压缩导致的不合适的比例r。循环的终止条件为整个系统中所有接触点处的相对速度均大于零，此时总的接触冲量等于每次循环所得的接触冲量之和。若用

代表总接触冲量，用

代表从第s次循环中得出的接触冲量，则：

其中num为循环的次数。在循环过程中，如果一个接触点(例如j点)依然有接触势能，那么它就将持续地为两边的物体施加接触应力。在其松弛阶段结束后，其接触势能降为零并且该点将不再为两边的物体施加力。j点接触应力的消失将对其周围的接触点造成影响。从公式10可以看出，当dP_j变为零时，Rⁱ中的第j个元素将从原本的常数r_ji变为零。这必然将对方程11的求解造成影响。因此，在为j点解出其接触冲量P_j以后，所有与其属于同一个物体的接触点需要重新构建它们的能量方程。

4)二次接触能量模型的系统完整求解

方程11描述了每个接触点处的接触势能，由此将方程4改写为如下形式：

E(P)＝AοPοP+BοP+C＝0 (13)

A、B、C这三个向量由所有接触点二次能量函数的系数构成。公式中的ο运算符为Hadamard积，即将两个具有相同维度的矩阵或向量中位置相同的元素两两相乘，并得出一个新的具有同样维度的矩阵或向量。从方程11中直接计算A、B、C并不容易，因为根据公式10，不同的接触点拥有不同的比例向量Rⁱ。为了能用统一的形式来进行表示与计算，并且为了避免在上文提出的由于分母上可能出现零所导致的错误，为所有接触点使用同一的向量F来代替其各自的Rⁱ。F中的每一个元素是对应接触点碰撞前的相对速度，因此Rⁱ中的第j个元素可以表示为F_j/F_i。因此对于任意接触点i，可以得到F＝F_iRⁱ。在方程11的等号两侧同时乘上F_i，可以得到：

这次乘法对于结果并没有任何影响，但是它使得A、B、C能够以一种更简单的形式表达出来。因为

因此A、B、C的初值可以被依此计算出：

A＝1/2JM^-1J^TF

C＝0 (16)

而能量方程变为：

FοE(P)＝0 (17)

在碰撞的过程中，所有接触点的碰撞在同一时刻开始，但是结束的时候并不相同。一旦一个点结束了它的碰撞过程，它周围的点的能量方程都需要重新构建。因此，碰撞终止的次序对结果有着决定性的影响。在公式3中时间项已经被消去了，因此需要找到另一种衡量碰撞终止次序的方法。对于任意一个拥有非零碰撞前速度的接触点i，根据公式10可以得出P＝P_iRⁱ，而根据F的定义又可以得到F＝F_iRⁱ，因此

P＝(P_i/F_i)F (18)

F_i是定值，而P_i在接触的过程中逐渐增加。因此，这个稳定随时间增加的(P_i/F_i)就是所要寻找的类似于时间的自变量。用s来抽象表示，并将公式18写成：

P＝sF (19)

对于任意一个接触点k，它拥有对应的使其接触势能重新降为零的冲量P_k。对于每一个P_k，都存在一个确定的s_k与之对应，类似于标志着碰撞结束的时间

P_k＝s_kF_k (20)

s_k越小，接触过程就越早结束。通过求解方程15获得所有的s_k。然后算法将不断从中挑选出最小的值，即对应于最早结束碰撞过程的点，将其碰撞过终止并重新构建其周围接触点的能量方程。重建过程在算法伪码中的19到24行。需要注意的是，算法中的一元二次方程总会有两个根r₁和r₂，r₁≤s_k≤r₂。本发明选择较大的根r₂作为实际的解，因为s是单调递增的。伪码：二次接触能量模型的仿真计算过程

为了加速整个算法，本发明采用最小堆而不是数组来作为S的数据结构。

5)具有恢复系数的碰撞模拟

对于完全弹性碰撞而言，使用二次接触能量模型来分析每个接触点的接触势能变化是非常直观并且有效的。在非完全弹性碰撞情况下，碰撞恢复过程可以以势能的变化来建模，而势能变化正是本发明的模型所着重关注的。因此，将恢复系数直接添加到本发明的模型中似乎并无不妥。然而，当采取了与GR类似的循环迭代策略后，本发明的模型也与其一样无法避免非弹性坍塌的问题，这一问题将导致在非完全弹性碰撞情况下，算法迭代次数迅速增加，甚至无法达到收敛。这在仿真领域是一个一直存在的问题，并且没有很好的解决方法。

由于直接加入恢复系数会导致非弹性坍塌的问题，本发明提出了一种新的方法，将本发明的二次接触能量模型与LCP模型结合到一起。二次接触能量模型能够在完全弹性碰撞的情况下得出可靠的结果，而LCP能够正确求解完全非弹性碰撞[ANITESCU,M.,ANDPOTRA,F.A.1997.Formulating dynamic multi-rigid-body contact problems withfriction as solvable linear complementarity problems.Nonlinear Dynamics 14,231–247.STEWART,D.2000.Rigid-body dynamics with friction and impact.SIAMReview 42,1,3–39.]。如在上文中所讨论的，对于碰撞后速度的错误预测是导致LCP不正确结果的主要原因。这个问题无法轻易解决，因为多刚体碰撞条件下整个系统过于复杂，无法轻易计算碰撞后的预期速度。然而在采用了二次接触能量模型之后，完全弹性碰撞后的广义速度

成为了已知。因此，在带恢复系数的LCP方程中

可以替代作为对碰撞后速度的预期。本发明的二次接触能量模型被称为QCE，将结合后的模型称为QCE-LCP模型，它的方程为如下形式：

在完全弹性碰撞中，λ_c是QCE-LCP方程的解。当恢复系数c_r＝0时，方程将退化到标准LCP的形式。这意味着在完全非弹性碰撞的条件下QCE-LCP能够与LCP得到相同的结果。当恢复系数c_r＝1时，根据QCE-LCP将得到与二次接触能量模型相同的结果。当0≤c_r≤1时，令

表示QCE-LCP的结果，则

等于

(LCP得出的结果)与

(二次接触能量模型得出的结果)之间的插值：

QCE-LCP并不仅仅是一个单纯的在二次能量模型与LCP之间的线性混合。其中的差异在处理摩擦时可以很明显地被观察到。如果仅仅采用简单的线性混合，就必须首先忽略摩擦用线性插值求解出法向接触冲量，然后在基于这个冲量计算切向摩擦。但是这一过程破坏了法向接触冲量与切向摩擦冲量之间的耦合关系。与之相反，在方程21中，QCE-LCP将两个模型结合到了一起。

图2～图4是采用本发明方法得到的基于精确波效应的多刚体碰撞仿真的效果图。其中图2是本发明模拟的经典的牛顿摆碰撞试验效果图，图3是本发明模拟的大量复杂物体相互碰撞的场景，图4是本发明仿真得到的震动台沙砾震动的奇特规则模板效果图。

以上实施例仅用以说明本发明的技术方案而非对其进行限制，本领域的普通技术人员可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明的精神和范围，本发明的保护范围应以权利要求书所述为准。

Claims (5)

1.一种基于精确波效应的多刚体碰撞仿真方法，其步骤包括：

1)在多体碰撞系统的基本物理准则约束基础上，增加波效应约束，所述波效应约束是指在碰撞过程中，应力波携带着动量和能量在相互接触的物体间传播；

2)根据多体碰撞系统中接触点处能量的瞬时变化过程，采用冲量的积分代替仿真时间的积分，并为接触势能构建关于接触冲量的二次能量方程，通过求解该方程计算出碰撞后刚体的可靠的运动状态，从而实现多体碰撞过程的仿真；

所述步骤2)将接触碰撞过程中的势能函数E_i函数的自变量由时间t替换为接触点的接触冲量P_i，由于冲量P_i可以表示为接触应力与时间的函数，即dP_i＝F_idt，并且相对速度v_i可以表示为广义速度与雅各布矩阵的乘积，即

因此接触势能的公式变换为：

P_i是接触点i处产生的接触冲量，它在碰撞的过程中不断增长，J_i是接触点i的雅各布矩阵，

是广义速度；

步骤2)所述关于接触冲量的二次能量方程为：

其中，E_i是接触点i处的接触势能，P_i是接触点i处产生的接触冲量，J_i是接触点i的雅各布矩阵，J是系统中所有接触点的雅各布矩阵，

是初始的广义速度，M是广义质量矩阵，Rⁱ是表示接触冲量P中每一个元素与P_i关系的比例向量；

通过求解所述二次能量方程，得到使接触势能向量中的所有元素归为零的每个接触点的接触冲量，从而得到整个系统碰撞后的运动状态。

2.如权利要求1所述的方法，其特征在于：对于接触点经历多次压缩-松弛阶段的复杂情况，采用循环迭代的方式将多次压缩-松弛阶段分离开来，在每一次循环中，一个接触点只能经历一次压缩-松弛过程；如果某个点在计算完成后，其相对速度再次变为负数，即需要进入再压缩，该点将在该次循环中被暂时忽略，其压缩阶段将被推迟到下次循环中执行。

3.如权利要求2所述的方法，其特征在于：在为一个接触点解出其接触冲量以后，重新构建所有与其属于同一个物体的接触点的能量方程。

4.如权利要求3所述的方法，其特征在于，采用s_k作为接触碰撞持续的时间来决定碰撞终止次序，s_k的计算方法是：为所有接触点使用同一的向量F来代替其各自的Rⁱ，F中的每一个元素是对应接触点碰撞前的相对速度；对于任意一个接触点k，它拥有对应的使其接触势能重新降为零的冲量P_k，对于每一个P_k，都存在一个确定的s_k与之对应，即

P_k＝s_kF_k，

其中，F_k表示接触点k碰撞前的相对速度；s_k越小则接触过程就越早结束，获得所有的s_k后，不断从中挑选出最小的值，即对应于最早结束碰撞过程的点，将其碰撞过终止并重新构建其周围接触点的能量方程。

5.如权利要求1所述的方法，其特征在于：对于非完全弹性碰撞的情况，通过与LCP模型结合来解决由于直接加入恢复系数而导致的非弹性坍塌的问题，结合后的方程为：

其中，J是系统中所有碰撞点的雅各布矩阵，M是广义质量矩阵，λ是包含所有冲量大小的向量，c_r是恢复系数，

是广义速度，

是完全弹性碰撞后的广义速度。

Images