VTI介质中多波多分量数值模拟及PML吸收边界条件
VTI介质中多波多分量数值模拟及PML吸收边界条件张为彪;洪宇;郑洁;夏弋峻;郭飞;邹清文【摘要】为了充分认识多波多分量地震勘探技术中各种波及其分量的传播特点,以便深入挖掘来自地下复杂介质信号中的有效信息,解决实际勘探中遇到的复杂问题,选择了一种特殊的各向异性介质——具有垂直对称轴的横向各向同性介质(VTI)中传播的波及其分量作为研究对象,对其波场进行了数值模拟,分析总结了其传播特点.为了解决边界反射的干扰问题,构建并应用了相应的PML吸收边界条件,模拟结果表明,吸收效果较好.【期刊名称】《石家庄经济学院学报》【年(卷),期】2017(040)005【总页数】6页(P7-12)【关键词】多波;多分量;各向异性;吸收边界条件;数值模拟【作者】张为彪;洪宇;郑洁;夏弋峻;郭飞;邹清文【作者单位】中海石油(中国)有限公司深圳分公司,广东深圳 518054;中海石油(中国)有限公司深圳分公司,广东深圳 518054;中海石油(中国)有限公司深圳分公司,广东深圳 518054;中海油能源发展股份有限公司工程技术分公司,广东深圳518054;中海石油(中国)有限公司深圳分公司,广东深圳 518054;中海石油(中国)有限公司深圳分公司,广东深圳 518054【正文语种】中文【中图分类】TE132随着油气勘探程度的增加，勘探难度不断加大，出现了很多具有挑战性的地质勘探问题，目前最为常用、理论和技术发展也最为成熟的地震勘探方法—纵波勘探已无法完全解决这些问题，而多波多分量地震勘探技术在解决这些问题上有其独特的优势，如：消除模糊带、识别真假亮点、探测裂缝、获得与横波速度有关的各种岩性参数等。
波场数值模拟是人们理解波在不同介质中传播规律，认识和正确解释采集波信号携带的重要信息的重要手段[1-2]。
为了理解多波多分量在地下介质中的传播特性，需对其波场进行数值模拟，模拟时，为了避免人工边界产生的反射干扰介质内的波场，必须引入吸收边界条件来处理边界。
Cerjan等人[3]通过在人工边界附近引入损耗介质来衰减向外传播的波，该方法的缺点是增加了较大的计算量。
另一种比较著名的吸收边界条件是Clayton吸收边界条件[4]，该吸收边界条件在旁轴近似理论的基础上导出，在特定的入射角和频率范围内，具有较好的吸收效果。
Berenger[5]针对电磁波传播情况，给出了一种高效的完全匹配层吸收边界条件，并在理论上证明该方法可以完全吸收来自各个方向、各种频率的电磁波，而不产生任何反射。
这种方法很快被引入到地震波场的模拟中，Frank（1996）把PML应用到一阶弹性波方程波场的模拟中[6]， Liu Q H（1999）把PML应用到球面坐标系和圆柱坐标系中的弹性波方程的波场的模拟中[7]，王守东（2003）针对声波方程，给出了完全匹配层吸收边界方法[8]。
为了充分认识多波多分量地震勘探技术中各种波及其分量的传播特点，以便深入挖掘来自地下复杂介质信号中的有效信息，解决实际勘探中遇到的复杂问题，需要对地下复杂介质中的各种波及其分量进行数值模拟，本文选择了一种特殊的各向异性介质—具有垂直对称轴的横向各向同性介质（VTI）中传播的波及其分量作为研究对象，对其波场进行了数值模拟和分析。
为了解决边界反射的干扰问题,构建并应用了相应的PML吸收边界条件。
均匀各向异性介质中的弹性波方程[9]为：其中，x1，x2，x3是质点在空间坐标系中的三个坐标分量，u1，u2，u3是质点沿空间坐标系三个坐标轴方向的位移，Cijkl是弹性系数。
观察（1）式可知，共有81个弹性参数。
但是，这81个弹性参数并不是完全独立的，其实，完全独立的弹性参数只有21个，即便如此，要求解该方程也是极其复杂的。
实际上，在地震勘探中，经常遇到的是一种特殊的各向异性介质，即横向各向同性（TI）介质，这种介质的各向异性主要是由裂隙和薄互层引起的，具有垂直对称轴的TI介质称为VTI介质，具有水平对称轴的TI介质称为HTI介质，HTI介质可以看成是VTI介质的对称轴旋转90o得到的。
对于VTI介质，在x-z平面内，（1）式可简化为：其中，ρ是密度，C11，C13，C33，C44，C66是弹性系数。
上面的方程组可以分解为两部分，一部分为方程（3），它为VTI介质中的SH波的方程，另一部分为方程（2）和方程（4）组成的方程组，它们为VTI介质中的准P波和准SV波的波动方程，这两种波是耦合在一起的，即在VTI介质中，P波传播可以引起SV波，SV波传播也可以引起P波。
（一）VTI介质中SH波波场模拟1．VTI介质中SH波方程差分格式根据1知，VTI介质中SH波的方程为：定义C=( C66/C44)1/2，V=( C44/ρ)1/2，则上述方程化为：式中，V为VTI介质中的SH波沿垂直于层面方向的速度，C是该介质的各向异性参数，表示VTI介质中地震波沿水平方向和沿垂直方向的速度比。
对于方程（6），时间导数采用二阶中心差分，空间导数采用2N阶差分精度的交错网格差分，可得到如下的差分格式[10-17]：式中，[u2(i,j,k)]、[u2(i,j,k)]同上。
2．VTI介质中SH波方程的PML吸收边界条件通过引入中间变量u2x，u2z，A2x和A2z，可以把式（6）表示为下面的等价方程组：很容易证明（8）与（6）是等价的。
由（8）式来构建新的方程组，即引入函数d1(x)和d2(z)，把方程组（8）变为：（9）的解与（8）的解相比是逐渐衰减的，函数d1(x)和d2(z)分别起到了x和z方向的衰减系数的作用。
这样可以利用（9）式采用下面的方式进行数值模拟：在模拟区域的四周引入完全匹配层，在模拟区域内，令d1(x)和d2(z)为0（此时式（9）等价于式（6）），波不受到衰减，在模拟区域外的完全匹配层内，d1(x)和d2(z)不为0，传播到完全匹配层内的波受到衰减，当波传播到完全匹配层的边界时，能量已经很弱，边界反射几乎为0，这样就可以得到正确的模拟结果。
（9）式即为VTI介质中SH波方程的PML吸收边界条件。
实际计算的时候，衰减系数d1(x) （d2(z)类似）采用下面的形式：其中，w为PML层的厚度，x为PML层内网格节点到PML内侧边界的距离，从上式中可看出，在PML内侧边界x为0处, 衰减系数d1(x)=0，在PML外侧边界x为w处，衰减系数d1(x)达到最大，为dmax，dmax可根据波传播一个网格的距离其振幅衰减为原来的十分之一来计算。
3．波场的数值模拟首先建立宽度和深度均为2000m的VTI介质模型，模型中，密度ρ=2.44g/cm3，弹性系数C44=7.0×109 N/m2，C66=21.0×109 N/m2，震源位于点（1000m，1000m）处，采用垂直于x-z平面的方向的集中力源作为震源[9]，选用30Hz的Ricker子波作为震源子波。
模拟波场时，x方向和z方向的步长均取为8m，t的步长取为1ms，空间导数采用10阶差分精度的差分格式。
波场的模拟结果见图1，其中(a)为t 等于250ms时SH波波场的等值线，(b)、(c)和(d)为t =530ms时SH波的波场等值线，(b)没加边界吸收条件，(c)加了阻尼衰减吸收边界条件，(d)加了上面推导出的PML吸收边界条件。
从模拟结果可以得出以下结论：（1）SH 波独立传播，且只存在y分量；（2）SH波沿横向的传播速度大于沿纵向的传播速度；（3）不加吸收边界条件时，边界反射干扰很严重，加了阻尼衰减吸收边界条件后，边界反射得到了一些压制，但压制效果不好，采用上面推导出的PML边界吸收条件后，边界反射得到了很好的压制，几乎完全消失了。
（二）VTI介质中准P-SV波的波场模拟1．VTI介质中准P-SV波方程的差分格式根据1知，VTI介质中准P-SV波的方程为：推导出（11）式的交错网格高阶有限差分格式[10-17]如下：其中，[ u( i, j, k) ]、 L2[ u( i, j, k )]同上， y2．PML吸收边界条件引入中间变量u1x，u1z，u3x，u3z，A1，A1x，A1z，A2，A2x，A2z，A3，A3x和A3z，把（11）式变为下面的等价方程组：容易证明方程组（13）和（11）是等价的。
引入函数d1(x)和d2(z)，把方程组（13）变为：该方程组的解与式（13）的解相比也是逐渐衰减的，用式（14）按上文中同样思路进行数值模拟，可有效压制边界反射。
（14）为VTI介质中准P-SV波方程的PML吸收边界条件。
3．准P-SV波波场的数值模拟建立图2所示的模型，模型中，ρ=2.44g/cm3，C11=25.5×109 N/m2，C13=1.0×109 N/m2，C33=18.4×109 N/m2，C44=5.6×109 N/m2，震源位于（1500m，1500m）处，采用倾斜作用集中力源（与x轴和z轴成45°）作为震源[9]，选用20Hz的Ricker子波作为震源子波。
模拟波场时，x方向和z方向的步长均取为10m，t的步长取为1ms，空间导数采用12阶差分精度的差分格式。
波场的模拟结果见图3和图4。
图3为模拟得到的VTI介质中准P-SV波波场的x分量，(a)为 t为 350ms时波场，(b)、(c)和 (d)都是 t为700ms时波场，(b)没加边界吸收条件，(c)加了阻尼衰减吸收边界条件，(d)加了上面推导出的PML吸收边界条件。
图4为模拟得到的 VTI介质中准P-SV 波波场的z分量， (a)为t为350ms时的波场，(b)、(c)和(d)都是t为700ms时的波场，(b)没加边界吸收条件，(c)加了阻尼衰减吸收边界条件，(d)加了上面推导出的PML吸收边界条件。
从这两组图中可以看出：（1）qP波和qSV波是耦合在一起的，且都存在x和z两个分量，波场较复杂；（2）qP波的传播速度大于qSV波的传播速度；（3）qSV波在与集中力源倾斜的方向正交的对角线方向上会出现三分叉现象；（3）不加吸收边界条件时，边界反射很严重，阻尼衰减吸收边界条件对边界反射有一些压制作用，但压制效果不理想，而PML边界吸收条件对边界反射压制效果较好。
1．纵波勘探基于声波方程，只考虑单一类型的纵波，而多波多分量地震勘探充分利用地下复杂介质中传播的各种波及其分量，以解决纵波勘探难以解决的一些复杂问题。
本文的模拟表明，在VTI介质中同时存在SH波、qP波和qSV波三类波，这三类波具有如下特点：（1）SH波是独立传播的，只存在y分量，且沿横向的传播速度大于沿纵向的传播速度；（2）qP波和qSV波是耦合在一起的，且都存在x和z两个分量，波场较复杂；（3）qP波的传播速度大于qSV波的传播速度；（4）qSV波在与集中力源倾斜的方向正交的对角线方向上会出现三分叉现象。
PML吸收边界条件下TM波的时域有限差分法分析
2) ∆ t 可由 Courant 稳定性条件得到。即：
对于二维情形:
∆t≤
1
(9)
c0
1+1 (∆x)2 (∆y)2
若 ∆ x= ∆ y= δ ，上式化为:
∆t≤ δ
(10)
c0 2
在实际计算中，我们一般取:
∆ t= δ
(11)
2c0
-2-
3. PML 吸收边界条件
4. 3 仿真分析
1)T=20 步的情况
-4-
2)T=40 步的情况
图 4 T=20 步时场传播情况
3)T=80 步的情况
图 5 T=40 步时场传播情况
图 6 T=80 步时场传播情况
图 6 看上去出现了反射，但是我们应该看到，这些反射都是出现在 PML 区域（共有八 层），而在场区没有出现任何反射，右侧图可以证明这一点。 3) T=120 步的情况
时域有限差分法的主要思想是把 Maxwell 方程在空间、时间上离散化，用差分方程代替 一阶偏微分方程，求解差分方程组，从而得出各网格单元的场值。FDTD 空间网格单元上电 场和磁场各分量的分布如图 1 所示。
图 1 FDTD 网格中电磁场分量的分布
电场和磁场被交叉放置，电场分量位于网格单元每条棱的中心，磁场分量位于网格单元 每个面的中心，每个磁场（电场）分量都有 4 个电场（磁场）分量环绕。这样不仅保证了介 质分界面上切向场分量的连续性条件得到自然满足，而且还允许旋度方程在空间上进行中心 差分运算，同时也满足了法拉第电磁感应定律和安培环路积分定律，也可以很恰当地模拟电 磁波的实际传播过程[2]。
于 1，这并不奇怪，如入射方向很接近左右侧分界面的波将最终几乎垂直地传到上或下侧边 界，在那里被吸收掉。
利用高阶交错网格有限差分法数值模拟VTI 介质井孔声场
声波测井是一种利用接收器接收井中声源激发 的弹性波信号获取井周围岩石介质性质的一种常用 的井中地球物理方法。它在煤田测井、石油勘探、 工程地质勘察等领域具有广泛的应用。井孔声场的
摘要: 横向各向同性(TI)介质是岩石地球物理中常见的一种现象，研究其井孔声场传播特征对声波
测井理论以及为声波测井解释提供依据具有重要意义。针对具有垂直对称轴的横向各向同性(VTI) 介质，根据柱坐标系条件下的弹性波波动方程，推导了速度—应力交错有限差分公式，采用时间 二阶、空间十阶的交错有限差分算法对 VTI 介质中的井孔声场进行数值模拟。给出了在均匀介质 中井孔声场不同时刻的波场快照，以及不同各向异性系数的 VTI 介质中的波场快照，计算了井轴 上声源激发出的声波全波列波形。结果表明，在其他条件不变的条件下，VTI 地层的各向异性系 数的增大对横波的传播影响不大，但会使得纵波在纵向上的传播速度相对变小，径向上变化不大。 各向异性系数的增大会使声波测井全波列首波信号时差变大，声波幅度略变小。
数值模拟是研究井中声场传播特征的有效途径，通 过井孔声场模拟可以了解井中各种声场的传播模 式，丰富测井方法，为声波测井解释提供依据。为 此，很多研究人员通过各种途径研究声场传播特征，
收稿日期: 2015-12-23 基金项目: 国家青年自然科学基金项目(41504090)；中央高校基金基础研究项目(310826161021, 310826161014) Foundation item：National Youth Natural Science Foundation of China(41504090); The Basic Research Project for the Central Universities of China(310826161021, 310826161014) 作者简介: 岳崇旺(1981—)，男，河南通许人，博士，从事地球物理测井方法相关研究. E-mail：ycw1981@ 引用格式: 岳崇旺，王飞. 利用高阶交错网格有限差分法数值模拟 VTI 介质井孔声场[J]. 煤田地质与勘探，2016，44(4)：125–131 YUE Chongwang，WANG Fei. The simulation of acoustic wave propagation in the borehole surrounded by vertical transversely isotropic (VTI) media using staggered-grid high-order finite–difference method[J]. Coal Geology & Exploration，2016，44(4)：125–131.
VTI介质多波速度分析方法及应用
VTI介质多波速度分析方法及应用摘要：在地层存在广泛各向异性的基础上，常规的速度分析已经不再适用。
本文以VTI介质为研究对象，对影响速度分析的各向异性参数进行了综合分析，且在提高速度分析的精度方面有了自己的见解。
首先研究了各向异性介质的基础理论和地震波传播特点，其次研究了各向异性介质下时距曲线的拟合公式，得出地层参数与速度的关系，再次研究转换波速度分析引入单平方根方程法。
利用交互式速度分析方法对模型进行了计算处理，结果表明该方法的正确性。
关键词：各向异性；VTI介质；速度分析；非双曲线时距曲线随着地震勘探技术的发展需求，各向异性逐渐成为地震勘探的前沿课题之一。
地震波在地下介质中的传播速度是地震数据处理和解释中非常重要的参数。
在早期求取这一参数的过程中，人们把地球介质假想成为各向同性，推动了石油地震勘探的发展。
但从理论和实际角度出发，各向异性广泛存在于地球地层中。
基于各向同性的地震勘探处理方法就不再适用于各向异性介质。
因此，发展各向异性的地震勘探处理方法是非常必要的。
由于地层存在各向异性，速度的提取就与地层各向异性有关，这个特点就是速度各向异性。
如果忽略速度各向异性就会造成地层速度提取不准确，进而影响以地层速度为基础的后续处理步骤[3]。
综上所述，各向异性地震勘探处理方法不仅是对理论研究的完善，更是实际应用迫切的需要[4]。
本文以VTI介质为理论为基础，研究了各向异性介质非零偏多波时距曲线公式，在此基础上将速度谱技术引入到非零偏多波速度分析中，研究了各向异性多波速度分析方法和步骤，建立了各向异性介质叠加速度模型。
一、各向异性介质速度分析原理1.各向异性介质基本定义在各向同性弹性体中，介质内任意一点沿不同方向弹性性质都一样，但也有许多物体其弹性性质在各个方向上不尽相同，因而是各向异性的，这种不同是材料本身的性质决定的，由介质的本构方程来描述。
各向同性材料中地震波的传播规律和各向异性材料中地震波的传播规律是不同的。
各向异性PML吸收层
ε 2 = ε1s ;
μ2 = μ1s ;
⎡ sx −1 ⎢ s =⎢ 0 ⎢ 0 ⎣
0 sx 0
0⎤ ⎥ 0⎥ sx ⎥ ⎦
(7.47)
无反射特性与入射角,极化特性,入射波的频率无关.且由(7.39)和(7.40),TE 和 TM 极化波的波 动特性是相同的.按照其单轴各向异性和完美匹配性质(横向特征阻抗完全匹配)可称之为”单 轴 PML”. 与 Berenger 的 PML 类似,2 中 UPML 的无反射特性对任意 sx 均成立.例如,选择
PML : − js β x − j β y − jβ x− jβ y H z = H 0 e x 1 x 1 y = H 0 e 1 x 1 y e −σ x xη1 cosθ − jβ x − jβ y Ex = − H 0η1 sin θ e 1 x 1 y e −σ x xη1 cosθ − jβ x− jβ y E y = H 0η1 cos θ e 1 x 1 y e −σ x xη1 cosθ
θi
kr 反射
ˆβ2 x + y ˆ β2 y 为区域 2 中的波矢,可推得波方程 β2 = x
β 2 × ( ε 2 −1 β 2 ) × H + ω 2 μ 2 H = 0
ki 入射
(7.37)
用矩阵算子表示叉积,波动方程可表示为矩阵形式
⎡ k 2 c − β 2 b −1 2y ⎢ 2 ⎢ −1 ⎢ β2 x β2 y b ⎢ ⎢ 0 ⎢ ⎣
β 2 = sx β1 = (1 − jσ x / ωε1 ) β1
x x
x
(7.48)
注意到其实部与 β1x 相等.联系(7.42),可知透射波的相速与入射角无关,即对任意角度来说都 是一样的.区域 2 中特征波阻抗和区域 1 中相等,这说明媒质是完美匹配的. 最后,代入(7.44)和(7.48)到(7.42)中,得到 TEz 入射波在 UPML 的场为
二维DGTD方法中UPML吸收边界的实现
二维DGTD方法中UPML吸收边界的实现李林茜;魏兵;杨谦;葛德彪【摘要】针对二维情形单轴各向异性完全匹配层吸收边界条件,研究了横磁波情形时域离散伽略金算法单轴各向异性完全匹配层吸收边界的理论基础和实现方案.借鉴时域有限差分方法——时域离散伽略金算法中吸收边界阻抗匹配、各向异性介质参数设置和匹配矩阵等思想,结合时域离散伽略金算法空间离散网格的非结构特性和离散单元之间场量传递的特点,给出了在单轴各向异性完全匹配层中电磁场量时域迭代公式,进一步离散成为离散时域迭代计算式.由于时域离散伽略金算法网格的非结构特性,一阶SM吸收边界条件一般对入射电磁波的反射率在-24 dB左右.仿真算例说明,给出的时域离散伽略金算法单轴各向异性完全匹配层吸收边界对电磁波的双重衰减达130dB,表明单轴各向异性完全匹配层具有良好的吸收效果.【期刊名称】《西安电子科技大学学报（自然科学版）》【年(卷),期】2016(043)006【总页数】5页(P86-90)【关键词】时域离散伽略金方法;单轴各向异性完全匹配层;横磁波【作者】李林茜;魏兵;杨谦;葛德彪【作者单位】西安电子科技大学物理与光电工程学院,陕西西安710071;西安电子科技大学信息感知协同创新中心,陕西西安710071;西安电子科技大学物理与光电工程学院,陕西西安710071;西安电子科技大学信息感知协同创新中心,陕西西安710071;西安电子科技大学物理与光电工程学院,陕西西安710071;西安电子科技大学信息感知协同创新中心,陕西西安710071;西安电子科技大学物理与光电工程学院,陕西西安710071;西安电子科技大学信息感知协同创新中心,陕西西安710071【正文语种】中文【中图分类】O441.4截断边界条件是许多电磁场数值方法精度保证的关键．在过去几十年，学着们提出许多种吸收边界条件，并成功应用于时域有限差分(Finite Difference Time Domain， FDTD)等方法中[1-2] ．文献[3]中提出了完全匹配层(Perfectly Matched Layer， PML)吸收边界，该边界是在计算域最外面增加一层非物理吸收层，该吸收层与相邻区域阻抗匹配，进入PML层的电磁波因为层内介质参数渐变而迅速衰减．后来有人提出的坐标伸缩完全匹配层(Coordinate stretched Perfectly Matched Layer， CPML)以及各向异性完全匹配层(Uniaxial anisotropy Perfectly Matched Layer， UPML)，均具有良好效果．时域离散伽略金(Discontinuous Galerkin Time Domain， DGTD)方法[4-5] 是基于时域体积元发展起来的一种兼备时域有限元(Finite Element Time Domian，FETD)方法网格剖分的灵活性和FDTD显式迭代特点的算法．文献[6]中将PML技术应用于DGTD，其后，文献[7]将PML技术应用于共形的三维DGTD计算．文献[8]提出一种DGTD提高PML稳定性的方案．文献[9-10]中讨论了高阶结点DGTD方法中的PML问题．上述文献主要讨论三维问题，然而在很多的情形，例如金属搭接缝中的电波传播问题、层状介质中的波传播问题、地面上方长传输线的电磁场分布、地下长线缆的电磁场分布问题和光纤里的光传输问题等都可以简化成二维问题来处理．考虑到二维方案在解决许多问题中方便快捷的特性，系统地研究DGTD算法的二维UPML吸收边界具有重要的意义．文中从各向异性介质中的Maxwell方程组和UPML基本理论出发，得到UPML介质层中的支配方程和辅助方程．考虑到DGTD空间离散网格的非结构特性，结合离散单元之间场量传递的特点，将方程加权积分得到单元矩阵方程．进一步，考虑相邻单元之间电磁场矢量切向分量突变的通量边界条件，通过基函数展开和矩阵单元计算，得到具有通量形式的时域公式，进而得到横磁(Transverse Magnetic，TM)波情形UPML中的时域迭代计算式．数值结果说明，文中UPML吸收边界对电磁波的衰减达 130 dB，吸收效果良好，满足二维情形电磁场精确仿真的需求．UPML是指在计算域外部增加一层非物理区域，该区域中为各向异性介质．实际操作中适当选择各向异性介质的本构参数，PML层波阻抗与相邻区域匹配，使入射到分界面上的电磁波无反射地进入该区域．在PML层中，电场和磁场满足的含辅助变量的支配方程为[11]其中，μ1和ε1为主要计算域中的介电系数和磁导率；a、b和c均为对角张量(具体表示为 a= diag(axx，ayy，azz)，其余类似)；辅助变量Ph和Pe与电场E和磁场H之间满足的微分方程为其中，d和κ同样为对角阵．UPML实现过程中，需要获得式(1)和式(2)的分量式并离散．对于TM波，关注的电磁场分量为Ez、Hx和Hy．因此，式(1)可转化为其中，辅助函数Ph和Pe的分量式为将式(4)代入式(3)，得TM波情形支配方程式(1)所满足的分量方程，即方程式(2)的分量式为分量方程式(5)和式(6)在具体实现数值计算时，需要改写为离散的时域递推公式．离散方案依赖于空间离散单元和基函数的选取等．在FDTD方法中常采用六面体离散计算域，并采用点匹配法，而DGTD与FETD类似常采用三维四面体或二维三角形离散单元和伽辽金法．由于采用非结构性网格离散，DGTD的UPML实现更为复杂．UPML区域划分如图1所示．在包含目标的计算域外部，吸收边界被划分成为4个棱边区和4个平面区．不同于FETD中单元交界面处采用电磁场连续的强制边界条件，DGTD采用数值通量的方式实现网格之间的场量交互．TM波情形在边界上交互时数值通量定义为[5]其中，m和m+分别表示当前元和相邻单元；和为通量系数；nx和ny分别表示单元中棱边的外法向单位矢量在x和y方向的投影．在式(7)通量定义的基础上，通过加权积分[12]和基函数展开，可得分量方程式(5)对应的矩阵方程，即其中，上标m为单元编号; Mm为单元质量矩阵；Sx和Sy为刚性矩阵；Fe，z、Fh，x和Fh，y分别指电场e和磁场h的通量变化量函数．将式(8)和式(9)在t=nΔt处离散，式(10)在t=(n+1/2) Δt处离散，可得其中，αgx、βgx、αgy、βgy、αgz和βgz为常系数．辅助方程式(6)无需加权积分，直接离散可得其中，ch、ci、cj和ck为常系数．式(11)～(16)就是二维TM波情形DGTD算法中UPML区域的递推计算公式．设二维计算域为8 m×8 m的矩形，区域离散为 10 902 个结点，共 21 474 个三角形单元，UPML层厚为 1 m．线源在计算域中心处．UPML层中电导率σ为非均匀分布，其渐变公式为其中，d为UPML区中单元的重心坐标位置；d0为UPML内侧界面位置；Δd为UPML层厚度；m为整数，本算例中 m=4．计算时，时间离散间隔Δt=0.16×10-10 s，当辐射电磁波的频率 f=0.3 GHz 时，文中方法的计算结果与利用Hankel函数得到的解析解的比较如图2所示．图2中三角形和实线分别表示文中算法和解析解的结果，由图可见，两种算法的计算结果一致，这说明了文中吸收边界算法的有效性．为更好地说明电磁波进入UPML层内的衰减情况，下面给出在辐射线源和计算域最外端连线上场值的分布情形．图3(a)是线元辐射的场值快照，线源位于A点处．线段AB上的电场如图3(b)所示，由图可见，电磁波在刚进入UPML层的阶段衰减不大(这是由于界面附近参数接近于内部计算域，便于电磁波无反射地进入吸收层)．随着电磁波深入UPML层，衰减量逐步加大(吸收层参数设置渐变所致)．吸收层界面上场值约为 -25 dB，而到吸收层外侧已下降到 -90 dB 以下．场值衰减为入射场的 1/500，考虑到返回计算域的电场值还需经过相同的路径，这一衰减的效果将加倍．文中从二维UPML的基本理论出发，根据空间离散网格的非结构特性，结合DGTD中离散单元之间场量传递的特点，在TM波情形下推导了DGTD各向异性完全匹配层吸收边界中的迭代计算式，实现了波在UPML层中的衰减．仿真算例说明，文中给出的DGTD UPML吸收边界对电磁波的双重衰减达 130 dB，吸收效果良好，满足二维情形电磁场精确仿真的需求．【相关文献】[1] 王林年, 褚庆昕, 梁昌洪. 交替方向隐式时域有限差分法中的Berenger理想匹配层[J]. 西安电子科技大学学报, 2006, 33(3): 458-461.WANG Linnian, CHU Qing xin, LIANG Changhong. Berenger’s Perfectly Matched Layer for the Alternating Direction Implicit Finite-difference Time-domain Method[J]. Journal of Xidian University, 2006, 33(3): 458-461.[2]姜彦南,葛德彪,杨利霞.卫星模型散射FDTD计算的共形边界研究[J]. 西安电子科技大学学报, 2007, 34(4): 587-589.JIANG Yannan, GE Debiao, YANG Lixia. Conformal Boundary Scheme in FDTD Computation of Satellite Model Scattering[J]. Journal of Xidian University, 2007, 34(4): 587-589.[3]BERENGER J P. A Perfectly Matched Layer for the Absorption of ElectromagneticWaves[J]. Journal of Computational Physics, 1994, 114(2): 185-200.[4]CHEN J, LIU Q H. Discontinuous Galerkin Time-domain Methods for Multiscale Electromagnetic Simulations: a Review[J]. Proceedings of the IEEE, 2013, 101(2): 242-254.[5]ALVAREZ J. A Discontinuous Galerkinfinite Element Method for the Time-domain Solution of Maxwell Equations[D]. Granada: University of Granada, 2013: 31-39.[6]XIAO T, LIU Q H. Three-dimensional Unstructured-grid Discontinuous Galerkin Method for Maxwell’s Equations with Wel l-posed Perfectly Matched Layer[J]. Microwave and Optical Technology Letters, 2005, 46(5): 459-463.[7]STYLIANOS D, LEE J F. Interior Penalty Discontinuous Galerkin Finite ElementMethod for the Time-dependent First Order Maxwell’s Equations[J]. IEEE Transa ctions on Antennas and Propagation, 2010, 58(12): 4085-4090.[8]PENG D, CHEN L, YIN W L, et al. An Efficient DGTD Implementation of the Uniaxial Perfectly Matched Layer[C]//Proceedings of the 2013 International Symposium on Antennas & Propagation: 2. Piscataway: IEEE, 2013: 1272-1275.[9]JI X, CAI W, ZHANG P. Reflection/Transmission Characteristics of a Discontinuous Galerkin Method for Maxwell’s Equations in Dispersive Inhomogeneous Media[J]. Journal of Computational Mathematics-international Edition, 2008, 26(3): 347-364. [10]NIEGEMANN J, KÖNIG M, STANNIGEL K, et al. Higher-order Time-domain Methods for the Analysis of Nano-photonic Systems [J]. Photonics and Nanostructures-Fundamentals and Applications, 2009, 7(1): 2-11.[11]GEDNEY S D, KRAMER T, LUO C, et al. The Discontinuous Galerkin Finite Element Time Domain Method (DGFETD)[C]//IEEE International Symposium on Electromagnetic Compatibility. Piscataway: IEEE, 2008: 1-4.[12]JIN J M. The Finite Element Method in Electromagnetics [M]. 3rd Edition. New York: John Wiley &Sons, 2014: 37.。
各向异性介质中二阶弹性波方程模拟PML吸收边界条件
瑞利面波数值模拟中的PML吸收边界条件
瑞利面波数值模拟中的PML吸收边界条件熊章强;唐圣松;张大洲【摘要】建立了弹性介质情况下完全匹配层(PML)吸收边界的2×12阶速度-应力交错网格有限差分算法,讨论了PML吸收边界条件的构建及其有限差分算法实现.通过与未加吸收边界及加常规指数衰减吸收边界3种情况下比较的波场模拟计算表明, PML吸收边界具有吸收更干净且能够吸收各种角度的边界反射等优点,其吸收率(吸收能量与未吸收能量之比)达到99.99%,很好地消除了周期折叠效应,使得所要计算的波场特征变得非常清晰,瑞利面波清楚地显示在波形记录上.【期刊名称】《物探与化探》【年(卷),期】2009(033)004【总页数】5页(P453-457)【关键词】瑞利面波;数值模拟;完全匹配层吸收边界;有限差分【作者】熊章强;唐圣松;张大洲【作者单位】中南大学,信息物理工程学院,湖南,长沙,410083;中南大学,信息物理工程学院,湖南,长沙,410083;中国地质大学,地球物理与空间信息学院,湖北,武汉,430074【正文语种】中文【中图分类】P631.4瑞利面波在地下复杂介质中的传播规律是近年来工程地球物理勘探中的重要研究课题［1－4］。
在利用计算机进行瑞利面波传播数值模拟计算过程中，总是对一个有限的空间区域进行计算，这样就自然引入了人为介质边界，从而使计算结果中不可避免地会产生一些假象，出现无实际意义的边界强反射，严重地干扰了对正常瑞利面波波场的认识。
假信号大都是在进行数值模拟计算时由四周边界的反射造成的，为了减弱或要消除这种假信号，就必须在人工边界上附加一定的边界条件，此条件应使向区域外传播的波不在边界上产生反射，因而它们一般被称为无反射边界条件或吸收边界条件。
为此，已发展了多种方法来削弱人为边界反射，一种好的边界条件不仅要吸收效果好，并且边界要能尽量贴近目标区域，这样可以节约计算机存储空间和计算时间。
计算吸收边界的方法有许多种，最早的吸收边界条件是由Lysmer和Kuhlemeyer 提出的［5］，其条件常被称为黏性边界条件，它的缺点是吸收人工反射波的效果较差。
PML边界条件下二维粘弹性介质波场模拟
的粘 弹性 体 ， 种粘 滞 性 表 现 为 基 质 内部存 在 “ 这 蠕
变 ” “ 弛” “ 、松 、 迟滞 ” 等现 象 。地震 波在 大地 中传 播
时， 质点 的振 动能 量 转 化 为 其 它各 种 形 式 的能 量 ，
目前 较 为普 遍 的说 法 即介 质 内部 的 热 效 应 现 象 。 这种 能量 的转 化使 地震 波在 粘弹 性介 质 中传播 时 ， 其 高频成 分 很 容 易 被 吸 收 ， 幅 近 似 指 数 规 律 衰 振
位 置发 生偏 离 。
的二 阶有 限差 分法 的一种 推 广 ［ 。此 外 ， ５ ］ 同有 限差 分法 和有 限元 法相 比 ， 伪谱 法不仅 克 服 了它们 对 高
频成 分 的限制 ， 以实 现 全 频 带 地震 波 模 拟 ， 且 可 而 节省 了前 两种 方 法所需 要 的大量 内存 。 我们 使 用伪谱 法模 拟粘 弹性 介 质 中的波 场 ， 并 引入 最 佳匹 配 层 （ ＭＬ） 理 边 界 ， 得 原 本 对 空 Ｐ 处 使 间做 的二维傅 里 叶变换 变成 了一维 傅里 叶变 换 ， 大 大 提高 了伪谱 法 的计 算效 率 。通过 模 型试算 ， 析 分 了粘 性 系数 对地 震波 的 吸收 和衰减 。
伪 谱 法 是 对空 间坐 标 通 过快 速 傅 里 叶 变换 在 时间 域作 差分运 算 ， 避免 了求 偏 导数 。由于伪谱 法 可 以看 成 是高 阶 网格 差 分 法 当其 空 间 差分 精 度 的
阶数 达到无 穷 时 的极 限 隋况 ， 故伪 谱 法可视 为 传统
减， 波形 和相位 失 真 ， 得分 辨率 降低 ， 使 反射 界 面 的
Ｐ Ｌ边 界 条 件 下 二 维 粘 弹 性 介 质 波 场 模 拟 Ｍ
多孔介质渗流传输特性数值模拟
多孔介质渗流传输特性数值模拟多孔介质是指由固体颗粒和孔隙组成的材料，其孔隙内充满了气体、液体或两者的混合物。
多孔介质在许多领域中发挥着重要作用，如地质工程、岩土力学、石油工程、环境科学等。
为了理解多孔介质中的渗流传输特性，数值模拟成为一种常用的工具。
数值模拟是通过建立多孔介质的物理模型和数学模型，运用计算机技术求解模型方程，从而获得多孔介质中渗流传输的各种参数和特性。
在数值模拟中，通常采用有限元法、有限差分法、边界元法等数值计算方法。
这些方法基于牛顿第二定律、达西定律、孔隙率、渗透率等物理规律，通过离散化和迭代求解，可以得到较为准确的渗流传输结果。
在进行多孔介质渗流传输特性的数值模拟时，首先需要建立二维或三维的几何模型。
几何模型可以根据实际多孔介质的形态进行构建，或者根据经验公式进行简化。
模型的精细程度对模拟结果的准确性有重要影响，因此需要根据研究目的和可用计算资源合理选择模型的细化程度。
接下来，需要确定多孔介质的物理性质参数，如孔隙度、孔径大小分布、渗透率等。
这些参数可以通过实验测量获得，也可以根据文献中的数据进行设定。
物理性质参数是决定多孔介质渗流传输特性的关键因素，因此选择合适的参数非常重要。
在模型建立和参数设定完成后，需要确定边界条件和初始条件。
边界条件包括入口流量和出口压力等，初始条件则指模拟开始时多孔介质内物理量的分布情况。
合理设定边界条件和初始条件可以更好地模拟多孔介质中的流体传输过程。
然后，通过数值计算方法，对模型进行离散化处理，并使用迭代算法求解模型方程。
在模拟过程中，需要考虑对流项、扩散项和源项等物理量的计算。
这些计算过程可以通过编程语言和计算软件实现，如MATLAB、Python、COMSOL等。
最后，根据模拟结果进行分析和评估。
分析包括流场分布、渗流速度、压力分布等多个方面。
这些结果可以帮助我们理解多孔介质中的渗流传输特性，指导实际工程的设计和优化。
评估模拟结果的准确性可以通过与实验数据的对比来进行，如果两者吻合较好，则说明模拟结果是可信的。
VTI介质交错网格FCT有限差分数值模拟
VTI介质交错网格FCT有限差分数值模拟张省;何兵寿;王玉凤【摘要】波动方程有限差分数值模拟是研究地震波在地下介质中的波场特征和传播机理的重要手段.对于常规有限差分技术,当采用大网格对计算空间进行差分离散时会出现严重的数值频散问题,降低了计算精度.通量校正(Flux- corrected transport method,FCT)技术能够有效压制粗网格情况下有限差分的数值频散.本文研究了具有垂直对称轴横向各向同性(Vertical Transverse Isotropy,VTI)介质的交错网格FCT有限差分技术.首先从一阶速度一应力弹性波方程出发,在交错网格空间中给出了该方程的高阶有限差分法格式及稳定性条件,在此基础上研究了波动方程正演过程中的数值频散FCT压制技术,二者结合实现了该方程的高精度有限差分数值模拟.同常规算法相比,本文算法不额外增加内存需求,少量增加计算量,但可有效压制VTI介质中弹性波动方程正演的数值频散现象.当采用大网格进行数值模拟时,本文方法明显提高了波场模拟精度.【期刊名称】《工程地球物理学报》【年(卷),期】2012(009)005【总页数】7页(P565-571)【关键词】VTI介质;数值模拟;FCT;有限差分;交错网格;数值频散;PML;Marmousi 【作者】张省;何兵寿;王玉凤【作者单位】中国海洋大学海底科学与探测技术教育部重点实验室,山东青岛266100;中国海洋大学海底科学与探测技术教育部重点实验室,山东青岛266100;中国海洋大学海底科学与探测技术教育部重点实验室,山东青岛266100【正文语种】中文【中图分类】P631.41 引言早期地震勘探的目标是寻找大的构造油气藏，都假设地球介质是完全弹性和各向同性的。
然而各向异性是地球介质中广泛存在的，如果采用基于各向同性介质的地震数据处理方法来处理各向异性介质的地震资料，就会导致地震分辨率降低和成像误差。
因此，开展地震各向异性研究是地震学家认知地球介质、进行地震学研究的发展趋势[1]，也是当今地震勘探与开发、地震预报、地球深部结构探测等研究领域中的重要课题之一。