趙 鐘, 何 磊, 張 健, 徐慶新, 張來平

(中國空氣動力研究與發(fā)展中心 計算空氣動力研究所, 四川 綿陽 621000)

0 引 言

隨著計算機科學技術的不斷發(fā)展，計算流體力學(Computational Fluid Dynamics，CFD)已經(jīng)成為與風洞試驗、飛行試驗互為補充的三大空氣動力學研究手段之一，被廣泛應用于航空航天飛行器設計過程之中。

長期以來，作為經(jīng)典物理中的世紀性難題，湍流被認為是一個重大的基礎科學問題，其流動機理至今仍是一個未解之謎[1]。利用CFD模擬層流流動已經(jīng)較為成熟，但是面對實際流動中的湍流問題，由于其物理本身的高度復雜性，只能做簡化處理。當前，湍流模擬主要有三類方法：雷諾平均N-S方程(RANS)、大渦模擬(LES)、直接數(shù)值模擬(DNS)。DNS能直接進行流動的數(shù)值模擬，但是由于需要極其龐大的計算資源，在現(xiàn)在和可預見的未來都難以用于模擬真實外形；LES可精確模擬流動分離，但是其計算資源仍然耗費巨大。RANS是目前唯一能夠直接用于實際工程計算的方法，且對常規(guī)流動問題的模擬能取得令人滿意的結果，是當前CFD工程應用中大量采用的方法。

目前的工程應用中大部分RANS數(shù)值模擬都是采用SA和k-ωSST兩類湍流模型，這兩類模型方程中均涉及壁面距離，即空間點到物體表面的最近距離。在計算前的預處理過程中，需要首先計算壁面距離。壁面距離計算方法主要有兩類：一類是通過幾何方法直接計算，另一類是通過求解偏微分方程的形式求解[2-4]。對于求解偏微分方程的方法，在簡單幾何外形上速度較快，但是在真實復雜外形上難以應用，主要原因有：(1)CFD代碼開發(fā)者除了要開發(fā)流動方程的求解器外，還要開發(fā)一個額外的求解器，帶來額外的工作負擔；(2)只要是偏微分方程求解類問題，都存在穩(wěn)定性問題，而在復雜外形上常常難以收斂，容易算發(fā)散；(3)偏微分方程求解得到的是近似解，而壁面距離的精確性直接影響到CFD求解器的穩(wěn)定性、收斂性和計算結果的精度(下文將詳細闡述)。因此，目前實際應用中基本上都是用幾何方法直接搜索計算。

直接搜索法，即對于空間中的任意一個點，搜索壁面上的離散單元(如三角形、四邊形)離其最近的距離。搜索過程中，需要一一比較與壁面上的每個單元的距離以找到最近點。對于N個空間點、M個壁面面元，搜索比較次數(shù)為N×M次，對于二維問題來說，這種搜索方式尚且能忍受，但是對于三維真實航空航天飛行器外形，N值一般為千萬、上億量級，M值一般為百萬、千萬量級，在復雜應用中甚至到達數(shù)十億，直接搜索將耗費大量時間，是整個CFD模擬所無法承受的。針對上述效率問題，研究者提出了幾種高效的計算方法。

NSMB[5]、Spalart[6]等介紹了一種沿網(wǎng)格線找到離其最近的壁面點的方法，但是這種方法只適應于結構網(wǎng)格計算，無法應用到非結構網(wǎng)格中，算法不具備通用性。王剛等[7]提出了一種逐層推進的方法，基本原理是從壁面開始，向流場逐層推進搜索最近距離，該方法適用于結構網(wǎng)格和非結構網(wǎng)格，但是實現(xiàn)過程比較繁瑣，編程困難，而且難以并行化。

