葉 靚, 張 穎, 楊 碩, 董 軍

(中國(guó)航空工業(yè)空氣動(dòng)力研究院 高速高雷諾數(shù)氣動(dòng)力航空科技重點(diǎn)實(shí)驗(yàn)室, 遼寧 沈陽(yáng) 110034)

0 引 言

旋翼流場(chǎng)是復(fù)雜非定常、渦占主導(dǎo)地位的流場(chǎng)。為了進(jìn)行流場(chǎng)細(xì)節(jié)(特別是槳尖渦)捕捉，相應(yīng)CFD計(jì)算使用的網(wǎng)格數(shù)量較多，為了獲得周期解，數(shù)值模擬需要進(jìn)行多個(gè)旋翼旋轉(zhuǎn)周期，而每個(gè)旋轉(zhuǎn)周期上又要進(jìn)行多個(gè)物理時(shí)間計(jì)算。一般情況下，由于槳葉/槳葉之間及槳葉/機(jī)身之間存在大幅相對(duì)運(yùn)動(dòng)，計(jì)算通常基于嵌套網(wǎng)格系統(tǒng)，這進(jìn)一步提高了計(jì)算的復(fù)雜度。旋翼流場(chǎng)數(shù)值模擬通常占用較大的計(jì)算資源且十分耗時(shí)[1]。

并行計(jì)算技術(shù)是有效縮短計(jì)算時(shí)間消耗的重要途徑。近年來(lái)，流場(chǎng)計(jì)算方面的并行技術(shù)已經(jīng)發(fā)展得相對(duì)成熟，如文獻(xiàn)[2-4]?？紤]到基于嵌套網(wǎng)格的旋翼流場(chǎng)數(shù)值模擬在每個(gè)計(jì)算物理時(shí)間上都要進(jìn)行網(wǎng)格裝配，較好的方法將會(huì)有效縮減網(wǎng)格裝配時(shí)間對(duì)流場(chǎng)計(jì)算時(shí)間的占比，從而進(jìn)一步提高整體并行的計(jì)算效率。

關(guān)于網(wǎng)格裝配，先前已經(jīng)開(kāi)展過(guò)不少工作，發(fā)展的方法可以分為兩大類(lèi)，一是顯式方法，即首先進(jìn)行網(wǎng)格的“洞切割”，圍繞生成的洞邊界，構(gòu)建插值邊界，進(jìn)行插值邊界的貢獻(xiàn)單元搜索，建立插值單元和貢獻(xiàn)單元的對(duì)應(yīng)關(guān)系。二是隱式方法，即對(duì)所有網(wǎng)格單元，均進(jìn)行貢獻(xiàn)單元搜索工作，通過(guò)比較單元及其貢獻(xiàn)單元之間的幾何關(guān)系(體積、長(zhǎng)寬比等)來(lái)自動(dòng)確定數(shù)據(jù)交換的邊界。顯式方法由于需要搜索的貢獻(xiàn)單元數(shù)目較少，因而具備較快的計(jì)算速度，但由于需要算法給定洞切割范圍，也需要進(jìn)行必須的幾何交叉問(wèn)題判斷及付出相應(yīng)的計(jì)算程序復(fù)雜度代價(jià)。比較經(jīng)典的方法如洞映射方法[5]、目標(biāo)X射線(xiàn)方法[6]和純幾何切割的方法等[7]。對(duì)于隱式方法[8]，由于需要全局貢獻(xiàn)單元搜索，時(shí)間消耗較為嚴(yán)重?？梢钥闯觯瑹o(wú)論是基于顯式還是隱式網(wǎng)格裝配方法，貢獻(xiàn)單元搜索都是不可缺少的步驟，目前發(fā)展的貢獻(xiàn)單元搜索方法包括窮舉法(Brute Force)、最近鄰居法Nearest Neighbor(NN)[9]，Neighbor to Neighbor(N2N)方法[10]，Alternating Digital Tree(ADT)方法[11]、結(jié)構(gòu)輔助網(wǎng)格方法Structured Auxiliary Mesh(SAM)[12]、近似反映射方法Approximate Inverse Map (AIM)[13]和精確反映射Exact Inverse Map(EIM)[14]等。其中ADT方法是一種較為經(jīng)典的方法，魯棒性較好，而使用結(jié)構(gòu)輔助網(wǎng)格方法能夠提供比ADT方法快一階的計(jì)算速度[15]，是一種高效的計(jì)算方法。N2N方法的問(wèn)題是搜索路徑可能遇到邊界面，需要附加計(jì)算重啟過(guò)程。精確反映射方法速度很快，但仍涉及很多復(fù)雜的幾何交叉判斷問(wèn)題。國(guó)內(nèi)研究者們的研究成果包括 “透視圖”方法[16]、“距離縮減”[17]和一些其他方法[18-20]等。

考慮到旋翼流場(chǎng)計(jì)算的多周期多物理時(shí)間步特性，網(wǎng)格裝配對(duì)整體計(jì)算效率影響較大。為縮短網(wǎng)格裝配時(shí)間，降低其對(duì)于整體流場(chǎng)計(jì)算時(shí)間的占比，本文進(jìn)行算法方面的改進(jìn)，目的是提高網(wǎng)格裝配效率的同時(shí)盡量減少幾何判斷和計(jì)算代碼的復(fù)雜程度，網(wǎng)格裝配算法具有一般性，不受物體形狀和網(wǎng)格拓?fù)涞南拗?，并使得?jì)算的內(nèi)存消耗在可容忍的范圍之內(nèi)。具體工作如下：一是采用基于構(gòu)造曲面外形的壓縮存儲(chǔ)的輔助結(jié)構(gòu)網(wǎng)格進(jìn)行并行洞切割，方法速度快且能夠適應(yīng)任意物面形狀。二是基于可自適應(yīng)的結(jié)構(gòu)輔助網(wǎng)格Adaptive Structured Auxiliary Mesh(ASAM)，完成貢獻(xiàn)單元初始位置定位，并結(jié)合N2N搜索和局部窮舉法，形成了無(wú)判斷邊界的并行貢獻(xiàn)單元確定方法，稱(chēng)為Inverse Map No Boundary Adjudication(IM_NBA)。通過(guò)以上措施，實(shí)現(xiàn)了網(wǎng)格裝配過(guò)程的最大限度并行化。

