水面舰船在作战中，会受到水雷、鱼雷等水下武器的攻击，对舰船生命力构成威胁.国内外学者对水下爆炸载荷及结构响应方面开展了广泛的研究.试验研究是最有效的手段，国内外开展了一些实船^[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×10⁸ Pa；γ=7.15；A=10⁵ Pa；ρ₀=1 000 kg/m³.

1.3 空间与时间离散

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

(5)

式中： =[n₁, n₂, n₃]^T为单元的外法向量；；Γ_j、Ω_j分别为在单元线上积分和单元面上积分.

对于二维问题，使用如下近似^[10]

(6)

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

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

(7)

式中U⁽⁰⁾为t_n时刻流场物理特性.

1.4 界面追踪

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

(8)

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

为了保证距离函数的性质，需要对Level Set方程时刻更新，使|▽Φ| =1，初始化方程如下

(9)

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

1.5 虚拟流动方法

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

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

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

(10)

式中|l₁|、|l₂|分别为线段l₁和l₂的长度.

设x_L和x_R处的状态参数为

(11)

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

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

(12)

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

(13)

式中，下标I表示黎曼解得到的中间状态.x_L和x_R之间界面的数值通量通过V_I^{L, R}计算得到.

2 算例验证 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.

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

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

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

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

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

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.3 船体板架结构响应分析

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

5 结论

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

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