基于徑向基函數(shù)的高效網(wǎng)格變形算法研究

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

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

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

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

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

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

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

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

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

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

表示Nsp個(gè)物面插值節(jié)點(diǎn)上的位移，

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

4 算例

4.1 NACA0012翼型變形

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

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

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

4.2 ONERA－M6機(jī)翼變形

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

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

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

5 結(jié)論

本文對(duì)于采用徑向基函數(shù)進(jìn)行動(dòng)態(tài)網(wǎng)格變形技術(shù)中如何減少計(jì)算量的問(wèn)題進(jìn)行了研究，與目前多數(shù)文獻(xiàn)的著眼于表面插值節(jié)點(diǎn)數(shù)不同，本文在減少表面插值節(jié)點(diǎn)數(shù)的時(shí)候也試圖減少空間待插值節(jié)點(diǎn)數(shù)，采用函數(shù)子空間逐級(jí)逼近的基本思想，每一次網(wǎng)格變形后，下一次網(wǎng)格變形的目標(biāo)更改為上一次產(chǎn)生的表面誤差，且每一次網(wǎng)格變形時(shí)的限定區(qū)域隨著誤差的減小而逐漸減小，由此使得空間待插值節(jié)點(diǎn)數(shù)較大時(shí)，表面插值節(jié)點(diǎn)數(shù)較小;表面插值節(jié)點(diǎn)數(shù)較大時(shí)，空間待插值節(jié)點(diǎn)數(shù)較小，從而節(jié)省了總的計(jì)算時(shí)間。算例結(jié)果表明:本文提出的網(wǎng)格變形算法在保證精度的同時(shí)，可以顯著的節(jié)省計(jì)算時(shí)間，并能支持大變形運(yùn)動(dòng)，且變形后的網(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)動(dòng)網(wǎng)格生成方法［J］.南京航空航天大學(xué)學(xué)報(bào)，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)動(dòng)網(wǎng)格計(jì)算方法-及其在包含運(yùn)動(dòng)界面的流場(chǎng)模擬中的應(yīng)用［M］.長(zhǎng)沙:國(guó)防科技大學(xué)出版社，2009.

［4］霍世慧，王富生，岳珠峰.彈簧近似法在二維非結(jié)構(gòu)動(dòng)網(wǎng)格生成技術(shù)中的應(yīng)用［J］.振動(dòng)與沖擊，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］伍貽兆，田書(shū)鈴，夏 鍵.基于非結(jié)構(gòu)動(dòng)網(wǎng)格的非定常數(shù)值模擬方法［J］.航空學(xué)報(bào).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é)報(bào)，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.