趙慧勇等[8]介紹了一種基于二分法將壁面循環(huán)劃分為若干個盒子的方法，即循環(huán)盒子法，首先計算離空間點最近的若干個盒子，再在這些盒子中搜索最近的壁面點。循環(huán)盒子法的本質是將整個壁面用八叉樹劃分，每個方盒子相當于一個八叉樹節(jié)點，從而實現(xiàn)了串行直接搜索方法的大幅加速。這種方法適用于任意網(wǎng)格，但是這種方法仍存在兩方面的不足：(1)算法難以并行化;(2)當網(wǎng)格量大幅增加后，效率將大幅下降。這是因為，如果盒子數(shù)過少，則每個盒子中的壁面面元數(shù)將大幅增加，每個盒子中的搜索相當于直接搜索，導致效率大幅下降；相反，如果盒子數(shù)過多，在計算離空間點最近的盒子過程中需要首先進行盒子距離排序，這一排序過程將耗費大量時間，甚至與盒子內單元直接搜索時間相當。上述兩個問題導致循環(huán)盒子法的擴展性受限。楊永國等[9]、李廣寧等[4]對這一方法中的盒子劃分方法、搜索方法進行了改進提高，但算法上并無本質區(qū)別。

Boger[10]采用了一種高效的ADT數(shù)據(jù)結構，利用三角不等式，設計了一種壁面距離的快速計算方法，可以有效利用空間中點的坐標信息，實現(xiàn)對直接搜索法的大幅加速。但是正如李廣寧[4]所述，這種方法存在兩個缺點，一是搜索半徑難以確定，二是不平衡的ADT樹導致外流計算問題中壁面距離計算效率大幅下降。

上述方法在幾何直接搜索法上進行了不同的改進，相對于原始方法而言，計算效率有了大幅的提升，對于一般工程應用中的百萬、千萬量級計算網(wǎng)格是適用的。但是，隨著設計人員對模擬分辨率要求的不斷提高、以及如DES等CFD計算方法不斷進入工程應用，目前工程應用中的計算網(wǎng)格量不斷增加，航空航天飛行器設計型號任務的網(wǎng)格量接近上億量級，甚至更多。另外，以“神威·太湖之光”、“天河二號”等為代表的高性能計算機硬件不斷發(fā)展，目前CFD軟件也在不斷進行硬件適配，以適應于超大規(guī)模并行計算，計算網(wǎng)格量可達十億、百億量級。使用大規(guī)模網(wǎng)格帶來的問題是，上述壁面距離計算方法已完全無法滿足大規(guī)模并行計算的預處理需求。

針對上述問題，本文提出了一種基于ADT數(shù)據(jù)結構的湍流模擬壁面距離并行計算方法，適用于十億、百億量級網(wǎng)格的大規(guī)模數(shù)值模擬問題。

1 壁面距離精確計算的重要性

1.1 湍流模型方程基礎

CFD工程應用中常用的湍流模型方程有一方程SA模型和兩方程k-ωSST模型。這里以SA模型為例，說明湍流壁面距離對計算結果的影響。SA湍流模型方程為[11-12]：

(1)

上式寫成非守恒物質導數(shù)形式為：

(2)

(3)

這里，下標wall、farfield分別指代固體壁面(即物面)和遠場。式(3)求在壁面上的湍流變量為0。SA方程是一個精巧的數(shù)學模型，通過式(2)中右端三項相互調節(jié)，使整個流場內的湍流黏性達到平衡狀態(tài)。尤其是，為了與式(3)“壁面上湍流變量為0”這一邊界條件相匹配，該方程通過右端三項驅動，使得空間湍流變量沿壁面法向，從空間到壁面逐漸衰減至0，在緊鄰壁面的一層網(wǎng)格的湍流變量值是一個接近0的小值。此外，定常狀態(tài)計算收斂時，要求湍流變量不再隨時間變化，即要求式(1)右端幾項之和為0。在未收斂前，右端幾項之和往往不等于0，將其定義為“殘差”。

1.2 壁面距離精確性的影響

三維情況下，壁面被離散為三角形或四邊形網(wǎng)格集，空間點的精確壁面距離是點到三角形面片或四邊形面片集的最近距離。一般情況下，為了編程方便、提高計算速度，會采用近似方法計算壁面距離，例如上文提到的求解偏微分方程的方法，或者幾何上的近似計算，即將點到三角形/四邊形面心的最近距離來代替壁面距離，從而帶來計算誤差。

