文章编号 : 1672 - 7649 ( 2006 ) 05 - 0022 - 07
余晓菲 ，刘土光 ，张 涛
(华中科技大学 交通科学与工程学院 ，湖北 武汉 430074 )
摘 要 : 应用商业有限元软件 M SC . D YTR AN ，比较了 2种耦合方法 (AL E and Gene r a l Co u p lin g) ，分析计算了
球形药包在近似无限水域中爆炸产生的冲击波 。对比分析了在不同参数设置下冲击波的传播过程 。通过将计算结 果与经验公式相比较 ，证明了 M S C . D Y TR AN 软件是一种合适的模拟水下爆炸冲击波传播的计算分析软件 。
关键词 : 爆炸力学 ；水下爆炸 ；耦合 ；冲击波 O 382 + . 1 中图分类号 : 文献标识码 : A
C om pu t a t i on of the b l a st loa d i n g stren g th of un derw a ter exp l o s i on shock w a ve s
Y U X i ao 2fe i, L I U Tu 2guang, ZHAN G Tao
环球经贸( C o ll ege of Traffi c Sc i ence and Engi nee r i ng, H u az hong U n i ve r sity
of Sc i ence and Techno l o gy, W uhan 430074 , Ch i na )
A b s tra c t :
B l a s t wave p r oduced by sp h e r e T N T cha r ge ex p l o d i ng unde r si m il a r i nfi n ite wa t e r regi on is d iscu ssed and comp u ted by u si ng comm e rc i a l FE M p r ogram M S
C .
D YTR AN. The b l a st wave s a re ana l y z ed when d iffe ren t rea sonab l e p a ram e te rs a re adap ted. B y comp a ri ng the si m u l a ti o n re su lt and the e m p i ri c a l r e 2 su l t , we can p r ove tha t the code M S C . D YTR AN is fit f o r s i m u l a t i ng unde r w a t e r ex p l o s i o n shock wave s .
Key word s :
m e chan i c s of ex p l o s i o n ; unde r w a t e r ex p l o s i o n; coup li ng; shock wave
Coup li ng ) 、快 速 耦 合 法 ( Fa s t G ene ra l Coup li ng ) 、考
虑失效的多重耦合 面法 、AL E 法 ( A rb itra ry L agrange
Eu l e r )及全 Eu l e r 方法 。它 采用 高效 的显 式 积分 技
术 ，支持广泛的材料模型和高度组合非线性分析及 流体 —结构的全耦合 。而且 ，该程序是基于三维非
线性动力学发展起来的计算程序 ，尤其擅长对爆炸
与冲击 ，如水下爆炸等对结构的影响及破坏的设计
优化分析 。 本文首先对比讨论了用 2种常用的耦合法 ，即一 般耦合法和 AL E 耦合法来进行水中爆炸冲击波数值
模拟的区别 。并分析了网格密度 ，边界条件对冲击波 的影响 。同时与经验公式对比 ，讨论了不同药包半径 冲击波的传播过程 。计算结果表明 ，利用正确的计算 参数和有限元模型 , M SC . D YTRAN 软件能够较好地 模拟水中及结构体附近水域不同位置爆炸冲击波的 传播 。
0 引 言
基于在海战中的重要作用 ，自二战以后人们就 开始对水下爆炸开展了系统的研 究 [ 1 - 7 ] 。最初 , 在 缺乏实验检测数据的情况下 , 大多运用库尔 ·P 公 式来计算 ，但要计算出符合实际情况的精确解却是 相当繁琐和困 难 的 。随 后 人 们 又 根 据 实 验 数 据 得 到
了一些半经验半理论公式 ，却也难以揭示水下爆 炸现象的本质 。随 着 计 算 机 技 术 和 计 算 方 法 的 进 步 ，通过数值计算方法求解方程组的数值解已经成 为现实 ，一系列针对水下爆炸模拟计算的软件相继 出现 ，并且通过计算能够得到较好的模拟结果 。其 中 , M S C. D YTRAN 以无可争辩 的 实力 , 赢得 世 界顶 尖高度非线性 、流 —固耦合及瞬态动力响应仿真软 件的地位 。M SC . D Y TRAN 软件一共提供了 5 种处 理流 —固耦合的 方 法 [ 8 ] , 包 括 一 般 耦 合 法 ( Gene ra l 收稿日期 : 2005 - 08 - 29
联系起来, 从而能够分析复杂的流固耦合问题。
1 水下爆炸冲击波数值计算
1. 1 欧拉方法
欧拉方法主要用于流体流动问题的分析以及固体材料发生很大变形的情况[ 8 ] 。它计算的是材料在体积恒定的网格中的运动，材料的质量、动量以及能量从一个元素流向另一个元素，这种方法对模拟计
算水下爆炸非常适合。M SC. D YTR AN 中的欧拉方法在空间域的离散上采用控制容积法，在时间域的离散上采用时间积分法，单元可以是八节点的六面体单元, 单元中可以充满材料，也可以是空的，也可以有一部分空间有材料。一个单元中可以同时有几种材料。材料可以是理想流体，也可以是非理想流体。
1. 2 耦合
水下爆炸是指在极短的时间内，在水下的极小体积内或面积上发生极大能量转换的过程[ 1 ] 。由于爆炸产物高速向外膨胀，在水中就形成了冲击波。结构或潜体承受水下爆炸载荷，使得其在很短的时间内, 在巨大冲击载荷作用下产生一种复杂的非线性动态响应过程。同时，这种爆炸还涉及到冲击波和结构体相互耦合作用的问题。流固耦合问题一直是水下爆炸研究工作中的一个难点[ 9 ] , 目前水下爆炸流固耦合问题的解决方法主要有3 种：一是忽略流固耦合的影响，这种方法不可避免的将产生较大计算误差; 二是采用二阶双重渐进近似法(DAA 法) ，延时法，将压缩流场边界积分法或附加质量法等近似处理，这些方法大多将流体产生的力作为“预先确定”的载荷作用在结构上进行分析; 三是本文中要详细阐述的M SC.
D YTR AN 的精确流固耦合计算方法。这种方法通过直接耦合结构网格(L agrange 网格) 和流体材料网格( Eu l e r网格)间的响应，自动精确地计算出流固界面处的物理性质。在这个过程中, Eu l e r材料
流动引起的压力载荷通过耦合算法，自动作用到结构的有限元网络上，在这种压力作用下，结构的有限元网格将发生变形; 另一方面结构的变形又反过来影响Eu l e r材料的流动和压力值。这种结构变形和流体载荷的相互影响，可以得到完全耦合的流体。
耦合的目的就是为了让欧拉网格和拉格朗日网格发生相互作用。最初，欧拉求解器与拉格朗日求解器是分开的，如果不定义耦合，即使拉格朗日单元恰好处在欧拉网格范围内，也不会对欧拉材料的流动产生任何影响，同时自身也不会受到任何来自欧拉材料D YTRAN 中的耦合有 2 种方式: 一般耦合和AL E 耦合。
一般耦合，根据M SC. D ytran内部的要求，耦合面应当是封闭的，而且必须具有正体积。这就要求所有面段的法线方向指向外面，也就是要将耦合面上所有单元的法线方向调整一致指向外面。拉格朗日网格和欧拉网格可以用在同一个界面相互耦合。该界面是欧拉网格中的材料的流场边界，同时欧拉网格中的材料对界面产生作用力，使拉格朗日网格变形。要将模型的欧拉部分和拉格朗日部分之间建立起耦合关系，首先要在拉格朗日模型上定义一层耦合面。该面是欧拉网格与拉格朗日网格之间相互作用的传递面，对于欧拉网格该面充当流场边界。同时，欧拉单元内的应力使得耦合界面上有力作用，引起拉格朗日单元发生变形。
另外，在实际建模的过程中，运用一般耦合建模时，如果计算不要求渗漏的发生，欧拉网格的密度最好大于拉格朗日网格的密度，以防止发生渗漏，影响计算。
任意拉格朗日—欧拉耦合(AL E) ，如同一般耦合一样，也可用于实现欧拉网格与拉格朗日网格之间产生相互作用。2个网格之间的耦合通过一个AL E 界面实现。与一般耦合不同的是，拉格朗日网格和欧拉网格的节点在界面上必须一一对应，位置重合，但是却是2个不同的节点，有2个不同的节点编号。也就是说，在程序允许的T o r l e rance 范围内, 在耦合界面上同一几何位置存在2 个节点。该界面一方面充当欧拉材料的流场边界; 另一方面，欧拉材料对该界面产生压力，从而使模型的拉格朗日部分受到载荷的作用。在AL E耦合的情况下欧拉网格与其他情况有所不同，一方面材料在欧拉网格中移动，另一方面，欧拉网格节点本身也在运动，使欧拉网格的位置和形状在不断调整。
1. 3 状态方程
郭庚茂简历M SC. D YTRAN 通过欧拉型的有限体积法求解流体动量方程、能量方程、连续性方程和状态方程。对于不同的流体介质(，水和空气) ，前3 种方程是相同的，但状态方程是不同的，所以要分别考虑。水和T N T材料的状态方程及参数设置如下:
( 1) 水的状态方程
水下爆炸冲击波压力一般在1 000～25 000个大
·24· 舰 船 科 学 技 术 第 28卷
介质后熵值变化很小 ，因此水的动力绝热方程为 [ 3 ]
( P 1
+ B ) ( P o
+ B ) ( 1 ) x
爱玛先声夺人x
1 o
2
2
4
式中 : P o = 1 kg / c m ;ρo = 102 kg · s /m ; 函数 x 表示
物业市场与系统的熵有关的系数 , x = 7. 15, B = 3 045 kg / c m 2 ,
水的声速为
[ 3 ]
x ( P + B
)
d P 2
。 ( 2 )
c = ρ
d ρ
因为 P o = B , , 所以式 ( 1 )可改写为水的等熵过程的状 态方程 [ 3 ]
:
图 1 水域几何模型
x
ρ P = B
- 1 ,
( 3 )
在实际建模过程中 , 把整个水域包括同时定
义为多材料的三维欧拉体 , 其相互覆盖部分的材料属
性用分级来区别 。公共区域的级别高 , 也就是说 T N T 的级别为 2, 而水单元的级别为 1。 水下冲击波
在深水部分传播时基本上不受自由
界面反射的稀疏波的影响 , 冲击波仍保持球对称峰 面 。因此为了简化建模 , 忽略水自由面对爆炸冲击的
影响 , 流场边界均定义为水介质可流出的无反射边 界 。圆柱壳的两端圆周上节点定义为简支 。 爆炸产
物以极大的速度向外扩散 , 并强烈压缩周
围的水介质 , 使水的压力呈现突跃式的升高 , 这就形
成了初始冲击波 , 随后冲击波将以极快的速度在水中 向外传播 。图 2分别显示了起爆后在 3个不同时刻 , ρo
其他参数则可按强冲击波的方法计算 。
根据 M S C . D YTR AN 程序可知 , 其对水的多项式 状态方程应为 [ 8 ]
2 3 2 3 p = α1μ +α2μ +α3μ + ( b 0 + b 1μ + b 2μ + b 3μ )ρ0 e , ( 4 )
其中 ,α1 ,α2 ,α3 , b 0 , b 1 , b 2 , b 3 均为常数 ;μ = (ρ/ρ0 ) - 1; e 为单位质量内能 ，本文中取值为 83. 950 kJ / kg [ 7 ]
。
( 2 ) TN T 的状态方程
的 引 爆 过 程 根 据 Chapm a n 2Jouge t 条 件 和 R a nk i ne 2H u gon i o t 关系来模拟 。通过起爆卡定义 起爆的时间和起爆点 ，某一单元的起爆时间则根 据该单 元 与 起 爆 点 的 距 离 和 爆 速 确 定 。在 M S C . D Y TRAN 中 ，全部引爆以后 ，其爆轰产物的压力 根据标准 Jone s - W l k i ng - L e e ( J W L )方程 [ 10 ]得到 :
T N T 球形药包爆炸生成的冲击波在水中的传播过程 。
由图 2可以看出球形药包爆炸后在水中的扩散云图 仍呈球形 。
η e η - R 1 η e η - R 2
+ωηρ e 。 ( 5 )
p = A 1 - + B 1 - 0 R 1
R 2
11
式中 : A , B , R 1 , R 2 , ω均为常数 , A = 3. 712 E Pa , B = 3. 231 E 9
Pa , R = 4. 15, R = 0. 95,ω = 0. 3; e 为单位质
1 2 量内能 , 对于 T N T 取 4. 19 MJ / k g ,η =ρ/ρ0 , T N T 3
的密度 ρ0 取 1 580 kg /m ,ρ为总体材料密度 , 爆轰传
播速度取 6 930 m / s 。
2 爆炸冲击波数值计算讨论
模拟水域长 3. 5 m ，宽 2. 0 m ，深 4 m ，采用正六
面体欧拉单元划分网格 。环向加筋圆柱壳长 2. 5 m , 半径为 0. 5 m ，壳体上有 19根肋骨 ，肋距为 0. 125 m , 采用拉格朗日四边形壳单元划分网格 , 材料为普通 钢 。圆柱壳轴线水平处于距水面 1 m 的水域中央 。
T N T 质量为球形 ，半径取 0. 14 m ，质量为 18. 16 kg，处于距圆柱壳轴线中点 2. 8 m 的水域正下方中
央 。图 1显示了在水域中剖面上水域模型及圆柱壳 和的相对几何位置 ，爆心处于坐标原点 。
图 2 冲击波传播图
2. 1 数值计算与理论计算比较
大量试验表明 , 在水深与装药半径之比大于 10 ～20 时 , 水中冲击波峰值压力不受或基本不受自由 水面和水底反射的影响 , 装药与无限水介质中爆炸时 相同或相近 , 可视为装药在无限水介质中爆炸 。根据 “爆炸相似律的分析 ”, 对于集中 T N T 球形装药在无 限水介质 中 爆 炸 水 中 冲 击 波 波 阵 面 上 的 压 力 峰 值
p w f , Co l e 经验计算公式为
ω ω
图 4 不同网格密度下冲击波时域曲线
1 / 3
α
= 53. 3 (W
/ R ) 。
( 7 )
的分布与经验值更加接近 , 吻合较好 。
2. 2 网格划分对水下爆炸冲击波的影响
理想情况下 , 网格划分的密度越大 , 单元的体积 越小 , 模型就越能真实反映物质结构 , 冲击波压力在
单元之间的传播振荡越小 , 其计算结果越精确 。但是
大量的网格单元划分 , 运行对计算机的要求很高 , 计
算时间也很漫长 。在实际计算中 , 当网格密度达到一
定值后 , 网格密度对计算结果的影响就不大了 。图 4
显示的是距爆心 0. 9 m 处 , 不同网格密度的冲击波时域曲线 。根据经验公式 , 在 0. 9 m 处冲击波的峰值应
为 179 M Pa , 从图 5中可以看出 , 网格密度越大 , 其计
课程教育研究
算结果越接近经验值 。沿圆柱壳长方向 , 网格划分密
度达到 40以后 , 加大划分密度对计算结果的影响就不大了 。网格密度取 35 时 , 计算峰值为 144 M Pa , 误差为 19. 6 % , 网 格 密 度 为 40 时 , 计 算 峰 值 为 163
M P a , 误差为 9 % , 网格密度为 50 时 , 计算峰值为 166
M Pa , 误差为 7. 3 % 。还有一个特点 , 随着网格划分密
度的增大 , 峰值压力后的压力脉动更加规则化 , 衰减 也更趋于稳定 。据经验公式可知 , 水中冲击波随时间
的衰减规律为指数形式 , 其值也是阶跃式的从 0一跃 达到峰值 , 然后以指数形式衰减 。而计算得到的曲线 在达到峰值之前有一个上升过程 , 峰值过后又以围绕
平滑的经验曲线上下波动的形式缓慢衰减 。各种外 界条件并非理想状态 , 计算冲击波曲线的这种衰减形 式更真实的反映实际情况 。直接利用 D YTR AN 计算 远场水下爆炸目前仍存在一定不足 。连续体的有限
元模型无法在数值上反映不连续区域的传播 , 当采
p w f 式中 : W 为装药量 , kg ; R 为测点波阵面距爆心的距 离 , m;α = 1. 13, 可看作是一个常数 。为得到一点处 的冲击波随时间的衰减曲线 , 则须引入以下经验公 式 :
R
R
g
σ0 p ( t) = p w f ex 。 ( 8 ) -t - t -
c
c w o w o
1 / 3
- 0. 24
1 / 3
式中 : p w f 的 单 位 为 M P a ; θ = ( W
/ R )
W
×
10 - 4
, s ; c
1 500 m / s ;
w o R
σ0 t -
c w o
( 9 )
随着 R 的增加 , 常数 θ也增加 , 冲击波的波长也
变大 。
图 3分别显示了 2 种网格划分时冲击波压力峰 值在距药包中心 0. 4 ～2. 0 m 距离内的分布情况 , 并 与经验值进行了比较 。图中实线是依据经验公式绘 出的经验曲线 。图中虚线分别显示的是沿圆柱壳长 方向划分 40 个单元和 20 个单元时的计算结果 。可 以看出 , M S C . D YTR AN 程序在整体上可以较好地模 拟爆炸冲击波压力峰值在空间上的分布 。当网格划 分密度增大时 , 冲击波压力峰值在 0. 4 ～2. 0 空间上
p
·26· 舰 船 科 学 技 术 第 28卷
用一种不具有算法阻尼的时间积分法来模拟冲击波 时 , 由于有限元网格只能产生有限的频率谱 , 因而会 产生严重的振荡伴随波阵 。为了控制冲击波的形成 , D Y TRAN 引入了人工粘性以减小伴随波阵的振荡 , 同时它使波阵处的压力随应变率增加 , 这样冲击波
扑在大约 5 个单元上 [ 8 ]
。按照这种计算方法计算 , 如果水单元划得不够密的话 , 则冲击波在水单元中 传播时峰值压力衰减会很快 。根据计算经验 , 单元 的尺寸在毫米数量级上 。这样对于远场水下爆炸而 言 , 不仅由于水单元的数目巨大 , 而在实际计算当 中几乎不可能实施 。同时计算结果也不稳定 。从理 论上来说 , 网格的划分应该是越细计算的结果就会越 准确 。但是实际计算考虑到建模计算时间和计算成 本 , 不可能满足网格足够小的要求 , 而在实际的工程 运算中 , 只要让计算结果保持在一定的精度范围内就 可以满足要求 。因此 , 在实际操作中 , 可以先建一个 有参考价值的简单模型 , 比较两种网格划分密度的计 算结果 , 如果两者之间差别不大 , 建议取较小的网格 密度再进行模型计算 。
2. 3 水域边界问题
由于计算可简化为在无限水介质中的爆炸 , 冲击 波的传播基本上不受自由界面反射的稀疏波的影响 , 而且为避免冲击波在欧拉网格内产生积压 , 忽略水自
由面对爆炸冲击的影响 , 流场边界均定义为水介质可 流出的无反射边界 。
从图 6可以看出 , 当将水域边界分别定义为有反 射边界和无反射边界两种情况时 , 计算所得的冲击波
时域曲线差别不大 , 2 条曲线几乎在相同时刻到达峰 值 ，而且峰值的数值相等 ，在爆炸初期 2条曲线几乎一 致 。但是在计算时间上 , 有反射边界设置的计算比无 反射边界设置的计算耗时多近一倍 。因此 , 为优化计 算 ，主要选用无反射边界作为计算中水域的边界条件 。
图 7 加筋圆柱壳不同位置附近冲击波时域曲线
爆面 、侧爆面和背爆面附近取点计算的结果 。由于 3 个测点距爆心的距离 不同 , 因 此达 到峰 值 的时 间 也就不同 。迎爆面距爆心最近 , 为 2. 25 m , 冲击 波
最先达到峰值 , 侧爆面和背爆面距爆心的距离分别 是 2. 94 m 和 3. 32 m , 侧爆面也比背爆面先达到峰
值 。从图可看出 , 3 条 冲击 波压 力曲 线随 时 间的 衰
减并 不是 很规 律 , 还出 现了 很 明显 的振 荡 , 这 除 了
二次加载的原因以外 , 还有加筋圆柱壳受到冲击波
作用 以后 , 由 于水 节点 与壳 节 点的 速度 不同 步 , 导
系统检测致水和壳之间出现反复的碰撞分离 , 再由于应力波 在壳结构中的传播等原因 , 致使壳附近的点冲击波
压力曲线的不规则振动 。
2. 5 一般耦合与 AL E 耦合
图 8 显示了当沿圆柱壳长方向划分 40 个单元 时距爆心 0. 9 m 处 运 用 2 种 耦 合 计 算 出 的 压 力 — 时间曲线与 经 验 公 式 计 算 结 果 的 比 较 。从 比 较 结
果可以看 出 , 当 网 格 的 划 分 密 度 和 单 元 形 式 相 同
时 , AL E 耦合和一般耦合的计算结果基本一致 。从
图中可以看 出 计 算 压 力 曲 线 值 并 不 是 突 跃 式 的 升 高 , 而 是 有一 个较 趋缓 的爬 升 过程 , 这 跟有 限 元 计 算的性质有关 , 它不能准确地反映冲击波波阵面上 的这种强间断 。在冲击波峰值以后 , 计算压力曲线 值会有较大幅度的波动 , 在有的计算中甚至会有几
个峰值 , 这可 能 是 由 于 爆 炸 余 能 的 扩 散 。当 t ≤7
θ 图 6 不同边界条件下冲击波时域曲线
2. 4 迎面冲击波与背面冲击波
图 7 中的 3 条 曲 线 分 别 是 在 加 筋 圆 柱 壳 的 迎 8
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议