唐 燁，解利軍，桂立業(yè)，何麗莎，鄭 耀

浙江大學(xué) 航空航天學(xué)院，杭州 310027

1 引言

流場可視化是科學(xué)可視化的經(jīng)典分支之一。通過直觀的圖像來表示抽象的流場數(shù)據(jù)[1]，既可以加快從數(shù)據(jù)中獲取信息的速度，又能讓不具備領(lǐng)域?qū)I(yè)知識的人理解這些信息。如今在航空航天、氣象、自動化等領(lǐng)域的科學(xué)研究過程中，流場可視化已經(jīng)成為不可或缺的技術(shù)。

流面可視化是流場可視化的一種方式，近幾年來越來越流行[2]，相比于常用的流線法，流面法更適合表現(xiàn)復(fù)雜的流場結(jié)構(gòu)[3]，并且展現(xiàn)更豐富的流場局部信息[4]。

1.1 相關(guān)概念

流場是一種矢量場。在空間坐標(biāo)系(x,y,z)中，任意點(diǎn)(x0,y0,z0)都有唯一確定的向量(Vx,Vy,Vz)與之對應(yīng)，就說在空間(x,y,z)范圍內(nèi)存在矢量場V。在流場中，該矢量場即速度場。

流線是一個(gè)無質(zhì)量的質(zhì)點(diǎn)（種子點(diǎn)）在某一穩(wěn)定流場中以某點(diǎn)為起點(diǎn)，在該流場中運(yùn)行的軌跡。流線上的每一點(diǎn)都與流場中該點(diǎn)對應(yīng)的速度相切。假設(shè)V(x)是一個(gè)三維向量場，在時(shí)刻t經(jīng)過點(diǎn)x0的流線為L(x0,t)，其中L(x0,t)滿足公式（1）：

類似于流線的概念，流面是某條種子線C的所有種子點(diǎn)以C中的位置為起點(diǎn)在穩(wěn)定流場中運(yùn)行的軌跡的集合，設(shè)該流面為S，S滿足公式（2），流面S中的點(diǎn)處處與流場中該點(diǎn)處的速度相切[5]。

其中S(s,t)是種子點(diǎn)C中s處的種子點(diǎn)在流場中的運(yùn)動軌跡，也就是由C中s處的種子點(diǎn)生成的流線，S為C中所有點(diǎn)生成的流線的集合，如圖1所示。

圖1 流面示意圖

1.2 相關(guān)工作

Hultquist首先在1992年提出了一種高效的流面生成算法——前沿推進(jìn)算法。其基本思想為相鄰的兩條流線之間形成流帶，每條流帶需要記錄兩個(gè)前沿位置點(diǎn)，算法通過推進(jìn)每條流帶的前沿位置的方式來構(gòu)造流面，其特別之處為前沿的個(gè)數(shù)并不是固定的，當(dāng)某條流帶前沿距離過大時(shí)可以通過加入新點(diǎn)分裂流帶來提升流面精確度，當(dāng)前沿過窄時(shí)可以考慮合并相鄰流帶來減少計(jì)算與存儲資源的消耗[5]。

前沿推進(jìn)法非常適用于流場中存在大量的流線分散與流線聚攏的情況，即流場中速度大小基本穩(wěn)定，但方向出現(xiàn)偏差的情況，并且也較易于實(shí)現(xiàn)。但當(dāng)流場中的速度大小分布不勻稱時(shí)，前沿推進(jìn)法就難以適用。為此Garth在2004年提出了一個(gè)基于Hultquist算法的新算法，其基本思想為，在前沿推進(jìn)時(shí)前進(jìn)的距離以弧長為參數(shù)而不是固定的時(shí)間參數(shù)，同時(shí)在分裂與合并流帶時(shí)考慮曲率，使得最后得出的圖像能夠更精確[6]。

Scheuermann提出了一種在四面體網(wǎng)格中生成流面的算法[7]，其基本思想為通過重心坐標(biāo)在每個(gè)四面體內(nèi)計(jì)算出規(guī)則表面，通過面上的線段組成最終流面。這種方法能夠有效地利用流場中的拓?fù)浣Y(jié)構(gòu)信息，但用這種方法最終生成的流面質(zhì)量十分依賴于網(wǎng)格的密度[7]。

