高效非結(jié)構(gòu)網(wǎng)格變形與流場(chǎng)插值方法

郭中州，何志強(qiáng)，趙文文，陳偉芳

浙江大學(xué) 航空航天學(xué)院，杭州 310027

含有動(dòng)邊界的非定常流動(dòng)問(wèn)題在航空航天領(lǐng)域十分常見(jiàn)，如火箭的級(jí)間分離、飛行器的機(jī)動(dòng)、外掛物投放等。這類問(wèn)題通常伴隨著復(fù)雜的氣動(dòng)干擾和強(qiáng)烈的非定常效應(yīng)。對(duì)這類非定常流動(dòng)問(wèn)題的精確和高效模擬已經(jīng)成為了當(dāng)前計(jì)算流體力學(xué)(Computational Fluid Dynamics, CFD)和航空航天領(lǐng)域的研究熱點(diǎn)之一。在此基礎(chǔ)上發(fā)展的“數(shù)值虛擬飛行”技術(shù)，更是被稱為航空領(lǐng)域CFD應(yīng)用的最后一個(gè)大挑戰(zhàn)[1]。

動(dòng)網(wǎng)格方法由于計(jì)算效率高，適合處理邊界變形問(wèn)題等特點(diǎn)，成為解決此類問(wèn)題的一種重要方法，得到了廣泛的關(guān)注。國(guó)外從20世紀(jì)八九十年代就開(kāi)展了利用動(dòng)網(wǎng)格技術(shù)求解動(dòng)邊界問(wèn)題的研究，Batina[2]首先使用彈簧動(dòng)網(wǎng)格模擬了NACA0012翼型的俯仰運(yùn)動(dòng)。Cavallo等[3]采用非結(jié)構(gòu)動(dòng)網(wǎng)格以及網(wǎng)格自適應(yīng)技術(shù)模擬了外掛物投放和火箭級(jí)間分離。但在模擬過(guò)程中進(jìn)行網(wǎng)格調(diào)整時(shí)需要進(jìn)行人工干預(yù)。Panagiotopoulos和Kyparissis[4]采用非結(jié)構(gòu)動(dòng)網(wǎng)格和歐拉方程，模擬了外掛物投放，并與試驗(yàn)結(jié)果進(jìn)行了對(duì)比。Luo和Baysal[5]模擬了高超聲速飛行器從助推器上的分離過(guò)程。Bartels[6]模擬了帶有擾流板的機(jī)翼的氣動(dòng)彈性問(wèn)題。Baum等[7]研究了F16C/D戰(zhàn)斗機(jī)拋擲油箱的動(dòng)態(tài)過(guò)程，并得到了與風(fēng)洞試驗(yàn)和飛行數(shù)據(jù)吻合良好的結(jié)果。

在國(guó)內(nèi)，劉君等[8]采用非結(jié)構(gòu)動(dòng)網(wǎng)格方法和歐拉方程，對(duì)飛行器的配平迎角計(jì)算、多體分離和結(jié)構(gòu)動(dòng)力學(xué)等問(wèn)題進(jìn)行了模擬。白曉征[9]結(jié)合網(wǎng)格變形和網(wǎng)格重構(gòu)并耦合多體動(dòng)力學(xué)，模擬了內(nèi)嵌式助推器的分離過(guò)程。陳堅(jiān)強(qiáng)等[10]采用結(jié)構(gòu)動(dòng)網(wǎng)格技術(shù)，研究了方形彈舵面操縱下導(dǎo)彈的俯仰運(yùn)動(dòng)響應(yīng)過(guò)程。張來(lái)平等[11]介紹了氣動(dòng)/運(yùn)動(dòng)/控制耦合一體化計(jì)算方法，探討了數(shù)值虛擬飛行中的一些挑戰(zhàn)性問(wèn)題。常興華[12]使用非結(jié)構(gòu)動(dòng)網(wǎng)格模擬了魚(yú)體的巡游和昆蟲(chóng)翼的拍動(dòng)。伍貽兆等[13]分別采用非結(jié)構(gòu)動(dòng)網(wǎng)格和非結(jié)構(gòu)重疊網(wǎng)格模擬了三維外掛物投放、直升機(jī)機(jī)身+旋翼前飛以及機(jī)翼氣動(dòng)彈性變形等問(wèn)題。王剛等[14]提出了基于貪心法逐級(jí)選擇徑向基函數(shù)空間子集的方法，提高了徑向基函數(shù)方法計(jì)算網(wǎng)格變形時(shí)的效率。陳龍等[15]發(fā)展了Delaunay圖映射彈簧混合動(dòng)網(wǎng)格方法，并應(yīng)用到機(jī)翼氣動(dòng)彈性計(jì)算中。

目前，國(guó)外成熟的商業(yè)CFD軟件中大多集成了重疊網(wǎng)格方法或者動(dòng)網(wǎng)格方法。如Fastran[16]采用了重疊網(wǎng)格方法。Fluent[17]中采用了非結(jié)構(gòu)動(dòng)網(wǎng)格方法。一些國(guó)外的開(kāi)源軟件，如OpenFOAM[18]，也集成了多種非結(jié)構(gòu)動(dòng)網(wǎng)格方法。由于邊界運(yùn)動(dòng)的形態(tài)多種多樣，商業(yè)軟件中大多具備一定的擴(kuò)展功能，如Fluent的用戶自定義函數(shù)(User Defined Function, UDF)功能模塊，以方便用戶處理各種不同類型的動(dòng)網(wǎng)格問(wèn)題。

