基于細分粒子的剛體破裂動畫

陳沸鑌，謝步瀛，冉修遠

(同濟大學建筑工程系，上海 200092)

日常生活中，脆性破裂是一種很常見的現象，如瓷器的破碎，混凝土磚墻的倒塌。如何對這種現象進行模擬是計算機圖形學領域的一個重要研究方向，在計算機游戲、影視特效、虛擬現實等領域中有著廣泛的應用。

在計算力學和材料力學領域，對于固體脆性破裂數值模擬的研究已經進行了數十年，研究人員提出了大量的數值計算方法。目前大致有兩種數值離散方法：基于網格的方法和無網格方法。有限元法是(finite element method，FEM)一種典型的基于網格的方法，它是通過將連續體模型剖分成一個個網格單元，根據各個網格單元的特性來計算網格節點上的屬性。有限元法的主要缺陷是對網格剖分的要求較高，當由外力產生的網格形變太大的時候，對控制方程求解的精度會產生很大的影響。相對而言，近年來引起研究人員極大興趣的無網格(粒子)方法由于在計算時不受到背景網格的約束可以有效地解決上述基于網格方法的缺陷。然而，在圖形學領域，使用基于物理的方法對脆性破裂進行動畫生成的研究卻為數不多?；谖锢淼墓腆w脆性破裂動畫主要是以斷裂力學以及材料力學為基礎，采用數值離散方法，來分析固體在外力作用下產生的內部應力。

與數值分析方法不同的是，在計算機圖形學中，為了達到相對快速的模擬速度，適當的模型簡化以及對數值方法的進行優化是非常重要的。與此同時，如何增強模擬的細節也是引起研究人員廣泛關注的另一個問題。單方面的追求計算速度往往導致模擬的結果過于粗糙，然而使用分辨率較好的原始模型又會帶來計算代價較大的問題。如何能在提高模擬框架的計算速度的前提下又能夠給出相對較好的模擬精度是一個具有挑戰性的問題。針對這種問題所提出的局部細分技術能夠很好的地解決這一矛盾，局部細分技術的主要思想是以一個相對較大的尺度進行初始離散建模，然后通過在需要增強細節的部分進行局部的細分方式來重采樣生成較小尺度的局部模型以此來提高模擬的細節。

本文提出一種基于細分粒子的剛體脆性破裂的模擬框架。首先將四面體化的固體網格模型綁定到一系列粒子上。然后使用光滑粒子流體動力學(smoothed particle hydrodynamics，SPH)方法來求解固體碰撞時產生的內部應力。通過更新SPH粒子所綁定網格節點的狀態，能夠快速的對開裂斷面進行維護從而減少在網格重構和拓撲更新上的計算代價。最后提出一種細分粒子的算法，從而在提高系統的計算效率的同時又增強了模擬效果的細節。

1 相關工作

與脆性破裂動畫研究有關的工作主要可以從數值計算模型和幾何模型兩個方面來進行闡述。

從數值計算模型的角度來說，在圖形學領域最早進行固體破裂動畫的研究可以追溯到上世紀80年代末。Terzopoulos等[1-2]最早使用基于物理的方法來進行彈塑性變形和破裂模擬，他們使用的數值模型是有限差分方法。同一時期研究較多的模型還有質點-彈簧模型[3]或者質點-約束模型[4]。上述模型的主要缺點是不能精確的計算開裂點的位置和方向，以此模擬的結果并不十分逼真。在破裂模擬中，更為精確的數值模型是基于連續介質力學[5]。作為一種典型的基于網格的數值離散方法，有限元法首次被O'Brien等[6-7]用來求解線彈性力學方程進行脆性和塑形破裂的模擬，從而得到了較為真實的效果。在此基礎上，許多研究人員針對有限元方法進行了各種改進，如對物理模型的簡化[8]，提高模擬的速度[9]以及加強計算的精度[10]。基于網格的數值方法的缺點在于得到計算精度高的結果需要高質量的網格，而無網格方法不存在這種問題。Pauly等[11]使用移動最小二乘法(moving least square，MLS)來求解應變張量從而進行脆性和塑性破裂的模擬。其不足之處在于為了保證計算的正確，每個粒子對其鄰接粒子的位置和數量有較高的要求。Liu等[12]采用了無網格局部 Petrov-Galerkin方法(meshless local petrov-Galerkin，MLPG)解決了這個問題，但是在求解的時候需要計算一個大規模的線性方程組，因此計算的開銷較大。

