水下爆炸船体结构响应间断伽辽金法数值模拟

于福临，郭 君，姚熊亮，任少飞

(哈尔滨工程大学 船舶工程学院, 150001 哈尔滨)



水下爆炸船体结构响应间断伽辽金法数值模拟

于福临，郭 君，姚熊亮，任少飞

(哈尔滨工程大学 船舶工程学院, 150001 哈尔滨)

为求解水下爆炸强间断流场，采用Level Set方法定位多相流界面位置，应用虚拟流体方法处理邻近界面两侧物理量，并用RKDG方法进行空间和时间的离散，求解流场的Euler方程，并进行一维、二维评价，计算结果能够较好地反映水下爆炸冲击波产生、传播、反射和爆炸产物的膨胀等现象.最后，结合大型非线性有限元软件ABAQUS，模拟了船体板在水下爆炸载荷作用下的变形和响应特征.模拟结果表明,间断迦辽金法能够实现对水下爆炸船体结构响应精确模拟，板架结构响应与爆心距离成反比.

水下爆炸;间断伽辽金;载荷计算;结构响应;有限元;数值模拟

水面舰船在作战中，会受到水雷、鱼雷等水下武器的攻击，对舰船生命力构成威胁.国内外学者对水下爆炸载荷及结构响应方面开展了广泛的研究.试验研究是最有效的手段，国内外开展了一些实船[1]、舱段以及模型试验[2]，但水下爆炸试验耗资巨大.解析法只能对简单边界的水下爆炸载荷进行求解，文献[3-4]给出了水下爆炸载荷的经验公式，但是经验公式无法对过程进行描述，而且难以对结构响应进行求解.随着计算机性能的提高，数值技术成为研究水下爆炸的一种有效手段，发展出了有限元、有限体积、边界元和无网格等数值方法.间断伽辽金方法是介于有限元法和有限体积法之间的一种算法[5]，具备两者的优点.一方面，本质上是一种高精度的微分方程空间离散方法，在离散过程中使用弱形式积分，具有有限维解，在单元内使用多项式逼近物理量的真实值，通过简单的增加多项式次数，即可直接提高单元阶数，构造高阶格式，与有限元方法相似;另一方面，间断伽辽金方法引入数值通量的概念，在单元与单元间计入流场间断，使其可用于冲击波等强间断现象的捕捉，同时还具备多种非线性的限制器，在保证其精度的前提下，抑制非物理振荡，这与有限体积法相似.间断伽辽金法本身就具有高分辨率捕捉的特性，能有效处理水下爆炸中出现的压力、密度、速度等物理量的间断.

本文基于间断伽辽金方法，结合虚拟流体方法(ghost fluid method)和Level Set方法，编制计算程序，对一维和二维算例进行了计算，并模拟了近壁面条件下的水下爆炸流场特性，最后结合ABAQUS，对水下爆炸载荷作用下船体结构的变形和响应特征进行了模拟.

1 数值模拟方法

1.1 控制方程

对于无黏可压缩流场，欧拉方程组可表示[6]为

(1)

对于二维计算区域，式(1)可表示为

(2)

式中：ρ为流体密度；p为压力； u为沿x轴速度；v为沿y轴速度；E为单位体积总能量.

1.2 状态方程

气体采用理想气体状态方程描述[7]

(3)

对于水而言，采用Tait状态方程描述[8]为

(4)

式中：ρ为密度；B取3.31×108Pa；γ=7.15；A=105Pa；ρ0=1 000 kg/m3.

1.3 空间与时间离散

对于水下爆炸无黏可压缩性流场，在空间上采用间断伽辽金法进行离散，设流体计算区域记为Ω，在式(1)两侧同时乘以试探函数w(x)，引入流通量得[9]

(5)

(6)

式中,Pm(x)、Pn(y)分别为Legendre多项式，将式(6)带入式(5)即得到离散方程.

