基于支持向量回歸的動網(wǎng)格技術研究

高 翔，廖海翔，徐傳福

(國防科技大學 計算機學院，量子信息研究所兼高性能計算國家重點實驗室，長沙 410073)

0 引 言

動網(wǎng)格生成技術是指網(wǎng)格做相應調整以適應模擬過程中計算區(qū)域的改變，是計算流體力學模擬運動界面流動問題的主要瓶頸之一。動網(wǎng)格方法在氣動外形優(yōu)化設計、流固耦合模擬和氣動彈性計算等領域有著廣泛的應用。為提升計算精度和結果可靠性，亟需一種高效健壯的動網(wǎng)格方法重新生成或移動網(wǎng)格以滿足計算域的變化。第一種方法是直接在各時間步內重構整體或局部的計算網(wǎng)格，但復雜區(qū)域的網(wǎng)格生成十分耗時，且新舊網(wǎng)格間的插值傳遞會引入額外的誤差[1]。第二種方法只變化網(wǎng)格節(jié)點坐標，但不改變連接關系，從而避免流場的插值傳遞，并能處理各種類型的網(wǎng)格單元，因此得到了廣泛使用，但這種方法很難直接用于多體分離等大尺度運動問題。

在上述第二種方法中，超限插值法（transfinite interpolation，TFI）按網(wǎng)格方向順序根據(jù)邊界點與內部點的相對位置更新網(wǎng)格，該方法快速有效但只適用于結構網(wǎng)格[2]。對于沒有規(guī)則方向的非結構網(wǎng)格運動，一般采用以下兩類方法：第一類將網(wǎng)格邊類比為彈簧[3]或彈性實體[4]，根據(jù)網(wǎng)格拓撲構建力平衡系統(tǒng)，因此需反復求解規(guī)模為網(wǎng)格節(jié)點數(shù)目的線性方程組，計算成本很大。第二類只基于與邊界點的距離計算各網(wǎng)格內部點的新坐標，由于計算量少且易于并行，近年來得到了廣泛應用。基于Delaunay 背景網(wǎng)格的插值方法先建立計算網(wǎng)格節(jié)點與背景單元的映射關系，邊界運動后網(wǎng)格節(jié)點即可根據(jù)相應背景單元的邊界點快速計算得到新坐標[5-7]。該方法效率較高且能推廣至三維空間，但目前還只能用于凸區(qū)域的網(wǎng)格運動。逆距加權插值法（inverse distance weighting，IDW）通過顯式函數(shù)移動網(wǎng)格內點，不需要求解隱式方程，但對于不同的運動形式和邊界類型需采用不同的指數(shù)函數(shù)，健壯性較差[8-11]。徑向基函數(shù)插值法（radial basis function，RBF）根據(jù)邊界點坐標和運動位移擬合一個徑向基函數(shù)來表示整個計算域內網(wǎng)格點位移的分布。由此需求解一個規(guī)模為邊界點數(shù)目的線性方程組得到函數(shù)的待定系數(shù)[12]。為進一步提高其計算效率，可只選擇部分邊界點作為控制點來計算分布函數(shù)。主要思路是基于貪心思想，給定初始的控制點集合，由此得到擬合函數(shù)后計算剩余邊界點的位移，然后再與其已知位移比較并將差異最大的若干個邊界點加入該控制點集合，繼而繼續(xù)迭代直至滿足設置的誤差條件[13-14]。這種方法能夠大幅度降低動網(wǎng)格的計算和存儲成本，因此被廣泛應用于各類問題。

