分層嵌套重疊網(wǎng)格自適應(yīng)樹結(jié)構(gòu)動(dòng)態(tài)組裝方法

李曉東， 屈崑， 蔡晉生

西北工業(yè)大學(xué) 航空學(xué)院， 西安 710072

分層嵌套重疊網(wǎng)格自適應(yīng)樹結(jié)構(gòu)動(dòng)態(tài)組裝方法

李曉東， 屈崑*， 蔡晉生

西北工業(yè)大學(xué) 航空學(xué)院， 西安 710072

采用重疊網(wǎng)格可以有效地進(jìn)行復(fù)雜流動(dòng)的大規(guī)模數(shù)值模擬，特別是包含運(yùn)動(dòng)部件(如旋翼、投彈)的動(dòng)態(tài)模擬。本文將樹結(jié)構(gòu)的自適應(yīng)直角網(wǎng)格用于重疊網(wǎng)格組裝過程中的切割和貢獻(xiàn)單元的搜索，大大加快重疊網(wǎng)格的組裝速度。通過二叉樹自適應(yīng)直角網(wǎng)格對(duì)物體外形進(jìn)行離散，實(shí)現(xiàn)切割過程的快速定位；采用八叉樹自適應(yīng)直角網(wǎng)格對(duì)流場(chǎng)區(qū)域進(jìn)行離散，高效地搜索貢獻(xiàn)單元。使用基于壁面距離準(zhǔn)則的重疊區(qū)域最小化方法和分層嵌套重疊策略，能提高重疊網(wǎng)格組裝的效率和質(zhì)量。對(duì)于具有運(yùn)動(dòng)部件的動(dòng)態(tài)重疊網(wǎng)格問題，采用多個(gè)二/八叉樹減少組裝過程中信息更新的冗余計(jì)算，從而大幅度減少重疊網(wǎng)格組裝的時(shí)間消耗。實(shí)際算例的重疊網(wǎng)格組裝結(jié)果說明本文發(fā)展的重疊網(wǎng)格組裝方法具有很高的計(jì)算效率,可以滿足運(yùn)動(dòng)邊界復(fù)雜流動(dòng)問題的動(dòng)態(tài)計(jì)算要求。

嵌套重疊網(wǎng)格； 二叉樹； 八叉樹； 自適應(yīng)直角網(wǎng)格； 動(dòng)態(tài)網(wǎng)格； 非定常流動(dòng)； 計(jì)算流體力學(xué)

隨著復(fù)雜流動(dòng)問題研究的進(jìn)一步發(fā)展，數(shù)值模擬對(duì)網(wǎng)格生成效率和質(zhì)量提出了更高的要求，高質(zhì)量網(wǎng)格是影響CFD計(jì)算的關(guān)鍵因素[1]。目前，許多圖形界面化的商業(yè)軟件采用多塊對(duì)接技術(shù)能處理復(fù)雜幾何形狀的網(wǎng)格生成問題，但對(duì)復(fù)雜外形生成統(tǒng)一的結(jié)構(gòu)化貼體網(wǎng)格仍是一項(xiàng)艱巨的工作。

采用重疊網(wǎng)格技術(shù)[2]，將流動(dòng)區(qū)域分成若干相互重疊的子區(qū)域，對(duì)各不同區(qū)域可單獨(dú)生成拓?fù)漭^簡單的結(jié)構(gòu)化網(wǎng)格。因此，重疊網(wǎng)格技術(shù)有效地降低了復(fù)雜外形網(wǎng)格的生成難度，大大拓展了結(jié)構(gòu)化網(wǎng)格求解技術(shù)的應(yīng)用范圍。同時(shí)，由于在處理運(yùn)動(dòng)部件方面的優(yōu)勢(shì)，重疊網(wǎng)格技術(shù)在載荷分離、旋翼等[3-4]領(lǐng)域得到廣泛應(yīng)用。

重疊網(wǎng)格技術(shù)的關(guān)鍵在于重疊網(wǎng)格的組裝，所謂組裝，包含了切割(也被稱為挖洞)和貢獻(xiàn)單元搜索2個(gè)過程。在切割過程中，需要將與固體區(qū)域重疊的網(wǎng)格單元(格點(diǎn)計(jì)算方式是網(wǎng)格點(diǎn))以及其他不必參與計(jì)算的網(wǎng)格單元標(biāo)記出來。而在貢獻(xiàn)單元搜索過程中，對(duì)于每個(gè)需要插值的網(wǎng)格單元，需要尋找一個(gè)包含它且能提供最佳插值模板的單元作為貢獻(xiàn)單元，并計(jì)算出相對(duì)應(yīng)的插值系數(shù)。然而重疊網(wǎng)格切割和貢獻(xiàn)單元搜索都是費(fèi)時(shí)的尋點(diǎn)過程。

國內(nèi)外學(xué)者對(duì)重疊網(wǎng)格的組裝方法進(jìn)行了大量的研究。Chiu和Meakin[5]提出采用均勻笛卡兒網(wǎng)格的洞映射(Hole Map)方法提高了挖洞效率，但洞映射方法在識(shí)別縫隙處的物面時(shí)需要較大的內(nèi)存。Meakin[6]提出了目標(biāo)X射線法進(jìn)行挖洞，該方法具有較好的魯棒性。PEGASUS 5[7]將洞映射方法和視線法結(jié)合來提高挖洞的效率。Wang等[8]采用ADT樹結(jié)構(gòu)實(shí)現(xiàn)重疊網(wǎng)格組裝的加速。Noack[9]采用二叉樹方法的自適應(yīng)直角網(wǎng)格加速了重疊網(wǎng)格的尋點(diǎn)。此外，為有效減小網(wǎng)格的重疊區(qū)域，提出了擴(kuò)張收縮法[5]、割補(bǔ)法[10]等；同時(shí)，SUGGAR++[11]采用網(wǎng)格體積、壁面距離等標(biāo)準(zhǔn)進(jìn)行了重疊區(qū)域最小化。在國內(nèi)，袁武[12-13]、王文[14-15]、淮洋[16]等改進(jìn)了洞映射方法和割補(bǔ)法，提高了洞映射法的分辨能力和重疊網(wǎng)格的組裝質(zhì)量。