與以上這些顯示計(jì)算曲面的算法不同的是，Wijk在1993年提出過一種通過定義網(wǎng)格邊緣的標(biāo)量函數(shù)的方式來隱式地計(jì)算曲面的算法[8]，但由于它無法將所有可能的曲面都計(jì)算出來，因此應(yīng)用范圍并不廣。

Schneider等人在2009年提出在構(gòu)造曲面的時(shí)候用施密特插值法取代線性插值法[9]，其計(jì)算所需的額外數(shù)據(jù)都可以在流線生成過程中通過計(jì)算得出，可以將曲面的精確度提升到4階。

Schulze等人提出過一種強(qiáng)制流面前沿線與流場中的速度垂直的方法，它可以在湍流區(qū)生成高質(zhì)量的流面[10]。

以上都是確定種子線情況下所用的流面生成算法，但是當(dāng)拿到流場數(shù)據(jù)時(shí)，無法確定在哪里布線能夠得到較好的結(jié)果。因此本文介紹一種種子線自動布局以及流面生成的算法。

本文的主要思想：流場中不同區(qū)域的流動情況并不一致。如圖2所示，以流線法表示，可能有直線區(qū)域、曲線區(qū)域、渦流區(qū)域?？梢酝ㄟ^聚類算法將這些區(qū)域進(jìn)行劃分。劃分后的每個(gè)聚類可以用一個(gè)流面表示。生成流面的方法為，首先從分區(qū)中心點(diǎn)出發(fā)，沿著流場的曲率場在該分區(qū)范圍內(nèi)積分得到種子線。選擇在曲率場中積分是因?yàn)椋海?）曲率場中任意點(diǎn)的方向與流場中該點(diǎn)處的速度方向垂直，而垂直于流場速度方向的種子線生成的流面能夠覆蓋較大的范圍。（2）沿著曲率場生成的種子線所生成的流面更能表現(xiàn)出流場的流動結(jié)構(gòu)。得到種子線后，通過前沿推進(jìn)法配合施密特插值法生成光滑流面。

圖2 流場不同區(qū)域示意圖

2 流場的分區(qū)和種子線布局

2.1 程序流程

該算法處理的對象是流場數(shù)據(jù)V(X)∈R3，其中V∈R3，X∈R3，X代表流場空間中點(diǎn)的坐標(biāo)，V(X)代表點(diǎn)X處的速度。

圖3為程序的流程圖，對于流場中每個(gè)點(diǎn)P，首先算出它的曲率|C|以及速度模的梯度G用于之后的聚類以及種子線生成，詳見2.2節(jié)。隨后可以用改進(jìn)版的K-means算法將流場中的點(diǎn)集劃分為K個(gè)類，詳見2.3節(jié)，其中K的值由用戶自己輸入。然后對于每一個(gè)類的聚類中心點(diǎn)，用四階龍格庫塔法在曲率場中進(jìn)行積分得到K條種子線，詳見第3章。對每條種子線，通過帶施密特插值的前沿推進(jìn)法生成并繪制流面，詳見第4章。如果需要，用戶可迭代該過程，修改參數(shù)再次進(jìn)行聚類與流面生成工作。

圖3 程序流程圖

2.2 流場中數(shù)據(jù)點(diǎn)間相似度計(jì)算

在三維空間點(diǎn)集中，兩點(diǎn)間的相似度一般用它們的歐式距離來表示，但在流場數(shù)據(jù)集中，每個(gè)點(diǎn)除了坐標(biāo)信息外，還有對應(yīng)的速度信息。速度有兩個(gè)基本屬性：大小和方向。但由于關(guān)注的是潛在的種子線位置，因此希望把種子線布置在能夠代表周圍流線（速度）發(fā)生變化的區(qū)域。因此，并不是直接對流場中不同點(diǎn)間的速度進(jìn)行比較，而要比較速度的變化情況。

同時(shí)，渦流位置也可以通過曲率值的大小來反映，而流場的結(jié)構(gòu)也與渦流的位置有關(guān)[11]。

在流場空間中某處布置種子點(diǎn)形成流線S，則通過計(jì)算可以得到S中任意點(diǎn)P處的曲率|C|[12]，曲率也可以量化流線的彎折程度，也就是速度方向的變化程度。其中