IDW 方法和RBF 方法最初都是加權處理所有邊界點對各內部點的影響，因此效率相對較低。而Delaunay 背景網(wǎng)格插值法只基于背景單元的幾個邊界點影響相應內部點，變形能力又相對較弱。近年來，也出現(xiàn)了一些結合背景網(wǎng)格方法與IDW 方法或RBF 方法的研究[11,15]，主要思路是先將IDW 或RBF方法用于稀疏的背景網(wǎng)格運動變形，然后也維持計算網(wǎng)格節(jié)點與背景網(wǎng)格單元的相對位置不變，從而插值計算出網(wǎng)格節(jié)點的新坐標。該方法綜合了兩類算法的優(yōu)點，但需要另外再生成一套計算域的稀疏背景網(wǎng)格?；诓逯祷驍M合函數(shù)，嘗試選擇部分邊界點與所有內部點、甚至不同邊界點與不同內部點建立相應映射關系，從而優(yōu)化和加速動網(wǎng)格方法，這是一個值得繼續(xù)研究的方向。

本文基于上述思想，進一步放寬貪心選點的徑向基函數(shù)插值法的約束條件，提出并完善了一種基于支持向量回歸（support vector regression，SVR）的動網(wǎng)格方法，該方法將貪心選點的插值問題轉化為自適應選點的優(yōu)化問題。在前期的工作中，已初步基于高斯核函數(shù)與貪心RBF 方法進行了對比，在保持網(wǎng)格運動健壯性的同時SVR 方法具有更高的計算效率。本文進一步通過系統(tǒng)性測試分析，以高斯核函數(shù)的SVR 方法為基準，針對SVR 動網(wǎng)格方法給出了幾個性質良好的核函數(shù)及其參數(shù)設置策略，并在若干典型二維/三維動網(wǎng)格案例上進行了成功應用，進一步優(yōu)化和完善了SVR 動網(wǎng)格算法。

1 方法描述

1.1 貪心選點RBF 方法約束條件的放寬

RBF 動網(wǎng)格方法在各坐標方向分別采用一個RBF 函數(shù)，以描述整個計算域內任意點在該方向將要移動的位移量，一般的表達形式如下：

求解上述線性方程組確定RBF 函數(shù)的系數(shù)后，從而得到了網(wǎng)格的運動函數(shù)。

為進一步提高RBF 方法在求解大規(guī)模問題時的計算效率，貪心選點的思想是逐次新增若干邊界點作為控制點，由此重新計算得到RBF 插值函數(shù)，最終將其他非控制點的邊界點位移量誤差控制在一定范圍內，從而減少待求方程組的規(guī)模。其數(shù)學內涵可以描述如下：

若進一步將公式（6）的約束條件放寬，使得對所有邊界點都允許計算得到的位移與已知位移存在一定誤差，即：

并在此基礎上求解某個最值問題，則將貪心選點的RBF 插值方法轉換成了一個約束優(yōu)化問題。而支持向量回歸模型正是解決這類問題的高效方法，且其回歸函數(shù)的表達形式也是由若干徑向基函數(shù)線性組合而成。因此，完全可以嘗試采用支持向量回歸作為網(wǎng)格運動的位移分布函數(shù)。

1.2 支持向量回歸

支持向量回歸的機器學習方法發(fā)展于統(tǒng)計學習的理論，提出近三十年來已廣泛應用于諸多領域的模式識別、分類以及回歸分析。如人臉識別、文本分類以及復雜工程分析的近似等[16-18]。SVR 的核心原理是基于一個非線性核函數(shù)將樣本空間映射到高維特征空間，并將原問題轉化為在此高維空間求解一個線性擬合問題[19]?；谥С窒蛄康母拍?，SVR 所得到的擬合函數(shù)即要求所有樣本點都落在給定的誤差范圍內，并使得函數(shù)在特征空間的斜率達到最小。

將SVR 應用于求解網(wǎng)格運動時，本文提出的方法把所有邊界點坐標及其某一坐標方向的已知位移量作為訓練樣本集，令其為：

其中假設線性情況下的擬合函數(shù)為：

如圖1 所示的樣本集擬合函數(shù)，所有訓練樣本點都落在函數(shù)向兩側平移 ε而形成的管道內。從中可知該函數(shù)實際上僅取決于部分樣本點，這些實質決定回歸函數(shù)的訓練樣本被稱為支持向量，也就相當于貪心選點RBF 方法中選擇出來的控制點。

