波动谱元模拟中透射边界稳定性分析

章旭斌，谢志南

(中国地震局工程力学研究所，中国地震局地震工程与工程振动重点实验室，哈尔滨 150080)

对于大多数只关心广义结构及其邻近介质中波动的问题，其本质是无限域波动问题。对无限域波动数值模拟的关键在于将其转换成有限域的数值计算，因此需引入人工边界截取有限计算区，并设置人工边界条件模拟外部无限域对外行波的能量辐射。无限域波动数值模拟涉及内域数值算法和人工边界条件，迄今常用的内域算法有差分法、有限元法和谱元法[1]等，而人工边界有粘弹性边界[2-3]、连分式边界[4-6]、透射边界[7-8]和完美匹配层[9-10]等。由内域和边界算法两者组合，将形成多种多样的波动数值模拟方案。

谱元法(Spectral element method, SEM)是一种节点非均匀分布的高阶有限元，其必然存在虚假模态。虚假模态也即高阶离散格式中存在的虚假解，描述了离散网格中允许的虚假运动形式，可由频散分析得到。MULDER[11]分析了谱元法中的虚假模态对数值模拟精度的影响并不明显。DE BASABE 等[12]分析了谱元频散和稳定条件，阐述了4 阶以上的谱元一波长只要4 个～5 个节点，就基本消除了数值频散和各向异性。KOMATITSCH和TROMP[13]将谱元法推广到三维波动数值模拟，并编写了SPECFEM3D 程序，推动了谱元法在地震波动数值模拟中的应用[14-16]。由廖振鹏等[7-8]提出的透射边界亦称多次透射公式(Multi-transmitting formula, MTF)，是一类典型的具有高阶精度、低计算量、且易于实现的人工边界条件，其已在地震工程、电磁学、流体力学等数值模拟领域得到应用[17-19]。已有研究将MTF 应用于波动谱元模拟中，戴志军等[20]在谱元法中应用了时间插值和空间插值形式的MTF，显示后者具有很好的稳定性。于彦彦等[21]对比了粘性边界和二阶MTF，MTF 明显提高了数值模拟精度。XING 等[22]基于数值实验讨论了MTF 的精度和稳定性，针对不同入射角设置不同人工波速的MTF，其吸收精度可接近完美匹配层的精度。上述研究显示了MTF 具有较好的模拟精度和稳定性，然而，在波动谱元模拟中MTF 仍可能引发数值失稳现象，其失稳机理和稳定条件尚不明确。下面梳理MTF 稳定性研究结果及其稳定性分析方法。

依据已有的研究结果，由MTF 引发的失稳形式表现为零频飘移失稳、反射放大失稳和高频振荡失稳，后两者为高频失稳。LIAO 和LIU[23]针对一维波动集中质量有限元模拟，基于边界反射系数揭示了MTF 引发的反射放大失稳机理，即在有限区域内人工边界对高频行波的反复反射放大所致。景立平等[24]针对二维SH 波动集中质量有限元模拟，基于GKS 定理的群速度解释[25]揭示了高频振荡失稳的机理，即边界和内域格式支持群速度指向内域的高频谐波。依据对高频失稳的认识，已提出诸多稳定措施[26-28]。零频飘移失稳是多种数值格式(如粘性边界和Higdon 边界等)常见的失稳现象，通常采用小量修正人工边界格式[29]或降阶方法[30]抑制飘移失稳。

以上研究工作基于反射系数分析和GKS 定理的群速度解释，揭示了波动有限元模拟中MTF 引发的反射放大失稳和高频振荡失稳机理。而关于稳定性分析方法，应提及BEAM 等[31]提出的P 稳定(Practical stability)，即在保证GKS 稳定条件下，仍需保证数值解不允许随时间增长。TREFETHEN[25]论证了P 失稳(P-instability)往往表现为反射放大失稳，并进一步阐明了GKS 定理群速度解释所针对的失稳模态的反射系数为无穷大，即对失稳模态而言，其仅有反射波而无入射波。因此，MTF引发的2 种高频失稳可统一采用反射系数来分析。

