本發(fā)明涉及計(jì)算機(jī)圖形學(xué)領(lǐng)域，具體涉及一種參數(shù)優(yōu)化的淺水方程模型水體建模方法。

背景技術(shù)：

復(fù)雜流體模擬技術(shù)一直以來都是計(jì)算機(jī)圖形學(xué)中一個熱門的話題。在計(jì)算機(jī)流體模擬的眾多方法中，基于物理的流體模擬使用流體力學(xué)方程來描述水體內(nèi)部速度、壓力、密度等物理量隨時間的變化，從而能夠真實(shí)呈現(xiàn)水體波動以及其他細(xì)節(jié)。

在描述流體動力學(xué)方程方面有著名的納維-斯托克斯方程(參考文件1：m.griebel,t.dornsheifer,t.neunhoeffer.numericalsimulationinfluiddynamics:apracticalintroduction[m].siam:societyforindustrialandappliedmathematics,1997.)，簡稱n-s方程。但n-s方程是一個非線性偏微分方程，在求解過程中需要大量的計(jì)算節(jié)點(diǎn)，導(dǎo)致計(jì)算復(fù)雜度非常高，因此使用傳統(tǒng)的n-s方程來對水體進(jìn)行建模，從計(jì)算機(jī)應(yīng)用的角度出發(fā)是不可行的。但是可以通過數(shù)值推導(dǎo)來簡化n-s方程的復(fù)雜度，從而進(jìn)行相應(yīng)的求解。從這個角度出發(fā)，將水體占據(jù)的計(jì)算空間按一定規(guī)則分為網(wǎng)格，然后分析空間中每一個網(wǎng)格位置上的水體速度、壓強(qiáng)、密度等參數(shù)隨時間和空間的變化，從而衍生出基于歐拉網(wǎng)格的淺水方程方法(shallowwaterequation)，簡稱swe方法。

