高效非結構網格變形與流場插值方法

郭中州，何志強，趙文文，陳偉芳

浙江大學 航空航天學院，杭州 310027

含有動邊界的非定常流動問題在航空航天領域十分常見，如火箭的級間分離、飛行器的機動、外掛物投放等。這類問題通常伴隨著復雜的氣動干擾和強烈的非定常效應。對這類非定常流動問題的精確和高效模擬已經成為了當前計算流體力學(Computational Fluid Dynamics, CFD)和航空航天領域的研究熱點之一。在此基礎上發展的“數值虛擬飛行”技術，更是被稱為航空領域CFD應用的最后一個大挑戰[1]。

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

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

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

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

網格變形、網格重構以及流場插值方法是非結構動網格的核心技術。此方面的研究仍在不斷發展之中。本文基于KD(K-Dimension)樹數據結構，提出了高效的網格變形方法和流場插值方法以及相應的并行化方法。發展了靈活設置邊界運動的UDF功能，并通過幾個算例的模擬對方法進行了驗證。

1 高效網格變形方法

彈簧近似方法[2，25-26]和徑向基函數(Radial Basis Functions, RBF)方法[27-30]是目前應用最為廣泛的兩類非結構網格變形方法，彈簧近似方法的物理意義清晰，實現簡單，但是變形能力較弱，只能處理較小幅度的網格變形。相比之下，RBF方法具有變形能力強，不受網格拓撲限制等優點。但是原始的RBF方法需要求解大型線性方程組，方程組的階數與動邊界的節點數目相同，當動邊界上的網格較密時，求解方程組的計算量很大。導致計算效率下降。

此外，RBF方法需要指定支撐半徑，網格的變形都局限在支撐半徑之內，無法傳遞到更遠的網格區域。對于動邊界與固定邊界相距較近的問題，難以選取合適的支撐半徑，因而需要頻繁進行網格重構，影響計算精度和計算效率。

1.1 KNN-RBF網格變形方法

為了高效地處理非結構網格的大變形問題，提出了K近鄰-徑向基函數(K Nearest Neighbor-Radial Basis Functions, KNN-RBF)網格變形方法。在計算網格變形時，首先分別查找網格點到動邊界的K個最近點及到其余邊界的K個最近點，然后通過式(1)計算網格點的位移。

(1)

式(1)中K的值可根據計算效率和網格來選取，代表了動邊界上對網格點的位移有影響的區域范圍，通常取K≤10即可。選取較大的K值時，表明網格點位移由動邊界上多個最近點的位移決定，可以更好地將動邊界的位移傳遞到網格點處，但計算量也會相應增大。對動邊界上網格較密或者外形較為簡單的算例，取K=1即可。

徑向基函數的形式多種多樣，本文采用了在網格變形中廣泛應用的Wendland C2[31]型基函數。其形式為

(2)

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

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

可以看出，RBF方法的主要瓶頸在于求解線性方程組。相比于RBF方法，KNN-RBF方法解決了支撐半徑的選取問題，同時大大提高了網格變形過程的計算效率。

圖1 KD樹的數據結構Fig.1 Data structure of KD tree

1.2 網格變形的并行化方法

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

在計算之前，首先使用多級K路圖算法[34]對網格進行分區，將原始網格分成多個數目均衡的網格塊，實現并行計算時各進程的負載平衡。

網格變形的每一步都需要使網格的分區邊界上對應點的坐標在變形前后保持一致，避免變形后出現網格重疊或撕裂。這是網格變形方法并行化的關鍵步驟。基于紅黑樹[35]數據結構，設計了交界面節點的一致性方法。

讀入網格后，為每個網格點設置一個唯一的全局編號，記為Uuid。進行網格分區時，交界面上的節點被分到了不同的網格塊中，如圖2所示。再為每個節點設置一個局部的編號，這樣，相鄰網格塊的分區邊界上的對應節點具有相同的全局編號。在讀入各個分區的網格時，建立分區對接點局部編號和全局編號之間的映射關系并采用紅黑樹存儲?？梢钥焖俑咝У夭檎揖植烤幪柡腿志幪柕膶P系。

圖2 交界面上對接節點的對應關系Fig.2 Corresponding relationships of patch nodes on zonal boundaries

在進行網格變形時，首先分別對各個網格塊采用網格變形方法計算各自網格點的位移，然后各進程進行通信，各相鄰的網格塊互相傳遞交界面上的對接點的全局編號和變形之后的新坐標值。各分區根據接收的節點全局編號找到在本分區的對應點，將接收的節點坐標與本分區對應點的坐標進行加權平均，作為交界面上點的新坐標。這樣就保證了分區交界面上網格點坐標的一致性。