基于反射系数分析得到的线性有限元中MTF稳定条件[23,32]并不适用于谱元法。谱元法是高阶有限元，其单元内多个非等间距节点的组合呈现周期延拓，故而存在多条包括真实模态和虚假模态在内的频散曲线，而线性有限元是单一节点的周期延拓，仅存在一条逼近波动方程的频散曲线。因此，在谱元法中分析人工边界反射系数时需考虑这一特点。本文将依据谱元节点周期延拓性质，通过构建谱元和MTF 数值格式的向量形式，进而基于向量形式推导MTF 反射系数，从而分析人工边界稳定性，揭示边界引发失稳的机理，并给出稳定条件，为实际波动模拟参数选取提供参考。同时其研究思路对于其他类型波动和人工边界稳定性分析具有参考价值。

1 数值计算格式

进而结合时间中心差分法，可得时空解耦的数值格式：

2 透射边界反射系数

在分析离散格式中人工边界的反射系数时，需要知道内域格式的频散关系，也就是内域格式在无限域上支持的谐波频率和波数之间的关系，可通过频散分析得到。而频散分析是将内域节点看作在无限域上的周期延拓。由于谱元是一种高阶有限元，其以单元内多个节点的组合呈现周期延拓，如四阶谱单元具有5 个节点，以1 个边界节点和3 个内节点进行周期延拓(如图1 所示)，其有别于线性有限元以单一节点的周期延拓。因此，分析谱元频散需考虑节点周期延拓这一性质。

图1 四阶谱单元节点周期延拓示意图Fig. 1 Periodic extension of four nodes in fourth-order SEM

同样的，在分析边界反射系数时也应考虑节点周期延拓性质，故而构建数值格式的向量形式，即将内域谱元周期延拓的节点组装成向量形式的内域格式，将边界节点与邻近谱单元内的节点组装成向量形式的边界格式，从而构建了与原格式等价的向量形式格式，同时这一向量形式规避谱元非等间距节点分布给分析造成的困难。实际波动数值模拟通常采用四阶谱元和二阶MTF，图2 为边界节点附近谱单元节点分布。依据MTF格式(8)和内域计算格式(3)，组装后的边界格式为：

图2 人工边界节点和四阶谱单元节点组合示意图Fig. 2 Combination of nodes in artificial boundary and fourth-order SEM

人工边界对入射波的吸收精度主要体现在反射系数的模，下文所指反射系数均表示反射系数的模。值得说明的是，式(10)中的入射波和反射波由其群速度方向确定，群速度向外的为入射波，群速度向内的为反射波[25]。因此，若谐波的相速度向外而群速度向内(由频散关系确定)，依据式(12)计算的反射系数取其倒数。

依据谱单元节点周期延拓性质(图1)和内域节点格式(3)，组装后的内域格式为：

图3(a)为 Δτ=0.1时四阶谱元的频散，其存在4 条频散曲线，仅其中1 条为具有物理意义的真实模态频散(实线)，即对应式(14)的最小特征值[12]。其余3 条均为虚假模态的频散，这是高阶格式所固有的特性，频散曲线的数量由周期延拓的节点个数所决定。图3(b)为MTF 的反射系数，包括边界对真实模态和虚假模态的反射系数。从图中可以看出，虚假模态的反射系数较大。虽然虚假模态对谱元模拟精度影响很小[11]，然而其对人工边界稳定性的影响却不可忽略。图3(c)为不同边界参数取值下，二阶MTF 对真实模态的反射系数，可以看出当S=Δτ时 (即ca=c，人工波速取为物理波速)，MTF 的吸收精度最好。

图3 四阶谱元的频散曲线和二阶MTF 的反射系数Fig. 3 The dispersion curves of fourth-order SEM and the reflection coefficients of second-order MTF

下面采用数值实验验证本文得到的反射系数。计算模型长1000 m，密度为2.2 g/cm3，波速为500 m/s。距离左端250 m 处施加Ricker 子波荷载，主频5 Hz。左端设置MTF，右端为自由边界，观测点设置在距离左边界0 m、50 m 和100 m 的位置。采用四阶谱元离散，Δx=10 m ，Δt=0.002 s，此时 Δτ=0.1。MTF 采用二阶格式，边界参数S=0.05。