當(dāng)邊界變形導(dǎo)致網(wǎng)格質(zhì)量過(guò)差時(shí)，為了改善網(wǎng)格質(zhì)量，需要進(jìn)行網(wǎng)格重構(gòu)，造成計(jì)算效率下降。重構(gòu)后的流場(chǎng)插值過(guò)程還會(huì)引入額外的耗散，甚至導(dǎo)致流場(chǎng)產(chǎn)生波動(dòng)，發(fā)展高效、健壯的插值方法至關(guān)重要。Smith和Quon[19]發(fā)展了基于點(diǎn)云和徑向基函數(shù)的重疊網(wǎng)格流場(chǎng)插值方法。Hner[20]總結(jié)了非結(jié)構(gòu)網(wǎng)格上的插值方法的不足。Olsson和Petersson[21]研究了一維重疊網(wǎng)格插值的穩(wěn)定性問(wèn)題。王永健[22]提出了一類適用于任意網(wǎng)格的守恒重映方法。董海波等[23]提出了基于格點(diǎn)格式有限體積法的流場(chǎng)數(shù)據(jù)傳遞方法。利用非結(jié)構(gòu)動(dòng)網(wǎng)格技術(shù)，將舊網(wǎng)格單元移動(dòng)到新網(wǎng)格單元，同時(shí)求解流場(chǎng)控制方程，以此來(lái)實(shí)現(xiàn)兩套網(wǎng)格間的信息傳遞。徐喜華等[24]提出了一種三維非結(jié)構(gòu)多面體二階保界全局重映算法，在舊網(wǎng)格上選取模板利用最小二乘構(gòu)造插值多項(xiàng)式，采用凸包算法計(jì)算多面體相交，最后使用局部保界修正技術(shù)修補(bǔ)重映后的越界量。

網(wǎng)格變形、網(wǎng)格重構(gòu)以及流場(chǎng)插值方法是非結(jié)構(gòu)動(dòng)網(wǎng)格的核心技術(shù)。此方面的研究仍在不斷發(fā)展之中。本文基于KD(K-Dimension)樹(shù)數(shù)據(jù)結(jié)構(gòu)，提出了高效的網(wǎng)格變形方法和流場(chǎng)插值方法以及相應(yīng)的并行化方法。發(fā)展了靈活設(shè)置邊界運(yùn)動(dòng)的UDF功能，并通過(guò)幾個(gè)算例的模擬對(duì)方法進(jìn)行了驗(yàn)證。

1 高效網(wǎng)格變形方法

彈簧近似方法[2，25-26]和徑向基函數(shù)(Radial Basis Functions, RBF)方法[27-30]是目前應(yīng)用最為廣泛的兩類非結(jié)構(gòu)網(wǎng)格變形方法，彈簧近似方法的物理意義清晰，實(shí)現(xiàn)簡(jiǎn)單，但是變形能力較弱，只能處理較小幅度的網(wǎng)格變形。相比之下，RBF方法具有變形能力強(qiáng)，不受網(wǎng)格拓?fù)湎拗频葍?yōu)點(diǎn)。但是原始的RBF方法需要求解大型線性方程組，方程組的階數(shù)與動(dòng)邊界的節(jié)點(diǎn)數(shù)目相同，當(dāng)動(dòng)邊界上的網(wǎng)格較密時(shí)，求解方程組的計(jì)算量很大。導(dǎo)致計(jì)算效率下降。

此外，RBF方法需要指定支撐半徑，網(wǎng)格的變形都局限在支撐半徑之內(nèi)，無(wú)法傳遞到更遠(yuǎn)的網(wǎng)格區(qū)域。對(duì)于動(dòng)邊界與固定邊界相距較近的問(wèn)題，難以選取合適的支撐半徑，因而需要頻繁進(jìn)行網(wǎng)格重構(gòu)，影響計(jì)算精度和計(jì)算效率。

1.1 KNN-RBF網(wǎng)格變形方法

為了高效地處理非結(jié)構(gòu)網(wǎng)格的大變形問(wèn)題，提出了K近鄰-徑向基函數(shù)(K Nearest Neighbor-Radial Basis Functions, KNN-RBF)網(wǎng)格變形方法。在計(jì)算網(wǎng)格變形時(shí)，首先分別查找網(wǎng)格點(diǎn)到動(dòng)邊界的K個(gè)最近點(diǎn)及到其余邊界的K個(gè)最近點(diǎn)，然后通過(guò)式(1)計(jì)算網(wǎng)格點(diǎn)的位移。

(1)