圖1 SVR 在樣本空間的示意圖Fig. 1 Schematic of SVR in the sample space

公式（9）的二次規(guī)劃問題一般采用拉格朗日乘數(shù)法（Lagrange Multiplier）轉化為對偶問題再更有效地求解。對偶后的函數(shù)可表示為如下形式：

此外，采用非線性核函數(shù)k(x,xi)代 替內積 〈x,xi〉則可求解更高維特征空間的線性回歸。這等價于將樣本數(shù)據(jù)高維擴展后再做內積，其結果等同于將樣本代入核函數(shù)得到的值[20]。通過這種方式從而極大提高了SVR 的擬合能力，擴展后的回歸函數(shù)為：

關于支持向量回歸對偶轉化和非線性擴展更詳細的推導，可參閱文獻[21]。

1.3 SVR 針對網(wǎng)格運動的適配

針對網(wǎng)格運動，本方法與貪心選點的RBF 插值法類似，也需要尋找合適的核函數(shù)k(x,xi)和誤差閾值ε，使之具有較強的擬合能力和變形能力。此外，適應各類運動的核函數(shù)和通用的參數(shù)設置方法對于動網(wǎng)格算法的推廣和應用也非常重要。前期相關研究初步采用了應用廣泛的高斯核函數(shù)處理網(wǎng)格運動，相比于貪心選點的RBF 方法，在保持相當網(wǎng)格質量的基礎上計算效率提高了數(shù)倍[21]。本文的一個重要研究內容是通過系統(tǒng)全面地測試，分析各類核函數(shù)在SVR 動網(wǎng)格方法框架下的性能表現(xiàn)，具體內容見第2 節(jié)。

對于SVR 方法的誤差閾值ε，為了減少累積誤差并避免網(wǎng)格發(fā)生交錯，閾值 ε一般應不大于網(wǎng)格點間最短距離Δmin與 變形迭代次數(shù)s的比值，即有：

在此基礎上，誤差閾值 ε越小，SVR 方法的擬合能力越強，但相應的建模效率可能會降低。一般可根據(jù)網(wǎng)格運動的復雜程度，在0 到1 之間粗略設置閾值系數(shù) λ即可。

通常而言，一般的SVR 應用會采用歸一化方法將樣本各個維度的數(shù)據(jù)縮放至[0, 1]范圍內。由于網(wǎng)格坐標各維度采用同一度量單位，且歸一化還會降低數(shù)據(jù)的精度，因此針對網(wǎng)格運動不需要進行這一操作。除此之外，因為所有網(wǎng)格邊界點都是有效數(shù)據(jù)，并且邊界在運動形變中通常是連續(xù)的，故而不存在所謂的噪聲數(shù)據(jù)，也就無需在算法模型中采取松弛變量處理例外[18]。

最后，針對SVR 動網(wǎng)格方法對應的特殊二次規(guī)劃問題，可采用序列最小最優(yōu)化（sequential minimal optimization，SMO）算法[22]進行高效求解。SMO 算法是在每次迭代中只選擇兩個變量而令其他為常數(shù)，從而將原問題分解為若干盡可能小的二次規(guī)劃問題，再采用解析方法直接求解。因此在每個時間步內進行大規(guī)模網(wǎng)格運動時，相對于迭代求解規(guī)模遞增方程組的貪心選點RBF 方法，SVR 方法有著很大的效率優(yōu)勢。此外，SMO 算法從訓練樣本集全局考慮，自適應選擇最合適的邊界點作為支持向量，且無需選擇若干邊界點作為初始的控制點集合，進一步減少了人工干預。

2 基準測試

SVR 動網(wǎng)格方法總體分為兩步，一是利用已知位移的邊界點建立SVR 模型，二是利用該模型預測計算內部點的位移。其中模型建立的時間開銷通常更大，且直接決定第二步的網(wǎng)格質量，因此本節(jié)主要關注第一步的計算效率和對邊界點位移值的擬合效果。對于內部點運動后的最終網(wǎng)格質量評估將在第3 節(jié)通過案例應用測試分析。

