并行化非結(jié)構(gòu)重疊網(wǎng)格隱式裝配技術(shù)

常興華，馬戎，張來平1,

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

重疊網(wǎng)格技術(shù)最早用于解決復(fù)雜構(gòu)型的結(jié)構(gòu)網(wǎng)格生成難題[1]，后來逐漸推廣應(yīng)用于非結(jié)構(gòu)網(wǎng)格以及運動邊界問題。目前已經(jīng)成為處理運動邊界問題的有效方法之一，在多體分離、旋翼等問題中應(yīng)用廣泛[2]，并已經(jīng)形成了一系列有代表性的軟件[3]。

網(wǎng)格裝配是重疊網(wǎng)格的核心流程，其目的是將網(wǎng)格中不參與計算的網(wǎng)格點(單元)刪除掉(稱為“挖洞”過程)，并且確定出重疊區(qū)域內(nèi)單元的插值關(guān)系。常見的重疊網(wǎng)格裝配過程是將挖洞過程與確定插值關(guān)系的過程分開處理。在挖洞過程中，需要根據(jù)挖洞曲面判斷點(單元)是否在曲面內(nèi)部，常見的判斷方法有矢量判別法[4]、射線求交法[5]、洞映射方法[6]、目標(biāo)X射線法[7]等。在確定初始的洞邊界之后，一般還要采取一些優(yōu)化措施如割補(bǔ)法[8-9]對洞邊界進(jìn)行優(yōu)化，使其遠(yuǎn)離壁面。“挖洞”結(jié)束之后，再通過查詢算法確定出插值點(單元)的貢獻(xiàn)單元。

理論上非結(jié)構(gòu)重疊網(wǎng)格可以沿用結(jié)構(gòu)網(wǎng)格上成熟的網(wǎng)格裝配技術(shù)，且由于非結(jié)構(gòu)網(wǎng)格更加適用于復(fù)雜構(gòu)型，在處理復(fù)雜構(gòu)型時不需要像結(jié)構(gòu)網(wǎng)格那樣生成很多的網(wǎng)格塊，因此相比結(jié)構(gòu)網(wǎng)格自動化程度更高。Nakahashi等[10]最早開展了非結(jié)構(gòu)重疊網(wǎng)格相關(guān)的研究工作，他們通過查詢節(jié)點的宿主單元并通過對比其壁面距離判斷節(jié)點的屬性(簡稱“壁面距離準(zhǔn)則”)，其實質(zhì)上是一種“隱式”挖洞方法(Implict Hole-Cutting, IHC，其概念之后由Lee和Baeder[11]在2003年提出)。Togashi等[12]進(jìn)一步將該方法推廣應(yīng)用于復(fù)雜的多體相對運動問題，并可以處理多體分離過程中物體相交的情況。Loehner[13]、Luo[14]等在Nakahashi方法的基礎(chǔ)上，采用單元的尺度和壁面距離的組合量作為單元屬性的判斷標(biāo)準(zhǔn)，使插值單元和貢獻(xiàn)單元的大小匹配，有助于提高插值穩(wěn)定性并減少插值誤差。國內(nèi)，田書玲[15]較早開展了非結(jié)構(gòu)重疊網(wǎng)格技術(shù)研究，通過改進(jìn)單元查詢算法提高了重疊網(wǎng)格的裝配效率，并將其應(yīng)用于多體分離問題的數(shù)值模擬。

這種基于單元品質(zhì)的非結(jié)構(gòu)重疊網(wǎng)格隱式裝配方法具有非常高的自動化程度，裝配效果較好，但是其仍然存在如下幾個瓶頸問題：

1) 非結(jié)構(gòu)網(wǎng)格本身占用內(nèi)存較大，由于在重疊網(wǎng)格裝配過程中需要存儲額外的拓?fù)湫畔?，因此?dāng)網(wǎng)格規(guī)模增加后內(nèi)存消耗量會急劇增加。

2) 查詢計算量較大。盡管可以采用交替數(shù)字樹(Alternation Digital Tree，ADT)數(shù)據(jù)結(jié)構(gòu)[16]或者一些基于網(wǎng)格拓?fù)湫畔⒌目焖俨樵兎椒ǎ怯嬎阈屎惋@式挖洞技術(shù)相比仍然有較大差距，并且查詢的計算量會隨著網(wǎng)格規(guī)模的增加而急劇增大。