Lee和Baeder[17]提出了一種重疊網(wǎng)格隱式切割方法，該方法采用統(tǒng)一的切割準(zhǔn)則進(jìn)行網(wǎng)格切割，簡化了重疊網(wǎng)格“切割”的概念。Landmann和Montagnac[18]提出采用ADT技術(shù)和并行切割的方法加速隱式切割進(jìn)程。劉秋洪等[19-20]提出了多層嵌套重疊隱式切割方法，將重疊網(wǎng)格組裝分為嵌套切割和重疊切割2個(gè)過程，簡化了貢獻(xiàn)單元選擇的復(fù)雜度。

盡管重疊網(wǎng)格的應(yīng)用在一定程度上減輕了結(jié)構(gòu)化網(wǎng)格生成的難度，然而對(duì)于工程中的復(fù)雜形狀建立重疊網(wǎng)格仍然不是一項(xiàng)容易的任務(wù)。在實(shí)際的工程應(yīng)用中，組裝過程需要較多的人工干預(yù)，并且重疊網(wǎng)格生成往往需要反復(fù)的迭代修改才能完成。尤其對(duì)于多體運(yùn)動(dòng)這類動(dòng)邊界問題，緩慢的組裝速度成為重疊網(wǎng)格在工業(yè)中應(yīng)用的瓶頸。因而研究高效、自動(dòng)的重疊網(wǎng)格組裝方法[21]是提高重疊網(wǎng)格應(yīng)用能力的關(guān)鍵途徑。

由于自適應(yīng)樹結(jié)構(gòu)具有訪問速度快、方便數(shù)據(jù)查詢的特點(diǎn)，所以本文采用二叉樹自適應(yīng)直角網(wǎng)格進(jìn)行切割，采用八叉樹自適應(yīng)直角網(wǎng)格進(jìn)行貢獻(xiàn)單元搜索，采用壁面距離準(zhǔn)則實(shí)現(xiàn)重疊區(qū)域最小化和分層嵌套重疊策略，采用多個(gè)二/八叉樹處理動(dòng)態(tài)問題?；谝陨戏椒?，本文發(fā)展了一種適用于動(dòng)態(tài)問題的快速重疊網(wǎng)格組裝方法。

1 重疊網(wǎng)格區(qū)域切割方法

由于二叉樹的數(shù)據(jù)結(jié)構(gòu)具有高效的尋址能力，因此在二叉樹直角網(wǎng)格中給定坐標(biāo)點(diǎn)能夠快速確定包圍該點(diǎn)的二叉樹葉子節(jié)點(diǎn)。同時(shí)二叉樹的各向異性自適應(yīng)特點(diǎn)特別適合對(duì)幾何形狀的近似解析。因此將其應(yīng)用于物體表面的近似解析和重疊網(wǎng)格切割時(shí)的尋點(diǎn)操作是非常適合的。

在下文中，對(duì)用于切割的二叉樹自適應(yīng)直角網(wǎng)格簡稱為網(wǎng)格切割二叉樹；該直角網(wǎng)格中的每個(gè)單元對(duì)應(yīng)著二叉樹的一個(gè)節(jié)點(diǎn)，簡稱為節(jié)點(diǎn)；根據(jù)節(jié)點(diǎn)之間的關(guān)系，有父節(jié)點(diǎn)、子節(jié)點(diǎn)以及同屬一個(gè)父節(jié)點(diǎn)的2個(gè)兄弟節(jié)點(diǎn)；二叉樹的葉子節(jié)點(diǎn)簡稱為葉子節(jié)點(diǎn)。

1.1 網(wǎng)格切割二叉樹構(gòu)造

在切割階段，首先構(gòu)造二叉樹根節(jié)點(diǎn)覆蓋原始計(jì)算網(wǎng)格區(qū)域，此時(shí)根節(jié)點(diǎn)包含了所有的壁面單元。通過比較節(jié)點(diǎn)網(wǎng)格尺度與所包含的壁面單元的包圍盒尺度，判斷是否要對(duì)該節(jié)點(diǎn)作二分加密和確定二分方向。這樣不斷加密，直到每個(gè)二叉樹葉子節(jié)點(diǎn)的尺度不大于其所包含的壁面單元尺度，或者葉子節(jié)點(diǎn)的尺度已經(jīng)達(dá)到預(yù)先指定的下限，亦或者樹的深度達(dá)到預(yù)設(shè)的上限。這樣就完成了網(wǎng)格切割二叉樹的構(gòu)造，此時(shí)每個(gè)葉子節(jié)點(diǎn)都維護(hù)了一個(gè)與其有相交關(guān)系的壁面單元列表。遍歷每個(gè)葉子節(jié)點(diǎn)，只要它包含了壁面單元，則該葉子節(jié)點(diǎn)被標(biāo)記為WALL節(jié)點(diǎn)。當(dāng)所有WALL節(jié)點(diǎn)被標(biāo)記后，即可采用淹沒填充法分別將固體區(qū)域和流場(chǎng)區(qū)域的葉子節(jié)點(diǎn)標(biāo)記為SOLID和FLOW。

