可控性的雙層粒子剛體脆性破裂模擬動畫

陳沸鑌， 王長波， 謝步瀛

(1. 同濟大學建筑工程系，上海 200092；2. 華東師范大學軟件學院，上海 200062)

可控性的雙層粒子剛體脆性破裂模擬動畫

陳沸鑌1， 王長波2， 謝步瀛1

(1. 同濟大學建筑工程系，上海 200092；2. 華東師范大學軟件學院，上海 200062)

提出一種可控性的雙層粒子模型的剛體脆性破裂動畫的模擬框架。首先，使用接觸力學分析撞擊時固體內部產生的彈性位移并采用光滑粒子流體動力學進行離散求解，在此基礎上建立一個可控性的內力分析模型。其次，提出一種雙層粒子模型來對剛體進行建模，在保證內力計算精度的同時提高了剛體動力學中碰撞檢測的效率。最后實現了多個剛體脆性模型場景的動畫。該算法可適用于剛體脆性破裂模擬的動畫應用。

可控性模型；雙層粒子；剛體脆性破裂模擬；光滑粒子流體動力學

固體的脆性破裂是一種常見的自然現象，如摔碎的瓷器、玻璃的碎裂等。表現這些現象并進行真實感繪制，在游戲動畫、影視特效等領域有著重要的應用。在計算機圖形學領域，脆性破裂的模擬已經引起了許多研究人員的注意，并做了許多工作。然而，現有的脆性破裂模擬方法仍然有值得改進的地方。首先，對于脆性破裂的力學分析往往直接采用了連續介質力學的方法，計算的結果雖然有較高的精度然而模擬的結果比較單一，缺乏對動畫模擬的可控性；其次，在大型場景的模擬時，大量的剛體運動和碰撞受力計算對于模擬系統會造成比較大的計算量，如何提高系統的運行效率也是值得改進的一個方面。

針對固體脆性破裂的特點，本文提出了一種具有可控性的雙層粒子的剛體脆性破裂模擬算法，能夠有效地對模擬結果進行控制并提高計算的效率。該方法能有效處理剛體受到撞擊后發生的破裂現象，使其有效地應用于動畫以及影視特效模擬中。

1 相關工作

對于脆性破裂的模擬最早可以追溯到上世紀80年代。Terzopoulos等[1-2]將固體表示成網格，將開裂的模擬假定為彈性受力問題并使用了有限差分法進行計算，類似的研究工作還有Norton等[3]。這些早期的固體脆性破裂模型的共同點是采用了彈簧-質點模型對固體進行建模，將固體的質量集中在各個質點上，并使用彈簧來模擬質點間的受力。這類模型簡單且計算速度快，但計算開裂點時其精度較低，無法準確地確定裂紋尖端，對于各項異性材料的模擬也有較大的難度。為了克服彈簧-質點模型的缺點，研究人員從力學上借鑒了很多更為精確的模型。在工程領域作為一種基于網格的方法，有限元(finite element method，FEM)在模擬變形和破裂方面已經被證明是一種具有較高精度的數值模型。因此在圖形學領域，有限元方法最早由O＜ Brien等人所引進用來模擬脆性破裂[4]和塑性破裂[5]，是將固體離散成一系列的節點和單元，由外部撞擊產生的內力通過對一個由所有節點的屬性組成的大型線性方程組進行求解計算而得到的。這種直接求解線性方程組的計算量是比較高的，而在圖形學領域中的模擬往往不需要如此高的精度，因此可以進行一定程度地簡化，如將撞擊問題假定為準靜態問題[6-7]。

在力學及工程領域，相對于基于網格的力學模型，近年來人們更關注的是無網格方法，即穩定且精度較高的數值方法。無網格方法使用粒子對求解區域進行離散建模，其求解不受網格的約束因而有較大的靈活性，對于開裂這一類求解區域不連續的問題有著較好地應用。在圖形學領域，最早使用無網格方法進行固體模擬的是Desbrun和Cani[8]。光滑 粒 子 流 體 動 力 學 (smoothed particle hydrodynamics，SPH)作為一種無網格方法，已經應用在流體模擬[9]，顆粒介質模擬[10]和固體變形模擬[11]方面。在固體破碎模擬方面，目前采用無網格方法進行破碎模擬的有無網格局部 petrov-galerkin法(meshless local petrov-galerkin，MLPG)[12]和移動最小二乘法(moving least squares，MLS)等方法[13]，而SPH模型尚未應用到固體破碎模擬上。