3) 由于挖洞過程是在查詢和判斷中自動完成的，有可能會產(chǎn)生一些 “孤島”從而導(dǎo)致計算域不是單連通域。

隨著CFD技術(shù)的不斷發(fā)展，工程領(lǐng)域?qū)FD的期望也越來越高。在航空航天以及民用領(lǐng)域，越來越多的工程應(yīng)用希望CFD能夠采用更大規(guī)模的計算網(wǎng)格，模擬更真實的復(fù)雜外形，并能夠更好地模擬分離流動、湍流流動現(xiàn)象。為了適應(yīng)超大規(guī)模的計算網(wǎng)格，并行化裝配是一個必然的選擇。

在并行化重疊網(wǎng)格裝配技術(shù)方面，Zhang和Owens[17]發(fā)展了并行化的非結(jié)構(gòu)重疊網(wǎng)格自動裝配技術(shù)，通過判斷單元的棱線是否和壁面相交確定出洞內(nèi)點，裝配過程分為并行挖洞和并行查詢兩個過程。Sitaraman等[18]開發(fā)了并行化的非結(jié)構(gòu)重疊網(wǎng)格隱式裝配軟件PUNDIT，通過將查詢區(qū)域分割成均勻的六面體區(qū)域，提高了查詢效率。Roget和Sitaraman[19]進(jìn)一步改進(jìn)了PUNDIT的裝配算法，提高了其自動化程度和魯棒性。目前PUNDIT已經(jīng)在美國CREAT-AV項目[20]的Helios系統(tǒng)中進(jìn)行了集成。Mavriplis研究團(tuán)隊[21]則針對非結(jié)構(gòu)網(wǎng)格的間斷-伽遼金格式，開發(fā)了適用于高階格式的非結(jié)構(gòu)重疊網(wǎng)格并行裝配軟件TIOGA。Zagaris[22]、Landmann[23]等針對結(jié)構(gòu)重疊網(wǎng)格，亦開展了并行裝配方面的研究工作。國內(nèi)，梁姍等[24]發(fā)展了一種并行化且具有局部數(shù)據(jù)特性的結(jié)構(gòu)網(wǎng)格挖洞技術(shù)，具有較好的并行計算效率。

從國內(nèi)外的相關(guān)研究工作來看，并行化重疊網(wǎng)格裝配技術(shù)研究剛起步不久，且由于重疊網(wǎng)格裝配方法自身的多樣性，其采用的并行化策略也不盡相同。然而，各種并行化重疊網(wǎng)格裝配技術(shù)的最終目標(biāo)是一致的，那就是自動化、魯棒、高效。就隱式裝配方法而言，其具有自動化程度高、裝配效果好的特點，與非結(jié)構(gòu)網(wǎng)格自身的特點非常吻合。但是其計算量較大，效率較低，從國外相關(guān)研究工作來看，相應(yīng)的并行計算技術(shù)仍在不斷發(fā)展完善之中，國內(nèi)更是缺乏相關(guān)的研究工作。為了適應(yīng)未來超大規(guī)模網(wǎng)格計算的需求，非常必要針對并行化的隱式裝配技術(shù)開展研究。

本文針對格心型的有限體積算法，發(fā)展了一種并行化的非結(jié)構(gòu)重疊網(wǎng)格隱式裝配技術(shù)。該技術(shù)采用壁面距離準(zhǔn)則判斷點的屬性，并通過物理邊界推進(jìn)的方法填充活躍區(qū)域，不需要進(jìn)行是否在物體內(nèi)部的判斷，并且能夠避免出現(xiàn)“孤島”問題。采用網(wǎng)格分區(qū)的方式實現(xiàn)了裝配技術(shù)的并行化，在每一個分區(qū)內(nèi)建立單獨的ADT子樹用于查詢，減小了內(nèi)存消耗量。二維以及三維的測試算例證明了該方法的有效性，并且展示了該方法在超大規(guī)模重疊網(wǎng)格計算中的應(yīng)用前景。

1 非結(jié)構(gòu)重疊網(wǎng)格裝配技術(shù)

