笛卡爾網格下精確高效的壁面距離計算方法

孟 爽，周 丹，李雪亮，畢 林，3，*

(1.中南大學 軌道交通安全教育部重點實驗室，長沙 410075；2.空氣動力學國家重點實驗室，綿陽 621000；3.中國空氣動力研究與發展中心 計算空氣動力研究所，綿陽 621000)

符號說明

0 引言

笛卡爾網格可以不依賴物面網格直接生成，因具有自動化程度高、復雜外形適應性好、非定常/多尺度流動結構捕捉能力強等優勢而廣受關注[1]。

壁面距離的精確高效計算對實現笛卡爾網格方法具有重要意義：一方面，由于笛卡爾網格的非貼體特性，通常采用虛擬單元法處理物面邊界[2-6]，需要根據物面確定虛擬單元的壁面距離，尋找參考點并插值獲得其物理量，壁面距離精確計算對虛擬單元物面處理精準度有較大影響；另一方面，對于非定常流動或運動物體，笛卡爾網格會根據流場的變化或物體的運動進行自適應[7-8]，自適應后的網格壁面距離需重新計算，壁面距離計算效率成為制約流動計算效率的關鍵因素之一。

目前，壁面距離的計算方式主要分為求解偏微分方程[9-12]和采用幾何方法直接計算[13-14]兩類。第一類方法，將最小距離轉換為關于波傳播問題的偏微分方程數值求解[15]，需要額外的計算花費，壁面距離計算精度受數值離散精度的影響，且對于復雜外形適應性較差。第二類方法是根據空間點與物面離散網格的幾何關系計算壁面距離。通常為簡單起見，在物面網格足夠密的情況下，求解空間網格點的壁面距離時，一般采用空間點到物面網格中心的距離代替最小距離。趙鐘等[16]發現，這種近似距離會引起計算誤差，不僅影響計算結果的精度，而且對計算穩定性造成影響。為提高計算精度，有學者發展了根據空間點關于物面網格所在平面投影點與物面網格的幾何關系計算壁面距離的3D 投影法[16]、通過坐標變換轉化為二維問題的2D 法[17]等，但這些方法涉及大量向量計算和關系判斷，計算量較大。

在幾何方法求解壁面距離時，由于物面網格數目較大，計算效率的核心問題是如何快速定位到最短距離相關的物面網格，物面網格數據結構是關鍵。

最常見的方法是將物面網格順序存入數組，采用遍歷法搜索物面點得到壁面距離。網格規模較小時，計算量尚在承受范圍之內。但是對于三維復雜幾何外形的流場計算，物面網格數量可達到O(105)～O(106)量級，空間網格點可達到O(108)～O(109)量級，采用遍歷法計算量達到O(1013)～O(1015)量級。此種規模的計算量對整個計算周期的影響是巨大的。Boger[18]提出采用ADT (alternating digital tree)二叉樹數據結構存儲物面網格，計算效率大幅提高，成為目前最常用的物面網格數據存儲方式。但對于復雜外形，ADT 平衡性下降，計算效率仍然差強人意。郭中州等[14]采用KDT (K-dimensional tree)存儲物面網格，解決了ADT 平衡性問題，計算效率進一步提升，但對于離物面較遠的空間點，在運用KDT 最近鄰搜索算法時，超球面與分割面位置關系判斷失準，需要反復回溯，計算量大。

本文針對上述問題，引入三角形參數化方法，將空間點到三角形物面離散網格的最小距離問題轉換為約束條件下一維極值問題；同時發展嵌套包圍盒概念的KDT 物面網格數據存儲結構，優化KDT 最近鄰搜索算法中距物面較遠數據點回溯過程，實現最小距離對應的三角形的快速定位。

1 點到三角形最小距離精確算法

點到空間三角形最小距離的問題[19]可描述為點P和三角形T(s,t)=B+sE0+tE1之間的最小距離，其中(s,t)∈D={(s,t):s≥0,t≥0,s+t≤1}，B為三角形的一個頂點，E0和E1分別為此頂點對應的三角形兩條邊。點P和三角形內任一點T(s,t)距離的平方為橢圓函數：

其 中，(s,t)∈D，a=E0·E0，b=E0·E1，c=E1·E1，d=E0·(B-P)，e=E1·(B-P)，f=(B-P)·(B-P)。點P和三角形之間的關系如圖1 所示。

圖1 點P 和三角形之間的關系Fig.1 The relationship between P and the triangle

