基于SPH的水體交互三維虛擬仿真

方貴盛，鄭高安

(浙江水利水電學(xué)院 機械與汽車工程學(xué)院，浙江 杭州 310018)

基于SPH的水體交互三維虛擬仿真

方貴盛，鄭高安

(浙江水利水電學(xué)院 機械與汽車工程學(xué)院，浙江 杭州 310018)

水體與復(fù)雜邊界物體相互作用時的運動情況非常復(fù)雜.目前用基于網(wǎng)格的方法很難準(zhǔn)確有效地實現(xiàn)三維模擬.采用基于無網(wǎng)格的光滑粒子動力學(xué)方法，通過求解納維—斯托克斯方程，在VC++和OpenGL編程環(huán)境下，實時動態(tài)模擬了水體相互碰撞、水流通過圓形橋墩、水流通過丁壩時的三維動作效果.通過與實際水流模擬實驗觀察效果相比較，驗證了所實現(xiàn)算法的有效性，說明SPH方法可廣泛應(yīng)用于水體復(fù)雜運動過程模擬，并可用于解決工程實際問題.

納維—斯托克斯方程；光滑粒子動力學(xué)；水體；虛擬仿真

水作為自然界中一種常見的物質(zhì)，對其研究和仿真具有重要的理論意義和實用價值.由于水在形態(tài)上是不規(guī)則的，同時也是動態(tài)變化的，因此水體的建模與仿真一直是計算機圖形學(xué)中的一個研究難點，同時也是真實感圖形學(xué)領(lǐng)域的一個研究熱點[1-2].

水體虛擬仿真建模方法包括了基于幾何曲面的方法[3]、基于波形函數(shù)的方法[4]、基于海浪譜的方法[5]、基于粒子系統(tǒng)的方法[6]、基于物理模型的方法[7-9]等.近年來，隨著計算機軟硬件技術(shù)的發(fā)展，基于流體動力學(xué)的方法被廣泛應(yīng)用于流體的仿真，取得了良好的效果.本文采用了光滑粒子動力學(xué)方法，通過求解納維—斯托克斯方程，在VC++和OpenGL編程環(huán)境下，動態(tài)模擬了水體相互碰撞、水流通過圓形橋墩、水流通過丁壩時的動作效果.

1 基于SPH的水流虛擬仿真基本原理

1.1 水流控制方程

描述水流運動的控制方程有圣維南方程、二維淺水動力方程、歐拉方程、納維—斯托克斯方程等.為了更好地模擬水流的運行，本文采用了納維—斯托克斯方程組(Navier-Stokes Equations，簡稱NSE)，它由一系列的偏微分方程來表示，包含了質(zhì)量守恒方程和動量守恒方程.

質(zhì)量守恒方程：

▽·U=0

(1)

動量守恒方程：

(2)

式中：U—速度矢量；▽—向量梯度算子；

▽2—拉普拉斯算子；t—時間變量；

p—流體的壓強；ρ—流體密度；

μ—流體動力粘滯系數(shù)；f—單位質(zhì)量力.

1.2 SPH離散化求解方法

光滑粒子動力學(xué)方法(Smoothed Particle Hydrodynamics，簡稱SPH)是近二十多年來逐步發(fā)展起來的一種無網(wǎng)格方法.SPH的基本思想是將連續(xù)的流體用相互作用的粒子組來描述，每個粒子點上攜帶獨立的物理信息，包括質(zhì)量、速度、密度等.通過求解粒子組的動力學(xué)方程和跟蹤每個粒子的運動軌跡，便可求得整個流體系統(tǒng)的力學(xué)行為.采用SPH方法求解NSE方程的關(guān)鍵是采用兩次逼近：一是核函數(shù)逼近，即將描述場的函數(shù)近似表示為該函數(shù)和核函數(shù)的乘積的積分表示；二是粒子逼近，即將場函數(shù)的積分表示用一系列粒子離散化來描述，從而降低了NSE方程求解的難度.根據(jù)上述思想，水流粒子的物理量，如密度、壓力和速度等，其SPH求解公式分別表示如下：

(3)

(4)

(5)

式中：mi，mj—粒子i，j的質(zhì)量；

