基于精確波效應(yīng)的多剛體碰撞仿真方法
【專利摘要】本發(fā)明涉及一種基于精確波效應(yīng)的多剛體碰撞仿真方法。該方法在多體碰撞系統(tǒng)的基本物理準(zhǔn)則約束基礎(chǔ)上，增加波效應(yīng)約束，所述波效應(yīng)約束是指在碰撞過程中，應(yīng)力波攜帶著動量和能量在相互接觸的物體間傳播；然后根據(jù)多體碰撞系統(tǒng)中接觸點(diǎn)處能量的瞬時(shí)變化過程，采用沖量的積分代替仿真時(shí)間的積分，并為接觸勢能構(gòu)建關(guān)于接觸沖量的二次能量方程，通過求解該方程計(jì)算出碰撞后剛體的可靠的運(yùn)動狀態(tài)，從而實(shí)現(xiàn)多體碰撞過程的仿真。本發(fā)明能夠精確模擬剛體碰撞過程中蘊(yùn)含的應(yīng)力波效應(yīng)，從而仿真出非常精細(xì)的剛體多體碰撞現(xiàn)象。
【專利說明】
基于精確波效應(yīng)的多剛體碰撞仿真方法
技術(shù)領(lǐng)域
[0001] 本發(fā)明屬于計(jì)算機(jī)圖形學(xué)、仿真動畫技術(shù)領(lǐng)域，具體涉及一種基于精確波效應(yīng)的 多剛體碰撞仿真計(jì)算及其動畫方法。
【背景技術(shù)】
[0002] 為多剛體碰撞問題進(jìn)行建模與計(jì)算求解是剛體仿真領(lǐng)域中的核心問題。從物理學(xué) 角度觀察，物體碰撞的過程可以依次分為壓縮和松弛這兩個(gè)階段。然而在計(jì)算機(jī)仿真應(yīng)用 中，仿真算法所能獲得的往往只有離散的物體位姿信息，用離散的剛體信息仿真這兩個(gè)連 續(xù)的階段是一件難以完成的事情。因此，在大多數(shù)圖形與動畫領(lǐng)域的仿真模型中，剛體之間 的碰撞被當(dāng)作瞬間發(fā)生的事件來處理。這些模型假設(shè)碰撞持續(xù)的時(shí)間為無窮小，僅關(guān)注碰 撞前與碰撞后這兩個(gè)瞬間的狀態(tài)而忽略其內(nèi)在的碰撞過程與階段轉(zhuǎn)換。然而這樣的假設(shè)會 導(dǎo)致一系列不適定的問題，比如碰撞發(fā)生次序的問題。同時(shí)模型和傳播模型是針對這一問 題的兩種處理方式。同時(shí)模型認(rèn)為所有的碰撞是在同一瞬間發(fā)生的，因此應(yīng)當(dāng)統(tǒng)一地進(jìn)行 求解。傳播模型認(rèn)為碰撞是次序發(fā)生的，因此需要逐個(gè)碰撞點(diǎn)按次序求解。
[0003] 同時(shí)和傳播模型有不同的適用場合，也都各自存在缺陷。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·]提出了五個(gè)物理 準(zhǔn)則約束，作為衡量剛體仿真算法模型正確性的準(zhǔn)則。這五個(gè)約束是：分離約束、對稱性約 束、能量約束、動量約束以及單邊沖量約束。
[0004] 在圖形與動畫領(lǐng)域，線性互補(bǔ)規(guī)劃（Linear Complementary Problem，LCP)模型 是最常用的剛體仿真模型。然而它無法正確處理牛頓擺的仿真，因此不滿足分離約束。LCP 錯(cuò)誤地采用了碰撞前的相對速度來預(yù)測碰撞后的相對速度，這會導(dǎo)致一系列問題。例如當(dāng) 碰撞前的相對速度為零，即兩個(gè)物體相對靜止地接觸在一起時(shí)，LCP計(jì)算得出的碰撞后相 對速度依然為零，即使在這兩個(gè)物體間有碰撞沖擊發(fā)生。這使得LCP模型無法滿足分離約 束，并且并不適合用于多體碰撞中。為了使其滿足分離約束，近年來有一些針對LCP的改 進(jìn)。Smith等人（參考文獻(xiàn)同上）在LCP模型的基礎(chǔ)上進(jìn)行了改進(jìn)，提出了 Generalized Reflection (GR)模型。GR模型可以在保證其他約束不被破壞的前提下，滿足分離約束。
[0005] 然而歸根溯源，分離現(xiàn)象產(chǎn)生的根本原因是在碰撞物體中應(yīng)力波的傳遞。為了從 更本質(zhì)的角度描述這一現(xiàn)象，本發(fā)明提出了一個(gè)新的物理準(zhǔn)則約束：波效應(yīng)約束。波效應(yīng)約 束是指：在碰撞過程中，應(yīng)力波攜帶著動量和能量在相互接觸的物體間傳播。牛頓擺中的分 離現(xiàn)象就是一種典型的由于波效應(yīng)造成的現(xiàn)象。波效應(yīng)也在實(shí)際的物理實(shí)驗(yàn)中被觀察到， 比如在圖1中的三球碰撞的場景。理想條件下，碰撞后兩側(cè)的球交換速度而中間的球依然 保持靜止。
[0006] 原則上來說，分離約束是波效應(yīng)約束的一種特例。由于其對于碰撞后速度的錯(cuò)誤 預(yù)測，LCP無法滿足分離約束，更無法滿足波效應(yīng)約束。GR在對于碰撞后速度的預(yù)測上采用 了與LCP基本相同的策略，因此在處理圖1中的情況時(shí)，GR會得到與LCP-樣的結(jié)果。這 個(gè)結(jié)果是與現(xiàn)實(shí)不符的。因此，GR也無法滿足波效應(yīng)約束。
[0007] 綜上所述，多剛體碰撞仿真問題是圖形仿真與動畫領(lǐng)域中一個(gè)富有挑戰(zhàn)性的問 題，其主要難點(diǎn)在于為同時(shí)多點(diǎn)接觸的情況提出一個(gè)精確適用的算法模型和仿真計(jì)算方 法。理想的用于描述多體碰撞的算法模型應(yīng)當(dāng)滿足一系列基本的物理準(zhǔn)則約束，例如能量 約束、對稱性約束等。然而現(xiàn)有的算法模型均無法同時(shí)滿足這些約束。

