孫曉旺，章 杰，王肖鈞，李永池，趙 凱，2

（1.中國科學(xué)技術(shù)大學(xué)近代力學(xué)系，安徽合肥230026；2.北京理工大學(xué)科學(xué)與技術(shù)國家重點(diǎn)實(shí)驗(yàn)室，北京100081）

應(yīng)力波數(shù)值計(jì)算中的SPH方法＊

孫曉旺1，章 杰1，王肖鈞1，李永池1，趙 凱1，2

（1.中國科學(xué)技術(shù)大學(xué)近代力學(xué)系，安徽合肥230026；2.北京理工大學(xué)科學(xué)與技術(shù)國家重點(diǎn)實(shí)驗(yàn)室，北京100081）

對(duì)一維波動(dòng)方程的SPH（smoothed particle hydrodynamics）格式和有限差分格式進(jìn)行比較，并采用SPH法模擬了一維應(yīng)力／應(yīng)變波，獲得1個(gè)可衡量SPH法模擬應(yīng)力波準(zhǔn)確性的重要指標(biāo)。結(jié)果表明，SPH法模擬應(yīng)力波傳播中采用的光滑長度必須不小于粒子間距；采用B－樣條核函數(shù)和高斯型核函數(shù)能夠獲得良好的應(yīng)力波圖像，而二次型核函數(shù)不能，因此二次型核函數(shù)不適用于沖擊動(dòng)力學(xué)的數(shù)值計(jì)算。

爆炸力學(xué)；核函數(shù)；光滑長度；SPH方法；應(yīng)力波

光滑粒子法（smoothed particle hydrodynamics，SPH）［1］是一種純Lagrange的無網(wǎng)格方法，通過引入光滑粒子和核函數(shù)，將空間和方程離散。作為一種Lagrange型粒子方法，SPH法不需要網(wǎng)格，比有限元或有限差分更適用于模擬沖擊動(dòng)力學(xué)問題。在應(yīng)力波數(shù)值計(jì)算方面，王肖鈞等［2］采用SPH法模擬了一維彈塑性應(yīng)變波，證明SPH在應(yīng)力波模擬中的適用性；卞梁等［3］模擬了鋁鋰合金材料中的應(yīng)變波和層裂，顯示了比有限差更高的精確度；章杰等［4］采用改進(jìn)的SPH法模擬了陶瓷材料層裂。

核函數(shù)確定了SPH方法中的加權(quán)方法，學(xué)者們提出了不同的核函數(shù)以適應(yīng)不同問題的模擬。J.J.Monaghan等［5］在三次樣條函數(shù)的基礎(chǔ)上提出了B－樣條核函數(shù)，它是文獻(xiàn)中應(yīng)用最廣泛的核函數(shù)；R.A.Gingold等［6］最早采用高斯型核函數(shù)模擬非球形星體；G.R.Johnson等［7］首次采用二次型核函數(shù)模擬高速?zèng)_擊問題。這3種核函數(shù)是目前沖擊動(dòng)力學(xué)中最常用的核函數(shù)。核函數(shù)及其光滑長度確定了粒子的加權(quán)方式及支持域大小，在SPH方法中占有重要的地位。但是在大多數(shù)研究中，核函數(shù)及光滑長度的選擇具有很大的隨意性，它們對(duì)模擬精度的影響沒有被提及。

本文中通過比較一維波動(dòng)方程的SPH格式和有限差格式，并采用SPH方法在不同核函數(shù)和不同光滑長度條件下模擬一維應(yīng)力／應(yīng)變波，研究了核函數(shù)和光滑長度對(duì)應(yīng)力波模擬精度的影響，發(fā)現(xiàn)了一個(gè)影響模擬精度的重要參數(shù)，可為提高應(yīng)力波模擬的精度提供參考。

1 一維波動(dòng)方程的SPH離散格式

在Lagrange觀點(diǎn)下，連續(xù)介質(zhì)的一維波動(dòng)方程（包括一維應(yīng)力／應(yīng)變波）如下：

式中：ρ表示密度，u表示質(zhì)點(diǎn)速度，t表示時(shí)間，σx表示x方向的應(yīng)力分量。

采用SPH法對(duì)式（1）進(jìn)行離散，得到一維波動(dòng)方程的SPH離散格式：

式中：下標(biāo)i、j為粒子編號(hào)，N是粒子i支持域內(nèi)的粒子總數(shù)，Wij為定義在粒子i、j上的核函數(shù)，Δxj＝mj／ρj表示粒子j在一維空間中的大小，mj為粒子j的質(zhì)量。

2 核函數(shù)

主要研究3種常用核函數(shù)及其光滑長度對(duì)應(yīng)力波模擬的影響。這3種核函數(shù)分別是B－樣條函數(shù)、高斯型核函數(shù)和二次型核函數(shù)。B－樣條核函數(shù)是J.J.Monaghan等［5］在三次樣條函數(shù)的基礎(chǔ)上提出的，是文獻(xiàn)中應(yīng)用最廣泛的核函數(shù)。其一維形式如下：