最小距離問題轉換為在D內求連續可微函數Q(s,t)的極值問題。令：

圖2 用三角域D 劃分s-t 平面以及不同等高線Fig.2 Partition of a s-t plane by triangle domain D and distance contours

圖3 以區域③為中心的等高線與三角形接觸的兩種情況Fig.3 Two contour levels with centers in region ③touching the triangle

在等值線上的任何一點，-Q的方向均朝向橢圓內部。為簡單起見，我們根據-Q沿頂點(0,1)的兩條邊的方向導數的符號進行判斷：

1）對于邊s+t=1，頂點(0,1)處的方向導數為如果值為負，最小值出現在邊s+t=1上。最小值的情形對應圖3 紅色曲線，處理方式和s∈[0,1]在區域②類似。

2）對于邊s=0，頂 點(0,1)處的方向導數為(0,-1)·?Q(0,1)。如果值為負，最小值出現在邊s=0上。最小值的情形對應圖3 綠色曲線，處理方式和 (,)在區域②類似。

2 物面單元的存儲

2.1 嵌套包圍盒概念的KDT 參數定義

與文獻[14]中只存儲物面網格中心點不同，本文發展了基于嵌套包圍盒概念的KDT，不僅按逆時針順序存儲物面網格頂點信息（方便笛卡爾網格物面處理），還同時包含KDT 節點對應所有物面網格的最小包圍盒，KDT 節點從上到下形成嵌套包圍盒形式。嵌套包圍盒的引入保證KDT 每個節點包含當前節點對應三角形的所有信息，且在最小距離查詢時，僅需要通過與剖分面（剖分面的大小由包圍盒確定）坐標對比，就可以快速排除與目標點距離大于當前最小距離的三角形，效率大幅提升。

2.2 嵌套包圍盒概念的KDT 建立步驟

步驟1：在正式創建KDT 之前，首先需要確定目標數據集G中三角形元素的數量N。然后，需要根據這N個元素確定G的最小包圍盒，包圍盒的兩個頂點分為 (Gimin,Gimax)，其中i=1，2，3，其定義如圖4 所示。

圖4 幾何外形包圍盒Fig.4 Geometry bounding box

步驟2：創建一個空的KDT 備用。

步驟3：若N為1，則將此元素插入到當前KDT的節點中去，同時，將此節點的左右子樹均置空。

步驟4：若N為2，則將第一個元素插入到當前KDT 的節點中去，同時規定第一維為劃分維。此時，只有第二個元素尚未被插入到KDT 中。將第二個元素中心點與第一個元素的中心點在第一維的值進行比較。若第二個元素大于第一個元素，則對當前節點的右子樹進行步驟3 的操作，同時當前節點的左子樹置空；否則，對當前節點的左子樹進行步驟3 的操作，同時當前節點的右子樹置空。此步驟中d加1。

步驟5：若N大于2，則首先確定劃分維。為了保證KDT 的平衡，劃分維的確定是計算N個三角形的中心點Tcenter在各個維度的方差，然后找出方差最大的維度即為劃分維。然后將此N個元素中心點在方差最大的維度上按從小到大的順序進行排序。將位于中間的元素插入到當前KDT 的節點中去。小于此元素的將全部位于當前節點的左子樹，大于此元素的將全部位于當前節點的右子樹。對當前節點左子樹和右子樹的元素遞歸執行步驟3 或步驟4 或步驟5，直至所有元素全部插入到KDT 中。

值得注意的是，在向KDT 中插入節點的同時，節點對應的包圍盒也在進行同步劃分。例如KDT 根節點對應的包圍盒即為數據集G的包圍盒 (Gimin,Gimax)。假設節點K對應的包圍盒為(Kimin,Kimax)，則K的左子節點KL的包圍盒(KLimin,KLimax)為：

由于節點包圍盒表示此空間區域內所有三角形元素集合的最小包圍盒，而在進行空間區域剖分時的依據是三角形中心點確定的平面，所以節點包圍盒與包圍盒之間，剖分面與剖分面之間會存在相交的部分。這種相交只會出現在三維及以上的情況。

2.3 創建KDT 時間復雜度分析

由于KDT 是一個二叉樹，所以對于給定的N個元素建成一個深度為d的最優KDT 有：

3 最小距離查詢

計算空間某點到物面的最小距離需要兩個步驟。第一，得到此點到物面最小距離的點所在的物面三角形；第二，采用點到三角形最小距離精確算法得到此點到物面的最小距離。本節基于KDT 最近鄰算法對第一步進行設計。