采用3阶TVDRunge-Kutta方法对流体力学方程组进行时间离散[11]，具体形式为

(7)

式中U(0)为tn时刻流场物理特性.

1.4 界面追踪

应用Level Set方法定位多相流界面位置，Level Set方程如下[12]

(8)

式中：V=(u,v,w)为运动速度；Φ为距离函数.

(9)

式中：Φ0为初始时刻Φ值;S(Φ0)为光滑函数.

1.5 虚拟流动方法

原始的虚拟流动方法简单的采用虚拟流体物理场的特性对应真实流体的物理场特性，方法较为简单，在处理强间断问题时，无法给出精确的界面状态[13]，误差较大，甚至无法计算.

本文使用修正的真实虚拟流动方法重构界面状态，使用两相流黎曼问题的求解方式，沿扫掠方向拾取状态参数.为了方便描述，沿x向描述这种算法.界面位于xL和xR之间，线段l1和l2代表界面.

界面法向n使用以下近似[14]:

(10)

设xL和xR处的状态参数为

(11)

速度分解成法向和切向分量为

使用投影态给出的参量求解黎曼问题为

(12)

最后，通过重新结合黎曼解得到的中间状态和切向速度得到了界面状态参数为

(13)

2 算例验证

2.1 一维算例

采用RKDG方法编制计算程序，模拟自由场空中爆炸流场，并与Henrych经验公式进行对比.表1和图1,2分别为不同爆距条件下冲击波超压峰值计算值与经验值对比结果.

表1 冲击波超压峰值计算结果

从表1可以看出，RKDG方法计算得到的冲击波超压峰值与经验公式吻合很好，最小误差仅为1.02%，最大误差不超过13%，说明本文提出的RKDG计算程序能够有效模拟空中爆炸冲击波载荷特性.

2.2 二维激波管算例

选取文献[15]中的算例，计算区域区取为[0,325]×[-44.5,44.5]，一Mach数为1.22的激波通过一个气泡，气泡中心位于(175,0)，半径25，激波位于x=225，向左传播，网格为320×80.无量纲化初始条件:气泡内ρ=0.138，u=0，v=0，p=1;激波前ρ=1，u=0，v=0，p=1;激波后ρ=1.376 4，u=-0.394，v=0，p=1.569 8.

图1 流场不同时刻密度

图2 流场不同时刻压力

由图1可以看出，平面冲击波与气泡相互作用后，气泡向左运动，当冲击波扫过气泡后，气泡的迎波面(右侧)向内凹陷，而背波面(左侧)形态基本保持不变.随着时间推移，迎波面继续向内凹陷，迎波面形成涡环状结构，这是由于冲击波诱导气泡发生射流而形成的漩涡.由图2可以看出，冲击波传播至气泡处时，气泡内部压力升高，并在气泡上下靠近壁面处形成局部高压，在t=100时，冲击波在传播方向的最前方形成高压区，随着时间推移，冲击波诱导气泡发生射流并在射流后方形成高压区.

由一维和二维算例可以看出本文的数值方法计算冲击波压力值较精确，能够实现对流场间断的高分辨率捕捉，对激波具有较高的分辨率，通过在间断附近构造Riemann问题的方式，能有效的抑制界面附近非物理震荡，并能消除间断附近的数值耗散.

3 近壁面水下爆炸冲击波载荷数值模拟

初始时刻，爆炸流场内密度为1 630 kg/m3，压力为7.8×109Pa，外部流场水的密度为1 000 kg/m3，压力为1×105Pa，爆炸初始流场半径为1 m，上边界为固壁边界条件.