圖1給出了二維三段翼30P30N的網(wǎng)格切割二叉樹，其中每個(gè)矩形區(qū)域代表二叉樹的一個(gè)葉子節(jié)點(diǎn)。圖1(a)中給出了計(jì)算域內(nèi)的葉子節(jié)點(diǎn)(NODE_IN)和固體內(nèi)葉子節(jié)點(diǎn)(NODE_OUT)；圖1(b)給出了壁面葉子節(jié)點(diǎn)(NODE_BND)。從圖中可以看出在壁面處，網(wǎng)格切割二叉樹具有較高的分辨率；而在離開壁面區(qū)域之后，會(huì)采用較粗的樹節(jié)點(diǎn)對(duì)區(qū)域進(jìn)行離散，體現(xiàn)了網(wǎng)格切割二叉樹具有較強(qiáng)的適應(yīng)能力。

1.2 重疊區(qū)域二叉樹切割

但是僅對(duì)二叉樹葉子節(jié)點(diǎn)標(biāo)記FLOW、WALL和SOLID可能造成切割錯(cuò)誤。如圖2所示，圖中給出了物體A的網(wǎng)格(藍(lán)色網(wǎng)格線)及其二叉樹的壁面葉子節(jié)點(diǎn)(紅色矩形)。圖中綠色網(wǎng)格線是另一物體的網(wǎng)格，其網(wǎng)格點(diǎn)P位于物體A的壁面葉子節(jié)點(diǎn)中，所以點(diǎn)P應(yīng)該被標(biāo)記為OUT。由于采用包圍盒對(duì)固體壁面進(jìn)行近似解析處理，造成物體A的網(wǎng)格點(diǎn)Q也位于其壁面葉子節(jié)點(diǎn)中，但是對(duì)點(diǎn)Q進(jìn)行切割顯然是錯(cuò)誤的。為了解決這個(gè)問題，WALL和SOLID葉子節(jié)點(diǎn)還必需記錄該節(jié)點(diǎn)所屬的物體。只有當(dāng)某一物體的網(wǎng)格點(diǎn)落在屬于其他物體的SOLID或者WALL葉子節(jié)點(diǎn)中時(shí)，才能被切割。

本文采用位域變量BBA(Bytes of Body Association)來處理這種問題。首先對(duì)于每個(gè)固體壁面邊界進(jìn)行編號(hào)i=1,2,…,Nw。對(duì)于每個(gè)葉子節(jié)點(diǎn)，使用一個(gè)至少Nw位的BBA來標(biāo)識(shí)它與每個(gè)壁面邊界的面元是否相交。而每個(gè)物體也使用這樣一個(gè)BBA變量來表示該物體的表面由哪些壁面邊界組成。要檢查一個(gè)葉子節(jié)點(diǎn)是否與另一個(gè)物體相交，只需要檢查兩者的BBA位域變量進(jìn)行按位與操作的結(jié)果是否為0即可。如果不為0，說明兩者屬于同一個(gè)物體，不應(yīng)進(jìn)行網(wǎng)格切割。

由于二叉樹的搜索速度僅與其樹深度成正比，所以可以快速地判斷網(wǎng)格點(diǎn)與固體區(qū)域的關(guān)系。根據(jù)葉子節(jié)點(diǎn)的標(biāo)記，將位于FLOW葉子節(jié)點(diǎn)的網(wǎng)格點(diǎn)標(biāo)記為IN，而位于WALL和SOLID葉子節(jié)點(diǎn)中的網(wǎng)格點(diǎn)標(biāo)記為OUT。然后進(jìn)一步完成網(wǎng)格單元的標(biāo)記，將包含OUT網(wǎng)格點(diǎn)的網(wǎng)格單元標(biāo)記為OUT，其余的單元標(biāo)為IN，這樣便完成了重疊區(qū)域網(wǎng)格切割。

2 貢獻(xiàn)單元搜索

當(dāng)網(wǎng)格單元被分類為IN和OUT之后，需要將IN類的網(wǎng)格單元中的插值單元標(biāo)記為FRINGE。沿著預(yù)先定義的重疊邊界和切割確定的洞邊界，直接將邊界內(nèi)兩層單元(因?yàn)椴捎枚A格式計(jì)算)標(biāo)記為FRINGE，即可確定所有插值單元。

在貢獻(xiàn)單元搜索階段，面對(duì)的問題是任意給定物體A中一個(gè)插值單元的中心點(diǎn)p，找到另外一個(gè)物體上包含點(diǎn)p的網(wǎng)格單元D。本文采用最速下降搜索，但是這種方法要求初始單元需接近單元D。所以本文采取等邊八叉樹自適應(yīng)直角網(wǎng)格輔助確定初始單元，以下簡稱八叉樹。

首先構(gòu)造一個(gè)比整個(gè)計(jì)算網(wǎng)格包圍盒稍大的立方體(二維情況下是正方形)作為八叉樹的根節(jié)點(diǎn)。然后根據(jù)其所包含的網(wǎng)格單元尺度不斷將節(jié)點(diǎn)八分加密，形成新的葉子節(jié)點(diǎn)。對(duì)于每個(gè)葉子節(jié)點(diǎn)，記錄其所覆蓋的網(wǎng)格單元的編號(hào)。當(dāng)葉子節(jié)點(diǎn)的尺度小于其覆蓋的網(wǎng)格單元尺度，或者節(jié)點(diǎn)尺度達(dá)到下限，或者包含的網(wǎng)格單元個(gè)數(shù)達(dá)到下限，或者樹的深度達(dá)到上限，則停止對(duì)該節(jié)點(diǎn)加密。此時(shí)每個(gè)葉子節(jié)點(diǎn)包含了一組網(wǎng)格單元編號(hào)。在這些單元中，對(duì)每個(gè)網(wǎng)格塊只保留一個(gè)作為初始單元。但是考慮到如圖3所示的情況，一個(gè)八叉樹節(jié)點(diǎn)(用實(shí)線矩形表示)跨越了一個(gè)薄物體(如翼型后緣)的固體區(qū)域。對(duì)于這樣的節(jié)點(diǎn)，它覆蓋區(qū)域內(nèi)的所有壁面邊界單元都需要被保留作為搜索初始單元。這樣就構(gòu)造了葉子節(jié)點(diǎn)的搜索初始單元列表。完成所有葉子節(jié)點(diǎn)的處理后，八叉樹構(gòu)造完成。

