应用光滑粒子流体动力学方法仿真引信撞击雨滴过程

刘　鹏，王雨时，闻　泉，张志彪

(南京理工大学机械工程学院，江苏 南京 210094)



应用光滑粒子流体动力学方法仿真引信撞击雨滴过程

刘鹏，王雨时，闻泉，张志彪

(南京理工大学机械工程学院，江苏 南京 210094)

摘要:针对引信防雨机构优化设计需要清晰了解高速撞击雨滴过程的问题，应用光滑粒子流体动力学方法，建立了其与有限元方法耦合数值模型，仿真了弹头引信高速撞击球形雨滴过程。仿真结果与已有试验结果和理论较为吻合，说明该数值分析方法可行、准确。以美国M594引信头部外凸刚性剪切式防雨触发机构为例，仿真了引信高速撞击雨滴过程并得到了相应的机构动态特性。

关键词:弹头触发引信；光滑粒子流体动力学；动态特性；防雨机构；固液碰撞；雨滴

0引言

为了确保引信全天候作战性能，研制时需考核防雨性能[1-2]。因自然雨场受气候影响很大，故用自然雨场进行引信防雨试验是不切实际的。人工雨场因成本昂贵、试验复杂，未能推广。在防雨性能考核时常采用等效靶板[3]替代雨滴，但理论结果与试验结果有较大差异，其准确性值得商榷。引信撞击雨滴速度极大、作用时间极短，高速录像难以捕捉，用人工雨场试验时，无法观测碰撞过程。固液撞击是一个十分复杂的过程，液体具有部分自由边界，与固体界面没有固定的接触面积，撞击过程中还伴随着液体飞溅，属于自由边界流体动力学问题。同时撞击过程还显现出高度非线性特征[4-5]，这些都为该问题数值模拟带来困难，因此数值分析结果与试验结果往往会存在较大差异。随着计算机技术和数值计算方法的发展，使采用数值手段分析液固高速撞击过程中流固耦合动力学响应过程成为可能。黄晓毛等人[6]曾利用反弹道原理进行了引信碰雨仿真试验，用高速模拟雨滴碰撞M739引信防雨杆，测得防雨杆的变形情况，对M739引信受雨滴冲击的动态特性有了初步认识；谢永慧等人[7]建立了光滑粒子流体动力学方法与有限元方法耦合模型，通过数值模拟获得了材料在高速液体撞击下的表面损伤过程，所得结果与试验结果较为吻合；王安文等人[8]对液滴分别采用光滑粒子流体动力学和任意拉格朗日-欧拉两种方法模拟了液滴高速冲击有机玻璃过程中液滴变形和固体结构损伤。结果表明：两种算法中结构表面损伤情况与Brunton实验结果吻合，证明了数值分析方法的可行性和精确性，还发现了光滑粒子流体动力学方法模拟液固撞击过程更具有优势。本文建立光滑粒子流体动力学方法与有限元方法耦合数值模型仿真引信高速撞击雨滴过程，并结合已有的理论和试验结果，与仿真结果进行对比分析。

1光滑粒子流体动力学仿真方法介绍

光滑粒子流体动力学方法即SPH是(Smooth Particle Hydrodynamic)方法，该方法由Lucy和Gingold等人于1977年提出，最初是应用于天文物理。由于光滑粒子流体动力学方法不用网格，没有网格畸变问题，所以能在拉格朗日格式下处理大变形问题。光滑粒子流体动力学方法已成功地应用于水下爆炸模拟、高速碰撞等材料动态响应的数值模拟等领域，可以模拟连续体结构的解体、破碎、成坑以及固体层裂和脆性断裂等现象。光滑粒子流体动力学粒子的离散特性，使得对液体飞溅的模拟成为可能。相对于液体的大变形，固体材料在撞击过程中的变形量是比较小的，利用有限元方法配合适当的材料模型便可以模拟。

光滑粒子流体动力学方法的基础是插值理论。在光滑粒子流体动力学方法中任一宏观变量(如密度、压力和温度等)都能方便地借助于一组无序点上的值表示成积分插值计算得到。各质点的相互作用借助于插值函数来描述。利用插值函数给出量场在一点处的核心估计值，将连续介质动力学的守恒定律由微分方程形式转换成积分形式，进而转换为求和。

