SSPH核近似的降元算法研究?

熱合買提江·依明，買買提明·艾尼

(1.新疆大學機械工程學院，新疆烏魯木齊 830049；2.新疆大學 數學與系統科學學院，新疆 烏魯木齊830046)

用基于網格的（有限差分，有限元、有限體積法等）方法進行數值計算時，在求解區域內的一些離散點進行離散化，每個離散點與相鄰離散點用網格來連接起來，這樣計算過程中避免不了網格重構，這樣不僅計算費用昂貴，而且會使計算精度受損.近年來發展較快的無網格方法之一：光滑粒子流體動力學(Smoothed Particle Hydrodynamics，簡稱SPH)方法，將求解區域離散成一系列粒子，每一個粒子攜帶相應的物理變量（速度、壓力、溫度、應力和應變等）并時時跟蹤求解域中物體的運動，可以完全拋開網格重構、容易編程實現而且保證良好的計算精度.

SPH方法是一種新的無網格純Lagrange粒子方法.1977年，由Lucy[1]和Monaghan[2]等人提出，最初用于解決天體物理學和宇宙學中模擬三維無邊界空間天體的演化問題.Gingold和Monaghan等人相繼提出了適用于SPH方法的人為粘性和人為熱流項后[3,4]，又延伸到了計算流體力學的各個領域.Libersky和Petschek[5]在SPH中引入了材料強度模型，模擬了動態斷裂和破片運動等固體的動態響應，此后SPH被迅速的應用到固體力學的研究當中.到今天，SPH方法已廣泛地應用于計算科學和模擬的各個領域[6].

隨著SPH方法應用范圍的擴展和對其研究的深入，它的一些缺點被人們發現，例如計算精度比較低，在邊界缺乏插值一致性，穩定性不高等等，于是便出現對SPH核近似方法的很多不同類型的修正，彌補傳統SPH核近似方法的一些缺點，這些修正和彌補導致了SPH核近似方法公式的多樣性.如重構（Reproducing）SPH核近似方法[7]，在SPH方法中引入修正項，滿足再生性條件，提高算法的精度和穩定性；正交化（Normalized）SPH核近似方法[8]推導出了新的離散公式，對于邊界附近粒子分布不對稱情況和同物質的界面附近在一定程度上提高了計算精度；改進(Corrective)核近似方法[9]不但能解決拉伸不穩定性問題，還能夠提高求解區域內及邊界區域處的計算精度；基于泰勒級數展開構造的修正（Modified）SPH核近似方法[10～12]和對稱（Symmetric）SPH核近似方法[13]也能夠有效的處理張力不穩定和邊界處缺陷問題等問題.

本文在研究SSPH核近似方法的基礎上，提出了SSPH核近似的降元算法并對其精度進行分析，主要包括以下幾個方面的內容.首先介紹了SPH方法的基本原理，其次給出了SSPH核近似方法和降元算法推導過程，然后對計算精度進行了比較，最后提出的算法用于非穩態熱傳導問題，進行了數值計算和模擬，驗證了算法的實際應用的有效性.

1 傳統SPH核近似方法

傳統SPH核近似方法基本原理比較簡單，將任意函數f(x)、一階導數和二階導數在某一點x處的SPH核近似式，用函數f(x)與核函數乘積的積分來表示.這里直接給出其表達式，函數f(x)的SPH核近似表達式為

一階導數?f(x)的可表示為

二階導數?2f(x)的可表示為

當粒子對稱分布時，SPH方法中函數、函數一階導數和函數二階導數的SPH核近似均有O(h2)階精度[14]，但當粒子不對稱分布或接近于邊界附近時，該SPH核近似方法精度會降低，直接會影響計算精度和穩定性.

