透射边界条件在波动谱元模拟中的实现：二维波动1)

邢浩洁 李鸿晶

(南京工业大学土木工程学院，南京211816)

固体力学

透射边界条件在波动谱元模拟中的实现：二维波动1)

邢浩洁 李鸿晶2)

(南京工业大学土木工程学院，南京211816)

将邢浩洁和李鸿晶提出的多次透射公式(multi-transmitting formula,MTF)的谱元格式应用于均匀介质中线弹性SH波动问题的谱元模拟.假定紧邻人工边界的一层谱单元为具有直线边界的四边形单元，以保证每个人工边界节点都唯一对应一条指向内域的离散网格线.人工边界节点在某时刻的位移由该离散网格线上的节点在前若干时刻的位移确定，按照MTF谱元格式进行计算.通过平面波以一定角度传播的外源问题算例和点源脉冲自由扩散的内源问题算例，验证了方法的可行性以及对实际复杂波动问题的适用性.通过不同类型初值问题算例，在时域内分析了插值多项式阶次、人工波速和透射阶次三个参数对反射误差的影响.结果表明：插值多项式阶次较高的格式会表现出更好的精度，但总体上对反射误差的影响较小；人工波速对反射误差具有显著影响，当人工波速小于介质物理波速时反射误差较大，而当人工波速等于或稍大于介质物理波速时反射误差处于较低水平；透射阶次对反射误差具有决定性影响，表现在不失稳的情形下提高透射阶次能够迅速降低反射误差，但内源问题从三阶MTF开始出现飘移失稳，外源问题从二阶MTF开始出现轻微的飘移失稳.

多次透射公式,谱元法,SH波动,MTF谱元格式,飘移失稳

引言

局部场地的地震波动数值模拟是地震学和地震工程等领域的重要课题，这类复杂介质中的大规模开放系统问题，对计算模型的空间离散以及人工边界条件都有很高要求.谱元法[2](spectral element method,SEM)兼具了有限元法的灵活性和谱方法的高精度，具有高效处理复杂几何或物理模型的能力，成为近年来备受关注的空间离散方法.从20世纪90年代开始，谱元法被用于求解地震波传播问题，在地球介质中的声波和弹性波模拟方面取得了进展[310].近十年来对谱元法的研究逐渐增多，在谱元法基本理论[11]、大规模地震动模拟[1214]、三角网格划分[1518]、并行计算[19]、非线性分析[20]以及与有限元法的对比[21]等方面都有成果发表.可以说，在处理大规模复杂波动问题方面，谱元法比传统有限差分法[22]或有限元法[23]具有显著优势.但由于谱元法在有限单元上采用了基于不等距节点的高阶插值格式，传统的以有限元法等采用低阶插值格式构造的人工边界条件计算公式在波动谱元模型中无法被直接使用，因此需要专门研究适合谱元模型的人工边界条件的具体实现方式.

Clayton-Engquist(CE)边界[24]是早期在谱元法地震波动模拟中使用最多的人工边界条件，它与谱元法结合的方式有两种，一种需将边界表达式写成沿边界的等效积分并纳入到谱元法的等效积分方程中[4]，另一种则将边界表达式写成关于边界上黏性阻尼力以及法向、切向应力的形式[79].作为一种时、空局部化的近似边界，它仅能吸收传播方向与边界法线夹角较小的外行波动.全局化的DtN边界在谱元法中有一定应用[2527]，该边界由波传播的远场解得到，是一种精确的人工边界.但全局边界的计算量很大，许多实际波动问题的远场解比较复杂，而在数值计算中需要对级数表达式进行截断[28]，这些都限制了该边界的推广.完美匹配层(perfectly matched layer,PML)边界是目前波动谱元模拟中比较热门的人工边界[2933]，其优点在于外行波在内域和吸收层的交界面上不发生反射，而且能够迅速地被吸收层内的衰减函数和阻尼作用消减掉.但无论是经典PML，还是复频移PML，形式都比较复杂，而且数值离散后的PML仍然存在精度和稳定性问题.最近，多次透射公式(multi-transmitting formula,MTF)[3435]被尝试引入波动谱元模拟中，文献[36]从时间回退和空间回退两个角度，给出了在谱元法中使用多次透射公式的方法并探讨了其稳定性问题.

MTF通过透射阶次和人工波速两个参数控制透射精度，谱元法则通过单元阶次和单元尺度来模拟复杂波场，兼具高精度和灵活性是二者的相似之处.因此，在波动谱元模拟中使用MTF作为人工边界条件，可能是提高波动模拟精度和计算效率的重要途径之一.MTF易于与内域离散格式结合，但关于这方面的研究不多.目前在波动有限差分或有限元模拟中广泛使用的是一种基于三点抛物线内插和齐次内插递推过程的MTF数值格式[37].在文献[1]中，我们详细探讨了在谱元离散网格中施加MTF时所面临的新问题，提出了一组基于拉格朗日多项式插值和简单内插方案的MTF谱元格式，并在一维波动问题中研究了这组格式的精度和稳定性特点.本文将根据文献[1]提出的MTF谱元格式，进一步探讨在二维波动谱元模拟中实现MTF的具体方法，分析方法的可行性以及复杂波动问题的物理特征与MTF计算参数之间的联系，并初步探讨与之相关的稳定性问题.

1 二维波动谱元模型

以SH型波动问题说明谱元法的基本公式以及MTF的实现过程.SH波动方程为如下的二维标量波方程

其中,u(x,y,t)为出平面波动位移，f(x,y,t)为外荷载，c为介质物理波速.模型所在区域为(x,y)∈Ω，计算时间t∈(0,T].

采用具有集中质量矩阵的高精度勒让德谱元法(Legendre spectralelementmethod,LSEM)对波动方程进行空间离散，通过以下3个步骤：(1)导出原方程的等效积分“弱”形式；(2)将整个计算区域的积分划分为各个单元积分的叠加；(3)将各个单元的积分和微分转换到标准的参考单元上.得到波动方程的等价形式如下