式(1)中K的值可根據(jù)計(jì)算效率和網(wǎng)格來(lái)選取，代表了動(dòng)邊界上對(duì)網(wǎng)格點(diǎn)的位移有影響的區(qū)域范圍，通常取K≤10即可。選取較大的K值時(shí)，表明網(wǎng)格點(diǎn)位移由動(dòng)邊界上多個(gè)最近點(diǎn)的位移決定，可以更好地將動(dòng)邊界的位移傳遞到網(wǎng)格點(diǎn)處，但計(jì)算量也會(huì)相應(yīng)增大。對(duì)動(dòng)邊界上網(wǎng)格較密或者外形較為簡(jiǎn)單的算例，取K=1即可。

徑向基函數(shù)的形式多種多樣，本文采用了在網(wǎng)格變形中廣泛應(yīng)用的Wendland C2[31]型基函數(shù)。其形式為

(2)

可以看出，KNN-RBF與其他RBF方法有較大區(qū)別：① KNN-RBF方法不需要設(shè)置支撐半徑。② RBF方法需要通過(guò)求解線性方程組來(lái)得到邊界點(diǎn)的權(quán)值，再根據(jù)權(quán)值計(jì)算網(wǎng)格點(diǎn)的位移。而KNN-RBF方法只需要進(jìn)行K近鄰查找。K近鄰的查找采用KD樹(shù)[32]實(shí)現(xiàn)。KD樹(shù)是一種平衡二叉樹(shù)，如圖1所示，它將多維數(shù)據(jù)按照各個(gè)維度進(jìn)行分割，然后形成一棵平衡二叉樹(shù)。利用二叉樹(shù)實(shí)現(xiàn)鄰近數(shù)據(jù)的快速查找。

設(shè)網(wǎng)格總的節(jié)點(diǎn)數(shù)目為M，動(dòng)邊界的節(jié)點(diǎn)數(shù)目為N，固定邊界的節(jié)點(diǎn)數(shù)目也為O(N)量級(jí)。使用KNN-RBF方法計(jì)算網(wǎng)格變形時(shí)，每個(gè)網(wǎng)格點(diǎn)需要查找2K個(gè)鄰近點(diǎn)。將動(dòng)邊界節(jié)點(diǎn)構(gòu)建成KD樹(shù)，利用KD樹(shù)查找一個(gè)鄰近點(diǎn)的平均時(shí)間復(fù)雜度為O(lbN)。則查找所有網(wǎng)格節(jié)點(diǎn)鄰近點(diǎn)的平均時(shí)間復(fù)雜度為O(2KM(lbN))，網(wǎng)格變形過(guò)程的時(shí)間復(fù)雜度為O(KM)。而RBF方法需要求解3個(gè)線性方程組，采用直接求解法，總時(shí)間復(fù)雜度為O(3N3)，網(wǎng)格變形過(guò)程的時(shí)間復(fù)雜度為O(MN)。對(duì)于網(wǎng)格總的節(jié)點(diǎn)數(shù)目和邊界節(jié)點(diǎn)數(shù)目的關(guān)系，可認(rèn)為O(M)=O(N1.5)。則KNN-RBF方法總時(shí)間復(fù)雜度為O(2KN1.5(lbN))，RBF方法總的時(shí)間復(fù)雜度為O(3N3)。

可以看出，RBF方法的主要瓶頸在于求解線性方程組。相比于RBF方法，KNN-RBF方法解決了支撐半徑的選取問(wèn)題，同時(shí)大大提高了網(wǎng)格變形過(guò)程的計(jì)算效率。

圖1 KD樹(shù)的數(shù)據(jù)結(jié)構(gòu)Fig.1 Data structure of KD tree

1.2 網(wǎng)格變形的并行化方法

為了適應(yīng)復(fù)雜外形的計(jì)算精度和效率的要求，使用消息傳遞接口[33](Message Passing Interface, MPI)技術(shù)在程序中實(shí)現(xiàn)了KNN-RBF方法的并行化。

在計(jì)算之前，首先使用多級(jí)K路圖算法[34]對(duì)網(wǎng)格進(jìn)行分區(qū)，將原始網(wǎng)格分成多個(gè)數(shù)目均衡的網(wǎng)格塊，實(shí)現(xiàn)并行計(jì)算時(shí)各進(jìn)程的負(fù)載平衡。

網(wǎng)格變形的每一步都需要使網(wǎng)格的分區(qū)邊界上對(duì)應(yīng)點(diǎn)的坐標(biāo)在變形前后保持一致，避免變形后出現(xiàn)網(wǎng)格重疊或撕裂。這是網(wǎng)格變形方法并行化的關(guān)鍵步驟。基于紅黑樹(shù)[35]數(shù)據(jù)結(jié)構(gòu)，設(shè)計(jì)了交界面節(jié)點(diǎn)的一致性方法。

讀入網(wǎng)格后，為每個(gè)網(wǎng)格點(diǎn)設(shè)置一個(gè)唯一的全局編號(hào)，記為Uuid。進(jìn)行網(wǎng)格分區(qū)時(shí)，交界面上的節(jié)點(diǎn)被分到了不同的網(wǎng)格塊中，如圖2所示。再為每個(gè)節(jié)點(diǎn)設(shè)置一個(gè)局部的編號(hào)，這樣，相鄰網(wǎng)格塊的分區(qū)邊界上的對(duì)應(yīng)節(jié)點(diǎn)具有相同的全局編號(hào)。在讀入各個(gè)分區(qū)的網(wǎng)格時(shí)，建立分區(qū)對(duì)接點(diǎn)局部編號(hào)和全局編號(hào)之間的映射關(guān)系并采用紅黑樹(shù)存儲(chǔ)?？梢钥焖俑咝У夭檎揖植烤幪?hào)和全局編號(hào)的對(duì)應(yīng)關(guān)系。

