波动数值模拟中的外推型人工边界条件1)

邢浩洁 李小军 刘爱文 李鸿晶 周正华 陈 苏

∗(中国地震局地球物理研究所,北京 100081)

†(北京工业大学建筑工程学院,北京 100124)

∗∗(南京工业大学土木工程学院,南京 211816)

††(南京工业大学交通运输工程学院,南京 211816)

引言

波动数值模拟是当前地球物理、地震学、地震工程、海洋声学、电磁波等多个学科领域广泛应用和研究的前沿热点技术.人工边界条件为该技术的基础问题之一,它表示有限计算模型边界上的辐射条件,其精度和稳定性对于模拟结果能否很好地反映实际无限介质中的波动过程至关重要.

自20 世纪60 年代以来,已有大量人工边界条件(artificial boundary condition,ABC)被学者和工程师们所提出,文献[1-2] 对它们作了较为系统的介绍.对于数量众多的ABC,如何有效掌握并在不同问题中进行选用已成为一个突出的实际问题.一直以来关于各个特定ABC 的应用及探讨较多,而从整体上针对ABC 理论与公式体系的系统性研究[3]则相对较少.这不仅导致较为复杂的新近研究成果不易被理解及吸收,甚至对于一些简单且有着广泛应用的经典人工边界条件,也仍然存在认识方面的不足.

廖氏透射边界[4-6]是中国学者在人工边界领域做出的重要成果.它采用等距离散点的外推公式(Newton-Gregory 公式) 来推算任意人工边界节点在各个时刻的运动.廖振鹏等将该边界条件与“误差波多次透射”的物理机制联系起来,据此命名为“多次透射公式”(multi-transmitting formula,MTF),也可称为廖氏透射边界.MTF 以简洁的离散表达式来描述任意外行波的一般传播过程,完美地将严格的数学公式和清楚的物理机制融合在一起,其基本思想和表达形式均具有普适性.在廖振鹏等的工作基础上,作者经过理论分析和波动模拟实践发现,MTF 边界的上述特征实际上奠定了一大类人工边界条件的思想和公式基础.本文将这类边界称为外推型人工边界条件,其特点是人工边界节点在每一时刻的运动,直接由一组临近离散网格点在前若干时刻的运动进行时空外推得到.具有时空外推特点的人工边界条件除了MTF 之外,主要包括经典的旁轴近似边界[7-8]、Higdon 边界[9-10]、Givoli-Neta 等[11-12]的辅助变量高阶边界、Hagstrom 等[13-14]的辅助变量高阶边界、Guddati 等[15-16]的任意大角度波动方程(arbitrarily wide-angle wave equations,AWWE)、Lindman 数值优化边界[17]、Peng-Toks¨oz 数值优化边界[18],以及Randall 势分解方法[19]、Liu-Sen 混合边界方法[20-21]等.这些边界使用相同的控制参数(边界阶次和一组计算波速),通过某种等价的中间形式相联系,并且在数值模拟中表现出相似的精度和稳定性.

现有研究虽然早已指出 MTF 边界和经典的Clayton-Engquist 旁轴近似边界以及Higdon 边界存在密切联系[9,22],但是在讨论和应用中仍然将它们视作不同的人工边界条件.由于上述边界至今尚未被归并到一个理论框架当中,导致对于它们的理论联系以及与精度和稳定性相关的诸多应用问题还缺乏系统性认识,如:表达形式与波动类型、内域离散格式、几何构型的无关性,边界阶次和计算波速参数对边界精度的控制原理,独立的边界离散格式所导致的稳定性问题等.这不利于该类人工边界条件的研究成果的相互借鉴和进一步发展与应用.另外值得指出的是,真正与上述边界在表达形式、数值离散格式、精度控制原理以及稳定性各个方面都存在显著不同的是另外两大类人工边界条件,分别为以黏性边界、黏弹性边界等为代表的应力型边界[2,23-24]和以函数衰减层、完美匹配层等为代表的衰减层型边界[1,25-26].

本文工作以新提出的两个基本边界公式作为出发点,通过公式推导证明它们的等价性并论证其与具有时空外推特点的上述各种ABC 的理论联系,建立外推型人工边界条件理论.深入研究外推型ABC的精度控制原理,分析并强调可调的计算波速参数对于处理复杂波动问题的价值.最后给出基本边界公式的数值离散格式,并通过算例验证外推型ABC理论的合理性.

1 外推型人工边界条件的基本公式

1.1 两个基本边界公式

本文提出两个基本公式作为应用和讨论外推型人工边界条件的出发点,即

这两个公式本身是最简单有效的外推型ABC,它们奠定了以时间-空间外推(简称时空外推)方式计算任意时刻单个人工边界节点运动的理论和方法基础.其他外推型ABC 大多可以从这两个基本边界公式转化得到,或者通过某种等价的中间形式与之相关联.式(1)为MTF 引入多个人工波速caj(j=1,2,···,N,这里N为边界阶次)得到的优化形式,本文记为优化的MTF 或caj-MTF 边界.它为离散表达式,人工边界节点的波场值由一系列参考点的波场值外推得到.式(2) 为Higdon 边界在一个统一局部坐标s0t上定义并使用人工波速caj作为边界参数的一般形式,记为统一的Higdon 或caj-Higdon 边界.它为连续表达式,是由多个一阶单向波动算子相乘形成的波动微分方程.

