對幾種SPH 方法的對比研究①

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

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

0 引言

光滑粒子流體動力學(Smoothed Particle Hydrodynamics，簡稱SPH)方法是一種新的無網格純Lagrange粒子方法．1977年，由 Lucy和 Monaghan[1]等人提出，最初用于解決天體物理學和宇宙學中模擬三維無界空間天體的演化問題．Monaghan等人相繼提出了適用于SPH方法的人為粘性和人為熱流項后[2]，又延伸到了計算流體力學的各個領域．Libersky和 Petschek[3]在 SPH 中引入了材料強度模型，模擬了動態(tài)斷裂和破片運動等固體的動態(tài)響應，此后SPH被迅速的應用到固體力學的研究當中．由于SPH方法計算過程不需要劃分網格、容易編程實現、能夠有效處理多種材料界面等特點，發(fā)展到今天，SPH方法已廣泛地應用于科學工程計算和模擬的各個領域．

SPH方法在許多問題上的應用極大地促進了傳統SPH方法的發(fā)展和改進．經典的SPH方法在模擬計算中，存在很多問題．例如，精度比較低，在邊界缺乏插值一致性，不穩(wěn)定性等等，于是，便出現了SPH的各種修正方法，以彌補傳統SPH方法的一些缺點，這些修正導致了SPH方法表達公式的多樣性．

J．K．Chen 等[4]，將函數的 Taylor級數展開用在函數及其空間導數近似上，提出了修正(Corrective)后的光滑粒子流體動力學方法(簡稱CSPH)．R．C．Batra 等[5，6]，通過將核函數及其空間導數乘到函數的Taylor級數展開式兩邊后組成方程組的思路，得到改進(Modified)后的光滑粒子流體動力學方法(簡稱 MSPH)．G．M．Zhang等[7]，通過將核函數及其冪函數序列乘到函數的Taylor級數展開式兩邊后組成方程組的思路，提出了對稱(Symmetric)的光滑粒子流體動力學方法(簡稱SSPH)．本文給出了每種SPH方法的推導過程和表達式，用一維和二維數值實例對這幾種方法進行了詳細的對比和誤差分析．

1 傳統SPH方法的理論基礎

1．1 核近似

在SPH方法中，將任意函數f(X)在空間某一點X處的值和各階(一階和二階)導數值用核函數與f(X)乘積的積分來表示，然后用粒子近似方法將積分轉化成離散形式的粒子求和式．

設函數f(X)，定義在區(qū)域Ω的任一連續(xù)函數，如果X∈Ω，則在Ω中將函數f(X)可以表示成如下積分形式

