石 磊,杜修力,樊 鑫

(北京工业大学建筑工程学院,北京 100124)

空气中爆炸冲击波峰值受空气和炸药的有限元网格尺寸影响,不同尺寸的网格会导致计算结果的巨大差异[1-6].网格尺寸按照 5mm,10mm……这样的确定数值进行划分,当网格尺寸不能被炸药体边长整除时,会带来一定的划分误差.当炸药的当量比较大时,会导致单元数量过大;针对某一具体当量炸药进行划分,当炸药当量不同时其使用存在一定的局限性;爆炸冲击波峰值比例距离越大,网格尺寸敏感性越差.

本文针对现存的问题,仔细研究了网格划分对爆炸冲击波峰值压力影响的规律,并且比较了计算精度和计算效率,给出了不同比例距离网格的优化划分方法.

1 炸药及空气单元划分方法

1.1 爆炸冲击波峰值经验公式

Josef[7]用实验的方法提出了无限域空气中炸药爆炸峰值压力经验计算公式

Baker[8]归纳了 TNT炸药包的 2个超压计算公式

式中,PSO为超压峰值(0.1MPa);Z为比例距离(m/kg1/3).

美国 TM5-1300手册[9]根据试验结果绘出了爆炸冲击波各参数之间的关系曲线以供设计使用.

1.2 空气及炸药状态方程及材料参数

炸药采用 LS-DYNA中的*MAT-HIGH-EXPLOSIVE-BURN材料模型[10],爆炸压力 P,单位体积内能及相对体积 V的关系采用 JWL状态方程[11]

式中,V为相对体积;E为单位体积炸药内能;A1,B1,R1,R2,ω为 JWL状态方程参数.JWL状态方程参数及炸药材料参数如表 1所示.

表1 相关计算参数表Tab le 1 Relative ca lcu lation parameters

空气采用 LS-DYNA中的*MAT-NULL材料模型[10]和线性多项式状态*EOS-LINEARPOLYNOMIAL来描述

其中,C0、C1、C2、C3、C4、C5、C6为用户定义的常量;E为初始能量;μ=1/V-1,V为初始相对体积.具体参数如表 2所示.

表2 相关计算参数表Tab le 2 Relative ca lcu lation parameters

1.3 单元网格尺寸对爆炸空气冲击波峰值的影响

分别取 1、20、40kg炸药作为计算模型,利用计算模型的对称性取1/8模型进行计算,对称面采用对称边界,非对称面采用透射边界[12].其单元划分分别采用 1/4,1/6,1/8,1/16倍的炸药体(整个炸药体,非1/8炸药体)边长作为网格尺寸,采用均匀划分的方法,如图 1所示.当比例距离小于 1.3时,网格尺寸对计算结果影响较大,这里计算并对比了比例距离在 0.1～1.3的空气冲击波峰值压力.

图1 1/8计算模型Fig.1 The one-eight symmetricalmodel

图2 (a)、(b)、(c)分别是 1、20、40 kg炸药 ,分别按照 1/4、1/6、1/8、1/16倍炸药体边长进行单元划分的峰值压力图,1)爆炸空气冲击波峰值压力随着单元尺寸的减小而越来越接近经验公式值;2)Henrych经验公式在比例距离 z为 0.4～0.7时,与 Bcker公式及 TM5-1300图表值相差较大,按以上单元尺寸划分的模型,在每一比例距离处都稍稍低于 Bcker公式及 TM 5-1300图表值,但是却在比例距离为 0.4～0.7稍稍大于 Henrych经验公式;3)在比例距离为 0.1～0.2时,数值模拟值比经验公式值大.

如图 3所示,对炸药体及 0.1～0.2比例距离内的空气体作进一步的划分,单元尺寸分别取 1/4,1/6,1/8,1/16,1/20,3/80倍炸药体边长.发现当进一步划分网格至 1/20,3/80倍炸药边长时,0.1～0.2比例距离的冲击波峰值压力很接近经验公式,没有了峰值压力明显大于经验公式值的情况,相比之下,按照 3/80倍炸药体边长进行网格划分效果更好.导致这种情况的主要原因可能是当单元比较大时,计算峰值提取点距爆心的比例距离总是小于 0.1,这导致了计算数值提取点与相应比例距离的偏差,而在 0.1～0.2比例距离内冲击波高频分量作用显著,受爆心距影响非常明显,这种微小的偏差足可以导致计算结果的巨大差异.当单元比较小时,这种偏差也就减小了,计算结果将更为精确.