從幾何模型的角度來說，破裂模擬實際上是對一個整體的網格進行拓撲的切割操作，因此，如何選擇一個合適的固體表面網格模型對模擬速度以及仿真結果有著重要的影響。固體的曲面模型大致可以分成2種：顯式的網格模型與隱式的網格模型。顯式模型在模擬過程中不需要在每個時間步長內反復進行表面網格的生成，只需要處理節點間的關系，因此處理代價低、計算速度高，但由于網格是初始計算時就生成好的，因此進行幾何切割操作較為復雜。如Müller等[9,13]使用的四面體網格模型和規則八面體網格模型、Molino等[14]提出的虛節點算法等，均是對初始的網格做了一定的約束條件之后才進行切割操作的。隱式曲面模型是使用隱式曲面方程來生成固體的表面，對復雜曲面變化的適用性較好而且不存在切割操作的問題，但是由于需要在每個時間步長內進行網格的計算，因此計算代價較高。如文獻[11]中提出的一種基于點的隱式跟蹤曲面模型，其開裂面的生成是通過構造型立體幾何表達法(constructive solid geometry，CSG)進行生成的。

綜上所述，本文脆性破裂模擬的框架提出了一種基于粒子和網格的混合模型。在數值離散方面使用了基于粒子SPH方法進行求解。在幾何方面使用了一種將網格綁定于粒子上幾何模型，能夠快速對開裂固體曲面進行維護操作。

2 基于粒子的固體模型

2.1 固體的離散

本文的數值計算模型是基于粒子的，因此首先需要將初始模型離散成一系列粒子。同時考慮到開裂斷面的生成，一個合適的固體幾何模型對于開裂時幾何拓撲的修改和模擬結果的渲染也是極其重要的?；谶@兩點的考慮，本文使用了一種將離散四面體綁定到粒子上的顯式曲面幾何模型。

在數值計算中，網格生成是一個關鍵性的技術，如何將一個連續的計算區域離散成一系列的四面體是一個重要的研究課題，作為數值計算和計算幾何的交叉研究內容，四面體網格生成技術經過近20年的發展已逐漸成熟并衍生出大量的網格生成技術，在這之中，Delaunay三角剖分由于其網格高質量特征和對復雜區域良好的逼近性受到研究人員的廣泛關注，是目前最流行的全自動網格生成方法之一。因此，本文選用Delaunay三角剖分來進行初始的四面體網格生成，詳細的生成算法可以參見文獻[15]。

首先，將整個固體區域通過Delaunay三角化剖分成一系列的四面體，在每個四面體的形心位置分配一個粒子。粒子的半徑r可根據四面體的體積V進行計算：為了進行后續的拓撲維護操作為每個粒子保存2種信息：原四面體的頂點坐標和指向相鄰粒子的指針(圖1中的“鏈接”)。四面體頂點坐標(圖1中藍色的點)主要有2個作用：維護開裂的曲面的生成以及計算內部應力。當粒子的鄰接粒子數小于4時將該粒子定義為固體表面粒子，其作用是進行碰撞處理和計算虛位移狀態。在初始化時，所有表面粒子所對應的三角形面都將加入到一個三角形網格中作為固體的曲面模型。 SPH粒子所保存的信息可以快速的進行開裂面的維護，將在第4節中具體介紹。

圖1 固體離散建模過程

2.2 線彈性固體力學模型

本文破裂模擬框架的力學模型是基于線彈性固體力學。在線彈性固體力學中，固體的應力與應變是與固體內部某個點的變形狀態有關的。設分別為固體中某個點的初始位置和變形后的位置。由u=[u,v,w]T定義了一個位置的向量場：xt=x0+u。該點的應變可以用位移的梯度場?u來進行計算。令J代表由x0到tx映射的雅克比矩陣。采用線性的柯西-格林應變張量，可以將應變張量表示為：

如果將固體近似考慮成線彈性各項同性材料，應力與應變張量的關系可以定義為：

對于各項同性材料，C是一個四階各項同性張量[5]。

其中E與v分別是材料的彈性模量和泊松比。

2.3 SPH離散求解

