基于CUDA的弱可壓SPH流體建模與仿真

段興鋒,任鴻翔,神和龍

(1.大連海事大學航海動態(tài)仿真和控制交通行業(yè)重點實驗室,遼寧 大連 116026；2.集美大學航海學院,福建 廈門 361021)

1 引言

流體在自然界中普遍存在，如水、煙等,是具有高度形變的物理實體，其仿真在近年來已成為計算機圖形學研究的熱點。水體的仿真主要采用幾何模型、波浪譜和物理模型等方法進行建模，近年來研究的重點主要是基于物理模型的方法，而基于物理模型的方法主要以Navier-Stokes(N-S)方程為基礎，可分為基于網(wǎng)格的歐拉法和基于粒子的拉格朗日法[1]。基于粒子的拉格朗日方法適用于各種規(guī)模場景的流體模擬，并能實現(xiàn)翻卷等逼真效果，因此正越來越受到流體模擬研究者的關注。目前，拉格朗日法模擬流體主要有MPS(Moving Particle Semi-implicit)[2]、LBM(Lattice Boltzmann Model)[3]和SPH(Smoothed Particle Hydrodynamics)[4]等方法。

目前，SPH方法被廣泛應用于流體動力學等領域，且該方法更加真實、穩(wěn)定且高效，已成為計算機圖形學領域流體仿真的主流方法之一，具有較好的應用前景[5]。Müller等人[4]首次將SPH方法應用于水面的繪制，其后SPH在水面繪制得到了長足的發(fā)展；Premo?e等人[6]指出，由于SPH方法基于理想氣體狀態(tài)方程，可能會產(chǎn)生很明顯的流體體積壓縮，不適合模擬水體等不可壓縮流體。Becker等人[7]對SPH方法進行了改進，提出了弱可壓SPH WCSPH(Weekly Compressible SPH)方法，用Tait方程代替理想氣體狀態(tài)方程，在時間步長較小的條件下，得到了近似體積守恒的效果。如要得到不可壓縮的流體，就需要對復雜的壓力泊松方程求解，非常費時。Solenthaler等人[8]提出了PCISPH(Predictive-Corrective Incompressible SPH)方法，通過對壓力不斷地修正、迭代，利用預測出的粒子密度來迭代更新當前每個粒子上的壓力，實現(xiàn)不可壓縮性流體的模擬，耗費大量的計算資源，而且不適合GPU并行計算。GPU的發(fā)展經(jīng)歷了固定管線、可編程GPU、GPU通用計算、CUDA等多個階段。Harada等人[9]提出利用空間網(wǎng)格進行鄰域粒子的搜索，大大提高了搜索的效率，將時間復雜度降為O(nlogn)，其缺陷是每個網(wǎng)格最多只能包含4個粒子，且不能進行并行化。溫嬋娟等人[10]提出了一種完全基于GPU的SPH流體模擬方法，采用通用GPU GPGPU(General Purpose GPU)技術明顯提高了SPH流體模擬的速度，但是GPGPU程序均使用AMD公司的流計算開發(fā)工具實現(xiàn)，不便于后續(xù)對流體場景的控制。陳曦等人[11]則采用CUDA技術對SPH流體物理計算進行了加速，模擬出更為豐富的細節(jié)效果，但流體表面渲染使用了POV Ray作為繪制工具，不便于采用可編程流水管線實現(xiàn)流體渲染。Ji等人[12]用GPU對SPH計算進行了并行加速，并將該框架擴展至多GPU計算平臺，4個GPU協(xié)同并行計算860萬個粒子較單個GPU得到了3.96倍的加速比。

傳統(tǒng)的SPH方法用于模擬水體時不能保證水體的不可壓縮性，基于WCSPH對流體進行建模，相比較于PCISPH方法，能滿足實時性的要求。為了進一步提高計算效率，本文提出了CPU-GPU混合架構計算方法，能較大地提高流體計算效率，采用三維均勻網(wǎng)格、并行前綴求和和并行計算排序算法對鄰域粒子搜索算法進行優(yōu)化。

