应力波数值计算中的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。