光滑粒子流体动力学方法中，质点近似函数定义为：

(1)

式(1)中，W是核函数(插值核)。

核函数W使用辅助函数θ定义：

(2)

式(2)中，d是空间维数；h是光滑长度，随时间和空间变化。

W(x,h)是尖峰函数，光滑粒子流体动力学方法中最常用的光滑核是三次B样条，定义为：

(3)

式(3)中，C是由空间维度确定的标准化常量，由归一化条件确定。在二维、三维情况下C分别等于5/(7πh2)、1/(πh3)，x为空间点与所求场函数值点之间的距离。

2仿真模型验证

文献[9]中给出了球形雨滴与静止平板碰撞产生的压力(单位面积上)为βρvc。这里ρ为雨滴密度、v为碰击速度，而c则是雨滴中的声速，β为系数，当速度v大于1 524 m/s即5 000 ft/s时，β近似等于0.5，而当速度v小于305 m/s即1 000 ft/s时，β近似等于0.7。为对照文献[3]中的射击试验，取雨滴直径为5.6 mm，初速为645 m/s，插值估取β等于0.644，此时P=βρvc=0.644×1 000×645×1 480=0.615 GPa。

利用LS-DYNA软件对雨滴撞击平板进行仿真建模，如图1所示。平板材料为30CrMnSiA钢，采用Johnson-Cook材料模型，平板圆柱面和底面设置为无反射边界；雨滴为SPH粒子，采用NULL材料模型，对材料模型赋予Grueisen状态方程；仿真参数如表1和表2所列。

图1　雨滴撞击力仿真模型Fig.1　The simulation model of raindrop impact force

密度ρ/(g/cm3)剪切模量G/GPa弹性模量E/GPa泊松比μ屈服应力A/MPa应变硬化模量B/MPa7.83772060.30792510应变硬化指数N应变率相关系数C温度相关系数M熔化温度Tmelt/K比热Cp/(J/(kg·K))失效应力Pc/(g·cm-1·μs-2)0.260.0141.031793477-9

表2　水模型参数[11]

图2给出了雨滴撞击平板时的压力分布云图；图3给出了平板中心附近点的压力变化，A到G为图2中沿径向依次远离中心的单元。撞击中心处的A单元压力最大，B单元到F单元的最大压力沿径向依次减小并最大值约为0.648 5 GPa，与上述理论计算结果0.615 GPa已很接近，误差约为5.4%。

图2　平板表面压力分布云图Fig.2　The pressure contours on plate surface

图3　平板中心附近径向各单元表面压力Fig.3　The surface pressure of each near plate center unit

文献[3]借鉴铜柱测压法在人工雨场中进行射击试验来测量雨滴对引信的作用力。为进一步与试验结果对比，建立有限元简化模型如图4所示，活塞直径为11.2 mm，材料为45钢，采用Rigid材料模型，铜柱为标准锥形铜柱，规格为Φ6×9.8，材料为高导无氧铜(OFHC)，各材料参数取自文献[12]。根据文献[3]给出的射击试验结果3 043 N，先对铜柱进行预压，预压力大小为3 000 N，持续时间为150 μs，预压结束后再使雨滴撞击活塞，仿真结果如图5所示。

图4　铜柱测压仿真模型Fig.4　The simulation model of pressure-measuring

图5　铜柱高度变化Fig.5　The height variation of copper cylinders with copper cylinders

由图5可以看出铜柱在预压力作用下高度降低到7.75 mm，随后在雨滴撞击力作用下，铜柱高度变为7.71 mm。对照标准铜柱压力-变形表可以查得铜柱高度7.75 mm时对应压力3 040 N，铜柱高度7.71 mm时对应压力为3 138 N，则仿真误差为：

3引信高速撞击雨滴过程仿真

以美国M594弹头引信和德国DM301结构为例，利用上述方法对引信高速撞击雨滴过程进行仿真。头部触杆直径为5.5 mm，凸起部直径为7 mm，凸缘厚度0.5 mm，外凸长度为3.5 mm，如图6所示。作为对比，图6另给出头部触杆内凹结构。利用LS-DYNA建立有限元模型，取雨滴直径为5.6 mm，初速为1 067 m/s，引信体材料取45钢，触杆的材料分别取铝合金2A12、45钢，其主要仿真参数如表3所列。通过比较圆台处的最大剪切应力来分析雨滴对引信的撞击作用。