其中V為點(diǎn)P處的速度，

點(diǎn)P處速度大小的變化情況|G|可以用速度模的梯度來表示，其中

某點(diǎn)處速度的變化情況可以通過|C|與|G|的值來反映，則流場中點(diǎn)P1、P2間的距離Dis(P1,P2)可以定義為Dis(P1,P2)=|C1-C2|+|G1-G2|。

當(dāng)P1、P2兩點(diǎn)坐標(biāo)間的歐式距離較大，而它們速度的變化量又相同的時(shí)候，Dis(P1,P2)為0，在聚類時(shí)會被歸為一類，如果這兩點(diǎn)的速度變化量較大，則它們應(yīng)該代表兩個(gè)潛在的布線區(qū)域。因此判定相似度時(shí)需要考慮兩點(diǎn)間的歐式距離：Dis(P1,P2)=|C1-C2|+|G1-G2|+|X1-X2|。

流場數(shù)據(jù)中空間坐標(biāo)與速度值的取值范圍并不相同，因此相比于點(diǎn)坐標(biāo)間的距離值，曲率值或梯度值間的距離值可能會顯得很小。此時(shí)兩點(diǎn)間的相似度無異于兩點(diǎn)坐標(biāo)間的歐氏距離，因此需要將它們?nèi)繗w一化處理。

在對不同的流場數(shù)據(jù)進(jìn)行聚類時(shí)，三個(gè)參數(shù)所占的權(quán)重未必一樣，因此需要對三個(gè)參數(shù)進(jìn)行加權(quán)，最終兩點(diǎn)間相似度Dis(P1,P2)表示為：

其中，α為兩點(diǎn)間相似度中二者的歐式距離所占的份額，β為速度變化量的兩個(gè)參數(shù)中曲率所占的份額，0<α<1，0<β<1。

2.3 基于改進(jìn)的K-means算法進(jìn)行流場聚類

聚類是將對象的集合按照一定的標(biāo)準(zhǔn)分割成多個(gè)相似子對象集合的過程，同一類中的對象間相似度較高，不同類間對象間相似度較低。

聚類算法可以分為兩大類：基于層次與基于劃分。在此基礎(chǔ)上還有三小類[13]：基于密度、基于模型、基于網(wǎng)格。

K-means算法屬于基于劃分的聚類法，與其他方法相比，有著實(shí)現(xiàn)簡潔與運(yùn)行快速兩大優(yōu)點(diǎn)。因此本文選擇該方法進(jìn)行流場聚類。

傳統(tǒng)K-means算法中初始聚類中心的選取完全是隨機(jī)的，而經(jīng)過實(shí)驗(yàn)發(fā)現(xiàn)K-means算法的效果會受到初始聚類中心的影響。選取不當(dāng)?shù)某跏季垲愔行臅咕垲愡^程更加漫長[14]。

因此，針對K-means算法在初始聚類中心選取上的盲區(qū)進(jìn)行一定的改進(jìn)。如果K個(gè)初始聚類中心互相之間具有較高的相似度，那么它們最后很有可能屬于同一個(gè)聚類，這會影響后續(xù)聚類的時(shí)間，于是希望K個(gè)初始聚類中心間的相似度盡可能得低。其過程如下：

（1）在數(shù)據(jù)集中隨機(jī)選擇一個(gè)初始種子點(diǎn)。

（2）對于數(shù)據(jù)集中所有點(diǎn)P，記錄其與聚類中心集合中的點(diǎn)Ci的距離的最小值D(P)。

（3）選擇一個(gè)新點(diǎn)加入到聚類中心集合中，選取的原則為D(P)值越大，被選中的概率也越大。

（4）如果得到K個(gè)聚類中心，執(zhí)行下一步，否則，返回到第（2）步。

（5）對于數(shù)據(jù)集中所有點(diǎn)，將其劃分到與它相似度最高（距離最?。┑木垲愔行乃鶎俚木垲愔?。

（6）對于每一個(gè)聚類，選擇到類中所有點(diǎn)的距離和最低的點(diǎn)作為新的聚類中心。

（7）如果本次迭代得到的K個(gè)聚類中心與上次迭代得到的K個(gè)聚類中心相同，則聚類情況已經(jīng)穩(wěn)定，結(jié)束程序，否則，返回到第（6）步。

