張靜，劉文龍，盧文虎

（1.天津市博盈科技發(fā)展有限公司天津市300180；2.中國(guó)地震局第一監(jiān)測(cè)中心天津市300180；3.國(guó)家海洋信息中心天津市300171）

基于SPH方法的海浪動(dòng)態(tài)演變模擬仿真

張靜1，劉文龍2，盧文虎3

（1.天津市博盈科技發(fā)展有限公司天津市300180；2.中國(guó)地震局第一監(jiān)測(cè)中心天津市300180；3.國(guó)家海洋信息中心天津市300171）

結(jié)合風(fēng)浪場(chǎng)復(fù)雜的動(dòng)力學(xué)、時(shí)空特性，針對(duì)海浪仿真模擬中存在的真實(shí)感不足，難以模擬海浪破碎等問(wèn)題，根據(jù)流體的物理特性，采用光滑粒子流體動(dòng)力學(xué)（SPH）的方法，從粒子初始化分配、粒子受力分析及狀態(tài)計(jì)算以及粒子運(yùn)動(dòng)過(guò)程中的碰撞檢測(cè)3個(gè)方面進(jìn)行了研究，實(shí)現(xiàn)了海浪的動(dòng)態(tài)演變仿真。

海浪模擬；粒子系統(tǒng)；光滑粒子流體動(dòng)力學(xué)；碰撞檢測(cè)

海浪仿真是近年來(lái)計(jì)算機(jī)圖形學(xué)領(lǐng)域研究的熱點(diǎn)之一，它在虛擬現(xiàn)實(shí)應(yīng)用、計(jì)算機(jī)游戲以及電影制作中起著重要的作用。另外，海浪的仿真模擬是構(gòu)建海洋環(huán)境信息可視化平臺(tái)的一個(gè)重要組成部分，平臺(tái)的建立可顯著提高我國(guó)海洋信息的可視化水平。

海浪具有復(fù)雜的動(dòng)力學(xué)、時(shí)空特性，對(duì)其三維模擬有著特殊的要求，這使得海浪建模及繪制成為海浪仿真模擬中的難點(diǎn)。海浪建模的本質(zhì)是海浪的真實(shí)感繪制，它是海洋環(huán)境仿真的基礎(chǔ)，在很大程度上決定了仿真效果的好壞[1]。針對(duì)海浪復(fù)雜的特性，本文采用一種能夠處理自由表面流動(dòng)的基于粒子的拉格朗日法（SPH方法），將粒子系統(tǒng)與物理模型結(jié)合進(jìn)行海浪的模擬。

1 SPH基本原理

SPH（Smoothed Particle Hydrodynamics，光滑粒子流體動(dòng)力學(xué)）方法是一種由一組粒子代替流體獲得流體動(dòng)力學(xué)公式的近似數(shù)學(xué)解決方案的方法[2]。它的基本思想是將連粒子質(zhì)點(diǎn)上承載著各種物理量，通過(guò)求解質(zhì)點(diǎn)組的動(dòng)力學(xué)方程和跟蹤每個(gè)粒子的運(yùn)動(dòng)軌跡，求得整個(gè)系統(tǒng)的力學(xué)行為。對(duì)于每個(gè)單獨(dú)的流體粒子，仍然遵循最基本的牛頓第二定律，由于在SPH算法中，流體的密度決定了流體單元的質(zhì)量，所以這里遵循牛頓第二定律公式可以表示為公式1.1的形式。

圖1 光滑核函數(shù)影影響域內(nèi)的粒子示意圖

SPH方法是一種無(wú)網(wǎng)格拉格朗日算法，它通過(guò)一系列粒子質(zhì)點(diǎn)的“核函數(shù)估值”將流體力學(xué)基本方程組轉(zhuǎn)換成數(shù)值計(jì)算用的公式。其核心是核函數(shù)，它表示在一定的光滑長(zhǎng)度范圍內(nèi)其臨近粒子質(zhì)點(diǎn)對(duì)研究粒子影響程度的權(quán)函數(shù)。假設(shè)流體中某點(diǎn)r，在光滑核半徑h范圍內(nèi)受數(shù)個(gè)粒子影響，其位置分別是r0，r1，r2…，rj，如圖1所示，該點(diǎn)位置處某項(xiàng)屬性A的累加值就可以用公式1.2來(lái)表示。

其中Aj是要累加的某種屬性，mj和ρj是周圍粒子的質(zhì)量和密度，rj是該粒子的位置，h是光滑核半徑，函數(shù)W為光滑核函數(shù)。