本文的重疊網(wǎng)格技術(shù)以格心型的非結(jié)構(gòu)有限體積算法為背景，裝配流程如圖1所示。首先將網(wǎng)格點的類型劃分為活躍(1屬性)和非活躍(0屬性)，然后再根據(jù)活躍點判斷出活躍單元(1屬性)，最后搜索出插值單元(-1屬性)并查詢其貢獻(xiàn)單元。在進(jìn)行網(wǎng)格生成時，根據(jù)實際問題，圍繞每一個物體或者部件生成獨立的網(wǎng)格塊(block)。對于嵌入在大網(wǎng)格內(nèi)的子網(wǎng)格塊，需要將其外邊界設(shè)定為重疊邊界條件。在網(wǎng)格裝配的預(yù)處理階段， 計算出每一個block內(nèi)網(wǎng)格的單元體積、單元中心點、面的面積、面的法向、面的中心點等幾何信息，重構(gòu)出計算所需的“點-體”和“體-點”拓?fù)湫畔?。之后求出所有網(wǎng)格點在該block內(nèi)的壁面距離，并在每一個block內(nèi)建立后續(xù)查詢所需要的ADT數(shù)據(jù)結(jié)構(gòu)。

圖1 非結(jié)構(gòu)重疊網(wǎng)格裝配流程Fig.1 Assembling procedure of unstructured overset grids

1.1 判斷點的類型

本文采用了和文獻(xiàn)[12]類似的基于單元(點)品質(zhì)的“隱式”挖洞方法。單元的品質(zhì)有多種取法，為了提高裝配的效果，文獻(xiàn)[12-13]采用壁面距離作為判斷準(zhǔn)則。文獻(xiàn)[14]則采用了單元尺度和壁面距離的組合量。文獻(xiàn)[15]的研究表明，以網(wǎng)格點的最小壁面距離作為標(biāo)準(zhǔn)，通過結(jié)合一定的洞邊界優(yōu)化措施，可以達(dá)到相當(dāng)不錯的裝配效果，因此本文選擇點的壁面距離作為判斷標(biāo)準(zhǔn)。

在預(yù)處理階段，已經(jīng)求出了各個block內(nèi)的網(wǎng)格點距離該block內(nèi)固壁邊界的距離。對于沒有固壁邊界的block，則直接給定一個固定的距離(這個給定的參數(shù)會影響該block內(nèi)活躍區(qū)域的大小)。之后搜索網(wǎng)格點在鄰居block內(nèi)的宿主單元，并通過插值求出其在鄰居block內(nèi)的壁面距離。對于多個block相互重疊的情況，某網(wǎng)格點可能存在多個宿主單元，此時應(yīng)當(dāng)選擇其中壁面距離的最小值。最后，根據(jù)網(wǎng)格點在本block內(nèi)的壁面距離dist_1以及其在鄰居塊內(nèi)的最小壁面距離dist_2確定其類型：如果dist_1大于dist_2，則該點屬性為0，否則屬性為1。

通過以上操作，可以確定出一部分網(wǎng)格點的屬性，但是仍然有一些網(wǎng)格點由于查找不到宿主單元而無法確定其屬性。這些剩余點有兩種可能：一種是處于非重疊區(qū)內(nèi)，其屬性應(yīng)當(dāng)為1；另一種是處在另外一個物體內(nèi)部，屬性應(yīng)當(dāng)為0。此時文獻(xiàn)中常用的方法是進(jìn)行是否在物體內(nèi)部的邏輯判斷，以便進(jìn)行進(jìn)一步的區(qū)分。本文則采用了一種物理邊界推進(jìn)的方法確定這些剩余點的類型：物理邊界上的點必定為活躍點，其屬性設(shè)置為1，然后根據(jù)“點-體”和“體-點”的拓?fù)潢P(guān)系逐步向計算域內(nèi)推進(jìn)并擴(kuò)大活躍點的范圍，直到推進(jìn)到0類型的點為止。在具體的推進(jìn)過程中采用了如下的循環(huán)迭代方式：針對所有單元進(jìn)行循環(huán)判斷，如果單元中的一個節(jié)點為活躍點，則單元屬性設(shè)置為活躍，其所有節(jié)點也設(shè)置為活躍，并立即生效且應(yīng)用到下一步循環(huán)。采用這種循環(huán)方式可以將物理邊界點的屬性快速拓展至計算域內(nèi)場，并且可以保證計算區(qū)域的單連通特性，避免了重疊網(wǎng)格隱式裝配方法中常見的“孤島”問題。

