SPH方法在聚能装药射流三维数值模拟中的应用＊

李 磊，沈兆武，李学岭，倪小军

（1.中国科学技术大学近代力学系，安徽 合肥 230027；2.安徽方圆机电股份有限公司，安徽 蚌埠 233010）

聚能装药技术广泛应用于弹药武器设计和工业生产行业，例如破甲弹、聚能切割、石油开采等领域。聚能装药起爆后形成聚能射流过程的数值模拟对聚能装药技术的研究及设计具有重要的意义。聚能射流形成过程是弹塑性流体动力学问题，对于整个过程进行数值模拟是一个难点。当采用Lagrange有限元法模拟射流形成过程时，必然遇到网格发生大畸变和滑移面处理等一系列问题，最终导致计算精度降低甚至计算终断。Euler方法不存在网格发生大畸变问题，但难以准确描述各类界面。光滑粒子动力学（Smoothed particle hydrodynamics，简称SPH）方法是一种Lagrange无网格数值方法，它采用带质量、动量、能量的粒子构成离散计算域，不同材料的粒子自然地构成界面，材料间的相互作用可以由粒子间的相互作用来自然地模拟，材料的变形不依赖于网格而通过粒子的运动来描述，因此在理论上，SPH方法能够“自然地”模拟侵彻贯穿、高速碰撞等物理现象。SPH方法最早由L.B.Lucy［1］和R.A.Gingold等［2］于1977年提出并应用于天体物理学，后来被广泛地应用于力学研究的各个领域，J.W.Swegle等［3］首先用SPH方法模拟爆炸问题，M.B.Liu［4－6］等应用SPH方法模拟了聚能装药的爆轰过程。

本文中，首先介绍SPH方法的基本理论，其次采用有限元程序LS－DYNA中的归一化粒子近似公式对聚能装药射流形成过程进行三维数值模拟，并对射流形成过程和机理进行分析；然后将SPH数值模拟结果与二维有限元数值模拟（FEM方法）结果及理论计算结果相比较，讨论SPH方法的精确性、有效性以及相对于FEM方法的优越性；最后分析和比较不同初始相邻粒子数对于计算结果的影响。

1 光滑粒子法

SPH方法是通过构造一个核函数对离散质点位置进行核估计来计算梯度的相关项，因此不需要网格来求解偏微分的差分，只需将微分形式的守恒方程转化为积分方程形式，计算出在任意一点上的各个场变量的核估计。函数f（x）在空间某一点x上的核估计可以通过函数f（x）在域Ω中的积分获得

式中：W（x－x′，h）为核函数或权函数，x为估计点的空间坐标，x′为对估计点有贡献作用的空间点坐标，h为紧支域的度量尺寸即光滑长度。

函数导数的核估计可通过将式（1）里的函数f（x）视为导数∂f（x）／∂x，并利用分步积分和核函数W在积分域Ω边界上为零的条件求得

由此可见光滑粒子法的基本思想之一是将函数导数的核估计转换成核函数的导数，核函数是预先设定的已知函数。将解域Ω划分为M个子域粒子，每个子域粒子j的质量和密度分别为mj＝m（xj），ρj＝ρ（xj），设f（x）在粒子i、j上的值分别为fi＝f（xi），fj＝f（xj），则f（x）及其导数在粒子i上的核估计式（1）～（2）的离散式为

式中：Wij＝W（xj－xi，h）为离散核函数，积分微体元dx′＝mj／ρj，α表示空间维序数，下标j为编号为i粒子的临近粒子编号。

考虑人工粘性Πij和人工热流Hi的影响，具有材料强度的全应力张量空间中的SPH插值公式可表示为［7－8］

式中：u、σ分别表示速度和应力张量，上标α、β表示空间方向；E为比内能；Γ、ε分别为偏应力和应变张量为狄拉克函数，当α＝β时为1，当α≠β时为0。

核函数是SPH方法中核心组成部分，常用的有B－spline核函数、Gauss核函数、二次核函数及指数核函数等，其中B－spline核函数是目前应用最广泛的核函数，其形式为［8］

2 数值模拟与分析

2.1 计算模型

图1 四分之一聚能装药模型（剖面图）Fig.1The quarter model of the shaped charge（profile map）

首先利用ANSYS有限元软件建立聚能装药三维模型并划分为Lagrange单元网格，根据对称性只建立四分之一模型。其次，生成PART以及最终求解文件（K文件）。然后，修改K文件：删除Lagrange单元信息，添加SPH质点及其相关属性设置，图1所示为修改K文件之后得到的光滑粒子计算模型。最后，利用LS＿DYNA 971求解器求解K文件。然而，采用这种方式建立SPH模型时需要注意实际长度和建模长度的区别［9］。