壁面距離計算的精確性對湍流模擬的健壯性、計算結果精度有重要的影響。這里以二維湍流平板算例為例說明。圖1(a)是算例邊界條件，下側x>0為壁面、x<0為對稱面，左右側為壓力入/出口，上側為遠場。圖1(b)是計算網(wǎng)格，流向和法向網(wǎng)格維度為35×25，為了說明壁面距離精確性，將緊鄰壁面的第一層網(wǎng)格劃分為各向異性三角形。要注意的是，為了觀察方便，圖1(b)中的x和y方向采用不同的長寬比可視化，真實的第一層三角形網(wǎng)格的長寬比為1×107量級的高度拉伸的各向異性三角形。

(a)邊界條件[12]

(b)計算網(wǎng)格

對于圖1(b)中的每一個三角形或四邊形計算網(wǎng)格，精確壁面距離是其y坐標。這里，首先采用精確的壁面距離將該算例計算到收斂狀態(tài)，使得湍流方程(1)的右端項即殘差接近于0，然后在第一層的三角形網(wǎng)格中人為引入壁面距離誤差，觀察對湍流方程的影響。圖2是給定第一層三角形網(wǎng)格不同的誤差后，對右端各項以及殘差的影響，可以看到，當誤差增加后，與壁面距離相關的破壞項迅速變化，殘差從0迅速增加到1×103量級。從局部放大圖看，即使誤差有小幅增加，也將使得殘差被快速放大。

(a)最大4倍誤差

(b)最大50%誤差

SA方程中通過壁面距離來判斷空間網(wǎng)格離壁面的相對遠近，第一層網(wǎng)格的壁面距離本應該相對小，殘差應該是接近于0的小值，其壁面距離值被人為增大后，使得其本來具備的“近壁面”屬性喪失而被“誤認為”是“遠壁面”單元，從而使得湍流模型方程中的壁面破壞效應迅速喪失，將嚴重影響CFD計算的收斂性，甚至導致計算發(fā)散。為了解決這一問題，在一些CFD代碼中，通過強行將壁面上的湍流方程通量修正為0的方式來提高計算魯棒性，然而，這樣可能會帶來另一個嚴重的問題：計算無法完全收斂。圖3中顯示了人為修正壁面通量對平板湍流模擬的影響，圖3(a)是x=1.05處邊界層內沿y向的湍流方程右端項分布。如果從整體看分布圖沒有任何異常，但是，當僅觀察壁面附近的幾層網(wǎng)格(圖3b)時，可以看到，即使對已經(jīng)收斂后的流場進行殘差修正，也會使得壁面附近的1～2層網(wǎng)格上的湍流項失真，湍流場在迭代計算過程中不斷呈現(xiàn)“收斂-修正-收斂-修正”的循環(huán)，使得湍流方程求解無法完全收斂。而當去掉殘差修正后(圖3c)，湍流方程收斂時其右端三項之和為0，符合SA模型方程的控制機理。

(a)邊界層區(qū)域

(b)局部有通量修正

(c)局部無通量修正

總的來說，壁面距離的精確與否，將直接影響計算的穩(wěn)定性和計算精度。對于不精確的壁面距離，如果在湍流方程求解過程中采用類似殘差修正等手段，雖然能保證計算順利進行，但是會影響計算收斂性，從而影響計算結果的準確性。

2 計算幾何基礎

一般情況下，真實的壁面外形在非間斷處是光滑連續(xù)的，在CAE或者計算機輔助幾何設計里，一般采用參數(shù)曲線、曲面描述。在CFD計算中生成計算網(wǎng)格時，會將這些連續(xù)的參數(shù)化曲面離散為三角形/四邊形面片，對于壁面距離計算問題，將轉化為計算空間點到三角形/四邊形面片最近距離的求解問題。一些CFD代碼為了降低實現(xiàn)難度，提高計算效率，往往會采用近似的幾何方法來計算壁面距離，最常用的是直接用空間單元中心到壁面面元面心的距離來求得最近壁面距離。本文嘗試采用較為精確的計算方法。幾何方法計算壁面距離的基礎是一些基本的幾何代數(shù)運算，這里先簡要介紹本文所述方法中采用的計算幾何基礎[13]。

