賈 斌，馬志濤，張 偉，龐寶君

(1.哈爾濱工業(yè)大學(xué)航天工程系，哈爾濱150080，jiabin@hit.edu.cn;2.清華大學(xué)航天航空學(xué)院，北京100084)

光滑粒子流體動(dòng)力學(xué)(Smoothed Particle Hydrodynamics，簡稱SPH)法最初是由Lucy L B［1］于1977年提出的一種純Lagrange無網(wǎng)格方法，并在天體物理領(lǐng)域里得到了廣泛應(yīng)用，成功地模擬了超新星爆發(fā)、銀河系的形成和演化、云團(tuán)的碰撞等超大尺度的復(fù)雜問題.該方法的基本思想是將整個(gè)流場的物質(zhì)離散為一系列具有質(zhì)量、速度和能量的“粒子”，然后通過一個(gè)稱為“核函數(shù)”的積分進(jìn)行“核函數(shù)估值”，從而求得流場中不同位置在不同時(shí)刻的各種動(dòng)力學(xué)量.1993年Libersky L D等［2］首先將材料強(qiáng)度效應(yīng)引入SPH方法，成功進(jìn)行了高速碰撞數(shù)值模擬.

但SPH方法中存在著一些內(nèi)在的缺點(diǎn)，其中一個(gè)主要問題便是1995年Swegle J W等人［3］指出的在具有材料強(qiáng)度的問題中存在著拉伸不穩(wěn)定性.為改進(jìn)這一問題，Morris J P［4］提出了應(yīng)用特殊核函數(shù)，因?yàn)槔觳环€(wěn)定與核函數(shù)的二階導(dǎo)數(shù)緊密相關(guān)，盡管此方法在一些情況中非常成功，但對(duì)于一般情況卻不能總得到令人滿意的結(jié)果.Wen Y等［5］及Balsara D S［6］提出了守恒光滑的方法來消除拉伸不穩(wěn)定，但該方法并不能解決所有的情況，而且還增加了計(jì)算時(shí)間.Dyka C T等［7－8］首次在一維算法中引入了應(yīng)力點(diǎn)法，隨后Randles P W等［9］將該方法擴(kuò)展到了多維空間，其中SPH粒子與應(yīng)力粒子交錯(cuò)排列，該方法雖然能夠解決拉伸不穩(wěn)定，但增加了大量的計(jì)算時(shí)間和內(nèi)存.Monaghan J J等［10－11］提出了用于穩(wěn)定計(jì)算的人工力，該力與具體的核函數(shù)有關(guān)，雖然在一定程度上解決了拉伸不穩(wěn)定的問題，但依然存在一定的震蕩.本文提出了一種改進(jìn)形式的人工力，并在人工粘性力的基礎(chǔ)上疊加抗拉伸不穩(wěn)定的人工力，建立了統(tǒng)一的張量形式的人工力.

1 SPH基本方程

考慮材料的彈塑性效應(yīng)，全應(yīng)力張量空間中的SPH插值公式的質(zhì)量守恒方程、動(dòng)量守恒方程、能量守恒方程分別表示為［12－13］

其中:ρ為密度;m為粒子質(zhì)量;v為粒子速度;σαβ為應(yīng)力張量(上角標(biāo)α和β表示張量坐標(biāo));e為比內(nèi)能;x為位置矢量;W為核函數(shù);下標(biāo)i為粒子編號(hào)，下標(biāo)j為i粒子的近鄰粒子編號(hào).核函數(shù)取最常用的三次B-Spline函數(shù)

其中:s=|x－x'|/h，h是光滑長度;在一維、二維、三維中，αd分別為1/h，15/(7πh2)，3/(2πh3);函數(shù)的影響域|x|=2h.

2 拉伸不穩(wěn)定性

2.1 壓縮失穩(wěn)與拉伸失穩(wěn)

所謂拉伸不穩(wěn)定性，是指當(dāng)粒子處于拉伸應(yīng)力狀態(tài)時(shí)，粒子的運(yùn)動(dòng)變得不穩(wěn)定，從而出現(xiàn)粒子凝集和較大的空洞.Swegle J W等［3］曾對(duì)其進(jìn)行了比較細(xì)致的研究，指出這種不穩(wěn)定性是由核函數(shù)的二階導(dǎo)數(shù)和應(yīng)力狀態(tài)的乘積決定的，并提出了拉伸失穩(wěn)產(chǎn)生的充分條件，如下式所示:

式中:W″是核函數(shù)的二階導(dǎo)數(shù).

式(4)的三次B-Spline函數(shù)曲線和一階導(dǎo)數(shù)曲線如圖1所示，圖中s=r/h(r為粒子i和粒子j的間距)，核函數(shù)一階導(dǎo)數(shù)的絕對(duì)值在s=2/3處有極大值.根據(jù)式(5)，如果核函數(shù)的二階導(dǎo)數(shù)為正，則SPH方法在拉伸應(yīng)力狀態(tài)下是不穩(wěn)定的，會(huì)出現(xiàn)拉伸失穩(wěn)現(xiàn)象;反之，則是穩(wěn)定的.因此，在B-Spline函數(shù)中，可以將極值點(diǎn)右側(cè)的區(qū)域稱為拉伸失穩(wěn)區(qū).當(dāng)s＞2/3時(shí)，核函數(shù)一階導(dǎo)數(shù)的絕對(duì)值隨著兩點(diǎn)距離的增大而減小，即鄰近點(diǎn)j對(duì)計(jì)算點(diǎn)i的貢獻(xiàn)隨著距離的增大而減小.當(dāng)s＜2/3時(shí)，核函數(shù)一階導(dǎo)數(shù)的絕對(duì)值隨著兩點(diǎn)距離的減小而減小，即j對(duì)i的貢獻(xiàn)隨著距離的減小反而減小，說明j對(duì)i的貢獻(xiàn)要小于離i更遠(yuǎn)的粒子對(duì)i的貢獻(xiàn)，并且這種趨勢隨兩點(diǎn)距離的減小而加劇，這便是壓縮失穩(wěn)產(chǎn)生的根本原因.因此，可將極值點(diǎn)左側(cè)的區(qū)域稱為壓縮失穩(wěn)區(qū).

圖1 三次B-Spline函數(shù)及其一階導(dǎo)數(shù)曲線

2.2 消除壓縮失穩(wěn)的人工粘性

壓縮失穩(wěn)可通過引入人工粘性來解決.人工粘性的形式很多，但通常采用的形式包含兩種不同類型的粘性，即適用于消除波后振蕩的線性粘性Π =－αρlc▽·v和抑制沖擊波強(qiáng)間斷面數(shù)值振蕩的二次粘性Π=βρl2(▽·v)2.其中:α和β是無量綱系數(shù);l是激波長度傳播尺度;c是聲速.

將上述兩種粘性考慮在一起，Monaghan J J等［14］給出了適合于SPH法的混合型人工粘性

引入人工粘性后，式(2)和式(3)分別變?yōu)?/p>

2.3 改進(jìn)拉伸失穩(wěn)的人工力

Swegle J W等［3］指出拉伸不穩(wěn)定既不依賴于人工粘性，也不依賴于時(shí)間積分項(xiàng)，它是典型的與狀態(tài)方程無關(guān)的問題，因?yàn)樵跔顟B(tài)方程中不產(chǎn)生拉伸應(yīng)力.Monaghan J J等［10－11］提出了用于穩(wěn)定計(jì)算的人工力，該力與具體的核函數(shù)有關(guān)，對(duì)于三階B-Spline函數(shù)，當(dāng)h＜rij≤2h時(shí)，

式中:rij為粒子i和粒子j的間距，h為兩粒子平均光滑長度，σ和ρ分別為應(yīng)力和密度;ε和n無量綱，通常為0.2和4.式(2)和式(3)分別變?yōu)?/p>

當(dāng)rij=2h時(shí)式(8)和(9)不連續(xù)，而且當(dāng)粒子間距rij增大時(shí)Tij大幅減小，難以抵消拉伸不穩(wěn)定的影響.因此雖然在一定程度上解決了拉伸不穩(wěn)定的問題，但依然存在一定的震蕩.為了進(jìn)一步解決拉伸不穩(wěn)定問題，考慮到物體處于彈性拉伸時(shí)彈性模量乘以應(yīng)變即是應(yīng)力的現(xiàn)象，本文提出了如下形式的人工力:

式中，κ通常取0.1～0.3;E為彈性模量;α和β的取值范圍通常為0.2～1.0和1.0～5.0.該人工力是連續(xù)的，其第一項(xiàng)隨粒子間距rij增大而增大，與式(7)相比可更有效地抵消拉伸不穩(wěn)定的影響;第二項(xiàng)為粒子分離時(shí)的人工粘性力［15］.該人工力起作用的條件為

1)粒子間距:當(dāng)h＜rij≤2h時(shí);

2)粒子反向運(yùn)動(dòng):當(dāng)vij·xij＞0時(shí);

3)粒子狀態(tài):相互作用的粒子是固體，并且尚未失效;

4)粒子歸屬:相互作用的粒子屬于同一物體.

對(duì)于二維和三維情況，Owen J M［16］提出了張量形式的人工粘性力，對(duì)消除壓縮失穩(wěn)得到了較好的結(jié)果.本文在該人工粘性力的基礎(chǔ)上疊加抗拉伸不穩(wěn)定的人工力，建立了統(tǒng)一形式的人工力，此時(shí)式(2)和式(3)分別變?yōu)?/p>

式中:統(tǒng)一張量形式的人工力Λαβij=Ταβij－Παβij，張量形式的人工粘性力Παβij有如下定義［16］:

3 算例

3.1 算例1

設(shè)一維鋁桿長 100 mm，左半部分具有－100 m/s的初速度，右半部分靜止.由于拉力作用，應(yīng)力波從中間向兩側(cè)傳播.將該鋁桿進(jìn)行均勻離散，光滑長度取為0.1 mm，采用三階B-Spline核函數(shù)和所提出的人工力，計(jì)算至5 μs.鋁桿材料為 Al 2024，采用 Mie-Grüneisen狀態(tài)方程和Johnson-Cook強(qiáng)度模型進(jìn)行模擬，材料參數(shù)見表1和表2.

表1 Mie-Grüneisen狀態(tài)方程參數(shù)

表2 Johnson-Cook強(qiáng)度模型參數(shù)

Mie-Grüneisen狀態(tài)方程為

其中:pH=ρ0cμ(1+μ)/［1－(s－1)μ］2，eH= (pHμ)/(2ρ0(1+μ))，μ=(ρ/ρ0)－1，Γρ= Γ0ρ0=const，U=c0+sup.式中:p和e分別為靜水壓力和比內(nèi)能;pH和eH分別為沖擊Hugoniot曲線上靜水壓力和比內(nèi)能的參考值;Γ和ρ分別為Grüneisun參數(shù)和密度，相應(yīng)地 Γ0和 ρ0為初始Grüneisun參數(shù)和初始密度;U和up分別為沖擊波波速和波后質(zhì)點(diǎn)速度，c0為體積聲速，s為U和up之間線性關(guān)系的斜率;μ為壓縮比.

Johnson-Cook強(qiáng)度模型為

式中:σy為屈服應(yīng)力;εp為等效塑性應(yīng)變;˙εp為等效塑性應(yīng)變率;參考應(yīng)變率=1 s－1;若T為溫度，TRoom為室溫，TMelt為熔點(diǎn)，則T*=(T－TRoom)/ (TMelt－TRoom);A，B，n，C和m是材料常數(shù).

圖2 5 μs時(shí)桿中應(yīng)力－位置曲線

計(jì)算得到的5 μs時(shí)桿中應(yīng)力 －位置曲線如圖2所示.其中，圖(a)為通過Ansys Autodyn v6.1利用有限元法的計(jì)算結(jié)果，不會(huì)出現(xiàn)拉伸不穩(wěn)定現(xiàn)象;圖(b)中顯示了沒有施加人工力時(shí)SPH法產(chǎn)生了極度不穩(wěn)定現(xiàn)象.當(dāng)加入Monaghan J J等提出的人工力進(jìn)行模擬后，結(jié)果如圖(c)所示，其曲線整體趨勢與圖(a)相同，在一定程度上解決了拉伸不穩(wěn)定的問題，但依然存在一定的振蕩.使用本文的人工力后計(jì)算結(jié)果如圖(d)所示，可以看出得到的曲線平滑且?guī)缀醪划a(chǎn)生振蕩，與有限元結(jié)果非常接近，拉伸不穩(wěn)定性得到了較大改善，其效果優(yōu)于Monaghan J J等提出的人工力.

