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

基于精确波效应的多刚体碰撞仿真方法
【专利摘要】本发明涉及一种基于精确波效应的多刚体碰撞仿真方法。该方法在多体碰撞系统的基本物理准则约束基础上，增加波效应约束，所述波效应约束是指在碰撞过程中，应力波携带着动量和能量在相互接触的物体间传播；然后根据多体碰撞系统中接触点处能量的瞬时变化过程，采用冲量的积分代替仿真时间的积分，并为接触势能构建关于接触冲量的二次能量方程，通过求解该方程计算出碰撞后刚体的可靠的运动状态，从而实现多体碰撞过程的仿真。本发明能够精确模拟刚体碰撞过程中蕴含的应力波效应，从而仿真出非常精细的刚体多体碰撞现象。
【专利说明】
基于精确波效应的多刚体碰撞仿真方法
技术领域
[0001] 本发明属于计算机图形学、仿真动画技术领域，具体涉及一种基于精确波效应的 多刚体碰撞仿真计算及其动画方法。
【背景技术】
[0002] 为多刚体碰撞问题进行建模与计算求解是刚体仿真领域中的核心问题。从物理学 角度观察，物体碰撞的过程可以依次分为压缩和松弛这两个阶段。然而在计算机仿真应用 中，仿真算法所能获得的往往只有离散的物体位姿信息，用离散的刚体信息仿真这两个连 续的阶段是一件难以完成的事情。因此，在大多数图形与动画领域的仿真模型中，刚体之间 的碰撞被当作瞬间发生的事件来处理。这些模型假设碰撞持续的时间为无穷小，仅关注碰 撞前与碰撞后这两个瞬间的状态而忽略其内在的碰撞过程与阶段转换。然而这样的假设会 导致一系列不适定的问题，比如碰撞发生次序的问题。同时模型和传播模型是针对这一问 题的两种处理方式。同时模型认为所有的碰撞是在同一瞬间发生的，因此应当统一地进行 求解。传播模型认为碰撞是次序发生的，因此需要逐个碰撞点按次序求解。
[0003] 同时和传播模型有不同的适用场合，也都各自存在缺陷。Smith等人 [SMITH, B. , KAUFMAN, D. Μ., V0UGA, Ε. , TAMSTORF, R. , AND GRINSPUN, Ε. 2012. Reflections on simultaneous impact. ACM Trans. Graph. 31, 4 (July) ,106:1 - 106:12·]提出了五个物理 准则约束，作为衡量刚体仿真算法模型正确性的准则。这五个约束是：分离约束、对称性约 束、能量约束、动量约束以及单边冲量约束。
[0004] 在图形与动画领域，线性互补规划（Linear Complementary Problem，LCP)模型 是最常用的刚体仿真模型。然而它无法正确处理牛顿摆的仿真，因此不满足分离约束。LCP 错误地采用了碰撞前的相对速度来预测碰撞后的相对速度，这会导致一系列问题。例如当 碰撞前的相对速度为零，即两个物体相对静止地接触在一起时，LCP计算得出的碰撞后相 对速度依然为零，即使在这两个物体间有碰撞冲击发生。这使得LCP模型无法满足分离约 束，并且并不适合用于多体碰撞中。为了使其满足分离约束，近年来有一些针对LCP的改 进。Smith等人（参考文献同上）在LCP模型的基础上进行了改进，提出了 Generalized Reflection (GR)模型。GR模型可以在保证其他约束不被破坏的前提下，满足分离约束。
[0005] 然而归根溯源，分离现象产生的根本原因是在碰撞物体中应力波的传递。为了从 更本质的角度描述这一现象，本发明提出了一个新的物理准则约束：波效应约束。波效应约 束是指：在碰撞过程中，应力波携带着动量和能量在相互接触的物体间传播。牛顿摆中的分 离现象就是一种典型的由于波效应造成的现象。波效应也在实际的物理实验中被观察到， 比如在图1中的三球碰撞的场景。理想条件下，碰撞后两侧的球交换速度而中间的球依然 保持静止。
[0006] 原则上来说，分离约束是波效应约束的一种特例。由于其对于碰撞后速度的错误 预测，LCP无法满足分离约束，更无法满足波效应约束。GR在对于碰撞后速度的预测上采用 了与LCP基本相同的策略，因此在处理图1中的情况时，GR会得到与LCP-样的结果。这 个结果是与现实不符的。因此，GR也无法满足波效应约束。
[0007] 综上所述，多刚体碰撞仿真问题是图形仿真与动画领域中一个富有挑战性的问 题，其主要难点在于为同时多点接触的情况提出一个精确适用的算法模型和仿真计算方 法。理想的用于描述多体碰撞的算法模型应当满足一系列基本的物理准则约束，例如能量 约束、对称性约束等。然而现有的算法模型均无法同时满足这些约束。