由于粒子質(zhì)點(diǎn)之間不存在網(wǎng)格關(guān)系，它可以避免因變形引起的網(wǎng)格扭曲從而造成精度破壞等問(wèn)題，而且也能方便地處理不同介質(zhì)的交界面，比較適合求解碰撞等動(dòng)態(tài)變形問(wèn)題。另外由于SPH方法使用粒子系統(tǒng)，通過(guò)描述粒子來(lái)代替對(duì)整個(gè)流體的描述，簡(jiǎn)化問(wèn)題的同時(shí)保證了質(zhì)量和動(dòng)量的守恒，減少了很多復(fù)雜的計(jì)算。因此SPH方法在許多領(lǐng)域有著重要的應(yīng)用。

2 海浪建模

SPH方法是一種基于粒子的插值方法，海浪建模時(shí)必須首先將其粒子化。這些海浪粒子在仿真初始的時(shí)候給定，所有的粒子都使用統(tǒng)一的光滑核半徑，而且都具有質(zhì)量、密度、速度、加速度、位置等屬性，這些屬性在模擬仿真過(guò)程中會(huì)不斷地發(fā)生變化。

2.1 基于空間網(wǎng)格的粒子分配與查找

由于海浪模擬過(guò)程中粒子數(shù)量非常大，造成了在后面粒子的屬性計(jì)算過(guò)程中，粒子搜索量較大，效率低等問(wèn)題。為了快速定位某空間位置周圍的粒子，本文使用空間網(wǎng)格的數(shù)據(jù)結(jié)構(gòu)，將粒子的分布空間劃分為網(wǎng)格。對(duì)于網(wǎng)格中每個(gè)粒子的各屬性值以及受力分析，可以通過(guò)該粒子支持域內(nèi)的其它粒子的插值計(jì)算得到，粒子支持域由核函數(shù)的光滑核半徑?jīng)Q定。

采用空間網(wǎng)格劃分的方式，對(duì)于一個(gè)粒子來(lái)說(shuō)，當(dāng)計(jì)算周圍核半徑內(nèi)粒子對(duì)其影響時(shí)，不必查找整個(gè)海浪范圍內(nèi)的粒子，只需查找周圍3*3*3個(gè)網(wǎng)格中所包含的粒子即可，這樣在很大程度上減少了粒子的搜索量與計(jì)算量。由于粒子具有時(shí)空變換的特征，隨著時(shí)間的變換其在空間中的分布范圍也將不斷變化，所以在每幀的開(kāi)始都需要重新劃分空間網(wǎng)格，并根據(jù)粒子的位置將粒子分配到相應(yīng)的網(wǎng)格單元中，圖2是粒子在空間網(wǎng)格中的分布示意圖。

圖2 空間網(wǎng)格中粒子的分布示意圖

在網(wǎng)格結(jié)構(gòu)中，對(duì)于粒子屬性的計(jì)算，需要搜索周圍27個(gè)網(wǎng)格中在光滑核半徑區(qū)域內(nèi)的粒子。本文采用鏈表結(jié)構(gòu)來(lái)存儲(chǔ)粒子，每個(gè)網(wǎng)格內(nèi)的粒子構(gòu)成一個(gè)單向鏈表，通過(guò)指針指向與其相連接的粒子。

2.2 海浪粒子受力分析及狀態(tài)計(jì)算

其中，ρ表示水粒子的密度，u表示海浪粒子的運(yùn)動(dòng)速度，p表示海浪粒子受到的壓力，μ表示水的粘滯系數(shù)。

根據(jù)SPH的核函數(shù)插值以及粒子遵循的拉格朗日流體控制方程1.2與2.1式可以計(jì)算得到海浪粒子所受力以及各個(gè)力產(chǎn)生的加速度變化，計(jì)算這些物理量需要的某位置處海浪粒子的密度計(jì)算公式可以由公式1.2推導(dǎo)得到：

根據(jù)公式1.2，用密度ρ代替式中的A，可以得到密度的計(jì)算公式：

其中密度計(jì)算使用的光滑核函數(shù)為Poly6函數(shù)，由于所有粒子的質(zhì)量相同且都為m，所以在三維情況下，ri處的密度計(jì)算公式最終可表示為：

海浪粒子密度已知后，其受力以及由受力引起的加速度變化可分別表示為：

（1）壓力及壓力加速度

