郭永恒，江雄，肖中云，王子維，盧風(fēng)順

中國空氣動(dòng)力研究與發(fā)展中心 計(jì)算空氣動(dòng)力研究所，綿陽 621000

隨著大規(guī)模并行計(jì)算技術(shù)的不斷進(jìn)步，計(jì)算流體力學(xué)(CFD)能夠處理的非定常數(shù)值模擬問題更加豐富[1]；對(duì)多體相對(duì)運(yùn)動(dòng)的情形，比如共軸雙旋翼直升機(jī)的氣動(dòng)干擾[2-3]、超聲速子母彈的播撒[4]、運(yùn)載火箭低空大動(dòng)壓頭罩分離[5]、水陸兩棲飛機(jī)靜水滑行[6]、外掛式導(dǎo)彈的發(fā)射[7]等，單純的結(jié)構(gòu)或非結(jié)構(gòu)網(wǎng)格已不能滿足分析需求。憑借較高的靈活性，重疊網(wǎng)格裝配已發(fā)展成為處理上述問題的關(guān)鍵技術(shù)之一。

重疊網(wǎng)格裝配有顯式與隱式兩類方法。在網(wǎng)格生成階段，顯式方法[8]需要在固壁附近構(gòu)造出挖洞邊界以輔助劃分單元類型，同時(shí)還要人工指定多組部件網(wǎng)格的重疊關(guān)系矩陣，從而導(dǎo)致裝配自動(dòng)化水平偏低，難以適應(yīng)網(wǎng)格運(yùn)動(dòng)時(shí)的分析需求。相比之下，隱式方法摒棄了人工洞邊界，結(jié)合網(wǎng)格拓?fù)浣Y(jié)構(gòu)以及網(wǎng)格組之間的度量關(guān)系進(jìn)行幾何計(jì)算，能夠極大地提升重疊裝配的自動(dòng)化水平。然而，在隱式裝配方法并行化過程中，當(dāng)網(wǎng)格規(guī)模和重疊區(qū)域增大時(shí)，算法涉及的查詢運(yùn)算量和進(jìn)程間的通信壓力急劇攀升，成為制約隱式裝配技術(shù)工程應(yīng)用的瓶頸問題[9]。

為了提高重疊網(wǎng)格隱式裝配效率、拓展其工程應(yīng)用范圍，大量專業(yè)化的重疊分析軟件得以發(fā)展起來。Overflow[10]與Beggar[11]是基于重疊網(wǎng)格的流場模擬軟件，已在工程實(shí)踐中得到廣泛應(yīng)用，但其重疊分析模塊與流場求解器緊密耦合，難以被第三方求解器使用。PEGASUS[12]能夠并行裝配運(yùn)動(dòng)部件網(wǎng)格，但其分析結(jié)果只能通過計(jì)算機(jī)文件系統(tǒng)進(jìn)行傳輸，且大量參數(shù)的人工控制導(dǎo)致其自動(dòng)化水平較低。SUGGAR裝配軟件與插值工具庫DiRTlib結(jié)合[13]能夠很好地解決結(jié)構(gòu)/非結(jié)構(gòu)重疊網(wǎng)格的裝配問題，并且可以嵌入到多種流場解算器中；但進(jìn)行復(fù)雜構(gòu)型網(wǎng)格挖洞時(shí)，控制參數(shù)需要經(jīng)大量嘗試才能達(dá)到合適的水平。SUGGAR++[14]在自動(dòng)化程度、魯棒性以及效率方面有了較大改進(jìn)，但是其直接挖洞算法無法適應(yīng)幾何構(gòu)型不封閉的情況。重疊網(wǎng)格裝配中間件TIOGA[15]將快速搜索算法和動(dòng)態(tài)負(fù)載平衡技術(shù)相結(jié)合，大幅提升了重疊網(wǎng)格的并行裝配效率；但算法的內(nèi)在屬性導(dǎo)致隱式挖洞邊界光滑性欠佳，進(jìn)而影響流動(dòng)模擬品質(zhì)的提升。

綜合來看，常見的重疊網(wǎng)格隱式裝配方法皆起源于Lee和Baeder提出的壁面距離準(zhǔn)則[16-17]。在其并行化實(shí)現(xiàn)過程中，當(dāng)前部件網(wǎng)格節(jié)點(diǎn)到自身所屬壁面距離的計(jì)算、當(dāng)前部件網(wǎng)格節(jié)點(diǎn)到其他部件壁面距離的插值以及貢獻(xiàn)單元的存在性判斷共同構(gòu)成了隱式裝配的自動(dòng)挖洞過程。正是上述3個(gè)步驟的高度耦合使幾何分析過程呈現(xiàn)出復(fù)雜和低效性質(zhì)，初始階段節(jié)點(diǎn)無差別查找操作影響并行裝配效率。

針對(duì)上述問題，本文根據(jù)壁面距離準(zhǔn)則，構(gòu)造出一種高度自動(dòng)化的重疊網(wǎng)格隱式裝配方法，其基本思路是：首先，結(jié)合協(xié)方差分析、切割盒子等快速算法，將壁面距離計(jì)算與貢獻(xiàn)單元存在性判斷解耦，實(shí)現(xiàn)網(wǎng)格組動(dòng)態(tài)重疊關(guān)系的自動(dòng)化識(shí)別；其次，結(jié)合集合分析，設(shè)計(jì)出并行化的自動(dòng)挖洞算法；最后，通過快速查詢方法建立重疊單元與貢獻(xiàn)單元的插值關(guān)系。其中，第1節(jié)使用集合論語言對(duì)重疊網(wǎng)格隱式裝配算法進(jìn)行描述；第2節(jié)詳細(xì)介紹算法的并行實(shí)現(xiàn)細(xì)節(jié)，包括提升效率和自動(dòng)化水平的各種關(guān)鍵技術(shù)；第3節(jié)使用多部件重疊的網(wǎng)格模型測試算法的并行效率，并結(jié)合流場解算器測試重疊插值結(jié)果的準(zhǔn)確性；最后給出結(jié)論與展望。

1 重疊網(wǎng)格隱式裝配算法的集合論表示

現(xiàn)有的重疊網(wǎng)格隱式裝配方法都是依據(jù)“網(wǎng)格節(jié)點(diǎn)到固體壁面的距離”進(jìn)行構(gòu)造的[18]，但是在實(shí)現(xiàn)過程中，流場計(jì)算單元、重疊插值單元以及洞內(nèi)單元的劃分準(zhǔn)則具有很大的差異；例如，很多研究者[18-20]根據(jù)網(wǎng)格圖形范例進(jìn)行算法描述。對(duì)于復(fù)雜幾何構(gòu)型、多部件網(wǎng)格重疊、部件網(wǎng)格外緣邊界與挖洞邊界相交等復(fù)雜情形，圖示法的嚴(yán)謹(jǐn)性便顯得不夠充分。因此，需要使用集合論語言對(duì)重疊網(wǎng)格隱式裝配原理進(jìn)行嚴(yán)格地描述。

設(shè)整個(gè)求解域上包含M個(gè)部件的集合為CompSet，對(duì)于任意部件compi∈CompSet都存在網(wǎng)格組gpi與之構(gòu)成一一對(duì)應(yīng)關(guān)系，其中i∈{1,2,…,M}。設(shè)網(wǎng)格組gpi的節(jié)點(diǎn)集合、面單元集合、體單元集合依次為NodeSeti、FaceElemSeti及VolElemSeti，那么根據(jù)計(jì)算流體力學(xué)通用符號(hào)系統(tǒng)(CFD General Notation System, CGNS)定義的拓?fù)浣Y(jié)構(gòu)，三者之間建立起空間鏈接關(guān)系；設(shè)固壁邊界面上的網(wǎng)格節(jié)點(diǎn)集合為SolidNodeSeti，它與NodeSeti的關(guān)系為

