柱状装药爆炸条件下厚壁圆筒爆室内流场及结构动力响应分析＊

姚哲芳，任辉启，沈兆武

（1.中国科学技术大学近代力学系，安徽 合肥 230027；2.西南科技大学土木工程与建筑学院，四川 绵阳 621010；3.总参工程兵科研三所，河南 洛阳 471023）

爆炸波模拟装置是模拟爆炸冲击波载荷环境的重要设备，对于防护工程与武器装备毁伤评估方面的研究，具有十分重要的意义。爆室为爆炸波模拟装置的重要部件，直接遭受炸药爆炸，受力环境十分恶劣。由于工程的实际需要，爆室设计为一端封闭一端开口的厚壁圆筒。为使厚壁圆筒能够在弹性范围内安全、重复地工作，需要深入研究它在柱状装药内爆炸条件下的强度问题。

对于上述问题，可参考火炮身管的设计经验以及爆炸容器的设计方法，而火炮发射药为火药，火药与高能炸药有较大区别。爆炸容器的研究始于1945年，中国于1963年建成了第一个爆炸容器，有关人员做了大量的研究工作［1－4］。尽管如此，由于爆炸容器设计涉及多个学科，动力学响应过程非常复杂，迄今为止还没有完整的爆炸容器设计标准。多数爆炸容器为两端封闭的柱形壳体，研究主要集中在两端封闭的柱形壳体（薄壁）在集团装药内爆炸载荷作用下的强度问题。而对一端封闭一端开口厚壁圆筒在柱状装药爆炸条件下强度问题的研究未见报道，应用于柱形壳体的一些研究方法（如等效静载法等），不再适用于上述问题。

由于理论分析存在较大困难，实验研究又受到场地条件、测试手段及经费等方面的限制，而计算机技术及数值算法有了很大发展，数值模拟成了上述问题研究的一个重要手段。本文中，主要采用数值模拟的方法，分析厚壁圆筒爆室内柱状装药爆炸非定常流场的演化过程以及筒体的动力响应规律，同时与实验结果及理论解进行对比，得出一些有价值的结论，为爆室的设计提供参考。

图1 厚壁圆筒爆室及装药剖面图Fig.1 Cutaway view of the blast chamber and the charge

1 有限元计算模型

厚壁圆筒爆室及装药结构如图1所示。圆筒内半径RI＝170mm，外半径RO＝340mm，筒体长度L1＝2 490mm，膛深L2＝2 198mm。TNT柱状装药直径68mm，长2 000mm，居中置放，左端位于筒底球心处，起爆点设置在爆室开口一端。

1.1 有限元模型及算法

由于爆室及装药结构是对称的，为了减少计算量，只取四分之一实体建模，在相关边界上施加对称边界条件。

有限元模型共有三部分，流体包括炸药和空气，结构为厚壁圆筒爆室。流体部分使用多物质ALE算法，建模时流体域取为足够大的规则六面体，全部以空气材料划分网格，网格尺寸为5mm×5mm×5mm，相关边界添加无反射边界条件，炸药材料利用在K文件中添加关键字＊INITIAL＿VOLUME＿FRACTION＿GEOMETRY方式在求解过程开始前自动填充到网格中。结构部分采用Lagrange算法，由于筒体球底结构的特殊性，划分网格时经过处理使筒体单元均为计算精度较高的八节点六面体单元，厚壁圆筒爆室有限元模型如图2所示。流固的耦合方式对模拟结果的准确性有重要影响，本文中采用基于罚函数的耦合算法完成流体与结构间的相互作用。有关算法可参见文献［5］。

1.2 测点布置

如图3所示，在靠近圆筒内壁处布置爆炸压力测点A～Y，以考察筒体内壁受到的爆炸载荷。为了全面考察筒体的动力响应，筒体上沿径向布置6排测点，建立原点位于筒底球面顶点的坐标系，6排测点的r坐标分别为0、RI／3、2RI／3、RI、（RI＋RO）／2和RO。筒体中部（与爆炸压力测点O 的z 坐标相等）三个测点编号为O1、O2和O3。

图2 厚壁圆筒爆室有限元模型Fig.2 FEM of the thick－walled cylinder

图3 测点布置Fig.3 Layout of the calculated points

1.3 材料模型及参数

TNT炸药密度1.63g／cm3，CJ爆速6 970m／s，爆压21GPa，采用JWL状态方程［5］，压力p 为比体积V与比内能E的函数

式中：A、B、R1、R2、ω为实验确定的常数。计算中，取A＝371GPa，B＝3.23GPa，R1＝4.15，R2＝0.95，ω＝0.3，V0＝1.0，E0＝7.0GPa。