hagen等人(參考文件2：t.r.hagen,j.m.hjelmervik,k.a.lieetal.visualsimulationofshallow-waterwaves[j].simulationmodellingpracticeandtheory,2005,13(8):716-726.)第一次將swe引入到計(jì)算機(jī)圖形學(xué)中，利用有限元方法在gpu(graphicsprocessingunit，圖形處理器)上求解swe，進(jìn)而模擬河流、海浪沖刷不規(guī)則地形的效果。利用高度場的方法，無法處理破碎的波浪效果，因此，thury和müller等人(參考文件3：n.thurey,m.muller-fischer,s.schirmetal.real-timebreakingwavesforshallowwatersimulations[a].computergraphicsandapplications,2007.pg'07.15thpacificconferenceon[c],2007:39-46)提出了一種能夠自動生成細(xì)化三角網(wǎng)格補(bǔ)丁的方法來模擬海水沖擊岸邊后產(chǎn)生的破碎效果，但是該算法不能很好地進(jìn)行并行計(jì)算，從而難以達(dá)到實(shí)時模擬。chentanez和müller(參考文件4：n.chentanez,m.müller.real-timesimulationoflargebodiesofwaterwithsmallscaledetails[a].proceedingsofthe2010acmsiggraph/eurographicssymposiumoncomputeranimation[c],2010:197-206)調(diào)研整合了之前的swe算法，并描述了兩種水固交互以及干濕跟蹤的方法，這兩種方法也在swe的網(wǎng)格模型的效率和穩(wěn)定性上進(jìn)行了進(jìn)一步的加強(qiáng)，其中還詳細(xì)介紹了如何使用swe構(gòu)建一個復(fù)雜水體場景，包括對海浪拍打島嶼岸邊、浪花飛濺、泡沫、漩渦等效果，并在其中實(shí)現(xiàn)了瀑布效果。但是，上述基于傳統(tǒng)淺水方程的水體建模都是建立在高度域的坐 標(biāo)軸方向固定為重力方向的基礎(chǔ)上的，因此當(dāng)與水體交互的固體表面過于陡峭時便會出現(xiàn)失真的現(xiàn)象，無法準(zhǔn)確描述水體特征。其次，傳統(tǒng)的淺水方程并沒有描述外作用力之間是如何共同影響水流變化的，包括水流表面的張力和人為的控制力，同時對于任意的外作用力并沒有給出隱式的描述，這樣會導(dǎo)致不能逼真地對水體進(jìn)行建模。

技術(shù)實(shí)現(xiàn)要素：

本發(fā)明的目的是為了解決上述現(xiàn)有傳統(tǒng)淺水方程模型水體建模方法的無法逼真模擬水體，以及計(jì)算過程中高數(shù)值耗散帶來的計(jì)算速度慢等缺點(diǎn)。本發(fā)明提出一種參數(shù)優(yōu)化的淺水方程模型水體建模方法，旨在提高水體建模過程中的計(jì)算效率和水體模擬效果的逼真度。

本發(fā)明提出的一種參數(shù)優(yōu)化的淺水方程模型水體建模方法，具體包括如下步驟：

步驟一：對n-s方程進(jìn)行數(shù)值推導(dǎo)和條件限制，推導(dǎo)出參數(shù)優(yōu)化的淺水方程模型。

淺水方程中的坐標(biāo)系包括，高度域空間坐標(biāo)軸以及平面空間坐標(biāo)軸。首先，沿著水體與固體接觸面的法向量方向建立高度域的空間坐標(biāo)軸，即h方向；然后，以垂直于該坐標(biāo)軸的平面建立平面的空間坐標(biāo)軸，即位置x方向，從而構(gòu)建淺水方程中使用的坐標(biāo)系統(tǒng)。參數(shù)優(yōu)化的淺水方程由式(1)和式(2)確定：

其中，u(x,t)是t時刻水體速度在平面x處的水平分量，ρ為水體的密度，p為外作用力(包括空氣壓力、表面張力以及重力沿h軸的分量)產(chǎn)生的壓強(qiáng)沿h軸的分量，a為外作用力在x平面方向產(chǎn)生的水平加速度，h(x,t)為t時刻位置x處的高度域。d＝h-b，其中d為水體的深度，b為地勢的高度。u和v為速度u(x,t)在x平面處的分量。符號為空間偏導(dǎo)，在二維曲面中運(yùn)算符“·”表示向量相乘。式中速度域u和高度域h為待求參數(shù)。

上述淺水方程模型的參數(shù)優(yōu)化主要體現(xiàn)在對n-s方程依據(jù)動量守恒和質(zhì)量守恒進(jìn)行數(shù)值推導(dǎo)時的過程簡化和形式變換。根據(jù)動量守恒定理優(yōu)化水體速度域偏微分方程時，將式(1)等號右邊部分優(yōu)化為三個組成部分：水體在x平面方向的平流項(xiàng)外作用力在h軸方向?qū)λw的影響以及外作用力在x平面方向?qū)λw的影響。其中，外作用力在h軸方向?qū)λw的影響由項(xiàng)決定，外作用力在x平面方向?qū)λw的影響由x平面方向的水流加速度a決定。

根據(jù)質(zhì)量守恒定理優(yōu)化水體高度域偏微方程時，將式(2)等號右邊部分優(yōu)化為兩個部分：水體在h軸方向的平流項(xiàng)和水流速度對高度域的影響。其中，水流速度對高度域的影 響由(u+v)d項(xiàng)決定。在水體運(yùn)動過程中，物理量如速度u和高度h是隨著時間和空間變化的，即速度和高度是時間t和空間x的函數(shù)，在式(1)和式(2)中分別表現(xiàn)為和統(tǒng)稱為平流項(xiàng)。平流項(xiàng)的存在，使得淺水方程是非線性的，因此在求解過程中比較復(fù)雜。本發(fā)明對傳統(tǒng)淺水方程進(jìn)行上述參數(shù)優(yōu)化有利于建立求解模型，提高方程的求解效率。

