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

1.中國空氣動力研究與發(fā)展中心 空氣動力學國家重點實驗室，綿陽 621000

2.中國空氣動力研究與發(fā)展中心 計算空氣動力研究所，綿陽 621000

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

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

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

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

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

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

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

自適應笛卡爾網(wǎng)格的加密方式主要分為兩種：各向同性和各向異性，本文討論的是各向同性的方式，即在三維情況下，網(wǎng)格被分割成8個相同的子單元。為了提高網(wǎng)格的生成效率，本文發(fā)展或應用了很多技術(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ù)外形的幾何尺寸設定計算域，并設定網(wǎng)格步長，劃分均勻的初始離散網(wǎng)格。

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

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

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

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

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

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

進行一定步數(shù)的計算，根據(jù)流場解特征進一步自適應，并重復步驟8。

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

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

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

采用叉樹數(shù)據(jù)結(jié)構(gòu)的優(yōu)勢之一是自適應的易操作性，但也有不足之處，就是網(wǎng)格信息查詢和調(diào)用的計算量偏大。例如，網(wǎng)格單元的相鄰單元沒有預先存儲，需要在調(diào)用的時候，通過查詢算法獲得。此外，叉樹結(jié)構(gòu)的存儲也存在一定的內(nèi)存浪費，因為葉子節(jié)點的網(wǎng)格單元，其分配的子單元指針并沒有指向，實際上是閑置的，在網(wǎng)格規(guī)模很大時，這部分浪費的內(nèi)存不容忽視。針對上述不足，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)的改進版本,它將叉樹數(shù)據(jù)結(jié)構(gòu)中葉子節(jié)點閑置的指針利用起來，指向鄰居單元的父單元，從而構(gòu)建快速訪問鄰居單元的“線程”。

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

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

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

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

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

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

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

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

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

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

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

(,,)=(,,)

(1)

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

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

(2)

式中：、、分別是沿各個坐標軸方向計算域的長度。

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

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

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

4.1.1 射線相交方法

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

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

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

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

(3)

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

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

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

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

(4)

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

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

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

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

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

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

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

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

下面以一些物面網(wǎng)格點為例，比較兩種數(shù)據(jù)結(jié)構(gòu)的操作復雜度。如圖8所示，這些點在數(shù)組中存儲為：/* 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 隨機排布的網(wǎng)格點Fig.8 Randomly distributed grids

對于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

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

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

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

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

(5)

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

構(gòu)建叉樹節(jié)點與空間域的對應關(guān)系。取三角形的中心點c(=1,2,…,)排序后的中值點作為剖分面經(jīng)過的點，比較各軸向方向上中心點坐標的方差，最大者對應的方向為剖面的法向，由上面信息即可以確定區(qū)域的劃分面。

節(jié)點的左、右節(jié)點對應的空間區(qū)域分別記作(,)、(,)，具體為

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

(6)

(7)

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

按照步驟2中的對應關(guān)系，將物面包圍盒逐個插入節(jié)點。

(8)

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

(9)

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

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

將射線與根節(jié)點對應的包圍盒進行相交判斷，滿足相交條件則進行步驟2，否則確定目標點為幾何外部的點。

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

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

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

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

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

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

1) 三角形整個包含于笛卡爾網(wǎng)格內(nèi)，這種情況的判斷只需比較三角形的包圍盒與笛卡爾網(wǎng)格的空間關(guān)系，計算量很小，可以很快獲得。用(LC, RC)表示笛卡爾網(wǎ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é)中的方式，則通過射線相交方法對網(wǎng)格中心點位置進行判斷，進而獲得網(wǎng)格是在內(nèi)部還是外部的信息。對于封閉的幾何外形，可以通過染色算法進一步提高與物面不相交的網(wǎng)格單元類型判斷的效率。這種算法是通過查詢目標網(wǎng)格的相鄰單元(非物面相交單元)的類型，直接給予其相同的網(wǎng)格類型，看起來就像被鄰居單元 “染色”了，故稱為染色算法。在本文采用的FTT數(shù)據(jù)結(jié)構(gòu)下，染色算法具有很好的適用性，因為鄰居單元查找很便捷。這相對于傳統(tǒng)的射線相交方法，計算量會大大減少。

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

圖13 染色算法Fig.13 Painting algorithm

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

5 幾何自適應加密

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

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

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

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

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

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

6 流場解自適應

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

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

(11)

(12)

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

7 結(jié)果和討論

7.1 圓球算例

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

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

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

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

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

激波距離球體駐點的距離計算公式為

(13)

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

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

=

(14)

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

(15)

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

圖17 激波位置對比Fig.17 Comparison of shock wave positions

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

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

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

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

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

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

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

8 結(jié) 論

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

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