基于上述的交界面節點一致性方法，實現了彈簧近似方法[2]、貪婪RBF方法[28]和KNN-RBF等多種網格變形方法的并行化。對于彈簧近似方法，并行策略為先計算各分區的網格點位移，然后設置交界面節點位移。對于RBF方法和KNN-RBF方法，需要各進程先收集所有邊界點，然后再計算各分區網格點的位移，最后設置交界面節點的位移。

2 邊界運動的UDF功能實現

對含有動邊界的非定常流動問題，流場中可能有多個運動物體，如圖3所示，每個物體又由多個部件組成。例如多體分離、子母彈拋撒等問題，且不同物體可能以不同的規律運動，如變形、強迫運動以及六自由度運動等等。為了能夠靈活設置邊界的運動，在自主發展的非結構混合網格程序GUO-CFD (General Unstructured Oriented-Computational Fluid Dynamics)中實現了UDF功能，可以方便地設置動邊界的各種運動。

使用UDF設置邊界的運動時，首先在邊界條件文件中將運動規律相同的邊界劃分為一個運動組，然后為每個運動組編寫相應的UDF，在UDF中設置運動參數的變化規律，計算時CFD程序通過動態調用UDF來改變運動組的運動參數。

圖3 運動組的定義Fig.3 Definition of motion groups

為了使CFD程序與UDF之間的數據交互通用化，采用了JSON (JavaScript Object Notation)格式的文件來描述每個運動組的運動參數。JSON是一種標準化的數據交換格式，可以靈活地表達不同邊界的運動。以某個運動組的六自由度運動為例，對應的JSON格式的運動參數如圖4所示。

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

解釋型UDF采用Python語言腳本實現，可以定義動邊界的質心坐標、質量、慣性矩、速度、角速度、外力及姿態角等參數的變化規律，通過返回JSON格式的運動參數與CFD程序進行交互，實現對邊界運動的控制。無需編譯即可在不同系統上運行。但是相比于編譯型UDF功能有限，需要安裝Python運行環境，且運行速度比編譯型UDF慢。

圖4 JSON表示的動邊界運動參數Fig.4 JSON represented motion parameters of moving boundary

3 網格重構

根據重構過程涉及的網格范圍，可以將重構策略分為局部重構和全局重構。

局部重構只在質量較差的網格單元附近區域重新生成網格，網格生成階段的計算量較小。但網格生成后需要與舊網格進行合并，消除網格的復邊。網格數據結構和流場變量都需要進行動態調整。過程比較復雜，只能改善局部的網格質量，通用性較差。且并行情況下分區交界面上的網格質量也可能變差，難以進行并行化處理。

相比而言，全局重構在網格生成階段的計算量較大，但可以有效地改善整體網格質量，對網格數據的處理也比較方便，可以方便地調用網格生成接口。并且能夠保持重構后并行計算的負載平衡，因此后續采用了全局網格重構的策略。在程序中實現了調用其他網格生成程序的接口，可以直接調用自編網格生成程序或者商業軟件。

3.1 網格質量準則

對于四面體網格，采用單元的半徑比來衡量網格質量，其表達式為

(3)

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

對于混合網格，q的表達式為

(4)

式中：Li為單元各頂點到單元中心的距離；E為Li的平均值；Cn為單元的節點數目。

在計算時需要設置q的閾值，當網格的最差質量低于此值時則認為需要進行網格重構。

3.2 網格重構流程

按照前述的網格質量準則，每步網格變形之后，對各進程的網格進行檢測，當發現網格質量較差時，進行網格重構。網格重構的步驟如下：

步驟1備份當前的計算結果，將各個進程的邊界網格合并成一個整體，作為網格重構的輸入。

步驟2阻塞其他進程，主進程調用面網格光順程序，對合并后的邊界網格進行調整，保證邊界網格的質量。然后調用網格生成子程序或者商業軟件，生成新網格并進行光順。保證重構后的網格質量滿足計算要求。

步驟3對新網格進行分區，各進程讀入對應的新網格塊并替換舊網格塊，網格的數據結構如圖5所示，需要相應更新網格的單元、節點和面元等信息，同時保持原有邊界的邊界條件、運動組和UDF等不變。

整個網格變形和重構過程如圖6所示，重構過程中不需要人工干預。

圖5 非結構網格的數據結構Fig.5 Data structure of unstructured mesh

圖6 網格的并行變形和重構過程Fig.6 Parallel mesh deformation and regeneration procedure

4 并行流場插值方法

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

新舊網格的分區通?；ゲ恢睾?。如圖7所示，若采用順序查找，需要在所有網格塊中查找貢獻單元。當網格量很大時，整個插值過程將耗時巨大。為了定位貢獻單元，要進行大量的幾何關系判斷，如果新網格點由于誤差等原因查找失敗，還要進行額外處理，健壯性難以保證。