其中，ue(ξ,η),ve(ξ,η)分别为参考单元上的波动位移函数和测试函数，顶标“¨”表示对时间的二阶偏导数，上标T表示向量或矩阵的转置.ne为模型的单元数；Λ为参考单元，ξ,η为参考单元坐标；|J|和J-1为雅可比行列式以及雅可比矩阵的逆阵.

参考单元上的波动位移函数，即LSEM的单元位移模式由下式定义

式中，单元的总体节点编号i以及ξ,η方向的局部节点编号j,k如图1所示.单元节点数ng=nξ×nη，这里nξ,nη分别为ξ,η方向的节点数.

图1 参考单元的节点编号Fig.1 Nodenumbering of a referenceelement

勒让德基函数为定义在GLL(Gauss Lobatto Legendre)节点上的拉格朗日插值函数，具有表达式

式中，GLL节点ξi,ξj是数值分析和工程计算中的一类重要节点，很多文献都给出了其数值，这里不再赘述.η方向的勒让德基函数定义与上式相同.

方程(2)再经过如下步骤：利用伽辽金原理，选取与波动位移函数相同的测试函数；将波动位移函数、测试函数、雅可比行列式以及雅可比逆矩阵代入方程(2)，根据二维GLL数值积分计算各个积分项，得到各个单元的运动方程；将单元方程组装为总体方程.得到空间离散后的谱元节点运动方程

式中，M为总体质量矩阵,K为总体刚度矩阵，F为外载荷向量，u为被求的节点位移.

这个运动方程是关于模型所有内域节点和自由边界节点的(谱元法基本公式已将自由边界条件包含在内)，再引入物理边界条件和人工边界条件，就能通过时域逐步积分对上述方程进行求解.其中，采用适当的人工边界条件，降低或消除人工边界反射波的干扰，是保证波动模拟精度的关键.

2 多次透射公式

有限模型人工边界节点的位移由MTF进行计算.MTF是一种直接以离散形式给出的一维化的时空外推公式，它由下式定义

由于MTF计算点与内域离散点位置不一定重合，上式无法直接用于波动数值计算.只需通过适当的插值方法，由临近内域离散点位移插值得到各个MTF计算点位移，就能得到适用于波动数值模拟的MTF数值格式.

2.1 人工边界附近的谱单元类型

在波动有限元模拟中使用MTF时，通常要求每个边界节点的内法线上有一系列等距分布的有限元节点，然后利用这些节点进行三点抛物线内插和齐次内插递推方案，得到实用的MTF有限元格式[3435].其中包含两个要点，一是边界法线，二是等距分布的有限元节点.然而研究发现：“法线”的条件不必严格满足，但它必须是一条将人工边界节点和一组内域有限元节点串连起来的直线，这条直线通常是从边界节点指向内域的离散网格线；等距分布的有限元节点是保证齐次内插过程得以实现的条件，因此在谱元不等距节点分布的条件下无法使用齐次内插方案.

为了在谱元离散网格中实现一维化的MTF，首先需要确保每个人工边界节点都对应一条上述“直线”.如图2所示的三角形谱单元和曲边谱单元显然不满足这一要求，因为前者有部分边界节点上没有指向内域的离散网格线，而后者有大量的离散网格线为曲线.

如图3所示的具有直线边界的四边形谱单元能够满足上述要求.同时由于谱单元尺寸较大，对于常用阶次的MTF，其计算点一般不会超过一个谱单元范围.于是对边界附近的谱单元类型作如下假定：紧邻人工边界的一层谱单元为具有直线边界的四边形单元.

图2 三角形谱单元和曲边谱单元Fig.2 Triangleand curvilinear spectralelements

图3 直边谱单元中的MTF计算点Fig.3 Computation nodesofMTF in a spectralelementw ith straightedges

基于该假定，每个人工边界节点的MTF计算点都位于一条指向内域的离散网格线上，该点位移可由文献[1]给出的MTF谱元格式进行计算.

2.2 边界位移计算公式

以图3中A点为例来说明在二维谱元离散网格中施加MTF的具体过程，其他人工边界节点与之类似.A点指向内域的离散网格线为一条直线，线上有不均匀分布的谱元节点A,A1,A2,···,A′，以及均匀分布的MTF计算点1a,2a,···,Na.各个MTF计算点分别位于不同的时刻，如图4所示.

根据文献[1]提出的MTF谱元格式，对各个MTF计算点的位移均采用同样的拉格朗日多项式插值公式进行计算，得到A点在p+1时刻的位移计算公式为

图4 计算A点位移的MTF谱元格式Fig.4 SpectralelementschemeofMTF for computing the displacement of node A

式中，上标T表示向量转置，M为插值多项式阶次.插值系数tj,0,tj,1,tj,2,···,tj,M与MTF计算点ja(j=1,2,···,N)和谱元网格点0,1,2,···,M一一对应.s0,s1,···,sM分别为节点A与A,A与A1,···,A与AM之间的距离.

插值多项式阶次M的取值范围为2,3,···,NEn；NEn为MTF计算方向的谱单元阶次，因此式(8)～式(11)表示的是一系列MTF谱元格式.不同格式的区别在于采用的插值多项式阶次M不同，即涉及的谱元离散点数目不同.其中，M=2时为三点抛物线内插格式，M=NEn时为与谱单元位移模式一致的格式.式(8)～式(11)表示的MTF谱元格式能够实现的MTF阶次范围，由保证各个插值均为“内插”的条件确定，为

可以看出，插值多项式阶次M越高，能够实现的MTF阶次N越高.对于最低的插值多项式阶次M=2，通常N取为2和3，所以上式对常用的MTF阶次一般不构成限制.只有在少数极端情形下，如人工波速大大超过物理波速或透射阶次很高时，才需要验证.