空气密度1.29mg／cm3，采用简化的线性多项式状态方程［5］

式中：γ为空气的比热比，E为比内能，计算中取初始比内能0.25MPa。

厚壁圆筒材料为炮钢，密度7.85g／cm3，弹性模量200GPa，泊松比0.33。采用简化的Johnson－Cook材料模型及Gruneisen状态方程。Johnson－Cook材料模型［5］定义流动应力为

2 爆炸流场分析

计算中，TNT药量为11.833kg，图4给出了不同时刻爆炸流场的压力云图。可以清晰地看出：由于爆炸产物的侧向飞散，先导冲击波拖着斜激波向筒底运动，当t＝50μs时，稳定的爆轰波结构已经形成。爆轰波沿药柱向筒底传播，从药柱边缘发出的斜激波（激波角约35°）打在筒壁而发生规则反射，反射激波与筒壁间形成高压区。同时，高压的爆轰产物向筒外方向运动。当t＝210μs时，反射激波已在轴线处相遇形成了新的高压区，该高压区随着透射激波的运动而逐渐扩大。当t≈300μs时，柱状炸药爆轰完毕，先导冲击波透射入筒底的空气中，由于没有了爆炸反应区提供的能量，其后的压力迅速降低。接近t＝330μs时，先导冲击波碰到筒底发生反射，因为筒底是球面，反射激波汇聚形成t＝430μs时左端的高压区。但该高压气团压力低于筒壁面反射激波在筒轴线汇聚形成的高压气团，二者的运动速度相反。后者向侧向扩散的同时向筒底运动，形成了强大的高压气体“射流”，最终作用于筒底部，见t＝430～500μs时的压力云图，这一点可以从筒底处所受的压力时程曲线（见图5）看出。图5中，第一个峰值为先导冲击波作用于筒底产生的最大压力，第二个峰值为反射激波在轴线处汇聚后作用于筒底产生的最大压力。当t＝500μs后筒底的波系结构已变得十分复杂，难于分辨。

图6为圆筒内壁测点Y、T、O、J和E的压力时程曲线。从图上可以看出，靠近筒口的部位（测点Y）受到的爆炸载荷较小。由于稳定的爆轰波传播，斜激波打在筒内壁后发生反射，反射点匀速地从筒口扫到筒底，筒中部（测点T、O和J）的爆炸载荷峰值压力几乎一致，峰值压力在306MPa左右。由于装药爆轰完毕后斜激波的强度降低，所以靠近筒底的部位（测点E）的峰值压力小于筒中部的峰值压力。从测点O的压力时程曲线可以看出，在1000μs内，筒内爆炸波在该测点位置发生了4次反射，分别对应图中的4个波峰。

图4 爆炸流场压力云图Fig.4 Pressure contour of blast flow field

图5 筒底（测点A）压力时程曲线Fig.5 Pressure－time curve at the bottom of the cylinder

图6 圆筒内壁压力时程曲线Fig.6 Pressure－time curves on inner wall of the cylinder

3 筒体动力响应分析

图7 厚壁圆筒的von－Mises等效应力云图Fig.7 Von－Mises stress contour of the thick－walled cylinder

图8 圆筒内壁面上爆炸压力峰值及冲量的分布Fig.8 The distribution of pressure peak and impulse on the inner wall of the cylinder

厚壁圆筒爆室采用Johnson－Cook材料模型和断裂准则，在整个计算过程中，所有单元均没有失效，并且没有产生塑性应变，说明在该药量下厚壁圆筒爆室始终在弹性范围内工作。图7给出了不同时刻的厚壁圆筒爆室的von－Mises等效应力云图，所有单元的von－Mises等效应力最大值均小于或等于706MPa，而实际厚壁圆筒爆室所采用的材料许可应力在1000MPa以上，由第四强度理论知，厚壁圆筒爆室是安全的。

图9 测点O1的应力时程曲线Fig.9 Stress histories at point O1

图8给出了筒内壁爆炸压力峰值及冲量沿轴线的分布情况。从图上可以看出，自筒口向筒底一定距离范围内峰值压力逐渐增加，稳定爆轰形成后，筒内壁压力峰值相等，为306MPa。爆轰完毕后，筒壁压力峰值逐渐减小，由于壁面反射激波在轴线汇聚后形成高压，筒底球顶处受到的压力峰值最大，为411MPa。