模型中炸药配置了151 827个粒子，药型罩配置了23 816个粒子。粒子间距为0.05～0.06cm，初始光滑长度为0.08cm。炸药粒子的最小、最大光滑长度分别为0.5倍和2.0倍初始光滑长度。药型罩粒子的最小、最大光滑长度分别为0.01倍和1.20倍初始光滑长度。每个粒子初始相邻粒子数设为1 200，时间步长系数为0.015。

炸药采用 MAT＿HIGH＿EXPLOSIVE＿BURN高能燃烧模型，状态方程为JWL［9－10］

式中：p为爆轰压力，E为炸药的比内能，V为相对比容，状态方程参数A＝371.2GPa，B＝3.21GPa，R1＝4.15，R2＝0.95，ω＝0.3，炸药密度ρ＝1.60kg／m3，C－J爆速vD＝7.74km／s。

药型罩材料为紫铜，采用 MAT＿JOHNSON＿COOK材料模型，状态方程为 Mie－Grüneisen［9－10］

式中：vs－vp曲线截距C＝0.394；vs－vp曲线斜率的相关因数S1＝1.49，S2＝0，S3＝0，vs为冲击波速度，vp为质点速度；Grüneisen因数γ0＝2.02，Grüneisen因数的一阶体积修正因数α′＝0.47，E′为材料单位体积内能，体积变化率分别为紫铜的瞬时和初始密度，＝8.93kg／m3，剪切模量γ＝47.7GPa，厚度d0＝0.18cm，外直径D＝4.0cm，顶角ξ＝53°，熔化温度T＝1 360K。

2.2 射流形成过程和机理分析

以药型罩内、外壁上编号分别为154 202和164 459的节点（粒子）为例，在不同时刻粒子速度大小的SPH数值模拟结果与FEM方法模拟结果的对比如图2～3中曲线所示。由图2～3可知，这2种方法的数值模拟结果十分接近。下面以SPH数值模拟结果为例简单地分析射流形成过程和机理。

图2 粒子（154 202）速度与FEM方法计算结果比较图Fig.2The velocity of the particle（154 202）calculated by SPH （3D）and finite element method（2D）

图3 粒子（164 459）速度与FEM方法计算结果比较图Fig.3The velocity of the particle（164 459）calculated by SPH and finite element method（2D）

当炸药起爆后产生球形爆轰波，爆轰波首先到达药型罩的顶端，当爆轰波前沿的冲击波传入药型罩后，驱使药型罩向轴线方向闭合，并在轴线处产生高压、高温、高速射流以及速度相对较低的杵体。粒子154 202于t＝5μs时刻开始向轴线方向闭合，由于该粒子在向轴线闭合过程中，不可避免地与其他轴线方向运动的粒子相互碰撞挤压，所以该粒子加速度是变化的，碰撞后粒子与罩内壁的其他粒子形成能量密度很高的射流。由图2中粒子154 202的速度曲线可以看出，当t＝8.8μs时刻，冲击波将越过粒子，其加速度逐渐降为0，而速度达到最大值为5 448m／s，然后在应力波的作用下粒子的速度发生小幅度震荡，并以约5 210m／s的平均速度沿着轴线方向运动。由图3可知，粒子164 459在t＝6.6μs时刻向轴线方向闭合运动，并与相邻粒子相互碰撞和挤压变形；当t＝12.0μs时速度达到最大值为1 184m／s，之后粒子速度逐渐降低并继续向轴线方向闭合运动；当t＝25.0μs时闭合过程完成，与药型罩外壁的其它粒子形成能量密度低的杵体，并以相对稳定的速度（约429m／s）跟随射流向前运动。由于射流头部至尾部存在速度梯度，且射流头部速度远大于杵体速度，所以射流在运动过程中不断拉伸，并在拉伸过程中产生颈缩与断裂，同时与杵体逐渐分离。当采用SPH法计算时，计算过程十分稳定，不存在网格畸变问题。而当采用Lagrange有限元方法计算时，在t＝100.0μs时刻网格发生了严重的扭曲缠绕变形，并导致计算中断。此外，采用FEM方法计算多介质问题时，常使用接触算法，然而该算法中参数较多，很难把握，经常出现穿透现象，如图4（b）所示，进而导致计算终断，而采用SPH方法计算多介质问题时，2种不同材料之间同样计算质点近似，交界面处不需要特别的处理，所以穿透现象很少发生，如图4（a）所示。因此，SPH方法比FEM方法更适合模拟此类的多介质、大变形问题。