當(dāng)八叉樹構(gòu)造完畢后，便可進(jìn)行包含單元搜索。在給定物體A的網(wǎng)格中，任意插值單元中心點(diǎn)為p，在八叉樹中能夠快速找到p所在的葉子節(jié)點(diǎn)N。選擇搜索初始單元列表中所有的初始單元，進(jìn)行最速下降搜索[22]。通過以上搜索可得到一組非OUT且非FRINGE的包含單元，采用物面距離作為貢獻(xiàn)單元的評(píng)價(jià)標(biāo)準(zhǔn)，在所有的包含單元中選擇物面距離最小者做貢獻(xiàn)單元。

3 重疊區(qū)域最小化方法

前文已經(jīng)指出，切割邊界內(nèi)兩層單元被標(biāo)記為FRINGE。但是這樣做可能導(dǎo)致重疊區(qū)域過大，浪費(fèi)計(jì)算資源，而且影響后處理效果。因此需要調(diào)整重疊邊界以有效地縮小重疊區(qū)域。

重疊區(qū)域進(jìn)行最小化的一個(gè)準(zhǔn)則是被插值單元和貢獻(xiàn)單元的壁面距離應(yīng)該盡量接近。例如在圖4中，兩個(gè)圓柱網(wǎng)格相互重疊，圖4(a)為未進(jìn)行重疊區(qū)域最小化時(shí)的網(wǎng)格組裝情況，具有較大的重疊區(qū)域。顯然在圓柱網(wǎng)格A的壁面附近，網(wǎng)格A中的貢獻(xiàn)單元的壁面距離要比圓柱網(wǎng)格B中作為插值單元的壁面距離值明顯偏??；反之，在圓柱網(wǎng)格B的壁面附近，網(wǎng)格B中的貢獻(xiàn)單元的壁面距離要比圓柱網(wǎng)格A中的插值單元的壁面距離偏小。但是在兩個(gè)圓柱的中間位置附近，網(wǎng)格A和網(wǎng)格B的單元各自的壁面距離是接近的。因此兩圓柱網(wǎng)格的重疊邊界應(yīng)在此位置附近。

在調(diào)整重疊區(qū)域時(shí)，從洞邊界開始在插值區(qū)域中做陣面推進(jìn)，以格心方式進(jìn)行說明。首先引入一個(gè)概念，偽固體(PSOLID)單元，是重疊區(qū)域最小化過程中對(duì)要進(jìn)行最小化的單元的一種標(biāo)記形式，最終被標(biāo)記的單元在經(jīng)過重疊區(qū)域最小化之后將被設(shè)置為OUT。具體來說，如果一個(gè)單元的壁面距離比貢獻(xiàn)單元的壁面距離大，并且貢獻(xiàn)單元插值模版中不包含PSOLID單元，則將這個(gè)插值單元標(biāo)記為PSOLID單元。如果推進(jìn)到某單元的貢獻(xiàn)單元插值模版中包含PSOLID單元的情況，說明在該處兩個(gè)重疊邊界發(fā)生接觸，不應(yīng)該繼續(xù)縮小重疊區(qū)域，并且該單元也不能作為插值單元，不進(jìn)行PSOLID單元標(biāo)記。此時(shí)可以根據(jù)所有的PSOLID單元確定出最小重疊區(qū)域以及重疊邊界，即將與計(jì)算單元鄰近的兩層PSOLID網(wǎng)格單元標(biāo)記為插值單元，并將剩余的PSOLID單元設(shè)置為OUT，即可完成重疊區(qū)域最小化。最后可以得到圖4(b)的結(jié)果。

4 分層嵌套重疊策略

對(duì)于復(fù)雜外形物體進(jìn)行重疊網(wǎng)格組裝時(shí)，多網(wǎng)格重疊的情況會(huì)導(dǎo)致貢獻(xiàn)單元選擇變得復(fù)雜。因此本文采用嵌套重疊網(wǎng)格[19,23]的思想，將各部件的網(wǎng)格分為不同的層次，按照網(wǎng)格間的層次關(guān)系進(jìn)行組裝，從而簡化組裝過程，提高貢獻(xiàn)單元搜索效率和網(wǎng)格組裝質(zhì)量。

4.1 分層嵌套策略

圖5給出了重疊網(wǎng)格的分層嵌套策略圖。在劃分網(wǎng)格時(shí)將較密的第n層網(wǎng)格嵌套在較粗的n-1 層網(wǎng)格內(nèi)部，將n-1層中被第n層覆蓋的網(wǎng)格區(qū)域進(jìn)行切割。而當(dāng)同層的網(wǎng)格出現(xiàn)相互切割時(shí)，需要根據(jù)單元壁面距離選擇距離較小者做計(jì)算單元。采用這種分層嵌套策略，可以按照部件的從屬關(guān)系進(jìn)行分層嵌套或同層重疊，降低對(duì)復(fù)雜外形網(wǎng)格的生成難度，簡化不同層之間貢獻(xiàn)單元的搜索，可生成高質(zhì)量的結(jié)構(gòu)化貼體重疊網(wǎng)格。

