基于ＯｐｅｎＧＬ的流體交互式仿真

摘 要：為避免二維波方程中的大量計(jì)算，使用二維波動(dòng)方程的中心差分近似來模擬流體運(yùn)動(dòng)。在該方法上實(shí)現(xiàn)了基于物理的浮力和阻力模型，模擬物體在流體表面的漂浮。在OpenGL上實(shí)現(xiàn)上述方法，實(shí)驗(yàn)證明該方法實(shí)現(xiàn)簡單，大大降低了運(yùn)算代價(jià)。

關(guān)鍵詞：虛擬現(xiàn)實(shí)；流體；波動(dòng)；柵格； 仿真

中圖法分類號(hào)：TP391.9文獻(xiàn)標(biāo)識(shí)碼：A

文章編號(hào)：1001—3695(2007)02—0186—03

近幾年來，隨著計(jì)算機(jī)硬件技術(shù)和計(jì)算機(jī)圖形學(xué)、三維仿真、虛擬現(xiàn)實(shí)等技術(shù)的快速發(fā)展，實(shí)時(shí)模擬真實(shí)的室外環(huán)境成為可能，而動(dòng)態(tài)的流體對(duì)任何室外游戲場景都可以增加極大的美感。國內(nèi)外的學(xué)者對(duì)流體模擬進(jìn)行了許多研究，如Trim [1]將水面看作一個(gè)繃緊伸展的彈性膜，忽略其重力，采用偏微分和積分來模擬水面，但這種方法計(jì)算量大。Press將流體表面分隔為大小平均的柵格，將上述積分變化為離散操作，利用快速傅里葉變換(FFT)進(jìn)行計(jì)算，成功地簡化了上述方法，但在交互式應(yīng)用中，如果柵格較大，其計(jì)算代價(jià)也高得驚人[2]。采用紋理映射中對(duì)紋理矩陣操作的方法可以非常簡單地實(shí)現(xiàn)浮動(dòng)的海面[3]，但逼真度低且無法實(shí)時(shí)模擬物體在水面上的移動(dòng)情況。通過計(jì)算級(jí)數(shù)中重要的項(xiàng)也可以得到近似解，雖然該方法直接且穩(wěn)定，但效率很低。因?yàn)榇蠖鄶?shù)的CPU計(jì)算一個(gè)三角函數(shù)要用大約30個(gè)周期，即使沿x，y方向僅取三個(gè)頻率最大的值，每個(gè)時(shí)間間隔每一點(diǎn)也要計(jì)算九個(gè)正弦函數(shù)，對(duì)于交互應(yīng)用中較大的柵格其代價(jià)更高。為減少運(yùn)行的代價(jià)且方便實(shí)現(xiàn)與水面漂浮物體的交互，本文通過使用二維波動(dòng)方程的中心差分近似，使每個(gè)點(diǎn)僅用幾個(gè)算術(shù)運(yùn)算就能模擬流體表面的水平運(yùn)動(dòng)。在這種表面的基礎(chǔ)上，對(duì)漂浮物體進(jìn)行受力分析，實(shí)現(xiàn)了物體在水面中的交互式模擬。

1 流體表面的模擬實(shí)現(xiàn)

流體的仿真一般采用如下的二維波方程:

其中常數(shù)c代表每單位時(shí)間位移的距離，即流體的速度；μ為非負(fù)常數(shù)，代表流體的黏性，該值大體上決定了波在一種表面上傳播的距離，也就是波消失的時(shí)間。較小的μ值可以讓波存在時(shí)間較長，適合模擬水波；較大的μ值會(huì)導(dǎo)致波的快速消失，如黏稠的油波。

本文假設(shè)邊界條件是齊次的（即邊緣不作上下運(yùn)動(dòng)），且流體表面的初始速度為0，對(duì)于一個(gè)L×L大小的正方形流體表面用三角形網(wǎng)格連接，網(wǎng)格的頂點(diǎn)排列成平均分隔的柵格，如圖1所示。式(1)的二維波方程可以用分離變量的解析方法求解，但其求解過程十分復(fù)雜、效率低，在實(shí)時(shí)仿真時(shí)需要相當(dāng)大的計(jì)算量。為簡化運(yùn)算，采用近似導(dǎo)數(shù)方法模擬波在流體表面的傳播過程。