2.3 数值模拟结果比较与分析

图4 药型罩顶部与炸药交界面处数值模拟图Fig.4The simulation images in the interface between the liner’s top and the main charge

2.3.1 射流头部速度和长度的理论计算

根据PER理论可以计算药型罩射流头部速度［11］

式中：vj为射流头部速度，θ为药型罩半锥角，φ、φ＋分别为动态压合角、稳态压合角，x为罩微元在轴线方向坐标（罩顶点为原点），v为指数压合速度，vx′为v对x的偏导数，t0为罩微元初始压垮时刻，t为罩微元压垮至轴线的时刻，v0为常数压合速度，τ为指数加速模型中的特征时间常数［12］

式中：m为单位面积药型罩的初始质量，pCJ＝28GPa，为C－J爆轰压力，c1、c2为经验参数。

对v0采用格尼模型计算［12］

式中：Q＝1.5（3χ2＋4χ＋1）／（4χ2＋4χ＋1），χ＝R0／Rr，μ＝mL／me，κ≈C0pCJτ（R0－Rr）［（R0／Rr）－1］，Rr、R0分别为罩微元x处炸药内表面和外表面半径，me、mL分别药型罩微元x处的炸药和罩的质量，C0为经验参数，E为炸药比内能。

射流长度的计算采用P.C.Chou等［11］提出的公式

式中：r（x，t）为坐标为x的罩微元在t时刻的位置，ζ为抛射角，t′为x罩微元刚到达轴线上的时间，L为罩顶部炸药高度，η（x，t）为x罩微元在t时刻的轴向应变。

2.3.2 SPH数值模拟结果、FEM方法模拟结果及理论计算结果的比较和分析

［11－12］，将x＝1.076 9cm处，θ＝26.5°，R0＝1.345 1cm，Rr＝0.560 6，μ＝0.629，m＝9.26g／cm2，me＝1.196g／cm2，φ＋＝43.3°［13］，φ＝55°［13］，L＝3cm 等计算参数代入式（11）～（18），计算射流头部速度vj（即粒子154 202速度）和射流长度lj，vj的理论计算结果为5 332m／s，FEM 方法计算结果为5 197m／s，SPH 方法的计算结果为5 210m／s，可知分别采用SPH方法和FEM方法计算的射流头部速度的误差分别为2.29%、2.53%；分别计算的射流长度的比较结果，如表1所示，2种方法计算的射流长度的平均误差分别为0.143%、4.425%，与理论计算结果的误差均小于5%，因此可认为SPH方法模拟聚能装药问题是可行的且具有更好的精度。

表1 不同时刻射流长度计算值Table 1Length of the jet at different time calculated by different method

2.4 光滑长度和初始相邻光滑粒子数对计算结果的影响分析

光滑长度h对于计算精度和效率至关重要。在聚能装药爆炸形成聚能射流过程中，炸药及药型罩物质变形速度和变形量很大。若采用固定光滑长度，在炸药爆轰膨胀时，核函数影响域内的粒子数量会不断减少，可能出现数值断裂等问题，而当药型罩在爆轰波作用下产生压缩变形时，核函数影响域内的粒子数会急剧增加，这将显着地降低计算效率，而且可能产生压缩不稳定性［14］。因此采用变光滑长度的时间积分格式［10］

式中：h（t）为与时间相关的光滑长度，d表示空间维度，vp为粒子速度。

图5 初始相邻粒子数为200和300时药型罩塑性应变云图Fig.5Plastic－strain contour of the liner when Nis 200or 300

初始相邻光滑粒子数量N是由初始光滑长度h决定的，由于初始时刻粒子间距的非均匀性，导致每个粒子的紧支域内的初始相邻粒子数不相同，即每个粒子的紧支域内包含不同的质量，这将造成非一致的核近似精度。因此在采用LS＿DYNA程序模拟时，需要设定一个一致的初始相邻粒子数量N用来调整初始化阶段的内存分配。在计算过程中发现变量N的设置除了影响计算效率外，对于计算精度也会产生一定影响。从能量角度分析并比较不同初始相邻粒子数量（N＝200～3 500）对计算结果的影响，计算模型与2.1节相同。当N设为200和300时，药型罩的边界自由面出现了“粒子飞散”的现象，如图5所示，而且在计算十几微秒之后，时间步长逐渐减小，最终导致计算困难。这可能是由于初始相邻粒子数量太少，药型罩边缘粒子的球形紧支域内没有足够的其他近似粒子，即光滑函数在该粒子的紧支域内呈现出严重的非正则化性质，进而影响了计算精度。