步驟二：為步驟一構(gòu)建的參數(shù)優(yōu)化的淺水方程模型建立求解模型。

步驟2.1，離散化步驟一中的待解參數(shù)。設(shè)淺水方程的空間計(jì)算網(wǎng)格為△τ，網(wǎng)格大小為△x×△x，時間步長為△t，在網(wǎng)格點(diǎn)(i,j)處，高度域?yàn)閔i,j，水平速度分量分別為ui,j和vi,j，作如下形式離散化：

步驟2.2，初始化淺水方程中的各項(xiàng)參數(shù)，包括初始高度場h，速度場u＝(u,v)以及地形高度b等。本發(fā)明中各參數(shù)的初始化如下：u＝0，v＝0，h＝2，b＝1。這樣，水體的初始深度為：d＝h-b＝1。實(shí)際操作過程中可以根據(jù)需要修改各項(xiàng)參數(shù)的初始值。

步驟2.3，計(jì)算淺水方程中的平流項(xiàng)。采用反向路徑追蹤和特征線法來計(jì)算平流項(xiàng)和反向路徑跟蹤的思想是，要計(jì)算t時刻位置x處的粒子在新的時間間隔△t后的速度，首先定義一條由t時刻到t+△t時刻的路徑p(x,s)，其中s為時間間隔[t,t+△t]中的一個時間點(diǎn)，新的速度項(xiàng)wnew由△t時間前舊的速度項(xiàng)wold通過計(jì)算得來：

wnew(x)＝wold(p(x,-△t))(6)

特征線法是一種解偏微分方程比較成熟的方法，其基本思想是通過把偏微分方程轉(zhuǎn)化為兩組常微分方程，再對常微分方程進(jìn)行求解。兩組常微分方程中的一組用于定義特征線，另一組用以描述解沿給定特征線的變化。這樣，結(jié)合反向路徑跟蹤和特征線法便可計(jì)算出平流項(xiàng)，詳見具體實(shí)施方式中步驟2.3。

步驟2.4，計(jì)算方程中外作用力對流體的作用，更新流體的速度場和高度場。設(shè)初始速度項(xiàng)為w0，計(jì)算外作用力對流體作用后的速度項(xiàng)，首先計(jì)算外作用力的直接作用效果：

w1＝w0(x)+△tf(x,t)(7)

其中，f(x,t)為外作用力在平面x處的分量，由加速度項(xiàng)a決定；w1為外作用力作用于水體后水體的速度。

然后對求得的w1進(jìn)行傅立葉變換：

對變換后的項(xiàng)進(jìn)行投影：

對投影后的項(xiàng)進(jìn)行傅立葉反變換，得到最后的結(jié)果：

上述計(jì)算所得結(jié)果w2即為新的速度域，將該新的速度域代入式(2)中并結(jié)合式(3)可以求出新的高度域。最后，用新的速度域和高度域更新舊的速度域和高度域，準(zhǔn)備下一輪的迭代求解。

步驟2.5，進(jìn)行邊界條件檢查，越過閾值處采用邊界截?cái)喾?。在?jì)算平流項(xiàng)時，反向路徑追蹤到的網(wǎng)格點(diǎn)可能超過定義的網(wǎng)格閾值范圍產(chǎn)生越界。在處理這種情況時，當(dāng)粒子距離邊界的距離小于半個網(wǎng)格，甚至超過了網(wǎng)格，設(shè)置其位置為距離邊界半個網(wǎng)格距離的位置。并根據(jù)上述步驟中計(jì)算出的速度域和高度域的值，調(diào)用繪圖程序繪制在當(dāng)前時間t和空間位置x處的水體。在模擬時間內(nèi)跳轉(zhuǎn)到步驟2.3，處理下一時間片。