2.1 測試案例設計

為測試不同核函數(shù)針對網(wǎng)格邊界運動規(guī)律的擬合能力、計算效率以及不同疏密網(wǎng)格下的健壯性，設計如圖2 所示的三套不同密度的二維笛卡爾網(wǎng)格。這些網(wǎng)格x、y方向的坐標范圍為[-2, 2]，z方向坐標為0，網(wǎng)格單元數(shù)分別為 2 0×20、 4 0×40和 8 0×80。

圖2 三套不同密度的二維笛卡爾網(wǎng)格Fig. 2 Three 2D Cartesian grids with different mesh densities

在實際動網(wǎng)格案例中SVR 方法是分別在三個方向獨立建模，因此研究其擬合能力時只需針對某一方向測試即可。根據(jù)不同類型的網(wǎng)格運動，本節(jié)采用四個相對復雜的解析函數(shù)指導測試網(wǎng)格在z方向的變形，其z坐標的位移值等于用該點坐標代入相應函數(shù)所得。正值表示網(wǎng)格點沿z軸正向運動，負值表示沿z軸反向運動。4 個運動函數(shù)如下：

變形后的網(wǎng)格如圖3 所示，經(jīng)F1 函數(shù)變形后的網(wǎng)格特點是平緩延展、曲面光滑；F2 函數(shù)變形后如對稱的渦流；F3 函數(shù)變形后網(wǎng)格呈現(xiàn)周期性上下起伏；F4 函數(shù)變形則同時具備周期性、激變性、對稱性等多種運動特點。因此，這些網(wǎng)格運動具有很好的代表性。

圖3 四種運動函數(shù)對應變形后的網(wǎng)格Fig. 3 Deformed grids for four motion functions

將初始網(wǎng)格的坐標點和運動變形后的位移量作為不同核函數(shù)下SVR 方法的訓練樣本，采用SMO 算法訓練得到相應的SVR 模型。再將初始網(wǎng)格點代入這些SVR 模型計算得到預測的位移量，從而可通過比較預測值和真實值的誤差、選擇的支持向量（控制點）考察不同核函數(shù)的擬合能力及訓練開銷等性能。

本文采用均方差（mean square error，MSE）來衡量回歸擬合的精確度，其計算公式如下：

其中，zi表 示網(wǎng)格點i的 真實位移值，～zi為預測位移值，Nb為邊界樣本點數(shù)目。MSE 值越大，說明對應SVR模型對網(wǎng)格運動的表示能力越差。

2.2 核函數(shù)類型

對于某個確定的回歸問題，Vapnik 等證明必然存在可行的核函數(shù)，且滿足Mercer 定理的函數(shù)均能作為核函數(shù)[16]。由于滿足Mercer 定理的函數(shù)較多，對于特定類型的問題，通常需要根據(jù)經(jīng)驗或大量實驗對比才能給出更適合的核函數(shù)及其參數(shù)，從而能顯著提高SVR 的建模效率和質量。徑向基核函數(shù)是處理無先驗知識數(shù)據(jù)集首選的核函數(shù)類型，其以空間中兩點的歐氏距離為自變量，影響周圍區(qū)域的樣本取值。

按照作用范圍的不同，徑向基函數(shù)可分為全支撐（Global Support）基函數(shù)和緊支撐（Compact Support）基函數(shù)兩類。全支撐基函數(shù)在整個定義域內均為非零值，其常見的函數(shù)類型如表1 所示，其中MQB 和IMQB 函數(shù)使用參數(shù)a控制基函數(shù)的形狀，本文設置其常用值a=10-3。參考標準高斯核函數(shù)（表1 中的Gauss 函數(shù)），本文進一步將全支撐基函數(shù)的自變量通過系數(shù)r進行縮放，從而控制核函數(shù)對樣本差異的敏感程度，即：

