基于全線程樹數(shù)據(jù)結(jié)構(gòu)的笛卡爾網(wǎng)格高效生成技術(shù)

陳浩，畢林，華如豪，周清清，唐志共，袁先旭,*

1.中國(guó)空氣動(dòng)力研究與發(fā)展中心 空氣動(dòng)力學(xué)國(guó)家重點(diǎn)實(shí)驗(yàn)室，綿陽(yáng) 621000

2.中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所，綿陽(yáng) 621000

隨著CFD在工程實(shí)際中的深入應(yīng)用，模擬中所面臨的的幾何外形也越來(lái)越復(fù)雜。在CFD計(jì)算中，網(wǎng)格生成需要頻繁的人機(jī)交互，并且生成的網(wǎng)格質(zhì)量嚴(yán)重依賴用戶經(jīng)驗(yàn)。有研究表明，網(wǎng)格生成過(guò)程占據(jù)了CFD周期中總用時(shí)的70%以上，是制約CFD效率的主要瓶頸。如何減少網(wǎng)格生成人工干預(yù)，提高高質(zhì)量網(wǎng)格生成自動(dòng)化水平，成為當(dāng)前網(wǎng)格技術(shù)領(lǐng)域的研究熱點(diǎn)之一。

通常使用的網(wǎng)格大體上分為3類：貼體結(jié)構(gòu)和非結(jié)構(gòu)，以及非貼體的笛卡爾網(wǎng)格。貼體結(jié)構(gòu)網(wǎng)格利于附面層特征的模擬，具有較高的質(zhì)量，但針對(duì)復(fù)雜外形生成要分多塊進(jìn)行，要考慮不同塊之間的過(guò)渡是否光滑。這在拓?fù)浣Y(jié)構(gòu)復(fù)雜時(shí)不易生成，即使是對(duì)網(wǎng)格生成有經(jīng)驗(yàn)的人員，仍然需要大量的時(shí)間。非結(jié)構(gòu)網(wǎng)格與之相比，生成難度降低，自動(dòng)化程度提高，但是黏性特征的捕捉能力卻不足，并且非結(jié)構(gòu)網(wǎng)格還存在存儲(chǔ)量較大、不利于高階高精度格式應(yīng)用的限制。笛卡爾網(wǎng)格方法相對(duì)前兩種方法，可以完全自動(dòng)化生成，而且其自適應(yīng)加密和粗化十分方便，使得網(wǎng)格的質(zhì)量較好，能夠較準(zhǔn)確地捕捉流動(dòng)結(jié)構(gòu)特征。因此，自適應(yīng)笛卡爾網(wǎng)格方法吸引了很多學(xué)者關(guān)注。

但是笛卡爾網(wǎng)格方法也有其不足，它的非貼體特性導(dǎo)致在與壁面相交處會(huì)產(chǎn)生大小、形狀差異比較大的單元，從而為壁面附近的高精度模擬帶來(lái)不便。為了降低這種不足帶來(lái)的負(fù)面影響，尤其是對(duì)于黏性流動(dòng)模擬，一般采用混合網(wǎng)格方法、切割或融合網(wǎng)格方法以及浸入邊界方法等壁面處理方式。其中，前兩者實(shí)際上是把非貼體的壁面網(wǎng)格轉(zhuǎn)換成貼體的形式，提供了一種可行的選擇，但是其壁面處理的難度也大大增加，降低了笛卡爾網(wǎng)格的靈活性。很多學(xué)者采用具備簡(jiǎn)單性和靈活性的浸入邊界方法，并通過(guò)提高壁面網(wǎng)格加密程度或借助于壁面函數(shù)的方式，提高黏性流動(dòng)的模擬精度。

此外，為了方便自適應(yīng)，笛卡爾網(wǎng)格通常采用四叉樹或八叉樹數(shù)據(jù)結(jié)構(gòu)。雖然這種數(shù)據(jù)結(jié)構(gòu)存儲(chǔ)量較少，但很多信息需要通過(guò)臨時(shí)計(jì)算的方式獲取，從而導(dǎo)致計(jì)算量較大。如果采用非結(jié)構(gòu)網(wǎng)格的數(shù)據(jù)結(jié)構(gòu)，雖然方便調(diào)用網(wǎng)格信息，但存儲(chǔ)量會(huì)較大。Khokhlov提出的全線程樹數(shù)據(jù)結(jié)構(gòu)給出了一種較好的選擇，其在叉樹數(shù)據(jù)結(jié)構(gòu)的基礎(chǔ)上充分利用閑置的內(nèi)存節(jié)點(diǎn)，重新組織網(wǎng)格的存儲(chǔ)信息和連接方式，方便相鄰單元的檢索，以減少計(jì)算量。但目前其尚未得到推廣應(yīng)用，需要進(jìn)一步發(fā)展。

通常，在相同的模擬精度時(shí)，笛卡爾網(wǎng)格量會(huì)明顯高于貼體結(jié)構(gòu)網(wǎng)格的數(shù)目。對(duì)于復(fù)雜外形的流動(dòng)模擬，尤其是在非定常情況下，網(wǎng)格量會(huì)進(jìn)一步增加，從而更加凸顯笛卡爾網(wǎng)格面臨的計(jì)算量和存儲(chǔ)量的問(wèn)題，使得網(wǎng)格生成耗時(shí)明顯增加，不可忽略。很多學(xué)者通過(guò)大規(guī)模并行技術(shù)來(lái)應(yīng)對(duì)上述難題，但是高性能計(jì)算平臺(tái)并未普及，研究人員仍然需要在個(gè)人計(jì)算機(jī)上完成CFD數(shù)值模擬工作。

