張志毅，呂之華

(西北農(nóng)林科技大學(xué)信息工程學(xué)院，陜西 楊凌 712100)

基于物理的波浪實(shí)時模擬方法

張志毅，呂之華

(西北農(nóng)林科技大學(xué)信息工程學(xué)院，陜西 楊凌 712100)

為解決廣域范圍內(nèi)波浪狀態(tài)的實(shí)時計(jì)算和可視化問題，結(jié)合無粘滯流體力學(xué)的物理模型，提出了一種以三角形為控制單元的有限體積簡化算法。該算法的優(yōu)點(diǎn)在于：以不規(guī)則邊界區(qū)域的非規(guī)則三角網(wǎng)格為基礎(chǔ)，通過簡化通量向量分裂方法獲得三角形控制單元的邊界數(shù)值通量，能快速逼近二維淺水方程的解進(jìn)而模擬非規(guī)則邊界淺水的實(shí)時流動。實(shí)驗(yàn)結(jié)果顯示，所提方法能在符合現(xiàn)實(shí)世界物理規(guī)律的前提下較好地實(shí)現(xiàn)大規(guī)模波浪的實(shí)時可視化模擬。

淺水方程;限定Delaunay三角剖分;有限體積法;通量向量分裂法

1 引言

水以其在自然與社會中的不可或缺及其本身所擁有的豐富特性和行為，一直以來都是人們的重點(diǎn)研究對象。近年來，因?qū)崟r水流仿真在現(xiàn)實(shí)應(yīng)用中的迫切需求，已吸引了眾多學(xué)者對其展開相關(guān)研究。Navier-Stokes方程是描述流體運(yùn)動的基本方程，但其數(shù)學(xué)表達(dá)和求解復(fù)雜且困難，在通常情況下甚至不可能求得解析解。當(dāng)求其數(shù)值解時，常因所選擇的方法導(dǎo)致耗時極長。對于海嘯模擬演化、洪水破壞預(yù)測、計(jì)算機(jī)動畫等實(shí)時性較強(qiáng)的應(yīng)用程序來說，模擬水流需要極快的速度。嚴(yán)格來說，水的流動是三維的，但為滿足快速計(jì)算的需要可將其降為二維問題，即將水體看成在某一個平面上的高度域，并使用“淺水”[1]條件簡化Navier-Stokes方程獲得雙曲型偏微分方程組，用于描述水流的表面信息。對滿足“淺水”條件的流體，可通過計(jì)算流體力學(xué)方法求得淺水方程的數(shù)值解。主要方法有三種：有限差分法、有限元法和有限體積法。有限差分法FDM(Finite Difference Method)[2]簡單易于理解，速度快，但精度不高。有限元法FEM(Finite Element Method)常用于固體力學(xué)，用于流體力學(xué)方面還有待改善。有限體積法FVM(Finite Volume Method)[3～5]的基本思路是將計(jì)算區(qū)域劃分為一系列不重復(fù)的控制體，并使每個網(wǎng)格點(diǎn)周圍有一個控制體，將待解的微分方程對每一個控制體積分，得出一組離散方程。使用FVM是對推導(dǎo)原始微分方程所用控制體途徑的回歸，其數(shù)值逼近的物理意義更直接明晰。

綜上所述，F(xiàn)VM能像FEM一樣適用于任意不規(guī)則網(wǎng)格，著眼于控制體上的逼近，具有守恒性，其處理效率與FDM相近，遠(yuǎn)高于FEM。也就是說，F(xiàn)VM集FEM的幾何靈活性和FDM的高效率于一體，因此常用于解決流體方面的問題。

2 有限體積法解淺水方程

淺水方程[6,7]的守恒形式可表示為：

(1)

其中，

式中，h是水深，u是x方向的速度，v是y方向的速度，g是重力加速度。B(U)是源項(xiàng)。淺水方程是一個非線性偏微分方程組，得不到解析解，需數(shù)值方法求解。本文使用有限體積法解淺水方程。

2.1 控制體的選擇與生成

