非結構嵌套網格中的一種改進型徑向基函數插值方法

靳晨暉，王剛，陳鑫，周豪

(西北工業大學 航空學院，西安 710072)

0 引 言

嵌套網格技術在航空航天領域有著廣泛的應用，針對一些復雜的幾何構型，使用嵌套網格技術可以分區生成網格，減少網格生成的難度。對于飛機副油箱分離[1]，座艙蓋分離/座椅彈射[2]，多級火箭助推器分離[3]等多體分離問題,不同物面間的相對位置會發生變化，同時伴有強烈的非定常效應。常規研究方法如風洞和飛行實驗常帶有一定的局限性，使用嵌套網格技術進行非定常計算可以對這些過程進行精確數值模擬。

在嵌套網格算法中，一般通過插值實現嵌套區域的流場信息傳遞，尤其是在解決帶有動邊界的問題，嵌套區域的網格在每個非定常步內都需要更新，插值模板也需要重新選取。因此，使用正確的插值方法是保證計算結果準確的關鍵。按照插值方法的不同可分為守恒型插值與非守恒型插值。對于守恒型插值而言，網格尺度對其插值精度的影響較小，但插值過程十分復雜；而非守恒型插值的實現則較為簡單，但需要控制嵌套區域網格的匹配程度，否則會損失計算精度。為了減少非守恒型插值過程中的精度損失，一般需要從兩方面入手：一方面要改善嵌套插值區域網格匹配程度，如在網格生成時，需要對背景網格中有相對運動的區域進行局部加密，同時使用高效的嵌套裝配算法[4]劃分嵌套區域，以此來保證物面運動的整個過程中插值單元與貢獻單元的匹配程度。另一方面需要研究不同插值方法的插值效果，C.Mastin等[5]研究了雙線性、三線性插值方法在嵌套區域的插值效果；田書玲等[6]研究了一種基于Lagrange插值的非結構嵌套網格插值方法；周乃春等[7]使用了逆向距離權方法對嵌套區域進行插值；黃宇等[8]使用了具有二階精度的插值方法，有效降低了由于在嵌套區域進行插值而引入的誤差。僅僅增加嵌套區域的網格量有時也難以達到較好的精度，網格量增加會導致仿真周期成倍增加，占用過多計算資源；同時，線性插值方法在處理復雜問題時的精度相對較低，選擇合適的插值方法對于處理不同的問題至關重要。

徑向基函數插值方法因其存儲十分方便和數據結構相對簡單的特點，在近二十年內迅猛發展。研究者們在不同的領域對徑向基函數的應用進行大量研究，如R.Franke[9]發現對散亂的數據進行擬合時，使用徑向基函數進行插值的精度要高于其他插值方法；林言中等[10]將徑向基函數插值方法應用于網格變形技術中的界面數據傳遞，對AGARD 445.6機翼顫振問題進行非定常計算；劉溢浪等[11]提出一種增量形式的RBF插值方法，并將其應用于有限體積方法的流場重構中；Wang G等[12-13]在RBF網格變形技術的基礎上，提出將拉普拉斯光順與網格變形相結合，有效改善了變形后的網格質量。但是，徑向基函數插值半徑的選取帶有一定的經驗性，傳統徑向基函數前述插值半徑選取不合適會嚴重影響插值的精度與效率。

為了解決嵌套區域的插值精度問題，本文對不同插值函數進行研究，同時采用一種改進的非結構嵌套網格徑向基函數(RBF)插值方法。在嵌套網格流場信息傳遞的過程中，除了使用線性插值(LINE)和以距離倒數作為權重插值(WA)兩種傳統插值方法，還對逆二次徑向基函數(IQ)和Wendland’s C2徑向基函數(C2)的插值效果進行研究。針對嵌套邊界區域出現的數值奇性現象，通過控制徑向基函數的作用半徑來調整插值矩陣條件數的方式進行改善。

1 非結構嵌套網格求解技術

1.1 控制方程