步驟三：利用步驟二所建立的參數(shù)優(yōu)化的淺水方程模型構(gòu)建水體場景，輸出結(jié)果并顯示。步驟二中建立的水體模型，僅僅是基于數(shù)學(xué)物理求解得出的結(jié)果，從視覺上描述，是由一系列網(wǎng)格和粒子構(gòu)成，并不能逼真地呈現(xiàn)現(xiàn)實(shí)中的水體場景，因此還需要對上述的模型進(jìn)行相應(yīng)的渲染操作，才能真實(shí)呈現(xiàn)水體模型。本發(fā)明中采用opengl提供的庫函數(shù)進(jìn)行編程，對水體進(jìn)行渲染，從而展現(xiàn)逼真的水體效果。

本發(fā)明方法的優(yōu)點(diǎn)和積極效果在于：本發(fā)明在建立參數(shù)優(yōu)化的淺水方程模型過程中，采用動態(tài)變化的坐標(biāo)系統(tǒng)，即h坐標(biāo)的方向不是固定為重力方向，而是根據(jù)水體與固體接觸面的法向量來決定，提高了模型的穩(wěn)定性，減輕了由于接觸面過于陡峭帶來的破碎效果；然后在求解方程的過程中采用向量投影與傅立葉變換相結(jié)合的方式計(jì)算外作用力的影響，以及采用反向路徑跟蹤和特征線法相結(jié)合的方式計(jì)算平流項(xiàng)，提高了方程的求解速率，同時更加逼真地模擬了水體場景。

附圖說明

圖1是本發(fā)明參數(shù)優(yōu)化的淺水方程模型水體建模方法的整體步驟流程圖；

圖2是本發(fā)明三維動態(tài)坐標(biāo)系統(tǒng)圖；

圖3是本發(fā)明物理模型示意圖；

圖4是本發(fā)明反向路徑跟蹤示意圖。

具體實(shí)施方式

下面將結(jié)合附圖和實(shí)施例對本發(fā)明的技術(shù)方案作進(jìn)一步的詳細(xì)說明。

本發(fā)明提出一種參數(shù)優(yōu)化的淺水方程模型水體建模方法，通過在建立淺水方程模型過程中，采用動態(tài)坐標(biāo)系統(tǒng)建模，提高了模型的穩(wěn)定性；然后在求解方程的過程中，對于平流項(xiàng)的計(jì)算采用反向路徑跟蹤和特征線法相結(jié)合的方式，對于外作用力的計(jì)算采用向量投影與傅立葉變換相結(jié)合的方式，提高了模型的計(jì)算效率和性能。

本發(fā)明提出一種參數(shù)優(yōu)化的淺水方程模型水體建模方法，旨在提高傳統(tǒng)淺水方程模型對水體進(jìn)行建模的性能，下面具體對實(shí)現(xiàn)步驟進(jìn)行說明。

步驟一：對n-s方程進(jìn)行數(shù)值推導(dǎo)和條件限制，推導(dǎo)出參數(shù)優(yōu)化的淺水方程模型。

步驟1.1，沿著水體與固體表面接觸面的法向量方向建立一維高度域坐標(biāo)系h，再以該高度域坐標(biāo)系為基礎(chǔ)建立二維平面坐標(biāo)系x，從而構(gòu)建本模型的三維動態(tài)變化的坐標(biāo)系統(tǒng)。

步驟1.2，對水體建模過程中的所有場元素，包括高度場、速度場以及外作用力等進(jìn)行物理建模。圖2三維動態(tài)坐標(biāo)系統(tǒng)圖中展示了各個場元素的物理模型。為了避免水體之間的相互碰撞，首先，在一個低分辨率接觸面處計(jì)算其平均曲面法線的值；然后，定義地勢b(x)為接觸面原始曲面法線與平均曲面法線的差值。這樣根據(jù)水體的高度域h(x,t)函數(shù)就可以計(jì)算出水體深度為d＝h-b。定義u(x,t)為水平速度(這里的水平并不是傳統(tǒng)意義上的水平面，而是平面x方向)，垂直速度由隱式給出。最后對外作用力進(jìn)行物理建模，外作用力包括重力、大氣壓力、表面張力以及人為控制力等。將外作用力分解為h軸方向的一維作用力p和x平面方向的二維作用力f(x,t)，其所產(chǎn)生的水平加速度為a。