pi，pj—粒子i，j所受的壓力；

ρi，ρj—粒子i，j的密度；

h—光滑核半徑.

2 基于SPH的水流虛擬仿真算法實現(xiàn)

2.1 算法的數(shù)據(jù)結(jié)構(gòu)設(shè)計

通過構(gòu)造粒子類和粒子系統(tǒng)類來確定算法的數(shù)據(jù)結(jié)構(gòu).單個粒子的屬性包括了粒子的位置、速度、預(yù)估速度、顏色、所受的壓力、密度、力(加速度)等參數(shù).粒子系統(tǒng)的屬性，包括了容納粒子系統(tǒng)的水箱容器大小、水體的大小、光滑核長度、時間步長、粒子的質(zhì)量、粘度、靜態(tài)密度、粒子半徑、粒子間距、光滑核半徑、內(nèi)部硬度、外部硬度、阻尼系數(shù)等參數(shù).

2.2 算法流程設(shè)計

采用SPH方法對NSE進行求解，以模擬水體動態(tài)流動過程，其總體設(shè)計思路(見圖1).

圖1 算法的總體設(shè)計思路

2.2.1 系統(tǒng)初始化

系統(tǒng)初始化過程需要構(gòu)造水體仿真場景，設(shè)置初始參數(shù)，包括水容器的參數(shù)、初始水體的參數(shù)、基本粒子參數(shù)、SPH參數(shù)，以及水粒子在容器中的排列方式等.為了提高粒子的最近鄰粒子的搜索速度，本文用了空間網(wǎng)格結(jié)構(gòu)，網(wǎng)格初始化構(gòu)建算法步驟如下.

步驟1：設(shè)置網(wǎng)格單元的大小(初始值為核半徑的兩倍)；

步驟2：確定網(wǎng)格的分辨率，即確定水箱容器X、Y、Z軸各排列有多少個網(wǎng)格單元；

步驟3：計算水容器總共包含有多少個網(wǎng)格單元.

步驟4：為每個網(wǎng)格單元分配指針；

步驟5：設(shè)置每個網(wǎng)格單元內(nèi)粒子的數(shù)目，初始值為0.

步驟6：繪制空間網(wǎng)格.

網(wǎng)格構(gòu)建完成以后，需要將一個個粒子加入到對應(yīng)的網(wǎng)格中，然后求出每個粒子所在網(wǎng)格的單元索引號，其實現(xiàn)算法如下所示：

步驟1：依次取出一個個水粒子i；

步驟2：確定粒子i的位置，并與空間網(wǎng)格的最小值進行比較，以確定粒子i所在的網(wǎng)格單元索引號；

步驟3：將粒子i插入到對應(yīng)的網(wǎng)格單元中；

步驟4：粒子i所在的網(wǎng)格單元索引中粒子數(shù)加1.

2.2.2 最近鄰粒子搜索

常用的最近鄰粒子搜索方法有直接粒子搜索技術(shù)、關(guān)聯(lián)鏈表搜索技術(shù)、樹形結(jié)構(gòu)搜索技術(shù)等.本文采用空間網(wǎng)格結(jié)構(gòu)和鏈表結(jié)構(gòu)來提高最近鄰粒子的搜索速度，其算法流程如下所示.

步驟1：取出粒子i，確定粒子i的位置；

步驟2：確定水容器X、Y、Z軸各排列有多少個網(wǎng)格單元；

步驟3：計算最小網(wǎng)格單元的位置；

步驟4：計算粒子i所在的網(wǎng)格單元索引號；

步驟5：按X、Y、Z軸正方向依次計算其余七個鄰居粒子網(wǎng)格單元索引號.

2.2.3 核函數(shù)選擇與計算

在SPH計算過程中，核函數(shù)的選擇至關(guān)重要，它不僅決定了函數(shù)近似式的形式，而且還決定核近似和粒子近似的一致性和精度.核函數(shù)一般應(yīng)包括歸一性、緊支性、對稱性、非負(fù)性等.常用的核函數(shù)有高斯型核函數(shù)、B樣條核函數(shù)、四次樣條核函數(shù)等，本文選用了Matthias Muller所采用的三種光滑核函數(shù)，分別用來求解粒子的密度、壓力和粘滯力等[10].