SPH是一種無網格粒子法，最早用于天體物理現象的模擬[16]，此后廣泛應用于流體力學和固體力學的數值分析中。在SPH中，每個粒子在連續的標量場f(x)的值是通過核近似法轉化為支持域內所有離散的粒子數據的疊加來計算得到：

相關SPH的詳細內容及公式推導可參見文獻[15]。

通過使用SPH方法可以對2.2節中線彈性固體力學方程來進行離散求解，在2.1節中已經將計算區域離散成與四面體有關的粒子。因此，對固體內部的每個粒子，位移梯度場的SPH近似表達式為：

其中向量uji代表的鄰接粒子uj與當前粒子iu的位移差：

2.4 虛位移計算

采用上述SPH方法可以對固體在外力作用下所產生的內部應力進行分析。在模擬過程中，外力主要是由于固體間的相互碰撞產生的。因此首先要進行的是碰撞的處理。將剛體離散成一系列粒子之后，剛體間的碰撞檢測實際上是由粒子來完成的。在模擬脆性破裂時，可將固體考慮成具有無限大剛性的剛體。為了能夠使用2.2節的線彈性固體力學模型進行分析，需要為發生碰撞的固體粒子計算一個虛位移的狀態來得到碰撞時發生的變形以此來計算應變和應力。對每個發生碰撞的粒子，本文使用粒子的碰撞信息來計算虛位移狀態，如圖2所示，碰撞粒子的虛位移的狀態是由穿刺深度所決定的。當每個碰撞粒子的虛位移狀態被確定之后，應力和應變張量就可以使用2.3節所述的公式進行計算。

圖2 碰撞粒子虛位移狀態的確定，虛位移狀態(黑色)是由碰撞的穿刺深度所決定

3 基于粒子細分的開裂處理

3.1 開裂準則

若將每個SPH粒子的所保存四面體頂點坐標作為內力分析的初始計算位置。當某個點的應力超出材料允許的值時就發生開裂現象。對于脆性破裂的模擬，本文采用經典的Rankine[17]條件：當固體內部某點應力σ的主應力所對應的特征值pmax超過了給定的材料閾值pm:pmax＞pm時，材料達到其破壞條件，裂縫由該點產生，初始開裂斷面的法向為主應力的方向。裂縫的延伸通過引入了一個開裂半徑的參數來模擬。

3.2 開裂斷面的快速生成

第2.1節提出的固體的表面模型可以快速地進行開裂面的更新。當初始裂縫的位置和方向被確定之后，首先在以初始點為圓心，開裂半徑的球體區域內進行查詢，任何與開裂法相所在平面相交的“鏈接”都會被移除。與該鏈接相交的三角形將被加入到一個三角形集合中用來作為渲染的幾何面。圖3是二維下開裂斷面處理的過程。

圖3 二維下破裂計算與處理

使用上述方法生成的開裂面能夠有效地維護幾何拓撲關系并且生成開裂斷面。但是，在模型初始離散的尺度較大的情況下其模擬的斷裂面效果是較為粗糙的，并且對原始計算開裂面的逼近程度不夠，然而如果提高初始模型的離散尺度又會增加系統的計算負擔。為了解決這個問題，這里提出了一種粒子細化的模型來增強開裂細節。首先，在進行開裂處理之前，對所有屬于移除“鏈接”的粒子進行標記。然后對每個標記的粒子，使用一種粒子細分方法將原始粒子進行分解以增強細節。

由于初始的每個粒子實際上是代表著一個四面體，因此，每個粒子細分后的粒子所代表的總四面體區域應該與初始四面體區域相等，粒子的細分實際上是對原有四面體區域進行細分。為了減少細分處理的復雜性，這里選擇了一種較為統一的四面體細分計算模型，將每個需要細分的粒子分解成12個粒子。

首先，對于每個需要細分的粒子，在其所代表的四面體各邊中點，將其分解為4個角的四面體和中心的一個八面體，對中心的八面體需要作進一步分割，一種普遍的方法是在八面體中插入一條對角線，使八面體分割成4個小四面體，如圖4(a)～(b)。由于可選擇的八面體對角線有兩條，在實際細分時需要選擇對角線進行操作，而如何選擇對角線的標準不好確定。為了避免選擇對角線的操作，本文使用一種新的細分操作，見圖4(c)。對于每個八面體單元，在八面體的重心處(八面體6個頂點坐標的平均值)插入一個新的頂點, 然后將八面體原有的6個頂點與新頂點連接, 形成8個四面體，如圖4(d)。