SolidNodeSeti?NodeSeti

(1)

1.1 網(wǎng)格組節(jié)點(diǎn)集合劃分

任取一個(gè)網(wǎng)格節(jié)點(diǎn)P∈NodeSeti，設(shè)它到第j個(gè)部件壁面的距離為

d(P,SolidNodeSetj):=

MinVal{d(P,Q)|Q∈SolidNodeSetj}

(2)

特別地，當(dāng)i=j時(shí)，d(P,SolidNodeSetj)是節(jié)點(diǎn)到自身網(wǎng)格組壁面的距離。如果存在一個(gè)網(wǎng)格組gpj(j≠i)同時(shí)滿足如下條件：

(3)

每個(gè)網(wǎng)格組節(jié)點(diǎn)集合NodeSeti都被劃分成活躍節(jié)點(diǎn)集合PosNodeSeti與非活躍節(jié)點(diǎn)集合NegNodeSeti，即

NodeSeti=PosNodeSeti∪NegNodeSeti

(4)

當(dāng)網(wǎng)格節(jié)點(diǎn)集合完成式(4)的分類后，便得到重疊網(wǎng)格的挖洞邊界。通過拓展挖洞邊界處的體單元，即可實(shí)現(xiàn)洞內(nèi)單元、重疊單元以及場單元的自動(dòng)識(shí)別。

1.2 網(wǎng)格組單元?jiǎng)澐?/h3>

每個(gè)網(wǎng)格組gpi的單元?jiǎng)澐职?個(gè)步驟(在不引起歧義的情況下，文中的集合符號(hào)均省略下標(biāo)i)：

1) 在部件網(wǎng)格組gp上，由于與延拓型物理邊界面相鄰的體單元不可能在計(jì)算過程中扮演場單元的角色，進(jìn)而也不可能成為其他網(wǎng)格組上重疊單元對(duì)應(yīng)的貢獻(xiàn)單元，因此VolElemSet可劃分為一般單元集合GeneralElemSet與延拓型單元集合ExtendElemSet，即

VolElemSet=GeneralElemSet∪ExtendElemSet

(5)

2) 如果體單元velem∈VolElemSet的所有節(jié)點(diǎn)皆為非活躍點(diǎn)集合NegNodeSet中的元素，那么定義velem為非活躍單元。如果體單元

velem∈GeneralElemSet=VolElemSet∩

(6)

的所有節(jié)點(diǎn)皆為活躍點(diǎn)集合PosNodeSet中的元素，那么定義velem為活躍單元，否則定義為過渡單元。不妨記活躍單元、非活躍單元以及過渡單元的集合分別為PosElemSet、NegElemSet和NeuElemSet，將

(7)

定義為延拓型重疊單元的集合。易知

VolElemSet=PosElemSet∪NegElemSet∪

NeuElemSet∪ExtendOvsSet

(8)

其中：等號(hào)右側(cè)4個(gè)元素之間互不相交。

式(7)表明，經(jīng)過壁面距離準(zhǔn)則篩查后，如果延拓型單元不屬于非活躍單元集合，那么它將被識(shí)別為重疊單元。

3) 任取體單元neg∈NegElemSet，如果存在體單元neu∈NeuElemSet，使

NodeList(neg)∩NodeList(neu)≠?

(9)

那么將neg歸入一般性重疊單元集合GeneralOvsSet，否則將其歸入洞內(nèi)單元集合HoleElemSet。NodeList(neg)和NodeList(neu)分別表示體單元neg和neu的節(jié)點(diǎn)元素集合。此處單元分類的幾何意義是：在洞邊界附近，與過渡單元相鄰的非活躍單元轉(zhuǎn)化為重疊單元；與過渡單元不相鄰的非活躍單元轉(zhuǎn)化為洞內(nèi)單元。

4) 場單元集合的定義為

FieldElemSet:=PosElemSet∪NeuElemSet

(10)

需要說明的是，過渡單元是同時(shí)具有活躍節(jié)點(diǎn)和非活躍節(jié)點(diǎn)的體單元，因?yàn)樗『每缭街丿B挖洞邊界，所以將與另一網(wǎng)格組上的過渡單元交織在一起；為避免不同網(wǎng)格組上的重疊單元互為對(duì)方的貢獻(xiàn)單元，在式(10)中，過渡單元最終轉(zhuǎn)化為場單元。重疊單元集合定義為

OvsElemSet:=GeneralOvsSet∪ExtendOvsSet

(11)

至此，結(jié)合第3)步的洞內(nèi)單元HoleElemSet，最終完成了單元分類

VolElemSet=FieldElemSet∪HoleElemSet∪

OvsElemSet

(12)

從本文重疊網(wǎng)格裝配算法的構(gòu)造過程可以看出，與圖示法相比，基于集合論語言的描述方式具有如下4個(gè)方面的優(yōu)勢：

1) 具有普適性，不受網(wǎng)格幾何形狀復(fù)雜程度或面、體單元類型多樣性的影響。

2) 能夠清晰表述多組網(wǎng)格重疊時(shí)的復(fù)雜映射關(guān)系。

3) 能夠精確判斷延拓型外緣邊界與距離挖洞邊界相交等復(fù)雜情形。

4) 可作為算法分析的有力工具，輔助提升并行裝配效率。

2 自動(dòng)化并行重疊網(wǎng)格隱式裝配方法

圖1顯示重疊網(wǎng)格的并行裝配流程，包含幾何信息輸入、網(wǎng)格組自動(dòng)配對(duì)、重疊單元隱式裝配、插值信息輸出等4個(gè)步驟。

圖1 重疊網(wǎng)格并行化裝配流程

步驟1左側(cè)(藍(lán)色)部分是算法主線，依次為分區(qū)網(wǎng)格幾何信息的輸入、部件網(wǎng)格組之間重疊關(guān)系的自動(dòng)識(shí)別、重疊單元的隱式裝配以及插值信息的輸出。

步驟2右上(綠色)部分是幾何數(shù)據(jù)的基本構(gòu)成，包含邊界映射信息、部件網(wǎng)格進(jìn)程分組信息、分區(qū)網(wǎng)格幾何數(shù)據(jù)以及網(wǎng)格塊對(duì)接信息。其中，邊界映射信息主要用來識(shí)別固體邊界面以及部件網(wǎng)格的外延拓邊界面；在挖洞過程完成后，網(wǎng)格塊對(duì)接關(guān)系用來跨進(jìn)程識(shí)別重疊單元。

步驟3右中(黃色)部分是重疊裝配所需的主要數(shù)據(jù)結(jié)構(gòu)，它們根據(jù)原始幾何信息構(gòu)建。其中，部件網(wǎng)格極小內(nèi)部盒子、外部盒子用來輔助網(wǎng)格組自動(dòng)配對(duì)；覆蓋固壁節(jié)點(diǎn)的遞歸盒子用來對(duì)查詢節(jié)點(diǎn)進(jìn)行初步篩查；空間體單元的二叉樹結(jié)構(gòu)用于單元節(jié)點(diǎn)或格心貢獻(xiàn)單元的快速查詢。

步驟4右下(橙色)部分顯示具體的重疊單元隱式裝配過程，包括并行計(jì)算節(jié)點(diǎn)到自身部件和鄰居部件距離、初篩待查詢節(jié)點(diǎn)、確認(rèn)節(jié)點(diǎn)激活屬性、識(shí)別挖洞邊界、識(shí)別重疊單元以及查找貢獻(xiàn)單元等7個(gè)步驟，(綠色)虛線條展示了各步驟所需的數(shù)據(jù)結(jié)構(gòu)。