上述MTF谱元格式是在较为稀疏的不等距谱元节点上进行插值，用于计算各个MTF计算点1a,2a,···,Na位移的谱元节点数均为M+1个；而传统MTF有限元格式是在相对密集的等距有限元节点上进行插值，用于计算各个MTF计算点1a,2a,···,Na位移的有限元节点数分别为3,5,7,···,2N+1个.二者的主要区别在于实现高阶MTF的方式不同，以及插值多项式阶次不同.MTF有限元格式实现高阶MTF的方式为齐次内插递推过程，这样的好处是反射系数表达式为R=-fN，f为一阶MTF的反射因子.此时各阶MTF的稳定界限相同，且在不失稳的情况下，随着透射阶次N的增加，反射系数能够稳步地减小.MTF谱元格式使用简单内插方案实现高阶MTF，其模拟效果比齐次内插方案有所降低，但是谱单元的高精度特性能够较好地弥补这方面不足，而高阶谱单元下插值多项式阶次M的取值则是MTF谱元格式面临的新问题，对此文献[1]已有初步讨论，本文还将进一步研究.

需要指出的是，上述MTF谱元格式和传统MTF有限元格式都来源于MTF定义式(7)，因此它们的模拟效果首先是由人工波速ca和透射阶次N这两个基本参数决定的.

3 方法可行性的验证

分别以无限均匀介质中的外源和内源SH波动问题为算例，验证本文方法的可行性.计算模型如图5所示，模型所在区域为宽600m、高400m的矩形，介质物理波速c=400m/s，模型四边均为MTF人工边界.外源问题的输入波为从模型左下方斜入射的平面波，斜入射方向与竖轴夹角为α.内源问题的输入波为作用在模型内部某点的一个脉冲波源，波源位置距模型左侧、下侧边界各100m.

输入波位移时程采用Ricker型地震子波，具有表达式

式中，f为中心频率，t0为脉冲波形的起始时间.由此不难得出脉冲宽度为2/f，峰值频率约为3f.计算时取f=10Hz，对应峰值频率约为30Hz.对于内源问题，上式为点波源的输入波位移时程；对于外源问题，上式仅为模型左下角点的输入波位移时程，模型左侧、下侧边界上其他节点的输入波位移时程，需在上式的基础上考虑相应的滞后时间.

图5 外源和内源SH波动问题计算模型Fig.5 Computationmodelof SH-typewavemotion problem w ith outer or innerwave source

采用四阶勒让德谱单元对模型进行空间离散，单元节点数为25，模型的单元数为1247(43×29=1247).两个方向的谱单元尺寸分别为13.95m和13.79m，接近于λmin=400/30≈13.33m，符合谱元离散的精度要求.时间步长取为0.003s，此时，与时域积分稳定性有关的库朗数(Courantnumber)为q=c∆t/∆xmin=400×0.003/2.38≈0.50，小于临界值(略小于1的某个常数)，符合稳定性要求.

3.1 外源问题

首先进行外源问题的模拟.平面波入射角度设置为α=30°，此时外行波在上边界和右边界的透射角度(相对于边界法线)分别为θ=30°和θ=60°.边界条件采用式(8)～式(11)表示的MTF谱元格式，插值多项式阶次取为M=4(与谱单元阶次一致).采用不同边界参数进行数值模拟，选取其中比较有代表性的两组结果，它们在0.3s,0.6s,1.1s,1.4s的波场快照如图6所示.

图6(a)为第一组模拟结果，使用的边界参数为：二阶MTF，人工波速ca取为介质物理波速c，这是MTF最常用的参数设置.前两幅子图显示，在模型的左侧和底部顺利地实现了平面波的输入，此时没有外行波，边界上不会发生反射.后两幅子图显示，当入射波到达右边界和上边界时，MTF边界能够将绝大部分外行波动能量透射出去，仅有很小一部分能量被反射回计算区域.从图中还可以看出，右边界的反射明显大于上边界的反射，其原因在于右边界为大角度透射(θ=60°)，而上边界的透射角度相对要小得多(θ=30°).这组模拟表明了本文提出的在二维谱元离散网格中使用MTF作为人工边界条件的方法是可行的.对于透射角度不大的外行波，常用的MTF参数设置就能取得较好的模拟效果；对于透射角度较大的外行波，模拟效果有所降低.

图6外源问题波场快照Fig.6 Wave fiel snapshotsof outer-source problem

图6 (b)为第二组模拟结果，使用的边界参数为：一阶MTF，人工波速ca取为各边的法向透射速度c/cosθ，即左、右边界ca=2c，上、下边界前两幅子图显示平面波能够在模型的左侧和底部顺利地输入计算区域，与上一组相同.后两幅子图则显示当入射波到达模型右边界和上边界时，能够完全透射出去，不发生任何反射.这种近乎完美的模拟结果表明,通过调整人工波速取值，可以很好地克服仅由透射角度带来的影响.这组模拟不仅验证了本文方法的可行性，还进一步显示了它在处理简单波动问题时的精确性.

平面波传播问题是一种理想化的波动模型，波动能量沿单一方向传播的特性，使得每个人工边界上外行波的视传播速度都可以确定，这为设置人工边界带来了方便.实际波动问题则要复杂得多，其复杂性表现为由透射角度或波动能量空间扩散的衰减效应导致的在人工边界不同位置处，外行波的视传播速度各不相同.对于这种复杂的实际波动问题，本文MTF谱元格式同样能够很好地模拟.

3.2 内源问题

接下来进行内源问题的模拟.为显示本文方法适应复杂波场的能力，将波源设置在模型左下角以增加边界附近波场的复杂性.此时对于左侧和下侧边界，在靠近波源的位置会受到较为明显的波场几何衰减效应的影响，远离波源的位置则存在外行波透射角度较大的问题；上侧和右侧边界受到的影响相对较小.采用式(8)～式(11)给出的MTF谱元格式，插值多项式阶次取为M=4.由两组不同边界参数得到的0.3s,0.6s,0.9s,1.5s的波场快照如图7所示.

图7内源问题波场快照Fig.7 Wave fiel snapshotsof inner-source problem

图7 (a)使用的是MTF常用的边界参数：二阶MTF，ca=c，它对弱衰减、小角度透射的外行波模拟效果较好.图中结果显示在人工边界各个位置处，大部分波动能量都能够顺利透出，表明了本文MTF谱元格式模拟复杂波场问题的有效性.仅有少量反射出现在左侧、下侧边界距波源较近的波场衰减区和远端的大角度透射区.如：0.6s子图显示了由于波场衰减效应导致的未被完全透出的波；0.9s子图中A点和B点反射波强度对比，表明本例的边界参数设置对大角度外行波动的模拟效果不如小角度外行波动；1.5s子图中下边界的反射波强度大于上边界，因为下边界的透射角度更大.