高斯型核函數(shù)的一維形式如下：

二次型核函數(shù)的一維形式為：

顯然，B－樣條核函數(shù)和二次型核函數(shù)的支持域是r≤2h，高斯型核函數(shù)的支持域?yàn)檎麄€(gè)計(jì)算域，因?yàn)槠渲笖?shù)衰減的性質(zhì)，本文中取r≤3h。

3 不同核函數(shù)及其光滑長度對(duì)應(yīng)力波模擬的影響

采用B－樣條核函數(shù)時(shí)，對(duì)粒子i起到作用的粒子在r≤2h的范圍內(nèi)。設(shè)初始時(shí)刻粒子均勻分布，粒子間距為Δx，取γ＝h／Δx。當(dāng)光滑長度取值范圍為0.5＜γ≤1時(shí)，對(duì)式（2）有貢獻(xiàn)的粒子只有j＝i－1，i，i＋1。注意，當(dāng)γ＝1時(shí)，粒子j＝i±2正好在i支持域邊界上，核函數(shù)及其導(dǎo)數(shù)都為零。將上面3個(gè)粒子代入式（2），并代入B－樣條核函數(shù)的導(dǎo)數(shù)，整理后可得：

而采用中心差分方法在空間上對(duì)一維波動(dòng)方程進(jìn)行離散，得到其有限差分格式：

當(dāng)光滑長度取值范圍為1＜γ≤1.5時(shí)，粒子i的支持域內(nèi)有5個(gè)粒子，采用相同的分析方法，可以得到：

首先模擬端部受峰值為800MPa、歷時(shí)12μs的半正弦脈沖作用的細(xì)長桿，得到15μs后的一維彈性應(yīng)力波圖像，如圖1所示。又模擬了0.1m厚的飛片以100m／s的速度撞擊同樣厚為0.1m的另一平板產(chǎn)生的一維彈塑性應(yīng)變波，得到15μs后的計(jì)算結(jié)果，如圖2所示。2次模擬的材料相同，作為理想彈塑性處理，性質(zhì)參數(shù)為：材料密度ρ＝7 800kg／m3；體積模量K＝222.5GPa；剪切模量G＝85.3GPa；簡單拉伸條件下的屈服極限Y0＝1.0GPa；側(cè)限彈性極限Y1＝1.97GPa。

圖1 不同γ值下一維應(yīng)力波在15μs時(shí)的波形Fig.1 Waveforms of one dimensional stress wave at 15μs under differentγ

圖2 不同γ值下一維應(yīng)變波在15μs時(shí)的波形Fig.2 Waveforms of one dimensional strain wave at 15μs under differentγ

表1列出了不同光滑長度下模擬得到的波速，還給出了由式（6）和式（8）確定的系數(shù)α，表中c表示模擬得到的波速，c0表示理論波速，c／c0就是量綱一波速。對(duì)于B－樣條核函數(shù)而言，光滑長度是模擬應(yīng)力波的重要指標(biāo)，當(dāng)且僅當(dāng)光滑長度大于等于粒子間距時(shí)，SPH方法給出的應(yīng)力波計(jì)算結(jié)果才是準(zhǔn)確的；彈性應(yīng)力波和彈塑性應(yīng)變波的量綱一波速都與α非常接近，所以α是SPH方法模擬應(yīng)力波的一個(gè)指標(biāo)，只有α接近1時(shí)，SPH方法給出的應(yīng)力波傳播結(jié)果才與理論值相符。

表1 采用B－樣條核函數(shù)在不同γ值下獲得的一維應(yīng)力波／應(yīng)變波波速Table 1 One dimensional stress／strain wave velocity obtained by B－spline kernel function using differentγ

采用高斯型核函數(shù)對(duì)相同的算例進(jìn)行模擬，結(jié)果見表2?？梢钥吹?，當(dāng)γ≥0.9時(shí)，采用高斯型核函數(shù)計(jì)算得到的波速結(jié)果與理論值符合良好；α作為波速模擬指標(biāo)，同樣適用于高斯型核函數(shù)。

表2 采用高斯型核函數(shù)在不同γ值下獲得的一維應(yīng)力波／應(yīng)變波波速Table 2 One dimensional stress／strain wave velocity obtained by Gaussian kernel function using differentγ

采用相同的方法對(duì)二次型核函數(shù)進(jìn)行分析，得到二次型核函數(shù)γ和α關(guān)系式如下：

采用二次型核函數(shù)模擬相同的算例，結(jié)果見表3。從表3可以看出，在所考察的γ區(qū)間，采用二次型核函數(shù)得到的波速會(huì)明顯小于理論波速，二次型核函數(shù)不適用于應(yīng)力波計(jì)算。同時(shí)，從表1～3可以看出，α作為衡量模擬應(yīng)力波精度的指標(biāo)的普適性。

表3 采用二次型核函數(shù)在不同γ值下獲得的一維應(yīng)力波／應(yīng)變波波速Table 3 One dimensional stress／strain wave velocity obtained by quadratic kernel function using differentγ

