三角网格间断有限元法弹性波模拟精度分析

2021-08-18韩德超刘卫华司文朋

石油地球物理勘探 2021年4期

韩德超刘卫华司文朋*

(①中国石化石油物探技术研究院，江苏南京 211103；②中国石化地球物理重点实验室，江苏南京 211103)

0 引言

地震波数值模拟对认识地震波传播规律和分析干扰波形成机理具有重要意义，因此常被用于优化观测系统和指导处理解释[1-5]。此外，地震波数值模拟也是逆时偏移和波形反演的重要内核[6]。目前已有许多基于波动方程的数值模拟方法，如：有限差分法(FDM)[7-10]、伪谱法(PSM)[11-12]、有限体积法(FVM)[13-14]、无网格法[15]、有限元法(FEM)[16-18]，谱元法(SEM)[19-23]和间断Galerkin有限元方法(DGFEM)[24-27]。FDM因容易编程实现、易于并行等特点，目前应用最广泛。但是FDM对起伏地表的适应性有限，且准确模拟地表自由边界较为困难。PSM同样也存在边界条件难以准确模拟的问题，且三维计算需要全局通信，难于实现并行。FVM和FEM都可以使用非结构化的三角形(三维是四面体)网格以适应复杂的地表结构，且天然满足自由边界条件，但是FVM的高阶计算难于实现，较少应用于地震波模拟。FEM则需要存储大型稀疏矩阵并求逆，对计算资源要求极高。SEM和DGFEM是对FEM的改进。SEM使用GLL(Gauss Lobatto Legendre)插值点得到对角形式的质量矩阵，从而避免了矩阵求逆过程。但是经典的SEM使用的是四边形(三维是六面体)网格，对复杂地表的刻画缺乏灵活性。DGFEM在单元内部使用有限元方法处理，而在单元边界上使用FVM的数值流通量的方式。DGFEM集合了有限元方法的网格灵活性和高精度，同时又可以逐元求解，降低了对计算资源的消耗且便于并行计算。

数值稳定性、数值频散以及耗散问题是研究波动方程数值解的重要内容。有限元方法因使用的网格形态比较灵活，所以一般需构建周期性网格进行分析。四边形网格的模拟精度分析因为周期网格构建容易、基函数可以解耦等原因，已有大量研究[16，28-32]。而对于三角形网格的数值精度研究较少，刘少林等[33]和曹丹平等[34]分析了经典FEM在不同形态的周期性三角形网格中的声波方程频散和稳定性条件，也有少量的关于三角形网格的SEM频散分析研究[35-37]。对于四边形网格的DGFEM，De Basabe等[38]分析了不同基函数对数值精度的影响；贺茜君等[39]分析了Runge-Kutta(RK)时间格式的频散和稳定性；Meng等[40]分析了不同基函数对频散和稳定性的影响。对于三角形网格的DGFEM，Hu等[41]分析了半离散声波方程的频散和耗散，但未考虑时间格式；Käser等[42]从实验角度分析了弹性波任意高阶导数法(Arbitrary high-order derivatives，ADER)时间格式的数值收敛率，但未分析频散和耗散；Delcourte等[43]分析了蛙跳(Frog-leap，FL)时间格式的基于中心数值流的收敛率；He等[44]分析了声波方程RK和ADER时间格式的数值精度。

本文针对弹性波的数值模拟精度问题，通过构造任意等腰三角形单元组成的周期性网格，可以方便地分析各种形态的三角形单元对数值模拟精度的影响。从理论和数值计算上分析了基于局部Lax-Friedrichs数值流和3阶总变差不增的Runge-Kutta(简称TVD RK)时间格式的三角形网格DGFEM弹性波模拟的数值频散、耗散和稳定性，可为基于三角形网格的DGFEM弹性波模拟提供参数选取指导。

1 数值模拟方法

二维各向同性介质中的弹性波速度—应力方程为[42]

(1)

式中：λ、μ为拉梅常数；ρ是密度；σxx和σzz为正应力，σxz为剪切应力；vx和vz分别为质点速度的x和z分量；S1～S5表示外力震源。将式(1)改写成一阶线性双曲偏微分方程矩阵形式，有

(2)

式中：u=(σxx，σzz，σxz，vx，vz)T；P、B是式(1)中弹性参数的矩阵形式。

1.1 空间离散

DGFEM的空间离散需要将计算空间划分成不重叠的子区域，即单元。本文讨论二维三角形网格，直接给出式(2)的DGFEM空间离散格式[36]

(3)

(4)

1.2 时间离散

