水下声波传播过程的时域间断Galerkin有限元方法

李瑞敏, 郭 攀*, 张景飞, 武文华, 王 飞

(1.郑州大学 力学与安全工程学院，郑州 450000；2.大连理工大学 运载工程与力学学部 工程力学系，工业装备结构分析国家重点实验室，大连 116024)

1 引 言

由于声压力波水下良好的传播性能，声波技术广泛应用于海洋环境监测、海洋军事、海洋航运、海洋科学等领域[1]。尤其是在水下发射、水下爆炸、水下爆破和水下气泡射流等问题中，近场强间断、高梯度复杂多变的冲击声压力波往往引起结构和舰船的毁伤[2]，使得冲击条件下声压力场结构的准确描述问题成为海洋工程和军事等领域一直以来的研究热点[3,4]；另外，梯度变化声速的海洋声波传播问题一直以来是海洋声学研究的重点[5-7]；复杂的海底地形边界条件下，海底地质勘探、障碍物的探测识别和目标定位等技术领域的研究近年来从未间断[8-12]。不难发现，这些问题往往伴随着空间或时域高梯度、强间断变化的特征。

目前，随着计算机软硬件的发展，在声场的研究领域出现了诸多数值计算方法，如有限体积方法[13]、有限元法、有限差分法[14]、边界元法[15]以及相关的混合方法等。其中有限元方法在工程领域中应用最多。以往研究表明，在模拟冲击作用下间断压力波的传播过程时，通常的时域连续Galerkin有限元方法如Newmark法会在压力波波后出现虚假的数值振荡，误差的累积使得问题求解的精度大幅降低[16]。时域间断Galerkin有限元方法 DGFEM(Discontinuous Galerkin Finite Element Method)方法在冲击问题中体现了区别于传统有限元法的优越性，其可以自动引入数值耗散，滤去虚假的高阶模式和数值振荡。Li等[17,18]在原有DGFEM基础上，分别采用三次Hermit插值函数和线性插值函数对基本未知量及其导数进行插值，将基本未知量及其导数解耦，大幅提高了计算效率。但是在高频振动的波传播问题仿真中，这类方法在波前仍有较小的数值振荡现象，Guo等[19]针对这种现象进行了进一步的研究，引入高频波动比较敏感的比例刚度阻尼，并将改进的DGFEM运用到了诸多领域，在层合压电材料冲击问题、非傅里叶热弹性问题及弹性固体冲击问题等算例的求解中均表现出了该方法的优越性[20-22]。

本文在以往研究基础上，基于声压波动传播方程，构建出水下声波传播问题的改进DGFEM(以下简称DGFEM)，针对几类典型的水下高梯度和强间断问题进行了探讨。

2 声传播控制方程及时域间Galerkin有限元方法

声传播方程为

(1)

式中p为流体压力，c0为流体中的声速。

对于声波传播方程，流体内单元分布可表示为

(2)

式中mf为流体单元的节点数，pe为单元的结构点压力向量，Ni为对应节点i的插值函数。

(3)

(4)

(5)

考虑到声源和水中阻尼的作用，式(5)可以略去下标f，进而推广为

(6)

式中C为阻尼矩阵，Q(t)为各个节点上的声源向量。声波传播方程的时域可以离散为

0<…

对任一时刻tn，其时刻处的阶跃函数可以表示为

(7)

(8)

(9)

(10)

基本未知函数的导数v在时间步内采用线性插值：

(11)

省略上标，式(11)可简化为

v=vnλ1+vn + 1λ2

(12)

(13)

式(13)分别对基本未知量及其时域一阶导数进行变分，可得

(14)

式(14)解耦后，可得

(15)

改进的时域间断Galerkin有限元方法是为了消除DGFEM方法波前虚假数值振荡。其基于高频振荡敏感于比例刚度阻尼和低频振荡敏感于比例质量阻尼的原理，在式(15)引入比例刚度阻尼Ck=βK，对原有DGFEM方法进行进一步发展。

其中β为阻尼系数:

β=2ξ/ωn

(16)