重疊網(wǎng)格分層嵌套策略中采用了簇和層2個(gè)物理概念。簇(Cluster)是包含一個(gè)部件或物體的所有網(wǎng)格塊，由一個(gè)或者多個(gè)具有對(duì)接(Match)邊界的網(wǎng)格塊組成，不同的網(wǎng)格簇之間可相互重疊(Overlap)。圖6中的兩個(gè)圓柱網(wǎng)格簇分別由4個(gè)網(wǎng)格塊(Block)組成，兩個(gè)網(wǎng)格簇之間相互重疊。層(Layer)由一個(gè)或幾個(gè)相互重疊的網(wǎng)格簇組成，不同的網(wǎng)格層之間可以相互嵌套(Embed)，圖7為2個(gè)背景網(wǎng)格層和一個(gè)NACA4412翼型網(wǎng)格層之間的相互嵌套關(guān)系。

4.2 分層嵌套重疊策略的統(tǒng)一實(shí)現(xiàn)

本文提出采用基于壁面距離的切割準(zhǔn)則將不同層間的嵌套切割和同層間的重疊切割2個(gè)過程統(tǒng)一于一個(gè)重疊網(wǎng)格切割過程，從而將嵌套重疊的概念用于重疊網(wǎng)格的切割中。整個(gè)過程遵循的基本原則是在網(wǎng)格發(fā)生重疊或嵌套時(shí)，將壁面距離小的網(wǎng)格單元用作計(jì)算單元，而壁面距離大的網(wǎng)格單元被切割。

具體來說，根據(jù)網(wǎng)格簇的概念對(duì)帶有物面的網(wǎng)格直接進(jìn)行壁面距離的計(jì)算；對(duì)于沒有壁面的背景網(wǎng)格，給定一個(gè)固定的壁面距離值。當(dāng)存在多個(gè)背景網(wǎng)格時(shí)，根據(jù)背景網(wǎng)格位于的層數(shù)給定背景網(wǎng)格壁面距離值，較大的壁面距離值表示較低的網(wǎng)格層。這樣通過壁面距離大小的比較就可以快速地將被高層網(wǎng)格區(qū)域覆蓋的低層網(wǎng)格設(shè)置為OUT，實(shí)現(xiàn)分層嵌套切割的效果。同時(shí)由于采用壁面距離實(shí)現(xiàn)了同層重疊區(qū)域的最小化，因而實(shí)現(xiàn)了重疊切割。因此，通過采用壁面距離準(zhǔn)則，可以將同層的重疊切割和不同層的嵌套切割統(tǒng)一地在一個(gè)切割過程中實(shí)現(xiàn)，簡化了重疊網(wǎng)格的切割過程。

5 運(yùn)動(dòng)部件的處理

在這里只考慮剛體運(yùn)動(dòng)的情況。每次物體位置和姿態(tài)變化后，都需要重新進(jìn)行切割和搜索貢獻(xiàn)單元，還需要最小化重疊區(qū)域。盡管基于樹結(jié)構(gòu)的切割與搜索速度很快，但是如果每次重新生成用于切割和貢獻(xiàn)單元搜索的二叉樹和八叉樹，必然導(dǎo)致效率不佳。因此，根據(jù)運(yùn)動(dòng)規(guī)律的不同將所有物體分為若干動(dòng)態(tài)組，而所有靜態(tài)物體被納入唯一的靜態(tài)組。各組分別構(gòu)造獨(dú)立的二叉樹和八叉樹。當(dāng)物體運(yùn)動(dòng)時(shí)，并不對(duì)物體的網(wǎng)格以及二/八叉樹的坐標(biāo)進(jìn)行變換，即各個(gè)動(dòng)態(tài)分組仍然用各自原有的局部坐標(biāo)系，但要記錄下各動(dòng)態(tài)組6個(gè)自由度的位移。在更新切割信息時(shí)，要重復(fù)將被切割的坐標(biāo)點(diǎn)變換到切割物體所在的局部坐標(biāo)內(nèi)。對(duì)于同一組內(nèi)的切割，由于相對(duì)位置不變，只需進(jìn)行一次切割，后續(xù)不再改變。

對(duì)于貢獻(xiàn)單元的搜索，也采取上述方式處理。這樣二叉樹和八叉樹不必重新構(gòu)造，節(jié)省了大量的計(jì)算時(shí)間。

6 分層嵌套重疊網(wǎng)格靜/動(dòng)態(tài)組裝算例

依據(jù)以上提出的重疊網(wǎng)格組裝方法開發(fā)了一套重疊網(wǎng)格組裝程序，采用二維三段翼30P30N算例[24]和ADEC三維投彈算例[25]分別驗(yàn)證該方法的靜態(tài)和動(dòng)態(tài)網(wǎng)格組裝能力。

6.1 二維三段翼30P30N算例

本文采用二維三段翼30P30N來驗(yàn)證重疊網(wǎng)格靜態(tài)組裝的可靠性。依據(jù)其幾何外形的結(jié)構(gòu)關(guān)系，劃分各部件不同層次的流場(chǎng)區(qū)域網(wǎng)格，采用由6個(gè)簇組成的4個(gè)網(wǎng)格層進(jìn)行重疊網(wǎng)格組裝。第4層(最高層)為翼型網(wǎng)格，其中包含三段翼前緣縫翼、主翼和后緣襟翼3個(gè)網(wǎng)格簇，對(duì)每個(gè)網(wǎng)格簇分別生成正交的結(jié)構(gòu)化網(wǎng)格；其余3層均為直角背景網(wǎng)格，第3層為翼型附近流場(chǎng)的網(wǎng)格，第2層為保證背景網(wǎng)格之間良好過度而采用的一層網(wǎng)格，第1層為包含遠(yuǎn)場(chǎng)邊界的背景網(wǎng)格。