使用有限體積法時，需先選擇適合的控制體作為計(jì)算區(qū)域中的最小單元[8]。本文選擇三角形作為基本控制體。對于非規(guī)則邊界的區(qū)域，通?？捎孟薅―elaunay 三角剖分CDT(Constrained Delaunay Triangulation)方法[9～11]來生成高質(zhì)量的非規(guī)則網(wǎng)格。

為實(shí)現(xiàn)數(shù)值離散方法[12～14]和通量向量分裂方法的快速精確求解，通常有兩種方法來選擇控制體： 格子中心式CC(Cell-Centered)和格子頂點(diǎn)式VC(Vertex-Centered)。以三角形網(wǎng)格為例來表達(dá)這兩種格式的區(qū)別是：CC 的每一個控制體是一個三角形，三角形中各物理量被認(rèn)為是常數(shù)；VC 的每一個控制體是以某一頂點(diǎn)為中心，以其周圍的三角形的形心/外心連線為邊界構(gòu)成的多邊形。本文選用三角形的CC作為控制體。

2.2 數(shù)值離散

當(dāng)選用格子中心式的三角形為基本控制體時，如圖1所示，n是邊的法向量，θ是n與x軸的夾角。對式(1)積分，可得式(2)。

Figure 1 Schematic diagram of control unit圖1 控制單元示意圖

(2)

每個控制體單元中的U以常數(shù)近似，即當(dāng)U的值均相同時，左端第一項(xiàng)可表示為：

(3)

其中，A是控制體的面積。

對于式(2)中左端第二項(xiàng)，由格林公式可得：

(4)

式(2)右端可表達(dá)為：

(5)

(6)

對于時間項(xiàng)，采用前向差分格式：

(7)

對于式(6)的左端第二項(xiàng)，若令：

E=Fcosθ+Gsinθ

(8)

則(6)式的左端第二項(xiàng)可表示為：

(9)

將式(7)和式(9)代入式(6)，可得到式(10)和式(11)：

(10)

(11)

因歐拉方程具有旋轉(zhuǎn)不變性[15]，而淺水方程的數(shù)學(xué)形式和歐拉方程是相同的，于是有：

F(U)cosθ+G(U)sinθ=T(θ)-1F(T(θ)U)

(12)

其中，

(13)

由式(12)和式(13)可得：

(14)

其中，

(15)

則式(11)可化為式(16)：

(16)

2.3 通量向量分裂法解黎曼問題

當(dāng)水流在空間或時間域產(chǎn)生不連續(xù)變化時，凡在單元交界面上發(fā)生的下述現(xiàn)象均屬于黎曼問題：由初始間斷的左右狀態(tài)，確定t>0時的波態(tài)、波強(qiáng)度及波間的流動特性的問題。黎曼問題是解決不連續(xù)物理問題的數(shù)學(xué)描述方程，是設(shè)計(jì)高解析度數(shù)值方法的核心。其微分方程形式為雙曲守恒律方程組：

(17)

其初值條件為：

(18)

(19)

其中J是雅克比矩陣。J可以表示為：

J=PΛP-1

(20)

其中Λ=diag(λi)，P是J的右特征向量構(gòu)成的矩陣。根據(jù)λi的符號不同，J可表示為：

J=PΛP-1=P(Λ++Λ-)P-1=J++J-

(21)

Figure 2 Schematic diagram of Riemann problem圖2 黎曼問題示意圖

由式(19)～式(21)可得：

(22)

(23)

(24)

其中，

(25)

(26)

通過式(16)、式(26)和給出的邊界條件，即可求解淺水方程。

2.4 邊界條件

水域的邊界分為兩種情況：開邊界和閉邊界。如圖2所示，對于邊界情況，L是存在的，R是不存在的?，F(xiàn)假設(shè)一個虛擬的R，然后給它賦值。

對于開邊界，法向的值可表示為：

(27)

對于閉邊界，法向的值可表示為：

(28)

