石油物探  2019, Vol. 58 Issue (4): 499-508  DOI: 10.3969/j.issn.1000-1441.2019.04.003
0
文章快速检索     高级检索

引用本文 

覃发兵, 高志伟, 解皓楠, 等. 完全匹配层在时域有限元弹性波数值模拟中的应用[J]. 石油物探, 2019, 58(4): 499-508. DOI: 10.3969/j.issn.1000-1441.2019.04.003.
QIN Fabing, GAO Zhiwei, XIE Haonan, et al. Application of perfectly matched layer in time-domain finite-element elastic wave modeling[J]. Geophysical Prospecting for Petroleum, 2019, 58(4): 499-508. DOI: 10.3969/j.issn.1000-1441.2019.04.003.

基金项目

湖北省教育厅科技研究计划青年人才项目(Q20161304)及油气资源与勘探技术教育部重点实验室(长江大学)开放基金(Q20161304)共同资助

作者简介

覃发兵(1979—), 男, 博士, 讲师, 现从事地震波场正演模拟、测井数据处理与软件开发等研究工作。Email:qinfabing@aliyun.com

文章历史

收稿日期:2018-05-29
改回日期:2018-11-02
完全匹配层在时域有限元弹性波数值模拟中的应用
覃发兵1,2 , 高志伟3 , 解皓楠3 , 徐振旺4     
1. 长江大学地球物理与石油资源学院, 湖北武汉 430100;
2. 长江大学管理学院, 湖北荆州 434023;
3. 中国石油天然气股份有限公司新疆油田分公司采油一厂, 新疆克拉玛依 834000;
4. 中国石油天然气股份有限公司辽河油田分公司勘探开发研究院, 辽宁盘锦 124010
摘要:完全匹配层(perfectly matched layer, PML)边界条件是消除人工边界虚假反射的经典方法之一, 但不易在时域有限元方法中实现, 尤其是求解二阶弹性波方程。为此, 详细推导了PML在时域有限元法求解二阶位移弹性波方程中的加载过程, 得到了含PML的有限元控制方程; 通过数值算例, 讨论了PML衰减参数中理论反射系数对PML吸收效果的影响以及在PML吸收层最外层加载狄利克雷边界条件对PML数值稳定性的影响。数值模拟结果表明:当PML吸收层厚度一定时, 理论反射系数越小, PML吸收效果越好; 当PML吸收层厚度为半个最大主波长时, 理论反射系数小于(等于)10-5, PML吸收效果最优; 虽然在PML吸收层的最外层加载狄利克雷边界条件可增强PML的数值稳定性, 但对处于自由表面上的PML吸收层最外层部分, 不可加载狄利克雷边界条件, 否则会产生严重的虚假反射。
关键词完全匹配层    弹性波    数值模拟    边界条件    狄利克雷条件    
Application of perfectly matched layer in time-domain finite-element elastic wave modeling
QIN Fabing1,2, GAO Zhiwei3, XIE Haonan3, XU Zhenwang4     
1. Geophysics and Oil Resource Institute, Yangtze University, Wuhan 430100, China;
2. School of Management, Yangtze University, Jingzhou 434023, China;
3. No.1 Production Plant, Xinjiang Oilfield Company, PetroChina, Karamay 834000, China;
4. Research Institute of Petroleum Exploration and Development, Liaohe Oilfield Company, CNPC, Panjin 124010, China
Foundation item: This research is financially supported by the Science and Technology Research Program Youth Talent Project of Hubei Provincial Department of Education (Grant No.Q20161304), and Open Fund (Yangtze University) from Key Laboratory of Oil and Gas Resources and Exploration Technology of Ministry of Education (Grant No.Q20161304)
Abstract: The perfectly matched layer (PML) is traditionally used to suppress artificial boundary reflections.It is difficult to apply the PML using a time-domain finite-element method, especially when solving second-order elastic wave equations.We proposed a detailed derivation that can be used to embed the PML into time-domain finite-element formulas and solve second-order displacement elastic wave equations; the final objective would be to obtain finite-element governing equations that include the PML.We also discussed the influence of a theoretical reflection coefficient on the PML, and the effects of Dirichlet boundary conditions (which were imposed on the outermost PML) on the PML numerical stability.The results showed that when the thickness of the absorption layer was constant, the theoretical reflection coefficient was small and the absorption effect of the PML was improved.On the other hand, when the thickness of the absorption layer was equal to half of the maximum dominant wavelength, the PML presented optimal absorbing properties only if the theoretical reflection coefficient was less than or equal to 10-5.The numerical stability of the PML could be enhanced imposing Dirichlet boundary conditions on its outermost layer; however, these conditions could not be imposed on the free surface of the outermost part of the PML without running the risk of serious false reflection.
Keywords: perfectly matched layer (PML)    elastic waves    numerical modeling    boundary condition    Dirichlet condition    

高精度弹性波数值模拟是弹性波动方程成像和反演的核心步骤之一。弹性波数值模拟时, 从计算可行性考虑, 无界的地下介质需加人工边界截断为有界区域, 但在人工边界处会出现虚假反射, 从而严重影响数值模拟的精度。因此, 需要在人工边界处进行特殊处理以消除虚假反射对数值模拟结果的影响。

根据地震波场数值模拟时是否在求解区域最外层附加额外的几何空间, 可将近四十年发展起来的人工边界处理技术分为两类:吸收边界方法和吸收衰减层方法。吸收边界方法是在计算区域最外层网格节点上引入额外变量压制虚假反射。由于没有增加计算区域的几何空间, 吸收边界方法产生的额外计算量较少。经典的吸收边界方法为CLAYTON等[1]于1977年提出的旁轴近似法, 该方法能有效吸收垂直入射的波场能量, 但对于大角度入射能量吸收效果较差, 难以满足精度要求。针对该问题, HIGDON[2]将旁轴近似法推广为高阶形式, 使其能吸收更大范围角度入射的波场能量, 但高阶形式的旁轴近似实现过程较复杂, 且引入的额外计算量较大。不同于吸收边界方法, 吸收衰减层方法是在实际物理求解区域外额外设置多个衰减层以逐渐吸收传入到衰减层中的波场能量。该方法显然会引入额外的计算量, 但往往能吸收较大范围入射角的波场能量。1994年, BERENGER[3]在求解Maxwell方程组时提出了目前最经典的吸收衰减层方法——完全匹配层(PML)。该方法理论上能使任意角度任意频率成分的边界入射波的反射系数为零, 并因此而得名。1996年, CHEW等[4]首次将PML推广到弹性波方程的求解。早期PML需要分裂波场来求解, 这对于大规模数值模拟是不容忽视的计算量。2007年, KOMATITSCH和MARTIN在地震波场模拟中发展了不需要分裂波场的卷积完全匹配层(CPML)[5], 显著提高了PML的计算效率及对近掠射的边界入射波的吸收效果。随后, CPML技术被迅速推广到声波方程模拟及更为复杂介质的数值模拟中[6-9]。近期, 高正辉等[10]在2.5维声波近似方程数值模拟中提出了一种改进的PML方法, 该方法可使边界入射波在PML吸收层中进行两次能量衰减, 从而进一步提高了PML的吸收效果。

然而, 关于PML的应用研究大多在有限差分法数值模拟方面[10-15], 这是因为PML在一阶形式的波动方程求解中较容易加载。有限元法数值模拟主要基于二阶波动方程, 对于PML的应用存在一定的障碍。2003年, KOMATITSCH等[16]提出了一组适用于二阶波动方程的PML方程组, 并将其成功应用于谱元法地震波场数值模拟。尽管如此, 关于PML在有限元法数值模拟中的应用研究仍然较少[17-18]

本文基于LIU等[19]提出的PML在谱元法中的加载方法, 推导了分裂式PML在二维弹性波时间域有限元数值模拟中的加载过程, 说明了理论反射系数选取对PML吸收效果的影响, 讨论了狄利克雷边界在PML外边界层和自由边界处的正确使用方式及其对PML吸收效果的影响。

1 PML的基本原理

频率域各向异性弹性波动方程可表示为:

$ -\rho \omega^{2} \boldsymbol{U}=\nabla \cdot(\boldsymbol{c} : \nabla \boldsymbol{U}) $ (1)

