﻿ 透射边界高频失稳机理及稳定实现——P-SV波动
 地球物理学报  2021, Vol. 64 Issue (10): 3646-3656 PDF

Mechanism of high frequency instability and stable implementation for transmitting boundary—P-SV wave motion
ZHANG XuBin, LIAO ZhenPeng, XIE ZhiNan
Institute of Engineering Mechanics, China Earthquake Administration, Key Laboratory of Earthquake Engineering and Engineering Vibration of China Earthquake Administration, Harbin 150080, China
Abstract: The stability of wave numerical simulation is the premise of obtaining reliable results. Transmitting boundary is a kind of artificial boundary with high order and high efficiency, but there exists high frequency instability induced by the improper coupling between the interior scheme and artificial boundary scheme. Aiming at the instability problem caused by transmitting boundary in the finite element simulation of P-SV wave motion, based on the group velocity interpretation of GKS theorem, the mechanism of instability is further clarified by the dispersion analysis of finite element and transmitting boundary scheme. The mechanism is that the interior scheme and artificial boundary scheme support high frequency P wave or SV wave whose group velocity points to interior domain, thus wave energy is allowed to enter interior domain spontaneously from artificial boundary to cause instability. At the same time, compared with the dispersion of the continuous model, it shows that the high frequency wave components leading to instability are introduced by the finite element discretization. In order to eliminate the above wave components, modified integration rules are used to change the stiffness of the finite element. Theory analysis and numerical experiments show that the proposed stabilization measures can ensure transmitting boundary to be implemented stably.
Keywords: Wave simulation    High frequency instability    GKS theorem    Dispersion analysis    Numerical integration
0 引言

P-SV弹性矢量波动方程是地震波动模拟常用的计算模型，在地形效应、盆地效应以及全球地震动模拟(Li et al., 2014)都存在应用.MTF在矢量波动模拟中相比标量波动更易引发失稳，而针对标量波动模拟中得到的MTF稳定条件(章旭斌等, 2015)并不能用于更为复杂的矢量波动模拟.另外，GKS定理的群速度解释虽然是针对一维和标量波动，但其失稳模态能自然地推广到P-SV波动中的P波和SV波.本文将在以往研究基础上，将GKS定理及其群速度解释应用到P-SV波动模拟中MTF的稳定性分析，揭示MTF引发的失稳机理，并采用修改的数值积分方法(Yue and Guddati, 2005)修改有限元刚度来调整内域格式的频散，消除引发MTF失稳的波动成分，从而稳定实现MTF.同时利用长持时数值实验加以验证.

1 数值格式

 (1)

 (2)

 (3)

 (4)

 (5)

 (6)

 (7)

 (8)

 (9)

 (10)

 (11)

 (12)

 (13)

 (14)

 (15)

 图 1 长方形单元和局部节点系 Fig. 1 Rectangle element and local grid system

 (16)

 (17)

2 失稳机理及稳定实现

 (18)

z=eiωΔtκ=e－ikxΔx时，正则模态退化为傅里叶模态.数值稳定也就是要求半无限模型不支持随时间步数无限增长的正则模态(|z|＞1, |κ|＞1)，即Godunov-Ryabenkii稳定条件(Gustafsson et al., 1972).其中|κ|＞1的含义为：在内域稳定条件下，数值格式已不再支持|z|＞1，|κ|=1的模态，并考虑到正则模态在无穷远处应为有限值，因此对右边界而言，|κ|＞1.图 2给出了正则模态示意图.

 图 2 正则模态示意图 Fig. 2 Sketch of normal model

Gustafsson等(1972)指出Godunov-Ryabenkii稳定条件未考虑|z|=1，|κ|=1时的失稳谐波模态，即z由|z|＞1趋近于|z|=1，相应的κ由|κ|＞1趋近于|κ|=1的极限状态，因此要求半无限模型不支持这一失稳谐波模态，并证明了Godunov-Ryabenkii稳定条件加上这一补充条件为稳定的充要条件，这个结果称为GKS定理.Trefethen(1983)基于扰动分析，对这一失稳模态给出了具有物理直观的群速度解释.由于

 (19)