基于上述理論，水面的高低變化可精確計(jì)算。尚存在的問題如圖3所示：C1至C5的高度能計(jì)算，但V點(diǎn)的高度不能。這是因?yàn)檫x取的控制體為三角形，而V并不是屬于某一個控制體內(nèi)的物理量。然而為了模擬水流的邊界情況，V點(diǎn)是必需的。考慮到剖分結(jié)果中的三角形都趨近于正三角形，當(dāng)其尺寸很小時，V點(diǎn)高度可近似用式(29)計(jì)算。

(29)

其中，n是V點(diǎn)周圍三角形的數(shù)量。當(dāng)三角形尺寸很小時，用平均法計(jì)算V點(diǎn)高度，對于建立實(shí)時水流動畫來說，可以取得速度和精度的平衡。

Figure 3 Obtain the water height in point V圖3 獲取V點(diǎn)的水高

3 實(shí)驗(yàn)結(jié)果

實(shí)驗(yàn)結(jié)果如圖4和圖5所示。本實(shí)驗(yàn)程序在Visual Studio 2010下開發(fā)，圖形庫使用OpenGL。程序運(yùn)行于Intel Pentium(R) Dual-Core CPU T4300，2.10 GHz，內(nèi)存為2 GB，顯卡為NVIDIA GeForce G210M。

圖4顯示對一組真實(shí)地形數(shù)據(jù)進(jìn)行三角剖分的結(jié)果。該數(shù)據(jù)來源于陜西省某天然湖面的邊界線，湖面長約4.8 km，寬約2.8 km，非規(guī)則邊界線約長17.2 km，原始數(shù)據(jù)為精度50 m的DEM數(shù)據(jù)。圖4a中每個三角形外接圓半徑都小于14 m。圖4b中每個三角形外接圓半徑都小于7 m。結(jié)果顯示，本三角剖分算法效果良好，且易于控制三角形尺寸和質(zhì)量。在實(shí)際使用中，三角形尺寸要根據(jù)情況選定，考慮到精度，建議外接圓半徑小于3 m。

圖5顯示如圖4所示地形的模擬結(jié)果。圖5中，三角形外接圓半徑取為1 m，三角剖分產(chǎn)生的三角形為56 046個，運(yùn)行時狀態(tài)間的時間間隔為Δt=30 ms，即每秒可實(shí)現(xiàn)至少33幀精細(xì)圖像，達(dá)到了實(shí)時模擬計(jì)算。圖5a～圖5d分別為激波后2 280 ms、5 460 ms、7 950 ms和12 810 ms時刻湖面狀態(tài)的模擬結(jié)果。因地形中有“灣”存在，數(shù)值計(jì)算結(jié)果和可視化結(jié)果均顯示：當(dāng)水流進(jìn)入“灣”后，因流量固定，受地形擠壓，其水面高度會高于寬闊“湖面”上的波高，且水流速度也較快。這真實(shí)地反映了波浪運(yùn)動遵循伯努利流體方程的實(shí)際情況。

Figure 4 Triangulation results of a real terrain boundary region圖4 真實(shí)地形邊界區(qū)域的三角剖分結(jié)果

Figure 5 Wave simulation results corresponding to Figure 4 terrain圖5 對應(yīng)于圖4中地形的波浪模擬結(jié)果

表1顯示了本方法與其它三種方法耗時對比結(jié)果。比較時所用的數(shù)據(jù)均為如圖4所示的相同數(shù)據(jù)，所生成的圖像分辨率均為1 024×768像素。其中差分逼近法是采用文獻(xiàn)[17]的方法實(shí)現(xiàn)；余弦波疊加法采用文獻(xiàn)[18]的方法實(shí)現(xiàn)；Perlin 噪聲法是采用文獻(xiàn)[19]的方法實(shí)現(xiàn)。本文所用方法和差分逼近法及余弦波疊加法均是以物理模型為基礎(chǔ)的仿真計(jì)算，而Perlin噪聲法僅是為滿足視覺效果而做的虛擬現(xiàn)實(shí)處理。從實(shí)驗(yàn)結(jié)果可以看出，本文所提的方法能真正滿足較大水域內(nèi)波浪的實(shí)時模擬。

Table 1 Time-consuming comparison among the present method and other methods

