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

(西北工業(yè)大學(xué) 航天學(xué)院，西安 710072)

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

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

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

1 采用徑向基函數(shù)進行網(wǎng)格變形

徑向基函數(shù)的基本形式是［8］:

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

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

式中，下標(biāo)s代表物面插值節(jié)點，

表示Nsp個物面插值節(jié)點上的位移，

是與每一插值節(jié)點相關(guān)的權(quán)重系數(shù)，未知。矩陣中每一元素是以表面插值節(jié)點中任意兩節(jié)點之間距離為參數(shù)的徑向基函數(shù)值，即:

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

2 表面插值節(jié)點的選擇方法

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

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

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

3 空間待插值節(jié)點的減縮方法

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

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

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

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

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

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

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

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

4 算例

4.1 NACA0012翼型變形

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

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

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

4.2 ONERA－M6機翼變形

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

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

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

5 結(jié)論

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

［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］張 軍，譚俊杰，褚 江，等.一種新的非結(jié)構(gòu)動網(wǎng)格生成方法［J］.南京航空航天大學(xué)學(xué)報，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］劉 君，白曉征，郭 正.非結(jié)構(gòu)動網(wǎng)格計算方法-及其在包含運動界面的流場模擬中的應(yīng)用［M］.長沙:國防科技大學(xué)出版社，2009.

［4］霍世慧，王富生，岳珠峰.彈簧近似法在二維非結(jié)構(gòu)動網(wǎng)格生成技術(shù)中的應(yīng)用［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］伍貽兆，田書鈴，夏 鍵.基于非結(jié)構(gòu)動網(wǎng)格的非定常數(shù)值模擬方法［J］.航空學(xué)報.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］王 剛，雷博琪，葉正寅.一種基于徑向基函數(shù)的非結(jié)構(gòu)混合網(wǎng)格變形技術(shù)［J］.西北工業(yè)大學(xué)學(xué)報，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.