針對(duì)以上問(wèn)題，有必要構(gòu)建一種高效的三維自適應(yīng)笛卡爾網(wǎng)格生成方法，減少大規(guī)模網(wǎng)格量時(shí)的網(wǎng)格生成時(shí)間，以提高笛卡爾網(wǎng)格方法對(duì)復(fù)雜外形和復(fù)雜流動(dòng)問(wèn)題的實(shí)用性。

1 自適應(yīng)笛卡爾網(wǎng)格生成流程

自適應(yīng)笛卡爾網(wǎng)格的加密方式主要分為兩種：各向同性和各向異性，本文討論的是各向同性的方式，即在三維情況下，網(wǎng)格被分割成8個(gè)相同的子單元。為了提高網(wǎng)格的生成效率，本文發(fā)展或應(yīng)用了很多技術(shù)方法，如采用全線程樹的笛卡爾網(wǎng)格數(shù)據(jù)結(jié)構(gòu)、基于KD樹的物面數(shù)據(jù)結(jié)構(gòu)、輔助判斷網(wǎng)格類型的染色算法等，具體的生成流程如下：

讀取基于STL(STereo Lithography)格式的外形數(shù)據(jù)文件，構(gòu)造物面離散三角形網(wǎng)格的KD樹，方便后續(xù)步驟快速檢索。

根據(jù)外形的幾何尺寸設(shè)定計(jì)算域，并設(shè)定網(wǎng)格步長(zhǎng)，劃分均勻的初始離散網(wǎng)格。

通過(guò)FTT數(shù)據(jù)結(jié)構(gòu)建立笛卡爾網(wǎng)格的信息存儲(chǔ)和網(wǎng)格間的連接方式。

判斷笛卡爾網(wǎng)格與物面的相交關(guān)系。

通過(guò)射線相交方法，結(jié)合染色算法，判斷笛卡爾網(wǎng)格相對(duì)物面的位置關(guān)系(內(nèi)部或者外部)。

將與物面相交的笛卡爾網(wǎng)格單元進(jìn)行一定次數(shù)的加密，并將曲率變化較大的位置進(jìn)一步加密，直至網(wǎng)格的密度達(dá)到了能夠準(zhǔn)確描述幾何外形特征的程度。

對(duì)于分割得到的新出現(xiàn)的網(wǎng)格，通過(guò)步驟4和步驟5判斷網(wǎng)格與物面相交與否，以及網(wǎng)格中心點(diǎn)相對(duì)物面的位置。

判斷新出現(xiàn)的網(wǎng)格與相鄰單元的層次差是否小于1，對(duì)不滿足的鄰居單元進(jìn)行加密，確保網(wǎng)格的質(zhì)量。

進(jìn)行一定步數(shù)的計(jì)算，根據(jù)流場(chǎng)解特征進(jìn)一步自適應(yīng)，并重復(fù)步驟8。

2 笛卡爾網(wǎng)格數(shù)據(jù)結(jié)構(gòu)

笛卡爾網(wǎng)格的數(shù)據(jù)結(jié)構(gòu)主要分為3種：結(jié)構(gòu)式、非結(jié)構(gòu)式和叉樹式，為了便于自適應(yīng)操作和節(jié)省內(nèi)存存儲(chǔ)，一般采用叉樹結(jié)構(gòu)居多。叉樹數(shù)據(jù)結(jié)構(gòu)本質(zhì)上是一種分層式數(shù)據(jù)結(jié)構(gòu)，不同層之間的劃分依據(jù)是網(wǎng)格的大小，而決定網(wǎng)格位置的是其所處的上一層單元，也就是父單元，同樣的，該單元也可以作為父單元?jiǎng)澐值玫綌?shù)個(gè)子單元,叉樹不同層級(jí)的父子單元通過(guò)指針建立訪問(wèn)和聯(lián)系的方式。由于本文工作針對(duì)三維外形，所以以八叉樹為例，如圖1所示。

圖1 基于八叉樹數(shù)據(jù)結(jié)構(gòu)的網(wǎng)格分割Fig.1 Cell subdivision using octree data structure

采用叉樹數(shù)據(jù)結(jié)構(gòu)的優(yōu)勢(shì)之一是自適應(yīng)的易操作性，但也有不足之處，就是網(wǎng)格信息查詢和調(diào)用的計(jì)算量偏大。例如，網(wǎng)格單元的相鄰單元沒(méi)有預(yù)先存儲(chǔ)，需要在調(diào)用的時(shí)候，通過(guò)查詢算法獲得。此外，叉樹結(jié)構(gòu)的存儲(chǔ)也存在一定的內(nèi)存浪費(fèi)，因?yàn)槿~子節(jié)點(diǎn)的網(wǎng)格單元，其分配的子單元指針并沒(méi)有指向，實(shí)際上是閑置的，在網(wǎng)格規(guī)模很大時(shí)，這部分浪費(fèi)的內(nèi)存不容忽視。針對(duì)上述不足，Khokhlov提出了全線程樹數(shù)據(jù)結(jié)構(gòu)(Fully Threaded Tree，F(xiàn)TT)。這種數(shù)據(jù)結(jié)構(gòu)是傳統(tǒng)的叉樹數(shù)據(jù)結(jié)構(gòu)的改進(jìn)版本,它將叉樹數(shù)據(jù)結(jié)構(gòu)中葉子節(jié)點(diǎn)閑置的指針利用起來(lái)，指向鄰居單元的父單元，從而構(gòu)建快速訪問(wèn)鄰居單元的“線程”。

