基于DYTRAN软件的水下爆炸数值计算

第17卷第6期2003年12月　华　东　船　舶　工　业　学　院　学　报(自然科学版)
Journal of East China Shipbuilding Institute(Natural Science Edition)
Vo1117No16
Dec.2003
文章编号:1006-1088(2003)06-0011-06
基于DYTRAN软件的水下爆炸数值计算
邓文彬,王国治
(华东船舶工业学院机械系,江苏镇江212003)
摘　要:采用瞬态动力学有限元软件D YTRAN,对舰船设备在水下爆炸冲击波作用后的动态响应进行
了研究。

介绍了如何利用D YTRAN软件模拟水下爆炸过程,并提出了优化计算文件来节约计算时间提
高解题效率的方法。

最后,对某舰船主机的隔振装置的抗冲击性能进行了分析。

关键词:D YTRAN软件;水下爆炸;冲击
中图分类号:TB535 文献标识码:A
Numerical Analysis of Explosion U nder w ater with DYTRAN Soft w are
D EN G Wen2bi n;W A N G Guo2z hi
(Dept.of Mechanical Eng.,East China Shipbuilding Institute,Zhenjiang Jiangsu212003,China)
Abstract:Principles and methods for dynamic response computation of ship to underwater explosion are in2 troduced in this paper.First,how to simulate explosion underwater by D YTRAN software.Second,The cost-effectiveness of the analysis can be increased by modifying document of computation.At last,an ex2 ample about the computation of shock response of marine diesel engine resilient mounting system is given in order to show the validity of the method.
K ey w ords:D YTRAN software;explosion underwater;shock
0　引　言
舰船设备抗冲击分析,对于优化舰船以及设备设计,改善船舶的工作性能,具有重要的指导意义。

由于水下爆炸的复杂性及计算方法和软件的缺乏,目前国内对舰船设备抗冲击分析还很不完善。

本文介绍了利用国外有限元计算软件D YTRAN对舰船主机隔振装置的抗冲击性能进行了研究。

1　水下爆炸问题计算分析的数值方法及软件
由于实船水下爆炸试验不容易实现,所以目前分析水下结构受爆炸冲击波作用后的动态响应主要局限于经验法、解析法和数值法3大类。

经验法是根据一些标准冲击曲线,如半正弦波,梯形波等,进行收稿日期:2003-09-16
基金项目:企业协作技术攻关课题(院编2001206)
作者简介:邓文彬(1978-),男,湖南祁东人,华东船舶工业学院硕士研究生。

冲击试验,得到一些经验公式及图表来估算结构的冲击加速度响应,此法简单,但精度差,不适合要求较高的场合;解析法是详细分析水下结构受爆炸冲击后的表面吸收压力,再分析结构的动态响应;此法精度高,但只适用于结构及边界几何形状较规则简单的问题,如球壳、无限长圆柱壳和椭球壳等。

用经验法和解析法来较精确求解水下复杂结构受爆炸冲击波作用后的动态响应是非常困难的,而采用数值方法能很好地解决此类问题。

D YTRAN利用先进的显式中心差分算法来求解结构动力响应[1],而大多数有限元程序采用隐式时间积分法来求解动力问题。

隐式算法利用系统在t+dt时刻的运动微分方程来求解t+dt时刻的响应[2],其等效刚度矩阵中含有系统刚度矩阵的贡献,积分时间步长只与结构的固有频率有关,可以做到条件稳定和无条件稳定;显式算法利用系统在t时刻的运动微分方程来求解t+dt时刻的响应,其等效刚度矩阵与系统刚度矩阵无关,这种算法是条件稳定的,它要求时间积分步必须取得很小。

对于解决瞬态、高度材料非线性、高度几何非线性、模型规模较大的问题时,采用隐式方法在积分的每一步上得到的线性代数方程组的系数矩阵(即等效刚度矩阵)是变化的,因而需要重新生成;而采用显式算法时省去了等效刚度矩阵地重新生成的工作,计算量相对比隐式法要小。

所以在求解像水下爆炸这类高度非线性问题时,采用显式法较隐式要优越得多。

D YTRAN具有丰富的材料模式、状态方程(J WL炸药方程)及各种起爆条件,能够用于模拟爆炸波的传播和爆炸产物的运动,以及爆炸冲击波对结构的响应。

而且在众多的商业有限元分析软件中, D YTRAN是提供精确的或称完全的流固耦合功能的优秀软件之一。

通常的软件在处理流体-结构相互作用问题时,将流体产生的力作为“预先确定”的载荷作用到结构上进行分析,而D YTRAN则不然,它是通过直接耦合结构网格(Lagrange网格)和流体材料网格(Euler网格)间的响应自动地、精确地算出每一时间步流固界面处的物理性质。