1.1 波動(dòng)的模擬偏導(dǎo)數(shù)實(shí)現(xiàn)

1.2 不穩(wěn)定性的限制

1.3 邊界的反射

水域通常不是正方形的，河流、湖泊和海洋有著不同斜度的不規(guī)則海岸線，而且島嶼也可能存在于水域內(nèi)部。如果水岸非常陡峭甚至是垂直的，波浪從海岸線反射回來就幾乎不會(huì)失去能量；然而如果水岸是傾斜的，波浪的反射波可能會(huì)非常弱，或根本沒有。如果波浪不是垂直擊向岸邊，那么就會(huì)以某個(gè)角度反射回來。可以通過一個(gè)本地阻尼系數(shù)d_{i， j}調(diào)節(jié)zⁿ⁺¹的值來模擬。系數(shù)為1表明允許在高度值上自由運(yùn)動(dòng)而不損失任何能量；系數(shù)為0代表限制在那個(gè)位置水域的所有運(yùn)動(dòng)。如果這些系數(shù)是根據(jù)地形特征來分配和調(diào)節(jié)的，波浪對(duì)海岸線做出的反應(yīng)就更自然。例如，如果海岸是陡峭的，阻尼系數(shù)就應(yīng)該從1（水）到0（陸地）作一個(gè)快速的轉(zhuǎn)變；如果海岸是逐漸傾斜的，阻尼系數(shù)就應(yīng)該從1漸變到0。實(shí)際上，將略小于1的阻尼系數(shù)用于“濕的”柵格比較好，可以產(chǎn)生較少的能量損耗；否則波動(dòng)會(huì)無限地延續(xù)下去。

1.4 算法的具體實(shí)現(xiàn)

為避免在時(shí)刻t_-1，t₀，t₁的Z值存儲(chǔ)三個(gè)柵格，將z_n-1被zn+1適當(dāng)代換，只需要兩個(gè)柵格， 通過在傳遞結(jié)束時(shí)，指向z_n+1和z_n的值的內(nèi)存指針進(jìn)行交換實(shí)現(xiàn)。因此開辟兩個(gè)緩沖區(qū)，每一個(gè)為包含n×m個(gè)頂點(diǎn)位置的數(shù)組。在顯示每一幀時(shí)，包含當(dāng)前頂點(diǎn)位置的緩沖區(qū)稱為當(dāng)前緩沖區(qū)，包含上一時(shí)間段的頂點(diǎn)位置的緩沖區(qū)稱為舊緩沖區(qū)。在計(jì)算頂點(diǎn)的新位移時(shí)，需要用新的頂點(diǎn)位移更新舊緩沖區(qū)中的頂點(diǎn)位置，則當(dāng)前緩沖區(qū)成為下一幀中的舊緩沖區(qū)，而舊緩沖區(qū)則成為下一幀中的當(dāng)前緩沖區(qū)。這樣，交替使用這兩個(gè)緩沖區(qū)進(jìn)行流體表面的顯示。

2 與水面交互

為真實(shí)模擬水面波動(dòng)，必須模擬水面波動(dòng)的開始（即水面的飛濺效果）和漂浮物體隨水面浮動(dòng)的效果。

2.1 飛濺

飛濺的效果可以通過在特定位置瞬時(shí)移動(dòng)一個(gè)或幾個(gè)Z值來創(chuàng)建。當(dāng)一個(gè)物體與流體表面接觸時(shí)（如投擲一塊磚頭）會(huì)引起流體表面的振動(dòng)，這時(shí)流體表面的位置可以通過修改接觸點(diǎn)周圍頂點(diǎn)的位移直接得到。將離接觸點(diǎn)最近的點(diǎn)壓縮為一點(diǎn)可以獲得滿意的模擬結(jié)果。實(shí)驗(yàn)證明，將離接觸點(diǎn)最近的八個(gè)相鄰頂點(diǎn)壓縮為一點(diǎn)較為合適。

2.2 物體受力分析