FTT在最初提出時(shí)，限于當(dāng)時(shí)的硬件水平，存儲(chǔ)信息較少。在本文工作中，通過(guò)增加關(guān)鍵信息的存儲(chǔ)，對(duì)該數(shù)據(jù)結(jié)構(gòu)進(jìn)行改進(jìn)，以減少網(wǎng)格生成過(guò)程中的計(jì)算量，進(jìn)而提高效率。

本文選擇存儲(chǔ)的信息是與空間網(wǎng)格單元相交的各個(gè)物面網(wǎng)格單元的序號(hào)，其對(duì)網(wǎng)格生成的效率有著重要影響。自適應(yīng)網(wǎng)格生成過(guò)程中一個(gè)計(jì)算量很大的操作是判斷網(wǎng)格是否與物面相交。傳統(tǒng)的方式中，網(wǎng)格類型判斷過(guò)程中父和子單元是獨(dú)立判斷的，使得空間網(wǎng)格單元與各個(gè)物面單元間相交信息需要重復(fù)計(jì)算。事實(shí)上，與子單元相交的物面網(wǎng)格肯定與父單元相交，遵循這個(gè)規(guī)律，如果能把與父單元相交物面進(jìn)行動(dòng)態(tài)存儲(chǔ)，那么子單元再進(jìn)行相交判斷時(shí)候只需要在與父單元相交的網(wǎng)格面這樣一個(gè)較小的范圍內(nèi)進(jìn)行搜索，能夠極大地降低計(jì)算量，從而提高網(wǎng)格單元類型判斷和自適應(yīng)生成的效率。

基于上述思路，利用FTT按照如圖2所示方式構(gòu)建數(shù)據(jù)結(jié)構(gòu)。圖中，父單元-子單元、子單元-父單元以及鄰居單元-鄰居單元之間的指針?lè)謩e用不同顏色的箭頭表示，通過(guò)Oct來(lái)組織，每個(gè)Oct包含8個(gè)單元。每個(gè)單元有一個(gè)與之關(guān)聯(lián)的物理狀態(tài)矢量，以及一個(gè)指向包含其子單元的Oct指針，如果沒(méi)有的話，指向NULL。每個(gè)Oct包含網(wǎng)格單元的層次OctLv、網(wǎng)格中心點(diǎn)坐標(biāo)(,,)、指向父單元的指針OctPr以及指向各個(gè)方向鄰居的父單元的指針數(shù)組OctNb(6)。本文增加了存儲(chǔ)與網(wǎng)格單元相交的各個(gè)物面網(wǎng)格單元的序號(hào)的動(dòng)態(tài)數(shù)組OctTC(∶)。所有單元存在目標(biāo)空間網(wǎng)格單元相交的物面網(wǎng)格單元時(shí)，動(dòng)態(tài)數(shù)組就會(huì)被分配內(nèi)存，并將所有相交物面網(wǎng)格單元的序號(hào)逐個(gè)存儲(chǔ)。在完成子單元的計(jì)算后，父單元的存儲(chǔ)數(shù)組會(huì)被釋放內(nèi)存。

圖2 全線程數(shù)據(jù)結(jié)構(gòu)中單元和Oct的關(guān)系Fig.2 Relationship between cells and Oct in FTT

大體上來(lái)說(shuō)，采用傳統(tǒng)的鄰居單元查找方式平均遍歷的樹層數(shù)約為2，而如果采用FTT結(jié)構(gòu)，平均遍歷層數(shù)則約為1/2，在搜索鄰居的時(shí)間上可以節(jié)省大約75%。在內(nèi)存占用方面，采用傳統(tǒng)的叉樹數(shù)據(jù)結(jié)構(gòu)，要存儲(chǔ)完上述信息，每個(gè)單元需要8 words的內(nèi)存，采用FTT方式，每個(gè)單元需要2 words的內(nèi)存，相對(duì)于傳統(tǒng)方式明顯降低。

由于增加了相交信息的存儲(chǔ)，每個(gè)單元內(nèi)存占用會(huì)有所增加。當(dāng)物面離散網(wǎng)格數(shù)目為，空間網(wǎng)格數(shù)目為，在考慮相鄰單元對(duì)于相交物面的重復(fù)記數(shù)的情況下，由相交信息存儲(chǔ)所引起的平均每個(gè)Oct的內(nèi)存增加量約為2/words，或者每個(gè)單元的內(nèi)存增加量約為/(4) words。排除極特殊的情況，一般空間網(wǎng)格量明顯多于物面網(wǎng)格數(shù)，即/，所以平均每個(gè)單元的內(nèi)存增加不大于1/4 words，整體上的內(nèi)存占用仍然在一個(gè)較低的水平。

3 初始網(wǎng)格生成

在進(jìn)行流動(dòng)計(jì)算前，首先需要根據(jù)幾何外形的尺寸設(shè)定計(jì)算域。根據(jù)數(shù)據(jù)文件里幾何構(gòu)型各個(gè)離散面元的信息，可以確定所有面元各個(gè)頂點(diǎn)坐標(biāo)的最小和最大值，從而建立幾何最小的方形包圍盒[(,,), (,,)]，如圖3所示。計(jì)算域[(,,), (,,)]應(yīng)完全包含整個(gè)幾何外形。

圖3 計(jì)算域和幾何包圍盒Fig.3 Computational domain and geometry bounding box

計(jì)算域的大小設(shè)定除了需要考慮幾何尺寸外，還有來(lái)流的情形和流場(chǎng)的結(jié)構(gòu)特點(diǎn)，這就得根據(jù)具體情況設(shè)定。需要說(shuō)明的是，本文工作是針對(duì)封閉的幾何外形計(jì)算外部流動(dòng)的問(wèn)題，其他情況(如圓管流動(dòng))的網(wǎng)格生成方法與本文的工作原理類似。