1.2 判斷單元的類型

在判斷出點的類型之后，再根據(jù)“體-點”的網(wǎng)格拓?fù)湫畔⑴袛喑鰡卧愋?。本文采用了如下判斷?zhǔn)則：一個單元的所有點為活躍點(1屬性)，則該單元為活躍單元(1屬性)；如果該單元的所有點為非活躍點，則該單元為非活躍單元(0屬性)；剩下的單元為插值單元(-1屬性)。

最后搜索插值單元的宿主網(wǎng)格塊和宿主單元。對于存在多個block相互重疊的情況，有可能會查找到多個宿主單元，此時仍然采用和1.1節(jié)中相同的做法，選取壁面距離最近的那個宿主單元建立最終的插值關(guān)系。

1.3 提高算法魯棒性的措施

為了提高非結(jié)構(gòu)重疊網(wǎng)格應(yīng)對復(fù)雜工程問題的能力，本文采用了一些提高網(wǎng)格裝配魯棒性和效率的措施。

1) 對插值邊界進(jìn)行優(yōu)化。采用非結(jié)構(gòu)網(wǎng)格中“點-體”和“體-點”的拓?fù)潢P(guān)系拓展1～2排活躍點，適當(dāng)擴(kuò)大活躍區(qū)域，使得插值單元完全落在鄰居zone的活躍區(qū)域內(nèi)。

2) 對于嵌入在大網(wǎng)格內(nèi)的子網(wǎng)格塊，需要在其重疊邊界上預(yù)留一排單元和一排點作為緩沖區(qū)，緩沖區(qū)內(nèi)的點和單元不參與宿主單元查詢。

3) 在重疊網(wǎng)格裝配的過程中引入容差ε，其大小應(yīng)當(dāng)小于最小的網(wǎng)格點間距(如最小網(wǎng)格間距的1/10)。在計算單元的“最小盒子”以及判斷點是否在單元內(nèi)時，留一個大小為ε的裕量，以防止由于某些極端原因而查找不到宿主單元的情況(例如點正好處在某個單元的邊界上，或者對于存在對稱面的情況)。

4) 為了提高查詢效率，采用ADT數(shù)據(jù)結(jié)構(gòu)用于宿主單元查詢。對于每一個網(wǎng)格塊，以圍繞每一個網(wǎng)格單元的最小盒子作為建立ADT樹的依據(jù)。

2 基于網(wǎng)格分區(qū)的并行通訊策略及裝配效率

基于第1節(jié)的非結(jié)構(gòu)重疊網(wǎng)格自動裝配方法，采用分區(qū)并行的策略，將每一個block分為若干個網(wǎng)格區(qū)(zone)，每一個進(jìn)程負(fù)責(zé)一個或者若干個zone的計算。和未分區(qū)情況下的裝配過程類似：對于每一個zone內(nèi)的網(wǎng)格點，在其他所有zone內(nèi)進(jìn)行查詢并尋找最小壁面距離，從而判斷點的屬性。

并行裝配的核心內(nèi)容是查詢所有zone內(nèi)網(wǎng)格點的所有宿主單元，根據(jù)查詢方式不同有兩種并行策略：第1種策略是每個進(jìn)程僅對其所負(fù)責(zé)的zone內(nèi)的網(wǎng)格點進(jìn)行查詢，然后進(jìn)行點和單元類型的區(qū)分。這種方法比較直觀，操作流程和第1節(jié)相同，易于實現(xiàn)。但是這種并行模式需要在本進(jìn)程內(nèi)存儲一個全局的計算網(wǎng)格，并且同時需要全局網(wǎng)格的幾何信息和拓?fù)湫畔?，?nèi)存占用較大。此外，搜集這個全局計算網(wǎng)格的過程需要較多的數(shù)據(jù)通訊，必然會影響并行效率。特別對于運動網(wǎng)格的問題，每一步重構(gòu)都需要對全局計算網(wǎng)格的信息進(jìn)行通訊，會影響到整個非定常計算的效率。因此，本文采用了第2種策略：每一個進(jìn)程都對所有計算點進(jìn)行查詢，但是只在本進(jìn)程所負(fù)責(zé)的zone內(nèi)查詢貢獻(xiàn)單元，最終的查詢結(jié)果通過各個進(jìn)程之間的并行通訊來確定。圖2給出了兩種并行模式的對比示意圖。

