基于SPH方法的快速逼真流體表面張力仿真

基于SPH方法的快速逼真流體表面張力仿真摘要隨著計算機仿真技術的發(fā)展，流體的模擬仿真被廣泛地應用于工程計算、游戲及電影特效中.逼真的流體仿真需要對流體的物理特性進行建模，基于物理的流體仿真是仿真領域最具挑戰(zhàn)的研究內容之一.基于光滑粒子流體動力學（SPH）的流體仿真方法是目前最為常用的流體仿真無網格法，其粒子特性特別適合處理流體的建模問題，同時能夠魯棒穩(wěn)定地解決流體仿真中的計算問題.但是,SPH方法無法實現(xiàn)精確的表面張力計算，無法獲得逼真的流體仿真效果.SPH表面計算中的數(shù)值不穩(wěn)定問題，引起的表面張力不真實、自由表面處粒子聚集、細節(jié)斷裂等現(xiàn)象，嚴重影響流體自由表面處仿真的視覺效果.為解決該問題，提出一種面向SPH流體的快速逼真流體表面張力仿真方法.首先快速精確找到流體表面粒子，精準計算其表面粒子法向量，其次通過構建水平集網格用于邊界處密度修正以更加準確地處理邊界處的計算，最后提出近表面粒子間作用力模型以計算流體表面張力.實驗結果表明該方法能實現(xiàn)逼真的流體表面張力仿真效果.關鍵詞光滑粒子流體動力學；流體仿真；表面張力；近表面粒子間作用力；密度修正FastandPreciseSurfaceTensionforSPH-BasedFluidSimulationAbstractWiththecontinuousprogressofcomputersimulationtechnology，fluidsimulationhasbeendevelopedrapidly，andithasbeenwidelyusedinmanyfields，especiallyinengineeringcomputation，gameandspecialeffectsinmovie.Nevertheless，withtheincreasingrequirementsoffidelityofthemotionandinteractionofsubjects，fluidsimulationrequrresmorerealisticsimulationresults.However，themotionofrealfluidandobjectsisoftendifficulttobeexpressedbysimpleprocess，whichcanonlybereaiizedbythephysicallawsintherealworld.Physics-basedfluidsimulationisoneofthemostchallengingresearchissuesduetoitscomplexityintheoryandcomputation.Amongnumerousfluidsimulationmodels，smoothedparticlehydrodynamics#SPH)basedfluidsimulationisthemostwidelyusedmeshlessmethodtosimulatetherealisticfluidmotion.Asasimpleandflexiblemethod，SPHparticlecharacteristicsmakeitnaturallysuitabeforf8uidsimuationmodeingandfeasibetodea8withmanyfuidsimuationcomputationprobemsinarobustandstabeway.Withtheimprovementofthereaisticeffectoffuidsimuation，thedetailedsimulationoffluidhasattractedmoreandmoreattention.But，currentSPHmethodcannotguaranteetheaccuratesurfacetensoncomputaton，whchreducesthefdeltyofthefludsimulation,failtoachieverealisticfluidmotionandmayleadtothedistortionoffluidsimulationinrealsimulationapplicationscenarios.Surfacetensonisthemostcommonandimportantphysicalpropertythatexhibitsthemicroscopicpropertiesoffluidsandisofgreatsignificanceinrealizingfluidsimulationwithhighfidelity.ItisstlllachallengetoaccuratelycalculatethesurfacetensionfortheSPHbasedfluidsimulation.Thenumerica1instabiiityintheSPHbasedfree—surfaceliquidsimulationcausesunrealsurfacetension,particlesclumping,andthinfeaturesrupture,whichseriouslyaffectsthevisualeffectsofsimulationatthefreesurfaceofthefluid.Inthefreesurfaceofthefluid,thelackofneighborhoodparticlescauseunderestimatedsurfaceparticlesdensitywhichiscacuatedbytheSPHmethod,resutinginnegativepressureandthus8eadingtothephenomenonofnon—physica!particescumping.Totackethisissue,weproposeafastandreaiisticmethodforfluidsurfacetensioncomputationinSPHbasedfluidsimulaion.Firstly,thefludsurfacepartclesareeffcentlydetectedandthesurfacenormalvectorofthepartcles9sprecselydefned.Then,basedonthedetectedsurfacepartcles,wedefnehelevelsetgrdbyusngthefreesurfacedefntonandsurfacenormalvector.Toallevatethedffculty9nhandlesurfacepartcledensty,weproposethelevel—setbaseddenstycorrectonmethodfornear—surfacepartclestomakethecalculatonattheboundarymoreaccurate.Fnally,nordertoachevefastandeffectvesurfacetenson,weproposeamodelwhchcombnesthenear—surfacepartclenteractonforcewththemnmumsurfaceareatocalculatetheaccuratesurfacetensonofthefluid.Theexperimentalresultsshowthatourmethodcanallevatetheparticleclumpingcausedbythenumericainstabiityatthefuidfreesurface,andourmethodcanachievefastandreaisticsurfacetensionsimuationinSPHfuid.Keywordssmoothedpariclehydrodynamics；fluidsimulaion；surfacetension；neai4surfaceinter4particleforce;densrtycorrection1引言隨著計算機硬件和圖形圖像技術的不斷進步,計算機仿真技術得以迅速發(fā)展，在游戲、電影特效、軍事仿真、工程仿真等領域，計算機仿真技術得到廣泛的應用.然而隨著各方面需求的日益增加，人們對計算機仿真中物體運動與交互的真實感提出更高的要求，希望得到更加真實的仿真效果.但是，真實的物體運動往往難以用過程表述，所以此時只有借助真實世界物理規(guī)律得以實現(xiàn).流體仿真作為計算機圖形學和虛擬現(xiàn)實的研究熱點之一，在動畫、游戲及電影特效等領域具有重要應用.由于流體仿真的復雜性，基于物理的流體仿真是最具挑戰(zhàn)性的研究問題之一.目前流體仿真有多種方法，其中SPH（SmoothedParticleHydrodynamics,光滑粒子流體動力學）方法是近年來應用最廣泛的無網格法.SPH作為一種簡單而靈活的方法，其粒子特性使其天然適合流體的仿真建模，能夠以魯棒穩(wěn)定的方式處理很多流體模擬問題.SPH方法將流體劃分為一系列離散的粒子，然后應用光滑核估計粒子動態(tài)物理量，每個粒子的空間導數(shù)可以從其周圍粒子推斷得到，最終粒子按其計算所得速度進行運動相較于傳統(tǒng)歐拉方法,SPH方法不需要計算復雜的網格，可更方便地處理復雜的固體邊界問題.但是SPH方法在流體仿真中不能保證流體的不可壓縮性，導致模擬場景中的流體仿真效果失真.隨著對流體模擬真實感效果要求的提升，流體細節(jié)仿真越來越受到人們的關注.表面張力是展現(xiàn)流體微觀特性最常見和最重要的物理特性，表面張力的模擬成為當前流體仿真中的研究熱點之一.從微觀角度來看,流體分子受到各個方向其他分子的作用力，該作用力隨著分子間距離的遠近分別表現(xiàn)引力和斥力作用.在流體內部，該作用力合力為零，流體分子處于平衡狀態(tài)；而在流體表面，該作用力合力不為零，方向指向流體內部，從而形成表面張力.為了實現(xiàn)逼真的流體仿真效果，當前基于SPH算法的流體自由表面處仿真仍舊是一個挑戰(zhàn)，其主要包括非真實準確人工表面張力、粒子聚集以及細節(jié)斷裂等.在流體自由表面處，鄰域粒子的缺失會引發(fā)SPH算法計算得到的表面粒子密度偏小，這種低估的密度值導致負壓力從而產生粒子的非物理聚集現(xiàn)象.為了解決該問題，Schechter和Bridson⑵提出GhostSPH方法，在自由表面創(chuàng)建一層幽靈粒子,并將該層幽靈粒子納入各物理量計算中.但是，幽靈粒子地不斷采樣較大提升了算法的時間復雜度.除了表面粒子密度計算問題,SPH表面張力計算也是處理流體表面問題的一個重要方面.當前基于SPH的表面張力仿真計算主要分兩類：宏觀模型和微觀模型.宏觀模型中，基于顏色場(ColorField)表面粒子檢測、表面法向及張力計算存在一定的不準確性,自由表面處的數(shù)值不穩(wěn)定導致表面粒子聚集，破壞了流體仿真中的細節(jié)之處，造成仿真效果不逼真的結果.微觀模型中，通過計算粒子間的聚斥力來仿真流體分子間作用力.如Hartakovsky等人⑶和Ak-mci等人也通過構造聚斥力模型模擬分子間作用力以對表面張力進行仿真.近年來，有關流體表面的仿真應用研究上還出現(xiàn)一些關于多相流的建模研究.如Ren等人通過體積分數(shù)表示模擬不同相流體的空間分布，由此改進混合模型實現(xiàn)對多相流流動的、通用的、魯棒的SPH流體仿真.而后Ren等人#句又通過將前述多相流仿真獲得的體積分數(shù)數(shù)據(jù)與簡單有效的二次氣泡模擬相結合，提出一種靈活的氣泡建模方法.Ren等人#0$通過補償可見氣體粒子密度和力的計算提出一種氣態(tài)流體的快速SPH仿真方法.Yan等人#)$通過拓展現(xiàn)有多相流SPH框架，實現(xiàn)了流體與可變形物體、顆粒狀物質的仿真.本文在Marrone等人##$關于表面粒子檢測及Akinci等人由關于表面張力定義的工作基礎上，提出一種快速逼真的流體表面張力仿真方法.通過該方法，自由表面處數(shù)值不穩(wěn)定引起的粒子聚集現(xiàn)象得到緩解，最終實現(xiàn)了快速逼真流體表面張力仿真效果.本文主要貢獻如下：通過對流體進行快速準確地表面檢測及法向量計算，提出一種結合近表面粒子間作用力與最小化表面積的模型計算流體的表面張力，所得表面張力更為真實，并緩解了自由表面處的數(shù)值不穩(wěn)定性問題；提出一種基于水平集網格的近表面粒子密度修正方法，以更加準確的處理邊界處計算，進一步緩解自由表面處的數(shù)值不穩(wěn)定性.2相關工作國內外研究人員提出了很多嵌入到流體仿真框架中的表面張力計算方法.早期Morris?提出從連續(xù)表面力模型出發(fā)模擬表面張力.后續(xù)Miiller等人#%%$在此工作基礎上提出基于顏色場方法進行表面張力求解.他們提出的表面張力模型側重于最小化表面曲率，為每一個粒子計算法向量，以此確定力的方向，而力的大小取決于粒子位置處曲率大小，而且對于流體內部的粒子，其法向量方向為任意的.強洪夫等人#12$在Morris提出的表面張力SPH方法基礎上，通過引入邊界核函數(shù)插值方法對邊界法向和曲率的計算進行修正，得到表面張力修正方程組，以實現(xiàn)穩(wěn)定的SPH表面張力計算.近些年，Yu等人#%'$提出在流體粒子表面重建網格，然后將網格頂點的曲率應用于其鄰域流體粒子的表面張力計算,證明在表面網格稠密的情況下可準確計算流體表面張力.吳恩華等人#1!$根據(jù)流體密度場提出一種基于平均曲率的表面張力模型，能夠有效地模擬水滴的表面張力現(xiàn)象.Yang等人#%&$通過體積保持的平滑網格方式對流體表面張力進行計算，實現(xiàn)穩(wěn)定流體模擬，特別是具有較強表面張力的場景仿真.這些工作的流體表面張力效果基本都依賴于準確的流體邊界檢測以及精確的表面法向量計算，同時曲率的計算涉及二次求導，這對粒子的分布情況要求比較高.為避免對表面曲率及表面法向量的依賴，研究人員提出了通過鄰域粒子間的相互作用力定義新的表面張力模型#"%*$.該方法可避免錯誤的法向量以及曲率計算，但是僅僅用粒子間作用力不能保證表面積最小化，無法實現(xiàn)準確表面張力計算.Clavet等人#10$提出一種改進的SPH模型對表面張力進行仿真，并提出一種基于距離的引力作用，獲得了較好的視覺效果.劉棟等人#%)$提出基于SPH的粒子間相互作用力模型以計算表面張力，獲得了與基于連續(xù)表面張力模型的SPH方法類似的表面張力效果.近年來，He等人#%#$提出最小化表面能量的表面張力模型，獲得較好的視覺效果.Akmci等人由提出一種基于SPH流體的通用表面張力模型，該模型通過結合內聚力項和最小化表面積項模擬表面張力.Wang等人#"$$在Akinci方法基礎上提出基于II-SPH的流體表面張力和吸附力方法，能夠較好地模擬流體的表面微觀特性，改善了數(shù)值穩(wěn)定性和提高