壁面距離計算的通用情況是計算空間點到凸N邊形的距離，為了提高計算精度，將凸N邊形的每條邊和多邊形中心連接，形成N個三角形，最終轉換為計算點到N個三角形的距離。分兩步計算點到三角形的距離：首先將空間點投影到三角形所在平面，然后判斷投影點是否在三角形內部。

圖4是點投影示意圖，將空間點P投影到△V0V1V2所在的平面P得到投影點PI。圖中，投影距離d=r·n，投影點PI=P-d·n，即將空間點沿逆法向方向投影到平面P。

圖4 點投影到三角形

圖5中，采用參數(shù)化坐標方法判斷投影點PI是否在三角形內[14]。參數(shù)化的平面P可以表示為式(4)。

V(s,t)=V0+s(V1-V0)+t(V2-V0)

=V0+su+tv

(4)

上式中，u=V1-V0和v=V2-V0分別是三角形的兩條邊向量。s和t是參數(shù)化坐標，對于任意點PI=V(s,t)，若位于三角形中，必滿足下式：

s≥0,t≥0,s+t≤1.0

(5)

因此，如要判斷投影點PI是否在三角形內，只需要找到其參數(shù)化坐標s和t即可。設w=V1-V0，PI的參數(shù)化坐標s和t可通過下式計算：

(6)

通過式(6)計算得到參數(shù)化坐標，然后通過式(5)判斷投影點PI是否在三角形內，若在三角形內，則投影距離d是空間點P投影到△V0V1V2的距離。特別地，若投影點不在所有物面三角形內部，則計算空間單元中心到物面點/三角形中心的最近距離來作為壁面距離，這種情況極少。

圖5 判斷點是否在三角形內部

3 大規(guī)模壁面距離快速計算

本文的大規(guī)模壁面距離計算方法幾個基本要素包括并行計算方法、基于ADT數(shù)據(jù)結構的搜索、MPI/OpenMP混合并行方法。

3.1 并行計算壁面距離

現(xiàn)有的壁面距離計算方法基本上都是串行的，一些學者提出的并行計算方法只是在某個步驟、分任務式的并行計算，其本質還是串行方式，可擴展性不強。隨著計算網(wǎng)格規(guī)模的增加，當達到上億甚至十億、百億量級后，現(xiàn)有串行計算方法都會失效，必須采用并行計算。

CFD中的并行計算基本上都是基于區(qū)域分解的，即將離散后整個塊計算區(qū)域分解為不同的空間子分區(qū)，這些子分區(qū)再被分配給不同的進程。非結構網(wǎng)格上的一個分區(qū)是指一個網(wǎng)格塊，而結構網(wǎng)格的分區(qū)可能包括多個網(wǎng)格塊。本文的分區(qū)包括兩個概念——“空間網(wǎng)格分區(qū)”和“壁面網(wǎng)格分區(qū)”?？臻g網(wǎng)格分區(qū)，是采用開源分區(qū)軟件METIS將原始的整個計算域劃分后的子分區(qū)，分區(qū)內的單元是四面體/六面體/金字塔等三維拓撲單元。壁面網(wǎng)格分區(qū)，是指空間網(wǎng)格分區(qū)中包含的壁面面元構成的子分區(qū)，分區(qū)內的單元是三角形/四邊形等二維拓撲單元。圖6是球體外形計算網(wǎng)格的空間網(wǎng)格分區(qū)和壁面網(wǎng)格分區(qū)?？臻g網(wǎng)格分區(qū)用于CFD模擬過程中分配給每個進程，而壁面網(wǎng)格分區(qū)天然地成為并行壁面計算時的分區(qū)。