为了避免以上出现的问题，增加初始相邻粒子数量，当N≥400时，“粒子飞散”的现象不再发生，射流成型质量较好，计算过程稳定。药型罩的内能e和动能ek时程曲线可知如图6～7所示，炸药爆轰能量主要以动能形式传递给药性罩，而药型罩的内能增加很小。当N取不同值时，计算得到的药型罩内能和动能在任意时刻的最大误差分别小于5%和1%，因此当N＝400～3 500时，计算结果之间的误差较小。将N取不同值时计算得到的射流头部速度vj与理论计算结果相比较如图8所示，从图中可以看出，当N≤300时误差最大，大于7.5% ；当400≤N≤3 500时误差小于4.5%。值得注意的是随N值的增加，计算量也会增加，当N＞1 800时计算效率大大地降低，而且从图8的曲线走势可以看出当N≥1 800时误差呈上升趋势，因此N的取值不宜太大。

图6 药型罩（1／4模型）内能时程曲线Fig.6Internal energy curve of the liner

图7 药型罩（1／4模型）动能时程曲线Fig.7Kinetic energy curve of the liner

图8 不同N值计算的射流头部速度与理论结果的误差Fig.8Error between jet tip’s velocity calculated with different Nand theoretical results

3 结 论

（1）SPH方法与FEM方法相比，避免了网格扭曲、缠绕和物质穿透等问题，计算过程更稳定，更适合诸如聚能装药等多介质、大变形问题。

（2）采用SPH法和FEM方法计算的射流头部速度的误差分别为2.29%、2.53%，射流长度的平均误差分别为0.143%、4.425%，SPH 方法的计算结果更接近理论计算值。

（3）初始相邻粒子数N影响计算精度和计算效率：当N≤300时，位于边界自由面上的粒子出现无物理意义的飞散现象，导致时间步长逐渐减小，与理论结果的误差大于7.5%；当400≤N≤3 500时，与理论结果的误差小于4.5%；当N≥1 800时计算效率大大降低，且误差呈上升趋势，因此N的最佳取值范围为N＝600～1 800，在该范围内误差小于3%。

参考文献：

［1］Lucy L B.A numerical approach to the testing of the fission hypothesis［J］.Astronomical Journal，1977，82（12）：1013－1024.

［2］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）：375－389.

［3］Swegle J W，Attaway S W.On the feasibility of using smoothed particle hydrodynamics for underwater explosion calculations［J］.Computational Mechanics，1995，17（3）：151－168.

［4］Liu M B，Liu G R，Zong Z，et al.Computer simulation of high explosive explosion using smoothed particle hydrodynamics methodology［J］.Computers ＆ Fluids，2003，32（3）：305－322.

［5］Liu M B，Liu G R，Lam K Y，et al.Meshfree particle simulation of the detonation process for high explosives in shaped charge unlined cavity configurations［J］.Shock Waves，2003，12（6）：509－520.

［6］Liu G R，Liu M B.Smoothed particle hydrodynamics：A meshfree particle method［M］.Singapore：World Scientif－ic，2003.

［7］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.

［8］Petschek A G，Libersky L D.Cylindrical smoothed particle hydrodynamics［J］.Journal of Computational Physics，1993，109（1）：76－83.

［9］李裕春，时党勇，赵远.ANSYS11.0／LS－DYNA基础理论与工程实践［M］.北京：中国水利水电出版社，2008.

［10］Hallquist J O.LS－DYNA keyword user’s manual［M］.Livermore，CA：Livermore Software Technology Corporation，2007.

［11］隋树元，王树山.终点效应学［M］.北京：国防工业出版社，2000.

［12］Chou P C，Carleone J，Flis W J，et al.Improved formulas for velocity，acceleration，and projection angle of explosively driven liners［J］.Propellants，Explosives，Pyrotechnics，1983，8（6）：175－183.

［13］谭多望，孙承纬.大锥角罩聚能装药射流理论计算方法［J］.高压物理学报，2006，20（3）：270－275.

TAN Duo－wang，SUN Cheng－wei.Analytical model for jet formation in shaped charge with wide cone angle［J］.Chinese Journal of High Pressure Physics，2006，20（3）：270－275.

［14］Swegle J W，Hicks D L，Attaway S W.Smoothed particle hydrodynamics stability analysis［J］.Journal of Computational Physics，1995，116（1）：123－134.