图6　引信头部结构与尺寸Fig.6　The structure and dimension of Fuze’s head

材料密度ρ/(g/cm3)剪切模量G/GPa泊松比μ屈服应力A/MPa应变硬化模量B/MPa应变硬化指数N应变率相关系数C温度相关系数M熔化温度Tmelt/K比热Cp/(J/(kg·K))2A12铝合金2.8260.334004240.350.0011.42686392145钢7.852000.325073200.280.0641.061793440

图7和图8为截锥结构触杆采用2A12材料时，外凸和内凹放置两种方式的仿真效果图。比较图7与图8可以明显看出，雨滴在撞击过程中存在明显的飞溅现象，触杆外凸放置时，由于周围是空气没有约束，所以雨滴会迅速飞溅开来，触杆虽然有大变形但未剪断；而当触杆内凹放置时，由于引信内壁的束缚，所以雨滴会集聚成束作用在触杆上，触杆被剪断。

图7　外凸截锥结构触杆采用2A12材料的仿真结果Fig.7　The simulation results of outer truncated cone structure trigger rod with 2A12 material

图8　内凹截锥结构触杆采用2A12材料的仿真结果Fig.8　The simulation results of inset truncated cone structure trigger rod with 2A12 material

对比采用45钢材料时，截锥形触杆两种放置方式圆台处的最大剪切应力，其仿真结果如图9所示。图中曲线起点不同是由于雨滴与触杆之间的距离差异造成的。从图9可以明显看出，内凹放置时触杆受到的作用力大于外凸放置，且作用持续时间更长。

图9　触杆为45钢时其凸缘处最大剪切力仿真结果Fig.9　The simulation results of 45 steel trigger rod flange’s maximum shear stress

4雨滴撞击力仿真

在LS-DYNA中提供方便、快捷的撞击力生成命令语句*DATABASE_RCFORC，其作用是记录各个单元上的接触力，同时自动求和得到撞击力的合力。在关键字文件中输入相应参数，通过ASC II选项卡即可提取撞击力时程曲线[15]。为了分析不同直径雨滴的撞击力大小给引信防雨设计提供参考，对三种直径的球形雨滴以不同速度撞击平板进行仿真，结果如表4所示，表中撞击持续时间为F≥50 N的时间段。表4给出直径Φ5.6 mm雨滴以撞击速度645 m/s撞击无限大平面时的最大撞击力为12 941 N,与理论计算结果14 307 N[3]接近，但是撞击面积接近雨滴最大横截面积条件下射击试验锥形铜柱测试结果3 043 N的4.25倍，估计此差异与机构响应特性即机械测试原理有关。雨滴撞击力随时间变化的曲线如图10所示，由图10可以看出雨滴撞击时的撞击时间极短，机械原理感知测试难以响应，因而理论得出的峰值测试不到。测试得到的撞击力3 043 N对应的脉冲时间宽度约为8 μs。引信防雨机构设计时建议按表4给出最大撞击力的四分之一估取后评价引信机械零件的碰雨强度。具体的引信防雨机构设计，应尽可能采用相应的仿真评价其防雨特性。

图10　直径Φ5.6 mm的球形雨滴以645 m/s速度撞击平板时的撞击力Fig.10　The force of 5.6 mm diameter spherical raindrop impact plate with 645 m/s speed

撞击速度v/(m/s)雨滴直径Φ2.5mm雨滴直径Φ4mm雨滴直径Φ5.6mm最大撞击力Fmax/N撞击持续时间t/μs最大撞击力Fmax/N撞击持续时间t/μs最大撞击力Fmax/N撞击持续时间t/μs20023811.269023.5149839.63005239.5141719.4277631.64009438.9248516.1498226.750014638.4395914.6769423.660019807.1568213.71124921.764522446.9667112.91294121.570026206.8761912.41509021.180033406.3945612.21941020.590041445.91212911.72354519.7100049155.81448111.32907819.3107656185.61552411.23292418.0110058995.41686111.13484617.2120070075.31920611.04195816.8130078695.22215310.54791916.2140091714.9253959.45462215.61500103094.6286499.26137114.8

5结论