表1 常見的全支撐基函數(shù)Table 1 Common global support basis functions

表2 常見的緊支撐基函數(shù)Table 2 Common compact support basis functions

新的SVR 動網(wǎng)格方法將基于以上核函數(shù)進行測試，更多有關徑向基函數(shù)的信息及其誤差和收斂特性可參閱文獻[23]。

2.3 結果分析

除了前期研究采用的高斯核函數(shù)外，本節(jié)通過大范圍調整核參數(shù)，對其余13 個核函數(shù)進行了大量擬合測試，得到了7 個能較好擬合公式（14）描述的四類網(wǎng)格運動的核函數(shù)及其核參數(shù)設置范圍，結果如表3 所示。其余核函數(shù)均不適用于SVR 動網(wǎng)格方法，后續(xù)不再討論。

表3 成功擬合網(wǎng)格運動的核函數(shù)Table 3 Kernel functions for successfully fitting mesh motions

當SVR 方法采用原來的高斯核函數(shù)時，基于測試案例 40×40的中等網(wǎng)格擬合四類運動，求解回歸方程占總的時間開銷均在85%以上。這應該是由于其為指數(shù)函數(shù)，運算復雜度較高。為評估不同核函數(shù)的擬合效率，本文進一步測試了它們在 40×40網(wǎng)格上的計算開銷，結果如圖4 所示?？梢钥闯觯咚购撕瘮?shù)的計算速度最慢，尤其是在擬合運動特性復雜的F4 函數(shù)時，表現(xiàn)明顯落后于除 CTPSC0以外的其他核函數(shù)。此外，IMQB、IQB 和 C PC2的效率最高。

圖4 不同核函數(shù)在網(wǎng)格上的時間開銷Fig. 4 Time cost on the grid for different kernel functions

為考察邊界點規(guī)模對SVR 方法擬合性能的影響，采用運動特性較全面的F4 函數(shù)測試了核函數(shù)在不同密度網(wǎng)格下（分別用20、40 和80 表示）完成變形的時間開銷、控制點選擇和均方差。圖5 給出了各類核函數(shù)在不同密度網(wǎng)格上的時間開銷，從中可知邊界點規(guī)模增大會明顯降低建模的速度，但不同核函數(shù)的影響差異顯著。高斯核函數(shù)和 CPC6函數(shù)在邊界點增加時效率下降嚴重，而其他核函數(shù)的效率下降相對緩慢。

圖5 不同密度網(wǎng)格上各類核函數(shù)時間開銷Fig. 5 Time cost on different density grids for various kernel functions

圖6 展示了不同核函數(shù)選擇控制點占邊界點數(shù)的比例，可以看出在處理較密兩套網(wǎng)格時高斯核函數(shù)和 CPC6函數(shù)依賴的控制點數(shù)更多，由于SVR 建模速度主要取決于選擇的支持向量數(shù)目，這也解釋了為何這兩個函數(shù)的時間開銷更大。其他核函數(shù)則都能在不改變核參數(shù)的情況下將控制點數(shù)量維持在一定范圍內，這一特性在處理高密度的大規(guī)模網(wǎng)格變形時將具有極大優(yōu)勢。

圖6 不同密度網(wǎng)格上各類核函數(shù)選擇的控制點數(shù)量Fig. 6 Number of control points selected by various kernel functions on different density grids

按照徑向基函數(shù)分類，全支撐組和緊支撐組選擇的控制點分布分別如圖7（a）和圖7（b）所示。從圖中可明顯看出，對于所有核函數(shù)，控制點的總體分布不會隨網(wǎng)格密度增加而有明顯改變。此外，不同核函數(shù)所選擇的控制點雖然分布半徑和密度有所區(qū)別，但基本都吻合同心圓樣式的排列，這正是圖3（d）中F4 函數(shù)運動變化最劇烈的區(qū)域。由此反映出SVR 方法能夠自適應地選擇位移最突出的邊界點來表示邊界運動。