计算中得到了筒体各测点的应力历史（主要考察等效应力），测点的应力变化范围应该可以反映整个筒体在动力响应过程中的应力历史。图9为筒体中部测点O1单元的应力历史，可以看出，由于筒壁体内应力波的传播很复杂，应力出现增长现象，即第一个应力峰值不是应力的最大值。

为了直观比较筒体各部位的等效应力峰值，图10给出了筒体上6排测点等效应力峰值的分布情况。可以看出，危险点在r＝0，z＝0，即筒底球面顶点，等效应力峰值达到706MPa。筒体材料许可应力取1000MPa，由第四强度理论得到爆室总体动态安全系数为1000／706≈1.4。筒底封堵部分三排测点（r＝0，RI／3，2RI／3）的等效应力水平相对整个筒体较高，均在550～706MPa范围内，说明筒底封堵部分是最危险的部位，实际使用时应适当保护，加入缓冲材料可以降低作用在筒底的峰值超压，降低该处的应力。另外，应考虑尽可能将柱状装药外移，增加装药与筒底的距离。

图10 筒体等效应力峰值的分布Fig.10 Von－Mises stress distribution of the cylinder

图11 测点O1、O2和O3环向应变时程曲线Fig.11 Hoop strain－time curves at points O1，O2and O3

图12 筒壁环向应变最大值的分布Fig.12 The distribution of hoop strain maximum

图13 测点O1环向应变的频谱Fig.13 Frequency spectrum of hoop strain at point O1

图11为筒中部测点O1、O2和 O3的环向应变波形，图12给出了筒壁环向应变最大值沿圆筒轴线的分布情况。对图11中的三个应变波形进行快速傅里叶变换分析（图13为内壁处测点O1应变波形的频谱），发现三者的频谱基本一致，径向运动的能量均主要集中在f＝3.516kHz附近的一个很小的区域内。

4 计算结果的验证

4.1 柱状装药的爆速

图4中的爆轰波结构类似于有限直径柱状装药的定常二维爆轰波结构，此时的斜激波运动速度与柱状装药的爆速相同。由激波关系式可知，爆速越高，斜激波波后的压力越高，打在筒内壁上形成的反射激波波后的压力也越高，即筒内壁受到的压力越高。另外，爆速实际上反映了炸药能量释放的快慢，对流场的非定常演化过程有重要影响。因此，为了准确地计算爆炸载荷与筒体的相互作用，必须准确地模拟柱状装药的爆速。柱状装药爆速D与直径d的经验关系式［7］为

式中：Di为理论爆速，A和df为拟合系数。由上式计算得到柱状装药爆速的经验值6942.4m／s。筒内壁压力测点等间距布设，由压力峰值的到达时间差可得到数值计算的柱状装药爆速6945.3m／s，与经验值吻合较好，说明数值计算结果是可信的。

4.2 计算结果与实测结果的比较

实验布置如图14所示，厚壁圆筒能够在支撑墙套管内滑动，封闭端顶有反后座装置防止后座行程过大，装药形式、尺寸及药量与数值模拟的情况一致。采用中北大学放入式电子测压器DT－FDCYQ对筒内底部靠近内壁面处的膛压进行测试，沿筒体轴线距离筒口不同位置上安放美国PCB公司137A系列传感器测试爆炸冲击波自由场压力。同时，在圆筒外壁面上距离开口端200mm处以及距离封闭端500mm处各对称粘贴一对应变片，测试环向应变。

图14 实验布置图Fig.14 Experimental layout

共进行了10炮次，发现筒体内外壁面无裂纹等损伤现象，内外径无变化，测得的数据重复性很好，说明在此工况下该厚壁圆筒爆室能够在弹性范围内可安全、重复地工作。筒内底部靠近内壁面处典型的膛压波形如图15（a）所示，图15（b）～（c）分别给出了沿筒体轴线距离筒口2.0和5.0m处的压力时程曲线。图16（a）～（b）分别给出了圆筒外壁面上距离开口端200mm处及距离封闭端500mm处环向应变，其中实测应变曲线经过了20kHz低通滤波处理。由图15～16可见，实测结果与数值计算结果吻合较好，证明本文中的有限元计算模型与方法是合理的。

图15 爆炸压力计算结果与实测结果的比较Fig.15 The comparison of blast pressure between simulation and experiment

图16 环向应变计算结果与实测结果的比较Fig.16 The comparison of hoop strain between simulation and experiment

4.3 轴对称平面应变假设条件下动力响应的理论解

在柱状装药内爆炸作用下，厚壁圆筒在弹性范围内工作，实际上是一个瞬态振动问题。圆筒长度比径向尺寸大得多，假设位于圆筒轴线上的柱状装药瞬间爆轰完毕，则筒壁所受的爆炸载荷沿长度方向不变化，问题简化为轴对称平面应变问题。