为了保证水利工程的顺利招标，我国出台了各种法律法规，以此来确保招标工作开展的科学性和合理性。但是，在实际的招标过程中，仍然存在着不按照国家相应的法律法规进行招投标的现象。尤其是有些小型企业在招标的过程中进行相应的暗箱操作，这样就在一定程度上打破了整个行业的竞争平衡。而有些企业则是依靠人脉、金钱等进行非法的招标，这样不仅给公平、公正以及公开的招标环境造成了影响，而且还会影响招标工作的顺利展开。

图7(b)使用的边界参数为：二阶MTF，ca=c/cos45°，此时它对透射角度为45°左右的外行波模拟效果最好.由于本例中外行波在大部分人工边界区域都是以一定角度透射的，所以不难预测采用45°角为透射中心，能够显著改善模型4个角点附近，以及左侧、下侧边界上远离波源部分的透射效果.图中结果显示这组边界参数的模拟效果总体上要优于上组边界参数.如：0.9s子图中B点的法向反射波依稀可见，而A点附近接近45°角的反射波则几乎被完全吸收；1.5s子图中C点附近区域与上一组模拟结果相比，小角度反射波几乎看不出差别，但是大角度反射波被吸收得更加充分.

内源问题算例不仅再次证明了本文方法的可行性，而且表明了它对复杂波动问题的适应性，这种适应性主要通过选取适当的边界参数来实现.就工程应用而言，常用的MTF参数配置(N=2，ca=c)对很多实际波动问题已具有足够的模拟精度.然而从研究角度看，深入研究MTF计算参数与波动物理特征之间的联系，根据外行波的视传播规律有针对性地设定MTF计算参数，对更好地求解复杂波动问题具有重要价值.

4 边界反射误差的时域分析

式(8)～式(11)描述的MTF谱元格式对外行波的透射效果受到插值多项式阶次M、人工波速ca和透射阶次N三个参数的影响，通过一组精心设计的SH波动数值试验，在时域内定量地研究不同MTF参数取值下，本文MTF谱元格式在人工边界上的反射误差.

如图8所示，采用两个计算模型：半空间模型Ω2，左侧使用式(8)～式(11)的MTF谱元格式，其余三侧为自由边界；全空间模型Ω3，四边均为自由边界.反射误差的分析区域为Ω1，这样半空间模型的计算结果为受到左侧人工边界影响的数值解，全空间模型计算结果为不受任何边界影响的大区域数值解(即精确解).介质物理波速为1m/s，计算时间为2s，这组参数能够保证计算时间内自由边界的反射波不会返回分析区域.

图8 人工边界反射误差分析模型Fig.8 Numericalmodel for analyzing the reflectio errorof artificia boundary

为便于分析，本文模拟初始波场的自由扩散问题.由于没有外部波动能量输入以及介质阻尼等因素的干扰，模型中波动能量的变化只受边界影响.初始波场所在区域为图中圆形区域：圆心为(0.5,0)，半径为0.45m.分别采用三种具有不同传播特征的初始波场：圆形扩散(circularspread)、法向透射(normal transmission)和45°角透射(45°transmission)，其中圆形扩散初始波场的表达式为如下高斯型脉冲

其中，r2=(x-0.5)2+y2.法向透射的初始波场由在上式右端乘以cos(8πx)得到，45°角透射的初始波场则通过在上式右端乘以cos(8πx-8πy)得到.根据公式λmin=1/kmax来确定3种初始波场中最短波长成分，其中kmax为波数的上限值.由于初始波场是二维的，为简化分析，取波场变化最剧烈的方向进行一维傅里叶分析，圆形扩散情形计算exp（-30r2）的傅里叶谱，后两种情形计算exp（-30r2）·cos(8πr)的傅里叶谱.得到圆形扩散情形kmax的值约为4m-1，后两种情形约为8m-1.3种情形的初始速度都为零.

x和y方向的谱单元阶次均取为5，即采用36节点的单元.对于圆形扩散初始波场，单元尺寸取为1/4，对于法向透射和45°角透射初始波场，单元尺寸取为1/7.区域Ω1内的误差波场定义为由半空间模型Ω2计算的数值解减去由全空间模型Ω3计算的精确解.衡量各个时刻MTF边界反射误差的指标为：区域Ω1内误差波场的弗罗贝尼乌斯范数与初始波场的弗罗贝尼乌斯范数的比值.其中，各个时刻波动位移矩阵U的弗罗贝尼乌斯范数的计算公式为

4.1 插值多项式阶次的影响

3种初始波场情形下，MTF边界上的反射误差随插值多项式阶次M的变化规律如图9所示(MTF计算参数：N=2，ca=c，M=2,3,4,5).

图9 插值多项式阶次M对反射误差的影响Fig.9 Reflectio errorofMTFa ff ected by interpolation polynomial order M

图9插值多项式阶次M对反射误差的影响(续)Fig.9 Reflectio errorofMTFa ff ected by interpolation polynom ial order M(continued)

图9 结果显示，人工边界反射误差随着插值多项式阶次M增加而逐步降低，其中，当M=2时反射误差最大，3种情形的峰值误差均接近4%;当M=5时反射误差最小，3种情形下的峰值误差分别为3.5%,0.9%和2.7%.显然，当插值多项式阶次与谱单元阶次一致时，模拟效果最好.另外图中结果还有两点值得重视：一是所有反射误差均在5%以内，这说明如果从工程精度考虑，由插值多项式阶次不同引起的误差变化几乎可以忽略不计；二是不同插值多项式阶次对透射方向比较单一的波动问题，如法向透射或45°角透射，具有层次较为分明的影响，而对透射方向比较丰富的波动问题，如圆形扩散，其影响程度较小.

