• ISSN 1001-1455  CN 51-1148/O3
  • EI Compendex、CA收录
  • 力学类中文核心期刊
  • 中国科技核心期刊、CSCD统计源期刊

应力波数值计算中的SPH方法

孙晓旺 章杰 王肖钧 李永池 赵凯

引用本文:
Citation:

应力波数值计算中的SPH方法

    作者简介: 孙晓旺(1987—),男,博士研究生,xiaowang@mail.ustc.edu.cn;
  • 基金项目: 爆炸科学与技术国家重点实验室开放基金项目 KFJJ13-9M
    爆炸冲击防灾减灾国家重点实验室开放基金项目 DPMEIKF201401
    国家自然科学基金项目 11202206
    国家自然科学基金项目 11472008
    国家自然科学基金项目 11402266

  • 中图分类号: O383

Application of SPH in stress wave simulation

  • CLC number: O383

  • 摘要: 对一维波动方程的SPH(smoothed particle hydrodynamics)格式和有限差分格式进行比较,并采用SPH法模拟了一维应力/应变波, 获得1个可衡量SPH法模拟应力波准确性的重要指标。结果表明,SPH法模拟应力波传播中采用的光滑长度必须不小于粒子间距;采用B-样条核函数和高斯型核函数能够获得良好的应力波图像,而二次型核函数不能,因此二次型核函数不适用于冲击动力学的数值计算。
  • 图 1  不同γ值下一维应力波在15 μs时的波形

    Figure 1.  Waveforms of one dimensional stress wave at 15 μs under different γ

    图 2  不同γ值下一维应变波在15 μs时的波形

    Figure 2.  Waveforms of one dimensional strain wave at 15 μs under different γ

    表 1  采用B-样条核函数在不同γ值下获得的一维应力波/应变波波速

    Table 1.  One dimensional stress/strain wave velocity obtained by B-spline kernel function using different γ

    γ α c/c0
    应力波 弹性应变波 塑性应变波
    0.6 0.309 0.332 0.324 0.367
    0.7 0.666 0.674 0.670 0.690
    0.8 0.879 0.881 0.878 0.888
    0.9 0.976 0.978 0.976 0.978
    1.0 1.000 1.000 1.000 1.000
    1.1 1.011 1.009 1.009 1.001
    1.2 1.022 1.015 1.020 1.002
    1.3 1.022 1.021 1.021 1.020
    1.4 1.010 1.008 1.007 1.005
    1.5 0.988 0.985 0.985 0.984
    下载: 导出CSV

    表 2  采用高斯型核函数在不同γ值下获得的一维应力波/应变波波速

    Table 2.  One dimensional stress/strain wave velocity obtained by Gaussian kernel function using different γ

    γ α c/c0
    应力波 弹性应变波 塑性应变波
    0.4 0.068 0.075 0.080 0.081
    0.5 0.330 0.339 0.337 0.341
    0.6 0.650 0.656 0.654 0.653
    0.7 0.862 0.867 0.866 0.870
    0.8 0.958 0.962 0.961 0.960
    0.9 0.990 1.000 1.000 1.000
    1.0 0.996 1.001 1.001 1.001
    1.1 1.000 1.001 1.001 1.001
    1.2 1.000 1.001 1.001 1.001
    下载: 导出CSV

    表 3  采用二次型核函数在不同γ值下获得的一维应力波/应变波波速

    Table 3.  One dimensional stress/strain wave velocity obtained by quadratic kernel function using different γ

    γ α c/c0
    应力波 弹性应变波 塑性应变波
    0.6 0.694 0.700 0.701 0.710
    0.7 0.875 0.878 0.878 0.880
    0.8 0.879 0.882 0.882 0.880
    0.9 0.823 0.825 0.826 0.823
    1.0 0.750 0.756 0.756 0.757
    1.1 0.902 0.903 0.904 0.902
    1.2 0.955 0.955 0.956 0.956
    1.3 0.956 0.956 0.956 0.956
    1.4 0.929 0.929 0.930 0.930
    1.5 0.889 0.891 0.891 0.890
    下载: 导出CSV
  • [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. doi: 10.3321/j.issn:1001-1455.2002.02.001
    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. doi: 10.3321/j.issn:1001-1455.2002.02.001
    [3] 卞梁, 王肖钧, 肖卫国, 等.应力波和层裂计算中的光滑粒子法[J].中国科学技术大学学报, 2007, 37(7):706-710, 723. doi: 10.3969/j.issn.0253-2778.2007.07.003
    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. doi: 10.3969/j.issn.0253-2778.2007.07.003
    [4] 章杰, 苏少卿, 郑宇, 等.改进SPH方法在陶瓷材料层裂数值模拟中的应用[J].爆炸与冲击, 2013, 33(4):401-407. doi: 10.3969/j.issn.1001-1455.2013.04.011
    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. doi: 10.3969/j.issn.1001-1455.2013.04.011
    [5] Monaghan J J. Particle methods for hydrodynamics[J]. Computer Physics Reports, 1985, 3(2):71-124. doi: 10.1016/0167-7977(85)90010-3
    [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.
  • 加载中
图(2)表(3)
计量
  • 文章访问数:  928
  • HTML全文浏览量:  172
  • PDF下载量:  654
  • 被引次数: 0
出版历程
  • 收稿日期:  2015-06-16
  • 录用日期:  2015-08-24
  • 刊出日期:  2017-01-25

应力波数值计算中的SPH方法

    作者简介:孙晓旺(1987—),男,博士研究生,xiaowang@mail.ustc.edu.cn
  • 1. 中国科学技术大学近代力学系, 安徽 合肥 230026
  • 2. 北京理工大学科学与技术国家重点实验室, 北京 100081
基金项目:  爆炸科学与技术国家重点实验室开放基金项目 KFJJ13-9M爆炸冲击防灾减灾国家重点实验室开放基金项目 DPMEIKF201401国家自然科学基金项目 11202206国家自然科学基金项目 11472008国家自然科学基金项目 11402266

摘要: 对一维波动方程的SPH(smoothed particle hydrodynamics)格式和有限差分格式进行比较,并采用SPH法模拟了一维应力/应变波, 获得1个可衡量SPH法模拟应力波准确性的重要指标。结果表明,SPH法模拟应力波传播中采用的光滑长度必须不小于粒子间距;采用B-样条核函数和高斯型核函数能够获得良好的应力波图像,而二次型核函数不能,因此二次型核函数不适用于冲击动力学的数值计算。

English Abstract

  • 光滑粒子法(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方法在不同核函数和不同光滑长度条件下模拟一维应力/应变波,研究了核函数和光滑长度对应力波模拟精度的影响,发现了一个影响模拟精度的重要参数,可为提高应力波模拟的精度提供参考。

    • 在Lagrange观点下,连续介质的一维波动方程(包括一维应力/应变波)如下:

      $ \rho \frac{{{\rm{d}}\mathit{u}}}{{{\rm{d}}\mathit{t}}} = \frac{{\partial {\sigma _x}}}{{\partial x}} $

      式中:ρ表示密度,u表示质点速度,t表示时间,σx表示x方向的应力分量。

      采用SPH法对式(1)进行离散,得到一维波动方程的SPH离散格式:

      $ \frac{{{\rm{d}}{\mathit{u}_i}}}{{{\rm{d}}\mathit{t}}} = \frac{1}{{{\rho _i}}}\sum\limits_{j = 1}^N {\left( {{\sigma _{x, j}} - {\sigma _{x, \mathit{i}}}} \right)} \frac{{\partial {W_{\mathit{ij}}}}}{{\partial {x_i}}}\Delta {x_j} $

      式中:下标ij为粒子编号,N是粒子i支持域内的粒子总数,Wij为定义在粒子ij上的核函数,Δxj=mj/ρj表示粒子j在一维空间中的大小,mj为粒子j的质量。

    • 主要研究3种常用核函数及其光滑长度对应力波模拟的影响。这3种核函数分别是B-样条函数、高斯型核函数和二次型核函数。B-样条核函数是J.J.Monaghan等[5]在三次样条函数的基础上提出的,是文献中应用最广泛的核函数。其一维形式如下:

      $ W = \frac{2}{{3\mathit{h}}}\left\{ \begin{array}{l} 1 - \frac{3}{2}{\left( {\frac{r}{h}} \right)^2} + \frac{3}{4}{\left( {\frac{r}{h}} \right)^3}\;\;\;\;\;\;0 \le \frac{r}{h} \le 1\\ \frac{1}{4}{\left( {2 - \frac{r}{h}} \right)^{\rm{3}}}\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;{\rm{1}} < \frac{r}{h} \le 2\\ 0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{r}{h} > 2 \end{array} \right. $

      式中:h为光滑长度,r=|x-xi|是粒子之间的距离。

      高斯型核函数的一维形式如下:

      $ W{\rm{ = }}\frac{1}{{h\sqrt {\rm{ \mathsf{ π} }} }}{{\rm{e}}^{ - {{\left( {\frac{r}{h}} \right)}^2}}} $

      二次型核函数的一维形式为:

      $ W{\rm{ = }}\frac{1}{{2h}}\left\{ \begin{array}{l} \frac{3}{8}{\left( {\frac{r}{h}} \right)^2} - \frac{3}{2}\left( {\frac{r}{h}} \right) + \frac{3}{2}\;\;\;\;\;\;\;\;\;0 \le \frac{r}{h} \le 2\\ 0\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{r}{h} > 2 \end{array} \right. $

      显然,B-样条核函数和二次型核函数的支持域是r≤2h,高斯型核函数的支持域为整个计算域,因为其指数衰减的性质,本文中取r≤3h

    • 采用B-样条核函数时,对粒子i起到作用的粒子在r≤2h的范围内。设初始时刻粒子均匀分布,粒子间距为Δx,取γ=hx。当光滑长度取值范围为0.5<γ≤1时,对式(2)有贡献的粒子只有j=i-1, i, i+1。注意,当γ=1时,粒子j=i±2正好在i支持域边界上,核函数及其导数都为零。将上面3个粒子代入式(2),并代入B-样条核函数的导数,整理后可得:

      $ \frac{{{\rm{d}}{\mathit{u}_i}}}{{{\rm{d}}\mathit{t}}} = {\gamma ^2}{\left( {2 - \frac{1}{\gamma }} \right)^2}\frac{1}{{2{\rho _i}\Delta x}}\left( {{\sigma _{x, i + 1}} - {\sigma _{x, i - 1}}} \right) $

      而采用中心差分方法在空间上对一维波动方程进行离散,得到其有限差分格式:

      $ \frac{{{\rm{d}}{\mathit{u}_i}}}{{{\rm{d}}\mathit{t}}} = \frac{1}{{2{\rho _i}\Delta x}}\left( {{\sigma _{x, i + 1}} - {\sigma _{x, i - 1}}} \right) $

      对比可知式(6)和式(7)之间差一个系数${\gamma ^2}{\left( {2 - \frac{1}{\gamma }} \right)^2}$,记为α

      当光滑长度取值范围为1<γ≤1.5时,粒子i的支持域内有5个粒子,采用相同的分析方法,可以得到:

      $ \frac{{{\rm{d}}{\mathit{u}_i}}}{{{\rm{d}}\mathit{t}}} = \frac{4}{{3{\gamma ^2}}}\left[ {\frac{3}{\gamma } - \frac{9}{{4{\gamma ^2}}} + 6{{\left( {1 - \frac{1}{\gamma }} \right)}^2}} \right]\frac{1}{{2{\rho _i}\Delta x}}\left( {{\sigma _{x, i + 1}} - {\sigma _{x, i - 1}}} \right) $

      对比可知式(8)和式(7)之间差的系数为$\frac{4}{{3{\gamma ^2}}}\left[ {\frac{3}{\gamma } - \frac{9}{{4{\gamma ^2}}} + 6{{\left( {1 - \frac{1}{\gamma }} \right)}^2}} \right]$,若同样记为α,式(6)和式(8)则在形式上完全一致,只是系数α有区别。由此可推断,由式(6)和式(8)获得的结果会受到系数α的影响,因而与有限差分结果有区别。通过对一维应力波和应变波的模拟来具体研究。

      首先模拟端部受峰值为800 MPa、历时12 μs的半正弦脉冲作用的细长杆,得到15 μs后的一维弹性应力波图像,如图 1所示。又模拟了0.1 m厚的飞片以100 m/s的速度撞击同样厚为0.1 m的另一平板产生的一维弹塑性应变波,得到15 μs后的计算结果,如图 2所示。2次模拟的材料相同,作为理想弹塑性处理,性质参数为:材料密度ρ=7800 kg/m3;体积模量K=222.5 GPa;剪切模量G=85.3 GPa;简单拉伸条件下的屈服极限Y0=1.0 GPa;侧限弹性极限Y1=1.97 GPa。

      图  1  不同γ值下一维应力波在15 μs时的波形

      Figure 1.  Waveforms of one dimensional stress wave at 15 μs under different γ

      图  2  不同γ值下一维应变波在15 μs时的波形

      Figure 2.  Waveforms of one dimensional strain wave at 15 μs under different γ

      表 1列出了不同光滑长度下模拟得到的波速,还给出了由式(6)和式(8)确定的系数α,表中c表示模拟得到的波速,c0表示理论波速,c/c0就是量纲一波速。对于B-样条核函数而言,光滑长度是模拟应力波的重要指标,当且仅当光滑长度大于等于粒子间距时,SPH方法给出的应力波计算结果才是准确的;弹性应力波和弹塑性应变波的量纲一波速都与α非常接近,所以α是SPH方法模拟应力波的一个指标,只有α接近1时,SPH方法给出的应力波传播结果才与理论值相符。

      γ α c/c0
      应力波 弹性应变波 塑性应变波
      0.6 0.309 0.332 0.324 0.367
      0.7 0.666 0.674 0.670 0.690
      0.8 0.879 0.881 0.878 0.888
      0.9 0.976 0.978 0.976 0.978
      1.0 1.000 1.000 1.000 1.000
      1.1 1.011 1.009 1.009 1.001
      1.2 1.022 1.015 1.020 1.002
      1.3 1.022 1.021 1.021 1.020
      1.4 1.010 1.008 1.007 1.005
      1.5 0.988 0.985 0.985 0.984

      表 1  采用B-样条核函数在不同γ值下获得的一维应力波/应变波波速

      Table 1.  One dimensional stress/strain wave velocity obtained by B-spline kernel function using different γ

      采用高斯型核函数,取$\frac{r}{h} \le 3$为支持域时,得到的γα关系式如下:

      $ \alpha = \left\{ \begin{array}{l} \frac{4}{{{\gamma ^3}\sqrt {\rm{ \mathsf{ π} }} }}\exp \left( { - \frac{1}{{{\gamma ^2}}}} \right)\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{1}{3} < \gamma \le \frac{2}{3}\\ \frac{4}{{{\gamma ^3}\sqrt {\rm{ \mathsf{ π} }} }}\left[ {\exp \left( { - \frac{1}{{{\gamma ^2}}}} \right) + 4\exp \left( { - \frac{4}{{{\gamma ^2}}}} \right)} \right]\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\;\frac{2}{3} < \gamma \le 1\\ \frac{4}{{{\gamma ^3}\sqrt {\rm{ \mathsf{ π} }} }}\left[ {\exp \left( { - \frac{1}{{{\gamma ^2}}}} \right) + 4\exp \left( { - \frac{4}{{{\gamma ^2}}}} \right) + 9\exp \left( { - \frac{9}{{{\gamma ^2}}}} \right)} \right]\;\;\;\;\;\;\;\;\;{\rm{1}} < \gamma \le \frac{{\rm{4}}}{{\rm{3}}} \end{array} \right. $

      采用高斯型核函数对相同的算例进行模拟,结果见表 2。可以看到,当γ≥0.9时,采用高斯型核函数计算得到的波速结果与理论值符合良好;α作为波速模拟指标,同样适用于高斯型核函数。

      γ α c/c0
      应力波 弹性应变波 塑性应变波
      0.4 0.068 0.075 0.080 0.081
      0.5 0.330 0.339 0.337 0.341
      0.6 0.650 0.656 0.654 0.653
      0.7 0.862 0.867 0.866 0.870
      0.8 0.958 0.962 0.961 0.960
      0.9 0.990 1.000 1.000 1.000
      1.0 0.996 1.001 1.001 1.001
      1.1 1.000 1.001 1.001 1.001
      1.2 1.000 1.001 1.001 1.001

      表 2  采用高斯型核函数在不同γ值下获得的一维应力波/应变波波速

      Table 2.  One dimensional stress/strain wave velocity obtained by Gaussian kernel function using different γ

      采用相同的方法对二次型核函数进行分析,得到二次型核函数γα关系式如下:

      $ \alpha = \left\{ \begin{array}{l} \frac{{\rm{3}}}{{{\rm{2}}{\gamma ^{\rm{2}}}}}\left( {{\rm{1}} - \frac{1}{{2\gamma }}} \right)\;\;\;\;\;\;\;\;\;\;\;\frac{1}{2} < \gamma \le 1\\ \frac{{\rm{3}}}{{{\rm{2}}{\gamma ^{\rm{2}}}}}\left( {3 - \frac{5}{{2\gamma }}} \right)\;\;\;\;\;\;\;\;\;\;\;1 < \gamma \le \frac{3}{2} \end{array} \right.\;\;\; $

      采用二次型核函数模拟相同的算例,结果见表 3。从表 3可以看出,在所考察的γ区间,采用二次型核函数得到的波速会明显小于理论波速,二次型核函数不适用于应力波计算。同时,从表 1~3可以看出,α作为衡量模拟应力波精度的指标的普适性。

      γ α c/c0
      应力波 弹性应变波 塑性应变波
      0.6 0.694 0.700 0.701 0.710
      0.7 0.875 0.878 0.878 0.880
      0.8 0.879 0.882 0.882 0.880
      0.9 0.823 0.825 0.826 0.823
      1.0 0.750 0.756 0.756 0.757
      1.1 0.902 0.903 0.904 0.902
      1.2 0.955 0.955 0.956 0.956
      1.3 0.956 0.956 0.956 0.956
      1.4 0.929 0.929 0.930 0.930
      1.5 0.889 0.891 0.891 0.890

      表 3  采用二次型核函数在不同γ值下获得的一维应力波/应变波波速

      Table 3.  One dimensional stress/strain wave velocity obtained by quadratic kernel function using different γ

    • 从波动方程的SPH离散格式出发,讨论了核函数和光滑长度在应力波数值计算中的作用;通过对一维波的具体计算和分析,获得以下结论:

      (1) 光滑长度是应力波计算中的重要参数,SPH方法模拟应力波时光滑长度不得小于粒子间距;

      (2) B-样条核函数和高斯型核函数都可以得到满意的应力波图像,但是二次型核函数不能,因此二次型函数不适用于应力波计算;

      (3) 另外获得了1个衡量SPH方法模拟应力波的重要指标。

参考文献 (7)

目录

    /

    返回文章
    返回