步驟1.3，根據(jù)步驟1.2中所建立的物理模型，將三維的n-s方程限制到二維曲面，來構(gòu)建本發(fā)明的參數(shù)優(yōu)化的淺水方程，如上述所構(gòu)建的式(1)和式(2)。其中式(1)由動量守恒推導(dǎo)而來，外作用力對水體的作用以加速度的形式體現(xiàn)，通過式(1)來更新水體的速度u；式(2)由質(zhì)量守恒推導(dǎo)而來，通過式(2)來更新水體的高度域h。

本發(fā)明步驟一中的三維動態(tài)坐標(biāo)系統(tǒng)與傳統(tǒng)的淺水方程坐標(biāo)系統(tǒng)相比，可以避免因水體與固體接觸面過于陡峭帶來的水體破碎等不穩(wěn)定的弊端，具有更好的穩(wěn)定性；同時，對水體的物理建模以及從三維n-s方程推導(dǎo)至二維的淺水方程，與傳統(tǒng)淺水方程相比具有更高的求解效率。

步驟二：為步驟一構(gòu)建的參數(shù)優(yōu)化的淺水方程模型建立求解模型。

步驟2.1，離散化步驟一中建立的物理模型，轉(zhuǎn)化為易于求解的數(shù)學(xué)模型。設(shè)該物理模型所使用的空間網(wǎng)格為△τ，每個網(wǎng)格的大小為△x×△x，每一個時間片的大小為△t，在網(wǎng)格在網(wǎng)格點(diǎn)(i,j)處，高度域?yàn)閔i,j，水平速度ui,j分量分別為ui,j和vi,j。在空間上進(jìn)行離散，將連續(xù)變化的流體離散為由一個個小網(wǎng)格構(gòu)成，對任意場的空間微分，只需求得與相鄰值的中央差分即可。在時間上進(jìn)行離散，將連續(xù)的的時間離散為由一個個時間片構(gòu)成，只需在每個時間片上計(jì)算流體運(yùn)動即可。因此，離散結(jié)果由式(3)～(5)給出。

步驟2.2，初始化方程中的各項(xiàng)參數(shù)，包括初始高度場，速度場以及地形高度等。本發(fā)明中各參數(shù)的初始化如下：u＝0，v＝0，h＝2，b＝1。這樣，水體的初始深度為：d＝h-b＝1。實(shí)際操作過程中可以根據(jù)需要修改各項(xiàng)參數(shù)的初始值。

步驟2.3，計(jì)算方程中的平流項(xiàng)。采用反向路徑追蹤和特征線法來計(jì)算平流項(xiàng)和圖3反向路徑跟蹤示意圖表明了這個過程。特征線法可以求解類似如下形式的平流問題：

其中為標(biāo)量，初始條件為定義向量u的特征線n(x0,t)是由位置x0處的粒子從時間t＝0開始平移而來，則有：

其中，初始條件為n(x0,0)＝x0。

然后，定義矢量是由特征線從初始狀態(tài)x0，t＝0平移而來。這樣可以通過下面的微分方程進(jìn)行計(jì)算：

其中特別的有所以在計(jì)算t時刻位置x處的值時，首先沿著特征線反向追蹤位置x處先前的位置，直到位置x0為止，然后初始化該位置的值：

根據(jù)上述特征線法可以計(jì)算本模型中的平流項(xiàng)，如上述式(6)。

