帶有漂浮對象的水波模擬

韓林呈， 陳喜春(機械化步兵學院 作戰訓練實驗中心， 石家莊 050083)

帶有漂浮對象的水波模擬

韓林呈， 陳喜春
(機械化步兵學院 作戰訓練實驗中心， 石家莊 050083)

針對實時水波模擬中二維波動方程解法復雜度高、編程求解困難的缺點，基于網格對水面建模，在時間和空間軸上對水波運動離散化處理，并使用數值方法對微分方程近似求解。且為滿足場景中水面與物體交互需要，對漂浮對象進行物理建模，通過計算漂浮對象在水中浮力和水波對其推力來模擬水面上的漂浮狀態。運行結果表明，該方法可對水波及漂浮物進行實時仿真，模擬速度快，視覺效果好，漂浮效果真實自然。

水波; 物理模型; 波動方程; 離散化; 浮力

0 引言

動態的水可以為室外環境的模擬增加較大的美感，一直是計算機圖形學領域研究的重點。作為一種非固定形態的流體，實時地對其進行精確描述較為困難，國內外許多學者均對此有過相關研究。

通常來說，水波模擬主要分為兩類：一類方法是基于數學函數模擬繪制水波形狀。如文獻[1]基于傅立葉合成的方法來模擬海洋波浪；文獻[2]采用Perlin噪聲源的預計算模擬水面的持續抖動來生成水波效果；文獻[3]對文獻[1]進行改進，基于逆傅立葉變換生成若干線性函數，將這些函數相互疊加后生成水波；文獻[4]將不同相位，不同頻率的正弦波的冪函數相互疊加模擬水面波浪；文獻[5]受文獻[4]啟發，但改用余弦的自然指數函數疊加的方法模擬波浪。另一類方法是基于物理模型構建水波。通過求解流體方程，獲得流體質點在各個時間的坐標。比較典型的如文獻[6]通過解Navies Stokes方程來模擬水面；文獻[7]在文獻[6]的基礎上，采取鄰域傳播的思想，離散求解Navies Stokes方程，提高了構造水波形狀的速度；文獻[8]引入物理上用于海水建模的Gerstner波模型對波浪進行構造；文獻[9]在充分考慮液體粘性和擴散特性的基礎上，基于Lattice Boltzmann流體模型實現了兩種液體的混合模擬；文獻[10]基于二維波動方程描述水波，使用3D引擎OGRE模擬了水面動態蕩漾的效果。

在上述兩類方法中，基于數學函數模擬繪制的方法實現簡單，在相位、振幅等方面易于控制，但生成的波形規律性較強，不夠自然；基于物理方程的方法產生的效果較為真實，但求解復雜、計算量大，實時性較差。為提高求解流體物理方程的效率，本文對水面質點坐標及運動時間分別進行了離散化處理，通過數值方法來近似求解二維波動方程，簡化了計算過程，提高了求解速度。此外，大多數學者對水波本身的構成和模擬進行了研究，無論在真實性和實時性上均取得了很好的效果，但對水波與水面漂浮物的互動仿真較少有涉足。本文通過計算物體浮力和水波推力，對物體在水面漂浮滑行的運動過程進行了仿真。

1 二維波動方程

波動方程是一種表示波動現象的偏微分方程，可對自然界的聲波、水波等現象進行描述。當水面上無限小的部分移動時，水粒子的直接鄰近點會施加線性“彈力”(表面張力)來最小化粒子之間的距離。由于水平方向的力是相等的，粒子僅在z方向運動。關于時間和水粒子空間的位置可以由二維波動方程來描述，如式(1)。

(1)

這里c是波越過水面傳播的速度。當邊界條件為齊次且水面的初始z方向速度為0時，對于一個L×L大小的正方形水域的通解為式(2)、(3)。

(2)

(3)

系數Amn由計算下面的積分得到式(4)。

(4)

這里f(x,y)是水面的初始形態[11]。如果將水面離散化為以z為高度域且平均分隔的柵格(圖1)，則式(4)可以用FFT變換解出[12]。但在實時應用中計算大量的三角函數會大幅降低程序的運行效率。為克服該缺點，本文采用數值方法求近似式(1)。

圖1 一個用來近似表示水面的L×L網格，每條邊有N個離散點

2 近似求解方法

將流體表面離散化表示為一個頂點排列成表面積為L×L，且每條邊含有N個點的規則柵格網絡，如圖2所示。

圖2 計算該點x方向切線斜率

