鈕倩倩,林綠開,李 毅,2*
(1.溫州大學 計算機與人工智能學院,浙江 溫州 325000;2.浙江大學 計算機科學與技術學院,浙江 杭州 310058)
對于流體模擬研究經歷了較長時間并不斷衍生出新的方法。從宏觀上來說,流體的模擬可分為基于物理的方法和非物理的方法。非物理的方法往往依賴于純粹的數學模型,需要通過觀測數據對實際數據進行建模,從而達到模擬仿真的效果,該類方法沒有引入物理控制方程,所以結果通常存在動量、能量不守恒的問題。
基于物理的方法往往通過離散的數值方法對控制方程來進行求解。根據采樣點的不同分布可以將其細分為歐拉方法和拉格朗日方法。歐拉方法是一種網格方法,典型的代表是有限差分方法[1],這種方法放棄了描述運動對象本身,轉而描述運動對象所處空間場的變化。由于該方法相當于模擬了一種空間內固定的背景網格,通過運動對象對網格的影響來近似描述運動對象,所以不必考慮對象復雜的運動。拉格朗日粒子法是將客觀世界中的流體看成由許多粒子組成的整體,該方法的基本思路是首先依照流體運動的物理定律,計算出粒子受到的外力和粒子之間的相互作用力,其次根據牛頓第二定律(Newton's Second Law of Motion-Force and Acceleration),計算出每個粒子在時間步長內的位置變化量,最后通過位置變化量來模擬出一段流體運動軌跡。主流的拉格朗日[2]粒子法包括 SPH(Smoothed Particle Hydrodynamics,光滑粒子流體動力學)[3-5]和PBF (Position Based Fluids,基于位置的流體模擬方法)[6]等。該文主要針對SPH進行研究。
SPH是一種基于積分插值理論的插值方法,通過離散采樣點(粒子)來近似連續場量的值和導數,這些粒子自身攜帶某些場量,粒子處的場量通過從相鄰粒子的值加權平均得到。有限差分法要求粒子排列在規則網格上,SPH則可以通過任意排布的粒子來求得近似。每個離散粒子都占據問題域的一小部分。
SPH源于計算天體物理學[7]領域,用于模擬天體物理現象。Desbrun等人[8]將SPH引入計算機圖形學領域,用于模擬可變形體。Müller[9]將SPH引入流體模擬,并在此后擴展到各種問題的模擬。
在計算機圖形學中,Solenthaler和Pajarola提出了一種有效的SPH變體,即PCISPH(Predictive-Corrective Incompressible SPH,預測校正不可壓縮SPH)[10]。在這種SPH變化中,通過迭代預測和校正密度來增強不可壓縮性。壓力的確定應使壓力能夠減小密度波動。預測校正不可壓縮SPH(PCISPH)允許更大的時間步長,如ISPH(隱式SPH)。此外,PCISPH解決了粒子級的不可壓縮性,因此不需要求解泊松方程。預測校正不可壓縮SPH(PCISPH)在一個模型中同時具有WCSPH(Weakly Compressible SPH,弱可壓縮SPH)[11]和ISPH[12]的優點,即每次物理更新的計算成本低,時間步長大。通過流體傳播估計的密度值,并以實現不可壓縮性的方式更新壓力。一旦達到每個單獨粒子的先前用戶定義的密度變化限制,傳播就會停止。其中WCSPH是最直接的顯式求解算法,基于連續性方程、動量方程和狀態方程,直接計算作用在各個顆粒上的壓力、體力、粘聚力和表面張力的合力,繼而求出加速度進而結合邊界條件進行速度位置的更新。
由于飛濺場景中流體的密度和壓力的速率變化非常大,因此模擬對離散化解的精度要求很高。傳統的不可壓縮SPH[13]在許多場景中都有很好的性能,例如電影效果等。但由于數值耗散,通常不可能長期精確模擬流體的運動變化,在對數值精度要求較高的場景中仍有很大的改進空間。
該文提出一種基于PCISPH的流體粒子飛濺的改進方法。使用Taichi[14]框架實現流體仿真,太極編程語言大大提高了并行編程的生產力。
二維場景中自由表面流體的方程由質量守恒方程和動量守恒方程組成[15]。該文研究的是基于粒子的流體模擬,因此將控制方程轉換為拉格朗日形式,如公式1和公式2所示。

(1)

(2)
其中,ρ表示流體的密度,t是時間,u是流體的速度,P是壓力,g是重力加速度,μ是剪切應力(包括粘滯力)。由于粘滯力的保留,式(1)和式(2)可以應用于牛頓流體和非牛頓流體。
在納維-斯托克斯(Navier-Stokes Equation)方程中,如果只考慮應力和重力,則從以下公式可以得到粒子的速度和位置。公式如下:

(3)
u*=ut+Δu*
(4)
r*=r+u*Δt
(5)
其中,ut和r表示粒子在時間t上的速度和位置。u*和r*表示預測步驟結束時粒子的中間速度和位置。Δu*表示預測期間粒子速度的變化量,Δt表示時間步長。
方程1中的計算不能保證流體的不可壓縮性,在X*位置處的密度ρ*由SPH得出:
(6)
其中,ri是起始粒子位置,rj是起始粒子的鄰居粒子位置,mj是鄰居粒子質量,h是光滑核函數半徑。在密度ρ*和恒定的不可壓縮流體密度ρ0之間存在偏差。
因此,校正步驟旨在將流體密度調整到初始值。在計算過程中,壓力項將用于更新中間步驟中計算的粒子速度。公式如下:

(7)
ut+1=u*+Δu**
(8)
其中,Δu**表示校正過程中粒子速度的變化量,ρ*表示預測步驟之后和校正步驟開始之前的中間密度,ut+1和Pt+1表示在時間t+1時刻的速度和壓力值。
強制執行不可壓縮性所需的壓力項來自方程1,如下:

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

(10)
壓力校正。
由于預測校正不可壓縮SPH(PCISPH)是從預測的密度變化中得出壓力變化,因此計算密度變化是整個過程的核心。在SPH中密度計算如下:
(11)
其中,ri表示目標粒子位置,rj表示鄰居粒子位置。
該文研究流體的飛濺場景。由于該場景中流體運動強烈,壓力變化率和變化范圍大,容易受到誤差累積的影響。因此,提出了改進的泊松壓力方程源項,如下所示:
(12)
其中,上標符號*表示在預測中計算所得。
將提出的密度變化解代入傳統不可壓縮流體的泊松壓力方程,得到了改進的泊松壓力方程式求解模型:

(13)
其中,ρ*表示預測中的流體密度。
公式12簡化如下:

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

(15)

(16)
又在預測校正不可壓縮SPH(PCISPH)中通過預測密度變化來校正壓力,公式如下:

由公式16得:

(18)
(19)
該文利用 Python環境與編輯器,采用 Anaconda3集成Python3和Taichi環境,使用 PyCharm 代碼編輯平臺進行數據預處理以及后續網絡框架的搭建與模型訓練。
Taichi起步于 MIT 的計算機科學與人工智能實驗室(CSAIL),設計初衷是便利計算機圖形學研究人員的日常工作,幫助他們快速實現適用于 GPU 的視覺計算和物理模擬算法。Taichi 是嵌入于 Python的,使用即時編譯(JIT)架構(如 LLVM,SPIR-V),將Python源代碼轉化為GPU或CPU的原生指令,在開發時和運行時均提供優越性能。Taichi的一大優勢在于可移植性,支持多種后端。
由于飛濺場景中流體密度和壓力的變化率非常大,因此模擬對離散解的精度要求很高。針對基于粒子的流體仿真在模擬流體飛濺時的數值不穩定性,提出改進的泊松壓力方程解,改進流體飛濺現象。將上文推導結果應用到模型中,得到改進后的流體模擬。下圖是按時間順序基于預測校正不可壓縮SPH的粒子飛濺的傳統方法和改進方法的效果對比。圖1為8 000粒子數的飛濺效果對比,圖2為4 000粒子數的飛濺效果對比。

圖1 8 000粒子飛濺效果對比

圖2 4 000粒子飛濺效果對比
表1是在8 000個粒子的三維流體模擬實驗中,在GPU上進行了9個獨立實驗。

表1 在GPU上運行的實驗數據 fps
基于NVIDIA和Windows平臺,記錄實驗結果,如表1所示(fps表示程序運行幀率)。
圖3是根據表1中的數據繪制出的折線圖。其中實驗編號表示時間先后順序。從圖中可以看出,在文中規定的測試環境下,傳統方法和改進方法的幀率基本都處于45到60之間,但改進方法的幀率比傳統方法的幀率略小一點,但是兩種方法均達到實時性仿真效果的要求。

圖3 基于GPU的實驗程序性能比較
由上面效果比較可以看出,改進后的流體粒子飛濺較之前的傳統流體粒子飛濺來說有所變化,改進后的飛濺的粒子更有規律性,符合現實中在重力作用下的效果。流體仿真模擬即要求仿真模擬結果更接近現實,改進后的更符合。從表1與圖3可以看出模擬均達到實時仿真要求,但改進后的性能較之改進前的較差一點。
實驗核心部分即為校正密度壓力值,表2為傳統方法中的密度誤差和改進方法中的密度誤差(數據選擇來自傳統方法和改進方法飛濺的粒子數量相差較多的情況)。

表2 校正密度誤差值
圖4是由表2數據繪制出的折線圖,其中實驗編號表示運行時間先后順序。

圖4 校正的密度誤差值對比
雖然改進后的幀率校之傳統方法有所下降,但是基于PCISPH的流體仿真是通過預測校正密度壓力值從而得到仿真效果。實驗中通過改進密度誤差從而改進壓力值實現流體仿真,從表2和圖4中可以看出,改進后的密度誤差明顯小于傳統方法的密度誤差,并且基本上在千分之五左右。改進后密度誤差變小,增強了流體的不可壓縮性,減小了壓力噪聲,得到了更真實、更穩定的飛濺仿真效果。從而驗證了改進方法的有效性。
該文主要針對預測校正不可壓縮SPH(PCISPH)流體仿真進行粒子飛濺改進,針對改進減小密度誤差,增強流體的不可壓縮性,減小了壓力噪聲,從而得到了優于原模型的仿真效果。但是該模型也存在相應的問題,模型不夠穩定,該文添加較多粒子來減小不穩定帶給實驗結果的影響,在未來將對這些方面進行深入研究,以達到穩定且較大時間步長的流體模擬。同時,該文沒有對該模型進行渲染,未來將對模型進行渲染。