在SPH方法中把整個求解域離散成一系列任意分布的有質量和容量的有限多個粒子如圖1所示，如果核近似式中粒子j處微元體的體積為?Vj、質量為mj、密度為ρj、位置矢量為xj，記函數f(x)在粒子j處的值為假設粒子i為中心的核函數影響范圍內的粒子總數為N，則函數f(x)在粒子i處的函數和導數的SPH核近似公式(1)～(3)寫成如下形式的SPH粒子求和式.

圖1 SPH核近似計算示意

2 對稱SPH核近似方法（簡稱SSPH方法）

2.1 一維空間中SSPH方法的推導

函數f(x)在x處泰勒級數展開得

由于線性方程組（9）的系數矩陣為對稱矩陣，因此方程組（9）確定的f(x),fx(x),fxx(x)稱為它們的對稱SPH核近似式.

方程組(9)在xi處寫成粒子求和式得解線性方程組（10）可得f(x),fx(x),fxx(x)在xi處用SSPH方法計算的近似值.

2.2 SSPH方法的通用公式的推導

式（11）中忽略二階以上的導數項，兩邊分別乘核函數W和ξ是空間維數有關；并在核函數影響區域? 上取積分得

當二維時，式(12a)～(12c)中函數f、各一階偏導數fx,fy和各二階偏導數fxx,fyyfxy為未知數可以聯立六元線性方程組.當三維時，式(12a)～(12c)中函數f、各一階偏導數和各二階偏導數為未知數可以聯立十元線性方程組.

3 SSPH核近似方法的降元算法

用式（8）和式（12）確定的函數近似值對整個方程組的計算結果帶來一定的影響，因此當f的初始值已知時，無需對它進行核近似并可以省略.本文提出的SSPH方法的降元算法是式（8）和式（12）中的函數核近似過程降元處理，函數f的值作為已知，函數的導數項作為未知數聯立線性方程組.

當一維空間時式(8b)～(8c)中函數一階導數fx(x)和二階導數fxx(x)作為未知數整理可組成二元方程組（13）

當二維空間時，式(12a)～(12b)中函數各一階偏導數fx,fy和各二階偏導數fxx,fyyfxy為未知數整理可以組成一個五元線性方程組.

當三維空間時，式(12a)～(12b)中函數各一階偏導數和各二階偏導數為未知數整理可以組成一個九元線性方程組.

SSPH方法的降元算法不僅保留SSPH方法已有的優點，而且在計算過程中未知數的數量減少一個.雖然降低一個未知數，但是它降低計算量的同時能節省內存空間和縮短計算時間.比如對二維空間而言把原系數矩陣中的36元素降低到25個，使計算量降低30%，同時節省內存空間30%以上.這對大規模數值計算具有一定的意義.

4 算法精度分析

下面以一維空間為例對SSPH核近似方法和降元算法的精度進行分析.

4.1 SSPH核近似方法的精度分析

定理1粒子對稱分布時，SSPH核近似式（9）確定的函數值對核函數光滑長度具有四階精度，函數一階和二階導數對核函數光滑長度具有二階精度；非對稱分布時，函數值對核函數光滑長度具有三階精度，函數一階導數對核函數光滑長度具有二階精度，函數二階導數對核函數光滑長度具有一階精度.

證明因為

當粒子對稱均勻分布時，a0=1，a1=0，a3=0；類此推導，可得f～O(h4)，f x～O(h2)，fxx～O(h2).

4.2 SSPH核近似的降元算法精度分析

定理2 粒子對稱分布時，SSPH核近似的降元算法式（13）確定的函數一階及其二階導數都對核函數光滑長度具有二階精度；非對稱分布時，函數一階導數對核函數光滑長度具有二階精度，函數二階導數對核函數光滑長度具有一階精度.

證明 因為

當粒子對稱分布時，a3=0；類此推導，可得fx～O(h2)，fxx～O(h2).

通過上述分析可知，無論粒子對稱分布還是不對稱分布，本文提出的降元算法與SSPH核近似具有相同的精度，另一方面沒有對函數值進行核近似計算，而直接使用函數正確值，避免了函數核近似計算帶來的新誤差，在計算精度上有一定的提高.