2.2.4 粒子的密度和所受的壓力計算

本文利用空間網(wǎng)格結(jié)構(gòu)來計算粒子所受的壓強力，其算法流程如下所示：

步驟1：依次取出一個個水粒子i，采用最近鄰粒子搜索方法，查找它的鄰居粒子所在的八個網(wǎng)格單元號；

步驟2：從鄰居網(wǎng)格單元中取出一個個粒子j，如果取出的粒子是i粒子，則取下一個粒子，否則計算粒子i與j之間的距離；

步驟3：如果粒子i與j之間的距離小于光滑核半徑，則計算式子(6)的值；

(6)

步驟4：判斷粒子i的鄰居粒子數(shù)是否小于40，如果滿足條件，則創(chuàng)建鄰居對列表，并計算鄰居粒子對的間距；

步驟5：根據(jù)公式(7)計算粒子i的密度；

(7)

步驟6：根據(jù)公式(8)計算粒子i所受的壓強力，式中K為常數(shù)，ρ0為初始密度；

pi=K(ρi-ρ0)

(8)

步驟7：根據(jù)壓強力的大小設(shè)置粒子的顏色值.

2.2.5 粒子所受的合力和加速度計算

本文采用空間網(wǎng)格結(jié)構(gòu)和鄰居列表結(jié)構(gòu)計算粒子的加速度，其算法流程如下所示：

步驟1：依次取出一個個水粒子i；

步驟2：從水粒子i的鄰居列表對中依次取出水粒子i的一個個鄰居粒子j；

步驟3：計算水粒子i與鄰居粒子j間的距離；

步驟4：根據(jù)公式(9)計算粒子i由于所受的壓力差而引起的加速度；

(9)

步驟5：根據(jù)公式(10)計算粒子i由于所受的粘滯力而引起的加速度.

(10)

2.2.6 碰撞檢測處理，粒子的加速度更新

在邊界處，粒子除了受到重力外，還受到固壁邊界的作用，此時粒子將改變原有的運動方向、速度和加速度.為了防止粒子穿越固壁邊界，本文采用了Matthias Muller所介紹的施加邊界力的方法[10].如果粒子與邊界間的距離小于2r(r為粒子的半徑)，則更新粒子的運動方向和加速度.

2.2.7 粒子的密度、速度和位置更新

SPH的時間積分方法，廣泛采用蛙跳法、預(yù)測校正法、龍格庫塔法等，均具有二階時間精度.由于在實際運用過程中，蛙跳法對存儲的需求量低，而且計算效率高，因此本文選用蛙跳積分方法來計算粒子演進的速度和位置.

2.2.8 可視化結(jié)果輸出

在每個計算步結(jié)束后，可以得到每個粒子的密度、加速度、速度、位置和所受的壓力等參數(shù).根據(jù)這些數(shù)據(jù)，用戶可以觀察到任何時刻的粒子位置分布、速度分布和壓力分布等，可以輸出不同時刻的速度、壓力等值線圖，也可以直觀地看到水體的產(chǎn)生、演進和變形的全過程，這對于評價或修改工程設(shè)計方案有很大幫助.

3 水體仿真實例

仿真實驗在Windows 7平臺上進行，采用VC++.net和OpenGL編程工具分別實現(xiàn)了兩側(cè)水體的相互碰撞效果仿真、水流通過丁壩時的動態(tài)效果仿真、水流通過圓形橋墩時的動態(tài)效果仿真等.實驗的硬件環(huán)境為Intel酷睿八核處理器，CPU 2.4 GHz，4 G內(nèi)存，1 G獨立顯存.

3.1 水體碰撞效果仿真

水箱容器長0.56 m，寬0.16 m，高0.32 m；左側(cè)水體長0.12 m，寬0.16 m，高0.045 m，粒子數(shù)為3 059個；右側(cè)水體長0.12 m，寬0.16 m，高0.097 m，粒子數(shù)為6 555個.粒子初始間距為0.006 m，按長方體規(guī)律分布，動力粘滯系數(shù)為0.2 kg/(m·s)，計算時間步長為0.004 s.初始時刻，兩側(cè)水體在重力作用下沿著壁面邊界流動(見圖2).