在三段翼30P30N中，縫翼與主翼以及主翼與襟翼之間的復(fù)雜流動(dòng)是整個(gè)數(shù)值模擬的關(guān)鍵，因此對(duì)這個(gè)區(qū)域的網(wǎng)格進(jìn)行了加密。網(wǎng)格組裝結(jié)果如圖8所示，其中圖8(a)為三層背景網(wǎng)格的嵌套切割結(jié)果，圖8(b)為前緣縫翼、主翼和后緣襟翼網(wǎng)格重疊切割結(jié)果。圖9和圖10分別給出了前緣縫翼、后緣襟翼與主翼縫隙處重疊網(wǎng)格的切割結(jié)果。重疊組裝后的網(wǎng)格中無網(wǎng)格孤島，說明了本文發(fā)展的重疊網(wǎng)格靜態(tài)組裝算法可靠。

6.2 ADEC三維投彈算例

為了進(jìn)一步驗(yàn)證本文算法對(duì)動(dòng)態(tài)重疊網(wǎng)格的組裝能力，選取ADEC三維投彈算例進(jìn)行驗(yàn)證。為了便于檢查動(dòng)態(tài)網(wǎng)格組裝結(jié)果，本文并沒有進(jìn)行CFD計(jì)算，而是通過python腳本向組裝程序發(fā)送位移數(shù)據(jù)?？紤]到投彈過程主要發(fā)生在豎直方向上，因此以導(dǎo)彈豎直下落過程來測(cè)試重疊網(wǎng)格的動(dòng)態(tài)組裝。圖11給出了機(jī)翼-掛架和導(dǎo)彈的物面網(wǎng)格。該算例中機(jī)翼的網(wǎng)格量約為365萬，導(dǎo)彈的網(wǎng)格量約為241萬。圖12給出了該算例的網(wǎng)格切割二叉樹的壁面葉子節(jié)點(diǎn)，從圖中可以看出采用網(wǎng)格切割二叉樹能有效地分辨出算例中掛架和導(dǎo)彈間的狹小縫隙。

該算例在一臺(tái)E5-2660處理器，內(nèi)存64 G，主頻為2.60 GHz的工作站上運(yùn)行，并采用單線程方式組裝網(wǎng)格。每步更新切割和貢獻(xiàn)單元搜索，只輸出組裝信息文件，以節(jié)省磁盤操作時(shí)間。

圖13給出了初始時(shí)刻投彈模型網(wǎng)格組裝結(jié)果，圖中藍(lán)色網(wǎng)格表示導(dǎo)彈的計(jì)算網(wǎng)格單元，紅色網(wǎng)格是機(jī)翼的計(jì)算網(wǎng)格單元。圖13(a)是xOy視圖上位于掛架中心處的一個(gè)切面，從圖中可以看出對(duì)兩者之間縫隙處的網(wǎng)格進(jìn)行了良好的組裝；圖13(b)給出了yOz方向上機(jī)翼和導(dǎo)彈之間縫隙處網(wǎng)格組裝的詳細(xì)結(jié)果，從圖中可以看出兩者之間的縫隙處有足夠多的網(wǎng)格單元，可以保證插值的精度，能有效地計(jì)算縫隙中的復(fù)雜流動(dòng)。

圖14為運(yùn)動(dòng)過程中4個(gè)不同時(shí)刻的網(wǎng)格組裝結(jié)果，圖中給出的是yOz和xOy方向上的切面, dt表示時(shí)間步長。從下落過程的組裝結(jié)果可以看出，在初始時(shí)刻機(jī)翼和導(dǎo)彈相距較近，兩者的網(wǎng)格經(jīng)過重疊最小化之后縫隙處距離壁面較近的網(wǎng)格都會(huì)被切割；隨著導(dǎo)彈遠(yuǎn)離掛架之后，導(dǎo)彈附近的網(wǎng)格都會(huì)被保留，這樣組裝的網(wǎng)格可以保證導(dǎo)彈周圍流場(chǎng)的精確計(jì)算。

按照是否進(jìn)行重疊區(qū)域最小化操作分別進(jìn)行時(shí)間統(tǒng)計(jì)，如表1所示。整個(gè)網(wǎng)格組裝過程總共30個(gè)時(shí)間步，在不進(jìn)行重疊最小化操作時(shí)，總共耗時(shí)437 s，平均每步用時(shí)約14.6 s；進(jìn)行重疊區(qū)域最小化操作的總共耗時(shí)為1 168 s，平均每步用時(shí)約39 s。其中，由于在初始時(shí)間步上會(huì)對(duì)網(wǎng)格數(shù)據(jù)進(jìn)行初始化和計(jì)算網(wǎng)格單元的壁面距離，所以初始時(shí)刻網(wǎng)格組裝的用時(shí)相對(duì)較多，當(dāng)導(dǎo)彈離開掛架一定的距離后，每步網(wǎng)格組裝時(shí)間會(huì)顯著下降。

Table 1 Time consumption of assembling overset grids for moving bodies

TypeofassemblingmethodTimeconsumption/sTotalAverageThefirststepWithoutoverlappingminimization43714.621.6Withoverlappingminimization116839108