從幾何學角度看，開裂的模擬實際上是一個網格的分解，這涉及到對一系列三角形或者四面體進行切割，即一個曲面模型的分解與合并的計算，在計算代價上，只在網格的節點上進行分解操作比直接切割網格的計算效率要高，因此許多模型均采用在節點上進行分解操作，如四面體網格模型[6]、規則八面體網格模型[14]、Molino等[15]提出的虛節點算法等，均是在對原始的幾何模型采用Voronoi圖的網格剖分后進行了一定的約束限制再進行網格的分解操作。最近，Müller等[16]提出了一種凸包的近似分解算法(volumetric approximate convex decomposion，VACD)可以在實時的程度上進行三角形的切割操作。總而言之，一種有效的幾何模型在破碎模擬上有著很重要地應用。

在計算機圖形學領域，對于一個模擬算法而言，如何在模擬規模比較大的時候提高計算效率是一個很重要的問題。對于一個粒子系統而言，在求解問題規模較大的時候效率往往是系統的一個瓶頸。針對這個問題，近些年研究人員也提出了許多方法，其中包括雙層/多層粒子建模技術。雙層/多層粒子的核心思想是在一個框架中通過采用不同尺度粒子同時對同一區域進行離散建模，從而產生不同層次的粒子，各層次的粒子在框架中的作用各不相同，該模型已經應用到固體形變模擬[17]、流體模擬[18]、顆粒物質(沙子)模擬[19]等，而在固體破裂模擬方面卻很少見到，從雙層/多層粒子模型的特性來看，其應用于固體脆性模擬具有很大的優勢。

2 固體脆性模擬框架

2.1 彈性接觸模型

為了使用基于物理的方式進行固體脆性模擬，必須建立一種合適的力學模型來分析撞擊產生的應變和應力，與以往采用的基于連續介質力學模型不同，為了提高模擬結果的可控性，本文采用了基于接觸力學的物理模型[20]。假定兩個物體之間產生了碰撞，首先將碰撞力分解為兩個方向的力：與接觸面所垂直的法向力P和與接觸面平行的切向力Q。以碰撞點O為坐標中心，建立一個坐標軸Oxyz，其中Oz的方向是與接觸面垂直的法向。根據接觸力學[20]中的勢能理論(由Love提出)，由于碰撞力P在碰撞點O產生彈性位移可以表示為：

同時，將切向力在碰撞點產生的彈性位移tU在 Ox和Oy軸上進行分解，在Ox軸方向上的彈性位移可以表示為：

式(1)～(3)給出了總彈性位移U，總接觸力P Q+ 以及固體內部任意一點到接觸點的距離ρ之間的關系，從上述公式中可以看出，在物體之間產生撞擊時，其接觸力越大，產生的彈性位移越大，離接觸點的距離越遠，則彈性位移越小。使用上述公式能在碰撞時對接觸域內任意點的位移進行求解計算，從而得到整個物體內部的位移場，在得到位移場后，根據green-saint-venant應變ε的定義，可以得到：

其中，I是單位矩陣，TJ u=? ， u? 是位移的梯度場，根據線彈性材料的虎克定律[21]，應力σ可表示為：

其中，C是一個由材料特性決定的6階方陣。

2.2 光滑粒子流體動力學的數值離散

SPH是一種無網格數值方法[22]，其核心是用一系列分散的粒子對問題區域中的某一個物理量進行求解，如位移、速度、應力等。求解區域中任意一個點的物理量值是由其支持域中所有鄰接粒子所對應的值來進行加權求和得到。

在SPH中，任意連續的物理場〈 〉f 在某個位置i的物理量值〈 fi〉可以由其支持半徑 h所定義的支持域中所有鄰接粒子對應的屬性值使用核函數 W( xij,h)來進行加權求和得到：

其中， xij是該求解點與鄰接粒子之間的距離。核函數 W( xij,h)通常取為具有歸一、對稱、非負的光滑函數?？梢岳蒙鲜綄ι瞎澋奈灰铺荻葓?? u進行離散求解：

其中， mj和 ρj分別是粒子j質量和密度，uji= uj- ui，ui和 uj分別是由式(3)計算得到的粒子i和j位置處的位移。