式(1) 和式(2) 定义在如图1 所示的以各个待求解的人工边界节点为原点的统一局部坐标s0t中,坐标空间轴s位于经过人工边界节点的一条指向内域的离散网格线上(正方向向内),时间轴t与整体模型一致.式中u=u(s,t)为局部坐标s0t中的波场函数;N为边界阶次,在式(1) 中对应于离散时空外推算子(I-B(caj))的个数,在式(2)中为集成的一阶单向波动微分算子(∂/∂t-caj∂/∂s)的个数;caj(j=1,2,···,N) 为一组计算波速参数称之为人工波速,其取值具有任意性,但较优的选择通常为等于或稍大于介质物理波速的值(后文将具体探讨);Δt为时间步长,取值同样具有任意性,不过为便于计算,通常取为与内域格式一致;I为单位算子,在乘积中可以略去.在局部坐标s0t中,离散算子B(caj) 定义为:B(caj)u(s,t)=u(s+cajΔt,t-Δt),即空间上向计算区域内部、时间上向以前时刻移动(cajΔt,-Δt),它用于定出caj-MTF 边界的各个参考点的位置(相对于人工边界节点).离散算子的求和、乘积运算法则为:[B(ca1)+B(ca2)]u(s,t)=u(s+ca1Δt,t-Δt)+u(s+ca2Δt,t-Δt),B(ca1)B(ca2)u(s,t)=u(s+ca1Δt+ca2Δt,t-2Δt),依次类推.边界表达式(1)和式(2)可直接用于处理单分量波动问题,如声波、SH 波动,也可以分别应用于多分量波动问题的各个分量,如弹性波、电磁波、固液两项介质波动等.

图1 用于定义基本边界公式(1)和(2)的统一局部坐标示意图Fig.1 Sketch map of the unified local coordinate used to define the basic boundary formulas(1)and(2)

图2 给出了二阶、三阶廖氏透射边界(MTF)和本文外推型边界公式(1) (caj-MTF) 所涉及的离散参考点示意图.基于单一人工波速的廖氏透射边界是边界(1) 的一个特例.利用上述离散算子廖氏透射边界可以表示为(I-B(ca))Nu=0.如图2(a)和图2(c) 所示,MTF 的参考点为等距离分布且在每个t-jΔt时刻只有一个点,各个参考点的算子表示依次为:u(s+caΔt,t-Δt)=B(ca)u(s,t),u(s+2caΔt,t-2Δt)=B2(ca)u(s,t),u(s+3caΔt,t-3Δt)=B3(ca)u(s,t),其余类推.二阶MTF 采用前两个参考点,其系数分别为2 和-1;三阶MTF 采用前3 个参考点,其系数分别为3,-3 和1.而从图2(b) 可以看到,二阶caj-MTF 与二阶MTF 相比,2B(ca)u(s,t)被优化成[B(ca1) +B(ca2)]u(s,t),-B2(ca)u(s,t) 被优化成-B(ca1)B(ca2)u(s,t).类似地,根据图2(d),三阶caj-MTF 与三阶MTF 相比,3B(ca)u(s,t) 被优化成[B(ca1)+B(ca2)+B(ca3)]u(s,t),-3B2(ca)u(s,t) 被优化成-[B(ca1)B(ca2)+B(ca2)B(ca3)+B(ca3)B(ca1)]u(s,t),B3(ca)u(s,t) 被优化成B(ca1)B(ca2)B(ca3)u(s,t).因此,caj-MTF 边界在各个时刻可能有一个或多个参考点,它对传统MTF 的优化是通过在各个离散时刻上由多个不同参考点的波场值代替原来单个参考点的波场值而实现的.如果从MTF“误差波多次透射”的物理机制角度来分析caj-MTF,那么后者可以被看作是保留了“多次透射”过程,但优化了对各阶“误差波”的描述.

图2 二阶、三阶MTF 和caj-MTF 边界的参考点示意图Fig.2 Sketch map of the reference points in 2nd-and 3rd-order MTF and caj-MTF boundaries

Higdon 边界对于声波和弹性波分别具有不同形式,它们基于整体坐标x,y或z给出,且根据外行波沿坐标的正、负方向传播而有所不同.式(2)可以视为这些现有Higdon 边界公式的统一形式.理由如下:考虑到边界参数取值的任意性,现有Higdon 边界的声波和弹性波形式实际上是一致的[9-10];用统一局部坐标的空间轴s代替整体坐标x,y或z,在推算人工边界节点运动时结果并无差异;同样,统一考虑s轴指向内域,与考虑s轴指向外域亦具有相同效果,Higdon[27]对此已有结论.用统一局部坐标的空间轴s代替整体坐标x,y或z还有一个优势,就是s轴可以不与整体坐标轴平行.式(2) 使得现有Higdon 边界针对不同方向边界分别给出表达式再进行数值离散的做法得到大大简化,它略去了不必要的冗余,回归到时空外推的本质.

