水下爆炸船体结构响应间断伽辽金法数值模拟
于福临;郭君;姚熊亮;任少飞
【摘 要】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.%为求解水下爆炸强间断流场，采用Level Set方法定位多相流界面位置，应用虚拟流体方法处理邻近界面两
侧物理量，并用RKDG方法进行空间和时间的离散，求解流场的Euler方程，并进行一维、二维评价，计算结果能够较好地反映水下爆炸冲击波产生、传播、反射和爆炸产物的膨胀等现象。最后，结合大型非线性有限元软件ABAQUS，模拟了船体板在水下爆炸载荷作用下的变形和响应特征。模拟结果表明，间断迦辽金法能够实现对水下爆炸船体结构响应精确模拟，板架结构响应与爆心距离成反比。
【期刊名称】《哈尔滨工业大学学报》
【年(卷),期】2016(048)004
【总页数】5页(P139-143)
【关键词】水下爆炸;间断伽辽金;载荷计算;结构响应;有限元;数值模拟
【作 者】于福临;郭君;姚熊亮;任少飞
【作者单位】哈尔滨工程大学 船舶工程学院，150001 哈尔滨;哈尔滨工程大学 船舶工程学院，150001 哈尔滨;哈尔滨工程大学 船舶工程学院，150001 哈尔滨;哈尔滨工程大学 船舶工程学院，150001 哈尔滨
【正文语种】中 文
【中图分类】U663.2
水面舰船在作战中，会受到、等水下武器的攻击，对舰船生命力构成威胁.国内外学者对水下爆炸载荷及结构响应方面开展了广泛的研究.试验研究是最有效的手段，国内外开展了一些实船[1]、舱段以及模型试验[2]，但水下爆炸试验耗资巨大.解析法只能对简单边界的水下爆炸载荷进行求解，文献[3-4]给出了水下爆炸载荷的经验公式，但是经验公式无法对过程进行描述，而且难以对结构响应进行求解.随着计算机性能的提高，数值技术成为研究水下爆炸的一种有效手段，发展出了有限元、有限体积、边界元和无网格等数值方法.间断伽辽金方法是介于有限元法和有限体积法之间的一种算法[5]，具备两者的优点.一方面，本质上是一种高精度的微分方程空间离散方法，在离散过程中使用弱形式积分，具有有限维解，在单元内使用多项式逼近物理量的真实值，通过简单的增加多项式次数，即可直接提高单元阶数，构造高阶格式，与有限元方法相似;另一方面，间断伽辽金方法引入数值通量的概念，在单元与单元间计入流场间断，使其可用于冲击波等强间断现象的捕捉，同时还具备多种非线性的限制器，在保证其精度的前提下，抑制非物理振荡，这与有限体积法相似.间断伽辽金法本身
就具有高分辨率捕捉的特性，能有效处理水下爆炸中出现的压力、密度、速度等物理量的间断.
本文基于间断伽辽金方法，结合虚拟流体方法(ghost fluid method)和Level Set方法，编制计算程序，对一维和二维算例进行了计算，并模拟了近壁面条件下的水下爆炸流场特性，最后结合ABAQUS，对水下爆炸载荷作用下船体结构的变形和响应特征进行了模拟.
1.1 控制方程
对于无黏可压缩流场，欧拉方程组可表示[6]为
对于二维计算区域，式(1)可表示为
式中：ρ为流体密度；p为压力； u为沿x轴速度；v为沿y轴速度；E为单位体积总能量.
1.2 状态方程
气体采用理想气体状态方程描述[7]
对于水而言，采用Tait状态方程描述[8]为
式中：ρ为密度；B取3.31×108 Pa；γ=7.15；A=105 Pa；ρ0=1 000 kg/m3.
1.3 空间与时间离散
对于水下爆炸无黏可压缩性流场，在空间上采用间断伽辽金法进行离散，设流体计算区域记为Ω，在式(1)两侧同时乘以试探函数w(x)，引入流通量得[9]
式中：=[n1,n2,n3]T为单元的外法向量；；Γj、Ωj分别为在单元线上积分和单元面上积分.
对于二维问题，i,j(x,y,t)使用如下近似[10]
式中,Pm(x)、Pn(y)分别为Legendre多项式，将式(6)带入式(5)即得到离散方程.
采用3阶TVD Runge-Kutta方法对流体力学方程组进行时间离散[11]，具体形式为
式中U(0)为tn时刻流场物理特性.
1.4 界面追踪
应用Level Set方法定位多相流界面位置，Level Set方程如下[12]
式中：V=(u,v,w)为运动速度；Φ为距离函数.
为了保证距离函数的性质，需要对Level Set方程时刻更新，使=1，初始化方程如下
式中：Φ0为初始时刻Φ值;S(Φ0)为光滑函数.
1.5 虚拟流动方法
原始的虚拟流动方法简单的采用虚拟流体物理场的特性对应真实流体的物理场特性，方法较为简单，在处理强间断问题时，无法给出精确的界面状态[13]，误差较大，甚至无法计算.
本文使用修正的真实虚拟流动方法重构界面状态，使用两相流黎曼问题的求解方式，沿扫掠方向拾取状态参数.为了方便描述，沿x向描述这种算法.界面位于xL和xR之间，线段l1和l2代表界面.
界面法向n使用以下近似[14]:
式中分别为线段l1和l2的长度.
设xL和xR处的状态参数为
速度分解成法向和切向分量为
使用投影态给出的参量求解黎曼问题为
最后，通过重新结合黎曼解得到的中间状态和切向速度得到了界面状态参数为
式中，下标I表示黎曼解得到的中间状态.xL和xR之间界面的数值通量通过计算得到.
2.1 一维算例
采用RKDG方法编制计算程序，模拟自由场空中爆炸流场，并与Henrych经验公式进行对比.表1和图1,2分别为不同爆距条件下冲击波超压峰值计算值与经验值对比结果.
从表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.reactive近场区
由图1可以看出，平面冲击波与气泡相互作用后，气泡向左运动，当冲击波扫过气泡后，气泡的迎波面(右侧)向内凹陷，而背波面(左侧)形态基本保持不变.随着时间推移，迎波面继续向内凹陷，迎波面形成涡环状结构，这是由于冲击波诱导气泡发生射流而形成的漩涡.由图2可以看出，冲击波传播至气泡处时，气泡内部压力升高，并在气泡上下靠近壁面处形成局部高压，在t=100时，冲击波在传播方向的最前方形成高压区，随着时间推移，冲击波诱导气泡发生射流并在射流后方形成高压区.
由一维和二维算例可以看出本文的数值方法计算冲击波压力值较精确，能够实现对流场间断的高分辨率捕捉，对激波具有较高的分辨率，通过在间断附近构造Riemann问题的方式，能有效的抑制界面附近非物理震荡，并能消除间断附近的数值耗散.
初始时刻，爆炸流场内密度为1 630 kg/m3，压力为7.8×109 Pa，外部流场水的密度为1 000
kg/m3，压力为1×105 Pa，爆炸初始流场半径为1 m，上边界为固壁边界条件.
初始时刻，冲击波和材料界面都迅速向外膨胀，并向周围流场中释放一道球面冲击波，称为主波，同时形成一道朝向气泡中心汇聚的内聚冲击波，如图3(a).当冲击波传播至壁面后，被刚性壁面完全反射，并仍以冲击波的形式朝向气泡传播.当冲击波传播至气泡表面后，由于气泡外部介质的阻抗高于气泡内部，因此一部分冲击波将以稀疏波的形式反射回流场中，而另一部分将以冲击波的形式透射入气泡内部，如图3(b)、3(c).由于主波相对壁面来说是以大角度入射的斜冲击波，因此开始出现马赫反射现象.主波脱离壁面表面，形成马赫杆结构，与壁面反射波一起交汇于“三点波”处，并在三波结构后出现一道滑移线.主波、壁面反射波以及马赫杆均为冲击波，波前波后的密度、压力和Ma均间断，而滑移线内侧的压力连续，密度和Ma间断，如图3(d)、(e).当稀疏波传播至壁面，即从低阻抗介质入射至高阻抗介质，又将以稀疏波的形式被反射，传播方向相反的两道稀疏波相互叠加，最终使得该区域内的压力低于临界压力以下，并形成空化，如图3(f)所示.
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).

版权声明：本站内容均来自互联网，仅供演示用，请勿用于商业和其他非法用途。如果侵犯了您的权益请与我们联系QQ:729038198，我们将在24小时内删除。