间断有限元求解有多种时间离散方法，如ADER、LF和RK法。ADER法将时间导数转化到空间导数的计算上，可以达到时间和空间任意精度，但是实现复杂[42]；LF时间格式通常和中心数值流一起使用[43，46]，且可以使用更大的时间步长从而提高计算效率，但是常用的LF时间格式精度只有二阶。本文使用RK法。

式(3)可简写为关于时间的常微分方程

(5)

式中忽略了震源项，L是空间离散的线性系统。常用的3阶TVD RK时间离散格式[25，39]为

(6)

2 数值稳定性及精度分析

2.1 基于三角形网格的理论分析方法

稳定性和数值精度的理论分析通常使用Von Neumann方法：假设研究区域为无界各向同性的均匀介质且无外力震源，分析平面波在研究区域的传播问题。

首先假设计算空间剖分为以顶角为θ、腰长为h的等腰三角形为基本单元的周期性网格(图1)，其中单元2-i是由单元1-i旋转所得，除此之外，单元完全相同。因此将单元1-1和单元2-1称为母单元，其他单元称为子单元。为了简化，将相邻的两类三角形组合在一起分析[41]，即图1中单元1-1和单元2-1作为整体进行分析。假设母单元中的边界顺序如图1b所示，子单元边界顺序同母单元。

图1 周期性三角形网格剖分示意图

考虑平面简谐波

Aexp(-iωt)exp(ikx)=A(t)exp(ikx)

(7)

则各子单元内的平面波与母单元平面波的相位关系如表1所示。假设波的传播方向与x轴夹角为γ，则kxh=2πδcosγ，kzh=2πδsinγ，其中k为波数矢量，δ=h/W是采样率，W是波长。

表1 各子单元与母单元的平面波相位关系

忽略震源项，单元1-2和单元2-1的间断Galerkin有限元法的半离散方程(式(3))可以写为

(8)

式中

(9)

(10)

其中的子矩阵Γ11、Γ12、Γ21、Γ22表达式见附录A。

2.2 稳定性分析

模拟时在网格尺寸h和介质速度cP已定的情况下，时间步长需满足Δt≤αmaxh/cP，其中α为Courant数，αmax=ΔtmaxcP/h是最大Courant数。DGFEM的3阶TVD RK方法可以写成[33]

(11)

(12)

式中ωh是与三角形腰长h有关的数值圆频率。可以看出，β=exp(-iωhΔt)是矩阵G的特征值。为了使模拟稳定，矩阵G的谱半径必须小于等于1，即|β|≤1[34]。需要注意的是，求解式(13)的特征值，系统会得到20N个特征值，要求每个特征值都满足|β|≤1，可用迭代方法求得最大Courant数αmax[40]。αmax越小说明模拟过程的稳定性越差，反之稳定性越强。

αmax随θ的变化曲线如图2a所示，可见，当单元过于狭长或扁平时模拟的稳定性较差。图2b为一、二、三阶单元αmax随r/h(r为内切圆半径)变化曲线，可知，各阶单元的αmax与r/h满足线性关系，其斜率分别为0.6137、0.3535、0.2683，约等于2/(2p+1)，其中p为单元阶次。由Δt=αh/cP可知，Δt与三角形内切圆半径r也满足线性关系，即Δt∝2/(2p+1)r/cP。

图2 不同阶单元αmax随θ(a)和r/h(b)的变化曲线

2.3 频散和耗散分析

(13)

将ωh代入式(12)左侧，可得

exp(ωiΔt)[cos(ωrΔt)-i sin(ωrΔt)]=βr+iβi

(14)

进一步可得

(15)

则数值频散系数d和耗散系数d′可分别表示为

(16)

(17)

当d大于1时，相速度大于介质速度，相反则说明相速度小于介质速度。d越接近1，模拟精度越高。由于P波和S波具有相同的频率、不同的空间采样率，因此分别对P波和S波进行分析。

当Courant数α=0.03、平面波传播方向角γ=0°时，计算θ分别为30°、60°、90°、120°时的二阶和三阶单元的P波和S波频散系数d随δ的变化曲线(图3)。由图3可以看出，二阶单元在δ小于0.4时，P、S波频散系数的相对误差小于0.9%，三阶单元的相对误差则小于0.05%，体现了DGFEM的弱频散特性。二阶单元一个波长内有2～3个网格即可满足频散要求。此外，由二阶单元的频散结果可以看出，在γ=0°方向上传播的平面波，θ为60°时频散最弱，其次是30°、90°、120°；S波的频散强弱相对关系与P波一致。不同形态(顶角θ不同)三阶单元的纵波频散相对关系与二阶单元一致。

图3 不同形态二阶(a)和三阶(b)单元的的P波和S波频散曲线