3.2 算例2

設(shè)一維鋁桿長 100 mm，左半部分具有－100 m/s的初速度，右半部分具有100 m/s的初速度.材料及離散方式與算例1相同，計(jì)算得到的5 μs時(shí)桿中應(yīng)力－位置曲線如圖3所示.

其中，通過Ansys Autodyn v6.1利用有限元法的計(jì)算結(jié)果如圖(a)所示;圖(b)顯示了沒有施加人工力時(shí)SPH法產(chǎn)生的極度不穩(wěn)定現(xiàn)象.當(dāng)加入Monaghan J J等提出的人工力進(jìn)行模擬后，結(jié)果如圖(c)所示.使用本文的人工力后計(jì)算結(jié)果如圖(d)所示.從圖中可以看出，通過施加本文人工力得到的曲線更加平滑且振蕩較小，與有限元結(jié)果接近，其效果優(yōu)于Monaghan J J等提出的人工力.

圖3 5 μs時(shí)桿中應(yīng)力－位置曲線

［1］LUCY L B.A numerical approach to the testing of the fission hypothesis［J］.Astronomical Journal，1977，82: 1013－1020.

［2］LIBERSKY L D，PETSCHEK A G，CARNEY T C，et al. High strain Lagrangian hydrodynamics:a three-dimensional SPH code for dynamic material response［J］.Journal of Computational Physics，1993，109(1):67－75.

［3］SWEGLE J W，HICKS D L，ATTAWAY S W.Smoothed Particle Hydrodynamics stability analysis［J］.Journal of Computational Physics，1995，116(1):123－134.

［4］MORRIS J P.Analysis of Smoothed Particle Hydrodynamics with applications［D］.Melbourne:Monash University，1996.

［5］WEN Y，HICKS D L，SWEGLE J W.Stabilizing SPH with conservative smoothing，SAND94－1932［R］.Albuquerque:Sandia National Laboratories，1994.

［6］BALSARA D S.Von Neumann stability analysis of Smoothed Particle Hydrodynamics－suggestions for optimal algorithms［J］.Journal of Computational Physics，1995，121(2):357－372.

［7］DYKA C T，INGEL R P.An approach for tension instability in Smoothed Particle Hydrodynamics［J］.Computers and Structures，1995，57(4):573－580.

［8］DYKA C T，RANDLES P W，INGEL R P.Stress points for tension instability in Smoothed Particle Hydrodynamics［J］.International Journal for Numerical Methods in Engineering，1997，40(13):2325－2341.

［9］RANDLES P W，LIBERSKY L D.Normalized SPH with stress points［J］.International Journal for Numerical Methods in Engineering，2000，47(10):1445－1462.

［10］MONAGHAN J J.SPH without a tensile instability［J］.Journal of Computational Physics，2000，159(2):290－311.

［11］GRAY J P，MONAGHAN J J，SWIFT R P.SPH elastic dynamics［J］.Computer Methods in Applied Mechanics and Engineering，2001，190(49－50):6641－6662.

［12］CUMMINS S J，RUDMAN M.An SPH projection method［J］.Journal of Computational Physics，1999，152 (2):584－607.

［13］MONAGHAN J J.Extrapolating B-Splines for interpolation［J］.Journal of Computational Physics，1985，60 (2):253－262.

［14］MONAGHAN J J，GINGOLD R A.Shock simulation by the particle method SPH［J］.Journal of Computational Physics，1983，52(2):374－389.

［15］SWEGLE J W，ATTAWAY S W，HEINSTEIN M W.An analysis ofSmoothed Particle Hydrodynamics，SAND93－2513［R］.Albuquerque:Sandia National Laboratories，1994.

［16］OWEN J M.A tensor artificial viscosity for SPH［J］.Journal of Computational Physics，2004，201(2):601－629.