2.1 重疊網(wǎng)格組并行分區(qū)策略

重疊網(wǎng)格并行裝配算法是與網(wǎng)格組分區(qū)策略緊密結(jié)合在一起的。假設(shè)存在N個(gè)進(jìn)程參與流場計(jì)算，那么上述M個(gè)網(wǎng)格組存在兩種分區(qū)策略。

策略T1把全體網(wǎng)格組包含的網(wǎng)格塊視為一個(gè)整體，根據(jù)體單元數(shù)負(fù)載均衡原則，使用Greedy算法或METIS算法進(jìn)行統(tǒng)一分配，如圖2(a)所示。

策略T2結(jié)合網(wǎng)格組體單元數(shù)在整體網(wǎng)格中所占的比重，計(jì)算其所需的進(jìn)程數(shù)，然后將當(dāng)前網(wǎng)格組均勻分配到相應(yīng)的進(jìn)程組中，如圖2(b)所示。

圖2 兩種剖分策略對(duì)比

文獻(xiàn)[21]中的重疊網(wǎng)格裝配方法使用T1策略能夠保證體單元分配結(jié)果的高度均衡性，但也容易使每個(gè)網(wǎng)格組上形成過多的碎片化網(wǎng)格塊，從而影響并行計(jì)算效率。在重疊網(wǎng)格并行化裝配算法中，兩個(gè)網(wǎng)格組涉及的進(jìn)程組間需要頻繁進(jìn)行數(shù)據(jù)通信；此時(shí)使用較少的進(jìn)程數(shù)容納當(dāng)前網(wǎng)格組的各個(gè)網(wǎng)格塊，能夠有效提高算法的并行效率。因此，采用T2策略進(jìn)行網(wǎng)格分區(qū)。

2.2 網(wǎng)格組自動(dòng)配對(duì)技術(shù)

在傳統(tǒng)的重疊網(wǎng)格裝配算法中，需要在圖形界面下人工指定網(wǎng)格組間的重疊關(guān)系，比如使用邏輯矩陣來描述。對(duì)于多體相對(duì)運(yùn)動(dòng)的情形，各網(wǎng)格組之間的相對(duì)位置與重疊關(guān)系通常是動(dòng)態(tài)變化的，這些特點(diǎn)決定了人工配對(duì)操作難以滿足重疊網(wǎng)格自動(dòng)化裝配的需求。常興華等[22]為提升大規(guī)模重疊網(wǎng)格的裝配效率與自動(dòng)化水平，使用包含每個(gè)網(wǎng)格塊的長方體輔助判斷重疊關(guān)系，其中長方體表面分別與坐標(biāo)軸垂直。這種方法本質(zhì)上是使用大量的簡單長方體實(shí)現(xiàn)網(wǎng)格組區(qū)域的幾何逼近，將網(wǎng)格組重疊關(guān)系判斷問題轉(zhuǎn)化為相應(yīng)長方體集合的求交邏輯運(yùn)算；但是，這種方法導(dǎo)致過于碎片化的分區(qū)結(jié)果，根據(jù)Navier-Stokes (N-S) 方程的數(shù)值迭代理論，它將使收斂歷程大大延長，從而降低流場求解的效率。受文獻(xiàn)[15]中方向適應(yīng)型邊界盒子的啟發(fā)，本文建立了關(guān)于部件表面網(wǎng)格點(diǎn)坐標(biāo)的協(xié)方差矩陣，結(jié)合譜分析技術(shù)，構(gòu)造了部件網(wǎng)格組的極小內(nèi)部盒子與極小外部盒子，以此輔助網(wǎng)格組區(qū)域重疊關(guān)系的快速識(shí)別。

在第i個(gè)部件compi∈CompSet上，設(shè)固壁節(jié)點(diǎn)集合SolidNodeSeti的元素個(gè)數(shù)為ms，那么可以把所有節(jié)點(diǎn)的坐標(biāo)依次存儲(chǔ)于3個(gè)ms維向量中：

(13)

式中：(xms,yms,zms)為編號(hào)s的固壁節(jié)點(diǎn)空間坐標(biāo)。協(xié)方差矩陣為

(14)

根據(jù)協(xié)方差運(yùn)算的可交換性可知：D為對(duì)稱矩陣，其特征系統(tǒng)中3個(gè)特征向量兩兩正交。通過Householder變換可將其約化成三對(duì)角陣，進(jìn)一步使用QR方法即可得到D的全部特征值及其特征向量。設(shè)特征向量的規(guī)范化結(jié)果為

ns=[As,Bs,Cs]Ts=0,1,2

(15)

給定方向序號(hào)s，以ns為基準(zhǔn)法向量建立空間平面參數(shù)方程：

Asx+Bsy+Csz+λs=0

(16)

依次遍歷SolidNodeSeti中的元素并取參數(shù)λs的上下確界，即得方向部件的2個(gè)臨界方程。3個(gè)特征方向的臨界方程使部件表面網(wǎng)格包含在一個(gè)極小長方體內(nèi)，稱其為極小內(nèi)部盒子MIBi，如圖3紅色線條所示。

如果保持特征方向不變，將上述運(yùn)算中的集合SolidNodeSeti替換成整體網(wǎng)格組節(jié)點(diǎn)集合NodeSeti，那么整體網(wǎng)格組將包含在新的極小長方體內(nèi)，稱其為極小外部盒子MEBi，如圖3藍(lán)色線條所示。對(duì)網(wǎng)格生成有如下附加要求：極小內(nèi)部盒子覆蓋的空間區(qū)域包含在整個(gè)網(wǎng)格組區(qū)域中，這在網(wǎng)格生成過程中是容易實(shí)現(xiàn)的；同時(shí)，網(wǎng)格組外緣表面適當(dāng)遠(yuǎn)離邊界層網(wǎng)格也是提升重疊插值精度的要求。圖3展示的是螺旋槳單葉片網(wǎng)格上極小內(nèi)部盒子與極小外部盒子的包含關(guān)系，而圖4展示了整個(gè)螺旋槳各部件網(wǎng)格相互重疊關(guān)系。

圖3 螺旋槳單葉片極小內(nèi)部/外部盒子

圖4 螺旋槳葉片各部件網(wǎng)格重疊關(guān)系

任意給定網(wǎng)格組gpi，設(shè)其覆蓋區(qū)域?yàn)镈i，那么存在集合關(guān)系式

MIBi?Di?MEBi

(17)

對(duì)于當(dāng)前網(wǎng)格組gpi(稱為源網(wǎng)格組SG)，如果對(duì)另一組網(wǎng)格gpj，關(guān)系式：

MEBi∩MEBj=?

(18)

成立，那么

Di∩Dj=?

(19)

即當(dāng)前狀態(tài)下gpi與gpj不存在相互重疊關(guān)系；否則，稱gpj為與gpi存在可能重疊關(guān)系的目標(biāo)網(wǎng)格組TG。式(18)的邏輯判斷可以根據(jù)圖5所示的分離平面法快速完成：依次遍歷源網(wǎng)格組SG和目標(biāo)網(wǎng)格組TG的3個(gè)特征方向，那么該式成立的充要條件為：存在以特征方向?yàn)榉ㄏ蛄康目臻g平面SP，使MEBi與MEBj分居其兩側(cè)。圖6描述了源網(wǎng)格組與目標(biāo)網(wǎng)格組(為目標(biāo)網(wǎng)格組的全局編號(hào))之間的重疊關(guān)系，其中，與當(dāng)前源網(wǎng)格組SG發(fā)生重疊關(guān)系的目標(biāo)網(wǎng)格組TG的總數(shù)為L；在并行計(jì)算環(huán)境下，源網(wǎng)格組的各個(gè)子塊grid_s分布在編號(hào)連續(xù)的NG個(gè)進(jìn)程中，而它的某個(gè)目標(biāo)網(wǎng)格組(以TG[2]為例)的各個(gè)子塊grid_t分布在編號(hào)連續(xù)的MG個(gè)進(jìn)程中。