為便于對固壁邊界的六自由度運動進行描述，采用ALE(Arbitrary Lagrangeian-Eulerian)方法對流動控制方程進行描述。ALE方法允許計算網格的進行剛體運動和變形，通過在流動控制方程中加入網格運動速度，將流體力學中的Lagrange方法和Euler方法進行統一描述。使用ALE方法描述的三維非定常雷諾平均Navier-Stokes方程(URANS)的積分形式為

(1)

式中：Ω為控制體；?Ω為控制體單元邊界；Q為守恒變量，Q=[ρρuρvρwρE]T，其中ρ為流體密度,u、v、w分別為旋轉彈在體軸系下三個軸向的運動速度E為總內能;n為控制體單元邊界上沿外法線方向的分量；F(Q)和G(Q)分別為N-S方程中的無粘通量項和有粘通量項。

對于式(1)，使用格心有限體積法(FVM)進行空間離散后，得到半離散形式為

(2)

(3)

式中：R為計算殘差；S為網格單元的面積；n為時間步；Vgrid為網格進行剛體運動的速度。為了保證足夠的時間精度，使用二階向后的歐拉隱式方法[14]求解式(2)，如式(4)所示

(4)

式中：Δt為時間步長；上標n為真實時間迭代的步數。直接求解式(4)較為復雜，一般采用偽時間迭代的方法對其進行求解，為此在方程的左端加入一個守恒變量對虛擬時間τ的導數，從而方程的解可以等價為求解一個關于虛擬時間τ的一階常微分方程組在虛擬時間τ趨近無窮時的漸進解：

(5)

(6)

由于式(5)的精度與偽時間迭代方法無關，因此，偽時間迭代過程可以使用一些加速收斂方法來提高流場的計算效率。具體方法可以參照文獻[15-16]。采用改進的LU-SGS[17]算法對式(5)進行向后的歐拉隱式時間迭代，同時使用基于OpenMP并行算法來提高計算的效率。

1.2 嵌套區域劃分與貢獻單元查找

嵌套網格裝配技術的兩個核心問題是如何確定嵌套邊界和如何查找貢獻單元。確定嵌套邊界也叫 “挖洞”，當計算網格量較大時，會顯著增加預處理的時間，因此高效并且合理地劃分出嵌套區域將是嵌套網格裝配技術的關鍵。本文在嵌套區域劃分時采取的策略是根據壁面距離來確定嵌套邊界區域，同時計算子網格到自身物面和其他物面的最短距離，通過比較到不同物面的最短距離來確定嵌套邊界，相比于傳統的枚舉法，有效減少了前處理時間。節點屬性判定如圖1所示，i點是屬于A網格的網格節點，在求解距離時，距離A網格更近，則i點的屬性標記為激活的狀態；而j點距離B網格更近，則j點的屬性標記為非激活的狀態。

圖1 節點屬性判定示意圖Fig.1 Schematic diagram of node attribute determination

當確定了嵌套邊界后，需要在嵌套區域查找貢獻單元。對于復雜模型來講，網格量較大，查找貢獻單元往往要花費較多時間，因此如何快速的查找到貢獻單元是提高嵌套效率的關鍵。本文根據文獻[4]采用了一種基于非結構嵌套網格的“蛇形查找方法”(如圖2所示)，在空間中以近乎最佳路徑的方式查找到了貢獻單元，有效縮短了查找貢獻單元的時間與距離。假設要尋找i點的貢獻單元，該方法的具體實現過程如下：

(1) 在空間中賦予i點一個初始的出發單元，這里假設是1號網格單元；

(2) 連接當前的出發單元的網格格心和i點，判斷其穿過當前出發單元的哪一條邊(三維為面)；

(3) 若穿過某一邊(面)，則將該邊(面)的另一側單元當做出發單元，循環第二步；

(4) 若沒有穿過當前的出發單元所有的邊(面)，則當前的出發單元就是i點的貢獻單元，結束查找。

圖2 蛇形查找示意圖Fig.2 Schematic diagram of snake search