在完成上述工作之后，就可以劃分初始的計(jì)算網(wǎng)格。首先需要確定大小合適的初始網(wǎng)格尺寸，既需要保證基本的流動(dòng)分辨率，又不能太小導(dǎo)致初始網(wǎng)格量太大。這里對(duì)計(jì)算域均勻離散：

(,,)=(,,)

(1)

式中：、、為3個(gè)坐標(biāo)軸方向的離散單元尺寸，為設(shè)定的初始網(wǎng)格邊長(zhǎng)。

沿笛卡爾坐標(biāo)軸各個(gè)方向的網(wǎng)格數(shù)目分別為

(2)

式中：、、分別是沿各個(gè)坐標(biāo)軸方向計(jì)算域的長(zhǎng)度。

4 網(wǎng)格類型判斷

對(duì)于本文的笛卡爾網(wǎng)格生成技術(shù)而言，很重要的一項(xiàng)工作是判斷網(wǎng)格類型。按照網(wǎng)格中心點(diǎn)對(duì)于封閉的物面的位置(內(nèi)部或外部)，分為物體外部網(wǎng)格點(diǎn)和物體內(nèi)部網(wǎng)格點(diǎn)。按照相交關(guān)系，又分為物面相交單元、物面內(nèi)部單元和物面外部單元。這是后續(xù)進(jìn)行網(wǎng)格自適應(yīng)生成和浸入邊界處理的基礎(chǔ)。

4.1 判斷網(wǎng)格點(diǎn)類型

4.1.1 射線相交方法

在建立幾何包圍盒后，不難發(fā)現(xiàn)，網(wǎng)格中心點(diǎn)在包圍盒外時(shí)，肯定是物面外部的點(diǎn)。該算法的計(jì)算量很小，可以幫助快速確定包圍盒外部的網(wǎng)格點(diǎn)的類型，但是包圍盒內(nèi)部的網(wǎng)格點(diǎn)類型仍然需要通用的判斷方法。很多學(xué)者通過(guò)簡(jiǎn)單而且高效的射線相交方法進(jìn)行上述工作。通常情況下，在要判斷的目標(biāo)點(diǎn)選定任意方向的射線之后，就可以通過(guò)計(jì)數(shù)射線與封閉的物面相交的點(diǎn)數(shù)，來(lái)判斷網(wǎng)格點(diǎn)相對(duì)于物面的位置，即相交數(shù)目為奇數(shù)時(shí)，點(diǎn)在內(nèi)部，偶數(shù)時(shí)則在外部。然而在使用射線相交方法時(shí)，可能存在射線恰好經(jīng)過(guò)外形拐點(diǎn)的情況，導(dǎo)致判斷的結(jié)果出現(xiàn)錯(cuò)誤，如圖4所示。圖中，和代表目標(biāo)點(diǎn)，和分別代表經(jīng)過(guò)目標(biāo)點(diǎn)的射線，代表幾何外形，?代表幾何邊界，1和2即為經(jīng)過(guò)拐點(diǎn)的射線。

圖4 基于射線相交方法判斷點(diǎn)在多邊形內(nèi)外Fig.4 Illustration of point inside or outside polygon using ray-casting method

針對(duì)上述問(wèn)題，提出了一種新的魯棒的判斷算法，與射線相交方法相配合，可以準(zhǔn)確識(shí)別射線經(jīng)過(guò)外形拐點(diǎn)的不同情形。這種方法是通過(guò)判斷這些線段或面是在射線的同側(cè)還是異側(cè)，若是前者的話，在射線與物面相交點(diǎn)計(jì)數(shù)時(shí)需要在拐點(diǎn)處計(jì)數(shù)增加1，若是后者的話，則正常計(jì)數(shù)，本文將這種算法稱為奇異性判斷算法。

首先以二維情況為例，通過(guò)圖5和圖6說(shuō)明。在圖5中，射線以待判斷的目標(biāo)點(diǎn)為起始點(diǎn)，經(jīng)過(guò)離散線段之間的拐點(diǎn)，兩條以為端點(diǎn)的線段分別為和。當(dāng)滿足：

(3)

圖5 射線經(jīng)過(guò)拐點(diǎn)示意圖：A和B在射線同側(cè)Fig.5 Illustration of ray passing inflection point: A and B on same side of ray

圖6 射線經(jīng)過(guò)拐點(diǎn)示意圖：A和B在射線異側(cè)Fig.6 Illustration of ray passing inflection point: A and B on different sides of ray

和均在射線的同一側(cè)。其中為以為起點(diǎn)，方向與相同的矢量，則為與之方向相反的矢量。如圖5所示。射線與線段對(duì)應(yīng)的向量以及射線的反向矢量與線段對(duì)應(yīng)的向量的夾角分別為和，顯然，兩條線段分別對(duì)應(yīng)的外法向?yàn)?span id="j5i0abt0b" class="emphasis_italic">和。若=0°或=0°，則射線與其中一條線段重合一定的長(zhǎng)度，本文采用重新選擇射線方向的處理方式，避免上述重合情況出現(xiàn)。

在圖6中，拐點(diǎn)處的線段和在射線的不同側(cè)。這種情形應(yīng)滿足的公式為

(4)

若〈,〉或〈,〉=0°,π,則同樣的,射線與其中一條線段重合一定的長(zhǎng)度,需要通過(guò)重新選擇射線方向的方式進(jìn)一步判斷和處理。

在三維外形的情形下，可以通過(guò)選取特定切割平面的方式，通過(guò)平面與離散三角形的相交線，轉(zhuǎn)換成二維情況進(jìn)行計(jì)算。

