韓林呈， 陳喜春(機(jī)械化步兵學(xué)院 作戰(zhàn)訓(xùn)練實(shí)驗(yàn)中心， 石家莊 050083)

帶有漂浮對(duì)象的水波模擬

韓林呈， 陳喜春
(機(jī)械化步兵學(xué)院 作戰(zhàn)訓(xùn)練實(shí)驗(yàn)中心， 石家莊 050083)

針對(duì)實(shí)時(shí)水波模擬中二維波動(dòng)方程解法復(fù)雜度高、編程求解困難的缺點(diǎn)，基于網(wǎng)格對(duì)水面建模，在時(shí)間和空間軸上對(duì)水波運(yùn)動(dòng)離散化處理，并使用數(shù)值方法對(duì)微分方程近似求解。且為滿足場(chǎng)景中水面與物體交互需要，對(duì)漂浮對(duì)象進(jìn)行物理建模，通過(guò)計(jì)算漂浮對(duì)象在水中浮力和水波對(duì)其推力來(lái)模擬水面上的漂浮狀態(tài)。運(yùn)行結(jié)果表明，該方法可對(duì)水波及漂浮物進(jìn)行實(shí)時(shí)仿真，模擬速度快，視覺(jué)效果好，漂浮效果真實(shí)自然。

水波; 物理模型; 波動(dòng)方程; 離散化; 浮力

0 引言

動(dòng)態(tài)的水可以為室外環(huán)境的模擬增加較大的美感，一直是計(jì)算機(jī)圖形學(xué)領(lǐng)域研究的重點(diǎn)。作為一種非固定形態(tài)的流體，實(shí)時(shí)地對(duì)其進(jìn)行精確描述較為困難，國(guó)內(nèi)外許多學(xué)者均對(duì)此有過(guò)相關(guān)研究。

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

在上述兩類方法中，基于數(shù)學(xué)函數(shù)模擬繪制的方法實(shí)現(xiàn)簡(jiǎn)單，在相位、振幅等方面易于控制，但生成的波形規(guī)律性較強(qiáng)，不夠自然；基于物理方程的方法產(chǎn)生的效果較為真實(shí)，但求解復(fù)雜、計(jì)算量大，實(shí)時(shí)性較差。為提高求解流體物理方程的效率，本文對(duì)水面質(zhì)點(diǎn)坐標(biāo)及運(yùn)動(dòng)時(shí)間分別進(jìn)行了離散化處理，通過(guò)數(shù)值方法來(lái)近似求解二維波動(dòng)方程，簡(jiǎn)化了計(jì)算過(guò)程，提高了求解速度。此外，大多數(shù)學(xué)者對(duì)水波本身的構(gòu)成和模擬進(jìn)行了研究，無(wú)論在真實(shí)性和實(shí)時(shí)性上均取得了很好的效果，但對(duì)水波與水面漂浮物的互動(dòng)仿真較少有涉足。本文通過(guò)計(jì)算物體浮力和水波推力，對(duì)物體在水面漂浮滑行的運(yùn)動(dòng)過(guò)程進(jìn)行了仿真。

1 二維波動(dòng)方程

波動(dòng)方程是一種表示波動(dòng)現(xiàn)象的偏微分方程，可對(duì)自然界的聲波、水波等現(xiàn)象進(jìn)行描述。當(dāng)水面上無(wú)限小的部分移動(dòng)時(shí)，水粒子的直接鄰近點(diǎn)會(huì)施加線性“彈力”(表面張力)來(lái)最小化粒子之間的距離。由于水平方向的力是相等的，粒子僅在z方向運(yùn)動(dòng)。關(guān)于時(shí)間和水粒子空間的位置可以由二維波動(dòng)方程來(lái)描述，如式(1)。

(1)

這里c是波越過(guò)水面?zhèn)鞑サ乃俣?。?dāng)邊界條件為齊次且水面的初始z方向速度為0時(shí)，對(duì)于一個(gè)L×L大小的正方形水域的通解為式(2)、(3)。

(2)

(3)

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

(4)

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

圖1 一個(gè)用來(lái)近似表示水面的L×L網(wǎng)格，每條邊有N個(gè)離散點(diǎn)

2 近似求解方法