总之，插值多项式阶次对透射效果的影响是在较高精度水平上的微小变化，它并不是决定式(8)～式(11)的MTF谱元格式模拟效果的关键因素.从实际使用角度出发，考虑到边界附近的谱单元阶次会受到特定波动问题的模拟精度和计算效率的影响而有所不同，若存在多种不同插值多项式阶次的边界格式可供选择，会大大增加问题的复杂性，于是笔者建议统一采用插值多项式阶次与谱单元阶次一致的边界格式.

4.2 人工波速的影响

3种初始波场情形下，MTF边界上的反射误差随人工波速ca的变化规律如图10所示(MTF计算参数：N=2，M=5，α=ca/c，取值为0.5,0.75,1.0,1.25,1.5).

图10 人工波速ca对反射误差的影响Fig.10 Refection errorofMTFa ff ected by artificia wave speed ca

对比图9和图10结果可知，不同人工波速ca取值导致的边界反射误差变化幅度，大大超过不同插值多项式阶次M导致的变化.这点不难理解，因为人工波速是MTF定义式中的一个基本参数.较大反射误差出现在人工波速取值小于物理波速的情形，如ca=0.5c时，峰值误差分别达到10%,6.2%和13%，ca=0.75c时，峰值误差分别为5.3%,1.6%和6%.而当人工波速取值等于或大于物理波速时，反射误差均处于较低水平，其峰值分别不超过3.5%,1.7%和2.7%.反射误差越小，表明人工波速取值越接近外行波沿边界法向的视传播速度.

波动问题的复杂性在于受到透射角度、波幅衰减或多个子波叠加等因素的影响，确定外行波的视传播速度常常比较困难，而且很多时候难以用一个值来精确地描述.但这部分研究能够帮助我们认识不同因素对视传播速度的具体影响，进而总结一些有助于优化人工波速取值的结论.图10(c)结果表明：外行波透射角度对视传播速度的影响是较为确定的，以θ角透射的外行波的视传播速度为c/cosθ.图中模拟的是45°角透射波，其视传播速度为α=1.5曲线与此最接近，因此反射误差最小，其他α取值偏离越远，反射误差越大.图10(b)模拟的是法向透射波，但反射误差最小值没有出现在α=1时，这表明了初始波峰附近波动能量几何扩散导致的波幅衰减效应带来的影响.从图中结果不难看出，采用稍大于介质物理波速的人工波速能够有效地弱化这一影响，而本例中最优人工波速取值在c和1.25c之间.图10(a)结果显示了单一人工波速的MTF谱元格式对实际复杂波动问题的模拟效果，它与不同时刻波动能量在边界上的透射区域和透射方向有关.在1.0 s之前，法向和小角度透射波先到达边界，此时人工波速等于c的模拟效果最好；在1.0s之后，大角度透射波成为主要成分，此时人工波速等于1.25c或1.5c的模拟效果要好于前者；而对于人工波速等于0.5c或0.75c的格式，其模拟效果自始至终都不理想.

以上分析表明，实际波动问题与理想的平面波法向透射问题相比，存在透射角度或波幅衰减，二者都会导致视传播速度增大，因此采用大于介质物理波速的人工波速能够达到更好的模拟效果.这一现象可理解为：从人工边界节点的视角来看，存在透射角度或波幅衰减加速了波动能量向人工边界的“积聚”过程.由此可以推断对于本例未涉及的多个子波叠加效应，其作用与前两者相似.

至此，可将均匀介质中标量波问题的人工波速取值规律总结为：不应取小于介质物理波速的值；对于大多数实际波动问题，取等于或稍大于介质物理波速的值就能满足要求；只有在少数特殊情形下，如透射角度较大、边界距离点波源很近或较多子波场在边界附近叠加时，才可能需要采用大大超过介质物理波速的值，但这往往会带来诸如稳定性等方面的问题.

4.3 透射阶次的影响

3种初始波场情形下，MTF边界上的反射误差随透射阶次N的变化规律如图11所示(MTF计算参数：M=5，ca=c，N=1,2,3).

图11 透射阶次N对反射误差的影响Fig.11 Reflectio errora ff ected by transm itting order N

将图11和图10、图9对比之后可以看出，透射阶次N对边界反射误差的影响是决定性的，它大大超过了人工波速ca或插值多项式阶次M的影响.其中，一阶MTF的反射误差较大，如圆形扩散和45°角透射情形，峰值误差分别为8.7%和10.2%，不满足工程精度要求，说明一阶MTF难以用来模拟复杂波动问题.二阶MTF的反射误差显著降低，三种情形下的峰值误差分别为3.5%,0.9%和2.7%，都小于5%，满足工程精度要求，说明通常采用二阶MTF进行波动模拟的做法是可取的.三阶MTF的反射误差没有像预想的那样进一步降低，反而出现了不合理的增长，如圆形扩散和法向透射情形，三阶MTF的反射误差大大超过一阶MTF，45°角透射情形，三阶MTF的反射误差与二阶MTF相当，但前者随时间增长的趋势一直在延续.我们还采用四阶MTF进行了模拟，其反射误差已经大到无法在图中显示的地步，更高阶次的MTF同样如此.这种现象表明N≥3时MTF边界发生了某种失稳.

观察三阶MTF模拟的各个时刻波场变化，我们发现这是一种飘移失稳[38]，如图12所示.飘移现象从MTF边界中部开始，逐渐向两端扩展，形成中间大、两端小的拱形位移，拱的高度随时间不断上升，而实际边界位移在主波形透过之后应当趋近于零.显然，当边界出现这种严重的失稳问题时，波动模拟过程失败.从实用角度看，不使用N≥3的MTF谱元格式就可以回避这个问题，因为这里二阶MTF没有出现失稳，而且通常具有足够的模拟精度.关于本文MTF谱元格式的飘移失稳机理及其抑制措施，我们进行了初步研究，已得到一些基本认识.

图12 区域Ω1在2s时的波场Fig.12 Wave fiel of domainΩ1at2s