在这个过程中,一方面,Euler材料流动引起的压力载荷通过耦合算法自动作用到结构的有限元网格上。

在这种压力作用下,结构的有限元网格将发生变形,结构的变形也反过来影响Euler材料的流动和压力值,这种结构变形和流体载荷间的相互影响使得我们可以得到完全耦合的流体-结构响应。

2　有限元法舰船水下非接触爆炸冲击环境计算
水下非接触爆炸是海军舰船主要的冲击环境能源。

众多的海上实战及实船水下爆炸冲击试验结果表明,水下非接触爆炸虽对舰船整个生命力来说并不形成致命的威胁,但可导致舰船上各类重要设备广泛的冲击破坏。

因此,舰船设备,尤其舰船动力装置的抗冲击性能研究是个重要的课题。

2.1　爆炸冲击波的理论分析
球形炸药在水中爆炸后形成的爆炸产物,其压力、温度远远大于周围水介质,所以爆炸产物迅速而剧烈地压缩周围介质,使水的压力呈现突跃式升高,这就是初始的冲击波。

流过爆炸反应区的物质质量,内能和动能是守恒的,可以用Rankin2Hugoniot方程来表示:
质量和动能守恒p-p0=D 2
V20
(V0-V)(1)
内能守恒e-e0=1
2
(p+p0)(V0-V)+q0(2)式中:p0,V0,e0,q0为冲击波前锋的初始压力,初始单位体积,初始单位内能和初始内能;p,V,e,q为冲击波后端的终时压力,终时单位体积,终时单位内能和终时内能;D为爆轰速度。

D YTRAN可以利用Jones2Wilkins2Lee(J WL)炸药方程式来模拟各种炸药的爆炸过程。

ρ=A1-ωη
R1
e-
R
1
η+B1-
`η
R2
e-
R
2
η+`ηρ0e(3)
式中:e为单位质量内能;ρ0为材料参考密度;η为密度比值;A、B、ω、R1、R2为常数。

对于密度ρ0为1717 kg/m3黑索金炸药,A为5.24229×1011;B为0.07678×1011;R1为4.2;R2为1.1;ω为0.34;e为4.95×106 21华东船舶工业学院学报(自然科学版)2003年
J/kg 。

D YTRAN 也可以利用理想气体状态方程来模拟各种炸药的爆炸过程,它是直接模拟爆炸后冲击波传播过程和气泡膨胀过程。

根据能量守恒定律保证炸药和爆炸气体的能量相同。

当气体比热γ为常量时,气体爆炸前后的状态参数可以用以下方程来表示:
p =2(γ-1)q 0p 0
(4)V =γ
γ+1V 0(5)
冲击波在深水部分中传播时基本上保持球对称波峰面。

在无限水域爆炸条件下,对于一定炸药当量,波峰面压力峰值随传播距离之间的规律可由库尔公式确定:
p m =W 13
R a ・k (6)
式中:p m 为压力峰值;W 为炸药当量;R 为距爆心距离。

对于TN T 炸药而言,α与k 分别为1.13和5.33×107。

对于其它炸药,可根据能量等效原理来估算冲击波的压力峰值。

爆炸产生的冲击波扫过流场中的某点时,该点的瞬时压力与时间的关系为:
p l =p m ・e -1φ 0<t <φ
(7)p l =0.368p m φ
t
φ≤t <150～100φ(8)式中:φ为指数衰减的时间常数;e 为自然对数;对于TN T 炸药而言,指数衰减的时间常数φ与装药量W 及爆距R 的关系如下:
φ=W
-13R -0.24・W -13×10-4(9)
但在浅水部分,由于受自由界面的影响,以及当水中存在固体结构时,受流固耦合的影响,流场中的压力规律自然不同于库尔公式所表示的。

此时流场中的压力等于入射压力(即冲击波压力)、绕射压力、辐射压力和稀疏波压力之和,从而形成不规则的反射区。

图1为50kg TN T 爆炸后耦合流场中的压力峰值与库尔公式的比较。

从图中可以看出当爆距接近30m 时(爆距30m 处存在船体结构和自由水界面
,因而在附近水域形成不规则区),不规则区的压力峰值要低于规则区的压力峰值,从而消弱了爆炸冲击波的威力。

a )库尔公式
b )耦合流场
图1　爆炸后冲击波压力峰值库尔公式与耦合流场之比较
Fig.1　Comparison of pressure
2.2　冲击波的数值模拟
图2所示算列为40m ×15m ×42m 的有限水域,采用欧拉实体单元,单元数201601个,节点数213436个。