漂浮在流體表面的物體被分解成M個(gè)相互連接的部分，由于每個(gè)部分對(duì)于該物體的任何其他部分都是不可平移或旋轉(zhuǎn)的，并且物體的質(zhì)心固定，因此我們將該物體當(dāng)作剛體對(duì)待。漂浮物體所受到的浮力等于它排出的水的重量。浮力的方向?qū)嶋H上是在氣壓剃度的方向上，在大多數(shù)情況下，取水面的法線方向就可以了。

在本系統(tǒng)中，將漂浮的物體外形看作由一組離散的點(diǎn)近似構(gòu)成，相鄰的一組點(diǎn)構(gòu)成一個(gè)區(qū)域片ΔAk，區(qū)域片的法向量為k，z。浮力通過對(duì)物體水下部分的體積積分來計(jì)算。水下部分物體排出的水的體積為

其中zwater是在PK的雙線性插值水面高度（可避免初級(jí)或一級(jí)間斷），雙線性插值可以避免可能產(chǎn)生的初級(jí)或一級(jí)間斷。在該位置的浮力Fk和轉(zhuǎn)動(dòng)力矩Nk分別為

這里rk是從質(zhì)心到pk的矢量，ρ為流體的密度。合力與轉(zhuǎn)動(dòng)力矩可以通過對(duì)每一個(gè)物體頂點(diǎn)分量的求和來獲得。頂點(diǎn)的數(shù)目決定物體的形狀和浮力的精確度。一般情況下，一個(gè)立方體可能需要20或30個(gè)點(diǎn)，而一棵有分支的樹可能需要數(shù)百個(gè)點(diǎn)來真實(shí)表現(xiàn)。漂浮物體受力分析如圖3所示。

只受到重力和浮力的作用，即物體受力總和為浮力減去重力。為了保持物體在水面上的滑行，可以通過求每一個(gè)頂點(diǎn)的分量和來計(jì)算拉力：

速度項(xiàng)vk，rel是在rk處船體相對(duì)水的速度，它通過計(jì)算頂點(diǎn)b的水流速度vwater和剛體總速度（包括線速度vcm和角速度ωcm）的差值獲得。因此除高度值外，如果一個(gè)三維速度與每一個(gè)柵格區(qū)域相關(guān)聯(lián)，水流就載著物體漂浮移動(dòng)。

2.3 流體表面法向量計(jì)算

為了正確地反映流體的流動(dòng)情況，要計(jì)算流體表面法向量，需要知道在每個(gè)頂點(diǎn)上的準(zhǔn)確法向量和切向量。在頂點(diǎn)(i， j)處，與x軸平行的切向量T和與y軸平行的切向量B可以表示為

法向量N（非規(guī)則化向量）可通過計(jì)算N=T×B得到，適當(dāng)放大（不會(huì)改變向量的方向）以消除分母，得到的最終表達(dá)式為

3 實(shí)驗(yàn)及其結(jié)論

根據(jù)上述算法，我們?cè)赩isual C++6.0平臺(tái)上，采用OpenGL進(jìn)行流體的模擬，實(shí)驗(yàn)?zāi)M對(duì)象為水面，大小為128×128，水面的阻力系數(shù)值μ=0.001，以保證波的震蕩時(shí)間較長。圖4為開始時(shí)流體平靜的狀態(tài)，物體靜止在水面上。圖5 為水面中剛引入飛濺時(shí)的效果。圖6、圖7為飛濺后不同時(shí)期的水面震蕩效果對(duì)比。其中圖6表明水波剛傳播到漂浮的物體處，引起物體隨著水波上下移動(dòng)；圖7表明物體被水波推動(dòng)，引起水平的移動(dòng)，并隔斷水波的傳遞，使其無法傳播到物體身后的區(qū)域。

本系統(tǒng)中假設(shè)水域的邊界是垂直的，波浪從邊界反射回來幾乎不損失能量。圖8是飛濺后，水波經(jīng)過邊界多次反射和疊加后的整個(gè)場景效果。

本文中所涉及到的圖表、注解、公式等內(nèi)容請(qǐng)以PDF格式閱讀原文。