并行計算壁面距離的流程如圖7。對于N個進程并行的情況，進程i的分區(qū)如果含有壁面面元，則將這些壁面面元壓縮到數(shù)據(jù)容器后發(fā)送到服務器進程(一般為0號進程)。服務器進程將所有進程的壁面面元收集后，將全局壁面面元發(fā)送給其他的所有進程，從而使得每個進程中都有全局壁面面元信息。圖8是寬體客機的壁面信息圖。該飛行器的計算網(wǎng)格被分為2048個子區(qū)，其中的若干個分區(qū)含有壁面面元，將這些壁面面元由服務器進程收集后發(fā)送給每個進程，使得每個進程都含有圖中的壁面面元信息的備份。最后，每個進程中的計算網(wǎng)格根據(jù)其上的全局壁面面元信息，計算每個空間網(wǎng)格單元中心到全局壁面的最近距離。

(a)空間網(wǎng)格分區(qū)

(b)壁面網(wǎng)格分區(qū)

圖7 并行計算流程

圖8 飛行器壁面信息及搜索過程

這種并行計算流程的優(yōu)點是，每個進程在進行壁面距離計算時，由于其中含有全局的壁面信息，因此可以采用幾乎現(xiàn)有的所有方法計算壁面距離，如直接搜索法、循環(huán)盒子法等，改造簡單，適應性和擴展性好。

3.2 基于ADT數(shù)據(jù)結構的搜索算法

圖7的并行壁面距離計算流程中，每個進程需要計算其中的單元中心到壁面面元的最近距離，如果采用直接搜索法，計算時間將取決于網(wǎng)格單元數(shù)最多的進程，大規(guī)模計算的每個分區(qū)單元數(shù)在數(shù)十萬量級，而壁面面元數(shù)可能達到千萬甚至億量級，直接搜索法的遍歷次數(shù)是二者的乘積，計算時間可能長達數(shù)天，效率存在瓶頸。本文采用了一種基于ADT(Alternating Digital Tree)數(shù)據(jù)結構的搜索算法。

ADT數(shù)據(jù)結構實際上是一種特殊的二叉樹[15]，其基本搜索原理是：對于每一個空間多邊形，將其幾何坐標的最小盒子，依次按照笛卡爾坐標三個方向插入到空間中。搜索時，在空間點的幾何坐標盒子中查找含有的空間多邊形，這種搜索方式是二分法在多維情況下的擴展。該方法由于效率極高而在計算幾何方法中得到廣泛的應用。例如，在網(wǎng)格生成、幾何求交等算法中，能實現(xiàn)給定空間范圍中元素的快速查找。

如果直接利用ADT數(shù)據(jù)結構搜索，會出現(xiàn)類似文獻[4]中提到的問題——當搜索對象離搜索目標距離較遠時(如外流問題的遠場區(qū)域)，采用ADT搜索將面臨極為尷尬的問題，即在給定的搜索半徑下要么搜索不到目標，要么搜索到的目標過多而退化為直接搜索法。針對這一問題，提出一種迭代搜索方法。首先，建立全局壁面面元的ADT樹，建立過程是將壁面面元的坐標方盒逐個插入到ADT樹的過程，詳細可參考文獻[15]。然后，遍歷當前進程所含分區(qū)中的所有單元，基于ADT查找計算最近的壁面距離。圖9是ADT搜索過程，其中的Nmax為人為給定的最大ADT節(jié)點數(shù)(其中存儲的數(shù)據(jù)是壁面面元)，其目的是限制搜索到的節(jié)點數(shù)，使其不至于過多而影響效率。

圖9 ADT搜索過程

表1是ADT迭代搜索算法描述。ADT查找的核心是確定搜索半徑，對于復雜外形，很難通過幾何外形直接計算得到搜索半徑。本文采用的這種迭代搜索方法，只需要在最開始給定一個初始搜索半徑，迭代過程中根據(jù)搜索到的ADT節(jié)點數(shù)量不斷地二分調整搜索半徑，直到搜索到理想的節(jié)點數(shù)為止。以圖8的搜索為例，第1次搜索時由于搜索半徑過小，沒有搜索到壁面面元；將其搜索半徑擴大2倍后進行第2次搜索，但導致搜索的壁面面元數(shù)超過了節(jié)點限制數(shù)Nmax；第3次搜索時，搜索半徑設定為前兩次搜索半徑之二分之一后，搜索到的壁面面元數(shù)滿足需求，然后在這些搜索到的節(jié)點內計算最近的壁面距離。