式(1) 给出的多人工波速离散边界表达式无法单独从传统廖氏透射边界理论导出,其借鉴了Higdon 边界的计算波速设置.式(2)给出的简洁而统一的连续边界表达式同样无法单独从现有Higdon 边界导出,其得益于廖氏透射边界思想和参数的引入.Liao[28]曾指出过廖氏透射边界与Higdon 边界的这种互补性.

1.2 添加小量修正的边界公式

在基本边界公式(1)和式(2)中添加一个小量对其进行修正所得到的形式,在消除高阶外推型ABC的计算漂移失稳[29]和处理强衰减外行波方面具有重要价值.添加小量修正的基本边界公式为

式(3) 中的小量γ 为一个小正数,如0.001～0.05 之间.式(4)中的小量εj(j=1,2,...,N)为非负数,其中至少一个不为零.

周正华和廖振鹏[30]指出添加小量在物理上考虑了介质的几何扩散特性或引入了介质的阻尼特性.可认为这两种物理特性对应于同一种数学机制,即考虑了波场的衰减性.由此不难得出,陈少林和廖振鹏[31]所提出的衰减波场MTF、以及Bayliss 和Turkel[32]针对柱面波或球面波扩散问题(即外行波含有显著的几何衰减)而提出的边界表达式,分别与式(3)、式(4)在形式上具有一致性.

2 caj-MTF 与caj-Higdon 边界的等价性

本节证明外推型ABC 理论与公式体系中的一个关键结论:离散的caj-MTF 边界与连续的caj-Higdon边界是等价的.有两种证明方法:一是利用二元函数泰勒展开法则对caj-MTF 中各个参考点的波场函数在点(s,t) 处进行展开并保留至所需精度阶,化简可得caj-Higdon;二是将caj-MTF 中的各个外推算子改写成某种与之等效的差分算子,使之形成类似于caj-Higdon 的结构,对差分算子取极限即得到后者.

证明 1:对于一阶caj-MTF,其具体形式为u(s,t)=u(s+ca1Δt,t-Δt).将等号右端项进行二元泰勒展开并保留至一阶项,得到

化简式(5)即得一阶caj-Higdon 表达式.对于二阶caj-MTF,先写出其具体形式,再将等号右端各项进行二元泰勒展开并保留至二阶项,得到

化简式(6)就能够得到二阶caj-Higdon 表达式.对于更高阶次边界,可用归纳法依次类推.

证明2:考虑caj-MTF 的单个算子形式(I-B(caj))u=0.分别定义空间移动算子K和时间移动算子Z-1:Ku(s,t)=u(s+cajΔt,t),Z-1u(s,t)=u(s,t-Δt),则式(1)的单个算子形式可以表示为(I-Z-1K)u=0.该式可改写为

将式(7) 中第一个Δs用cajΔt替换然后对每一项乘以caj,得到

对式(8) 中的差分算子取极限,即可得到(∂/∂t-∂/∂s)u=0.这表明caj-MTF 中每个离散形式的时空外推算子的极限形式为caj-Higdon 的各个单向波动微分算子,所以这两个边界表达式等价.

由于证明1 和证明2 都建立在基本的数学原理之上,证明过程中并未涉及任何波动方程,所以caj-MTF 与caj-Higdon 边界的等价性是普遍存在的.这远远强于已有研究所指出的Higdon 边界或MTF 边界与Clayton-Engquist 声波边界之间的相似性[9,22].这种等价性将离散caj-MTF 边界与连续caj-Higdon 边界紧密联系在一起,在讨论和应用中二者通常可以互相替换.

3 外推型ABC 的理论联系

这一部分将通过公式推导来论证其他外推型ABC 大多可以从基本边界公式(1) 和(2) 转化得到,或通过某种等价的中间形式与之相关联.

3.1 Reynolds 边界、Keys 边界

Reynolds[33]和Keys[34]各自独立提出了一种ABC.以外法向为x轴负向的边界为例,其表达式分别为

将式(9)、式(10) 与式(2) 相比较可知,它们都是式(2)的特例.因此,caj-Higdon 边界已将这两种边界概括在内.

3.2 旁轴近似边界

旁轴近似边界主要为经典的 Clayton-Engquist声波和弹性波边界(CE 边界)[7]、Halpern-Trefethen 大角度单向波动方程[8]、以及Fuyuki-Matsumoto(FM)[35]、Emerman-Stephen (ES)[36]、Stacey[37]对CE弹性波边界的改进.

对于声波旁轴近似边界,Higdon[9]已证明它们可以由Higdon 边界结合声波方程来导出.以外法向为z轴负向的边界为例,考虑二阶情形,CE 和Halpern-Trefethen 所提出的各种二阶声波边界表达形式为

