陳沸鑌,謝步瀛,冉修遠(yuǎn)
(同濟(jì)大學(xué)建筑工程系,上海 200092)
日常生活中,脆性破裂是一種很常見(jiàn)的現(xiàn)象,如瓷器的破碎,混凝土磚墻的倒塌。如何對(duì)這種現(xiàn)象進(jìn)行模擬是計(jì)算機(jī)圖形學(xué)領(lǐng)域的一個(gè)重要研究方向,在計(jì)算機(jī)游戲、影視特效、虛擬現(xiàn)實(shí)等領(lǐng)域中有著廣泛的應(yīng)用。
在計(jì)算力學(xué)和材料力學(xué)領(lǐng)域,對(duì)于固體脆性破裂數(shù)值模擬的研究已經(jīng)進(jìn)行了數(shù)十年,研究人員提出了大量的數(shù)值計(jì)算方法。目前大致有兩種數(shù)值離散方法:基于網(wǎng)格的方法和無(wú)網(wǎng)格方法。有限元法是(finite element method,F(xiàn)EM)一種典型的基于網(wǎng)格的方法,它是通過(guò)將連續(xù)體模型剖分成一個(gè)個(gè)網(wǎng)格單元,根據(jù)各個(gè)網(wǎng)格單元的特性來(lái)計(jì)算網(wǎng)格節(jié)點(diǎn)上的屬性。有限元法的主要缺陷是對(duì)網(wǎng)格剖分的要求較高,當(dāng)由外力產(chǎn)生的網(wǎng)格形變太大的時(shí)候,對(duì)控制方程求解的精度會(huì)產(chǎn)生很大的影響。相對(duì)而言,近年來(lái)引起研究人員極大興趣的無(wú)網(wǎng)格(粒子)方法由于在計(jì)算時(shí)不受到背景網(wǎng)格的約束可以有效地解決上述基于網(wǎng)格方法的缺陷。然而,在圖形學(xué)領(lǐng)域,使用基于物理的方法對(duì)脆性破裂進(jìn)行動(dòng)畫生成的研究卻為數(shù)不多?;谖锢淼墓腆w脆性破裂動(dòng)畫主要是以斷裂力學(xué)以及材料力學(xué)為基礎(chǔ),采用數(shù)值離散方法,來(lái)分析固體在外力作用下產(chǎn)生的內(nèi)部應(yīng)力。
與數(shù)值分析方法不同的是,在計(jì)算機(jī)圖形學(xué)中,為了達(dá)到相對(duì)快速的模擬速度,適當(dāng)?shù)哪P秃?jiǎn)化以及對(duì)數(shù)值方法的進(jìn)行優(yōu)化是非常重要的。與此同時(shí),如何增強(qiáng)模擬的細(xì)節(jié)也是引起研究人員廣泛關(guān)注的另一個(gè)問(wèn)題。單方面的追求計(jì)算速度往往導(dǎo)致模擬的結(jié)果過(guò)于粗糙,然而使用分辨率較好的原始模型又會(huì)帶來(lái)計(jì)算代價(jià)較大的問(wèn)題。如何能在提高模擬框架的計(jì)算速度的前提下又能夠給出相對(duì)較好的模擬精度是一個(gè)具有挑戰(zhàn)性的問(wèn)題。針對(duì)這種問(wèn)題所提出的局部細(xì)分技術(shù)能夠很好的地解決這一矛盾,局部細(xì)分技術(shù)的主要思想是以一個(gè)相對(duì)較大的尺度進(jìn)行初始離散建模,然后通過(guò)在需要增強(qiáng)細(xì)節(jié)的部分進(jìn)行局部的細(xì)分方式來(lái)重采樣生成較小尺度的局部模型以此來(lái)提高模擬的細(xì)節(jié)。
本文提出一種基于細(xì)分粒子的剛體脆性破裂的模擬框架。首先將四面體化的固體網(wǎng)格模型綁定到一系列粒子上。然后使用光滑粒子流體動(dòng)力學(xué)(smoothed particle hydrodynamics,SPH)方法來(lái)求解固體碰撞時(shí)產(chǎn)生的內(nèi)部應(yīng)力。通過(guò)更新SPH粒子所綁定網(wǎng)格節(jié)點(diǎn)的狀態(tài),能夠快速的對(duì)開(kāi)裂斷面進(jìn)行維護(hù)從而減少在網(wǎng)格重構(gòu)和拓?fù)涓律系挠?jì)算代價(jià)。最后提出一種細(xì)分粒子的算法,從而在提高系統(tǒng)的計(jì)算效率的同時(shí)又增強(qiáng)了模擬效果的細(xì)節(jié)。
與脆性破裂動(dòng)畫研究有關(guān)的工作主要可以從數(shù)值計(jì)算模型和幾何模型兩個(gè)方面來(lái)進(jìn)行闡述。
從數(shù)值計(jì)算模型的角度來(lái)說(shuō),在圖形學(xué)領(lǐng)域最早進(jìn)行固體破裂動(dòng)畫的研究可以追溯到上世紀(jì)80年代末。Terzopoulos等[1-2]最早使用基于物理的方法來(lái)進(jìn)行彈塑性變形和破裂模擬,他們使用的數(shù)值模型是有限差分方法。同一時(shí)期研究較多的模型還有質(zhì)點(diǎn)-彈簧模型[3]或者質(zhì)點(diǎn)-約束模型[4]。上述模型的主要缺點(diǎn)是不能精確的計(jì)算開(kāi)裂點(diǎn)的位置和方向,以此模擬的結(jié)果并不十分逼真。在破裂模擬中,更為精確的數(shù)值模型是基于連續(xù)介質(zhì)力學(xué)[5]。作為一種典型的基于網(wǎng)格的數(shù)值離散方法,有限元法首次被O'Brien等[6-7]用來(lái)求解線彈性力學(xué)方程進(jìn)行脆性和塑形破裂的模擬,從而得到了較為真實(shí)的效果。在此基礎(chǔ)上,許多研究人員針對(duì)有限元方法進(jìn)行了各種改進(jìn),如對(duì)物理模型的簡(jiǎn)化[8],提高模擬的速度[9]以及加強(qiáng)計(jì)算的精度[10]?;诰W(wǎng)格的數(shù)值方法的缺點(diǎn)在于得到計(jì)算精度高的結(jié)果需要高質(zhì)量的網(wǎng)格,而無(wú)網(wǎng)格方法不存在這種問(wèn)題。Pauly等[11]使用移動(dòng)最小二乘法(moving least square,MLS)來(lái)求解應(yīng)變張量從而進(jìn)行脆性和塑性破裂的模擬。其不足之處在于為了保證計(jì)算的正確,每個(gè)粒子對(duì)其鄰接粒子的位置和數(shù)量有較高的要求。Liu等[12]采用了無(wú)網(wǎng)格局部 Petrov-Galerkin方法(meshless local petrov-Galerkin,MLPG)解決了這個(gè)問(wèn)題,但是在求解的時(shí)候需要計(jì)算一個(gè)大規(guī)模的線性方程組,因此計(jì)算的開(kāi)銷較大。
從幾何模型的角度來(lái)說(shuō),破裂模擬實(shí)際上是對(duì)一個(gè)整體的網(wǎng)格進(jìn)行拓?fù)涞那懈畈僮?,因此,如何選擇一個(gè)合適的固體表面網(wǎng)格模型對(duì)模擬速度以及仿真結(jié)果有著重要的影響。固體的曲面模型大致可以分成2種:顯式的網(wǎng)格模型與隱式的網(wǎng)格模型。顯式模型在模擬過(guò)程中不需要在每個(gè)時(shí)間步長(zhǎng)內(nèi)反復(fù)進(jìn)行表面網(wǎng)格的生成,只需要處理節(jié)點(diǎn)間的關(guān)系,因此處理代價(jià)低、計(jì)算速度高,但由于網(wǎng)格是初始計(jì)算時(shí)就生成好的,因此進(jìn)行幾何切割操作較為復(fù)雜。如Müller等[9,13]使用的四面體網(wǎng)格模型和規(guī)則八面體網(wǎng)格模型、Molino等[14]提出的虛節(jié)點(diǎn)算法等,均是對(duì)初始的網(wǎng)格做了一定的約束條件之后才進(jìn)行切割操作的。隱式曲面模型是使用隱式曲面方程來(lái)生成固體的表面,對(duì)復(fù)雜曲面變化的適用性較好而且不存在切割操作的問(wèn)題,但是由于需要在每個(gè)時(shí)間步長(zhǎng)內(nèi)進(jìn)行網(wǎng)格的計(jì)算,因此計(jì)算代價(jià)較高。如文獻(xiàn)[11]中提出的一種基于點(diǎn)的隱式跟蹤曲面模型,其開(kāi)裂面的生成是通過(guò)構(gòu)造型立體幾何表達(dá)法(constructive solid geometry,CSG)進(jìn)行生成的。
綜上所述,本文脆性破裂模擬的框架提出了一種基于粒子和網(wǎng)格的混合模型。在數(shù)值離散方面使用了基于粒子SPH方法進(jìn)行求解。在幾何方面使用了一種將網(wǎng)格綁定于粒子上幾何模型,能夠快速對(duì)開(kāi)裂固體曲面進(jìn)行維護(hù)操作。
本文的數(shù)值計(jì)算模型是基于粒子的,因此首先需要將初始模型離散成一系列粒子。同時(shí)考慮到開(kāi)裂斷面的生成,一個(gè)合適的固體幾何模型對(duì)于開(kāi)裂時(shí)幾何拓?fù)涞男薷暮湍M結(jié)果的渲染也是極其重要的?;谶@兩點(diǎn)的考慮,本文使用了一種將離散四面體綁定到粒子上的顯式曲面幾何模型。
在數(shù)值計(jì)算中,網(wǎng)格生成是一個(gè)關(guān)鍵性的技術(shù),如何將一個(gè)連續(xù)的計(jì)算區(qū)域離散成一系列的四面體是一個(gè)重要的研究課題,作為數(shù)值計(jì)算和計(jì)算幾何的交叉研究?jī)?nèi)容,四面體網(wǎng)格生成技術(shù)經(jīng)過(guò)近20年的發(fā)展已逐漸成熟并衍生出大量的網(wǎng)格生成技術(shù),在這之中,Delaunay三角剖分由于其網(wǎng)格高質(zhì)量特征和對(duì)復(fù)雜區(qū)域良好的逼近性受到研究人員的廣泛關(guān)注,是目前最流行的全自動(dòng)網(wǎng)格生成方法之一。因此,本文選用Delaunay三角剖分來(lái)進(jìn)行初始的四面體網(wǎng)格生成,詳細(xì)的生成算法可以參見(jiàn)文獻(xiàn)[15]。
首先,將整個(gè)固體區(qū)域通過(guò)Delaunay三角化剖分成一系列的四面體,在每個(gè)四面體的形心位置分配一個(gè)粒子。粒子的半徑r可根據(jù)四面體的體積V進(jìn)行計(jì)算:為了進(jìn)行后續(xù)的拓?fù)渚S護(hù)操作為每個(gè)粒子保存2種信息:原四面體的頂點(diǎn)坐標(biāo)和指向相鄰粒子的指針(圖1中的“鏈接”)。四面體頂點(diǎn)坐標(biāo)(圖1中藍(lán)色的點(diǎn))主要有2個(gè)作用:維護(hù)開(kāi)裂的曲面的生成以及計(jì)算內(nèi)部應(yīng)力。當(dāng)粒子的鄰接粒子數(shù)小于4時(shí)將該粒子定義為固體表面粒子,其作用是進(jìn)行碰撞處理和計(jì)算虛位移狀態(tài)。在初始化時(shí),所有表面粒子所對(duì)應(yīng)的三角形面都將加入到一個(gè)三角形網(wǎng)格中作為固體的曲面模型。 SPH粒子所保存的信息可以快速的進(jìn)行開(kāi)裂面的維護(hù),將在第4節(jié)中具體介紹。
圖1 固體離散建模過(guò)程
本文破裂模擬框架的力學(xué)模型是基于線彈性固體力學(xué)。在線彈性固體力學(xué)中,固體的應(yīng)力與應(yīng)變是與固體內(nèi)部某個(gè)點(diǎn)的變形狀態(tài)有關(guān)的。設(shè)分別為固體中某個(gè)點(diǎn)的初始位置和變形后的位置。由u=[u,v,w]T定義了一個(gè)位置的向量場(chǎng):xt=x0+u。該點(diǎn)的應(yīng)變可以用位移的梯度場(chǎng)?u來(lái)進(jìn)行計(jì)算。令J代表由x0到tx映射的雅克比矩陣。采用線性的柯西-格林應(yīng)變張量,可以將應(yīng)變張量表示為:
如果將固體近似考慮成線彈性各項(xiàng)同性材料,應(yīng)力與應(yīng)變張量的關(guān)系可以定義為:
對(duì)于各項(xiàng)同性材料,C是一個(gè)四階各項(xiàng)同性張量[5]。
其中E與v分別是材料的彈性模量和泊松比。
SPH是一種無(wú)網(wǎng)格粒子法,最早用于天體物理現(xiàn)象的模擬[16],此后廣泛應(yīng)用于流體力學(xué)和固體力學(xué)的數(shù)值分析中。在SPH中,每個(gè)粒子在連續(xù)的標(biāo)量場(chǎng)f(x)的值是通過(guò)核近似法轉(zhuǎn)化為支持域內(nèi)所有離散的粒子數(shù)據(jù)的疊加來(lái)計(jì)算得到:
相關(guān)SPH的詳細(xì)內(nèi)容及公式推導(dǎo)可參見(jiàn)文獻(xiàn)[15]。
通過(guò)使用SPH方法可以對(duì)2.2節(jié)中線彈性固體力學(xué)方程來(lái)進(jìn)行離散求解,在2.1節(jié)中已經(jīng)將計(jì)算區(qū)域離散成與四面體有關(guān)的粒子。因此,對(duì)固體內(nèi)部的每個(gè)粒子,位移梯度場(chǎng)的SPH近似表達(dá)式為:
其中向量uji代表的鄰接粒子uj與當(dāng)前粒子iu的位移差:
采用上述SPH方法可以對(duì)固體在外力作用下所產(chǎn)生的內(nèi)部應(yīng)力進(jìn)行分析。在模擬過(guò)程中,外力主要是由于固體間的相互碰撞產(chǎn)生的。因此首先要進(jìn)行的是碰撞的處理。將剛體離散成一系列粒子之后,剛體間的碰撞檢測(cè)實(shí)際上是由粒子來(lái)完成的。在模擬脆性破裂時(shí),可將固體考慮成具有無(wú)限大剛性的剛體。為了能夠使用2.2節(jié)的線彈性固體力學(xué)模型進(jìn)行分析,需要為發(fā)生碰撞的固體粒子計(jì)算一個(gè)虛位移的狀態(tài)來(lái)得到碰撞時(shí)發(fā)生的變形以此來(lái)計(jì)算應(yīng)變和應(yīng)力。對(duì)每個(gè)發(fā)生碰撞的粒子,本文使用粒子的碰撞信息來(lái)計(jì)算虛位移狀態(tài),如圖2所示,碰撞粒子的虛位移的狀態(tài)是由穿刺深度所決定的。當(dāng)每個(gè)碰撞粒子的虛位移狀態(tài)被確定之后,應(yīng)力和應(yīng)變張量就可以使用2.3節(jié)所述的公式進(jìn)行計(jì)算。
圖2 碰撞粒子虛位移狀態(tài)的確定,虛位移狀態(tài)(黑色)是由碰撞的穿刺深度所決定
若將每個(gè)SPH粒子的所保存四面體頂點(diǎn)坐標(biāo)作為內(nèi)力分析的初始計(jì)算位置。當(dāng)某個(gè)點(diǎn)的應(yīng)力超出材料允許的值時(shí)就發(fā)生開(kāi)裂現(xiàn)象。對(duì)于脆性破裂的模擬,本文采用經(jīng)典的Rankine[17]條件:當(dāng)固體內(nèi)部某點(diǎn)應(yīng)力σ的主應(yīng)力所對(duì)應(yīng)的特征值pmax超過(guò)了給定的材料閾值pm:pmax>pm時(shí),材料達(dá)到其破壞條件,裂縫由該點(diǎn)產(chǎn)生,初始開(kāi)裂斷面的法向?yàn)橹鲬?yīng)力的方向。裂縫的延伸通過(guò)引入了一個(gè)開(kāi)裂半徑的參數(shù)來(lái)模擬。
第2.1節(jié)提出的固體的表面模型可以快速地進(jìn)行開(kāi)裂面的更新。當(dāng)初始裂縫的位置和方向被確定之后,首先在以初始點(diǎn)為圓心,開(kāi)裂半徑的球體區(qū)域內(nèi)進(jìn)行查詢,任何與開(kāi)裂法相所在平面相交的“鏈接”都會(huì)被移除。與該鏈接相交的三角形將被加入到一個(gè)三角形集合中用來(lái)作為渲染的幾何面。圖3是二維下開(kāi)裂斷面處理的過(guò)程。
圖3 二維下破裂計(jì)算與處理
使用上述方法生成的開(kāi)裂面能夠有效地維護(hù)幾何拓?fù)潢P(guān)系并且生成開(kāi)裂斷面。但是,在模型初始離散的尺度較大的情況下其模擬的斷裂面效果是較為粗糙的,并且對(duì)原始計(jì)算開(kāi)裂面的逼近程度不夠,然而如果提高初始模型的離散尺度又會(huì)增加系統(tǒng)的計(jì)算負(fù)擔(dān)。為了解決這個(gè)問(wèn)題,這里提出了一種粒子細(xì)化的模型來(lái)增強(qiáng)開(kāi)裂細(xì)節(jié)。首先,在進(jìn)行開(kāi)裂處理之前,對(duì)所有屬于移除“鏈接”的粒子進(jìn)行標(biāo)記。然后對(duì)每個(gè)標(biāo)記的粒子,使用一種粒子細(xì)分方法將原始粒子進(jìn)行分解以增強(qiáng)細(xì)節(jié)。
由于初始的每個(gè)粒子實(shí)際上是代表著一個(gè)四面體,因此,每個(gè)粒子細(xì)分后的粒子所代表的總四面體區(qū)域應(yīng)該與初始四面體區(qū)域相等,粒子的細(xì)分實(shí)際上是對(duì)原有四面體區(qū)域進(jìn)行細(xì)分。為了減少細(xì)分處理的復(fù)雜性,這里選擇了一種較為統(tǒng)一的四面體細(xì)分計(jì)算模型,將每個(gè)需要細(xì)分的粒子分解成12個(gè)粒子。
首先,對(duì)于每個(gè)需要細(xì)分的粒子,在其所代表的四面體各邊中點(diǎn),將其分解為4個(gè)角的四面體和中心的一個(gè)八面體,對(duì)中心的八面體需要作進(jìn)一步分割,一種普遍的方法是在八面體中插入一條對(duì)角線,使八面體分割成4個(gè)小四面體,如圖4(a)~(b)。由于可選擇的八面體對(duì)角線有兩條,在實(shí)際細(xì)分時(shí)需要選擇對(duì)角線進(jìn)行操作,而如何選擇對(duì)角線的標(biāo)準(zhǔn)不好確定。為了避免選擇對(duì)角線的操作,本文使用一種新的細(xì)分操作,見(jiàn)圖4(c)。對(duì)于每個(gè)八面體單元,在八面體的重心處(八面體6個(gè)頂點(diǎn)坐標(biāo)的平均值)插入一個(gè)新的頂點(diǎn), 然后將八面體原有的6個(gè)頂點(diǎn)與新頂點(diǎn)連接, 形成8個(gè)四面體,如圖4(d)。
圖4 一個(gè)八面體細(xì)分成四面體的算法
這樣做的優(yōu)點(diǎn)是每次劃分一個(gè)四面體時(shí)只需要計(jì)算原四面體六條邊的中點(diǎn)和形心,避免了分解方式的不確定性。圖5是一個(gè)四面體分解的過(guò)程。首先將粒子所代表的四面體的各邊中心分解為一個(gè)4個(gè)小四面體和一個(gè)八面體,再將八面體的分解為8個(gè)四面體,從而將原粒子細(xì)分為12個(gè)小粒子。在粒子細(xì)分完成之后再根據(jù)細(xì)分后的粒子進(jìn)行開(kāi)裂面的生成。圖6是使用粒子細(xì)分進(jìn)行開(kāi)裂處理的過(guò)程。
圖5 單個(gè)粒子細(xì)分的過(guò)程
使用上述粒子細(xì)分的模型,僅僅需要在開(kāi)裂處理時(shí)計(jì)算細(xì)分的粒子,因此對(duì)于模擬框架的整體計(jì)算效率影響不大,卻又能得到細(xì)節(jié)較好的結(jié)果。
圖6 二維下使用粒子細(xì)分進(jìn)行開(kāi)裂處理
綜上所述,整個(gè)開(kāi)裂模擬過(guò)程可以分成以下4個(gè)步驟:
(1) 進(jìn)行剛體動(dòng)力學(xué)計(jì)算和碰撞檢測(cè),發(fā)生碰撞時(shí),使用SPH方法對(duì)所有粒子的坐標(biāo)頂點(diǎn)進(jìn)行應(yīng)力分析。若某個(gè)點(diǎn)主應(yīng)力所對(duì)應(yīng)的特征值大于給出的材料參數(shù)值;將其加入到開(kāi)裂點(diǎn)集合中。
(2) 對(duì)每個(gè)開(kāi)裂點(diǎn),在用戶定義的開(kāi)裂半徑的區(qū)域內(nèi)進(jìn)行搜索,將與開(kāi)裂斷面相交的所有“鏈接”的粒子進(jìn)行標(biāo)記
(3) 對(duì)每個(gè)粒子進(jìn)行細(xì)分處理,將細(xì)分后的粒子進(jìn)行開(kāi)裂斷面的生成
(4) 根據(jù)粒子間關(guān)系來(lái)進(jìn)行碎片的生成。
圖7是整個(gè)模擬算法的流程,在每個(gè)時(shí)間步長(zhǎng)內(nèi),首先使用剛體動(dòng)力學(xué)來(lái)模擬固體的運(yùn)動(dòng)。碰撞檢測(cè)是基于固體表面粒子的,我們使K-D樹(shù)進(jìn)行碰撞檢測(cè)的加速,碰撞反應(yīng)采用了文獻(xiàn)[18]中基于沖量的方法。分屬于不同固體的表面粒子在發(fā)生碰撞時(shí)產(chǎn)生虛位移并使用SPH進(jìn)行內(nèi)力分析,將與計(jì)算開(kāi)裂法相關(guān)的粒子進(jìn)行細(xì)分后進(jìn)行開(kāi)裂斷面的生成。
圖7 整體框架算法流程圖
以上方法在Intel (R) Core i3 2100 (R) 3.2 GHz CPU,4 GB 內(nèi)存, NVIDIA GeForce GTX460 圖形卡,1 GB顯存的臺(tái)式機(jī)上進(jìn)行編程。使用Intel? TBB庫(kù)進(jìn)行了在CPU上的并行加速。模擬的結(jié)果采用V-Ray(www.chaosgroup.com)進(jìn)行真實(shí)感渲染。
圖8是模擬9個(gè)下落的磚塊在撞擊地面產(chǎn)生的破裂效果,九個(gè)磚塊從不同角度掉落地面,從而產(chǎn)生不同的破碎效果。圖9是一個(gè)實(shí)心球撞擊磚墻的模擬結(jié)果,圖9(a)是不采用粒子細(xì)分算法和采用粒子細(xì)分算法進(jìn)行開(kāi)裂面生成的模擬結(jié)果,圖9(b)是采用粒子細(xì)分算法進(jìn)行開(kāi)裂模擬的放大圖,通過(guò)對(duì)比可以看出在采用粒子細(xì)分的算法進(jìn)行開(kāi)裂面生成對(duì)細(xì)節(jié)有明顯的提高。而總的計(jì)算速度卻沒(méi)有大幅度下降,僅僅在進(jìn)行斷面生成時(shí)需要額外的計(jì)算時(shí)間。
表1是本文算法進(jìn)行模擬的性能統(tǒng)計(jì)。表2是本文方法與相關(guān)文獻(xiàn)方法的比較。文獻(xiàn)[8]、[11]、[13]是近年來(lái)圖形學(xué)領(lǐng)域中進(jìn)行破碎模擬的幾種典型方法。相對(duì)于文獻(xiàn)[8]基于網(wǎng)格的有限元方法,由于采用了將四面體網(wǎng)格幾何信息與形心粒子相結(jié)合的方法,本文方法在維護(hù)幾何拓?fù)涞男噬嫌袃?yōu)勢(shì)。相對(duì)于同樣使用粒子的MLS(文獻(xiàn)[11])及MLPG(文獻(xiàn)[13])方法,由于MLS及MLPG在內(nèi)力分析時(shí)對(duì)鄰接粒子位置的要求較高,因此離散時(shí)需要較小的粒子,而本文的數(shù)值方法是基于SPH,其數(shù)值穩(wěn)定性較高,而且采用了細(xì)分粒子的方法,能夠在較小的計(jì)算代價(jià)下得到更為豐富的細(xì)節(jié),既保證計(jì)算效率,又不丟失細(xì)節(jié)。
圖8 磚塊落地破碎的模擬效果
圖9 磚墻受沖擊倒塌的模擬效果(左側(cè):使用原始粒子;右側(cè):使用細(xì)分粒子)
表1 不同實(shí)驗(yàn)的性能對(duì)比
表2 本文方法與相關(guān)文獻(xiàn)方法的比較
本文提出了一種基于細(xì)分粒子的剛體破裂快速模擬算法。該算法具有如下特點(diǎn):
(1) 該算法使用粒子固體進(jìn)行離散建模。通過(guò)對(duì)每個(gè)碰撞粒子估計(jì)一個(gè)虛位移狀態(tài)并使用來(lái)SPH方法對(duì)碰撞時(shí)產(chǎn)生的局部?jī)?nèi)力進(jìn)行了分析求解,能有效地提高計(jì)算效率。
(2) 將固體四面體化后的網(wǎng)格綁定到粒子上,并提出一種細(xì)化粒子的算法來(lái)進(jìn)行開(kāi)裂面的生成,同時(shí)兼顧計(jì)算速度和模擬細(xì)節(jié)。
本文固體模型是基于粒子的,規(guī)則固體在進(jìn)行碰撞時(shí)需要的檢測(cè)的粒子數(shù)量較多,使之成為快速模擬的瓶頸,下一步可以考慮多尺度的粒子框架。
未來(lái)的工作包括:采用GPU技術(shù)加速計(jì)算,將模型擴(kuò)展到塑形破裂的模擬,考慮更為復(fù)雜大型場(chǎng)景的模擬。
[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.