在流場(chǎng)計(jì)算方面，本文采用并行化的隱式雙時(shí)間方法數(shù)值求解Navier-Stokes方程，無(wú)黏通量采用ROE格式和梯度重構(gòu)方法計(jì)算，湍流項(xiàng)用Spalart-Allmaras模型模擬。

為了便于演示，文中首先給出了一個(gè)二維嵌套網(wǎng)格洞切割和貢獻(xiàn)單元搜索的算例。然后基于以上方法，對(duì)不同的旋翼模型網(wǎng)格系統(tǒng)，在工作站上，通過(guò)計(jì)算測(cè)試給出了采用改進(jìn)的網(wǎng)格裝配方法時(shí)的裝配過(guò)程并行加速效果及裝配時(shí)間對(duì)計(jì)算時(shí)間的占比。

1 并行網(wǎng)格裝配方法

1.1 網(wǎng)格分區(qū)

對(duì)于網(wǎng)格數(shù)目和拓?fù)洳蛔兦闆r下，網(wǎng)格分區(qū)僅需在計(jì)算開(kāi)始的時(shí)候進(jìn)行一次，即把網(wǎng)格劃分成多個(gè)區(qū)域，并把網(wǎng)格信息傳遞到對(duì)應(yīng)處理器上。對(duì)于多塊嵌套的旋翼網(wǎng)格系統(tǒng)，設(shè)背景網(wǎng)格塊編號(hào)為B1，各個(gè)小網(wǎng)格塊編號(hào)分別為B2，B3，…，BN，其中N為網(wǎng)格塊總數(shù)，網(wǎng)格分區(qū)采用Metis工具[21]對(duì)每塊網(wǎng)格分別進(jìn)行。分區(qū)完成后的網(wǎng)格為BN_1，BN_2，…，BN_np，其中np為處理器個(gè)數(shù)。分布在第一個(gè)處理器上的子網(wǎng)格塊為B1_1，B2_1，…，BN_1；第np個(gè)處理器上的子網(wǎng)格塊為B1_np，B2_np，…，BN_np。網(wǎng)格分區(qū)的邏輯視圖如圖1所示。

(a)

(b)

需要注意的是，進(jìn)行網(wǎng)格劃分后要自行生成網(wǎng)格分區(qū)之間的交疊邊界，以實(shí)現(xiàn)處理器之間的數(shù)據(jù)交互。由于本文計(jì)算采用二階格式故構(gòu)建兩層單元的數(shù)據(jù)交換邊界。如圖2所示，其中有陰影部分為兩塊分區(qū)網(wǎng)格向外擴(kuò)展的的重疊區(qū)域。

圖2 網(wǎng)格子塊擴(kuò)展邊界示意圖Fig.2 Schematic of halo cells for grid sub-blocks

1.2 網(wǎng)格裝配

1.2.1 并行洞切割

并行洞切割過(guò)程包括“洞構(gòu)造面”(Scons)構(gòu)建，輔助結(jié)構(gòu)網(wǎng)格生成、洞切割及插值邊界確定等。

1.2.1.1 “洞構(gòu)造面”構(gòu)建

由主處理器讀入整體網(wǎng)格，對(duì)每塊網(wǎng)格上的所有物面，向外推出幾層，構(gòu)建出一層網(wǎng)格面，用來(lái)作為洞邊界生成的判據(jù)。該過(guò)程只需在計(jì)算開(kāi)始的時(shí)候進(jìn)行一次。把該構(gòu)造面的幾何信息廣播給所有處理器。

在所有處理器上均儲(chǔ)存該構(gòu)造面信息。由于構(gòu)造面只是一層網(wǎng)格面，因而對(duì)計(jì)算內(nèi)存的需求較小。

1.2.1.2 洞切割輔助結(jié)構(gòu)網(wǎng)格生成和壓縮存儲(chǔ)

以構(gòu)造面為基準(zhǔn)，生成“輔助結(jié)構(gòu)網(wǎng)格立柱”，根據(jù)某種規(guī)則(如每塊網(wǎng)格上的單元數(shù))確定輔助結(jié)構(gòu)網(wǎng)格的“分辨率”(Nresolution)。輔助結(jié)構(gòu)網(wǎng)格的空間范圍由Scons的最大最小坐標(biāo)確定。設(shè)與坐標(biāo)軸y軸平行的方向?yàn)椤傲⒅较颉?Direction)，該方向作為壓縮存儲(chǔ)方向，其余方向?yàn)椤傲⒅矫娣较颉?，?duì)三維問(wèn)題，空間三個(gè)方向(x，y，z)排除“立柱方向”后剩余(x，z)兩個(gè)方向，申請(qǐng)數(shù)組大小Nresolution×Nresolution；對(duì)于二維問(wèn)題，剩余一個(gè)方向(x)，申請(qǐng)數(shù)組大小為Nresolution。每個(gè)輔助結(jié)構(gòu)網(wǎng)格“立柱單元”(以平面方向?yàn)榛?，以立柱方向?yàn)楦?坐標(biāo)范圍由Scons和“分辨率”確定，根據(jù)每塊構(gòu)造曲面網(wǎng)格面在非“立柱方向”的最大最小坐標(biāo)，即可確定該網(wǎng)格面所交疊的輔助結(jié)構(gòu)網(wǎng)格立柱。