飘移失稳现象不是本文MTF谱元格式所特有的，传统MTF有限元格式可能出现飘移失稳已是众所周知[3841].而其他类型的人工边界，如文献[42]提出的几种Higdon边界形式，同样可能出现飘移失稳.关于飘移失稳机理目前还没有统一的解释，如：李小军等[3839]认为飘移失稳是由于外源问题采用连续介质中的波动理论解作为输入波场，而它与有限元离散解之间存在差异所造成的.周正华等[40-41]认为是由于MTF不满足GKS准则导致数值解中的零频和零波数成分可以通过边界进入计算区，从而引起飘移失稳.Higdon[42]则将其解释为边界格式允许波动模拟中出现关于零频成分的广义特征值.我们认为这可能是由于人工边界条件作为计算模型的支承条件本身就存在一定缺陷所致.

数值试验结果表明，本文MTF谱元格式可能出现飘移失稳的透射阶次，对于内源问题为N≥3，对于外源问题为N≥2.其中，前者与文献[42]的结果一致，后者则与文献[38]结论相同.由此不难看出，飘移失稳可能是边界高阶形式固有问题，而当边界上存在外部波动输入时，这一问题更为严重.已有研究提出了几种抑制MTF飘移失稳的方法，如降低MTF阶次[38]、在边界区添加弹簧和阻尼元件[39]，以及在MTF公式中增加一个微小正数[40].我们认为这些方法都可以应用到本文的MTF谱元格式当中，这部分内容值得深入探讨，但已超出本文的范畴.

最后，对本文MTF谱元格式的高频振荡失稳特性做一简要说明.研究过程中我们进行了大量SH波动模拟，边界均未出现高频振荡失稳.结合文献[1]关于一维波动问题的研究，我们将其解释为：本文MTF谱元格式的稳定条件比内域离散格式宽松，所以一般情况下不会发生高频振荡失稳，只有在人工波速大大超过介质物理波速或透射阶次很高时，才有可能超出稳定界限，在边界上出现高频振荡失稳.

5 结论

(1)一维化的MTF公式需要通过一条直线上的离散网格点来实现，本文通过假定紧邻人工边界的一层谱单元为直线边界的四边形单元，保证每个人工边界节点都唯一地对应一条指向内域的离散网格线，用于实现文献[1]所提出的MTF谱元格式.外源和内源SH波动问题算例证实了方法的可行性以及对复杂波动问题的适应性.

(2)插值多项式阶次M对边界反射误差的影响较小，M取值较高的MTF谱元格式精度略好于M取值较低的格式.从使用方便的角度，我们建议M取值统一采用MTF计算方向的谱单元阶次，这同时也是精度最好的选择.

(3)人工波速ca取值对边界反射误差具有显著影响.就本文研究的均匀介质中SH波动问题而言，外行波存在透射角度、波幅衰减或多个子波场叠加，都会导致波动能量更快地向边界“积聚”，因此采用等于或稍大于介质物理波速的人工波速能得到较好的模拟效果.采用小于或等于介质物理波速的人工波速，推算边界节点位移的准确性较差.其他类型波动问题，如涉及到成层场地、P-SV波动，或Rayleigh面波等，情况较为复杂，可在此基础上进行专门讨论.

(4)透射阶次N对本文MTF谱元格式的模拟效果具有决定性影响.不发生失稳时，提高透射阶次对降低反射误差的效果最为明显.但是，高阶MTF容易出现飘移失稳，其中，内源问题和外源问题开始出现飘移失稳的MTF阶次分别为3和2.

1邢浩洁,李鸿晶.透射边界条件在波动谱元模拟中的实现：一维波动.力学学报,2017,49(2)：367-379(Xing Haojie,LiHongjing.Implementation ofmulti-transm itting boundary condition for thewave motion simulation by spectralelementmethod：onedimension case.Chinese JournalofTheoreticaland Applied Mechanics,2017,49(2)：367-379(in Chinese))

2 Patera AT.A spectral elementmethod for flui dynam ics：Lam inar fl w in a channelexpansion.JournalofComputational Physics,1984,54(3)：468-488

3 Priolo E,SerianiG.A numerical investigation of Chebyshev spectral elementmethod for acoustic wave propagation//Proceedings of the13th IMACSConferenceon Comparative Applied Mathematics,Dublin,Ireland,1991,2：551-556

4 SerianiG,Priolo E.Spectralelementmethod foracousticwavesimulation in heterogeneousmedia.Finite Elements in Analysis and Design,1994,16(3)：337-348

5 Seriani G.A parallel spectral elementmethod for acoustic wave modeling.JournalofComputationalAcoustics,1997,5(1)：53-69

6 SerianiG,Su Chang.Wave propagationmodeling in highly heterogeneousmedia by a ploy-grid Chebyshev spectralelementmethod.JournalofComputationalphysics,2012,20(2)：345-351

7 Komatitsch D,Vilotte JP.The spectralelementmethod：an e ffi cient tool to simulate theseism ic responseof 2D and 3D geologicalstructures.Bulletin ofthe SeismologicalSociety ofAmerica,1998,88(2)：368-392

8 Komatitsch D,Vilotte JP,VaiR,etal.The spectralelementmethod for elastic wave equations—application to 2-D and 3-D seism ic problems.International Journal for Numerical Methods in Engineering,1999,45(9)：1139-1164

9 Komatitsch D,Barnes C,Tromp J.Wave propagation near a flui solid interface：a spectral-element approach.Geophysics,2000,65(2)：623-631

10 Komatitsch D,Liu Qinya,Tromp J,et al.Simulations of ground motion in the Los Angeles basin based upon the spectral-element method.Bulletin of the Seismological Society of America,2004,94(1)：187-206

11王秀明,SerianiG,林伟军.利用谱元法计算弹性波场的若干理论问题.中国科学：G辑,2007,37(1)：1-19(Wang Xium ing,SerianiG,Lin Weijun.Several theoretic aspects for the calculation of elasticwave fiel using spectralelementmethod.Science China(G series),2007,37(1)：1-19(in Chinese))

12严珍珍,张怀,杨长春等.汶川大地震地震波传播的谱元法数值模拟研究.中国科学：D辑,2009,39(4)：393-402(Yan Zhenzhen,Zhang Huai,Yang Changchun,et al.Numerical investigation of seismic wave propagation of Wenchuan Earthquake using spectral elementmethod.Science China(D series),2009,39(4)：393-402(in Chinese))