圖5 空間分離平面法

圖6 源網(wǎng)格組與目標(biāo)網(wǎng)格組的重疊關(guān)系

使用圖6所示的映射關(guān)系，可以對(duì)多組網(wǎng)格間的復(fù)雜重疊關(guān)系進(jìn)行分解。從整體上來說，優(yōu)先遍歷當(dāng)前的源網(wǎng)格組SG；對(duì)于特定SG，局部遍歷與之重疊的多個(gè)目標(biāo)網(wǎng)格組TG。

2.3 自動(dòng)挖洞的并行化實(shí)現(xiàn)

在第1節(jié)中，描述重疊網(wǎng)格自動(dòng)挖洞的壁面距離準(zhǔn)則包含兩個(gè)方面內(nèi)容：一是網(wǎng)格節(jié)點(diǎn)壁面距離的計(jì)算；二是網(wǎng)格節(jié)點(diǎn)與體單元之間的覆蓋關(guān)系。在當(dāng)前網(wǎng)格組gpi上，如果直接遍歷集合NodeSeti并對(duì)每個(gè)網(wǎng)格點(diǎn)元素P∈NodeSeti同時(shí)執(zhí)行上述分析，那么不僅會(huì)使程序編寫難度急劇提高，而且會(huì)使源網(wǎng)格組與目標(biāo)網(wǎng)格組之間的通訊壓力迅速增長，進(jìn)而限制該算法在大規(guī)模并行計(jì)算中的應(yīng)用。反之，如果首先對(duì)NodeSeti中所有元素基于距離參數(shù)比較法計(jì)算距離參數(shù)，就可以大大降低判斷節(jié)點(diǎn)與體單元覆蓋關(guān)系時(shí)的計(jì)算量與通訊壓力，其原理描述如下。

對(duì)于?gpi所對(duì)應(yīng)的目標(biāo)網(wǎng)格組gpj(j=j0,j1,…,jL-1)，如果?P∈NodeSeti，條件C1：

d(P,SolidNodeSetj)>d(P,SolidNodeSeti)

(20)

成立，那么P∈PosNodeSeti。作為式(3)的互斥條件之一，式(20)使用單純的距離判據(jù)排除了大量待查詢的節(jié)點(diǎn)。

此外，結(jié)合2.2節(jié)關(guān)于極小外部盒子的定義可知，如果?P∈NodeSeti，條件C2：

(21)

成立，那么可以推出

(22)

作為式(3)的另一個(gè)互斥條件，式(22)根據(jù)網(wǎng)格組的極小外部盒子的空間相對(duì)位置進(jìn)行間接篩查，同樣可知P∈PosNodeSeti。C1與C2稱為距離初篩條件與覆蓋初篩條件，據(jù)此可提升重疊網(wǎng)格并行化挖洞算法的效率。

圖7描述了當(dāng)前源網(wǎng)格組SG壁面幾何信息的收集以及與各目標(biāo)網(wǎng)格組TG共享的過程：非阻塞通訊模式下，在SG包含的各個(gè)進(jìn)程組內(nèi)，對(duì)相應(yīng)網(wǎng)格塊grid_si(i=0,1,…,NG-1，NG為分區(qū)后SG包含的網(wǎng)格塊總數(shù))上的壁面節(jié)點(diǎn)坐標(biāo)序列進(jìn)行數(shù)據(jù)壓碼操作，形成NG個(gè)二進(jìn)制數(shù)據(jù)流solid_s，在主進(jìn)程中通過數(shù)據(jù)容器solid_node_collector加以匯總后，分別在SG與TGj(j=0,1,…,L-1，L為SG對(duì)應(yīng)的目標(biāo)網(wǎng)格組個(gè)數(shù))所在的進(jìn)程組內(nèi)共享并解碼，形成完整的SolidNodeSet信息后，對(duì)節(jié)點(diǎn)壁面距離信息進(jìn)行更新。

圖7 源網(wǎng)格組壁面幾何信息收集與共享

圖8描述了節(jié)點(diǎn)壁面距離參數(shù)的定義：對(duì)源網(wǎng)格組SG(gpi)上任意給定網(wǎng)格節(jié)點(diǎn)P∈NodeSeti，令ds(P)∶=d(P,SolidNodeSeti)為P到SG壁面的距離，dtj∶=d(P,SolidNodeSetj)為P到TGj壁面的距離，其中j取0～L-1范圍內(nèi)的整數(shù)；令dt(P)∶=MinVal{dt[0],dt[1],…,dt[L-1]}，那么，ds與dt就構(gòu)成了距離比較判據(jù)的兩個(gè)重要參數(shù)。在并行計(jì)算框架下，若直接對(duì)ds與dt進(jìn)行計(jì)算，則存在兩種方式：① 把P的信息發(fā)送L次，分別與SolidNodeSetj計(jì)算出dtj，然后通過全局規(guī)約操作得到dt值；② 把SolidNodeSetj的信息發(fā)送到P所在的進(jìn)程，對(duì)其依次遍歷，按最小值準(zhǔn)則不斷更新dt值。可見，方式①的通訊壓力較大，全局規(guī)約操作將影響計(jì)算效率的提升，因此采用后一種方式計(jì)算ds與dt。

圖8 節(jié)點(diǎn)壁面距離參數(shù)定義

空間節(jié)點(diǎn)的壁面距離是重疊網(wǎng)格隱式裝配方法的關(guān)鍵參數(shù)之一，如果采用線性查找方法，遍歷所有的壁面點(diǎn)后取最小值，那么，當(dāng)壁面節(jié)點(diǎn)總數(shù)很大時(shí)，計(jì)算成本將急劇升高；同時(shí)，這種簡單查找法也難以適應(yīng)網(wǎng)格運(yùn)動(dòng)、網(wǎng)格自適應(yīng)時(shí)壁面距離頻繁更新的情形。相比之下，基于點(diǎn)集樹形排序思想的方盒切割算法[23]能夠很好地處理這一問題，它包括固壁點(diǎn)集臨界長方體(方盒)序列的生成和最近距離點(diǎn)的查找兩部分內(nèi)容，其中，針對(duì)臨界長方體序列的輔助分析，可以在點(diǎn)集遍歷過程中動(dòng)態(tài)排除無關(guān)的子集，從而大幅提升壁面距離的計(jì)算效率。

臨界長方體序列的構(gòu)造分為如下4個(gè)步驟：

步驟1取網(wǎng)格組gpi共享的壁面節(jié)點(diǎn)集合SolidNodeSet[i]，沿著平行于坐標(biāo)平面的方向，根據(jù)式(23)分析壁面節(jié)點(diǎn)坐標(biāo)分量的上下界，構(gòu)成臨界長方體序列根結(jié)點(diǎn)的6個(gè)基本參數(shù)，它們分別是

(23)

步驟2在節(jié)點(diǎn)集坐標(biāo)分量最大的變化區(qū)間上進(jìn)行排序與二分法操作，構(gòu)造當(dāng)前節(jié)點(diǎn)集的左右子集以及相應(yīng)的臨界長方體參數(shù)。分別計(jì)算