由弹性力学基本方程，易得厚壁圆筒径向振动微分方程为

式中：ρ、E和μ分别为材料的密度、弹性模量和泊松比。

问题的边界条件与初始条件分别为

文献［8］中给出了位移解u（r，t），再利用弹性方程求得径向应力σr（r，t）、环向应力σθ（r，t）和纵向应力σz（r，t），进而得到von－Mises等效应力σe（r，t）。

从图10可以看出，厚壁圆筒的管体部分（z／L2＞0）最大等效应力均发生在内壁（r＝RI），中部一段（0.3＜z／L2＜0.7）内壁的等效应力峰值基本相等，为570MPa左右。因此，将数值计算得到的圆筒中部测点O的爆炸压力曲线代入文献［8］中的动力响应理论公式，借助Maple数学计算软件，最终计算得到等效应力时程曲线，最大等效应力亦发生在内壁，为688MPa。可见，在轴对称平面应变假设前提下，厚壁圆筒动力响应的理论计算值偏大，可以作为保守估算。

由文献［8］中式（2.12）计算得到厚壁圆筒第一振型自由振动频率f＝3.540kHz。而由前面对应变波形的频谱分析，筒中部径向运动的能量主要集中在f＝3.516kHz附近的一个很小区域内。二者基本一致，说明在内柱状装药爆炸条件下，筒中部的径向运动因受端部效应的影响小以第一振形为主。

5 结 论

ANSYS／LS－DYNA有限元程序能够较好地模拟厚壁圆筒爆室在柱状装药爆炸作用下的动力响应过程，可以得到爆炸流场以及动力响应较全面的描述，弥补了实验的不足。本文中的研究方法可为工程提供参考。

在本文算例装药条件下，厚壁圆筒爆室在弹性范围内工作，由于先导冲击波以及筒壁面反射激波在轴线汇聚形成高压气团的作用，筒底球顶处的等效应力最大，为706MPa，总体动态安全系数为1.4，强度设计是安全的。

爆速实际上反映了炸药能量释放的快慢，对流场的非定常演化过程以及爆炸载荷与筒体的相互作用有重要影响。本文中数值计算得出柱状装药的爆速为6 945.3m／s，经验公式得出的为6 942.4m／s，二者吻合得较好。

筒体中部环向应变波形的频谱分析结果证明，此处的径向运动受端部效应影响小以第一振形为主。

轴对称平面应变假设下厚壁圆筒动力响应的理论解可以作为保守估算。

［1］Baker W E，Hu W C L，Jackson T R.Elastic response of thin spherical shells to axisymmetric blast loading［J］.Journal of Applied Mechanics，1966，33（4）：800－806.

［2］Duffey T A，Romero C.Strain growth in spherical explosive chambers subjected to internal blast loading［J］.International Journal of Impact Engineering，2003，28（9）：967－983.

［3］胡八一，柏劲松，张明，等.真实爆炸容器壳体动力响应的强度分析［J］.应用力学学报，2001，18（3）：136－139.HU Ba－yi，BAI Jin－song，ZHANG Ming，et al.The dynamic response analysis of a real explosion－container vessel［J］.Chinese Journal of Applied Mechanics，2001，18（3）：136－139.

［4］马圆圆，郑津洋，陈勇军，等.椭圆封头圆柱形爆炸容器动力响应的数值模拟［J］.爆炸与冲击，2009，29（3）：249－254.MA Yuan－yuan，ZHENG Jin－yang，CHEN Yong－jun，et al.Numerical simulation on dynamic response of the cylindrical explosion containment vessel with an elliptical cover［J］.Explosion and Shock Waves，2009，29（3）：249－254.

［5］Hallquist J O.LS－DYNA theory manual［Z］.Livermore：Livermore Software Technology Corporation，2006.

［6］Rusinek A，Rodríguez－Martínez J A，Arias A，et al.Influence of conical projectile diameter on perpendicular impact of thin steel plate［J］.Engineering Fracture Mechanics，2008，75（10）：2946－2967.

［7］孙承纬，卫玉章，周之奎.应用爆轰物理［M］.北京：国防工业出版社，2000.

［8］林祖森.受冲击内压作用的厚壁圆筒的动力分析［J］.兵工学报，1986，（1）：57－64.LIN Zu－sen.Dynamic analysis of thick cylinders subjected to internal impulsive pressure［J］.Acta Armamentarii，1986（1）：57－64.