通過SPH物理計算得到粒子的密度、位置、速度等屬性值后，需要進行流體粒子的表面提取，并實現(xiàn)流體渲染，主要有Marching Cubes[13,14]、光線跟蹤[15]、屏幕空間繪制[16,17]等算法。其中，光線跟蹤算法是最經(jīng)典的三維數(shù)據(jù)體繪制方法之一，生成的圖像質(zhì)量較高，但大量光線與景物求交計算時間較長，隨著場景越來越復雜，其中的景物數(shù)量會增加，算法的計算量也不斷地增加，繪制一幅復雜場景往往耗時很長，根本滿足不了人們的交互性、實時性要求，只能做為一種離線的繪制方法。屏幕空間的方法是先將粒子透視到屏幕上獲取深度圖，然后采用后處理的方法對深度圖進行平滑處理，最后根據(jù)平滑后的深度圖求取法向量，從而渲染場景，提高了繪制效率，逼真度依賴于屏幕空間網(wǎng)格的分辨率，是一種以逼真度換速率的方法。Marching Cubes是一種最流行的三維表面提取算法，算法簡潔但計算量大。因此，本文為了提高繪制速度，采用GPU加速算法，尤其是流體表面的標量場的計算和空體素的剔除部分，最后利用環(huán)境貼圖表現(xiàn)流體的反射和折射效果，采用openGL著色語言GLSL(openGL Shading Language)實現(xiàn)了流體表面的著色。

2 基于CUDA的WCSPH的流體建模

2.1 SPH建模

SPH是一種無網(wǎng)格的拉格朗日粒子法，其核心是一種插值，引入光滑核函數(shù)表示粒子的影響域，通過影響域的積分插值計算粒子的密度、壓力、粘性力，從而得到位置、速度和加速度等場變量。對不可壓縮流體的N-S方程進行簡化，其計算式如下：

(1)

其中，ρ為流體的密度，u為流體的速度，p為流體壓強，f為外力，μ為流體粘度。

對于液體(如水體)而言，內(nèi)力可以分為壓力和粘性力，外力則包括重力、表面張力等。密度、壓力等宏觀變量均可通過SPH粒子近似法進行求解，該方法通過粒子近似將SPH的核近似有關的連續(xù)積分表達式離散化為由支持域內(nèi)粒子之和的形式：

,h)

(2)

其中,mj，ρj，rj分別代表粒子j的質(zhì)量、密度和位置；Aj表示粒子j的支持域內(nèi)粒子的物理量；W為光滑核函數(shù)。

Figure 1 CPU-GPU mixed architecture flow chart of SPH fluids modeling using CUDA圖1 基于CUDA的SPH流體模擬CPU-GPU混合架構流程圖

2.2 并行計算流程

為了能夠更好地模擬流體的細節(jié)，整個系統(tǒng)需要的粒子數(shù)至少幾萬個，且SPH流體粒子的位置、壓力、粘性力、加速度、速度等物理量的計算具有高度的并行性，而GPU具有強大的并行計算能力，其浮點計算能力已經(jīng)達到了P級。CUDA是NVIDIA專門為GPU通用計算設計的一種全新的架構，在CUDA架構下，開發(fā)人員可以通過CUDA C對GPU編程，從而避免了需要具備圖形學知識，通過OpenGL或者DirectX等API來訪問GPU的缺點。因此，采用CUDA對流體物理引擎計算，使得在一定數(shù)量的流體粒子的模擬能達到實時性，大大提高了程序的執(zhí)行效率。整個SPH流體模擬的CPU-GPU混合架構流程如圖1所示。

基于WCSPH的流體物理引擎計算和模擬的步驟如下：

(1)初始化并設置WCSPH、粒子系統(tǒng)參數(shù)和緩存。

(2)將有關參數(shù)和緩存數(shù)據(jù)轉(zhuǎn)移至GPU內(nèi)存。