【发明内容】

[0008] 本发明针对多刚体碰撞仿真问题，提出了一种基于精确波效应的多刚体碰撞仿真 方法，采用二次接触能量模型（QCE，Quadratic Contact Energy)描述碰撞发生时物体与物 体之间表面形成的接触所产生的能量变化，能够精确模拟刚体碰撞过程中蕴含的应力波效 应，从而仿真出非常精细的刚体多体碰撞现象。
[0009] 为实现上述目的，本发明采用如下技术方案：
[0010] 一种基于精确波效应的多刚体碰撞仿真方法，其步骤包括：
[0011] 1)在多体碰撞系统的基本物理准则约束基础上，增加波效应约束，所述波效应约 束是指在碰撞过程中，应力波携带着动量和能量在相互接触的物体间传播；
[0012] 2)根据多体碰撞系统中接触点处能量的瞬时变化过程，采用冲量的积分代替仿真 时间的积分，并为接触势能构建关于接触冲量的二次能量方程，通过求解该方程计算出碰 撞后刚体的可靠的运动状态，从而实现多体碰撞过程的仿真。
[0013] 进一步地，步骤2)将接触碰撞过程中的势能函数E1函数的自变量由时间t替换 为该点的接触冲量Pi，由于冲量Pi可以表示为接触应力与时间的函数，即dP^Fid，并且相 对速度Vl可以表示为广义速度与雅各布矩阵的乘积，即vJ lCl，因此接触势能的公式变换 为：
[0015] ?1是接触点i处产生的接触冲量，它在碰撞的过程中不断增长，1是接触点i的雅 各布矩阵，是广义速度；通过这次变换，得以用接触冲量而非接触时间作为自变量来分析 接触势能的变化，从而避免由于对远小于仿真时间步长的碰撞接触持续时间做微分而引入 的较大误差。
[0016] 进一步地，步骤2)所述关于接触冲量的二次能量方程为：
[0018] 其中，E1是接触点i处的接触势能，P 1是接触点i处产生的接触冲量，J 1是接触点 i的雅各布矩阵，J是系统中所有接触点的雅各布矩阵，爲:是初始的广义速度，Μ是广义质量 矩阵，R1是表示接触冲量Ρ中每一个元素与P i关系的比例向量；
[0019] 通过求解上述一元二次方程，得到使接触势能向量中的所有元素归为零的每个接 触点的接触冲量，从而得到整个系统碰撞后的运动状态。
[0020] 进一步地，对于接触点经历多次压缩-松弛阶段的复杂情况，采用循环迭代的方 式将多次压缩-松弛阶段分离开来，在每一次循环中，一个接触点只能经历一次压缩-松弛 过程；如果某个点在计算完成后，其相对速度再次变为负数，即需要进入再压缩，该点将在 该次循环中被暂时忽略，其压缩阶段将被推迟到下次循环中执行。在为一个接触点解出其 接触冲量以后，重新构建所有与其属于同一个物体的接触点的能量方程。
[0021] 与现有技术相比，本发明的有益效果如下：
[0022] 本发明精确刻画了多体碰撞系统中接触点处能量的变化过程，并为接触势能构建 了关于接触冲量的二次方程，通过求解该方程，计算出可靠的碰撞后刚体的运动状态。该动 力学模型和仿真方法可以满足所有的基本物理约束，并且可以仿真出多种复杂碰撞现象， 特别是刚体碰撞中产生的波效应。除此之外，本发明的模型还具有高兼容性，可以较为容易 地与其它传统模型如LCP (线性互补规划）模型结合起来，结合后的模型可以为带恢复系数 和摩擦系数的刚体仿真提供更准确的结果。本发明中二次能量方程可以解析求解，避免了 数值求解带来的误差和庞大时间开销，可以快速并且准确地对同时多点接触的情况进行仿 真，仿真结果与现实情况高度吻合。
【附图说明】
[0023] 图1是刚体碰撞中由于波效应导致的速度交换现象示意图。A，B，C是三个质量 相同的球体，碰撞前，B球静止，A球和C球以完全相反的方向运动，A的速率大于C的速率 (速度的绝对值）；A球和C球完全同时正碰上B球，碰撞的瞬间状态分成两部分（接触前、 接触后），在接触后，A球和C球的运动状态发生了改变，即A和C的运动方向完全调换，同 时A和C的运动速率也发生了调换，也就是速度交换现象。
[0024] 图2是本发明模拟的经典的牛顿摆碰撞试验效果图。
[0025] 图3是本发明模拟的大量复杂物体相互碰撞的场景。
[0026] 图4是本发明仿真得到的震动台沙砾震动的奇特规则模板效果图。
【具体实施方式】
[0027] 为使本发明的上述目的、特征和优点能够更加明显易懂，下面通过具体实施例对 本发明做进一步说明。
[0028] 1.波效应理论基础
[0029] 波效应是刚体碰撞中的一种常见效应，然而对其进行仿真并不容易。在牛顿摆这 一例子中，对碰撞后速度的不正确预估使得LCP未能将原本贴在一起的物体分开。考虑一 个有η个刚体的系统，其广义坐标可以表示为q s ,广义质量矩阵为M S 1^^·。在此系 统中有m个接触点，这些接触点组成了一个活跃集Α。系统中物体的恢复系数为。对于此系 统，其典型的LCP方程为：
[0031] 公式里K 足系统中所有碰撞点的雅各布矩阵，λ是包含所有冲量大小 的向量。If-:是用于预测碰撞后相对速度的向量，在LCP中它被简单地用-/4,即碰撞前相 对速度的反向速度来代替。这一替换与物理实验中观察到的现象并不吻合，比如牛顿摆实 验。因此，LCP无法满足分离约束和波效应约束。
[0032] 为 了弥补 LCP 的缺陷，Smith 等人提出了 GR 模型[5]\〇1'!1，8.，1^耶]\^，0· M. , VOUGA, E. , TAMSTORF, R. , AND GRINSPUN, E. 2012. Reflections on simultaneous impact. ACM Trans. Graph. 31，4 (July)，106:1 - 106:12.]。在 GR 使用一个新集合 V 代替活 跃集A，只有广义法向与广义速度方向相反的接触点才会被添加进V中。直观地来 说，V只包含相对速度为负值（即相互接近中）的接触点。在牛顿摆实验中，在应力波的传 递过程中，任意时刻只会有一个接触点拥有负相对速度。因此在GR的每次迭代中，V只包 含一个接触点，而其他的接触点将在构建LCP方程是暂时被忽略掉。GR的算法核心就是利 用集合V进行迭代，不断构建LCP方程并求解。通过这个方法，GR能够满足分离约束。
[0033] 通过暂时忽略相对速度为零的接触点，GR得以仿真出物体分离的相关现象，然而 GR并无法满足波效应约束。举例来说，在图1的场景中，两个接触点都拥有负相对速度，因 此这两个接触点都将被添加到集合V中。在这种情况下，V与A相同，因此对于这个例子GR 得到的结果与LCP的结果完全一致，两者都无法仿真出由波效应导致的速度交换现象。
[0034] 2.本发明技术方案的具体实施方法
[0035] 本发明的贡献在于：提出了一种新的基于精确波效应的二次接触能量模型与仿真 计算方法，用以生成准确的碰撞后物体运动状态；该模型满足以上所有的物理约束，包括 Smith等人提出的五个约束以及本发明提出的波效应约束。由于多刚体碰撞问题的复杂性， 在现有的算法模型中，准确的碰撞后运动状态求解往往被认为是非常耗时并且不必要的。 因此一些物理约束，比如波效应约束就被忽略了。这一问题可以通过引入一个二次接触能 量模型来解决。本发明通过将接触势能精确建模为关于接触冲量的二次函数，详细分析并 计算了碰撞的细节过程，使得仿真结果与现实情况高度吻合。模型中的二次能量方程可以 解析求解，避免了数值求解带来的误差和庞大时间开销。下面进行具体说明。
[0036] 1)刚体接触势能分析与建模
[0037] 为了分析波的传递，碰撞通常被认为是一个拥有过程的事件而非瞬时的事件，物 体之间的碰撞必然伴随着物体表面的接触。在局部形变的作用下，物体的速度随着时间而 变化，应力波也藉此在物体间传播。因此波效应可以通过分析物体速度的变化而计算得出。 受此启发，本发明从分析碰撞细节着手来进行建模。
[0038] 根据碰撞点的相对速度，碰撞过程可以被分为两个阶段：压缩阶段和松弛阶段。在 压缩阶段中，两个物体相互挤压，它们之间的相对速度从一个负值变为零，动能转化为接触 势能。在松弛阶段中，两个物体相互弹开，相对速度从零变为一个正值，弹性势能转化为动 能。在这两个阶段中，局部形变是将动能与弹性势能相互转换并提供接触应力的关键。接 触点处的接触势能与其局部形变存在正相关关系。在碰撞开始时，接触势能为零。其在压 缩阶段逐渐增加，并在松弛阶段逐渐减少。让接触势能重新将为零时，该接触点的碰撞过程 结束。因此，本发明首先将每个接触点的接触势能变化作为求解多刚体碰撞问题的关键。
[0039] 对于任意一个接触点i，它的接触势能Ei可以表示为接触应力所做功的负数：
[0040]
[0041] 其中，EiS接触点i的接触势能，F 接触点i的接触应力，v 接触点i的相对 速度，T为接触的时间。对上式作进一步调整，将E1函数的自变量替换为该点的接触冲量 Ρι。由于冲量Ρι可以表示为接触应力与时间的函数，即dPiiFidt，并且相对速度 Vl可以表 示为广义速度与雅各布矩阵的乘积，即^ =人4,因此接触势能的公式变换为：
[0043] ?1是接触点i处产生的接触冲量，它在碰撞的过程中不断增长。通过这次变换，得 以用接触冲量而非接触时间（碰撞时间或者碰撞接触时间）作为自变量来分析接触势能的 变化。因为碰撞接触持续的时间通常远小于仿真时间步长，对这样微小的量再做微分将不 可避免地引入较大误差，而采用接触冲量来描述这一过程可以避免该问题。当整个系统没 有能量损耗时，接触势能E 1在碰撞结束时降回至零。因此使E i降回至零的接触冲量P i，即 是在碰撞过程中该点最终产生的接触冲量。解多体碰撞系统的关键即是求解出冲量向量P， 使得接触势能向量E中的所有元素归为零，SP
[0044] E(P) = 0 (4)
[0045] 2)接触势能转换成二次接触能量模型及每个接触点的求解
[0046] 为了求解这一势能方程，本发明提出二次接触能量模型。从公式3开始，进行一系 列的推导和变换以进行求解。首先，广义速度#可以分为两部分：初始的广义速度為 1以及由 于接触冲量P引起的速度增量，即
[0047] f ^fs 4 ||~siTF (5)
[0048] 因此接触点i处的接触势能表示为：
[0051] 在公式6中，除去P和外的所有量在碰撞过程中都是常量。为了求出这个非 常量积分P中的每一个元素与关系应当被求解出。
[0052] 为了求解出这个关系，从接触力学中接触应力-穿透深度的方程出发，推导出法 向冲量之间的比例关系：
[0054] 其中γ jl= k 4和k汾别表示第j个点和第i个点处的接触刚度。η表示 接触的类型
[0055] (η = 3/2代表赫兹接触模型，η = 1代表线性弹簧模型）。^是接触点j和i 处接触势能的比例，即Εη= E j (Pj) /?的）
[0056] 通过该比例关系可以完成积分.因此方程6可以被改写成入公式8所示的 隐式方程的形式：
[0057] E1(P1，E)=0 (8)
[0058] 然而，直接求解该隐式方程是非常耗时间的。作为一个复杂的一阶的差分方程，它 只能通过迭代的方式近似求解。在图形学的仿真中，时间开销是非常关键的衡量标准。
[0059] 为了能快速求解方程6,本发明提出一种新的用于描述P中每一个元素与？1关系 的方式。与冲量相关比例（impulse correlation ratio, ICR)类似，假设dPydPjli例是一 个常量。ICR是用于描述冲量之间比例的一个常量，
[0061 ] ICR的值仅仅依赖于系统的固有属性以及碰撞前的相对速度。本发明对ICR进行 了扩展，假设ru= dP /ch这一比例是一个常量。与ICR类似，r u是碰撞前相对速度的比 例。这使得P中每一个元素与Pi的关系变得非常清楚，并可以用一个比例向量R 1来表示， 即：
[0063] 这个比例关系反映了一个事实，即碰撞前相对速度的比例影响了碰撞的结果。现 实物理世界中，具有更高碰撞前速度的接触点能够产生更高的接触应力，因此其接触冲量 的微分dPi也更高，这是与本发明提出的这一比例关系相符合的。因此，用碰撞前相对速度 的比例来估计P中每一个元素与？ 1的关系不仅仅是一种对于ICR的扩充，更是符合直觉、并 且有理可循的。需要注意的是，接触点i处的碰撞前相对速度有可能为零。然而作为比例， r]1并不允许其分母上出现零。为了解决这一问题，本发明实际使用的是碰撞前相对速度， 而非它们的比例。这一替换将在下文中详细讨论，这里依然使用比例向量R 1来解释本发明 的模型。
[0064] 有了这一个比例向量，方程6中的E1就可以被转换为一种更加简单的形式：
[0068] 可以看出，i点的接触势能被化成了一个关于其接触冲量？1的二次函数。找出让 EJ3为零的P 再是一个复杂的隐式方程求解问题，而是一个简单的一元二次方程求解问 题。这个一元二次方程可以很容易地进行解析求解，避免了数值方法带来的昂贵的时间开 销和数值误差。对方程11进行求解可以得到每个点的接触冲量，因此整个系统的碰撞后运 动状态就可以被计算得出。
[0069] 3)多次压缩-松弛的处理与计算
[0070] 某些复杂的情况下，接触点将会经历多次压缩-松弛阶段。这将为计算该点相关 的比例r带来一些困难。对于一个刚刚进入再压缩阶段的接触点，其相对速度V。仅仅略大 于零，因此，当将这个相对速度代入公式10来计算相应的比例时，算出的比例将会过大（当 其位于分母时）或过小（当其位于分子时）。这将为接下来的计算带来误差。
[0071] 为了解决这一问题，本发明采用循环迭代的方式来将多次压缩-松弛阶段分离开 来。在每一次循环中，一个接触点只能经历一次压缩-松弛过程。如果某个点在计算完成 后，其相对速度再次变为负数，即需要进入再压缩，该点将在该次循环中被暂时忽略，其压 缩阶段将被推迟到下次循环中执行。通过这个方法，避免了由于再压缩导致的不合适的比 例r。循环的终止条件为整个系统中所有接触点处的相对速度均大于零，此时总的接触冲 量等于每次循环所得的接触冲量之和。若用矣S 代表总接触冲量，用丨# s 代表从第s 次循环中得出的接触冲量，贝1J :
[0073] 其中num为循环的次数。在循环过程中，如果一个接触点（例如j点）依然有接 触势能，那么它就将持续地为两边的物体施加接触应力。在其松弛阶段结束后，其接触势能 降为零并且该点将不再为两边的物体施加力。j点接触应力的消失将对其周围的接触点造 成影响。从公式10可以看出，当dP/变为零时，R 1中的第j个元素将从原本的常数r μ变为 零。这必然将对方程11的求解造成影响。因此，在为j点解出其接触冲量匕以后，所有与 其属于同一个物体的接触点需要重新构建它们的能量方程。
[0074] 4)二次接触能量模型的系统完整求解
[0075] 方程11描述了每个接触点处的接触势能，由此将方程4改写为如下形式：
[0076] E (P) = A 〇 P 〇 P+B 〇 P+C = 0 (13)
[0077] A、B、C这三个向量由所有接触点二次能量函数的系数构成。公式中的。运算符 为Hadamard积，即将两个具有相同维度的矩阵或向量中位置相同的元素两两相乘，并得出 一个新的具有同样维度的矩阵或向量。从方程11中直接计算A、B、C并不容易，因为根据公 式10,不同的接触点拥有不同的比例向量R 1。为了能用统一的形式来进行表示与计算，并 且为了避免在上文提出的由于分母上可能出现零所导致的错误，为所有接触点使用同一的 向量F来代替其各自的R 1。F中的每一个元素是对应接触点碰撞前的相对速度，因此R1中 的第j个元素可以表示为F/G。因此对于任意接触点i，可以得到FiRR 1。在方程11的 等号两侧同时乘上G，可以得到：
[0079] 这次乘法对于结果并没有任何影响，但是它使得A、B、C能够以一种更简单的形式 表达出来。因为
[0081] 因此A、B、C的初值可以被依此计算出：
[0082] A = 1/2JM YF
[0083]
[0084] C = 0 (16)
[0085] 而能量方程变为：
[0086] F 〇 E (P) = 0 (17)
[0087] 在碰撞的过程中，所有接触点的碰撞在同一时刻开始，但是结束的时候并不相同。 一旦一个点结束了它的碰撞过程，它周围的点的能量方程都需要重新构建。因此，碰撞终止 的次序对结果有着决定性的影响。在公式3中时间项已经被消去了，因此需要找到另一种 衡量碰撞终止次序的方法。对于任意一个拥有非零碰撞前速度的接触点i，根据公式10可 以得出P = PJ1，而根据F的定义又可以得到F = RR1，因此
[0088] P = (P./F^F (18)
[0089] R是定值，而Pi在接触的过程中逐渐增加。因此，这个稳定随时间增加的（PyFj 就是所要寻找的类似于时间的自变量。用S来抽象表示，并将公式18写成：
[0090] P = sF (19)
[0091] 对于任意一个接触点k，它拥有对应的使其接触势能重新降为零的冲量Pk。对于 每一个p k，都存在一个确定的&与之对应，类似于标志着碰撞结束的时间
[0092] Pk= s kFk (20)
[0093] sk越小，接触过程就越早结束。通过求解方程15获得所有的s k。然后算法将不断 从中挑选出最小的值，即对应于最早结束碰撞过程的点，将其碰撞过终止并重新构建其周 围接触点的能量方程。重建过程在算法伪码中的19到24行。需要注意的是，算法中的一 元二次方程总会有两个根6和r 2, s r 2。本发明选择较大的根r2作为实际的解，因 为s是单调递增的。伪码：二次接触能量模型的仿真计算过程

[0096] 为了加速整个算法，本发明采用最小堆而不是数组来作为S的数据结构。
[0097] 5)具有恢复系数的碰撞模拟
[0098] 对于完全弹性碰撞而言，使用二次接触能量模型来分析每个接触点的接触势能变 化是非常直观并且有效的。在非完全弹性碰撞情况下，碰撞恢复过程可以以势能的变化来 建模，而势能变化正是本发明的模型所着重关注的。因此，将恢复系数直接添加到本发明的 模型中似乎并无不妥。然而，当采取了与GR类似的循环迭代策略后，本发明的模型也与其 一样无法避免非弹性坍塌的问题，这一问题将导致在非完全弹性碰撞情况下，算法迭代次 数迅速增加，甚至无法达到收敛。这在仿真领域是一个一直存在的问题，并且没有很好的解 决方法。
[0099] 由于直接加入恢复系数会导致非弹性坍塌的问题，本发明提出了一种新的方法， 将本发明的二次接触能量模型与LCP模型结合到一起。二次接触能量模型能够在完全弹 性碰撞的情况下得出可靠的结果，而LCP能够正确求解完全非弹性碰撞[ANITES⑶，M.，AND POTRA, F. A. 1997.Formulating dynamic multi-rigid-body contact problems with friction as solvable linear complementarity problems. Nonlinear Dynamics 14, 231 - 247. STEWART, D. 2000. Rigid-body dynamics with friction and impact. SIAM Review 42, 1，3 - 39.]。如在上文中所讨论的，对于碰撞后速度的错误预测是导致LCP不正 确结果的主要原因。这个问题无法轻易解决，因为多刚体碰撞条件下整个系统过于复杂，无 法轻易计算碰撞后的预期速度。然而在采用了二次接触能量模型之后，完全弹性碰撞后的 广义速度成为了已知。因此，在带恢复系数的LCP方程中餐/可以替代作为对 碰撞后速度的预期。本发明的二次接触能量模型被称为QCE，将结合后的模型称为QCE-LCP 模型，它的方程为如下形式：
[0101] 在完全弹性碰撞中，λ。是QCE-LCP方程的解。当恢复系数c 0时，方程将退化 到标准LCP的形式。这意味着在完全非弹性碰撞的条件下QCE-LCP能够与LCP得到相同 的结果。当恢复系数4=1时，根据QCE-LCP将得到与二次接触能量模型相同的结果。当 0彡cr彡1时，令表示QCE-LCP的结果，则等于#/ (LCP得出的结果）与(二次 接触能量模型得出的结果）之间的插值：
[0103] QCE-LCP并不仅仅是一个单纯的在二次能量模型与LCP之间的线性混合。其中的 差异在处理摩擦时可以很明显地被观察到。如果仅仅采用简单的线性混合，就必须首先忽 略摩擦用线性插值求解出法向接触冲量，然后在基于这个冲量计算切向摩擦。但是这一过 程破坏了法向接触冲量与切向摩擦冲量之间的親合关系。与之相反，在方程21中，QCE-LCP 将两个模型结合到了一起。
[0104] 图2~图4是采用本发明方法得到的基于精确波效应的多刚体碰撞仿真的效果 图。其中图2是本发明模拟的经典的牛顿摆碰撞试验效果图，图3是本发明模拟的大量复杂 物体相互碰撞的场景，图4是本发明仿真得到的震动台沙砾震动的奇特规则模板效果图。
[0105] 以上实施例仅用以说明本发明的技术方案而非对其进行限制，本领域的普通技术 人员可以对本发明的技术方案进行修改或者等同替换，而不脱离本发明的精神和范围，本 发明的保护范围应以权利要求书所述为准。
【主权项】
1. 一种基于精确波效应的多刚体碰撞仿真方法，其步骤包括： 1) 在多体碰撞系统的基本物理准则约束基础上，增加波效应约束，所述波效应约束是 指在碰撞过程中，应力波携带着动量和能量在相互接触的物体间传播； 2) 根据多体碰撞系统中接触点处能量的瞬时变化过程，采用冲量的积分代替仿真时间 的积分，并为接触势能构建关于接触冲量的二次能量方程，通过求解该方程计算出碰撞后 刚体的可靠的运动状态，从而实现多体碰撞过程的仿真。2. 如权利要求1所述的方法，其特征在于，步骤2)将接触碰撞过程中的势能函数E 1函 数的自变量由时间t替换为该点的接触冲量？1，由于冲量？1可以表示为接触应力与时间的 函数，即F lClt，并且相对速度\可以表示为广义速度与雅各布矩阵的乘积，SP ^ 因此接触势能的公式变换为：?1是接触点i处产生的接触冲量，它在碰撞的过程中不断增长，J 1是接触点i的雅各布 矩阵，t;是广义速度；通过这次变换，得以用接触冲量而非接触时间作为自变量来分析接触 势能的变化，从而避免由于对远小于仿真时间步长的碰撞接触持续时间做微分而引入的较 大误差。3. 如权利要求1或2所述的方法，其特征在于，步骤2)所述关于接触冲量的二次能量 方程为：其中，E1是接触点i处的接触势能，P 1是接触点i处产生的接触冲量，J 1是接触点i的 雅各布矩阵，J是系统中所有接触点的雅各布矩阵，知:是初始的广义速度，Μ是广义质量矩 阵，R1是表示接触冲量Ρ中每一个元素与P i关系的比例向量； 通过求解上述一元二次方程，得到使接触势能向量中的所有元素归为零的每个接触点 的接触冲量，从而得到整个系统碰撞后的运动状态。4. 如权利要求3所述的方法，其特征在于：对于接触点经历多次压缩-松弛阶段的复 杂情况，采用循环迭代的方式将多次压缩-松弛阶段分离开来，在每一次循环中，一个接触 点只能经历一次压缩-松弛过程；如果某个点在计算完成后，其相对速度再次变为负数，即 需要进入再压缩，该点将在该次循环中被暂时忽略，其压缩阶段将被推迟到下次循环中执 行。5. 如权利要求4所述的方法，其特征在于：在为一个接触点解出其接触冲量以后，重新 构建所有与其属于同一个物体的接触点的能量方程。6. 如权利要求5所述的方法，其特征在于，采用s k作为接触碰撞持续的时间来决定碰 撞终止次序，sk的计算方法是：为所有接触点使用同一的向量F来代替其各自的R \ F中的 每一个元素是对应接触点碰撞前的相对速度；对于任意一个接触点k，它拥有对应的使其 接触势能重新降为零的冲量P k，对于每一个Pk，都存在一个确定的81<与之对应，即 Pk= s kFk， 其中，Fk表示接触点k碰撞前的相对速度；s k越小则接触过程就越早结束，获得所有的 %后，不断从中挑选出最小的值，即对应于最早结束碰撞过程的点，将其碰撞过终止并重新 构建其周围接触点的能量方程。7.如权利要求1所述的方法，其特征在于：对于非完全弹性碰撞的情况，通过与LCP模 型结合来解决由于直接加入恢复系数而导致的非弹性坍塌的问题，结合后的方程为：其中，J是系统中所有碰撞点的雅各布矩阵，Μ是广义质量矩阵，λ是包含所有冲量大 小的向量，Α是恢复系数，％是广义速度，是完全弹性碰撞后的广义速度。
【文档编号】G06F17/50GK105868425SQ201510025999
【公开日】2016年8月17日
【申请日】2015年1月19日
【发明人】张天翔, 李胜, 汪国平
【申请人】北京大学