5 數值驗證與比較

5.1 一維函數算例誤差比較

5.1.1 算法誤差比較準則

對不同類型數值方法進行誤差比較時，可選擇的類型有多種多樣，本文對誤差的定義為

5.1.2 誤差比較采用SSPH方法和本文提出的降元算法計算函數f(x)=sin(πx)在區間[0,1]內的一階、二階導數.區間[0,1]上分別均勻分布11、21、41、81和161個粒子,所得的一階、二階導數誤差進行比較，如圖2所示.計算過程都用三次樣條（Cubic B-spline）函數，光滑長度取粒子間距的1.2倍.

圖2 誤差比較

結果表明，每種方法的誤差都隨粒子數的增加而減少；計算一階導數時，SSPH方法和本文提出算法的誤差沒有明顯的區別，但降元算法的誤差小一點，計算二階導數時，本文提出的降元算法誤差小于SSPH方法的誤差.

5.2 熱傳導實例

5.2.1 一維熱傳導實例

一維熱傳導問題

如果初始和邊界條件為T(x,0)=T0,T(0,t)=T0,T(1,t)=T1時，這個問題的解析解是[15]

圖3表示一維熱傳導問題解析解與數值解的對比.同樣的初始和邊界條件下、用本文提出的方法進行計算所得結果與解析解吻合的非常好.

5.2.2 二維熱傳導實例

二維熱傳導問題可用下式來描述，

如果初始和邊界條件分別為T(x,y,0)=10sin(πx)sin(πy)，T(x,0,t)=0，T(x,1,t)=0，T(0,y,t)=0，T(1,y,t)=0時，這個問題的解析解是[15]

圖3 一維熱傳導問題解析解與數值解的比較

圖4 是求解域上均勻分布441個粒子,時間步長取t=5×10?4s，選用三次樣條（B-spline）核函數，光滑長度取粒子間距的1.2倍時，t=0.1s時刻的SSPH降元算法與解析法計算結果對比圖.對比這兩個圖可以發現SSPH降元算法與解析法的計算結果吻合非常好.

圖4 二維熱傳導問題解析解與降元SSPH數值解

圖5 (a)是初始條件T(x,y,0)=0和恒溫邊界條件T(x,0,t)=0,T(x,1,t)=1,T(0,y,t)=0,T(1,y,t)=1，當t=0.1s時刻的計算結果.圖5(b)是初始條件T(x,y,0)=?1和含對流邊界條件T(x,0,t)=1,T(x,1,t)=時刻的計算結果，兩個問題計算結果符合實際熱傳導規律.

圖5 不同邊界條件下二維熱傳導問題的SSPH降元算法數值解

在求解熱傳導控制方程中，二階導數的計算精度非常重要，因此本文提出的降元算法解這類問題時，與其他方法對比具有一定的優勢.

6 結論與展望

本文對SSPH核近似方法的算法進行理論研究和精度對比分析，提出了SSPH核近似方法的降元算法.通過誤差分析和熱傳導問題進行數值計算，得出如下結論.

（1）雖然計算一階和二階導數時，本文提出的SSPH降元算法與SSPH方法具有相同的精度.但是在相等的粒子數，選用相同的核函數和光滑長度的條件下，對同樣問題進行數值計算，發現本方法的誤差小于SSPH方法的誤差.由于本方法沒有對函數值進行核近似計算，而直接使用函數值，因此避免了函數核近似計算所帶來的新誤差，在計算精度上有一定的提高.

（2）計算過程中，雖然僅降低一個未知數，但是它可以節省內存空間，同時降低計算量，可以提高計算效率.這在大規模數值計算和模擬中具有一定的意義和價值.

（3）用本方法對非穩態熱傳導問題進行了數值計算和對比分析，驗證了本方法在實際問題應用中的有效性，有望推廣到其它實際問題.