圖2 交界面上對(duì)接節(jié)點(diǎn)的對(duì)應(yīng)關(guān)系Fig.2 Corresponding relationships of patch nodes on zonal boundaries

在進(jìn)行網(wǎng)格變形時(shí)，首先分別對(duì)各個(gè)網(wǎng)格塊采用網(wǎng)格變形方法計(jì)算各自網(wǎng)格點(diǎn)的位移，然后各進(jìn)程進(jìn)行通信，各相鄰的網(wǎng)格塊互相傳遞交界面上的對(duì)接點(diǎn)的全局編號(hào)和變形之后的新坐標(biāo)值。各分區(qū)根據(jù)接收的節(jié)點(diǎn)全局編號(hào)找到在本分區(qū)的對(duì)應(yīng)點(diǎn)，將接收的節(jié)點(diǎn)坐標(biāo)與本分區(qū)對(duì)應(yīng)點(diǎn)的坐標(biāo)進(jìn)行加權(quán)平均，作為交界面上點(diǎn)的新坐標(biāo)。這樣就保證了分區(qū)交界面上網(wǎng)格點(diǎn)坐標(biāo)的一致性。

基于上述的交界面節(jié)點(diǎn)一致性方法，實(shí)現(xiàn)了彈簧近似方法[2]、貪婪RBF方法[28]和KNN-RBF等多種網(wǎng)格變形方法的并行化。對(duì)于彈簧近似方法，并行策略為先計(jì)算各分區(qū)的網(wǎng)格點(diǎn)位移，然后設(shè)置交界面節(jié)點(diǎn)位移。對(duì)于RBF方法和KNN-RBF方法，需要各進(jìn)程先收集所有邊界點(diǎn)，然后再計(jì)算各分區(qū)網(wǎng)格點(diǎn)的位移，最后設(shè)置交界面節(jié)點(diǎn)的位移。

2 邊界運(yùn)動(dòng)的UDF功能實(shí)現(xiàn)

對(duì)含有動(dòng)邊界的非定常流動(dòng)問(wèn)題，流場(chǎng)中可能有多個(gè)運(yùn)動(dòng)物體，如圖3所示，每個(gè)物體又由多個(gè)部件組成。例如多體分離、子母彈拋撒等問(wèn)題，且不同物體可能以不同的規(guī)律運(yùn)動(dòng)，如變形、強(qiáng)迫運(yùn)動(dòng)以及六自由度運(yùn)動(dòng)等等。為了能夠靈活設(shè)置邊界的運(yùn)動(dòng)，在自主發(fā)展的非結(jié)構(gòu)混合網(wǎng)格程序GUO-CFD (General Unstructured Oriented-Computational Fluid Dynamics)中實(shí)現(xiàn)了UDF功能，可以方便地設(shè)置動(dòng)邊界的各種運(yùn)動(dòng)。

使用UDF設(shè)置邊界的運(yùn)動(dòng)時(shí)，首先在邊界條件文件中將運(yùn)動(dòng)規(guī)律相同的邊界劃分為一個(gè)運(yùn)動(dòng)組，然后為每個(gè)運(yùn)動(dòng)組編寫(xiě)相應(yīng)的UDF，在UDF中設(shè)置運(yùn)動(dòng)參數(shù)的變化規(guī)律，計(jì)算時(shí)CFD程序通過(guò)動(dòng)態(tài)調(diào)用UDF來(lái)改變運(yùn)動(dòng)組的運(yùn)動(dòng)參數(shù)。

圖3 運(yùn)動(dòng)組的定義Fig.3 Definition of motion groups

為了使CFD程序與UDF之間的數(shù)據(jù)交互通用化，采用了JSON (JavaScript Object Notation)格式的文件來(lái)描述每個(gè)運(yùn)動(dòng)組的運(yùn)動(dòng)參數(shù)。JSON是一種標(biāo)準(zhǔn)化的數(shù)據(jù)交換格式，可以靈活地表達(dá)不同邊界的運(yùn)動(dòng)。以某個(gè)運(yùn)動(dòng)組的六自由度運(yùn)動(dòng)為例，對(duì)應(yīng)的JSON格式的運(yùn)動(dòng)參數(shù)如圖4所示。

為了方便不同情況下的使用，UDF分為了編譯型UDF和解釋型UDF。編譯型UDF采用C++ 語(yǔ)言實(shí)現(xiàn)，在程序首次執(zhí)行時(shí)將UDF文件編譯為動(dòng)態(tài)鏈接庫(kù)，通過(guò)運(yùn)行時(shí)加載其中的函數(shù)實(shí)現(xiàn)自定義的邊界運(yùn)動(dòng)。編譯之后可以在Windows和Linux平臺(tái)上執(zhí)行，運(yùn)行速度快?？梢灾苯釉O(shè)置動(dòng)邊界的各個(gè)運(yùn)動(dòng)參數(shù)，也可以通過(guò)返回JSON格式的運(yùn)動(dòng)參數(shù)與CFD程序進(jìn)行交互。