图2 不同网格尺寸数值计算值与经验公式值对比Fig.2 Comparison between numerical results using differentmesh sizes and empirical predictions

图3 0.1～0.2比例距离内的对比Fig.3 Comparison in scaled-distance 0.1-0.2

1.4 不同当量炸药相同网格划分方法时冲击波峰值压力对比

当网格划分采用相同的炸药边长比例系数时,在各个比例距离上,不同当量的炸药,相同比例距离空气冲击波峰值非常接近,为此作了进一步的对比.图 4分别是 1、20、40 kg炸药,采用相同炸药边长比例系数进行网格划分时的峰值压力对比图,当采用相同的炸药边长比例系数时,不同当量的炸药在相同的比例距离下峰值压力基本相同.由此认为,当按照炸药边长的比例划分不同当量炸药单元及空气单元的时候,采用相同的比例系数,在比例距离相同时,会得到相同的爆炸冲击波峰值,本文的网格划分方法对于各种当量的炸药都有很强的适用性.

2 单元网格的渐进划分方法

计算模型还是 1、20、40 kg炸药.炸药单元的网格尺寸始终采用 3/80倍的炸药体边长,因为在比例距离为 0.2之内,爆炸荷载峰值对单元网格尺寸反应非常敏感,本文计算模型的渐进开始端分别取比例距离为 0.2处和炸药单元边处,相邻单元的渐进放大比例系数分别取 1.01,1.02,1.03进行网格划分.计算结果如图 5所示,当渐进比例系数为 1.01时计算值最接近经验公式,而为 1.03时计算误差最大,1.02介于两者之间.此外,渐进划分方法爆炸冲击波峰值的衰减规律更符合常理及经验,衰减曲线也更平滑,均匀.在比例距离为 0.1～0.6,无论采用那种渐进比例,都可以得到比较精确的峰值压力计算结果.但是在比例距离为 0.6以外,渐进比例大的计算误差也相应的增大.本文建议,在进行计算时,如果所关心的爆距在比例距离 0.1～0.6,可以采用 1.03倍的渐进划分方法;当所关心的爆距为比例距离在 0.6～1.4,衡量计算消耗与计算精度之后,建议采用渐进系数 1.02进行单元网格划分.

图4 不同当量炸药相同网格尺寸的数值计算值与经验公式值对比Fig.4 Comparisonbetweendifferentchargeweightsnumericalresultsusingthesamemeshsizes andempiricalpredictions

图5 渐进网格划分方法数值计算值与经验公式值对比Fig.5 Comparisonbetweennumericalresultsusingprogressivemeshmethodsandempiricalpredictions

由于渐进开始端的取法对单元数量及计算消耗有很大的影响,本文比较了不同渐进开始端的爆炸冲击波峰值压力.图 6是 1、20、40kg炸药分别以 1.01、1.02、1.03的渐进比例划分网格,起始端分别为炸药单元边处和比例距离为 0.2处的峰值压力对比图,当起始端为炸药单元边时,峰值压力在比例距离 0.1～ 0.2稍稍低于起始边为比例距离 0.2处的,这是因为,由于在比例距离 0.1～0.2内存在渐进划分,单元网格普遍大于 3/80倍的炸药单元边长.虽然峰值压力略微变小,但是仍然很贴近经验公式,所以计算结果仍然比较理想.而在比例距离 0.2～1.3,这两种网格划分方法无论哪个渐进比例系数,与经验公式相比都可以得到比较精确的计算结果.

图6 不同渐进起始点数值计算结果与经验公式值对比Fig.6 Comparisonbetweennumericalresultswithdifferentprogressivestartingpointandempiricalpredictions

