靳晨暉,王剛,陳鑫,周豪
(西北工業(yè)大學(xué) 航空學(xué)院,西安 710072)
嵌套網(wǎng)格技術(shù)在航空航天領(lǐng)域有著廣泛的應(yīng)用,針對一些復(fù)雜的幾何構(gòu)型,使用嵌套網(wǎng)格技術(shù)可以分區(qū)生成網(wǎng)格,減少網(wǎng)格生成的難度。對于飛機副油箱分離[1],座艙蓋分離/座椅彈射[2],多級火箭助推器分離[3]等多體分離問題,不同物面間的相對位置會發(fā)生變化,同時伴有強烈的非定常效應(yīng)。常規(guī)研究方法如風(fēng)洞和飛行實驗常帶有一定的局限性,使用嵌套網(wǎng)格技術(shù)進(jìn)行非定常計算可以對這些過程進(jìn)行精確數(shù)值模擬。
在嵌套網(wǎng)格算法中,一般通過插值實現(xiàn)嵌套區(qū)域的流場信息傳遞,尤其是在解決帶有動邊界的問題,嵌套區(qū)域的網(wǎng)格在每個非定常步內(nèi)都需要更新,插值模板也需要重新選取。因此,使用正確的插值方法是保證計算結(jié)果準(zhǔn)確的關(guān)鍵。按照插值方法的不同可分為守恒型插值與非守恒型插值。對于守恒型插值而言,網(wǎng)格尺度對其插值精度的影響較小,但插值過程十分復(fù)雜;而非守恒型插值的實現(xiàn)則較為簡單,但需要控制嵌套區(qū)域網(wǎng)格的匹配程度,否則會損失計算精度。為了減少非守恒型插值過程中的精度損失,一般需要從兩方面入手:一方面要改善嵌套插值區(qū)域網(wǎng)格匹配程度,如在網(wǎng)格生成時,需要對背景網(wǎng)格中有相對運動的區(qū)域進(jìn)行局部加密,同時使用高效的嵌套裝配算法[4]劃分嵌套區(qū)域,以此來保證物面運動的整個過程中插值單元與貢獻(xiàn)單元的匹配程度。另一方面需要研究不同插值方法的插值效果,C.Mastin等[5]研究了雙線性、三線性插值方法在嵌套區(qū)域的插值效果;田書玲等[6]研究了一種基于Lagrange插值的非結(jié)構(gòu)嵌套網(wǎng)格插值方法;周乃春等[7]使用了逆向距離權(quán)方法對嵌套區(qū)域進(jìn)行插值;黃宇等[8]使用了具有二階精度的插值方法,有效降低了由于在嵌套區(qū)域進(jìn)行插值而引入的誤差。僅僅增加嵌套區(qū)域的網(wǎng)格量有時也難以達(dá)到較好的精度,網(wǎng)格量增加會導(dǎo)致仿真周期成倍增加,占用過多計算資源;同時,線性插值方法在處理復(fù)雜問題時的精度相對較低,選擇合適的插值方法對于處理不同的問題至關(guān)重要。
徑向基函數(shù)插值方法因其存儲十分方便和數(shù)據(jù)結(jié)構(gòu)相對簡單的特點,在近二十年內(nèi)迅猛發(fā)展。研究者們在不同的領(lǐng)域?qū)较蚧瘮?shù)的應(yīng)用進(jìn)行大量研究,如R.Franke[9]發(fā)現(xiàn)對散亂的數(shù)據(jù)進(jìn)行擬合時,使用徑向基函數(shù)進(jìn)行插值的精度要高于其他插值方法;林言中等[10]將徑向基函數(shù)插值方法應(yīng)用于網(wǎng)格變形技術(shù)中的界面數(shù)據(jù)傳遞,對AGARD 445.6機翼顫振問題進(jìn)行非定常計算;劉溢浪等[11]提出一種增量形式的RBF插值方法,并將其應(yīng)用于有限體積方法的流場重構(gòu)中;Wang G等[12-13]在RBF網(wǎng)格變形技術(shù)的基礎(chǔ)上,提出將拉普拉斯光順與網(wǎng)格變形相結(jié)合,有效改善了變形后的網(wǎng)格質(zhì)量。但是,徑向基函數(shù)插值半徑的選取帶有一定的經(jīng)驗性,傳統(tǒng)徑向基函數(shù)前述插值半徑選取不合適會嚴(yán)重影響插值的精度與效率。
為了解決嵌套區(qū)域的插值精度問題,本文對不同插值函數(shù)進(jìn)行研究,同時采用一種改進(jìn)的非結(jié)構(gòu)嵌套網(wǎng)格徑向基函數(shù)(RBF)插值方法。在嵌套網(wǎng)格流場信息傳遞的過程中,除了使用線性插值(LINE)和以距離倒數(shù)作為權(quán)重插值(WA)兩種傳統(tǒng)插值方法,還對逆二次徑向基函數(shù)(IQ)和Wendland’s C2徑向基函數(shù)(C2)的插值效果進(jìn)行研究。針對嵌套邊界區(qū)域出現(xiàn)的數(shù)值奇性現(xiàn)象,通過控制徑向基函數(shù)的作用半徑來調(diào)整插值矩陣條件數(shù)的方式進(jìn)行改善。
為便于對固壁邊界的六自由度運動進(jìn)行描述,采用ALE(Arbitrary Lagrangeian-Eulerian)方法對流動控制方程進(jìn)行描述。ALE方法允許計算網(wǎng)格的進(jìn)行剛體運動和變形,通過在流動控制方程中加入網(wǎng)格運動速度,將流體力學(xué)中的Lagrange方法和Euler方法進(jìn)行統(tǒng)一描述。使用ALE方法描述的三維非定常雷諾平均Navier-Stokes方程(URANS)的積分形式為
(1)
式中:Ω為控制體;?Ω為控制體單元邊界;Q為守恒變量,Q=[ρρuρvρwρE]T,其中ρ為流體密度,u、v、w分別為旋轉(zhuǎn)彈在體軸系下三個軸向的運動速度E為總內(nèi)能;n為控制體單元邊界上沿外法線方向的分量;F(Q)和G(Q)分別為N-S方程中的無粘通量項和有粘通量項。
對于式(1),使用格心有限體積法(FVM)進(jìn)行空間離散后,得到半離散形式為
(2)
(3)
式中:R為計算殘差;S為網(wǎng)格單元的面積;n為時間步;Vgrid為網(wǎng)格進(jìn)行剛體運動的速度。為了保證足夠的時間精度,使用二階向后的歐拉隱式方法[14]求解式(2),如式(4)所示
(4)
式中:Δt為時間步長;上標(biāo)n為真實時間迭代的步數(shù)。直接求解式(4)較為復(fù)雜,一般采用偽時間迭代的方法對其進(jìn)行求解,為此在方程的左端加入一個守恒變量對虛擬時間τ的導(dǎo)數(shù),從而方程的解可以等價為求解一個關(guān)于虛擬時間τ的一階常微分方程組在虛擬時間τ趨近無窮時的漸進(jìn)解:
(5)
(6)
由于式(5)的精度與偽時間迭代方法無關(guān),因此,偽時間迭代過程可以使用一些加速收斂方法來提高流場的計算效率。具體方法可以參照文獻(xiàn)[15-16]。采用改進(jìn)的LU-SGS[17]算法對式(5)進(jìn)行向后的歐拉隱式時間迭代,同時使用基于OpenMP并行算法來提高計算的效率。
嵌套網(wǎng)格裝配技術(shù)的兩個核心問題是如何確定嵌套邊界和如何查找貢獻(xiàn)單元。確定嵌套邊界也叫 “挖洞”,當(dāng)計算網(wǎng)格量較大時,會顯著增加預(yù)處理的時間,因此高效并且合理地劃分出嵌套區(qū)域?qū)⑹乔短拙W(wǎng)格裝配技術(shù)的關(guān)鍵。本文在嵌套區(qū)域劃分時采取的策略是根據(jù)壁面距離來確定嵌套邊界區(qū)域,同時計算子網(wǎng)格到自身物面和其他物面的最短距離,通過比較到不同物面的最短距離來確定嵌套邊界,相比于傳統(tǒng)的枚舉法,有效減少了前處理時間。節(jié)點屬性判定如圖1所示,i點是屬于A網(wǎng)格的網(wǎng)格節(jié)點,在求解距離時,距離A網(wǎng)格更近,則i點的屬性標(biāo)記為激活的狀態(tài);而j點距離B網(wǎng)格更近,則j點的屬性標(biāo)記為非激活的狀態(tài)。
圖1 節(jié)點屬性判定示意圖Fig.1 Schematic diagram of node attribute determination
當(dāng)確定了嵌套邊界后,需要在嵌套區(qū)域查找貢獻(xiàn)單元。對于復(fù)雜模型來講,網(wǎng)格量較大,查找貢獻(xiàn)單元往往要花費較多時間,因此如何快速的查找到貢獻(xiàn)單元是提高嵌套效率的關(guān)鍵。本文根據(jù)文獻(xiàn)[4]采用了一種基于非結(jié)構(gòu)嵌套網(wǎng)格的“蛇形查找方法”(如圖2所示),在空間中以近乎最佳路徑的方式查找到了貢獻(xiàn)單元,有效縮短了查找貢獻(xiàn)單元的時間與距離。假設(shè)要尋找i點的貢獻(xiàn)單元,該方法的具體實現(xiàn)過程如下:
(1) 在空間中賦予i點一個初始的出發(fā)單元,這里假設(shè)是1號網(wǎng)格單元;
(2) 連接當(dāng)前的出發(fā)單元的網(wǎng)格格心和i點,判斷其穿過當(dāng)前出發(fā)單元的哪一條邊(三維為面);
(3) 若穿過某一邊(面),則將該邊(面)的另一側(cè)單元當(dāng)做出發(fā)單元,循環(huán)第二步;
(4) 若沒有穿過當(dāng)前的出發(fā)單元所有的邊(面),則當(dāng)前的出發(fā)單元就是i點的貢獻(xiàn)單元,結(jié)束查找。
圖2 蛇形查找示意圖Fig.2 Schematic diagram of snake search
通過以上步驟,嵌套區(qū)域的網(wǎng)格裝配完畢。在進(jìn)行流場計算時,網(wǎng)格實際被劃分為三個區(qū)域,分別是計算區(qū)域、關(guān)閉區(qū)域和嵌套插值區(qū)域。在計算區(qū)域,進(jìn)行正常的CFD求解,當(dāng)接觸到邊界時,面外側(cè)位于嵌套區(qū)域的網(wǎng)格信息由貢獻(xiàn)單元插值得到;在關(guān)閉區(qū)域,只計算網(wǎng)格單元體積和網(wǎng)格單元面積等預(yù)處理操作;在嵌套插值區(qū)域,在每個偽時間迭代結(jié)束后進(jìn)行插值,嵌套插值區(qū)域的網(wǎng)格不進(jìn)行計算,而是被當(dāng)作嵌套邊界來處理,這類邊界具有插值后獲得的流場信息。
流場信息的插值傳遞方法較多,本文主要考察了線性插值、距離倒數(shù)權(quán)重插值兩種傳統(tǒng)插值方法,以及逆二次徑向基函數(shù)、Wendland’s C2徑向基函數(shù)插值方法。由于線性插值結(jié)構(gòu)較為簡單,此處不再贅述,下面重點對其余三種插值方法進(jìn)行簡要介紹。
以距離倒數(shù)作為權(quán)重來進(jìn)行插值的原理是:對待插值點附近點進(jìn)行加權(quán)平均,其中權(quán)重的大小與需要插值的點與相鄰點之間的距離密切相關(guān),與距離呈反比,因此也叫做距離反比插值方法。其具體形式為
(7)
本文所使用的徑向基函數(shù)的基本形式為
(8)
式中:下標(biāo)i為支撐點;r為該點的位置矢量;N為網(wǎng)格節(jié)點的總數(shù);φ(‖rp-ri‖)為所采用的基函數(shù)的形式;ωi為與第i個支撐點相關(guān)的權(quán)重系數(shù)。徑向基函數(shù)大致可分為以下三類:全局型(Global)徑向基函數(shù)、局部型(Local)徑向基函數(shù)、緊支型(Compact)徑向基函數(shù),這里簡單介紹一下本文所使用的全局型/局部型(IQ)與緊支型(C2)兩種徑向基函數(shù)。
(1) 逆二次徑向基函數(shù)方法
(9)
在對嵌套區(qū)域的P點進(jìn)行流場信息傳遞時,根據(jù)之前查找得到的貢獻(xiàn)單元,將貢獻(xiàn)單元以及它的相鄰單元一起作為徑向基函數(shù)的插值模版構(gòu)建徑向基函數(shù)的系數(shù)矩陣。插值模版示意圖如圖3所示,i網(wǎng)格是P的貢獻(xiàn)單元,則使用i、j、k、l四個網(wǎng)格構(gòu)建徑向基函數(shù)系數(shù)矩陣。
圖3 插值模版示意圖Fig.3 Interpolation template schematic diagram
為了選取比較光滑的模版值,在選取模版時,每個物理量減去他們的平均值,最后再加上被減去的值。將各個點的流場信息做差后的ΔZ帶入到RBF插值函數(shù)之中,得到的線性方程組如下:
ΦW=ΔZ
(10)
(11)
W=[ωi,ωj,ωk,ωl]T
(12)
式中:‖ri-rj‖為歐式距離,即每個基網(wǎng)格格心之間的距離,求解方程后得到每個基網(wǎng)格的權(quán)重系數(shù)ωi,將需要插值的節(jié)點P坐標(biāo)rP帶入到徑向基函數(shù)中就可以重構(gòu)出該點的值。
(2) Wendland’s C2型徑向基函數(shù)
基本形式如下所示:
(13)
插值半徑的選取對于徑向基函數(shù)插值方法的效率有著很重大的影響。較大的插值半徑會擴大物面插值的影響域,與此同時,也會使得基矩陣的維度增加并引起計算量的增加。選擇較小的插值半徑可以減小矩陣維度,降低計算量。
徑向基函數(shù)插值精度高,靈活性很強,使用不同類型的插值函數(shù)的插值效果也有所差異。對于同一個徑向基函數(shù),選取不同的參數(shù)d也會得到不同的插值效果。本文所使用的C2徑向基函數(shù),作用半徑d的初始值為插值模版中格心和P點(嵌套區(qū)域需要進(jìn)行插值的點)的最大距離的50倍。矩陣的條件數(shù)不僅是判斷矩陣病態(tài)與否的一種度量,也可以表示矩陣計算對誤差的敏感度強弱,條件數(shù)越大,數(shù)值穩(wěn)定性越差。對于C2徑向基函數(shù)的系數(shù)矩陣Φ,取二范數(shù)后的條件數(shù)可以表示為:
(14)
為了改善在一些網(wǎng)格上的條件數(shù)過大的情況,本文采用了一種基于控制條件數(shù)的徑向基函數(shù)插值算法,通過人為的改變C2徑向基函數(shù)的作用半徑d來調(diào)整插值矩陣的條件數(shù),從而在條件數(shù)過大的插值模版上進(jìn)行有效的控制,有效改善了由于條件數(shù)過大使得穩(wěn)定性較差以及出現(xiàn)數(shù)值奇性的現(xiàn)象。具體實施的步驟為:當(dāng)流場求解開始時,人為輸入一個期望達(dá)到的條件數(shù),當(dāng)使用徑向基函數(shù)插值進(jìn)行流場信息傳遞時,每次迭代對貢獻(xiàn)單元插值矩陣的插值半徑自動乘以0.9,計算所有貢獻(xiàn)單元插值矩陣的條件數(shù),不滿足要求則返回上一步繼續(xù)縮小插值矩陣半徑,直到所有貢獻(xiàn)單元插值矩陣的條件數(shù)都滿足要求。
為了分析本文所使用的插值函數(shù)進(jìn)行解析的精度,使用解析函數(shù)作為插值精度驗證算例。該函數(shù)的數(shù)學(xué)形式為:
(15)
嵌套網(wǎng)格組裝過程示意圖如圖4所示,圓心在(0,0)處,半徑為1。
(a) 嵌套組裝前網(wǎng)格
(b) 嵌套完成后網(wǎng)格圖4 嵌套網(wǎng)格組裝過程示意圖Fig.4 Diagram of overset grid assembly process
使用線性插值,以距離倒數(shù)作為權(quán)重插值,逆二次徑向基函數(shù)插值、Wendland’s C2徑向基函數(shù)插值四種不同的方法進(jìn)行插值,四種插值方法的誤差統(tǒng)計如表1所示。
表1 四種插值方式在嵌套區(qū)域的誤差統(tǒng)計Table 1 Error statistics of four interpolation modes in nest area
從表1可以看出:使用C2徑向基函數(shù)進(jìn)行插值具有較小的誤差,得到的結(jié)果和解析解最為接近。在整個插值的區(qū)域內(nèi),C2徑向基函數(shù)插值的平均絕對誤差僅為0.007 254 9%,而IQ徑向基函數(shù)插值為0.036 369%,距離反比函數(shù)插值(WA)為0.197 52%,線性插值(LINE)為2.361 2%,說明使用C2徑向基函數(shù)進(jìn)行插值具有更高的計算精度,本文將優(yōu)先選擇使用C2徑向基函數(shù)插值進(jìn)行流場信息傳遞。
當(dāng)直接使用C2徑向基函數(shù)用于嵌套網(wǎng)格插值過程中時,在某些點可能存在數(shù)值原因引起的振蕩或插值誤差過大的點。以MD30P/30N多段翼為例,計算馬赫數(shù)Ma=0.2,雷諾數(shù)Re=9.0×106,迎角α=8.010 9°,湍流模型選擇Spalart-Allmaras(SA)。圍繞前緣縫翼、主翼、后緣襟翼分別生成計算網(wǎng)格,其中前緣縫翼網(wǎng)格量約為6萬,后緣襟翼網(wǎng)格量約為11萬,主翼網(wǎng)格量約為22萬。三套網(wǎng)格按照之前距離判定的準(zhǔn)則自動進(jìn)行網(wǎng)格嵌套組裝。三套網(wǎng)格均按照預(yù)期完成組裝示意圖如圖5所示。
圖5 MD30P/30N嵌套網(wǎng)格組裝示意圖Fig.5 Assembly schematic of MD30P/30N overset grid
未采用控制條件數(shù)與采用控制條件數(shù)后的C2徑向基函數(shù)插值得到的流場圖如圖6所示。從流場計算結(jié)果來看,直接使用C2函數(shù)進(jìn)行插值得到了圖6(a)所示的結(jié)果,可以看出:在嵌套區(qū)域存在有數(shù)值奇性。根據(jù)之前提出的條件數(shù)控制策略,對于不同的插值模板,通過調(diào)整他們的插值半徑來控制條件數(shù)。經(jīng)過條件數(shù)控制為10,即cond=10時,得到了圖6(b)所示的計算結(jié)果,整個流場的求解更加光滑,消除了之前出現(xiàn)的數(shù)值奇性現(xiàn)象。
(a) 未采用控制條件數(shù)
(b) 采用控制條件數(shù)圖6 MD30P/30N流場計算壓力分布云圖Fig.6 Pressure distribution cloud diagram of MD30P/30N flow field
MD30P/30N表面Cp分布與實驗對比示意圖如圖7所示,可以看出,當(dāng)直接使用C2徑向基函數(shù)插值進(jìn)行流場信息傳遞時,表面Cp分布在主翼上表面與后緣襟翼下表面與實驗數(shù)據(jù)差別較大,當(dāng)采取條件數(shù)控制策略的C2徑向基函數(shù)插值策略將矩陣條件數(shù)控制在10后,所得到的翼型表面Cp分布與實驗數(shù)據(jù)吻合度較高,說明該插值方法能夠很好地解決直接使用C2徑向基函數(shù)插值帶來的數(shù)值奇性問題,且計算精度較高。MD30P/30N殘差收斂歷程對比示意圖如圖8所示,當(dāng)條件數(shù)控制在10時,殘差下降的速率較快,說明采取條件數(shù)控制策略的C2徑向基函數(shù)插值策略能有效改善殘差的收斂效果。
圖7 MD30P/30N表面Cp分布示意圖Fig.7 Schematic diagram of MD30P/30N surface Cp distribution
圖8 MD30P/30N殘差收斂歷程對比示意圖Fig.8 Comparison diagram of MD30P/30N residual convergence process
AEDC標(biāo)準(zhǔn)外掛物分離模型的計算參數(shù)[18-19]如圖9所示。計算馬赫數(shù)Ma=0.95,雷諾數(shù)Re=7.874×106,迎角α=0°,湍流模型選擇Spalart-Allmaras(SA)。圍繞機翼和外掛物分別生成計算網(wǎng)格,網(wǎng)格總數(shù)約450萬,兩套網(wǎng)格按照之前的距離判定的準(zhǔn)則自動進(jìn)行網(wǎng)格嵌套組裝。
圖9 標(biāo)準(zhǔn)外掛物分離模型的計算參數(shù)Fig.9 The calculation parameters of standard separation model
為了考察使用控制條件數(shù)策略的C2徑向基函數(shù)的插值效果,將插值矩陣的條件數(shù)控制在100,即cond=100時,如圖10所示。從流場計算結(jié)果來看,可以看出:未采用控制條件數(shù)時,在嵌套區(qū)域的邊界上,出現(xiàn)了數(shù)值奇性,這些異常導(dǎo)致在尾跡區(qū)出現(xiàn)了和周圍顯著不同的區(qū)域;當(dāng)采用控制條件數(shù)后,尾跡區(qū)流場過渡的更加光滑和連續(xù),這說明采用條件數(shù)控制策略的插值方法可以有效消除嵌套邊界區(qū)域出現(xiàn)的數(shù)值奇性的問題。
(a) 無條件數(shù)控制的z-x截面
(b) 采用條件數(shù)控制的z-x截面
(c) 無條件數(shù)控制的x-y截面
(d) 采用條件數(shù)控制的x-y截面圖10 外掛物截面流場對比示意圖Fig.10 Flow field comparison diagram of the cross section
外掛物表面Cp分布與實驗值對比示意圖如圖11所示,可以看出:直接使用C2徑向基函數(shù)和未使用控制條件數(shù)策略的C2徑向基函數(shù)進(jìn)行插值得到的外掛物周向Cp分布均與實驗值較為吻合,當(dāng)加入控制插值矩陣條件數(shù)后,外掛物在彈尾處的Cp分布與實驗值吻合度更高一些。
(a) 外掛物兩側(cè)表面Cp分布
(b) 外掛物上下表面Cp分布圖11 外掛物表面Cp分布與實驗對比示意圖Fig.11 Schematic diagram of Cp distribution on the surface of the plug-in and comparison with the experiment
AEDC外掛物殘差收斂歷程對比曲線如圖12所示,可以看出:未采用控制條件數(shù)時,殘差收斂的歷程波動較為劇烈;采用控制插值矩陣條件數(shù)后,殘差收斂效果較好,呈現(xiàn)不斷下降的趨勢,說明采用條件數(shù)控制策略的C2徑向基函數(shù)插值策略能提高殘差收斂的效果。
圖12 殘差收斂歷程對比Fig.12 Comparison of residual convergence process
非定常數(shù)值模擬得到的外掛物質(zhì)心處的位移Xg、Yg、Zg與實驗值的對比結(jié)果如圖13所示,可以看出:采用控制條件數(shù)策略的C2徑向基函數(shù)插方法數(shù)值模擬得到的水平位移Xg整體與實驗吻合度更高一些。兩種不同的插值方式得到外掛物沿豎直方向的位移與實驗結(jié)果吻合較好。
圖13 質(zhì)心位移與實驗值的對比圖Fig.13 Comparison diagram of center of mass displacement and experimental value
外掛物分離過程中的姿態(tài)角變化曲線如圖14所示。
圖14 姿態(tài)角與實驗值的對比圖Fig.14 Comparison diagram of attitude angle and experimental value
從圖14可以看出:使用控制條件數(shù)策略的C2徑向基函數(shù)進(jìn)行插值得到的俯仰角θ變化曲線在分離中后期與實驗吻合更好一些,而偏航角ψ與滾轉(zhuǎn)角φ使用兩種插值方法數(shù)值模擬得到的結(jié)果較為接近,均與實驗結(jié)果吻合度較高。同時可以看出彈射力矩與氣動力的共同作用,使得外掛物的俯仰角θ會呈現(xiàn)先增加后減小的趨勢;外掛物受到三角翼下洗與側(cè)洗的影響,偏航角也會逐漸增加。
(1) 使用解析函數(shù)作為插值模板,對線性插值、距離倒數(shù)作為權(quán)重插值,逆二次徑向基函數(shù)插值、Wendland’s C2徑向基函數(shù)插值四種不同的嵌套網(wǎng)格插值方法的精度進(jìn)行驗證。計算結(jié)果表明C2徑向基函數(shù)的插值精度較高,同時收斂較快。
(2) 改進(jìn)后的C2徑向基函數(shù)插值方法能夠有效消除嵌套插值區(qū)域數(shù)值奇性,同時計算收斂的速率更快,魯棒性較好。