船体总长22.2m ;型宽3.25m ;型深1.33m ;设计吃水0.5m 。

炸药水深30m ,处于舯部正龙骨下方。

图2给出了50kg TN T 爆炸后,存在船体结构耦合时流场中的冲击波扩散时历云图。

31第6期 邓文彬等:基于D YTRAN 软件的水下爆炸数值计算
从图3可以看出,爆炸冲击波在水中传播的速度和声速非常接近,大约1500m/s 。

冲击波经过某点时,该点的瞬时压力与时间的关系呈指数分布,与理论趋势非常接近。

在离爆炸中心距离较近时,冲击波出现几个峰值,这是因为爆炸冲击波余能的扩散,使周围流体的压力一次或多次反复迅速上升。

a )13ms
b )20ms
图2　存在船体结构耦合时,50kgTN T 爆炸后流场中球形冲击波扩散时历图
Fig.2　Shock wave pervasion
sequence
a )爆距5m
b )爆距10
m
c )爆距20m
d )爆距30m
图3　50kgTN T 爆炸后纯流场中某点冲击波压力时历曲线
Fig.3　The cuves of pressure 2time
3　带浮筏装置的舰船机舱段在爆炸冲击波作用下的动态响应
水下爆炸问题涉及到水下冲击波和船体结构的相互耦合作用,具有较强的非线性因素。

由于一般模型总的单元数在6位数以上,计算周期长。

因此,数值计算工作量非常巨大,在不改变硬件设备的情况下,寻找可行的方法以提高计算效率,并且使计算精度在许可的范围内,是非常重要的。

3.1　质量因子法
影响计算的关键因素是积分时间步长。

系统默认,积分时间步长由D YTRAN 自动计算,用户只能定义初始时间步长,这样确保最小积分时间步长小于应力波穿过最小单元的时间,使得计算趋于稳定。

在每个积分过程中,时间步长将重新计算。

自动时间步长适合于大部分计算分析,一般不要修改它,除
41华东船舶工业学院学报(自然科学版)2003年
非你有更好的原因。

修改时间步长,可以定义参数STEPFCT 值,内部计算得到的时间步长乘以参数值所得到的时间就是实际所使用的时间步长。

系统默认STEPFCT 值为0.666,对于大部分计算模型,STEPFCT 值可以设置到0.9,但不能超过1,否则计算不稳定。

除了定义参数STEPFCT 值以外,D Y 2TRAN 还提供了一种质量因子法(mass scale )来增大时间步长,提高解题效率。

所谓质量因子法就是在每个单元上增加计算质量,使积分时间步长总是大于所定义的最小时间步长。

尤其在以下2类问题中,质量因子法可以提高解题效率。

①如果有限元网格中包含有一些非常小的单元,而时间步长由最小的单元决定的。

②如果在分析过程中,一些单元发生了严重的变形而失效,而时间步长由变形最严重的单元决定的。

时间步长太小会终止计算的进行。

出现时间步长太小的原因之一是建模过程中没有考虑结构的失效,但是,这些失效单元有时只存在相对较小的区域,没有必要建立全局的失效模型。

因此,质量因子法是个很好的方法来人为地提高计算效率。

质量因子法可以通过在计算文件3.dat 定义参数SCAL EMAS 来激活。

PARAM ,SCAL EMAS ,D TM IN ,MAXPERC ,STEPS
PARAM ,STEPFCT ,VAL U E
此时,最小时间步长dt =STEPFCT 3D TM IN ,积分时间步长总是大于dt 。

3.2　计算文件的优化
除了增加积分时间步长提高计算效率外,通过修改D YTRAN 的计算文件3.dat 也可以提高每个积分步的计算速度:①将Analysis 菜单中Output Controls 子菜单下的Write COMM2’s Summary 设置为NO 。

这个参数定义是否输出集中质量以及它们的内能和动量。

设置NO 可以节省内存和计算时间。

系统默认为YES 。

②在Parameter Section 域中定义COSUBMAX 参数,这个参数定义更新Euler/Lagrange 耦合面的最大子循环数目。

而在每一个子循环中,将不更新耦合面。

系统默认在每个子循环中耦合面将被更新一次,更新一次耦合面占有大量的计算时间,参数值越大,计算时间将越少。

这使得在微机上进行大型计算变得可能,但同时会导致精度降低。

一般,如果耦合面变化比较缓慢,可以采用较大的参数值。

③在Parameter Section 域中定义COSUBCYC 参数,这个参数定义在耦合计算中,子循环的最大增长步。

系统默认为1,可以设置更大值。

