·， ·江

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

1 引 言

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

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

2 傳統(tǒng)SPH 方法

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

(1)

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

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

(2)

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

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

(3)

(4)

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

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

(5)

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

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

(6)

同樣的方法可以得到

(7)

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

(8)

以此類(lèi)推，

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

(9)

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

(10)

(11)

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

3 函數(shù)導(dǎo)數(shù)計(jì)算的SODF-SPH方法

3.1 一元函數(shù)導(dǎo)數(shù)計(jì)算的SODF-SPH方法

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

(12)

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

(13)

很顯然，當(dāng)核函數(shù)的影響區(qū)域沒(méi)有受到求解區(qū)域的邊界截?cái)鄷r(shí)，用式(13)得到的f″(x)的核近似式具有o(h)精度；否則不能保證o(h)精度。

3.2 多元函數(shù)偏導(dǎo)數(shù)計(jì)算的SODF-SPH方法

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

U=x-x′，則

(14)

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

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

Mlk2f(x)]

(15)

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

4 函數(shù)導(dǎo)數(shù)計(jì)算的SODF-SPH方法的修正

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

(16)

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

(17)

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

(18)

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

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

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

(19)

式中

(20)

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

Bl kf(x)]+o(h)

(21)

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

5 數(shù)值驗(yàn)證和誤差對(duì)比分析

定義在數(shù)值驗(yàn)證和誤差對(duì)比分析過(guò)程中將要用到的誤差的無(wú)窮范數(shù)、二范數(shù)和絕對(duì)誤差為

(22)

(23)

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

(24)

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

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

5.1 函數(shù)二階導(dǎo)數(shù)的計(jì)算

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

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

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

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

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

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

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

5.2 熱傳導(dǎo)問(wèn)題的模擬

在直角坐標(biāo)系下，瞬態(tài)熱傳導(dǎo)問(wèn)題的控制方程為

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

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

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

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

(25)

初始條件為T(mén)(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) 一維熱傳導(dǎo)問(wèn)題的模擬

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

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

(2) 二維熱傳導(dǎo)問(wèn)題的模擬

二維熱傳導(dǎo)問(wèn)題中邊界和初始條件分別為T(mén)(0,y,t)=0，T(π,y,t)=π2，T(x,0,t)=x2，T(x,π,t)=x2，T(x,y,0)=x2+sinxsiny，熱源強(qiáng)度Q(x,t)=-2，此時(shí)該問(wèn)題對(duì)應(yīng)的解析解為

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

(26)

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

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

Fig.3 Compare of temperature curve in different time steps

圖4 溫度隨時(shí)間的變化分布

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é) 論

本文基于泰勒級(jí)數(shù)展開(kāi)對(duì)已提出的計(jì)算一元函數(shù)函數(shù)二階導(dǎo)數(shù)和多元函數(shù)的二階偏導(dǎo)數(shù)的SODF-SPH方法進(jìn)行了修正。將修正公式應(yīng)用到一元函數(shù)的二階導(dǎo)數(shù)和多元函數(shù)的二階偏導(dǎo)數(shù)的計(jì)算和各種熱傳導(dǎo)問(wèn)題的模擬，得出如下結(jié)論。

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

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