其中c为声波波速,p0和p2为不同边界的控制参数,具体见文献[8].容易验证,在式(2)的二阶形式中引入声波方程c2(∂2u/∂x2+∂2u/∂z2)=∂2u/∂t2,将关于z的二阶偏导用关于x和t的二阶偏导所替换,就导出了式(11).式(2) 中参数ca1,ca2与式(11) 中参数p0,p2存在如下对应关系:p0=(1+α1α2)/(α1+α2),p2=α1α2/(α1+α2),其中α1=ca1/c,α2=ca2/c.

对于弹性波旁轴近似边界,作者分析发现它们缺乏严格的理论基础,难以达到较好的精度.由于弹性波方程不能由其频散关系唯一地确定,旁轴近似技术实际上无法用于推导弹性波边界.Clayton 和Engquist[7]仅通过粗略地模仿声波边界给出了如下弹性波边界(z轴负向)

然而,式(12)既无法从旁轴近似角度进行解释,也无法由caj-Higdon 边界结合弹性波方程导出.Fuyuki 和Matsumoto[35]、Emerman 和Stephen[36]、Stacey[37]在式(12) 基础上所做的改进并未解决这一问题.数值试验表明这些边界精度较差.为了从理论上揭示所谓的旁轴近似弹性波边界的不合理性,在波数平面上绘出了CE,FM,ES 和Stacey 弹性波边界的频散曲线,如图3 所示.作为对比,图中还给出了声波边界(11)的频散曲线,A1,A2,A3 分别表示一阶、二阶、三阶声波边界.

图3 中各曲线与上半圆的接近程度反映了人工边界条件的精度.各个弹性波边界频散曲线与上半圆的吻合度很差,这证实了旁轴近似弹性波边界的不合理性.与之相比,旁轴近似声波边界精度较好,且随着阶次上升迅速提高,与理论结果相符.

图3 弹性波、声波旁轴近似边界的频散曲线Fig.3 Dispersion curves of acoustic-and elastic-wave paraxial-approximation boundaries

3.3 辅助变量高阶边界

针对高阶Higdon 边界数值离散比较困难的问题,Givoli-Neta 等[11-12]、Hagstrom-Warburton 等[13-14]分别提出一种可顺利离散至任意阶次的辅助变量高阶边界(G-N 边界和H-W 边界).这两种边界都是先给出一个与Higdon 边界等价的辅助变量递推关系,再引入声波方程将递推关系中的边界法向导数转换为切向导数,得到只含有边界切向导数和时间导数的辅助变量递推关系,形成实用的高阶边界.由于声波方程的引入,这两种为声波边界.

G-N 边界表达式见文献[11-12],它所使用的辅助变量递推关系为

逐步消去辅助变量,可得式(13)的合并形式为

将式(14) 与式(2) 进行比较可知它们是一致的.式(14)中的计算波速参数c1,c2,···,cN对应于式(2)中的人工波速参数ca1,ca2,···,caN.所以在声波问题中,相同参数的G-N 边界与caj-Higdon 边界具有相同的精度.

H-W 边界表达式见文献[13-14],它基于的辅助变量递推关系为

B´ecache 等[38]经过推导,证明式(15) 的合并形式为

式 (16) 中计算波速由声波波速c和任意系数a0,a1,a2,···,aN组合而成,为c/a0,c/a1,c/a2,···,c/aN.比较式(16) 和式(2) 可知,在声波问题中,1,2,3,··· 阶H-W 边界分别对应于3,5,7,··· 阶caj-Higdon 边界,即存在2N+1 倍关系.

3.4 任意大角度波动方程

Guddati 等[15-16]提出一种矩阵形式的辅助变量高阶边界AWWE,他们通过两种方式导出这种边界:一是将声波方程的旁轴近似递推关系按偶数阶提取出来,将其写成矩阵形式,再返回到时域成为人工边界条件;二是将半无限介质分层,利用各层动力刚度的有限单元近似,在频域内导出矩阵关系式,进而返回时域得到人工边界条件.根据本文第3.2 节所指出的声波旁轴近似边界可以由caj-Higdon 边界与声波方程相结合得到,可以推断出1,2,3,··· 阶AWWE边界分别对应于2,4,6,··· 阶caj-Higdon 边界,即存在2N倍关系.AWWE 边界与caj-Higdon 边界的等价中间形式为

3.5 其他相关边界方法

Lindman[17]提出一种数值优化边界,表达式见文献[17].Kausel[39]证明了该边界的连续形式为一种高阶旁轴近似边界,与Clayton-Engquist 声波边界属于一个类型.

Peng-Toks¨oz[18]提出另一种数值优化边界,表达式为

式中un(i,j,k)=u(iΔx,jΔy,kΔz,nΔt)表示离散网格中波场的节点值.编号0 为边界节点,1 和2 为临近的内域节点.a01,a02,···,a22为控制参数,由一系列限定条件建立的方程组来确定.显然,式(18)借鉴了二阶Higdon 边界数值离散格式所采用的3×3 节点组合,只是确定计算系数的方式有所不同.

