基于徑向基函數的高效網格變形算法研究

謝 亮，徐 敏，張 斌，安效民

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

在非定常流動模擬、流固耦合時域仿真、氣動外形優化中，為了使得計算網格適應物面邊界的變化，可以采用網格重構的方法，然而更廣泛的是應用網格變形技術，這樣可以保留原始網格的拓撲結構與密度、正交性等特性，不會在求解器中引入額外的誤差，同時計算量較小。目前存在多種網格變形方法，無限代數插值方法［1］計算量較小，然而卻僅適用于結構網格，且難于處理復雜的拓撲結構，通用性不好。對于非結構網格、混合網格，可以應用彈性體方法［2］、彈簧比擬法［3-4］以及Delaunay圖映射［5］方法。彈性體方法將網格所占據的區域比擬為一個彈性體，通過采用有限元或者邊界元方法求解彈性體的變形來實現動態網格變形技術，通用性較好，效果優良，變形后的網格質量極好，尤其適用于大變形情況，然而計算量非常龐大。彈簧比擬法將網格線比擬為彈簧，通過求解彈簧變形后的平衡位置實現網格變形，其通用性較好，適用于任意網格，網格變形后的質量也較好，然而其計算量仍然較大，且難于處理大變形問題。Delauny圖映射方法魯棒性與效率較好，但是對于大變形問題仍然難以處理。

近些年來，由于具備諸多優良的特性，基于徑向基函數的網格變形技術得到了極大的發展［6-9］與廣泛的應用［10-12］。其主要思想是利用已知的物面網格變形，構造一個徑向基函數序列，再使用此徑向基函數序列將物面變形光滑地插值到空間氣動網格上去。此方法計算過程中不需要利用網格節點之間的聯系，各個計算結點的變形計算是完全不相干擾的，因此非常便于并行化;可支持復雜外形的大變形運動，變形后的網格質量較好，通用性與魯棒性極好;編程也較方便，易于實現;物面與空間網格節點的處理采用一樣的方式，因而不需要額外的數據插值方法;既可用于結構網格，也可用于非結構與混合網格;可支持任意形式的變形，剛體運動、彈性變形甚至塑性變形都可支持。然而，采用徑向基函數進行網格變形時的計算量與Nvp×Nsp成正比，其中Nvp是待插值的氣動空間網格點數，其量級一般在105～106，Nsp是用于插值的表面節點數，其量級一般在102～103，相較于代數插值方法，其計算量是相當大的。在過去的數年間，為了減少計算量，不少研究員［6-8］做出了相當多的努力。其中，Rendall等［6-7］采用貪心算法，在物面插值過程中將出現最大誤差值的點納入插值節點集合中，從而在保證精度的前提下選擇得到盡可能少的物面插值節點，以此來減少總的計算時間。最近，王剛等［8］將Rendall提供的算法加以改進，引入子空間逐級逼近的思想，每一次插值過程中選擇一插值節點集合，下一次插值時以上一次插值過程中在物面產生的誤差作為插值對象，從而將選擇插值節點的時間大為減少。然而目前來看，為了達到減少計算量，幾乎所有的文獻都集中在減少表面插值節點的數目上，而少有人考慮在插值過程中減少空間待插值節點的方法，僅見有Andreas［9］為了分離飛行器部件運動對其它部件的影響，提出過限制徑向基函數插值區域的算法。

本文受王剛在選擇表面插值節點過程中引入子空間逼近的思想的啟發，同時將Andreas提出的限制徑向基函數插值區域的方法加以改進，兩者結合，提出了在計算過程中逐級減少空間待插值節點的方案，此方案不損失插值精度，且支持大變形。算法的關鍵在于在選擇表面插值節點個數的時候，隨著插值節點個數的逐步增加，插值誤差會逐步減少，當以前一步插值過程中產生的插值誤差為插值對象的時候，此時的插值限制區域可以取得較小，這樣就可以使得當Nvp比較大時，Nsp比較小，而Nvp比較大時，Nsp比較小，從而減少總的計算量。算例結果表明本文提出的方案可以顯著地提高采用徑向基函數進行網格變形時的計算效率，同時保留對大變形問題的支持及插值精度。

1 采用徑向基函數進行網格變形

徑向基函數的基本形式是［8］:

式中:F(r)是插值函數，在網格變形問題中，它代表網格變形量;Nsp代表插值問題所使用的徑向基函數的總數目，在網格變形問題中，它等于用于插值的物面節點數;φ(‖r-ri‖)是徑向基函數的一般形式，ri是第i個物面插值節點的位置，r是空間任意一點的位置矢量，進行網格變形時，它就是CFD網格空間節點的位置矢量，‖r-ri‖是空間任意一點到第i個物面插值節點的距離;ωi是與第i個插值節點相對應的權重系數。徑向基函數有許多種，文獻［6］列出了數種形式，經比較認為Wendland's C2函數計算效率與網格變形的質量都較好，故而本文采用它作為徑向基函數來實施網格變形，其形式如下［6，8］:

給定物面插值節點及該點上的位移后，采用式(1)、(2)求解空間任意一點的變形時，r，ri已知，且徑向基函數φ(η)已經給定，故而唯一未知的是與各個物面插值節點相關的權重系數ωi，該值可以使用物面上的插值結果必須與給定位移相吻合這一條件來求得，即求解下述方程即可:

式中，下標s代表物面插值節點，

表示Nsp個物面插值節點上的位移，

是與每一插值節點相關的權重系數，未知。矩陣中每一元素是以表面插值節點中任意兩節點之間距離為參數的徑向基函數值，即:

求解方程(3)～(5)后可得到權重系數，然后根據式(1)可求得空間任意位置上的網格變形。

2 表面插值節點的選擇方法

采用徑向基函數進行動態網格變形時，計算量與Nvp×Nsp成正比，其中Nvp是待插值的氣動空間網格點數，Nsp是用于插值的表面節點數，在流固耦合時域仿真時，它一般取自結構模態。為了減少總的計算量，當前主流做法是減少Nsp的值，Rendall等［6-7］采用貪心法根據最大插值誤差位置逐步添加插值節點的方法來實現徑向基函數序列的精簡。其基本過程是:

首先，任意選擇Ⅰ(一般可取Ⅰ=3)個物面節點形成初始節點集合P0=(p1，p2，…，pI)，采用此集合進行徑向基函數插值，通過求解方程得到相應的權重系數，然后求得所有物面節點上的網格變形，顯然這樣建立的初始插值函數對于P0中的所有節點是精確的，但是對于不屬于P0中的物面節點將產生誤差，確定出現最大誤差的物面節點，根據貪心法的原則，將此節點納入P0中形成下一個節點集合P1，再次采用此節點集合P1進行徑向基函數插值。反復執行此過程直到物面節點上的插值誤差滿足事先給定的誤差限。在選擇物面節點過程中亦可以每插值一次，將誤差大于平均誤差或者給定的一個限值的所有物面節點都選擇到插值節點集合中;或者混合上述兩種方法，每n步交替使用此兩種方法，n步內選擇誤差最大的節點，第n+1步選擇誤差較大的一些節點，以加快選擇節點的速度。

最近，為了加快數據精簡的效率，王剛等［8］提出了一套新的數據精簡方法，其基本思想是按照貪心法的原則，引入函數空間子集逐級逼近的基本思想，具體做法是先選擇N0個節點進行徑向基函數插值，得到所有物面節點上的誤差ΔS(0)，然后將徑向基函數的插值對象由最初的網格變形更改為當前物面節點的誤差ΔS(0)，再次運用貪心法選取N1個物面節點進行插值。重復此步驟直到殘差滿足要求，最后將此n步選擇得到的插值節點和權重系數迭加，得到最終的徑向基函數插值系數。采用此種方法的目的在于加快數據精簡的效率，因為在計算插值系數的過程中，每一次都要求解一個Ni×Ni的線性代數方程組，通過此種算法，可以將每一次計算插值系數過程中的Ni值控制在一個較小量上，由此加快數據精簡過程。然而我們將看到，將函數空間子集逐級逼近的思想加以擴展，可以得到一種減縮空間待插值節點的方法。

3 空間待插值節點的減縮方法

目前，為了減少采用徑向基函數進行網格變形的計算量，研究人員的著眼點都集中在減少Nsp即表面插值節點上，但是對于在實施徑向基函數過程中如何減少Nsp，卻少有人提及。目前僅見有 Andreas［9］為了分離單個部件運動對其它部件的影響提出了限制徑向基函數插值區域的算法。其算法的基本思想是將部件的運動限制區域以長方體描述，該長方體包括了運動部件區域，運動部件上的變形以給定為要求值，長方體表面上的變形給定為零，在長方體表面與運動部件上選擇插值節點后進行徑向基插值，長方體外的空間網格不作插值。

將函數空間子集逐級逼近加以擴展，并將Andreas提出的限制插值區域的方法加以改進，兩者相結合，得到了一種新的效果優良的減縮Nvp的方法，可支持大變形運動。方案主要步驟為:

(1)確定插值限定區域Ω0，不同于文獻中的做法，本文選擇到物面距離小于一限制值R0=k·ΔSmax的點的集合作為Ω0，其中k為一系數，一般取k≥5。并確定Ω0中除去物面的邊界，可按下列方法確定此邊界，掃描Ω0內各點，如果某點周邊單元上的節點中有不屬于Ω0的節點，則此點為邊界，否則是內點，將邊界上的位移值給定為0。給定一個初始誤差限E0，此誤差可取得較大，比如1.0E-2，選擇N(0)個插值節點使得表面插值誤差的最大值小于此給定誤差限，并記錄表面插值誤差ΔS0。在區域Ω0內進行插值，此時，Nsp較大，Nsp=N(0)，由于誤差限取得較大，求得的N(0)會比較少，即總計算量會比較少。

(2)再以上一步的表面插值誤差ΔS0為插值對象，確定插值限定區域Ω1。由于此時表面插值量最大位移才不到E0，同樣按到物面距離小于R1=k·ΔS0，max為標準，由于此時 ΔS0，max≤E0，故而此時的 Ω1可以取得較小，從而使得Nvp較小，這時可以給定一較小的誤差限E1來選擇N(1)個插值點，由于誤差較小，故選擇得到的N(1)較大。然而總的計算量仍然會比較小。

(3)重復上述步驟直到誤差限滿足要求。

需要注意的是，當Ri較小的時候，插值區域限定過小，區域外邊界與物面距離過近，使得選擇的插值節點中有些節點之間距離過近，會造成求解插值系數時系數矩陣十分病態，因此，實際上對于Ri有一限制，本文限定為Ri=k1·max(Ei-1，dmax)，其中k1為一保險系數，一般取為2，Ei-1為上一步的誤差限，而dmax為與物面相連的最長網格線的長度。加上此限制之后，一般第二次插值時Ri即達到限制值了，因此，第三步就用不著了。同時應當注意的是，在插值區域的外邊界上選擇插值節點時，并不要求在此邊界上的誤差與物面上保持同樣小的誤差，僅需此邊界上的誤差不造成該邊界附近網格質量出現較大下降，因此，該邊界上插值節點的選擇以保證該邊界上最大誤差小于附近單元尺寸最小值的k倍，k可取0.1。這樣可以避免插值過程中Nsp較大。

(4)如果最后一次插值中誤差限已經衰減到遠小于空間網格點到物面的最小距離了，此時可選擇更多的插值點，但是僅在物面進行插值，進一步減小物面插值的誤差，因為此時物面的變形不會造成空間網格的扭曲。這時Nvp就等于氣動網格的物面網格數了，這一步計算量會比較小，然而可以較大程度地減小整個物面的插值誤差。

在上述過程中，確定插值限定區域的計算量較大，但是此過程僅實施一次，之后每一次網格變形時的限定區域都與初始確定的區域一致，因此這一步的計算時間并不增加多次執行的網格變形所需要的時間。

4 算例

4.1 NACA0012翼型變形

取二維NACA0012結構網格作為初始網格，變形位移方程為 Δy=0.03sin(4πx)［8］。網格大小為 323 ×81，取誤差限為1E-6，分別采用單次插值與本文所提出的多次插值方法進行計算，在計算過程中，多次插值中第一次插值由于本來變形比較大，且遠場網格本身比較稀疏，簡單起見，直接將全場網格當做Ω0，在第二次插值中減縮空間網格點數，選擇Ω1中包含3 553個網格節點，第二次插值計算完畢后，物面插值誤差限已經足夠小，因此第三次插值只計算物面網格變形，不用計算Ω2。插值效率的比較見表1，可見本文提出的方法確實較好地提高了徑向基函數的插值效果，對于多次插值中的第二次插值，雖然表面插值節點較多，但是由于空間待插值節點大幅減少，效率仍然得到提高，而第三步由于空間插值節點數退化為表面插值節點數，效率提高得更加明顯，從而印證了本算法思路的正確性。網格變形結果見圖1。

圖1 NACA0012翼型變形前后對比Fig.1 Comparison of deforming grid for naca0012 foil

表1 完成NACA0012翼型網格變形所需時間的比較Tab.1 Comparison of efficiency between one interpolation and multi-level interpolation for grid-deforming of naca0012 foil

4.2 ONERA－M6機翼變形

取ONERA-M6機翼非結構網格為初始網格，機翼幾何形狀與參數見文獻［13］，網格單元數為1 883 014，結點數為336 684，物面節點數為33 206，給定以下變形模式:其中:b為機翼展長，z為z方向坐標，Δy為y方向變形。取最終誤差限為1.0E-6，分別采用單次插值與本文所提出的多次插值方法進行計算，計算過程與上例相同。計算過程所花費的時間與單次插值效率的比較見表2，由表可知本文所提出的多次插值方法可以較大程度上減少空間待插值網格數，從而顯著地提高了網格變形效率。初始網格與變形后網格的比較如圖2所示。圖2(a)，(b)對比了初始的與變形后的表面網格，可見式(6)所給出的變形幅度是相當大的。圖2(c)，(d)分別給出了機翼附近空間網格變形前后的情況，圖2(e)，(f)分別給出了翼梢附近網格變形前后的細節，可見變形后的網格基本保持了變形前的密度分布，且質量較好。經流場求解程序的檢驗，變形后的網格質量沒出現問題。