式中:ρ为密度, ω表示圆频率; U=U(x, ω)为频率域位移矢量, x= $x \boldsymbol{\boldsymbol{\hat{x}}}+y \boldsymbol{\hat{y}}+z \boldsymbol{\hat{z}}$为位置矢量; c为具有21个独立参数的弹性张量; “:”表示张量双点积运算, “ ·”为散度算子。在均匀各向同性条件下, (1)式具有简谐平面波解, 其一般形式为A exp{-i(k·xωt}。其中A表示振幅和平面波极化的方向, k= $k_{x} \boldsymbol{\hat{x}}+k_{y} \boldsymbol{\hat{y}}+k_{z} \boldsymbol{\hat{z}}$为具有笛卡尔分量的波数矢量。平面P波满足A×k= 0, 且波数k=(kx2+ky2+kz2)1/2=ω/α, 这里α表示P波波速; 平面S波满足A·k=0, 且k=ω/β, 这里β表示S波波速。

图 1所示, PML边界层区域位于x>0的范围, 而求解的目标区域位于x≤0的范围。首先需要定义阻尼函数dx(x), 在目标区域内, dx(x)=0;在PML区域内, dx(x)>0。注意这里的下标x指的是x坐标轴, 括号内的x为实变量。PML的核心思想是依据阻尼函数dx(x)构建一个复数坐标${\tilde{x}}$, 位置矢量变为x= $\tilde{x} \boldsymbol{\hat{x}}+y \boldsymbol{\hat{y}}+z \boldsymbol{\hat{z}}$。复数坐标的表达式[5]为:

图 1 目标区域和PML边界层区域图示
$ \tilde{x}(x)=x-\frac{\mathrm{i}}{\omega} \int_{0}^{\overline{x}} d_{x}(s) \mathrm{d} s $ (2)

依据链式法有:

$ \partial_{x}=\partial_{\tilde{x}} \frac{\mathrm{d} \tilde{x}}{\mathrm{d} x}=\frac{\mathrm{i} \omega+d_{x}}{\mathrm{i} \omega} \partial_{\tilde{x}} $ (3)

x表示对x坐标的偏导(下文同)。则:

$ \partial_{\tilde{x}}=\frac{i \omega}{\mathrm{i} \omega+d_{x}} \partial_{x}=\frac{1}{s_{x}} \partial_{x} $ (4)

其中

$ s_{x}=\frac{\mathrm{i} \omega+d_{x}}{\mathrm{i} \omega}=1+\frac{d_{x}}{\mathrm{i} \omega} $ (5)

PML从纯数学上讲, 是将对实坐标的偏导做如(4)式所示对复数坐标偏导的复数变换。对应地, PML的物理意义可借助平面波解来很好地理解。将(2)式代入平面波解:

$ \boldsymbol{A} \mathrm{e}^{-\mathrm{i}\left(k_{x} \tilde{x}+k_{y} y+k_{z} z-\omega t\right)}=\boldsymbol{A} \mathrm{e}^{-(k \cdot x-\omega t)} \mathrm{e}^{-k_{x} / \omega \int_{0}^{\overline{x}} d_{x}(s) \mathrm{d} s} $ (6)

(6) 式可理解为PML吸收人工边界反射的基本原理。在目标区域, 由于dx=0, 解的形式与原来一致。但在PML区域, 从PML内边界到PML外边界(x≥0方向)dx逐渐增大, 衰减系数exp{kx/ω·∫ 0xdx(s)ds}逐渐增大, 从而逐渐衰减了向人工边界传播的波场。这里的衰减系数与平面波圆频率ω成反比, 即低频成分比高频成分衰减更快。

2 含PML的有限元控制方程 2.1 含PML的弹性波方程

以二维各向同性完全弹性介质为例, 其本构方程为:

$ \left[\begin{array}{c}{\sigma_{x x}} \\ {\sigma_{z z}} \\ {\sigma_{x z}}\end{array}\right]=\left[\begin{array}{ccc}{C_{11}} & {C_{12}} & {0} \\ {C_{12}} & {C_{22}} & {0} \\ {0} & {0} & {C_{66}}\end{array}\right]\left[\begin{array}{c}{\partial_{x} U_{x}} \\ {\partial_{z} U_{z}} \\ {\partial_{x} U_{z}+\partial_{z} U_{x}}\end{array}\right] $ (7)

式中:σxxσzzσxz为应力分量; 弹性参数C11=C22=λ+2μ, C12=λ, C66=μ, λμ为拉梅常数; Uj(j=x, z)约定为频率域位移场量。纵波速度vP和横波速度vS可分别表示为vP= $\sqrt{(\lambda+2 \mu) / \rho}$vS= $\sqrt{\mu / \rho}$。频率域运动方程可表示为:

$ \left\{\begin{array}{l}{\rho(\mathrm{i} \omega)^{2} U_{x}=\partial_{x} \sigma_{x x}+\partial_{z} \sigma_{x z}} \\ {\rho(\mathrm{i} \omega)^{2} U_{z}=\partial_{x} \sigma_{x z}+\partial_{z} \sigma_{z z}}\end{array}\right. $ (8)

将(7)式代入(8)式, 有:

$ \left\{\begin{array}{cl}{\rho(\mathrm{i} \omega)^{2} U_{x}=} & {\partial_{x}\left(C_{11} \partial_{x} U_{x}+C_{12} \partial_{z} U_{z}\right)+} \\ {} & {\partial_{z}\left[C_{66}\left(\partial_{x} U_{z}+\partial_{z} U_{x}\right)\right]} \\ {\rho(\mathrm{i} \omega)^{2} U_{z}=} & {\partial_{x}\left[C_{66}\left(\partial_{x} U_{z}+\partial_{z} U_{x}\right)\right]+} \\ {} & {\partial_{z}\left(C_{12} \partial_{x} U_{x}+C_{22} \partial_{z} U_{z}\right)}\end{array}\right. $ (9)

根据前文PML原理, 加载PML的数学本质是将实坐标系变换到复坐标系, 即有对应关系∂x$\partial_{\tilde{x}}$和∂z$\partial_{\tilde{z}}$, 则:

$ \left\{ \begin{array}{l} \rho {\left( {{\rm{i}}\omega } \right)^2}{U_x} = {\partial _{\tilde x}}\left( {{C_{11}}{\partial _{\tilde x}}{U_x} + {C_{12}}{\partial _{\tilde z}}{U_z}} \right) + \\ \;\;\;\;\;\;\;\;\;\;{\partial _{\tilde z}}\left[ {{C_{66}}\left( {{\partial _{\tilde x}}{U_z} + {\partial _{\tilde z}}{U_x}} \right)} \right]\\ \rho {\left( {{\rm{i}}\omega } \right)^2}{U_z} = {\partial _{\tilde x}}\left[ {{C_{66}}\left( {{\partial _{\tilde x}}{U_z} + {\partial _{\tilde z}}{U_x}} \right)} \right] + \\ \;\;\;\;\;\;\;\;\;\;{\partial _{\tilde z}}\left( {{C_{12}}{\partial _{\tilde x}}{U_x} + {C_{22}}{\partial _{\tilde z}}{U_z}} \right) \end{array} \right. $ (10)

又知:

$ \left\{\begin{array}{l}{\partial_{\tilde{x}}=\frac{\mathrm{i} \omega}{\mathrm{i} \omega+d_{x}} \partial_{x}} \\ {\partial_{\tilde{z}}=\frac{\mathrm{i} \omega}{\mathrm{i} \omega+d_{z}} \partial_{z}}\end{array}\right. $ (11)

将(11)式代入(10)式得:

$ \left\{ \begin{array}{l} \rho {\left( {{\rm{i}}\omega } \right)^2}{U_x} = \frac{{{\rm{i}}\omega }}{{{\rm{i}}\omega + {d_x}}}{\partial _x}\left( {{C_{11}}\frac{{{\rm{i}}\omega }}{{{\rm{i}}\omega + {d_x}}}{\partial _x}{U_x} + } \right.\\ \;\;\;\;\;\;\left. {{C_{12}}\frac{{{\rm{i}}\omega }}{{{\rm{i}}\omega + {d_z}}}{\partial _z}{U_z}} \right) + \frac{{{\rm{i}}\omega }}{{{\rm{i}}\omega + {d_z}}}{\partial _z}\left[ {{C_{66}}\left( {\frac{{{\rm{i}}\omega }}{{{\rm{i}}\omega + {d_x}}} \cdot } \right.} \right.\\ \;\;\;\;\;\;\left. {\left. {{\partial _x}{U_z} + \frac{{{\rm{i}}\omega }}{{{\rm{i}}\omega + {d_z}}}{\partial _z}{U_x}} \right)} \right]\\ \rho {\left( {{\rm{i}}\omega } \right)^2}{U_z} = \frac{{{\rm{i}}\omega }}{{{\rm{i}}\omega + {d_x}}}{\partial _x}\left[ {{C_{66}}\left( {\frac{{{\rm{i}}\omega }}{{{\rm{i}}\omega + {d_x}}}{\partial _x}{U_z} + } \right.} \right.\\ \;\;\;\;\;\;\left. {\left. {\frac{{{\rm{i}}\omega }}{{{\rm{i}}\omega + {d_z}}}{\partial _z}{U_x}} \right)} \right] + \frac{{{\rm{i}}\omega }}{{{\rm{i}}\omega + {d_z}}}{\partial _z}\left( {{C_{12}}\frac{{{\rm{i}}\omega }}{{{\rm{i}}\omega + {d_x}}} \cdot } \right.\\ \;\;\;\;\;\;\left. {{\partial _x}{U_x} + {C_{22}}\frac{{{\rm{i}}\omega }}{{{\rm{i}}\omega + {d_z}}}{\partial _z}{U_z}} \right) \end{array} \right. $ (12)

将位移场x分量Ux分裂为Ux=Ux1+Ux2+Ux3, 并将其代入(12)式中的第一个方程, 整理后有:

$ \left\{ \begin{array}{l} \rho {U_{x1}} = \frac{1}{{{{\left( {{\rm{i}}\omega + {d_x}} \right)}^2}}}{\partial _x}\left( {{C_{11}}{\partial _x}{U_x}} \right) - \frac{{{{d'}_x}}}{{{{\left( {{\rm{i}}\omega + {d_x}} \right)}^3}}} \cdot \\ \;\;\;\;\;\;\;\;\;{C_{11}}{\partial _x}{U_x}\\ \rho {U_{x2}} = \frac{1}{{\left( {{\rm{i}}\omega + {d_x}} \right)\left( {{\rm{i}}\omega + {d_z}} \right)}}{\partial _x}\left( {{C_{12}}{\partial _z}{U_z}} \right) + \\ \;\;\;\;\;\;\;\;\;\;\;\frac{1}{{\left( {{\rm{i}}\omega + {d_x}} \right)\left( {{\rm{i}}\omega + {d_z}} \right)}}{\partial _z}\left( {{C_{66}}{\partial _x}{U_z}} \right)\\ \rho {U_{x3}} = \frac{1}{{{{\left( {{\rm{i}}\omega + {d_z}} \right)}^2}}}{\partial _z}\left( {{C_{66}}{\partial _z}{U_x}} \right) - \frac{{{{d'}_z}}}{{{{\left( {{\rm{i}}\omega + {d_z}} \right)}^3}}} \cdot \\ \;\;\;\;\;\;\;\;\;{C_{66}}{\partial _z}{U_x} \end{array} \right. $ (13)

同理, 将Uz分裂为Uz=Uz1+Uz2+Uz3, 并将其代入(12)式中第二个方程, 整理后有:

$ \left\{ \begin{array}{l} \rho {U_{z1}} = \frac{1}{{{{\left( {{\rm{i}}\omega + {d_x}} \right)}^2}}}{\partial _x}\left( {{C_{66}}{\partial _x}{U_z}} \right) - \frac{{d_x^\prime }}{{{{\left( {{\rm{i}}\omega + {d_x}} \right)}^3}}}{C_{66}}{\partial _x}{U_z}\\ \rho {U_{z2}} = \frac{1}{{\left( {{\rm{i}}\omega + {d_x}} \right)\left( {{\rm{i}}\omega + {d_z}} \right)}}{\partial _x}\left( {{C_{66}}{\partial _z}{U_x}} \right) + \\ \;\;\;\;\;\;\;\;\;\;\frac{1}{{\left( {{\rm{i}}\omega + {d_x}} \right)\left( {{\rm{i}}\omega + {d_z}} \right)}}{\partial _z}\left( {{C_{12}}{\partial _x}{U_x}} \right)\\ \rho {U_{z3}} = \frac{1}{{{{\left( {{\rm{i}}\omega + {d_z}} \right)}^2}}}{\partial _z}\left( {{C_{22}}{\partial _z}{U_z}} \right) - \frac{{d_z^\prime }}{{{{\left( {{\rm{i}}\omega + {d_z}} \right)}^3}}}{C_{22}}{\partial _z}{U_z} \end{array} \right. $ (14)

(13) 式和(14)式中均含有高次项, 故引入以下4个中间变量:

$ \left\{ \begin{array}{l} {P_{xx}} = - \frac{1}{\rho }\frac{1}{{{\rm{i}}\omega + {d_x}}}d_x^\prime {C_{11}}{\partial _x}{U_x}\\ {P_{xz}} = - \frac{1}{\rho }\frac{1}{{{\rm{i}}\omega + {d_z}}}d_z^\prime {C_{66}}{\partial _z}{U_x}\\ {P_{zx}} = - \frac{1}{\rho }\frac{1}{{{\rm{i}}\omega + {d_x}}}d_x^\prime {C_{66}}{\partial _x}{U_z}\\ {P_{zz}} = - \frac{1}{\rho }\frac{1}{{{\rm{i}}\omega + {d_z}}}d_z^\prime {C_{22}}{\partial _z}{U_z} \end{array} \right. $ (15)

将(15)式代入(13)式和(14)式中, 有:

$ \left\{ \begin{array}{l} \rho {\left( {{\rm{i}}\omega + {d_x}} \right)^2}{U_{x1}} = {\partial _x}\left( {{C_{11}}{\partial _x}{U_x}} \right) + \rho {P_{xx}}\\ \rho \left( {{\rm{i}}\omega + {d_x}} \right)\left( {{\rm{i}}\omega + {d_z}} \right){U_{xz}} = {\partial _x}\left( {{C_{12}}{\partial _z}{U_z}} \right) + \\ \;\;\;\;\;\;\;\;\;\;{\partial _z}\left( {{C_{66}}{\partial _x}{U_z}} \right)\\ \rho {\left( {{\rm{i}}\omega + {d_z}} \right)^2}{U_{x3}} = {\partial _z}\left( {{C_{66}}{\partial _z}{U_x}} \right) + \rho {P_{xz}}\\ {P_{xx}} = - \frac{1}{\rho }\frac{1}{{{\rm{i}}\omega + {d_x}}}d_x^\prime {C_{11}}{\partial _x}{U_x}\\ {P_{xz}} = - \frac{1}{\rho }\frac{1}{{{\rm{i}}\omega + {d_z}}}d_z^\prime {C_{66}}{\partial _z}{U_x} \end{array} \right. $ (16)

$ \left\{ {\begin{array}{*{20}{l}} {\rho {{\left( {{\rm{i}}\omega + {d_x}} \right)}^2}{U_{z1}} = {\partial _x}\left( {{C_{66}}{\partial _x}{U_z}} \right) + \rho {P_{zx}}}\\ {\rho \left( {{\rm{i}}\omega + {d_x}} \right)\left( {{\rm{i}}\omega + {d_z}} \right){U_{z2}} = {\partial _x}\left( {{C_{66}}{\partial _z}{U_x}} \right) + }\\ \;\;\;\;\;\;{{\partial _z}\left( {{C_{12}}{\partial _x}{U_x}} \right)}\\ {\rho {{\left( {{\rm{i}}\omega + {d_z}} \right)}^2}{U_{z3}} = {\partial _z}\left( {{C_{22}}{\partial _z}{U_z}} \right) + \rho {P_{zz}}}\\ {{P_{zx}} = - \frac{1}{\rho }\frac{1}{{{\rm{i}}\omega + {d_x}}}d_x^\prime {C_{66}}{\partial _x}{U_z}}\\ {{P_{zz}} = - \frac{1}{\rho }\frac{1}{{{\rm{i}}\omega + {d_z}}}{d^\prime }z{C_{22}}{\partial _z}{U_z}} \end{array}} \right. $ (17)

对(16)式和(17)式进行时间傅里叶反变换, 可得包含了PML的弹性波动方程:

$ \left\{\begin{array}{l}{\rho\left(\partial_{t}+d_{x}\right)^{2} u_{x 1}=\partial_{x}\left(C_{11} \partial_{x} u_{x}\right)+\rho p_{x x}} \\ {\rho\left(\partial_{t}+d_{x}\right)\left(\partial_{t}+d_{z}\right) u_{x 2}=\partial_{x}\left(C_{12} \partial_{z} u_{z}\right)+} \\ \;\;\;\;\;\;{\partial_{z}\left(C_{66} \partial_{x} u_{z}\right)} \\ {\rho\left(\partial_{t}+d_{z}\right)^{2} u_{x 3}=\partial_{z}\left(C_{66} \partial_{z} u_{x}\right)+\rho p_{x z}} \\ {\rho\left(\partial_{t}+d_{x}\right) p_{x x}=-d_{x}^{\prime} C_{11} \partial_{x} u_{x}} \\ {\rho\left(\partial_{t}+d_{z}\right) p_{x z}=-d_{z}^{\prime} C_{66} \partial_{z} u_{x}}\end{array}\right. $ (18)

$ \left\{\begin{array}{l}{\rho\left(\partial_{t}+d_{x}\right)^{2} u_{z 1}=\partial_{x}\left(C_{66} \partial_{x} u_{z}\right)+\rho p_{z x}} \\ {\rho\left(\partial_{t}+d_{x}\right)\left(\partial_{t}+d_{z}\right) u_{z 2}=\partial_{x}\left(C_{66} \partial_{z} u_{x}\right)+} \\ \;\;\;\;\;\;{\partial_{z}\left(C_{12} \partial_{x} u_{x}\right)} \\ {\rho\left(\partial_{t}+d_{z}\right)^{2} u_{z 3}=\partial_{z}\left(C_{22} \partial_{z} u_{z}\right)+\rho p_{z z}} \\ {\rho\left(\partial_{t}+d_{x}\right) p_{z x}=-d_{x}^{\prime} C_{66} \partial_{x} u_{z}} \\ {\rho\left(\partial_{t}+d_{z}\right) p_{z z}=-d^{\prime}_{z} C_{22} \partial_{z} u_{z}}\end{array}\right. $ (19)

式中:uj(j=x, z)表示对应Uj(j=x, z)的时间域位移场量, pxxpxzpzxpzz为时间域中间变量。

2.2 弹性波方程的等效积分弱形式

Ω为全区域, Γ为PML外边界, nx为单位外法向量的x分量, nz为单位外法向量的z分量。采用Galerkin法, 取试探函数φxφz分别为位移x分量和z分量的变分, 将(18)式和(19)式分别乘以φxφz并求其加权积分。根据分部积分性质有:

$ \begin{array}{l} \int_\mathit{\Omega } \rho \left( {{\partial _{tt}}{u_{x1}} + 2{d_x}{\partial _t}{u_{x1}}} \right){\varphi _x}{\rm{d}}\mathit{\Omega } = \int_\mathit{\Gamma } {{C_{11}}} {\partial _x}{u_x}{\varphi _x}{n_x}{\rm{d}}\mathit{\Gamma } - \\ \int_\mathit{\Omega } {\left( {{C_{11}}{\partial _x}{u_x}{\partial _x}{\varphi _x} - \rho {p_{xx}}{\varphi _x}} \right)} {\rm{d}}\mathit{\Omega } - \int_\mathit{\Omega } \rho d_x^2{u_{x1}}{\varphi _x}{\rm{d}}\mathit{\Omega }\\ \int_\mathit{\Omega } \rho \left[ {{\partial _{tt}}{u_{x2}} + \left( {{d_x} + {d_z}} \right){\partial _t}{u_{x2}}} \right]{\varphi _x}{\rm{d}}\mathit{\Omega } = \int_\mathit{\Gamma } {\left( {{C_{12}}{\partial _z}{u_z}} \right.} \\ {n_x}{\varphi _x} + {C_{66}}{\partial _x}{u_z}{n_z}{\varphi _x}){\rm{d}}\mathit{\Gamma } - \\ \int_\mathit{\Omega } {\left( {{C_{12}}{\partial _z}{u_z}{\partial _x}{\varphi _x} + {C_{66}}{\partial _x}{u_z}{\partial _z}{\varphi _x}} \right)} {\rm{d}}\mathit{\Omega } - \int_\mathit{\Omega } \rho {d_x}{d_z}\\ {u_{x2}}{\varphi _x}{\rm{d}}\mathit{\Omega }\\ \int_\mathit{\Omega } \rho \left( {{\partial _{tt}}{u_{x3}} + 2{d_z}{\partial _t}{u_{x3}}} \right){\varphi _x}{\rm{d}}\mathit{\Omega } = \int_\mathit{\Gamma } {{C_{66}}} {\partial _z}{u_x}{\varphi _x}{n_x}{\rm{d}}\mathit{\Gamma } - \\ \int_\mathit{\Omega } {\left( {{C_{66}}{\partial _z}{u_x}{\partial _z}{\varphi _x} - \rho {p_{xz}}{\varphi _x}} \right)} {\rm{d}}\mathit{\Omega } - \int_\mathit{\Omega } \rho d_z^2{u_{x3}}{\varphi _x}{\rm{d}}\mathit{\Omega }\\ \int_\mathit{\Omega } \rho \left( {{\partial _t}{p_{xx}} + {d_x}{p_{xx}}} \right){\varphi _x}{\rm{d}}\mathit{\Omega } = - \int_\mathit{\Omega } {{C_{11}}} d_x^\prime {\partial _x}{u_x}{\varphi _x}{\rm{d}}\mathit{\Omega }\\ \int_\mathit{\Omega } \rho \left( {{\partial _t}{p_{xz}} + {d_z}{p_{xz}}} \right){\varphi _x}{\rm{d}}\mathit{\Omega } = - \int_\mathit{\Omega } {{C_{66}}} d_z^\prime {\partial _z}{u_x}{\varphi _x}{\rm{d}}\mathit{\Omega } \end{array} $ (20)

$ \begin{array}{l} \begin{array}{*{20}{c}} {\int_\mathit{\Omega } \rho \left( {{\partial _{tt}}{u_{z1}} + 2{d_x}{\partial _t}{u_{z1}}} \right){\varphi _z}{\rm{d}}\mathit{\Omega } = \int_\mathit{\Gamma } {{C_{66}}} {\partial _x}{u_z}{\varphi _x}{n_x}{\rm{d}}\mathit{\Gamma } - }\\ {\int_\mathit{\Omega } {\left( {{C_{66}}{\partial _x}{u_z}{\partial _x}{\varphi _z} - \rho {p_{zx}}{\varphi _z}} \right)} {\rm{d}}\mathit{\Omega } - \int_\mathit{\Omega } \rho d_x^2{u_{z1}}{\varphi _z}{\rm{d}}\mathit{\Omega }} \end{array}\\ \begin{array}{*{20}{l}} {\int_\mathit{\Omega } \rho \left[ {{\partial _{tt}}{u_{z2}} + \left( {{d_x} + {d_z}} \right){\partial _t}{u_{z2}}} \right]{\varphi _z}{\rm{d}}\mathit{\Omega } = \int_\mathit{\Gamma } {\left( {{C_{66}}{\partial _z}{u_x}{n_x}{n_x}} \right.} }\\ {{\varphi _z} + {C_{12}}{\partial _x}{u_x}{n_z}{\varphi _z}){\rm{d}}\mathit{\Gamma } - }\\ {\int_\mathit{\Omega } {\left( {{C_{66}}{\partial _z}{u_x}{\partial _x}{\varphi _z} + {C_{12}}{\partial _x}{u_x}{\partial _z}{\varphi _z}} \right)} {\rm{d}}\mathit{\Omega } - \int_\mathit{\Omega } \rho {d_x}}\\ {{d_z}{u_{z2}}{\varphi _z}{\rm{d}}\mathit{\Omega }} \end{array}\\ \begin{array}{*{20}{c}} {\int_\mathit{\Omega } \rho \left( {{\partial _{tt}}{u_{z3}} + 2{d_z}{\partial _t}{u_{z3}}} \right){\varphi _x}{\rm{d}}\mathit{\Omega } = \int_\mathit{\Gamma } {{C_{22}}} {\partial _z}{u_z}{\varphi _z}{n_z}{\rm{d}}\mathit{\Gamma } - }\\ {\int_\mathit{\Omega } {\left( {{C_{22}}{\partial _z}{u_z}{\partial _z}{\varphi _z} - \rho {p_{zz}}{\varphi _z}} \right)} {\rm{d}}\mathit{\Omega } - \int_\mathit{\Omega } \rho d_z^2{u_{z3}}{\varphi _z}{\rm{d}}\mathit{\Omega }} \end{array}\\ \int_\mathit{\Omega } \rho \left( {{\partial _t}{p_{zx}} + {d_x}{p_{zx}}} \right){\varphi _z}{\rm{d}}\mathit{\Omega } = - \int_\mathit{\Omega } {{C_{66}}} d_x^\prime {\partial _x}{u_z}{\varphi _z}{\rm{d}}\mathit{\Omega }\\ \int_\mathit{\Omega } \rho \left( {{\partial _t}{p_{zz}} + {d_z}{p_{zz}}} \right){\varphi _z}{\rm{d}}\mathit{\Omega } = - \int_\mathit{\Omega } {{C_{22}}} d_z^\prime {\partial _z}{u_z}{\varphi _z}{\rm{d}}\mathit{\Omega } \end{array} $ (21)

(20) 式和(21)式即为包含PML的弹性波方程等效积分弱形式。

2.3 空间离散及有限元控制方程

有限元的主要思想是将全区域Ω的求解转变为一系列子区域Ωe(或单元区域)上的局部求解, 将这些局部解组合即得到全局解。这里将全区域分解为若干线性三角形单元, 即每个单元内有3个节点, 进行双线性插值, 则单元内部的位移场可近似为:

$ {u^i} \approx \mathit{\boldsymbol{N}} \cdot \mathit{\boldsymbol{u}}_e^{\rm{T}} $ (22)

注意此处上标i表示位移xz分量, N=[N1   N2   N3]为形函数向量, uei=[ue1i    ue2i   ue3i]为局部解向量, 下标e表示第e个三角形单元, 则试探函数可表示为φi=δuei·NT。同理, 单元内部阻尼函数也可近似为:

$ \left\{\begin{array}{l}{d_{x} \approx \boldsymbol{N} \cdot\left(\boldsymbol{d}_{x}^{e}\right)^{\mathrm{T}}} \\ {d_{z} \approx \boldsymbol{N} \cdot\left(\boldsymbol{d}_{z}^{e}\right)^{\mathrm{T}}}\end{array}\right. $ (23)

其中

$ \left\{\begin{array}{lll}{\boldsymbol{d}_{x}^{e}} & {=\left[\begin{array}{ccc}{d_{x 1}^{e}} & {d_{x 2}^{e}} & {d_{x 3}^{e}}\end{array}\right]} \\ {\boldsymbol{d}_{z}^{e}} & {=\left[\begin{array}{lll}{d_{z 1}^{e}} & {d_{z 2}^{e}} & {d_{z 3}^{e}}\end{array}\right]}\end{array}\right. $ (24)

那么, 某一单元内dxdzdxdzdx2dz2可分别近似为:

$ \left\{\begin{array}{l}{d_{x}' \approx\left(\partial_{x} \boldsymbol{N}\right) \cdot\left(\boldsymbol{d}_{x}^{e}\right)^{\mathrm{T}}} \\ {d_{z}' \approx\left(\partial_{z} \boldsymbol{N}\right) \cdot\left(\boldsymbol{d}_{z}^{e}\right)^{\mathrm{T}}} \\ {d_{x}^{2} \approx \boldsymbol{N} \cdot\left(\boldsymbol{d}_{x}^{e} \otimes \boldsymbol{d}_{x}^{e \mathrm{T}}\right.} \\ {d_{z}^{2} \approx \boldsymbol{N} \cdot\left(\boldsymbol{d}_{z}^{e} \otimes \boldsymbol{d}_{z}^{e}\right)^{\mathrm{T}}} \\ {d_{x} d_{z} \approx \boldsymbol{N} \cdot\left(\boldsymbol{d}_{x}^{e} \otimes \boldsymbol{d}_{z}^{e}\right)^{\mathrm{T}}}\end{array}\right. $ (25)

其中

$ \left\{\begin{array}{lll}{\partial_{x} \boldsymbol{N}} & {=\left[\begin{array}{lll}{\partial_{x} N_{1}} & {\partial_{x} N_{2}} & {\partial_{x} N_{3}}\end{array}\right]} \\ {\partial_{z} \boldsymbol{N}} & {=\left[\begin{array}{lll}{\partial_{z} N_{1}} & {\partial_{z} N_{2}} & {\partial_{z} N_{3}}\end{array}\right]}\end{array}\right. $ (26)

符号“⊗”表示两个向量之间对应元素相乘。综上所述, 可将(20)式和(21)式写成有限元控制方程形式:

$ \left\{\begin{array}{l}{\boldsymbol{M} \partial_{t t} u_{x 1}+2 \boldsymbol{C}_{x} \partial_{t} u_{x 1}=\boldsymbol{M} p_{x x}-\boldsymbol{K}_{x 1} u_{x}-\boldsymbol{C}_{x x} u_{x 1}} \\ {\boldsymbol{M} \partial_{t t} u_{x 2}+\left(\boldsymbol{C}_{x}+\boldsymbol{C}_{z}\right) \partial_{t} u_{x 2}=-\boldsymbol{K}_{x 2} u_{z}-\boldsymbol{C}_{x z} u_{x 2}} \\ {\boldsymbol{M} \partial_{t t} u_{x 3}+2 \boldsymbol{C}_{z} \partial_{t} u_{x 3}=\boldsymbol{M} p_{x z}-\boldsymbol{K}_{x 3} u_{x}-\boldsymbol{C}_{z z} u_{x 3}} \\ {\boldsymbol{M} \partial_{t} p_{x x}+\boldsymbol{C}_{x} p_{x x}=-\boldsymbol{K}_{x 4} u_{x}} \\ {\boldsymbol{M} \partial_{t} p_{x z}+\boldsymbol{C}_{z} p_{x z}=-\boldsymbol{K}_{x 5} u_{x}}\end{array}\right. $ (27)

$ \left\{\begin{array}{l}{\boldsymbol{M} \partial_{t t} u_{z 1}+2 \boldsymbol{C}_{x} \partial_{t} u_{z 1}=\boldsymbol{M} p_{z x}-\boldsymbol{K}_{z 1} u_{z}-\boldsymbol{C}_{x x} u_{z 1}} \\ {\boldsymbol{M} \partial_{t t} u_{z 2}+\left(\boldsymbol{C}_{x}+\boldsymbol{C}_{z}\right) \partial_{t} u_{z 2}=-\boldsymbol{K}_{z 2} u_{x}-\boldsymbol{C}_{x z} u_{z 2}} \\ {\boldsymbol{M} \partial_{t t} u_{z 3}+2 \boldsymbol{C}_{z} \partial_{t} u_{z 3}=\boldsymbol{M} p_{z z}-\boldsymbol{K}_{z 3} u_{z}-\boldsymbol{C}_{z z} u_{z 3}} \\ {\boldsymbol{M} \partial_{t} p_{z x}+\boldsymbol{C}_{x} p_{z x}=-\boldsymbol{K}_{z 4} u_{z}} \\ {\boldsymbol{M} \partial_{t} p_{z z}+\boldsymbol{C}_{z} p_{z z}=-\boldsymbol{K}_{z 5} u_{z}}\end{array}\right. $ (28)

式中:M为全局质量矩阵; Kx1Kx2Kx3Kx4Kx5Kz1Kz2Kz3Kz4Kz5均为全局刚度矩阵; CxCzCxzCxxCzz为全局阻尼矩阵。

3 有限元控制方程的时间离散

观察(27)式和(28)式可见, 加载PML的弹性波有限元控制方程有如下一般形式:

$ \boldsymbol{M} \partial_{t t} \boldsymbol{u}+\boldsymbol{C} \partial_{t} \boldsymbol{u}+\boldsymbol{K} \boldsymbol{u}=\boldsymbol{f} $ (29)
$ \partial_{t} \boldsymbol{u}=f(\boldsymbol{u}) $ (30)

采用Euler方法进行时间离散求解类似(30)式的方程, 可得到四个中间变量pxxpxzpzxpzz的时间递推格式:

$ \left\{\begin{array}{l}{p_{x x}^{i+1}=\left(1-\Delta t \cdot \frac{\boldsymbol{C}_{x}}{\boldsymbol{M}}\right) p_{x x}^{i}-\Delta t \cdot \frac{\boldsymbol{K}_{x 4}}{\boldsymbol{M}} u_{x}^{i}} \\ {p_{x z}^{i+1}=\left(1-\Delta t \cdot \frac{\boldsymbol{C}_{z}}{\boldsymbol{M}}\right) p_{x z}^{i}-\Delta t \cdot \frac{\boldsymbol{K}_{x 5}}{\boldsymbol{M}} u_{x}^{i}} \\ {p_{z x}^{i+1}=\left(1-\Delta t \cdot \frac{\boldsymbol{C}_{x}}{\boldsymbol{M}}\right) p_{z x}^{i}-\Delta t \cdot \frac{\boldsymbol{K}_{z 4}}{\boldsymbol{M}} u_{z}^{i}} \\ {p_{z z}^{i+1}=\left(1-\Delta t \cdot \frac{\boldsymbol{C}_{z}}{\boldsymbol{M}}\right) p_{z z}^{i}-\Delta t \cdot \frac{\boldsymbol{K}_{z 5}}{\boldsymbol{M}} u_{z}^{i}}\end{array}\right. $ (31)

其中, Δt为时间步长。对于类似(29)式的方程, 可用中心差分法得到:

$ \left\{ \begin{array}{l} \left( {\frac{\mathit{\boldsymbol{M}}}{{\Delta {t^2}}} + \frac{{{\mathit{\boldsymbol{C}}_x}}}{{\Delta t}}} \right)u_{x1}^{i + 1} = \left( {\frac{{\mathit{\boldsymbol{2M}}}}{{\Delta {t^2}}} + {\mathit{\boldsymbol{C}}_{xx}}} \right)u_{x1}^i - {\mathit{\boldsymbol{K}}_{x1}}u_x^i + \\ \;\;\;\;\;\;\;\;\;\;\left( {\frac{{{\mathit{\boldsymbol{C}}_x}}}{{\Delta t}} - \frac{\mathit{\boldsymbol{M}}}{{\Delta {t^2}}}} \right)u_{x1}^{i - 1} + \mathit{\boldsymbol{M}}p_{xx}^i\\ \begin{array}{*{20}{c}} {\left( {\frac{\mathit{\boldsymbol{M}}}{{\Delta {t^2}}} + \frac{{{\mathit{\boldsymbol{C}}_x} + {\mathit{\boldsymbol{C}}_z}}}{{2\Delta t}}} \right)u_{x2}^{i + 1} = \left( {\frac{{2\mathit{\boldsymbol{M}}}}{{\Delta {t^2}}} - {\mathit{\boldsymbol{C}}_{xz}}} \right)u_{x2}^i + }\\ {\left( {\frac{{{\mathit{\boldsymbol{C}}_x} + {\mathit{\boldsymbol{C}}_z}}}{{2\Delta t}} - \frac{\mathit{\boldsymbol{M}}}{{\Delta {t^2}}}} \right)u_{x2}^{i - 1} - {\mathit{\boldsymbol{K}}_{x2}}u_z^i} \end{array}\\ \begin{array}{*{20}{c}} {\left( {\frac{\mathit{\boldsymbol{M}}}{{\Delta {t^2}}} + \frac{{{\mathit{\boldsymbol{C}}_z}}}{{\Delta t}}} \right)u_{x3}^{i + 1} = \left( {\frac{{2\mathit{\boldsymbol{M}}}}{{\Delta {t^2}}} - {\mathit{\boldsymbol{C}}_{zz}}} \right)u_{x3}^i - {\mathit{\boldsymbol{K}}_{x3}}u_x^i + }\\ {\left( {\frac{{{\mathit{\boldsymbol{C}}_z}}}{{\Delta t}} - \frac{\mathit{\boldsymbol{M}}}{{\Delta {t^2}}}} \right)u_{x3}^{i - 1} + \mathit{\boldsymbol{M}}p_{xz}^i} \end{array} \end{array} \right. $ (32)

$ \begin{array}{*{20}{c}} {\left( {\frac{M}{{\Delta {t^2}}} + \frac{{{\mathit{\boldsymbol{C}}_x}}}{{\Delta t}}} \right)u_{z1}^{i + 1} = \left( {\frac{{2\mathit{\boldsymbol{M}}}}{{\Delta {t^2}}} - {\mathit{\boldsymbol{C}}_{xx}}} \right)u_{z1}^i - {\mathit{\boldsymbol{K}}_{z1}}u_z^i + }\\ {\left( {\frac{{{\mathit{\boldsymbol{C}}_x}}}{{\Delta t}} - \frac{\mathit{\boldsymbol{M}}}{{\Delta {t^2}}}} \right)u_{z1}^{i - 1} + \mathit{\boldsymbol{M}}p_{zx}^i} \end{array} $ (33a)
$ \begin{array}{*{20}{c}} {\left( {\frac{\mathit{\boldsymbol{M}}}{{\Delta {t^2}}} + \frac{{{\mathit{\boldsymbol{C}}_x} + {\mathit{\boldsymbol{C}}_z}}}{{2\Delta t}}} \right)u_{z2}^{i + 1} = \left( {\frac{{2\mathit{\boldsymbol{M}}}}{{\Delta {t^2}}} - {\mathit{\boldsymbol{C}}_{xz}}} \right)u_{z2}^i + }\\ {\left( {\frac{{{\mathit{\boldsymbol{C}}_x} + {\mathit{\boldsymbol{C}}_z}}}{{2\Delta t}} - \frac{\mathit{\boldsymbol{M}}}{{\Delta {t^2}}}} \right)u_{z2}^{i - 1} - {\mathit{\boldsymbol{K}}_{zz}}u_x^i} \end{array} $ (33b)
$ \begin{array}{*{20}{c}} {\left( {\frac{\mathit{\boldsymbol{M}}}{{\Delta {t^2}}} + \frac{{{\mathit{\boldsymbol{C}}_z}}}{{\Delta t}}} \right)u_{z3}^{i + 1} = \left( {\frac{{2\mathit{\boldsymbol{M}}}}{{\Delta {t^2}}} - {\mathit{\boldsymbol{C}}_{zz}}} \right)u_{z3}^i - {\mathit{\boldsymbol{K}}_{z3}}u_z^i + }\\ {\left( {\frac{{{\mathit{\boldsymbol{C}}_z}}}{{\Delta t}} - \frac{\mathit{\boldsymbol{M}}}{{\Delta {t^2}}}} \right)u_{z3}^{i - 1} + \mathit{\boldsymbol{M}}p_{zz}^i} \end{array} $ (33c)

(31)、(32)和(33)式即为求取位移场各分量所采用的时间递推格式。第i+1时刻的位移场分量可表示为uxi+1=ux1i+1+ux2i+1+ux3i+1uzi+1=uz1i+1+uz2i+1+uz3i+1, 第i时刻速度场分量可由中心差分表示为:

$ \left\{ {\begin{array}{*{20}{l}} {v_x^i = \frac{{u_x^{i + 1} - u_x^{i - 1}}}{{2\Delta t}}}\\ {v_z^i = \frac{{u_z^{i + 1} - u_z^{i - 1}}}{{2\Delta t}}} \end{array}} \right. $ (34)

于是, 第i时刻总速度场为$v^{i}=\sqrt{\left(v_{x}^{i}\right)^{2}+\left(v_{z}^{i}\right)^{2}}$

4 模型算例 4.1 全空间均匀介质模型

为了说明理论反射系数的设定对PML吸收效果的影响, 进行了二维弹性波均匀全空间数值模拟。图 2a为均匀各向同性弹性介质模型, 密度为2 000 kg/m3, 纵波速度为1 800 m/s, 横波速度为1 100 m/s, 泊松比约为0.202。震源为正z方向激发的方向力源, 子波函数为雷克子波, 主频为5 Hz, 位于(2 000 m, 1 000 m)点处。采用对称结构化线性三角形网格(图 2b), 先将全区域用边长为10 m的正方形单元剖分, 然后将每个正方形单元划分成4个等腰直角三角形。在全空间条件下(即无面波情形), 最短波长为1 100/5/2.5=88.0 m, 三角形最大边长为10 m, 则一个波长内的采样点数至少为8.8。此外, 三角形单元的高最小为5 m, 线性三角形单元CFL(courant-friedrichs-lewy)数约为0.67, 时间步长应满足稳定性条件Δt≤0.67×5/1 800=1.86 ms, 本算例中取时间步长为1 ms。

图 2 均匀各向同性弹性介质模型(a)与对称结构化线性三角形网格(b)

根据前人研究[5], 阻尼函数dx可表示为:

$ d_{x}(\overline{x})=-(n+1) \cdot \lg R \cdot \frac{v_{\mathrm{P}}}{2 L} \cdot\left(\frac{\overline{x}}{L}\right)^{n} $ (35)

式中:n为常数, 控制阻尼系数变化的速率, 通常设置为2;L为PML吸收层的厚度。本算例中最大主波长为1 800/5=3 60 m, PML吸收层厚度L取为200 m, 对应于0.56倍的最大主波长(约10个网格点); x为PML区域内一点到PML内边界x坐标方向的垂直距离, 所以0≤ x/L≤1;R为理论反射系数, 一般在10-2~10-6之间取值。图 3展示了理论反射系数R取10-2、10-3、10-4、10-5对PML吸收效果的影响, 其中波场振幅强度均放大了10倍(下文作了相同处理)。从图 3可以看出, 理论反射系数为10-5时已观察不到明显的人工边界反射(图 3d), 说明在PML吸收层厚度一定时, 理论反射系数取值越小, PML吸收效果越好, 这与周凤玺等[20]对PML衰减因子参数优化的分析结果一致。为了更清晰地展示理论反射系数对边界反射的影响, 提取了(1 700 m, 1 300 m)点处的单道速度场记录, 如图 4所示。图 4中的参考解通过扩边模拟获得。可以看出, 随着理论反射系数的减小, 边界反射误差越来越小; 当R=10-5时, 加载PML模拟的速度场记录与参考解几乎一致, 说明PML吸收层厚度为半个最大主波长时, 理论反射系数应设置为R≤10-5

图 3 不同理论反射系数1.5 s时刻的总速度场波场快照对比 a理论反射系数R=10-2; b理论反射系数R=10-3; c理论反射系数R=10-4; d理论反射系数R=10-5
图 4 不同理论反射系数条件下(1 700 m, 1 300 m)点处地震记录 a速度场x分量; b速度场z分量
4.2 半空间均匀介质模型

狄利克雷边界条件不同于自由边界条件, 是一种刚性边界条件, 即波传播至此会被完全反射回去。假设PML吸收层最外层区域为Ωu, 那么狄利克雷边界条件可表示为:

$ {u^i} = 0\;在\;{\mathit{\Omega }_u}\quad i = x,z $ (36)

下面讨论半空间均匀介质情况下, 在PML吸收层最外层施加狄利克雷边界条件对PML吸收效果的影响。震源位置设为(2 000 m, 50 m), z=0处为自由表面。图 5图 6分别为未加载狄利克雷边界条件和不包括地表自由表面部分的PML最外吸收层上加载了狄利克雷边界条件的波场快照。从图 5a图 6a均可看到自由表面处有明显的面波产生。对比图 5b图 6b可见, PML最外吸收层上不加载狄利克雷边界条件会降低PML边界的稳定性, 产生可见的虚假波场(图 5b黄色虚线圈内)。因此, 在半空间均匀介质情形下, PML吸收层最外层加载狄利克雷边界条件能增强PML的数值稳定性。

图 5 未加狄利克雷边界条件的总速度场波场快照 a 2.0 s时刻; b 2.5 s时刻
图 6 加载了狄利克雷边界条件的总速度场波场快照 a 2.0 s时刻; b 2.5 s时刻

需要指出的是, 在PML吸收层最外层加载狄利克雷边界条件存在一个严重的数值陷阱。图 7为在包括地表自由表面部分的所有PML最外吸收层上加载了狄利克雷边界条件模拟的波场快照。可以看出, 在自由表面左右两端产生了非常严重的面波虚假反射, 严重干扰了目标区域内的波场求解, 说明处于地表自由表面的PML吸收层最外层部分不可加载狄利克雷边界条件。

图 7 在处于自由表面上的PML吸收层最外层上加载狄利克雷边界条件的总速度场波场快照 a 2.0 s时刻; b 2.5 s时刻
4.3 起伏地表非均匀介质模型

为了检验PML吸收层中包含狄利克雷边界在复杂介质条件下的吸收效果, 采用如图 8所示的起伏地表非均匀介质模型进行了有限元数值模拟测试。模型x方向3 000 m, 最大深度为1 000 m, 地表最大高程为210 m。上层介质密度为2 100 kg/m3, 纵波速度2 800 m/s, 横波速度1 500 m/s; 下层介质密度为2 250 kg/m3, 纵波速度3 400 m/s, 横波速度1 800 m/s。震源函数为雷克子波, 主频20 Hz, 位于(1 500 m, 100 m)点处。模型全区域被剖分为686 025个非结构化线性三角形单元, 其最大边长约为5.4 m, 高最小为3.1 m。由于稳定性要求, 时间采样间隔取为0.5 ms。PML吸收层厚度设置为170 m, 对应于一个最大主波长(约20个网格点)。理论反射系数取10-6, 以保证复杂模型中PML吸收层能充分吸收不同方向不同波长的能量。

图 8 起伏地表非均匀介质模型

图 9为0.4 s、1.0 s和2.0 s时刻的总速度场波场快照。其中, 图 9a图 9c图 9e为PML吸收层最外层不加狄利克雷边界时的模拟结果; 图 9b图 9d图 9f为PML吸收层最外层加载了狄利克雷边界条件时的模拟结果。从图 9可以看出, 在0.4 s和1.0 s时刻, PML吸收层最外层加载和不加载狄利克雷边界条件得到的波场信息一致; 但在2.0 s时刻, PML吸收层最外层不加狄利克雷边界条件时, 模拟结果出现了不稳定现象(见图 9e黄色虚线框部分), 这进一步说明了在PML吸收层最外层加载狄利克雷边界条件能增强PML的数值稳定性。

图 9 不同时刻总速度场波场快照(从上到下波场快照时刻依次为0.4 s、1.0 s和2.0 s) a, c, e PML吸收层最外层未加狄利克雷边界条件; b, d, f PML吸收层最外层加载狄利克雷边界条件
4.4 讨论

在全空间均匀介质算例中, PML吸收层厚度一定时, 理论反射系数越小, PML吸收效果越好。由于PML是一种衰减性质的边界条件, 这类加衰减层的边界条件经验上遵循衰减层越厚吸收效果越好的规律, 因此推断理论反射系数一定时, PML厚度越大, PML吸收效果越好。总之, 理论反射系数越小, PML吸收层越厚, PML吸收效果越好。针对不同复杂程度的速度模型, 理论反射系数和PML吸收层厚度的选择可根据数值试验确定。观察分析PML吸收层最外层加载狄利克雷边界条件对PML吸收效果的影响, 我们认为其与狄利克雷边界本身性质有关。狄利克雷边界条件是一种刚性边界条件, 在有限区域内向外传播的波场经过PML衰减层逐步衰减, 传播到PML最外层时遇到狄利克雷边界, 未被PML层完全吸收的波场会被完全反射回PML衰减层, 进行第二次指数衰减。因此, 理论分析认为, 在PML最外层添加狄利克雷边界条件可使传播到人工边界处的波场经历两倍厚度的PML吸收层的指数衰减, 减小了时间递推过程中计算误差的传递, 从而增强了PML的数值稳定性。但在自由表面PML吸收层区域加载狄利克雷边界条件时, 产生了严重的面波虚假反射, 这是因为自由表面的Rayleigh面波沿自由表面传播至PML吸收层时, 遇到狄利克雷边界产生了全反射, 继续传播至PML外边界自由表面角落时, 如同产生了二次震源, 最终造成强虚假反射。

5 结束语

本文推导和实现了PML吸收边界条件在二阶弹性波方程有限元数值解法中的变分形式, 通过全空间均匀介质数值实验, 测试了理论反射系数对PML吸收效果的影响, 发现当PML吸收层厚度约为半个最大主波长、理论反射系数设置小于(或等于)10-5时, PML吸收效果达到最优。

半空间均匀介质数值实验结果表明, 在PML吸收层最外层加载狄利克雷边界条件可提高PML的数值稳定性, 但是处于自由表面的PML吸收层最外层部分不可加载狄利克雷边界条件, 否则会产生严重的虚假反射。

起伏地表复杂介质数值模拟结果进一步表明, 在PML吸收层最外层加载狄利克雷边界条件对PML数值稳定性改善效果显著, 说明狄利克雷边界条件的引入可作为提高PML数值稳定性的有效手段之一。

参考文献
[1]
CLAYTON R, ENGQUIST B. Absorbing boundary conditions for acoustic and elastic wave equations[J]. Bulletin of the Seismological Society of America, 1977, 67(6): 1529-1540.
[2]
HIGDON R L. Absorbing boundary conditions for elastic waves[J]. Geophysics, 1991, 56(2): 231-241. DOI:10.1190/1.1443035
[3]
BERENGER J P. A perfectly matched layer for the absorption of electromagnetic waves[J]. Journal of Computational Physics, 1994, 114(2): 185-200. DOI:10.1006/jcph.1994.1159
[4]
CHEW W C, LIU Q. Perfectly matched layers for elastodynamics:A new absorbing boundary condition[J]. Journal of Computational Acoustics, 1996, 4(4): 341-359. DOI:10.1142/S0218396X96000118
[5]
KOMATITSCH D, MARTIN R. An unsplit convolutional perfectly matched layer improved at grazing incidence for the seismic wave equation[J]. Geophysics, 2007, 72(5): M155-M167. DOI:10.1190/1.2757586
[6]
CHEN H, ZHOU H, LIN H, et al. Application of perfectly matched layer for scalar arbitrarily wide-angle wave equations[J]. Geophysics, 2013, 78(1): T29-T39.
[7]
CHEN H, ZHOU H, LI Y. Application of unsplit convolutional perfectly matched layer for scalar arbitrarily wide-angle wave equation[J]. Geophysics, 2014, 79(6): T313-T321. DOI:10.1190/geo2014-0103.1
[8]
张衡, 刘洪, 李博, 等. VTI介质声波方程非分裂式PML吸收边界条件研究[J]. 石油物探, 2016, 55(6): 781-792.
ZHANG H, LIU H, LI B, et al. The research on unsplit PML absorbing boundary conditions of acoustic equation for VTI media[J]. Geophysical Prospecting for Petroleum, 2016, 55(6): 781-792. DOI:10.3969/j.issn.1000-1441.2016.06.002
[9]
柯璇, 石颖, 宋利伟, 等. 基于褶积完全匹配吸收边界的声波方程数值模拟[J]. 石油物探, 2017, 56(5): 637-643.
KE X, SHI Y, SONG L W, et al. Numerical modeling of acoustic wave equations based on convolutional perfectly matched layer absorbing boundary condition[J]. Geophysical Prospecting for Petroleum, 2017, 56(5): 637-643. DOI:10.3969/j.issn.1000-1441.2017.05.003
[10]
高正辉, 孙建国, 孙章庆, 等. 基于完全匹配层构建新方法的2.5维声波近似方程数值模拟[J]. 石油地球物理勘探, 2016, 51(6): 1128-1133.
GAO Z H, SUN J G, SUN Z Q, et al. 2.5D acoustic approximate equation numerical simulation with a new construction method of the perfectly matched layer[J]. Oil Geophysical Prospecting, 2016, 51(6): 1128-1133.
[11]
陈可洋. 声波完全匹配层吸收边界条件的改进算法[J]. 石油物探, 2009, 48(1): 76-79.
CHEN K Y. Improved algorithm for absorbing boundary condition of acoustic perfectly matched layer[J]. Geophysical Prospecting for Petroleum, 2009, 48(1): 76-79. DOI:10.3969/j.issn.1000-1441.2009.01.014
[12]
陈可洋. 完全匹配层吸收边界条件研究[J]. 石油物探, 2010, 49(5): 472-477.
CHEN K Y. Study on perfectly matched layer absorbing boundary condition[J]. Geophysical Prospecting for Petroleum, 2010, 49(5): 472-477. DOI:10.3969/j.issn.1000-1441.2010.05.006
[13]
高刚, 贺振华, 黄德济, 等. 完全匹配层人工边界条件中的衰减因子分析[J]. 石油物探, 2011, 50(5): 430-433.
GAO G, HE Z H, HUANG D J, et al. Analysis on attenuation factor in the processing of artificial boundary conditions of PML[J]. Geophysical Prospecting for Petroleum, 2011, 50(5): 430-433. DOI:10.3969/j.issn.1000-1441.2011.05.002
[14]
张衡, 刘洪, 李博, 等. TTI介质声波方程分裂式PML吸收边界条件研究[J]. 石油物探, 2017, 56(3): 349-361.
ZHANG H, LIU H, LI B, et al. The research on split PML absorbing boundary conditions of acoustic equation for TTI media[J]. Geophysical Prospecting for Petroleum, 2017, 56(3): 349-361. DOI:10.3969/j.issn.1000-1441.2017.03.005
[15]
马锐, 邹志辉, 芮拥军, 等. 基于SPML和海绵边界的伪谱法弹性波模拟复合吸收边界条件[J]. 石油物探, 2018, 57(1): 94-103.
MA R, ZOU Z H, RUI Y J, et al. A composite absorbing boundary based on the SPML and sponge absorbing boundary for pseudo-spectral elastic wave modeling[J]. Geophysical Prospecting for Petroleum, 2018, 57(1): 94-103. DOI:10.3969/j.issn.1000-1441.2018.01.013
[16]
KOMATITSCH D, TROMP J. A perfectly matched layer absorbing boundary condition for the second-order seismic wave equation[J]. Geophysical Journal International, 2003, 154(1): 146-153. DOI:10.1046/j.1365-246X.2003.01950.x
[17]
ZHAO J, SHI R. Perfectly matched layer-absorbing boundary condition for finite-element time-domain modeling of elastic wave equations[J]. Applied Geophysics, 2013, 10(3): 323-336. DOI:10.1007/s11770-013-0388-y
[18]
刘有山, 滕吉文, 刘少林, 等. 稀疏存储的显式有限元三角网格地震波数值模拟及其PML吸收边界条件[J]. 地球物理学报, 2013, 56(9): 3085-3099.
LIU Y S, TENG J W, LIU S L, et al. Explicit finite element method with triangle meshes stored by sparse format and its perfectly matched layers absorbing boundary condition[J]. Chinese Journal of Geophysics, 2013, 56(9): 3085-3099.
[19]
LIU Y S, TENG J W, LAN H Q, et al. A comparative study of finite element and spectral element methods in seismic wavefield modeling[J]. Geophysics, 2014, 79(2): T91-T104. DOI:10.1190/geo2013-0018.1
[20]
周凤玺, 张家齐, 张海威. 完全匹配层中衰减函数的参数优化分析[J]. 石油物探, 2016, 55(6): 793-799.
ZHOU F X, ZHANG J Q, ZHANG H W. The parameter optimization analysis for the attenuation function in the perfectly layer[J]. Geophysical Prospecting for Petroleum, 2016, 55(6): 793-799. DOI:10.3969/j.issn.1000-1441.2016.06.003