13李孝波.基于谱元法的玉田震害异常研究.[博士论文].哈尔滨：中国地震局工程力学研究所,2014(Li Xiaobo.The investigation of seism ic damage anomaly in Yutian based on spectral element method.[PhD Thesis].Harbin：Institute of Engineering Mechanics,Chinese Earthquake Adm inistration,2014(in Chinese))

14李孝波,薄景山,齐文浩等.地震动模拟中的谱元法.地球物理学进展,2014,29(5)：2029-2039(LiXiaobo,Bo Jingshan,QiWenhao,etal.Spectralelementmethod in seism ic groundmotion simulation.Progress in Geophysics,2014,29(5)：2029-2039(in Chinese))

15刘有山,滕吉文,徐涛等.三角网格谱元法地震波场数值模拟.地球物理学进展,2014,29(4)：1715-1726(Liu Youshan,Teng Jiwen,Xu Tao,etal.Numericalmodeling of seism ic wave fiel w ith the SEM based on triangles.Progress in Geophysics,2014,29(4)：1715-1726(in Chinese))

16李琳,刘韬,胡天跃.三角谱元法及其在地震正演模拟中的应用.地球物理学报,2014,57(4)：1224-1234(Li Lin,Liu Tao,Hu Tianyue.Spectral elementmethod w ith triangularmesh and its application in seismicmodeling.Chinese Journal ofGeophysics,2014,57(4)：1224-1234(in Chinese))

17李洪建,韩立国,巩向博.复杂构造网格化及高精度地震波场谱元法数值模拟.石油物探,2014,53(4)：375-383(LiHongjian,Han Liguo,Gong Xiangbo.High precision spectral elementmethod based on grid discretization of complicated structure for seism ic wavefiel numerical simulation.Geophysical Prospecting for Petroleum,2014,53(4)：375-383(in Chinese))

18曹丹平,周建科,印兴耀.三角网格有限元法波动模拟的数值频散及稳定性研究.地球物理学报,2015(5)：1717-1730(Cao Danping,Zhou Jianke,Yin Xingyao.The study for numerical dispersion and stability ofwavemotion w ith triangle-based finit element algorithm.Chinese JournalofGeophysics,2015(5)：1717-1730(in Chinese))

19 He Chunhui,Wang Jinting,Zhang Chuhan.Nonlinear spectralelementmethod for 3D seismic-wave propagation.Bulletin of the Seismological Society ofAmerica,2016,106(3)：1074-1087

20林灯,崔涛,冷伟等.一种求解地震波方程的高效并行谱元格式.计算机研究与发展,2016,53(5)：1147-1155(Lin Deng,Cui Tao,Leng Wei,et al.An e ffi cient parallel spectral element scheme for solving seism ic wave equation.JournalofComputer Research and Development,2016,53(5)：1147-1155(in Chinese))

21 Liu Youshan,Teng Jiwen,Lan Haiqiang,etal.A comparativestudy of finit elementand spectralelementmethods in seism ic wavefiel modeling.Geophysics,2014,79(2)：T91-T104

22 Chen Yushu,Yang Guangwen,Ma Xiao,et al.A time-space domain stereo finit di ff erencemethod for3D scalarwavepropagation.Computers&Geosciences,2016,96：218-235

23 Liu Shaolin,LiXiaofan,WangWenshuai,etal.Am ixed-grid finit elementmethod w ith PML absorbing boundary conditions for seism ic wavemodelling.Journal ofGeophysics&Engineering,2014,11(5)：1-13

24 Clayton R,EngquistB.Absorbing boundary conditions for acoustic and elastic waveequations.Bulletin ofthe Seismological Society of America,1977,67(6)：1529-1540

25 Keller JB,Dan G.Exactnon-reflectin boundary conditions.JournalofComputational Physics,1989,82(1)：172-192

26 Chaljub E,Komatitsch D,Vilotte JP,etal.Spectral-elementanalysis in seismology.Advances in Geophysics,2007,48：365-419

27 He Ying,M in M isun,Nicholls DP.A spectralelementmethod w ith transparentboundary condition for periodic layered media scattering.JournalofScientifi Computing,2016,68(2)：772-802

28赵密.近场波动有限元模拟的应力型时域人工边界条件及其应用.[博士论文].北京：北京工业大学,2009(Zhao M i.Stress-type time-domain artificia boundary condition for finite-elemen simulation ofnear-fiel wavemotion and itsengineering application.[PhD Thesis].Beijing：Beijing University of Technology,2009(in Chinese))

29 Berenger JP.A perfectlymatched layer for theabsorption ofelectromagnetic waves.Journal ofComputational Physics,1994,114(2)：185-200

30 Komatitsch D,Tromp J.A perfectlymatched layerabsorbingboundary condition for thesecond-order seism icwaveequation.Geophysical Journal International,2003,154(1)：146-153

31 Martin R,Komatitsch D,Gedney SD.A variational formulation of a stabilized unsplit convolutional perfectly matched layer for the isotropic or anisotropic seism ic wave equation.Computer Modeling in Engineering&Sciences,2008,37(3)：274-304

32廉西猛,单联瑜,隋志强.地震正演数值模拟完全匹配层吸收边界条件研究综述.地球物理学进展,2015,30(4)：1725-1733(Kang Ximeng,Shan Lianyu,Sui Zhiqiang.An overview of research on perfectly matched layers absorbing boundary condition of seismic forward numericalsimulation.Progress in Geophysics,2015,30(4)：1725-1733(in Chinese))

33 Xie Zhinan,Matzen R,Cristini P,etal.A perfectlymatched layer for fluid-soli problems：Application to ocean-acousticssimulations w ith solid oceanbottoms.Journalofthe AcousticalSociety ofAmerica,2016,140(1)：165-175

34廖振鹏,黄孔亮,杨柏坡等.暂态波透射边界.中国科学：A辑,1984,(6)：556-564(Liao Zhenpeng,Huang Kongliang,Yang Baipo,et al.A transm itting boundary for transientwave.Scientia Sinica(SeriesA),1984,(6)：556-564(in Chinese))

