基于SPH方法的海浪動態演變模擬仿真

張靜，劉文龍，盧文虎

（1.天津市博盈科技發展有限公司天津市300180；2.中國地震局第一監測中心天津市300180；3.國家海洋信息中心天津市300171）

基于SPH方法的海浪動態演變模擬仿真

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

（1.天津市博盈科技發展有限公司天津市300180；2.中國地震局第一監測中心天津市300180；3.國家海洋信息中心天津市300171）

結合風浪場復雜的動力學、時空特性，針對海浪仿真模擬中存在的真實感不足，難以模擬海浪破碎等問題，根據流體的物理特性，采用光滑粒子流體動力學（SPH）的方法，從粒子初始化分配、粒子受力分析及狀態計算以及粒子運動過程中的碰撞檢測3個方面進行了研究，實現了海浪的動態演變仿真。

海浪模擬；粒子系統；光滑粒子流體動力學；碰撞檢測

海浪仿真是近年來計算機圖形學領域研究的熱點之一，它在虛擬現實應用、計算機游戲以及電影制作中起著重要的作用。另外，海浪的仿真模擬是構建海洋環境信息可視化平臺的一個重要組成部分，平臺的建立可顯著提高我國海洋信息的可視化水平。

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

1 SPH基本原理

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

圖1 光滑核函數影影響域內的粒子示意圖

SPH方法是一種無網格拉格朗日算法，它通過一系列粒子質點的“核函數估值”將流體力學基本方程組轉換成數值計算用的公式。其核心是核函數，它表示在一定的光滑長度范圍內其臨近粒子質點對研究粒子影響程度的權函數。假設流體中某點r，在光滑核半徑h范圍內受數個粒子影響，其位置分別是r0，r1，r2…，rj，如圖1所示，該點位置處某項屬性A的累加值就可以用公式1.2來表示。

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

由于粒子質點之間不存在網格關系，它可以避免因變形引起的網格扭曲從而造成精度破壞等問題，而且也能方便地處理不同介質的交界面，比較適合求解碰撞等動態變形問題。另外由于SPH方法使用粒子系統，通過描述粒子來代替對整個流體的描述，簡化問題的同時保證了質量和動量的守恒，減少了很多復雜的計算。因此SPH方法在許多領域有著重要的應用。

2 海浪建模

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

2.1 基于空間網格的粒子分配與查找

由于海浪模擬過程中粒子數量非常大，造成了在后面粒子的屬性計算過程中，粒子搜索量較大，效率低等問題。為了快速定位某空間位置周圍的粒子，本文使用空間網格的數據結構，將粒子的分布空間劃分為網格。對于網格中每個粒子的各屬性值以及受力分析，可以通過該粒子支持域內的其它粒子的插值計算得到，粒子支持域由核函數的光滑核半徑決定。

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

圖2 空間網格中粒子的分布示意圖

在網格結構中，對于粒子屬性的計算，需要搜索周圍27個網格中在光滑核半徑區域內的粒子。本文采用鏈表結構來存儲粒子，每個網格內的粒子構成一個單向鏈表，通過指針指向與其相連接的粒子。

2.2 海浪粒子受力分析及狀態計算

其中，ρ表示水粒子的密度，u表示海浪粒子的運動速度，p表示海浪粒子受到的壓力，μ表示水的粘滯系數。

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

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

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

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

（1）壓力及壓力加速度

依據公式1.2，壓力場-▽p可以表示為：

兩個水粒子在相互作用的時候，兩者之間的作用力是不對稱的。通常兩個粒子位置上的壓力是不對稱的、不等的。這里提出一個簡單的解決辦法是使用兩個粒子壓強的算術平均值來代替單個粒子的壓力，這樣保證了粒子的速度和穩定性。

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

結合公式2.5、2.6得到壓力產生的加速度：

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