在圖7中，射線起點(diǎn)為，經(jīng)過(guò)位于3個(gè)離散面上的拐點(diǎn)。此時(shí)選取不經(jīng)過(guò)點(diǎn)的任意一條邊上任意一點(diǎn)(除了端點(diǎn))，本文選取的是點(diǎn)，將其與射線上的任意一點(diǎn)(除了起點(diǎn))以及起點(diǎn)組成切割平面，該平面會(huì)切割這些離散三角形得到切割線和。在該切割平面內(nèi)，就可以按照前面二維的方式來(lái)判斷射線與離散面之間的位置關(guān)系。

圖7 通過(guò)切割面將三維問(wèn)題轉(zhuǎn)化為二維Fig.7 Conversion of 3D problem to 2D by cutting plane passing ray

4.1.2 物面網(wǎng)格快速查找方法

在采用射線相交方法進(jìn)行判斷時(shí)，需要計(jì)算多個(gè)點(diǎn)乘和叉乘，當(dāng)空間網(wǎng)格和物面網(wǎng)格量比較大時(shí)，整體計(jì)算量還是很大的。實(shí)際上網(wǎng)格類型的判斷占據(jù)了整個(gè)網(wǎng)格生成過(guò)程的大部分時(shí)間。雖然前面采用幾何包圍盒進(jìn)行了過(guò)濾，可以縮短一定的判斷時(shí)間，但十分有限，有必要發(fā)展高效的網(wǎng)格類型判斷方法。為此，一些學(xué)者采用ADT(Alternating Digital Tree)數(shù)據(jù)結(jié)構(gòu)來(lái)存儲(chǔ)物面網(wǎng)格，以縮小可能相交的物面單元的檢測(cè)范圍，從而提高射線相交判斷的效率。

設(shè)空間網(wǎng)格數(shù)目為，物面離散單元(線段或三角形)數(shù)目為，對(duì)于三維構(gòu)型，是(10)～(10)，而一般是(10)～(10)。在不采用任何優(yōu)化算法的情況下，即射線相交判斷時(shí)是逐個(gè)物面單元進(jìn)行計(jì)算，此時(shí)的時(shí)間復(fù)雜度為(·)。在采用ADT結(jié)構(gòu)后，可以將相交計(jì)算的時(shí)間復(fù)雜度降低到(·lg)。然而傳統(tǒng)的ADT所構(gòu)建的二叉樹的平衡性往往受物面單元的分布密度影響，導(dǎo)致查詢效率不穩(wěn)定。為此，采用其改進(jìn)版本—KDT結(jié)構(gòu)，其更容易構(gòu)造平衡優(yōu)質(zhì)的二叉樹。

KDT是一種多維的二叉樹結(jié)構(gòu)，樹中的每個(gè)節(jié)點(diǎn)都是一個(gè)-維節(jié)點(diǎn)，且每一層都和某個(gè)維度相關(guān)。其在構(gòu)造時(shí)，需要將當(dāng)前層的相關(guān)點(diǎn)進(jìn)行排序，然后取中值點(diǎn)，接著確定經(jīng)過(guò)該點(diǎn)與相關(guān)維度垂直的平面，以剖分空間搜索區(qū)域。由于是以每次排序后得到的中值點(diǎn)作為剖分位置，而不是ADT中的在區(qū)域等分點(diǎn)進(jìn)行劃分，因此可以保證多層叉樹的平衡性。

下面以一些物面網(wǎng)格點(diǎn)為例，比較兩種數(shù)據(jù)結(jié)構(gòu)的操作復(fù)雜度。如圖8所示，這些點(diǎn)在數(shù)組中存儲(chǔ)為：/* A, B, C, D, E, F, G, H, I, J, K, L, M, N, O, P, Q */。如果遍歷操作的話，需要1+2+3+4+5+6+7+8+9+10+11+12+13+14+15+16+17=153次。

圖8 隨機(jī)排布的網(wǎng)格點(diǎn)Fig.8 Randomly distributed grids

對(duì)于ADT叉樹結(jié)構(gòu)而言，如圖9所示，遍歷操作需要1+2+2+3+3+3+4+4+4+4+5+5+5+6+6+7+8=72次。

圖9 ADT結(jié)構(gòu)示意圖Fig.9 Illustration of ADT

而對(duì)于KDT叉樹結(jié)構(gòu)而言，如圖10所示，遍歷操作需要1+2+2+3+3+3+3+4+4+4+4+4+4+4+4+5+5=59次，相對(duì)ADT構(gòu)建了更為平衡的叉樹結(jié)構(gòu)，操作的復(fù)雜度明顯降低。

圖10 KDT結(jié)構(gòu)示意圖Fig.10 Illustration of KDT

對(duì)于物面離散三角形網(wǎng)格構(gòu)建KD樹的主要流程如下：

建立三角形物面網(wǎng)格的空間包圍盒，包圍盒與坐標(biāo)軸方向一致，對(duì)應(yīng)著三角形的最大分布空間，可以表示為(,)，其中：

(5)

為三角形網(wǎng)格單元的序號(hào);為空間維度。

構(gòu)建叉樹節(jié)點(diǎn)與空間域的對(duì)應(yīng)關(guān)系。取三角形的中心點(diǎn)c(=1,2,…,)排序后的中值點(diǎn)作為剖分面經(jīng)過(guò)的點(diǎn)，比較各軸向方向上中心點(diǎn)坐標(biāo)的方差，最大者對(duì)應(yīng)的方向?yàn)槠拭娴姆ㄏ颍缮厦嫘畔⒓纯梢源_定區(qū)域的劃分面。

