馬志強，多 佳，李 曼

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

1 引言

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

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

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

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

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

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

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

3 兩種判斷外邊界的算法

3.1 傳統(tǒng)算法

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

步驟1:

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

步驟2:

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

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

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

3.2 改進算法

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

步驟1:

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

步驟2:

逐個對比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中點的個數(shù)。

4 兩種算法的實例比較

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

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

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

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

表1 各算法用時統(tǒng)計表

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

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

5 總結(jié)

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

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

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