解釋型UDF采用Python語(yǔ)言腳本實(shí)現(xiàn)，可以定義動(dòng)邊界的質(zhì)心坐標(biāo)、質(zhì)量、慣性矩、速度、角速度、外力及姿態(tài)角等參數(shù)的變化規(guī)律，通過(guò)返回JSON格式的運(yùn)動(dòng)參數(shù)與CFD程序進(jìn)行交互，實(shí)現(xiàn)對(duì)邊界運(yùn)動(dòng)的控制。無(wú)需編譯即可在不同系統(tǒng)上運(yùn)行。但是相比于編譯型UDF功能有限，需要安裝Python運(yùn)行環(huán)境，且運(yùn)行速度比編譯型UDF慢。

圖4 JSON表示的動(dòng)邊界運(yùn)動(dòng)參數(shù)Fig.4 JSON represented motion parameters of moving boundary

3 網(wǎng)格重構(gòu)

根據(jù)重構(gòu)過(guò)程涉及的網(wǎng)格范圍，可以將重構(gòu)策略分為局部重構(gòu)和全局重構(gòu)。

局部重構(gòu)只在質(zhì)量較差的網(wǎng)格單元附近區(qū)域重新生成網(wǎng)格，網(wǎng)格生成階段的計(jì)算量較小。但網(wǎng)格生成后需要與舊網(wǎng)格進(jìn)行合并，消除網(wǎng)格的復(fù)邊。網(wǎng)格數(shù)據(jù)結(jié)構(gòu)和流場(chǎng)變量都需要進(jìn)行動(dòng)態(tài)調(diào)整。過(guò)程比較復(fù)雜，只能改善局部的網(wǎng)格質(zhì)量，通用性較差。且并行情況下分區(qū)交界面上的網(wǎng)格質(zhì)量也可能變差，難以進(jìn)行并行化處理。

相比而言，全局重構(gòu)在網(wǎng)格生成階段的計(jì)算量較大，但可以有效地改善整體網(wǎng)格質(zhì)量，對(duì)網(wǎng)格數(shù)據(jù)的處理也比較方便，可以方便地調(diào)用網(wǎng)格生成接口。并且能夠保持重構(gòu)后并行計(jì)算的負(fù)載平衡，因此后續(xù)采用了全局網(wǎng)格重構(gòu)的策略。在程序中實(shí)現(xiàn)了調(diào)用其他網(wǎng)格生成程序的接口，可以直接調(diào)用自編網(wǎng)格生成程序或者商業(yè)軟件。

3.1 網(wǎng)格質(zhì)量準(zhǔn)則

對(duì)于四面體網(wǎng)格，采用單元的半徑比來(lái)衡量網(wǎng)格質(zhì)量，其表達(dá)式為

(3)

式中：Rs和rs分別為四面體外接球和內(nèi)切球的半徑。

對(duì)于混合網(wǎng)格，q的表達(dá)式為

(4)

式中：Li為單元各頂點(diǎn)到單元中心的距離；E為L(zhǎng)i的平均值；Cn為單元的節(jié)點(diǎn)數(shù)目。

在計(jì)算時(shí)需要設(shè)置q的閾值，當(dāng)網(wǎng)格的最差質(zhì)量低于此值時(shí)則認(rèn)為需要進(jìn)行網(wǎng)格重構(gòu)。

3.2 網(wǎng)格重構(gòu)流程

按照前述的網(wǎng)格質(zhì)量準(zhǔn)則，每步網(wǎng)格變形之后，對(duì)各進(jìn)程的網(wǎng)格進(jìn)行檢測(cè)，當(dāng)發(fā)現(xiàn)網(wǎng)格質(zhì)量較差時(shí)，進(jìn)行網(wǎng)格重構(gòu)。網(wǎng)格重構(gòu)的步驟如下：

步驟1備份當(dāng)前的計(jì)算結(jié)果，將各個(gè)進(jìn)程的邊界網(wǎng)格合并成一個(gè)整體，作為網(wǎng)格重構(gòu)的輸入。

步驟2阻塞其他進(jìn)程，主進(jìn)程調(diào)用面網(wǎng)格光順程序，對(duì)合并后的邊界網(wǎng)格進(jìn)行調(diào)整，保證邊界網(wǎng)格的質(zhì)量。然后調(diào)用網(wǎng)格生成子程序或者商業(yè)軟件，生成新網(wǎng)格并進(jìn)行光順。保證重構(gòu)后的網(wǎng)格質(zhì)量滿足計(jì)算要求。

步驟3對(duì)新網(wǎng)格進(jìn)行分區(qū)，各進(jìn)程讀入對(duì)應(yīng)的新網(wǎng)格塊并替換舊網(wǎng)格塊，網(wǎng)格的數(shù)據(jù)結(jié)構(gòu)如圖5所示，需要相應(yīng)更新網(wǎng)格的單元、節(jié)點(diǎn)和面元等信息，同時(shí)保持原有邊界的邊界條件、運(yùn)動(dòng)組和UDF等不變。