節(jié)點(diǎn)的左、右節(jié)點(diǎn)對(duì)應(yīng)的空間區(qū)域分別記作(,)、(,)，具體為

圖11 定義三角面元的坐標(biāo)范圍Fig.11 Definition of coordinate limits for triangular elements

(6)

(7)

從根節(jié)點(diǎn)開始，通過(guò)循環(huán)劃分，可以將二叉樹節(jié)點(diǎn)與劃分的空間區(qū)域逐一對(duì)應(yīng)。

按照步驟2中的對(duì)應(yīng)關(guān)系，將物面包圍盒逐個(gè)插入節(jié)點(diǎn)。

(8)

則插入父親節(jié)點(diǎn)的左側(cè)，成為左子節(jié)點(diǎn)。如果滿足：

(9)

則插入父親節(jié)點(diǎn)的右側(cè)，成為右子節(jié)點(diǎn)。

完成上述工作之后，進(jìn)行射線相交判斷時(shí)就可以利用構(gòu)建的KD樹，快速查找相交的物面網(wǎng)格，主要步驟如下：

將射線與根節(jié)點(diǎn)對(duì)應(yīng)的包圍盒進(jìn)行相交判斷，滿足相交條件則進(jìn)行步驟2，否則確定目標(biāo)點(diǎn)為幾何外部的點(diǎn)。

基于KD樹進(jìn)行篩選，逐層縮小檢索范圍，直至到葉子結(jié)點(diǎn)那一層。

判斷射線與相交的葉子節(jié)點(diǎn)對(duì)應(yīng)的包圍盒內(nèi)的三角形相交情況。

4.2 網(wǎng)格與物面相交判斷方法

前面的工作是基于射線相交方法對(duì)網(wǎng)格中心點(diǎn)的類型(內(nèi)部點(diǎn)或外部點(diǎn))進(jìn)行判斷，這不足以開展網(wǎng)格自適應(yīng)加密，還需要確定空間網(wǎng)格是否與物面相交。

最原始的方法是對(duì)各個(gè)網(wǎng)格頂點(diǎn)分別判斷網(wǎng)格點(diǎn)的類型，如果同時(shí)在內(nèi)部(或者外部)，則與物面不相交，否則與物面相交。上述方法雖然簡(jiǎn)單，但是計(jì)算量較大。一些學(xué)者通過(guò)ADT或KD方法直接確定網(wǎng)格與物面的相交情況，可以極大提高網(wǎng)格單元相交判斷的效率。其思路是將笛卡爾網(wǎng)格對(duì)應(yīng)的空間區(qū)域作為目標(biāo)區(qū)域，在叉樹中快速檢索與之相交的物面三角形包圍盒，進(jìn)而判斷笛卡爾網(wǎng)格與包圍盒內(nèi)的三角形相交與否。這與前面基于KD樹的射線相交的快速判斷類似。

對(duì)于網(wǎng)格與三角形的相交判斷，可以歸為兩類：

1) 三角形整個(gè)包含于笛卡爾網(wǎng)格內(nèi)，這種情況的判斷只需比較三角形的包圍盒與笛卡爾網(wǎng)格的空間關(guān)系，計(jì)算量很小，可以很快獲得。用(LC, RC)表示笛卡爾網(wǎng)格對(duì)應(yīng)的空間域，當(dāng)滿足以下關(guān)系：

LC≤,LC≤

=1,2,…,

(10)

三角形在笛卡爾網(wǎng)格內(nèi)部,其中為待判斷的笛卡爾網(wǎng)格總數(shù)。

2) 三角形與笛卡爾網(wǎng)格面相交，如圖12所示，可以分解成線段與笛卡爾網(wǎng)格面的相交判斷。

圖12 笛卡爾網(wǎng)格與三角面元相交示意圖Fig.12 Illustration of intersection of Cartesian grids and triangles

4.3 染色算法輔助網(wǎng)格類型判斷

在確定笛卡爾網(wǎng)格與物面相交與否后，需要判斷非相交網(wǎng)格單元是在物面內(nèi)部還是外部。若按照4.1節(jié)中的方式，則通過(guò)射線相交方法對(duì)網(wǎng)格中心點(diǎn)位置進(jìn)行判斷，進(jìn)而獲得網(wǎng)格是在內(nèi)部還是外部的信息。對(duì)于封閉的幾何外形，可以通過(guò)染色算法進(jìn)一步提高與物面不相交的網(wǎng)格單元類型判斷的效率。這種算法是通過(guò)查詢目標(biāo)網(wǎng)格的相鄰單元(非物面相交單元)的類型，直接給予其相同的網(wǎng)格類型，看起來(lái)就像被鄰居單元 “染色”了，故稱為染色算法。在本文采用的FTT數(shù)據(jù)結(jié)構(gòu)下，染色算法具有很好的適用性，因?yàn)猷従訂卧檎液鼙憬?。這相對(duì)于傳統(tǒng)的射線相交方法，計(jì)算量會(huì)大大減少。

如圖13所示，在使用染色方法時(shí)，首先要確定染色源，可以取入口邊界的第1個(gè)網(wǎng)格單元，其必然為物體外部網(wǎng)格，然后與其相鄰的單元只要不與物面相交，則同樣為外部網(wǎng)格。這種判斷方式是面向與物面不相交的網(wǎng)格單元，可以在很快的時(shí)間內(nèi)確定這些單元是在物面內(nèi)部還是外部，見表1。

圖13 染色算法Fig.13 Painting algorithm

表1 染色算法與傳統(tǒng)方法確定網(wǎng)格單元類型對(duì)比Table 1 Performance of painting algorithm vs traditional means for determining cell types