(3)為每個粒子分配空間網(wǎng)格。

(4)對于每個粒子，進行三維空間鄰域粒子搜索，然后：

①計算壓強和密度；

②計算壓力和粘性力；

③計算其他外力,包括重力和表面張力；

④計算粒子的加速度和速度；

⑤計算粒子的位置；

⑥更新粒子的位置。

(5)直接采用GPU紋理緩存數(shù)據(jù)如位置、速度和顏色場進行三維可視化渲染。

(6)重復步驟(4)。

其中,SPH和粒子系統(tǒng)參數(shù)的初始化在主機端執(zhí)行，并傳送至設備端的常量存儲器中，利用常量存儲器的緩存機制加快訪問速度，最快訪問速度為一個GPU時鐘周期。而流體粒子的位置、速度和顏色則采用VBO的方式進行存儲，以便流體粒子繪制時直接訪問VBO緩沖區(qū)(位置、速度和顏色)，而不需要將數(shù)據(jù)從GPU傳輸至CPU端，然后再進行繪制，這將大大減少GPU和CPU之間的數(shù)據(jù)傳輸。另外，由于繪制流體粒子時不需要壓強、密度、壓力、粘性力等數(shù)據(jù)，因此可將這些數(shù)據(jù)直接存儲于GPU全局存儲器，但由于GPU全局存儲器無緩存功能，故訪問速度慢。

在CUDA實現(xiàn)SPH物理引擎計算時，首先通過調(diào)用cudaGraphicsGLRegisterBuffer，告訴CUDA運行時希望在OpenGL和CUDA中使用OpenGL VBO緩沖區(qū)，該函數(shù)返回一個指向VBO緩沖區(qū)的句柄。然后，在隨后對CUDA運行時的調(diào)用中，將通過這個句柄來訪問VBO緩沖區(qū)。粒子的插入、排序，粒子的壓強和密度、壓力和粘性力、表面張力的計算，粒子的更新均采用CUDA核函數(shù)在GPU端并行計算，實現(xiàn)時調(diào)用cudaGraphicsMapResources和cudaGraphicsResourceGetMappedPointer映射到VBO緩沖區(qū)，并獲得該緩沖區(qū)的GPU指針。

2.3 kernel函數(shù)

第一個kernel函數(shù)根據(jù)粒子的位置將粒子插入至三維空間網(wǎng)格中，構成粒子-網(wǎng)格對。第二個kernel函數(shù)為并行計數(shù)排序函數(shù)，使粒子-網(wǎng)格對按照網(wǎng)格編號進行排序，形成網(wǎng)格-粒子對。第三個kernel函數(shù)為密度和壓強的計算函數(shù)，將它們放到同一內(nèi)核中求解，以提高計算效率。根據(jù)SPH方法，當前粒子i的密度計算公式為：

,h)

(3)

壓強可以由當前粒子的密度直接計算得到，計算公式采用WCSPH方法中推薦的Tait方程，如下：

(4)

其中，k1，k2為常數(shù)項，k1=1119 kPa，k2=7。

第四個kernel函數(shù)為壓力和粘性力的計算函數(shù)，同樣也將它們放入同一內(nèi)核中求解，同理，由SPH方法，當前粒子i的壓力和粘性力的計算式為：

▽spikyW(ri-rj,h)

(5)

(6)

第五個kernel函數(shù)為粒子的更新，包括邊界處理、碰撞檢測、外力計算。求取各流體粒子內(nèi)力、外力后，便可求得最終的加速度ai，時間積分采用跳蛙法，更新速度ui和位置ri，其優(yōu)點是計算時所需要的存儲量低，每次計算中只需要進行一次優(yōu)化估值。

3 關鍵問題

3.1 弱可壓縮性

很多工作采用理想氣體狀態(tài)方程計算壓強，該方程中壓強和密度呈簡單線性關系，從而可能出現(xiàn)高度的可壓縮性，適用于氣體，而水體具有不可壓縮性。