35 Liao Zhenpeng,Wong Hongliang.A transm itting boundary for the numerical simulation of elastic wave propagation.International JournalofSoilDynamicsand Earthquake Engineering,1984,3(4)：174-183

36戴志军,李小军,侯春林.谱元法与透射边界的配合使用及其稳定性研究.工程力学,2015,32(11)：40-50(Dai Zhijun,Li Xiaojun,Hou Chunlin.A combination usageof transm itting formulaand spectral elementmethod and the study of its stability.Engineering Mechanics,2015,32(11)：40-50(in Chinese))

37陈少林,赵宇昕.一种三维饱和土--基础--结构动力相互作用分析方法.力学学报,2016,48(6)：1362-1371(Chen Shaolin,Zhao Yuxin.A method for three-dimensional saturated soil-foundationstructure dynam ic interaction analysis.Chinese JournalofTheoreticaland Applied Mechanics,2016,48(6)：1362-1371(in Chinese))

38李小军,廖振鹏.时域局部透射边界的计算飘移失稳.力学学报,1996,28(5)：627-632(LiXiaojun,Liao Zhenpeng.Thedriftinstability of local transmitting boundary in time domain.Acta Mechanica Scinica,1996,28(5)：627-632(in Chinese))

39李小军,杨宇.透射边界稳定性控制措施探讨.岩土工程学报,2012,34(4)：641-645(LiXiaojun,Yang Yu.Measures for stability controlof transmitting boundary.Chinese Journal ofGeotechnical Engineering,2012,34(4)：641-645(in Chinese))

40周正华,廖振鹏.消除多次透射公式飘移失稳的措施.力学学报,2001,33(4)：550-554(Zhou Zhenghua,Liao Zhenpeng.A measure for eliminating drift instability of the multi-transmitting formula.Acta Mechanica Scinica,2001,33(4)：550-554(in Chinese))

41廖振鹏,周正华,张艳红.波动数值模拟中透射边界的稳定实现.地球物理学报,2002,45(4)：533-545(Liao Zhenpeng,Zhou Zhenghua,Zhang Yanhong.Stable implementation of transm itting boundary in numericalsimulation ofwavemotion.Chinese Journal ofGeophysics,2002,45(4)：533-545(in Chinese))

42 Higdon RL.Absorbing boundary conditions for di ff erence approximations to themulti-dimensionalwave equation.Mathematics of Computation,1986,47(176)：437-459

43景立平.多次透射公式实用形式稳定性分析.地震工程与工程振动,2004,24(4)：20-24(Jing Liping.Stability analysis of practical formula formulti-transm itting boundary.Earthquake Engineering and Engineering Vibration,2004,24(4)：20-24(in Chinese))

44章旭斌,廖振鹏,谢志南.透射边界高频耦合失稳机理及稳定实现——SH波动.地球物理学报,2015,58(10)：3639-3648(Zhang Xubin,Liao Zhenpeng,Xie Zhinan.Mechanism ofhigh frequency coupling instability and stable implementation for transmitting boundary——SH wavemotion.Chinese Journal ofGeophysics,2015,58(10)：3639-3648(in Chinese))

IMPLEMENTATIONOFMULTI-TRANSM ITTING BOUNDARY CONDITION FOR WAVEMOTION SIMULATION BY SPECTRAL ELEMENTMETHOD:TWO DIMENSION CASE1)

Xing Haojie LiHongjing2)
(College ofCivilEngineering，Nanjing Tech University，Nanjing 211816，China)

Thespectralelementschemeofmulti-transm itting formula(MTF)which proposed by Xing and Liisdeveloped and expanded to numerical simulation of SH-typewavemotion in homogeneous linearmedium.The spectral elements adjacent to the artificia boundariesare supposed to be rectilinear quadrilateralelements,in order thateach node on the artificia boundary locates on a unique grid line pointing to the interior domain,hence the displacementof a particular artificia boundary node at a certain time can be inferred from displacements of nodes on the grid line atearlier times.Two numericalexamples,the planewave propagationw ith a certain angle(outer-source problem)and the free spread ofa pulsewaveoriginated from a point(inner-source problem),areemployed for verificatio of thiswavemotion simulationprocedurewhich combines the spectral elementmethod w ith MTF boundary.Themain parameters a ff ecting reflectio error of the MTF scheme,such as order of interpolation polynom ial,artificia wave speed and transm itting order,are investigated in time domain via a seriesof initial-value problems.The results show that the order of interpolation polynomial influence the reflectio error very little,while higher interpolation ordermay lead to better accuracy.For the choice of artificia wave speed,it is preferable to choose values equal or slightly greater than the physicalwave speed,otherw ise bigger reflectio errors come about.Transm itting order exerts themost significan impacton reflectio error,which would be reduced greatly with the increase of transm itting order of MTF,but the undesired drift instabilitymay arisewhen transmitting order reaches three or two for inner-source and outer-source problems,respectively.Although thework combining MTF boundary w ith spectralelementmethod is conducted on simulation of SH wave propagation in infinit homogeneousmedium in thispaper,ithas laid the foundation of research onmore complicated situations,such as issuesof layered soilsites,propagation of P-SV waveorRayleighwave,etc.

multi-transm itting formula,spectral elementmethod,SH-type wavemotion,spectral element scheme of MTF,drift instability

P315

A

10.6052/0459-1879-16-393

2016-12-25收稿，2017-03-01录用,2017-03-06网络版发表.

1)国家自然科学基金资助项目(51278245).

2)李鸿晶，教授，主要研究方向：地震工程学.E-mail：hjing@njtech.edu.cn

邢浩洁,李鸿晶.透射边界条件在波动谱元模拟中的实现：二维波动.力学学报,2017,49(4)：894-906

Xing Haojie,LiHongjing.Implementation ofmulti-transm itting boundary condition forwavemotion simulation by spectral element method：two dimension case.Chinese JournalofTheoreticaland Applied Mechanics,2017,49(4)：894-906