應力波數值計算中的SPH方法＊

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

（1.中國科學技術大學近代力學系，安徽合肥230026；2.北京理工大學科學與技術國家重點實驗室，北京100081）

應力波數值計算中的SPH方法＊

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

（1.中國科學技術大學近代力學系，安徽合肥230026；2.北京理工大學科學與技術國家重點實驗室，北京100081）

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

爆炸力學；核函數；光滑長度；SPH方法；應力波

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

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

本文中通過比較一維波動方程的SPH格式和有限差格式，并采用SPH方法在不同核函數和不同光滑長度條件下模擬一維應力／應變波，研究了核函數和光滑長度對應力波模擬精度的影響，發現了一個影響模擬精度的重要參數，可為提高應力波模擬的精度提供參考。

1 一維波動方程的SPH離散格式

在Lagrange觀點下，連續介質的一維波動方程（包括一維應力／應變波）如下：

式中：ρ表示密度，u表示質點速度，t表示時間，σx表示x方向的應力分量。

采用SPH法對式（1）進行離散，得到一維波動方程的SPH離散格式：

式中：下標i、j為粒子編號，N是粒子i支持域內的粒子總數，Wij為定義在粒子i、j上的核函數，Δxj＝mj／ρj表示粒子j在一維空間中的大小，mj為粒子j的質量。

2 核函數

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

高斯型核函數的一維形式如下：

二次型核函數的一維形式為：

顯然，B－樣條核函數和二次型核函數的支持域是r≤2h，高斯型核函數的支持域為整個計算域，因為其指數衰減的性質，本文中取r≤3h。

3 不同核函數及其光滑長度對應力波模擬的影響

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

而采用中心差分方法在空間上對一維波動方程進行離散，得到其有限差分格式：

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

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

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

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

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

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

采用高斯型核函數對相同的算例進行模擬，結果見表2?？梢钥吹?，當γ≥0.9時，采用高斯型核函數計算得到的波速結果與理論值符合良好；α作為波速模擬指標，同樣適用于高斯型核函數。

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

采用相同的方法對二次型核函數進行分析，得到二次型核函數γ和α關系式如下：

采用二次型核函數模擬相同的算例，結果見表3。從表3可以看出，在所考察的γ區間，采用二次型核函數得到的波速會明顯小于理論波速，二次型核函數不適用于應力波計算。同時，從表1～3可以看出，α作為衡量模擬應力波精度的指標的普適性。

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

4 結 論

從波動方程的SPH離散格式出發，討論了核函數和光滑長度在應力波數值計算中的作用；通過對一維波的具體計算和分析，獲得以下結論：

（1）光滑長度是應力波計算中的重要參數，SPH方法模擬應力波時光滑長度不得小于粒子間距；

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

（3）另外獲得了1個衡量SPH方法模擬應力波的重要指標。

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

［2］王肖鈞，張剛明，劉文韜，等.彈塑性波計算中的光滑粒子法［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］卞梁，王肖鈞，肖衛國，等.應力波和層裂計算中的光滑粒子法［J］.中國科學技術大學學報，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］章杰，蘇少卿，鄭宇，等.改進SPH方法在陶瓷材料層裂數值模擬中的應用［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國標學科代碼：1303520

A

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

（責任編輯 王易難）

2015－06－16；

2015－08－24

國家自然科學基金項目（11402266，11202206，11472008）；

爆炸科學與技術國家重點實驗室開放基金項目（KFJJ13－9M）；

爆炸沖擊防災減災國家重點實驗室開放基金項目（DPMEIKF201401）

孫曉旺（1987— ），男，博士研究生，xiaowang＠mail.ustc.edu.cn。