輔助結(jié)構(gòu)網(wǎng)格采用壓縮存儲(chǔ)方式，即對(duì)于每個(gè)“立柱單元”，僅標(biāo)記不連續(xù)的包含構(gòu)造面的單元位置，并存儲(chǔ)該單元下一單元是否落洞的信息，如此可以節(jié)約大量?jī)?nèi)存。下面以二維翼型圖的方式進(jìn)行說(shuō)明。

圖3(a)為根據(jù)構(gòu)造網(wǎng)格面標(biāo)記的輔助結(jié)構(gòu)網(wǎng)格中的單元(標(biāo)記了M)，圖3(b)給出了不使用壓縮存儲(chǔ)方式時(shí)，需要記錄的結(jié)構(gòu)網(wǎng)格單元位置(標(biāo)記了M)。圖3(c)給出了y方向立柱壓縮存儲(chǔ)時(shí)需要記錄的輔助結(jié)構(gòu)網(wǎng)格單元位置，圖3(d)給出了x方向立柱壓縮存儲(chǔ)時(shí)需要記錄的輔助結(jié)構(gòu)網(wǎng)格單元位置(M(0)表示本輔助網(wǎng)格被記錄，且在該立柱中該單元后面到下一個(gè)記錄單元之間的網(wǎng)格為落洞，M(1)表示不落洞)。如此可以節(jié)約大量的存儲(chǔ)空間。注意采用y方向立柱壓縮存儲(chǔ)的圖(c)與(a)比較去除了連續(xù)重復(fù)位置(M(1，2)，M(1，3)，M(1，4)，M(1，5)，M(4，3)，M(4，4)，M(5，2))，X方向壓縮存儲(chǔ)的(d)與(a)比較中去除了連續(xù)重復(fù)位置(M(2，3)，M(3，3)，M(4，3)，M(2，6))。進(jìn)行壓縮存儲(chǔ)的存儲(chǔ)量在二維情況下是O(Nresolution)，三維情況下是O(Nresolution×Nresolution)。

結(jié)構(gòu)網(wǎng)格單元后方的單元是否落洞的判斷方法是從后面結(jié)構(gòu)網(wǎng)格的格心引-Direction方向射線(xiàn)，按照與該立柱內(nèi)構(gòu)造面交點(diǎn)的個(gè)數(shù)奇數(shù)還是偶數(shù)來(lái)判定，僅需判斷有限的幾個(gè)網(wǎng)格面，計(jì)算量很小。

可以發(fā)現(xiàn)本文的壓縮存儲(chǔ)法對(duì)于構(gòu)造面是“凹面”的情況仍然有效，具有一般性。

(a) Structured auxiliary mesh marked by scons

(b) Theoretical storage

(c) Actual storage(Y column)

(d) Actual storage(X column)

1.2.1.3 “洞切割”和插值邊界元interpolate point(IP)確定

在每個(gè)處理器上，對(duì)該處理器上的所有網(wǎng)格循環(huán)，通過(guò)調(diào)用輔助結(jié)構(gòu)網(wǎng)格壓縮存儲(chǔ)數(shù)組及其檢索數(shù)組，根據(jù)簡(jiǎn)單坐標(biāo)比較，即可確定單元是否是洞單元，而后把所有與洞單元鄰接的非洞單元標(biāo)記成插值邊界元，進(jìn)行后續(xù)的貢獻(xiàn)單元搜索工作。

1.2.2 并行貢獻(xiàn)單元搜索

確定了插值邊界單元后，需要尋找這些單元的貢獻(xiàn)單元。這里提出了一種可自適應(yīng)的輔助結(jié)構(gòu)網(wǎng)格方法Adaptive Structured Auxiliary Mesh(ASAM)，用輔助結(jié)構(gòu)網(wǎng)格方法給定貢獻(xiàn)單元搜索的起始點(diǎn)，給定搜索開(kāi)始點(diǎn)后，采用N2N方法進(jìn)行貢獻(xiàn)單元確定，對(duì)于搜索路徑遇到邊界面的問(wèn)題，采用了單個(gè)結(jié)構(gòu)輔助網(wǎng)格內(nèi)的單元局部窮舉法，避免了復(fù)雜幾何交叉的判定(輔助結(jié)構(gòu)網(wǎng)格自適應(yīng)過(guò)程也可以有效排除邊界面)。稱(chēng)之為不判斷邊界的反映射方法Inverse Map No Boundary Adjudication(IM_NBA)。具體步驟如下：

1) 網(wǎng)格塊外邊界范圍信息廣播。在各個(gè)處理器之間交換信息，使得每個(gè)處理器均知道其他處理器上所有網(wǎng)格塊的最大最小坐標(biāo)范圍。

2) 插值邊界單元交互。對(duì)任意處理器，對(duì)其上面所有要進(jìn)行貢獻(xiàn)單元搜索的單元IP，判斷這些待尋找貢獻(xiàn)單元的坐標(biāo)是否在其他處理器上網(wǎng)格邊界范圍內(nèi)，如果在，打包發(fā)送給對(duì)應(yīng)的處理器。