2.3 可控性的內力分析模型

在動畫和電影特效的制作過程中，制作人員往往需要根據所繪場景的需要對模擬結果進行調整，如控制模擬物體運動的速度，調整物體運動的位移等。在破裂模擬中，模擬模型是否能對破碎的結果進行一定控制，并使模擬的結果具有可控性，是其應用價值所在。合適的模型動畫制作人員通過幾個參數就可以對模擬的結果進行人為干預，并得到預期的結果。針對這個問題，本節提出一種具有可控性的內力分析模型。

為了模擬固體的開裂，需要選擇一定的材料屈服準則來確定開裂點以及開裂法向。對于脆性模擬，可以選擇經典的Rankine準則[23]：當固體內部某個點的主應力超出了某個指定的閾值時，材料達到屈服，發生開裂，其開裂面的法向為主應力的方向。

計算主應力的方向和大小時，先根據2.1節中的公式通過接觸力大小(這里的接觸力是由剛體動力學計算得到的)來計算位移場，再使用2.2節中的SPH式(7)計算位移的梯度場，最后根據式(4)和(5)可以得到應力的大小和方向。

本文使用 SPH數值離散方法來計算位移的梯度場。SPH方法的一個重要特性是某個位置的值是由其鄰接粒子相應的值所決定的。從式(7)可以看出，某個位置位移的梯度是通過該點的鄰接粒子位移值的加權求和計算的。如能選擇性地使用鄰接粒子的位移來計算該點的位移梯度值，所得的應力主方向必然是不相同的，在模擬的結果上表現為同一點產生多條不同的裂縫。盡管從數值方法的角度上說此方法無法保證精度，然而從圖形學的角度來看，卻能達到人為干預裂縫生成的要求。圖1是內力分析的過程，具體的內力分析流程如下：

對某個應力分析點：

獲得在SPH光滑半徑內的鄰接粒子對所有的鄰接粒子：

計算與分析點間的距離

根據人為設定的裂紋數量n進行分組對于每個鄰接粒子組：

使用式(3)～(5)計算應力大小得到n個應力

對于計算的每個應力：

計算應力張量的大小maxd 及方向

如果： max threshold

d d＞

材料達到屈服極限，產生裂紋

圖1 內力分析示意圖

3 雙層粒子模型

如圖2所示，本文提出了一種雙層粒子的模型來對固體的運動過程進行模擬，通過將其離散成兩個不同層次的粒子(物理層粒子邊界層粒子)來分別進行內力分析和動力學計算。從而在保證模擬細節的同時提高了模擬的效率。

圖2 雙層粒子模型示意圖

3.1 物理層粒子

力學分析模型是采用 SPH方法進行內力分析計算的，因此，需將原始的固體模型離散成一系列的粒子。使用文獻[23-24]中的方法對固體進行離散建模，首先將固體離散成一系列的四面體，再根據四面體生成相應的粒子。在四面體離散算法上，采用了 Delaunay網格[25]計算對初始的固體模型進行四面體網格化。具體模型的粒子化過程可見圖3。

圖3 物理層粒子生成過程示意圖

初始的固體模型通過 Delaunay三角化算法得到一個整體的四面體網格；在每個四面體的形心分配一個粒子，其質量等于四面體的質量，粒子的半徑 r是由相應四面體體積 V計算得到的：為了避免開裂處理時網格維護幾何拓撲的計算開銷過大，使用了文獻[24]中基于約束的方法進行幾何網格的切割操作，從而能夠有效地減少計算代價。

在剛體動力學計算時，直接采用上述物理層的粒子進行碰撞檢測是一個簡單而有效的方式。然而，為了保證內力分析時的計算正確性，物理層的粒子尺度往往需要比較小，在碰撞檢測過程中會增加計算代價。當模擬的場景規模比較大時，大量的粒子碰撞檢測會降低模擬計算的速度。為了解決這個問題，在物理層粒子的基礎上額外建立一個邊界層的粒子用來進行剛體動力學檢測的運算。

3.2 邊界層粒子

將物理層的粒子分為兩類：內部粒子和表面粒子，內部粒子是位于固體內部，而表面粒子是位于固體表面。在碰撞檢測時，可使用表面粒子進行碰撞檢測計算。但是直接采用表面粒子進行碰撞檢測計算會受制于物理層粒子的離散尺寸。為了提高檢測效率，在固體表面的物理層粒子的基礎上，需建立一個邊界層的粒子用來進行碰撞檢測。