Randall[19]针对弹性波边界精度不高的问题,提出一种将弹性波场进行势分解,对分解出来的波势函数分别使用高精度的Lindman 声波边界,再将结果合成为弹性波场的方法.弹性波场向势函数的分解与合成分别表示为

式中,u为弹性波位移场向量,Φ和Ψ分别为P 波和S 波势函数.∇为向量微分算子,cp和cs分别为P 波和S 波波速.Randall 势分解方法是提高外推型ABC 模拟弹性波的精度的一个有益尝试,但由于势分解与合成过程以及Lindman 边界所带来的双重复杂度,作者未见其他学者应用或探讨该边界方法.

Liu 和Sen[20-21]将外推型ABC 与过渡区技术相结合,发展出一套混合边界方法,如图4 所示.该方法在计算模型的外边界附近设置一定宽度的过渡区,过渡区波场由内域离散格式和边界数值格式分别计算的波场进行加权和得到.二者权系数之和始终为1,其分配比例沿着过渡区宽度的变化如图中三角形所示.

图4 Liu-Sen 混合边界示意图Fig.4 Schematic of Liu-Sen hybrid ABC

大量波动模拟实践表明Liu-Sen 混合边界能够显著提高人工边界的精度,并有效降低外推型ABC可能面临的数值失稳风险.过渡区的引入在相当程度上阻止了反射误差向内域传播,并延缓了不稳定结果的积累.在数值试验中观察到过渡区内波场表现出类似完美匹配层内波场的快速衰减特征.我们认为本文外推型ABC 与Liu-Sen 混合边界方法相结合,应是解决波动数值模拟中人工边界问题一个值得重视的方向.

4 外推型ABC 的精度控制原理

4.1 离散边界公式(1)的时空外推原理

离散边界公式(1)的单人工波速形式(即廖氏透射边界)存在一种明确的几何解释,如图5 所示.外行波沿人工边界节点所在的局部坐标空间轴s的传播过程,可以在sot时空坐标平面上展开成一个二维“静态”波场.在这个二维静态波场中,沿着由单一人工波速和离散时间步距所确定的时空外推方向,可以截取出一维波场曲线g.人工边界节点运动的实际值由连续曲线g决定,但在数值模拟中只能利用曲线g上若干个参考点的值来推算其端部人工边界节点的值.这个推算过程的精度即为MTF 边界的精度.

从图5 中可以看出,边界阶次N(即参考点数目)代表基于N-1 阶多项式的外推精度.曲线g的弯曲程度则由外行波复杂程度和所选取的时空外推方向(它与人工波速有关) 共同决定.于是这种几何解释可以用一句话概括为:N阶MTF 是以基于N-1 阶多项式的外推精度对由人工波速所确定的时空外推方向上的外行波曲线的一种近似.显然,提高边界阶次N,或者选取与实际视波速相接近的人工波速以获得较为平缓的曲线g,都是提高MTF 边界精度的有效措施.

图5 MTF 边界的时空外推原理Fig.5 Time-space extrapolation of MTF boundary

MTF 边界这种明确的几何解释,加之其简洁的数学形式和“多次透射”的物理机制,奠定了外推型人工边界理论的思想基础.

4.2 连续边界公式(2)的反射系数乘积原理

连续边界公式(2) 是Higdon 边界的统一形式.Higdon 边界所采用的多个单向波动微分算子乘积形式是从满足边界要求的一系列离散表达式中总结而来的[9].这表明连续边界公式(2)的理论源头可以追溯到离散边界公式(1),即上述MTF 边界的时空外推原理.因此,外推型人工边界条件具有统一的思想基础,它来源于离散边界.各种微分方程形式的连续边界主要为外推型ABC 提供了丰富的变化形式.

连续边界公式(2) 的精度可以利用反射系数加以分析.该边界对于入射角为θ 的外行平面波的反射系数为

式中c为介质物理波速,caj为边界的人工波速参数.式(21)为多个因子的乘积,这对应了式(2)的单向波动微分算子乘积.图6 绘出了几种不同N和caj取值的反射系数R(θ).

在图6 中,根据一阶边界反射系数R1 可知单个反射因子的幅值总是小于或等于1,所以高阶边界所具有的多个反射因子的乘积能够迅速提高边界精度.另外,多个反射因子如果具有不同的零反射角度(由不同人工波速取值决定),能够在更大透射角范围内提高边界精度.

图6 连续边界公式的反射系数乘积原理Fig.6 Product of reflection coefficients of multiple first-order wave-propagation operators

图7 绘出了caj-Higdon 边界(Hig)、Hagstrom-Warburton 边界(HW)、AWWE 边界的反射系数.根据第3.3 节和第3.4 节内容可知,Givoli-Neta 边界(GN)反射系数与caj-Higdon 边界相同,故未在图中列出;N阶HW 边界、AWWE 边界分别对应于2N+1 阶、2N阶caj-Higdon 边界.