通過以上步驟，嵌套區域的網格裝配完畢。在進行流場計算時，網格實際被劃分為三個區域，分別是計算區域、關閉區域和嵌套插值區域。在計算區域，進行正常的CFD求解，當接觸到邊界時，面外側位于嵌套區域的網格信息由貢獻單元插值得到；在關閉區域，只計算網格單元體積和網格單元面積等預處理操作；在嵌套插值區域，在每個偽時間迭代結束后進行插值，嵌套插值區域的網格不進行計算，而是被當作嵌套邊界來處理，這類邊界具有插值后獲得的流場信息。

1.3 插值方法簡介

流場信息的插值傳遞方法較多，本文主要考察了線性插值、距離倒數權重插值兩種傳統插值方法，以及逆二次徑向基函數、Wendland’s C2徑向基函數插值方法。由于線性插值結構較為簡單，此處不再贅述，下面重點對其余三種插值方法進行簡要介紹。

以距離倒數作為權重來進行插值的原理是：對待插值點附近點進行加權平均，其中權重的大小與需要插值的點與相鄰點之間的距離密切相關，與距離呈反比，因此也叫做距離反比插值方法。其具體形式為

(7)

本文所使用的徑向基函數的基本形式為

(8)

式中：下標i為支撐點；r為該點的位置矢量；N為網格節點的總數；φ(‖rp-ri‖)為所采用的基函數的形式；ωi為與第i個支撐點相關的權重系數。徑向基函數大致可分為以下三類：全局型(Global)徑向基函數、局部型(Local)徑向基函數、緊支型(Compact)徑向基函數，這里簡單介紹一下本文所使用的全局型/局部型(IQ)與緊支型(C2)兩種徑向基函數。

(1) 逆二次徑向基函數方法

(9)

在對嵌套區域的P點進行流場信息傳遞時，根據之前查找得到的貢獻單元，將貢獻單元以及它的相鄰單元一起作為徑向基函數的插值模版構建徑向基函數的系數矩陣。插值模版示意圖如圖3所示，i網格是P的貢獻單元，則使用i、j、k、l四個網格構建徑向基函數系數矩陣。

圖3 插值模版示意圖Fig.3 Interpolation template schematic diagram

為了選取比較光滑的模版值，在選取模版時，每個物理量減去他們的平均值，最后再加上被減去的值。將各個點的流場信息做差后的ΔZ帶入到RBF插值函數之中，得到的線性方程組如下：

ΦW=ΔZ

(10)

(11)

W=[ωi,ωj,ωk,ωl]T

(12)

式中：‖ri-rj‖為歐式距離，即每個基網格格心之間的距離，求解方程后得到每個基網格的權重系數ωi，將需要插值的節點P坐標rP帶入到徑向基函數中就可以重構出該點的值。

(2) Wendland’s C2型徑向基函數

基本形式如下所示：

(13)

插值半徑的選取對于徑向基函數插值方法的效率有著很重大的影響。較大的插值半徑會擴大物面插值的影響域，與此同時，也會使得基矩陣的維度增加并引起計算量的增加。選擇較小的插值半徑可以減小矩陣維度，降低計算量。

1.4 改進的徑向基函數插值策略

徑向基函數插值精度高，靈活性很強，使用不同類型的插值函數的插值效果也有所差異。對于同一個徑向基函數，選取不同的參數d也會得到不同的插值效果。本文所使用的C2徑向基函數，作用半徑d的初始值為插值模版中格心和P點(嵌套區域需要進行插值的點)的最大距離的50倍。矩陣的條件數不僅是判斷矩陣病態與否的一種度量，也可以表示矩陣計算對誤差的敏感度強弱，條件數越大，數值穩定性越差。對于C2徑向基函數的系數矩陣Φ，取二范數后的條件數可以表示為：

(14)

為了改善在一些網格上的條件數過大的情況，本文采用了一種基于控制條件數的徑向基函數插值算法，通過人為的改變C2徑向基函數的作用半徑d來調整插值矩陣的條件數，從而在條件數過大的插值模版上進行有效的控制，有效改善了由于條件數過大使得穩定性較差以及出現數值奇性的現象。具體實施的步驟為：當流場求解開始時，人為輸入一個期望達到的條件數，當使用徑向基函數插值進行流場信息傳遞時，每次迭代對貢獻單元插值矩陣的插值半徑自動乘以0.9，計算所有貢獻單元插值矩陣的條件數，不滿足要求則返回上一步繼續縮小插值矩陣半徑，直到所有貢獻單元插值矩陣的條件數都滿足要求。