對于圖2(b)所示的第2種策略，其流程和圖1 所示的串行情況下的流程是類似的，首先各個進(jìn)程要對其所負(fù)責(zé)的zone進(jìn)行預(yù)處理操作，之后再分別進(jìn)行點和單元類型的判斷。

圖2 兩種并行化的裝配策略Fig.2 Two approaches for parallel assembling

在判斷點的類型時，需要先搜集各個zone的網(wǎng)格點坐標(biāo)并進(jìn)行廣播通訊，在每個進(jìn)程內(nèi)形成整體的待查詢點集，然后在本進(jìn)程所負(fù)責(zé)的zone內(nèi)查詢該點集的宿主單元并計算壁面距離。點在所有鄰居zone內(nèi)的最小壁面距離最終通過各個進(jìn)程間的并行通訊來完成。

并行情況下同樣通過物理邊界推進(jìn)來填充活躍區(qū)域。但由于每一個block被分解成了若干個zone，每推進(jìn)一次，就需要在分區(qū)邊界上進(jìn)行單元屬性的通訊。圖3所示為分區(qū)并行情況下物理邊界推進(jìn)過程示意圖，其仍然借助了“點-體”和“體-點”的拓?fù)湫畔?，其中從點到體的過程與串行情況一致，唯一不同的是在從體到點的推進(jìn)過程：分區(qū)并行情況下需要首先對分區(qū)邊界面上的單元屬性進(jìn)行通訊，將鄰居zone的單元屬性儲存在本zone分區(qū)邊界上的虛擬單元內(nèi)，之后借助分區(qū)邊界面的“面-點”關(guān)系使屬性順利拓展到邊界點。在插值邊界優(yōu)化以及判斷單元屬性的操作中，需要用到圖3所示的通訊過程，因此將該通訊過程進(jìn)行了封裝，以方便程序的調(diào)用。

在查找插值單元的宿主單元過程中，采用和圖2(b)類似的策略：首先搜集所有zone里面需要進(jìn)行查詢的插值單元，通過廣播并行通訊在每個進(jìn)程內(nèi)形成整體的待查詢單元。然后各個進(jìn)程在其所負(fù)責(zé)的zone內(nèi)進(jìn)行查詢，最后各個進(jìn)程再通過并行通訊來確定其插值單元的宿主區(qū)和宿主單元。

不失一般性，假設(shè)兩套規(guī)模為N(網(wǎng)格點和單元數(shù))的block進(jìn)行重疊裝配，則在未分區(qū)串行裝配的情況下，裝配所需的計算量約為2N個網(wǎng)格點在N個單元內(nèi)查詢所需的時間。在采用ADT快速查詢算法的前提下，裝配所需計算量約為2N(log2N)(一個點在規(guī)模為N的block內(nèi)的查詢計算量約為log2N)。

圖3 分區(qū)并行情況下物理邊界推進(jìn)流程Fig.3 Procedure of boundary advancing in parallel environment

在本文的分區(qū)并行策略下，假設(shè)每一個block被均勻分成了M個網(wǎng)格區(qū)(zone)，則兩套規(guī)模為N的block之間的查詢轉(zhuǎn)變?yōu)?M套規(guī)模為N/M的網(wǎng)格間的查詢。顯然屬于同一個block的zone之間不需要進(jìn)行查詢，因此對于每一個zone而言，其查詢計算量為Nlog2(N/M)。

在實際查詢過程中，可以引入一些預(yù)判方法加速查詢過程：① 計算出包圍每一個zone的最小盒子，根據(jù)該盒子信息可以對兩個zone是否存在重疊區(qū)域進(jìn)行快速的預(yù)判斷，如果兩個盒子不重疊，則兩個zone之間不進(jìn)行查詢；② 查詢一個點在某zone內(nèi)的宿主單元時，先判斷其是否落在該最小盒子之內(nèi)，否則不進(jìn)行查詢。加入上述預(yù)判斷之后，對于某zone而言僅需要對其他一部分zone內(nèi)的一部分點進(jìn)行判斷即可。假設(shè)網(wǎng)格分布均勻，分區(qū)大小均勻，在完全理想的情況下，處于非重疊區(qū)域內(nèi)的zone其查詢計算量為0，因為沒有點落在其盒子內(nèi)；而處于重疊區(qū)內(nèi)的zone，大約僅需要查詢N/M個點(落在其最小盒子之內(nèi)的其他block的網(wǎng)格點數(shù))。因此，在完全理想的情況下，單個zone查詢的計算量最多為N/M(log2(N/M))，這也決定了并行情況下進(jìn)行網(wǎng)格裝配所需的時間。