3) 輔助結(jié)構(gòu)網(wǎng)格定位。對(duì)任一處理器n，子網(wǎng)格塊為B1_n，B2_n，……BN_n，對(duì)其每個(gè)子網(wǎng)格塊，按照某一基準(zhǔn)(如子網(wǎng)格塊的單元數(shù))建立一定分辨率的輔助結(jié)構(gòu)網(wǎng)格(注意：該輔助結(jié)構(gòu)網(wǎng)格的最大最小坐標(biāo)由其他處理器發(fā)送過(guò)來(lái)的插值邊界元IP的最大最小坐標(biāo)動(dòng)態(tài)確定，而不是該塊網(wǎng)格在該處理器上的坐標(biāo)范圍確定，如此大大縮小了輔助結(jié)構(gòu)網(wǎng)格的范圍，排除了大量的不可能待搜索單元。)。通過(guò)對(duì)待搜索插值單元及網(wǎng)格塊上的網(wǎng)格單元進(jìn)行一次大循環(huán)，建立壓縮存儲(chǔ)數(shù)組，輸入輔助結(jié)構(gòu)網(wǎng)格單元號(hào)，即可取出與該輔助結(jié)構(gòu)網(wǎng)格單元有空間交疊的候選貢獻(xiàn)單元網(wǎng)格及落入該輔助結(jié)構(gòu)網(wǎng)格空間范圍內(nèi)的待搜索單元。實(shí)際存儲(chǔ)時(shí)，不存儲(chǔ)所有的輔助結(jié)構(gòu)網(wǎng)格，而是僅對(duì)用到的輔助結(jié)構(gòu)網(wǎng)格單元(輔助結(jié)構(gòu)網(wǎng)格內(nèi)有待搜索插值點(diǎn)及候選貢獻(xiàn)單元)進(jìn)行存儲(chǔ)，這樣可以降低存儲(chǔ)容量，壓縮存儲(chǔ)數(shù)組通過(guò)建立的索引數(shù)組來(lái)查詢(xún)。

由以上操作可知，選擇合適的輔助結(jié)構(gòu)網(wǎng)格分辨率是較為重要的，太大的分辨率導(dǎo)致每個(gè)輔助結(jié)構(gòu)網(wǎng)格內(nèi)存儲(chǔ)了太多可能的候選貢獻(xiàn)單元，使得搜索路徑增長(zhǎng)，搜索時(shí)間效率下降。太小的分辨率將會(huì)帶來(lái)索引數(shù)組的容量增加，也無(wú)疑增加了存儲(chǔ)負(fù)擔(dān)。

4) 無(wú)邊界交叉判斷的N2N貢獻(xiàn)單元搜索和輔助結(jié)構(gòu)網(wǎng)格細(xì)化。對(duì)任一輔助結(jié)構(gòu)網(wǎng)格單元，取出其包含的待搜索單元及候選貢獻(xiàn)單元，對(duì)每個(gè)待搜索的插值單元，以該輔助網(wǎng)格內(nèi)的候選貢獻(xiàn)單元范圍為邊界，使用N2N方法，進(jìn)行貢獻(xiàn)單元搜索。搜索時(shí)，把每一個(gè)途經(jīng)的單元進(jìn)行記錄，當(dāng)搜索路徑遭遇邊界面(外場(chǎng)面、網(wǎng)格分區(qū)面等)，搜索路徑指向邊界面外時(shí)，需要在該輔助結(jié)構(gòu)網(wǎng)格單元內(nèi)的候選貢獻(xiàn)單元中取出路徑未經(jīng)歷過(guò)的單元，重新啟動(dòng)N2N搜索，此過(guò)程為局部的窮舉法，如果N2N搜索指向一個(gè)已經(jīng)途徑過(guò)的候選單元時(shí)，再次重啟N2N搜索，直至找到需要的貢獻(xiàn)單元，或者遍歷經(jīng)過(guò)了該輔助結(jié)構(gòu)網(wǎng)格中的所有貢獻(xiàn)單元，仍然沒(méi)有找到貢獻(xiàn)單元時(shí)，認(rèn)為待搜索點(diǎn)在該處理器的該網(wǎng)格塊上沒(méi)有貢獻(xiàn)單元。通過(guò)以上步驟，完全避免了繁瑣的邊界判斷問(wèn)題。還需注意的是，對(duì)于待搜索點(diǎn)的位置位于兩個(gè)或多個(gè)網(wǎng)格的交線(xiàn)、交點(diǎn)很近時(shí)，N2N搜索的路徑將會(huì)在可能的貢獻(xiàn)單元之間跳躍(僅對(duì)三維情況，四邊形網(wǎng)格面四點(diǎn)不共面)，這時(shí)啟用模糊判斷規(guī)則，即搜索路徑出現(xiàn)重復(fù)時(shí)，即認(rèn)為鎖定了貢獻(xiàn)單元。圖4給出了某個(gè)輔助結(jié)構(gòu)網(wǎng)格單元內(nèi)的待搜索貢獻(xiàn)單元(I1，I2,I3)及其候選貢獻(xiàn)單元，其中I1在邊界外，I2、I3有貢獻(xiàn)單元，并分別給出了三個(gè)點(diǎn)的搜索路徑。

(a)

(b)

對(duì)于一個(gè)輔助結(jié)構(gòu)網(wǎng)格內(nèi)包含較多候選貢獻(xiàn)單元的情況，搜索路徑會(huì)比較長(zhǎng)，時(shí)間效率下降，因而可以對(duì)輔助結(jié)構(gòu)網(wǎng)格再次細(xì)分。對(duì)細(xì)分化后的結(jié)構(gòu)網(wǎng)格單元，仍需做一次前處理，把屬于各個(gè)細(xì)分化輔助結(jié)構(gòu)網(wǎng)格單元中的待搜索點(diǎn)和候選貢獻(xiàn)單元放入對(duì)應(yīng)位置，取出細(xì)分化結(jié)構(gòu)網(wǎng)格，即可取出其對(duì)應(yīng)的待搜索點(diǎn)和候選貢獻(xiàn)單元。進(jìn)行輔助結(jié)構(gòu)網(wǎng)格細(xì)分化后，貢獻(xiàn)單元的搜索路徑可能會(huì)大幅縮短，如圖5所示。注意，附加的前處理過(guò)程仍需要耗費(fèi)計(jì)算時(shí)間，在實(shí)際計(jì)算時(shí)，僅對(duì)輔助結(jié)構(gòu)網(wǎng)格內(nèi)包含較多數(shù)量的待搜索點(diǎn)和貢獻(xiàn)單元的情況進(jìn)行網(wǎng)格細(xì)分，如此在儲(chǔ)存容量和搜索速度之間進(jìn)行了平衡。細(xì)化過(guò)程還會(huì)去除一些輔助結(jié)構(gòu)網(wǎng)格(輔助結(jié)構(gòu)網(wǎng)格體積縮小后，其內(nèi)不再包含待搜索點(diǎn)或候選貢獻(xiàn)單元)。