圖7 不同密度網(wǎng)格上各類核函數(shù)的控制點分布Fig. 7 Distribution of control points for various kernel functions on different density grids

圖8 進一步給出了不同核函數(shù)擬合的均方差，說明所有被測核函數(shù)均能有效擬合邊界的運動規(guī)律，但高斯核函數(shù)的擬合能力在密度最大的網(wǎng)格上有所下降，說明其健壯性較差。同時注意到 C PC4和 C PC6相比緊支撐組的其他核函數(shù)，其特性更接近全支撐函數(shù)，且自身階數(shù)較高，較難體現(xiàn)緊支撐核函數(shù)高計算效率的優(yōu)勢。

圖8 不同密度網(wǎng)格上各類核函數(shù)擬合的均方差Fig. 8 MSE of various kernel functions fitting on different density grids

綜合以上分析，針對SVR 動網(wǎng)格方法，全支撐核函 數(shù)IQB、IMQB 和 緊 支 撐 核 函 數(shù) CPC0、 CPC2、CTPSC0在建模效率和健壯性方面都要明顯優(yōu)于高斯核函數(shù)，且 CPC2核函數(shù)的效率最高。這些核函數(shù)對網(wǎng)格質量的影響將在第3 節(jié)做進一步評估。

3 案例應用

本節(jié)通過三個典型案例測試分析核函數(shù)及其參數(shù)設置對SVR 動網(wǎng)格方法的全面影響。對于變形運動后的網(wǎng)格質量，本文采用文獻[12]提出的方法對三角形/四面體單元進行評估，并統(tǒng)計網(wǎng)格平均質量和最差單元質量。該方法通過函數(shù)fss(0 ≤fss≤1)綜合相對尺寸度量和形狀斜交度量評估變形后的網(wǎng)格單元，當函數(shù)值為0 時說明單元已退化、為1 時說明單元為正多邊形/正多面體且其面積/體積保持不變，且對應函數(shù)值越大說明單元質量越高。

3.1 二維矩形塊旋轉平移

本案例主要考察在較大形變運動中采用SVR 方法進行多步變形后的網(wǎng)格質量以及核函數(shù)的影響。初始網(wǎng)格如圖9 所示，其中外邊界尺寸為 25×25，正中的矩形邊界尺寸為 5×1。實驗設置矩形塊經(jīng)由多步沿x和y軸負方向平移5 個單位長度，并繞中心逆時針旋轉 60°。考慮網(wǎng)格變形程度較為嚴重，設置較小的回歸閾值系數(shù) λ =0.4。

圖9 二維矩形塊旋轉平移案例的初始網(wǎng)格Fig. 9 Initial mesh for the 2D case of rectangular block rotation and translation

在進行實際動邊界問題數(shù)值模擬時，從初始到最終狀態(tài)往往需要計算多個時間步，因此也需要經(jīng)過多步將初始網(wǎng)格變形到最終狀態(tài)。首先設總的變形步數(shù)為20，研究不同核參數(shù)對網(wǎng)格變形后質量的影響。圖10 給出了最終網(wǎng)格的最差單元質量度量值，從中可知 CPC0和 CTPSC0函數(shù)的變形質量始終較差，予以排除。三組核函數(shù)各自的核參數(shù)可行范圍并不類似，其中高斯核可行取值范圍最小，全支撐組的范圍最大，此外 CPC2函數(shù)的最差質量在可行核參數(shù)范圍內最穩(wěn)定。

圖10 分組對比不同核參數(shù)對最差單元質量的影響Fig. 10 Group comparison of the worst cell quality by different kernel parameters