在實際問題中，網(wǎng)格單元大小、網(wǎng)格分區(qū)大小、分區(qū)形狀肯定不均勻。例如遠(yuǎn)場區(qū)域的計算網(wǎng)格尺度必然較大，此處的分區(qū)區(qū)域也較大；而靠近壁面的網(wǎng)格尺度較小，分區(qū)區(qū)域也較小。因此，對于某zone而言，落在其最小盒子之內(nèi)的其他block的點數(shù)肯定不可能為理想狀態(tài)下的N/M。假設(shè)實際情況下落在某zone的最小盒子內(nèi)的點數(shù)為k(N/M)，查詢的計算量變?yōu)閗(N/M)log2(N/M),k∈(1,M),k和計算網(wǎng)格的均勻度、分區(qū)數(shù)、分區(qū)的均勻度相關(guān)，理想情況下k為最小值1。而在實際情況中，由于在計算域遠(yuǎn)場存在較大的網(wǎng)格分區(qū)，其最小盒子甚至有可能覆蓋了其他整個block，因此k的最大值有可能為M。

所以，對于實際問題的多區(qū)并行裝配，其裝配時間是和分區(qū)數(shù)、網(wǎng)格均勻度以及分區(qū)均勻度密切相關(guān)的。為了減少針對實際問題的裝配時間，最有效的辦法是增加分區(qū)數(shù)，這樣一方面會減少每個zone內(nèi)用于查詢的網(wǎng)格單元規(guī)模，另一方面則會顯著減少其“最小盒子”的大小，減少實際進(jìn)行查詢的點數(shù)。

除此之外，本文還采用了另外一種加速策略：對重疊區(qū)進(jìn)行預(yù)判斷。首先求出覆蓋每一個block的最小盒子，然后對于某一個block，判斷其中的點是否落在其他的block內(nèi)，如果是則標(biāo)注為重疊區(qū)域點，否則標(biāo)注為非重疊區(qū)域點。對于非重疊區(qū)域的點，不參與建立整體查詢點集，不參與宿主單元搜索，同時在計算zone的最小盒子時也不予考慮；對于非重疊區(qū)域的單元，則不參與建立ADT。經(jīng)過上述重疊/非重疊區(qū)的預(yù)判斷，可以有效提高網(wǎng)格裝配效率，原因在于：① 將遠(yuǎn)場非重疊區(qū)的zone進(jìn)行了進(jìn)一步細(xì)分，減少了其中的有效點，從而縮小了其“最小盒子”，也就減少了實際進(jìn)入該zone內(nèi)進(jìn)行查詢的點數(shù)；② 處于非重疊區(qū)內(nèi)的點不參與整體查詢點集的建立，因此減少了并行通訊量，減少了每個進(jìn)程中的內(nèi)存開銷；③ 減少了非重疊區(qū)域內(nèi)zone的有效單元數(shù)，也就減小了用于查詢的ADT的規(guī)模。

在網(wǎng)格生成時著重考慮網(wǎng)格的均勻性，并通過上述這些優(yōu)化措施，可以有效減少k值的大小，從而減少并行裝配的計算時間。

3 非結(jié)構(gòu)重疊網(wǎng)格裝配實例

3.1 二維裝配實例

圖4 二維非結(jié)構(gòu)重疊網(wǎng)格裝配實例(3個圓柱)Fig.4 Example of assembling of 2D unstructured overset grids (three cylinders)