實(shí)驗(yàn)表明，采用改進(jìn)版的K-means算法能夠在參數(shù)一致的前提下，減少聚類過程中的迭代次數(shù)，減少聚類所需時(shí)間，如表1所示。

表1 兩種K-means算法比較

聚類效果如圖4所示，其中圖4（a）為 α=0.5、β=0.3、k=2時(shí)的聚類效果。圖4（b）為 α=0.8、β =0.8、k=4時(shí)的聚類效果，圖4（c）為 α=0.5、β =0.8、k=6時(shí)的聚類效果。

圖4 聚類效果圖

3 種子線生成

得到了流場的聚類與中心數(shù)據(jù)后，就可以著手種子線的生成工作，流場的聚類中心雖然能代表流場的聚類，但單憑一個(gè)點(diǎn)難以確定種子線，還需確定種子線生成的方向與長度。

通過觀察可以得知，與速度方向平行的種子線無法生成流面，與速度方向垂直的種子線生成的流面能覆蓋較大范圍[15]，如圖5所示，圖中藍(lán)色帶箭頭曲線代表流線方向，x為垂直流場方向種子線，紅色曲線為不與流場垂直的種子線。二者生成同一流面，但種子線X更短。

圖5 種子線示意圖

種子線的朝向需要考慮到流場局部位置的結(jié)構(gòu)信息，否則，生成流面效果將大打折扣。如果流面處處與曲率方向相切并且平行于旋轉(zhuǎn)軸，那么它可以表現(xiàn)出螺旋流的特性。如果該面垂直于曲率方向和旋轉(zhuǎn)軸，該流面就難以傳達(dá)出足夠的有效信息。如圖6所示，圖6（a）為以垂直于速度與曲率方向生成的種子線所形成的流面，圖6（b）為垂直于速度方向平行于曲率向量方向積分形成的種子線生成的流面，可以看到圖6（b）中的內(nèi)容包含了圖6（a）的主要內(nèi)容。

圖6 種子線垂直曲率與平行曲率效果圖

根據(jù)流線中的曲率定義，曲率向量C由流線的速度向量V與加速度向量a通過叉積得出，因此C向量的方向一定是與V垂直，即給定一個(gè)初始點(diǎn)P，在常流場生成曲率場中的運(yùn)動軌跡必定處處與流場中該點(diǎn)處的速度垂直。因此，最終種子線生成方法是：以聚類中心點(diǎn)P為初始點(diǎn)，在曲率場中以生成流線的方法通過積分生成種子線。

由于每個(gè)聚類中心只能代表各自所屬的聚類。因此，種子線沿著曲率場積分到達(dá)分區(qū)邊界時(shí)，停止積分。

積分使用龍格庫塔法（Runge_Kutta）來實(shí)現(xiàn)。

四階龍格庫塔法為：

種子線效果如圖7所示，其中圖7（a）為α=0.5、β=0.3、k=2時(shí)生成的種子線，圖7（b）為 α =0.8、β =0.8、k=4時(shí)生成的種子線，圖7（c）為 α=0.5、β =0.8、k=6時(shí)生成的種子線。

圖7 種子線效果示意圖

4 流面生成

得到種子線后就可以開展流面的生成工作，基本的思想采用前沿推進(jìn)法，如圖8（a）所示。前沿推進(jìn)法將流面劃分為若干流帶，相鄰兩條流線之間構(gòu)成一條流帶，記每條流線當(dāng)前最新采樣點(diǎn)為前沿，每個(gè)流帶對應(yīng)一個(gè)追蹤器記錄其兩個(gè)前沿。追蹤器通過單鏈表的方式連接，如圖8（b）所示，流帶通過前沿在流場中推進(jìn)的過程來生成。當(dāng)流場中前沿處速度為0，或者前沿位置超出流場空間時(shí)，該流帶結(jié)束生成。

圖8 三角化法基本思想示意圖

流帶通過三角形面片法表示，即在相鄰的兩條流線的采樣點(diǎn)上生成三角形條帶來近似化表現(xiàn)流帶，三角形的形狀需要盡量規(guī)則，角度不宜過大，這可以通過局部貪心法，即每次新增對角線較短的三角形，來實(shí)現(xiàn)。兩條流線的前沿在流場中各生成一個(gè)新的采樣點(diǎn)后，與追蹤器中記錄的前沿點(diǎn)一起就可以構(gòu)成一個(gè)四邊形，選擇較短對角線與兩前沿連線構(gòu)成三角形面片加入流帶中，并更新追蹤器中的前沿信息。其過程可用如下偽代碼表示。