图4(a)为观测点位移时程，100 m 处观测点入射波和反射波完全分离且间隔时间较短，可由其入射波和反射波的傅里叶谱比得到MTF 反射系数。从图4(b)可以看出，数值实验和理论计算的真实模态的反射系数在0 Hz～10 Hz 范围内基本完全一致。

图4 三个观测点的位移时程以及基于理论计算和数值实验的二阶MTF 反射系数Fig. 4 The displacement time history of three receivers and the reflection coefficient of second-order MTF computed by theoretical calculation and numerical experiment

3 透射边界稳定条件

波动有限元模拟中MTF 引发反射放大失稳和高频振荡失稳。前者为在有限的计算区内边界对外行波的反复发射放大所致，后者为边界和内域共同支持了群速度指向内域的平面谐波。实际上这两类失稳可统一采用人工边界反射系数来分析，两者仅在反射系数上表现差异，前者为反射系数大于1，后者为反射系数无穷大，即仅有反射波，而无入射波[25]。因此，MTF 引发的两种高频失稳可统一采用反射系数来分析。

另外，谱元法中存在虚假模态，虽然对数值模拟的精度影响很小，但对数值稳定性不可忽略。因此，在分析边界稳定性时，必须同时避免虚假模态和真实模态的反射放大。由反射系数式(12)，可以得到MTF 在谱元法中的稳定条件为：

由于四阶谱元并未有解析形式的频散关系，故而通过数值方法给出谱元法中MTF 的稳定条件。图5(a)为一阶MTF 稳定条件S≤0.195，与Δτ无关，其中 Δτ取值限定在四阶谱元的稳定条件上限0.147。图5(b)为二阶MTF 稳定条件，其表现为无量纲边界参数S=caΔt/Δx和谱元参数Δτ=cΔt/Δx之 间的关系，当 Δτ取 值变大时，S取值上限变小。一维线性有限元中MTF 的稳定条件[23]为S≤1.5，MTF 在高阶谱元中的稳定性较线性有限元更为复杂。从图5 中S=Δτ曲线可以发现，对于一阶MTF 而言，这一取值方式在整个区间内都是稳定的，而对于二阶MTF 而言，这一取值方式存在失稳的情况，即当 Δτ取值接近四阶谱元稳定条件上限0.147 时。

图5 四阶谱元法中MTF 的稳定条件Fig. 5 The stability condition of MTF in SEM

下面通过第2 节的数值算例加以验证，取S=0.4，其他参数不变。图6(a)中3 个观测点的位移时程在较短时间内出现数值失稳，其中100 m处观测点的失稳时程与自由面反射波叠加，因此显得不对称。图6(b)为模拟参数取值下二阶MTF 反射系数，有两条反射系数曲线均存在大于1 的值。

图6 三个观测点的位移时程和模拟参数条件下二阶MTF 的反射系数Fig. 6 The displacement time history of three receivers and the reflection coefficient of second-order MTF under the given simulation parameters

实际波动数值模拟中MTF 人工波速通常取为物理波速，此时S=Δτ。图5(b)显示二阶格式S=Δτ=0.147时，处于失稳区域。同样采用第2 节算例加以验证。图7(a)中的观测点位移模拟结果显示失稳。图7(b)为截取的20 s 数值结果显示了失稳形式，其失稳峰值在逐步放大。图7(c)为模拟参数取值下二阶MTF 的反射系数，其中一条虚假模态的反射系数大于1。因此，数值失稳机理为MTF 对虚假模态的反射放大及在有限区域内的反复反射放大所致。该数值实验中失稳出现的时间很晚，其原因在于上述取值情况下虚假模态对应的反射系数没有远大于1，并且数值模拟中的虚假模态成分很小，以致其失稳过程缓慢。保持其他参数取值不变，取S=0.137，处于稳定范围。图8(a)数值模拟结果未出现失稳，图8(b)为计算了500 万步的最后20 s 位移时程，未出现微小扰动，没有失稳迹象。图8(c)为模拟参数取值下二阶MTF 的反射系数，其值均小于等于1。