(24)

不妨假設(shè)

δz<δy<δx

(25)

然后，對(duì)壁面節(jié)點(diǎn)序列按照x坐標(biāo)大小排序，并采用二分法生成SolidNodeSeti的2個(gè)子集SubNodeSeti0與SubNodeSeti1。

步驟4當(dāng)臨界長方體構(gòu)造完成后，僅收集所有的葉子結(jié)點(diǎn)(不可分的臨界長方體)，形成壁面距離計(jì)算的基礎(chǔ)查詢數(shù)據(jù)集。經(jīng)過若干層級(jí)的分解，原始的SolidNodeSeti可以寫成一系列節(jié)點(diǎn)子集的并集，即

SolidNodeSeti=SubNodeSeti0∪SubNodeSeti1∪

…∪SubNodeSeti(NSub-1)

此時(shí)，對(duì)于每個(gè)SubNodeSetis，分別以空間平面

(26)

構(gòu)成臨界長方體Capsules，以此覆蓋集合SubNodeSetis，如圖9和圖10所示。

圖9 定點(diǎn)和臨界長方體的距離

圖10 壁面距離計(jì)算的臨界長方體算法

任意給定空間網(wǎng)格節(jié)點(diǎn)P，可以按照如下4個(gè)步驟在集合SolidNodeSeti中快速查找出最近的固壁節(jié)點(diǎn)：

步驟1依次計(jì)算當(dāng)前網(wǎng)格節(jié)點(diǎn)到基礎(chǔ)查詢數(shù)據(jù)集中每個(gè)臨界長方體的距離，其定義需要滿足這樣的條件：無論網(wǎng)格點(diǎn)在臨界長方體內(nèi)部還是外部，它與臨界長方體的距離將小于到其中任意固壁節(jié)點(diǎn)的距離。記P到每個(gè)Capsules區(qū)域的距離為d(P,Capsules)。根據(jù)圖9，設(shè)P的坐標(biāo)為(x,y,z)、Capsules的格心坐標(biāo)為(xC,yC,zC)，那么，對(duì)于任意相對(duì)位置，d(P,Capsules)具有統(tǒng)一形式：

d(P,Capsules)=

(27)

式中：

f(ξ,η,ζ):=

(28)

為每個(gè)Capsules對(duì)象設(shè)置邏輯變量Logic表征其激活屬性，初始值為TRUE；同時(shí)，初始化

d(P,SolidNodeSeti)=LARGE_DISTANCE

(29)

式中：LARGE_DISTANCE=1.0×1040。

步驟2在當(dāng)前所有激活的臨界長方體中進(jìn)行遍歷，選取使d(P,Capsules)達(dá)到最小值的Capsules。如果符合條件的Capsules存在，那么記錄其編號(hào)s以及d(P,Capsules)，此時(shí)最近的固壁點(diǎn)位于其中的可能性很大，為進(jìn)一步確認(rèn)，通過步驟3對(duì)臨界長方體內(nèi)的固壁點(diǎn)集進(jìn)行更加詳細(xì)的判斷；否則，在激活狀態(tài)的臨界長方體不存在時(shí)，轉(zhuǎn)向步驟4即可。

步驟3一邊計(jì)算網(wǎng)格節(jié)點(diǎn)與固壁節(jié)點(diǎn)的距離，一邊將其與處于激活狀態(tài)的各個(gè)臨界長方體的距離進(jìn)行比較，如果后者大于當(dāng)前的距離，那么其中的任意固壁節(jié)點(diǎn)皆不可能成為最近距離點(diǎn)，所以不再予以考慮，具體過程描述為：在Capsules覆蓋的子集SubNodeSetis上，計(jì)算

d(P,SubNodeSetis)∶=

MinVal(d(P,Q)|Q∈SubNodeSetis)

(30)

令Capsules的邏輯變量Logic=FALSE且更新d(P,SubNodeSetis)，即

d(P,SolidNodeSeti)=

min(d(P,SolidNodeSeti),d(P,SubNodeSetis))

(31)

遍歷全部臨界長方體，如果臨界長方體Capsules滿足條件

d(P,SolidNodeSeti)

(32)

那么，為其邏輯變量賦值Logic=FALSE；最后，返回步驟2。

步驟4當(dāng)步驟2和3的循環(huán)判斷結(jié)束后，即可得到最終的壁面距離d(P,SolidNodeSeti)。

在網(wǎng)格組gpi上，結(jié)合距離初篩條件與覆蓋初篩條件，定義查詢節(jié)點(diǎn)的集合

QueryNodeSeti:={P|dt(P)

Dj0∪Dj1∪…∪DjL-1}

(33)

且QueryNodeSeti滿足

NegNodeSeti?QueryNodeSeti?NodeSeti

(34)

式中：源網(wǎng)格組gpi與目標(biāo)網(wǎng)格組gpjs極小外部盒子參數(shù)可以在并行編程框架下實(shí)現(xiàn)共享，js∈{j0,j1,…,jT-1}為目標(biāo)網(wǎng)格組的編號(hào)。

對(duì)于集合關(guān)系式(34)，根據(jù)覆蓋關(guān)系式(17)從QueryNodeSeti篩選出NegNodeSeti，在并行編程框架下的查詢算法以及數(shù)據(jù)流如圖11所示。

圖11 節(jié)點(diǎn)查詢算法及插值信息數(shù)據(jù)流

設(shè)當(dāng)前源網(wǎng)格組SG各子網(wǎng)格塊分布在編號(hào)連續(xù)的S個(gè)進(jìn)程中，在每個(gè)進(jìn)程上，子網(wǎng)格塊將按照六維超子空間分割法[24]生成獨(dú)立的二叉樹結(jié)構(gòu)，為任意外部節(jié)點(diǎn)提供快速查詢服務(wù)；同時(shí)，按照式(5)，根據(jù)延拓型物理邊界，SG在每個(gè)進(jìn)程內(nèi)子網(wǎng)格塊體單元被初步劃分為一般單元與延拓單元。設(shè)SG的任意目標(biāo)網(wǎng)格組為TG，因?yàn)椴樵児?jié)點(diǎn)分布在編號(hào)連續(xù)的T個(gè)進(jìn)程中，每個(gè)進(jìn)程對(duì)應(yīng)查詢節(jié)點(diǎn)的子集QueryNodeListjs。在非阻塞通訊模式下，主進(jìn)程對(duì)QueryNodeListjs進(jìn)行收集，即

QueryNodeSetj=QueryNodeListj0∪

QueryNodeListj1∪…∪QueryNodeListjT-1

(35)

然后將其共享到SG所在的進(jìn)程組內(nèi)進(jìn)行查詢操作，其步驟如下：

步驟1貢獻(xiàn)單元的并行查找與插值信息的計(jì)算。在上述共享機(jī)制作用下，SG涉及的每個(gè)進(jìn)程內(nèi)皆包含目標(biāo)網(wǎng)格組TG上查詢節(jié)點(diǎn)集合QueryNodeSetj的完整信息。任意給定查詢節(jié)點(diǎn)對(duì)象P∈QueryNodeSetj，通過局部二叉樹查詢算法，如果存在體單元velem∈GenElemSeti覆蓋P，那么在對(duì)象P中記錄如下3種信息：① velem所在的當(dāng)前進(jìn)程編號(hào)以及局部體單元序號(hào)；②P在velem中的相對(duì)坐標(biāo)；③P到SG壁面的距離d(P,SolidNodeSeti)。由于已計(jì)算出velem各個(gè)節(jié)點(diǎn)到自身網(wǎng)格組的壁面距離參數(shù)ds，因此結(jié)合②中的相對(duì)坐標(biāo)以及距離函數(shù)的空間連續(xù)性，即對(duì)d(P,SolidNodeSeti)進(jìn)行三線性插值。如果當(dāng)前進(jìn)程內(nèi)不存在體單元velem∈GenElemSeti覆蓋P，那么擦除對(duì)象P在當(dāng)前進(jìn)程的全部信息。