輸入：以單鏈表形式相互連接的追蹤器Tracer，流場數(shù)據(jù)V，三角形面片集合Output

輸出：無

Advance_Font（Tracer）

1.Left_advanced=false；

2.while（1）

3.L0=Tracer->L0，R0=Tracer->R0；

4.在V中通過龍格庫塔法以L0、R0計(jì)算L1、R1

5.L_diag=|L1R0|，R_diag=|L0R1|

6.M_diag=min（L_diag，R_diag）；

7.if（M_diag==L_diag）

8. ad_left=true；

9.else ad_left=false；

10.if（Left_advanced&&ad_left）return；

11.if（ad_left）

12.Output->push_triangle（L0，R0，L1）；

13.Tracer->L0=L0=L1；

14.Left_advanced=true；

15.else

16.Output->push_triangle（L0，R0，R1）；

17.Advance_Font（Tracer->next）；

18.Tracer->R0=R0=R1；

隨著流場中流線的分散與聚攏，流帶的寬度也會產(chǎn)生變化，當(dāng)流線不斷聚攏時(shí)，流帶的寬度會變窄，置之不理的話會生成多余的三角形，影響到繪制的時(shí)間，此時(shí)可以考慮將其與相鄰流帶合并。當(dāng)流線不斷分散時(shí)，流帶的寬度會增大，置之不理的話，會影響到生成曲面的精確度，此時(shí)可以考慮將該流帶分裂。

判斷合并的標(biāo)準(zhǔn)：如圖9（a）所示，流帶Ribbon1當(dāng)前的兩前沿分別為L0、M0。其相鄰流帶Ribbon2的兩前沿為M0、R0。通過L0、M0、R0三點(diǎn)積分得到的采樣點(diǎn)為L1、M1、R1。定義四邊形L0L1R1R0的高H為L1與R1兩點(diǎn)到線段L0R0所屬直線的距離的較小值。當(dāng)點(diǎn)L1與點(diǎn) R1間的距離小于 H，且 L0、L1、M0、M1、R0、R1六點(diǎn)近似共面時(shí)，將流帶Ribbon1與流帶Ribbon2合并。

圖9 流帶合并分裂示意圖

合并流帶的方法：如圖9（a）所示，在流帶Ribbon2中，選擇對角線M0R1，將?M0R0R1加入到流帶中，更新右前沿為 R1，在流帶Ribbon1中，選擇對角線L1M0，將?L0M0L1加入到流帶中，更新追蹤器1左前沿為L1，將?L1M0R1加入到流帶中，將追蹤器1右前沿更新為R1，并將追蹤器1右邊的追蹤器更新為追蹤器2右邊的追蹤器。其過程如下偽代碼所示：

輸入：非隊(duì)尾追蹤器Tracer，流場數(shù)據(jù)V，三角形面片集合Output

輸出：無

Merge_Ribbon（Tracer）

1.L0=Tracer->L0，M0=Tracer->R0；

2.R0=Tracer->next->R0；

3.在V中通過龍格庫塔法以L0、M0、R0計(jì)算L1，M1，R1

4.Output->push_triangle（L0，L1，M0）；

5.Output->push_triangle（M0，R0，R1）；

6.Output->push_triangle（L1，M0，R1）；

7.Tracer->L0=L1，Tracer->R0=R0；

8.Tmp=Tracer->next

9.Tracer->next=Tracer->next->next；

10.Delete tmp；

判斷分裂的標(biāo)準(zhǔn)：當(dāng)新采樣點(diǎn)L1與點(diǎn)R1的距離|L1R1|大于四邊形 L0L1R1R0的高 H的兩倍，即|L1R1|≥2H時(shí)，流帶需要分裂。

分裂流帶的方法：如圖9（b）所示，在流帶中用適當(dāng)?shù)姆椒尤胄曼c(diǎn)M1，將?L0M1R0加入流帶中，生成新流帶Ribbon2，其左前沿為M1，右前沿為Ribbon1的右前沿R0，更新Ribbon1的右前沿為M1。將Ribbon2右邊的流帶設(shè)置為Ribbon1右邊的流帶，Ribbon1右邊的流帶設(shè)置為Ribbon2。其過程如下偽代碼所示：