另一方面，雖然采用重疊區(qū)域最小化會(huì)一定程度增加網(wǎng)格組裝的時(shí)間，但提高了網(wǎng)格的組裝質(zhì)量。從網(wǎng)格量、時(shí)間消耗和網(wǎng)格組裝質(zhì)量上來看，兩種方式的動(dòng)態(tài)重疊網(wǎng)格組裝效率均可以滿足運(yùn)動(dòng)物體非定常計(jì)算的要求。

7 結(jié) 論

本文采用二/八樹結(jié)構(gòu)的自適應(yīng)直角網(wǎng)格方法，發(fā)展了一種動(dòng)態(tài)重疊網(wǎng)格快速組裝方法，該方法只需要初始的組裝網(wǎng)格和網(wǎng)格邊界條件，無需人工干預(yù)便可自動(dòng)實(shí)現(xiàn)網(wǎng)格組裝。

1) 通過采用二/八樹結(jié)構(gòu)的自適應(yīng)直角網(wǎng)格方法，提高了切割的效率并節(jié)約內(nèi)存空間，實(shí)現(xiàn)了高效搜索和定位貢獻(xiàn)單元。提出的重疊區(qū)域最小化方法顯著減小了網(wǎng)格重疊區(qū)域，提高了網(wǎng)格組裝質(zhì)量，有利于提高流動(dòng)計(jì)算的效率。

2) 采用壁面距離準(zhǔn)則的分層嵌套重疊切割策略，在一個(gè)過程中同時(shí)實(shí)現(xiàn)重疊切割和嵌套切割功能，簡化了貢獻(xiàn)單元的搜索，提高了復(fù)雜外形網(wǎng)格切割效率和網(wǎng)格組裝質(zhì)量。

3) 對(duì)不同運(yùn)動(dòng)部件分別采用隨部件一起運(yùn)動(dòng)的獨(dú)立的二/八叉樹結(jié)構(gòu)，能夠有效地降低計(jì)算冗余，加快動(dòng)態(tài)重疊網(wǎng)格的組裝速度。

4) 算例說明本文提出的重疊網(wǎng)格組裝方法快速高效，提升了網(wǎng)格組裝效率，可解決帶有運(yùn)動(dòng)邊界的非定常流動(dòng)數(shù)值模擬過程中的動(dòng)網(wǎng)格生成問題。

[1] LEVY D W, LAFLIN K R, TINOCO E N, et al. Summary of data from the fifth computational fluid dynamics drag prediction workshop[J]. Journal of Aircraft， 2014, 51(4): 1194-1213.

[2] BENEK J A, STEGER J L, DOUGHERTY F C. A flexible grid embedding technique with application to the Euler equations：AIAA-1983-1944[R]. Reston: AIAA, 1983.

[3] NOACK R W, SLOTNICK J P. A summary of the 2004 overset symposium on composite grids and solution technology：AIAA-2005-0921[R]. Reston: AIAA, 2005.

[4] CHAN W M. Overset grid technology development at NASA Ames Research Center[J]. Computers & Fluids， 2009, 38(3): 496-503.

[5] CHIU I T, MEAKIN R L. On automating domain connectivity for overset grids：AIAA-1995-0854[R]. Reston: AIAA, 1995.

[6] MEAKIN R L. Object X-rays for cutting holes in composite overset structured grids：AIAA-2001-2537[R]. Reston: AIAA, 2001.

[7] ROGERS S E, SUHS N E, DIETZ W E. PEGASUS 5: An automated preprocessor for overset-grid computational fluid dynamics[J]. AIAA Journal， 2003, 41(6): 1037-1045.

[8] WANG Z J, PARTHASARATHY V, HARIHARAN N. A fully automated Chimera methodology for multiple moving body problems：AIAA-1998-0217[R]. Reston: AIAA, 1998.

[9] NOACK R W. SUGGAR: A general capability for moving body overset grid assembly：AIAA-2005-5117[R]. Reston: AIAA, 2005.

[10] CHO K W, KWON J H, LEE S. Development of a fully systemized chimera methodology for steady/unsteady problems[J]. Journal of Aircraft， 1999, 36(6): 973-980.

[11] NOACK R W, BOGER D A, KUNZ R F, et al. Suggar++: An improved general overset grid assembly capability：AIAA-2009-3992[R]. Reston: AIAA, 2009.

[12] 袁武， 閻超， 席柯. 洞映射方法的研究和改進(jìn)[J]. 北京航空航天大學(xué)學(xué)報(bào), 2012, 38(4): 563-568. YUAN W， YAN C， XI K. Investigation and enhancement of hole mapping method[J]. Journal of Beijing University of Aeronautics and Astronautics, 2012, 38(4): 563-568 (in Chinese).

[13] 袁武， 閻超， 楊威. 割補(bǔ)法的改進(jìn)和應(yīng)用[J]. 空氣動(dòng)力學(xué)學(xué)報(bào), 2013, 31(6): 704-709. YUAN W， YAN C， YANG W. Enhancement and application of cut-paste algorithm[J]. Acta Aerodynamica Sinica, 2013, 31(6): 704-709 (in Chinese).

[14] 王文， 閻超， 袁武， 等. 新型重疊網(wǎng)格洞面優(yōu)化方法及其應(yīng)用[J]. 航空學(xué)報(bào), 2016, 27(3): 826-835. WANG W， YAN C， YUAN W， et al. A new kind of overlap optimization algorithm and its applica-tions[J]. Acta Aeronautica et Astronautica Sinica, 2016, 27(3): 826-835 (in Chinese).

[15] 王文， 閻超， 袁武， 等. 魯棒的結(jié)構(gòu)網(wǎng)格自動(dòng)化重疊方法[J]. 航空學(xué)報(bào), 2016， 37(10): 2980-2991. WANG W， YAN C， YUAN W， et al. A robust and automatic structured overlapping grid approach[J]. Acta Aeronautica et Astronautica Sinica, 2016， 37(10): 2980-2991 (in Chinese).

