避免核函數(shù)導數(shù)的二階導數(shù)核近似方法的修正

·， ·江

(新疆大學 數(shù)學與系統(tǒng)科學學院,烏魯木齊 830046)

1 引 言

SPH方法是Monaghan等[1,2]提出的發(fā)展比較快、應用研究最為廣泛的無網(wǎng)格方法之一，最初用于解決天體物理學和宇宙學中模擬三維無界空間天體的演化問題。現(xiàn)在廣泛地應用于計算固體力學[3]、流體力學[4]和多相流[5]等各個領(lǐng)域[6]。用經(jīng)典的SPH方法進行計算模擬時，在邊界處和不連續(xù)處有較大誤差，其主要是由于此處的粒子缺失和粒子不均勻分布等因素造成的，于是便出現(xiàn)很多不同類型的SPH 核近似方法的修正[7-12]，以彌補傳統(tǒng)SPH核近似方法的一些缺點，導致了SPH核近似方法表達公式的多樣性[13,14]。

本文對已提出的避免核函數(shù)導數(shù)的二階導數(shù)核近似方法SODF-SPH方法[15]進行了修正。首先，簡單介紹SPH核近似方法的基本原理和計算函數(shù)導數(shù)的SODF-SPH方法；其次，在一元函數(shù)和多元函數(shù)的情形下給出計算函數(shù)二階導數(shù)SODF-SPH方法的修正公式；最后，通過函數(shù)二階導數(shù)的數(shù)值計算和熱傳導問題的模擬進行誤差對比分析，從而驗證該方法對提高計算精度和加快收斂速度的有效性。

2 傳統(tǒng)SPH 方法

傳統(tǒng)SPH方法的基本原理見式(1)，在區(qū)域Ω上定義的任何函數(shù)f(x)表示為

(1)

式中δ(x-x′)為狄拉克函數(shù)，定義為

式(1)的狄拉克函數(shù)用光滑核函數(shù)W(x-x′,h)=(1/hd)W[(x-x′)/h]來代替，得到函數(shù)f(x)的核近似表達式

(2)

式中d為函數(shù)自變量的個數(shù)，光滑核函數(shù)依賴于兩點之間的距離|x-x′|和定義核函數(shù)影響區(qū)域ΩW的光滑長度h,核函數(shù)影響區(qū)域在一維、二維和三維情形下分別為以x為中心的區(qū)間、圓和球，半徑為κh，不同的核函數(shù)中κ取不同的值[16]。

考慮函數(shù)乘積的導數(shù)、格林公式及核函數(shù)的性質(zhì)，n階偏導數(shù)的核近似可表示為

(3)

(4)

很顯然M(k)n為常量，并不依賴于h。

將式(2)的f(x)在x處的泰勒級數(shù)展開式替換后得

(5)

由核函數(shù)的對稱性和歸一性得M0=1和M1=0，因此

〈f(x)〉=f(x)+o(h2)

(6)

同樣的方法可以得到

(7)

〈f′(x)〉=f′(x)+o(h2)

(8)

以此類推，

〈f″(x)〉=f″(x)+o(h)

(9)

很顯然核函數(shù)影響范圍在求解區(qū)域內(nèi)時，核近似式(2)有o(h2)精度。式(2)中，當n=1時，有o(h2)精度；當n=2時，有o(h)精度。核函數(shù)影響范圍不在求解區(qū)域內(nèi)或接近于邊界時，該SPH 核近似方法精度會降低，直接影響計算精度和穩(wěn)定性。

(10)

(11)

粒子和近似式的精度不僅與核函數(shù)光滑長度h有關(guān)，而且與微元體的體積ΔV有關(guān)[16]。

3 函數(shù)導數(shù)計算的SODF-SPH方法

3.1 一元函數(shù)導數(shù)計算的SODF-SPH方法

由第2節(jié)討論的函數(shù)核近似表達的基本思路可得

(12)

式中 如果核函數(shù)的影響區(qū)域未受求解區(qū)域的邊界截斷，由核函數(shù)對稱性質(zhì)有M3=0，此時解出f″(x)為

(13)

很顯然，當核函數(shù)的影響區(qū)域沒有受到求解區(qū)域的邊界截斷時，用式(13)得到的f″(x)的核近似式具有o(h)精度；否則不能保證o(h)精度。

3.2 多元函數(shù)偏導數(shù)計算的SODF-SPH方法

