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

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

符號說明

0 引言

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

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

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

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

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

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

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和三角形內(nèi)任一點T(s,t)距離的平方為橢圓函數(shù)：

其 中，(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

最小距離問題轉(zhuǎn)換為在D內(nèi)求連續(xù)可微函數(shù)Q(s,t)的極值問題。令：

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

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

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

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

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

2 物面單元的存儲

2.1 嵌套包圍盒概念的KDT 參數(shù)定義

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

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

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

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

步驟2：創(chuàng)建一個空的KDT 備用。

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

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

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

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

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

2.3 創(chuàng)建KDT 時間復雜度分析

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

3 最小距離查詢

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

3.1 算法設計

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

2）根據(jù)1）確定的子樹，對此子樹的節(jié)點進行步驟1）的操作。

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

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

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

3.2 算法優(yōu)化

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

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

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

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

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

4 算例測試

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

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

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

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

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

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

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

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

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

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

5 結論

在本文工作中，我們基于笛卡爾網(wǎng)格發(fā)展了一種壁面距離計算方法，具備精確性、高效性以及通用性?；诒疚膬?nèi)容，主要有以下3 個方面結論：

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

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

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

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