邊界層粒子的計算過程見圖4。一個由用戶定義的采樣數量n被引入，本文實驗將采樣數量定為4。在所有的邊界粒子中，將每n個物理層表面粒子合并為一個邊界層粒子，每個邊界層粒子的質量與其所對應的物理層表面粒子的總質量相等，其位置為n個物理表面層粒子。在邊界層粒子的生成過程中，可以發現有些物理層表面粒子(如位于固體表面比較尖銳處的粒子)無法被合并成一個邊界層粒子，這個問題可以有兩種解決方法：直接使用這些粒子作為邊界層粒子，或者將這些粒子忽略。

當使用邊界層粒子進行碰撞檢測時，在動力學計算的每個時刻，所有邊界層粒子的位置都會隨著剛體的運動而改變，因此便于下一步碰撞檢測。只有在進行內力分析計算時，才會更新物理層粒子的位置來進行內力的分析計算，通過以上方法可以有效地提高計算的效率。

圖4 邊界層粒子采樣示意圖

4 算法的實現與結果分析

根據用上述章節提出的模擬框架，在Intel (R) Core i3-2100 (R) 3.2 GHz CPU、4 GB內存、NVIDIA GeForce GTX460 圖形卡、1 GB顯存的微機上進行模擬系統的算法編制，實現了多個剛體脆性模擬的場景，并采用3dmaxV-Ray對模擬結果進行離線渲染，以期達到真實感繪制的效果。固體模型的四面體 離 散 化 采 用 了 NETGEN 工 具(sourceforge.net/projects/netgen-mesher)進行生成。剛體動力學部分的算法在第三方物理引擎“Newton Game Dynamics” (newtondynamics.com)的基礎上進行了一定的改進從而實現。

圖5是一個圓環形瓷器下落摔碎的模擬場景，圖 5(b)～(d)是使用本文內力分析模型進行開裂控制的結果，通過使用不同的鄰接粒子組在同一個碰撞點進行開裂法向的計算可以得到不同的結果。

圖6是一個較為復雜的場景的模擬，同時有13個圓環瓷器下落并破碎，對于每一個圓環，在內力分析時均采用了不同參數對裂紋生成進行控制。圖7是3個玻璃墻受到撞擊后連續破碎的效果。表1是圖5～7不同場景模擬的復雜度統計表。

表1 不同場景模擬的復雜度統計

表2是本文方法在不同場景的物理計算時間的統計，由于使用了雙層粒子的模型，在剛體動力學以及碰撞檢測上的時間開銷有了很大降低。

圖5 圓環摔碎的模擬場景通過不同鄰接粒子組來進行開裂法向的計算

圖6 多個圓環摔碎的模擬場景

圖7 3個玻璃墻被撞碎的模擬場景

表2 不同場景的物理計算時間的統計(ms)

5 總結與展望

本文提出了一種剛體脆性破碎的模擬算法。該算法具有如下特點：

(1) 該算法在接觸力學和SPH數值方法內力分析的基礎上進行了改進，提出了一種具有可控性的內力計算模型，從而能夠達到對模擬結果進行用戶交互控制的目的

(2) 針對粒子系統模擬效率的不足，提出了一種雙層粒子的模型對固體進行建模，在保證模擬的精度上又提高了系統運行的效率。

本文的提出模擬框架還有一定的局限性，首先，無法模擬一些復雜的裂紋，如環形的裂紋，需要在現有的內力分析模型上做進一步地改進。其次，雙層粒子的采樣算法需要進行進一步的優化使得能夠支持復雜的幾何模型。

未來的工作包括：采用 GPU技術應用到文中的算法上，將算法擴展到塑性破裂模擬上，進行更為復雜場景的模擬。

[1]Terzopoulos D, Platt J, Barr A, et al. 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]// Proceedingsof the 15th Annual Conference on Computer Graphics and Interactive Techniques, Atlanta, Georgia, USA, 1988： 269-278.

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

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

[5]O＇ Brien J F, Bargteil A, Hodgins J K. Graphical modeling and animation of ductile fracture [J]. ACM Transactions on Graphics, 200, 21(3)： 291-294.

[6]Müller M, McMillan L, Dorsey J, et al. 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.