將流體表面離散化表示為一個(gè)頂點(diǎn)排列成表面積為L(zhǎng)×L，且每條邊含有N個(gè)點(diǎn)的規(guī)則柵格網(wǎng)絡(luò)，如圖2所示。

圖2 計(jì)算該點(diǎn)x方向切線斜率

(5)

(6)

(7)

二階導(dǎo)數(shù)可通過(guò)計(jì)算一階導(dǎo)數(shù)的近似值的方法得出，即計(jì)算一階導(dǎo)數(shù)差的平均值來(lái)近似求出二階導(dǎo)數(shù)。該點(diǎn)處一階導(dǎo)數(shù)的平均差為式(8)。

(8)

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

(9)

(10)

可以看出，在某個(gè)頂點(diǎn)處二階導(dǎo)數(shù)的計(jì)算過(guò)程中要用到與該頂點(diǎn)距離兩個(gè)頂點(diǎn)間距處的頂點(diǎn)位移。將坐標(biāo)系以(i,j)頂點(diǎn)為中心縮小1/2，同時(shí)Δx也縮小1/2，就可以獲得對(duì)于x的二階導(dǎo)數(shù)的等價(jià)近似方程，即式(11)。

(11)

同理可得z相對(duì)于y和t的二階導(dǎo)數(shù)近似值為式(12)、(13)。

(12)

(13)

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

(14)

(15)

上式說(shuō)明了點(diǎn)zi,j僅受其相鄰點(diǎn)的影響如圖3所示。

圖3 點(diǎn)zi,j的運(yùn)動(dòng)僅受相鄰點(diǎn)影響

3 浮力對(duì)象

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

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

(18)

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

(19)

且轉(zhuǎn)動(dòng)力矩為式(20)。

Nk=rkFk

(20)

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

(21)

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

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

(22)

用本文第2節(jié)所述方法近似(22)式的導(dǎo)數(shù)得式(23)、(24)。

(23)

(24)

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

(25)

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

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

(26)

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

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

(27)

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

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

4.1 計(jì)算z值緩沖區(qū)

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

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

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

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

floatd[N][N];//阻尼系數(shù)

…

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;//計(jì)算粘性阻力

z1[i][j] -=u[i][j]; //應(yīng)用粘性阻力

//z1[i][j] *=d[i][j]; //或者應(yīng)用阻尼系數(shù)

}

}

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

swap(z,z1 );

4.2 增加波源

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

const longN= 128 + 1;//N為水面邊界的離散點(diǎn)數(shù)目

…

z[N/3][N/3] =z1[N/3][N/3] = 10; //創(chuàng)建波浪

5 測(cè)試

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

圖5 網(wǎng)格水平面

圖6 創(chuàng)建水波，波浪從該點(diǎn)發(fā)出

圖7 水波開(kāi)始擴(kuò)散

圖8 物體漂浮在水面

圖9 水波推動(dòng)漂浮物體滑行

圖10 加入光照及進(jìn)行渲染后效果1

圖11 加入光照及進(jìn)行渲染后效果2

圖12 加入光照及進(jìn)行渲染后效果3

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

7 總結(jié)

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

[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] 聶衛(wèi)東, 康鳳舉, 褚彥軍, 等. 基于線性海浪理論的海浪數(shù)值模擬[J]. 系統(tǒng)仿真學(xué)報(bào), 2005, 5(17): 1037-1044.

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

[5] 劉曉平, 謝文軍. 實(shí)時(shí)水面模擬方法研究[J]. 工程圖學(xué)學(xué), 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] 吳獻(xiàn), 董蘭芳, 盧德唐. 一種基于鄰域傳播的水波模擬方法[J]. 中國(guó)科技大學(xué)學(xué)報(bào), 2010, 40(3): 278-282.

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

[9] 朱紅斌，劉學(xué)慧，柳有權(quán)，等．基于Lattice Boltzmann模型的液-液混合流模擬[J]． 計(jì)算機(jī)學(xué)報(bào), 2006, 29(12): 2071-2079.

[10] 孫曉鵬, 李翠芳. 三維游戲中基于OGRE的動(dòng)態(tài)水面模擬算法[J]. 計(jì)算機(jī)工程與設(shè)計(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-)，男，石家莊市人，碩士，講師，研究方向?yàn)橛?jì)算機(jī)圖形學(xué)。

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

TP391

A

2016.10.29)