計(jì)算流體力學(xué)三維非結(jié)構(gòu)網(wǎng)格可視化的高效算法研究

馬志強(qiáng)，多 佳，李 曼

(1.大連理工大學(xué)土木水利學(xué)院，遼寧 大連 116024;2.黑龍江科技學(xué)院數(shù)力系，哈爾濱 150027;3.黑龍江省水利水電勘測設(shè)計(jì)研究院，哈爾濱 150080)

1 引言

隨著有限元(FEM)和有限體積法(FVM)等離散方法的發(fā)展和計(jì)算機(jī)自身運(yùn)算能力的提高，三維流體力學(xué)數(shù)值模擬越來越多的得到應(yīng)用，由于對復(fù)雜地形的適應(yīng)性較強(qiáng)，因此大部分三維流體力學(xué)數(shù)值模擬都采用非結(jié)構(gòu)網(wǎng)格。針對這種數(shù)據(jù)結(jié)構(gòu)，如何把計(jì)算結(jié)果高效快速的動態(tài)顯示出來是計(jì)算流體力學(xué)科學(xué)可視化的面臨的一個(gè)問題[1]。

由于在一般的三維可視化過程中需要顯示的只是所選擇的計(jì)算區(qū)域外表面的信息，因此三維非結(jié)構(gòu)化網(wǎng)格數(shù)值模擬計(jì)算結(jié)果可視化快速實(shí)現(xiàn)問題就轉(zhuǎn)化為快速的判斷在三維網(wǎng)格中哪些組成單元的面是計(jì)算域外表面的問題。

圖1 常見的網(wǎng)格形式

2 非結(jié)構(gòu)化網(wǎng)格的基本定義

在說明算法以前我們先有如下的定義:

定義1:F為在區(qū)域D中面的集合，對于任意一個(gè)面f∈F，如果f同屬于兩個(gè)不同的單元，則f屬于內(nèi)部面，記為I。如果f只屬于唯一單元，則f屬于外部面，記為Г。I∪Г=F。

定義 2:f∈F，v1，v2，……，vL是面 f的頂點(diǎn)，記為 f=如果f的頂點(diǎn)的排列符合順時(shí)針排列，f稱為順時(shí)針面。

3 兩種判斷外邊界的算法

3.1 傳統(tǒng)算法

傳統(tǒng)搜索外邊界面的算法是對所有單元的面進(jìn)行提取，進(jìn)而按所提取的面進(jìn)行循環(huán)，搜索出循環(huán)到的面是否在已提取的集合中有相同的面，若有重復(fù)則說明此面是內(nèi)部面，否則此面為外部面[2]。

步驟1:

從每個(gè)單元中提取面元，使每個(gè)面頂點(diǎn)的排序滿足定義2，Sj為單元ej所包含的面的個(gè)數(shù)，通過對一個(gè)單元的面元的提取，形成單元ej中面的集合Fj={f(1)j，f(2)j，……fsjj}，而每個(gè)面記為 f(s)j={njm1，njm2，……，njnks}。其中 j=1，2，……，Me;njmk∈ej?Np(Np為所有點(diǎn)的集合)，s=1，2，……，S，k=1，2，……，Ks。Ks為平面f(s)j的頂點(diǎn)數(shù)。非結(jié)構(gòu)化網(wǎng)格中所有單元所包含的面記為F={Fj|j=1，2，……，Me}。

步驟2:

逐個(gè)對比F中的元素，區(qū)分出外部面集合Г和內(nèi)部面集合I。因?yàn)橥粋€(gè)單元沒有重復(fù)面，若一個(gè)面在集合F中出現(xiàn)兩次，則說明為兩個(gè)單元共有，為內(nèi)部面。若只出現(xiàn)一次，則為外部面。此過程可用以下方程表示:

對于 1≤j≤Me-1，1≤k≤sj，如果，那么就存在 i0＞ j，1≤s0≤si0，滿足，因此反之，如果，那么

當(dāng)?shù)诙浇Y(jié)束后，在非空集合F中就提取出外表面的集合Г。

3.2 改進(jìn)算法