依據(jù)公式1.2，壓力場(chǎng)-▽p可以表示為：

兩個(gè)水粒子在相互作用的時(shí)候，兩者之間的作用力是不對(duì)稱的。通常兩個(gè)粒子位置上的壓力是不對(duì)稱的、不等的。這里提出一個(gè)簡(jiǎn)單的解決辦法是使用兩個(gè)粒子壓強(qiáng)的算術(shù)平均值來(lái)代替單個(gè)粒子的壓力，這樣保證了粒子的速度和穩(wěn)定性。

由于單個(gè)粒子的壓力P，可以用理想氣體狀態(tài)方程P=K（ρ-ρ0）來(lái)表示，其中ρ0表示的是流體的靜態(tài)密度，K是與流體相關(guān)的常數(shù)，只和溫度有關(guān)。上式中壓力計(jì)算用的光滑核函數(shù)為Spiky函數(shù)，在三維情況下，

結(jié)合公式2.5、2.6得到壓力產(chǎn)生的加速度：

（2）粘滯力及粘滯力加速度

這里同樣存在不平衡的問(wèn)題，由于每個(gè)粒子具有不同的速度，粘滯力取決于粒子的相對(duì)速度而不是絕對(duì)速度，所以應(yīng)該使用相對(duì)速度來(lái)計(jì)算粘滯力。

在三維情況下，結(jié)合粘滯力采用的核函數(shù)：

然后結(jié)合公式2.9與2.10便可以得到粘滯力產(chǎn)生的加速度：

（3）總加速度

根據(jù)拉格朗日流體控制方程—公式2.1，粒子i在位置ri處的加速度a（ri）的推導(dǎo)公式為：

最后將上面得到的壓力產(chǎn)生的加速度（公式2.7）與粘滯力產(chǎn)生的加速度（公式2.11）代入到上面的計(jì)算公式中，便可以得到粒子運(yùn)動(dòng)過(guò)程中最后的總加速度。

2.3 海浪粒子與障礙物的碰撞檢測(cè)

海浪由大量粒子組成，而障礙物的邊界可以用多邊形集合來(lái)描述，即障礙物的面可由多個(gè)三角形或者四邊形構(gòu)成。因此，海浪與障礙物的碰撞檢測(cè)可以轉(zhuǎn)化為粒子與多邊形的碰撞檢測(cè)。也就是在粒子的運(yùn)動(dòng)過(guò)程中，檢測(cè)其是否與多邊形發(fā)生碰撞。所以結(jié)合海浪粒子實(shí)際運(yùn)動(dòng)情況，最終可以轉(zhuǎn)化為檢測(cè)粒子在一定時(shí)間段內(nèi)所經(jīng)過(guò)的路徑是否會(huì)與多邊形發(fā)生碰撞。

在各個(gè)時(shí)刻，海浪粒子的運(yùn)動(dòng)路徑可以用粒子的初始位置、速度向量以及時(shí)間步長(zhǎng)來(lái)表示。由于時(shí)間步長(zhǎng)較短，本研究試驗(yàn)采用的時(shí)間步長(zhǎng)為0.004 s，因此粒子的運(yùn)動(dòng)軌跡可近似看作是一條直線段。于是上面海浪粒子與環(huán)境障礙物的碰撞檢測(cè)問(wèn)題便可歸結(jié)為線段與多邊形面片的相交檢測(cè)問(wèn)題。

在整個(gè)碰撞檢測(cè)算法流程中，需要計(jì)算的前提條件有：在某時(shí)刻海浪粒子的速度、加速度以及粒子位置的變化量，即粒子在某時(shí)間間隔內(nèi)的運(yùn)動(dòng)軌跡，此運(yùn)動(dòng)軌跡可看作一條線段來(lái)處理。根據(jù)粒子加速度計(jì)算方法，得到海浪粒子的加速度，進(jìn)一步計(jì)算得到海浪粒子在下一時(shí)刻的速度。在速度、加速度以及時(shí)間間隔已知的情況下，粒子在該時(shí)刻位置的變化量便可很容易得到。另外還要計(jì)算構(gòu)成障礙物三角網(wǎng)的邊向量與法向量。

在上述計(jì)算條件準(zhǔn)備充分之后就可以判斷海浪粒子與障礙物某平面是否會(huì)發(fā)生碰撞，判斷過(guò)程如下：

