1 引言
网格划分方法对空气爆炸的冲击波峰值数值模拟精度有较大影响，且不同炸药当量下其网格尺寸取值亦有较大区别[2-5]。爆炸冲击波峰值的数值精度随着网格尺寸的减小而增大，然而网格尺寸过小会导致单元数目呈指数倍增长，极大程度地降低了计算效率。因此，需要提出一种既综合考虑计算精度和效率，又对不同当量TNT均有较强适用性的网格划分方法，为大型结构采用流固耦合方法进行抗爆分析提供空气域网格划分方法的参考。
炸药在空气中爆炸后瞬间产生强烈的爆炸冲击波，是爆炸对周围目标产生破坏效应的重要因素之一，自20世纪40年代以来学者们便开展了广泛的研究[1]。柱状炸药作为主要的装药形式之一，其在空气中的传播规律吸引了众多学者的研究目光。索强等[2]利用试验数据，比较了不同网格划分比对爆炸冲击波数值结果的影响，结果表明当网格划分比为炸药体边长的1/32时，数值结果与试验结果的误差在5%以内；张社荣等[3]采用AUTODYN软件对不同网格划分比的水下爆炸传播规律进行了研究，结果表明采用炸药半径的1/3作为网格尺寸可满足工程精度要求，并给出了不同网格划分比下的误差拟合公式；石磊等[4]采用LS-DYNA软件研究了网格划分方法对爆炸冲击波的影响，结果表明网格尺寸的渐进比例和渐进比例起始位置对计算结果有较大影响；姚成宝等[5]利用LS-DYNA软件研究了无反射边界条件对冲击波参量的影响，结果表明无反射边界尺寸为考察区域2倍时数值结果与实验结果相似；陈铭等[6]采用LS-DYNA软件对圆柱装药的径高比和起爆方式进行了研究，结果表明径高比对比冲量影响较大；高轩能等[7]采用LS-DYNA软件分析了炸药材料参数、TNT药量、网格尺寸、空气域形状和炸药形状等对冲击波参量的影响；Luccioni等[8]利用流体力学软件研究了爆炸冲击波的网格尺寸效应，结果表明10 cm的网格尺寸便可较为精确的模拟冲击波的传播规律。
由此可见，学者们对爆炸冲击波的数值模拟研究较为全面，但是不同学者在研究中提出的网格划分方法存在较大差别，而网格划分方法对计算效率和计算精度的影响较大。本文采用LS-DYNA软件，对比分析不同炸药当量下网格划分方法对计算精度及效率的影响，提出了一种对不同当量均有较强适用性的网格划分方法。
2 基于经验公式的超压峰值
不同类型的炸药爆炸产生的能量可用等效当量TNT炸药来衡量。当爆炸发生在自由空气时，其传播规律由峰值超压、冲量、持续时间等参量来描述。目前，最常用的实验方法是采用Hopkinson比例定律[10]，以小当量爆炸来研究大当量爆炸的特性，比例距离Z的表达式如下：
Z=R/W1/3
(1)
式中，R为测点与爆心之间的距离，W为等效TNT药量。
目前，具有较高权威性的实验成果为UFC 3-340-02规范[9]中给出的设计图表；除此之外，其他学者以及我国的人防规范亦给出了相应的经验公式[11]，如图1所示，不同学者提出的经验公式结果存在较大差别，尤其是在靠近爆心差别越显著。图2所示，Henrych和baker提出的经验公式与UFC 3-340-02规范中提出的图表较为一致，当比例距离大于1时，各个公式预测的结果比较接近。文献[12]表明：当冲击波传播距离r>L时(r为测点距轴心的距离，L为柱形炸药的高度)，可近似看成球形装药的爆炸。为此，本研究采用Henrych公式对大当量柱状炸药的超压数值结果进行验证。
图1 不同经验公式下的冲击波超压峰值
Fig.1 Comparison of peak overpressure of shock wave
图2 Henrych、Baker与UFC 3-340-02超压峰值
Fig.2 Comparison of peak overpressure of Henrych，Baker and UFC 3-340-02
3 柱形炸药的数值计算及对比分析
3.1 材料的状态方程及参数
空气采用*MAT_NULL模型[7]，应用线性多项式状态方式*EOS_LINEAR_POLUNOMIAL描述，表达式如下：
(2)
μ=1/V
(3)
式中： C0、C1、C2、C3、C4、C5、C6 为常数，E为单位初始体积内能，V为初始相对体积。参数取值如表1所示。
表1 空气的材料参数
Table 1 Material parameters of air
C0C1C2C3C4-0.10000.4C5C6VE/(J·m-3)ρ/(kg·m-3)0.4012.51.29
炸药采用*MAT_HIGH_EXPLOSIVE_BURN材料模型[7]，定义压力P为相对体积V和初始能量E的函数，表达式如下：
(4)
式中: V为炸药的相对体积，E0为单位体积内能，A、B、R1、R2为表征炸药材料特性的常数，上式的三项由左至右分别代表峰值超压的高压段、中压段及低压段[13]。文献[14-16]采用的TNT炸药参数在CJ状态下满足约束守恒方程组，且能较好地模拟结构在爆炸工况下的动态响应过程，其参数取值如表2所示。
表2 TNT炸药的材料参数
Table 2 Material parameters of TNT explosive
A/GPaB/GPaR1R2ωV03713.2314.150.950.31.0E0/(GJ·m-3)ρ/(kg·m-3)D/(m·s-1)Pcj/GPa7.01 6306 93021.0
3.2 网格划分方法对数值结果的影响
文献[17]已验证了建立1/8模型进行爆炸冲击波研究的可行性。文献[5]通过定义相对厚度T=L/D，其中L为数值模型的最短边尺寸，D为柱状装药的半径，研究发现当T>2时，无反射边界对计算结果影响较小。文献[3]研究表明网格尺寸小于炸药半径3倍以上时，其水中爆炸数值结果满足工程精度要求。
本研究采用LS-DYNA软件建立TNT爆炸的数值计算模型。为了提高计算效率，本研究采用1/8长方体模型，单元类型取8节点SOLID164单元，空气域尺寸取为4 m×4 m×13m(高×宽×长)，通过关键字*CONSTRAINED_GLOBAL及*BOUNDARY_NON_REFLECTING定义模型的对称边界及无反射边界，材料及状态方程参数见表1和表2，以研究网格划分方法对爆点在原点、长径比为3的800kgTNT圆柱炸药的峰值超压的影响；通过关键字*DATABASE_ELOUT和*DATABASE_HISTORY_SOLID输出以爆炸中心为起点，沿空气域长度方向比例距离分别为0.1～1.1 m/kg1/3、时间间隔为1的峰值超压曲线，计算总时间为0.01 s。
3.2.1 网格划分比对数值结果的影响
网格划分比指“网格尺寸与炸药半径之比”。考虑爆炸近区的冲击波的衰减性较大，对比例距离小于0.3的爆炸区域进行加密(以爆炸中心为角点的3 m×3 m×3 m长方体区域)，加密区段网格尺寸分别为4 cm、5 cm、7.5 cm、10 cm,非加密区段网格尺寸均为10 cm[8]。
图3表明：当比例距离小于0.4时，数值结果与Henrych经验公式相差较大；图4表明：当比例距离大于0.4时，4 cm与5 cm网格尺寸计算结果与经验公式吻合度较好。不同网格尺寸下的峰值超压与Henrych经验公式的误差分析如表3所示，当比例距离大于0.4时，4 cm、5 cm、7.5 cm和10 cm网格尺寸的最大误差分别为10.1%、-12.5%、-31.8%和-31%，单元数目分别为1 264 375、735 000、325 000和208 000。结果表明：当比例距离大于0.4时，峰值超压的计算精度随着网格尺寸的减小而增加，4 cm与5 cm网格尺寸的数值结果精度相差不大，误差均在10%左右，但4 cm网格的单元数目5 cm网格的1.7倍，综合考虑计算效率和计算精度的关系，加密区采用5 cm的网格尺寸。
图3 网格尺寸对比例距离小于0.4的超压峰值曲线
Fig.3 Influence of grid size on peak overpressure with proportional distance less than 0.4
图4 网格尺寸对比例距离大于0.4的超压峰值曲线
Fig.4 Influence of grid size on peak overpressure with proportional distance greater than 0.4
表3 网格尺寸对超压峰值计算误差的影响(%)
Table 3 Influence of grid size on calculation error of overpressure peak(%)
网格尺寸/cm比例距离/m/kg1/30.20.30.40.50.60.70.80.91.01.14-49.5-33.4-17.40.66.09.09.88.310.18.05-57.3-40.3-35.4-10.7-12.5-10.5-7.7-5.4-3.2-0.77.5-72.4-61.6-44.7-31.8-23.9-21.3-18.3-16.1-8.8-15.310-59.3-66.8-44.0-31.0-25.4-20.9-18.3-18.1-16.5-24.9
3.2.2 加密范围对数值结果的影响
以3.2.1节的结果为依据，进一步探讨加密范围分别为4 m、5 m、6 m对峰值超压结果的影响。网格划分方法如下：以爆炸中心为角点，加密区域分别为4/5/6 m×4/5/6 m×4/5/6 m的长方体区域。若空气域某一方向小于该加密长度，取空气域该方向的最大长度作为加密长度，加密区网格尺寸取为5 cm，非加密区网格尺寸均为10 cm[8]。
如图5所示，图例中如1-0.05-3表示加密范围为3 m且加密范围为5 cm的等比例网格。图5表明：当比例距离小于0.4时，增大加密范围对数值结果与Henrych经验公式的误差影响不大；图6表明：当比例距离大于0.4时，增大加密范围，能在一定程度上提高数值结果的精度。如表4所示，当比例距离大于0.4时，峰值超压随着加密范围的增大而增大，当加密区长度由3 m增加至4 m时，最大误差由-12.5%减少至-4.4%，误差在5%以内，精度有了明显提高；加密区长度进一步由4 m增大至5 m和6 m时，对计算精度的提高影响不明显。由此可见，800 kg TNT当量下，加密区长度在4 m以内发生改变对精度有明显的影响，据此推断出其加密界限为0.4 m/kg1/3。
图5 加密范围对比例距离小于0.4的超压峰值曲线
Fig.5 Influence of encryption range on peak overpressure with proportional distance less than 0.4
图6 加密范围对比例距离大于0.4的超压峰值曲线
Fig.6 Influence of encryption range on peak overpressure with proportional distance greater than 0.4
表4 加密范围对超压峰值计算误差的影响(%)
Table 4 Influence of encryption range on calculation error of overpressure peak(%)
加密长度/m比例距离/(m·kg-1/3)0.20.30.40.50.60.70.80.91.01.13-57.3-40.3-35.4-10.7-12.5-10.5-7.7-5.4-3.2-0.74-57.3-40.3-11.4-4.41.23.54.02.33.7-0.35-57.3-40.3-11.47.02.06.06.87.06.12.26-57.3-40.3-11.47.07.36.37.16.98.75.3
3.2.3 渐进比例对数值结果的影响
以3.2.2节研究结果为依据，为了提高计算效率、节省计算时间，进一步探讨加密范围内渐进比例网格对峰值超压结果的影响。网格划分方法如下：加密范围取4 m，加密区网格尺寸取5 cm，加密范围网格尺寸分别以1.01、1.02、1.03的渐进比例系数沿x/y/z三个方向过渡到非加密范围，非加密网格尺寸取10 cm[8]。为了控制非加密区网格尺寸不超过10 cm和加密范围为4 m，渐进系数为1.01的渐进网格由爆心开始，渐进系数为1.02和1.03的渐进网格分别从距离爆心x/y/z三个方向的2.3 m和2.6 m开始，其单元数目分别为532 800、630 784、1 000 936。
如图7、图8所示，图例中如“1.01-0.05”代表加密区渐进比例系数为1.01，网格尺寸为5 cm的数值模型。图7表明当比例距离小于0.4时，渐进比例系数对近区冲击波超压影响不大；图8表明当比例距离大于0.4时，3种渐进比例系数都与经验公式有较好的吻合效果。表5展示了不同渐进比例系数下数值结果与Henruch的计算误差，可以看出，3种渐进比例系数下，当比例距离大于0.4时，误差都在10%左右，且当比例距离大于0.5时，误差在5%以内，计算精度较高。此外，当渐进系数为1.01时，单元数目较其他2种渐进比例少，约为等比例网格划分方式的0.5倍，在计算机处理器为Intel(R) Xeon(R) CPU E5-2620 v4 @ 2.10 GHz、内存为16 GB、LS-DYNA计算版本为ls971_s_R5.1.1等条件下，计算时间为23 min，较2.2.2节的加密区等比例网格划分方法缩短了约1/2计算时间。
图7 渐进比例对比例距离小于0.4的超压峰值曲线
Fig.7 Influence of progressive proportion on peak overpressure with proportional distance less than 0.4
图8 渐进比例对比例距离大于0.4的超压峰值曲线
Fig.8 Influence of progressive proportion on peak overpressure with proportional distance greater than 0.4
表5 渐进比例对超压峰值计算误差的影响(%)
Table 5 Influence of progressive proportion on calculation error of peak overpressure(%)
渐进比例比例距离/(m·kg-1/3)0.20.30.40.50.60.70.80.91.01.11-57.3-40.3-11.4-4.41.23.54.02.33.7-0.31.01-58.3-50.2-26.5-10.5-4.9-1.6-0.40.30.14.81.02-59.4-46.0-23.3-7.1-1.81.42.51.32.93.01.03-57.3-42.7-21.5-6.0-0.92.13.01.73.33.4
3.3 网格划分方法对不同当量TNT的适用性
网格划分方法的结果为：加密范围不小于0.4 m/kg1/3，加密区网格划分比不大于2/15、非加密区网格划分比不大于1/4或10 cm[8](取两者的较小值)，加密区渐进比例为1.01。根据以上研究成果，分别验证该网格划分方法对500 kg和1 500 kg当量，长径比为3的柱状TNT炸药的适用性。500 kg当量TNT和1 500 kg当量TNT的网格划分方法如下：
1) 500 kg当量。空气域尺寸为3.5 m×3.5 m×11 m(高×宽×长)，加密范围为3.5 m，加密区网格取4 cm，由爆心开始以1.01渐进比例沿x/y/z三个方向过渡到非加密区，非加密区网格尺寸均为8 cm。
2) 1 500 kg当量。空气域尺寸为4.5 m×4.5 m×15 m(高×宽×长)，加密范围为4.5 m，加密区网格取6 cm，以1.01渐进比例系数沿x/y/z三个方向过渡到非加密区，非加密区网格尺寸均为10 cm。
如图9所示，当比例距离小于0.4时，该网格划分方法所得的峰值超压与Henrych经验公式的误差较大；当比例距离大于0.4时，500 kg当量和1 500 kg当量TNT均与Henrych经验公式吻合度较好。由此验证了当比例距离大于0.4时，该网格划分方法对不同当量TNT均有较强的适用性。
图9 500 kg/1 500 kg当量的超压峰值曲线
Fig.9 Peak overpressure curve of 500 kg/1 500 kg equivalent
4 与文献提出的网格划分方法的对比
文献[2-4]提出的网格划分方法表现为：在全局范围内划分等比例网格，网格划分比分别为1/32、3/80和1/3。如表6所示，索强[2]和石磊[4]提出的网格划分方法误差在5%以内，精度较高，但在空气域尺寸为4 m×4 m×13 m(高×宽×长)的几何模型下，单元数目在千万级以上，计算资源要求高，容易导致在大当量(800 kg)爆炸工况下由于计算时间过长、内存较大等原因而难以进行结构的爆炸动态响应分析。张社荣[3]提出的网格划分方法单元数目少，但是计算精度较低。对比文献[2-4]提出的网格划分方法，本文所提出的网格划分方法在比例距离大于0.4 m/kg1/3时，误差在10%以内，虽然较文献[2，4]提出的方法精度低，适用的比例距离范围较小，但该方法在同等的空气域尺寸下，划分的单元数目低于1百万、计算内存占比低、计算效率较高，且在比例距离大于0.4 m/kg1/3时精度在10%以内，符合工程误差的标准，可以更好地应用于几何模型较大的结构在大当量TNT柱状炸药工况下的数值计算。
表6 本文的网格划分方法与文献[2-4]方法的精度及效率对比
Table 6 Comparison of the accuracy and efficiency between the meshing method in this paper and the method in literature[2-4]
爆炸分类炸药形状适用范围/(m·kg-1/3)单元数目/个炸药当量/kg计算误差/%文献[2]空气爆炸长方体0.1～1.350 528 448<1 000<5文献[4]空气爆炸长方体0.1～1.430 459 375<40<5文献[3]水中爆炸圆柱体0.1～1.6133 947<1 000<35本文方法空气爆炸圆柱体>0.4532 800<1 500<10
注：误差的计算以Henrych经验公式为基准；以4 m×4 m×13 m的空气域尺寸进行单元数目计算
5 结论
本文使用LS-DYNA软件，对大当量柱状TNT炸药在自由空气中的冲击波传播过程进行了数值模拟，通过分析网格划分比、网格加密范围和爆炸近区渐进系数对不同当量下长径比为3的TNT柱状炸药的计算精度及效率的影响，提出了一种“加密区网格划分比不大于2/15，非加密区网格划分比不大于1/4或10 cm(取两者的较小值)，加密范围不小于0.4 m/kg1/3，加密区渐进比例系数为1.01”的网格划分方法，研究结果表明：
1)院该网格划分方法对不同当量长径比为3的TNT柱状炸药均有较强的适用性，当比例距离大于0.4时，与Henrych经验公式相比，超压峰值误差均在10%以内。
2) 该网格划分方法单元数目较少、计算内存占比低、计算效率较高，当比例距离大于0.4 m/kg1/3时，精度在10%以内，适用于几何模型较大的结构在大当量TNT柱状炸药工况下的数值计算。
[1] KREHL P O K.History of shock waves，explosions and impact:a chronological and biographical reference[M].Springer Science and Business Media，2008:1-9.
[2] 索强，徐鹏，尤文斌.网格划分对冲击波波形的影响分析[J].兵器装备工程学报，2020，41(02):198-203.
Suo Q，Xu P，You W B.Analysis of the influence of meshing on shock wave waveform[J].Journal of Ordnance Engineering，2020，41(02):198-203.
[3] 张社荣，李宏壁，王高辉，等.水下冲击波数值模拟的网格尺寸确定方法[J].振动与冲击，2015，34(08):93-100.
Zhang S R，Li H B，Wang G H，et al.Mesh size determination method for underwater shock wave numerical simulation[J].Vibration and Shock，2015，34(08):93-100.
[4] 石磊，杜修力，樊鑫，等.爆炸冲击波数值计算网格划分方法研究[J].北京工业大学学报，2010，36(11):1465-1470.
Shi L，Du X L，Fan X，et al.Research on meshing method for blast wave numerical calculation[J].Journal of Beijing University of Technology，2010，36(11):1465-1470.
[5] 姚成宝，王宏亮，张柏华，等.TNT空中爆炸冲击波传播数值模拟及数值影响因素分析[J].现代应用物理，2014，5(01):39-44.
Yao C B，Wang H L，Zhang B H，et al.Numerical Simulation of Shock Wave Propagation in TNT Air Explosion and Analysis of Numerical Influential Factors[J].Modern Applied Physics，2014，5(01):39-44.
[6] 陈铭，周云波，孙晓旺.TNT径高比对爆炸冲击波比冲量的影响研究[J].兵器装备工程学报，2018，39(06):62-66.
Chen M，Zhou Y B，Sun X W.Influence of diameter to height ratio of TNT on specific impulse of blast wave[J].Journal of Ordnance Engineering，2018，39(06):62-66.
[7] 高轩能，吴彦捷.TNT爆炸数值计算及其影响因素[J].火炸药学报，2015，38(03):32-39.
Gao X N，Wu Y J.Numerical calculation of TNT explosion and its influencing factors[J].Journal of Explosives，2015，38(03):32-39.
[8] Luccioni B，Ambrosini D，Danesi R.Blast load assessment using hydrocodes[J].Engineering Structures，2006，28(12):1736-1744.
[9] UFC 3-340-02.Structures to resist the effects of accidental explosions[S].Washington DC，USA:USA Army Corps of Engineering，2008.
[10] Hopkinson B.British ordnance board minutes[R].British Ordnance Office，London，UK，1915.
[11] 李忠献，任其武，师燕超，等.重要建筑结构抗恐怖爆炸设计爆炸荷载取值探讨[J].建筑结构学报，2016，37(03):51-58.
Li Z X，Ren Q W，Shi Y C，et al.Discussion on the value of explosion load in anti-terrorist explosion design of important building structures[J].Journal of Building Structures，2016，37(03):51-58.
[12] 周旭.导弹毁伤效能试验与评估[M].北京:国防工业出版社，2014.
Zhou X.Missile damage effectiveness test and evaluation[M].Beijing:National Defense Industry Press，2014.
[13] 胡兆颖，唐德高.空气中TNT爆炸的数值模拟[J].爆破，2014，31(04):41-45.
Hu Z Y，Tang D G.Numerical simulation of TNT explosion in air[J].Blasting，2014，31(04):41-45.
[14] 汪维，刘光昆，刘瑞朝，等.近爆作用下方形板表面爆炸载荷分布函数研究[J].中国科学(物理学·力学·天文学)，2020，50(02):144-152.
Wang W，Liu G K，Liu R Z，et al.Study on the distribution function of blast load on square plate surface under near detonation[J].Science in China(Physics，Mechanics，Astronomy)，2020，50(02):144-152.
[15] Yan B，Liu F，Song D Y，et al.Numerical study on damage mechanism of RC beams under close-in blast loading[J].Engineering Failure Analysis，2015，51(11):9-19.
[16] 杨赞，韩国振，严波，等.内爆荷载作用下PC箱梁桥的动态响应过程[J].高压物理学报，2021，35(01):1-11.
Yang Z，Han G Z，Yan B，et al.Dynamic response process of PC box girder bridge under implosion load[J].Journal of High Pressure Physics，2021，35(01):1-11.
[17] 李晓勇，崔村燕，陈景鹏，等.LS-DYNA软件开展爆炸冲击波计算时需考虑的问题[J].装备学院学报，2014，25(04):79-84.
Li X Y，Cui C Y，Chen J P，et al.Problems to be considered when LS-DYNA software is used to calculate blast wave[J].Journal of Equipment Institute，2014，25(04):79-84.