圖2 水體碰撞仿真效果

從圖2可以看出，初始時刻，兩側(cè)水體在重力的作用下流動，由于受水箱側(cè)壁的影響，左側(cè)水體向右流動，右側(cè)水體向左流動.由于右側(cè)水體比左側(cè)水體大，右側(cè)水體的流動速度比左側(cè)水體快，經(jīng)過了0.232 s后兩側(cè)水體相互碰撞，壓力增大.左側(cè)水體受右側(cè)水體的沖擊，水體流動方向發(fā)生改變，向上躍起并折回，同時濺起水花.接著兩側(cè)水體融合在一起，繼續(xù)向左側(cè)流動，當(dāng)碰到左側(cè)壁時，左側(cè)的水體升高，右側(cè)降低，然后又在重力和慣性力的作用下，向右側(cè)流動，當(dāng)碰到右側(cè)壁時，又向左側(cè)移動，如此反復(fù)，經(jīng)過5 s后水體趨向平衡靜止.圖2(7)、(8)中反映了兩側(cè)水體的流動速度狀況.剛開始時，水平流動速度較小，然后逐漸增大，當(dāng)兩側(cè)水體發(fā)生碰撞時，又有一個速度變小的過程，與實際觀察的結(jié)果相吻合.

3.2 水流通過丁壩時的動態(tài)效果仿真

水箱容器長0.8 m，寬0.16 m，高0.32 m；梯形丁壩的上底為0.04 m，下底為0.08 m，高度為0.024 m，長度為0.08 m，丁壩最左側(cè)位置與水箱左側(cè)壁間的距離為0.4 m；水體長0.16 m，寬0.16 m，高0.1 m，粒子數(shù)為9 928個(見圖3).粒子初始間距為0.006 m，按長方體規(guī)律分布，動力粘滯系數(shù)為0.2 kg/(m·s)，計算時間步長為0.004 s.

圖3 水流通過丁壩時的動態(tài)效果仿真

初始時刻，水體在重力作用下沿著壁面邊界流動.當(dāng)水流到達梯形丁壩時，由于受到丁壩坡面的影響，一部分水體沿著坡面向上運動并越過丁壩，一部分水體受丁壩的阻擋被折回，還有一部分水體繞過丁壩向前推進.當(dāng)水流到達水箱右側(cè)壁時，在壁面的影響下，水體被折回，然后向左流動，并再次通過丁壩.由于水流通過丁壩時，速度不一致，在右側(cè)壁的影響下，將會產(chǎn)生渦流效果.圖3(1)～(6)反映的是水流通過梯形丁壩時的位置及壓力分布效果，圖3(7—8)為水流速度分布狀況，青色為低速部分粒子，黃色為高速部分粒子.

3.3 水流通過圓形橋墩時的動態(tài)效果仿真

水箱容器長0.8 m，寬0.16 m，高0.32 m；橋墩半徑為0.032 m，高度為0.32 m，橋墩軸心與水箱左側(cè)壁間的距離為0.48 m；水體長0.16 m，寬0.16 m，高0.1 m，粒子數(shù)為9 936個(見圖4).粒子初始間距為0.006 m，按長方體規(guī)律分布，動力粘滯系數(shù)為0.2 kg/m·s，計算時間步長為0.004 s.

圖4 水流通過圓形橋墩時的動態(tài)效果仿真

初始時刻，水體在重力作用下沿著壁面邊界流動，速度逐漸增大.當(dāng)水流到達圓形橋墩時，由于受到橋墩的阻力，靠近橋墩表面水粒子壓力增大，速度減慢，水位升高.橋墩兩側(cè)的水粒子兵分兩路繞過橋墩向前推進，速度增大，當(dāng)達到水箱右側(cè)壁時，由于受到阻力的影響，水流向后折回，同時濺起水花.當(dāng)5 s后，水位趨于平穩(wěn).圖4(1)～(7)反映的是水流通過圓形橋墩時的位置及壓力效果，圖4(8)、(9)為水流速度分布效果，青色為低速部分粒子，黃色為高速部分粒子.

