地屁处理2处理24非矩形不规则网格有限差分方法弹性波模拟孙卫涛杨慧珠(清华大学工程力学系)摘要本文提出了一种空间非矩形不规则网格有限Falk[Tess‘2〕等给出了交错]网格上的可变网格差分方法和R’5uu差分方法模拟具有弯曲界面和低速弯曲薄层的非mer[‘3Hestohlm4〕d[’用坐标变换方法均匀各向异性介质中弹性波传播过程该方法通过ft模拟曲线边界Mo[i]研究了三维模型复杂边界具有二阶时间精度和二阶空间精度的任意四边形不规则交错网格差分算子导出一阶弹性波动差分方程’6的处理方法oPrsa[l〕给出了非均匀介质中的二阶Ptark波动方程的矩形不规则网格差分方法ia[’7〕提在边界上应用吸收边界条件消除人工边界反出了各向同性介质中矩形不规则网格的有限差分方法射波该方法无需在精细网格和粗糙网格间插值本文给出了一种高阶精度非均匀各向异性波所有网格点上的计算在同一次空间迭代中完成在动方程非矩形不规则交错网格有限差分方法这种方法允许在弯曲几何界面附近和物质属性间断处使用任意四边形网格无需坐标变换和网格间的插值不规则曲线界面附近使用非矩形不规则网格能够精确的模拟模型几何构造;在几何尺度变化很大的模型中无孺网格加密就可以达到很高的计算精度简单易行而且占用内存少计算量小适合于复杂介理论分析和数值算例表明该方法十分适合处理抵质模型的弹性波正演问题文中推导出的有限差分糙界面断层和弯曲薄层等复杂几何构造不但节省了大t内存和计算时间而且具有令人满意的稳定性和精度更高方程具有一般形式其退化形式与规则网格方法完全吻合本文给出了该差分格式的数值频散分析在模拟复杂地层构造中地震波传播时并导出了差分稳定性条件数值算例表明这种新该方法比规则网格和矩形不规则网格差分方法效率方法具有很高的计算效率和精度己l,l绪侣「J理论公式笛卡儿坐标系斜方晶系各向异性介质弹性波动有限差分方法是求解波动方程的常用数值方法之一地震波正演模拟中有限差分方法的较早研究方程速度应力公式为PV=:包括Aly[4temranord[’]和Karal川玫幻re[2〕Alf[’6KelDT(1)(2)一]和viiSv波模拟reux]等学者的著作Levander[7〕在Pt二CDT.Vv中引入了四阶空间差分算子Dblaina[“〕其中应力矢量pT(x)为介质密度质点速度矢量(x(xt)aev[提出了高阶差分算子的方法Grs9t〕给出了等效普)微分算子矩阵D和斜方晶系介质:物质参数的三维四阶速度应力有限差分方法弹性系数矩阵VTTC=二(9个独立系数)定义如下v二通差分方法基于笛卡儿坐标系中的规则正方形或矩形网格在模拟曲线界面时出现阶梯状边界必须v,v(3)(4)T(Tx采用精细网格才能减少虚假绕射波致计算量的大大增加shor另外局部物理参数的剧烈变化也会要求加密整个模型网格导D二aeeLapl(5)ey[’0]tl首先研究了方程中的不规essasra则网格有限差分方法Jtm和Tmer[“]处理2c124e1300Y。(a一=月一imjo(yn)为c22c23c33o(y为1)一1)22o(y一湘一1()xx一Px;)11000一(凡一为)儿)%)o(y一为)o(y为)(一xC=(6)刃(湘1一(为(湘,一儿),儿)2(为(为,一儿)(x%)(x,,,一x))l(为2一2一2一2一x1在交错网格上离散一阶偏微分算子笛卡儿坐标系x聊metric(15)y平面内二维非矩形不规则网格见图l,一。二一日「全!川2孟,士二‘,,刁价二丁刁广不厂川以不泛异口X小刀一:癸全十二、‘X,一x,2;y1…夕十2。)价(16)这里N言是8个相关节点(m2一lmm十lmn;一1nn十1n+)空间坐标的函数可以通29)的矩阵求逆得到过公式(,。价是空间节点q处波价~刁协得和举场值同理容到“易“一举a夕az~~~~式(9)可退化为矩形不规则网格公式,目~*。二。a:图1非矩形不规则网格有限差分示意图{{!一Dxx,0必Y沉节点(mi汁和(m十n1)l的中点命名为节点]DY1\\LO(17):节点(m..二n)和(m十)的中点命名为节点l:交错网格一阶时间导数差分算子定义为Ja价~相关节点处的波场值小可以写成级数形式==J价1一笋jatxx二nx++v,vnv++o(△3)o(留)3(7)(8)一△t(18)nxDv数值频散和稳定性考虑平面波e略去空间步长高阶小量O(△)后空间导数可以表示为价的线性组合:(‘i“)r传播方向与xy:s坐ly岭):(熟m(价XX二YmY‘of,0(9)标轴成=二k/yly:ey37:夹角=,k/这些角度可以表示成co}kJso!k}和二soc,)3=k/Ik{其中k二k,和k是波数(xk:(kP二kk)的分量空间位:1,二价,。价1,,笋2)T(10))T移矢量rT=y一)波频散关系式.=(尹m1,价二‘笋二1‘笋m+2}李}其中‘三十‘互十‘已。(19)。(11)了为P波速度’,是数值角频率。由交错网格:DX=刁价Jxa笋a夕22刁价ZaZ价2!Jx2!日x刁夕tn笋j汀|订l产时间中心差分算子得到数值角频率=一的表达式(12)(13)量aZDY=x,笋2刁尹!刁少2n2!刁夕Jxm必,(20)\\2一由非矩形不规则网格差分算子得到数值波数分△t一“鱼。;n{竺。(a(xo月一=mj)二却)一…一今l)即)(xo一x(2一却一1(%)o(y)一Py)l1川门Jw卜we片l(几一(几(今(今一今)2(xo(,(却lse孟代花y工一名(N言qlnisl(+‘。一、k二,一动一+l夕。一夕lk,}z。一221k))(21)一(,1一几)1一几)几)l一x)(为,%)儿)=一(,2一几)2一22一几)(为2一习(叫qlsni(lx+。一x,lkxlk+}夕。一夕lk,12。一z,))(22)处理2zk4一+烈}为n‘,X‘代一i一少,一x二,xkyly:y3:都小于1得到P波稳定性条件x.2(*1一x二)卜k,+}z。一z‘卜k))(23)一引进无量纲量氛Hi(i二xy)后得到了与z二又矩形不规则网格相似的频散公式y,几和又闪局一“汤“a是为Lev咒汤,a欠(31)z坐标方向波长a二a,口是xyz坐标方z当网格为正方形时这个稳定性条件可以退化7d[near]给出的公式;采用二阶差分时该稳定向P波速度下标xyz分别表示沿xy坐标方向的量乞二乞二ieux性条件退化为vir[61的公式△ta一\"、汤t+~跳产_工l一毛毯(夕一少△t、1一夕)算例口乞l+△t1一1一yl(24)y考察弯曲薄地层海底模型检验非矩形不规则网格有限差分方法的有效性200mxH从=(x二,x二+1一x孟x模型尺寸Zomx=(z,十,一z,)在深度80一140m范围内有一个弯曲薄地层又‘(25)P交界面(见图2)界面层上方为海水界面层下方为将氛和H代人数值频散关系式得到波数高速地层步长为高斯爆炸源位于(IOOm6Om)震源中l值速度与真实速度之比称、,一心频率为40H:000025模型密度和弹性常数见表计算时间为035时间众51一}f1使用8]Higdon[’公乞一吸收边界条件消除边界反射分别采用规则网格(26)同理得到S波数值速度与真实速度之比凡为x矩形不规则网格和本文的非矩形不规则网格差分方法计算弹性波传播非矩形不规则网格划分示意图见图3q方向上axS1波速度。,1,2222距离0(100)20凡妥尹毖二Slnl=屯二IAa口C}生习一J+a刁〕一任‘,拼A一吠H吮y+-aa于‘,拼AJlO扣~.~曰‘~山~‘~~.~‘~一泛H吮二(27)其中A二:.0一十}!y郭I(之分釜…、snixH,oc一21转侧司加~+1一y}lvclIso,2+z之一21l一l+从zcHso82月\\功们,一一一一一一曰扛图2A,一名!衅sin·;是是一六笼去…卜…、。(…老弋…一…、弯曲薄地层海底模型衰1育曲薄地层海底徽型介质密度和弹性常教8O计算参数介质I1000a介质n26506介质m221112e10(29)Apcl;c12(掩/me22cl3c33c233一(象为y+N;·‘·训X二几之勿+1xq一x份一Xm,…zxHC叱“(P)(P)a225函1因225e9039e50e9e。。‘55。66一yl一y…、。。s,2(Pa)029e9359e一21一乞l之l+1sHOczy非矩形不规则网格方法能够根据模型几何构造的形状划分网格具有很大的自由度克服了传统差(30)令P波数值频散中的si-n‘项在任意人射角分方法均匀正方形网格不适应复杂几何模型的弱处理点24同时在曲线边界附近使用变形的网格无需网“”削弱了虚假绕射(图5);本文非矩形不规则网格方格加密就可以消除阶梯状界面产生的非物理绕射波与矩形不规则网格方法相比节省了计算量提高了计算精度法用变形的网格模拟曲线界面在未增加网格密度的情况下很好的消除了虚假绕射波(图6)如果在相同网格规模(75x75)下规则网取得相同的计算精度规则网格方法要在整个模型中采用精细网格网格规模将达到(300格矩形不规则网格和非矩形不规则网格差分方法x300);矩形不规则网格方法在界面附近加密网格网格规模达50到(1X220)两种方法的计算量急剧增加具有几乎相同精度的非矩形不规则网格有限差分方法明显降低了计算时间(见图7)和计算机内存使用量但是网格分布的对称性会受到一定影响图3弯曲薄地层海底模型非矩形不规则网格有限差分网格划分图6非矩形不规则网格有限差分水平应5力分量、在0102时刻波场快照图4规则网格有限差分水平应力分量5时刻波场快照等值线、在0102图7相同计算精度规则网格矩形不规则网格和非矩形不规则网格有限差分计算时间对比结论本文提出一种空间非矩形不规则网格有限差分图5矩形不规则网格有限差分水平应5力分t、在0102时刻波场快照方法模拟地震波在复杂几何构造各向异性弹性介质中的传播该方法无需坐标变换或者粗细网格之间的插值比规则网格和矩形不规则网格有限差分模拟弹性波场快照见图4一图6图4中规则网格方法占用内存少大大降低了计算量;适合于弹性常方法在薄地层界面附近产生虚假绕射;矩形不规则网格方法在弯曲界面附近加密网格在一定程度上数密度以及几何构造剧烈变化的大规模数值模拟理论分析和数值算例表明该方法具有良好的稳定处理24性和精度同时这种方法具有网格划分灵活的特9GravseRWSimulati吃seismiewaveProPagationin3De差lastiemediausingstaggerdegridfinitedifferene点为自适应分方法提供了良好基础seBulletinfotheeSismolgoiealS〕eietyofAmeriea199686:1091一今考文献110610eShortlyGHandWellerRNumeriealoslutionofuai月terman2KaralFCJrPropagationofelastiewapvseaLlaees叫tnoJPhys19389:334一348astraerPAplinlayeredmediabyfinitediffereneemeth5u11mCand浏BlletinfoJTEElsatiem浏elingnoagrideSismoeaer:witivarthegO1ilSx二ietyofAmiea196858(l)367一hverteallymseyingsgGeoPhy哪cinsP哪P199434982:357一370a3刀reMrimiewFlkJTeranevsTubewavem记eDFinitedifereneemethl记5ofsesave12propanbesmEdGajkiDingythefinite-differenees,tioninheterogenoeusmateiralsIn:Meth记iBoltmeth司withva斗inggridPsaeBAEdC冶mputatinoalph邓icscAadmeiePrseIneingPageoPh199614:77一93197213TsemerEKoslfofDandBehleAElastiewaveProPagationsiRx〕ermulat月肠记MKellyKRandBDMynoiinthepneeaeetopographyGoeat:finitedifferneeemdoelingoftheascAeuarefo因uatinoPhySJIntern1992serusd108621一63214HestGeOPhaysc197439:834一842uoctiewaveholm50andRuudB02DFinitediffereneeetavesuKelliyKRWardRWadnTreitelSetalS”thetielsaiewm司elingineludingfraeetoopgraPhyGeosepismogarmsafinitediftereneeapPeG印physicsh邓Psorp199442:371一390197641:2一27aorhMufti151RaLgrecsalethree-dimensionalseismiemdoelsantiinteretive1sesVirieuxJSHwaveProPagationinheterogenoeusmedia:dher甲5,ificaneeG印phyi199055:vecolitystrsefinitediffereneemethdoGeophysies1984(9)1166一1182an49:1933一194216IvoOPrasldJiriaZhradnikElastiefinitediffereneeetrirreVirieuxJPSVwaveProP昭atinoinheterogensemh记fo,largirdsG印physics199964(1):avecolitystrsefiniteets4differeneemh记Guoemdi印phyics20一250198651:889一90117PitarkaA3Delsatiefinit-ediffereneemdeelingofseimsie路er记swithunb四g朋der月anRFb以tr卜以derfiint已diflerenee-PSV妇sromtionsuinsta幻g乎别比曰(尧扣h邓isc198853:1425一1436BulleSigmshgdS二(Anl199989(onifo1):54一mrsPaein68aDlbainMATheaPPlieatinoofhighorderdiftereneing18HidonRLAb沁rbing比unda职seneeaPPxtototheacslarwaveuationGeophysies198651:54一roimatoisnhemutltidimnocditionofrdiferensi叫noalwave闪uat66ionMathC冶mP198647:437一459161