5 幾何自適應(yīng)加密

很多學(xué)者都探索和建立了可靠的笛卡爾網(wǎng)格幾何自適應(yīng)加密方法。其基本思路是：

1) 將物面相交笛卡爾網(wǎng)格和過(guò)渡層網(wǎng)格初步加密特定的次數(shù)，具體次數(shù)需要根據(jù)初始網(wǎng)格和幾何外形的尺寸確定。

2) 將物面曲率變化較大的位置進(jìn)一步加密，主要是通過(guò)比較與目標(biāo)網(wǎng)格相交的相鄰物面網(wǎng)格單元法向矢量變化角度，當(dāng)大于一定閾值時(shí)即進(jìn)行加密。

本文采用上述方式開展幾何自適應(yīng)加密。在加密過(guò)程中，需要注意的是，有些區(qū)域的網(wǎng)格可能并不合理或質(zhì)量不高，影響網(wǎng)格對(duì)應(yīng)叉樹結(jié)構(gòu)的平衡性以及流場(chǎng)計(jì)算的精度和穩(wěn)定性。例如，在加密一定次數(shù)后，相鄰網(wǎng)格單元在數(shù)據(jù)結(jié)構(gòu)中的層次之差大于1，如圖14所示。

圖14 相鄰單元層次差大于1Fig.14 Level difference between adjacent cells larger than 1

針對(duì)上述情況，在網(wǎng)格加密過(guò)程中需要建立一個(gè)檢測(cè)算法，識(shí)別上述情形，并對(duì)質(zhì)量差的網(wǎng)格進(jìn)行加密。本文通過(guò)查詢目標(biāo)單元在各個(gè)方向上的鄰居單元，比較各個(gè)鄰居單元與目標(biāo)單元的層次差，當(dāng)大于1時(shí)即對(duì)該鄰居單元加密，直至層次差等于1。

6 流場(chǎng)解自適應(yīng)

在完成幾何自適應(yīng)之后，就得到了初始計(jì)算網(wǎng)格，但其分辨率對(duì)于流動(dòng)結(jié)構(gòu)的精細(xì)捕捉還是不足的，需要基于流場(chǎng)解進(jìn)行自適應(yīng)加密，使網(wǎng)格的分布更合理、更有利于刻畫流動(dòng)特征結(jié)構(gòu)。

本文對(duì)流動(dòng)的數(shù)值求解基于NND格式對(duì)Navier-Stokes方程進(jìn)行空間離散，時(shí)間推進(jìn)則采用Runge-Kutta方法。為了保證對(duì)激波、剪切層、旋渦等多種流動(dòng)結(jié)構(gòu)的捕捉能力，采用一種綜合判據(jù)——速度的散度和旋度相結(jié)合方式。該判據(jù)既可以通過(guò)速度的散度捕捉激波，又可以基于速度的旋度捕捉渦旋和剪切層等特征，相對(duì)于單一判據(jù)具有明顯的優(yōu)勢(shì)。該判據(jù)可表述為

(11)

(12)

式中：、、、為經(jīng)驗(yàn)系數(shù)，在不同的流動(dòng)問(wèn)題中并不相同。由于各類流動(dòng)問(wèn)題主導(dǎo)特征不同，上述各量在決定網(wǎng)格加密或粗化時(shí)起到的作用不同，因而系數(shù)設(shè)置并不統(tǒng)一。后續(xù)研究工作中，將對(duì)笛卡爾網(wǎng)格流場(chǎng)自適應(yīng)判據(jù)系數(shù)的自動(dòng)化配置進(jìn)行研究，進(jìn)一步提高笛卡爾網(wǎng)格方法的自動(dòng)化水平。

7 結(jié)果和討論

7.1 圓球算例

表2顯示了在CPU為intel core i7-8700，3.2 GHz 臺(tái)式機(jī)上的圓球算例測(cè)試結(jié)果，網(wǎng)格如圖15所示。采用相同的初始網(wǎng)格，固定的幾何加密次數(shù)為7次，物面網(wǎng)格數(shù)目是變化的。雖然物面網(wǎng)格數(shù)目不同，但每個(gè)單元的CPU時(shí)間很接近。這得益于對(duì)幾何三角形面元構(gòu)建的KD樹結(jié)構(gòu)，使得物面網(wǎng)格密度對(duì)于每個(gè)單元平均的計(jì)算耗時(shí)影響不大。

表2 圓球算例網(wǎng)格單元生成類型對(duì)比Table 2 Mesh generation after refinement around sphere surface with variable facets

圖15 幾何自適應(yīng)網(wǎng)格Fig.15 Geometric adaptive refinement mesh

對(duì)本文流場(chǎng)自適應(yīng)技術(shù)進(jìn)行測(cè)試，自適應(yīng)網(wǎng)格和云圖如圖16所示。其中來(lái)流馬赫數(shù)=1.5，球半徑為=0.5。網(wǎng)格幾何自適應(yīng)加密次數(shù)AMR=5，流場(chǎng)自適應(yīng)加密次數(shù)為4。經(jīng)過(guò)幾何自適應(yīng)后笛卡爾網(wǎng)格數(shù)目為196 186，經(jīng)過(guò)流場(chǎng)解自適應(yīng)后網(wǎng)格數(shù)目為567 324。

圖16 流場(chǎng)自適應(yīng)網(wǎng)格和云圖Fig.16 Adaptive mesh and contour of flow features

激波距離球體駐點(diǎn)的距離計(jì)算公式為

(13)