(a)

(b)

在進(jìn)行輔助網(wǎng)格自適應(yīng)時(shí)，由于其是在各個(gè)單個(gè)輔助結(jié)構(gòu)網(wǎng)格內(nèi)進(jìn)行搜索時(shí)進(jìn)行，僅針對(duì)該輔助結(jié)構(gòu)網(wǎng)格內(nèi)的插值單元和候選貢獻(xiàn)單元進(jìn)行再劃分歸類(lèi)，過(guò)程針對(duì)某一個(gè)輔助結(jié)構(gòu)網(wǎng)格內(nèi)單元進(jìn)行動(dòng)態(tài)內(nèi)存申請(qǐng)和釋放；另一方面，由于網(wǎng)格細(xì)分化，導(dǎo)致一些細(xì)分的輔助結(jié)構(gòu)網(wǎng)格內(nèi)不再擁有待搜索點(diǎn)，其交疊的候選貢獻(xiàn)單元也可以不再存儲(chǔ)。因而只要限定細(xì)分化的上限為一有限數(shù)值，則對(duì)整體存儲(chǔ)容量不會(huì)產(chǎn)生較大影響，甚至能夠有效縮減存儲(chǔ)容量，如此可以采用分辨率較低的基礎(chǔ)輔助結(jié)構(gòu)網(wǎng)格，獲得較快速的貢獻(xiàn)單元搜索效果。

1.2.3 貢獻(xiàn)單元信息回收和確定

如果一個(gè)待搜索插值單元的貢獻(xiàn)單元或貢獻(xiàn)單元的鄰元是邊界元，則認(rèn)為該貢獻(xiàn)單元是不合理的(二階格式構(gòu)建)，認(rèn)為未找到貢獻(xiàn)單元。此方法保證了網(wǎng)格分區(qū)邊界的擴(kuò)展邊界元不能作為貢獻(xiàn)單元。一個(gè)待搜索插值單元在多個(gè)處理器網(wǎng)格塊上的貢獻(xiàn)單元通過(guò)MPI數(shù)據(jù)傳遞回收，如果該待搜索點(diǎn)在所有處理器上均沒(méi)有貢獻(xiàn)單元，認(rèn)為此點(diǎn)沒(méi)有貢獻(xiàn)單元。

圖6為2個(gè)進(jìn)程時(shí)的并行網(wǎng)格裝配過(guò)程的計(jì)算流程圖。

圖6 并行網(wǎng)格裝配流程Fig.6 Schematic of parallel grid assembly

2 并行流場(chǎng)計(jì)算方法

2.1 控制方程形式

基于有限體積法的慣性系下三維Navier-Stokes方程可以寫(xiě)為

F= [ρV·nρV(V·n)·n+Pn

ρE(V·n)+PV·n]T

其中，P、ρ、V分別為壓強(qiáng)、密度和速度；n為單位外法矢；體積為Ω，邊界為?Ω；W、H分別為守恒變量和對(duì)流矢量。

2.2 隱式雙時(shí)間Hybrid-LUSGS方法推進(jìn)

對(duì)于旋翼非定常流場(chǎng)，采用隱式雙時(shí)間方法來(lái)進(jìn)行時(shí)間推進(jìn)。引入二階時(shí)間精度的雙時(shí)間法[22]后，控制方程可以寫(xiě)為:

式中：m為偽時(shí)間；n為物理時(shí)間；Δτ為偽時(shí)間增量；Δt為物理時(shí)間增量；對(duì)于控制方程AΔW=-R，其中A為nelt×nelt矩陣(nelt為網(wǎng)格單元的總數(shù))。隱式無(wú)矩陣存儲(chǔ)的LU-SGS方法[23]如下。

AΔW=(L+D+U)ΔW=-R(3)

向前掃略步：

(L+D)ΔW*=-R(4)

向后掃略步：

(D+U)ΔW=DΔW*(5)

在并行計(jì)算時(shí)，由于網(wǎng)格分區(qū)，每個(gè)線(xiàn)程上的網(wǎng)格不是整體網(wǎng)格，因而原始的LU-SGS方法不再適用，本文采用以下步驟(Hybrid LUSGS)[24]：

1) 雅克比迭代方式計(jì)算每個(gè)單元上的守恒量增量：

3) 各個(gè)處理器上計(jì)算守恒變量的增量：

更新守恒變量(網(wǎng)格分區(qū)擴(kuò)展邊界元不更新)；

4) 不同處理器之間交換插值邊界的守恒變量。

2.3 通量計(jì)算

使用動(dòng)網(wǎng)格上的ROE格式[25]進(jìn)行無(wú)粘通量計(jì)算，為保證計(jì)算空間精度，交接面的通量確定采用數(shù)值重構(gòu)技術(shù)。用Venkatakrishnan限制器[26]避免插值過(guò)程中的過(guò)沖。

粘性通量計(jì)算采用中心格式，其中渦粘性系數(shù)的計(jì)算采用Spalart-Allmaras湍流模型[27]，如下：