圖4 一個八面體細分成四面體的算法

這樣做的優點是每次劃分一個四面體時只需要計算原四面體六條邊的中點和形心，避免了分解方式的不確定性。圖5是一個四面體分解的過程。首先將粒子所代表的四面體的各邊中心分解為一個4個小四面體和一個八面體，再將八面體的分解為8個四面體，從而將原粒子細分為12個小粒子。在粒子細分完成之后再根據細分后的粒子進行開裂面的生成。圖6是使用粒子細分進行開裂處理的過程。

圖5 單個粒子細分的過程

使用上述粒子細分的模型，僅僅需要在開裂處理時計算細分的粒子，因此對于模擬框架的整體計算效率影響不大，卻又能得到細節較好的結果。

圖6 二維下使用粒子細分進行開裂處理

綜上所述，整個開裂模擬過程可以分成以下4個步驟：

(1) 進行剛體動力學計算和碰撞檢測，發生碰撞時，使用SPH方法對所有粒子的坐標頂點進行應力分析。若某個點主應力所對應的特征值大于給出的材料參數值；將其加入到開裂點集合中。

(2) 對每個開裂點，在用戶定義的開裂半徑的區域內進行搜索，將與開裂斷面相交的所有“鏈接”的粒子進行標記

(3) 對每個粒子進行細分處理，將細分后的粒子進行開裂斷面的生成

(4) 根據粒子間關系來進行碎片的生成。

4 算法實現與結果分析

4.1 算法實現

圖7是整個模擬算法的流程，在每個時間步長內，首先使用剛體動力學來模擬固體的運動。碰撞檢測是基于固體表面粒子的，我們使K-D樹進行碰撞檢測的加速，碰撞反應采用了文獻[18]中基于沖量的方法。分屬于不同固體的表面粒子在發生碰撞時產生虛位移并使用SPH進行內力分析，將與計算開裂法相關的粒子進行細分后進行開裂斷面的生成。

圖7 整體框架算法流程圖

以上方法在Intel (R) Core i3 2100 (R) 3.2 GHz CPU,4 GB 內存, NVIDIA GeForce GTX460 圖形卡,1 GB顯存的臺式機上進行編程。使用Intel? TBB庫進行了在CPU上的并行加速。模擬的結果采用V-Ray(www.chaosgroup.com)進行真實感渲染。

4.2 繪制結果與分析

圖8是模擬9個下落的磚塊在撞擊地面產生的破裂效果，九個磚塊從不同角度掉落地面，從而產生不同的破碎效果。圖9是一個實心球撞擊磚墻的模擬結果，圖9(a)是不采用粒子細分算法和采用粒子細分算法進行開裂面生成的模擬結果，圖9(b)是采用粒子細分算法進行開裂模擬的放大圖，通過對比可以看出在采用粒子細分的算法進行開裂面生成對細節有明顯的提高。而總的計算速度卻沒有大幅度下降，僅僅在進行斷面生成時需要額外的計算時間。

表1是本文算法進行模擬的性能統計。表2是本文方法與相關文獻方法的比較。文獻[8]、[11]、[13]是近年來圖形學領域中進行破碎模擬的幾種典型方法。相對于文獻[8]基于網格的有限元方法，由于采用了將四面體網格幾何信息與形心粒子相結合的方法，本文方法在維護幾何拓撲的效率上有優勢。相對于同樣使用粒子的MLS(文獻[11])及MLPG(文獻[13])方法，由于MLS及MLPG在內力分析時對鄰接粒子位置的要求較高，因此離散時需要較小的粒子，而本文的數值方法是基于SPH，其數值穩定性較高，而且采用了細分粒子的方法，能夠在較小的計算代價下得到更為豐富的細節，既保證計算效率，又不丟失細節。

圖8 磚塊落地破碎的模擬效果

圖9 磚墻受沖擊倒塌的模擬效果(左側：使用原始粒子；右側：使用細分粒子)

表1 不同實驗的性能對比

表2 本文方法與相關文獻方法的比較

5 總結與展望

本文提出了一種基于細分粒子的剛體破裂快速模擬算法。該算法具有如下特點：