目前，WCSPH和PCISPH提出了解決水體不可壓縮性問題的方法，其中PCISPH方法通過對壓力不斷地修正、迭代，利用預測出的粒子密度來迭代更新當前每個粒子上的壓力，實現(xiàn)不可壓縮性流體的模擬。

但是，大量的SPH粒子本身非常耗費內(nèi)存，需要保存位置、速度、壓力、壓強、粘性力等多個分量，在PCISPH方法中，還需要保存預測的部分物理量，且不適合GPU并行計算。WCSPH方法采用Tait方程計算壓強，該方程是一個經(jīng)驗方程，它僅僅需要很小的計算量，能滿足小尺度范圍流體的模擬需要，保證流體的弱可壓縮性。可見，WCSPH方法具有比解泊松方程和PCISPH方法高得多的時間效率，而且實現(xiàn)方法簡單，類似于基本SPH方法。其缺點在于壓強隨著密度的增大而持續(xù)增長，不符合物理規(guī)律，需要采用較小的時間步長，否則系統(tǒng)在運行一段時間之后穩(wěn)定性會大大降低，產(chǎn)生數(shù)值耗散，其原因在于密度的微小變化都會導致壓力值產(chǎn)生較大的波動。實驗表明，部分粒子會做無規(guī)則的隨機運動，甚至飛出系統(tǒng)邊界。本文采用類似于延訶等人[18]的解決方法優(yōu)化Tait方程，采用類似于Lennard-Jones方程的壓力狀態(tài)方程，解決了過高的指數(shù)項帶來的過大數(shù)值耗散，其公式如下：

(7)

其中,k1調(diào)整為1 958 kPa，其取值主要保證密度的波動率不超過1%。

3.2 人工粘度

在應用人工壓縮性的時候，采用修正的SPH方法XSPH(Corrected SPH)[19],粒子按以下形式運動：

(8)

其中，ε為常數(shù)，且ε∈[0,1]，用于控制速度校正的強度，在許多情況下，ε=0.3是模擬不可壓縮流體的最佳選擇。通過人工增加動量方程中的粘性力，使得中心粒子以更接近于相鄰粒子的平均速度運動，因而粒子運動變得更加有序，顯得相當整齊。

3.3 并行的鄰域粒子搜索

SPH方法的實質(zhì)是從大量離散粒子的運動和相互作用中重構連續(xù)介質(zhì)，粒子間是否發(fā)生相互作用的前提是判斷它們是否鄰近。每個粒子一般有30～40個鄰域粒子，假定流體物理引擎中有100 000個粒子，那么每幀就有300 000～400 000個粒子參與計算。可見，鄰域粒子搜索算法對于程序的效率顯得尤為重要。本文首先采用三維空間網(wǎng)格對整個模擬區(qū)域進行均勻網(wǎng)格劃分，并將粒子插入至網(wǎng)格中，如圖2所示，此時粒子緩沖區(qū)是按粒子編號排序的；然后采用并行前綴求和算法求得網(wǎng)格中的粒子數(shù)量和偏移地址；接著采用并行計數(shù)排序算法按網(wǎng)格編號排列粒子緩沖區(qū)，如表1所示；最后在查找最近鄰粒子時僅查找鄰近網(wǎng)格內(nèi)的粒子，若近鄰粒子與計算粒子的距離小于光滑核半徑，則進行求和計算，否則，該近鄰粒子不參與物理計算。

Figure 2 3D uniform grid圖2 三維均勻網(wǎng)格

排序前粒子編號網(wǎng)格編號排序后粒子編號網(wǎng)格編號偏移地址粒子個數(shù)03072001122421112321222133103031421331435315316328317202327283163295110409110409511011152115211112531253121

鄰域粒子搜索算法采用三維空間網(wǎng)格對粒子進行網(wǎng)格劃分后，查找近鄰粒子時僅搜索粒子所在網(wǎng)格鄰域內(nèi)的27個網(wǎng)格，從全配對搜索時間復雜度O(n3)降低至O(n2)，主要步驟如下：

