班曉娟, 周靖, 王笑琨, 劉斯諾, 徐衍睿
(北京科技大學(xué) 計(jì)算機(jī)與通信工程學(xué)院, 北京 100083)
基于物理的流體模擬動(dòng)畫是計(jì)算機(jī)圖形學(xué)領(lǐng)域的重要方面之一,流體模擬與交互是流體動(dòng)畫制作的核心問題,該類研究在游戲動(dòng)畫、虛擬現(xiàn)實(shí)等人機(jī)交互領(lǐng)域有著廣泛的應(yīng)用和較強(qiáng)的落地價(jià)值.
湍流細(xì)節(jié)模擬一直是流體模擬中的熱點(diǎn)問題.SPH作為一種拉格朗日方法能夠高效地模擬具有極大變形特征的流體,十分適合于模擬混亂的湍流.目前,許多研究通過強(qiáng)制不可壓縮性[1-3],解決了壓強(qiáng)力求解的數(shù)值耗散問題. 然而,黏性力計(jì)算過程中的數(shù)值耗散問題仍然存在,使得諸如湍流等高頻細(xì)節(jié)被平滑掉,極大地影響了實(shí)驗(yàn)結(jié)果的逼真度[4-5].
渦旋是湍流的主要組成成分,渦度是渦旋的重要屬性之一,但傳統(tǒng)的SPH方法很少考慮渦度場(chǎng).物理學(xué)中將渦旋分為兩類:剛體渦旋和無旋渦旋,見圖1所示.在兩種渦旋上分別放置兩個(gè)剛性元件.圖1(a)為渦旋像剛體一樣旋轉(zhuǎn),具有非零渦度.圖1(b)為無旋渦旋,除中心外渦度處處為0.剛體渦旋的渦度是角速度的兩倍,而無旋渦旋除中心點(diǎn)外渦度處處為0.黏性力對(duì)流體粒子的旋轉(zhuǎn)運(yùn)動(dòng)起主要作用,而用于描述流場(chǎng)內(nèi)旋轉(zhuǎn)趨勢(shì)的渦度場(chǎng)并未被傳統(tǒng)SPH方法納入考慮.這就是在SPH方法中湍流細(xì)節(jié)損失的主要原因.
圖1 剛體渦旋和無旋渦旋Fig.1 Rigid body vortexes and irrotational vortex
為了保持流體表面上的復(fù)雜湍流和渦旋細(xì)節(jié),一些研究通過上采樣方法,對(duì)表面點(diǎn)進(jìn)行精細(xì)化處理[6],增加表面分辨率.基于網(wǎng)格的方法使用自適應(yīng)體積網(wǎng)格達(dá)成該目標(biāo)[7].然而,以上方法都僅在表面添加細(xì)節(jié),并沒有考慮到流體內(nèi)部運(yùn)動(dòng).目前,機(jī)器學(xué)習(xí)方法被引入到流體模擬領(lǐng)域[8],以實(shí)現(xiàn)高分辨率湍流的快速合成.但因?yàn)椴灰?guī)則的離散化,機(jī)器學(xué)習(xí)方法極難與基于粒子的方法結(jié)合.
Navier-Stokes方程在理論上能夠準(zhǔn)確地描述流體運(yùn)動(dòng),但離散化過程仍會(huì)導(dǎo)致計(jì)算誤差.在SPH方法中,N-S方程中的黏性力不能直接離散化并用于流體模擬,因此人工黏度被廣泛地使用.黏性計(jì)算過程中的數(shù)值耗散將導(dǎo)致渦旋和湍流細(xì)節(jié)損失,然而目前只有少數(shù)研究涉及該類型的數(shù)值耗散.
為了解決上述問題,本文引入了一種基于黏性的湍流模擬方法,通過黏度計(jì)算過程中的能量耗散率來校正渦度場(chǎng).借助流函數(shù),利用渦度場(chǎng)的增強(qiáng)來提高速度,從而恢復(fù)逼真且可控的渦旋和湍流效果.該方案不僅可以增加現(xiàn)有的渦旋,還可以在潛在位置產(chǎn)生新的湍流.此外,渦度場(chǎng)的增強(qiáng)和速度的校正可與大部分類型的SPH方法兼容,且不影響隨后的壓力計(jì)算步驟.
MONAGHAN用SPH方法解決了自由表面流動(dòng)的模擬問題[9],為流體模擬奠定了基礎(chǔ).隨后MULLER等[10]提出使用理想氣體狀態(tài)方程來模擬流體,但其壓縮性對(duì)結(jié)果造成了負(fù)面影響.BECKER和TESCHNER[11]提出了一種改進(jìn)的弱壓縮SPH (weakly compressible SPH)方法來降低壓縮性,然而其效率受到時(shí)間步長(zhǎng)的限制.HE等[12]和SOLENTHALER,PAJAROLA[1]提出了一種預(yù)測(cè)-校正方法,通過迭代步驟計(jì)算壓強(qiáng),保證流體的不可壓縮性.這種思想也被應(yīng)用在基于位置的流體模擬方法中(position based fluids, PBF)[13].IHMSEN等[2]遵循壓力預(yù)測(cè)策略從而提出了隱式不可壓縮SPH (implicit incompressible SPH, IISPH)方法.此外,BENDER和KOSCHIER[3]提出了一種強(qiáng)制零壓縮和無發(fā)散速度場(chǎng)的方法.本文使用IISPH方法作為實(shí)驗(yàn)中計(jì)算效率和穩(wěn)定性比較的基礎(chǔ).高頻細(xì)節(jié)恢復(fù)問題一直以來都是流體模擬中的重要課題[14-16].對(duì)于歐拉方法,STAMS方案[17]首次將流體模擬引入計(jì)算機(jī)圖形學(xué)中.然而,時(shí)間和空間的一階精度使得該方法飽受數(shù)值耗散的影響.KIM等[18]對(duì)達(dá)到更高階近似做出了貢獻(xiàn).
隨后研究的混合方法進(jìn)一步減少了數(shù)值耗散.ZHU和BRIDSONS[19]在不考慮某些旋轉(zhuǎn)耗散的情況下,顯著地消除了不可壓縮流體的FLIP (fluid-implicit-particle)平流中數(shù)值耗散.JIANG等[20]使用混合方法恢復(fù)了流體大部分的旋轉(zhuǎn)運(yùn)動(dòng).雖然上面提到的模擬方法能夠從宏觀角度上處理這個(gè)問題,但歐拉方法和拉格朗日方法都難以模擬湍流等高頻細(xì)節(jié).因此,出現(xiàn)了一些專門用于模擬湍流細(xì)節(jié)的方法.這些方法可分為3大類,即上采樣法,渦度限制法和拉格朗日渦旋法[21].
上采樣法旨在粗離散化下添加高頻率表面細(xì)節(jié).MERCIER等[6]提出了一種獨(dú)立的后處理方法: 在基于粒子的流體表面上施加細(xì)湍流.首先在曲率評(píng)估之后接種高分辨率表面點(diǎn),然后在粗粒子上產(chǎn)生精細(xì)的表面波.之后,EDWARDS等[7]使用自適應(yīng)間斷法提出了一種基于網(wǎng)格的流體的自適應(yīng)體積網(wǎng)格的方法.總之,上面提到的上采樣法只能改善表面效果.一些機(jī)器學(xué)習(xí)方法也可應(yīng)用于流體模擬領(lǐng)域,如卷積神經(jīng)網(wǎng)絡(luò)(CNN)[22],能夠在高分辨率源和粗糙模擬結(jié)果上合成精細(xì)湍流.
渦度限制方法旨在識(shí)別現(xiàn)有的渦旋并減少耗散.該方法通過添加新的強(qiáng)制項(xiàng)來增加目標(biāo)位置的速度,并強(qiáng)制進(jìn)行渦旋的旋轉(zhuǎn).LENTINE等[23]在中改進(jìn)了渦度限制,同時(shí)保持了能量和動(dòng)量的守恒.JANG等[24]使得渦度限制多層次化,從而獲得更好的結(jié)果.盡管渦度限制方法是保留渦旋的簡(jiǎn)單方法,但是它無法產(chǎn)生額外的湍流細(xì)節(jié).更重要的是,他們傾向于通過可調(diào)參數(shù)向系統(tǒng)添加過多的能量,這會(huì)違反能量守恒.
拉格朗日渦旋粒子法是建立在Navier-Stokes方程的渦度表示上的,理論上不受數(shù)值耗散的影響.該方法可通過粒子[25-26],曲線[27],細(xì)絲[28]和表面[29]實(shí)現(xiàn)湍流細(xì)節(jié)的捕捉與恢復(fù),但難以處理邊界問題,例如非剛性障礙物和自由表面等.ZHU等[30]提出了一種模擬運(yùn)動(dòng)物體周圍渦旋細(xì)節(jié)的方法,但該方法需要借助歐拉網(wǎng)格實(shí)現(xiàn).之后GOLAS等[31]處理歐拉網(wǎng)格上的邊界問題.以上這些方法的另一個(gè)缺點(diǎn)是需求解Biot-Savart積分或矢量值泊松方程來恢復(fù)速度場(chǎng).在2018年,BENDER等[21]提出了用于非黏性流體的微極元流體模型,該模型可以捕獲流體粒子的微旋轉(zhuǎn)并獲得較好的湍流視覺效果.ZHANG等[32]提出了一種對(duì)流運(yùn)動(dòng)學(xué)的集成渦度(integrated vorticity of convective Kinematics, IVOCK)方法,通過計(jì)算平流中的渦度損失來恢復(fù)能量耗散.這些方法都是采用優(yōu)化計(jì)算來恢復(fù)速度場(chǎng)的.受其流函數(shù)的啟發(fā),本文將其擴(kuò)展到基于拉格朗日的流體模擬,能夠從渦度場(chǎng)中得到速度的細(xì)化.
本文的方法可被視為拉格朗日渦旋方法.不僅可以增加現(xiàn)有的渦流,還可以在新渦流的潛在位置產(chǎn)生湍流.此外,方法不需要求解Biot-Savart積分或矢量值泊松方程.相反,借助流函數(shù)通過渦量場(chǎng)來細(xì)化速度,極大地提高了仿真效率.
在拉格朗日粒子描述下,流體粒子的加速度可表示為
(1)
(2)
(3)
在傳統(tǒng)的流體模擬中,N-S方程中的黏性項(xiàng)用于計(jì)算黏性力,可以按如下方式離散:
(4)
式中:d為維度;xij=xi-xj為粒子i和粒子j之間的距離;vij=vi-vj;h為支持半徑;Wij為光滑核.因?yàn)槔绽顾阕硬贿m用于粒子系統(tǒng),所以人工黏度被廣泛采用.
由于人工黏度造成的耗散,流體系統(tǒng)的能量不斷減小.某些區(qū)域的渦旋難以維持,丟失了許多細(xì)節(jié).隨著誤差的累積,一段時(shí)間后會(huì)出現(xiàn)更明顯的失真.為了保持這些細(xì)節(jié),本文提出了基于黏性力校正的SPH湍流模擬方法.
本文方法與拉格朗日渦旋方法有密切關(guān)系:通過渦度恢復(fù)速度場(chǎng).但本文通過流函數(shù)避免求解泊松方程,并利用黏度相關(guān)算法直接修改速度.在本文中,旋度場(chǎng)是通過粒子撞擊前后黏性力導(dǎo)致的動(dòng)能變化率來調(diào)整的.然后通過旋轉(zhuǎn)場(chǎng)校正每個(gè)粒子的速度.該方法不僅可以放大現(xiàn)有的渦流,還可以在新渦流的潛在位置產(chǎn)生湍流.在本文方法中,渦度場(chǎng)的求解和速度的校正屬于平流步驟,不影響隨后的壓力投影步驟,保證不可壓縮性.
在本文算法中,除了速度之外,每個(gè)粒子都有一個(gè)渦度屬性γ.渦度是用于描述局部位置的卷曲的矢量場(chǎng).通常,它由速度場(chǎng)的旋度表示,在粒子系統(tǒng)中,渦度表示粒子旋轉(zhuǎn).離散化形式為
(5)
式中:m為質(zhì)量;Wij為光滑核函數(shù),在本文中使用3次樣條函數(shù).
在流動(dòng)過程中,湍流和渦旋的細(xì)節(jié)由渦度場(chǎng)所展現(xiàn)的.無黏性流體只能在理想條件下存在.實(shí)際上,需要考慮黏性力對(duì)流體的影響.在傳統(tǒng)的SPH方法中,黏性力的計(jì)算導(dǎo)致數(shù)值耗散,這造成了速度和渦度的計(jì)算誤差.若不考慮黏性力,就無法準(zhǔn)確描述流體的運(yùn)動(dòng),在一些場(chǎng)景中許多細(xì)節(jié)如爆炸、攪動(dòng)等將丟失.
隨著時(shí)間的推移,流體系統(tǒng)中的能量將在流體粒子之間轉(zhuǎn)移.當(dāng)粒子的速度發(fā)生變化時(shí),渦度也發(fā)生變化并影響其運(yùn)動(dòng)軌跡.在此過程中,如果未計(jì)算渦度的特定值,或者未將其反饋應(yīng)用于速度的糾正,則會(huì)丟失某些重要細(xì)節(jié).在預(yù)測(cè)過程中,黏性力和外力影響著粒子的運(yùn)動(dòng).一般情況下,外力是重力.因?yàn)橹亓]有旋轉(zhuǎn),所以本文方法只需考慮由黏度造成渦度所帶來的速度變化,從而保持粒子的速度并且減少能量耗散.在流體粒子的運(yùn)動(dòng)過程中,它們受到相鄰粒子的力的影響并導(dǎo)致速度的變化.這可歸納為圖2兩種情況.
圖2 兩種能量轉(zhuǎn)換Fig.2 Two kinds of energy conversion
在平流-投影模型中,平流過程部分主要用于解決除壓力之外的黏性力和外力.此時(shí),由速度變化(黏性力造成)引起的能量增量和能量變化率平方根可表示為
(6)
(7)
對(duì)于粒子系統(tǒng),總能量可以表示為動(dòng)能和勢(shì)能之和,與熱量無關(guān).顯然,粒子的動(dòng)能是E=1/2mv2,其中一部分是可以用角速度表示,如圖3.由于力的作用,粒子之間的碰撞將改變角速度.因?yàn)榻撬俣饶芰渴莿?dòng)能的一部分,在獲得粒子的動(dòng)能變化率之后,若將角速度表示的能量以相同的比率變化,則其值仍然在合理的范圍內(nèi).對(duì)于由角速度表示的能量是E=Iω2,角速度的變化率可以表示為λ.
圖3 角速度所攜帶的能量Fig.3 The energy carnied by the angular velocity
在物理學(xué)中,理想的渦旋可以分為兩類: 剛體渦旋和無旋渦旋.其中,剛體的渦度是固定軸角速度的兩倍,無旋渦旋的渦度處處為0.因此,渦度與角速度有一定的關(guān)系,并且至多等于角速度的兩倍.當(dāng)角速度增加P倍時(shí),渦度也將增加f(P).渦度增量由以下函數(shù)表示
δγ=ηλγ
(8)
式中,δγ表示渦度校正,湍流水平增強(qiáng)參數(shù)η用于表示系統(tǒng)中渦旋的程度.參數(shù)越大,湍流效果越明顯.
如圖3所示,角速度所攜帶的能量是動(dòng)能的一部分.當(dāng)流體靜止時(shí),兩種速度的能量都為0; 當(dāng)只有線速度而沒有角速度.此時(shí),角速度能量為0,動(dòng)能為E=1/2mv2,當(dāng)流體旋轉(zhuǎn)時(shí),總動(dòng)能為E=1/2mv2=Iω2,其中由角速度表示的能量等于動(dòng)能.因此,角速度所表示的能量總是小于或等于動(dòng)能.
(9)
(10)
最后
(11)
現(xiàn)已在渦度場(chǎng)和速度場(chǎng)之間建立了聯(lián)系.通過式(8)的黏性影響程度限制式(11),避免增加的能量大于耗散能量.在三維空間,ψ和u可以分別設(shè)置為(ψx,ψy,ψz)和(u,v,w).
在傳統(tǒng)方法中,需要求解泊松方程以將δγ轉(zhuǎn)換為速度.Green函數(shù)提供了一種半解析解.在三維空間中表示為δφ=-ρ的泊松方程可以轉(zhuǎn)換為
(12)
式中:φ為計(jì)算通量;ρ為密度.
通過亥姆霍茲分解以數(shù)學(xué)方式推廣,流函數(shù)是速度場(chǎng)v的向量勢(shì).向量勢(shì)ψ可以根據(jù)上面Green函數(shù)結(jié)果定義:
(13)
考慮到標(biāo)量泊松問題:
ψx(p)=0,p→∞
(14)
所以可以使用以下離散化表達(dá)式進(jìn)行空間離散化[30]:
(15)
式中:γj渦度;Vj為采樣點(diǎn)體積;xi-xj代表第i個(gè)采樣點(diǎn)和第j個(gè)采樣點(diǎn)位置之間的距離.
(16)
在本節(jié)中,對(duì)本文方法在若干場(chǎng)景中進(jìn)行了測(cè)試,與傳統(tǒng)的IISPH方法和微極元方法(MP solver)[20]進(jìn)行了比較.所有上述模型都與IISPH方法集成,密度誤差控制在0.1%以內(nèi).此外,本文采用了AKINCI等[33]提出的邊界處理方法.通過C ++實(shí)現(xiàn)基于物理的仿真及重構(gòu)[34]框架,動(dòng)畫由Blender渲染.仿真平臺(tái)為一圖形工作站,采用Intel Xeon E5-2 637 v2 (15M高速緩存,3.50 GHz @ 4核心) CPU,80 GB RAM,NVIDIA Quadro M5 000 GPU.
圖4 左列分別顯示了不同湍流水平控制參數(shù)η=0,0.2,0.4,0.6的模擬結(jié)果.隨著η的增加,球體右側(cè)的湍流效果逐漸變得豐富且在實(shí)驗(yàn)場(chǎng)景左側(cè)的水面波紋顯得更加細(xì)致.這說明本文方法可以使用可控參數(shù)恢復(fù)或添加的湍流細(xì)節(jié).η作為一控制因子用于限制對(duì)于黏度造成的能量耗散的彌補(bǔ),確保其不會(huì)過分增加能量,這樣整個(gè)系統(tǒng)能量是趨于守恒的.較高的η值表示粗采樣和大時(shí)間步長(zhǎng)下的視覺效果.因此,在隨后的對(duì)比實(shí)驗(yàn)中,本文方法中的參數(shù)均以0.6確保實(shí)驗(yàn)質(zhì)量.實(shí)驗(yàn)還將本文方法與微極元方法進(jìn)行了比較.圖4 右列顯示具有不同傳遞系數(shù)νt=0,0.2,0.4,0.6的微極元方法下的實(shí)驗(yàn)效果,可見人為地增加傳遞系數(shù),并且在大參數(shù)情況下的湍流細(xì)節(jié)的豐富變化.但是,在相同的控制參數(shù)下,微極元方法的視覺效果并不像本文方法那么明顯.此外,隨著參數(shù)越高,數(shù)值不穩(wěn)定性問題開始出現(xiàn).當(dāng)參數(shù)值大于0.6時(shí),其數(shù)值計(jì)算結(jié)果不夠穩(wěn)定,開始出現(xiàn)無法控制的混沌湍流.因此,在隨后的對(duì)比實(shí)驗(yàn)中,微極元方法中的參數(shù)均為0.4以確保穩(wěn)定性.當(dāng)控制參數(shù)η=0或νt=0時(shí),算法可以被視為標(biāo)準(zhǔn)SPH方法,在本文中為IISPH方法.
圖4 本文方法和微極元方法Fig.4 This method and the MP solver
圖5展示了在1.072 M流體粒子情況下的傳統(tǒng)IISPH方法和本文方法的模擬結(jié)果.在IISPH的模擬結(jié)果中,湍流細(xì)節(jié)的數(shù)量很快消失.在本文方法中,湍流更明顯.特別是在接下來的兩排中,球的左側(cè)也有明顯的波浪.本文方法更好地捕捉了湍流和渦旋細(xì)節(jié).
圖5 IISPH下潰壩實(shí)驗(yàn),傳統(tǒng)方法和η=0.6的本文方法Fig.5 3D breakdam experiment under IISPH and our method with η=0.6
為了證明方法的穩(wěn)定性和效果,實(shí)驗(yàn)使用動(dòng)態(tài)邊界條件模擬了兩個(gè)復(fù)雜場(chǎng)景,并將本文方法與微極元方法進(jìn)行了比較.在圖6中,將木板緩緩放入水中然后橫向移動(dòng).在這里,可觀察到使用IISPH和微極元方法不會(huì)產(chǎn)生和保持復(fù)雜流動(dòng)和湍流效應(yīng).與微極元方法相比,本文方法以物理上合理的方式增加能量,并在自由表面上產(chǎn)生逼真的湍流細(xì)節(jié).此外,使用本文方法可以逼真地處理此場(chǎng)景中的高速的流體粒子和高度湍流運(yùn)動(dòng),這證明了本文方法出色的穩(wěn)定性.
圖6 木板撥水實(shí)驗(yàn),488 k流體粒子Fig.6 Board interacts with water containing 488 k particles
由于以前渦流或湍流方案傾向于增加過多的能量并且不穩(wěn)定,實(shí)驗(yàn)進(jìn)行了與2.9 M流體粒子交互的沉船實(shí)驗(yàn).在圖7中,船被放入水中然后逐漸下沉.船的巨大勢(shì)能轉(zhuǎn)化為流體粒子的巨大動(dòng)能.當(dāng)船落入水中時(shí),水首先猛烈地運(yùn)動(dòng),隨著船沉沒,水逐漸平靜下來,最終達(dá)到穩(wěn)定狀態(tài).在這個(gè)實(shí)驗(yàn)中,可發(fā)現(xiàn)使用IISPH和微極元方法產(chǎn)生了微弱的湍流細(xì)節(jié),并且因?yàn)閿?shù)值耗散,湍流效果很快丟失.相比之下,本文方法顯示出更直觀的運(yùn)動(dòng),模擬出具有逼真的湍流效應(yīng)的流體運(yùn)動(dòng).當(dāng)船沉沒時(shí),流體流動(dòng)不會(huì)受到明顯的阻礙.當(dāng)船完全在水下時(shí),有著明顯的湍流和渦流.此外,隨著時(shí)間的推移,流體逐漸平靜下來,這意味著本文方法不會(huì)添加過多的能量.
圖7 沉船實(shí)驗(yàn)Fig.7 Boat sinks into the water
這兩個(gè)場(chǎng)景表明,本文方法可以在處理激烈碰撞等極端條件時(shí)保持穩(wěn)定,并在物理上保留動(dòng)能.且本文方法不僅放大了現(xiàn)有的渦旋,還產(chǎn)生了新的渦旋.此外,湍流細(xì)節(jié)的計(jì)算開銷也保持在了合理范圍內(nèi),表1列出了IISPH、微極元方法和本文方法在兩個(gè)典型場(chǎng)景下的計(jì)算時(shí)間,與傳統(tǒng)IISPH方法相比,MP方法與本文方法都相應(yīng)增加了計(jì)算時(shí)間,部分場(chǎng)景下本文方法的計(jì)算效率優(yōu)于MP方法.
表1 不同方法時(shí)間對(duì)比
本文提出了一種基于黏性的SPH湍流模擬方法,該方法從能量耗散中推導(dǎo)出渦度變化,從而恢復(fù)速度場(chǎng).該方法可以顯著地減少數(shù)值耗散,增加現(xiàn)有渦流并在潛在位置產(chǎn)生新的湍流.使用控制參數(shù),可以很容易地控制湍流的精細(xì)度.實(shí)驗(yàn)證明,相比于傳統(tǒng)SPH方法和微極元SPH方法,本文方法能更好地增強(qiáng)湍流效果.此外,即使使用大時(shí)間步長(zhǎng)和大粒子半徑,本文方法仍可以較好地保持能量,可以在極端模擬條件和復(fù)雜的大規(guī)模場(chǎng)景下穩(wěn)健運(yùn)行.目前實(shí)驗(yàn)表明其適用于不可壓縮流體,但無數(shù)據(jù)證明該方法在其他種類流體中的穩(wěn)定性.
未來計(jì)劃將該方法擴(kuò)展到其他流體,進(jìn)一步研究微極元模型,并嘗試將本文的模型與其合并,因?yàn)槲O元模型在粗糙條件下顯示出巨大的作用并且與黏度密切相關(guān).此外將提高計(jì)算精度作為另一個(gè)研究方向.盡管離散流函數(shù)和渦量計(jì)算在某種程度上可以產(chǎn)生精細(xì)效果,但仍然可以找到具有有效且更精確的計(jì)算方法.