AUTODYN水下爆炸数值模拟研究.pdf

第2 6 卷第3 期 爆破 V 0 1 ．2 6N 。．3 2 0 0 9 年9 月 B L A S T I N G S e p ．2 0 0 9 - l l _ l _ o _ _ o I l _ - - _ l _ l _ - _ l _ _ I - l - _ _ l _ _ _ _ _ l - - l I _ _ _ _ _ - _ I _ - I l _ l _ l l - _ l I - _ I _ _ - _ _ _ _ _ 文章编号1 0 0 1 4 8 7 x 2 0 0 9 1 0 3 0 0 1 8 0 4 A U T O D Y N 水下爆炸数值模拟研究 刘科种1 ，徐更光1 ，辛春亮2 ，杨拯磊1 ，秦建1 ’3 1 ．北京理工大学爆炸科学与技术国家重点实验室，北京1 0 0 0 8 1 ；2 ．中国运载火箭技术研究院，北京1 0 0 0 7 6 ； 3 ．海军装备研究院，北京1 0 0 0 7 3 摘要采用A U T O D Y N 数值计算软件，对1 N T 水下爆炸冲击波传播过程进行了数值模拟。计算结果表 明，爆炸初期炸药中的爆轰波和水中冲击波是高频波，炸药及其附近水域需要足够细密的网格才能反映出这 种特性，建议按照峰压衰减时间常数内至少采样1 0 0 次的标准划分网格；同时应采用较小的人工粘性系数， 减少对冲击波峰值超压的抹平，并对计算结果进行低通滤波，这样计算得到的冲击波压力曲线与经验公式计 算曲线吻合得较好，不但可以准确捕捉冲击波峰值压力，而且能够比较准确地计算冲击波在水中传榆时的衰 减演化过程。 关键词爆炸力学；数值模拟；网格密度；人工粘性 中图分类号0 3 8 9文献标识码A R e s e a r c ho nN u m e r i c a lS i m u l a t i o ni nU n d e r w a t e r E x p l o s i o nb yA U T O D Y N L I UK e ．z h o n 9 1 ，X UG e n g ．g u a n 9 1 ，X I NC h u n ．1 i a n 9 2 ，Y A N GZ h e n g l e i l ，口，ⅣJ i a n L 3 1 ．S t a t eK e yL a b o r a t o r yo fE x p l o s i o nS c i e n c ea n dT e c h n o l o g y ，B e i j i n gI n s t i t u t eo fT e c h n o l o g y ， B e i j i n g1 0 0 0 8 1 ，C h i n a ；2 ．C h i n aA c a d e m yo fL a u n c hV e h i c l eT e c h n o l o g y ，B e i j i n g1 0 0 0 7 6 ，C h i n a ； 3 ．N a v a lA c a d e m yo fA r m a m e n t ，B e i j i n g1 0 0 0 7 3 ，C h i n a A b s t r a c t T h ep r o p a g a t i o no fu n d e r w a t e rs h o c kw a v eo fT N Tw a ss i m u l a t e db yA U T O D Y N ．T h er e s u l t ss h o wt h a t t h ed e t o n a t i o nw a v ei ne a r l ye x p l o s i o na n dw a t e rs h o c kw e r eh i g h f r e q u e n c yw a v e ，t h i sp h e n o m e n o nn e e d e df i n em e - s h e si ne x p l o s i v ea n dn e a r b yw a t e ra r e a ，W ea d v i s e dm e s h i n gb a s e do nt h es a m p l i n gf r e q u e n c yl a r g e rt h a n10 0d u r i n g d a m p i n gt i m eo fp e a kp r e s s u r e ，鹪w e l la f ts m a l lp s e u d o - v i s c o s i t yc o e f f i c i e n t s ，t or e d u c et h et o w e l i n go fp e a kp r e s - s u r e ．a n dl o w - p a s s i n gt h ec a l c u l a t i o n ．T h ep r e s s u r ec u r v e so fs h o c kw a v eb ys i m u l a t i o na r ea p p r o x i m a t e l yi na c c o r d - a n e ew i t ht h ee m p i r i c a le q u a t i o n ，a n dt h es i m u l a t i o nn o to n l yC a na c c u r a t e l yc a p t u r ep e a kp r e s s u r eo fs h o c kw a v e ，b u t a l s ob e t t e rc a l c u l a t et h ep r o p a g a t i o na n da t t e n u a t i o ni nw a t e r ． K e yw o r d s e x p l o s i o nm e c h a n i c s ；n u m e r i c a ls i m u l a t i o n ；m e s hd e n s i t y ；p 北u d o v i s c o s i t y 0引言 炸药水下爆炸产生的冲击波是对目标产生破坏 的重要因素，在进行舰船结构的抗爆设计和分析时， 水中冲击波作为输入载荷，其计算精度直接影响到 收稿日期2 0 0 9 0 6 2 6 作者简介刘科种 1 9 8 2 一 ．男，博士研究生，主要从事爆炸力学数 值计算及炸药应用研究，E - m a i l l i u k e z h o n g b i t ．e d u ．c n 。 整体分析的精度，因此有必要建立炸药水中爆炸冲 击波较精确的计算方法。张振华、方斌通过计算 T N T 水中爆炸冲击波传播情况，分析了网格密度和 人工粘性对D Y T R A N 数值模拟结果的影响，为了逼 近冲击波峰值压力远场经验公式计算值，他们还对 水的状态方程参数进行了大幅度修改引。 A U T O D Y N 软件是1 个显式非线性有限差分程 序，它拥有E u l e r F C T 、E u l e r ．G o d u n o v 、A L E 二维和三 万方数据 第2 6 卷第3 期 刘科种，徐更光，辛春亮，等A U T O D Y N 水下爆炸数值模拟研究 1 9 维求解器，可用来解决固体、流体、气体及相互作用 的高度非线性动力学问题，在水中兵器威力仿真和 舰船抗爆炸冲击数值模拟方面具有重要的应用价 值，采用A U T O D Y N 数值计算软件，进行了T N T 水 中爆炸的数值模拟，并对影响A U T O D Y N 计算精度 的因素进行了分析。 1 数值计算模型 数值计算用炸药为lk gT N T ，计算水域为1 0 m ，采用一维轴对称楔形网格计算模型，模型中x 坐 标的起点位置不是位于炸药的中心，而是距离炸药 中心不到0 ．4c m 处 不到炸药半径的1 0 ％，由此带 来的质量上的误差不到O ．1 ％，可以忽略不计 ，这 样做的目的是为了避免炸药中心处微小单元尺寸带 来计算上的不稳定。 1 ．1 炸药的状态方程 炸药爆轰产物的压力用ⅣL 状态方程来描述， 其表达形式为 p 州 卜南 e 埘，枷 1 _ 茄 e 啦～等 式中，E 为单位质量内能；V 为比容；A 、召、R 。、R 2 、∞ 为常数。其中，方程式右端第1 项在高压段起主要 作用，第2 项在中压段起主要作用，第3 项代表低压 段。在爆轰产物膨胀后期，方程式前2 项的作用可 以忽略，为了加快求解速度，将炸药爆轰产物J w L 状态方程转换为更为简单的理想气体状态方程 绝 热指数7 ∞ 1 。A U T O D Y N 材料数据库中T N T 的J W L 状态方程参数有2 套，分别来自L E EF I N G E R 1 9 7 3 和D O B R A T Z 1 9 8 1 ∞引，根据计算经验 表明，采用后者提供的参数的计算结果更贴近经验 计算值，其具体参数见表I 。 表1T N T 炸药J W L 状态方程参数 1 ．2 水的状态方程 水采用多项式状态方程，材料受压时 P A l u A 2 1 1 , 2 A 3 1 1 , 3 B o B 1 “ p o E 材料受拉时状态方程为 P A l “ A s “2 A 3 u 3 B o B 1 u p o E 式中，A 。、A 、A 3 、B ”B 1 、L 、疋为常数， Ⅱ p ／p 耐一1 ，E 为单位质量内能，计算方法为 E t p g h P o 、 ／p ／B o 这里p 是水的密度，h 是水深，P 。是大气压力， 如果忽略水深的影响，水的比内能为1 0 13 2 5 ／ l0 0 0 ／0 ．2 8 3 6 1 ．8 7 5J ／k g 。 水这种材料难以承受负压引起的拉伸。一些文 献认为纯净的水最多能够承受1 个大气压的负压， 而生活中的水含有气体杂质，这使之只能承受的负 压为零，在计算中将截止压力统一设成零，水的状态 方程参数如表2 所示。 表2 水的状态方程参数 1 ．3 网格密度的选取 容许网格尺寸为 爆炸初期炸药中的爆轰波和水中冲击波是高频 A d C O ／1 0 0 波，炸药及其附近水域需要足够细密的网格才能反 映出这种特性，否则计算出的峰值压力会小于实际 值。黄正平认为，在水下爆炸测试时，测量系统的采 样速率是按峰压衰减时间常数百采样点法确定 的M 】，也就是在峰压衰减肘间常数内至少采样1 0 0 次，否则难以捕捉峰值。根据文献[ 7 ] ，’I N T 近场爆 炸时间常数计算公式为 0 0 ．4 5 R n r ．0 4 5 1 0 - 1 0 ．4 5 R n4 5 ／R o o ‘5 5 1 0 _ 3 这样，为准确捕捉冲击波峰值超压，数值计算时 这里，尺。为装药初始半径，r 为爆心到测点距 离，尺是当地水中声速，对于远场水下爆炸，C 14 8 3r n ／s ，对于中近场，水中声速要大于此值。 从公式可以看出，由于炸药特征尺寸的影响，近 场衰减时间常数并不符合爆炸相似律，而且药量越 大，衰减时间常数越小，由此计算出的容许网格尺寸 也就越小。而且，在近场范围内，邻近爆心处的测点 虽然衰减时间常数小，但由于水中声速很大，对应的 容许网格尺寸可能会大于距离爆心稍远处测点的容 万方数据 2 0爆破 2 0 0 9 年9 月 许网格尺寸。 对于1k gT N T 水中爆炸，距离爆心0 ．3m 水中 声速为19 4 0m ／s ，按公式计算出的衰减时间常数为 5 ．2 1 0 ～s ，如果按照峰压衰减时间常数内采样1 0 0 次的标准划分网格，距离爆心0 ．3m 处网格大小不 应该大于l 1 0 一m 。如果计算水域较大，即使采 用一维模型，计算规模也将极其庞大。作者建议，首 先采用小模型进行计算，计算水域为2m ，网格划分 总数为40 0 0 ，即网格大小为0 ．5 1 0 一m 。冲击波 将要传播到边界时终止计算，然后采用大模型，计算 水域为1 0m ，网格划分总数为50 0 0 ，即网格大小为 2 1 0 3m ，并把小模型计算结果映射到大模型中继 续计算。 I ．4 人工粘性系数 对于冲击波传播这种强间断问题，在数值计算 中，往往通过引入人工粘性系数来制造1 个小的物 理量的过渡层，解决微分方程不连续问题。但人工 粘性在几个网格内抹平了冲击波，过大的人工粘性 系数使得峰值压力计算值小于真实值，人工粘性系 数越大，计算峰值压力越低，而且上游的相对误差会 累积到下游，导致下游冲击波峰值压力的计算值更 加偏低，但较小的人工粘性系数难以抑制峰值过后 曲线的剧烈伪振荡。A U T O D Y N 中人工粘性一次 项、二次项系数缺省值分别为O ．2 和1 ，参考文献 [ 8 ] ，将人工粘性系数调低为O ．0 2 和0 ．1 ，计算完毕 后，对计算曲线低通滤波，去掉伪振荡。 2 数值计算结果分析 2 ．1C O L E 和Z A M Y S H L Y A Y E V 经验公式 C O L E 1 9 4 8 很早给出了简单明了、应用广泛 的水下爆炸冲击波及其能量和冲量公式一】。 Z A M Y S H L Y A Y E V 1 9 7 3 在C O L E 研究成果的基础 上作了进一步的发展，将水下爆炸爆炸载荷分为5 个阶段‘7 | 指数衰减阶段、倒数衰减阶段、倒数衰减 后段、气泡膨胀收缩阶段和脉动压力段，对于球形 T N T 装药，前4 个阶段对应的计算公式为 P t P 。e 叫8 ，t 如 譬 竿 5 ．9 7 8 酽菩一 3 “ 6 5 ㈡一警蹦一， 。E 蔫W 1 ／31 3 ．1 ．5 喜≥ ，‘。 57 6 8 矿w l ／3 ．} 0 ．8 9 ； ‰- 8 ．1 4 1 0 4 ∥力 警 2 ∞； ‰警 譬； ”i 3 ．5R o ．厕 F ≥3 0 ； ”譬㈣毗 m 1 1 ．4 1 0 ．6 ／F 。1 3 1 ．5 1 ／F 1 2 6 F 瓦R ；f 挚r 观- ，黪； 7 2 瓦5‘2 - ￡； 2 。22 1 1 矛； △P 譬 56 3 5 产5 4 - 0 ．1 1 3 P 0 1 ．～2 ； ∥ 杀器52m；天f ． 一 ““’ 铲 豢一劳 m 譬； 百万知“9 1 0 。叩m 习意； P o P 。m P g H o ； F 一二旦 ‘o P ．m 。 式中，P 。为冲击波峰值压力，P a ；0 为冲击波时间常 数，s ；形为T N T 药球质量，k g ；R 为爆心到测点的距 离，m ；R 。为药球初始半径，m ；t a 为冲击波到达时间， s ；t 。为冲击波正压作用时间，s ；P 。为爆心处流体静 压，P a ；P 。为大气压，P a ；c 为水中声速，m ／s ；％为爆 心的初始深度，m 。 2 ．2 数值计算与经验公式对比 图1 一图3 是距离爆心不同距离处冲击波压力 历时曲线数值计算与经验公式计算结果对比情况， 表3 、表4 列出了这些距离处冲击波峰值超压和冲 量大小 积分时间取5 倍时间常数 对比情况。 万方数据 第2 6 卷第3 期刘科种，徐更光，辛春亮，等A U T O D Y N 水下爆炸数值模拟研究 2 l 小的粘性系数所致，但这不影响冲击波峰值压力和 冲量大小的数值。从表3 和表4 也可以看出，冲击 波峰值超压误差控制在5 ％以内，冲量大小误差控 制在1 0 ％以内，可以满足工程计算的需要。数值计 算值均低于经验公式计算结果。 通过模拟结果和经验公式结果的对比，可知模 拟结果和经验结果是相符合的，因此，可以认为在数 值计算中，通过选择合理的模型和计算参数，可以较 好地模拟冲击波的传播情况。 图1 距离爆心3m 处冲击波压力曲线 3 ，J 、结 图2 距离爆心5m 处冲击波压力曲线 C - 图3 距离爆心7m 处冲击波压力曲线 表3 峰值超压数值计算结果与经验公式对比 表4冲量大小数值计算结果与经验公式对比 从图中可以看出，在3m 、5i n 和7m 处，数值计 算冲击波压力衰减曲线与经验公式计算结果符合得 较好，基本重合。只是数值计算冲击波压力历时曲 线波形上存在一定的振荡现象，这是由于选择了较 计算了Ik g ，I N T 炸药水下爆炸冲击波的传播 过程，数值计算结果与经验公式结果吻合得较好，能 够较为准确的反应水中爆炸冲击波传输过程，通过 计算表明 1 爆炸初期炸药中的爆轰波和水中冲击波是 高频波，炸药及其附近水域需要足够细密的网格才 能反映出这种高频特性，建议按照峰压衰减时间常 数内至少采样1 0 0 次的标准划分网格。 2 采用较小的人工粘性系数能够减小对冲击 波峰值的抹平，这样计算出的冲击波压力曲线与经 验公式计算曲线吻合得较好，但过小的人工粘性系 数难以抑制峰值过后曲线的剧烈伪振荡，为此，需要 对计算曲线进行低通滤波。 3 在数值计算中，通过建立合理的模型并选取 合适的模型参数，A U T O D Y N 不但可以准确捕捉冲 击波峰值压力，而且能够比较准确地计算冲击波在 水中传输时的衰减演化过程。 参考文献 [ 1 ] 张振华，朱锡，白雪飞．水下爆炸冲击波的数值模拟 研究[ J ] ．爆炸与冲击，2 0 0 4 ，2 4 2 1 8 2 - 1 8 8 ． [ 2 ] 方斌，朱锡，张振华，等．水下爆炸冲击波数值模 拟中的参数影响[ J ] ．哈尔滨工程大学学报，2 0 0 5 ， 2 6 4 4 2 0 - 4 2 4 ． [ 3 ] A N S Y SA U T O D Y NM a n u a l sV e r s i o n1 1 ．0 [ M ] ．A N S Y S ，2 唧． [ 4 ] L E EF i n g e r ，C O L L I N S ．肌E q u a t i o n so fS t a t eC o e f f i c i e n t s l 矗强S hE x p l o s i v e s [ J ] ．U C I D - 1 6 1 8 9 ，J a n u a r y1 9 7 3 ． [ 5 ] D O B R A I ZBM ．L a w r e n c eL i v e r m o r eN a n o n a lL a b o r a t o r yE x - p l o s i v e sH a n d b o o k [ M ] ．U C R L - 5 2 9 9 1 ，L a w r e n c e ，C A ，1 9 8 1 ． [ 6 ] 黄正平．爆炸与冲击电测技术[ M ] ．北京国防工业出 版社，2 0 0 6 ． [ 7 ] Z A M Y S H L A Y E VBV ．D y n a m i cL o a d si nU n d e r w a t e rE x p l o s i o n [ J ] ．A D - 7 5 7 1 8 3 ，1 9 7 3 ． [ 8 ]辛春亮．高能炸药爆炸能量输出结构的数值分析[ D ] ． 北京理工大学博士学位论文，2 0 0 8 ． [ 9 ] C O L ERH ．U n d e r w a t e rE x p l o s i o n s 。P r i n c e t o nU n i v e r s i t y P r e s s [ M ] ．P r i n c e t o n ，N J ，1 9 4 8 ． 万方数据