（1）計(jì)算粒子位置變化量P0P1與平面法向量n的點(diǎn)積，以及粒子某時(shí)刻初始位置到三角面所屬平面上任一點(diǎn)的向量P0P與平面法向量n的點(diǎn)積，如圖3所示。通過(guò)兩個(gè)點(diǎn)積值乘積的正負(fù)判斷海浪粒子是否與三角面所在的平面相交，若乘積為正，則海浪粒子與三角面所在平面有交點(diǎn)，否則沒(méi)有。

圖3 粒子運(yùn)動(dòng)軌跡與平面相交示意圖

（2）若海浪粒子與三角面所在平面相交，然后運(yùn)用Moller提出的射線與三角形求交算法，判斷在設(shè)定的時(shí)間間隔內(nèi)，某時(shí)刻的海浪粒子所在射線是否與三角面相交。

（3）若粒子軌跡所在射線與三角面相交，則計(jì)算返回射線與三角面的交點(diǎn)O，見(jiàn)圖5。交點(diǎn)O的坐標(biāo)值可由上一步中得到的t值、粒子初始位置P0以及定義射線的向量求得。然后比較粒子初始位置到交點(diǎn)之間的距離P0O與粒子在此時(shí)刻短時(shí)間段內(nèi)的運(yùn)動(dòng)距離P0P1之間的關(guān)系，若P0O小于P0P1的值，則海浪粒子與三角面相交，從而完成了海浪粒子與三角面的相交檢測(cè)。

（4）計(jì)算海浪粒子與環(huán)境障礙物碰撞后的加速度、速度以位置等變量。

海浪粒子在與障礙物的碰撞過(guò)程中，將粒子看作是理想的剛體，與障礙物的碰撞后做反彈運(yùn)動(dòng)，如圖4所示。

圖4 粒子與平面碰撞示意圖

3 算法實(shí)現(xiàn)

本文在.NET環(huán)境下利用OpenGL設(shè)計(jì)并實(shí)現(xiàn)一個(gè)粒子系統(tǒng)，一個(gè)微觀尺度下海浪運(yùn)動(dòng)的演示平臺(tái)。在微觀尺度下，基于粒子系統(tǒng)的海浪運(yùn)動(dòng)演示系統(tǒng)由SPH海浪運(yùn)動(dòng)模擬模塊與海浪表面建模渲染模塊組成。其中SPH海浪運(yùn)動(dòng)模擬模塊主要實(shí)現(xiàn)的功能有：初始化海浪粒子，通過(guò)海浪粒子的受力分析與狀態(tài)計(jì)算，以及海浪粒子運(yùn)動(dòng)過(guò)程中與環(huán)境障礙物的碰撞檢測(cè)分析建立海浪運(yùn)動(dòng)模型。

整個(gè)演示系統(tǒng)的設(shè)計(jì)過(guò)程分為下面幾個(gè)步驟：

（1）初始化粒子的屬性，生成基本的海浪粒子類；

（2）分析海浪粒子受力情況，計(jì)算粒子速度、加速度等屬性的變化，模擬海浪的運(yùn)動(dòng)狀態(tài)；

（3）計(jì)算粒子位置，判斷是否會(huì)與周圍環(huán)境障礙物發(fā)生碰撞；若發(fā)生了，重新計(jì)算粒子加速度、速度等變量，以確定粒子碰撞后的位置。圖5為渲染效果圖。

圖5 SPH渲染效果

4 結(jié)語(yǔ)

結(jié)合海浪基本特征，利用物理模型與粒子系統(tǒng)相結(jié)合的方法—SPH方法完成了對(duì)海浪流動(dòng)的模擬，在基于SPH方法模擬海浪流動(dòng)的基礎(chǔ)上，添加環(huán)境障礙物，完成了海浪與障礙物碰撞過(guò)程的模擬。由于海浪運(yùn)動(dòng)場(chǎng)景比較小，不用使用大量的粒子來(lái)模擬，所以模擬的實(shí)時(shí)效果較好。但是當(dāng)場(chǎng)景較大，粒子數(shù)目增多后，模擬的實(shí)時(shí)性會(huì)較差。

[1]丁紹潔.虛擬海洋環(huán)境生成及場(chǎng)景特效研究[D].哈爾濱工程大學(xué)，2008.

[2]高峰.基于SPH的化學(xué)溶液傾倒過(guò)程仿真[D].東北師范大學(xué)，2010.

[3]張靜.基于粒子系統(tǒng)的海浪動(dòng)態(tài)演變模擬仿真研究[D].山東科技大學(xué)，2012.

2016-10-17