輸入：追蹤器Tracer，流場數(shù)據(jù)V，三角形面片集合Output

輸出：無

Split_Ribbon（Tracer）

1.L0=Tracer->L0，R0=Tracer->R0

2.通過龍格庫塔法在V中以L0、R0計(jì)算L1、R1

3.計(jì)算新加入點(diǎn)M1的坐標(biāo)

4.Output->push_triangle（L0，M1，R0）；

5.Tracer->R0=M1；

6.New Tracer2；

7.Tracer2->L0=M1，Tracer2->R0=R0；

8.Tracer2->next=Tracer->next；

9.Tracer->next=Tracer2；

加入新點(diǎn)M1的坐標(biāo)在取線段L1R1中點(diǎn)的基礎(chǔ)上配合施密特插值法使得最后生成的流面更平滑。其基本思想如圖10所示。

圖10 施密特插值法示意圖

在插入新點(diǎn)時(shí)考慮到流場沿著種子線方向的變化，決定點(diǎn)M1是位于中點(diǎn)偏上還是偏下，其計(jì)算方法為：

其中，，由于是基于中值法改進(jìn)，通常r取0.5。

流面生成效果如圖11所示，直接三角化法得出的圖像由于無法適應(yīng)流線變化會產(chǎn)生失真，如圖11（a）所示。前沿推進(jìn)法產(chǎn)生的圖像不夠平滑，如圖11（b）所示。配合施密特插值法則能得到平滑又精確的流面，如圖11（c）所示。

圖11 流面效果圖

5 實(shí)驗(yàn)結(jié)果分析

測試平臺為Win 10操作系統(tǒng)，Intel CORE i7，8 GB內(nèi)存，顯卡為NVIDIA GEFORCE GTX 860M。

圖12展現(xiàn)了通過自定義種子線與通過本文算法自動布線在Delta-wing流場數(shù)據(jù)中生成流面的情況：圖12（a）為連接端點(diǎn)（0.334 5，－200，40.404）與（0.334 5，200，40.404）為種子線生成的流面，圖12（b）為連接端點(diǎn)（201.003，－200，90.909）與（201.003，200，90.909）為種子線生成的流面。可以看到它們都屬于比較穩(wěn)定的流面，且值覆蓋了流場中一小部分，難以對流場整體結(jié)構(gòu)進(jìn)行描繪。圖12（c）是用本文方法以 K=2，α=0.5，β=0.3生成的流面，可以看到它畫出了兩塊近似對稱的區(qū)域。圖12（d）是用本文方法以 K=4，α=0.5，β =0.8生成的流面，可以看到在兩塊對稱區(qū)域的基礎(chǔ)上增加了兩塊平滑曲面，圖12（e）是用本文方法以 K=6，α=0.8，β=0.8生成的流面，可以看到相對于圖12（d）中圖像，新增了兩塊平滑曲面，它們找出了流場中變化較劇烈的部分，并生成了有代表性的圖像，并且由于圖12（e）中的圖像相比圖12（d）中的圖像并無劇烈的變化，不再追加聚類個(gè)數(shù)。另外，用Esturo的方法在Telda-wing上所得圖像[15]，只能生成一張流面，容易給人造成流場中只有該區(qū)域有流動的錯(cuò)覺，從而對流場整體的狀況產(chǎn)生誤判。

圖13展示了通過自定義種子線與本文方法在5Jet流場數(shù)據(jù)中所生成流面的情況。圖13（a）為連接端點(diǎn)（1，64，90）與端點(diǎn)（128，64，90）為種子線所生成的流面，圖13（b）為連接端點(diǎn)（1，64，64）與端點(diǎn)（128，64，64）為種子線所生成的流面?？梢钥闯?，這兩張流面都只能反映出流場局部的流動情況，無法對該流場整體的流動情況作出判斷。圖13（c）是本文方法以 K=2，α=0.3，β=0.3生成的流面。圖13（d）是本文方法以K=4，α=0.5，β=0.8生成的流面，保持整體結(jié)構(gòu)的同時(shí)增加了平滑曲面。圖13（e）是以 K=6，α=0.5，β =0.3生成的流面，保持圖13（d）中整體結(jié)構(gòu)的同時(shí)新增了兩張平滑曲面。至此，認(rèn)為5JET數(shù)據(jù)的整體流動結(jié)構(gòu)已得出，不再增加聚類個(gè)數(shù)。