GKS定理的群速度解释物理概念清晰，易于操作.边界和内域频散曲线的交点就是两者共同支持的谐波模态，其交点切线斜率就是谐波的群速度.对于右边界而言，如果群速度为负，则表示向内行进，边界将会引发失稳.值得说明的是，GKS定理的群速度解释虽是针对一维和SH波动的，但其失稳模态能自然地推广到P-SV波动中的P波和SV波.考虑到P波和SV波的频散是完全独立的，因此分别针对P波和SV波，验证边界和内域格式是否支持群速度指向内域的失稳谐波模态.

 (20)

 (21)

 (22)

 (23)

 (24)

 (25)

 (26)

 图 3 连续模型和传统有限元的频散曲线，及在垂直x方向人工边界上的MTF频散曲线，其中参数ε=, β=1, Δτ=0.4 (a)(b) 分别为连续模型的P波和SV波频散；(c)(d) 分别为传统有限元的P波和SV波频散. Fig. 3 The dispersion curves of continuous model, classic FEM and MTF scheme when artificial boundary is perpendicular to x-axis with these parameters, ε=, β=1, Δτ=0.4 (a) and (b) The dispersion curves of P wave and SV wave in continuous model; (c) and (d) The dispersion curves of P wave and SV wave in classic FEM.

 (27)

 (28)

 (29)

 (30)

 (31)

 图 4 修正有限元和在垂直x方向人工边界上的MTF频散曲线，其中参数ε=, β=1, Δτ=0.5 (a)(b) 分别为修正有限元()的P波和SV波频散；(c)(d) 分别为修正有限元(α=1)的P波和SV波频散. Fig. 4 The dispersion curves of modified FEM and MTF scheme when artificial boundary is perpendicular to x-axis with these parameters, ε=, β=1, Δτ=0.5 (a) and (b) The dispersion curves of P wave and SV wave in modified FEM (); (c) and (d) The dispersion curves of P wave and SV wave in modified FEM (α=1).

 (32)

 (33)

3 数值实验 3.1 基岩半空间

 (34)

 图 5 基岩半空间模型 Fig. 5 Rock half-space model

 图 6 传统有限元结合MTF计算的接收点位移时程 (a) x方向；(b) y方向. Fig. 6 The displacement time history of the receiver computed by classic FEM with MTF (a) x direction; (b) y direction.
 图 7 修正有限元结合MTF计算的接收点位移时程 (a) x方向；(b) y方向. Fig. 7 The displacement time history of the receiver computed by modified FEM with MTF (a) x direction; (b) y direction.
 图 8 (a) 计算模型参数条件下，传统有限元和MTF的频散曲线，红线为一条水平频散曲线，红色方框为引发MTF失稳的谐波频率范围；(b) 截取的接收点失稳增长数值解；(c) 失稳增长数值解的频谱 Fig. 8 (a) The dispersion curves of classic FEM and MTF at the given parameters of computational model, the red line is horizontal dispersion curve, and the red box indicates the harmonic frequency range that causes MTF instability; (b) Truncated numerical instability solution of receiver; (c) Spectrum of numerical instability solution
3.2 盆地-基岩半空间

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

 图 10 修正有限元结合MTF和黏弹性边界以及远置边界时计算的接收点位移时程 (a) x方向；(b) y方向. Fig. 10 The displacement time history of the receiver computed by modified FEM with MTF, viscous-spring boundary and extend boundary (a) x direction; (b) y direction.
4 结论

P-SV弹性波动有限元模拟中MTF引发的高频失稳机理可用GKS定理群速度解释加以阐明，即有限元离散引入了在连续方程中并不存在的群速度指向内域的高频P波或SV波，并得到了MTF的支持，从而导致数值失稳.可见数值模拟中由于引入人工边界所导致的失稳并非仅是边界的问题.本文进而采用修改的数值积分方法修改有限元刚度来调整内域格式的频散，消除了引发边界失稳的波动成分，稳定实现了MTF，从而构建了稳定的近场P-SV波动模拟方案.对于实际复杂问题，通过保证边界区域为规则的长方形单元且满足稳定条件，数值实验显示MTF同样能够稳定实现.本文研究思路可进一步推广到三维弹性波动数值模拟，对于其他类型数值格式中人工边界的稳定问题亦具有参考价值.

 (A1)

 (A2)

 (A3)

 (A4)

 (A5)

 (A6)

 (A7)