图4是一阶单元和二阶单元频散系数随传播方向的变化曲线，由图中可以看出，当单元为等边三角形时，频散系数随传播方向变化的周期为60°，且沿单元边界方向传播时频散最小，而垂直网格线传播时频散最大。而经典有限元方法在等边三角形网格中基本没有数值各向异性。二者的区别可能是由DGFEM的数值流形式引起的。其他形态三角形单元的频散系数随传播方向的变化以180°为周期。等边三角形单元的频散曲线更接近圆形，数值各向异性最弱，而其他形态的三角形单元的频散多呈条带状，数值各向异性更强。在图4b中，当θ=30°时频散曲线在γ=15°和195°处频散最大，即垂直于等腰三角形底边的方向，最小频散的方向为沿底边方向；当θ=90°时，最大频散的方向为γ=135°和225°，当θ=120°时，最大频散的方向为γ=150°和330°，即沿底边的传播方向频散最大，垂直于底边方向频散最小。

“人设”先行，架空了历史，扭曲了情节。热播的《延禧攻略》中，主角魏璎珞的“人设”可谓鲜明至极：她聪明果敢、为人刚强，见招拆招、快意恩仇，以披靡之势走向人生巅峰。“底层逆袭”“人人都爱我”等现实中无法实现的白日梦在其中得到代偿，使得这样的“人设”很是讨喜。网剧情节有所夸张无可厚非，但如果为了凸显“人设”而不顾历史常识，违背起码的生活逻辑，生造出在现实中不可能发生的情节，“人设”就失去了坚实的支撑与丰满的内里，虽满足了受众一时的感官刺激，却成了虚假而空洞的存在。

图4 不同形态单元纵波频散系数随传播方向角γ的变化曲线

θ分别为30°、60°、90°、120°时，计算二阶和三阶单元的网格P波和S波的耗散系数随采样率的变化曲线(图5)，其中α=0.03。为了明显区分曲线，图5是模拟了200个时间步长后的耗散曲线(即exp(200ωiΔt))。由图可见，随着δ的增加，数值耗散越强；对于相同的δ，S波的耗散略小于P波，通常情况下，同一介质内的S波采样率大于P波的采样率；另外，随着单元阶数的增大，耗散的影响减弱。通常地震波需要模拟几千甚至几万个时间步长，数值耗散将对模拟结果产生较大影响。

图5 不同形态二阶(a)和三阶(b)单元的的P波和S波耗散曲线

图6为θ不同时，一阶单元网格的exp(ωiΔt)随传播方向的变化曲线，其中α=0.03、δP=0.11、δS=0.20。由图可见，耗散曲线的对称性与频散曲线(图4)相同。单元顶角θ=60°时，耗散系数整体上更接近1。

图6 不同形态一阶单元的耗散系数随传播方向角γ的变化曲线

表2梳理了上述周期性网格的理论分析结果。表中数值各向异性是基于已比较过的单元形态。

表2 不同形态三角形单元的最大αmax及数值各向异性特征

3 数值实验

首先通过模拟二维均匀各向同性介质平面波单频成分分析不同形态单元的误差。不考虑震源项，均匀各向同性介质中弹性平面波应力和速度分量的解析解表达式分别为

(18)

(19)

上式每个分量都包含一个沿着(nx，nz)T方向传播的平面P波和S波。取t=0时刻的解析值作为数值模拟的初始条件。

分别用θ为30°、60°、90°的周期性网格进行均匀模型模拟，并对每个θ取其数值各向异性的两个对称轴的方向作为平面波传播方向，即θ=30°时，取γ=15°或105°；θ=60°时，取γ=0°或30°；θ=90°时，取γ=45°或135°。等腰三角形单元的腰长h分别设置为10.0、12.5、20.0、25.0m。均匀模型尺寸为5000m×5000m，cP=4000m/s、cS=2000m/s、ρ=2000kg/m3；模拟频率f0=40Hz，Δt=2ms，模拟时长tmax=250ms。

图7为θ=90°、γ=45°时，tmax时刻的模拟波场与解析解的对比。由于没有使用吸收和周期性边界，所以取模型中没有受到边界反射干扰的区域。由图7可以看出，使用一阶单元且h=25.0m时的模拟波场与解析解在能量上有很大差异，模拟误差主要来源于是数值耗散(纵波空间采样率为0.25，横波采样率为0.5)。使用二阶单元，h=10.0m时的模拟波场与解析解接近(纵、横波的空间采样率分别为0.1和0.2)。图8为图7深度为2000m处的波场，可见：二阶单元h=10.0m和h=12.5m两种网格的计算结果与解析解对应较好；当使用一阶单元时，即使空间采样率取0.1，也会产生很强的数值耗散。