為考察網(wǎng)格整體質量，測試了不同步數(shù)下變形到最終狀態(tài)的網(wǎng)格單元最差質量和平均質量，結果如圖11 所示。從中可知 CPC2函數(shù)和所有全支撐核函數(shù)的最差質量均隨著變形步數(shù)增加而增大，且超過一定步數(shù)后質量趨于穩(wěn)定。此外， CPC2函數(shù)的最差質量值是所有函數(shù)最高的，且在步數(shù)超過4 和8 后，IMQB 和IQB 函數(shù)的最差質量也分別超過了對照的高斯核函數(shù)。所有核函數(shù)變形后的網(wǎng)格平均質量度量值都較高，但需注意由于 CPC0和 CTPSC0等函數(shù)未能很好地將邊界運動擴散至周邊網(wǎng)格，導致矩形塊附近的網(wǎng)格單元質量差而其余單元受影響較小保持了原有高質量，所以反而其網(wǎng)格平均質量更高。

圖11 不同核函數(shù)下多步變形后最終的網(wǎng)格質量對比Fig. 11 Comparison of the final mesh quality after multi-step deformation by different kernel functions

圖12 直觀展示了采用IMQB 和高斯核函數(shù)通過4 步變形到最終狀態(tài)的網(wǎng)格及其質量度量分布。在矩形塊經(jīng)過大幅度旋轉平移后，其周圍的網(wǎng)格單元仍與之保持著較好的連接關系，且采用IMQB 核函數(shù)變形后的網(wǎng)格質量明顯更優(yōu)。

圖12 采用4 步變形的最終網(wǎng)格Fig. 12 Final mesh deformed after four steps by IMQB and Gaussian kernel function

3.2 二維多段翼型旋轉

本案例通過多段翼的多體相對運動，進一步研究不同密度的網(wǎng)格變形質量受核函數(shù)的影響。初始網(wǎng)格的遠場邊界為圓形，半徑長為25L，其中L為尾翼弦長，令內部邊界中三段翼型的尾翼沿順時針方向旋轉1 5°，旋轉中心為尾翼根部，其余部件保持靜止。兩套不同密度的初始網(wǎng)格如圖13 所示，其中尾翼與最近的邊界距離僅約0.1L。低密度網(wǎng)格的單元數(shù)為13 120，高密度網(wǎng)格的單元數(shù)為131 772。

圖13 兩套三段翼初始網(wǎng)格（局部）Fig. 13 Two initial grids of the three-segment wing (partial in detail)

本案例網(wǎng)格變形較緩和，取回歸閾值系數(shù) λ=1。此外，考慮到動靜邊界距離很近，應當取較小的核參數(shù)。結合前面實驗的分析，設置高斯核寬度 σ為0.125L，全支撐組的核參數(shù)r為150L，緊支撐組的支撐半徑R為5L。在兩套網(wǎng)格上應用相同的核參數(shù)，不同變形步數(shù)下最終網(wǎng)格的最差單元質量如圖14 所示。首先可以發(fā)現(xiàn)全支撐組函數(shù)和 CPC2的變形能力明顯優(yōu)于高斯核函數(shù)，且緊支撐組的其余兩個函數(shù)的表現(xiàn)明顯好于之前案例，可以推斷是由于動靜邊界距離較小的緣故。

圖14 兩套網(wǎng)格變形后的最差單元質量Fig. 14 Worst cell quality after deformation for two grids

低密度和高密度網(wǎng)格的初始最差單元質量度量值分別為0.71 和0.69，對比兩套網(wǎng)格變形后的質量，可知在同一核參數(shù)下，網(wǎng)格密度的增加令高斯核函數(shù)的變形能力明顯下降，但對其他核函數(shù)影響不大。圖15 分別展示了兩套網(wǎng)格最終狀態(tài)下動邊界附近的網(wǎng)格單元質量分布，結果表明采用 CPC2和IMQB 函數(shù)進行變形后的網(wǎng)格質量明顯優(yōu)于高斯核，且質量較差單元也更少。綜上所述，在多體相對運動的案例中，全支撐組和緊支撐組的核函數(shù)比對照的高斯核函數(shù)具有更好的變形能力和參數(shù)魯棒性。

