趙 鐘, 何 磊, 張 健, 徐慶新, 張來(lái)平

(中國(guó)空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所, 四川 綿陽(yáng) 621000)

0 引 言

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

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

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

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

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

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

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

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

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

1 壁面距離精確計(jì)算的重要性

1.1 湍流模型方程基礎(chǔ)

CFD工程應(yīng)用中常用的湍流模型方程有一方程SA模型和兩方程k-ωSST模型。這里以SA模型為例，說(shuō)明湍流壁面距離對(duì)計(jì)算結(jié)果的影響。SA湍流模型方程為[11-12]：

(1)

上式寫成非守恒物質(zhì)導(dǎo)數(shù)形式為：

(2)

(3)

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

1.2 壁面距離精確性的影響

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

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

(a)邊界條件[12]

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

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

(a)最大4倍誤差

(b)最大50%誤差

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

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

(b)局部有通量修正

(c)局部無(wú)通量修正

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

2 計(jì)算幾何基礎(chǔ)

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

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

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

圖4 點(diǎn)投影到三角形

圖5中，采用參數(shù)化坐標(biāo)方法判斷投影點(diǎn)PI是否在三角形內(nèi)[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ù)化坐標(biāo)，對(duì)于任意點(diǎn)PI=V(s,t)，若位于三角形中，必滿足下式：

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

(5)

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

(6)

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

圖5 判斷點(diǎn)是否在三角形內(nèi)部

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

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

3.1 并行計(jì)算壁面距離

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

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

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

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

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

圖7 并行計(jì)算流程

圖8 飛行器壁面信息及搜索過(guò)程

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

3.2 基于ADT數(shù)據(jù)結(jié)構(gòu)的搜索算法

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

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

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

圖9 ADT搜索過(guò)程

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

表1 ADT迭代搜索算法

3.3 MPI/OpenMP混合并行計(jì)算

上面介紹了基于區(qū)域分解的并行壁面距離計(jì)算方法。一般來(lái)說(shuō)，每個(gè)進(jìn)程上分配一個(gè)網(wǎng)格分區(qū)，在這種條件下，要求每個(gè)進(jìn)程中都必須包含全局壁面的數(shù)據(jù)備份。這里的“全局壁面數(shù)據(jù)”包含兩類數(shù)據(jù)，一是所有壁面上點(diǎn)的坐標(biāo)、單元與點(diǎn)的拓?fù)潢P(guān)系，二是壁面面元構(gòu)建的ADT樹數(shù)據(jù)。

當(dāng)進(jìn)行大規(guī)模網(wǎng)格計(jì)算時(shí)，壁面面元數(shù)可能達(dá)到上千萬(wàn)，每個(gè)進(jìn)程中都要存儲(chǔ)全局壁面數(shù)據(jù)將會(huì)導(dǎo)致計(jì)算節(jié)點(diǎn)的內(nèi)存不夠。可以采用的解決方法是，在一個(gè)節(jié)點(diǎn)內(nèi)只存儲(chǔ)一個(gè)全局壁面信息，節(jié)點(diǎn)中的網(wǎng)格分區(qū)逐個(gè)、串行地利用全局壁面信息計(jì)算壁面距離，但如此一來(lái)將導(dǎo)致一個(gè)節(jié)點(diǎn)內(nèi)只有一個(gè)進(jìn)程在進(jìn)行壁面距離計(jì)算，其他進(jìn)程處于空閑狀態(tài)。

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

表2 MPI/OpenMP混合并行算法

4 計(jì)算方法測(cè)試

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

4.1 計(jì)算精度測(cè)試

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

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

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

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

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

(a)六面體結(jié)構(gòu)網(wǎng)格

(b)各向異性拉伸非結(jié)構(gòu)網(wǎng)格

4.2 計(jì)算效率測(cè)試

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

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

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

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

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

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

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

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

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

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

5 結(jié) 論

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

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