波动方程数值模拟的边界条件研究.docx

波动方程数值模拟的边界条件研究陈国新 摘要：本文主要介绍了在利用有限差分法进行波动方程数值模拟时为消除人工边界产生的虚假反射而采用的旁轴近似法,多向透射边界条件，黏性边界，衰减吸收边界，完全匹配层法及其相应的改进。通过分析和比较这些吸收边界条件各自的优缺点来更直观的阐述吸收边界条件消除人工边界影响的基本原理，为今后实际工作中的吸收边界条件的应用和改进提供理论上帮助。关键词：数值模拟 边界条件 旁轴近似 海绵吸收 完全匹配层引言波动方程尤其是地震波场的数值模拟在研究地震波的传播规律，指导设计数据的观测系统，验证地震解释，为数据处理提供理论数据，制定试验处理流程，为反演提供正演方法等方面具有广泛的应用。但是由于计算模型的限制，在正演模拟中一般要人为的限定数值模拟范围，这样就不可避免的引入了了人工边界，而如何消除人工边界对地震波的虚假反射则成为地震波数值模拟中的重要问题。为消除人工边界带来的影响基于不同原理的吸收边界条件相应提出。按照其思想划分主要分为两种：一种是基于单程波近似的吸收边界条件，如旁轴近似法，多向透射边界条件。一种是基于加吸收层的衰减边界条件，如黏性边界，海绵吸收边界条件，完全匹配层吸收边界条件。.第一类方法不严格满足散射波向无穷远传播的波动性质， 即不满足向无穷远逸射的辐射定律， 而是近似的，具有计算简单，运算时间短的特点，但是其效果依赖于波的入射角度，而第二类方法相对于第一类方法计算较为复杂，但是吸收效果较好，实用性更强。一 基于单程波近似原理的边界条件旁轴近似法最早是由Robert Clayton 和 Bjorn Engquist于1977年提出[1],其主要思想是在于构造一个波只能沿一个方向传播的方程， 以代替波动方程， 并在小入射角前提下求替代方程的截断级数解.以二维的标量声波方程为例： + (1)对方程变换（1）做傅氏变换就可得到频率波数域方程 (2)我们考虑方程（1）的空间外推，比如说沿Z方向，那么我们可以将方程（2）改成如下的频散关系式： =(/v) (3)方程（3）中的正号表示波场沿Z轴正方向传播，即为“上行波”，负号表示波场沿Z轴负方向传播，即为“下行波”。为了保持算法的稳定性我们将不考虑（3）式中|/w|1所对应的损耗波成分。当波场沿Z轴正方向传播时有方程（3）可得 /w= （4）我们对方程（4）的右边做合理的近似，一般有三种近似方法： a1： /w=1+0() (5) a2：/w=1-+0() (6)a3：/w= +0() （7） 将(5),(6),(7)式做傅里叶反变换就可得到它们的微分形式的近似： A1: +（1/v）=0 (8) A2:+（1/v）-(v/2)=0 (9) A3:-(/4)+(1/v)-(3v/4)=0 (10)旁轴近似法吸收人工边界反射的原理就是通过对方程（3）做合理化近似，然后在下边界利用下行波方程，上边界利用上行波方程从而达到吸收边界反射的目的。将（3）式中与的位置互换就可得到相应的左行波方程与右行波方程。同样的道理，对左边界与右边界分别使用左行波方程和右行波方程吸收边界反射。对人工边界反射的吸收效果可以用边界处的有效反射系数r的大小来衡量，以一个沿Z轴正方向传播的平面波为例，平面波方程为： =exp(+i-iwt) (11)则在底边界引起的反射波： =r*exp(i-i-iwt) (12) 其中的r就是底边界上的有效反射系数。在边界附近，波场+既满足边界条件，也满足模拟区域内部的波动方程，并且波场值在边界上为零，以A2近似方法为例可得： w(-r)-(-)（-r）=0 （13） 求解方程（13）可得r: r= (14)设平面波入射方向与底边界法线方向夹角为，设/w=cos,结合方程（4）可以得到有效反射系数与角的函数关系式： r() = - (15)由[0,],可以看出|r()|关于是增函数，即表明随着角的增大有效反射系数变大，也就意味着吸收边界的吸收效果变差。从图一可以看出利用三种近似方法时有效反射系数r与角大小联系。可见旁轴近似法对边界反射的吸收效果依赖于入射角。 图一 采用旁轴近似法的入射波入射角度与有效反射系数关系图Clayton 和Engquist[2]对针上述算法做了一些改进，加强旁轴近似法的实用性 。仍以上文提到的沿Z轴正向的传播的平面波为例。在底边界反射波为，在这里引入一个线性算子B,使它作用于边界上的波场时，波场值为0，即： B(+)=0 (16)基于算子B是线性算子，所以对方程（16）整理可得： r=B/B() （17） 因为与是关于,,w的坡场函数，所以（17）式可以改写成： r=B(,w)/B(,w) （18） 从（18