一種sph與動態(tài)表面網(wǎng)格相結合的流體仿真方法【專利摘要】本發(fā)明涉及一種SPH與動態(tài)表面網(wǎng)格相結合的流體仿真方法。該方法使顯式表面追蹤用粒子系統(tǒng)采樣的流體，并在仿真過程中自動處理自交和拓撲變化，維護三角形網(wǎng)格的質量；該方法提出了新的基于顯式表面的自適應重采樣方法。SPH仿真不易展現(xiàn)出很細很薄的特征，尤其是在粒子分辨率不足的情況下，甚至還會形成不自然的空洞。在傳統(tǒng)SPH基礎上如果想體現(xiàn)出更好的細節(jié)，可以增大粒子的分辨率，但是增大分辨率會相應地增加時間和空間上的計算代價，本發(fā)明基于顯式表面的自適應重采樣方法既可以模仿出非常精微的流體表面細節(jié)從而提升仿真的效果，同時又可以提高仿真的效率和速度?！緦＠f明】一種SPH與動態(tài)表面網(wǎng)格相結合的流體仿真方法【
技術領域：
】[0001]本發(fā)明屬于計算機圖形與動畫【
技術領域：
】，涉及一種流體仿真數(shù)據(jù)表示方法以及流體運動和動畫的高效仿真計算方法?！?br>背景技術：
】[0002]在自然界以及日常生活中，流體是很常見的物質形態(tài)。在很多計算機圖形學應用中，流體的仿真也占據(jù)了很重要的地位，例如在動畫、游戲等娛樂領域，以及在科學可視化和醫(yī)療應用領域的應用等?，F(xiàn)在計算機圖形學上的流體仿真與動畫方法有很多種，主流的有基于歐拉觀點的有限差分方法(通稱歐拉法)和基于拉格朗日觀點的粒子法(無網(wǎng)格法，通稱拉格朗日法)。兩種方法各有優(yōu)缺點，各有應用。粒子法在描述流體的細節(jié)特征(波動以及渦旋)方面具有天然的優(yōu)勢。[0003]歐拉法(Euler法)最早由Stam在1999年引入圖形領域(JosStam.Stablefluids[C].Proceedingsofthe26thannualconferenceonComputergraphicsandinteractivetechniques.ACMPress/Addison-ffesleyPublishingC0.,1999,121-128)。在這篇文章中作者描述了使用有限差分法進行流體仿真的算法流程，給出了在Euler框架下求解Navier-stokes方程的方法，并給出了平流和投影的數(shù)值解法。這篇文章提供了在圖形學中進行基于物理的流體仿真的基本框架。[0004]SPH方法(光滑粒子動力學法)是一種基于Lagrange觀點的典型的無網(wǎng)格方法。2003年M--ulIer將該方法引入計算機圖形學中進行流體仿真(MatthiasM"uller,DavidCharypar,MarkusGross.Particle-basedfluidsimulationforinteractiveapplications[C]?Proceedingsofthe2003ACMSIGGRAPH/EurographicssymposiumonComputeranimation.EurographicsAssociation,2003，154-159)。SPH方法避免了基于網(wǎng)格方法中頻繁的重網(wǎng)格化操作，能夠增加程序的魯棒性，降低計算復雜度。SPH方法在連續(xù)介質特別是流體的仿真領域中有廣泛的應用。圖形學中的流體仿真有多種方法，其中比較常用的是基于Euler觀點的有限差分法以及基于Lagrange方法的無網(wǎng)格方法(SPH)。相比于Euler方法，SPH方法能更好地展現(xiàn)流體的細節(jié),例如流體的飛派效果。另外,通過調節(jié)粒子的規(guī)?？梢院苋菀椎氖褂肧PH方法達到實時仿真而不使效果有很大損失。[0005]動態(tài)顯式表面的主要思想是在仿真的初始階段給定一個滿足閉合流形條件的初始顯式表面(通常即為三角形網(wǎng)格)，隨著仿真的不斷進行，每個時間步都得到一個空間中的速度場，由該速度場驅動顯式表面的頂點進行運動并保持連接關系不變。之后得到了一個可能有相交以及退化三角形的網(wǎng)格，對此網(wǎng)格上的三角形進行操作以保證三角形的質量(由面積、邊長比、角度等指標來衡量)，同時修正已經(jīng)相交的部分使網(wǎng)格仍然保持閉合流形性質，在此過程中應能處理由于運動引起的拓撲變化。由于動態(tài)顯式表面的輸入僅需要初始網(wǎng)格和速度場，與底層的仿真無關，只要是能提供速度場的數(shù)值方法都可以使用。又由于動態(tài)顯式表面可以很好地處理形變以及拓撲變化，很適合應用在流體仿真中作為流體的表面。傳統(tǒng)的流體仿真使用隱式表面，即通過仿真模型生成一個距離場，再從距離場上抽取出一個等值面作為流體的表面。顯式表面與隱式表面相比有如下好處:[0006]?顯式表面比隱式表面更有利于保持住表面細節(jié)，這一點上要優(yōu)于動態(tài)隱式表面，后者常常會因為平流中的數(shù)值耗散而抹掉細節(jié)；[0007]?顯式表面能夠很好地表示幾何上退化的特征，例如三維場景下的一維或二維特征。而隱式表面則受限于空間離散化的分辨率大小，往往很難表現(xiàn)這類特征；[0008]?顯式表面能夠很好地保持表面上的參數(shù)化特征，例如紋理坐標等等。使用動態(tài)顯式表面意味著可以省去每步重建表面的過程，不僅可以節(jié)省大量的內存空間(以及通常會節(jié)省計算時間)，更重要的是可以天然地保持住表面上各點在時間序列上的演變關系，從而可以使參數(shù)化特征的傳遞有跡可循；[0009]?對于流體來講，引入顯式表面可以有助于以一種更直接的方式來計算表面物理，例如由于表面張力引起的表面壓差和毛細波等現(xiàn)象。此外，有了顯式表面，可以借助離散微分幾何上的大量知識來理解當前的仿真系統(tǒng)，例如可以便捷地計算表面積以及體積等等；[0010]?隨著仿真的進行，顯式表面有時間維度上的連續(xù)性。由于在每一幀中隱式表面是作為當時時間步的后處理操作來進行的，即隱式表面的形狀只取決于那一時間步生成隱式表面的等值面的形狀。有時在不同幀之間由于粒子分布的不同等值面的變化可能不連續(xù)，這體現(xiàn)為視覺效果上閃動。而顯式表面則是由上一時間步的表面運動而來，保留了上一時間步的狀態(tài)，很大程度上可以避免這種閃動的發(fā)生。[0011]在使用SPH進行仿真計算時，如果粒子系統(tǒng)在局部受到拉伸作用，會產(chǎn)生拉伸不穩(wěn)定(TensileInstability)現(xiàn)象。這種現(xiàn)象的產(chǎn)生是由SPH的計算模型所決定的,表現(xiàn)出來的結果就是粒子系統(tǒng)在比較細、薄的特征處的粒子運動會趨向于不穩(wěn)定，可能徹底分散或是聚成小團。我們希望在粒子變得過于稀疏之前可以對粒子系統(tǒng)進行自適應的重采樣，從而使仿真可以展現(xiàn)出一些細薄特征?！?br/>發(fā)明內容】[0012]由上述內容可知，基于SPH思想的仿真方法不易表示和展現(xiàn)出很細很薄的特征，尤其是在粒子分辨率不足的情況下。在傳統(tǒng)SPH方法基礎上如果想體現(xiàn)出更好的細節(jié)，可以增大粒子的分辨率，但是增大分辨率會相應地增加時間和空間上的計算代價。如果能夠實現(xiàn)自適應的重采樣則可以把有限的計算資源根據(jù)細節(jié)的粗細進行合理地分配。本發(fā)明針對該問題，提出了一種顯式表面與SPH流體相結合的流體仿真方法，可以使顯式表面追蹤用粒子系統(tǒng)采樣的流體，并在仿真過程中自動處理自交和拓撲變化，維護三角形網(wǎng)格的質量。[0013]具體來說，本發(fā)明采用的技術方案如下:[0014]一種SPH與動態(tài)表面網(wǎng)格相結合的流體仿真方法，其步驟包括:[0015]I)根據(jù)腳本在3D空間中對初始流體物質進行采樣以生成粒子系統(tǒng)，計算由粒子定義的隱式函數(shù)V并由等值面穸=0生成初始三角形網(wǎng)格表示流體的外表面形態(tài)，后續(xù)的流體運動仿真過程的計算將分為粒子管線和表面管線兩部分組成；[0016]2)對于粒子管線，首先進行初始粒子設置以表示流體所處的初始狀態(tài)，然后根據(jù)SPH方法計算粒子所受壓強力、粘滯力和表面張力，并移動粒子，其中表面張力由顯式表面上的信息所提供；然后對移動后的粒子系統(tǒng)重新計算隱函數(shù)V并得到其等值面-=Oi[0017]3)對于表面管線，首先初始化流體的顯式表面，該顯式表面為具有閉合流形性質的三角形網(wǎng)格，然后從粒子系統(tǒng)得到顯式表面上每一個頂點的插值速度，并相應地移動頂點至新位置；[0018]4)將步驟3)所得的運動后的網(wǎng)格表面投影到步驟2)重新計算后的等值面0上，在此過程中保證顯式表面仍具有閉合流形性質且無明顯的相交，同時處理有可能產(chǎn)生的拓撲變化事件，然后根據(jù)表面及兩級距離場對流體進行重采樣，從而完成當前時間步之內的流體仿真的所有計算，之后進入下一時間步的迭代，重復執(zhí)行上述步驟2)-4)；[0019]5)當仿真時間滿足一定條件或者仿真狀態(tài)滿足一定條件時，整個流體仿真過程結束。[0020]進一步的，步驟2)中通過粒子系統(tǒng)計算-并得到其等值面?=0,該過程首先定義在三維空間的標量函數(shù)勢對于空間中一點X，計算【權利要求】1.一種SPH與動態(tài)表面網(wǎng)格相結合的流體仿真方法，其步驟包括:1)根據(jù)腳本在3D空間中對初始流體物質進行采樣以生成粒子系統(tǒng)，計算由粒子定義的隱式函數(shù)妒并由等值面黌=O生成初始三角形網(wǎng)格表示流體的外表面形態(tài)，后續(xù)的流體運動仿真過程的計算分為粒子管線和表面管線兩部分；2)對于粒子管線，首先進行初始粒子設置以表示流體所處的初始狀態(tài)，然后根據(jù)SPH方法計算粒子所受壓強力、粘滯力和表面張力，并移動粒子，其中表面張力由顯式表面上的信息所提供；然后對移動后的粒子系統(tǒng)重新計算隱函數(shù)f并得到其等值面-=Oi3)對于表面管線，首先初始化流體的顯式表面，該顯式表面為具有閉合流形性質的三角形網(wǎng)格，然后從粒子系統(tǒng)得到顯式表面上每一個頂點的插值速度，并相應地移動頂點至新位置；4)將步驟3)所得的運動后的網(wǎng)格表面投影到步驟2)重新計算后的等值面-=0上，在此過程中保證顯式表面仍具有閉合流形性質且無明顯的相交，同時處理有可能產(chǎn)生的拓撲變化事件，然后根據(jù)表面及兩級距離場對流體進行重采樣，從而完成當前時間步之內的流體仿真的所有計算，之后進入下一時間步的迭代，重復執(zhí)行上述步驟2)-4)；5)當仿真時間滿足一定條件或者仿真狀態(tài)滿足一定條件時，整個流體仿真過程結束。2.如權利要求1所述的方法，其特征在于:步驟2)通過粒子系統(tǒng)計算-并得到其等值面黌=0,該過程首先定義在三維空間的標量函數(shù)料對于空間中一點X，計算其中i是粒子的下標，W是平滑函數(shù)，r是粒子i的影響域半徑；在一個規(guī)則格子上對史進行采樣，在進行采樣之后，將粒子沒有覆蓋到的格點上的爐都設成一個負值，以保證在有粒子覆蓋的地方妒值是正的，離粒子系統(tǒng)較遠的地方則呈負值，V=O等值面則可以用來定義粒子系統(tǒng)的表面。3.如權利要求1所述的方法，其特征在于:步驟3)從粒子系統(tǒng)得到顯式表面上每一個頂點的插值速度，該過程構造在空間任意一點X處的插值速度公式為4.如權利要求1所述的方法，其特征在于，步驟4)采用的投影算法包括如下步驟:a)對于顯式表面上的頂點V，其以插值速度運動之后空間位置為X，法向為nv，在標量函數(shù)場w上計算樹-r)和斜X—--,.)，其中C=Algfl(樹r是粒子系統(tǒng)的初始設定半徑，通常是b倍的粒子采樣間距，如果妒是精確的到某個閉合曲面的距離場，則使用e=sign{(p{x))?可得到更快的收斂速度；b)如果則在從X到X+env的線段上用二分法迭代查找黌的0值點；c)如果3x,使得錢X)嗔x+env)>0，表明此時頂點v的拓撲發(fā)生變化，此時如果I史(x+env)j<I斜x)|，則將v移至X+env,否則不移動v的位置，另外設置一個拓撲變化標志來提醒后續(xù)步驟的算法需要處理拓撲變化。5.如權利要求1所述的方法，其特征在于:步驟4)在處理顯式表面相交和拓撲事件時，不對顯式表面地整個表面進行重建，而只對發(fā)生拓撲變化或是自交的區(qū)域進行局部重建，具體處理步驟包括:a)在空間中建立一個正則格子；b)在表面附近的格點上計算有向距離場；c)識別發(fā)生拓撲事件以及自交的區(qū)域，標記為重網(wǎng)格化區(qū)域；d)對表及區(qū)域重網(wǎng)格化，并修補與外部區(qū)域網(wǎng)格的接縫。6.如權利要求5所述的方法，其特征在于，所述識別發(fā)生拓撲事件以及自交的區(qū)域的具體步驟包括:a)依次檢查一個單元是否含有復雜邊、復雜面以及完全包含在單元內部的閉合曲面分支，判斷是否滿足復雜單元的條件；b)距離隱式表面爐=O超過一定距離的單元，稱為深單元，從深單元開始向外擴展以包含周圍所有的復雜單元，直到該區(qū)域的邊界單元都是簡單的，這樣既能保證與拓撲變化相關的復雜單元都被標記為重網(wǎng)格化區(qū)域，同時盡可能保住了高分辨率的幾何特征；c)進行擴展的復雜單元的標記，即對MarchingCubes的15個模板中有4個是在一個格子面上允許有兩條三角形邊的情況標記為重網(wǎng)格化區(qū)域。7.如權利要求5所述的方法，其特征在于，所述重網(wǎng)格化標記區(qū)域的步驟是:a)處理標記區(qū)域的邊界，使內部和外部的三角形片段沿區(qū)域邊界分隔開；b)刪除標記區(qū)域內部三角形,使用MarchingCubes模板在區(qū)域內部重新建立三角形；c)處理邊界，將內部外部三角形片段縫合。8.如權利要求1所述的方法，其特征在于，步驟4)進行重采樣的步驟包括:a)檢測采樣不足處的粒子，標記重采樣區(qū)域；b)在標記區(qū)域內進行粒子的快速泊松盤重采樣并增加采樣點；c)重新分配重采樣區(qū)域的粒子質量和動量；d)合并符合質量條件、密度條件以及質量約束這三條合并條件的粒子。9.如權利要求8所述的方法，其特征在于，所述粒子的快速泊松盤重采樣方法的步驟包括:a)在每個格子單元內部，隨機生成t個探測位置，對每生成一個探測點，檢測其是否離現(xiàn)有的粒子以及表面足夠遠，即探測點的位置X必須滿足樹x}>0且I(X-Xi)|>kpdQ且Ix-VjlSkvCltl,其中樹4>0即X必須在流體內部，i和j分別是粒子與表面頂點下標，kp和kv是常數(shù)，d0是仿真開始時粒子的初始間距；b)如探測點滿足以上要求，則將探測點處插入一個新的粒子，并將粒子映射到格子上，參數(shù)kp和kv限制可以插入粒子的區(qū)域，即不僅在流體內部，還要離粒子和表面一定的距離；c)對新產(chǎn)生的粒子進行一步松弛操作以優(yōu)化新粒子的位置，對于每個新產(chǎn)生的粒子Pi,其以r為半徑的鄰域內的粒子下標集合為NPi,鄰域內的表面頂點下標集合為NVi,計算距離Pi與其周圍頂點或粒子的最小距離⑶:.r,.,v；之后，在Pi的鄰域范圍內隨機撒若干探測點，對于每個探測點‘，如果其計算出的Clnlaxnlin大于原來的值，則用ft的位置替換Pi;最后，Pi將被移動到其鄰域內與周圍頂點和粒子距離近似最大的位置。10.如權利要求8所述的方法，其特征在于，重采樣區(qū)域內的粒子的質量和動量的重新分配，其步驟包括:a)將重采樣區(qū)域進行聚類，并按照聚類結果進行分配，聚類的方法是設置一個隊列并將一個重采樣單元放入隊尾，之后進行迭代循環(huán)，將隊首單元彈出，加入集合S，并檢查這個單元的鄰接單元，如果也為重采樣單元則加入隊尾，直到隊列為空，此時S中的單元聚為一類，如所有重采樣單元都已被檢查過，則聚類完成，否則繼續(xù)將一個未檢查的單元放入隊尾進行以上操作；b)在質量的分配上，采用的方法是平分一類內部的粒子質量，設粒子數(shù)為n。，原粒子的下標集為P。，粒子Pi的質量為IV則重采樣后粒子上的采樣質量為【文檔編號】G06T17/00GK103714575SQ201310744592【公開日】2014年4月9日申請日期:2013年12月30日優(yōu)先權日:2013年12月30日【發(fā)明者】李勝,任毅,王衡,汪國平申請人:北京大學