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

段興鋒,任鴻翔,神和龍

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

1 引言

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

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

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

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

2 基于CUDA的WCSPH的流體建模

2.1 SPH建模

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

(1)

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

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

,h)

(2)

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

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

2.2 并行計算流程

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

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

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

(2)將有關參數和緩存數據轉移至GPU內存。

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

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

①計算壓強和密度；

②計算壓力和粘性力；

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

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

⑤計算粒子的位置；

⑥更新粒子的位置。

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

(6)重復步驟(4)。

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

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

2.3 kernel函數

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

,h)

(3)

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

(4)

其中，k1，k2為常數項，k1=1119 kPa，k2=7。

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

▽spikyW(ri-rj,h)

(5)

(6)

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

3 關鍵問題

3.1 弱可壓縮性

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

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

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

(7)

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

3.2 人工粘度

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

(8)

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

3.3 并行的鄰域粒子搜索

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

Figure 2 3D uniform grid圖2 三維均勻網格

排序前粒子編號網格編號排序后粒子編號網格編號偏移地址粒子個數03072001122421112321222133103031421331435315316328317202327283163295110409110409511011152115211112531253121

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

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

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

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

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

3.4 邊界處理

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

fij=njR(y)P(x)

(9)

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

(10)

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

(11)

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

4 SPH流體的繪制

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

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

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

5 實驗結果

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

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

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版本物理計算的鄰域粒子搜索算法分別為鄰近網格的粒子搜索法和CUDA并行鄰域粒子搜索法。顯然，采用CUDA架構的GPU并行計算大大節省了時間，加速比最大達到60.7。

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

Table 3 SPH physical computing timewith different numbers of particles

6 結束語

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