表1 ADT迭代搜索算法

3.3 MPI/OpenMP混合并行計算

上面介紹了基于區(qū)域分解的并行壁面距離計算方法。一般來說，每個進程上分配一個網(wǎng)格分區(qū)，在這種條件下，要求每個進程中都必須包含全局壁面的數(shù)據(jù)備份。這里的“全局壁面數(shù)據(jù)”包含兩類數(shù)據(jù)，一是所有壁面上點的坐標、單元與點的拓撲關系，二是壁面面元構建的ADT樹數(shù)據(jù)。

當進行大規(guī)模網(wǎng)格計算時，壁面面元數(shù)可能達到上千萬，每個進程中都要存儲全局壁面數(shù)據(jù)將會導致計算節(jié)點的內存不夠?？梢圆捎玫慕鉀Q方法是，在一個節(jié)點內只存儲一個全局壁面信息，節(jié)點中的網(wǎng)格分區(qū)逐個、串行地利用全局壁面信息計算壁面距離，但如此一來將導致一個節(jié)點內只有一個進程在進行壁面距離計算，其他進程處于空閑狀態(tài)。

為此，設計了一種MPI/OpenMP混合并行的壁面距離計算模式。CFD軟件通常運行在分布式集群上，采用MPI進行進程間的通信。對于多核處理器來說，節(jié)點內部的內存是共享的，如果在節(jié)點內部也采用MPI通信，可能會由于帶寬的限制帶來通信總量的增加，采用單一的MPI通信往往無法得到最好的效果。其中一種解決方案是采用MPI/OpenMP混合并行，即在節(jié)點內部用OpenMP并行，而在節(jié)點之間用MPI并行，實現(xiàn)節(jié)點間和節(jié)點內部的兩級并行，充分利用消息傳遞和共享內存兩種編程模型的優(yōu)點。每個進程在進行壁面距離計算時，各進程無需通信，因此可以在節(jié)點內部對分區(qū)的循環(huán)遍歷實施OpenMP并行。當節(jié)點內的各分區(qū)都計算完壁面距離后，在節(jié)點間采用MPI進行數(shù)據(jù)同步，從而實現(xiàn)MPI/OpenMP混合并行。表2是混合并行計算算法。

表2 MPI/OpenMP混合并行算法

4 計算方法測試

從精度和速度兩方面，分別在國產(chǎn)定制集群和“天河二號”集群上，對上述的方法進行了測試。測試平臺是國產(chǎn)通用CFD軟件平臺PHengLEI[16-17](前身為HyperFLOW)，該軟件是作者所在團隊研發(fā)的具有完全自主知識產(chǎn)權的通用CFD軟件平臺，主要面向航空航天飛行器設計應用部門。軟件采用C++語言編程、面向對象的軟件設計理念。設計的體系結構和數(shù)據(jù)結構具有良好的可擴展性和通用性，以適應不同類型網(wǎng)格(結構網(wǎng)格、非結構網(wǎng)格、混合網(wǎng)格、重疊網(wǎng)格)的計算。在該軟件平臺上同時開發(fā)了結構、非結構求解器，二者可獨立運行，也可同步耦合計算，即在流場中同時含有結構、非結構網(wǎng)格的情況下，在結構網(wǎng)格上調用結構求解器，在非結構網(wǎng)格上調用非結構求解器。

4.1 計算精度測試

為了測試不同計算方法的精度，采用了結構/非結構兩套高度拉伸的各向異性網(wǎng)格。(1)結構網(wǎng)格。將圖1中的二維四邊形網(wǎng)格直接在展向拉伸，得到三維六面體結構網(wǎng)格。(2)非結構網(wǎng)格。將圖1中的二維四邊形網(wǎng)格三角化，然后展向拉伸，形成圖10中高度拉伸、高度扭曲的各向異性三棱柱網(wǎng)格。這兩套網(wǎng)格都有很大的拉伸率，在壁面附近，長寬比達到2×107。二者的不同點是網(wǎng)格正交性，結構網(wǎng)格在法向具有很好的正交性，全場網(wǎng)格單元的夾角都幾乎為90°；而三棱柱化后的各向異性非結構網(wǎng)格具有很大的扭曲度，在壁面附近最小夾角僅為0.27°。采用一套正交性很好、一套正交性很差的網(wǎng)格，其目的是測試精確方法和近似方法的壁面距離計算精度。