【發(fā)明內(nèi)容】

[0008] 本發(fā)明針對多剛體碰撞仿真問題，提出了一種基于精確波效應(yīng)的多剛體碰撞仿真 方法，采用二次接觸能量模型（QCE，Quadratic Contact Energy)描述碰撞發(fā)生時(shí)物體與物 體之間表面形成的接觸所產(chǎn)生的能量變化，能夠精確模擬剛體碰撞過程中蘊(yùn)含的應(yīng)力波效 應(yīng)，從而仿真出非常精細(xì)的剛體多體碰撞現(xiàn)象。
[0009] 為實(shí)現(xiàn)上述目的，本發(fā)明采用如下技術(shù)方案：
[0010] 一種基于精確波效應(yīng)的多剛體碰撞仿真方法，其步驟包括：
[0011] 1)在多體碰撞系統(tǒng)的基本物理準(zhǔn)則約束基礎(chǔ)上，增加波效應(yīng)約束，所述波效應(yīng)約 束是指在碰撞過程中，應(yīng)力波攜帶著動量和能量在相互接觸的物體間傳播；
[0012] 2)根據(jù)多體碰撞系統(tǒng)中接觸點(diǎn)處能量的瞬時(shí)變化過程，采用沖量的積分代替仿真 時(shí)間的積分，并為接觸勢能構(gòu)建關(guān)于接觸沖量的二次能量方程，通過求解該方程計(jì)算出碰 撞后剛體的可靠的運(yùn)動狀態(tài)，從而實(shí)現(xiàn)多體碰撞過程的仿真。
[0013] 進(jìn)一步地，步驟2)將接觸碰撞過程中的勢能函數(shù)E1函數(shù)的自變量由時(shí)間t替換 為該點(diǎn)的接觸沖量Pi，由于沖量Pi可以表示為接觸應(yīng)力與時(shí)間的函數(shù)，即dP^Fid，并且相 對速度Vl可以表示為廣義速度與雅各布矩陣的乘積，即vJ lCl，因此接觸勢能的公式變換 為：
[0015] ?1是接觸點(diǎn)i處產(chǎn)生的接觸沖量，它在碰撞的過程中不斷增長，1是接觸點(diǎn)i的雅 各布矩陣，是廣義速度；通過這次變換，得以用接觸沖量而非接觸時(shí)間作為自變量來分析 接觸勢能的變化，從而避免由于對遠(yuǎn)小于仿真時(shí)間步長的碰撞接觸持續(xù)時(shí)間做微分而引入 的較大誤差。
[0016] 進(jìn)一步地，步驟2)所述關(guān)于接觸沖量的二次能量方程為：
[0018] 其中，E1是接觸點(diǎn)i處的接觸勢能，P 1是接觸點(diǎn)i處產(chǎn)生的接觸沖量，J 1是接觸點(diǎn) i的雅各布矩陣，J是系統(tǒng)中所有接觸點(diǎn)的雅各布矩陣，爲(wèi):是初始的廣義速度，Μ是廣義質(zhì)量 矩陣，R1是表示接觸沖量Ρ中每一個(gè)元素與P i關(guān)系的比例向量；
[0019] 通過求解上述一元二次方程，得到使接觸勢能向量中的所有元素歸為零的每個(gè)接 觸點(diǎn)的接觸沖量，從而得到整個(gè)系統(tǒng)碰撞后的運(yùn)動狀態(tài)。
[0020] 進(jìn)一步地，對于接觸點(diǎn)經(jīng)歷多次壓縮-松弛階段的復(fù)雜情況，采用循環(huán)迭代的方 式將多次壓縮-松弛階段分離開來，在每一次循環(huán)中，一個(gè)接觸點(diǎn)只能經(jīng)歷一次壓縮-松弛 過程；如果某個(gè)點(diǎn)在計(jì)算完成后，其相對速度再次變?yōu)樨?fù)數(shù)，即需要進(jìn)入再壓縮，該點(diǎn)將在 該次循環(huán)中被暫時(shí)忽略，其壓縮階段將被推遲到下次循環(huán)中執(zhí)行。在為一個(gè)接觸點(diǎn)解出其 接觸沖量以后，重新構(gòu)建所有與其屬于同一個(gè)物體的接觸點(diǎn)的能量方程。
[0021] 與現(xiàn)有技術(shù)相比，本發(fā)明的有益效果如下：
[0022] 本發(fā)明精確刻畫了多體碰撞系統(tǒng)中接觸點(diǎn)處能量的變化過程，并為接觸勢能構(gòu)建 了關(guān)于接觸沖量的二次方程，通過求解該方程，計(jì)算出可靠的碰撞后剛體的運(yùn)動狀態(tài)。該動 力學(xué)模型和仿真方法可以滿足所有的基本物理約束，并且可以仿真出多種復(fù)雜碰撞現(xiàn)象， 特別是剛體碰撞中產(chǎn)生的波效應(yīng)。除此之外，本發(fā)明的模型還具有高兼容性，可以較為容易 地與其它傳統(tǒng)模型如LCP (線性互補(bǔ)規(guī)劃）模型結(jié)合起來，結(jié)合后的模型可以為帶恢復(fù)系數(shù) 和摩擦系數(shù)的剛體仿真提供更準(zhǔn)確的結(jié)果。本發(fā)明中二次能量方程可以解析求解，避免了 數(shù)值求解帶來的誤差和龐大時(shí)間開銷，可以快速并且準(zhǔn)確地對同時(shí)多點(diǎn)接觸的情況進(jìn)行仿 真，仿真結(jié)果與現(xiàn)實(shí)情況高度吻合。
【附圖說明】
[0023] 圖1是剛體碰撞中由于波效應(yīng)導(dǎo)致的速度交換現(xiàn)象示意圖。A，B，C是三個(gè)質(zhì)量 相同的球體，碰撞前，B球靜止，A球和C球以完全相反的方向運(yùn)動，A的速率大于C的速率 (速度的絕對值）；A球和C球完全同時(shí)正碰上B球，碰撞的瞬間狀態(tài)分成兩部分（接觸前、 接觸后），在接觸后，A球和C球的運(yùn)動狀態(tài)發(fā)生了改變，即A和C的運(yùn)動方向完全調(diào)換，同 時(shí)A和C的運(yùn)動速率也發(fā)生了調(diào)換，也就是速度交換現(xiàn)象。
[0024] 圖2是本發(fā)明模擬的經(jīng)典的牛頓擺碰撞試驗(yàn)效果圖。
[0025] 圖3是本發(fā)明模擬的大量復(fù)雜物體相互碰撞的場景。
[0026] 圖4是本發(fā)明仿真得到的震動臺沙礫震動的奇特規(guī)則模板效果圖。
【具體實(shí)施方式】
[0027] 為使本發(fā)明的上述目的、特征和優(yōu)點(diǎn)能夠更加明顯易懂，下面通過具體實(shí)施例對 本發(fā)明做進(jìn)一步說明。
[0028] 1.波效應(yīng)理論基礎(chǔ)
[0029] 波效應(yīng)是剛體碰撞中的一種常見效應(yīng)，然而對其進(jìn)行仿真并不容易。在牛頓擺這 一例子中，對碰撞后速度的不正確預(yù)估使得LCP未能將原本貼在一起的物體分開?？紤]一 個(gè)有η個(gè)剛體的系統(tǒng)，其廣義坐標(biāo)可以表示為q s ,廣義質(zhì)量矩陣為M S 1^^·。在此系 統(tǒng)中有m個(gè)接觸點(diǎn)，這些接觸點(diǎn)組成了一個(gè)活躍集Α。系統(tǒng)中物體的恢復(fù)系數(shù)為。對于此系 統(tǒng)，其典型的LCP方程為：
[0031] 公式里K 足系統(tǒng)中所有碰撞點(diǎn)的雅各布矩陣，λ是包含所有沖量大小 的向量。If-:是用于預(yù)測碰撞后相對速度的向量，在LCP中它被簡單地用-/4,即碰撞前相 對速度的反向速度來代替。這一替換與物理實(shí)驗(yàn)中觀察到的現(xiàn)象并不吻合，比如牛頓擺實(shí) 驗(yàn)。因此，LCP無法滿足分離約束和波效應(yīng)約束。
[0032] 為 了彌補(bǔ) 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 使用一個(gè)新集合 V 代替活 躍集A，只有廣義法向與廣義速度方向相反的接觸點(diǎn)才會被添加進(jìn)V中。直觀地來 說，V只包含相對速度為負(fù)值（即相互接近中）的接觸點(diǎn)。在牛頓擺實(shí)驗(yàn)中，在應(yīng)力波的傳 遞過程中，任意時(shí)刻只會有一個(gè)接觸點(diǎn)擁有負(fù)相對速度。因此在GR的每次迭代中，V只包 含一個(gè)接觸點(diǎn)，而其他的接觸點(diǎn)將在構(gòu)建LCP方程是暫時(shí)被忽略掉。GR的算法核心就是利 用集合V進(jìn)行迭代，不斷構(gòu)建LCP方程并求解。通過這個(gè)方法，GR能夠滿足分離約束。
[0033] 通過暫時(shí)忽略相對速度為零的接觸點(diǎn)，GR得以仿真出物體分離的相關(guān)現(xiàn)象，然而 GR并無法滿足波效應(yīng)約束。舉例來說，在圖1的場景中，兩個(gè)接觸點(diǎn)都擁有負(fù)相對速度，因 此這兩個(gè)接觸點(diǎn)都將被添加到集合V中。在這種情況下，V與A相同，因此對于這個(gè)例子GR 得到的結(jié)果與LCP的結(jié)果完全一致，兩者都無法仿真出由波效應(yīng)導(dǎo)致的速度交換現(xiàn)象。
[0034] 2.本發(fā)明技術(shù)方案的具體實(shí)施方法
[0035] 本發(fā)明的貢獻(xiàn)在于：提出了一種新的基于精確波效應(yīng)的二次接觸能量模型與仿真 計(jì)算方法，用以生成準(zhǔn)確的碰撞后物體運(yùn)動狀態(tài)；該模型滿足以上所有的物理約束，包括 Smith等人提出的五個(gè)約束以及本發(fā)明提出的波效應(yīng)約束。由于多剛體碰撞問題的復(fù)雜性， 在現(xiàn)有的算法模型中，準(zhǔn)確的碰撞后運(yùn)動狀態(tài)求解往往被認(rèn)為是非常耗時(shí)并且不必要的。 因此一些物理約束，比如波效應(yīng)約束就被忽略了。這一問題可以通過引入一個(gè)二次接觸能 量模型來解決。本發(fā)明通過將接觸勢能精確建模為關(guān)于接觸沖量的二次函數(shù)，詳細(xì)分析并 計(jì)算了碰撞的細(xì)節(jié)過程，使得仿真結(jié)果與現(xiàn)實(shí)情況高度吻合。模型中的二次能量方程可以 解析求解，避免了數(shù)值求解帶來的誤差和龐大時(shí)間開銷。下面進(jìn)行具體說明。
[0036] 1)剛體接觸勢能分析與建模
[0037] 為了分析波的傳遞，碰撞通常被認(rèn)為是一個(gè)擁有過程的事件而非瞬時(shí)的事件，物 體之間的碰撞必然伴隨著物體表面的接觸。在局部形變的作用下，物體的速度隨著時(shí)間而 變化，應(yīng)力波也藉此在物體間傳播。因此波效應(yīng)可以通過分析物體速度的變化而計(jì)算得出。 受此啟發(fā)，本發(fā)明從分析碰撞細(xì)節(jié)著手來進(jìn)行建模。
[0038] 根據(jù)碰撞點(diǎn)的相對速度，碰撞過程可以被分為兩個(gè)階段：壓縮階段和松弛階段。在 壓縮階段中，兩個(gè)物體相互擠壓，它們之間的相對速度從一個(gè)負(fù)值變?yōu)榱?，動能轉(zhuǎn)化為接觸 勢能。在松弛階段中，兩個(gè)物體相互彈開，相對速度從零變?yōu)橐粋€(gè)正值，彈性勢能轉(zhuǎn)化為動 能。在這兩個(gè)階段中，局部形變是將動能與彈性勢能相互轉(zhuǎn)換并提供接觸應(yīng)力的關(guān)鍵。接 觸點(diǎn)處的接觸勢能與其局部形變存在正相關(guān)關(guān)系。在碰撞開始時(shí)，接觸勢能為零。其在壓 縮階段逐漸增加，并在松弛階段逐漸減少。讓接觸勢能重新將為零時(shí)，該接觸點(diǎn)的碰撞過程 結(jié)束。因此，本發(fā)明首先將每個(gè)接觸點(diǎn)的接觸勢能變化作為求解多剛體碰撞問題的關(guān)鍵。
[0039] 對于任意一個(gè)接觸點(diǎn)i，它的接觸勢能Ei可以表示為接觸應(yīng)力所做功的負(fù)數(shù)：
[0040]
[0041] 其中，EiS接觸點(diǎn)i的接觸勢能，F(xiàn) 接觸點(diǎn)i的接觸應(yīng)力，v 接觸點(diǎn)i的相對 速度，T為接觸的時(shí)間。對上式作進(jìn)一步調(diào)整，將E1函數(shù)的自變量替換為該點(diǎn)的接觸沖量 Ρι。由于沖量Ρι可以表示為接觸應(yīng)力與時(shí)間的函數(shù)，即dPiiFidt，并且相對速度 Vl可以表 示為廣義速度與雅各布矩陣的乘積，即^ =人4,因此接觸勢能的公式變換為：
[0043] ?1是接觸點(diǎn)i處產(chǎn)生的接觸沖量，它在碰撞的過程中不斷增長。通過這次變換，得 以用接觸沖量而非接觸時(shí)間（碰撞時(shí)間或者碰撞接觸時(shí)間）作為自變量來分析接觸勢能的 變化。因?yàn)榕鲎步佑|持續(xù)的時(shí)間通常遠(yuǎn)小于仿真時(shí)間步長，對這樣微小的量再做微分將不 可避免地引入較大誤差，而采用接觸沖量來描述這一過程可以避免該問題。當(dāng)整個(gè)系統(tǒng)沒 有能量損耗時(shí)，接觸勢能E 1在碰撞結(jié)束時(shí)降回至零。因此使E i降回至零的接觸沖量P i，即 是在碰撞過程中該點(diǎn)最終產(chǎn)生的接觸沖量。解多體碰撞系統(tǒng)的關(guān)鍵即是求解出沖量向量P， 使得接觸勢能向量E中的所有元素歸為零，SP
[0044] E(P) = 0 (4)
[0045] 2)接觸勢能轉(zhuǎn)換成二次接觸能量模型及每個(gè)接觸點(diǎn)的求解
[0046] 為了求解這一勢能方程，本發(fā)明提出二次接觸能量模型。從公式3開始，進(jìn)行一系 列的推導(dǎo)和變換以進(jìn)行求解。首先，廣義速度#可以分為兩部分：初始的廣義速度為 1以及由 于接觸沖量P引起的速度增量，即
[0047] f ^fs 4 ||~siTF (5)
[0048] 因此接觸點(diǎn)i處的接觸勢能表示為：
[0051] 在公式6中，除去P和外的所有量在碰撞過程中都是常量。為了求出這個(gè)非 常量積分P中的每一個(gè)元素與關(guān)系應(yīng)當(dāng)被求解出。
[0052] 為了求解出這個(gè)關(guān)系，從接觸力學(xué)中接觸應(yīng)力-穿透深度的方程出發(fā)，推導(dǎo)出法 向沖量之間的比例關(guān)系：
[0054] 其中γ jl= k 4和k汾別表示第j個(gè)點(diǎn)和第i個(gè)點(diǎn)處的接觸剛度。η表示 接觸的類型
[0055] (η = 3/2代表赫茲接觸模型，η = 1代表線性彈簧模型）。^是接觸點(diǎn)j和i 處接觸勢能的比例，即Εη= E j (Pj) /?的）
[0056] 通過該比例關(guān)系可以完成積分.因此方程6可以被改寫成入公式8所示的 隱式方程的形式：
[0057] E1(P1，E)=0 (8)
[0058] 然而，直接求解該隱式方程是非常耗時(shí)間的。作為一個(gè)復(fù)雜的一階的差分方程，它 只能通過迭代的方式近似求解。在圖形學(xué)的仿真中，時(shí)間開銷是非常關(guān)鍵的衡量標(biāo)準(zhǔn)。
[0059] 為了能快速求解方程6,本發(fā)明提出一種新的用于描述P中每一個(gè)元素與？1關(guān)系 的方式。與沖量相關(guān)比例（impulse correlation ratio, ICR)類似，假設(shè)dPydPjli例是一 個(gè)常量。ICR是用于描述沖量之間比例的一個(gè)常量，
[0061 ] ICR的值僅僅依賴于系統(tǒng)的固有屬性以及碰撞前的相對速度。本發(fā)明對ICR進(jìn)行 了擴(kuò)展，假設(shè)ru= dP /ch這一比例是一個(gè)常量。與ICR類似，r u是碰撞前相對速度的比 例。這使得P中每一個(gè)元素與Pi的關(guān)系變得非常清楚，并可以用一個(gè)比例向量R 1來表示， 即：
[0063] 這個(gè)比例關(guān)系反映了一個(gè)事實(shí)，即碰撞前相對速度的比例影響了碰撞的結(jié)果?，F(xiàn) 實(shí)物理世界中，具有更高碰撞前速度的接觸點(diǎn)能夠產(chǎn)生更高的接觸應(yīng)力，因此其接觸沖量 的微分dPi也更高，這是與本發(fā)明提出的這一比例關(guān)系相符合的。因此，用碰撞前相對速度 的比例來估計(jì)P中每一個(gè)元素與？ 1的關(guān)系不僅僅是一種對于ICR的擴(kuò)充，更是符合直覺、并 且有理可循的。需要注意的是，接觸點(diǎn)i處的碰撞前相對速度有可能為零。然而作為比例， r]1并不允許其分母上出現(xiàn)零。為了解決這一問題，本發(fā)明實(shí)際使用的是碰撞前相對速度， 而非它們的比例。這一替換將在下文中詳細(xì)討論，這里依然使用比例向量R 1來解釋本發(fā)明 的模型。
[0064] 有了這一個(gè)比例向量，方程6中的E1就可以被轉(zhuǎn)換為一種更加簡單的形式：
[0068] 可以看出，i點(diǎn)的接觸勢能被化成了一個(gè)關(guān)于其接觸沖量？1的二次函數(shù)。找出讓 EJ3為零的P 再是一個(gè)復(fù)雜的隱式方程求解問題，而是一個(gè)簡單的一元二次方程求解問 題。這個(gè)一元二次方程可以很容易地進(jìn)行解析求解，避免了數(shù)值方法帶來的昂貴的時(shí)間開 銷和數(shù)值誤差。對方程11進(jìn)行求解可以得到每個(gè)點(diǎn)的接觸沖量，因此整個(gè)系統(tǒng)的碰撞后運(yùn) 動狀態(tài)就可以被計(jì)算得出。
[0069] 3)多次壓縮-松弛的處理與計(jì)算
[0070] 某些復(fù)雜的情況下，接觸點(diǎn)將會經(jīng)歷多次壓縮-松弛階段。這將為計(jì)算該點(diǎn)相關(guān) 的比例r帶來一些困難。對于一個(gè)剛剛進(jìn)入再壓縮階段的接觸點(diǎn)，其相對速度V。僅僅略大 于零，因此，當(dāng)將這個(gè)相對速度代入公式10來計(jì)算相應(yīng)的比例時(shí)，算出的比例將會過大（當(dāng) 其位于分母時(shí)）或過小（當(dāng)其位于分子時(shí)）。這將為接下來的計(jì)算帶來誤差。
[0071] 為了解決這一問題，本發(fā)明采用循環(huán)迭代的方式來將多次壓縮-松弛階段分離開 來。在每一次循環(huán)中，一個(gè)接觸點(diǎn)只能經(jīng)歷一次壓縮-松弛過程。如果某個(gè)點(diǎn)在計(jì)算完成 后，其相對速度再次變?yōu)樨?fù)數(shù)，即需要進(jìn)入再壓縮，該點(diǎn)將在該次循環(huán)中被暫時(shí)忽略，其壓 縮階段將被推遲到下次循環(huán)中執(zhí)行。通過這個(gè)方法，避免了由于再壓縮導(dǎo)致的不合適的比 例r。循環(huán)的終止條件為整個(gè)系統(tǒng)中所有接觸點(diǎn)處的相對速度均大于零，此時(shí)總的接觸沖 量等于每次循環(huán)所得的接觸沖量之和。若用矣S 代表總接觸沖量，用丨# s 代表從第s 次循環(huán)中得出的接觸沖量，貝1J :
[0073] 其中num為循環(huán)的次數(shù)。在循環(huán)過程中，如果一個(gè)接觸點(diǎn)（例如j點(diǎn)）依然有接 觸勢能，那么它就將持續(xù)地為兩邊的物體施加接觸應(yīng)力。在其松弛階段結(jié)束后，其接觸勢能 降為零并且該點(diǎn)將不再為兩邊的物體施加力。j點(diǎn)接觸應(yīng)力的消失將對其周圍的接觸點(diǎn)造 成影響。從公式10可以看出，當(dāng)dP/變?yōu)榱銜r(shí)，R 1中的第j個(gè)元素將從原本的常數(shù)r μ變?yōu)?零。這必然將對方程11的求解造成影響。因此，在為j點(diǎn)解出其接觸沖量匕以后，所有與 其屬于同一個(gè)物體的接觸點(diǎn)需要重新構(gòu)建它們的能量方程。
[0074] 4)二次接觸能量模型的系統(tǒng)完整求解
[0075] 方程11描述了每個(gè)接觸點(diǎn)處的接觸勢能，由此將方程4改寫為如下形式：
[0076] E (P) = A 〇 P 〇 P+B 〇 P+C = 0 (13)
[0077] A、B、C這三個(gè)向量由所有接觸點(diǎn)二次能量函數(shù)的系數(shù)構(gòu)成。公式中的。運(yùn)算符 為Hadamard積，即將兩個(gè)具有相同維度的矩陣或向量中位置相同的元素兩兩相乘，并得出 一個(gè)新的具有同樣維度的矩陣或向量。從方程11中直接計(jì)算A、B、C并不容易，因?yàn)楦鶕?jù)公 式10,不同的接觸點(diǎn)擁有不同的比例向量R 1。為了能用統(tǒng)一的形式來進(jìn)行表示與計(jì)算，并 且為了避免在上文提出的由于分母上可能出現(xiàn)零所導(dǎo)致的錯(cuò)誤，為所有接觸點(diǎn)使用同一的 向量F來代替其各自的R 1。F中的每一個(gè)元素是對應(yīng)接觸點(diǎn)碰撞前的相對速度，因此R1中 的第j個(gè)元素可以表示為F/G。因此對于任意接觸點(diǎn)i，可以得到FiRR 1。在方程11的 等號兩側(cè)同時(shí)乘上G，可以得到：
[0079] 這次乘法對于結(jié)果并沒有任何影響，但是它使得A、B、C能夠以一種更簡單的形式 表達(dá)出來。因?yàn)?br>[0081] 因此A、B、C的初值可以被依此計(jì)算出：
[0082] A = 1/2JM YF
[0083]
[0084] C = 0 (16)
[0085] 而能量方程變?yōu)椋?br>[0086] F 〇 E (P) = 0 (17)
[0087] 在碰撞的過程中，所有接觸點(diǎn)的碰撞在同一時(shí)刻開始，但是結(jié)束的時(shí)候并不相同。 一旦一個(gè)點(diǎn)結(jié)束了它的碰撞過程，它周圍的點(diǎn)的能量方程都需要重新構(gòu)建。因此，碰撞終止 的次序?qū)Y(jié)果有著決定性的影響。在公式3中時(shí)間項(xiàng)已經(jīng)被消去了，因此需要找到另一種 衡量碰撞終止次序的方法。對于任意一個(gè)擁有非零碰撞前速度的接觸點(diǎn)i，根據(jù)公式10可 以得出P = PJ1，而根據(jù)F的定義又可以得到F = RR1，因此
[0088] P = (P./F^F (18)
[0089] R是定值，而Pi在接觸的過程中逐漸增加。因此，這個(gè)穩(wěn)定隨時(shí)間增加的（PyFj 就是所要尋找的類似于時(shí)間的自變量。用S來抽象表示，并將公式18寫成：
[0090] P = sF (19)
[0091] 對于任意一個(gè)接觸點(diǎn)k，它擁有對應(yīng)的使其接觸勢能重新降為零的沖量Pk。對于 每一個(gè)p k，都存在一個(gè)確定的&與之對應(yīng)，類似于標(biāo)志著碰撞結(jié)束的時(shí)間
[0092] Pk= s kFk (20)
[0093] sk越小，接觸過程就越早結(jié)束。通過求解方程15獲得所有的s k。然后算法將不斷 從中挑選出最小的值，即對應(yīng)于最早結(jié)束碰撞過程的點(diǎn)，將其碰撞過終止并重新構(gòu)建其周 圍接觸點(diǎn)的能量方程。重建過程在算法偽碼中的19到24行。需要注意的是，算法中的一 元二次方程總會有兩個(gè)根6和r 2, s r 2。本發(fā)明選擇較大的根r2作為實(shí)際的解，因 為s是單調(diào)遞增的。偽碼：二次接觸能量模型的仿真計(jì)算過程

[0096] 為了加速整個(gè)算法，本發(fā)明采用最小堆而不是數(shù)組來作為S的數(shù)據(jù)結(jié)構(gòu)。
[0097] 5)具有恢復(fù)系數(shù)的碰撞模擬
[0098] 對于完全彈性碰撞而言，使用二次接觸能量模型來分析每個(gè)接觸點(diǎn)的接觸勢能變 化是非常直觀并且有效的。在非完全彈性碰撞情況下，碰撞恢復(fù)過程可以以勢能的變化來 建模，而勢能變化正是本發(fā)明的模型所著重關(guān)注的。因此，將恢復(fù)系數(shù)直接添加到本發(fā)明的 模型中似乎并無不妥。然而，當(dāng)采取了與GR類似的循環(huán)迭代策略后，本發(fā)明的模型也與其 一樣無法避免非彈性坍塌的問題，這一問題將導(dǎo)致在非完全彈性碰撞情況下，算法迭代次 數(shù)迅速增加，甚至無法達(dá)到收斂。這在仿真領(lǐng)域是一個(gè)一直存在的問題，并且沒有很好的解 決方法。
[0099] 由于直接加入恢復(fù)系數(shù)會導(dǎo)致非彈性坍塌的問題，本發(fā)明提出了一種新的方法， 將本發(fā)明的二次接觸能量模型與LCP模型結(jié)合到一起。二次接觸能量模型能夠在完全彈 性碰撞的情況下得出可靠的結(jié)果，而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.]。如在上文中所討論的，對于碰撞后速度的錯(cuò)誤預(yù)測是導(dǎo)致LCP不正 確結(jié)果的主要原因。這個(gè)問題無法輕易解決，因?yàn)槎鄤傮w碰撞條件下整個(gè)系統(tǒng)過于復(fù)雜，無 法輕易計(jì)算碰撞后的預(yù)期速度。然而在采用了二次接觸能量模型之后，完全彈性碰撞后的 廣義速度成為了已知。因此，在帶恢復(fù)系數(shù)的LCP方程中餐/可以替代作為對 碰撞后速度的預(yù)期。本發(fā)明的二次接觸能量模型被稱為QCE，將結(jié)合后的模型稱為QCE-LCP 模型，它的方程為如下形式：
[0101] 在完全彈性碰撞中，λ。是QCE-LCP方程的解。當(dāng)恢復(fù)系數(shù)c 0時(shí)，方程將退化 到標(biāo)準(zhǔn)LCP的形式。這意味著在完全非彈性碰撞的條件下QCE-LCP能夠與LCP得到相同 的結(jié)果。當(dāng)恢復(fù)系數(shù)4=1時(shí)，根據(jù)QCE-LCP將得到與二次接觸能量模型相同的結(jié)果。當(dāng) 0彡cr彡1時(shí)，令表示QCE-LCP的結(jié)果，則等于#/ (LCP得出的結(jié)果）與(二次 接觸能量模型得出的結(jié)果）之間的插值：
[0103] QCE-LCP并不僅僅是一個(gè)單純的在二次能量模型與LCP之間的線性混合。其中的 差異在處理摩擦?xí)r可以很明顯地被觀察到。如果僅僅采用簡單的線性混合，就必須首先忽 略摩擦用線性插值求解出法向接觸沖量，然后在基于這個(gè)沖量計(jì)算切向摩擦。但是這一過 程破壞了法向接觸沖量與切向摩擦沖量之間的親合關(guān)系。與之相反，在方程21中，QCE-LCP 將兩個(gè)模型結(jié)合到了一起。
[0104] 圖2~圖4是采用本發(fā)明方法得到的基于精確波效應(yīng)的多剛體碰撞仿真的效果 圖。其中圖2是本發(fā)明模擬的經(jīng)典的牛頓擺碰撞試驗(yàn)效果圖，圖3是本發(fā)明模擬的大量復(fù)雜 物體相互碰撞的場景，圖4是本發(fā)明仿真得到的震動臺沙礫震動的奇特規(guī)則模板效果圖。
[0105] 以上實(shí)施例僅用以說明本發(fā)明的技術(shù)方案而非對其進(jìn)行限制，本領(lǐng)域的普通技術(shù) 人員可以對本發(fā)明的技術(shù)方案進(jìn)行修改或者等同替換，而不脫離本發(fā)明的精神和范圍，本 發(fā)明的保護(hù)范圍應(yīng)以權(quán)利要求書所述為準(zhǔn)。
【主權(quán)項(xiàng)】
1. 一種基于精確波效應(yīng)的多剛體碰撞仿真方法，其步驟包括： 1) 在多體碰撞系統(tǒng)的基本物理準(zhǔn)則約束基礎(chǔ)上，增加波效應(yīng)約束，所述波效應(yīng)約束是 指在碰撞過程中，應(yīng)力波攜帶著動量和能量在相互接觸的物體間傳播； 2) 根據(jù)多體碰撞系統(tǒng)中接觸點(diǎn)處能量的瞬時(shí)變化過程，采用沖量的積分代替仿真時(shí)間 的積分，并為接觸勢能構(gòu)建關(guān)于接觸沖量的二次能量方程，通過求解該方程計(jì)算出碰撞后 剛體的可靠的運(yùn)動狀態(tài)，從而實(shí)現(xiàn)多體碰撞過程的仿真。2. 如權(quán)利要求1所述的方法，其特征在于，步驟2)將接觸碰撞過程中的勢能函數(shù)E 1函 數(shù)的自變量由時(shí)間t替換為該點(diǎn)的接觸沖量？1，由于沖量？1可以表示為接觸應(yīng)力與時(shí)間的 函數(shù)，即F lClt，并且相對速度\可以表示為廣義速度與雅各布矩陣的乘積，SP ^ 因此接觸勢能的公式變換為：?1是接觸點(diǎn)i處產(chǎn)生的接觸沖量，它在碰撞的過程中不斷增長，J 1是接觸點(diǎn)i的雅各布 矩陣，t;是廣義速度；通過這次變換，得以用接觸沖量而非接觸時(shí)間作為自變量來分析接觸 勢能的變化，從而避免由于對遠(yuǎn)小于仿真時(shí)間步長的碰撞接觸持續(xù)時(shí)間做微分而引入的較 大誤差。3. 如權(quán)利要求1或2所述的方法，其特征在于，步驟2)所述關(guān)于接觸沖量的二次能量 方程為：其中，E1是接觸點(diǎn)i處的接觸勢能，P 1是接觸點(diǎn)i處產(chǎn)生的接觸沖量，J 1是接觸點(diǎn)i的 雅各布矩陣，J是系統(tǒng)中所有接觸點(diǎn)的雅各布矩陣，知:是初始的廣義速度，Μ是廣義質(zhì)量矩 陣，R1是表示接觸沖量Ρ中每一個(gè)元素與P i關(guān)系的比例向量； 通過求解上述一元二次方程，得到使接觸勢能向量中的所有元素歸為零的每個(gè)接觸點(diǎn) 的接觸沖量，從而得到整個(gè)系統(tǒng)碰撞后的運(yùn)動狀態(tài)。4. 如權(quán)利要求3所述的方法，其特征在于：對于接觸點(diǎn)經(jīng)歷多次壓縮-松弛階段的復(fù) 雜情況，采用循環(huán)迭代的方式將多次壓縮-松弛階段分離開來，在每一次循環(huán)中，一個(gè)接觸 點(diǎn)只能經(jīng)歷一次壓縮-松弛過程；如果某個(gè)點(diǎn)在計(jì)算完成后，其相對速度再次變?yōu)樨?fù)數(shù)，即 需要進(jìn)入再壓縮，該點(diǎn)將在該次循環(huán)中被暫時(shí)忽略，其壓縮階段將被推遲到下次循環(huán)中執(zhí) 行。5. 如權(quán)利要求4所述的方法，其特征在于：在為一個(gè)接觸點(diǎn)解出其接觸沖量以后，重新 構(gòu)建所有與其屬于同一個(gè)物體的接觸點(diǎn)的能量方程。6. 如權(quán)利要求5所述的方法，其特征在于，采用s k作為接觸碰撞持續(xù)的時(shí)間來決定碰 撞終止次序，sk的計(jì)算方法是：為所有接觸點(diǎn)使用同一的向量F來代替其各自的R \ F中的 每一個(gè)元素是對應(yīng)接觸點(diǎn)碰撞前的相對速度；對于任意一個(gè)接觸點(diǎn)k，它擁有對應(yīng)的使其 接觸勢能重新降為零的沖量P k，對于每一個(gè)Pk，都存在一個(gè)確定的81<與之對應(yīng)，即 Pk= s kFk， 其中，F(xiàn)k表示接觸點(diǎn)k碰撞前的相對速度；s k越小則接觸過程就越早結(jié)束，獲得所有的 %后，不斷從中挑選出最小的值，即對應(yīng)于最早結(jié)束碰撞過程的點(diǎn)，將其碰撞過終止并重新 構(gòu)建其周圍接觸點(diǎn)的能量方程。7.如權(quán)利要求1所述的方法，其特征在于：對于非完全彈性碰撞的情況，通過與LCP模 型結(jié)合來解決由于直接加入恢復(fù)系數(shù)而導(dǎo)致的非彈性坍塌的問題，結(jié)合后的方程為：其中，J是系統(tǒng)中所有碰撞點(diǎn)的雅各布矩陣，Μ是廣義質(zhì)量矩陣，λ是包含所有沖量大 小的向量，Α是恢復(fù)系數(shù)，％是廣義速度，是完全彈性碰撞后的廣義速度。
【文檔編號】G06F17/50GK105868425SQ201510025999
【公開日】2016年8月17日
【申請日】2015年1月19日
【發(fā)明人】張?zhí)煜? 李勝, 汪國平
【申請人】北京大學(xué)