水下爆炸冲击波的数值模拟研究

书书 ! 第 "# 卷! 第 " 期 爆 炸 与 冲 击 $%&’ "#，(%’ "! ! "))# 年 * 月 +,-./01/( 2(3 04/56 72$+0 89:;<，"))#! 文章编号：=))=>=#??（"))#）)">)=@">)A 水下爆炸冲击波的数值模拟研究 ! 张振华，朱! 锡，白雪飞 （海军工程大学船舶与海洋工程系，湖北 武汉 #*))**） ! ! 摘要：应用商业有限元程序 805’ 3BCD2(数值模拟了球形药包在无限水域中爆炸产生的冲击波。通过 和经验公式计算结果的比较，证明采用合理的计算参数和有限元模型能够较好地模拟水下爆炸冲击波的传 播过程。通过调整水的状态方程参数，达到了提升冲击波峰值应力的效果，从而有效地降低了单元数量和计 算时间。 ! ! 关键词：爆炸力学；水下爆炸；数值模拟；冲击波；805’ 3BCD2( ! ! 中图分类号：/*@"’ =；EFA#’ A)*! ! ! 国标学科代码：=*)·*?")! ! ! 文献标志码：2 !" 引" 言 ! ! 对水下爆炸冲击波的研究一直是舰船抗爆防护的重点，然而由于其复杂性，该领域一直主要以实验 研究为主［= G #］。近年来随着计算机技术和计算理论的快速发展，使得人们可以通过数值模拟的方法对 水下爆炸的各种现象进行预报。世界上相继开发了一些大型商业有限元程序（如 (20CD2(、231(2、 3B(2*3、2062、3BCD2(等），其中不少程序，尤其是 805’ 3BCD2( 取得了较为广泛的应用。该软件 具有流体>结构相互作用的不同的描述方式，流体可以用 +H&I:J9K 单元描述，固体则采用 .9L:9KLI 单元 模拟。美国、日本和意大利都曾考核了计算软件模拟水下爆炸的可行性，认为 805’ 3BCD2( 是一种合 适的计算分析软件［?］。可见 805’ 3BCD2(代表了当今世界上比较先进的爆炸力学计算程序，同时可 以利用商业有限元程序强大的前后处理功能，有可能真实地再现水下爆炸冲击波产生和传播的过程。 ! ! 水中爆炸冲击波的数值模拟是进一步研究结构在水下爆炸载荷作用下动态响应的基础。在本文 中，探讨了采用 805’ 3BCD2(进行水中爆炸冲击波数值模拟的相关问题，计算结果表明，采用合适的 计算参数和有限元模型，能够较好地模拟水中爆炸冲击波的各状态参量。同时研究了网格密度和药包 大小对计算结果的影响。最后通过调整计算参数提升冲击波峰值压力，大大减少了水单元数量和计算 所需时间，同时分析了这种方法的利弊。 #" 水下爆炸冲击波的数值模拟 #$ !" 欧拉求解器简介 ! ! 805’ 3BCD2(中包含拉格朗日求解器和欧拉求解器，其中欧拉求解器主要用于流体流动问题的分 析以及固体材料发生很大变形的情况。当采用欧拉方法时，节点固定在空间中，由相关节点连接而成的 单元仅仅是空间的划分，分析对象的材料在网格中流动，材料的质量、动量和能量从一个单元流向另一 个单元。因而对于水中爆炸现象，采用欧拉方法求解要合适一些。805’ 3BCD2(中的欧拉求解器在空 间上的离散采用控制容积法，在时间域上的离散采用显式积分法。采用显式积分法的优点是不需要作 矩阵分解，稳定性与积分的时间步长无关。但是对于显式积分法，要保持计算的稳定，积分步长必 须小于网格的最小固有周期，即时间步长必须小于应力波跨越网格中最小元素的时间。3BCD2(程 ! 收稿日期："))*>)?>"F；修回日期："))*>==>"@ ! ! ! 作者简介：张振华（=MAF—! ），男，博士研究生，讲师。 万方数据 书书书 序所采用的的单元均为一阶单元，基本单元为八节点的任意六边形单元，还有三菱柱单元及四节点的四 面体单元。 !" !# 状态方程的参数设置 ! ! 状态方程用来定义固体或流体在各种不同状态下的压力和密度以及比内能之间的函数关系，正确 选取状态方程中的参数对于计算结果至关重要。水和 "#"材料状态方程的参数设置如下述。 !" !" $# 水的状态方程 ! ! 在化爆情况下，水下爆炸冲击波压力一般小于 $%&&& 个大气压，属于中等强度冲击波。这时爆炸冲 击波通过介质后熵值变化很小，接近于等熵过程，由文献［’］得到水的等熵状态方程为 ! " # ! !( )& $ %[ ]( （(） 式中：# ) *&+ ,% -. / 0$；$ ) 1+ (% 。 ! ! 234+ 56789:中并无式（(）形式的状态方程，程序包含多项式状态方程如下所示 ! " &(" ’ &$" $ ’ &*" * ’（(& ’ ((" ’ ($" $ ’ (*" *）!&) （$） 式中：!& 为水在常温状态下的密度；)为单位质量的内能，" )（! / !&）; (。 ! ! 因此需要将式（(）式变换为多项式形式。将式（(）按 "96<=8展开有 ! " $* (> + (&?" ’ ’* ’? + (&?"$ ’ (* (% + (&(&"* （*） ! ! 对比式（$）可得状态方程中各参数的值。取水的初始密度 !& ) (&&& @. / 0 *，初始质量比内能 )& ) >*+ ?%& @A / @.。图 ( 比较了式（(）和式（*）的计算结果。从图中可以看到，当 " B &+ $ 时，两者吻合得很 好。即使是 " ) &+ $% 时，两者相差也不到 %C。由文献［$］可知当水的密度达到 ($%& @. / 0* 时，冲击波 压力约为 (+ % -D9。而 (&&& @. "#"水中爆炸时，在距离爆炸中心 ( 0处的压力也不过 (+ $’ -D9［1］。因 此这个范围能够满足一般化爆的要求。 !" !" !# "#"的状态方程 ! ! "#"的状态方程采用的 AEF方程［,］ ! " ,（( % #$-( ）)% -( $ ’ #（( % #$-$ ）)% -$ $ ’ #$!&) （,） 图 ( 指数形式状态方程和多项式状态方程的比较 GH.+ ( 4=0I98HJ=: KL7MLL: LNI=:L:7H9< J797L LOP97H=: 9:Q I=<6:=0H9< J797L LOP97H=: 式中：$ " ! . !&；,、#、#、-(、-$ 为常数。 ! ! 对于 "#"装药，, ) *1(+ $ -D9，R ) *+ $*( -D9， S( ) ,+ (%，S$ ) &+ ?%，" ) &+ *。以上参数取自 234+ 2TU3UV#材料库。"#" 的密度取 (%>& @. / 0*，初始 质量比内能为 ,+ (? 2A / @.，初始爆轰速度为 ’+ ?*& @0 / J。 !" %# 水下爆炸数值模拟实例 ! ! 模拟水域 $+ , 0 W$+ , 0，球形 "#"装药设置在 坐标原点，药包质量为 ,>+ % @.，半径为 &+ (?, 0。 由于球形药包爆炸时的球对称性，仅模拟了装药爆 炸时 ( / > 的水域部分。球形药包也只有 ( / > 与水单 元相交，相交部分的材料属性根据其级别来判断，水 单元的级别为 (，"#" 的级别为 $，所以相交部分为 "#"。采用正六面体单元划分单元。在 /，0，1三个方向上均匀划分 ,& 个单元，每个单元长 &+ &’ 0。为 避免冲击波在欧拉网格内产生积压，所有边界设置为流出边界。 ! ! 计算结果如图 $ 所示，图中显示了不同时刻水域中的压力分布即爆炸冲击波在水中的传播情况。 冲击波在 /，0，1三个方向上基本呈球对称分布。可以看到冲击波的波阵面约分布在沿波传播方向上 的五个单元之内。这是因为连续体的有限元模型无法在数值上模拟不连续状态的传播，而冲击波在水 *>(! 第 $ 期! ! ! ! ! ! ! ! ! ! ! ! 张振华等：水下爆炸冲击波的数值模拟研究 万方数据 中传播时会产生强间断，因而程序中引入了人工体积粘性。人工粘性使得波阵面处的压力随应变率增 加，它大约铺在 ! 个单元之上［"］。 图 # 不同时刻水中压力的分布 $%&’ # ()*++,)* -%+.)%/,.%01 02 34.*) 4. -%22*)*1. 505*1. 6 6 为了便于以后研究潜体结构在水下爆炸冲击波作用下的动态响应，应该重点对 789:;< 程序与 =0>*经验公式的冲击波峰值压力、比冲量以及水质点速度的计算结果进行比较。=0>* 的关于冲击波峰 值压力 !5 和比冲量 "的经验公式如下 !5 # $（% ? & @ * & ’）! " # (%? & @*（% ? & @ * & ’） { " （!） 式中：%* 为装药量；’为球面距爆心距离；!，"，$，( 为经验系数，对于 9<9 装药有：! A ?’ ?@，" A B’ "C，$ A !#’ #D E(4，( A !DF" <·+ G 5#。图 @ 和图 H 是 789:;<程序计算结果和经验公式计算结果的比较。 图 @ 冲击波压力时程曲线结果的比较 $%&’ @ =05I4)%+01 02 I)*++,)* .%5* J%+.0)K /*.3**1 *LI*)%*1M* 20)5,>4 41- 789:;< 图 H 比冲量时程曲线比较 $%&’ H =05I4)%+01 02 %5I,>+* J%+.0)K /*.3**1 *LI*)%*1M* 20)5,>4 41- 789:;< 6 6 图 @ 显示了距药包中心 ?’ " 5处的计算压力曲线与经验公式计算结果的比较情况。计算的峰值压 力为 ??@ E(4，经验公式的计算值为 ??F E(4，相差 #’ FN。经验公式计算的压力在 ? 5+时突然跃升到 峰值压力，然后按指数规律衰减。由于有限元不可能反映这种强间断，所以压力的爬升需要一定的时 间。在峰值压力过后，计算值还会有较大的扰动，出现压力的双峰或多峰现象。图 H 显示了 789:;< 与经验公式计算的比冲量的比较情况。在 #’ ! 5+之前，计算值和经验值吻合较好。这段时间实际上是 冲击波的高压作用时间，之后冲击波的压力与峰值压力相比已经很小。此时压力的经验值已经接近为 零，而计算的压力仍然存在小的脉动。这就是计算比冲量在 #’ ! 5+ 之后比经验公式计算值偏高的原 因。 6 6 图 ! 显示了 789:;<计算的水质点速度与理论公式计算结果的比较情况。根据爆炸力学的理论， 冲击波压力和水质点速度有如下关系 H"? 爆6 6 炸6 6 与6 6 冲6 6 击6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 第 #H 卷6 万方数据 !（ "）# $（ "） !% & ! !’! " " $（"）#" （$） % % 即在爆炸产生的球面冲击波的压力场中，水质点的速度不仅与该点的压力有关，而且与扰动到该点 后的压力变化过程有关，即所谓的滞后流。根据式（$）可知，当压力衰减到 " 时，右边的第一项变为 "， 第二项也停止增长，如图 & 中实线所示。然而 ’()*+, 计算的压力在冲击波尾部仍然存在扰动，使得 滞后流速度继续增加，因而数值计算速度比经验公式的结果大。图 $ 显示了冲击波压力峰值在距药包 中心 !- . / .- . 0之间的分布与经验公式的比较。可以看到，计算值和经验值整体吻合较好。 图 & 水质点速度曲线结果的比较 123- & 450678295: 5; <7=>8 678=2?@> A>@5?2=B C>=<>>: >D6>82>:?> ;580E@7 7:# ’()*+, 图 $ 冲击波压力峰值在空间场上的分布 123- $ 450678295: 5; 6>7F 68>99E8> A78B <2=G #29=7:?> C>=<>>: >D6>82>:?> ;580E@7 7:# ’()*+, !" 网格密度对计算结果的影响分析 % % 水域和药包大小与 .- H 节中相同，在水域上均匀划分 !"""，I""" 和 .J""" 个正六面体单元，如果定 义网格密度为水域长度与单元实际长度之比，则以上 H 个模型和 .- H 节中实例模型的网格密度分别为 !"，."，H" 和 K"。图 J 显示了不同网格密度下距离爆炸中心 !- I 0处的冲击波时程曲线的变化情况。 图 J 不同网格密度下冲击波时域曲线 123- J L8>9E8> =20> G29=58B 5; #2;;>8>:= 0>9G % % 随着网格密度的增大，冲击波压力时域曲线由平坦逐渐变陡，压力爬升到峰值的时间逐渐缩短，脉 冲宽度减小。同时峰值压力有了比较明显的提高，峰值压力的提高较网格密度的提高缓慢。另外一个 &I!% 第 . 期% % % % % % % % % % % % 张振华等：水下爆炸冲击波的数值模拟研究 万方数据 特点就是当网格密度大于 !" 后，冲击波峰值压力之后出现了明显的压力脉动，各波峰的压力值基本上 呈指数衰减。当网格密度达到 #" 时，冲击波的峰值压力波形与指数衰减函数已经很接近了。 !" 水状态方程参数调整对结果的影响 $ $ 根据大量的网格密度对计算结果的影响的数值试验得到，要使数值模拟的冲击波压力峰值与经验 公式一致，应该使单元尺寸划分到药包半径的 % & # ’ % & ( 之间。这样 %"" )* 装药在水中爆炸时，水单元 需要划分到 "+ ", -，即要达到厘米量级。对于近距离水中爆炸现象的数值模拟，并不构成什么问题。 但是如果要模拟较远距离的冲击波，或者是对整条舰艇进行水下抗冲击计算，需要划分的水单元数量是 很惊人的。因而希望能够通过调整计算参数，使得在划分较少水单元的情况下，能够达到和细化网格同 样的目的。我们在这方面作了初步的尝试。 $ $ 在 %" -见方的水域上均匀划分 #""" 个水 单元，!，"，# 三个方向的网格密度均为 !"。分 别使式（!）中的系数 $%，$! 和 $( 为式（(）中正 常参数的 %倍。比较距离爆心 .+ ./ - 处的第 %0 号水单元冲击波峰值压力和比冲量的变化 情况。根据 1234的经验公式计算，该处的峰值 压力为 !!+ /0 567。比冲量为 %"+ 80 )9· : & -!。各种状态参数下冲击波峰值压力和比冲 量列于表 % 中。图 , 显示了 ;<=>?9 计算的 冲击波峰值压力和比冲量与经验公式计算结果 之比随状态方程参数倍数的变化情况。 表 # 水的状态方程系数对冲击波峰值压力和比冲量的影响 $%&’( # $)( (**(+,-./ .* 0%1%2(,(1 ,. 0(%3 01(4451( %/6 -205’4( 状态方程参数 的倍数 % 峰值压力 &- & 567 比冲量 ’ &（)9·: & -!） % /+ ,#8 %!+ " (" !%+ #, ,+ #! (/ !!+ 0# ,+ %8 #" !(+ /0 .+ ., /" !/+ ". .+ /! %"" ("+ %( 0+ /. 图 , 水状态方程参数倍数与峰值压力的变化曲线图 @A*+ , 67B7-4C4B 2D :C7C4 4EF7CA2G 2D H7C4B I7BJ HACK L47) LB4::FB4 7GM A-LF3:4 图 8 冲击波波形随参数的变化情况 @A*+ 8 6B4::FB4 CA-4 KA:C2BJ I7BJ HACK L7B7-4C4B 2D :C7C4 4EF7CA2G 2D H7C4B $ $ 从图 , 中可以看到当水的状态方程系数为正常参数的 (/ 倍时，;<=>?9 计算的冲击波峰值压力 和 1234的经验公式很吻合，此时比冲量约为经验公式的 ./N。采用正常参数时，由于网格较粗峰值压 力不到经验公式的 ("N，但比冲量与经验公式较为吻合。当提高状态方程参数时，冲击波峰值压力有 明显提高，但冲击波比冲量随水状态方程参数的增加则逐渐降低。它们增加和降低的速度比参数增加 的速度要慢，从整体上看，比冲量的降低速度较峰值压力为低。因此采用这种方法模拟水下爆炸冲击波 时，需要以牺牲比冲量为代价来提高冲击波峰值压力。根据计算经验，如果网格划分的细一些，可以使 比冲量的降低速度减缓一些。根据冲击动力学理论，当结构响应时间和冲击波的作用时间在一个数量 级上时，结构的动态响应主要取决于冲击波的峰值压力；而当结构的响应时间远大于冲击波的作用时间 时，结构的动态响应主要由冲击波的冲量决定。因此在计算结构在水下爆炸冲击波作用下的动态响应 0,% 爆$ $ 炸$ $ 与$ $ 冲$ $ 击$ $ $ $ $ $ $ $ $ $ $ $ $ $ $ 第 !# 卷$ 万方数据 时，可根据峰值压力和冲量对结构响应的影响程度具体选择参数。 ! ! 图 " 显示了参数提高时冲击波波形的变化情况。与细化网格的效果一样，当参数提高时，峰值压力 提高，但作用时间减少，波形变陡。同时由于压力提高，冲击波波速明显增加，使得冲击波到达该点的时 间明显提前。总的来看，由于单元划分的不够细，没有出现图 # 中的峰值压力成指数衰减的现象，也没 有出现压力脉动，而是类似于三角脉冲载荷。 !" 结" 论 ! ! 从应用 $%&’ ()*+,-数值模拟水下爆炸冲击波的结果来看，按照所取的计算参数和建立模型的 方法，能够很好的模拟水下爆炸冲击波的各种状态量，并得到如下结论。 ! !（.）由于有限元不能反映爆炸冲击波这种强间断，所以当冲击波到达水中某点时，压力的爬升需要 一定的时间，爬升时间与网格密度有关。在峰值压力过后，()*+,-计算的压力值还会有较大的扰动， 出现压力的双峰或多峰现象，而经验公式计算的压力则是平滑下降。一般试验所测得水下爆炸冲击波 压力时域曲线一般是不太规则的强间断曲线，且带有多峰现象。经验公式中将压力时域曲线作为一个 指数衰减函数处理带有一定的近似性。从这一点上来说，数值模拟的这种脉动更符合实际情况。 ! !（/）在冲击波的高压作用时间，比冲量的计算值和经验值吻合较好。之后由于计算压力仍然存在 波动，使得计算比经验值偏高。由于同样的原因，在冲击波压力衰减到很小后，水质点速度的计算值比 经验值要高一些。而冲击波的峰值在空间上的分布，数值计算结果与经验公式基本吻合。 ! !（#）用有限元程序精确模拟水下爆炸现象，必须在一定距离上划分足够细的欧拉网格。当网格密 度较小时，冲击波压力时域曲线过于平坦，脉冲宽度过大。随着网格密度的增大，冲击波压力时域曲线 逐渐变陡，峰值压力有了比较明显的提高，压力爬升的峰值的时间逐渐缩短，脉冲宽度减小。 ! !（0）当模拟较远距离的冲击波或者是对整条舰艇进行水下抗冲击计算时，需要划分的水单元数量 是很大的。这种情况下可以通过调整水的状态方程中的参数来达到提升冲击波峰值压力的目的。当提 高状态方程参数时，冲击波峰值压力有明显提高，但是比冲量却有所降低。因而采用这种方法模拟水下 爆炸冲击波时，需要以牺牲比冲量为代价。计算结构在水下爆炸冲击波作用下的动态响应时，可根据冲 击波峰值压力和冲量对结构响应的影响程度具体选择参数。 参考文献： ［.］! 库尔’水下爆炸［$］’罗耀杰，等译’北京：国防工业出版社，."12’ ［/］! 张鹏翔，顾文斌，叶序双’浅层水中爆炸直达波压力峰值计算方法探讨［3］’解放军理工大学学报，/22/，#（/）：45— 4"’ 67,-8 9:;<=>?@;<，8A B:;=C?;，)D EF=GHF@;<’ (?GIFGG?J; I@KIFK@L?J; M:LHJN JO P:@Q PR:GGFR: GFCS:IL:N LJ GH@KKJT K@U:R T@L:R :>PKJGJ?;［3］’ 3JFR;@K JO 9V, A;?W:RG?LU JO %I?:;I: @;N *:IH;JKJPKJG?J; F;N:R GH@KKJT T@L:R［3］’ YK@GL?;<，."""，.1（.）："—.#’ ［0］! 顾文彬，叶序双，张鹏祥，等’浅层水中爆炸水底影响的试验研究［3］’解放军理工大学学报，/22.，/（/）：44—4Z’ 8A B:;=C?;，)D EF=GHF@;<，67,-8 9:;<=>?@;<，:L @K’ D>P:R?M:;L GLFN?:G JO CJLLJM ?;OKF:;I: ?; GH@KKJT K@U:R T@L:R :>= PKJG?J;［3］’ 3JFR;@K JO 9V, A;?W:RG?LU JO %I?:;I: @;N *:IHKJ’= ?@ 6:&’;’ >&+;A’ ’B:5<=.(A *(=’; .(1.(.7’ 8+7’; ;’A.<( .6 6.C*5+7’= ?@ *6.(A >.+5 DEF :;+5 ;’6*576I D.(+55@，7&’ :’+N :;’66*;’ <1 ?5+67 8+9’ .6 .(>;’+6’= ?@ +=O*67.(A 7&’ :+;+C’7’; .( :<5@(’= ’11’>7.9’5@I 6#( 2),’$：C’>&+(.>6 <1 ’B:5<6.<(；*(=’;8+7’; ’B:5<6.<(，(*C’;.>+5 6.C*5+7.<(；6&<>N 8+9’；FGHI JKL) M#$ QQR 爆S S 炸S S 与S S 冲S S 击S S S S S S S S S S S S S S S 第 T2 卷S " H<;;’6:<(=.(A +*7&<;：!"#$% !&’()&*+ S S S E)C+.5 +==;’66：TU6UV RW3I >