在來(lái)流馬赫數(shù)為1.5時(shí)，上述理論計(jì)算值Δ=0302，基于本文方法得到的計(jì)算值為 Δ=0310，與理論值相差2.65%，具有較高的吻合度，證明本文數(shù)值方法可以對(duì)激波位置進(jìn)行較好的預(yù)測(cè)。

在平面上，激波形狀公式為

=

(14)

是鈍體最高點(diǎn)曲率半徑，其通過(guò)與來(lái)流馬赫數(shù)的關(guān)系確定：

(15)

通過(guò)式(13)～式(15)可以確定激波形狀曲線，如圖17所示，與本文數(shù)值模擬云圖結(jié)果對(duì)比可以看出，本文結(jié)果和理論值能夠較好的吻合，對(duì)激波位置能夠較好的捕捉，證明了本文流場(chǎng)自適應(yīng)方法和數(shù)值格式的可靠性。注意到，本文結(jié)果相對(duì)理論值仍有細(xì)微的偏離，主要是因?yàn)榱鲌?chǎng)自適應(yīng)過(guò)程中，細(xì)化網(wǎng)格的流場(chǎng)值需要通過(guò)粗網(wǎng)格的值插值得到，粗網(wǎng)格分辨率不足和插值計(jì)算都會(huì)帶來(lái)一定的誤差。

圖17 激波位置對(duì)比Fig.17 Comparison of shock wave positions

7.2 復(fù)雜構(gòu)型算例

對(duì)三維復(fù)雜構(gòu)型的算例進(jìn)行測(cè)試和評(píng)估，如圖18～圖20所示。表3為網(wǎng)格生成算例,物面網(wǎng)格數(shù)目從翼身組合體構(gòu)型的43 128到導(dǎo)彈構(gòu)型的51 094。表3的第2列和第3列分別為物面三角形數(shù)目和加密后的笛卡爾網(wǎng)格數(shù)目，第3和第4列分別為在串行時(shí)總的CPU時(shí)間和每個(gè)單元在生成過(guò)程中花費(fèi)的CPU時(shí)間。

表3 網(wǎng)格生成復(fù)雜構(gòu)型算例Table 3 Mesh generation tests of complex configurations

圖18 導(dǎo)彈構(gòu)型自適應(yīng)網(wǎng)格Fig.18 Adaptive mesh of missile

圖19 機(jī)翼-副翼構(gòu)型自適應(yīng)網(wǎng)格Fig.19 Adaptive mesh of wing-aileron configuration

圖20 翼身組合體構(gòu)型自適應(yīng)網(wǎng)格Fig.20 Adaptive mesh of wing-body configuration

綜合各個(gè)算例的數(shù)據(jù)可以看出，本文構(gòu)建的三維自適應(yīng)笛卡爾網(wǎng)格生成技術(shù)效率很高，這主要?dú)w功于以下幾個(gè)方面：第一，采用的FFT數(shù)據(jù)結(jié)構(gòu)使得鄰居查詢遍歷層數(shù)少，更快捷；第二，通過(guò)構(gòu)建的物面網(wǎng)格單元的KD樹結(jié)構(gòu)，使得確定網(wǎng)格單元的位置所需要查找的物面網(wǎng)格的范圍極大地縮??；第三，染色算法可以幫助快速確定與物面未相交的單元是在物面內(nèi)部還是外部。值得注意的是，射線相交和網(wǎng)格與物面相交占據(jù)了主要的網(wǎng)格生成CPU時(shí)間。KDT結(jié)構(gòu)的應(yīng)用，可以將潛在的與射線或網(wǎng)格相交的面單元縮小到一個(gè)很小的數(shù)目。采用該方式，可以查詢的復(fù)雜度是(·lg)而不是(·)。因此，每個(gè)單元生成需要的平均CPU時(shí)間受物面網(wǎng)格分辨率和空間網(wǎng)格數(shù)目的影響不大。

8 結(jié) 論

本文發(fā)展了一個(gè)面向三維構(gòu)型的自適應(yīng)笛卡爾網(wǎng)格生成技術(shù)，具備通用性和高效性。為了提高網(wǎng)格生成的效率，采用了幾個(gè)有效的方法。首先，使用FTT數(shù)據(jù)結(jié)構(gòu)作為笛卡爾網(wǎng)格數(shù)據(jù)結(jié)構(gòu)，而不是傳統(tǒng)的八叉樹數(shù)據(jù)結(jié)構(gòu)，使得鄰居查詢更快捷，同時(shí)內(nèi)存利用率更高。對(duì)于物面網(wǎng)格的組織方式則是采用平衡性相對(duì)傳統(tǒng)ADT更好的KDT，盡可能的縮小物面單元查詢檢索的范圍，提高確定網(wǎng)格或射線與物面相交與否的計(jì)算效率。網(wǎng)格點(diǎn)在物面內(nèi)外的確定主要是采用射線相交方法，同時(shí)結(jié)合染色方法幫助快速確定與物面不相交的單元在物面內(nèi)部或外部。為了保證射線相交方法的魯棒性，本文構(gòu)建了識(shí)別經(jīng)過(guò)拐點(diǎn)的射線的判斷算法，幫助準(zhǔn)確的相交計(jì)數(shù)。此外，對(duì)于基于速度散度和旋度點(diǎn)的三維綜合判據(jù)進(jìn)行了發(fā)展，經(jīng)算例測(cè)試，顯示了較高的可靠性。

通過(guò)幾個(gè)三維的幾何構(gòu)型測(cè)試了本文方法，并且得到了魯棒高效的結(jié)果。在接下來(lái)工作中，可以通過(guò)在一些步驟引入基于GPU的并行加速進(jìn)程，進(jìn)一步提高計(jì)算效率，縮短網(wǎng)格生成時(shí)間。