3.1 算法設計

1）首先從KDT 的根節點開始查詢，根據目標點與根節點剖分面的位置關系，確定與目標點T潛在的最近點是在根節點的左子樹還是右子樹。

2）根據1）確定的子樹，對此子樹的節點進行步驟1）的操作。

3）當遞歸訪問至葉子節點時，記錄T與當前節點之間的距離作為當前最近距離lmin，當前節點作為與目標點T最近的節點記為Knearest。

4）進行回溯操作。回溯操作從當前最近的點Knearest開始向上進一步查找。構造以T為球心，lmin為半徑的球。若球面與Knearest父節點的剖分面相交，則查找Knearest的兄弟節點即Knearest父節點的另外子樹的節點，同時計算T與Knearest父節點之間的距離，若小于當前最小距離，則Knearest和lmin均進行更新；若球面與剖分面不相交，則說明Knearest的兄弟節點中不存在比Knearest與T更近的節點。

5）執行步驟4）直至根節點。得到的Knearest即為與T最近的點，lmin為Knearest與T的最小距離。至此最小距離查詢模塊結束。

3.2 算法優化

在實際流場計算中，存在空間網格點與物面距離較遠的情況，此時采用上述最小距離查詢方法將面臨回溯操作效率低下的問題。為了直觀說明問題，本小節以二維數據點為例進行說明，給定隨機分布的17 個點，對這些點創建KDT，如圖5 所示。給定目標單元T，正常查找得到的最鄰近單元是G。在回溯操作中，以T為圓心，T和G之間的距離lTG為半徑的圓與根節點J的剖分軸m相交。根據上述回溯步驟，此時需在J的左子樹進行查詢。然而，實際情況是圓與J的左子樹并不相交。因此，后續的回溯操作均為無效回溯，這將大大浪費計算資源。為此，本文在回溯操作時，首先判斷以T為圓心，T和G之間的距離lTG為半徑的圓是否與線段ab相交（ab由剖分軸和當前包圍盒確定）。若相交，則進行正常的回溯操作；若不相交，則直接排除對應節點另外子樹存在最小距離的可能性。采用改進前后的回溯法進行的回溯操作如圖5（b～c）中綠色箭頭所示。從圖中可以看出，改進后的回溯方法可極大減少不必要的回溯過程，進而提高最小距離查詢效率。

圖5 圓與剖分軸相交且回溯無效的情況Fig.5 Sketch of an optimized backtracking process when the circle intersects the split axis

圖6 為球、導彈和DPW6 三種幾何外形及空間笛卡爾網格。為了定量對比改進前后回溯算法的效率，本文對這三種不同幾何外形進行了壁面距離查詢測試。對幾何外形進行壁面距離查找時，存儲所有空間網格點的壁面距離，同時記錄得到壁面距離所需的查找次數，然后將查找次數的總和除以空間網格點數得到平均查找次數，結果如圖7 所示。從圖中可以看出，針對不同外形，本文改進的回溯方法均可有效減小查找次數。針對同一幾何外形，物面三角形數量越多，查找次數減少的越明顯。

圖6 三種幾何外形Fig.6 Three geometries

圖7 求最小距離所需查找次數Fig.7 Number of queries for finding the minimum distance

4 算例測試

本小節對本文發展的壁面距離算法的準確性和效率進行了測試。測試的幾何外形仍為球、導彈和DPW6。為了保證測試結果的合理性，本文所有工況均在AMD EPYC 7702 的處理器上測試，測試平臺為Window 10 系統及Visual Studio 2012 版本。每個工況測試10 次，統計平均獲得程序運行的時間。

首先對壁面距離結果的準確性進行對比測試。對壁面距離的兩種計算方式進行對比，第一種是近似壁面距離，即空間點到物面三角形中心的最小距離，第二種是精確壁面距離，即采用本文發展的算法計算得到的最小距離。計算結果如圖8 所示，圖8（a）為近似壁面距離計算的結果，圖8（b）為精確壁面距離計算所得結果。從圖8（a）中可以明顯看出，球的近似壁面距離云圖及等值線呈現波浪狀，與實際物理特征明顯不符。采用本文方法計算得出的精確壁面距離等值線光順平滑，符合物理特征。對于幾何外形較為復雜的DPW6，在遠離壁面區域，兩種計算方法得出的壁面距離結果接近；而在近壁區域，采用近似壁面距離計算方法得到的壁面距離為250 的等值線與實際幾何特征誤差較大，而精確計算結果的近壁面等值線可較好反映近壁物面特征。