整個(gè)網(wǎng)格變形和重構(gòu)過(guò)程如圖6所示，重構(gòu)過(guò)程中不需要人工干預(yù)。

圖5 非結(jié)構(gòu)網(wǎng)格的數(shù)據(jù)結(jié)構(gòu)Fig.5 Data structure of unstructured mesh

圖6 網(wǎng)格的并行變形和重構(gòu)過(guò)程Fig.6 Parallel mesh deformation and regeneration procedure

4 并行流場(chǎng)插值方法

重構(gòu)過(guò)程完成后，還需要將舊網(wǎng)格上的物理量插值到新網(wǎng)格上。為了保持流場(chǎng)的守恒性，理論上應(yīng)當(dāng)采用守恒的插值方法，但是目前還難以構(gòu)造一種高效、健壯而且通用的守恒性插值方法。目前在重疊網(wǎng)格和動(dòng)網(wǎng)格插值中都廣泛地采用了非守恒插值，并取得了良好效果。最常用的是線性插值，進(jìn)行插值時(shí)，首先要在舊網(wǎng)格上尋找貢獻(xiàn)單元，然后使用貢獻(xiàn)單元的物理量插值得到新網(wǎng)格上的物理量。

新舊網(wǎng)格的分區(qū)通常互不重合。如圖7所示，若采用順序查找，需要在所有網(wǎng)格塊中查找貢獻(xiàn)單元。當(dāng)網(wǎng)格量很大時(shí)，整個(gè)插值過(guò)程將耗時(shí)巨大。為了定位貢獻(xiàn)單元，要進(jìn)行大量的幾何關(guān)系判斷，如果新網(wǎng)格點(diǎn)由于誤差等原因查找失敗，還要進(jìn)行額外處理，健壯性難以保證。

圖7 舊網(wǎng)格和新網(wǎng)格的網(wǎng)格分區(qū)Fig.7 Mesh partitions of old mesh and regenerated mesh

4.1 兩級(jí)KD樹(shù)插值方法

為提高插值過(guò)程的計(jì)算效率和健壯性，基于KD樹(shù)數(shù)據(jù)結(jié)構(gòu)，提出了兩級(jí)KD樹(shù)方法，實(shí)現(xiàn)并行條件下對(duì)貢獻(xiàn)單元的快速查找。方法分為以下步驟：

步驟1查找貢獻(xiàn)單元所在的網(wǎng)格塊。分別將舊網(wǎng)格各網(wǎng)格塊的邊界點(diǎn)構(gòu)建成KD樹(shù)，作為第1級(jí)KD樹(shù)，在其中查找新網(wǎng)格單元的最近點(diǎn)，如果由新網(wǎng)格單元點(diǎn)指向最近點(diǎn)的矢量與舊網(wǎng)格邊界在最近點(diǎn)的外法向的點(diǎn)積為正，則認(rèn)為貢獻(xiàn)單元位于此網(wǎng)格塊內(nèi)。在靠近分區(qū)界面處，如果網(wǎng)格較為扭曲，可能有不止一個(gè)網(wǎng)格塊滿足條件，則認(rèn)為貢獻(xiàn)單元位于多個(gè)網(wǎng)格塊內(nèi)。

步驟2在對(duì)應(yīng)的網(wǎng)格塊中定位貢獻(xiàn)單元。將舊網(wǎng)格各個(gè)網(wǎng)格塊的單元中心點(diǎn)分別構(gòu)建成KD樹(shù)，作為第2級(jí)KD樹(shù)，在定位了新網(wǎng)格的點(diǎn)位于舊網(wǎng)格的哪個(gè)分區(qū)之后，再利用這個(gè)分區(qū)的單元點(diǎn)的KD樹(shù)，得到新網(wǎng)格點(diǎn)在舊網(wǎng)格上的鄰近點(diǎn)，即為插值的貢獻(xiàn)單元。

步驟3使用貢獻(xiàn)單元的物理量進(jìn)行插值。采用線性插值公式，根據(jù)鄰近點(diǎn)的空間位置關(guān)系將物理量插值到新網(wǎng)格點(diǎn)。

由于最近點(diǎn)必然存在，因此兩級(jí)KD樹(shù)插值方法的健壯性可以得到保證，其本質(zhì)上也屬于線性插值，具有一階精度。由于線性插值是基于泰勒展開(kāi)，到插值點(diǎn)的距離越近，插值引起的耗散越小，而兩級(jí)KD樹(shù)方法使用的是最近點(diǎn)進(jìn)行插值，因此相比其他線性插值方法具有更小的數(shù)值耗散。

4.2 時(shí)間復(fù)雜度分析

假設(shè)并行計(jì)算時(shí)有Np個(gè)網(wǎng)格塊，每個(gè)網(wǎng)格塊的單元數(shù)目為Nc，每個(gè)網(wǎng)格塊的邊界網(wǎng)格點(diǎn)數(shù)目為Nb，且新舊網(wǎng)格的單元數(shù)目為同一量級(jí)。如果采用線性查找，需要依次遍歷各個(gè)分區(qū)的網(wǎng)格以定位貢獻(xiàn)單元。定位一個(gè)貢獻(xiàn)單元的平均時(shí)間復(fù)雜度為O(NpNc)。