多元函數(shù)偏導數(shù)計算的SODF-SPH 方法與一元函數(shù)一樣的思路相似，為表述方便，先引入與核函數(shù)光滑長度h無關(guān)的常量

U=x-x′，則

(14)

式中i和j為與函數(shù)自變量個數(shù)有關(guān)的啞標，l,k=1,2,…,d。

與一元函數(shù)類似條件下，有Milk3=0，此時由式(15)解出

Mlk2f(x)]

(15)

利用第2節(jié)討論的SPH方法粒子和近似表示的基本思路，可以得出一元函數(shù)導數(shù)和多元函數(shù)偏導數(shù)計算的SODF-SPH方法核近似公式 (12，14)在粒子i處的粒子和近似表達式。

4 函數(shù)導數(shù)計算的SODF-SPH方法的修正

4.1 一元函數(shù)二階導數(shù)計算的SODF-SPH方法的修正

(16)

由式(16)將f′(x)表示為

(17)

將式(17)代入式(12)整理得

(18)

式中M3=0時為一元函數(shù)導數(shù)計算的SODF-SPH公式(13)，由式(18)的推導過程可知，不論M3=0還是M3≠0，式(18)得到的f″(x)的核近似式均具有o(h)精度。

4.2 多元函數(shù)二階偏導數(shù)計算的SODF-SPH方法的修正

對多元函數(shù)同樣有

(19)

式中

(20)

式(20)代入式(14)，整理得

Bl kf(x)]+o(h)

(21)

同樣可以得出一元函數(shù)導數(shù)和多元函數(shù)偏導數(shù)計算的SODF-SPH方法核近似修正公式(18,21)在粒子i處的粒子和近似表達式。

5 數(shù)值驗證和誤差對比分析

定義在數(shù)值驗證和誤差對比分析過程中將要用到的誤差的無窮范數(shù)、二范數(shù)和絕對誤差為

(22)

(23)

ei=|fExact(xi)-fComputed(xi)|

(24)

其中M是求解域內(nèi)的總粒子數(shù)，fExact(xi)和fComputed(xi)分別是xi處的精確值和SPH計算值。

本文所有數(shù)值算例，如果沒有特別說明，計算過程中核函數(shù)的光滑長度等于粒子間距的1.2倍，使用的核函數(shù)是三次B -樣條函數(shù)[16]。

5.1 函數(shù)二階導數(shù)的計算

(1) 一元函數(shù)算例?？紤]函數(shù)f1(x)=e- x2，f2(x)=sin(πx)，f3(x)=x3，x∈[-1,1]。

在區(qū)間[-1,1]上均勻分布不同個數(shù)的粒子時，用SODF-SPH方法和修正方法計算函數(shù)二階導數(shù)所產(chǎn)生誤差的二范數(shù)的比較列入表1。由表1可知，隨著粒子數(shù)的增大，用SODF-SPH方法計算函數(shù)二階導數(shù)的誤差越來越大，其主要原因是隨著粒子數(shù)的增加，在邊界附近誤差劇增。然而隨著粒子數(shù)的增大，本文提出的修正方法計算函數(shù)二階導數(shù)的誤差越來越小，修正方法的誤差遠小于SODF-SPH方法。

數(shù)控專業(yè)是一門實踐性很強的專業(yè)，其中涉及到很多抽象的理論知識，對于剛剛接觸這些知識的學生而言，在學習時極易遇到一些困難。所謂數(shù)字化課程，其是將網(wǎng)絡(luò)作為支持環(huán)境，運用數(shù)字化學習材料，使用數(shù)字化信息傳遞與人際交互方式的一種教學方式。在科學技術(shù)發(fā)展之下，數(shù)字化課程成為當前教育行業(yè)中一種新型輔助性教學方式，可以有效提升教學質(zhì)量與效率。

當Δx=0.05，并在區(qū)間[-1,1]上均勻分布41個粒子，對不同的函數(shù)進行計算二階導數(shù)時，用SODF-SPH方法和修正方法的誤差二范數(shù)隨著核函數(shù)光滑長度h的變化比較如圖1所示，可以看出，本文提出的修正方法選擇不同的光滑長度時，誤差的變化不太大，且比SODF-SPH方法的誤差小，即本文修正方法對提高精度起了很大的作用。

(2) 多元函數(shù)算例。考慮函數(shù)f1(x)=e- x2- y2，(x，y)∈[-1,1]×[-1,1]。

表1 一元函數(shù)二階導數(shù)誤差二范數(shù)的比較Tab.1 Comparison of error 2-norms in the second order derivative