圖7 舊網格和新網格的網格分區Fig.7 Mesh partitions of old mesh and regenerated mesh

4.1 兩級KD樹插值方法

為提高插值過程的計算效率和健壯性，基于KD樹數據結構，提出了兩級KD樹方法，實現并行條件下對貢獻單元的快速查找。方法分為以下步驟：

步驟1查找貢獻單元所在的網格塊。分別將舊網格各網格塊的邊界點構建成KD樹，作為第1級KD樹，在其中查找新網格單元的最近點，如果由新網格單元點指向最近點的矢量與舊網格邊界在最近點的外法向的點積為正，則認為貢獻單元位于此網格塊內。在靠近分區界面處，如果網格較為扭曲，可能有不止一個網格塊滿足條件，則認為貢獻單元位于多個網格塊內。

步驟2在對應的網格塊中定位貢獻單元。將舊網格各個網格塊的單元中心點分別構建成KD樹，作為第2級KD樹，在定位了新網格的點位于舊網格的哪個分區之后，再利用這個分區的單元點的KD樹，得到新網格點在舊網格上的鄰近點，即為插值的貢獻單元。

步驟3使用貢獻單元的物理量進行插值。采用線性插值公式，根據鄰近點的空間位置關系將物理量插值到新網格點。

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

4.2 時間復雜度分析

假設并行計算時有Np個網格塊，每個網格塊的單元數目為Nc，每個網格塊的邊界網格點數目為Nb，且新舊網格的單元數目為同一量級。如果采用線性查找，需要依次遍歷各個分區的網格以定位貢獻單元。定位一個貢獻單元的平均時間復雜度為O(NpNc)。

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

5 算例驗證

5.1 標量函數插值

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

圖8 給定函數在非結構網格上的插值誤差Fig.8 Interpolation errors of given function on unstructured meshes

5.2 NACA0012翼型流場插值

采用前述的并行流場插值方法，對NACA0012翼型的定常流場插值進行了驗證。計算采用的兩套網格如圖9所示。

計算狀態對應的來流馬赫數為0.755，迎角為2.526°。首先采用六面體網格計算出流場結果，然后將計算結果插值到非結構網格上，得到的壓力等值線如圖10所示。圖10中的粗線為截面與網格塊邊界的交線。可以看出，插值前后流場的壓力等值線云圖一致性符合較好，驗證了插值方法的有效性。

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

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

5.3 外掛物投放模擬

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

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

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

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

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

在整個過程中，彈簧近似方法進行了15次重構，貪婪RBF方法進行了11次重構，KNN-RBF方法進行了4次網格重構。

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

為了考察3種方法的計算效率，統計了計算過程中各進程的動網格模塊占用的CPU時間的平均值。以KNN-RBF方法的平均CPU時間為基準，用無量綱CPU時間來衡量3種方法的計算效率，結果如圖14所示。對本算例而言，只考慮網格變形子程序占用的時間時，彈簧近似方法的效率最高。KNN-RBF方法和貪婪RBF方法都需要收集各網格塊的邊界點，故并行計算效率比彈簧方法低。而計入網格重構和流場插值占用的時間后，由于KNN-RBF方法變形能力強，所需重構次數少，相比其他兩種方法具有明顯優勢。同時可看出網格重構和流場插值會對計算效率產生重要影響。

圖13 重構前后的網格分布對比Fig.13 Comparison of mesh distribution before and after regeneration

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

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

圖15 外掛物運動參數的CFD結果與試驗值的對比Fig.15 Comparison of CFD and test results of motion parameters of store

6 結 論

1) 發展了高效的KNN-RBF網格變形方法，解決了RBF方法計算效率低和支撐半徑的選取問題。通過采用高效的數據結構，發展了網格分區節點一致性方法，實現了網格變形的并行化。

2) 在自主開發的非結構CFD程序中實現了編譯型和解釋型UDF功能，能夠靈活設置動邊界的運動規律，增強了程序求解含動邊界非定常流動問題時的擴展能力。

3) 采用全局網格重構策略，解決了大變形情況下網格質量下降的問題。針對重構后流場的插值問題，提出了基于兩級KD樹的并行流場插值方法，具有高效、健壯和低耗散的特點，實現了并行條件下新舊網格流場的快速插值。

4) 通過模擬解析函數、翼型定常流動和外掛物投放等算例，驗證了網格變形和流場插值方法的適用性。此外，兩級KD樹插值方法只需要網格點的坐標信息，對網格拓撲沒有限制，不僅適用于非結構動網格，也可應用于重疊網格和滑移網格等網格類型的流場插值問題中。