了計算效率.Yang等人⑵"提出一種可實時模擬不同界面處交互作用的成對力SPH模型，該模型可在一次仿真中同時處理多種類型的交互作用，并且緩解了自由表面的粒子聚集現(xiàn)象.由于SPH方法對于所有粒子的場量都會通過核插值計算其空間導數(shù)，所以在流體自由表面處的計算中，存在非完整核半徑引發(fā)的密度計算偏小問題，導致數(shù)值不穩(wěn)定、表面處的粒子聚集現(xiàn)象，影響流體表面視覺效果.為了解決數(shù)值不穩(wěn)定的缺陷，Monaghan等人提出一種人工壓力項以緩解該問題，但導致產生較為虛假的表面張力.Oger等人!提出了使用對SPH場量空間導數(shù)的高階估計，然而較大地提升了計算復雜度.Kulasegaram等人提出了一種不使用邊界粒子的SPH密度方程以防止粒子聚集.它們通過支持半徑中固體邊界占據(jù)的面積引入修正項，從而避免在邊界處計算錯誤密度.Feldman等人!通過將球體中的密度函數(shù)與二維空間中多個線段表示的邊界相結合計算密度修正項.Ihmsen等人通過壓力項修正了粒子在邊界處的分布.Schechter和Bridson!""提出GhostSPH方法，通過在液氣、液固交接處采樣幽靈粒子實現(xiàn)密度修正解決非完整核半徑問題，然而對粒子地不斷采樣等步驟較大地增加了計算代價.He等人通過結合雙尺度壓力估計和各向異性壓力過濾，實現(xiàn)了魯棒的、具有較好視覺效果的流體模擬.近年來，F(xiàn)ujisava等人⑵"提出一種無需采樣任何邊界內部粒子的密度修正方法，該方法有效地減少了邊界粒子的內存消耗，但卻增加了流固耦合的計算時間.Ren等人通過補償可見氣體粒子密度和力的計算提出一種表面附近粒子密度修正的方法.3流體仿真計算基礎3.1SPH/WCSPHSPH算法將流體劃分為一系列離散的粒子，然后應用光滑核估計粒子動態(tài)物理量，最終粒子按照其本身速度進行運動，每個粒子的空間導數(shù)可以從其周圍粒子推斷得到.對于粒子z,其某一場量A可以表示為周圍支持半徑九內所有粒子的加權和：D(7)="D(7)竺W(r>6) (1)其中=7—7>,7?分別為粒子i,j的位置，皿表示周圍粒子>的質量，)>表示為周圍粒子j的密度,A(ri)是粒子i的場量.函數(shù)R(r>6)是核函數(shù).本文使用WCSPH!1*"作為基本框架.WCSPH提出了一種基于Hait方程的弱可壓縮SPH方法.通過對連續(xù)性方程求解,得到粒子密度更新如下：>VW7>6) (2)其中>>=>i—>j;>i>j分別表示粒子i打的速度.得到粒子密度后可得粒子運動方程如下:坐=—(8, DW(r6),s(3)1 // "12[Vi'’>#,/1*v°/dz > \piP>/其中8 分別表示粒子i,j的壓強.壓強通過Hait方程得到，定義如下:(4)其中，P$為初始密度，c,為流體中的音速，％為一可調系數(shù)，原文建議值為7.使用Hait方程計算壓強使得密度只能小范圍波動，保證了流體的弱可壓縮性.為了提高流體仿真中的數(shù)值穩(wěn)定性,WCSPH定義粘性力如下:d> >DW(r>,h), vHr><0Z=B> HH⑸$， v>r>,0/'j=T+(|r>|2+62/'j=T+(|r>|2+62(6)P+Pj6"為防止分母為0的修正項，e=0.01.本文方法以WCSPH為基本框架計算壓力與粘性力，而密度計算使用SPH密度求和方式得到：p,="m}W(r>6) (7)j在計算過程中，得到初始粒子密度后，再使用本文第5節(jié)提出的密度修正方法進行粒子的密度修正以緩解由自由表面粒子錯誤的密度計算導致的數(shù)值不穩(wěn)定.3.2快速表面粒子檢測與水平集網格構建3.2.1表面粒子檢測與表面法向量計算精確地定義SPH流體粒子的邊界是對流體表面張力進行仿真的基礎.不準確的表面粒子檢測將導致后續(xù)不準確的水平集網格,繼而導致不準確的表面張力計算.本文基于將文獻［9］的方法引入到SPH流體仿真領域，將其用于快速準確地檢測SPH流體表面粒子.首先定義重歸一化矩陣尋找自由表面周圍的粒子:B(r,)=$"DW(rt],h)Irt]-V]%41 (8)其中-V,是粒子，的體積.定義N為所有流體粒子集合,$8N為N的子集，表示自由表面粒子集合.Doring[Z8]證明矩陣64的最小特征值人取決于被計算點(周圍的粒子，的空間分布.理論上，當粒子遠離流體區(qū)域時，該特征值趨向于0,當越趨近于流體內部，該特征值越趨近于1.所以可以通過人值確定粒子處于流體區(qū)域的位置.通過計算B4的最小特征值人，可定義三種互補的子集：舊為飛濺的粒子集合，人值較?。?為遠離自由表面的內部粒子集合，人值較大；'為靠近自由表面的粒子或者呈現(xiàn)非均勻分布區(qū)域的粒子集合球值相對居中.本文通過選取人值的閾值來精準區(qū)分表面粒子，為了確定較為合適的人值閾值，進行了多個場景的仿真測試.其最終不同子集對應閾值如下:在實際實驗中，當某一粒子沒有鄰域粒子時，重歸一化矩陣是不存在的.此種情況下，本文通過設置該粒子的法向量為0進行處理$上述為流體表面粒子粗檢測并得到粒子的法向量.利用下式可以判斷粒子是否屬于表面粒子：若|r>>%Zh,\r1T<h kB#，若|r><%2h,arccos( K(L$\m4其他情況 K'$$(12)&EJA<0.2J'J0.2<A<0.75 (9)(&J*〉0.75在流體內部區(qū)域，鄰域粒子的梯度核函數(shù)之和接近于0;當粒子接近于自由表面時，鄰域粒子的梯度核函數(shù)之和接近該粒子的表面法向量$29%.再次使用重歸一化矩陣得到表面粒子的法向量:n(r()=〃(r()I&(r,)|(10)其中#(11)其中，$為表面粒子集合，T為粒子'法向方向且與'距離h的點,n表示粒子'的法向量,N'表示粒子(周圍的粒子集合$圖1所示為采用顏色場以及本文方法對流體表面進行檢測的視覺效果對比.其展示了本文所引用方法的有效性.對比圖1(a)、(c)與圖1(b)、(d)可見,采用顏色場的表面檢測結果不準確,導致表面不連續(xù),可導致后續(xù)基于表面粒子的法向量計算及表面張力計算不準確,而本方法能夠精確地檢測流體表面粒子.在該場景中，隨著粒子數(shù)的增加，本方法與顏色場方法檢測表面粒子用時如表1所示.從表中可知，隨著粒子數(shù)的增加,兩種方法耗時都隨之增加,并且本方法相對顏色場方法耗時也在增加,相對時間占比隨粒子數(shù)增多而增大.本方法與顏色場方法都能快速對流體表面粒子進行檢測，在實驗中，本方法雖比顏色場方法平均用時多19.4%，但其在接受范圍內，并且也取得了更加精確的結果.(a)顏色場方法&(r()=-B(r()"DW(r>,h)-V,(c)顏色場方法圖1流體粒子表面檢測結果表1表面粒子檢測用時對比粒子數(shù)顏色場耗時/ms本方法耗時/ms相對耗時增加占比/%100001.5131691118200003.949450614I130000!7.4148I73017I840000!11I90114I79824I350000!17I70222I87429I2平均19I43.2.2水平集網格構建基于檢測得到的SPH流體表面粒子,Marrone等人$9%通過計算level-set函數(shù)以劃分流體內部或流體外部區(qū)域的網格節(jié)點.對于每一個接近于自由表面粒子集$的節(jié)點N,檢測其最近的表面粒子Fn,則標量2nF%可以由下式計算：d%pN=3%—r%）#f% （13）此處#4表示Fn處的自由表面法向量.這樣對于每一個網格頂點：N，可以定義level-set函數(shù)：，_1 2NF%<_262（：n）=v2nFn%—"6—"6V2nFnV26（14）+1 2nFn,26該函數(shù)在流體內部為正，在流體外部為負，在自由表面處為0（因為2酒%=。）.D快速表面張力計算與密度修正自由表面的存在意味著SPH插值精度應當保持到接近流體域邊界處.當實現(xiàn)SPH插值計算，如壓力梯度、速度散度時，接近流體表面處被以與流體內部同樣的方式所處理.然而表面處的計算常常被人們所忽視,導致壓力梯度不能收斂到接近自由表面處的正確值.相反，本文主張通過快速、精確的表面粒子檢測，對沒有數(shù)值截斷的自由表面流動進行仿真.圖2為本文方法計算模式.首先，根據(jù)前述的快速、精確的表面粒子檢測方法獲得精確的自由表面定義與表面法向.根據(jù)獲得的自由表面定義和表面法向構建流體區(qū)域水平集網格，以區(qū)分流體內外區(qū)域的網格頂點.本文通過SPH密度求和方式得到流體密度，所以對于擁有非完整核半徑的流體粒子，由于鄰域粒子的缺失，其計算所得密度將偏小.通過前述構建的水平集網格，可對粒子密度進行修正，詳見4.3節(jié).實際情況下，不是所有的粒子密度計算都是錯誤的，只有靠近表面的粒子密度才是被低估的.本表面檢測 其他外力V水平集粒子運動WCSPH壓力、粘性力L與修正近表面最小化表面積近表面粒子間作用力圖2本文方法計算模式文在精確表面粒子檢測的基礎上，將擁有非完整核半徑的粒子定義為近表面粒子，即自由表面一個核半徑內的所有流體粒子，如圖3中實線和虛線（含線上）之間的流體粒子所示.所以，密度修正只需對近表面粒子操作，以減少密度修正計算量，提高計算效率.為了實現(xiàn)快速、有效的表面張力，本文提出了一種近表面粒子間作用力，詳見4.1節(jié).為了更好地模擬流體的表面微觀特性，還需加入近表面最小化表面積項保證流體表面積最小化，詳見4.2節(jié).最后，通過WCSPH框架求解壓力、粘性力及其他外力后，計算流體粒子的加速度驅動粒子運動.圖3近表面粒子4.1近表面粒子間作用力模型在毛細、水珠等現(xiàn)象中，表面張力扮演著十分重要的角色.微觀本質上，表面張力是分子間作用力的一種表現(xiàn).當分子間距離較近時，表現(xiàn)為斥力作用，當距離較遠時，則表現(xiàn)為引力作用.真實分子間作用力是十分復雜的.勒那德-瓊斯（Lennarc—Jones）模型是對分子間作用力的一種粗糙近似，盡管如此，其仍然具有一定的參考意義，其近似模型圖如圖4所示.根據(jù)古根海姆（Guggenheim）模型，液體表面層是一個有一定厚度的過渡區(qū)，既占有體積又具有物質.所圖4Lennard-Jones模型近似圖

以在表面層，分子較為稀疏，分子間距離較遠，其作用力表現(xiàn)為吸引力，表面粒子間的相互吸引作用形成了表面張力.而對于液體內部分子，其處于一種平衡狀態(tài)，所受分子間作用力合力為零.為了準確計算粒子間的表面張力作用，在Hartakovsky等人提出的表面張力計算模型基礎上，本文提出一種精確快速的近表面粒子間作用力模型用于計算表面張力.在Tartakovsky表面張力模型中，表面張力定義如下：坐=%F,(15)其中是粒子(的鄰域粒子間相互作用合力，定義(15)F,="f> (16)而粒子間相互作用力定義如下:=bcos(16T")；'“<(17)［o, |r>|>h其中，氣是粒子(和粒子，之間的作用力強度系數(shù).圖4為Lennard-Jones模型曲線大致形狀.圖5為本文模型及其他兩個對比模型粒子間相互作用力曲線形狀對比圖，其中，橫軸表示到所計算粒子的中心距離與核半徑的比值|r|/h,縱軸表示相互作用力大小，正值表示斥力作用，負值表示引力作用，并對最大值進行歸一化處理.對比圖4Lennard-Jones模型曲線與圖5Akinci方法"4#和Tartakovsky方法曲線可知，雖然后者具有引力、斥力作用，但是一方面該作用力隨粒子間距離的變化與Lennarc-Jones模型有較大出入；另一方面，考慮其在SPH模型的計算中，該模型也忽視了對核半徑的使用.很多文獻通過構造復雜的樣條函數(shù)以解決上述兩點問題，而事實上可以通過直接將文獻"3］模型引入到SPH模型中解決問題.令粒子間作用力強度系數(shù)為S,并引入圖5本文模型與其他模型曲線形狀比較

核半徑，定義粒子間相互作用如下:4=S Ir>I)7>,W(r>,h)(18)P> \h /|r>從Lennard-Jones模型曲線可以看出，引力作用與斥力作用的極值比值小于1,且在引力作用區(qū)域，引力大小從0到引力極值部分曲線的平均梯度絕對值是大于極值右側曲線的平均梯度絕對值的.對比圖4與圖5本文方法曲線，將原模型引入到SPH模型后，引力極值和斥力極值的比值、引力梯度變化程度都更加貼近于Lennard-Jones模型.同時，因為該模型的斥力作用，由非完整核半徑導致的表面數(shù)值不穩(wěn)定的粒子非物理聚集現(xiàn)象，也將得到緩解.除此之外，考慮到表面張力本質上是分子間作用力的結果，而該作用力在流體內部合力為零.所以為了提高計算效率，只需對表面粒子計算粒子間作用力即可.但是另一方面，考慮到在SPH模型中，只針對表面粒子計算將導致數(shù)值空間梯度的不連續(xù)，為了保證空間梯度的連續(xù)性，本文在計算粒子間作用力時針對所有具有非完整核半徑的近表面粒子.這樣的處理方式，在保證計算效率的基礎上，保證了數(shù)值穩(wěn)定性，使仿真效果更加自然.而表面粒子則可以通過前述3.2節(jié)的方法進行快速、精確表面檢測.4.2最小化表面積項近表面粒子間作用力模型，能夠構造類似于分子間作用力的效果，可以有效地解決液體表面粒子的非物理聚集問題.但是僅僅使用近表面粒子作用力模型模擬表面張力，將導致流體最終所處的形狀取決于粒子的初始分布，不能保證流體表面積最小化.所以為了更好地模擬流體的表面微觀特性，獲得更加真實的表面張力仿真效果，還需加入最小化表面積項保證流體表面積最小化.Akinci等人"4］提出一種最小化流體表面積方法定義如下：F>=—&m#> (19)其中，&表示可調的最小化表面積項系數(shù)=#—n>,n(表示表面粒子(和，的法向量.在SPH中，從粒子分布計算曲率是十分容易出錯的.為了避免明確的曲率計算,Aknci通過對平滑顏色場梯度的SPH估計，計算粒子法向量：mn=4"」DW(.>,h) (20)>P>與一般顏色場方法的不同在于該方法多一項獨立的、針對法向量的縮放系數(shù)4,表示粒子的法向量