因本文所提方法簡化了流體內(nèi)部的摩擦損耗，模擬結(jié)果與實(shí)際水流波浪相比較，運(yùn)動速度略快，且在波的干涉現(xiàn)象表現(xiàn)上尚有不足。

4 結(jié)束語

本文提供了一種模擬非規(guī)則邊界水流運(yùn)動的方法。簡單來說分為兩部分：三角剖分和淺水方程的數(shù)值計(jì)算。為獲得良好的計(jì)算區(qū)域，本文使用了高效的限定Delaunay三角剖分算法，并進(jìn)行尺寸控制和質(zhì)量控制，使得三角剖分足夠密且更趨近于正三角形。接著使用有限體積法解偏微分方程，為了解決控制體邊界上的數(shù)值通量問題，使用了通量向量分裂法(FVS)。實(shí)驗(yàn)結(jié)果表明，本方法簡單穩(wěn)定且速度很快，能滿足實(shí)時性要求。

本文主要側(cè)重于水流模擬的理論研究，因此在水流的渲染上還有待提高。此后的工作重點(diǎn)是水流的渲染，主要包括粒子生成和光照影響。

[1] Tan Wei-yan. Computational shallow water hydrodynamics—the application of the finite volume method[M]. Beijing:Tsinghua University Press, 1998.(in Chinese)

[2] Causon D M, Mingham C G. Introductory finite difference methods for PDEs[EB/OL].[2010-05-16]. http://bookboon.methods-for-pdes-ebook.

[3] Li Ren-xian.The basis of finite volume method[M]. Beijing:National Defence Industry Press, 2008.(in Chinese)

[4] Versyeeg H K,Malalasekera W.An introduction to computational fluid dynamics[M]. England:Longman Group Ltd, 1995.

[5] Bowyer A. Computing dirichlet tessellations[J]. The Computer Journal, 1980,24(2):162-166.

[6] Chippada S, Dawson C N，Martinet M L, et al. Compute methods appL[J]. Mech Eng，1998,151(1-2):105-129.

[7] Rogers B M, Fujihara M, Borthwick A G L. Adaptive Q-tree Godunov-type scheme for shadow water equation[J]. International Journal for Numerical Methods in Fluids, 2001, 35(3):247-280.

[8] Pan C, Dai S, Chen S. Numerical simulation for 2D shallow water equations by using Godunov-type scheme with unstructured mesh[J]. Journal of Hydrodynamics, 2005,18(4):475-480.

[9] Yang Qin. Constrained Delaunay triangulation mesh generation technology[M]. Beijing:Publishing House of Electronics Industry, 2005.(in Chinese)

[10] Du Min. Unstructured triangulation grid generation and its application in 2D hydrodynamics model[D]. Tianjin:Tianjin University, 2005.(in Chinese)

[11] Zhang Zhi-yi, Konno K, Tokuyama Y. Curve mesh reconstruction based on mountain contours [J]. Journal of the Institute of Image Information and Television Engineers, 2006, 60(11):1803-1810.

[12] Zhao Di-hua, Yao Qi, Jiang Yan, et al. FVS scheme in 2D depth-averaged flow-pollutants modeling[J]. Advances In Water Science, 2002, 13(6):701-705.(in Chinese)

[13] Zhang Li-qiong,Cui Guang-bai,Yang Jue.Flux vector splitting and finite volume method in the application of hydraulic modeling[J]. Water Resources and Hydropower Engineering, 2001,8(32):8-11.(in Chinese)

[14] Roshandel A,Hedayat N,Kiamanesh H.Simulation of dam break using finite volume method[J]. World Academy of Science, Engineering and Technology, 2010,48(4):112-115.

[15] Spekerijse S P. Second order accurate upwind solutions of the 2D steady Euler equations by the use of a defect correction method[C]∥Proc of the 2nd European Conference on Multigrid Method, 1986:285-300.

[16] Steger J L, Warming R. Flux vector splitting of the inviscid gasdynamic equations with application to finite difference methods[J]. Journal of Computational Physics, 1981, 40(2):263-293.