④在Parameter Section 域中定义EL SUBMAX 参数,这个参数定义单元子循环中的最大组群数。

系统默认为4,可以设置1<Value <21。

增大参数值可以使系统允许时间步的比率大于8。

但同时可能导致计算的不稳定以及精度的降低。

每个单元组的更新时间为:dtgroup =233(ngroup )3dtmin (式中:ngroup =1,2…,EL SUBMAX ,dtmin 为最小的稳
定时间步。

)⑤在Parameter Section 域中定义EL SUBCYC 参数,单元子循环运算法则把单元分成等时间步的各个单元组群。

每个单元组随组时间步长而更新。

特别地,当非常小的单元决定时间步长时。

可以节省大量的计算时间。

系统默认为NON E 。

笔者经过大量的实验证明,光是修改COSUBMAX 参数和EL SUBCYC 参数可以使计算速度提高5～10倍,但同时会引起计算上的偏差。

其偏差程度与计算模型和修改的参数有关。

对于像船体结构在水下爆炸冲击波作用下的动态响应这类问题,因为爆炸冲击波压力以及冲击波作用在结构上的力随时间变化非常快,为微秒量级。

所以COSUBMAX 参数值越小,每一积分时间步作用在船体结构上的冲击力越接近实际值。

随着COSUBMAX 参数值的增大,计算偏差越大,在爆炸冲击波压力增大到峰值之前,计算值要小于实际值;在冲击波衰减过程中,计算值要大于实际值,从整个爆炸过程上来讲,CO 2SUBMAX 参数值过大,更新流固耦合面的时间过长,就会错过爆炸冲击波最大峰值压力,使作用在结构上的最大冲击力小于实际值,计算结果偏小。

读者可以根据需要在误差允许的范围内设置其中几项或全部参数。

3.3　计算结果及分析
算例为某艇主机舱段[3],如图4所示,直径为10m ,总长3m ,主机湿重:4500kg ;几何尺寸:3250mm ×1350mm ×1650mm ;筏体几何尺寸:3500mm ×2700mm ×715mm ,筏体质量:2150kg ;浮筏
51第6期 邓文彬等:基于D YTRAN 软件的水下爆炸数值计算
上层和下层隔振器刚度和阻尼:k =7×106N/m ,c =5000N ・m/s ;筏体与船体的相连基础部分采用16Mn 钢;水域为30m ×30m ×40m.炸药当量为10kg TN T ;球形药包;艇深30m ,炸药爆距5m 。

本模型中节点数为191029,单元总数为130543。

本文算例中的部分计算文件3.dat 如下(其中最后两行参数用来定义质量因子法):
＄———Parameter Section ———
PARAM ,CONM2OU T ,NO
PARAM ,COSUBMAX ,10
PARAM ,COSUBCYC ,2
PARAM ,EL SUBMAX ,12
PARAM ,EL SUBCYC ,ALL
PARAM ,SCAL EMAS ,1E -6,100,1
PARAN ,STEPFCT ,0.9
图5为取某艇壳体、筏体下层、筏体上层以及机组上关键点的加速度响应值。

从图5可以看出,机组的垂向最大冲击加速度a 小于10g (机组无损坏最大加速度估算值),浮筏起到了良好的抗冲击作用。

图4　某艇主机舱段有限元模型图
Fig.4　FEM of the submarine ’s engine room 图5　主机舱段加速度时历曲线Fig.5　Acceleration response of engine room
4　结　论
1)利用D YTRAN 计算软件来模拟水下爆炸问题,球形炸药在水中爆炸后以球形冲击波的形式在水中扩散,速度和声速基本接近。

冲击波扫过流场中某点时,该点的压力随时间的衰减呈指数分布,模拟结果和理论基本一致。

2)质量因子法是个很好的方法,它使时间积分步长总是大于所定义的最小时间步长,提高了某些模型的计算效率。

3)通过修改D YTRAN 计算文件3.dat 中的参数,在误差许可的范围内可以大幅度提高模型的计算速度,使得大型计算在微机上可以实现。

参考文献:
[1]　D YTRAN User ’s Manual[Z].MSC ,2001.
[2]　柳贡民,刘志刚.以逐步积分法计算舰船柴油机隔振装置冲击响应[J ].哈尔滨工程大学学报,1997,18(4):19-24.
[3]　邓文彬,王国治.舰船浮筏的声学特性及水下爆炸冲击响应分析[A ].第九届船舶水下噪声学术讨论会[C].苏州:
《船舶力学》编辑部,2003.67-73.
(责任编辑:汪时美)
61华东船舶工业学院学报(自然科学版)2003年。