對(duì)于兩級(jí)KD樹(shù)方法，每個(gè)網(wǎng)格塊的邊界點(diǎn)都構(gòu)成一棵平衡二叉樹(shù)，在每棵樹(shù)中查找的平均時(shí)間復(fù)雜度為O(lbNb)，則查找貢獻(xiàn)單元所屬舊網(wǎng)格塊的平均時(shí)間復(fù)雜度為O(Np(lbNb))，在舊網(wǎng)格塊中查找貢獻(xiàn)單元的時(shí)間復(fù)雜度為O(lbNc)。整個(gè)過(guò)程總的時(shí)間復(fù)雜度為O(Np(lbNb)+lbNc)。通常情況下，有Np?Nb?Nc，則定位一個(gè)貢獻(xiàn)單元總的時(shí)間復(fù)雜度為O(lbNc)。因此，與線性查找相比，兩級(jí)KD樹(shù)方法具有更高的效率。

5 算例驗(yàn)證

5.1 標(biāo)量函數(shù)插值

采用了兩套非結(jié)構(gòu)網(wǎng)格。網(wǎng)格A為六面體單元，網(wǎng)格B為混合網(wǎng)格。函數(shù)的表達(dá)式為z=x2+y2，先根據(jù)函數(shù)直接得到網(wǎng)格格心處的函數(shù)值。然后插值到另外一個(gè)網(wǎng)格上，分別計(jì)算了由網(wǎng)格A插值到網(wǎng)格B和由網(wǎng)格B插值到網(wǎng)格A兩個(gè)算例，得到的截面上的插值誤差分布如圖8所示。最大插值誤差均在5%之內(nèi)，表明插值方法具有較高精度。

圖8 給定函數(shù)在非結(jié)構(gòu)網(wǎng)格上的插值誤差Fig.8 Interpolation errors of given function on unstructured meshes

5.2 NACA0012翼型流場(chǎng)插值

采用前述的并行流場(chǎng)插值方法，對(duì)NACA0012翼型的定常流場(chǎng)插值進(jìn)行了驗(yàn)證。計(jì)算采用的兩套網(wǎng)格如圖9所示。

計(jì)算狀態(tài)對(duì)應(yīng)的來(lái)流馬赫數(shù)為0.755，迎角為2.526°。首先采用六面體網(wǎng)格計(jì)算出流場(chǎng)結(jié)果，然后將計(jì)算結(jié)果插值到非結(jié)構(gòu)網(wǎng)格上，得到的壓力等值線如圖10所示。圖10中的粗線為截面與網(wǎng)格塊邊界的交線。可以看出，插值前后流場(chǎng)的壓力等值線云圖一致性符合較好，驗(yàn)證了插值方法的有效性。

圖9 插值采用的六面體網(wǎng)格和四面體網(wǎng)格Fig.9 Hexahedron and tetrahedron meshes used for interpolation

圖10 原始結(jié)果和插值得到的壓力云圖對(duì)比Fig.10 Comparison of pressure contours of original and interpolated results

5.3 外掛物投放模擬

算例外形為具有詳細(xì)試驗(yàn)數(shù)據(jù)的Wing/Pylon/Finned store外形[36]，該算例被廣泛用于驗(yàn)證動(dòng)網(wǎng)格和重疊網(wǎng)格程序的正確性。采用的網(wǎng)格為四面體網(wǎng)格，單元數(shù)目約156萬(wàn)，在掛架和外掛物之間的縫隙等細(xì)節(jié)位置對(duì)網(wǎng)格進(jìn)行了加密，表面網(wǎng)格的分布如圖11所示。

采用有限體積法求解ALE(Arbitrary Lagrangian-Euler)形式的控制方程，其中，梯度采用加權(quán)最小二乘法計(jì)算，對(duì)流通量格式采用考慮網(wǎng)格運(yùn)動(dòng)速度的Van Leer格式，限制器采用Venkatakrishnan限制函數(shù)。為提高計(jì)算效率，采用32個(gè) 計(jì)算核心并行計(jì)算。首先模擬了高度為10 km，馬赫數(shù)為1.2的定常流場(chǎng)作為初場(chǎng)。然后使用雙時(shí)間步長(zhǎng)方法模擬了外掛物投放后0.8 s內(nèi)的非定常流場(chǎng)。不同時(shí)刻外掛物位置如圖12所示。

圖11 Wing/Pylon/Finned store外形Fig.11 Configuration of Wing/Pylon/Finned store

圖12 不同時(shí)刻外掛物的位置Fig.12 Position of external store at different interval times

分別采用彈簧近似方法、貪婪RBF方法和KNN-RBF方法進(jìn)行了計(jì)算。其中，彈簧近似方法的迭代次數(shù)為50次，貪婪RBF方法的插值誤差取10-5m，支撐半徑為1 m，KNN-RBF方法取K=1，3種方法的最差網(wǎng)格質(zhì)量閾值均取qmin=0.001。

