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

高 翔，廖海翔，徐傳福

(國防科技大學(xué) 計(jì)算機(jī)學(xué)院，量子信息研究所兼高性能計(jì)算國家重點(diǎn)實(shí)驗(yàn)室，長沙 410073)

0 引 言

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

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

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

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

1 方法描述

1.1 貪心選點(diǎn)RBF 方法約束條件的放寬

RBF 動網(wǎng)格方法在各坐標(biāo)方向分別采用一個RBF 函數(shù)，以描述整個計(jì)算域內(nèi)任意點(diǎn)在該方向?qū)⒁苿拥奈灰屏?，一般的表達(dá)形式如下：

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

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

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

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

1.2 支持向量回歸

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

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

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

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

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

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

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

關(guān)于支持向量回歸對偶轉(zhuǎn)化和非線性擴(kuò)展更詳細(xì)的推導(dǎo)，可參閱文獻(xiàn)[21]。

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

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

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

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

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

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

2 基準(zhǔn)測試

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

2.1 測試案例設(shè)計(jì)

為測試不同核函數(shù)針對網(wǎng)格邊界運(yùn)動規(guī)律的擬合能力、計(jì)算效率以及不同疏密網(wǎng)格下的健壯性，設(shè)計(jì)如圖2 所示的三套不同密度的二維笛卡爾網(wǎng)格。這些網(wǎng)格x、y方向的坐標(biāo)范圍為[-2, 2]，z方向坐標(biāo)為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

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

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

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

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

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

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

2.2 核函數(shù)類型

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

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

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

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

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

2.3 結(jié)果分析

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

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

當(dāng)SVR 方法采用原來的高斯核函數(shù)時，基于測試案例 40×40的中等網(wǎng)格擬合四類運(yùn)動，求解回歸方程占總的時間開銷均在85%以上。這應(yīng)該是由于其為指數(shù)函數(shù)，運(yùn)算復(fù)雜度較高。為評估不同核函數(shù)的擬合效率，本文進(jìn)一步測試了它們在 40×40網(wǎng)格上的計(jì)算開銷，結(jié)果如圖4 所示?？梢钥闯觯咚购撕瘮?shù)的計(jì)算速度最慢，尤其是在擬合運(yùn)動特性復(fù)雜的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

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

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

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

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

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

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

圖8 進(jìn)一步給出了不同核函數(shù)擬合的均方差，說明所有被測核函數(shù)均能有效擬合邊界的運(yùn)動規(guī)律，但高斯核函數(shù)的擬合能力在密度最大的網(wǎng)格上有所下降，說明其健壯性較差。同時注意到 C PC4和 C PC6相比緊支撐組的其他核函數(shù)，其特性更接近全支撐函數(shù)，且自身階數(shù)較高，較難體現(xiàn)緊支撐核函數(shù)高計(jì)算效率的優(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)格質(zhì)量的影響將在第3 節(jié)做進(jìn)一步評估。

3 案例應(yīng)用

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

3.1 二維矩形塊旋轉(zhuǎn)平移

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

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

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

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

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

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

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

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

3.2 二維多段翼型旋轉(zhuǎn)

本案例通過多段翼的多體相對運(yùn)動，進(jìn)一步研究不同密度的網(wǎng)格變形質(zhì)量受核函數(shù)的影響。初始網(wǎng)格的遠(yuǎn)場邊界為圓形，半徑長為25L，其中L為尾翼弦長，令內(nèi)部邊界中三段翼型的尾翼沿順時針方向旋轉(zhuǎn)1 5°，旋轉(zhuǎn)中心為尾翼根部，其余部件保持靜止。兩套不同密度的初始網(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。此外，考慮到動靜邊界距離很近，應(yīng)當(dāng)取較小的核參數(shù)。結(jié)合前面實(shí)驗(yàn)的分析，設(shè)置高斯核寬度 σ為0.125L，全支撐組的核參數(shù)r為150L，緊支撐組的支撐半徑R為5L。在兩套網(wǎng)格上應(yīng)用相同的核參數(shù)，不同變形步數(shù)下最終網(wǎng)格的最差單元質(zhì)量如圖14 所示。首先可以發(fā)現(xiàn)全支撐組函數(shù)和 CPC2的變形能力明顯優(yōu)于高斯核函數(shù)，且緊支撐組的其余兩個函數(shù)的表現(xiàn)明顯好于之前案例，可以推斷是由于動靜邊界距離較小的緣故。

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

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

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

3.3 三維返回艙俯仰振蕩

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

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

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

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

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

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

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

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

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

4 結(jié) 論

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

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