步驟2.4，計(jì)算方程中外作用力對流體的作用，更新流體的速度場和高度場。該步驟主要針對式(2)動量守恒形式。

著名的helmholtz-hodege分解可以將一個任意的向量η分解成如下形式：

其中，向量ξ具有零散度，即a是標(biāo)量。這個分解表明任何一個向量可以由一個質(zhì)量守恒的域和一個標(biāo)量域的和組成。因此，可以定義一個算子k將任意向量η投影到一個零散度的向量ξ＝κη，于是有：

將上式代入到式(2)的等式兩邊可得：

其中ku＝u，并且根據(jù)該分解可以建立本發(fā)明中的投影模型。

將本模型中的邊界條件限制在周期邊界中，利用如下傅立葉變換進(jìn)行相應(yīng)的計(jì)算：

其中k＝|k|。

因此，在計(jì)算外作用力對水體的作用以及更新速度場和高度域時，首先計(jì)算力的直接作用效果如式(7)，然后對求得的項(xiàng)進(jìn)行傅立葉變換如式(8)，接著采用helmholtz-hodege分解投影如式(9)，最后對投影后的項(xiàng)進(jìn)行傅立葉反變換求得最終結(jié)果，并進(jìn)行相應(yīng)的域更新。

步驟2.5，進(jìn)行邊界條件檢查，越過閾值處采用邊界截?cái)喾?。在?jì)算平流項(xiàng)時，反向路徑追蹤到的網(wǎng)格點(diǎn)可能超過定義的網(wǎng)格閾值范圍產(chǎn)生越界。在處理這種情況時，當(dāng)粒子距離邊界的距離小于半個網(wǎng)格，甚至超過了網(wǎng)格，設(shè)置其位置為距離邊界半個網(wǎng)格距離的位置。并根據(jù)上述步驟中計(jì)算出的速度域和高度域的值，調(diào)用繪圖程序繪制在當(dāng)前時間t和空間位置x處的水體。在模擬時間內(nèi)跳轉(zhuǎn)到步驟2.3，處理下一時間片。

步驟二建立的計(jì)算模型通過將基于helmholtz-hodege分解的向量投影和傅立葉變換相結(jié)合來求解外作用力對水體的作用，同時將反向路徑跟蹤和特征線法相結(jié)合來求解平流項(xiàng)，能夠提高淺水方程的求解速度，從而提高水體建模的效率。

步驟三：利用步驟二所建立的參數(shù)優(yōu)化的淺水方程模型構(gòu)建水體場景，輸出結(jié)果并顯示。步驟二中建立的水體模型，僅僅是基于數(shù)學(xué)物理求解得出的結(jié)果，從視覺上描述，是由一系列網(wǎng)格和粒子構(gòu)成，并不能逼真地呈現(xiàn)現(xiàn)實(shí)中的水體場景，因此還需要對上述的模型進(jìn)行相應(yīng)的渲染操作，才能真實(shí)呈現(xiàn)水體模型。本發(fā)明中采用opengl提供的庫函數(shù)進(jìn)行編程，對水體進(jìn)行渲染，從而展現(xiàn)逼真的水體效果。

表1展示了傳統(tǒng)的淺水方程模型水體建模方法與本發(fā)明提出的參數(shù)優(yōu)化的淺水方程模型水體建模方法的模擬速率對比表?？梢钥吹?，本發(fā)明提出的方法平均模擬速率為24.9幀/秒，相對于傳統(tǒng)淺水方程方法平均模擬速率提高了29.016％。

表1模擬速率比較

本發(fā)明所有測試都是在聯(lián)想y460筆記本電腦上進(jìn)行的，cpu型號是英特爾corei3m350@2.27ghz，顯示卡型號是atimobilityradeonhd5650(1024mb)。平均模擬速率是通過記錄多次模擬實(shí)驗(yàn)的結(jié)果，求平均值得來的，不同硬件環(huán)境所得結(jié)果可能有所不同。