圖4所示為3個二維圓柱網(wǎng)格的裝配結(jié)果。圓柱直徑為1，背景網(wǎng)格的壁面距離設(shè)置為1，每個網(wǎng)格塊被分成為10個區(qū)。這種隱式裝配方法具有較好的重疊網(wǎng)格裝配效果，插值邊界都盡量遠(yuǎn)離壁面附近的大梯度流動區(qū)域，使插值單元和貢獻(xiàn)單元的大小更為匹配，有助于減少插值誤差并增加數(shù)值計算的穩(wěn)定性。圖5(a)為7301兩段翼外形的重疊網(wǎng)格裝配結(jié)果。共有兩個網(wǎng)格塊，每個網(wǎng)格塊分成8個區(qū)。圖5(b)是局部放大圖，圖中綠色的網(wǎng)格單元為插值單元。

圖5 二維非結(jié)構(gòu)重疊網(wǎng)格裝配實例(7301兩段翼形)Fig.5 Example of assembling of 2D unstructured overset grids (7301 airfoil)

3.2 三維裝配實例

圖6 三維非結(jié)構(gòu)重疊網(wǎng)格裝配算例 (5個圓球)Fig.6 Example of assembling of 3D unstructured overset grids (5 spheres)

圖6為三維情況下5個圓球的裝配結(jié)果。每個圓球網(wǎng)格單元數(shù)為120萬，背景網(wǎng)格單元數(shù)為448萬，總體網(wǎng)格單元總數(shù)為1 048萬。采用該算例對并行網(wǎng)格裝配的效率進(jìn)行了測試。計算機(jī)集群采用國產(chǎn)飛騰CPU，主頻為1.0～1.5 GHz。首先考慮分區(qū)數(shù)和進(jìn)程數(shù)相同的情況，表1給出了不同進(jìn)程情況下的裝配時間對比。隨著進(jìn)程數(shù)以及分區(qū)數(shù)的增加，查詢節(jié)點以及查詢插值單元所需的時間單調(diào)減少，但是當(dāng)進(jìn)程數(shù)增加以后，各個進(jìn)程之間需要通過廣播通訊建立整體的網(wǎng)格點集，搜集計算點所需的時間會有所增加。在1 152個分區(qū)和進(jìn)程情況下，完成裝配過程所需時間為4.78 s。表2為分區(qū)數(shù)和進(jìn)程數(shù)不同的情況，網(wǎng)格分成為1 152個分區(qū)，采用不同的進(jìn)程進(jìn)行裝配。查詢節(jié)點和單元所需時間隨著進(jìn)程數(shù)增加而單調(diào)減少，但是當(dāng)進(jìn)程數(shù)增加到一定程度之后，查詢所需時間趨于一個定值，而由于搜集計算點所需時間隨著進(jìn)程數(shù)增加而單調(diào)增加，最終導(dǎo)致并行效率降低，原因在于：① 算法中引入了一些優(yōu)化措施，通過對zone進(jìn)行預(yù)判斷，對于不在重疊區(qū)或者最小盒子與其他zone不重疊的zone，將直接跳過而不參與查詢，因此導(dǎo)致負(fù)載不平衡；② 根據(jù)第2節(jié)的效率分析，由于網(wǎng)格單元、分區(qū)區(qū)域大小的不均勻，每個zone里面查詢的計算量是不同的，計算最耗時的zone決定了并行裝配的最少時間。因此，重疊網(wǎng)格裝配算法本身是導(dǎo)致并行效率不高的主要原因，文獻(xiàn)[18,24]的并行裝配方法在進(jìn)程數(shù)增加后也會遇到相同的問題。因此，如何改進(jìn)并行裝配策略，提高其并行計算效率，仍然是一個需要進(jìn)一步研究的內(nèi)容。

在上述算例的基礎(chǔ)上，通過增加每個block的網(wǎng)格量以及圓球數(shù)量，對更大規(guī)模的計算網(wǎng)格進(jìn)行了并行裝配測試。圖7為裝配結(jié)果示意圖，網(wǎng)格單元總數(shù)為1.23億，采用了10 968個分區(qū)，采用1 024個進(jìn)程并行裝配，完成整個裝配所需時間為39.76 s，說明該方法已經(jīng)具備了大規(guī)模網(wǎng)格高效自動化并行裝配能力。

