1、2023年第47卷第4期20声 学 基 础声 学 基 础coustics FoundationA文献引用格式：翟宇文，马广彬，田小强，等.PML 在声波方程间断有限元计算中的应用 J.电声技术，2023，47（4）：20-22.ZHAI Y W，MA G B，TIAN X Q，et al.Application of perfectly matched layer in discontinuous galerkin method calculation for acoustic equationJ.Audio Engineering，2023，47（4）：20-22.中图分类号：P631.4
2、文献标识码：A DOI：10.16311/j.audioe.2023.04.006PML 在声波方程间断有限元计算中的应用翟宇文1，马广彬1*，田小强2，董兴蒙1（1.中国电子科技集团公司第二十二研究所，河南 新乡 453000；2.中国石油集团测井有限公司物资装备公司，重庆 401120）摘要：在模拟声波在无限空间中的传播时，标准完美匹配层（PerfectMatchedLayer，PML）对于一阶偏微分方程而言简单直接。将 PML 引入到二维声波方程的间断有限元计算中，在二维声波一阶方程组中引入 PML 边界条件，然后将带 PML 边界的方程组整理成适用于间断 Galerkin（Discon
3、tinuousGalerkin，DG）求解的双曲守恒形式，并用 DG 进行计算。数值结果表明，PML 能有效吸收介质中的声波，极大地衰减进入其中的声场能量，实现了良好的吸收效果。关键词：完美匹配层（PML）；间断有限元；声波方程Application of Perfectly Matched Layer in Discontinuous Galerkin Method Calculation for Acoustic EquationZHAIYuwen1,MAGuangbin1*,TIANXiaoqiang2,DONGXingmeng1(1.N0.22ResearchInstituteofCh
4、inaElectronicsTechnologyGroupCorporation,Xinxiang453000,China;2.CNLCMaterialsandEquipmentCompany,Chongqing401120,China)Abstract:Whensimulatingthepropagationofsoundwavesininfinitespace,thePerfectlyMatchedLayer(PML)issimpleanddirectforthefirst-orderpartialdifferentialequation.Inthispaper,PMLisintroduc
5、edintothediscontinuousfiniteelementcalculationoftwo-dimensionalacousticwaveequations,andPMLboundaryconditionsareintroducedintotwo-dimensionalacousticwavefirst-orderequations.Then,theequationswithPMLboundaryaresortedintohyperbolicconservationformssuitableDiscontinuousGalerkin(DG)solution,andcalculate
6、dbyDG.ThenumericalresultsshowthatPMLcaneffectivelyabsorbsoundwavesinthemedium,greatlyattenuatethesoundfieldenergyenteringit,andachievegoodabsorptioneffect.Keywords:PerfectMatchedLayer(PML);discontinuousfiniteelementmethod;acousticwaveequation0 引 言在数值求解声波在无限或者半无限空间中的传播问题时，通常需要截断计算区域。BERENGER1 提出了十分高效
7、的完美匹配层（Perfect Matched Layer，PML）吸收边界条件，PML 现已广泛应用在有限差分法或显式有限元法计算声波与弹性波的传播方面2。间断有限元方法最早由 Reed 等人在求解中子输运方程时提出，经过多年的发展，已经出现了多种多种多样的间断有限元方法，其中 Cockburn 和Shu 提出的求解双曲守恒方程的龙格库塔间断有限元方法的应用尤为广泛3。贺茜君将 DG 方法引入数值求解地震波方程中4，成功实现了复杂区域地震波传播数值模拟。本文推导了带 PML 声波方程的双曲守恒形式，并用 DG 方法完成声场计算。1 PML 声波方程的双曲守恒形式1.1 标准 PMLPML 吸收
8、层将波场分量在吸收边界区域分裂，分别对各个分裂的波场分量赋以不同的耗损。在计算区域截断边界外，PML 层是一种非物理的特殊吸收介质，该层的波阻抗与相邻介质的波阻抗完全匹配，因而入射波将无反射地穿过截断边界进入作者简介：翟宇文（1990），男，硕士，工程师，研究方向为波场正演、地球物理探测。通信作者：马广彬（1989），男，硕士，工程师，研究方向为地球物理探测、装备设计与研发。E-mail：。2023年第47卷第4期21Acoustics FoundatioN声 学 基 础声 学 基 础PML 层内。由于 PML 层被设置为吸收介质，进入PML 层内的入射波将被迅速衰减，最终实现削弱边界反射的效
9、果。PML 吸收边界可以看作是对复空间中实坐标的解析扩展，具体做法如图 1 所示。计算区域dx 0，dz=0dx 0，dz=0dx 0，dz=0 xz图 1 PML 吸收边界示意图在内部计算区域，采用一般的声波方程。而在PML 层区域，以 x 方向为例，定义一个新的复坐标系统5，即()()0dxxix xxdss=?（1）式中：为角频率，dx为 x 方向的阻尼因子。对上式求导，可得复拉伸函数为()()00020d11log11xzxzxxxxxnxix xxdssidxxsxxi+=|=|=+|=|=|?（2）对于 PML 吸收边界，只需在频率空间域将方程中对 x 方向的偏导数()()()()
10、00020d11log11xzxzxxxxxnxix xxdssidxxsxxixxsxhPvvctxzvPtxvPddHVdRtzH=+=|=|=+|=|=|?替换为对复坐标的偏导数()()()()00020d11log11xzxzxxxxxnxix xxdssidxxsxxixxsxhPvvctxzvPtxvPddHVdRtzH=+=|=|=+|=|=|?即可实现，即()()()()00020d11log11xzxzxxxxxnxix xxdssidxxsxxixxsx=|=|=+|=|=|?（3）定义衰减因子表达式为0020d11log11xzxzxxxxxnxhddHVdRH=+=|=
11、|=+|=|=|?（4）式中：h 为 x 从计算区域边界起进入 PML 的厚度，其取值范围为 0 H；H 为 PML 的总厚度；为微调系数，可以取 1 4；V 为声波速度；n 为阶数，通常取 2 4；R 为反射系数。当 PML 为 5 层时，R=0.01；当 PML 为 10 层时，R=0.001；当 PML为 20 层时，R=0.000 1。1.2 带 PML 的声波方程双曲守恒形式在各向同性介质中，二维声波的一阶速度-压力方程可表示为2011xzxzxxxxxnxix xxdssidxxsxxixxsxhPvvctxzvPtxvPtz=+=|=|=+|=|=|（5）式中：P 为声压，vx为
12、质点在 x 方向的振速，vz为质点在z方向的振速，c0为介质声速，为介质密度，t 为传播时间。将声压P从x方向和z方向拆分为Px和Pz，可得()()()()202020202020 xzxxzzxxzzxzxxzzxxxxxxxPPPPvctxPvctzPvctxPvctzPPvtxPPvtz=|=|+|=|+|=|（6）记=(x,z)，在 PML 层内，对式（6）进行复坐标变换，将其整理成双曲守恒形式，可以得到带 PML的声波方程式（7）。该方程可以直接进行 DG 求解。当 dx=dz=0 时式（7）退化为式（6）。TT2020()()00()00 xzxzxxzzxxzzxzxzxzftv
13、vPPfd vd vd Pd PPPPPc vc v+=|=|+|+|=|uF uu=F（7）2 数值算例建立二维数值计算模型，定义模型为 80 m 80 m 的正方形，x 方向和 z 方向的坐标范围均为-40,40。采用自由三角形划分网格，在 x 方向和 z方向上的最大网格尺寸（即 x，z）均为 2m，四周 PML 厚度为 10 m（划分 5 层网格）。PML 数值计算模型如图 2 所示。声源 S 位于模型中心，其坐标为(0,0)。选择 A点作为声场的接收测试点，A 点坐标为(-20,20)。A 点距离左边 PML 层和下边 PML 层的距离均为 2023年第47卷第4期22声 学 基 础声
14、 学 基 础coustics FoundationA10 m，虚线内部为计算区域，虚线以外为 PML 区域。当不使用 PML 边界时，虚线设定为常规的匹配阻抗边界。数值算例中，设置声速 c0=1 500 ms-1，密度=1 000 kgm-3，声源主频 f0=300 Hz。声源函数使用雷克子波，表达式为()()()TT202202220200002 1 exp()()00()001122xzxzxxzzxxzzxzxzxzftvvPPfd vd vd Pd PPPPPS tfttfc vc vvtt=+uF uu=F（8）式中：t0为雷克子波的延迟时间，通常选用一个周期的长度，即 t0=1/f
15、0。分别使用匹配阻抗边界和 PML 边界求解声场后，A 点接收的波形对比结果如图 3 所示。时间/s1.21.00.80.60.40.20-0.2-0.4-0.6-0.8-1.0PML 边界匹配阻抗边界0.020.060.040.080归一化幅度图 3 A 点接收波形对比当采用常规的匹配阻抗边界进行计算时，从 A点的接收波形中可以看出，A 点在 0.03 s、0.06 s均接收到来自边界的反射波；采用 PML 边界时，A点接收波形中没有来自边界的反射。定义声场能量为62220200002 1 exp()()00()001122vEcp=+（9）式中：v 为质点振速，p 为声压值。由式（9）可以
16、计算整个区域内的声场能量，计算结果如图 4 所示。时间/s10210110010-110-210-310-410-510-610-7PML 边界匹配阻抗边界0.020.060.040.080声场能量/J图 4 声场能量随时间的变化曲线由图 4 可知：0.02 s 时，声场还未进入 PML层，无论采用何种边界整个声场内能量均为 272 J；0.04 s 时，采用匹配阻抗边界的声场能量值为 0.4 J，而采用 PML 边界的声场能量值为 0.005 J，二者能量值相差数倍。3 结 语本文详细推导了带 PML 声波方程的双曲守恒形式，并用 DG 方法实现了对方程的数值求解。数值算例结果表明，PML 可以有效消除介质边界处发生的反射。对于 5 个网格厚度的 PML 层来说，声场传播进入 PML 层后，其能量从 272 J 迅速降低至 0.005 J。参考文献：1BERENGERJP.AperfectlymatchedlayerfortheabsorptionofelectromagneticwavesJ.JournalofComputationalPhysics，1994，114（2）：185