湍流模型方程與主控方程采用松耦合的方式進(jìn)行，時(shí)間離散方法與主控方程相同。需要注意的是由于槳葉等物體間存在相對(duì)運(yùn)動(dòng)，任意網(wǎng)格點(diǎn)的壁面距離需要在每個(gè)物理時(shí)間重新計(jì)算，該過(guò)程也可以完全并行進(jìn)行(各個(gè)處理器上都存儲(chǔ)壁面坐標(biāo)即可)。

2.4 嵌套網(wǎng)格間數(shù)據(jù)交換

插值單元物理量的值通過(guò)貢獻(xiàn)單元的鄰居上物理量的值插值得到。

3 算例研究

3.1 二維多段翼型示例

為了方便展示本文基于自適應(yīng)輔助結(jié)構(gòu)網(wǎng)格貢獻(xiàn)單元搜索方法的效果，針對(duì)二維多段翼型網(wǎng)格進(jìn)行了洞切割和貢獻(xiàn)單元搜索測(cè)試。其中主翼網(wǎng)格作為背景網(wǎng)格，拓?fù)洳捎没旌暇W(wǎng)格形式(四邊形/三角形)，襟翼和縫翼網(wǎng)格塊均采用四邊形網(wǎng)格。以物面為邊界，向外擴(kuò)展4層構(gòu)建構(gòu)造曲面，以構(gòu)造曲面為依據(jù)生成洞切割輔助結(jié)構(gòu)網(wǎng)格，圖7為用于洞切割的輔助結(jié)構(gòu)網(wǎng)格，未加粗部分為理論存儲(chǔ)(黑色)，加粗部分(藍(lán)色)為采用壓縮存儲(chǔ)時(shí)實(shí)際存貯的輔助網(wǎng)格?？梢?jiàn)采用壓縮存儲(chǔ)技術(shù)可以大大縮減內(nèi)存消耗，使得采用更加致密的輔助結(jié)構(gòu)網(wǎng)格以獲取更好的洞切割效果成為可能。圖8是洞切割后用于計(jì)算的三段翼型網(wǎng)格示意圖。

圖9給出了采用不同的輔助結(jié)構(gòu)網(wǎng)格進(jìn)行貢獻(xiàn)單元初始位置定位的示意圖(不同顏色的網(wǎng)格線(xiàn)代表

圖7 洞切割輔助結(jié)構(gòu)網(wǎng)格示意圖Fig.7 Schematic of hole cutting SAM

圖8 洞切割后的三段翼型網(wǎng)格Fig.8 Multi-wing grid after hole cutting

其屬于不同的處理器)。其中三種情況分別對(duì)應(yīng)著均勻輔助結(jié)構(gòu)網(wǎng)格(同一處理器，同一網(wǎng)格塊使用尺度相同的結(jié)構(gòu)網(wǎng)格輔助單元，劃分的尺度人工給定)(圖9(a))，自適應(yīng)輔助結(jié)構(gòu)網(wǎng)格(圖9(b))和極致密輔助結(jié)構(gòu)網(wǎng)格(圖9(c))。圖中還給出了輔助結(jié)構(gòu)網(wǎng)格中對(duì)應(yīng)的候選貢獻(xiàn)單元(黑色網(wǎng)格線(xiàn))?？梢?jiàn)，對(duì)于均勻的輔助結(jié)構(gòu)網(wǎng)格，每個(gè)結(jié)構(gòu)網(wǎng)格單元內(nèi)包含了較多的候選貢獻(xiàn)單元，使得每個(gè)結(jié)構(gòu)輔助網(wǎng)格內(nèi)的貢獻(xiàn)單元搜索過(guò)程較長(zhǎng)；圖9(b)是圖9(a)網(wǎng)格依據(jù)其中每個(gè)輔助結(jié)構(gòu)網(wǎng)格單元中的待搜索點(diǎn)數(shù)和候選貢獻(xiàn)單元數(shù)進(jìn)行自適應(yīng)的結(jié)果，在較小的代價(jià)下，縮小了貢獻(xiàn)單元搜索的范圍；圖9(c)采用了極致密的輔助結(jié)構(gòu)網(wǎng)格(人工給定極小的劃分尺度)，幾乎一次鎖定了需要的貢獻(xiàn)單元。需要注意的是，采用極致密的輔助結(jié)構(gòu)網(wǎng)格，對(duì)于三維情況，內(nèi)存開(kāi)銷(xiāo)將過(guò)于巨大，因而此處給出的僅是對(duì)于二維計(jì)算問(wèn)題的演示性說(shuō)明。

(a) uniform

(b) adaptive

(c) fine

3.2 不同背景網(wǎng)格下的“Caradonna & Tung旋翼”網(wǎng)格并行裝配計(jì)算效果比較

“Caradonna & Tung旋翼”[28]，槳葉2片，翼型NACA0012，展弦比6，無(wú)扭轉(zhuǎn)尖削。計(jì)算狀態(tài)為槳尖馬赫數(shù)0.877，安裝角8°。網(wǎng)格由背景和兩塊槳葉網(wǎng)格組成。為方便比較，使用兩種不同的背景網(wǎng)格(BA1,BA2)，網(wǎng)格數(shù)(百萬(wàn))和拓?fù)湟?jiàn)表1。

表1 計(jì)算網(wǎng)格Table 1 Grids for calculation

圖10為兩種不同背景網(wǎng)格時(shí)的網(wǎng)格示意圖。