4 結(jié) 論

從波動(dòng)方程的SPH離散格式出發(fā)，討論了核函數(shù)和光滑長度在應(yīng)力波數(shù)值計(jì)算中的作用；通過對(duì)一維波的具體計(jì)算和分析，獲得以下結(jié)論：

（1）光滑長度是應(yīng)力波計(jì)算中的重要參數(shù)，SPH方法模擬應(yīng)力波時(shí)光滑長度不得小于粒子間距；

（2）B－樣條核函數(shù)和高斯型核函數(shù)都可以得到滿意的應(yīng)力波圖像，但是二次型核函數(shù)不能，因此二次型函數(shù)不適用于應(yīng)力波計(jì)算；

（3）另外獲得了1個(gè)衡量SPH方法模擬應(yīng)力波的重要指標(biāo)。

［1］Liu G R，Liu M B.Smoothed particle hydrodynamics：A meshfree particle method［M］.Singapore：World Scientific，2003：309－341.

［2］王肖鈞，張剛明，劉文韜，等.彈塑性波計(jì)算中的光滑粒子法［J］.爆炸與沖擊，2002，22（2）：97－103.Wang Xiaojun，Zhang Gangming，Liu Wentao，et al.Computations of elastic－plastic waves by smoothed particle hydrodynamics［J］.Explosion and Shock Waves，2002，22（2）：97－103.

［3］卞梁，王肖鈞，肖衛(wèi)國，等.應(yīng)力波和層裂計(jì)算中的光滑粒子法［J］.中國科學(xué)技術(shù)大學(xué)學(xué)報(bào)，2007，37（7）：706－710，723.Bian Liang，Wang Xiaojun，Xiao Weiguo，et al.Numerical simulation of stress waves and spallation by smoothed particle hydrodynamics［J］.Journal of University of Science and Technology of China，2007，37（7）：706－710，723.

［4］章杰，蘇少卿，鄭宇，等.改進(jìn)SPH方法在陶瓷材料層裂數(shù)值模擬中的應(yīng)用［J］.爆炸與沖擊，2013，33（4）：401－407.Zhang Jie，Su Shaoqing，Zheng Yu，et al.Application of modified SPH method to numerical simulation of ceramic spallation［J］.Explosion and Shock Waves，2013，33（4）：401－407.

［5］Monaghan J J.Particle methods for hydrodynamics［J］.Computer Physics Reports，1985，3（2）：71－124.

［6］Gingold R A，Monaghan J J.Smoothed particle hydrodynamics：Theory and application to non－spherical stars［J］.Monthly Notices of the Royal Astronomical Society，1977，181（2）：375－389.

［7］Johnson G R，Stryk R A，Beissel S R.SPH for high velocity impact computations［J］.Computer Methods in Applied Mechanics and Engineering，1996，139（1／2／3／4）：347－373.

Application of SPH in stress wave simulation

Sun Xiaowang1，Zhang Jie1，Wang Xiaojun1，Li Yongchi1，Zhao Kai1，2
（1.Department of Modern Mechanics，University of Science and Technology of China，Hefei 230026，Anhui，China；2.State Key Laboratory of Explosion Science and Technology，Beijing Institute of Technology，Beijing100081，China）

Obtaining accurate waveforms is significant in impact mechanics numerical calculation.This paper is to analyze how the kernel functions and smooth length affect the result of stress wave simulation.The SPH（smoothed particle hydrodynamics）formulations with different kernel functions and smooth lengths of one dimensional wave equation was compared with the finite difference formulation，which was derived in this paper.One dimensional stress and strain waves were simulated using the SPH method with different kernel functions and smooth lengths，and waveforms were gained accurately by B－spline and Gaussian kernels when the smooth length was equal to or greater than the particle interval.The wave velocity obtained by the quadratic kernel is below the theoretical value，no matter what the smooth length is.A parameter was deduced in this paper as roughly equal to the dimensionless wave velocity.Several conclusions were drawn.Firstly，the smooth length is equal to or greater than the particle interval，which is the necessary prerequisite for accurate stress wave simulation with SPH.Then，the quadratic kernel is not suitable in impact mechanics numerical calculation.Finally，the parameter deduced in this paper is a significant index to evaluate the stress wave simulation result of SPH.

mechanics of explosion；kernel function；smooth length；SPH method；stress wave

O383國標(biāo)學(xué)科代碼：1303520

A

10.11883／1001－1455（2017）01－0010－05

（責(zé)任編輯 王易難）

2015－06－16；

2015－08－24

國家自然科學(xué)基金項(xiàng)目（11402266，11202206，11472008）；

爆炸科學(xué)與技術(shù)國家重點(diǎn)實(shí)驗(yàn)室開放基金項(xiàng)目（KFJJ13－9M）；

爆炸沖擊防災(zāi)減災(zāi)國家重點(diǎn)實(shí)驗(yàn)室開放基金項(xiàng)目（DPMEIKF201401）

孫曉旺（1987— ），男，博士研究生，xiaowang＠m(xù)ail.ustc.edu.cn。