在兩套網(wǎng)格中，y=0平面即為壁面，各棱柱單元中心的y向坐標即為其精確壁面距離，可以用于測試壁面距離計算的誤差。

圖10 平板各向異性三棱柱網(wǎng)格

圖11是不同方法計算的壁面距離比較，其中，“Y coordinate”是y方向的坐標也是壁面距離精確值，“Approximation”是直接用單元中心到壁面面元面心的距離來求壁面距離的近似方法，這種方法在很多代碼中都用到。從圖11(a)的結構網(wǎng)格計算結果看，采用直接搜索方法和ADT搜索方法計算得到的壁面距離和精確值幾乎完全重合，“Approximation”近似法雖然在遠離壁面處有一些誤差，但是在壁面附近(y接近0)和精確值符合很好。從圖11(b)的非結構網(wǎng)格計算結果看，同樣地，直接搜索方法和ADT搜索方法計算得到的壁面距離和精確值幾乎完全重合，但是“Approximation”方法計算的壁面距離在壁面附近誤差很大。

正如第1.1節(jié)所述，對于RANS方程湍流模擬而言，其魯棒性和收斂性主要受壁面附近壁面距離精確性的影響。從精度測試結果看，傳統(tǒng)方法中常用的“Approximation”近似計算方法，對于具有良好正交性的網(wǎng)格，計算的壁面距離在壁面附近能滿足湍流模擬的需要，而在網(wǎng)格正交性不好的情況下，計算精度差，將嚴重影響CFD模擬結果。因此，近似法受網(wǎng)格質量的影響很大，嚴重依賴于網(wǎng)格生成人員的經(jīng)驗。而對于直接搜索法和ADT搜索法，計算的壁面距離精度都是滿足需要的，受網(wǎng)格質量的影響較小。

(a)六面體結構網(wǎng)格

(b)各向異性拉伸非結構網(wǎng)格

4.2 計算效率測試

測試算例是第三屆美國AIAA高升力CFD預測活動中提供的JSM(JAXA Standard Model)標模外形。該外形是客機、運輸機等大型飛行器的簡化研究外形，主要用于評估CFD技術對高升力外形的預測能力，為工程應用提供經(jīng)驗指導。為了說明方法的可擴展性，分別采用了粗/密兩套不同規(guī)模的計算網(wǎng)格在不同的運行環(huán)境下進行了測試。圖12是該外形的幾何外形及計算網(wǎng)格。

圖12 JSM外形及粗計算網(wǎng)格

表3是計算網(wǎng)格、測試環(huán)境和計算效率對比。計算網(wǎng)格為棱柱/四面體/金字塔組成的混合網(wǎng)格，其中，粗網(wǎng)格單元為0.52億，共2048個空間網(wǎng)格分區(qū)，在國產(chǎn)定制集群上測試；密網(wǎng)格單元為33.2億，共12 288個空間網(wǎng)格分區(qū)，在“天河二號”上測試。對于兩套網(wǎng)格，分別采用直接搜索方法和ADT搜索方法計算了壁面距離，二者都采用并行方式計算。可以看到，在兩套網(wǎng)格、兩種測試環(huán)境下，并行ADT搜索方法都有大幅的效率提升。尤其是，直接搜索法時間將隨網(wǎng)格量的增加而指數(shù)增長，當網(wǎng)格量達到上億量級后，壁面距離計算很可能耗費與CFD模擬過程差不多的時間，這將是不可接受的。