表2為計(jì)算得到的不同背景網(wǎng)格、不同的線(xiàn)程數(shù)時(shí)，單次流場(chǎng)內(nèi)迭代計(jì)算情況下，網(wǎng)格裝配時(shí)間與計(jì)算時(shí)間的比率。其中“hc”表示洞切割，ds表示貢獻(xiàn)單元搜索，mpga表示網(wǎng)格裝配過(guò)程中用于處理器之間信息傳遞和比較確認(rèn)的時(shí)間，1 iter表示一次流場(chǎng)內(nèi)迭代。由表可以看出，在網(wǎng)格裝配過(guò)程中，貢獻(xiàn)單元搜索時(shí)間占網(wǎng)格裝配時(shí)間的較大比重，洞切割時(shí)間

與貢獻(xiàn)單元搜索時(shí)間在同一量級(jí)，但比重相對(duì)較小。而對(duì)于本文這種顯式構(gòu)建插值邊界情況，由于待搜索單元數(shù)目較少，處理器之間數(shù)據(jù)傳遞的用時(shí)很少，在目前測(cè)試的算例中幾乎可以忽略。在該測(cè)試中，進(jìn)行一次流場(chǎng)內(nèi)迭代情況下，網(wǎng)格裝配占整體時(shí)間的最大比例約為6%。對(duì)于采用雙時(shí)間方法計(jì)算的旋翼流場(chǎng)，網(wǎng)格裝配時(shí)間對(duì)全局計(jì)算時(shí)間的占比取決于內(nèi)迭代次數(shù)的選取，可見(jiàn)對(duì)于本文算例，內(nèi)迭代大于7次的情況下，網(wǎng)格裝配的時(shí)間則下降到流場(chǎng)計(jì)算時(shí)間的1%以下。由圖還可以看出，隨著線(xiàn)程數(shù)的增加，用于貢獻(xiàn)單元搜索的時(shí)間占比也在增加(網(wǎng)格裝配部分加速比低于流場(chǎng)計(jì)算加速比)。

(a) background1 (b) background2

表2 單次流場(chǎng)計(jì)算時(shí)各部分用時(shí)比率Table 2 Time consumption comparison for difference processes in calculation with 1 flow field iteration

圖11為兩套網(wǎng)格下的流場(chǎng)計(jì)算和網(wǎng)格裝配時(shí)間的加速比(均假定2進(jìn)程時(shí)的加速比為2)，可以發(fā)現(xiàn)，隨著線(xiàn)程數(shù)的增加，用于流場(chǎng)計(jì)算(圖中標(biāo)注“cal”)的加速比增長(zhǎng)率較好，而用于網(wǎng)格裝配(圖中標(biāo)注“ga”)部分的加速比效果相對(duì)較差?！癝tr”表示六面體背景網(wǎng)格，“Uns”表示四面體背景網(wǎng)格。究其原因是在進(jìn)行多線(xiàn)程計(jì)算時(shí)，網(wǎng)格切分過(guò)程保證了各個(gè)線(xiàn)程上擁有近似相同的計(jì)算單元數(shù)，但對(duì)于每個(gè)處理器，其進(jìn)行貢獻(xiàn)單元搜索的單元數(shù)并不平均且隨時(shí)間變化，圖12給出了8線(xiàn)程情況下各線(xiàn)程上面待進(jìn)行貢獻(xiàn)單元搜索單元的個(gè)數(shù)(背景網(wǎng)格1)，貢獻(xiàn)單元搜索完成時(shí)間決定于進(jìn)行貢獻(xiàn)單元搜索耗時(shí)最長(zhǎng)的線(xiàn)程的執(zhí)行時(shí)間。

圖11 計(jì)算和網(wǎng)格裝配的加速比Fig.11 Acceleration ratio of calculation and grid assembly

圖12 一旋翼旋轉(zhuǎn)周期內(nèi)各處理器上待搜索貢獻(xiàn)單元數(shù)Fig.12 Interpolation point numbers in different threads in 1 rotor revolution

圖13為計(jì)算得到的各個(gè)槳葉剖面壓強(qiáng)系數(shù)及與試驗(yàn)值的對(duì)比，在多數(shù)截面上，計(jì)算與試驗(yàn)的吻合度較好(r/R=0.68處有一定差異)，且采用兩套背景網(wǎng)格計(jì)算得到的結(jié)果實(shí)際差異很小，圖14為計(jì)算得到的特征截面渦量等值圖。

(a) background1

(b) background2

(a) background1 (b) background2

3.3 自適應(yīng)輔助結(jié)構(gòu)網(wǎng)格貢獻(xiàn)單元搜索影響

為討論自適應(yīng)輔助結(jié)構(gòu)網(wǎng)格定位貢獻(xiàn)單元初始位置對(duì)網(wǎng)格裝配過(guò)程的影響，這里采用4片槳葉的“HELISHAPE 7A”旋翼為計(jì)算模型，計(jì)算狀態(tài)是槳尖懸停馬赫數(shù)0.616，前進(jìn)比0.167，揮舞變距規(guī)律見(jiàn)文獻(xiàn)[29]。計(jì)算網(wǎng)格采用六面體網(wǎng)格，背景網(wǎng)格數(shù)約為816萬(wàn)，每片槳葉的小網(wǎng)格數(shù)約為228萬(wàn)。

為驗(yàn)證可自適應(yīng)輔助網(wǎng)格的計(jì)算加速效果，采用了不同的輔助結(jié)構(gòu)網(wǎng)格分辨率進(jìn)行了計(jì)算效率的測(cè)試，計(jì)算分別采用4個(gè)和8個(gè)線(xiàn)程并行進(jìn)行。采用4個(gè)線(xiàn)程時(shí)，每個(gè)線(xiàn)程上，背景網(wǎng)格單元約為204萬(wàn)，每塊槳葉網(wǎng)格約為57萬(wàn)。輔助結(jié)構(gòu)網(wǎng)格分辨率Nres按照以下公式確定:

其中：coefficient為系數(shù)，dimen是維度，nelt為計(jì)算單元數(shù)。計(jì)算采用輔助網(wǎng)格的分辨率(每個(gè)方向)近似值見(jiàn)表3(背景網(wǎng)格為BA，槳葉網(wǎng)格為BL，T代表線(xiàn)程數(shù))，計(jì)算時(shí)，限定最大的自適應(yīng)分辨率為3，則可以把計(jì)算內(nèi)存波動(dòng)影響減到很小。

表3 輔助結(jié)構(gòu)網(wǎng)格分辨率Table 3 SAM resolution of grid blocks

圖15為測(cè)試得到的不同分辨率時(shí)，采用與不采用自適應(yīng)輔助結(jié)構(gòu)網(wǎng)格方法情況下，從輔助結(jié)構(gòu)網(wǎng)格生成到貢獻(xiàn)單元搜索確定的時(shí)間長(zhǎng)度比較，圖中時(shí)間按照t0(采用網(wǎng)格分辨率系數(shù)0.1時(shí)的網(wǎng)格裝配時(shí)間)正則化處理。可以看出，隨著輔助結(jié)構(gòu)網(wǎng)格分辨率增加，網(wǎng)格裝配時(shí)間迅速縮短，過(guò)大的分辨率不會(huì)獲得計(jì)算時(shí)間優(yōu)勢(shì)。在采用系數(shù)0.2時(shí)，采用自適應(yīng)方法得到的網(wǎng)格裝配時(shí)間已經(jīng)接近網(wǎng)格裝配速度的極值，這表明采用自適應(yīng)輔助結(jié)構(gòu)網(wǎng)格時(shí)，可以采用更小的網(wǎng)格分辨率，獲得較高的網(wǎng)格裝配速度。對(duì)三維情況，內(nèi)存占用與分辨率三次方成比例，這表明采用自適應(yīng)輔助結(jié)構(gòu)網(wǎng)格方法可以帶來(lái)較大的內(nèi)存節(jié)約。同時(shí)也可以看出，在采用較粗的輔助網(wǎng)格時(shí)(如采用網(wǎng)格分辨率系數(shù)0.1的輔助網(wǎng)格)，改進(jìn)的方法可以提供近似十倍于與原方法的加速效果。這表明當(dāng)前方法在低內(nèi)存開(kāi)銷(xiāo)情況下，可以用來(lái)大幅度提高貢獻(xiàn)單元搜索速度。

(a)

(b)

圖16為采用系數(shù)0.2時(shí)的特征截面上輔助結(jié)構(gòu)和自適應(yīng)輔助結(jié)構(gòu)網(wǎng)格圖。

圖16 輔助結(jié)構(gòu)與自適應(yīng)輔助結(jié)構(gòu)網(wǎng)格示意圖Fig.16 Schematic of SAM and ASAM

圖17為計(jì)算得到的不同槳葉截面在方位角為90°時(shí)，計(jì)算的表面壓強(qiáng)與試驗(yàn)值對(duì)比。計(jì)算與試驗(yàn)的吻合程度較好。

圖17 槳葉各截面壓強(qiáng)(方位角90°)Fig.17 Sectional Cp of blade surface (azimuth=90°)

4 結(jié) 論

對(duì)用于旋翼流場(chǎng)計(jì)算的并行網(wǎng)格裝配方法進(jìn)行了改進(jìn)研究，并進(jìn)行了計(jì)算測(cè)試，得到結(jié)論如下：

1) 壓縮存儲(chǔ)的并行輔助結(jié)構(gòu)網(wǎng)格洞切割方法有效降低了計(jì)算內(nèi)存消耗，內(nèi)存消耗為輔助網(wǎng)格分辨率(維度-1)次方量級(jí)，使得采用更加致密的洞切割輔助結(jié)構(gòu)網(wǎng)格成為可能，計(jì)算方法對(duì)凹形邊界仍然有效，具有一般性。

2) 基于輔助結(jié)構(gòu)網(wǎng)格、N2N及局部窮舉法，構(gòu)建了新的無(wú)邊界交叉幾何判斷的貢獻(xiàn)單元搜索方法，采用輔助結(jié)構(gòu)網(wǎng)格自適應(yīng)技術(shù)，有效縮減了貢獻(xiàn)單元搜索的路徑長(zhǎng)度。在當(dāng)前測(cè)試情況下，采用較粗的輔助網(wǎng)格時(shí)改進(jìn)的方法(ASAM)可以提供近似十倍于原方法(SAM)的加速效果。這表明當(dāng)前方法在低內(nèi)存開(kāi)銷(xiāo)情況下，可以用來(lái)大幅度提高貢獻(xiàn)單元搜索速度。

3) 實(shí)現(xiàn)了完全并行化的網(wǎng)格裝配過(guò)程，并行網(wǎng)格裝配方法理論上不受CPU數(shù)的限制，盡管網(wǎng)格裝配部分的加速效果低于流場(chǎng)部分并行加速效果，但當(dāng)前計(jì)算測(cè)試情況下，方法仍有效地縮減了網(wǎng)格裝配時(shí)間，使得在多進(jìn)程并行的情況下，網(wǎng)格裝配時(shí)間占單次流場(chǎng)計(jì)算迭代時(shí)總時(shí)間長(zhǎng)度的6%以下，在使用雙時(shí)間方法且內(nèi)迭代大于7次時(shí)，網(wǎng)格裝配的時(shí)間則下降到流場(chǎng)計(jì)算時(shí)間的1%以下。對(duì)于需要多個(gè)物理時(shí)間迭代的旋翼非定常流場(chǎng)計(jì)算來(lái)說(shuō)，方法具備較好的計(jì)算效果。