在區(qū)域[-1,1]×[-1,1]內(nèi)均勻分布1681個粒子，且Δx=0.05，Δy=0.05時，用SODF-SPH方法和修正方法計算函數(shù)對x的二階偏導數(shù)，圖2給出角落區(qū)域[0.85,1]×[0.85,1]內(nèi)兩種方法計算偏導數(shù)的絕對誤差對比?？梢钥闯觯拚椒ㄔ谶吔绺浇恼`差小于SODF-SPH方法，對其他類型的函數(shù)f2(x)=sin(πxy)和f3(x)=(xy)3，(x，y)∈[-1,1]×[-1,1]進行計算時，結(jié)論一致。

5.2 熱傳導問題的模擬

在直角坐標系下，瞬態(tài)熱傳導問題的控制方程為

圖1 一元函數(shù)二階導數(shù)誤差比較

Fig.1 Comparation of error in the second order on one variete case

圖2 多元函數(shù)二階偏導數(shù)絕對誤差比較

Fig.2 Comparation of absolute error in the second order derivative on multivariate case

(25)

初始條件為T(x,0)=T0(x)，x∈Ω

Dirchlet邊界條件和Neumann邊界條件分別為

T(x,t)=T1(t) (x∈?ΓD)

T(x,t)·n+bT(x,t)=0 (x∈?ΓN)

式中n為區(qū)域Ω邊界外單位向量，bT(x,t)為邊界熱源，?ΓD∩?ΓN=Φ， ?ΓD∪?ΓN=?Ω， ?Ω為區(qū)域邊界。

(1) 一維熱傳導問題的模擬

初始和邊界條件分別為T(x,0)=T0，T(0,t)=T0，T(L,t)=T1，Q(x,t)=0。x∈[0,L]時，該問題有解析解[17]。

本文的模擬取L=1，T0=0 ℃，T1=1 ℃，并分別用SODF-SPH 方法和修正方法進行模擬，圖3 給出設(shè)立了均勻分布的21個粒子，Δt=1×10-3s時,在不同時刻的數(shù)值模擬解與解析解的對比。結(jié)果表明，用本文提出的修正方法可以得到更精確的數(shù)值解。

(2) 二維熱傳導問題的模擬

二維熱傳導問題中邊界和初始條件分別為T(0,y,t)=0，T(π,y,t)=π2，T(x,0,t)=x2，T(x,π,t)=x2，T(x,y,0)=x2+sinxsiny，熱源強度Q(x,t)=-2，此時該問題對應的解析解為

T(x,y,t)=x2+e-2tsinxsiny

(26)

求解區(qū)域內(nèi)均勻分布20×20粒子，時間步長Δt=1×10-3s時進行計算，圖4為沿x和y軸方向在不同時刻的修正方法數(shù)值解和解析解?？梢钥闯?，修正方法數(shù)值解與解析解吻合良好。表2列出了修正方法與SODF-SPH方法的誤差對比?？梢钥闯?修正方法的誤差比較小，提高了數(shù)值解的精度。

圖3 不同時刻溫度曲線數(shù)值解比較

Fig.3 Compare of temperature curve in different time steps

圖4 溫度隨時間的變化分布

Fig.4 Temperature distribution in different time steps

表2 修正方法與SODF-SPH方法誤差范數(shù)的比較

Tab.2 Comparison of error norms between the corrective method and the SODF-SPH method

t/sSODF-SPH方法修正方法L∞L2L∞L20.010.02189128240.00455192480.00011614460.00005520450.10.05614959190.01408896600.00097016550.00045266920.50.07057004720.02115713860.00215224800.00101575080.90.07510599540.02352524300.00173207510.0008219932

6 結(jié) 論

本文基于泰勒級數(shù)展開對已提出的計算一元函數(shù)函數(shù)二階導數(shù)和多元函數(shù)的二階偏導數(shù)的SODF-SPH方法進行了修正。將修正公式應用到一元函數(shù)的二階導數(shù)和多元函數(shù)的二階偏導數(shù)的計算和各種熱傳導問題的模擬，得出如下結(jié)論。

(1) 在粒子間距、核函數(shù)和光滑長度分別相同的情況下，對函數(shù)的二階導數(shù)進行計算，本文提出的修正方法在提高SODF-SPH方法的計算精度和收斂速度方面起了很大的作用，特別在邊界附近能有效地減小誤差。

(2) 本文修正方法的計算量比原方法大，但任何點均具有o(h)精度。