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

(1．新疆大學(xué)數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院，新疆烏魯木齊830046．2．新疆大學(xué)機(jī)械工程學(xué)院，新疆烏魯木齊830047)

0 引言

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

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

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

1 傳統(tǒng)SPH方法的理論基礎(chǔ)

1．1 核近似

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

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

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

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

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

核函數(shù)一般若干個(gè)性質(zhì)[8]:

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

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

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

1．2 粒子近似法

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

其中

2 其它格式的SPH方法

2．1 CSPH 方法

為了表示簡(jiǎn)便，假設(shè)函數(shù)f(X)在直角坐標(biāo)系中以 x1，x2和x3為自變量的函數(shù)，即X=(x1，x2，x3)，則函數(shù)f(X)在X'處展開(kāi)泰勒級(jí)數(shù)得．

其中，當(dāng) i=j時(shí) δij=1，當(dāng) i≠ j時(shí) δij=0．

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

因?yàn)椤姚?xi－xi')W(X －X'，h)dX=0并忽略二級(jí)及二階以上的導(dǎo)數(shù)項(xiàng)可得函數(shù)的CSPH核近似式為

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

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

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

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

2．2 MSPH 方法

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

其中

忽略二階以上的項(xiàng)函數(shù)f(X)在X'處的泰勒級(jí)數(shù)展開(kāi)式

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

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

2．3 SSPH 方法

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

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

3 數(shù)值算例

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

3．1 一維算例

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

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

圖1 函數(shù)核近似

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

其中N是整個(gè)計(jì)算區(qū)域內(nèi)的粒子數(shù)，fExact(xi)分別表示函數(shù)及其一階、二階導(dǎo)數(shù)在i粒子處的值，fCompute(xi)分別表示響應(yīng)函數(shù)在i粒子處的核近似值．

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

圖2 一維情況誤差比較

3．2 二維算例

圖3給出函數(shù) f(x，y)=sin(πx)+sin(πy)，在[0，1]×[0，1]上分別均勻分布 121，441，1681，6561，25921個(gè)粒子時(shí)，函數(shù)值和兩個(gè)偏導(dǎo)數(shù)值誤差范數(shù)進(jìn)行比較，結(jié)果表明SSPH方法誤差還是最?。畬?duì)其它類型的f(x，y)=(0．5 － xy)4，f(x，y)=x2y2等函數(shù)及其偏導(dǎo)數(shù)進(jìn)行了計(jì)算和比較，發(fā)現(xiàn)SSPH方法的精度最高．

圖3 二維情況誤差比較

4 結(jié)論

本文對(duì)幾種SPH方法進(jìn)行了研究，詳細(xì)介紹了每種方法在一維、二維和三維情況下核近似計(jì)算公式的推導(dǎo)過(guò)程．用每一種方法對(duì)一維和二維問(wèn)題進(jìn)行數(shù)值計(jì)算和誤差比較，得出以下幾個(gè)結(jié)論．

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

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

(3)用SSPH方法進(jìn)行計(jì)算時(shí)，避免了核函數(shù)導(dǎo)數(shù)的計(jì)算，這樣不僅提高了計(jì)算效率，而且放寬了對(duì)核函數(shù)導(dǎo)數(shù)的要求．

[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] 陳建設(shè)，徐飛，黃其青．光滑粒子流體動(dòng)力學(xué)方法研究[J]．機(jī)械強(qiáng)度，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．