图7 θ=90°、γ=45°时tmax时刻的不同网格模拟的40Hz单频波场与解析解的对比

图8 θ=90°、γ=45°时tmax时刻、深度2000m处不同网格模拟的波场与解析解的对比

(20)

式中Ω为选取的计算区域。具体地，当θ=30°时，Ω为4600m≤x≤4900m，1200m≤z≤1400m；当θ=60°时，Ω为3500m≤x≤4500m，2000m≤z≤3000m；当θ=90°时，Ω为2000m≤x≤3000m，2000m≤z≤3000m。

图9、表3给出了不同θ的一阶和二阶单元网格在不同传播方向上的模拟误差。可以看出，由于误差主要受数值耗散的影响，对于相同θ，两个对称轴方向上的误差相对关系与耗散随传播方向的变化曲线(图5和图6)一致；对于相同阶次的单元，其网格尺寸与误差在双对数坐标系下呈线性关系：一阶单元时，斜率约为2；二阶单元时，斜率约为4；单元从一阶提高到二阶，网格越小，误差减小越明显。

表3 不同形态的一阶和二阶单元在不同传播方向上的误差统计

同时，计算了在一般性非规则网格中的3阶TVD RK时间格式的DGFEM的误差(图9、表4)。图10是误差统计区域内的非结构性网格，呈无规律排列。图9中非规则网格误差曲线同样表现为线性，而且在45°和135°两个方向上没有出现明显的数值各向异性。

表4 非结构网格中一阶和二阶单元在45°和135° 方向上的误差

图9 不同θ的时一阶和二阶单元网格在不同传播方向上的模拟误差随网格尺寸的变化曲线蓝色、红色和紫色线分别代表顶角为60°、30°、90°的周期性网格结果，绿色为非规则网格结果

图10 模拟使用的非结构性网格

应用上述均匀各向同性介质模型进行点源激发模拟，以更加直观地说明网格对DGFEM的波场特征的影响。震源为主频40Hz的Ricker子波，时间步长Δt=8ms。图11为h取不同值时一阶直角三角形周期性网格与一般非结构化网格模拟得到的0.3s时刻的vz波场快照对比。图12为二阶直角三角形周期性网格与一般非结构化网格模拟得到的0.3s时刻的vz波场快照对比。由图11可以看出，当网格阶数为一时，随着网格尺寸的增加，两种网格的波场能量总体都有所下降。由h=20.0m的直角三角形周期性网格模拟的波场快照(图11a右和图12a右)可以看出，在135°和315°方向，纵波和横波的波形变宽，旁瓣更加明显，能量比45°方向弱很多，这是数值耗散所致。而在图11b右和图12b右所示的h=20.0m的一般非结构化网格模拟的性网格的波场快照上，观察不到明显的方向差异，说明非结构数值各向异性不明显。图12b右中有部分杂乱的波场，这是由于网格尺寸增大导致的点加载的震源噪声。通常可以将震源以高斯分布的方式加载在震源点附近的一定范围内以减弱这一噪声[27]。

图11 h取不同值时一阶直角三角形周期性网格(a)与一般非结构化网格(b)模拟的0.3s时刻的vz波场快照对比

图12 h取不同值时二阶直角三角形周期性网格(a)与一般非结构化网格(b)模拟的0.3s时刻的vz波场快照对比

4 结论

本文通过构造周期性三角形网格模型，分析了基于Runge-Kutta的间断Galerkin有限元方法的稳定性条件、数值频散以及数值耗散，得到以下结论和认识。

(1)DGFEM稳定模拟的最大时间步长Δtmax与单元内切圆半径存在线性关系。分析结果表明，等边三角形的稳定性条件最宽松。

(2)从理论和实验角度说明了基于局部Lax-Friedrichs数值流的DGFEM的弱频散和强耗散特性，且在不同传播方向上，二者都表现出方向各向异性。数值各向异性呈双轴对称性。等边三角形网格的数值各向异性最弱，而直角三角形和钝角三角形网格表现出的数值各向异性较强。一般的非结构性网格由于网格排列没有一定规律，因此没有表现出明显的数值各向异性，说明在网格生成时，应尽量避免钝角三角形和呈大面积排列的直角三角形，以减小数值各向异性的影响。

(3)误差分析显示，模拟误差与网格尺寸在双对数坐标系下呈线性关系。数值实验表明，模拟误差主要来源是基于局部Lax-Friedrichs数值流的数值耗散。实验中，使用一阶单元模拟的波场均有较强耗散。而选取二阶单元，一个横波波长内包含四个单元则耗散明显减弱。

关于其他数值流形式的频散和耗散特性，以及其对模拟精度的影响仍需进一步研究。