4 結(jié) 論

基于無網(wǎng)格的SPH水體運動仿真算法，能夠解決基于網(wǎng)格的水體模擬方法中存在的單元劃分、網(wǎng)格纏結(jié)和扭曲等問題，在模擬水體大變形及波浪破碎等方面具有獨特的優(yōu)勢，因此得到廣泛的應(yīng)用.本文采用基于光滑粒子動力學(xué)方法，通過求解納維—斯托克斯方程，在VC++和OpenGL編程環(huán)境下，實時動態(tài)地模擬了水體相互碰撞、水流通過圓形橋墩、水流通過丁壩時的三維動作效果，模擬結(jié)果比較理想，與實驗效果吻合較好.對于更復(fù)雜的水體交互問題，有待于進一步研究.

[1] 方貴盛,潘志庚.水體虛擬仿真與應(yīng)用綜述[J].計算機仿真，2012，29(10):30-33,361.

[2] 方貴盛，潘志庚.水體虛擬仿真建模技術(shù)研究進展[J].系統(tǒng)仿真學(xué)報，2013，26(9):1981-1989.

[3] JAY Allen FAULKNER. Beauty Waves: an Artistic Representation of Ocean Waves Using Bezier Curves [D]. Texas: Texas A & M University,2006.

[4] 趙 欣，李鳳霞，戰(zhàn)守義.基于小振幅波理論的淺海波浪建模及動態(tài)仿真[J].系統(tǒng)仿真學(xué)報,2008,20(2):281-284.

[5] 楊懷平，孫家廣.基于海浪譜的波浪模擬[J].系統(tǒng)仿真學(xué)報，2002，14(9):1175-1178.

[6] 張尚弘，陳 壘，趙登峰.基于粒子系統(tǒng)的流場實時模擬[J].水利水電技術(shù)，2004，35(9):47-50.

[7] TSUNEMI TAKAHASHI,HIROKO FUJII,ATSUSHI KUNIMATSU. Realistic Animation of Fluid with Splash and Foam [J]. Computer Graphics,2002,21(3):736-744.

[8] MARKUS IHMSEN,NADIRAKINCI,MARC GISSLER. Boundary Handling and Adaptive Time-Stepping for PCISPH[C]]// Kenny Erleben. Proceedings of the 7thWorkshop on Virtual Reality Interaction and Physical Simulation,Aire-la-Ville: The Eurographics Association,2010:79-88.

[9] NADIR AKINCI,MARKUS IHMSEN,GIZEM AKINCI. Versatile Rigid-Fluid Coupling for Incompressible SPH [J],Acm Tranctions on Graph,2012,31(4):1-62.

[10] MATTHIAS MULLER,DAVID CHARYPAR,MARKUS GROSS. Particle-based Fluid Simulation for Interactive Applications[C]]//D. Breen. Proceedings of the 2003 ACM SIGGRAPH/Eurographics symposium on Computer animation,Aire-la-Ville: The Eurographics Association,2003:154-159.

3DVirtualSimulationofWaterInteractionBasedonSPH

FANG Gui-sheng,ZHENG Gao-an

(College of Mechanical and Automotive Engineering,Zhejiang University of Water Resources and Electric Power,Hangzhou 310018,China)

When water interacts with complex boundary objects,the movement is very complicated. It is difficult to realize 3D simulation of this phenomenon by using the grid-based method. In this paper,a real-time virtual simulation method based on SPH is proposed. In the programming environment of VC++ and OpenGL,the Navier-Stokes equation is solved to simulate some complex water phenomena,such as interaction between water and water,interaction between water and pier,and interaction between water and spur dike. Experimental results show that the method is realistic and can be used to solve some project problems.

Navier-Stokes Equations; smoothed particle hydrodynamics; water; virtual simulation

2016-03-18

2013年度浙江省水利廳科研基金資助項目(RC1337)

方貴盛(1973-)，男，江西婺源人，博士，教授，主要研究方向為人機交互與虛擬仿真技術(shù).

TP391.9

A

1008-536X(2016)10-0068-06