局部近似特別解的無網(wǎng)格方法

在傳統(tǒng)的基于網(wǎng)格的微分方程數(shù)值方法（如有限差分法（fdm）和有限元法（fem）中，傳統(tǒng)的基于網(wǎng)格的數(shù)值法廣泛應(yīng)用于計算水流和固體力學(xué)領(lǐng)域。但這些方法仍然存在諸多不足之處,因?yàn)槠湓诤艽蟪潭壬先Q于網(wǎng)格劃分的質(zhì)量,當(dāng)劃分的網(wǎng)格出現(xiàn)畸變時會嚴(yán)重影響數(shù)值解的精度甚至導(dǎo)致計算失敗,從而使得這些方法在許多問題上的應(yīng)用受到了限制。于是無網(wǎng)格方法應(yīng)運(yùn)而生,并且在求解偏微分方程方面得到很大的關(guān)注。無網(wǎng)格方法包括:Lucy和Gingold提出的光滑粒子流體動力學(xué)方法(SmoothedParticleHydrodynamics,SPH),基本解方法(MFS),局部近似特別解方法(LMAPS)等。在眾多的無網(wǎng)格方法中基于徑向基函數(shù)(RadialBasisFunctions,RBFs)的無網(wǎng)格方法備受矚目。RBFs是一類以場點(diǎn)和源點(diǎn)之間的歐式距離為自變量的函數(shù),它具有形式簡單,各項(xiàng)同性等優(yōu)點(diǎn),已被成功應(yīng)用于多變量插值。但基于全局性質(zhì)的RBFs的無網(wǎng)格方法產(chǎn)生的系數(shù)矩陣常常是滿秩,有時甚至是病態(tài)的,適于解決有限插值點(diǎn)的問題,但應(yīng)用到解決大規(guī)模散亂點(diǎn)的數(shù)值計算問題時則受到限制。為了克服上述問題,文中采用局部性質(zhì)的近似特別解方法解決大規(guī)模問題,研究表明效果顯著。該方法只需創(chuàng)建局部區(qū)域并利用該區(qū)域的點(diǎn)構(gòu)造局部低階矩陣,然后將局部形式推廣至全局形式形成大規(guī)模稀疏矩陣,最后通過求解這個稀疏線性方程組就可以得到偏微分方程的近似解,從而可以成功避免上述問題。1局部化的map方法首先考慮如下的橢圓形偏微分方程:其中,Δ為二階線性Laplace算子;a(x),b(x),c(x),f(x)和g(x)為給定的函數(shù);B為邊界算子。MAPS的主要思想是將式(1)轉(zhuǎn)化為Poisson型方程,即其中,H(x,u,ux,uy)=-a(x)ux-b(x)uy-c(x)·u+f(x),并且H可以用近似表示運(yùn)用中常取MQ函數(shù)作為基函數(shù),其表達(dá)形式為其中,c為形參,決定基函數(shù)的形狀。當(dāng)c很小時,MQ函數(shù)成錐形形狀;當(dāng)c逐漸變大時,函數(shù)圖像則越來越平滑。由此可以看出,形參c對于近似式的精確性至關(guān)重要。影響c的取值因素主要有3個,即支撐域內(nèi)節(jié)點(diǎn)的數(shù)量、支撐域內(nèi)的節(jié)點(diǎn)分布以及支撐域的大小。一般情況下,c越大,近似誤差就越小。應(yīng)用MQ函數(shù)插值近似會出現(xiàn)這樣的情況:為提高近似精度增大c值時得到的系數(shù)矩陣往往會變得越來越病態(tài),因而舍入誤差會變大,增加了解的不穩(wěn)定性。此外,如何適當(dāng)選取c值并不是一件容易的事情,這也成為學(xué)界研究討論的熱點(diǎn)問題。根據(jù)式(5)和式(6)推得由式(3)、式(4)和式(7)可知式中:重組式(8),有其中并且式(2)可寫為下面具體介紹局部化的MAPS方法。由于Ωi={xki}kni=1中的每一個點(diǎn)對式(11)都成立,于是有在式(12)中右端的方陣可以用Φni表示,并且Φni是非奇異的,因此是可逆陣,所以得到的式(12)中未知系數(shù)為其中因此,在式(11)中的^u(xi)可以表示為ni個節(jié)點(diǎn)函數(shù)值的線性組合。其中類似于式(15)有對于邊界條件,在局部區(qū)域上,有于是在全局區(qū)域上有所以由式(17)和式(19)可以構(gòu)造一個N維稀疏線性方程組,通過求解該方程組可以得到求解該區(qū)域上所有結(jié)點(diǎn)上的近似解。2局部近似方法利用LMAPS進(jìn)行數(shù)值計算時,找到區(qū)域中每個計算結(jié)點(diǎn)相鄰近的ni個結(jié)點(diǎn)是很重要的,當(dāng)插值點(diǎn)為很龐大的數(shù)值時,有效的搜索算法非常重要。在這些算法中,kd-tree算法是一種k維數(shù)據(jù)空間的有效搜索算法。為了驗(yàn)證文中數(shù)值精度,分別定義均方根誤差(RMSE),最大絕對誤差(MAE),最大相對誤差(RAE)如下:式中:N為測試點(diǎn)的數(shù)目;^uk為近似解;uk為精確解。算例1考慮Dirichlet邊界型條件的Poisson方程:其中,Ω∪ue785Ω為如圖1所示的一個不規(guī)則區(qū)域。由圖1可以看到,內(nèi)部和邊界結(jié)點(diǎn)都是均勻分布的,在計算時取14408個內(nèi)部結(jié)點(diǎn)和300個邊界結(jié)點(diǎn)。該二維Poisson方程的精確解為u(x,y)=e-2x+3y。圖2為分別利用LMAPS方法和傳統(tǒng)的LDQ方法得到的絕對誤差。由圖2(a)和圖2(b)比較可以看出,LMAPS方法得到的誤差比LDQ方法更小,因此,LMAPS方法略優(yōu)于LDQ方法。同時還可以看出,誤差在不規(guī)則區(qū)域內(nèi)的大部分情況都是下降平穩(wěn)的,除了在邊界周圍的一部分點(diǎn)有一定的跳躍。但作為局部方法,它們與全局的徑向基函數(shù)方法比較,計算效率還是比較顯著的。算例2采用局部近似特別解方法求解方腔驅(qū)動流問題,以檢驗(yàn)這種方法的精度和準(zhǔn)確性。不可壓流體是由如下原始變量形式的Navier-Stokes方程控制。動量方程:連續(xù)方程:圖3為Re=50和Re=500時的方腔驅(qū)動流的流線圖。當(dāng)Re很小時,靠近腔體上壁有一個呈順時針方向的漩渦;隨著Re的增大,方腔內(nèi)部受到流體流動的影響也會變大,從而使得靠近壁面的速度梯度越來越大且腔內(nèi)主漩渦逐漸向腔體右下方移動。當(dāng)Re=50腔體右下角和左下角先后次生出一個小渦;隨著Re繼續(xù)增大,這兩個次生的小渦會越來越大,而主漩渦則逐漸向腔內(nèi)中心位置移動。圖4為不同雷諾數(shù)下方腔垂直中心線上速度u和水平中心線上速度v的分布情況。根據(jù)文獻(xiàn)將得到的數(shù)值解與Ghia和Sanyasiraju的結(jié)果進(jìn)行比較,得到的解和參考解具有相似的精度。如果增加區(qū)域內(nèi)結(jié)點(diǎn)分布的密度,局部近似特別解方法的精度還會有一定的提高。由圖3中可以看出,渦流開始在方腔中心軸線的上半部分且處于一個穩(wěn)定的狀態(tài),隨著雷諾數(shù)的增大渦流逐漸向方腔的中心移動。由圖4可以看出,雷諾數(shù)從50到500,在薄邊界層的附近,u,v的速度梯度在不斷地增加。3大力采用