ξ为结构的阻尼比，ωn为结构的固有频率。

式(15)可以改写为

(17)

3 算例及结果分析

算例1选取结构长度为50 m，截面积为0.5 m2的一维流体区域，如图1所示。在区域最左端为完全反射边界条件，最右端受到脉冲阶跃荷载的作用，如图2所示。流体介质密度为1000 kg/m3，声波在该流体中的传播速度为1500 m/s。结构划分为100个单元，时间步为0.000 1 s，分别采用 Newmark 方法与改进的DGFEM进行模拟。

图1 流体区域受冲击作用模型

Fig.1 Model of fluid region by impact action

图2 外载阶跃荷载

图3给出了传统的Newmark方法与改进的DGFEM方法在0.005 s，0.01 s和0.02 s时水中压力波分布的情况。可以看出，压力波从右向左传播，传统的Newmark法计算得到的压力波在波后和波前均出现了虚假的数值振荡，而且从波前和波后数值振荡的区域分布来看，这类数值振荡随着时间增长数值振荡范围不断增加，误差累积不断增大，这势必降低问题求解的准确性和稳定性。而改进的DGFEM方法可以很好地滤去不同时刻波前和波后虚假数值振荡，相对较好地保证了这类间断声波传播问题求解的精度和稳定性。

算例2为了验证新方法在复杂边界问题求解中的可靠性，考虑到水下障碍物的声探测技术，设计了如图4所示的模型。随机分布的孔表示水中的障碍物，设定为全反射边界。在数值模拟中，流体的密度为1000 kg/m3，声音在流体中的传播速度为1500 m/s，模拟区域的尺寸为20 m×20 m，将模型划分为1681个单元，时间步长为1.0×10-4s。声源位于点O(10.8 m,10.8 m)，声压力荷载如图5所示。

图3 不同时刻下压力波传播情况

图4 含障碍物二维声传播模型网格划分

图5 外载阶跃荷载

图6和图7是在t=0.002 s，0.005 s，0.008 s，0.01 s四个不同时刻，DGFEM和Newmark方法所模拟的流体区域中声压分布情况。可以看出，声源激发声波后，向周围传播，当它们接触障碍物时，声波发生反射并继续传播。当采用Newmark方法进行计算时，数值振荡较为明显，区域内压力碎片状分布相对无序。压力波传播过程势必会受到大型虚假数值振荡的干扰，而压力信号真实分布特征受到扰乱。使用改进的DGEFM进行模拟时，压力波整个时程清晰明确。从云图可以明显观察到压力波扩散传播过程中反射和叠加等特征。

图6 不同时刻声压分布情况(Newmark)

图7 不同时刻声压分布情况(DGFEM)

图8 层合液体网格划分

图9比较了Newmark方法和改进的DGFEM在t=0.003 s，t=0.006 s，t=0.01 s，t=0.013 s，t=0.016 s，t=0.02 s时的声波传播情况。可以看出，Newmark方法得到的结果随着压力波的传播出现了数值振荡，当压力波到达层合界面处时，层合界面处表现出更强的数值振荡，随着压力波的传播，这种数值振荡不断增强，误差累积增大。而改进的DGFEM可以很好滤除具有间断特征的压力波传播过程以及层合介质间断面两侧出现的数值振荡，三维压力波梯度结构明确，压力波间断特征较好描述。

4 结 论

本文构建了水下声波传播问题的改进时域间断Galerkin有限元方法，并针对具有强间断、高梯度变化特征的水下含多障碍物复杂边界问题以及层合液体问题进行了求解。结果表明，通常的时域连续Galerkin有限元方法如Newmark方法求解这类问题时会在波前波后出现虚假的数值振荡现象，随着波的传播，数值振荡范围会不断扩大，误差累积也不断增大，影响问题计算精度。而改进的DGFEM能较好地消除声波在不同时刻以及不同介质层合间断面处波前波后的数值振荡，使整个传播过程清晰明确，具有较高的精度和可靠性。问题的求解，为进一步流固声耦合问题研究奠定了基础。

图9 不同时刻声压分布情况