步驟2貢獻(xiàn)單元存在時(shí)節(jié)點(diǎn)的篩查。在SG涉及的每個(gè)進(jìn)程內(nèi)，剩余查詢節(jié)點(diǎn)信息通過圖11中的信息收集器collector匯總后，在TG包含的進(jìn)程組內(nèi)進(jìn)行信息共享。在TG各子網(wǎng)格塊所在的進(jìn)程內(nèi)，皆包含collector共享的有效查詢結(jié)果。如果查詢節(jié)點(diǎn)P原始的進(jìn)程編號(hào)與當(dāng)前進(jìn)程一致，那么結(jié)合其局部節(jié)點(diǎn)編號(hào)，在當(dāng)前進(jìn)程網(wǎng)格塊的相應(yīng)位置更新信息：如果dt(P)

步驟4剩余查詢節(jié)點(diǎn)的處理。在式(34)中，查詢節(jié)點(diǎn)集包含了非激活點(diǎn)集，在確認(rèn)階段，排除所有的非激活點(diǎn)集后，剩余的查詢節(jié)點(diǎn)P都需歸入一般節(jié)點(diǎn)集合，即P∈GenNodeSetj。

主要對(duì)4種非結(jié)構(gòu)體單元進(jìn)行處理，即四面體、六面體、三棱柱和金字塔。為確定給定節(jié)點(diǎn)與體單元的相對(duì)位置，需要把物理坐標(biāo)系下的體單元映射成為計(jì)算坐標(biāo)系下的標(biāo)準(zhǔn)計(jì)算單元。圖12展示了4種類型單元與標(biāo)準(zhǔn)單元之間的微分同胚關(guān)系，其中物理坐標(biāo)系下體單元的覆蓋區(qū)域可用形函數(shù)的線性組合加以描述：

圖12 4種類型單元與標(biāo)準(zhǔn)單元之間的微分同胚關(guān)系

(36)

式中：(x,y,z)為體單元內(nèi)部任意一點(diǎn)的物理坐標(biāo)；(ξ,η,ζ)為對(duì)應(yīng)的計(jì)算坐標(biāo)；M為單元角點(diǎn)總數(shù)；(xs,ys,zs)為角點(diǎn)物理坐標(biāo)。不同類型velem單元的形函數(shù)序列如表1所示。

表1 4種類型單元的插值形函數(shù)序列

為了確定給定節(jié)點(diǎn)P(x*,y*,z*)在標(biāo)準(zhǔn)計(jì)算單元中的相對(duì)坐標(biāo)(ξ,η,ζ)，建立非線性方程組

(37)

構(gòu)造Newton迭代序列

(38)

式中：下標(biāo)n與n+1用來區(qū)分當(dāng)前時(shí)間層與更新時(shí)間層的變量。設(shè)收斂結(jié)果(ξ,η,ζ)滿足條件：

那么，可以得到P∈velem的結(jié)論。

至此，已實(shí)現(xiàn)網(wǎng)格節(jié)點(diǎn)分類算法的并行化工作。在上述“初篩+確認(rèn)”模式中，基于兩種距離參數(shù)以及網(wǎng)格組極小外部盒子參數(shù)對(duì)可疑節(jié)點(diǎn)的篩查，使各進(jìn)程內(nèi)查詢節(jié)點(diǎn)集合的規(guī)模得到了有效的控制；同時(shí)，對(duì)于目標(biāo)網(wǎng)格組TG來說，伴隨著洞內(nèi)節(jié)點(diǎn)的不斷剔除，它在與另一組網(wǎng)格組SG進(jìn)行配對(duì)時(shí)，查詢節(jié)點(diǎn)的總量將進(jìn)一步的縮小。這些措施都極大降低了并行通訊的壓力。

2.4 重疊單元的識(shí)別及并行化裝配

根據(jù)第1節(jié)的介紹，在完成網(wǎng)格節(jié)點(diǎn)分類的基礎(chǔ)上，可以結(jié)合單元與節(jié)點(diǎn)的拓?fù)潢P(guān)系，完成場單元、洞內(nèi)單元以及重疊單元的劃分。在2.3節(jié)描述的單元分類過程中，第1)2)4)步可以在各自進(jìn)程內(nèi)通過“節(jié)點(diǎn)-單元染色”方式獨(dú)立進(jìn)行。然而，在第3)步中，與過渡單元相鄰的非活躍單元被定義為一般性重疊單元；如果過渡單元與非活躍單元的分界面與網(wǎng)格分區(qū)邊界面存在重合(如圖13所示)，那么需要通過分區(qū)邊界處點(diǎn)對(duì)點(diǎn)的并行通訊完成單元類型信息的傳遞，其算法流程如圖14所示。

圖13 過渡單元與非活躍單元分界面與網(wǎng)格分區(qū)邊界面重合情況

在圖14所示的各個(gè)進(jìn)程內(nèi)，沿著箭頭指向，重疊單元識(shí)別算法的并行化共分為5個(gè)步驟：① 根據(jù)網(wǎng)格節(jié)點(diǎn)的激活屬性將體單元?jiǎng)澐譃榛钴S體單元、非活躍體單元、過渡體單元與延拓型重疊單元4種類型；② 對(duì)過渡體單元包含的網(wǎng)格節(jié)點(diǎn)(稱為過渡節(jié)點(diǎn))進(jìn)行染色標(biāo)記；③ 在當(dāng)前進(jìn)程內(nèi)，根據(jù)已染色的過渡節(jié)點(diǎn)識(shí)別一般類型重疊單元與洞內(nèi)單元；④ 根據(jù)網(wǎng)格分區(qū)邊界處的并行通訊，識(shí)別一般型重疊單元；⑤把一般型重疊單元與延拓型重疊單元統(tǒng)一歸并到重疊單元集合中，活躍體單元與過渡體單元統(tǒng)一歸并到場單元集合中。

圖14 重疊單元的并行識(shí)別過程

圖15展示了使用二叉樹技術(shù)建立重疊單元格心與貢獻(xiàn)體單元之間“插值-覆蓋”映射的過程。圖15與圖11有較高的相似度，特別是在網(wǎng)格節(jié)點(diǎn)分類(即并行化自動(dòng)挖洞)時(shí)使用的二叉樹結(jié)構(gòu)未發(fā)生變化。此處主要有2點(diǎn)不同：① 目標(biāo)網(wǎng)格組上收集的節(jié)點(diǎn)信息變更為重疊單元格心信息；② 在源網(wǎng)格組上，二叉樹查詢算法對(duì)應(yīng)的有效貢獻(xiàn)單元附加約束屬性由一般型單元變更為場單元。有效貢獻(xiàn)單元的約束條件，有效避免了不同網(wǎng)格組重疊單元互相插值的異?，F(xiàn)象，從而使流場迭代過程順利進(jìn)行。

圖15 重疊單元的并行裝配

3 算例驗(yàn)證

本文采用面向?qū)ο缶幊趟枷耄褂肅++語言和MPI并行編程模型實(shí)現(xiàn)了重疊網(wǎng)格隱式裝配算法的并行化，從而建立了圖16所示的“重疊網(wǎng)格縫紉機(jī)”(Overlapping Grid Sewing Machine, OGSM)工具庫。

圖16 OGSM運(yùn)行流程