[17] Pang Ming-yong. The real time dynamic simulation method of two dimensional water waves based on the discrete model[J]. Shuili Xuebao, 2007,38(11):1358-1363.(in Chinese)

[18] Gao Zhi-yi, Yu Fu-jiang, Xu Fu-xiang. Preparing 3-D animation of sea wave forecasting[J]. Marine Science Bulletin, 2011,30(2):173-178.(in Chinese)

[19] Li Yun-fei,Cheng Tian-tian,He Wei.Simulation method for lake surface wave[J]. Journal of System Simulation,2009,21(23):7507-7510.(in Chinese)

附中文參考文獻(xiàn)：

[1] 譚維炎. 計(jì)算淺水動力學(xué)—有限體積法的應(yīng)用[M]．北京：清華大學(xué)出版社, 1998.

[3] 李人憲. 有限體積法基礎(chǔ) [M]. 北京：國防工業(yè)出版社, 2008.

[9] 楊欽. 限定Delaunay三角網(wǎng)格剖分技術(shù)[M]. 北京：電子工業(yè)出版社,2005.

[10] 杜敏. 非結(jié)構(gòu)三角網(wǎng)格生成及其在二維水動力學(xué)模型中的應(yīng)用[D]. 天津：天津大學(xué), 2005.

[12] 趙棣華, 姚琪, 蔣艷, 等. 通量向量分裂格式的二維水流-水質(zhì)模擬 [J]. 水科學(xué)進(jìn)展, 2002, 13(6):701-705.

[13] 張麗瓊, 崔廣柏, 楊玨. 通量向量分裂格式及有限體積法在水流模擬中的應(yīng)用 [J]. 水利水電技術(shù),2001,8(32):8-11.

[17] 龐明勇. 基于離散模型的二維水波實(shí)時動態(tài)模擬方法 [J]. 水利學(xué)報(bào), 2007,38(11):1358-1363.

[18] 高志一，于福江，許富祥. 海浪預(yù)報(bào)三維動畫計(jì)算原理與制作方法 [J]. 海洋通報(bào), 2011,30(2):173-178.

[19] 李云飛，程甜甜，何偉. 一種湖面波浪模擬的方法 [J]. 系統(tǒng)仿真學(xué)報(bào), 2009,21(23):7507-7510.

ZHANG Zhi-yi,born in 1974,PhD,associate professor,his research interest includes computer graphics.

A method of real time wave simulation based on physical

ZHANG Zhi-yi,Lü Zhi-hua
(College of Information Engineering,Northwest A&F University,Yangling 712100,China)

In order to solve the real-time computation and visualization of the wave state in the wide area, by combining the physical model of non-viscous fluid dynamics, a simplified finite volume algorithm is proposed, which uses triangle as control unit. The advantage of the proposed algorithm is that: Based on irregular triangular mesh of irregular boundary area, the boundaries' numerical flux of triangles can be obtained by simplifying flux-vector splitting method. Then shallow water flow with irregular boundary can be real-time simulated by fast approaching the solution of the two-dimensional shallow water equations. Experimental results show that the real-time visual simulation of the wave can be better achieved by using this method in the physical laws of the real world under the premise.

shallow water equations;constrained Delaunay triangulation;finite volume method;flux-vector splitting method

2012-09-20;

2012-12-19

教育部留學(xué)回國人員科研啟動基金資助項(xiàng)目(K314020901)；中央高?；究蒲袠I(yè)務(wù)費(fèi)專項(xiàng)基金資助項(xiàng)目(Z109021004)

1007-130X(2014)03-0508-05

TP391.41

A

10.3969/j.issn.1007-130X.2014.03.023

張志毅(1974-),男,山西夏縣人，博士，副教授，研究方向?yàn)橛?jì)算機(jī)圖形學(xué)。E-mail:815802490@qq.com

通信地址：712100 陜西省楊凌市西農(nóng)路22號西北農(nóng)林科技大學(xué)信息工程學(xué)院

Address:College of Information Engineering,Northwest A&F University,22 Xinong Rd,Yangling 712100,Shaanxi,P.R.China