式(1)中，X為位置矢量，δ(X－X')為狄拉克δ(Dirac)函數．

如果可以找到一個近似于狄拉克函數的核函數W(X － X'，h)，用W(X － X'，h)來代替(1)中的 δ(X－X')函數，則函數f(X)的積分表達式可寫成．

上式簡稱函數的核近似，式(2)中W(X－X'，h)稱為核函數，它依賴于兩點之間的距離|X－X'|和定義核函數影響區(qū)域的光滑長度h．由于核函數W(X －X'，h)不是δ(X －X')函數，因此，式(2)只能是近似式．在光滑粒子法中，用 ＜ ＞來表示函數的核近似式，于是，式(2)可以改寫為

核函數一般若干個性質[8]:

在式(3)中的函數f(X)用函數的空間導數▽f(X)替換，并用分部積分法、高斯定理和核函數性質，可得到函數梯度▽f(X)的核近似式(4)．

同理可得函數的二階梯度▽2f(X)的核近似式(5)．

可證函數核近似式(3)，(4)和(5)對于光滑長度具有二階精度[8]．

1．2 粒子近似法

在SPH方法中把整個求解域離散成一系列任意分布有質量和容量的有限多個粒子，如果核近似式中粒子j處微元體dX'的體積為△Vj、質量為mj、密度為ρj、位置矢量為Xj，記函數f(X)在粒子j處的值為，fj=f(Xj)粒子i為中心的核函數影響范圍內的粒子總數為N，則函數f(X)在粒子i處的函數及其一階、二階導數核近似式(3)、(4)、(5)的核近似粒子求和式分別為．

其中

2 其它格式的SPH方法

2．1 CSPH 方法

為了表示簡便，假設函數f(X)在直角坐標系中以 x1，x2和x3為自變量的函數，即X=(x1，x2，x3)，則函數f(X)在X'處展開泰勒級數得．

其中，當 i=j時 δij=1，當 i≠ j時 δij=0．

式(9)的兩邊乘核函數W(X－X'，h)并在核函數影響區(qū)域Ω內取積分得

因為∫Ω(xi－xi')W(X －X'，h)dX=0并忽略二級及二階以上的導數項可得函數的CSPH核近似式為

式(9)的兩邊乘核函數偏導數▽xjW(X－X'，h)(j=1，2，3)并取積分，忽略二階及二階以上的偏導數項后整理可得

當函數為一維時，整理式(12)直接可得函數一階導數的CSPH核近似式，當函數二維和三維時式(12)分別變成二元和三元方程組．

式(9)的兩邊乘核函數二階偏導數▽xixjW(X－ X'，h)(i，j=1，2，3)在Ω上取積分并忽略二階以上的偏導數項整理得

當函數為一維時，由(13)式直接可得函數的二階導數CSPH核近似式，當函數二維和三維時式(13)分別變成三元和六元方程組．

2．2 MSPH 方法

假設函數f(X)在直角坐標系中以x1，x2和x3為自變量的函數，

其中

忽略二階以上的項函數f(X)在X'處的泰勒級數展開式

式(14)的兩邊分別乘M并在Ω上取積分得

當函數為一維、函數二維和三維函數時，式(15)分別變?yōu)槿?、六元和十元方程組．

2．3 SSPH 方法

假設S=WPT，式(13)的兩邊乘S并在Ω上取積分得

當函數為一維、函數二維和三維函數時，式(17)分別變?yōu)槿?、六元和十元方程組．

3 數值算例

本文中所有數值算例中選用三次B－樣條函數，光滑長度為粒子間距的1．2倍．

3．1 一維算例

考慮函數 f(x)=(x － 0．5)4，x∈[0，1]．

區(qū)間[0，1]上均勻分布21個粒子時，圖1a～1c給出正確值與用SPH，CSPH，MSPH，SSPH方法計算得的函數及其一階、二階導數的核近似值對比圖．從圖1可以看出，在內部區(qū)域每一種方法的計算精度沒有太大的區(qū)別，在邊界處SPH方法和CPSPH方法的精度不如MSPH方法和SSPH方法，特別計算二階導數時，SPH方法和CPSPH方法在邊界處的誤差變的更大，而SSPH方法的誤差最?。?/p>

圖1 函數核近似

對幾種不同類型SPH方法的精度進行更好對別分析，本文引入誤差范數比較，誤差范數可選擇類型有多種多樣，本文中誤差范數的定義為

其中N是整個計算區(qū)域內的粒子數，fExact(xi)分別表示函數及其一階、二階導數在i粒子處的值，fCompute(xi)分別表示響應函數在i粒子處的核近似值．

圖2給出函數f(x)=(x－0．5)4，x∈[0，1]的區(qū)間上均勻分布11，21，41，81，161個粒子時，誤差范數比較，結果表明，4種方法的誤差都隨粒子數的增加而減少，其中SPH方法的誤差最大，SSPH方法的誤差最小，誤差由大到小依次排列為SPH，CSPH，MSPH，SSPH．計算一階和二階導數時，隨粒子數的增大SPH方法誤差的減小不太明顯，而MSPH和SSPH方法的誤差非常接近．對其它類型的函數如 f(x)=sin(πx)，x ∈ [0，1];f(x)=e－x2，x ∈[－ 1，1]進行計算發(fā)現，計算區(qū)間上粒子的分布和光滑長度如何改變，SSPH方法的誤差還是最?。?/p>

圖2 一維情況誤差比較

3．2 二維算例

圖3給出函數 f(x，y)=sin(πx)+sin(πy)，在[0，1]×[0，1]上分別均勻分布 121，441，1681，6561，25921個粒子時，函數值和兩個偏導數值誤差范數進行比較，結果表明SSPH方法誤差還是最?。畬ζ渌愋偷膄(x，y)=(0．5 － xy)4，f(x，y)=x2y2等函數及其偏導數進行了計算和比較，發(fā)現SSPH方法的精度最高．

圖3 二維情況誤差比較

4 結論

本文對幾種SPH方法進行了研究，詳細介紹了每種方法在一維、二維和三維情況下核近似計算公式的推導過程．用每一種方法對一維和二維問題進行數值計算和誤差比較，得出以下幾個結論．

(1)用SPH方法和CSPH進行數值計算時，在邊界處誤差比較大，特別是在計算二階導數時誤差變得更大．用MSPH方法和SSPH進行一階導數時，在邊界處誤差非常小，在邊界處計算二階導數時有明顯地誤差，但SSPH方法的誤差比較?。?/p>

(2)用每一種方法進行誤差分析，對比結果顯示，誤差由大到小依次排列為SPH，CSPH，MSPH，SSPH．

(3)用SSPH方法進行計算時，避免了核函數導數的計算，這樣不僅提高了計算效率，而且放寬了對核函數導數的要求．

[1] R．A．Gingold＆ J．J．Monaghan;Smoothed Particle Hydrodynamics:Theory and Application to Non － Spherical Stars[J]．Mon．Not．R．Astr．Soc．1977，181:375 －389．

[2] J．J．Monaghan;Shock Simu1ation by the Particle Method SPH[J]．Computer Physics，1983，52:374 －389．

[3] L．D．Leberssky，A．G．Petschek;Smoothed Particle Hydrodyna1nics with Strength of Materials[J]．Advances in the Free Lagrange Method，Lecture Notes in Physics，1990，395:248 －257．

[4] J．K．Chen，J．E．Beraun ＆ T．C．Carney;A Corrective Smoothed Particle Method for Boundary Value Problems in Heat Conduction[J]．Comput．Methods App．Mech．Engrg．1999a，23:279－287．

[5] G．M．Zhang＆ R．C．Betra;Modified Smoothed Particle Hydrodynamics Method and Its Application to Transient Problems[J]．Computational Mechanics，2004，34:137 －146．

[6] 陳建設，徐飛，黃其青．光滑粒子流體動力學方法研究[J]．機械強度，2008，30(2):78 －82．

[7] G．M．Zhang＆ R．C．Betra;Symmetric Smoothed Particle Hydrodynamics(SSPH)Method and Its Application to Elastic Problems[J]．Computational Mechanics，2009，43:321 － 340．

[8] G．R．Liu，M．B．Liu．Smoothed Particle Hydrodynamics:A Meshless Particle Method[M]．World Scientific Publishing，2003．