周诚尧;汪勇;蔡伟祥;桂志先
【摘 要】首先将迎风机制引入五点八阶超紧致有限差分(CCD8)格式,得到五点七阶迎风超紧致(UCCD7)格式,并对两种格式进行数值频散分析和精度分析;其次建立了二阶声波方程的位移场时间四阶离散格式,将五点CCD8格式和五点UCCD7格式分别应用于位移场空间导数的求取,并分析这两种格式的稳定性条件;最后基于PML边界条件,将上述两种格式分别应用于声波方程的均匀介质、水平层状介质及Marmousi模型的数值模拟和波场特征分析及对比.研究结果表明:相较于普通紧致差分,五点CCD8格式具有小截断误差、高模拟精度、低数值频散、高稳定性、所需网格点数少的优点;引入迎风机制后,声波方程的五点UCCD7格式稳定性得到进一步提高;模型试算的结果验证了五点CCD8格式适用于复杂介质的数值模拟,模拟精度和计算效率都高. 【期刊名称】《石油物探》
【年(卷),期】2019(058)002
【总页数】11页(P176-186)
【关键词】五点八阶超紧致差分;五点七阶迎风超紧致差分;数值模拟;数值频散;稳定性条件
【作 者】周诚尧;汪勇;蔡伟祥;桂志先
【作者单位】油气资源与勘探技术教育部重点实验室(长江大学),湖北武汉430100;油气资源与勘探技术教育部重点实验室(长江大学),湖北武汉430100;中石化重庆涪陵页岩气勘探开发有限公司,重庆408014;油气资源与勘探技术教育部重点实验室(长江大学),湖北武汉430100
【正文语种】中 文
【中图分类】P631
地震正演模拟是探索地震波在介质中传播规律的重要手段,基于波动方程的数值模拟方法是一种地震正演模拟方法[1]。目前波动方程法主要包括伪谱法[2]、有限元法[3]、边界元法[4]、谱元法[5]以及有限差分法。有限元法和有限差分方法相比,二者精度相同时,有限元法耗时长[6]。有限差分法因其具有简单灵活、计算效率高以及占用内存小等优势被广泛应用 于地震波场数值模拟[7]。近年来，随着人们对数值模拟结果精度的要求不断提升，传统有限差分数值模拟方法的不足逐步凸显:一是数值计算时所需的网格点数多;二是具有严重的数值频散[8]。相较传统的差分格式,紧致差分格式在相同的计算网格时,精度和分辨率高,且稳定性高[9-10]。紧致差分格式的发展可以追溯到1989年DENNIS等针对Navier-Stokes方程首次提出了空间四阶紧致格式[11];1992年LELE对Pade格式进行了研究,提出了求解一阶导数和二阶导数的对称型紧致差分格式[12],该格式精度可达到谱方法的精度;1998年CHU等提出了三点六阶超紧致差分(combined compact difference 6,CCD6)格式用于对流扩散方程[13];林东等在三点CCD6格式的基础上提出了组合型超紧致差分格式,其精度最高可达到十阶[14]。在紧致差分格式中引入迎风机制可有效抑制非物理振荡,理论上可以提高紧致差分格式的稳定性。1994年傅德薰在紧致差分格式中引入迎风机制[15];1997年FU等又提出了更适合于多尺度复杂流场计算的五点五阶精度的迎风紧致格式[16];2002年王书强等在地震波场数值模拟时,首次将常规紧致差分格式应用于弹性波方程的数值模拟[17],截至目前，超紧致差分格式的地震波场数值模拟应用未见文献报道。
本文首先讨论了五点八阶超紧致差分(combined compact difference 8,CCD8)格式,然后通过引入迎风机制,在五点八阶超紧致差分格式的基础上,构造了五点七阶迎风超紧致差分(up
wind combined compact difference 7,UCCD7)格式,并将这两种格式应用于位移场空间导数的求取,实现了二维均匀介质模型、水平层状介质模型以及Marmousi模型的地震波场数值模拟。
zssi
1 超紧致差分格式及其迎风型方法原理
早期的紧致差分格式是基于Hermite多项式构造形成的,1992年LELE对Hermite公式进行了扩展,构造了紧致差分格式[12]:
(1)
式中:h为网格间距;a,b,c,β1,β2均为差分系数;fi为i节点的函数值;和分别表示i节点的一阶和二阶导数;fi+1,fi+2,fi+3,fi-1,fi-2,fi-3分别表示i节点依次向前3个节点和依次向后3个节点的函数值;分别表示i节点依次向前两个节点和依次向后两个节点的一阶导数;分别表示i节点依次向前两个节点和依次向后两个节点的二阶导数。
对上述差分格式进行泰勒展开并求解差分系数,可以得到一、二阶导数的普通八阶的紧致差分(compact difference 8,CD8)格式。
一阶导数普通八阶紧致差分格式:
(2)
二阶导数普通八阶紧致差分格式:
(3)
从公式(2)和公式(3)可知,普通紧致差分格式若要达到八阶精度,需用到7个节点。此外,紧致差分格式也适用于交错网格,其中八阶紧致交错网格差分(compact staggered difference 8,CSD8)格式可表示为:橙子剥皮机
(4)
式中:a1,b1,b2,b3是CSD8格式的差分系数;xi,Δx分别为x方向上的第i个节点、x方向上的空间步长。为了达到2n+2阶精度,只需2n个节点函数值。
超紧致有限差分格式在紧致差分格式基础上发展而来,1998年KRISHNAN提出了五点CCD8格式[18]:
(5)
将迎风机制引入五点CCD8格式,可得到如下的五点七阶迎风超紧致差分(UCCD7)格式,其精度可达七阶[19]:
(6)
根据超紧致差分的原理,将超紧致差分方法用于求解声波方程初值问题[20]。在完全弹性介质中声波方程可表示为:
(7)
式中:u(x,z)为二维介质的质点在x方向、z方向的位移;∂2u/∂t2为u(x,z)对时间t的二阶导数;∂2u/∂x2,∂2u/∂z2分别为u(x,z)在x方向和z方向的二阶导数;v为声波速度。
对(7)式采用四阶中心差分对时间导数离散可得:
(8)
智能点菜系统式中:un+1,un,un-1分别为n+1,n和n-1时刻位移;Δt为时间步长。
以z方向为例,说明利用五点CCD8格式求解公式(8)中u对z方向的空间偏导数的过程,x方向的求解过程不再赘述。假设U为位移场,其大小为m×n,则求∂u/∂z和∂2u/∂z2的过程可以用矩阵表示为AF=BU。其中A为公式(5)左端的差分系数矩阵,大小为2m×2m;B为公式(5)右端的差分系数矩阵,大小为2m×m;F为待求位移场空间一阶和二阶导数矩阵,大小为2m×n。矩阵分别为:
B=
(11)
根据AF=BU,可得F=A-1BU,其中F的奇数行矩阵元素为∂u/∂z,偶数行矩阵元素为∂2u/∂z2。同理,将位移场转置,可求得∂u/∂x和∂2u/∂x2。利用UCCD7格式求公式(8)中的空间偏导数的方法与此相同,仅有差分系数矩阵不同,故不再赘述。
以声波方程为例,说明了五点CCD8格式或五点UCCD7格式在地震声波波场数值模拟中的应用方法,该方法同样可以用于弹性波方程的数值模拟。各向同性完全弹性介质中的二维弹性波方程可以表示为:
(12)
式中:u(x,z)和w(x,z)分别表示x方向和z方向的质点位移;vP和vS分别表示纵、横波速度。u和w的时间四阶精度差分格式为[21]:
(13)
(14)
由公式(13)和公式(14)可以看出,弹性波方程与声波方程的差分格式相比,增加了许多二阶和四阶混合偏导数。对于这些混合偏导数,均可以按不同方向进行多次导数求取,求取顺序对结果无影响,求取方法不再赘述。
2 数值频散、稳定性及精度分析
水汽分离器
对波动方程的时间和空间偏导数的离散化造成了数值频散,它使相速度变成了网格间距的函数。数值模拟时,如果空间网格长度过大,会造成较大的求解误差,产生数值频散,降低模拟精度[22-23]。压制数值频散是采用有限差分方法时面临的关键问题之一[24]。频散关系分析
既是判断数值模拟所用方法优劣的重要依据,也是确定空间网格大小的重要依据[25],五点CCD8的二阶导数数值波数可以由公式(15)表示:
(15)
式中:h为空间网格大小;k为真波数;kx为x方向的真波数;θ为波的传播方向与x轴的夹角;为x方向的数值波数。
同理可得z方向的数值波数将和代入声波方程的差分格式(8),利用欧拉公式,化简可得频散关系:
(16)
式中:v′为数值波速;v是真波速;vΔt/h称为courant数。
数值波速与真波速的比值γ可表示为:
(17)
将简谐波解代入公式(5)得到如下公式:
(18)
为方便表达,引入中间变量数值求解公式(18)得到UCCD7的二阶导数数值波数,按照CCD8数值波速与真波速的比值求取方法可求得UCCD7的数值波速与真波度的比值。
将多种有限差分方法的速度比进行比较,结果如图1所示。图1a,图1b和图1c展示了不同的θ条件下4种常见网格差分格式的速度比曲线,其中α=vΔt/h;图1d中增加了八阶紧致交错网格格式的速度比曲线;图1e增加了五点七阶迎风超紧致差分格式的速度比曲线。考虑到参与对比的方法的稳定性,此处courant数均为0.25,横坐标φ∈[0,π],为波数与空间步长的乘积,纵坐标γ为速度比。速度比为1表示数值波速与理论波速一致,说明该方法能很好地压制数值频散,反之,则说明该方法的数值频散严重。
由图1可以看出,五点CCD8格式的速度比曲线偏离1时,其横坐标的数值最大,所以五点CCD8格式压制数值频散的能力强。图1e说明将五点CCD8格式引入迎风机制后的五点UCCD7格式,压制频散效果并无明显降低。五点UCCD7,五点CCD8,三点CCD6,CD6,CD8和CSD8格式在偏离1时所对应的横坐标数值分别为1.57,1.69,1.12,0.98,1.55和1.61,这6种方法在保证无数值频散时所需的最低网格点数分别为4.00,3.70,5.60,6.40,4.05和3.90个。五点CCD8格
式的数值频散的压制效果最好,CD6格式的数值频散的压制效果最差。
稳定性条件是影响差分方法计算效率的重要因素。本文使用Fourier方法[26-27]进行稳定性分析,略去推导过程,直接给出声波方程时间四阶精度的五点CCD8和五点UCCD7格式的稳定性条件αCCD8和αUCCD7分别为:
(19)
式中:vmax是指正演模拟时,所建立模型的最大声波速度,与五点CCD8格式相比,引入迎风机制的五点UCCD7格式的稳定性有了一定的提高,可使用略大的时间网格,在一定程度上减少计算时间,提高计算效率。
截断误差决定了有限差分格式和地震波场数值模拟精度,将文中所提到的几种方法的二阶导数差分格式截断误差进行对比。由Taylor级数展开[28]可以计算得到CD6、三点CCD6、CD8和五点CCD8格式的截断误差分别为-4.2163×10-4,-4.9603×10-5,2.7000×10-5和2.2046×10-6,即五点CCD8格式的截断误差最小,数值模拟精度最大。nand闪存
图1 不同差分方法的速度比曲线a θ=0,α=0.25时4种常见差分格式; b θ=π/6,α=0.25时4种常
重油燃烧器
见差分格式; c θ=π/3,α=0.25时4种常见差分格式; d θ=0,α=0.25时八阶紧致交错网格格式与4种常见差分格式; e θ=π/3,α=0.25时五点七阶迎风超紧致差分格式与4种常见差分格式
3 PML边界条件
数值模拟时,由于模型大小的限制,必然会存在人工边界,而人工边界的处理直接影响到数值模拟的精度及计算效率,若不对其进行处理则会干扰正常的地震波场。因此,本文在模型试算时,采用完全匹配层(perfectly matched layer,简称PML)对人工边界进行处理,其基本思想是:在吸收边界区域匹配与计算区域相同的波阻抗,使入射波无反射的进入吸收层,吸收层为衰减介质,可使得入射波迅速衰减直至消除[29]。理论上,PML边界能够吸收任何入射角度和频率的入射波,实践也证明PML吸收边界条件应用效果良好,可应用于声波和弹性波的研究[30-31],下式为声波方程的PML边界条件的控制方程:
豫公网安备41160202000603
站长QQ:729038198
关于我们
投诉建议