并行化非結構重疊網格隱式裝配技術

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

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

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

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

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

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

1) 非結構網格本身占用內存較大，由于在重疊網格裝配過程中需要存儲額外的拓撲信息，因此當網格規模增加后內存消耗量會急劇增加。

2) 查詢計算量較大。盡管可以采用交替數字樹(Alternation Digital Tree，ADT)數據結構[16]或者一些基于網格拓撲信息的快速查詢方法，但是計算效率和顯式挖洞技術相比仍然有較大差距，并且查詢的計算量會隨著網格規模的增加而急劇增大。

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

隨著CFD技術的不斷發展，工程領域對CFD的期望也越來越高。在航空航天以及民用領域，越來越多的工程應用希望CFD能夠采用更大規模的計算網格，模擬更真實的復雜外形，并能夠更好地模擬分離流動、湍流流動現象。為了適應超大規模的計算網格，并行化裝配是一個必然的選擇。

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

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

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

1 非結構重疊網格裝配技術

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

圖1 非結構重疊網格裝配流程Fig.1 Assembling procedure of unstructured overset grids

1.1 判斷點的類型

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

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

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

1.2 判斷單元的類型

在判斷出點的類型之后，再根據“體-點”的網格拓撲信息判斷出單元類型。本文采用了如下判斷準則：一個單元的所有點為活躍點(1屬性)，則該單元為活躍單元(1屬性)；如果該單元的所有點為非活躍點，則該單元為非活躍單元(0屬性)；剩下的單元為插值單元(-1屬性)。

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

1.3 提高算法魯棒性的措施

為了提高非結構重疊網格應對復雜工程問題的能力，本文采用了一些提高網格裝配魯棒性和效率的措施。

1) 對插值邊界進行優化。采用非結構網格中“點-體”和“體-點”的拓撲關系拓展1～2排活躍點，適當擴大活躍區域，使得插值單元完全落在鄰居zone的活躍區域內。

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

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

4) 為了提高查詢效率，采用ADT數據結構用于宿主單元查詢。對于每一個網格塊，以圍繞每一個網格單元的最小盒子作為建立ADT樹的依據。

2 基于網格分區的并行通訊策略及裝配效率

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

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

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

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

在判斷點的類型時，需要先搜集各個zone的網格點坐標并進行廣播通訊，在每個進程內形成整體的待查詢點集，然后在本進程所負責的zone內查詢該點集的宿主單元并計算壁面距離。點在所有鄰居zone內的最小壁面距離最終通過各個進程間的并行通訊來完成。

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

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

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

圖3 分區并行情況下物理邊界推進流程Fig.3 Procedure of boundary advancing in parallel environment

在本文的分區并行策略下，假設每一個block被均勻分成了M個網格區(zone)，則兩套規模為N的block之間的查詢轉變為2M套規模為N/M的網格間的查詢。顯然屬于同一個block的zone之間不需要進行查詢，因此對于每一個zone而言，其查詢計算量為Nlog2(N/M)。

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

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

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

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

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

3 非結構重疊網格裝配實例

3.1 二維裝配實例

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

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

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

3.2 三維裝配實例

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

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

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

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

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

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

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

圖7 大規模非結構重疊網格裝配算例 (6個圓球)Fig.7 Example of assembling of large scale unstructured overset grids (6 spheres)

圖8 機翼以及外掛物的初始計算網格Fig.8 Initial computation grids over wing/store configuration

圖9 機翼及外掛物分離過程中的動態重疊網格Fig.9 Overset dynanic grids over wing/store configuration during store separation

圖10 分離過程中外掛物的運動學參數以及氣動力系數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 結 論

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

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