這里同樣存在不平衡的問題，由于每個粒子具有不同的速度，粘滯力取決于粒子的相對速度而不是絕對速度，所以應該使用相對速度來計算粘滯力。

在三維情況下，結合粘滯力采用的核函數：

然后結合公式2.9與2.10便可以得到粘滯力產生的加速度：

（3）總加速度

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

最后將上面得到的壓力產生的加速度（公式2.7）與粘滯力產生的加速度（公式2.11）代入到上面的計算公式中，便可以得到粒子運動過程中最后的總加速度。

2.3 海浪粒子與障礙物的碰撞檢測

海浪由大量粒子組成，而障礙物的邊界可以用多邊形集合來描述，即障礙物的面可由多個三角形或者四邊形構成。因此，海浪與障礙物的碰撞檢測可以轉化為粒子與多邊形的碰撞檢測。也就是在粒子的運動過程中，檢測其是否與多邊形發生碰撞。所以結合海浪粒子實際運動情況，最終可以轉化為檢測粒子在一定時間段內所經過的路徑是否會與多邊形發生碰撞。

在各個時刻，海浪粒子的運動路徑可以用粒子的初始位置、速度向量以及時間步長來表示。由于時間步長較短，本研究試驗采用的時間步長為0.004 s，因此粒子的運動軌跡可近似看作是一條直線段。于是上面海浪粒子與環境障礙物的碰撞檢測問題便可歸結為線段與多邊形面片的相交檢測問題。

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

在上述計算條件準備充分之后就可以判斷海浪粒子與障礙物某平面是否會發生碰撞，判斷過程如下：

（1）計算粒子位置變化量P0P1與平面法向量n的點積，以及粒子某時刻初始位置到三角面所屬平面上任一點的向量P0P與平面法向量n的點積，如圖3所示。通過兩個點積值乘積的正負判斷海浪粒子是否與三角面所在的平面相交，若乘積為正，則海浪粒子與三角面所在平面有交點，否則沒有。

圖3 粒子運動軌跡與平面相交示意圖

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

（3）若粒子軌跡所在射線與三角面相交，則計算返回射線與三角面的交點O，見圖5。交點O的坐標值可由上一步中得到的t值、粒子初始位置P0以及定義射線的向量求得。然后比較粒子初始位置到交點之間的距離P0O與粒子在此時刻短時間段內的運動距離P0P1之間的關系，若P0O小于P0P1的值，則海浪粒子與三角面相交，從而完成了海浪粒子與三角面的相交檢測。

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

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

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

3 算法實現

本文在.NET環境下利用OpenGL設計并實現一個粒子系統，一個微觀尺度下海浪運動的演示平臺。在微觀尺度下，基于粒子系統的海浪運動演示系統由SPH海浪運動模擬模塊與海浪表面建模渲染模塊組成。其中SPH海浪運動模擬模塊主要實現的功能有：初始化海浪粒子，通過海浪粒子的受力分析與狀態計算，以及海浪粒子運動過程中與環境障礙物的碰撞檢測分析建立海浪運動模型。

整個演示系統的設計過程分為下面幾個步驟：

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

（2）分析海浪粒子受力情況，計算粒子速度、加速度等屬性的變化，模擬海浪的運動狀態；

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

圖5 SPH渲染效果

4 結語

結合海浪基本特征，利用物理模型與粒子系統相結合的方法—SPH方法完成了對海浪流動的模擬，在基于SPH方法模擬海浪流動的基礎上，添加環境障礙物，完成了海浪與障礙物碰撞過程的模擬。由于海浪運動場景比較小，不用使用大量的粒子來模擬，所以模擬的實時效果較好。但是當場景較大，粒子數目增多后，模擬的實時性會較差。

[1]丁紹潔.虛擬海洋環境生成及場景特效研究[D].哈爾濱工程大學，2008.

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

[3]張靜.基于粒子系統的海浪動態演變模擬仿真研究[D].山東科技大學，2012.

2016-10-17