正如上文所言，由于ADT搜索過程中每個進程都需要有全局壁面信息的備份，可能帶來內存開銷的大幅增加，上文的MPI/OpenMP混合并行的目的正是為了解決內存開銷問題。表4中給出了粗網(wǎng)格在不同線程/進程數(shù)情況下的內存耗費和測試時間，可以看到，即使對于粗網(wǎng)格而言，如果用純MPI并行計算(2048×1)，每個節(jié)點內的16個CPU核都擁有一個全局壁面信息備份，單節(jié)點內存耗費將達到23.2%；而在MPI/OpenMP混合并行情況下，如果1個節(jié)點內只用2個MPI進程(256×8)，單節(jié)點內存耗費降到7.8%。對于密網(wǎng)格，如果采用純MPI并行計算，其內存耗費約為120 GB，將遠遠超過節(jié)點內存，必須采用混合并行計算。

表3 網(wǎng)格、環(huán)境及測試結果

對于直接搜索法和ADT搜索法，雖然都可以采用3.3節(jié)的MPI/OpenMP混合并行計算方案，但是二者在混合并行時有很大的差異。為了說明這一點，采用了ONERA-M6機翼外形(圖13)，在國產(chǎn)定制集群上進行了測試。計算網(wǎng)格為棱柱/四面體/金字塔組成的混合網(wǎng)格，共46萬單元，分為了64個空間網(wǎng)格分區(qū)，包括43個壁面網(wǎng)格分區(qū)，壁面三角形總數(shù)為1.89萬。之所以采用這個相對較小的算例來比較，是由于直接搜索法計算效率低，如果用大規(guī)模算例測試時間太長。表5是不同的MPI進程數(shù)、OpenMP線程數(shù)情況下兩種方法的測試結果。對于直接搜索法，隨著MPI進程數(shù)減少、OpenMP線程增加，計算時間大幅增加。而對于ADT搜索法，計算時間除了4進程×16線程外，其他情況下計算時間沒有明顯區(qū)別。

表4 MPI/OpenMP計算測試(JSM外形粗網(wǎng)格)

圖13 ONERA-M6機翼計算網(wǎng)格

表5 兩種方法的MPI/OpenMP混合計算時間對比(M6外形)

出現(xiàn)以上結果的原因是：對于直接搜索法，每個空間點都要遍歷所有的壁面面元，因此，如果多個線程訪問同一個內存中的壁面信息備份，必然會出現(xiàn)數(shù)據(jù)競爭；對于ADT搜索法，由于是根據(jù)其空間坐標關系遍歷其坐標周圍一定搜索半徑的區(qū)域，因此搜索時只和該區(qū)域中的其他對象產(chǎn)生競爭，而根據(jù)ADT數(shù)據(jù)結構，即使在同一個區(qū)域中的空間對象，也是依次在3個笛卡爾坐標方向輪換存儲于內存中的，又進一步減少了數(shù)據(jù)競爭的機會。對于4進程×16線程情況，意味著一個節(jié)點內只有一個壁面信息備份，該壁面信息被全部16個線程訪問，會增加數(shù)據(jù)競爭的概率，即使是ADT搜索方法，其計算效率也會受到影響。

5 結 論

針對超大規(guī)模網(wǎng)格下CFD湍流模擬中壁面距離計算效率低的問題，設計了一種基于ADT搜索的并行計算方法，采用了MPI/OpenMP混合并行方法解決內存開銷過大的問題，實現(xiàn)了復雜外形大規(guī)模網(wǎng)格壁面距離高精度、高效并行計算。通過33.2億網(wǎng)格大規(guī)模并行測試，結果表明，即使與并行化的直接搜索方法相比，新方法在效率上也有量級的提升，計算時間完全能夠滿足大規(guī)模CFD湍流模擬的需求，解決了面向下一代E級CFD計算的前處理問題。

該方法還存在待改進之處。ADT搜索實際上是一種三維的二分法，其搜索“方盒”的各個面天然地與Oxyz笛卡爾坐標面平行，若幾何外形的某處局部也與笛卡爾坐標面平行，將導致ADT搜索效率下降，即要么搜索不到任何點,要么搜索到該平面中所有的點。解決方法是，將幾何坐標旋轉一個角度(小于90°)，從而使得壁面面元不再與笛卡爾坐標面平行。下一步工作將繼續(xù)開展深化研究。