初始时刻，冲击波和材料界面都迅速向外膨胀，并向周围流场中释放一道球面冲击波，称为主波，同时形成一道朝向气泡中心汇聚的内聚冲击波，如图3(a).当冲击波传播至壁面后，被刚性壁面完全反射，并仍以冲击波的形式朝向气泡传播.当冲击波传播至气泡表面后，由于气泡外部介质的阻抗高于气泡内部，因此一部分冲击波将以稀疏波的形式反射回流场中，而另一部分将以冲击波的形式透射入气泡内部，如图3(b)、3(c).由于主波相对壁面来说是以大角度入射的斜冲击波，因此开始出现马赫反射现象.主波脱离壁面表面，形成马赫杆结构，与壁面反射波一起交汇于“三点波”处，并在三波结构后出现一道滑移线.主波、壁面反射波以及马赫杆均为冲击波，波前波后的密度、压力和Ma均间断，而滑移线内侧的压力连续，密度和Ma间断，如图3(d)、(e).当稀疏波传播至壁面，即从低阻抗介质入射至高阻抗介质，又将以稀疏波的形式被反射，传播方向相反的两道稀疏波相互叠加，最终使得该区域内的压力低于临界压力以下，并形成空化，如图3(f)所示.

图3 近壁面水下爆炸流场物理特性

4 船体板架结构响应数值模拟

4.1 结构模型及工况设置

建立船体板架结构有限元模型，模型尺寸为6 m×4 m，厚度20 mm，综合考虑计算精度和计算耗时，选取网格数量120×80，编制水下爆炸载荷计算程序并结合非线性有限元软件ABAQUS，对船体板架结构响应进行数值模拟.上边界为船体板结构，板四周刚性固定.应用RKDG方法计算整个流场载荷，并与ABAQUS软件结合，计算船体板架响应特性.设板架中心的坐标为O(0 m,0 m)，选取4个考核点，坐标依次为A(2.0 m,0 m)、B(1.5 m,0 m)、C(1.0 m,0 m)、D(0.5 m,0 m).

4.2 船体板架结构毁伤变形分析

当冲击波传播作用到板架结构时，板架承受压力从0.1 MPa增大GPa至级，达到材料屈服极限，结构出现塑性变形，计算得到的船体板架结构变形如图4所示.从图4可以看出，冲击波作用范围内船体板产生塑性变形，且冲击波作用的圆形区域内部结构应变较小，结构产生圆盘状拱起变形，与板格形成塑性铰.

图4 船体板架结构塑性变形

4.3 船体板架结构响应分析

图5给出了1.2 ms时船体板架结构响应云图.图6给出了A、B、C、D这4个考核点的板架结构响应计算结果.从图6可以看出，冲击波传播过程中，位移及加速度响应大小沿板格中心向外逐渐减小，离爆心最近的D点位移、加速度响应最大，而A点最小，即板架结构响应与距离爆心距离呈反比.

5 结 论

1)能够实现对流场强间断的高分辨率捕捉，对激波具有较高的分辨率，通过在间断附近构造Riemann问题的方式，能有效的抑制界面附近非物理震荡，并能消除间断附近的数值耗散.

2)RKDG-GFM-LSM方法能很好的模拟近壁面水下爆炸冲击波产生、传播、反射和爆炸产物的膨胀等现象，并观察到冲击波传播中的马赫杆以及气泡空化现象.船体板架结构水下爆炸模拟中发现，结构产生圆盘状拱起变形，与板格形成塑性铰，板架结构响应与距离爆心距离成反比.

图5 船体板架结构响应云图

图6 考核点位移、加速度响应时历曲线

[1] SHIN Y S, SCHNEIDER N A. Ship shock trial simulation of USS Winston S. Churchill (DDG 81): modeling and simulation strategy and surrounding fluid volume effects[C]//Proceedings of the 74thShock and Vibration Symposium. San Diego, CA, Oct. : [s.n.], 2003: 27-31.

[2] 崔杰, 张阿漫, 郭君, 等. 舱段结构在气泡射流作用下的毁伤效果[J]. 爆炸与冲击, 2012, 32(4): 355-361.

[3] 崔杰. 近场水下爆炸气泡载荷及对结构毁伤试验研究[D]. 哈尔滨: 哈尔滨工程大学, 2013.