在整個(gè)過(guò)程中，彈簧近似方法進(jìn)行了15次重構(gòu)，貪婪RBF方法進(jìn)行了11次重構(gòu)，KNN-RBF方法進(jìn)行了4次網(wǎng)格重構(gòu)。

圖13為KNN-RBF方法重構(gòu)前后的網(wǎng)格對(duì)比，可以看出，隨著外掛物的下落，掛架和外掛物之間的網(wǎng)格被不斷拉長(zhǎng)，質(zhì)量明顯變差，但緊貼物面處的網(wǎng)格質(zhì)量保持較好。重構(gòu)之后，中間區(qū)域的網(wǎng)格分布得到改善。

為了考察3種方法的計(jì)算效率，統(tǒng)計(jì)了計(jì)算過(guò)程中各進(jìn)程的動(dòng)網(wǎng)格模塊占用的CPU時(shí)間的平均值。以KNN-RBF方法的平均CPU時(shí)間為基準(zhǔn)，用無(wú)量綱CPU時(shí)間來(lái)衡量3種方法的計(jì)算效率，結(jié)果如圖14所示。對(duì)本算例而言，只考慮網(wǎng)格變形子程序占用的時(shí)間時(shí)，彈簧近似方法的效率最高。KNN-RBF方法和貪婪RBF方法都需要收集各網(wǎng)格塊的邊界點(diǎn)，故并行計(jì)算效率比彈簧方法低。而計(jì)入網(wǎng)格重構(gòu)和流場(chǎng)插值占用的時(shí)間后，由于KNN-RBF方法變形能力強(qiáng)，所需重構(gòu)次數(shù)少，相比其他兩種方法具有明顯優(yōu)勢(shì)。同時(shí)可看出網(wǎng)格重構(gòu)和流場(chǎng)插值會(huì)對(duì)計(jì)算效率產(chǎn)生重要影響。

圖13 重構(gòu)前后的網(wǎng)格分布對(duì)比Fig.13 Comparison of mesh distribution before and after regeneration

圖14 3種方法的CPU時(shí)間比較Fig.14 Comparison of CPU time of three methods

圖15展示了外掛物質(zhì)心、歐拉角的計(jì)算值與試驗(yàn)值的對(duì)比。偏航角在整個(gè)過(guò)程中與試驗(yàn)值都吻合較好，俯仰角在0.6 s之前符合較好，之后偏差逐漸增大。而滾轉(zhuǎn)角在0.4 s之后偏差逐漸增大，由于外掛物滾轉(zhuǎn)慣性矩的量級(jí)較小，對(duì)氣動(dòng)力矩的誤差十分敏感，再加上計(jì)算時(shí)未考慮黏性影響，故其計(jì)算值和試驗(yàn)值差異較大。質(zhì)心位移在整個(gè)過(guò)程中符合較好，但后期也出現(xiàn)偏差增大的趨勢(shì)。總體而言，隨著計(jì)算時(shí)間的增加，上述6個(gè)運(yùn)動(dòng)參數(shù)的偏差都有不斷增大的趨勢(shì)。

圖15 外掛物運(yùn)動(dòng)參數(shù)的CFD結(jié)果與試驗(yàn)值的對(duì)比Fig.15 Comparison of CFD and test results of motion parameters of store

6 結(jié) 論

1) 發(fā)展了高效的KNN-RBF網(wǎng)格變形方法，解決了RBF方法計(jì)算效率低和支撐半徑的選取問(wèn)題。通過(guò)采用高效的數(shù)據(jù)結(jié)構(gòu)，發(fā)展了網(wǎng)格分區(qū)節(jié)點(diǎn)一致性方法，實(shí)現(xiàn)了網(wǎng)格變形的并行化。

2) 在自主開(kāi)發(fā)的非結(jié)構(gòu)CFD程序中實(shí)現(xiàn)了編譯型和解釋型UDF功能，能夠靈活設(shè)置動(dòng)邊界的運(yùn)動(dòng)規(guī)律，增強(qiáng)了程序求解含動(dòng)邊界非定常流動(dòng)問(wèn)題時(shí)的擴(kuò)展能力。

3) 采用全局網(wǎng)格重構(gòu)策略，解決了大變形情況下網(wǎng)格質(zhì)量下降的問(wèn)題。針對(duì)重構(gòu)后流場(chǎng)的插值問(wèn)題，提出了基于兩級(jí)KD樹(shù)的并行流場(chǎng)插值方法，具有高效、健壯和低耗散的特點(diǎn)，實(shí)現(xiàn)了并行條件下新舊網(wǎng)格流場(chǎng)的快速插值。

4) 通過(guò)模擬解析函數(shù)、翼型定常流動(dòng)和外掛物投放等算例，驗(yàn)證了網(wǎng)格變形和流場(chǎng)插值方法的適用性。此外，兩級(jí)KD樹(shù)插值方法只需要網(wǎng)格點(diǎn)的坐標(biāo)信息，對(duì)網(wǎng)格拓?fù)錄](méi)有限制，不僅適用于非結(jié)構(gòu)動(dòng)網(wǎng)格，也可應(yīng)用于重疊網(wǎng)格和滑移網(wǎng)格等網(wǎng)格類型的流場(chǎng)插值問(wèn)題中。