圖12 Delta-Wing實(shí)驗(yàn)結(jié)果圖

圖13 5JET實(shí)驗(yàn)結(jié)果圖

6 結(jié)論

本文設(shè)計(jì)了一個(gè)流場種子線自動布局以及流面自動生成的算法，通過坐標(biāo)、曲率以及梯度參數(shù)對數(shù)據(jù)點(diǎn)進(jìn)行聚類，得到中心點(diǎn)，在曲率場中通過四階龍格庫塔法得到種子線，然后基于前沿推進(jìn)法得到流面，能夠有效地在未知流場數(shù)據(jù)中生成表現(xiàn)力較強(qiáng)的流面，也能為人工布線位置提供較好的參考。目前，該算法還需要用戶自己輸入聚類個(gè)數(shù)K，希望將來能夠開發(fā)出自動確定聚類個(gè)數(shù)的算法。該算法依然是串行執(zhí)行的版本，未來將實(shí)現(xiàn)并行計(jì)算的版本，以提高程序運(yùn)行的速度。此外，該算法所生成的流面依然存在流面間遮擋的問題，將來希望開發(fā)出能夠根據(jù)觀察者視角與流面遮擋情況自動調(diào)節(jié)局部透明度的版本。

[1]宋漢戈，劉世光.三維流場可視化綜述[J].系統(tǒng)仿真學(xué)報(bào)，2016（9）：1929-1936.

[2]Edmunds M，Laramee R S，Chen G，et al.Surface-based flow visualization[J].Computers&Graphics，2012，36（8）：974-990.

[3]Mcloughlin T，Laramee R S，Peikert R，et al.Over two decades of integration-based，geometric flow visualization[J].Computer Graphics Forum，2010，29（6）：1807-1829.

[4]Laramee R S，Garth C，Doleisch H，et al.Visual analysis and exploration of fluid flow in a cooling jacket[C]//Proceedings IEEE Visualization 2005，2005：623-630.

[5]Hultquist J P M.Constructing stream surfaces in steady 3D vector fields[C]//IEEE Conference on Visualization（Visualization’92），1992：171-178.

[6]Garth C，Tricoche X，Salzbrunn T，et al.Surface techniques for vortex visualization[C]//Symposium on Visualization，Konstanz（Vissym 2004），Germany，May 2004：155-164.

[7]Scheuermann G，Bobach T，Hagen H，et al.A tetrahedrabased stream surface algorithm[C]//IEEE Visualization 2001，San Diego，CA，USA，2001：151-553.

[8]Wijk J J V.Implicit stream surfaces[C]//IEEE Conference on Visualization（Visualization’93），1993：245-252.

[9]Schneider D，Wiebel A，Scheuermann G.Smooth stream surfaces of fourth order precision[J].Computer Graphics Forum，2009，28（3）：871-878.

[10]Schulze M，Germer T，R?ssl C，et al.Stream surface parametrization by flow-orthogonal front lines[C]//Computer Graphics Forum，2012：1725-1734.

[11]Banks D C，Singer B.Vortex tubes in turbulent flows：Identification，representation，reconstruction[C]//IEEE Conference on Visualization（Visualization’94），1994：132-139.

[12]Roth M.Automatic extraction of vortex core lines and other line type features for scientific visualization[microform][D].Swiss Federal Institute of Technology Zurich，2000.

[13]Rokach L.A survey of clustering algorithms[M]//Data Mining and Knowledge Discovery Handbook.Boston，MA：Springer，2009：269-298.

[14]Arthur D，Vassilvitskii S.k-means++：The advantages of careful seeding[C]//Eighteenth ACM-SIAM Symposium on Discrete Algorithms（SODA 2007），New Orleans，Louisiana，USA，January 2007：1027-1035.

[15]Esturo J M，Schulze M，R?ssl C，et al.Global selection of stream surfaces[J].Computer Graphics Forum，2013，32（2pt1）：113-122.