图7 曲线表明相同人工波速取值下,1,2,3 阶HW 反射系数与3,5,7 阶Higdon 相同,1,2,3 阶AWWE 反射系数与2,4,6 阶Higdon 相同.这证实了前文理论分析所给出的结论.由此可见,GN、HW和AWWE 3 种辅助变量高阶边界使用与caj-Higdon相同数量的一阶单向波动微分算子(∂/∂t-caj∂/∂s)实现了与之相同的精度,即它们提高边界精度的原理与caj-Higdon 并无二致.

图7 Higdon,G-N,H-W,AWWE 边界的反射系数Fig.7 Reflection coefficients of Higdon,G-N,H-W,AWWE boundaries

4.3 基于视波速的多人工波速分析

基本边界公式(1)和(2)所采用的多人工波速参数caj(j=1,2,···,N)并不仅仅是形式上的优化,它对于具有多种物理波速的复杂波动问题具有非常重要的实用价值.这些问题包括多层介质问题、弹性波、两相介质波动、频散波动等.上述基于透射角度θ 的反射系数式(21)未考虑物理波速的变化.为同时考虑介质物理波速和透射角度的变化,这里采用外行波的视波速cn=c/cos θ 作为基本变量,得出caj-Higdon边界反射系数的另一种表达形式

图8 所示为存在两种物理波速时caj-Higdon 边界的反射系数幅值与外行波视波速cn的关系.这里两种物理波速为c1=v和c2=3v,v表示某个可约去的波速,它的取值不影响分析结果.3 倍的物理波速差异将对人工边界性能提出挑战.

在图8 中,主要波动能量的视波速范围有两段,其中虽然有透射角度的影响,但差异较大的两种物理波速起了主要作用.理想的人工边界反射系数应当在所关心的两个区段中都处于低值.可以看出,采用单一人工波速的3 种边界R(v,v),R(2v,2v) 和R(3v,3v)都不够理想,只有采用多人工波速的二阶边界R(v,3v)和三阶边界R(v,3v,3v)较好地兼顾了对两段波动能量的吸收.因此,在具有多种差异较大的物理波速的波动问题中,使用以多人工波速为参数的外推型ABC 很有必要也非常实用.

图8 多人工波速分析Fig.8 Analysis of various artificial wave velocities

5 数值算例

5.1 声波大角度透射问题

以一个均匀介质声波问题算例来检验上述外推型ABC 的精度.模型尺度为200 m×500 m,介质波速为c=500 m/s.顶面为自由边界,左、右、底边界为人工边界.在顶部中点施加中心频率为10 Hz 的Ricker 子波激励.计算模型的空间、时间离散均采用标准的二阶中心差分格式[9].网格尺寸为2 m×2 m,时间步长满足稳定条件cΔt/Δx≤0.71.人工边界分别采用本文提出的外推型ABC 基本公式caj-MTF 和caj-Higdon,以及现有的Clayton-Engquist (CE) 旁轴近似边界[7]、辅助变量高阶边界Givoli-Neta(GN)[11]、Hagstrom-Warburton(HW)[13]和AWWE[16].本文ABC 的数值离散格式见附录.

图9 给出了0.9 s 时的波场快照.此时外行波在侧边界上的透射角度约75°,为大角度透射,这种情形对ABC 性能要求较高.图中括号内参数如(c,1.5c,2.5c)表示所采用的人工波速参数caj(j=1,2,···,N).

图9 采用不同外推型人工边界的声波模拟结果Fig.9 Modeling results of acoustic wave propagation using different ABCs

图9 中不同ABC 的反射误差对比证实了前文理论分析所给出的结论,主要包括:(1)caj-MTF 和caj-Higdon 等价(见第3 部分),这在本算例中表现为尽管用于实现caj-MTF 和caj-Higdon 的数值计算格式有所不同,但只要当二者所采用的边界参数N和caj(j=1,2,···,N)相同时,其模拟结果就会保持一致;(2)随着阶次N升高,边界精度提高,且多人工波速方案能够进一步改善精度(见图6),图9(a)～图9(e)中不同边界阶次结果对比和相同阶次下多人工波速方案与单一人工波速方案结果对比证实了这个结论;(3)CE 边界和辅助变量高阶边界GN,HW 和AWWE 具有与基本边界公式caj-Higdon 相同的精度控制原理(见第3.2～3.4 节和图7),这预示着相同边界参数下它们的模拟结果一致,如图9 中二阶CE 与二阶caj-Higdon在0.9 s 时的结果非常接近,二阶HW(相当于五阶caj-Higdon)和三阶AWWE(相当于六阶caj-Higdon)模拟结果与三阶caj-Higdon 结果相似.不同外推型ABC除了呈现出上述共同特征之外,它们各自还有一些不同的表现,如图中CE 边界存在角点反射问题(见0.5 s 快照图),GN 边界出现了失稳(该失稳首先出现在边界上距离震源较近的区域,为典型的人工边界条件失稳特征),HW 和AWWE 边界由于数值离散误差影响显著,其三阶以上的高阶精度难以实现.这些问题从侧面说明本文所提出的基本边界公式(1) 和式(2)最为简单实用.