圖8 不同幾何外形近似壁面距離與精確壁面計算結果云圖及等值線Fig.8 Contours of approximate and accurate calculation results for wall distances of different geometric shapes

此外，以球為研究對象，我們將本文方法計算得到的壁面距離與解析解進行對比。結果表明本文方法的計算結果與解析解差異在百萬分之一以內。因此，本文方法得到的壁面距離精確性較高。

在計算效率測試方面，對遍歷法+精確距離、ADT+精確距離、KDT+精確距離，KDT+近似距離4 種方法的計算效率進行對比。物面三角形數量、空間網格數量以及CPU 運行時間如表1 所示。從表中可以看出，采用KDT 方法的計算效率相比遍歷法和ADT 方法均有量級的提升。尤其對于復雜的DPW6 來說，采用遍歷法的時間為3 616.87 s，KDT 只需2.92 s，效率提升超過1 000 倍。而對于簡單的球來說，效率提升僅不到300 倍。因此，壁面距離查詢效率的提升程度與幾何外形密切相關。

表1 壁面距離查詢所需時間Table 1 Wall-distance querying time

精確距離的計算要比近似距離的計算多調用一個點到三角形距離的函數，此函數在回溯過程中頻繁調用會增加計算耗時。從定量結果來看，精確距離計算耗時為近似距離計算耗時的2 倍左右，但這對整個計算周期的耗時影響并不明顯。因此，通過增加少量計算時間來獲取計算精度收益是值得的。

為了進一步測試本文發展的方法在大規模網格計算時的效率，我們選用DPW6 外形進行物面網格離散并生成空間笛卡爾網格。物面三角形網格數量跨越萬、十萬以及百萬3 個量級；空間笛卡爾網格數量跨越千萬、億以及十億3 個量級。效率測試結果如表2 所示，結果表明，對于同一種物面離散方式，CPU運行時間與空間網格數量大致呈正比；對于同一數量空間網格（由于物面網格數量不同，加密相同層級生成的空間笛卡爾網格數量有略微差異），物面網格數量從30 021 增加到1 417 181 時，CPU 運行時間僅增加為原來的1～2 倍，這再次說明了本文基于KDT 的回溯方法的高效性。

表2 DPW6 壁面距離計算效率Table 2 Wall-distance calculation time for DPW6

為了與文獻[16]中計算效率進行對比，本文選用JSM（JAXA standard model）標模進行壁面距離查詢測試，網格規模為33 億。本文方法（無并行）進行壁面距離查詢所需時間大致為0.51 h，文獻[16]中采用并行方法對進行壁面距離查詢所需時間為0.42 h。因此，除去網格類型、測試環境帶來的誤差，總體來說，本文方法在大規模網格下仍表現出較高效率。

5 結論

在本文工作中，我們基于笛卡爾網格發展了一種壁面距離計算方法，具備精確性、高效性以及通用性。基于本文內容，主要有以下3 個方面結論：

1）在提高壁面距離計算的準確性方面，我們將物面三角形進行參數化（s和t）表示，根據s和t的取值范圍將s-t平面分為7 個區域，對不同區域進行分情況討論，從而將求最小距離問題巧妙地轉化為求約束條件下一維極值問題，計算量大幅減少。

2）在提高計算效率方面，采用平衡的KDT 存儲物面三角形所有頂點信息并建立相應的嵌套包圍盒，使得在最小距離查詢時盡可能快速排除在目標范圍之外的單元；同時優化了當前超球面與剖分平面的位置判斷法則，解決了在距離物面較遠處因無效回溯次數過多導致查詢效率下降的問題。

3）此外，本文的壁面距離查詢模塊是一個單獨的模塊，在計算壁面距離時，僅需知道空間網格點的位置信息以及物面三角形信息。因此，本文發展的壁面距離計算方法不僅對笛卡爾網格，對結構網格、非結構網格以及重疊網格等也可較好兼容，可方便程序的移植。

本文采用3 個三維幾何外形對本文方法進行測試，結果表明計算得到的壁面距離與解析值的誤差在百萬分之一以內，有較高的精確性；十億量級網格規模下的單核計算效率與已有文獻[16]的并行計算效率相當。雖然在大規模網格測試中表現出較高的計算效率，但仍有提升空間。下一步我們計劃添加OpenMP/MPI 并行方法進一步提升計算效率，以滿足下一代超大規模CFD 湍流模擬的需求。