OGSM描述的流程共分為8個(gè)節(jié)點(diǎn)，它們分別是：① 讀入當(dāng)前進(jìn)程子區(qū)域的網(wǎng)格拓?fù)湫畔ⅲ虎?建立網(wǎng)格分區(qū)邊界的并行通訊架構(gòu)；③ 自動(dòng)化注冊(cè)網(wǎng)格組重疊關(guān)系；④ 基于壁面距離準(zhǔn)則識(shí)別重疊單元；⑤ 在背景網(wǎng)格上查找重疊單元的最優(yōu)貢獻(xiàn)單元；⑥ 重疊裝配信息的跨進(jìn)程整理；⑦ 統(tǒng) 計(jì)當(dāng)前進(jìn)程內(nèi)子區(qū)域上貢獻(xiàn)單元的總數(shù)；⑧ 輸 出重疊插值以及單元分類信息。其中，左側(cè)(藍(lán)色)方框表示的是OGSM的外部接口，右側(cè)(綠色)方框內(nèi)是重疊網(wǎng)格并行裝配算法的內(nèi)部實(shí)現(xiàn)。

3.1 五球體部件

為了獨(dú)立測試OGSM的并行裝配能力，對(duì)包含5個(gè)球體部件的重疊網(wǎng)格進(jìn)行裝配。如圖17所示，每個(gè)球體部件獨(dú)立生成非結(jié)構(gòu)網(wǎng)格，體單元總量為1.2億。在z=0平面內(nèi)進(jìn)行切割，可以看到5個(gè)區(qū)域之間的自動(dòng)裝配結(jié)果。其中，藍(lán)色曲線為網(wǎng)格分區(qū)邊界，綠色部分為重疊單元。在該算例中，部件網(wǎng)格的外延拓邊界與挖洞邊界相交，球體之間的挖洞邊界遵循了壁面距離準(zhǔn)則，本文的自動(dòng)挖洞邏輯的準(zhǔn)確性得到了檢驗(yàn)。圖18展示了算法的并行效率曲線，隨著分區(qū)進(jìn)程數(shù)的增加，裝配時(shí)間單調(diào)下降。

圖17 重疊單元和場單元切片(截面方程：z=0)

圖18 不同進(jìn)程對(duì)應(yīng)的墻上時(shí)間

3.2 WPFS模型

為了進(jìn)一步檢驗(yàn)OGSM重疊單元與貢獻(xiàn)單元插值關(guān)系的準(zhǔn)確性，本文通過接口模式將其嵌入至流場求解器PMB3D[25]中，對(duì)WPFS(Wing-Pylon-Finned Store)投放算例的非定常流場進(jìn)行了數(shù)值模擬。

WPFS模型[26]是美國發(fā)起的用于CFD驗(yàn)證的機(jī)翼外掛物分離投放簡化模型，其中單外掛物分離投放試驗(yàn)由AEDC于1990年完成，該模型作為標(biāo)準(zhǔn)算例通常用于對(duì)外掛物分離數(shù)值方法進(jìn)行考核。如圖19所示，研究模型的機(jī)翼為45°后掠的三角翼，采用NACA64A010翼型，無扭轉(zhuǎn)。外掛物的頭部為尖拱形，中間為圓柱，共有4片尾翼，呈X形分布。尾翼采用NACA0008翼型。外掛物質(zhì)心位置距頭部頂點(diǎn)1.417 3 m，參考長度0.508 m，參考面積0.202 69 m2。實(shí)驗(yàn)?zāi)M研究兩個(gè)狀態(tài)下的投放分離，狀態(tài)1的環(huán)境高度為7 925 m(26 000 ft)，來流馬赫數(shù)0.95；狀態(tài)2的環(huán)境高度為11 582 m(38 000 ft)，來流馬赫數(shù)1.2，機(jī)翼無攻角和側(cè)滑角。為了確保外掛物安全分離，在投放初始瞬間，加入了兩點(diǎn)彈射力的作用，彈射力1的作用點(diǎn)距頭部尖點(diǎn)1.24 m，大小10 679.4 N，彈射力2距頭部尖點(diǎn)1.75 m，大小42 717.5 N。在外掛物離開載機(jī)0.1 m(分離時(shí)間0.055 s)以后，彈射力消失。圖19中給出了外掛物的質(zhì)量以及3個(gè)方向上的轉(zhuǎn)動(dòng)慣量，由于外掛物關(guān)于軸和軸對(duì)稱，它的慣性積Ixy=Iyz=Izx=0 kg·m2。

圖19 外掛分離試驗(yàn)?zāi)Ｐ蛥?shù)

WPFS算例中，在機(jī)翼和彈體部件周圍生成多塊結(jié)構(gòu)網(wǎng)格，前者對(duì)應(yīng)的節(jié)點(diǎn)總數(shù)和體單元總數(shù)分別為5 124 792和4 798 080，后者對(duì)應(yīng)的節(jié)點(diǎn)總數(shù)和體單元總數(shù)分別為1 722 242和1 585 664。整個(gè)計(jì)算區(qū)域劃分到64個(gè)進(jìn)程中。在每個(gè)進(jìn)程上，合并重合的節(jié)點(diǎn)單元和面單元，多塊結(jié)構(gòu)網(wǎng)格即轉(zhuǎn)化為獨(dú)立的非結(jié)構(gòu)網(wǎng)格(原有各個(gè)網(wǎng)格塊可以是不連通的)。使用OGSM進(jìn)行重疊裝配后，依舊采用多塊結(jié)構(gòu)網(wǎng)格求解器進(jìn)行RANS方程的非定常數(shù)值模擬，其中對(duì)流項(xiàng)和黏性項(xiàng)分別使用Roe格式和SA湍流模式進(jìn)行離散，時(shí)間方向上采用雙時(shí)間步推進(jìn)，結(jié)合minmod限制器抑制非物理解的生成。

在雙時(shí)間步推進(jìn)過程中，流場求解器PMB3D的每個(gè)物理迭代步劃分為300個(gè)子迭代步，單個(gè)子迭代步消耗的時(shí)間為1.49 s；因?yàn)閺楏w相對(duì)于機(jī)翼處于不同位置時(shí)，依據(jù)壁面距離準(zhǔn)則得到的挖洞邊界形態(tài)各不相同，從而導(dǎo)致重疊裝配時(shí)間發(fā)生一定的變化，所以本文僅選取初始時(shí)刻對(duì)OGSM模塊的裝配時(shí)間進(jìn)行統(tǒng)計(jì)，對(duì)應(yīng)的裝配時(shí)間為9.29 s，占每個(gè)物理迭代步消耗時(shí)間的比率為2.1%，能夠滿足當(dāng)前算例的數(shù)值模擬需求；同時(shí)，壁面距離的計(jì)算時(shí)間為3.12 s，可疑查詢節(jié)點(diǎn)的并行確認(rèn)時(shí)間為4.63 s，它們與總裝配時(shí)間的比率分別為33.5%和49.9%；相比之下，生成背景網(wǎng)格二叉樹的時(shí)間為0.22 s，重疊單元插值系數(shù)的并行計(jì)算時(shí)間為0.69 s，它們與總裝配時(shí)間的比率僅為2.3%和7.4%。結(jié)合本文的算法流程可以判斷，隨著進(jìn)程數(shù)的增加，切割盒子序列的結(jié)構(gòu)未發(fā)生明顯變化，為適應(yīng)壁面節(jié)點(diǎn)規(guī)模不斷增加的情形，今后需對(duì)壁面距離的高效并行算法做更加深入的研究；同時(shí)，基于部件網(wǎng)格的幾何特征，挖掘新的初篩判據(jù)，在早期排除更多的待查詢節(jié)點(diǎn)，也有助于進(jìn)一步降低并行通訊壓力，從而提升OGSM的裝配效率。