改進(jìn)后的算法是按節(jié)點(diǎn)對單元的面進(jìn)行提取，找出最小(或最大)頂點(diǎn)序號相同的一組面元進(jìn)行提取形成集合，并判斷這些集合中的面元是否有重復(fù)，若有重復(fù)則說明此面是內(nèi)部面，否則此面為外部面。

步驟1:

找出與結(jié)點(diǎn)ni(ni∈Np)相關(guān)聯(lián)的面，且ni是這個(gè)面結(jié)點(diǎn)序號最小的值，這些符合條件的面形成矩陣Eni，Eni中的每個(gè)元素是滿足以上條件的面除ni外其它結(jié)點(diǎn)序號的集合。

步驟2:

逐個(gè)對比Eni中的元素，定義，如果有，那么就存在 j0≠k(j0＞ k)，使，那么}一定是內(nèi)部面，保存在集合I中。如果且，那么平面{ni，nl1，nl2，……，nlk}一定是外部面，保存在集合Г中。

步驟3:

k=k+1，如果k＜kni，那么重新執(zhí)行第二步，并且把內(nèi)部面和外部面分別記錄在集合I和Г中。如果k=kni，那么執(zhí)行第四步。

步驟4:

i=i+1，如果i≤N，那么重復(fù)執(zhí)行1-3步，把新判斷的面加到集合I和Г。如果i＞N，則判斷結(jié)束，其中N為Np中點(diǎn)的個(gè)數(shù)。

4 兩種算法的實(shí)例比較

傳統(tǒng)算法中對外邊界判斷計(jì)算次數(shù)是O(M2e)，而利用改進(jìn)算法需要的計(jì)算次數(shù)是O(N)，一般來說Me和N是同一個(gè)量級，因此改進(jìn)算法的計(jì)算次數(shù)要遠(yuǎn)遠(yuǎn)小于傳統(tǒng)算法。

我們應(yīng)用一個(gè)實(shí)際的例子來比較兩種算法的運(yùn)算速度:

鐵成尾礦庫隸屬于河北省承德市灤平縣小營鄉(xiāng)管轄，尾礦庫庫區(qū)最高海拔高程900 m，最低海拔高程410 m，總面積約為109萬m2。

計(jì)算區(qū)域網(wǎng)格由ansys劃分，為四面體、六面體混合網(wǎng)格，共160523個(gè)節(jié)點(diǎn)，739987個(gè)單元，圖2為計(jì)算網(wǎng)格局部圖。計(jì)算平臺為CPU P4 2.4G，1G內(nèi)存，編程環(huán)境為VS.Net 2005，表1為計(jì)算結(jié)果，結(jié)合自行開發(fā)的基于OpenGL的流體力學(xué)可視化軟件，采用Dividing Cube等值面填充的方法顯示計(jì)算區(qū)域的含水度，將結(jié)果動態(tài)顯示出來，見圖3。通過表1可知通過算法的改進(jìn)大大提高了顯示效率。

表1 各算法用時(shí)統(tǒng)計(jì)表

圖2 鐵成尾礦庫三維計(jì)算局部網(wǎng)格圖

圖3 鐵成尾礦庫含水度三維動態(tài)演示圖

5 總結(jié)

綜上，改進(jìn)算法較之傳統(tǒng)算法在大型三維計(jì)算的可視化顯示中更具效率，可實(shí)現(xiàn)對外邊界的快速判斷，這種算法不但可以應(yīng)用到數(shù)值模擬可視化工作中來，也可以應(yīng)用到數(shù)值模擬的計(jì)算中，實(shí)現(xiàn)對邊界條件的快速判斷，還可擴(kuò)展到2D計(jì)算的工作中來，有著很高的應(yīng)用價(jià)值。

[1] 劉曉波，華祖林，何國建.計(jì)算流體力學(xué)的科學(xué)計(jì)算可視化研究進(jìn)展[J].水動力學(xué)研究與進(jìn)展，2004，(1):120-125.

[2] 唐澤圣.三維數(shù)據(jù)場可視化[M].北京:清華大學(xué)出版社，1999.