取炸药体边长的 3/80作为炸药单元划分尺寸.当计算所关心的爆距在比例距离 0～0.2时,建议采用从比例距离为 0.2处开始渐进的划分方法.相邻单元间的渐进比例系数取 1.02;当计算所关心的爆距在比例距为 0.2以外时,从炸药单元边处开始渐进划分网格,在保证计算精度的前提下,可以减少计算单元数,在计算效率上占有很大的优势,所以建议采用从炸药单元边处开始渐近划分.此时当所关心的爆距在比例距离 0.1～0.6,渐进比例系数可取 1.03,当所关心的爆距为比例距离在 0.6～1.4,建议采用渐进系数 1.02进行单元网格划分.

3 结束语

通过与常用经验公式的比较,研究了不同网格尺寸对爆炸空气冲击波峰值的影响,建议了采用均匀划分方法时的网格尺寸与炸药体边长的比例系数.并对使用相同比例系数进行网格划分的不同当量炸药爆炸冲击波峰值进行了对比,结果表明对各种当量的炸药爆炸,此比例系数都有很好的适用性及较高准确性.为了进一步提高计算效率,合理的利用冲击波高频分量对其峰值的影响随着爆心距的变大逐渐变弱这一特性,本文发展了渐进的网格划分方法,并给出了渐进的比例系数.使用该方法,对不同当量炸药爆炸进行计算的结果表明:此方法在各种当量炸药爆炸情况中,不但具有理想的计算精度和更高的计算效率,而且具备更为广泛的适用性.

[1]GANTESC J,PNEVMATIKOSN G.Elastic-plastic response spectra for exponential b last loading[J].International Journal of Impact Engineering,2004,30:323-343.

[2]SHIYan-chao,LIZhong-xian,HAOHong.Mesh size effect in numerical simulation ofblastwave propagation and interaction with structures[J].Transactions of Tianjin University,2008,14(6):396-402.

[3]GONG Shun-feng,LU Yong,JINWei-liang.Simulation of airblast load and its effect on RC Structures[J].Transactions of Tian jin University,2006,12(Suppl):165-170.

[4]廖维张.常规武器爆炸作用下地下结构的动力响应及智能隔震研究[D].北京:北京工业大学建筑工程学院,2007.LIAO Wei-zhang.Dynamic response and intelligent vibration isolation of underground structures under explosion of cohventional weapons[D].Beijing:The College of Architecture and Civil Engineering of Beijing University of Technology,2007.(in Chinese)

[5]杨鑫,石少卿,程鹏飞.空气中 TNT爆炸冲击波超压峰值的预测及数值模型[J].爆破,2008,25(1):15-18.YANG Xin,SHIShao-qing,CHENG Peng-fei.Forecastand simulation of peak overpressure of TNT explosion shock wave in the air[J].Blasting,2008,25(1):15-18.(in Chinese)

[6]李顺波,东兆星,齐燕军,等.爆炸冲击波在不同介质中传播衰减规律及数值模拟[J].振动与冲击,2009,28(7):115-117.LIShun-bo,DONG Zhao-xing,QIYan-jun,et al.Numerical simu lation for spread decay of blasting shock ware in different media[J].Journal of Ibration and Shock,2009,28(7):115-117.(in Chinese)

[7]JOSEFH.爆炸动力学及其应用[M].熊建国,译.北京:科学出版社,1987:126-128.

[8]BAKERW E.Explosions in air[M].University of Texas Press,Austin,TX,1973:7-15.

[9]TM5-1300.Structures to resist the effects of accidental explosion[S].America:Department of the Army TechnicalManual,Department of the Navy Publication NAVFAC P-397,Department of the Air Force Manual AFM 88-22,Department of the Army,the Navy,and the Air Force,1969.

[10]LS-DYNA Keyword User's Manual(Version 950)[M].Livemore Software Technology Corporation,1999:20,33-20.37.

[11]白金泽.LS-DYNA3D理论基础与实例分析[M].北京:科学出版社,2005:74-75.

[12]时党勇,李裕春,张胜民.基于ANSYS/LS-DYNA 8.1进行显式动力分析[M].北京:清华大学出版社,2005:216-235.