5.2 均匀介质弹性波问题

以一个均匀介质弹性波问题算例来进一步检验上述外推型ABC 的性能.模型尺度为2500 m ×2500 m,介质纵、横波速为cp=2000 m/s,cs=1000 m/s.四边均为人工边界.在模型中心施加中心频率为10 Hz 的Ricker 子波激励.计算模型的空间、时间离散均采用标准的二阶中心差分格式[35].网格尺寸为5 m×5 m,时间步长满足稳定条件Δt/Δx≤1.人工边界分别采用caj-MTF 和caj-Higdon,以及第3.2 节提及的几种弹性波旁轴近似边界.此时附录A 中的边界离散格式分别用于人工边界节点的水平和竖向波场分量.

图10 给出了0.9 s 和1.5 s 时水平分量的波场快照.此时存在两种差异较大的物理波速(cp等于2倍cs),这种情形能够凸显采用多人工波速方案的必要性及优势.caj-Higdon 与caj-MTF 结果相同,未在图中列出.考虑了3 种单一人工波速和一种多人工波速的二阶caj-MTF 边界,并以几种旁轴近似边界作为对比.

图10 中caj-MTF 边界模拟结果证明了对于多物理波速的复杂波动问题采用多人工波速外推型ABC的必要性及优势(第4.3 节).此时3 种单一人工波速的二阶caj-MTF 均不能很好地实现对P 波和S 波的同步透射,而多人工波速(cs,cp) 方案几乎完美地实现了这一目标.与基本边界caj-MTF 相比,所谓的弹性波旁轴近似边界反射量很大,这对应了图3 关于这几个弹性波边界理论上具有不合理性的分析结果.Clayton-Engquist 和Fuyuki-Matsumoto(它与前者高度相似) 边界出现失稳,Emerman-Stephen 和Stacey 边界(文献[36-37] 指出它们改善了CE 弹性波边界的稳定性)以及本文提出的caj-MTF 边界在20 s 长持时模拟结果中保持了稳定.因此,综合考虑精度与稳定性,模拟弹性波的外推型ABC 建议采用基本边界公式caj-MTF 或caj-Higdon.

图10 采用不同外推型人工边界的弹性波模拟结果Fig.10 Modeling results of elastic wave propagation by using different ABCs

5.3 纵、横向不均匀介质弹性波模拟

设计如图11 所示的纵、横向不均匀介质场地模型,检验本文边界条件的吸收效果.该模型为右下方整块基岩与左侧、上部成层土体形成的复杂场地,岩土体参数如图所示.在(x,z)=(300 m,100 m)处节点的竖向分量上施加中心频率为10 Hz 的Ricker 子波激励.上部为自由地表,其余3 边分别采用自由边界(作为对照) 或本文提出的人工边界,进行两组模拟.采用集中质量有限元与中心差分相结合的数值模拟方法[29-30].单元尺寸为2.5 m×2.5 m,时间步长满足稳定条件

图11 纵、横向不均匀介质模型Fig.11 Computational model of a longitudinal and transverse inhomogeneous media

这个模拟选用所提出的caj-MTF (或caj-Higdon,与前者结果一致)边界,边界参数为N=2,caj=cs,cp.不考虑其他外推型ABC,因为前文已指出:Givoli-Neta,Hagstrom-Warburton,AWWE 等高阶辅助变量边界和Lindman 数值优化边界只适用于声波;算例5.2 表明几种旁轴近似弹性波边界的精度和稳定性远不如caj-MTF 或caj-Higdon;引言中提及的Peng-Toks¨oz 数值优化边界和Randall 势分解边界十分复杂,少有应用和讨论,其是否有精度或稳定性方面的优势尚无定论;Liu-Sen 混合边界是外推型ABC 与过渡区技术相结合的综合性解决方案,已发表大量文献,而本文关注的仅是外推型ABC 本身的性能.模拟结果如图12 所示.

在图12 中,由于自由边界会完全反射波动能量,该模型像一个封闭波动能量的“盒子”,其模拟结果与实际问题中波可以无阻挡地穿过(因为原问题中此处无界面)人工边界再传向远处的过程相去甚远.caj-MTF(或caj-Higdon)边界则几乎看不到虚假反射波,说明它们所推算的边界节点运动与原来无限介质中该处运动非常接近(或基本一致),因此这个有限计算模型给出了与实际问题相符的模拟结果.

图12 不均匀介质弹性波模拟结果Fig.12 Simulation results of elastic waves in the inhomogeneous model

5.4 关于复杂介质模拟效果和稳定性的讨论