進(jìn)一步耦合求解六自由度動力學(xué)方程、非定常流動控制方程，采用機(jī)翼外掛物投放的標(biāo)準(zhǔn)算例對重疊網(wǎng)格技術(shù)以及耦合算法進(jìn)行測試。計算軟件為自主研制的通用CFD軟件平臺“HyperFLOW”[25-26]，其中的非定常計算方法見文獻(xiàn)[27-28]，這里不再贅述。圖8所示為機(jī)翼以及外掛物的初始計算網(wǎng)格，物面采用了三角形和四邊形單元離散，空間網(wǎng)格有三棱柱、金字塔、六面體、四面體等單元類型。機(jī)翼網(wǎng)格單元總數(shù)為320萬，外掛物單元總數(shù)為360萬。網(wǎng)格共分成為768個分區(qū)，采用128個進(jìn)程并行計算，在主頻為1.0～1.5 GHz的集群系統(tǒng)上，單步重疊網(wǎng)格裝配所需總時間約為7 s。

表1 裝配時間對比(分區(qū)數(shù)和進(jìn)程數(shù)相同)Table 1 Comparison of assembling time (same numbers of sub-zones and processors)

表2 裝配時間對比(分區(qū)數(shù)和進(jìn)程數(shù)不同)Table 2 Comparison of assembling time (different numbers of sub-zones and processors)

分離計算條件詳見文獻(xiàn)[29]，圖9所示為分離過程中典型時刻的動態(tài)重疊網(wǎng)格，其中綠色區(qū)域為外掛物網(wǎng)格的插值邊界；圖10為外掛物的運動學(xué)參數(shù)以及氣動力系數(shù)，圖中各變量定義如下：dx、dy、dz為質(zhì)心位移；Vx、Vy、Vz為質(zhì)心運動速度；CA、CY、Cn分別為軸向力、側(cè)向力、法向力系數(shù)；φ、θ、ψ分別滾轉(zhuǎn)角、俯仰角、偏航角；p、q、r分別為滾轉(zhuǎn)角速率、俯仰角速率、偏航角速率；Cl、Cm、Cn分別滾轉(zhuǎn)力矩、俯仰力矩、偏航力矩系數(shù)。圖中實線表示本文計算結(jié)果，離散點表示試驗結(jié)果，計算結(jié)果與試驗結(jié)果吻合較好。圖11所示為3個典型時刻的物面壓力p云圖(無量綱壓力，以來流的動壓為參考值)，激波干擾等現(xiàn)象捕捉較好。

圖7 大規(guī)模非結(jié)構(gòu)重疊網(wǎng)格裝配算例 (6個圓球)Fig.7 Example of assembling of large scale unstructured overset grids (6 spheres)

圖8 機(jī)翼以及外掛物的初始計算網(wǎng)格Fig.8 Initial computation grids over wing/store configuration

圖9 機(jī)翼及外掛物分離過程中的動態(tài)重疊網(wǎng)格Fig.9 Overset dynanic grids over wing/store configuration during store separation

圖10 分離過程中外掛物的運動學(xué)參數(shù)以及氣動力系數(shù)Fig.10 Kinematic parameters and aerodynamic coefficients during store separation

圖11 分離過程中3個典型時刻的物面壓力云圖Fig.11 Surface pressure contours during store separation at three typical time

4 結(jié) 論

建立了一種并行化的非結(jié)構(gòu)重疊網(wǎng)格隱式裝配技術(shù)。采用壁面距離作為判斷插值邊界的準(zhǔn)則，并采用物理邊界推進(jìn)的方式快速確定出活躍區(qū)域，避免了隱式裝配方法中的“孤島”問題。采用分區(qū)并行的思路實現(xiàn)了隱式裝配技術(shù)的并行化，所采用的并行計算策略有助于提高網(wǎng)格裝配效率并減少內(nèi)存消耗。效率分析以及針對實際問題的裝配結(jié)果表明，該技術(shù)能夠適用于多塊重疊網(wǎng)格，自動化程度較高，裝配效果較好，能夠推廣應(yīng)用于超大規(guī)模計算網(wǎng)格。

然而隨著進(jìn)程數(shù)的增多，查詢算法所占時間比重逐漸減少，此時并行通訊所占的比重增多并逐漸成為影響裝配效率的一個關(guān)鍵問題。本文的并行裝配方法在并行效率上仍有較大的改進(jìn)空間，下一步工作中將針對算法中影響并行效率的諸多細(xì)節(jié)問題如單元屬性通訊、節(jié)點搜集等進(jìn)行細(xì)致分析和改進(jìn)，進(jìn)一步提升該方法的適用性。