(1)對三維空間進行均勻網(wǎng)格劃分，并根據(jù)每個網(wǎng)格在空間中的位置給每個網(wǎng)格分配一個整型變量作為網(wǎng)格編號，分配若干網(wǎng)格緩沖區(qū)，包括粒子數(shù)、偏移地址。

(2)分配若干粒子緩沖區(qū)，包括位置、速度、顏色等，并按粒子的空間位置插入至網(wǎng)格中，建立粒子編號與網(wǎng)格編號之間的關系。

(3)采用并行前綴求和算法求得每個網(wǎng)格中的粒子數(shù)量和偏移地址，接著采用并行計數(shù)排序算法按網(wǎng)格編號排列粒子緩沖區(qū)。其中，并行計數(shù)排序算法和前綴求和算法均采用CUDA并行加速策略，以提高程序的性能。

(4)求取計算粒子的鄰近27個網(wǎng)格單元，對于網(wǎng)格單元中的每個近鄰粒子，判斷其是否在支持域內(nèi)，若近鄰粒子與計算粒子的距離小于光滑核半徑，則進行求和計算；否則，該近鄰粒子不參與物理計算。

3.4 邊界處理

SPH方法處理流體的自由表面時，其邊界條件是自由滿足的，不需要特殊處理。但是，在流固耦合的邊界位置，對于邊界處的流體粒子，由于核函數(shù)積分域被邊界截斷，支持域內(nèi)粒子數(shù)不足，計算的物理屬性將帶來一定的誤差。固壁邊界排斥力模型主要有Lennard-Jones與法向力模型。其中Lennard-Jones模型將固壁邊界離散成一排“邊界粒子”，當流體粒子運動到其影響范圍之內(nèi)時，沿著兩粒子的中心線對水粒子施加一個強排斥力，而法向力模型則認為邊界粒子對流體粒子只提供法向力作用，沒有切向力，本文采用法向力模型[20]。該法向力的表達式為：

fij=njR(y)P(x)

(9)

其中,nj為邊界粒子j的法向量，y為水粒子i到邊界粒子j的法向距離，x為切向距離。函數(shù)R(y)的表達式為：

(10)

函數(shù)P(x)可保證粒子平行于邊界運動時，排斥力保持不變，其表達式為：

(11)

其中,q=y/(2△p)，△p為粒子的初始間距，當q<1時，邊界粒子對水粒子施加排斥力。

4 SPH流體的繪制

本文采用基于CUDA并行加速的Marching Cubes算法實現(xiàn)流體表面提取(Isosurface Extraction)，利用環(huán)境貼圖表現(xiàn)流體的反射和折射效果，以實現(xiàn)流體表面的著色，采用折射效果可較好地提高場景繪制的真實感效果。Marching Cubes流體表面提取算法主要涉及網(wǎng)格節(jié)點處的密度計算、體素分類、交點坐標和法矢量計算。但是，Marching Cubes算法需耗費大量的計算資源，該算法最為耗時的操作在于數(shù)據(jù)處理階段，特別是網(wǎng)格節(jié)點密度的計算和大量的空體素占據(jù)了計算資源。因此，本文采用CUDA并行計算網(wǎng)格節(jié)點處的密度值和粒子網(wǎng)格節(jié)點密度值，計算思想由影響域改為支持域，即當前網(wǎng)格節(jié)點的密度值由鄰域范圍內(nèi)所有粒子計算所得，而不是計算一個粒子點對鄰域范圍內(nèi)產(chǎn)生的影響，如圖3所示，這將導致不同的密度計算方法，圖3a為基于支持域思想的計算方法，有利于CUDA并行、圖3b為基于影響域思想的計算方法，不利于CUDA并行。

Figure 3 Two methods for calculating scalar field of the fluid surface圖3 兩種不同的流體表面的標量場計算方法