圖15 三段翼動邊界附近的最終網(wǎng)格質量分布Fig. 15 Final mesh quality distribution near the moving boundary of the three-segment wing

3.3 三維返回艙俯仰振蕩

為進一步測試SVR 方法在三維動網(wǎng)格上的表現(xiàn)，本案例對三維周期性運動后的網(wǎng)格質量穩(wěn)定性進行研究。初始的四面體網(wǎng)格如圖16 所示，計算域的外邊界成圓柱形，底部半徑為5L，高約16L，L為返回艙中軸線長，返回艙最大半徑為0.25L。共有61 514個邊界點，303 927 個內部點。

圖16 返回艙初始網(wǎng)格Fig. 16 Initial mesh of the return capsule

圖17 給出了采用本文方法變形返回艙一個周期的運動示意圖，其俯仰中心位于中軸線尾部五分之一處，并在y軸方向以 20°的振幅做往復運動。為測試網(wǎng)格質量的穩(wěn)定性，令返回艙運動3 個周期，分120 步完成，每10 步完成1/4 個周期并運動至最大俯仰角。

本案例動邊界結構復雜，為減小潛在的插值累積誤差，取較小的回歸閾值系數(shù) λ=0.1，且由于是做周期運動，計算回歸閾值時應取運動到最大形變時的步數(shù)計算，即s=10。此外，由于動靜邊界距離大于運動部件長度L的10 倍，應取較大的核參數(shù)，這里將高斯核寬度 σ設為1.5L，全支撐組的核參數(shù)r為1 200L，緊支撐組的支撐半徑R為5L。

圖18 給出了每一變形步后四面體網(wǎng)格單元的最差質量和平均質量。當達到最大俯仰角時，高斯核函數(shù)相應的網(wǎng)格質量略好于其他核函數(shù)，但相差很小且網(wǎng)格質量均較好。另一方面從圖17 也可直觀看出運動過程中保持了良好的網(wǎng)格質量。此外，之前實驗表現(xiàn)最好的全支撐函數(shù)和 CPC2函數(shù)，其對應網(wǎng)格質量能很好地跟隨返回艙周期運動而表現(xiàn)出周期性變化。不僅各周期對應變形步質量相當，且在三個運動周期后最差網(wǎng)格質量均能回到與初始幾乎相等的水平，說明在整個變形過程中網(wǎng)格質量并未衰退。

圖17 返回艙的周期俯仰運動Fig. 17 Periodic pitching motion of the return capsule

圖18 返回艙各變形步完成后的網(wǎng)格質量Fig. 18 Mesh quality after each deformation step of the return capsule

表4 進一步給出了采用不同核函數(shù)進行變形的總時間開銷，并以高斯核函數(shù)為基準列出了其他函數(shù)所節(jié)省的時間比例。其中采用 CPC2核函數(shù)的計算效率比高斯核提高了68%，其余全支撐函數(shù)也能節(jié)省超過40%的時間。因此，在大規(guī)模網(wǎng)格變形時，應用 CPC2、IMQB 和IQB 核函數(shù)具有更好的性能優(yōu)勢。

表4 不同核函數(shù)針對返回艙案例的時間開銷Table 4 Time cost of different kernel functions for the return capsule case

4 結 論

本文通過進一步放寬貪心選點RBF 插值法的約束條件，提出了一種更高效的基于支持向量回歸的機器學習動網(wǎng)格技術，并以前期研究采用的高斯核函數(shù)為基準，總結了所提方法的核參數(shù)設置規(guī)律。對于常見運動形式，一般可將 CPC2函數(shù)支撐半徑設為3L～5L，IMQB 函數(shù)縮放系數(shù)設為300L～500L，其中L表示運動部件的特征長度。

在包含網(wǎng)格運動的大規(guī)模數(shù)值模擬中，動網(wǎng)格方法的計算效率始終是影響模擬周期的重要因素。在下一步的工作中，可以嘗試將本文算法進行不同模式的并行化改造，使之適用于大規(guī)模復雜問題的工程實踐。