[7]Bao Zhaosheng, Hong J M, Teran J, et al. Fracturing rigid materials [J]. IEEE Transactions on Visualization and Computer Graphics, 2007, 13(2)： 370-378.

[8]Desbrun M, Cani M P. Space-time adaptive simulation of highly deformable substances [R]. Rapport de Recherche RR-3829, INRIA (1999).

[9]Müller M, Charypar D, Gross M. Particle-based fluid simulation for interactive applications [C]//Proceedings of the 2003 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, San Diego, California, USA, 2003： 154-159.

[10]Alduan I, Otaduy M A. SPH granular flow with friction and cohesion [C]//Proceedings of the 2011 ACM SIGGRAPH/Eurographics Symposium on Computer Animation, Vancouver, British Columbia, Canada, 2011：25-32.

[11]Becker M, Ihmsen M, Teschner M. Corotated SPH for deformable solids [C]//Proceedings of the Fifth Eurographics conference on Natural Phenomena, Munich, Germany, 2009： 27-34.

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

[13]Pauly M, Keiser R, Adams B, et al. Meshless animation of fracturing solids [J]. ACM Transactions on Graphics, 2005, 24(3)： 957-964.

[14]Müller M, Teschner M, Gross M H. Physically-based simulationof objects represented by surface meshes [C]// Proceedingsof the Computer Graphics International, Hersonissos, Crete, Greece, 2004： 26-33.

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

[16]Müller M, Chentanez N, Kim T Y. Real time dynamic fracture with volumetric approximate convex decompositions [J]. ACM Transactions on Graphics, 2013, 32(4)： 115：1-115：10.

[17]Adams B, Pauly M, Keiser R, et al. Adaptively sampled particle fluids [J]. ACM Transactios on Graphics, 2007, 26(3)： 115-120.

[18]Solenthaler B, Gross M. Two-scale particle simulation [J]. ACM Transactios on Graphics, 2011, 30(4)： 81：1-81：8.

[19]Ihmsen M, Wahl A, Teschner M. High-resolution simulation of granular material with SPH [C]// Proceedings of Workshop on Virtual Reality Interaction and Physical Simulation, Darmstadt, Germany, 2012：53-60.

[20]Johnson K. Contact mechanics [M]. Cambirdge：Cambridge University Press, 1985： 108-112.

[21]Fung Y. A first course in continuum mechanics [M]. Englewood Cliffs, N.J： Prentice-Hall, 1969： 75-95.

[22]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.

[23]Gross D, Seelig T. Fracture mechanics： with an introduction to micromechanics [M]. Springer Verlag, 2006： 42-52.

[24]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.

[25]Chen Siuwing, Dey T K, Shewchuk J. Delaunay mesh generation [M]. Boca Raton, Florida： CRC Press, 2012：1-8.

Controllable Two Scale Particles for Fracture Animation of Brittle Materials

Chen Feibin1, Wang Changbo2, Xie Buying1
(1. Department of Architectural Engineering, Tongji University, Shanghai 200092, China; 2. Software Engineering Institute, East China Normal University, Shanghai 200062, China)

A controllable two scale particle framework for simulate the brittle fracture is presented. Firstly, the elastic displacement induced by the collision between two rigid bodies is analyzed based on contact mechanics and solved by using of smoothed particle hydrodynamics, and a stress analysis model is proposed to control the result of simulation. Secondly, a two scale particle model is introduced to present the rigid; this two scale particle system serves for the purposes of more accurate fracture modeling and more effective collision handling. Finally, rendering the fracture of brittle material is achieved in some scenes. Our method can apply for simulating the fracture of brittle material.

controllable model; two scale particles; simulation of brittle fracture; smoothed particle hydrodynamics

TP 391.1

A

2095-302X(2015)01-0111-06

2014-08-11；定稿日期：2014-08-20

國家自然科學基金面上資助項目(61272199)；教育部博士點基金資助項目(20130076110008)；上海市教委科研創新項目重點資助項目(12ZZ042)

陳沸鑌(1983-)，男，福建南平人，博士研究生。主要研究方向為基于物理的仿真、流體模擬。E-mail：feibin.chen@gmail.com

王長波(1976-)，男，湖北隨州人，教授，博士。主要研究方向為計算機圖形學、虛擬現實等。E-mail：cbwangcg@gmail.com