完全匹配层在时域有限元弹性波数值模拟中的应用

2019-08-12覃发兵高志伟解皓楠徐振旺

石油物探 2019年4期

覃发兵,高志伟,解皓楠,徐振旺

(1.长江大学地球物理与石油资源学院,湖北武汉430100;2.长江大学管理学院,湖北荆州434023;3.中国石油天然气股份有限公司新疆油田分公司采油一厂,新疆克拉玛依834000;4.中国石油天然气股份有限公司辽河油田分公司勘探开发研究院,辽宁盘锦124010)

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

根据地震波场数值模拟时是否在求解区域最外层附加额外的几何空间,可将近四十年发展起来的人工边界处理技术分为两类:吸收边界方法和吸收衰减层方法。吸收边界方法是在计算区域最外层网格节点上引入额外变量压制虚假反射。由于没有增加计算区域的几何空间,吸收边界方法产生的额外计算量较少。经典的吸收边界方法为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的基本原理

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

(1)

图1 目标区域和PML边界层区域图示

(2)

依据链式法有:

(3)

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

(4)

其中

(5)

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

(6)

2 含PML的有限元控制方程

2.1 含PML的弹性波方程

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

(7)

(8)

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

(9)

(10)

又知:

(11)

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

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

(13)

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

(14)

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

(15)

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

(16)

和

(17)

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

(18)

和

(19)

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

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

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

(20)

和

φz+C12∂xuxnzφz)dΓ-

dzuz2φzdΩ

(21)

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

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

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

(22)

(23)

其中

(24)

(25)

其中

(26)

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

(27)

和

(28)

式中:M为全局质量矩阵;Kx1、Kx2、Kx3、Kx4、Kx5、Kz1、Kz2、Kz3、Kz4和Kz5均为全局刚度矩阵;Cx、Cz、Cxz、Cxx和Czz为全局阻尼矩阵。

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

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

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

(31)

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

(32)

和

(33a)

(33b)

(33c)

(34)

4 模型算例

4.1 全空间均匀介质模型

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

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

(35)

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

4.2 半空间均匀介质模型

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

(36)

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

图3 不同理论反射系数1.5s时刻的总速度场波场快照对比a 理论反射系数R=10-2; b 理论反射系数R=10-3; c 理论反射系数R=10-4; d 理论反射系数R=10-5

图4 不同理论反射系数条件下(1700m,1300m)点处地震记录a 速度场x分量; b 速度场z分量

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

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

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

图7 在处于自由表面上的PML吸收层最外层上加载狄利克雷边界条件的总速度场波场快照a 2.0s时刻; b 2.5s时刻

4.3 起伏地表非均匀介质模型

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

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

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

图9 不同时刻总速度场波场快照(从上到下波场快照时刻依次为0.4s、1.0s和2.0s)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外边界自由表面角落时,如同产生了二次震源,最终造成强虚假反射。