2 計算結果與分析

2.1 插值精度分析

為了分析本文所使用的插值函數進行解析的精度，使用解析函數作為插值精度驗證算例。該函數的數學形式為：

(15)

嵌套網格組裝過程示意圖如圖4所示，圓心在(0,0)處,半徑為1。

(a) 嵌套組裝前網格

(b) 嵌套完成后網格圖4 嵌套網格組裝過程示意圖Fig.4 Diagram of overset grid assembly process

使用線性插值，以距離倒數作為權重插值，逆二次徑向基函數插值、Wendland’s C2徑向基函數插值四種不同的方法進行插值，四種插值方法的誤差統計如表1所示。

表1 四種插值方式在嵌套區域的誤差統計Table 1 Error statistics of four interpolation modes in nest area

從表1可以看出：使用C2徑向基函數進行插值具有較小的誤差,得到的結果和解析解最為接近。在整個插值的區域內，C2徑向基函數插值的平均絕對誤差僅為0.007 254 9%，而IQ徑向基函數插值為0.036 369%，距離反比函數插值(WA)為0.197 52%，線性插值(LINE)為2.361 2%，說明使用C2徑向基函數進行插值具有更高的計算精度，本文將優先選擇使用C2徑向基函數插值進行流場信息傳遞。

2.2 MD30P/30N多段翼流場數值模擬

當直接使用C2徑向基函數用于嵌套網格插值過程中時，在某些點可能存在數值原因引起的振蕩或插值誤差過大的點。以MD30P/30N多段翼為例，計算馬赫數Ma=0.2，雷諾數Re=9.0×106，迎角α=8.010 9°，湍流模型選擇Spalart-Allmaras(SA)。圍繞前緣縫翼、主翼、后緣襟翼分別生成計算網格，其中前緣縫翼網格量約為6萬，后緣襟翼網格量約為11萬，主翼網格量約為22萬。三套網格按照之前距離判定的準則自動進行網格嵌套組裝。三套網格均按照預期完成組裝示意圖如圖5所示。

圖5 MD30P/30N嵌套網格組裝示意圖Fig.5 Assembly schematic of MD30P/30N overset grid

未采用控制條件數與采用控制條件數后的C2徑向基函數插值得到的流場圖如圖6所示。從流場計算結果來看，直接使用C2函數進行插值得到了圖6(a)所示的結果，可以看出：在嵌套區域存在有數值奇性。根據之前提出的條件數控制策略，對于不同的插值模板，通過調整他們的插值半徑來控制條件數。經過條件數控制為10，即cond=10時，得到了圖6(b)所示的計算結果，整個流場的求解更加光滑，消除了之前出現的數值奇性現象。

(a) 未采用控制條件數

(b) 采用控制條件數圖6 MD30P/30N流場計算壓力分布云圖Fig.6 Pressure distribution cloud diagram of MD30P/30N flow field

MD30P/30N表面Cp分布與實驗對比示意圖如圖7所示，可以看出，當直接使用C2徑向基函數插值進行流場信息傳遞時，表面Cp分布在主翼上表面與后緣襟翼下表面與實驗數據差別較大，當采取條件數控制策略的C2徑向基函數插值策略將矩陣條件數控制在10后，所得到的翼型表面Cp分布與實驗數據吻合度較高，說明該插值方法能夠很好地解決直接使用C2徑向基函數插值帶來的數值奇性問題，且計算精度較高。MD30P/30N殘差收斂歷程對比示意圖如圖8所示，當條件數控制在10時，殘差下降的速率較快，說明采取條件數控制策略的C2徑向基函數插值策略能有效改善殘差的收斂效果。

圖7 MD30P/30N表面Cp分布示意圖Fig.7 Schematic diagram of MD30P/30N surface Cp distribution

圖8 MD30P/30N殘差收斂歷程對比示意圖Fig.8 Comparison diagram of MD30P/30N residual convergence process