(1) 該算法使用粒子固體進行離散建模。通過對每個碰撞粒子估計一個虛位移狀態并使用來SPH方法對碰撞時產生的局部內力進行了分析求解，能有效地提高計算效率。

(2) 將固體四面體化后的網格綁定到粒子上，并提出一種細化粒子的算法來進行開裂面的生成，同時兼顧計算速度和模擬細節。

本文固體模型是基于粒子的，規則固體在進行碰撞時需要的檢測的粒子數量較多，使之成為快速模擬的瓶頸，下一步可以考慮多尺度的粒子框架。

未來的工作包括：采用GPU技術加速計算，將模型擴展到塑形破裂的模擬，考慮更為復雜大型場景的模擬。

[1]Terzopoulos D, Platt J, Barr A, Fleischer K.Elastically deformable models [C]// Proceedings of the 14th Annual Conference on Computer Graphics and Interactive Techniques.Anaheim, California, USA, 1987: 205-214.

[2]Terzopoulos D, Fleischer K.Modeling inelastic deformation: viscolelasticity, plasticity, fracture [C]//Proceedings of the 15th Annual Conference on Computer Graphics, and Interactive Techniques.Atlanta, Georgia,USA, 1988: 269-278.

[3]Norton A, Turk G, Bacon B, Gerth J, Sweeney P.Animation of fracture by physical modeling [J].The Visual Computer, 1991, 7(4): 210-219.

[4]Smith J, Witkin A, Baraff D.Fast and controllable simulation of the shattering of brittle objects [J].Computer Graphics Forum, 2001, 20(2): 81-90.

[5]Fung Y C.A first course in continuum mechanics [M].Prentice-Hall, Englewood Cliffs, N.J, 1969: 10-15.

[6]O'Brien J F, Hodgins J K.Graphical modeling and animation of brittle fracture [C]// Proceedings of the 26th Annual Conference on Computer Graphics and Interactive Techniques.Los Angeles, California, USA,1999: 137-146.

[7]O'Brien J F, Bargteil A W, Hodgins J K.Graphical modeling and animation of ductile fracture [J].ACM Transactions on Graphics, 2002, 21(3): 291-294.

[8]Bao Zhaosheng, Hong J M, Teran J, Fedkiw R.Fracturing rigid materials [J].IEEE Transactions on Visualization and Computer Graphics, 2007, 13(2):370-378.

[9]Müller M, McMillan L, Dorsey J, Jagnow R.Real-time simulation of deformation and fracture of stiff materials [C]// Proceedings of the Eurographic Workshop on Computer Animation and Simulation.Manchester, UK, 2001: 113-124.

[10]Parker E G, O'Brien J F.Real-time deformation and fracture in a game environment [C]// Proceedings of the 2009 ACM SIGGRAPH/Eurographics Symposium on Computer Animation.New Orleans, Louisiana, USA,2009: 165-175.

[11]Pauly M, Keiser R, Adams B, Dutr′e P, Gross M, Guibas L J.Meshless animation of fracturing solids [J].ACM Transactions on Graphics , 2005, 24(3): 957-964.

[12]Liu Ning, He Xiaowei, Li Sheng, Wang Guoping.Meshless simulation of brittle fracture [J].Computer Animation and Virtual Worlds, 2011, 22(2-3): 115-124.

[13]Müller M, Teschner M, Gross M.Physically-based simulation of objects represented by surface meshes [C]// Proceedings of the Computer Graphics International, Hersonissos, Crete, Greece, 2004: 26-33.

[14]Molino N, Bao Zhaosheng, Fedkiw R.A virtual node algorithm for changingmesh topology during simulation [J].ACM Transactions on Graphics, 2004,23(3): 385-392.

[15]Cheng S W, Tamal Krishna D, Shewchuk J R.Delaunay mesh generation [M].Broken Sound Parkway NW: CRC Press, 2012: 410.

[16]Liu M B, Liu G R.Smoothed particle hydrodynamics(SPH): an overview and recent developments [J].Archives of Computational Methods in Engineering,2010, 17(1): 25-76.

[17]Rankine W J M.On the stability of loose earth [J].Philosophical Transactions of the Royal Society of London, 1857, 147: 9-27.

[18]Bender J.Impulse-based dynamic simulation in linear time [J].Computer Animation and Virtual Worlds, 2007,18(4-5): 225-233.