圖20給出的是狀態(tài)1(Ma=0.95)外掛物表面壓力系數(shù)Cp和空間流線，可以看到初始時(shí)刻外掛物受到來流外洗的作用，產(chǎn)生向右的偏航力矩。在該算例中，機(jī)翼為前緣后掠角45°的三角翼，從翼根開始，氣流在受到機(jī)翼的阻滯作用后，產(chǎn)生沿橫向流動(dòng)的速度分量。低于外掛物來說，頭部靠機(jī)翼內(nèi)側(cè)的地方迎著橫向流動(dòng)，壓力系數(shù)較另一側(cè)高，因此產(chǎn)生頭部向機(jī)翼外側(cè)方向偏轉(zhuǎn)的偏航力矩。

圖20 外掛物表面壓力系數(shù)和空間流線分布

圖21顯示的是外掛物處于掛載狀態(tài)時(shí)展向截面馬赫數(shù)云圖，機(jī)翼上下表面大部分區(qū)域都是超聲速流動(dòng)，靠近機(jī)翼后緣的地方，上下表面均形成了激波，波后流動(dòng)速度為亞聲速。圖21同時(shí)給出了截面上的網(wǎng)格分布，從中能夠清晰地看到重疊邊界的位置，在跨過重疊邊界的地方，壓力系數(shù)云圖分布連續(xù)，反映了在這些地方流場信息得到了正確的傳遞。

圖21 外掛物處于掛載狀態(tài)時(shí)展向截面馬赫數(shù)云圖

圖22給出的是外掛物分離軌跡示意圖，其中模擬時(shí)長t為0.5 s。從圖中可以看到，外掛物在分離過程中，產(chǎn)生了頭部向機(jī)翼外側(cè)的偏航運(yùn)動(dòng)；在俯仰方向上，分離初始階段產(chǎn)生的是抬頭運(yùn)動(dòng)，隨后在尾翼的穩(wěn)定作用下產(chǎn)生低頭運(yùn)動(dòng)。圖22反映了外掛物分離過程的總體運(yùn)動(dòng)趨勢，定性地描述了運(yùn)動(dòng)過程，下面將給出定量描述。

圖22 外掛物分離軌跡

圖23和圖24分別給出的是外掛物質(zhì)心位置(dx、dy、dz)和線速度(U、V、W)變化曲線，可以看到在重力和彈射力作用下，外掛物向機(jī)翼下方移動(dòng)，其運(yùn)動(dòng)幅度在3個(gè)方向上是最大的；同時(shí)在阻力作用下，外掛物向機(jī)翼的后方移動(dòng)；在橫向方向上(y方向)，外掛物先是向機(jī)翼的內(nèi)側(cè)移動(dòng)，然后在t=0.35 s以后，外掛物開始向機(jī)翼外側(cè)移動(dòng)，盡管移動(dòng)幅度較小，計(jì)算仍比較準(zhǔn)確地模擬這種變化趨勢?？偟膩碚f，在所顯示的時(shí)間范圍內(nèi)(t=0.5 s)，計(jì)算得到的質(zhì)心線速度和位移和實(shí)驗(yàn)值都比較吻合，稍顯不足的是由于黏性計(jì)算的阻力大于試驗(yàn)值，所以模擬得到的后向移動(dòng)速度略大。

圖23 外掛物質(zhì)心位置變化曲線

圖24 外掛物線速度變化曲線

圖25和圖26分別給出的是外掛物的姿態(tài)角和角速度的時(shí)間歷程，其中P、Q、R分別為體軸系的3個(gè)角速度分量，Φ、Θ、Ψ分別為滾轉(zhuǎn)角、俯仰角和偏航角，用來表示從地軸系到體軸系的變換。對(duì)于俯仰運(yùn)動(dòng)來說，在分離初始階段，彈射力產(chǎn)生較大的抬頭力矩，彈射力維持的時(shí)間到t=0.055 s， 之后彈射力作用消失，作用在外掛物上的氣動(dòng)力和力矩開始起主導(dǎo)作用，這時(shí)候氣動(dòng)力產(chǎn)生的是低頭力矩，在t=0.2 s以后，外掛物開始下俯運(yùn)動(dòng)，本文計(jì)算的最大俯仰角度和下俯開始位置都和試驗(yàn)值吻合良好。在偏航方向，從分離開始外掛物偏向機(jī)翼的外側(cè)，在t=0.5 s偏航角約為16°，在t≤0.4 s計(jì)算的偏航角與試驗(yàn)值吻合良好，t>0.4 s后開始出現(xiàn)偏差，影響到后續(xù)計(jì)算的正確性。在滾轉(zhuǎn)方向，分離后外掛彈體向機(jī)翼外側(cè)滾轉(zhuǎn)，在t=0.5 s之前滾轉(zhuǎn)角的計(jì)算和試驗(yàn)值比較一致，由于x方向的轉(zhuǎn)動(dòng)慣量較小(比其他兩個(gè)方向小1個(gè)量級(jí))，氣動(dòng)力估計(jì)的誤差很容易在這里得到放大，從而對(duì)計(jì)算結(jié)果造成污染。

圖25 外掛物姿態(tài)角隨時(shí)間變化情況

圖26 外掛物角速度隨時(shí)間變化情況

結(jié)合數(shù)值模擬方法，本文認(rèn)為空間離散和時(shí)間離散的誤差來源于如下3個(gè)方面：① 重疊網(wǎng)格隱式裝配算法中，為簡明計(jì)，沿著挖洞邊界僅延拓了一層重疊單元，各部件網(wǎng)格間的物理信息交換量略顯不充分；② 雙時(shí)間推進(jìn)過程中，限制器造成的數(shù)值解污染在一定程度上得以累積，影響了后期數(shù)值模擬的精度；③ SA湍流模型經(jīng)驗(yàn)參數(shù)的選取欠佳，影響了RANS方程黏性項(xiàng)的離散精度。今后需繼續(xù)研究重疊邊界附近的高精度插值技術(shù)；同時(shí)，不斷提高流場解算器的離散精度，從而在整體上降低CFD仿真過程中產(chǎn)生的誤差。

4 結(jié) 論

本文基于網(wǎng)格節(jié)點(diǎn)的壁面距離參數(shù)，構(gòu)造了一種自動(dòng)化的并行重疊網(wǎng)格隱式裝配算法，同時(shí)使用集合論對(duì)該算法的并行化過程進(jìn)行了詳細(xì)地描述。網(wǎng)格組極小內(nèi)部-外部盒子、盒子切割方法、初篩+確認(rèn)兩步分析法、二叉樹排序-查找等技術(shù)的綜合運(yùn)用，降低了重疊裝配算法的并行通訊壓力，使之能夠適應(yīng)多體相對(duì)運(yùn)動(dòng)時(shí)復(fù)雜非定常流場的數(shù)值模擬。數(shù)值計(jì)算結(jié)果表明，該重疊裝配算法具備了較高的計(jì)算效率與魯棒性，能夠輔助實(shí)現(xiàn)針對(duì)外掛物分離的非定常流場準(zhǔn)確模擬。

下一步工作內(nèi)容包括兩方面：一是將自適應(yīng)笛卡爾網(wǎng)格融入到當(dāng)前的重疊裝配框架中，進(jìn)一步拓展其工程應(yīng)用范圍；二是繼續(xù)優(yōu)化初篩技術(shù)，提升大規(guī)模重疊網(wǎng)格的并行裝配效率。