(5)

(6)

(7)

二階導數可通過計算一階導數的近似值的方法得出，即計算一階導數差的平均值來近似求出二階導數。該點處一階導數的平均差為式(8)。

(8)

將(5)式代入(8)式得式(9)。

(9)

(10)

可以看出，在某個頂點處二階導數的計算過程中要用到與該頂點距離兩個頂點間距處的頂點位移。將坐標系以(i,j)頂點為中心縮小1/2，同時Δx也縮小1/2，就可以獲得對于x的二階導數的等價近似方程，即式(11)。

(11)

同理可得z相對于y和t的二階導數近似值為式(12)、(13)。

(12)

(13)

將式(11)、(12)、(13)代入式(1)，可以得出式(14)。

(14)

(15)

上式說明了點zi,j僅受其相鄰點的影響如圖3所示。

圖3 點zi,j的運動僅受相鄰點影響

3 浮力對象

當物體的密度小于水密度時，就可浮于水面之上，物體上的浮力等于其排除的水的重量。該力的方向與氣壓梯度(pressure gradient)相同，在近似模擬中，可以取水面的法線方向。將船體按照離散的點建模，則其浸入水中的浮力方向，見圖4所示。

圖4 漂浮物的外形以一組pi近似

(18)

上式中zwater是在pk的雙線性差值水面高度。在這個位置的浮力為式(19)。

(19)

且轉動力矩為式(20)。

Nk=rkFk

(20)

上式中rk是從質心到pk的矢量。一個對任意參數化的右手螺旋三維表面的法向量可通過式(21)進行計算[13]為式(21)。

(21)

如果將x和y作為參數，則水面可以用下面的矢量描述為式(22)。

S(x,y)water=[x,y,z(x,y,t)]

(22)

用本文第2節所述方法近似(22)式的導數得式(23)、(24)。

(23)

(24)

則位于第i,j柵格處的法向量為式(25)。

(25)

將該向量放大2h倍，得到等價有效法向量為式(26)。

ni,j=[zi-1,j-zi+1,j,zi,j-1-zi,j+1,2h]

(26)

漂浮物在水面上的推力為漂浮物各點的分力之和為式(27)。

Fdrag=∑-bvk,rel=∑b[vwater-(vcm+ωcm×rk)]

(27)

其中，vk,rel是在rk處漂浮物與水的相對速度，b是相對速度與合力之比。

4 算法實現

4.1 計算z值緩沖區

為實現式(15)所示算法，似乎有必要為在時刻t-1,t0,t1開辟3塊內存區域。但是如果zn和zn+1適當代換，則只需要兩塊內存區域即可。在首次計算結束后，指向zn+1和zn的指針交換。在下一次迭代中，原本指向zn的指針指向了zn+1，而指向zn+1的指針則指向了zn+2。這樣只通過兩塊內存區域就可以實現本文算法。下面是實現該節省內存方法的主要代碼。

floatz[N][N];//z^n值

floatz1[N][N];//z^(n-1)和z^(n+1)值

floatu[N][N];//粘性阻力

floatd[N][N];//阻尼系數

…

const floatA= (c*dt/h)*(c*dt/h);

const floatB= 2 - 4*A;

longi,j;

for(i=1 ;ilt;N-1 ;i++ )

{

for(j=1 ;jlt;N-1 ;j++ )

{

u[i][j] =z1[i][j];

//以z^(n-1)和z^n更新z^(n+1)

z1[i][j] =A*(z[i-1][j] +z[i+1][j] +z[i][j-1] +z[i][j+1] )+B*z[i][j] -z1[i][j];

u[i][j] = (z1[i][j]-u[i][j])/2t;//計算粘性阻力

z1[i][j] -=u[i][j]; //應用粘性阻力

//z1[i][j] *=d[i][j]; //或者應用阻尼系數

}

}

//交換z^n和z^(n+1)

swap(z,z1 );

4.2 增加波源

為了形成水波，需在水中加入波源。程序運行初始，在某點增加該位置的z值，則波浪以該點為中心形成。以下是增加波源的主要代碼。

const longN= 128 + 1;//N為水面邊界的離散點數目

…

z[N/3][N/3] =z1[N/3][N/3] = 10; //創建波浪

5 測試

在CPU為i7-5500 2.40GHz，顯卡為AMD Radeon R5，內存4GB，32位Win7系統下，基于C++和OPENGL3.2對本文算法進行編程測試，效果如圖5-圖12所示。

