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

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

0 引 言

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

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

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

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

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

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

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

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

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

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

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

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

(1)

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

(2)

(3)

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

1.2 壁面距離精確性的影響

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

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

(a)邊界條件[12]

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

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

(a)最大4倍誤差

(b)最大50%誤差

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

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

(b)局部有通量修正

(c)局部無通量修正

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

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

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

壁面距離計(jì)算的通用情況是計(jì)算空間點(diǎn)到凸N邊形的距離，為了提高計(jì)算精度，將凸N邊形的每條邊和多邊形中心連接，形成N個三角形，最終轉(zhuǎn)換為計(jì)算點(diǎn)到N個三角形的距離。分兩步計(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)，對于任意點(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可通過下式計(jì)算：

(6)

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

圖9 ADT搜索過程

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

表1 ADT迭代搜索算法

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

上面介紹了基于區(qū)域分解的并行壁面距離計(jì)算方法。一般來說，每個進(jìn)程上分配一個網(wǎng)格分區(qū)，在這種條件下，要求每個進(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ù)可能達(dá)到上千萬，每個進(jìn)程中都要存儲全局壁面數(shù)據(jù)將會導(dǎo)致計(jì)算節(jié)點(diǎn)的內(nèi)存不夠。可以采用的解決方法是，在一個節(jié)點(diǎn)內(nèi)只存儲一個全局壁面信息，節(jié)點(diǎn)中的網(wǎng)格分區(qū)逐個、串行地利用全局壁面信息計(jì)算壁面距離，但如此一來將導(dǎo)致一個節(jié)點(diǎn)內(nèi)只有一個進(jìn)程在進(jìn)行壁面距離計(jì)算，其他進(jìn)程處于空閑狀態(tài)。

為此，設(shè)計(jì)了一種MPI/OpenMP混合并行的壁面距離計(jì)算模式。CFD軟件通常運(yùn)行在分布式集群上，采用MPI進(jìn)行進(jìn)程間的通信。對于多核處理器來說，節(jié)點(diǎn)內(nèi)部的內(nèi)存是共享的，如果在節(jié)點(diǎn)內(nèi)部也采用MPI通信，可能會由于帶寬的限制帶來通信總量的增加，采用單一的MPI通信往往無法得到最好的效果。其中一種解決方案是采用MPI/OpenMP混合并行，即在節(jié)點(diǎn)內(nèi)部用OpenMP并行，而在節(jié)點(diǎn)之間用MPI并行，實(shí)現(xiàn)節(jié)點(diǎn)間和節(jié)點(diǎn)內(nèi)部的兩級并行，充分利用消息傳遞和共享內(nèi)存兩種編程模型的優(yōu)點(diǎn)。每個進(jìn)程在進(jìn)行壁面距離計(jì)算時，各進(jìn)程無需通信，因此可以在節(jié)點(diǎn)內(nèi)部對分區(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ì)算方法測試

從精度和速度兩方面，分別在國產(chǎn)定制集群和“天河二號”集群上，對上述的方法進(jìn)行了測試。測試平臺是國產(chǎn)通用CFD軟件平臺PHengLEI[16-17](前身為HyperFLOW)，該軟件是作者所在團(tuán)隊(duì)研發(fā)的具有完全自主知識產(chǎn)權(quán)的通用CFD軟件平臺，主要面向航空航天飛行器設(shè)計(jì)應(yīng)用部門。軟件采用C++語言編程、面向?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ì)算。在該軟件平臺上同時開發(fā)了結(jié)構(gòu)、非結(jié)構(gòu)求解器，二者可獨(dú)立運(yùn)行，也可同步耦合計(jì)算，即在流場中同時含有結(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ì)算精度測試

為了測試不同計(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)格都有很大的拉伸率，在壁面附近，長寬比達(dá)到2×107。二者的不同點(diǎn)是網(wǎng)格正交性，結(jié)構(gòu)網(wǎng)格在法向具有很好的正交性，全場網(wǎng)格單元的夾角都幾乎為90°；而三棱柱化后的各向異性非結(jié)構(gòu)網(wǎng)格具有很大的扭曲度，在壁面附近最小夾角僅為0.27°。采用一套正交性很好、一套正交性很差的網(wǎng)格，其目的是測試精確方法和近似方法的壁面距離計(jì)算精度。

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

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

圖11是不同方法計(jì)算的壁面距離比較，其中，“Y coordinate”是y方向的坐標(biāo)也是壁面距離精確值，“Approximation”是直接用單元中心到壁面面元面心的距離來求壁面距離的近似方法，這種方法在很多代碼中都用到。從圖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é)所述，對于RANS方程湍流模擬而言，其魯棒性和收斂性主要受壁面附近壁面距離精確性的影響。從精度測試結(jié)果看，傳統(tǒng)方法中常用的“Approximation”近似計(jì)算方法，對于具有良好正交性的網(wǎng)格，計(jì)算的壁面距離在壁面附近能滿足湍流模擬的需要，而在網(wǎng)格正交性不好的情況下，計(jì)算精度差，將嚴(yán)重影響CFD模擬結(jié)果。因此，近似法受網(wǎng)格質(zhì)量的影響很大，嚴(yán)重依賴于網(wǎng)格生成人員的經(jīng)驗(yàn)。而對于直接搜索法和ADT搜索法，計(jì)算的壁面距離精度都是滿足需要的，受網(wǎng)格質(zhì)量的影響較小。

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

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

4.2 計(jì)算效率測試

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

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

表3是計(jì)算網(wǎng)格、測試環(huán)境和計(jì)算效率對比。計(jì)算網(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搜索方法計(jì)算了壁面距離，二者都采用并行方式計(jì)算?？梢钥吹?，在兩套網(wǎng)格、兩種測試環(huán)境下，并行ADT搜索方法都有大幅的效率提升。尤其是，直接搜索法時間將隨網(wǎng)格量的增加而指數(shù)增長，當(dāng)網(wǎng)格量達(dá)到上億量級后，壁面距離計(jì)算很可能耗費(fèi)與CFD模擬過程差不多的時間，這將是不可接受的。

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

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

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

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

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

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

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

5 結(jié) 論

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

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