[4] SUSANTO A, IVAN L, De STERCK H, et al. High-order central ENO finite-volume scheme for ideal MHD[J]. Journal of Computational Physics, 2013, 250: 141-164.

[5] JEONG W, SEONG J. Comparison of effects on technical variances of computational fluid dynamics (CFD) software based on finite element and finite volume methods[J]. International Journal of Mechanical Sciences, 2014, 78: 19-26.

[6] QIU Jianxian, LIU Tiegang, KHOO B C. Simulations of compressible two-medium flow by runge-kutta discontinuous galerkinmethods with the ghost fluid method[J]. Communications in Computational Physics, 2008, 3(3): 479-504.[7] 张学莹，赵宁，王春武. 可压缩多介质流体数值模拟中的Level-Set间断跟踪方法[J]. 计算物理,2006, 23(5): 518-524.[8] JOHNSEN E. Numerical simulations of non-spherical bubble collapse: with applications to shockwave lithotripsy[D]. Pasadena: California Institute of Technology, 2008.

[9] 任少飞. 基于间断迦辽金法的舰船水下及空中爆炸模拟[D]. 哈尔滨：哈尔滨工程大学, 2012.

[10]ZHU Jun, LIU Tiegang, QIU Jianxian, et al. RKDG methods with WENO limiters for unsteady cavitating flow[J]. Computers & Fluids, 2012, 57: 52-65.

[11]BILLET G, RYAN J, BORREL M. A Runge Kutta Discontinuous Galerkin approach to solve reactive flows on conforming hybrid grids: the parabolic and source operators[J]. Computers & Fluids, 2014, 95: 98-115.

[12]OSHER S, FEDKIW R P. Level set methods: an overview and some recent results[J]. Journal of Computational Physics, 2001, 169(2): 463-502.

[13]陈荣三. 大密度比和大压力比可压缩流的数值计算[J]. 应用数学和力学, 2008, 29(5): 609-617.

[14]BO W, GROVE J W. A volume of fluid method based ghost fluid method for compressible multi-fluid flows[J]. Computers & Fluids, 2014, 90: 113-122.

[15]HAAS J F, STURTEVANT B. Interaction of weak shock waves with cylindrical and spherical gas inhomogeneities[J]. Journal of Fluid Mechanics, 1987, 181: 41-76.

(编辑 张 红)

Numerical simulation of the response for Hull plates subjected to underwater explosion based on RKDG method

YU Fulin, GUO Jun, YAO Xiongliang, REN Shaofei

(College of Shipbuilding Engineering, Harbin Engineering University, 150001 Harbin，China)

In order to solve the underwater explosion flow field with large discontinuities, Level Set method was applied to track the interface position of the multi-medium flow, Ghost Fluid method was used to calculate the physical parameter of both sides of the interface, time and space were discretized by Runge-Kutta Discontinuous Galerkin Method, Euler equations of the flow field were solved. One-dimensional and two-dimensional assessments were conducted by RKDG approach. The results reflect the phenomena of underwater explosion shock wave generation, propagation, reflection and explosion products expansion. Finally, the shock responses and damage characteristics of hull plates under shock load were simulated with the nonlinear FEM softeware ABAQUS. The RKDG method can be applied to simulate the hull plates response with high accuracy. The response of hull plates is inversely proportional to the blast center distance.

underwater explosion; RKDG; shock load calculation; hull plates response; FEM; numerical simulation

10.11918/j.issn.0367-6234.2016.04.023

2014-11-26.

国家自然科学基金(51109042)；中国博士后科学基金(2012M520707)；黑龙江省自然科学基金 (E201124).

于福临(1988—)，男，博士研究生；

姚熊亮(1963—)，男，教授，博士生导师.

姚熊亮，yaoxiongliang@163.com.

U663.2

A

0367-6234(2016)04-0139-05