通过仿真结果看，用光滑粒子流体动力学方法仿真研究弹头触发引信防雨性能是可行的。该方法可以模拟雨滴与引信高速碰撞过程，并能得到试验时难以观测到的多种数据。例如，当球形雨滴直径为5.6 mm、引信外凸触杆直径为5.5 mm、弹丸速度为1 076 m/s时，作用在引信头部触杆剪切凸缘上最大剪切应力为615 MPa(外凸头部结构)至701.6 MPa(内凹头部结构)，约相当于雨滴最大冲击力为5 313 N(外凸头部结构)至6 061 N(内凹头部结构)，峰值持续时间为微秒级约为5～10 μs，且呈振荡规律衰减。对于其他头部结构的引信的雨滴撞击后特性，可参照此方法作详细仿真研究。

参考文献:

[1]张德智，王雨时，曹纯柱. 国外触发引信防雨装置研究和防雨性能现状[J]. 沈阳工业学院学报，1989，8(2)：62-73.

[2]王雨时，张德智，曹纯柱. 触发引信防雨装置研究进展[J]. 现代引信，1991，9(1)：23-31.

[3]王雨时，张德智，曹纯柱. 触发引信防雨试验的等效靶板模拟[J]. 沈阳工业学院学报，1994，13(1)：27-37.

[4]孙弼，鄢宇鹏，张荻，等. 液滴与可变形固面撞击的二维非线性激波模型[J]. 应用力学学报，1996，13(3)：33-38.

[5]张荻，周屈兰，谢永慧，等. 液固撞击的非线性波动模型的研究[J]. 西安交通大学学报，2002，36(11)：1138-1142.

[6]黄晓毛，程翔，王炅，等. M739引信防雨杆受雨滴冲击作用力实验研究[J]. 弹道学报，1999，11(4)：930-933.

[7]谢永慧，汪勇，陈建辉，等. 固体材料在高速液体撞击下的表面损伤研究[J]. 兵工学报，2009，30(S2)：139-144.

[8]王安文，徐绯，张岳青. SPH方法在液固撞击数值模拟中的应用[J]. 计算物理，2012，29(4)：525-532.

[9]Mullinix B R, Ranson W F, Swinson W F. The Impact Pressure Distribution of a Hypervelocity Raindrop[R]. US: ADA,1975..

[10]杨翔. 某航空火箭弹弹头触发引信相关技术研究——弹道环境与触发机构动态特性[D]. 南京：南京理工大学，2014.

[11]林晓东，卢义玉，汤积仁，等. 基于SPH-FEM耦合算法的磨料水射流破岩数值模拟[J]. 振动与冲击，2014， 33(18)：170-176.

[12]徐绍帅. 某枪榴弹弹底机械触发引信相关技术研究[D]. 南京：南京理工大学，2013.

[13]张伟，肖新科，魏刚. 2A12铝合金本构关系和失效模型[J].爆炸与冲击，2013，34(3)：276-282.

[14]陈刚，陈忠富，徐伟芳，等. 45钢的J-C损伤失效参量研究[J].爆炸与冲击，2007，27(2)：131-135.

[15] 刘思明，崔堃鹏，夏禾，等.基于LS-DYNA的车辆对桥墩撞击力仿真分析[J].铁道标准设计，2013(8)：70-74.

Approximate Analytical Expression of Supersonic Projectile Low Trajectory

LIU Peng, WANG Yushi, WEN Quan, ZHANG Zhibiao

(College of Mechanical Engineering, NUST, Nanjing 210094, China)

Abstract:Aiming at the problem of low trajectory parameters of supersonic projectile for fuze and ammunition test solution was complicated, rational empirical formulas for air resistance law of projectile at supersonic section was uesed to solve mass center motion equations of low trajectory. A analytic function of speed and time of low trajectory for supersonic projectile was present. The analytic function was simple to be programmed and convenient to use. Exam results showed that the approximate analytic expression of low trajectory for supersonic projectile had higher accuracy.

Key words:external ballistics; mathematical model; empirical formula; low trajectory; ballistic approximate solution; ballistic forecast

中图分类号:TJ430.1

文献标志码:A

文章编号：1008-1194(2016)01-0018-06

作者简介:刘鹏(1991—)，男，湖南长沙人，硕士研究生，研究方向:引信及弹药技术。E-mail：skiy906@163.com。

*收稿日期:2015-07-21