[16] 淮洋， 郝海兵， 姚冰， 等. 一種改進(jìn)型洞映射法[J]. 航空計(jì)算技術(shù), 2015, 45(2): 31-34. HUAI Y， HAO H B， YAO B， et al. An improved method of hole-map[J]. Aeronautical Computing Technique, 2015, 45(2): 31-34 (in Chinese).

[17] LEE Y, BAEDER J D. Implicit hole cutting-A new approach to overset grid connectivity: AIAA-2003-4128[R]. Reston: AIAA, 2003.

[18] LANDMANN B, MONTAGNAC M. A highly automated parallel Chimera method for overset grids based on the implicit hole cutting technique[J]. International Journal for Numerical Methods in Fluids, 2011, 66(6): 778-804.

[19] 劉秋洪， 屈崑， 蔡晉生， 等. 嵌套重疊網(wǎng)格的構(gòu)造策略及其隱式切割[J]. 中國科學(xué): 物理學(xué) 力學(xué) 天文學(xué), 2013, 42(2): 186-198. LIU Q H， QU K， CAI J S， et al. An automated Chimera method based on a hierarchical overset grid strategy and the implicit hole cutting technique[J]. SCIENTIA SINICA Physica, Mechanica & Astronomica, 2013, 42(2): 186-198 (in Chinese).

[20] 徐嘉， 劉秋洪， 蔡晉生， 等. 基于隱式嵌套重疊網(wǎng)格技術(shù)的阻力預(yù)測(cè)[J]. 航空學(xué)報(bào), 2013, 34(2): 208-217. XU J， LIU Q H， CAI J S， et al. Drag prediciton based on overset grids with implict hole cutting technique[J]. Acta Aeronautica et Astronautica Sinica, 2013, 34(2): 208-217 (in Chinese).

[21] DRUYOR JR C T, JONES W T, KARMAN JR S L. A survey of overset domain assembly methods: AIAA-2015-0910[R]. Reston: AIAA, 2015.

[22] DIETZ W E, JACOCKS J L, FOX J H. Application of domain decomposition to the analysis of complex aerodynamic configurations[C]//3rd International Symposium on Domain Decomposition Methods for Partial Differential Equations. Houston, Texas: SIAM, 1989: 428-450.

[23] CAI J, TSAI H M, LIU F. A parallel viscous flow solver on multi-block overset grids[J]. Computers & Fluids, 2006, 35(10): 1290-1301.

[24] CHIN V D, PETERS D W, SPAID F W, et al. Flowfield measurements about a multi-element airfoil at high Reynolds numbers: AIAA-1993-3137[R]. Reston: AIAA, 1993.

[25] HEIM E R. CFD wing/pylon/finned store mutual interference wind tunnel experiment:ADB152669[R]. Arnold: Arnold Engineering Development Center, 1991.

(責(zé)任編輯：李明敏)

＊Corresponding author. E-mail: kunqu@nwpu.edu.cn

A dynamic assembling method based on adaptive tree structure for hierarchical overset grids

LI Xiaodong, QU Kun*, CAI Jinsheng

SchoolofAeronautics,NorthwesternPolytechnicalUniversity,Xi’an710072,China

The overset or chimera grid is an effective method for large-scale numerical simulation of aerodynamic problems with complex geometry, especially those with bodies in relative motion (such as the rotorcraft, store/aircraft separation). The method of adaptive cartesian grids based on the tree structure is applied to perform fast hole cutting and donor searching in the overset grid assembly process. In this technique, the geometry is represented by adaptive cartesian grids of binary trees, and this can achieve quick locating in the hole cutting. Adaptive cartesian grids of octrees are used to decompose the flow field, making the donor searching efficient. By means of a minimal overlapping-domain approach and a hierarchical strategy for overset grids, which both employ the rule of wall distance, the efficiency and quality of overset grid assembly are improved. For dynamic problems of moving bodies, multiple binary trees and octrees are applied to decrease the redundant work of updating data in procedures of overset grid assembly, substantially reducing the time consumption. The test cases demonstrate that the technique in this paper is computationally efficient and can be successfully employed to problems of multiple moving bodies.

overset grid; binary trees; octrees; adaptive cartesian grid; dynamic grid; unsteady flow; computational fluid dynamics

2016-03-21； Revised：2016-04-12； Accepted：2016-04-27； Published online：2016-05-04 14:16

URL：www.cnki.net/kcms/detail/11.1929.V.20160504.1416.010.html

National Basic Research Program of China

http://hkxb.buaa.edu.cn hkxb@buaa.edu.cn

10.7527/S1000-6893.2016.0135

2016-03-21； 退修日期：2016-04-12； 錄用日期：2016-04-27； 網(wǎng)絡(luò)出版時(shí)間：2016-05-04 14:16

www.cnki.net/kcms/detail/11.1929.V.20160504.1416.010.html

國家“973”計(jì)劃

＊通訊作者.E-mail: kunqu@nwpu.edu.cn

李曉東， 屈崑， 蔡晉生． 分層嵌套重疊網(wǎng)格自適應(yīng)樹結(jié)構(gòu)動(dòng)態(tài)組裝方法[J]． 航空學(xué)報(bào), 2017, 38(3): 120243. LI X D, QU K, CAI J S. A dynamic assembling method based on adaptive tree structure for hierarchical overset grids[J]. Acta Aeronautica et Astronautica Sinica, 2017, 38(3): 120243.

V211.3

A

1000-6893(2017)03-120243-10