图7 失稳的观测点位移时程和模拟参数条件下二阶MTF 的反射系数Fig. 7 The unstable displacement time history of receiver and the reflection coefficient of second-order MTF under the given simulation parameters

图8 稳定的观测点位移时程和模拟参数条件下二阶MTF 的反射系数Fig. 8 The stable displacement time history of receiver and the reflection coefficient of second-order MTF under the given simulation parameters

对于实际波动模拟而言，本文给出的稳定条件过于严格。因为实际波动模拟通常在较短的时间内，而对于不满足稳定条件的S取值，若其反射系数并未远大于1，失稳误差在较短的时间内并不显著。因此，对于实际波动模拟，S取值可适当放宽。

本文给出了一维波动谱元模拟中MTF 的稳定条件，可作为多维波动稳定条件的参考。下面以二维SH 波动模拟为例，阐述MTF 的稳定性。对于二维平面波入射问题(如图9)，为了吸收大角度入射波，通常需要增大MTF 的人工波速取值。若平面波入射角为θ，SH 波速为cs，则右侧MTF 的人工波速可取为：

图9 平面波入射示意图Fig. 9 Diagram of plane wave incident

由式(16)可知，当入射角很大时，ca会取得很大，以致S=caΔt/Δx取值过大可能引发数值失稳。

图10 为盆地-基岩半空间模型，计算模型长度4 km、深度2 km，盆地长度2 km、深度0.3 km，基岩介质密度 ρ=2.5 g/cm3，cs=2.5 km/s，盆地介质密度 ρ=1.5 g/cm3，cs=0.5 km/s。采用四阶谱元离散，单元平均尺寸 Δx=Δy=50 m，时间步长Δt=0.0014 s ， 此时 Δτ=0.07。除自由地表外，其余三边均采用二阶MTF，左边和底边MTF 的人工波速取为基岩SH 波速，右边MTF 的人工波速按式(16)取值。3 个观测点(图10 中菱形)设置在右边边界上。入射平面波时间函数采用Ricker 子波，主频10 Hz。

图10 盆地-基岩半空间模型Fig. 10 Basin-rock half space model

图11(a)为 θ=80°时计算的观测点位移时程，MTF 引发数值失稳，其中MTF 人工波速依据式(16)计算为ca≈14.4 km/s ，此时S≈0.403。参考一维稳定条件S≤0.243 ， 可得ca＜8.6 km/s。图11(b)为 θ=80°时，取ca=8.0 km/s计算的观测点位移时程未出现数值失稳。图12 为各个时刻的波场图，从1.2 s 波场图可以看出右侧边界MTF 基本没有反射误差。对于大角度平面波入射，MTF 人工波速取值可参考一维稳定条件，可适当调大以提高对大角度入射波的吸收效果，但取值不宜过大，以免引发数值失稳。

图11 右边界上三个观测点的位移时程Fig. 11 The displacement time history of three receivers located at the right boundary

图12 谱元和MTF 模拟的各个时刻波场Fig. 12 Propagation of wave simulated by SEM with MTF

4 结论

本文分析了谱元法中MTF 的反射系数，揭示了MTF 引发失稳的机理，同时给出了MTF 稳定条件，并通过数值实验加以验证。主要结论有：

(1) 谱元法中虚假模态虽然对数值模拟的精度影响很小，但会引发人工边界数值失稳。谱元法中MTF 引发高频失稳的机理为在有限计算区内边界对虚假模态的反复反射放大所致；

(2) 一维波动谱元模拟中MTF 稳定条件表现为谱元参数和边界参数之间的关系，可作为多维波动稳定条件的参考；(3)通过构建高阶谱元格式和MTF 格式的向量形式，进而基于边界反射系数是分析MTF 稳定性的可行途径。关于多维弹性波动谱元模拟中MTF 的稳定性，以及其他类型内域格式和人工边界的稳定性仍需深入研究。