2.3 AEDC外掛物投放標模驗證

AEDC標準外掛物分離模型的計算參數[18-19]如圖9所示。計算馬赫數Ma=0.95，雷諾數Re=7.874×106，迎角α=0°，湍流模型選擇Spalart-Allmaras(SA)。圍繞機翼和外掛物分別生成計算網格，網格總數約450萬，兩套網格按照之前的距離判定的準則自動進行網格嵌套組裝。

圖9 標準外掛物分離模型的計算參數Fig.9 The calculation parameters of standard separation model

為了考察使用控制條件數策略的C2徑向基函數的插值效果，將插值矩陣的條件數控制在100，即cond=100時，如圖10所示。從流場計算結果來看，可以看出：未采用控制條件數時，在嵌套區域的邊界上，出現了數值奇性，這些異常導致在尾跡區出現了和周圍顯著不同的區域；當采用控制條件數后，尾跡區流場過渡的更加光滑和連續，這說明采用條件數控制策略的插值方法可以有效消除嵌套邊界區域出現的數值奇性的問題。

(a) 無條件數控制的z-x截面

(b) 采用條件數控制的z-x截面

(c) 無條件數控制的x-y截面

(d) 采用條件數控制的x-y截面圖10 外掛物截面流場對比示意圖Fig.10 Flow field comparison diagram of the cross section

外掛物表面Cp分布與實驗值對比示意圖如圖11所示，可以看出：直接使用C2徑向基函數和未使用控制條件數策略的C2徑向基函數進行插值得到的外掛物周向Cp分布均與實驗值較為吻合，當加入控制插值矩陣條件數后，外掛物在彈尾處的Cp分布與實驗值吻合度更高一些。

(a) 外掛物兩側表面Cp分布

(b) 外掛物上下表面Cp分布圖11 外掛物表面Cp分布與實驗對比示意圖Fig.11 Schematic diagram of Cp distribution on the surface of the plug-in and comparison with the experiment

AEDC外掛物殘差收斂歷程對比曲線如圖12所示，可以看出：未采用控制條件數時，殘差收斂的歷程波動較為劇烈；采用控制插值矩陣條件數后，殘差收斂效果較好，呈現不斷下降的趨勢，說明采用條件數控制策略的C2徑向基函數插值策略能提高殘差收斂的效果。

圖12 殘差收斂歷程對比Fig.12 Comparison of residual convergence process

非定常數值模擬得到的外掛物質心處的位移Xg、Yg、Zg與實驗值的對比結果如圖13所示，可以看出：采用控制條件數策略的C2徑向基函數插方法數值模擬得到的水平位移Xg整體與實驗吻合度更高一些。兩種不同的插值方式得到外掛物沿豎直方向的位移與實驗結果吻合較好。

圖13 質心位移與實驗值的對比圖Fig.13 Comparison diagram of center of mass displacement and experimental value

外掛物分離過程中的姿態角變化曲線如圖14所示。

圖14 姿態角與實驗值的對比圖Fig.14 Comparison diagram of attitude angle and experimental value

從圖14可以看出：使用控制條件數策略的C2徑向基函數進行插值得到的俯仰角θ變化曲線在分離中后期與實驗吻合更好一些，而偏航角ψ與滾轉角φ使用兩種插值方法數值模擬得到的結果較為接近，均與實驗結果吻合度較高。同時可以看出彈射力矩與氣動力的共同作用，使得外掛物的俯仰角θ會呈現先增加后減小的趨勢；外掛物受到三角翼下洗與側洗的影響，偏航角也會逐漸增加。

3 結 論

(1) 使用解析函數作為插值模板，對線性插值、距離倒數作為權重插值，逆二次徑向基函數插值、Wendland’s C2徑向基函數插值四種不同的嵌套網格插值方法的精度進行驗證。計算結果表明C2徑向基函數的插值精度較高，同時收斂較快。

(2) 改進后的C2徑向基函數插值方法能夠有效消除嵌套插值區域數值奇性，同時計算收斂的速率更快，魯棒性較好。