圖2 ONERA-M6機翼網格變形前后對比圖Fig.2 Comparison between original and deforming grid of ONERA-M6 wing

表2 完成ONERA－M6機翼網格變形所需時間的比較Tab.2 Comparison of efficiency betweenone interpolation and multi－level interpolation for grid-deforming of ONERA－M6 wing

5 結論

本文對于采用徑向基函數進行動態網格變形技術中如何減少計算量的問題進行了研究，與目前多數文獻的著眼于表面插值節點數不同，本文在減少表面插值節點數的時候也試圖減少空間待插值節點數，采用函數子空間逐級逼近的基本思想，每一次網格變形后，下一次網格變形的目標更改為上一次產生的表面誤差，且每一次網格變形時的限定區域隨著誤差的減小而逐漸減小，由此使得空間待插值節點數較大時，表面插值節點數較小;表面插值節點數較大時，空間待插值節點數較小，從而節省了總的計算時間。算例結果表明:本文提出的網格變形算法在保證精度的同時，可以顯著的節省計算時間，并能支持大變形運動，且變形后的網格質量仍然合乎要求。

［1］ Li J，Huang S Z，Jiang S Q，et al.Unsteady viscous flow simulations by a fully implicit method with deforming grid［M］.AIAA 2005-1221.

［2］張 軍，譚俊杰，褚 江，等.一種新的非結構動網格生成方法［J］.南京航空航天大學學報，2007，39(5):633-636.

ZHANG Jun，TAN Jun-jie，CHU Jiang，et al.New method for generating unstructured moving grids［J］.Journal of Nanjing University of Aeronautics ＆ Astronautics，2007，39(5):633-636.

［3］劉 君，白曉征，郭 正.非結構動網格計算方法-及其在包含運動界面的流場模擬中的應用［M］.長沙:國防科技大學出版社，2009.

［4］霍世慧，王富生，岳珠峰.彈簧近似法在二維非結構動網格生成技術中的應用［J］.振動與沖擊，2011，30(10):177-182.

HUO Shi-hui， WANG Fu-sheng， YUE Zhu-feng. Spring analogy method for generating of 2D unstructured dynamic meshes［J］.Journal of Vibration and Shock，2011，30(10):177-182.

［5］伍貽兆，田書鈴，夏 鍵.基于非結構動網格的非定常數值模擬方法［J］.航空學報.2011:32(1):15-26.

WU Yi-zhao，TIAN Shu-ling，XIA Jian.Unstructured grid methods for unstready flow simulation［J］.Acta Aeronautica et Astronautica Sinica，2011，32(1):15-26.

［6］ Rendall T C S，Allen C B.Efficient mesh motion using radial basis functions with data reduction algorithms［J］.Journal of Computational Physics，2009，229(7):6231-6249.

［7］Rendall T C S，Allen C B.Reduced surface point selection options for efficient mesh deformation using radial basis functions［J］.JournalofComputationalPhysics，2010，229(8):2810-2820.

［8］王 剛，雷博琪，葉正寅.一種基于徑向基函數的非結構混合網格變形技術［J］.西北工業大學學報，2011，29(5):783-788.

WANG Gang， LEI Bo-qi， YE Zheng-yin. An efficient deformation tachnique for hybrid unstructured grid using radialbasis functions［J］. Journal of Northwestern Polytechnical University，2011，29(5):783-788.

［9］Andreas K M.Aircraft control surface deflection using RBFBased Mesh deformation.

［10］ Beckert A，Wendland H.Multivariate interpolation for fluitstructure-interaction problems using radial basis functions.

［11］de Boer A，van der Schoot M S，Bijl H.Mesh deformation based on radial basis function interpolation［J］.Computers ＆Structures，2007，85:784-795.

［12］ Rendall T C S，Allen C B.Parallel efficient mesh motion using radial basis functions with application to multi-bladed rotors［J］.International Journal for Numerical Methods in Engineering.2010，81:89 ～105.

［13］Schmitt V，Charpin F.Pressure Distributions on the ONERAM6-Wing at Transonic Mach Numbers［R］.Experimental Data Base for Computer Program Assessment.Report of the Fluid Dynamics Panel Working Group 04，AGARD AR 138，May，1979.