基于PCISPH的流體粒子飛濺改進(jìn)方法

鈕倩倩,林綠開(kāi),李 毅,2*

(1.溫州大學(xué) 計(jì)算機(jī)與人工智能學(xué)院,浙江 溫州 325000;2.浙江大學(xué) 計(jì)算機(jī)科學(xué)與技術(shù)學(xué)院,浙江 杭州 310058)

0 引 言

對(duì)于流體模擬研究經(jīng)歷了較長(zhǎng)時(shí)間并不斷衍生出新的方法。從宏觀上來(lái)說(shuō),流體的模擬可分為基于物理的方法和非物理的方法。非物理的方法往往依賴(lài)于純粹的數(shù)學(xué)模型,需要通過(guò)觀測(cè)數(shù)據(jù)對(duì)實(shí)際數(shù)據(jù)進(jìn)行建模,從而達(dá)到模擬仿真的效果,該類(lèi)方法沒(méi)有引入物理控制方程,所以結(jié)果通常存在動(dòng)量、能量不守恒的問(wèn)題。

基于物理的方法往往通過(guò)離散的數(shù)值方法對(duì)控制方程來(lái)進(jìn)行求解。根據(jù)采樣點(diǎn)的不同分布可以將其細(xì)分為歐拉方法和拉格朗日方法。歐拉方法是一種網(wǎng)格方法,典型的代表是有限差分方法[1],這種方法放棄了描述運(yùn)動(dòng)對(duì)象本身,轉(zhuǎn)而描述運(yùn)動(dòng)對(duì)象所處空間場(chǎng)的變化。由于該方法相當(dāng)于模擬了一種空間內(nèi)固定的背景網(wǎng)格,通過(guò)運(yùn)動(dòng)對(duì)象對(duì)網(wǎng)格的影響來(lái)近似描述運(yùn)動(dòng)對(duì)象,所以不必考慮對(duì)象復(fù)雜的運(yùn)動(dòng)。拉格朗日粒子法是將客觀世界中的流體看成由許多粒子組成的整體,該方法的基本思路是首先依照流體運(yùn)動(dòng)的物理定律,計(jì)算出粒子受到的外力和粒子之間的相互作用力,其次根據(jù)牛頓第二定律(Newton's Second Law of Motion-Force and Acceleration),計(jì)算出每個(gè)粒子在時(shí)間步長(zhǎng)內(nèi)的位置變化量,最后通過(guò)位置變化量來(lái)模擬出一段流體運(yùn)動(dòng)軌跡。主流的拉格朗日[2]粒子法包括 SPH(Smoothed Particle Hydrodynamics,光滑粒子流體動(dòng)力學(xué))[3-5]和PBF (Position Based Fluids,基于位置的流體模擬方法)[6]等。該文主要針對(duì)SPH進(jìn)行研究。

SPH是一種基于積分插值理論的插值方法,通過(guò)離散采樣點(diǎn)(粒子)來(lái)近似連續(xù)場(chǎng)量的值和導(dǎo)數(shù),這些粒子自身攜帶某些場(chǎng)量,粒子處的場(chǎng)量通過(guò)從相鄰粒子的值加權(quán)平均得到。有限差分法要求粒子排列在規(guī)則網(wǎng)格上,SPH則可以通過(guò)任意排布的粒子來(lái)求得近似。每個(gè)離散粒子都占據(jù)問(wèn)題域的一小部分。

SPH源于計(jì)算天體物理學(xué)[7]領(lǐng)域,用于模擬天體物理現(xiàn)象。Desbrun等人[8]將SPH引入計(jì)算機(jī)圖形學(xué)領(lǐng)域,用于模擬可變形體。Müller[9]將SPH引入流體模擬,并在此后擴(kuò)展到各種問(wèn)題的模擬。

在計(jì)算機(jī)圖形學(xué)中,Solenthaler和Pajarola提出了一種有效的SPH變體,即PCISPH(Predictive-Corrective Incompressible SPH,預(yù)測(cè)校正不可壓縮SPH)[10]。在這種SPH變化中,通過(guò)迭代預(yù)測(cè)和校正密度來(lái)增強(qiáng)不可壓縮性。壓力的確定應(yīng)使壓力能夠減小密度波動(dòng)。預(yù)測(cè)校正不可壓縮SPH(PCISPH)允許更大的時(shí)間步長(zhǎng),如ISPH(隱式SPH)。此外,PCISPH解決了粒子級(jí)的不可壓縮性,因此不需要求解泊松方程。預(yù)測(cè)校正不可壓縮SPH(PCISPH)在一個(gè)模型中同時(shí)具有WCSPH(Weakly Compressible SPH,弱可壓縮SPH)[11]和ISPH[12]的優(yōu)點(diǎn),即每次物理更新的計(jì)算成本低,時(shí)間步長(zhǎng)大。通過(guò)流體傳播估計(jì)的密度值,并以實(shí)現(xiàn)不可壓縮性的方式更新壓力。一旦達(dá)到每個(gè)單獨(dú)粒子的先前用戶(hù)定義的密度變化限制,傳播就會(huì)停止。其中WCSPH是最直接的顯式求解算法,基于連續(xù)性方程、動(dòng)量方程和狀態(tài)方程,直接計(jì)算作用在各個(gè)顆粒上的壓力、體力、粘聚力和表面張力的合力,繼而求出加速度進(jìn)而結(jié)合邊界條件進(jìn)行速度位置的更新。

由于飛濺場(chǎng)景中流體的密度和壓力的速率變化非常大,因此模擬對(duì)離散化解的精度要求很高。傳統(tǒng)的不可壓縮SPH[13]在許多場(chǎng)景中都有很好的性能,例如電影效果等。但由于數(shù)值耗散,通常不可能長(zhǎng)期精確模擬流體的運(yùn)動(dòng)變化,在對(duì)數(shù)值精度要求較高的場(chǎng)景中仍有很大的改進(jìn)空間。

該文提出一種基于PCISPH的流體粒子飛濺的改進(jìn)方法。使用Taichi[14]框架實(shí)現(xiàn)流體仿真,太極編程語(yǔ)言大大提高了并行編程的生產(chǎn)力。

1 準(zhǔn)備工作

1.1 不可壓縮泊松方程模型

二維場(chǎng)景中自由表面流體的方程由質(zhì)量守恒方程和動(dòng)量守恒方程組成[15]。該文研究的是基于粒子的流體模擬,因此將控制方程轉(zhuǎn)換為拉格朗日形式,如公式1和公式2所示。

(1)

(2)

其中,ρ表示流體的密度,t是時(shí)間,u是流體的速度,P是壓力,g是重力加速度,μ是剪切應(yīng)力(包括粘滯力)。由于粘滯力的保留,式(1)和式(2)可以應(yīng)用于牛頓流體和非牛頓流體。

1.2 速度預(yù)測(cè)

在納維-斯托克斯(Navier-Stokes Equation)方程中,如果只考慮應(yīng)力和重力,則從以下公式可以得到粒子的速度和位置。公式如下:

(3)

u*=ut+Δu*

(4)

r*=r+u*Δt

(5)

其中,ut和r表示粒子在時(shí)間t上的速度和位置。u*和r*表示預(yù)測(cè)步驟結(jié)束時(shí)粒子的中間速度和位置。Δu*表示預(yù)測(cè)期間粒子速度的變化量,Δt表示時(shí)間步長(zhǎng)。

1.3 壓力密度校正

方程1中的計(jì)算不能保證流體的不可壓縮性,在X*位置處的密度ρ*由SPH得出:

(6)

其中,ri是起始粒子位置,rj是起始粒子的鄰居粒子位置,mj是鄰居粒子質(zhì)量,h是光滑核函數(shù)半徑。在密度ρ*和恒定的不可壓縮流體密度ρ0之間存在偏差。

因此,校正步驟旨在將流體密度調(diào)整到初始值。在計(jì)算過(guò)程中,壓力項(xiàng)將用于更新中間步驟中計(jì)算的粒子速度。公式如下:

(7)

ut+1=u*+Δu**

(8)

其中,Δu**表示校正過(guò)程中粒子速度的變化量,ρ*表示預(yù)測(cè)步驟之后和校正步驟開(kāi)始之前的中間密度,ut+1和Pt+1表示在時(shí)間t+1時(shí)刻的速度和壓力值。

強(qiáng)制執(zhí)行不可壓縮性所需的壓力項(xiàng)來(lái)自方程1,如下:

(9)

獲得壓力的泊松方程如下:

(10)

2 飛濺改進(jìn)

壓力校正。

由于預(yù)測(cè)校正不可壓縮SPH(PCISPH)是從預(yù)測(cè)的密度變化中得出壓力變化,因此計(jì)算密度變化是整個(gè)過(guò)程的核心。在SPH中密度計(jì)算如下:

(11)

其中,ri表示目標(biāo)粒子位置,rj表示鄰居粒子位置。

該文研究流體的飛濺場(chǎng)景。由于該場(chǎng)景中流體運(yùn)動(dòng)強(qiáng)烈,壓力變化率和變化范圍大,容易受到誤差累積的影響。因此,提出了改進(jìn)的泊松壓力方程源項(xiàng),如下所示:

(12)

其中,上標(biāo)符號(hào)*表示在預(yù)測(cè)中計(jì)算所得。

將提出的密度變化解代入傳統(tǒng)不可壓縮流體的泊松壓力方程,得到了改進(jìn)的泊松壓力方程式求解模型:

(13)

其中,ρ*表示預(yù)測(cè)中的流體密度。

公式12簡(jiǎn)化如下:

(14)

結(jié)合公式13和公式14得到:

(15)

(16)

又在預(yù)測(cè)校正不可壓縮SPH(PCISPH)中通過(guò)預(yù)測(cè)密度變化來(lái)校正壓力,公式如下:

由公式16得:

(18)

(19)

3 實(shí)驗(yàn)結(jié)果與分析

3.1 實(shí)驗(yàn)環(huán)境

該文利用 Python環(huán)境與編輯器,采用 Anaconda3集成Python3和Taichi環(huán)境,使用 PyCharm 代碼編輯平臺(tái)進(jìn)行數(shù)據(jù)預(yù)處理以及后續(xù)網(wǎng)絡(luò)框架的搭建與模型訓(xùn)練。

Taichi起步于 MIT 的計(jì)算機(jī)科學(xué)與人工智能實(shí)驗(yàn)室(CSAIL),設(shè)計(jì)初衷是便利計(jì)算機(jī)圖形學(xué)研究人員的日常工作,幫助他們快速實(shí)現(xiàn)適用于 GPU 的視覺(jué)計(jì)算和物理模擬算法。Taichi 是嵌入于 Python的,使用即時(shí)編譯(JIT)架構(gòu)(如 LLVM,SPIR-V),將Python源代碼轉(zhuǎn)化為GPU或CPU的原生指令,在開(kāi)發(fā)時(shí)和運(yùn)行時(shí)均提供優(yōu)越性能。Taichi的一大優(yōu)勢(shì)在于可移植性,支持多種后端。

3.2 實(shí)驗(yàn)對(duì)比

由于飛濺場(chǎng)景中流體密度和壓力的變化率非常大,因此模擬對(duì)離散解的精度要求很高。針對(duì)基于粒子的流體仿真在模擬流體飛濺時(shí)的數(shù)值不穩(wěn)定性,提出改進(jìn)的泊松壓力方程解,改進(jìn)流體飛濺現(xiàn)象。將上文推導(dǎo)結(jié)果應(yīng)用到模型中,得到改進(jìn)后的流體模擬。下圖是按時(shí)間順序基于預(yù)測(cè)校正不可壓縮SPH的粒子飛濺的傳統(tǒng)方法和改進(jìn)方法的效果對(duì)比。圖1為8 000粒子數(shù)的飛濺效果對(duì)比,圖2為4 000粒子數(shù)的飛濺效果對(duì)比。

圖1 8 000粒子飛濺效果對(duì)比

圖2 4 000粒子飛濺效果對(duì)比

表1是在8 000個(gè)粒子的三維流體模擬實(shí)驗(yàn)中,在GPU上進(jìn)行了9個(gè)獨(dú)立實(shí)驗(yàn)。

表1 在GPU上運(yùn)行的實(shí)驗(yàn)數(shù)據(jù) fps

基于NVIDIA和Windows平臺(tái),記錄實(shí)驗(yàn)結(jié)果,如表1所示(fps表示程序運(yùn)行幀率)。

圖3是根據(jù)表1中的數(shù)據(jù)繪制出的折線圖。其中實(shí)驗(yàn)編號(hào)表示時(shí)間先后順序。從圖中可以看出,在文中規(guī)定的測(cè)試環(huán)境下,傳統(tǒng)方法和改進(jìn)方法的幀率基本都處于45到60之間,但改進(jìn)方法的幀率比傳統(tǒng)方法的幀率略小一點(diǎn),但是兩種方法均達(dá)到實(shí)時(shí)性仿真效果的要求。

圖3 基于GPU的實(shí)驗(yàn)程序性能比較

由上面效果比較可以看出,改進(jìn)后的流體粒子飛濺較之前的傳統(tǒng)流體粒子飛濺來(lái)說(shuō)有所變化,改進(jìn)后的飛濺的粒子更有規(guī)律性,符合現(xiàn)實(shí)中在重力作用下的效果。流體仿真模擬即要求仿真模擬結(jié)果更接近現(xiàn)實(shí),改進(jìn)后的更符合。從表1與圖3可以看出模擬均達(dá)到實(shí)時(shí)仿真要求,但改進(jìn)后的性能較之改進(jìn)前的較差一點(diǎn)。

實(shí)驗(yàn)核心部分即為校正密度壓力值,表2為傳統(tǒng)方法中的密度誤差和改進(jìn)方法中的密度誤差(數(shù)據(jù)選擇來(lái)自傳統(tǒng)方法和改進(jìn)方法飛濺的粒子數(shù)量相差較多的情況)。

表2 校正密度誤差值

圖4是由表2數(shù)據(jù)繪制出的折線圖,其中實(shí)驗(yàn)編號(hào)表示運(yùn)行時(shí)間先后順序。

圖4 校正的密度誤差值對(duì)比

雖然改進(jìn)后的幀率校之傳統(tǒng)方法有所下降,但是基于PCISPH的流體仿真是通過(guò)預(yù)測(cè)校正密度壓力值從而得到仿真效果。實(shí)驗(yàn)中通過(guò)改進(jìn)密度誤差從而改進(jìn)壓力值實(shí)現(xiàn)流體仿真,從表2和圖4中可以看出,改進(jìn)后的密度誤差明顯小于傳統(tǒng)方法的密度誤差,并且基本上在千分之五左右。改進(jìn)后密度誤差變小,增強(qiáng)了流體的不可壓縮性,減小了壓力噪聲,得到了更真實(shí)、更穩(wěn)定的飛濺仿真效果。從而驗(yàn)證了改進(jìn)方法的有效性。

4 結(jié)束語(yǔ)

該文主要針對(duì)預(yù)測(cè)校正不可壓縮SPH(PCISPH)流體仿真進(jìn)行粒子飛濺改進(jìn),針對(duì)改進(jìn)減小密度誤差,增強(qiáng)流體的不可壓縮性,減小了壓力噪聲,從而得到了優(yōu)于原模型的仿真效果。但是該模型也存在相應(yīng)的問(wèn)題,模型不夠穩(wěn)定,該文添加較多粒子來(lái)減小不穩(wěn)定帶給實(shí)驗(yàn)結(jié)果的影響,在未來(lái)將對(duì)這些方面進(jìn)行深入研究,以達(dá)到穩(wěn)定且較大時(shí)間步長(zhǎng)的流體模擬。同時(shí),該文沒(méi)有對(duì)該模型進(jìn)行渲染,未來(lái)將對(duì)模型進(jìn)行渲染。