(A7) 式即为(29)式的前半部分，对于β≥1这一必要条件，可采用反证法证明.如果β＜1，取值s22=，即Q2=0.为避免Q3等于0以消除(A4)式中的Q3，因此再取值s1→1，即kxΔx→π，此时.在β＜1条件下，满足s22∈[0, 1]，因此上述取值方式是存在的，从而(A4)式简化为

 (A8)

s1→1可知，c1→0，(Q′1)2→0，(Q′3)2.由于s22为0到1之间的正常数，因此(Q′3)2趋向于一个正常数，而(Q′1)2却趋向于0，故(A8)式无法成立，从而可得

 (A9)

(A7) 和(A9)式的组合即(29)式，因此(A3)和(A4)式的必要条件为(29)式.

 (A10)

 (A11)

 (A12)

Qi, Q′i代入(A11)式，经整理可知，要证明(A11)式仅需证明不等式

 (A13)

 (A14)

 (A15)

References
 Bérenger J P. 2007. Perfectly Matched Layer (PML) for Computational Electromagnetics. Arcueil: Morgan & Claypool. Bowden D C, Tsai V C. 2017. Earthquake ground motion amplification for surface waves. Geophysical Research Letters, 44(1): 121-127. DOI:10.1002/2016GL071885 Chen H Q. 2006. Discussion on seismic input mechanism at dam site. Journal of Hydraulic Engineering (in Chinese), 37(12): 1417-1423. Duru K, Kozdon J E, Kreiss G. 2015. Boundary conditions and stability of a perfectly matched layer for the elastic wave equation in first order form. Journal of Computational Physics, 303: 372-395. DOI:10.1016/j.jcp.2015.09.048 Gustafsson B, Kreiss H O, Sundström A. 1972. Stability theory of difference approximations for mixed initial boundary value problems.Ⅱ. Mathematics of Computation, 26(119): 649-686. DOI:10.1090/S0025-5718-1972-0341888-3 Hagstrom T, Mar-Or A, Givoli D. 2008. High-order local absorbing conditions for the wave equation: Extensions and improvements. Journal of Computational Physics, 227(6): 3322-3357. DOI:10.1016/j.jcp.2007.11.040 Higdon R L. 1987. Numerical absorbing boundary conditions for the wave equation. Mathematics of Computation, 49(179): 65-90. DOI:10.1090/S0025-5718-1987-0890254-1 Jing L P, Liao Z P, Zou J X. 2002. A high-frequency instability mechanism in numerical realization of multi-transmitting formula. Earthquake Engineering and Engineering Vibration (in Chinese), 22(1): 7-13. Li D Z, Helmberger D, Clayton R W, et al. 2014. Global synthetic seismograms using a 2-D finite-difference method. Geophysical Journal International, 197(2): 1166-1183. DOI:10.1093/gji/ggu050 Li X J, Yang Y. 2012. Measures for stability control of transmitting boundary. Chinese Journal of Geotechnical Engineering (in Chinese), 34(4): 641-645. Liao Z P, Liu J B. 1992. Numerical instabilities of a local transmitting boundary. Earthquake Engineering and Structural Dynamics, 21(1): 65-77. DOI:10.1002/eqe.4290210105 Liao Z P. 2002. Introduction to Wave Motion Theories in Engineering (in Chinese). 2nd ed. Beijing: Science Press. Liao Z P, Zhou Z H, Zhang Y H. 2002. Stable implementation of transmitting boundary in numerical simulation of wave motion. Chinese Journal of Geophysics (in Chinese), 45(4): 533-545. DOI:10.3321/j.issn:0001-5733.2002.04.011 Liu Y S, Xu T, Wang Y H, et al. 2019. An efficient source wavefield reconstruction scheme using single boundary layer values for the spectral element Method. Earth and Planetary Physics, 3(4): 342-357. DOI:10.26464/epp2019035 Liu J B, Li B. 2005. A unified viscous-spring artificial boundary for 3-D static and dynamic applications. Science in China Series E Engineering & Materials Science, 48(5): 570-584. Liu Y, Sen M K. 2012. A hybrid absorbing boundary condition for elastic staggered-grid modelling. Geophysical Prospecting, 60(6): 1114-1132. DOI:10.1111/j.1365-2478.2011.01051.x Semblat J F, Lenti L, Gandomzadeh A. 2011. A simple multi-directional absorbing layer method to simulate elastic wave propagation in unbounded domains. International Journal for Numerical Methods in Engineering, 85(12): 1543-1563. DOI:10.1002/nme.3035 Strikwerda J C. 2004. Finite Difference Schemes and Partial Differential Equations. 2nd ed. Philadelphia: SIAM. Taflove A, Hagness S C. 2005. Computational Electrodynamics: The Finite-Difference Time-Domain Method. 3rd ed. Norwood, MA: Artech House, Inc. Trefethen L N. 1983. Group velocity interpretation of the stability theory of Gustafsson, Kreiss, and Sundström. Journal of Computational Physics, 49(2): 199-217. DOI:10.1016/0021-9991(83)90123-7 Xie Z N, Liao Z P. 2012. Mechanism of high frequency instability caused by transmitting boundary and method of its elimination——SH wave. Chinese Journal of Theoretical and Applied Mechanics (in Chinese), 44(4): 745-752. Xu G, Hamouda A M S, Khoo B C. 2016. Time-domain simulation of second-order irregular wave diffraction based on a hybrid water wave radiation condition. Applied Mathematical Modelling, 40(7-8): 4451-4467. DOI:10.1016/j.apm.2015.11.034 Yue B, Guddati M N. 2005. Dispersion-reducing finite elements for transient acoustics. The Journal of the Acoustical Society of America, 118(4): 2132-2141. DOI:10.1121/1.2011149 Zhang X B, Liao Z P, Xie Z N. 2015. Mechanism of high frequency coupling instability and stable implementation for transmitting boundary-SH wave motion. Chinese Journal of Geophysics (in Chinese), 58(10): 197-206. DOI:10.6038/cjg20151017 Zhang X B, Xie Z N, Liao Z P. 2019. Reflection amplification instability of the transmitting boundary in SH wave simulation. Journal of Harbin Engineering University (in Chinese), 40(6): 1031-1035. Zhao M, Du X L, Liu J B, et al. 2011. Explicit finite element artificial boundary scheme for transient scalar waves in two-dimensional unbounded waveguide. International Journal for Numerical Methods in Engineering, 87(11): 1074-1104. DOI:10.1002/nme.3147 Zhou Z H, Liao Z P. 2001. A measure for eliminating drift instability of the multi-transmitting formula. Acta Mechanica Sinica (in Chinese), 33(4): 550-554. 陈厚群. 2006. 坝址地震动输入机制探讨. 水利学报, 37(12): 1417-1423. DOI:10.3321/j.issn:0559-9350.2006.12.004 景立平, 廖振鹏, 邹经湘. 2002. 多次透射公式的一种高频失稳机制. 地震工程与工程振动, 22(1): 7-13. DOI:10.3969/j.issn.1000-1301.2002.01.002 李小军, 杨宇. 2012. 透射边界稳定性控制措施探讨. 岩土工程学报, 34(4): 641-645. 廖振鹏. 2002. 工程波动理论导论. 2版. 北京: 科学出版社. 廖振鹏, 周正华, 张艳红. 2002. 波动数值模拟中透射边界的稳定实现. 地球物理学报, 45(4): 533-545. DOI:10.3321/j.issn:0001-5733.2002.04.011 谢志南, 廖振鹏. 2012. 透射边界高频失稳机理及其消除方法——SH波动. 力学学报, 44(4): 745-752. DOI:10.6052/0459-1879-11-312 章旭斌, 廖振鹏, 谢志南. 2015. 透射边界高频耦合失稳机理及稳定实现——SH波动. 地球物理学报, 58(10): 197-206. DOI:10.6038/cjg20151017 章旭斌, 谢志南, 廖振鹏. 2019. SH波动模拟中透射边界反射放大失稳研究. 哈尔滨工程大学学报, 40(6): 1031-1035. 周正华, 廖振鹏. 2001. 消除多次透射公式飘移失稳的措施. 力学学报, 33(4): 550-554. DOI:10.3321/j.issn:0459-1879.2001.04.015