外推型ABC 模拟复杂介质的效果已被现有文献给出的各种波动问题模拟结果所证实,如文献[4]混凝土坝开裂模拟中的边界处理,文献[11] 中频散介质波动模拟,文献[16] 对力学性质渐变介质的反演,文献[21] 对三维VTI 介质的波动过程正演,文献[31]中适用于两相介质波动问题的边界,文献[35]对Rayleigh 波边界的探讨等.在本文建立的外推型人工边界条件理论与公式体系中,这些复杂波动过程的绝大多数特征已得到充分考虑,具体为:(1)多种人工波速的参数配置能够很好地考虑由多种物理波速(包括弹性波、两相介质波动、频散波动、介质交界面等)以及不同透射角度所决定的沿外推型ABC计算方向视传播速度的复杂性;(2)由几何扩散、阻尼效应或其他原因引起的波场的衰减性,可以由添加小量修正的边界表达式(3) 和式(4) 很好地模拟,或者只是单纯地将边界(1)和边界(2)的阶次提得足够高,也能够达到满意精度(其原理参考图5);(3)对于角点、介质交界面、曲线边界等不同几何构型的适应性,caj-MTF 和caj-Higdon 可以无差别地应用,因为它们只涉及到一条由边界节点指向内域的离散网格线上的信息,与边界形状无关,而其他外推型ABC若涉及到边界本身的信息,则会受到边界几何构型的制约.除了上述特征之外,在不同复杂波动问题中是否还存在其他影响边界精度的重要因素,如多种波动成分之间的耦合效应等,将在后续工作中具体探讨.

数值试验表明,本文基本边界公式caj-MTF 和caj-Higdon 的稳定性高度一致,且与现有MTF 和Higdon 边界基本没有区别.外推型ABC 的稳定性可以从所引文献及相关学者的研究成果中总结得到,这个研究课题目前仍然活跃.本文和相关文献数值试验表明了外推型ABC 的实用性,或者说它们一般能够在所关心的模拟时间内给出稳定结果,这对于解决大多数持时不太长的波动问题已经足够.对外推型ABC 稳定性问题完整而详细的探讨显然不是本文所能完成的任务,不过,关于稳定性、边界失稳特征及其解决方法,作者从已有研究中总结出以下几点主要认识:(1)不同于应力型ABC(如黏弹性边界)和衰减层型ABC(如完美匹配层边界)需要附着在内域离散模型之上,外推型ABC 有着独立于内域离散模型的数值计算格式,后者两种不同离散格式的耦合常常不如前两者的整体计算格式稳定.(2)当边界出现某种失稳时,其不稳定值的积累速度由边界条件的反射幅值(高阶边界大于低阶边界)和波在整个模型与所有失稳节点之间的反复反射放大共同决定,因此三维模型(边界节点多) 失稳放大最快,二维模型次之,一维模型最慢;同样维度下,小尺寸模型(来回反射距离短)比大尺寸模型失稳积累得更快;高波速问题(来回反射所需时间少)比低波速问题失稳积累得要快;复杂介质问题(多个界面之间反射次数多)比均匀介质问题失稳积累得快.(3)若边界出现漂移失稳(指边界节点运动向正值或负值单方向的不正常累积),通常采用修正的边界公式(3)或(4)来解决,亦可考虑将高阶边界与低价边界组合使用,或者在边界上附加弹簧和阻尼元件等.(4)若边界出现振荡失稳,可以考虑在整个计算模型中增加少量阻尼来解决,或者在与一些特定的内域离散格式结合时能够避免这一问题,此外,增大计算模型,或使用具有耗能特性的时间积分格式,或采用Liu-Sen 混合边界方法等都是值得考虑的选择.

6 结论

本文新提出离散形式和连续形式 (即微分方程形式) 两个基本边界公式,通过证明二者的等价性并阐明其与经典的廖氏透射边界、旁轴近似边界、Higdon 边界以及Givoli-Neta、Hagstrom-Warburton、AWWE 辅助变量边界等之间的理论联系,据此建立了一种外推型人工边界条件的理论与公式体系.主要研究结论如下:

(1) 上述经典人工边界条件可以统一地从外推型ABC 角度来讨论和应用.因为它们计算人工边界节点运动的方式都是利用边界附近一组节点的运动(不同边界使用的节点组合可能有一定差异)进行时空外推,使用相同的控制参数(边界阶次和一组计算波速),并且在数值模拟中表现出相似的精度和稳定性.

(2)本文提出的基本边界公式在精度和稳定性方面均更有优势,是该类边界最简单实用的形式.

(3) 真正与外推型ABC 不同类型的是以黏性边界、黏弹性边界等为代表的应力型边界[40-41]和以函数衰减层、完美匹配层等为代表的衰减层型边界[25-26].在后续研究和应用中,外推型ABC[42-43]的稳定性、高阶应力型边界或完美匹配层边界应用于不同波动问题的实现方式与计算效率、融合不同方法的人工边界问题综合性解决方案等,都是值得继续探讨并更好地解决的重要问题.

附录: caj-MTF 边界的数值离散格式

图A1 caj-MTF 边界数值计算格式的局部节点编号Fig.A1 A1 Local index numbers of the numerical scheme of caj-MTF boundaries