然后，本文改進了Marching Cubes算法，對體素分類后進行體素壓縮，剔除空體素，并采用CUDA并行加速。顯然，經(jīng)CUDA并行技術后，相對于CPU版本的Marching Cubes表面提取方式，該算法可以大幅節(jié)省數(shù)據(jù)處理時間，在流體粒子數(shù)量較多的情況下，效率提升更加明顯。

5 實驗結果

本文采用基于CUDA的SPH方法對流體進行建模，CUDA并行加速的Marching Cubes算法實現(xiàn)了流體的表面提取，利用環(huán)境貼圖表現(xiàn)流體的反射和折射效果，以實現(xiàn)流體表面的著色，實現(xiàn)了水的波動和木塊在水中晃動的效果，如圖4～圖6所示。其中，圖4為用點精靈渲染的水波動，圖5為本文方法渲染的水波動，圖6為本文方法渲染的木塊在水中晃動的效果圖。可見，本文方法較好地實現(xiàn)了流體的小尺度模擬，如水的波動、翻卷以及流固交互的效果。

Figure 4 Water wave rendering by point sprite圖4 用點精靈渲染的水波動圖

Figure 5 Simulation of water wave圖5 水的波動模擬效果圖

Figure 6 Simulation of wood shaking in the water圖6 木塊在水中搖晃模擬效果圖

實驗測試平臺CPU：Intel(R) Core(TM) i5-6400T CPU 2.2 GHz，內(nèi)存：8 GB，顯卡：NVIDIA GeForce GTX960，硬盤：1 TB。SPH計算中鄰域粒子搜索算法從全配對搜索算法改為基于網(wǎng)格的鄰域粒子搜索后，時間復雜度由O(n3)降低至O(n2)，在CPU平臺上，經(jīng)測試，插入粒子并排序、壓強和密度計算、壓力和粘性力計算和粒子更新所需要的時間詳見表2。可見采用三維空間網(wǎng)格后，除了粒子插入至網(wǎng)格并排序的算法時間耗費增大外，其他壓強、密度、壓力和粘性力的計算效率均得到了較大的增高，但仍不能滿足實時需要。

Table 2 SPH physical computing timeof different search methods

為此，我們采用CUDA架構對鄰域粒子搜索算法進行了GPU并行加速，對65 536、131 072、262 144、524 288和1 048 576個粒子的情況分別進行了CPU和GPU版本的程序測試，物理計算時間和加速比如表3所示。其中CPU和GPU版本物理計算的鄰域粒子搜索算法分別為鄰近網(wǎng)格的粒子搜索法和CUDA并行鄰域粒子搜索法。顯然，采用CUDA架構的GPU并行計算大大節(jié)省了時間，加速比最大達到60.7。

隨著粒子數(shù)目的增加，本文采用的基于CUDA的SPH流體計算能夠帶來更大的加速比，在100萬個粒子左右，流體仿真系統(tǒng)的物理計算能達到實時的要求。

Table 3 SPH physical computing timewith different numbers of particles

6 結束語

SPH、WCSPH、PCISPH等方法在流體模擬中得到了廣泛的應用，考慮到水體的不可壓縮性以及實時性的要求，采用WCSPH方法實現(xiàn)流體的建模，并用CUDA進行了GPU加速和渲染。該方法思想簡單，易于實現(xiàn)，能滿足小尺度流體場景實時模擬，實現(xiàn)了水的波動和木塊在水中晃動的效果。從實驗結果可知，對于100萬個粒子的WCSPH計算能滿足實時要求；采用CUDA并行加速后的Marching Cubes算法進行表面提取和繪制，連同物理計算也能滿足20萬個粒子實時計算和渲染的要求，大幅度提高流體物理計算和渲染的整體性能，但是Marching Cubes算法用于流體繪制時不能展現(xiàn)泡沫、浪花飛濺等特效。我們下一步的工作將實現(xiàn)破碎浪和流固交互的模擬，并對SPH物理計算進一步優(yōu)化，采用GPU集群系統(tǒng)實現(xiàn)大尺度復雜流體場景的模擬。