圖5 網格水平面

圖6 創建水波，波浪從該點發出

圖7 水波開始擴散

圖8 物體漂浮在水面

圖9 水波推動漂浮物體滑行

圖10 加入光照及進行渲染后效果1

圖11 加入光照及進行渲染后效果2

圖12 加入光照及進行渲染后效果3

在程序開始運行時，首先初始化所有的z值為0，即處于一個平靜的水平面(圖5)。然后選取任一點，設置一個較大的z值，從該點開始進行擴散(圖6、7)。該點的z值決定了一開始水波波動的高度。圖8、圖9為物體漂浮及在受水波推動滑行的效果。圖10-12為加入光照并進行渲染后的效果圖。文中場景均為128×128個網格，刷新率達到59幀/秒，運行流暢。

7 總結

本文基于二維波動方程，采用數值方法對該物理方程快速求解，同時對水面漂浮物所受浮力和推力進行計算，最后基于OPENGL對水面的波動效果和與漂浮物的交互效果進行了模擬。該方法程序設計簡單，計算量小，運算速度快，對漂浮物體的仿真效果真實感較強。但通過本文算法實現的水波還是較為規則，與真實自然環境中多變的情況仍有差距。在下一步工作中，還將參考更多流體物理方程(如N-S方程)的引入，并加入隨機因素及其他水面特效(如浪花)等。

[1] MASTIN G A，WATTERBERG P A，MAREDA J F．Fourier synthesis of ocean scenes[J]．IEEE Computer Graphics amp; Applications, 1987, 7(3): 16-23.

[2] Hugo Elias. Perlin noise [EB/OL]. http://free space.virgin.net/hugo.elias/models/m_perlin.htm.

[3] 聶衛東, 康鳳舉, 褚彥軍, 等. 基于線性海浪理論的海浪數值模擬[J]. 系統仿真學報, 2005, 5(17): 1037-1044.

[4] 方建文, 于金輝, 馬文龍. 圖形硬件加速的實時水面繪制[J]. 計算機工程與應用, 2006, 42(15): 86-88.

[5] 劉曉平, 謝文軍. 實時水面模擬方法研究[J]. 工程圖學學, 2010, 1: 79-83.

[6] FOSTER N, M ETEXAS D. Realistic animation of liquids [J]． Graphical Models and Image Processing, 1996, 58(5): 471-483.

[7] 吳獻, 董蘭芳, 盧德唐. 一種基于鄰域傳播的水波模擬方法[J]. 中國科技大學學報, 2010, 40(3): 278-282.

[8] 王海玲, 印桂生, 張菁, 等. 基于改進曲面熵的動態水面模擬方法[J]. 計算機工程, 2011, 37(6): 24-26.

[9] 朱紅斌，劉學慧，柳有權，等．基于Lattice Boltzmann模型的液-液混合流模擬[J]． 計算機學報, 2006, 29(12): 2071-2079.

[10] 孫曉鵬, 李翠芳. 三維游戲中基于OGRE的動態水面模擬算法[J]. 計算機工程與設計, 2011, 32(12):4122-4124.

[11] Trim D W. Applied Partial Differential Equations[M]. PWS-Kent, 1990.

[12] William H, Teukolsky Saul A, Vetterling,etal. Numerical Recipes in C (second edition)[M]. Cambridge. The Press Syndicate of the University of Cambridge, 1992.

[13] Davis, Harry F., Snider, Arthur David, Introduction to Vector Analysis (sixth edition)[M]. William C. Brown Publishers, 1991.

WaveSimulationwithaFloatingObject

Han Lincheng, Chen Xichun
(Combat Training Experiment Center, Mechanized Infantry Academy, Shijiazhuang 050083)

For the high complexity and programming difficulty in solving 2-dimensional wave equation in real-time wave simulation, a model is bullt based on water face grid, the wave motion is olescreted in time and space axis, and numerical methods are used to solv differential equation. In order to meet the need of the interaction between the surface and the object in the scene, the floating object is modeled by the physical modeling, and the floating state of the object is simulated by calculating the buoyancy and thrust. The results show that the method can simulate the wave and floating object in real time, the simulation speed is fast, the visual result is good, and the floating effect is real and natural.

Wave; Physical modeling; Wave equation; Discrete; Buoyancy

韓林呈(1984-)，男，石家莊市人，碩士，講師，研究方向為計算機圖形學。

1007-757X(2017)11-0069-05

TP391

A

2016.10.29)