地屁处理2处理
24
非矩形不规则网格有限差分方法弹性波模拟孙卫涛杨慧珠
(清华大学工程力学系)
摘要本文提出了一种空间非矩形不规则网格有限差分方法模拟具有弯曲界面和低速弯曲薄层的非均匀各向异性介质中弹性波传播过程该方法通过具有二阶时间精度和二阶空间精度的任意四边形不规则交错网格差分算子导出一阶弹性波动差分方程在边界上应用吸收边界条件消除人工边界反射波该方法无需在精细网格和粗糙网格间插值所有网格点上的计算在同一次空间迭代中完成在不规则曲线界面附近使用非矩形不规则网格能够精确的模拟模型几何构造;在几何尺度变化很大的模型中无孺网格加密就可以达到很高的计算精度理论分析和数值算例表明该方法十分适合处理抵糙界面断层和弯曲薄层等复杂几何构造不但节省了大t内存和计算时间而且具有令人满意的稳定性和精度在模拟复杂地层构造中地震波传播时该方法比规则网格和矩形不规则网格差分方法效率更高Falk[‘2〕等给出了交错网格上的可变网格差分方法Tessmer[‘3]Hestholm和Ruud[
’4〕用坐标变换方法
模拟曲线边界Mofti[’5]研究了三维模型复杂
边界
的处理方法oPrsal[’6〕给出了非均匀介质中的二阶
波动方程的矩形不规则网格差分方法Pitarka[’7〕提出了各向同性介质中矩形不规则网格的有限差分方法本文给出了一种高阶精度非均匀各向异性波
动方程非矩形不规则交错网格有限差分方法这种方法允许在弯曲几何界面附近和物质属性间断处使
用任意四边形网格无需坐标变换和网格间的插值
简单易行而且占用内存少计算量小适合于复杂介质模型的弹性波正演问题文中推导出的有限差分
方程具有一般形式其退化形式与规则网格方法完全吻合本文给出了该差分格式的数值频散分析并导出了差分稳定性条件数值算例表明这种新方法具有很高的计算效率和精度
己l侣绪,l「J
理论公式
有限差分方法是求解波动方程的常用数值方法之一地震波正演模拟中有限差分方法的较早研究包括Alterman和Karal川玫幻re[2
〕Alford[’]Kel
一y[4]和virieux[’6]等学者的著作Levander[7〕在
P
Sv波模拟中引入了四阶空间差分算子D
ablain[
“
〕
提出了高阶差分算子的方法Graves[9〕给出了等效物质参数的三维四阶速度应力有限差分方法普通差分方法基于笛卡儿坐标系中的规则正方形或矩形网格在模拟曲线界面时出现阶梯状边界必须
采用精细网格才能减少虚假绕射波另外局部物理参数的剧烈变化也会要求加密整个模型网格导致计算量的大大增加shortley[’0]首先研究了Laplaee
方程中的不规
则网格有限差分方法Jastram和Tessmer[“]
笛卡儿坐标系斜方晶系各向异性介质弹性波动方程速度应力公式为:
PV=DTt二CDT.V
(1)(2)
其中p(x)为介质密度质点速度矢量v(xt)应力矢量T(xt)微分算子矩阵D和斜方晶系介质
弹性系数矩阵C(9个独立系数)定义如下:
VT=
v二v,v
TT
二
(
Tx
D二
(3)(4)
(5)处理24c12e1300
Y。(a月=imjn
)
c22c23
聊C=c3300
0(6)
metr
ic
一(yo一为1)(yo一为1)2(yo一湘1)(x一xP;)1
一(凡一为)(yo一为)2(yo一为)(
x
一x
刃
1
(湘1一儿)(为,一儿),(为,
一儿
)(x
,,一x
)l
(为2一%)(湘2一儿)2(为2一%)(x,2一x)1
(15
)
在交错网格上离散一阶偏微分算子笛卡儿坐标系xy平面内二维非矩形不规则网格见图l,一。二,,刁价二士二‘
一日「全!川2孟,丁异刁广不厂一川以不泛小
刀:
口X
癸二全、‘X,一
x,
2
;y
1
…夕十
2)价
。
(16)这里N言是8个相关节点
(m一lmm十lm
十2;n一1nn十1n+2)
空间坐标的函数可以通
过公式(9)的矩阵求逆得到价。是空间节点q处波
,,目~*。二。a价~刁协
场值同理容易得到举和举
~~“~
“
~
一a夕az
式(9)可退化为矩形不规则网格公式:
图1非矩形不规则网格有限差分示意
图
{Dx{一
!x
,
\DY1LO(17)沉必0节点(m汁和(m十1l)的中点命名为节点i节点(mn)和(mn十l)的中点命名为节点]相关节点处的波场值小可以写成级数形式:.二=x二nx+v,nv+o(△3)(7).=xnx+vDv+o(留)(8)略去空间步长高阶小量O(△3)后空间导数可以表示为价的线性组合:Y
交错网格一阶时间导数差分算子定义为:a价J~价J1一笋
j
at
一△
t
(18)
数值频散和稳定性
X二
Ym
XY‘fo,
0(9)
(价m1,价二,价。1,笋,2
)T(10))T(11)(12)(13)考虑平面波ei(‘“r)传播方向与xy:坐标轴成yly:y3夹角这些角度可以表示成cosyl=k二/}kJeos7:=k,/!k}和cos)3=k/Ik{其中k二k,和k是波数k二(k二k,k)的分量空间位移矢量rT=(xy:)P
波频散关系式:
(:)熟岭
.=
(尹m1,价二‘笋二1‘笋
m+
2
}李}一‘三十‘互十‘
已
其中。为P波速度。是数值角频率(19)由交错网格
时间中心差分算子得到数值角频率。的表达式:
DX=
,’=鱼。;n{竺。
一△t一“\2一(20)
DY=
笋ntj必m,
由非矩形不规则网格差分算子得到数值波数分
x,(
a
月=
刁价刁2价aZ价Jx2!JxZ2!日x刁夕a笋aZ笋刁2尹a夕2!刁少22!刁夕Jxmjn)二量汀|产
订l了
…一(xo一今l)(xo一却)(x一却1)(%一yPl)一(几一即)(几一今)2(xo一,)(yo一动一(,1一几)(今1一几)2(,l一x)(为,一%)一(,2一几)(今2一几)2(却2一几)(为2一儿)一名(N言sin(l‘。一、ql+l夕。一夕lk,+}z。一221k二k))(21)=习(叫ql+}夕。一夕lsin(lx。一x,lkxk,+12。一z,lk))工y孟代花(22)门lleswJwe川卜1片处理24
kz一烈‘代一in‘,X,一二,kx+}为一少卜k,+}z。一z‘
卜k))(23)
引进无量纲量氛Hi(i二xyz)后得到了与
矩形不规则网格相似的频散公式又二几,和又是xyz坐标方向波长a二a,口是xyz
坐标方
向P波速度下标xyz分别表示沿xyz坐标方向的量yly:y3都小于1得到P波稳定性条件:2(x.*1一x二)一闪局一“汤“a咒汤,a欠(31)当网格为正方形时这个稳定性条件可以退化为Levander[7]给出的公式;采用二阶差分时该稳定性条件退化为virieux[61的公式△ta一_~、汤
"t+l一工跳产
算例少一毯毛乞二乞二口△t乞l+1一yl△t(夕、1
一夕)
(24
)
Hx=
从=
(x
二,1一x二
孟x(z,十,一z,)又
y+
1
一x
将氛和H‘代人数值频散关系式得到(25)P波数值速度与真实速度之
比称
乞一公
、,
一
众
51
一}
(26
)
同理得到S波数值速度与真实速度之比q
凡为x方向上S波速度ax凡二妥尹
毖
f。,1
2222
11口C,a刁〕aJl
Slnl=屯二IA+一任‘,拼A+-
于‘,拼A
}a
生习一J吠H吮一ya泛H吮一
二
考察弯曲薄地层海底模型检验非矩形不规则网格有限差分方法的有效性模型尺寸Zom
x
200m在深度80一140
m
范围内有一个弯曲薄地层
交界面(见图2)界面层上方为海水界面层下方为
高速地层高斯爆炸源位于(IOOm6Om)震源中
心频率为40H:模型密度和弹性常数见表l时间
步长为000025计算时间为035使用Higdon[
’8]
吸收边界条件消除边界反射分别采用规则网格矩形不规则网格和本文的非矩形不规则网格差分方法计算弹性波传播非矩形不规则网格划分示意图
见图3距离
()
010020
O扣~.~曰‘~山~‘~~.~‘~
(27)
加司侧
转
其中
:
A二一
郭、sin(I
之分釜
…
Hxco
一
.0
十}~}Ilvcos
,
2+
!y+1一ylz,一
21
之l+l一21
A,一
名!衅
sin
(…老弋……、一
…、。O8一卜笼是是;,扛一一一一一一曰图2弯曲薄地层海底模型衰1育曲薄地层海底徽型介质密度和弹性常教·…六去
A一象
(N;
·‘
·
(29)计算参数p(掩/m3功们28\月cosHz从xq一x份X勿+1一Xm训二之,几…、。。s,2…HxC叱“cl;e22c33(Pa)c12cl3c23(Pa)e。。‘55。66(Pa)介质I
1000225函225e90介质n265061因03e929e9介质m221112e1050e935e9yOs
Hzc
为一y
y+l一yz
一
21
之l+1一乞l(30)
令P波数值频散中的sin-‘项在任意人射角非矩形不规则网格方法能够根据模型几何构造的形状划分网格具有很大的自由度克服了传统差分方法均匀正方形网格不适应复杂几何模
型
的弱