數(shù)量級與曲率是成某一比例的.雖然該方法沒有顯式地計算曲率，但其與一般顏色場計算方法并無二致.而基于顏色場的粒子法向量是存在一定不準確性的.本文通過'2節(jié)重歸一化矩陣得到的精確表面粒子及表面法相量代替顏色場方法，可以得到比顏色場方法更加精確的粒子法向量.本文將式(10)所得法向量引入最小化表面積項中，將其推廣到SPH模型中并對稱化，得到最小化表面積項定義如下：——nt]W(rt]6) (21)>p對于上式容易得知，在流體平坦區(qū)域和流體內部區(qū)域，該項為零.因為對于平坦區(qū)域n(,得到#功=n(—n>20;對于內部粒子而言，由于n,20和n>20,故該項也為零.所以在計算該項時，為了提高計算效率，亦只針對表面粒子進行計算；而為了保證空間梯度的連續(xù)性，將所有具有非完整核半徑的粒子也加入計算之中.在提高效率的同時保證了數(shù)值穩(wěn)定性.4.3基于水平集網格的密度修正在SPH流體仿真中，計算所得的液氣、液固交界處粒子密度比其真實密度小,這種錯誤的密度值導致負壓力，引起粒子的聚集現(xiàn)象,即所謂的數(shù)值不穩(wěn)定.概括來講,流體表面粒子的非完整核半徑導致流體密度的錯誤計算是自由表面處數(shù)值不穩(wěn)定與粒子非物理聚集現(xiàn)象的原因.而自由表面處的粒子非物理聚集對近表面粒子間作用力的計算有較大的影響.近表面粒子間作用力模型只有在粒子分布較為均勻自然的情況下才能得到較好的仿真效果.考慮到式(14)所得的level-set函數(shù)值可以判斷某水平集網格是否屬于流體內部，本文通過定義流體體積占比系數(shù)R,對表面粒子進行密度修正.流體體積占比系數(shù)E定義如下："s(2))-V)W(r,,h)R,= (22)"-V1W(r,,h)其中，△匕表示粒子j的體積大小，對于4(2,),當2>40時4(2,)=1,否則4(2,)=0.密度修正公式定義如下："m,W(rj6)pp='—E (23)對于流體內部粒子E,21,無需進行密度修正,對流體近表面粒子,R,<1.故同樣只需針對具有非完整核半徑的粒子進行密度修正即可.圖6通過level-set函數(shù)進行密度修正算法1.流體仿真迭代.8.ENDFOR1.WHILEanimatingDO9.FOREACHpartice(DO2.updateputeviscosity圖6通過level-set函數(shù)進行密度修正算法1.流體仿真迭代.8.ENDFOR1.WHILEanimatingDO9.FOREACHpartice(DO2.updateputeviscosity,gravityandexterna1putedensity11.IFparticle,isneaLsurfaceparticleTHEN4.findfree-surfaceparticles12.!computenear-surface9nterpartclesforces5.updatelevel-setputetheforcewhichminimizesthesurface6.FOREACHneaLsurfaceparticle,DOarea7.!densitycorrection14.!ENDIFENDFORcomputevelocityup：＜te□。，讓仍口，1）iendwhile5實驗結果為了驗證本文所述方法的有效性!我們實驗對比了相同場景下Akinci方法等與本文方法的效果，并采用本方法實現(xiàn)了逼真的多液塊融合、水滴浸潤、滴水效果以及流體流動仿真效果.SPH的光滑核半徑6取值為'R,R為SPH采樣粒子的半徑,對于所有的實驗，本文還融入了Akmci等人提出的流固耦合方法.本文的實驗平臺為：（1）硬件：英特爾Intel5-4570M（3.20GHz）CPU,8GB內存，顯卡GeforceGT755M；（2）軟件：VisualStudio2012,Mitsuba,OpenGL.表2所示為本文方法對上述場景進行實驗所得的時間性能統(tǒng)計.表21000幀平均時間性能統(tǒng)計 （單位:ms）場景流體粒子數(shù)（最多）表面檢測水平集網格表面張力密度修正總計多液塊融合150003.6933.32!.356.19102.56水滴浸潤50001.0528.98滴水效果10000!8822.581.733.0862.1!液態(tài)巧克力流經兔子300008.3135.5!3.635.36127.15水流沖擊兔子!00001!.8286.6911.1222.16336.61圖7圖7多液塊融合仿真圖7為本文方法多液塊的融合仿真過程.對四團液塊在失重環(huán)境下的融合過程進行了逼真的模擬.液塊相互接觸吸弓I,在表面張力的作用下，最終融合為一個整體，形成了一個圓球形的水珠.多團液塊在無重力的作用下最終融合形成球形水珠的變化過程，說明了本文表面張力包括近表面作用力與近表面最小化表面積的有效性.對于此場景，本文還分別對GhostSPH、Akinci山方法、Yang等人#21$方法進行了模擬.仿真過程效果與本文方法相似，無實質不同，最終都形成了球形水珠.圖8為四種方法最終粒子分布效果圖.除了GhostSPH方法外，其他三種方法最終都形成了完美的球形水珠.說明本文方法與Akinci、Yang方法，在表面張力有效性方面相同.而GhostSPH最終粒子分布不是完美的球形,即沒有實現(xiàn)表面積最小化.本文實現(xiàn)的GhostSPH表面張力是基于顏色場的，說明基于顏色場的表面張力效果不真實，分析其主要原因在于其易出錯的表面曲率與表面法向計算.盡管后三種方法在最終效果上相似，但在實驗過程中發(fā)現(xiàn)液塊融合過程中形成穩(wěn)定球形水珠的速度不一致.在本文實驗環(huán)境下，仿真時間步長timestep=0.001s,Akinci方法與本文方法約在800幀